Derivatives of stochastic optimization are correct, new save function…

…ality implemented, more post-processing implemented

sources/Makefile

@@ -2,10 +2,10 @@
 
 ############################################################################################################
 
- ALLFILES= globals initial stochastic surface rdsurf rdknot rdcoils packdof bfield \
+ ALLFILES= globals initial saving diagnos stochastic surface rdsurf rdknot rdcoils packdof bfield \
  fdcheck datalloc solvers descent congrad lmalg \
  bmnharm bnormal torflux length surfsep nissin torsion coilsep curvature \
- saving diagnos specinp poinplot boozer wtmgrid focus 
+ specinp poinplot boozer wtmgrid focus 
  HFILES= $(ALLFILES:=.f90) # raw source files
  FFILES= $(ALLFILES:=_m.F90) # Fortran 90 files 
  PFILES= $(ALLFILES:=.pdf) # documentations

sources/bfield.f90

@@ -27,7 +27,7 @@ subroutine bfield0(icoil, x, y, z, tBx, tBy, tBz)
 ! Biot-Savart constant and currents are not included for later simplication. 
 ! Be careful if coils have different resolutions.
 !------------------------------------------------------------------------------------------------------ 
- use globals, only: dp, coil, surf, Ncoils, Nteta, Nzeta, cosnfp, sinnfp, &
+ use globals, only: dp, coil, surf, Ncoils, Nteta, Nzeta, cosnfp, sinnfp, machprec, &
  zero, myid, ounit, Nfp, pi2, half, two, one, bsconstant, MPI_COMM_FOCUS
  use mpi
  implicit none
@@ -89,6 +89,7 @@ subroutine bfield0(icoil, x, y, z, tBx, tBy, tBz)
  dlx = xx - coil(icoil)%xx(kseg)
  dly = yy - coil(icoil)%yy(kseg)
  dlz = zz - coil(icoil)%zz(kseg)
+ if ( dlx**2+dly**2+dlz**2 .lt. machprec ) cycle
  rm3 = (sqrt(dlx**2 + dly**2 + dlz**2))**(-3)
  ltx = coil(icoil)%xt(kseg)
  lty = coil(icoil)%yt(kseg)

sources/bnormal.f90

@@ -31,16 +31,17 @@ subroutine bnormal( ideriv )
 
  INTEGER, INTENT(in) :: ideriv
  !--------------------------------------------------------------------------------------------
- INTEGER :: astat, ierr
+ INTEGER :: astat, ierr, request
  INTEGER :: icoil, iteta, jzeta, idof, ND, NumGrid, isurf
  !--------------------------initialize and allocate arrays-------------------------------------
  isurf = plasma 
  NumGrid = Nteta*Nzeta
  ! reset to zero;
- bnorm = zero 
+ bnorm = zero
  surf(isurf)%Bx = zero; surf(isurf)%By = zero; surf(isurf)%Bz = zero; surf(isurf)%Bn = zero 
  dBx = zero; dBy = zero; dBz = zero; Bm = zero
  bn = zero
+ call MPI_BARRIER( MPI_COMM_FOCUS, ierr )
  !-------------------------------calculate Bn-------------------------------------------------- 
  if( ideriv >= 0 ) then
  do jzeta = 0, Nzeta - 1
@@ -78,6 +79,7 @@ subroutine bnormal( ideriv )
  call MPI_ALLREDUCE( MPI_IN_PLACE, surf(isurf)%Bz, NumGrid, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FOCUS, ierr )
  call MPI_ALLREDUCE( MPI_IN_PLACE, surf(isurf)%Bn, NumGrid, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FOCUS, ierr )
  call MPI_ALLREDUCE( MPI_IN_PLACE, bnorm, 1 , MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FOCUS, ierr )
