Easy helical

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

regcoil_init_plasma_mod.f90

@@ -497,12 +497,20 @@ 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)) - TF_coil_pol_current
+       ! subtract the current that is already taken care of by the TF coils?
+       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 = (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."

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 = 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.
   integer :: mpol_coil_filter = 24, ntor_coil_filter = 24
@@ -85,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.
@@ -169,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
+       fixed_norm_sensitivity_option, coil_plasma_dist_lse_p, helicity_ratio, TF_coil_pol_current
 
 end module regcoil_variables
 

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)