Merge pull request #42 from PrincetonUniversity/develop

Separate plasma surface with limiter surface for coil-surface separation AND stellarator symmetry

.gitignore

@@ -7,7 +7,7 @@
 
 *.[oa]
 *.mod
-*.F90
+*_m.F90
 *.pyc
 *.log
 *.h5

examples/limiter_surface/ellipse.boundary

@@ -0,0 +1,13 @@
+#bmn bNfp nbf
+4 2 0
+#plasma boundary
+# n m Rbc Rbs Zbc Zbs
+0 0 3.00 0.0 0.0 0.00
+0 1 0.30 0.0 0.0 -0.30
+1 0 0.00 0.0 0.0 -0.06
+1 1 -0.06 0.0 0.0 -0.06
+#Bn harmonics
+# n m bnc bns
+0 0 1.0 0.0
+0 1 0.5 0.25
+1 0 0.5 0.0

examples/limiter_surface/ellipse.limiter

@@ -0,0 +1,13 @@
+#bmn bNfp nbf
+4 2 0
+#plasma boundary
+# n m Rbc Rbs Zbc Zbs
+0 0 3.00 0.0 0.0 0.00
+0 1 0.50 0.0 0.0 -0.50
+1 0 0.00 0.0 0.0 -0.06
+1 1 -0.06 0.0 0.0 -0.06
+#Bn harmonics
+# n m bnc bns
+0 0 1.0 0.0
+0 1 0.5 0.25
+1 0 0.5 0.0

examples/limiter_surface/limiter.input

@@ -0,0 +1,75 @@
+ &FOCUSIN
+ ISQUIET = -1,
+ ISSYMMETRIC = 0,
+ INPUT_SURF = 'ellipse.boundary',
+ !LIMITER_SURF = 'ellipse.limiter',
+ INPUT_HARM = 'target.harmonics',
+ INPUT_COILS = 'none',
+ CASE_SURFACE = 0,
+ KNOTSURF = 0.200000000000000 ,
+ ELLIPTICITY = 0.000000000000000E+000,
+ NTETA = 64,
+ NZETA = 64,
+ CASE_INIT = 1,
+ CASE_COILS = 1,
+ NCOILS = 16,
+ INIT_CURRENT = 1000000.00000000 ,
+ INIT_RADIUS = 0.60000000000000 ,
+ ISVARYCURRENT = 0,
+ ISVARYGEOMETRY = 1,
+ NFCOIL = 4,
+ NSEG = 128,
+ ISNORMALIZE = 1,
+ ISNORMWEIGHT = 0,
+ CASE_BNORMAL = 0,
+ CASE_LENGTH = 1,
+ WEIGHT_BNORM = 1.00000000000000 ,
+ BHARM_JSURF = 0,
+ WEIGHT_BHARM = 0.000000000000000E+000,
+ WEIGHT_TFLUX = 0.000000000000000E+000,
+ TARGET_TFLUX = 0.000000000000000E+000,
+ WEIGHT_TTLEN = 1.000000000000000E-003,
+ TARGET_LENGTH = 5.000000000000000E+000,
+ WEIGHT_CSSEP = 0.000000000000000E-004,
+ CSSEP_FACTOR = 1.00000000000000 ,
+ WEIGHT_SPECW = 0.000000000000000E+000,
+ WEIGHT_CCSEP = 0.000000000000000E+000,
+ WEIGHT_INORM = 1.00000000000000 ,
+ WEIGHT_GNORM = 1.00000000000000 ,
+ WEIGHT_MNORM = 1.00000000000000 ,
+ CASE_OPTIMIZE = 1,
+ EXIT_TOL = 1.000000000000000E-004,
+ DF_MAXITER = 0,
+ DF_XTOL = 1.000000000000000E-008,
+ DF_TAUSTA = 0.000000000000000E+000,
+ DF_TAUEND = 1.00000000000000 ,
+ CG_MAXITER = 20,
+ CG_XTOL = 1.000000000000000E-008,
+ CG_WOLFE_C1 = 0.100000001490116 ,
+ CG_WOLFE_C2 = 0.899999976158142 ,
+ LM_MAXITER = 0,
+ LM_XTOL = 1.000000000000000E-008,
+ LM_FTOL = 1.000000000000000E-008,
+ LM_FACTOR = 100.000000000000 ,
+ HN_MAXITER = 0,
+ HN_XTOL = 1.000000000000000E-008,
+ HN_FACTOR = 100.000000000000 ,
+ TN_MAXITER = 0,
+ TN_REORDER = 0,
+ TN_XTOL = 1.000000000000000E-008,
+ TN_CR = 0.100000001490116 ,
+ CASE_POSTPROC = 3,
+ SAVE_FREQ = 1,
+ SAVE_COILS = 1,
+ SAVE_HARMONICS = 0,
+ SAVE_FILAMENTS = 0,
+ UPDATE_PLASMA = 0,
+ PP_PHI = 0.000000000000000E+000,
+ PP_RAXIS = 0.000000000000000E+000,
+ PP_ZAXIS = 0.000000000000000E+000,
+ PP_RMAX = 0.000000000000000E+000,
+ PP_ZMAX = 0.000000000000000E+000,
+ PP_NS = 10,
+ PP_MAXITER = 1000,
+ PP_XTOL = 1.000000000000000E-006
+ /

examples/rotating_ellipse/ellipse.input

@@ -1,33 +1,33 @@
 &focusin
  IsQuiet = -1 ! -2 verbose and including unconstrained cost functions; -1: verbose; 0: normal; 1: concise