+ call MPI_BARRIER( MPI_COMM_FOCUS, ierr )
  bnorm = bnorm * half * discretefactor
  bn = surf(isurf)%Bn + surf(isurf)%pb ! bn is B.n from coils
  ! bn = surf(isurf)%Bx * surf(isurf)%nx + surf(isurf)%By * surf(isurf)%ny + surf(isurf)%Bz * surf(isurf)%nz
@@ -198,6 +200,6 @@ subroutine bnormal( ideriv )
  endif
  endif
  !--------------------------------------------------------------------------------------------
- call MPI_barrier( MPI_COMM_FOCUS, ierr )
+ call MPI_BARRIER( MPI_COMM_FOCUS, ierr )
  return
 end subroutine bnormal

sources/curvature.f90

@@ -303,6 +303,13 @@ subroutine CurvDeriv1(icoil, derivs, ND, NF) !Calculate all derivatives for a co
 
  derivs = pi2*derivs/(NS*leng)
 
+ DALLOCATE(dxtdDoF)
+ DALLOCATE(dytdDoF)
+ DALLOCATE(dztdDoF)
+ DALLOCATE(dxadDoF)
+ DALLOCATE(dyadDoF)
+ DALLOCATE(dzadDoF)
+
  return
 
 end subroutine CurvDeriv1

sources/datalloc.f90

@@ -13,7 +13,7 @@ subroutine AllocData(type)
 
  INTEGER, intent(in) :: type
 
- INTEGER :: icoil, idof, ND, NF, icur, imag, isurf, NS, mm, iseg
+ INTEGER :: icoil, idof, ND, NF, icur, imag, isurf, NS, mm, iseg, n
  REAL :: xtmp, mtmp, tt
 
  isurf = plasma
@@ -33,6 +33,9 @@ subroutine AllocData(type)
  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)
+ SALLOCATE(DoF(icoil)%xtof, (0:coil(icoil)%NS, 1:ND), zero)
+ SALLOCATE(DoF(icoil)%ytof, (0:coil(icoil)%NS, 1:ND), zero)
+ SALLOCATE(DoF(icoil)%ztof, (0:coil(icoil)%NS, 1:ND), zero)
  ! allocate and calculate trignometric functions for re-use 
  SALLOCATE( FouCoil(icoil)%cmt, (0:NS, 0:NF), zero )
  SALLOCATE( FouCoil(icoil)%smt, (0:NS, 0:NF), zero ) 
@@ -67,6 +70,16 @@ subroutine AllocData(type)
  DoF(icoil)%yof(0:NS-1, 3*NF+3:4*NF+2) = FouCoil(icoil)%smt(0:NS-1, 1:NF) !y/ys
  DoF(icoil)%zof(0:NS-1, 4*NF+3:5*NF+3) = FouCoil(icoil)%cmt(0:NS-1, 0:NF) !z/zc
  DoF(icoil)%zof(0:NS-1, 5*NF+4:6*NF+3) = FouCoil(icoil)%smt(0:NS-1, 1:NF) !z/zs
+
+ do n = 1,NF
+ DoF(icoil)%xtof(0:NS,n+1) = -1.0*FouCoil(icoil)%smt(0:NS,n) * real(n) !xt/xc
+ DoF(icoil)%xtof(0:NS,n+NF+1) = FouCoil(icoil)%cmt(0:NS,n) * real(n) !xt/xs
+ DoF(icoil)%ytof(0:NS,n+2*NF+2) = -1.0*FouCoil(icoil)%smt(0:NS,n) * real(n) !yt/yc
+ DoF(icoil)%ytof(0:NS,n+3*NF+2) = FouCoil(icoil)%cmt(0:NS,n) * real(n) !yt/ys
+ DoF(icoil)%ztof(0:NS,n+4*NF+3) = -1.0*FouCoil(icoil)%smt(0:NS,n) * real(n) !zt/zc
+ DoF(icoil)%ztof(0:NS,n+5*NF+3) = FouCoil(icoil)%cmt(0:NS,n) * real(n) !zt/zs
+ enddo
+
  ! allocate xyz data
  SALLOCATE( coil(icoil)%xx, (0:coil(icoil)%NS), zero )
  SALLOCATE( coil(icoil)%yy, (0:coil(icoil)%NS), zero )
