Merge pull request #20 from PrincetonUniversity/develop

Develop

.gitignore

@@ -22,4 +22,4 @@
 bin/
 Old/
 docs/
-pyfocus/
+pyfocus/

README.md

@@ -1,6 +1,6 @@
 # <img alt="FOCUS" src="https://joshburns.net/blog/wp-content/uploads/2013/08/focus.jpg" height="80"> Flexible Optimized Coils Using Space curves
 
-**FOCUS**is a nonlinear optimization code for designing 3D coils. 
+**FOCUS** is a nonlinear optimization code for designing 3D coils. 
 
 - **Website (including documentation):** https://princetonuniversity.github.io/FOCUS/
 - **Source:** https://github.com/PrincetonUniversity/FOCUS
@@ -15,4 +15,6 @@ There are several branches available. Please use the correct one.
 - **gh-pages:** the branch hosting GitHub pages for the offcicial website.
 - **others:** task-oriented branches.
 
+
 *Photo credit: iStockphoto.com*
+

examples/d3d_RMP/d3d.input

@@ -12,7 +12,7 @@
  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
  init_current = 1.000D+06 ! initial coil currents (Amper); only valid when case_init = 1
- init_radius = 0.500D+00 ! 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
  IsVaryGeometry = 0 ! 0: all the geometries fixed; 1: geometries can be changed; overwritten by ext.focus
  NFcoil = 8 ! number of Fourier harmonics representing the coils; overwritten by ext.focus
@@ -27,11 +27,13 @@
  weight_tflux = 0.000D+00 ! weight for toroidal flux error
  target_tflux = 0.000D+00 ! target for the toroidal flux
  weight_ttlen = 0.000D+00 ! weight for coil length error
+ weight_cssep = 1.000D+00 ! weight for coil surface separation constraint
+ cssep_factor = 3.00D+00 ! exponential factor for cssep
  target_length = 0.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_ccsep = 0.000D+00 ! weight for coil-coil separation constraint
- weight_inorm = 1.000D+00 ! weight for normalization of current. Larger weight makes the derivatives more important.
- weight_gnorm = 1.000D+00 ! weight for normalization of geometric coefficients. Larger weight makes the derivatives more important.
+ weight_specw = 0.000D+00 ! weight for spectral condensation error (not ready)
+ weight_ccsep = 0.000D+00 ! weight for coil-coil separation constraint (not ready)
+ weight_inorm = 1.000D+00   ! weight for normalization of current. Larger weight makes the derivatives more important.
+ weight_gnorm = 1.000D+00  ! weight for normalization of geometric coefficients. Larger weight makes the derivatives more important.
 
  case_optimize = 1 ! -2: check the 2nd derivatives (not ready); -1: check the 1st derivatives; 0: no optimizations performed; 1: optimizing with algorithms using the gradient (DF and/or CG); 2: optimizing with algorithms using the Hessian (HT and/or NT)
  exit_tol = 1.000D-04 ! Exit the optimizer if the percent change in the cost function over the last 5 steps is below this threshold
@@ -60,4 +62,4 @@
  save_coils = 1 ! flag for indicating whether write example.focus and example.coils
  save_harmonics = 1 ! flag for indicating whether write example.harmonics
  save_filaments = 0 ! flag for indicating whether write .example.filaments.xxxxxx
-/
+/

examples/rotating_ellipse/ellipse.input

@@ -12,15 +12,15 @@
  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
  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
+ 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
  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+00 ! weight for real space Bn errors
  weight_bharm = 0.000D+00 ! weight for Bnm harmonic errors
@@ -29,19 +29,19 @@
  weight_ttlen = 0.100D+00 ! weight for coil length error
  target_length = 0.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_ccsep = 0.000D+00 ! weight for coil-coil separation constraint
+ 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.
  weight_gnorm = 1.000D+00 ! weight for normalization of geometric coefficients. Larger weight makes the derivatives more important.
 
- case_optimize = 1  ! -2: check the 2nd derivatives (not ready); -1: check the 1st derivatives; 0: no optimizations performed; 1: optimizing with algorithms using the gradient (DF and/or CG); 2: optimizing with algorithms using the Hessian (HT and/or NT)
+ case_optimize = 1 ! -2: check the 2nd derivatives (not ready); -1: check the 1st derivatives; 0: no optimizations performed; 1: optimizing with algorithms using the gradient (DF and/or CG); 2: optimizing with algorithms using the Hessian (HT and/or NT)
  exit_tol = 1.000D-04 ! Exit the optimizer if the percent change in the cost function over the last 5 steps is below this threshold
 
  DF_maxiter = 0 ! maximum iterations allowed for using Differential Flow (DF)
  DF_xtol = 1.000D-08 ! relative error for ODE solver
  DF_tausta = 0.000D+00 ! starting value of τ. Usually 0.0 is a good idea
  DF_tauend = 1.000D-00 ! ending value of τ. The larger value of τend − τsta, the more optimized
 
