Skip to content

Commit

Permalink
ThinCurr: Use Hurwitz and Landreman limiting form of self-inductance …
Browse files Browse the repository at this point in the history 
…for filamentary coils (#73)

- Update tests using filament self-inductance
  • Loading branch information
@hansec
hansec committed Aug 19, 2024
1 parent e89a5a9 commit 1789073
Show file tree
Hide file tree
Showing 2 changed files with 27 additions and 13 deletions.
20 changes: 17 additions & 3 deletions src/physics/thin_wall.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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<coil_min_rad))CALL oft_abort("Invalid radius for passive coil", &
"tw_setup", __FILE__)
END DO
WRITE(*,'(2A)')oft_indent,'Loading I(t) driver coils'
CALL xml_get_element(self%xml,"icoils",coil_element,error_flag)
IF(error_flag==0)CALL tw_load_coils(coil_element,self%n_icoils,self%icoils)
DO i=1,self%n_icoils
DO j=1,self%icoils(i)%ncoils
self%icoils(i)%radius(j)=MAX(coil_min_rad,self%icoils(i)%radius(j)) ! Remove dummy radius on Icoils
END DO
END DO
#endif
WRITE(*,*)
! WRITE(*,'(2A)')oft_indent,'Thin-wall model loaded:'
Expand Down Expand Up @@ -692,13 +700,18 @@ SUBROUTINE tw_compute_Ael2dr(tw_obj,save_file)
END SUBROUTINE tw_compute_Ael2dr
!------------------------------------------------------------------------------
!> 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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
20 changes: 10 additions & 10 deletions src/tests/physics/test_ThinCurr.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)),
Expand All @@ -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)),
Expand All @@ -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)),
Expand All @@ -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),),
Expand All @@ -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)),
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)))
Expand Down Expand Up @@ -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),),
Expand Down Expand Up @@ -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)),
Expand Down

0 comments on commit 1789073

Please sign in to comment.