@@ -82,7 +95,26 @@ subroutine AllocData(type)
  SALLOCATE( coil(icoil)%zb, (0:coil(icoil)%NS), zero )
  SALLOCATE( coil(icoil)%dl, (0:coil(icoil)%NS), zero )
  SALLOCATE( coil(icoil)%dd, (0:coil(icoil)%NS), zero )
- coil(icoil)%dd = pi2 / NS ! discretizing factor;
+ coil(icoil)%dd = pi2 / NS ! discretizing factor also set in rdcoils?
+
+ if ( filforce .eq. 1 ) then
+ SALLOCATE( coil(icoil)%Fx, (0:coil(icoil)%NS), 0.0 )
+ SALLOCATE( coil(icoil)%Fy, (0:coil(icoil)%NS), 0.0 )
+ SALLOCATE( coil(icoil)%Fz, (0:coil(icoil)%NS), 0.0 )
+ SALLOCATE( coil(icoil)%Bxx, (0:coil(icoil)%NS), 0.0 )
+ SALLOCATE( coil(icoil)%Byy, (0:coil(icoil)%NS), 0.0 )
+ SALLOCATE( coil(icoil)%Bzz, (0:coil(icoil)%NS), 0.0 )
+ endif
+
+ if ( calcfb .eq. 1 ) then
+ SALLOCATE( coil(icoil)%nfbx, (0:coil(icoil)%NS), 0.0 )
+ SALLOCATE( coil(icoil)%nfby, (0:coil(icoil)%NS), 0.0 )
+ SALLOCATE( coil(icoil)%nfbz, (0:coil(icoil)%NS), 0.0 )
+ SALLOCATE( coil(icoil)%bfbx, (0:coil(icoil)%NS), 0.0 )
+ SALLOCATE( coil(icoil)%bfby, (0:coil(icoil)%NS), 0.0 )
+ SALLOCATE( coil(icoil)%bfbz, (0:coil(icoil)%NS), 0.0 )
+ endif
+
  case(2)
 #ifdef dposition
  DoF(icoil)%ND = coil(icoil)%Lc * 5 ! number of DoF for permanent magnet
@@ -219,7 +251,7 @@ subroutine AllocData(type)
  if (type == 0 .or. type == 1) then ! 0-order cost functions related arrays;
 
  ! Bnorm and Bharm needed;
- if (weight_bnorm > sqrtmachprec .or. weight_bharm > sqrtmachprec .or. IsQuiet <= -2) then
+ if (weight_bnorm > sqrtmachprec .or. weight_bharm > sqrtmachprec .or. weight_sbnorm > sqrtmachprec .or. IsQuiet <= -2) then
  SALLOCATE( bn, (0:Nteta-1,0:Nzeta-1), zero ) ! Bn from coils; 
  SALLOCATE( surf(isurf)%bn, (0:Nteta-1,0:Nzeta-1), zero ) ! total Bn;
  SALLOCATE( surf(isurf)%Bx, (0:Nteta-1,0:Nzeta-1), zero ) ! Bx on the surface;
@@ -250,7 +282,7 @@ subroutine AllocData(type)
  SALLOCATE( deriv, (1:Ndof, 0:10), zero )
 
  ! Bnorm related;
- if (weight_bnorm > sqrtmachprec .or. weight_bharm > sqrtmachprec) then
+ if (weight_bnorm > sqrtmachprec .or. weight_bharm > sqrtmachprec .or. weight_sbnorm > sqrtmachprec ) then
  SALLOCATE( t1B, (1:Ndof), zero ) !total d bnorm / d x;
  SALLOCATE( dBn, (1:Ndof), zero ) !total d Bn / d x;
  SALLOCATE( dBm, (1:Ndof), zero ) !total d Bm / d x;