Merge pull request #89 from PrincetonUniversity/Itotal

Add a new cost function to constrain the total current

examples/check_deriv/deriv.focus

@@ -4,7 +4,7 @@
  #coil_type coil_symm coil_name
  1 2 1 
  #Nseg current Ifree Length Lfree target_length
- 128 -1.500725500000000E+05 0 2.033518907783376E+00 1 2.200000000000000E+00
+ 128 -1.500725500000000E+05 1 2.033518907783376E+00 1 2.200000000000000E+00
  #NFcoil
  4
  #Fourier harmonics for coils ( xc; xs; yc; ys; zc; zs) 
@@ -18,7 +18,7 @@
  #coil_type coil_symm coil_name
  1 2 2 
  #Nseg current Ifree Length Lfree target_length
- 128 -1.500725500000000E+05 0 2.095801591241713E+00 1 2.300000000000000E+00
+ 128 -1.500725500000000E+05 1 2.095801591241713E+00 1 2.300000000000000E+00
  #NFcoil
  4
  #Fourier harmonics for coils ( xc; xs; yc; ys; zc; zs) 

examples/check_deriv/deriv.input

@@ -23,11 +23,14 @@
  WEIGHT_GNORM= 1.0000000000000000 ,
  WEIGHT_MNORM= 1.0000000000000000 ,
  WEIGHT_BNORM= 1.0000000000000000 ,
+ NPERT = 5,
  WEIGHT_SBNORM = 1.0, 
  CASE_BNORMAL=0 ,
  WEIGHT_BHARM= 0.0000000000000000 ,
  BHARM_JSURF=0 ,
  INPUT_HARM="target.harmonics ",
+ WEIGHT_ISUM = 1.0
+ TARGET_ISUM = -3E6
  WEIGHT_TFLUX= 0.1000000000000000 ,
  TARGET_TFLUX= 0.0000000000000000 ,
  WEIGHT_TTLEN= 0.2000000000000000 ,

examples/lhd/lhd.input

@@ -60,4 +60,4 @@
  save_coils = 1 ! flag for indicating whether write example.focus and example.coils
  save_harmonics = 0 ! flag for indicating whether write example.harmonics
  save_filaments = 0 ! flag for indicating whether write .example.filaments.xxxxxx
-/
+/

sources/Makefile

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

sources/current.f90

@@ -0,0 +1,83 @@
+!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+subroutine current( ideriv )
+!------------------------------------------------------------------------------------------------------ 
+! Calculate the difference between the total FREE current and the target.
+! Formula is isum = ((SUM(I_i) - target_isum) / target_isum)**2
+! ideriv = 0 -> only calculate the total current constraint;
+! ideriv = 1 -> calculate the toroidal flux constraint and its first derivatives;
+! ideriv = 2 -> calculate the toroidal flux constraint and its first & second derivatives;
+!------------------------------------------------------------------------------------------------------ 
+ use globals, only: dp, zero, half, one, pi2, sqrtmachprec, ncpu, myid, ounit, &
+ coil, DoF, Ncoils, Cdof, Nfp, &
+ isum, t1I, Ndof, weight_isum, target_isum, &
+ iisum, misum, LM_fvec, LM_fjac, MPI_COMM_FOCUS
+ use mpi
+ implicit none
+
+ INTEGER, INTENT(in) :: ideriv
+ !--------------------------------------------------------------------------------------------
+ INTEGER :: astat, ierr
+ INTEGER :: icoil, idof, ND
+ REAL :: total_current
+ REAL, dimension(1:Ndof) :: current_array, symmetry
+ !--------------------------initialize and allocate arrays------------------------------------- 
+ isum = zero
+ total_current = zero
+ current_array = zero
+ symmetry = zero
+ !--------------------------function only ------------------------------------------------------ 
+ if( ideriv >= 0 ) then
+ do icoil = 1, Ncoils 
+ if(coil(icoil)%Ic /= 0) then ! exclude fixed currents
+ ! check if the coil is stellarator symmetric
+ select case (coil(icoil)%symm) 
+ case ( 0 )
+ symmetry(icoil) = 1
+ case ( 1 )
+ symmetry(icoil) = Nfp
+ case ( 2) 
+ symmetry(icoil) = 2*Nfp
+ end select
+ !if( myid.ne.modulo(icoil-1,ncpu) ) cycle ! parallelization loop;
+ current_array(icoil) = coil(icoil)%I
+ endif
+ enddo
+ total_current = SUM(current_array * symmetry)
+ isum = ((total_current - target_isum) / target_isum)**2
+
+ if (misum > 0) then ! L-M format of targets
+ LM_fvec(iisum+1) = weight_isum * ((total_current - target_isum) / target_isum)
+ endif
+ endif 
+
+ !--------------------------function & 1st deriv. -----------------------------------------------
+
+ if ( ideriv >= 1 ) then
+ t1I = zero
+ idof = 0 
+ do icoil = 1, Ncoils
+ ND = DoF(icoil)%ND
+ if ( coil(icoil)%Ic /= 0 ) then !if current is free;
+ idof = idof +1
+ t1I(idof) = 2 * ((total_current - target_isum) / target_isum) * symmetry(icoil) / target_isum
+
+ if (misum > 0) then ! L-M format of targets
+ LM_fjac(iisum + 1, idof) = weight_isum * symmetry(icoil) / target_isum
+ endif
+ endif
+ if ( coil(icoil)%Lc /= 0 ) then !if geometry is free;
+ idof = idof + ND
+ endif
+
+ enddo ! end icoil;
+ FATAL( current , idof .ne. Ndof, counting error in packing )
+ endif
+
+ !--------------------------------------------------------------------------------------------
+
+ call MPI_barrier( MPI_COMM_FOCUS, ierr )
+
+ return
+end subroutine current
+
+!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!