- CG_maxiter = 100 ! maximum iterations allowed for using Conjugate Gradient (CG)
+ CG_maxiter = 50  ! maximum iterations allowed for using Conjugate Gradient (CG)
  CG_xtol = 1.000D-08 ! the stopping criteria of finding minimum; if |dχ2/dX| < CG xtol, exit the optimization
  CG_wolfe_c1 = 1.000D-04 ! c1 value in the strong wolfe condition for line search; 
  CG_wolfe_c2 = 0.9 ! c2 value in the strong wolfe condition for line search; if one CG step takes too long, try to increase c2, but remember 0 < c1 < c2 < 1

sources/Makefile

@@ -3,7 +3,7 @@
 ############################################################################################################
 
  ALLFILES= globals initial datalloc rdsurf rdknot rdcoils packdof bfield bnormal bmnharm fdcheck \
- torflux length solvers descent congrad saving diagnos focus
+ torflux length surfsep solvers descent congrad saving diagnos focus
  HFILES= $(ALLFILES:=.h)
  FFILES= $(ALLFILES:=.F90)
  PFILES= $(ALLFILES:=.pdf)

sources/bfield.h

@@ -1,5 +1,5 @@
 
-!title (field) ! Computes magnetic field.
+!title (bfield) ! Computes magnetic field.
 
 !latex \briefly{Computes magnetic field given coil geometry.}
 
@@ -50,7 +50,7 @@ subroutine bfield0(icoil, iteta, jzeta, Bx, By, Bz)
  dly = zero; lty = zero; By = zero
  dlz = zero; ltz = zero; Bz = zero
 
- do kseg = 1, coil(icoil)%NS
+ do kseg = 0, coil(icoil)%NS-1
 
  dlx = surf(1)%xx(iteta,jzeta) - coil(icoil)%xx(kseg)
  dly = surf(1)%yy(iteta,jzeta) - coil(icoil)%yy(kseg)
@@ -92,7 +92,7 @@ subroutine bfield1(icoil, iteta, jzeta, Bx, By, Bz, ND)
 
  INTEGER :: ierr, astat, kseg, NS
  REAL :: dlx, dly, dlz, r, rm3, rm5, ltx, lty, ltz, rxp
- REAL, dimension(1:1, 1:coil(icoil)%NS) :: dBxx, dBxy, dBxz, dByx, dByy, dByz, dBzx, dBzy, dBzz
+ REAL, dimension(1:1, 0:coil(icoil)%NS-1) :: dBxx, dBxy, dBxz, dByx, dByy, dByz, dBzx, dBzy, dBzz
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
@@ -107,7 +107,7 @@ subroutine bfield1(icoil, iteta, jzeta, Bx, By, Bz, ND)
  dly = zero; lty = zero; By = zero
  dlz = zero; ltz = zero; Bz = zero
 
- do kseg = 1, NS
+ do kseg = 0, NS-1
 
  dlx = surf(1)%xx(iteta,jzeta) - coil(icoil)%xx(kseg)
  dly = surf(1)%yy(iteta,jzeta) - coil(icoil)%yy(kseg)

sources/congrad.h

@@ -49,8 +49,8 @@ SUBROUTINE congrad
  p(1:Ndof) = -gradk ! initial step direction;
  alpha = 1.0 ! initial step size;
 
- if (myid == 0) write(ounit, '("output : "A6" : "9(A12," ; "))') "iout", "time (s)", "chi", "dE_norm", &
- "Bnormal", "Bmn harmonics", "tor. flux", "coil length", "spectral", "c-c sep." 
+ if (myid == 0) write(ounit, '("output : "A6" : "8(A12," ; "))') "iout", "time (s)", "chi", "dE_norm", &
+ "Bnormal", "Bmn harmonics", "tor. flux", "coil length", "c-s sep." 
 
  do
 

sources/datalloc.h

@@ -28,9 +28,9 @@ subroutine AllocData(itype)
  ND = (6*NF + 3) ! total variables for geometry
  DoF(icoil)%ND = coil(icoil)%Lc * ND !# of DoF for icoil;
  SALLOCATE(DoF(icoil)%xdof, (1:DoF(icoil)%ND), zero)
