diff --git a/src/physics/thin_wall.F90 b/src/physics/thin_wall.F90 index 70318b9..1f21b68 100644 --- a/src/physics/thin_wall.F90 +++ b/src/physics/thin_wall.F90 @@ -171,6 +171,7 @@ MODULE thin_wall REAL(r8) :: quad_tols(3) = [0.75d0, 0.95d0, 0.995d0] !< Distance tolerances for quadrature order selection INTEGER(i4) :: quad_orders(3) = [18, 10, 6] !< Quadrature order for each tolerance REAL(r8), PARAMETER :: target_err = 1.d-8 +REAL(r8), PARAMETER :: coil_min_rad = 1.d-6 integer(i4), public, parameter :: tw_idx_ver=1 !< File version for array indexing character(LEN=16), public, parameter :: tw_idx_path="ThinCurr_Version" !< HDF5 field name CONTAINS @@ -205,10 +206,17 @@ SUBROUTINE tw_setup(self,hole_ns) DO i=1,self%n_vcoils IF(ANY(self%vcoils(i)%res_per_len<0.d0))CALL oft_abort("Invalid resistivity for passive coil", & "tw_setup", __FILE__) + IF(ANY(self%vcoils(i)%radius Compute coupling from thin-wall model elements to flux loop sensors +!! +!! @note The asymptotic form of self-inductance for thin circular coils derived +!! by Hurwitz and Landreman [arXiv:2310.09313 (2023)] is used for all inductance +!! calculations. Note that this is not strictly valid for mutual inductances, +!! but is instead used to avoid integration challenges with very-closely-spaced coils. !------------------------------------------------------------------------------ SUBROUTINE tw_compute_Lmat_coils(tw_obj) TYPE(tw_type), INTENT(inout) :: tw_obj !< Thin-wall model object LOGICAL :: exists INTEGER(4) :: i,ii,j,jj,k,kk,l,ik,jk,ih,ihp,ihc,file_counts(4),ierr,io_unit,iquad REAL(8) :: tmp,dl,pt_i(3),pt_j(3),evec_i(3,3),evec_j(3,3),pts_i(3,3),pt_i_last(3) -REAL(8) :: rvec_i(3),r1,rmag,rvec_j(3),z1,cvec(3),cpt(3),coil_thickness +REAL(8) :: rvec_i(3),r1,rmag,rvec_j(3),z1,cvec(3),cpt(3),coil_thickness,sqrt_e REAL(8) :: rgop(3,3),norm_i(3),area_i,f(3),dl_min,dl_max,pot_last,pot_tmp REAL(8), allocatable :: atmp(:,:),Acoil2sen_tmp(:,:),Acoil2coil_tmp(:,:) CLASS(oft_bmesh), POINTER :: bmesh @@ -731,6 +744,7 @@ SUBROUTINE tw_compute_Lmat_coils(tw_obj) Acoil2coil_tmp=0.d0 WRITE(*,*)'Building coil<->coil inductance matrix' IF(tw_obj%n_vcoils>0.AND.ncoils_tot>0)THEN + sqrt_e=SQRT(EXP(1.d0)) !$omp parallel private(i,ii,j,k,kk,tmp,pt_i,rvec_i,atmp,cvec,cpt,coil_thickness) ALLOCATE(atmp(ncoils_tot,1)) !$omp do @@ -742,12 +756,12 @@ SUBROUTINE tw_compute_Lmat_coils(tw_obj) pt_i = (coils_tot(l)%coils(i)%pts(:,ii)+coils_tot(l)%coils(i)%pts(:,ii-1))/2.d0 DO j=1,ncoils_tot DO k=1,coils_tot(j)%ncoils - coil_thickness=MAX(coils_tot(l)%radius(i),coils_tot(j)%radius(k)) + coil_thickness=(MAX(coils_tot(l)%radius(i),coils_tot(j)%radius(k))**2)/sqrt_e tmp=0.d0 DO kk=2,coils_tot(j)%coils(k)%npts cvec = coils_tot(j)%coils(k)%pts(:,kk)-coils_tot(j)%coils(k)%pts(:,kk-1) cpt = (coils_tot(j)%coils(k)%pts(:,kk)+coils_tot(j)%coils(k)%pts(:,kk-1))/2.d0 - tmp = tmp + DOT_PRODUCT(rvec_i,cvec)/MAX(coil_thickness,SQRT(SUM((pt_i-cpt)**2))) + tmp = tmp + DOT_PRODUCT(rvec_i,cvec)/SQRT(SUM((pt_i-cpt)**2) + coil_thickness) END DO atmp(j,1)=atmp(j,1)+coils_tot(j)%scales(k)*tmp END DO diff --git a/src/tests/physics/test_ThinCurr.py b/src/tests/physics/test_ThinCurr.py index fdda64e..76e308a 100644 --- a/src/tests/physics/test_ThinCurr.py +++ b/src/tests/physics/test_ThinCurr.py @@ -410,7 +410,7 @@ def test_td_plate(direct_flag,python): @pytest.mark.parametrize("direct_flag", ('F', 'T')) @pytest.mark.parametrize("python", (False, True)) def test_td_plate_volt(direct_flag,python): - sigs_final = (4.E-3, 4.766694E-4, 4.010512E-4) + sigs_final = (4.E-3, 4.580643E-4, 3.854292E-4) assert ThinCurr_setup("tw_test-plate.h5",1,direct_flag, vcoils=((0.5, 0.1),), floops=((0.5, -0.05), (0.5, -0.1)), @@ -432,7 +432,7 @@ def test_td_cyl(direct_flag,python): @pytest.mark.parametrize("direct_flag", ('F', 'T')) @pytest.mark.parametrize("python", (False, True)) def test_td_cyl_volt(direct_flag,python): - sigs_final = (4.E-3, 1.710824E-4, 1.453702E-4) + sigs_final = (4.E-3, 1.504279E-4, 1.276624E-4) assert ThinCurr_setup("tw_test-cyl.h5",1,direct_flag, vcoils=((1.1, 0.25), (1.1, -0.25)), floops=((0.9, 0.5), (0.9, 0.0)), @@ -457,7 +457,7 @@ def test_td_torus(direct_flag,python): @pytest.mark.parametrize("direct_flag", ('F', 'T')) @pytest.mark.parametrize("python", (False, True)) def test_td_torus_volt(direct_flag,python): - sigs_final = (4.E-3, 6.501287E-5, 4.640848E-6) + sigs_final = (4.E-3, 5.653338E-5, 4.035387E-6) assert ThinCurr_setup("tw_test-torus.h5",1,direct_flag, vcoils=((1.5, 0.5), (1.5, -0.5)), floops=((1.4, 0.0), (0.6, 0.0)), @@ -470,7 +470,7 @@ def test_td_torus_volt(direct_flag,python): @pytest.mark.parametrize("direct_flag", ('F', 'T')) @pytest.mark.parametrize("python", (False, True)) def test_td_passive(direct_flag,python): - sigs_final = (4.E-3, 7.706778E-4, 7.903190E-4) + sigs_final = (4.E-3, 8.349309E-4, 8.364054E-4) assert ThinCurr_setup("tw_test-passive.h5",1,direct_flag,eta=1.E4, icoils=((0.5, 0.1),), vcoils=((0.5, 0.0),), @@ -483,7 +483,7 @@ def test_td_passive(direct_flag,python): @pytest.mark.parametrize("direct_flag", ('F', 'T')) @pytest.mark.parametrize("python", (False, True)) def test_td_passive_volt(direct_flag,python): - sigs_final = (4.E-3, 4.228515E-4, 4.338835E-4) + sigs_final = (4.E-3, 4.379235E-4, 4.389248E-4) assert ThinCurr_setup("tw_test-passive.h5",1,direct_flag,eta=1.E4, vcoils=((0.5, 0.0), (0.5, 0.1)), floops=((0.5, -0.05), (0.5, -0.1)), @@ -519,7 +519,7 @@ def test_eig_torus(direct_flag,python): @pytest.mark.parametrize("direct_flag", ('F', 'T')) @pytest.mark.parametrize("python", (False, True)) def test_eig_passive(direct_flag,python): - eigs = (1.483589E-1, 6.207849E-2, 2.942791E-2, 2.693574E-2) + eigs = (1.504155E-1, 6.423383E-2, 3.190175E-2, 2.942398E-2) assert ThinCurr_setup("tw_test-passive.h5",2,direct_flag,eta=1.E4, vcoils=((0.5, 0.1), (0.5, 0.05), (0.5, -0.05), (0.5, -0.1)),python=python) @@ -549,7 +549,7 @@ def test_mred_torus(direct_flag): @pytest.mark.coverage @pytest.mark.parametrize("direct_flag", ('F', 'T')) def test_mred_passive(direct_flag): - eigs = (1.483589E-1, 6.207849E-2, 2.942791E-2, 2.693574E-2) + eigs = (1.504155E-1, 6.423383E-2, 3.190175E-2, 2.942398E-2) assert ThinCurr_setup("tw_test-passive.h5",4,direct_flag,eta=1.E4, vcoils=((0.5, 0.1), (0.5, 0.05), (0.5, -0.05), (0.5, -0.1))) @@ -595,8 +595,8 @@ def test_fr_torus(direct_flag,python): @pytest.mark.parametrize("direct_flag", ('F', 'T')) @pytest.mark.parametrize("python", (False, True)) def test_fr_passive(direct_flag,python): - fr_real = (1.777366E-1, 1.868689E-1) - fr_imag = (-2.331033E-4, -1.671967E-4) + fr_real = (1.947713E-1, 1.990873E-1) + fr_imag = (-2.174952E-4, -1.560016E-4) assert ThinCurr_setup("tw_test-passive.h5",3,direct_flag,eta=1.E4,freq=5.E3,fr_limit=0, icoils=((0.5, 0.1),), vcoils=((0.5, 0.0),), @@ -633,7 +633,7 @@ def test_td_aca(python): @pytest.mark.coverage @pytest.mark.parametrize("python", (False, True)) def test_td_volt_aca(python): - sigs_final = (4.E-3, 1.720867E-4, 1.470760E-4) + sigs_final = (4.E-3, 1.512679E-4, 1.291681E-4) assert ThinCurr_setup("tw_test-cyl_hr.h5",1,'F',use_aca=True, vcoils=((1.1, 0.25), (1.1, -0.25)), floops=((0.9, 0.5), (0.9, 0.0)),