- IsSymmetric = 0 ! 0: no stellarator symmetry enforced; 1: plasma periodicity enforced; 2: coil periodicity enforced
+ IsSymmetric = 2 ! 0: no stellarator symmetry enforced; 1: plasma periodicity enforced; 2: coil periodicity enforced
 
  case_surface = 0 ! 0: general VMEC-like format (Rbc, Rbs, Zbc, Zbs); 1: read axis for knots
  knotsurf = 0.200D-00 ! minor plasma radius for knototrans, only valid for case surface = 1
  ellipticity = 0.000D+00 ! ellipticity of plasma for knototrans, only valid for case surface = 1
- Nteta = 128 ! poloidal number for discretizing the surface
- Nzeta = 128 ! toroidal number for discretizing the surface
+ Nteta =  64 ! poloidal number for discretizing the surface
+ Nzeta =  64 ! toroidal number for discretizing the surface
 
  case_init = 1 ! -1: read coils.ext file; 0: read ext.focus file; 1: initialize with circular coils; 2: initialize with dipoles
  case_coils = 1 ! 0: using piecewise linear representation; (not ready); 1: using Fourier series representation
- Ncoils = 16 ! number of coils; only valid when case_init = 1
+ Ncoils = 4  ! number of coils; only valid when case_init = 1
  init_current = 1.000D+06 ! initial coil currents (Amper); only valid when case_init = 1
  init_radius = 0.500D+00 ! initial coil radius (meter); only valid when case_init = 1
- IsVaryCurrent = 1 ! 0: all the currents fixed; 1: currents can be changed; overwritten by ext.focus
+ IsVaryCurrent = 0 ! 0: all the currents fixed; 1: currents can be changed; overwritten by ext.focus
  IsVaryGeometry = 1 ! 0: all the geometries fixed; 1: geometries can be changed; overwritten by ext.focus
  NFcoil = 4 ! number of Fourier harmonics representing the coils; overwritten by ext.focus
  Nseg = 128 ! number of coil segments for discretizing; overwritten by ext.focus
 
  IsNormalize = 1 ! 0: do not normalize coil parameters; 1: normalize; I = I/I0, x = x/R0; I0 & R0 are quadrtic mean values.
  IsNormWeight = 1 ! 0: do not normalize the weights; 1: normalize the weights
- case_bnormal = 1 ! 0: keep raw Bn error; 1: Bn residue normalized to local |B|
+ case_bnormal = 0 ! 0: keep raw Bn error; 1: Bn residue normalized to local |B|
  case_length = 1 ! 1: quadratic format, converging the target length; 2: exponential format, as short as possible
  weight_bnorm = 1.000D+02 ! weight for real space Bn errors
  weight_bharm = 0.000D+00 ! weight for Bnm harmonic errors
  weight_tflux = 0.000D+00 ! weight for toroidal flux error
  target_tflux = 0.000D+00 ! target for the toroidal flux
- weight_ttlen = 0.100D-02 ! weight for coil length error
- target_length = 3.500D+00 ! target value (or for normalization) of the coils length, if zero, automatically set to initial actual length
+ weight_ttlen = 1.000D+01 ! weight for coil length error
+ target_length = 7.000D+00 ! target value (or for normalization) of the coils length, if zero, automatically set to initial actual length
  weight_specw = 0.000D+00 ! weight for spectral condensation error
  weight_cssep = 0.010D+00 ! weight for coil-surface separation constraint
  weight_inorm = 1.000D+00 ! weight for normalization of current. Larger weight makes the derivatives more important.
@@ -46,7 +46,7 @@
  CG_wolfe_c1 = 0.1 ! c1 value in the strong wolfe condition for line search, (0.0, 0.5) 
  CG_wolfe_c2 = 0.9 ! c2 value in the strong wolfe condition for line search; 0 < c1 < c2 < 1
 
- LM_maxiter = 0 ! maximum iterations allowed for using Levenberg-Marquard (LM)
+ LM_maxiter = 5 ! maximum iterations allowed for using Levenberg-Marquard (LM)
  LM_xtol = 1.000D-08 ! if the relative error between two consecutivec iterates is at most xtol, the optimization terminates
  LM_ftol = 1.000D-08 ! if both the actual and predicted relative reductions in the sum of squares are at most ftol, the optimization terminates;
  LM_factor = 100.0 ! the initial step bound, which is set to the product of factor and the euclidean norm of diag*x if nonzero
@@ -66,5 +66,17 @@
  pp_ns = 10 ! number of following fieldlines
  pp_maxiter = 1000 ! number of periods for each fieldline following
  pp_xtol = 1.000D-06 ! tolarence of ODE solver during fieldline fowllowing
-
+/
+&mgrid
+! mgrid file dimensions
+Rmin = 0.0
+Rmax = 0.0
+Zmin = 0.0
+Zmax = 0.0
+Pmin = 0.0
+Pmax = 6.283
+! resolutions
+NR = 101
+NZ = 101
+NP = 72
 /

sources/Makefile

@@ -2,7 +2,7 @@
 
 ############################################################################################################
 
- ALLFILES= globals initial rdsurf rdknot rdcoils packdof bfield bmnharm bnormal fdcheck \
+ ALLFILES= globals initial surface rdsurf rdknot rdcoils packdof bfield bmnharm bnormal fdcheck \
  torflux length surfsep datalloc solvers descent congrad lmalg saving diagnos \
  specinp poinplot boozer wtmgrid focus
  HFILES= $(ALLFILES:=.f90) # raw source files