sources/datalloc.f90

@@ -208,7 +208,7 @@ subroutine AllocData(type)
 
  SALLOCATE( xdof, (1:Ndof), zero ) ! dof vector;
  SALLOCATE( dofnorm, (1:Ndof), one ) ! dof normalized value vector;
- SALLOCATE( evolution, (1:Nouts+1, 0:13), zero ) !evolution array;
+ SALLOCATE( evolution, (1:Nouts+1, 0:14), zero ) !evolution array;
  if (Tdof >= 1) then
  SALLOCATE( coilspace, (1:Nouts+1, 1:Tdof), zero ) ! all the coil parameters;
  endif 
@@ -360,6 +360,11 @@ subroutine AllocData(type)
  SALLOCATE( t1F, (1:Ndof), zero )
  endif
 
+ ! isum needed;
+ if (weight_isum > sqrtmachprec) then
+ SALLOCATE( t1I, (1:Ndof), zero )
+ endif
+
  ! ttlen needed;
  if (weight_ttlen > sqrtmachprec) then
  SALLOCATE( t1L, (1:Ndof), zero )
@@ -420,6 +425,12 @@ subroutine AllocData(type)
  mtflux = Nzeta
  LM_mfvec = LM_mfvec + mtflux
  endif
+
+ if (weight_isum > sqrtmachprec) then
+ iisum = LM_mfvec
+ misum = 1 ! put together
+ LM_mfvec = LM_mfvec + misum
+ endif 
 
  if (weight_ttlen > sqrtmachprec) then
  ittlen = LM_mfvec

sources/diagnos.f90

@@ -36,6 +36,7 @@ SUBROUTINE diagnos
  write(ounit, 100) "Bnormal", bnorm
  write(ounit, 100) "B_mn harmonics", bharm
  write(ounit, 100) "Toroidal flux", tflux
+ write(ounit, 100) "Total current penalty", tflux
  write(ounit, 100) "Coil length", ttlen
  write(ounit, 100) "Coil curvature", curv
  write(ounit, 100) "Coil torsion", tors

sources/globals.f90

@@ -16,7 +16,7 @@ module globals
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
 
- CHARACTER(10), parameter :: version='v0.17.01' ! version number
+ CHARACTER(10), parameter :: version='v0.18.00' ! version number
 
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
@@ -113,6 +113,9 @@ module globals
  ! toroidal flux
  REAL :: weight_tflux = 0.000D+00
  REAL :: target_tflux = 0.000D+00
+ ! total coil current
+ REAL :: weight_isum = 0.000D+00
+ REAL :: target_isum = 1.000D+06 
  ! coil length
  REAL :: weight_ttlen = 0.000D+00
  INTEGER :: case_length = 1
@@ -269,6 +272,8 @@ module globals
  input_harm ,&
  weight_tflux ,&
  target_tflux ,&
+ weight_isum ,&
+ target_isum ,&
  weight_ttlen ,&
  case_length ,&
  target_length ,&