- SALLOCATE(DoF(icoil)%xof , (1:coil(icoil)%NS, 1:ND), zero)
- SALLOCATE(DoF(icoil)%yof , (1:coil(icoil)%NS, 1:ND), zero)
- SALLOCATE(DoF(icoil)%zof , (1:coil(icoil)%NS, 1:ND), zero)
+ SALLOCATE(DoF(icoil)%xof , (0:coil(icoil)%NS-1, 1:ND), zero)
+ SALLOCATE(DoF(icoil)%yof , (0:coil(icoil)%NS-1, 1:ND), zero)
+ SALLOCATE(DoF(icoil)%zof , (0:coil(icoil)%NS-1, 1:ND), zero)
  case default
  FATAL(AllocData, .true., not supported coil types)
  end select
@@ -43,7 +43,7 @@ subroutine AllocData(itype)
  Tdof = Tdof + 1 + 6*(FouCoil(icoil)%NF)+3
  if (DoF(icoil)%ND >= Cdof) Cdof = DoF(icoil)%ND ! find the largest ND for single coil;
 
- enddo
+  enddo
 
  if(Ndof == 0) then ! no DOF;
  Nouts = 0
@@ -52,7 +52,7 @@ subroutine AllocData(itype)
 
  SALLOCATE( xdof, (1:Ndof), zero ) ! dof vector;
  SALLOCATE( dofnorm, (1:Ndof), zero ) ! dof normalized value vector;
- SALLOCATE( evolution, (1:Nouts+1, 0:8), zero ) !evolution array;
+ SALLOCATE( evolution, (1:Nouts+1, 0:7), zero ) !evolution array;
  SALLOCATE( coilspace, (1:Nouts+1, 1:Tdof), zero ) ! all the coil parameters;
 
  idof = 0
@@ -125,6 +125,11 @@ subroutine AllocData(itype)
  SALLOCATE( t1L, (1:Ndof), zero )
  endif
 
+ ! cssep needed;
+ if (weight_cssep > sqrtmachprec) then
+ SALLOCATE( t1S, (1:Ndof), zero )
+ endif
+
  endif
  !--------------------------------------------------------------------------------------------- 
 

sources/descent.h

@@ -50,8 +50,8 @@ subroutine descent
  SALLOCATE( work, (1:100+21*Ndof), zero )
 
  call denergy(tau, lxdof, dE)
- if (myid == 0) write(ounit, '("output : "A6" : "9(A12," ; "))') "iout", "tau", "chi", "dE_norm", &
- "Bnormal", "Bmn harmonics", "tor. flux", "coil length", "spectral", "c-c sep." 
+ if (myid == 0) write(ounit, '("output : "A6" : "8(A12," ; "))') "iout", "tau", "chi", "dE_norm", &
+ "Bnormal", "Bmn harmonics", "tor. flux", "coil length", "c-s sep." 
  !call output(t0)
 
  do itau = 1, DF_maxiter

sources/diagnos.h

@@ -6,7 +6,8 @@ SUBROUTINE diagnos
 ! diagonose the coil performance
 !------------------------------------------------------------------------------------------------------ 
  use globals, only: zero, one, myid, ounit, sqrtmachprec, IsQuiet, case_optimize, coil, surf, Ncoils, &
- Nteta, Nzeta, bnorm, bharm, tflux, ttlen, specw, ccsep, coilspace, FouCoil, iout, Tdof, case_length
+ Nteta, Nzeta, bnorm, bharm, tflux, ttlen, specw, ccsep, coilspace, FouCoil, iout, Tdof, case_length, &
+ cssep
 
  implicit none
  include "mpif.h"
@@ -23,9 +24,9 @@ SUBROUTINE diagnos
  if (case_optimize == 0) call AllocData(0) ! if not allocate data;
  call costfun(0)
 
- if (myid == 0) write(ounit, '("diagnos : "6(A12," ; "))') , &
- "Bnormal", "Bmn harmonics", "tor. flux", "coil length", "spectral", "c-c sep." 
- if (myid == 0) write(ounit, '(" : "6(ES12.5," ; "))') bnorm, bharm, tflux, ttlen, specw, ccsep
+ if (myid == 0) write(ounit, '("diagnos : "5(A12," ; "))') , &
+ "Bnormal", "Bmn harmonics", "tor. flux", "coil length", "c-s sep." 
+ if (myid == 0) write(ounit, '(" : "6(ES12.5," ; "))') bnorm, bharm, tflux, ttlen, cssep
 
  !save all the coil parameters;
  if (allocated(coilspace)) then
