Merge pull request #16 from PrincetonUniversity/develop

Develop

New/Makefile

@@ -20,7 +20,7 @@ ifeq ($(CC),gfortran)
  RFLAGS=-O3 -Wall -fdefault-real-8 -ffixed-line-length-none -march=native -ffast-math
  DFLAGS=-g3 -Wextra -Wtarget-lifetime -fbacktrace -fbounds-check -ffpe-trap=zero -fcheck=all -DDEBUG
 else
- RFLAGS=-r8 -mcmodel=large -O3 -m64 -unroll0 -fno-alias -ip -traceback -vec_report0 #-ipo -xhost
+ RFLAGS=-r8 -mcmodel=large -O3 -m64 -unroll0 -fno-alias -ip -traceback #-vec_report0 #-ipo -xhost
  DFLAGS=-check all -check noarg_temp_created -debug full -D DEBUG
 endif
 
@@ -61,7 +61,7 @@ $(FFILES): %.F90: %.h
 ############################################################################################################
 
 clean:
- rm -f $(FFILES) ; rm -f $(ROBJS) ; rm -f $(DOBJS) ; rm -f *.mod ; rm -f *.pdf
+ rm -f $(FFILES) ; rm -f $(ROBJS) ; rm -f $(DOBJS) ; rm -f $(NUMOBJ) ; rm -f *.mod ; rm -f *.pdf
 
 ############################################################################################################
 

New/congrad.h

@@ -72,7 +72,7 @@ SUBROUTINE congrad
  if(myid .eq. 0) write(ounit, '("congrad : EXITING--------Local minimum reached!")')
  exit ! reach minimum
  endif
- 
+
  beta = sum(gradf**2) / sum( (gradf-gradk)*p )
  p = -gradf + beta*p ! direction for next step;
  gradk = gradf !save for the current step;

New/initial.h

@@ -420,7 +420,7 @@ subroutine initial
  FATAL( initial, TN_maxiter < 0, must be non-negative )
  endif
 
- write(ounit, '(8X,": IsNormalize = "I" ; IsNormWeight = "I)') IsNormalize, IsNormWeight
+ write(ounit, '(8X,": IsNormalize = "I1" ; IsNormWeight = "I1)') IsNormalize, IsNormWeight
 
  select case ( case_bnormal )
  case ( 0 )

New/rdcoils.h

@@ -84,7 +84,7 @@ subroutine rdcoils
  include "mpif.h"
 
  LOGICAL :: exist
- INTEGER :: icoil, maxnseg, ifirst, NF, itmp, ip, icoef
+ INTEGER :: icoil, maxnseg, ifirst, NF, itmp, ip, icoef, total_coef
  REAL :: Rmaj, zeta, totalcurrent, z0, r1, r2, z1, z2
  !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
@@ -383,17 +383,20 @@ subroutine rdcoils
  if (IsNormalize > 0) then
  Gnorm = 0
  Inorm = 0
+ total_coef = 0 ! total number of coefficients
  do icoil = 1, Ncoils
+ NF = FouCoil(icoil)%NF
+ total_coef = total_coef + (6*NF + 3)
  do icoef = 0, NF
  Gnorm = Gnorm + FouCoil(icoil)%xs(icoef)**2 + FouCoil(icoil)%xc(icoef)**2
  Gnorm = Gnorm + FouCoil(icoil)%ys(icoef)**2 + FouCoil(icoil)%yc(icoef)**2
  Gnorm = Gnorm + FouCoil(icoil)%zs(icoef)**2 + FouCoil(icoil)%zc(icoef)**2
  enddo
  Inorm = Inorm + coil(icoil)%I**2
  enddo
