From ce62ab8bf91d17a32d01aa66f4fee13307b6d1ae Mon Sep 17 00:00:00 2001 From: dpanici Date: Thu, 13 Apr 2023 10:24:58 -0400 Subject: [PATCH 1/4] make it easy to make the tor current a helicity multiple of poloidal current --- makefile | 4 ++-- regcoil_build_matrices.f90 | 2 +- regcoil_init_plasma_mod.f90 | 9 +++++++-- regcoil_variables.f90 | 4 +++- 4 files changed, 13 insertions(+), 6 deletions(-) diff --git a/makefile b/makefile index 3d2eb6f..77724ae 100644 --- a/makefile +++ b/makefile @@ -68,7 +68,7 @@ else ifeq ($(HOSTNAME),pppl_intel) LIBSTELL_DIR=$(STELLOPT_PATH)/LIBSTELL/Release LIBSTELL_FOR_REGCOIL=$(LIBSTELL_DIR)/libstell.a REGCOIL_COMMAND_TO_SUBMIT_JOB = srun -N 1 -n 1 -c 8 -q debug --mem 8G -else ifeq ($(HOSTNAME),stellar) +else ifeq ($(HOSTNAME),stellar-intel.princeton.edu) REGCOIL_HOST=stellar FC = ifort NETCDF_F = $(NETCDFDIR) @@ -78,7 +78,7 @@ else ifeq ($(HOSTNAME),stellar) -lpthread \ -I$(NETCDF_F)/include -I$(NETCDF_C)/include EXTRA_LINK_FLAGS = -fopenmp -L$(NETCDF_C)/lib -lnetcdf -L$(NETCDF_F)/lib -lnetcdff -lmkl_gf_lp64 -lmkl_core -lmkl_sequential -lhdf5_hl -lhdf5_fortran -lhdf5 -lpthread -lz -lm - REGCOIL_COMMAND_TO_SUBMIT_JOB = srun -N 1 -n 1 -c 8 -q debug --mem 8G + REGCOIL_COMMAND_TO_SUBMIT_JOB = srun -N 1 -n 1 -c 8 -q debug --mem 8G --time 00:10:00 LIBSTELL_FOR_REGCOIL=$(LIBSTELL_DIR)/mini_libstell.a else ifeq ($(HOSTNAME),osx_brew) REGCOIL_HOST=osx_brew diff --git a/regcoil_build_matrices.f90 b/regcoil_build_matrices.f90 index 0d62644..66973fd 100644 --- a/regcoil_build_matrices.f90 +++ b/regcoil_build_matrices.f90 @@ -42,7 +42,7 @@ subroutine regcoil_build_matrices() if (verbose) print *,"Initializing basis functions and f" ! Initialize Fourier arrays - call regcoil_init_Fourier_modes(mpol_potential, ntor_potential, mnmax_potential, xm_potential, xn_potential, .false.) + call regcoil_init_Fourier_modes(mpol_potential, ntor_potential, mnmax_potential, xm_potential, xn_potential, .false., helicity_ratio) xn_potential = xn_potential * nfp select case (symmetry_option) diff --git a/regcoil_init_plasma_mod.f90 b/regcoil_init_plasma_mod.f90 index 04a112a..510cef2 100644 --- a/regcoil_init_plasma_mod.f90 +++ b/regcoil_init_plasma_mod.f90 @@ -497,12 +497,17 @@ subroutine regcoil_init_plasma() case (2,3,4) ! A VMEC wout file is available ! VMEC stores the toroidal Boozer component B_zeta as "bvco", using the HALF mesh - net_poloidal_current_Amperes = 2*pi/mu0*(1.5_dp*bvco(ns)-0.5_dp*bvco(ns-1)) + net_poloidal_current_Amperes = 2*pi/mu0*(1.5_dp*bvco(ns)-0.5_dp*bvco(ns-1)) - 200.0*13264.873 + ! subtract the current that is already taken care of by the TF coils? + + net_toroidal_current_Amperes = net_poloidal_current_Amperes / helicity_ratio / nfp ! curpol is a number which multiplies the data in the bnorm file. - curpol = (2*pi/nfp)*(1.5_dp*bsubvmnc(1,ns) - 0.5_dp*bsubvmnc(1,ns-1)) + curpol = 1!I am supplying it in T already so should be ok... (2*pi/nfp)*(1.5_dp*bsubvmnc(1,ns) - 0.5_dp*bsubvmnc(1,ns-1)) if (verbose) print *,"Overriding net_poloidal_current_Amperes with value from the VMEC wout file." if (verbose) print *,"G = ", net_poloidal_current_Amperes, " ; curpol = ", curpol + if (verbose) print *,"I = ", net_toroidal_current_Amperes, " ; helicity_ratio = ", helicity_ratio + case default if (abs(net_poloidal_current_Amperes-1)<1e-12) then if (verbose) print *,"No VMEC file is available, and the default value of net_poloidal_current_Amperes (=1) will be used." diff --git a/regcoil_variables.f90 b/regcoil_variables.f90 index 8e49cdc..9d1dd2f 100644 --- a/regcoil_variables.f90 +++ b/regcoil_variables.f90 @@ -66,6 +66,8 @@ module regcoil_variables real(dp), dimension(:), allocatable :: rmns_plasma, zmnc_plasma, rmnc_plasma, zmns_plasma real(dp), dimension(:), allocatable :: rmns_coil, zmnc_coil, rmnc_coil, zmns_coil integer :: nfp + integer :: helicity_ratio ! to produce coils with this helicity i.e. contours of constant curr potential lie on curves of +!theta = 2*pi - helicity_ratio*Nfp*zeta logical :: lasym integer :: max_mpol_coil = 24, max_ntor_coil = 24 ! These variables are upper limits on the # of Fourier modes used to describe a uniform-offset coil surface. integer :: mpol_coil_filter = 24, ntor_coil_filter = 24 @@ -169,7 +171,7 @@ module regcoil_variables target_option, target_value, lambda_search_tolerance, & sensitivity_option, nmax_sensitivity, mmax_sensitivity, & sensitivity_symmetry_option, target_option_p, & - fixed_norm_sensitivity_option, coil_plasma_dist_lse_p + fixed_norm_sensitivity_option, coil_plasma_dist_lse_p, helicity_ratio end module regcoil_variables From d83b56e93ef91885a23d737935f0e680e2d2dfaa Mon Sep 17 00:00:00 2001 From: dpanici Date: Thu, 13 Apr 2023 10:27:40 -0400 Subject: [PATCH 2/4] add var for coil current --- regcoil_variables.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/regcoil_variables.f90 b/regcoil_variables.f90 index 9d1dd2f..699ef08 100644 --- a/regcoil_variables.f90 +++ b/regcoil_variables.f90 @@ -66,7 +66,7 @@ module regcoil_variables real(dp), dimension(:), allocatable :: rmns_plasma, zmnc_plasma, rmnc_plasma, zmns_plasma real(dp), dimension(:), allocatable :: rmns_coil, zmnc_coil, rmnc_coil, zmns_coil integer :: nfp - integer :: helicity_ratio ! to produce coils with this helicity i.e. contours of constant curr potential lie on curves of + integer :: helicity_ratio = 1 ! to produce coils with this helicity i.e. contours of constant curr potential lie on curves of !theta = 2*pi - helicity_ratio*Nfp*zeta logical :: lasym integer :: max_mpol_coil = 24, max_ntor_coil = 24 ! These variables are upper limits on the # of Fourier modes used to describe a uniform-offset coil surface. @@ -87,6 +87,7 @@ module regcoil_variables real(dp) :: net_poloidal_current_Amperes = 1 real(dp) :: net_toroidal_current_Amperes = 0 + real(dp) :: TF_coil_pol_current = 0 ! current in external TF coils logical :: load_bnorm = .false. character(len=200) :: bnorm_filename="" real(dp) :: curpol = 1 ! number which multiplies data in bnorm file. @@ -171,7 +172,7 @@ module regcoil_variables target_option, target_value, lambda_search_tolerance, & sensitivity_option, nmax_sensitivity, mmax_sensitivity, & sensitivity_symmetry_option, target_option_p, & - fixed_norm_sensitivity_option, coil_plasma_dist_lse_p, helicity_ratio + fixed_norm_sensitivity_option, coil_plasma_dist_lse_p, helicity_ratio, TF_coil_pol_current end module regcoil_variables From fa27fd907ec6f40b02b1f916f3ed26243caca7bb Mon Sep 17 00:00:00 2001 From: dpanici Date: Thu, 13 Apr 2023 10:30:27 -0400 Subject: [PATCH 3/4] fix error: --- regcoil_build_matrices.f90 | 2 +- regcoil_init_plasma_mod.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/regcoil_build_matrices.f90 b/regcoil_build_matrices.f90 index 66973fd..0d62644 100644 --- a/regcoil_build_matrices.f90 +++ b/regcoil_build_matrices.f90 @@ -42,7 +42,7 @@ subroutine regcoil_build_matrices() if (verbose) print *,"Initializing basis functions and f" ! Initialize Fourier arrays - call regcoil_init_Fourier_modes(mpol_potential, ntor_potential, mnmax_potential, xm_potential, xn_potential, .false., helicity_ratio) + call regcoil_init_Fourier_modes(mpol_potential, ntor_potential, mnmax_potential, xm_potential, xn_potential, .false.) xn_potential = xn_potential * nfp select case (symmetry_option) diff --git a/regcoil_init_plasma_mod.f90 b/regcoil_init_plasma_mod.f90 index 510cef2..d2d5554 100644 --- a/regcoil_init_plasma_mod.f90 +++ b/regcoil_init_plasma_mod.f90 @@ -497,7 +497,7 @@ subroutine regcoil_init_plasma() case (2,3,4) ! A VMEC wout file is available ! VMEC stores the toroidal Boozer component B_zeta as "bvco", using the HALF mesh - net_poloidal_current_Amperes = 2*pi/mu0*(1.5_dp*bvco(ns)-0.5_dp*bvco(ns-1)) - 200.0*13264.873 + net_poloidal_current_Amperes = 2*pi/mu0*(1.5_dp*bvco(ns)-0.5_dp*bvco(ns-1)) - TF_coil_pol_current ! subtract the current that is already taken care of by the TF coils? net_toroidal_current_Amperes = net_poloidal_current_Amperes / helicity_ratio / nfp From 75a08dd735a4e5786e30d8aa72255020e5addb74 Mon Sep 17 00:00:00 2001 From: dpanici Date: Sun, 21 May 2023 21:45:50 -0400 Subject: [PATCH 4/4] add if statement to allow for helicity of zero, and save the helicity ratio in the output --- regcoil_init_plasma_mod.f90 | 7 +++++-- regcoil_write_output.f90 | 7 ++++++- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/regcoil_init_plasma_mod.f90 b/regcoil_init_plasma_mod.f90 index d2d5554..fd9c532 100644 --- a/regcoil_init_plasma_mod.f90 +++ b/regcoil_init_plasma_mod.f90 @@ -499,8 +499,11 @@ subroutine regcoil_init_plasma() ! VMEC stores the toroidal Boozer component B_zeta as "bvco", using the HALF mesh net_poloidal_current_Amperes = 2*pi/mu0*(1.5_dp*bvco(ns)-0.5_dp*bvco(ns-1)) - TF_coil_pol_current ! subtract the current that is already taken care of by the TF coils? - - net_toroidal_current_Amperes = net_poloidal_current_Amperes / helicity_ratio / nfp + if (helicity_ratio .ne. 0) then + net_toroidal_current_Amperes = net_poloidal_current_Amperes / helicity_ratio / nfp + else + net_toroidal_current_Amperes = 0 + end if ! curpol is a number which multiplies the data in the bnorm file. curpol = 1!I am supplying it in T already so should be ok... (2*pi/nfp)*(1.5_dp*bsubvmnc(1,ns) - 0.5_dp*bsubvmnc(1,ns-1)) diff --git a/regcoil_write_output.f90 b/regcoil_write_output.f90 index bf7eada..46e7089 100644 --- a/regcoil_write_output.f90 +++ b/regcoil_write_output.f90 @@ -54,7 +54,8 @@ subroutine regcoil_write_output vn_coil_plasma_dist_max = "coil_plasma_dist_max", & vn_coil_plasma_dist_min_lse = "coil_plasma_dist_min_lse", & vn_coil_plasma_dist_max_lse = "coil_plasma_dist_max_lse", & - vn_coil_plasma_dist_lse_p = "coil_plasma_dist_lse_p" + vn_coil_plasma_dist_lse_p = "coil_plasma_dist_lse_p", & + vn_helicity_ratio = "helicity_ratio" ! Arrays with dimension 1 character(len=*), parameter :: & @@ -289,6 +290,9 @@ subroutine regcoil_write_output call cdf_setatt(ncid, vn_net_toroidal_current_Amperes, 'Net current (in Amperes) that flows on the coil winding surface in the toroidal direction. ' // & 'This quantity corresponds to I in the 2017 Nuclear Fusion paper. For modular coils, this quantity is 0.') + call cdf_define(ncid, vn_helicity_ratio, helicity_ratio) + call cdf_setatt(ncid, vn_helicity_ratio, 'Helicity of the helical coils (i.e. number of poloidal turns per field period).') + call cdf_define(ncid, vn_curpol, curpol) call cdf_define(ncid, vn_nlambda, nlambda) @@ -544,6 +548,7 @@ subroutine regcoil_write_output call cdf_write(ncid, vn_net_poloidal_current_Amperes, net_poloidal_current_Amperes) call cdf_write(ncid, vn_net_toroidal_current_Amperes, net_toroidal_current_Amperes) call cdf_write(ncid, vn_curpol, curpol) + call cdf_write(ncid, vn_helicity_ratio, helicity_ratio) call cdf_write(ncid, vn_nlambda, nlambda) call cdf_write(ncid, vn_total_time, total_time) call cdf_write(ncid, vn_exit_code, exit_code)