@@ -84,16 +85,16 @@ SUBROUTINE diagnos
  itmp = icoil + 1
  if(icoil .eq. Ncoils) itmp = 1
 
- SALLOCATE(Atmp, (1:3,1:coil(icoil)%NS), zero)
- SALLOCATE(Btmp, (1:3,1:coil(itmp )%NS), zero)
+ SALLOCATE(Atmp, (1:3,0:coil(icoil)%NS-1), zero)
+ SALLOCATE(Btmp, (1:3,0:coil(itmp )%NS-1), zero)
 
- Atmp(1, 1:coil(icoil)%NS) = coil(icoil)%xx(1:coil(icoil)%NS)
- Atmp(2, 1:coil(icoil)%NS) = coil(icoil)%yy(1:coil(icoil)%NS)
- Atmp(3, 1:coil(icoil)%NS) = coil(icoil)%zz(1:coil(icoil)%NS)
+ Atmp(1, 0:coil(icoil)%NS-1) = coil(icoil)%xx(0:coil(icoil)%NS-1)
+ Atmp(2, 0:coil(icoil)%NS-1) = coil(icoil)%yy(0:coil(icoil)%NS-1)
+ Atmp(3, 0:coil(icoil)%NS-1) = coil(icoil)%zz(0:coil(icoil)%NS-1)
 
- Btmp(1, 1:coil(itmp)%NS) = coil(itmp)%xx(1:coil(itmp)%NS)
- Btmp(2, 1:coil(itmp)%NS) = coil(itmp)%yy(1:coil(itmp)%NS)
- Btmp(3, 1:coil(itmp)%NS) = coil(itmp)%zz(1:coil(itmp)%NS)
+ Btmp(1, 0:coil(itmp )%NS-1) = coil(itmp)%xx(0:coil(itmp )%NS-1)
+ Btmp(2, 0:coil(itmp )%NS-1) = coil(itmp)%yy(0:coil(itmp )%NS-1)
+ Btmp(3, 0:coil(itmp )%NS-1) = coil(itmp)%zz(0:coil(itmp )%NS-1)
 
  call mindist(Atmp, coil(icoil)%NS, Btmp, coil(itmp)%NS, tmp_dist)
 
@@ -110,12 +111,12 @@ SUBROUTINE diagnos
  minCPdist = infmax
  do icoil = 1, Ncoils
 
- SALLOCATE(Atmp, (1:3,1:coil(icoil)%NS), zero)
+ SALLOCATE(Atmp, (1:3,0:coil(icoil)%NS-1), zero)
  SALLOCATE(Btmp, (1:3,1:(Nteta*Nzeta)), zero)
 
- Atmp(1, 1:coil(icoil)%NS) = coil(icoil)%xx(1:coil(icoil)%NS)
- Atmp(2, 1:coil(icoil)%NS) = coil(icoil)%yy(1:coil(icoil)%NS)
- Atmp(3, 1:coil(icoil)%NS) = coil(icoil)%zz(1:coil(icoil)%NS)
+ Atmp(1, 0:coil(icoil)%NS-1) = coil(icoil)%xx(0:coil(icoil)%NS-1)
+ Atmp(2, 0:coil(icoil)%NS-1) = coil(icoil)%yy(0:coil(icoil)%NS-1)
+ Atmp(3, 0:coil(icoil)%NS-1) = coil(icoil)%zz(0:coil(icoil)%NS-1)
 
  Btmp(1, 1:(Nteta*Nzeta)) = reshape(surf(1)%xx(0:Nteta-1, 0:Nzeta-1), (/Nteta*Nzeta/))
  Btmp(2, 1:(Nteta*Nzeta)) = reshape(surf(1)%yy(0:Nteta-1, 0:Nzeta-1), (/Nteta*Nzeta/))

sources/fdcheck.h

@@ -30,8 +30,9 @@ SUBROUTINE fdcheck( ideriv )
  !--------------------------------------------------------------------------------------------
 
  INTEGER :: astat, ierr, idof
- REAL :: tmp_xdof(1:Ndof), small=1.0E-6, fd, negvalue, posvalue, diff, rdiff
+ REAL :: tmp_xdof(1:Ndof), fd, negvalue, posvalue, diff, rdiff
  REAL :: start, finish
+ REAL, parameter :: small=1.0E-4
  !--------------------------------------------------------------------------------------------
 
  if(myid == 0) write(ounit, *) "-----------Checking derivatives------------------------------"