- Gnorm = sqrt(Gnorm) * weight_gnorm
- Inorm = sqrt(Inorm) * weight_inorm
- Inorm = Inorm * (NF + 1) * 6 ! compensate for the fact that there are so many more spatial variables
+ Gnorm = sqrt(Gnorm/total_coef) * weight_gnorm ! quadratic mean
+ Inorm = sqrt(Inorm/Ncoils) * weight_inorm ! quadratic mean
+ !Inorm = Inorm * 6 ! compensate for the fact that there are so many more spatial variables
 
  FATAL( rdcoils, abs(Gnorm) < machprec, cannot be zero )
  FATAL( rdcoils, abs(Inorm) < machprec, cannot be zero )

New/saving.h

@@ -29,7 +29,7 @@ subroutine saving
 
 
  INTEGER :: ii, jj, icoil, NF
- CHARACTER(LEN=10) :: version='v0.1.0'
+ CHARACTER(LEN=10) :: version='v0.1.02'
 
 
  ! the following are used by the macros HWRITEXX below; do not alter/remove;
@@ -209,10 +209,10 @@ subroutine saving
  if( save_coils == 1 ) then
 
  open(funit,file=trim(outcoils), status="unknown", form="formatted" )
- write(funit,'("periods 1")')
+ write(funit,'("periods "I3)') Nfp
  write(funit,'("begin filament")')
  write(funit,'("mirror NIL")')
- do icoil = 1, Ncoils
+ do icoil = 1, Ncoils*Npc
  do ii = 0, coil(icoil)%NS-1
  write(funit,1010) coil(icoil)%xx(ii), coil(icoil)%yy(ii), coil(icoil)%zz(ii), coil(icoil)%I
  enddo

New/solvers.h

@@ -311,13 +311,13 @@ end subroutine costfun
 subroutine normweight
  use globals, only : zero, one, machprec, ounit, myid, xdof, bnorm, bharm, tflux, ttlen, specw, ccsep, &
  weight_bnorm, weight_bharm, weight_tflux, weight_ttlen, weight_specw, weight_ccsep, target_tflux, &
- psi_avg, coil, Ncoils, case_length
+ psi_avg, coil, Ncoils, case_length, Bmnc, Bmns, tBmnc, tBmns
 
  implicit none 
  include "mpif.h"
 
  INTEGER :: ierr, icoil
- REAL :: tmp, cur_tflux
+ REAL :: tmp, cur_tflux, modBn, modtBn
 
  !----------------------------------------------------------------------------------------------------
 
@@ -345,23 +345,31 @@ subroutine normweight
 
  endif
 
- !-!-!-!-!-!-!-!-!-!-bnorm-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+ !-!-!-!-!-!-!-!-!-!-bharm-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
- if( weight_bnorm >= machprec ) then
+ if( weight_bharm >= machprec ) then
 
- call bnormal(0) 
- if (abs(bnorm) > machprec) weight_bnorm = weight_bnorm / bnorm
- if( myid == 0 ) write(ounit, 1000) "weight_bnorm", weight_bnorm
+ call bmnharm(0)
+ modBn = sqrt(sum(Bmnc**2 + Bmns**2))
+ modtBn = sqrt(sum(tBmnc**2 + tBmns**2))
+ do icoil = 1, Ncoils
+ coil(icoil)%I = coil(icoil)%I * modtBn / modBn
+ enddo
+ if(myid .eq. 0) write(ounit,'("solvers : rescale coil currents with a factor of "ES12.5)') &
+ modtBn / modBn
+ call bmnharm(0)
+ if (abs(bharm) > machprec) weight_bharm = weight_bharm / bharm
+ if( myid == 0 ) write(ounit, 1000) "weight_bharm", weight_bharm
 
  endif
 
  !-!-!-!-!-!-!-!-!-!-bnorm-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
- if( weight_bharm >= machprec ) then
+ if( weight_bnorm >= machprec ) then
 
- call bmnharm(0) 
- if (abs(bharm) > machprec) weight_bharm = weight_bharm / bharm
- if( myid == 0 ) write(ounit, 1000) "weight_bharm", weight_bharm
+ call bnormal(0) 
+ if (abs(bnorm) > machprec) weight_bnorm = weight_bnorm / bnorm
+ if( myid == 0 ) write(ounit, 1000) "weight_bnorm", weight_bnorm
 
  endif
 