@@ -447,8 +452,12 @@ module globals
 !latex \subsection{Optimization}
  ! General target functions;
  INTEGER :: iout, Nouts, LM_iter, LM_mfvec
- INTEGER :: ibnorm = 0, ibharm = 0, itflux = 0, ittlen = 0, icssep = 0, icurv = 0, istr=0, iccsep = 0, itors = 0, inissin = 0 ! starting number
- INTEGER :: mbnorm = 0, mbharm = 0, mtflux = 0, mttlen = 0, mcssep = 0, mcurv = 0, mstr=0, mccsep = 0, mtors = 0, mnissin = 0 ! numbers of targets
+ INTEGER :: ibnorm = 0, ibharm = 0, itflux = 0, iisum = 0, ittlen = 0, & 
+ & icssep = 0, icurv = 0, istr=0, iccsep = 0, itors = 0, &
+ & inissin = 0 ! starting number
+ INTEGER :: mbnorm = 0, mbharm = 0, mtflux = 0, misum = 0, mttlen = 0, &
+ & mcssep = 0, mcurv = 0, mstr=0, mccsep = 0, mtors = 0, &
+ & mnissin = 0 ! numbers of targets
  REAL :: chi, discretefactor, sumDE
  REAL , allocatable :: t1E(:), t2E(:,:), evolution(:,:), coilspace(:,:), deriv(:,:)
  REAL , allocatable :: LM_fvec(:), LM_fjac(:,:)
@@ -466,6 +475,9 @@ module globals
  INTEGER :: tflux_sign = -1 ! default theta : counter-clockwise
  REAL :: tflux, psi_avg
  REAL , allocatable :: t1F(:), t2F(:,:)
+ ! total current
+ REAL :: isum
+ REAL , allocatable :: t1I(:)
  ! Length constraint
  REAL :: ttlen
  REAL , allocatable :: t1L(:), t2L(:,:)

sources/initial.f90

@@ -668,11 +668,10 @@ subroutine check_input
  FATAL( initial, .true., selected case_stright is not supported )
  end select
 
-
-
  FATAL( initial, weight_bnorm < zero, illegal )
  FATAL( initial, weight_bharm < zero, illegal )
  FATAL( initial, weight_tflux < zero, illegal )
+ FATAL( initial, weight_isum < zero, illegal )
  FATAL( initial, weight_ttlen < zero, illegal )
  FATAL( initial, weight_specw < zero, illegal )
  FATAL( initial, weight_ccsep < zero, illegal )

sources/saving.f90

@@ -102,6 +102,8 @@ subroutine saving
  HWRITERV( 1 , weight_bharm , weight_bharm )
  HWRITERV( 1 , weight_tflux , weight_tflux )
  HWRITERV( 1 , target_tflux , target_tflux )
+ HWRITERV( 1 , weight_isum , weight_isum )
+ HWRITERV( 1 , target_isum , target_isum )
  HWRITERV( 1 , weight_ttlen , weight_ttlen )
  HWRITERV( 1 , target_length , target_length )
  HWRITERV( 1 , curv_k0 , curv_k0 ) 
@@ -289,7 +291,7 @@ subroutine saving
  HWRITERV( 1 , Gnorm , Gnorm )
  HWRITERV( 1 , Mnorm , Mnorm )
  HWRITERV( 1 , overlap , overlap )
- HWRITERA( iout, 14 , evolution , evolution(1:iout, 0:13) )
+ HWRITERA( iout, 15 , evolution , evolution(1:iout, 0:13) )
  if (allocated(coilspace)) then
  HWRITERA( iout, Tdof , coilspace , coilspace(1:iout, 1:Tdof) )
  endif
@@ -319,7 +321,9 @@ subroutine saving
  HWRITEIV( 1 , ibharm , ibharm )
  HWRITEIV( 1 , mbharm , mbharm ) 
  HWRITEIV( 1 , itflux , itflux )
- HWRITEIV( 1 , mtflux , mtflux ) 
+ HWRITEIV( 1 , mtflux , mtflux ) 
+ HWRITEIV( 1 , iisum , iisum )
+ HWRITEIV( 1 , misum , misum ) 
  HWRITEIV( 1 , ittlen , ittlen )
  HWRITEIV( 1 , mttlen , mttlen ) 
  HWRITEIV( 1 , icssep , icssep )