@@ -403,7 +411,7 @@ subroutine normweight
 !!$ 
 !!$ endif
 
-1000 format("solvers : "A12" is normalized to" ES23.15)
+1000 format("solvers : "A12" is normalized to " ES23.15)
 
  call packdof(xdof)
 

New/torflux.h

@@ -246,9 +246,12 @@ subroutine bpotential0(icoil, iteta, jzeta, Ax, Ay, Az)
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
- FATAL( bpotential0, icoil .lt. 1 .or. icoil .gt. Ncoils*Npc, icoil not in right range )
- FATAL( bpotential0, iteta .lt. 0 .or. iteta .gt. Nteta , iteta not in right range )
- FATAL( bpotential0, jzeta .lt. 0 .or. jzeta .gt. Nzeta , jzeta not in right range )
+ FATAL( bpotential0, icoil .lt. 1 .or. icoil .gt. Ncoils*Npc, &
+ icoil not in right range )
+ FATAL( bpotential0, iteta .lt. 0 .or. iteta .gt. Nteta , &
+ iteta not in right range )
+ FATAL( bpotential0, jzeta .lt. 0 .or. jzeta .gt. Nzeta , &
+ jzeta not in right range )
 
  dlx = zero; ltx = zero; Ax = zero
  dly = zero; lty = zero; Ay = zero
@@ -300,9 +303,12 @@ subroutine bpotential1(icoil, iteta, jzeta, Ax, Ay, Az, ND)
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
- FATAL( bpotential1, icoil .lt. 1 .or. icoil .gt. Ncoils*Npc, icoil not in right range )
- FATAL( bpotential1, iteta .lt. 0 .or. iteta .gt. Nteta , iteta not in right range )
- FATAL( bpotential1, jzeta .lt. 0 .or. jzeta .gt. Nzeta , jzeta not in right range )
+ FATAL( bpotential1, icoil .lt. 1 .or. icoil .gt. Ncoils*Npc, &
+ icoil not in right range )
+ FATAL( bpotential1, iteta .lt. 0 .or. iteta .gt. Nteta , &
+ iteta not in right range )
+ FATAL( bpotential1, jzeta .lt. 0 .or. jzeta .gt. Nzeta , &
+ jzeta not in right range )
  FATAL( bpotential1, ND <= 0, wrong inout dimension of ND )
 
  NS = coil(icoil)%NS

Old/Makefile

@@ -1,10 +1,14 @@
 #!/bin/sh
 
 ############################################################################################################
+# comment out on 20180228 for NAG incompative
+# ALLFILES= globals numrec initial surface rdknot rdcoils knotxx iccoil bfield bnormal torflux tlength equarcl \
+# specwid coilsep nlinear fdcheck denergy descent congrad truncnt restart diagnos exinput identfy \
+# bnftran hessian focus
 
- ALLFILES= globals numrec initial surface rdknot rdcoils knotxx iccoil bfield bnormal torflux tlength equarcl \
+ ALLFILES= globals numrec initial surface rdcoils iccoil bfield bnormal torflux tlength equarcl \
  specwid coilsep nlinear fdcheck denergy descent congrad truncnt restart diagnos exinput identfy \
- bnftran focus
+ bnftran hessian focus
 
 ############################################################################################################
 
@@ -19,12 +23,17 @@
  MACROS=macros
  CC=intel # if want to use gfortran; make CC=gfortran
  FC=mpif90
- FLAGS=-r8 -mcmodel=large -O3 -m64 -unroll0 -fno-alias -ip -traceback -vec_report0
+ FLAGS=-r8 -mcmodel=large -O3 -m64 -unroll0 -fno-alias -ip -traceback -D BNORM #-vec_report0
  DFLAGS=-check all -check noarg_temp_created -debug full -D DEBUG
 
 ############################################################################################################
 
- NAG=-L$(NAG_ROOT)/lib -lnag_nag
+ #NAG=-L$(NAG_ROOT)/lib -lnag
+ NAG= -I$(NAG_ROOT)/nag_interface_blocks $(NAG_ROOT)/lib/libnag_mkl.a \
+ -Wl,--start-group $(NAG_ROOT)/mkl_intel64_11.3.3/lib/libmkl_intel_lp64.a \
+ $(NAG_ROOT)/mkl_intel64_11.3.3/lib/libmkl_intel_thread.a \
+ $(NAG_ROOT)/mkl_intel64_11.3.3/lib/libmkl_core.a -Wl,--end-group \
+ -liomp5 -lpthread -lm -ldl -lstdc++
 
  HDF5=-I$(HDF5_HOME)/include -L$(HDF5_HOME)/lib -lhdf5hl_fortran -lhdf5_hl -lhdf5_fortran \
 -lhdf5 -lpthread -lz -lm
@@ -39,12 +48,12 @@
 
 ############################################################################################################
 
-xfocus: $(OBJS) tnpack.o hybrj.o
- $(FC) -o xfocus $(OBJS) tnpack.o hybrj.o $(OCULUSFLAGS) $(NAG) $(HDF5)
+xfocus: $(OBJS) tnpack.o hybrj.o modchl.o
+ $(FC) -o xfocus $(OBJS) tnpack.o hybrj.o modchl.o $(OCULUSFLAGS) $(NAG) $(HDF5) $(MKL)
  @echo "Compiling xfocus finished."
 
-dfocus: $(DOBJS) tnpack.o hybrj.o
- $(FC) -o dfocus $(DOBJS) tnpack.o hybrj.o $(OCULUSFLAGS) $(NAG) $(HDF5)
+dfocus: $(DOBJS) tnpack.o hybrj.o modchl.o
+ $(FC) -o dfocus $(DOBJS) tnpack.o hybrj.o modchl.o $(OCULUSFLAGS) $(NAG) $(HDF5) $(MKL)
  @echo "Compiling dfocus finished."
 
 ############################################################################################################
@@ -54,11 +63,14 @@ hybrj.o: hybrj.f
 tnpack.o: tnpack.f
  $(FC) -c $(FLAGS) $(DFLAGS) -o $@ $<
 
+modchl.o: modchl.f
+ $(FC) -c $(FLAGS) $(DFLAGS) -o $@ $<
+
 $(OBJS): %_r.o: %.F90
- $(FC) -c $(FLAGS) -o $@ $< $(OCULUSFLAGS) $(NAG) $(HDF5)
+ $(FC) -c $(FLAGS) -o $@ $< $(OCULUSFLAGS) $(NAG) $(HDF5) $(MKL)
 
 $(DOBJS): %_d.o: %.F90
- $(FC) -c $(FLAGS) $(DFLAGS) -o $@ $< $(OCULUSFLAGS) $(NAG) $(HDF5)
+ $(FC) -c $(FLAGS) $(DFLAGS) -o $@ $< $(OCULUSFLAGS) $(NAG) $(HDF5) $(MKL)
 
 $(FFILES): %.F90: %.h
  m4 -P $(MACROS) $< > $@

Old/bfield.h

@@ -15,68 +15,68 @@
 
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
-subroutine bfield( RpZ, itangent, dBRpZ, ibfield ) 
-
-! DATE : April 2016
-! using bs00aa in oculus ( actually NAG adaptive integration routines ) calculating magnetic field;
-! only be used at the beginning stage; may be deleted later.
-
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
-
- use kmodule, only : zero, one, pi2, vsmall, myid, ounit, Ncoils, coil, icoil, bsfield
-
- use oculus , only : bs00aa
-
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
-
- implicit none
-
- include "mpif.h"
-
- INTEGER :: itangent, ibfield
- REAL :: RpZ(1:3), dBRpZ(1:3,0:3)
-
- INTEGER :: ierr, ibs00aa
- REAL :: zeta, Bx, By, Bz, czeta, szeta
-
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
-
- dBRpZ(1:3,0:3) = zero ! set default intent(out) ; 11 Oct 15;
-
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
-
- zeta = modulo( RpZ(2), pi2)
-
- czeta = cos(zeta)
- szeta = sin(zeta)
-
- bsfield%x = RpZ(1) * czeta
- bsfield%y = RpZ(1) * szeta
- bsfield%z = RpZ(3)
-
- Bx = zero
- By = zero
- Bz = zero
-
- do icoil = 1, Ncoils ! icoil is a global variable which is passed through to auxiliary.h; 11 Oct 15;
-
- ibs00aa = 0 ; bsfield%LB = .true. ; bsfield%LA = .false. ; bsfield%LL = .false. ; call bs00aa( bsfield, ibs00aa )
-
- Bx = Bx + bsfield%Bx * coil(icoil)%I
- By = By + bsfield%By * coil(icoil)%I
- Bz = Bz + bsfield%Bz * coil(icoil)%I
-
- enddo ! end of do icoil; 11 Oct 15;
-
- dBRpZ(1,0) = ( Bx * czeta + By * szeta )
- dBRpZ(2,0) = ( - Bx * szeta + By * czeta ) / RpZ(1)
- dBRpZ(3,0) = Bz
-
- ibfield = 0
-
- return
-
-end subroutine bfield
+!!$subroutine bfield( RpZ, itangent, dBRpZ, ibfield ) 
+!!$
+!!$! DATE : April 2016
+!!$! using bs00aa in oculus ( actually NAG adaptive integration routines ) calculating magnetic field;
+!!$! only be used at the beginning stage; may be deleted later.
+!!$
+!!$!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+!!$ 
+!!$ use kmodule, only : zero, one, pi2, vsmall, myid, ounit, Ncoils, coil, icoil, bsfield
+!!$ 
+!!$ use oculus , only : bs00aa
+!!$ 
+!!$!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+!!$ 
+!!$ implicit none
+!!$ 
+!!$ include "mpif.h"
+!!$ 
+!!$ INTEGER :: itangent, ibfield
+!!$ REAL :: RpZ(1:3), dBRpZ(1:3,0:3)
+!!$ 
+!!$ INTEGER :: ierr, ibs00aa
+!!$ REAL :: zeta, Bx, By, Bz, czeta, szeta
+!!$ 
+!!$!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+!!$ 
+!!$ dBRpZ(1:3,0:3) = zero ! set default intent(out) ; 11 Oct 15;
+!!$ 
+!!$!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+!!$ 
+!!$ zeta = modulo( RpZ(2), pi2)
+!!$ 
+!!$ czeta = cos(zeta)
+!!$ szeta = sin(zeta)
+!!$ 
+!!$ bsfield%x = RpZ(1) * czeta
+!!$ bsfield%y = RpZ(1) * szeta
+!!$ bsfield%z = RpZ(3)
+!!$ 
+!!$ Bx = zero
+!!$ By = zero
+!!$ Bz = zero
+!!$ 
+!!$ do icoil = 1, Ncoils ! icoil is a global variable which is passed through to auxiliary.h; 11 Oct 15;
+!!$ 
+!!$ ibs00aa = 0 ; bsfield%LB = .true. ; bsfield%LA = .false. ; bsfield%LL = .false. ; call bs00aa( bsfield, ibs00aa )
+!!$ 
+!!$ Bx = Bx + bsfield%Bx * coil(icoil)%I
+!!$ By = By + bsfield%By * coil(icoil)%I
+!!$ Bz = Bz + bsfield%Bz * coil(icoil)%I
+!!$ 
+!!$ enddo ! end of do icoil; 11 Oct 15;
+!!$ 
+!!$ dBRpZ(1,0) = ( Bx * czeta + By * szeta )
+!!$ dBRpZ(2,0) = ( - Bx * szeta + By * czeta ) / RpZ(1)
+!!$ dBRpZ(3,0) = Bz
+!!$ 
+!!$ ibfield = 0
+!!$ 
+!!$ return
+!!$ 
+!!$end subroutine bfield
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 ! DATE: 06/15/2016