(old) Bugfix in Igeom=1 Lcons=3 and output of stability variables

Utilities/pythontools/py_spec/ci/test.py

@@ -23,6 +23,8 @@ def compare(data, reference, localtol=1e-6, action='ERR'):
     """
     global match
     for key, value in vars(data).items():
+        if key in ['Lconstraint', 'LBeltrami']:
+            continue
         if isinstance(value, SPECout):  # recurse data (csmiet: I'm not the biggest fan of this recursion...)
             print('------------------')
             print('Elements in '+key)

ci/G1V03L2Fi/G1V03L2Fi.001.sp

@@ -63,7 +63,7 @@ Vns(0,2)    =  0.000000000000000E+00 Bns(0,2)    =  0.000000000000000E+00 Vnc(0,
  Mregular    =        -1
 /
 &locallist
- LBeltrami   =         2
+ LBeltrami   =         4
  Linitgues   =         1
 /
 &globallist

src/dforce.f90

@@ -98,14 +98,14 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)
 
   use numerical, only : logtolerance
 
-  use fileunits, only : ounit
+  use fileunits, only : ounit, munit
 
   use inputlist, only : Wmacros, Wdforce, Nvol, Ntor, Lrad, Igeometry, &
                         epsilon, &
                         Lconstraint, Lcheck, dRZ, &
                         Lextrap, &
                         mupftol, &
-                        Lfreebound, LHmatrix
+                        Lfreebound, LHmatrix, gamma, pscale, adiabatic
 
   use cputiming, only : Tdforce
 
@@ -134,7 +134,8 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)
                         LocalConstraint, xoffset, &
                         solution, IPdtdPf, &
                         IsMyVolume, IsMyVolumeValue, WhichCpuID, &
-                        ext ! For outputing Lcheck = 6 test
+                        ext, & ! For outputing Lcheck = 6 test
+                        BetaTotal, vvolume
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
@@ -164,6 +165,9 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)
 
   LOGICAL              :: LComputeAxis, dfp100_logical
 
+  REAL                 :: press, voltotal
+  REAL                 :: betavol(1:Mvol)
+
 #ifdef DEBUG
   INTEGER              :: isymdiff
   REAL                 :: dvol(-1:+1), evolume, imupf(1:2,-2:2), lfactor
@@ -239,8 +243,13 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)
       ! Mvol-1 surface current plus 1 poloidal linking current constraints
       Ndofgl = Mvol
     else
-      ! Mvol-1 surface current constraints
-      Ndofgl = Mvol-1
+      ! add an additional constraint to make the total pflux = 0 
+      if(Igeometry.eq.1) then
+        Ndofgl = Mvol
+      else
+        ! Mvol-1 surface current constraints
+        Ndofgl = Mvol-1
+      endif
     endif
 
     SALLOCATE( Fvec, (1:Ndofgl), zero )
@@ -257,6 +266,11 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)
 
         ! one step Newton's method
         dpflux(2:Mvol) = dpflux(2:Mvol) - dpfluxout(1:Mvol-1)
+
+        if(Igeometry.eq.1) then
+          dpflux(1) = dpflux(1) - dpfluxout(Mvol)
+        endif
+
         if( Lfreebound.eq.1 ) then
           dtflux(Mvol) = dtflux(Mvol  ) - dpfluxout(Mvol    )
         endif
@@ -305,6 +319,38 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)
 
     enddo ! end of do vvol = 1, Mvol
 
+    !add an additional constraint to make the total pflux 0
+    if(Igeometry.eq.1) then
+
+      vvol = 1
+
+      WCALL(dforce, IsMyVolume, (vvol))
+
+      if( IsMyVolumeValue .EQ. 0 ) then
+
+      else if( IsMyVolumeValue .EQ. -1) then
+          FATAL(dforce, .true., Unassociated volume)
+      else
+        NN = NAdof(vvol)
+
+        SALLOCATE( solution, (1:NN, 0:2), zero)
+
+        ! Pack field and its derivatives
+        packorunpack = 'P'
+        WCALL( dforce, packab, ( packorunpack, vvol, NN, solution(1:NN,0), 0 ) ) ! packing;
+        WCALL( dforce, packab, ( packorunpack, vvol, NN, solution(1:NN,2), 2 ) ) ! packing;
+
+        ! compute the field with renewed dpflux via single Newton method step
+        solution(1:NN, 0) = solution(1:NN, 0) - dpfluxout(Mvol) * solution(1:NN, 2)
+
+        ! Unpack field in vector potential Fourier harmonics
+        packorunpack = 'U'
+        WCALL( dforce, packab, ( packorunpack, vvol, NN, solution(1:NN,0), 0 ) ) ! unpacking;
+
+        DALLOCATE( solution )
+      endif
+    endif
+
     DALLOCATE(Fvec)
     DALLOCATE(dpfluxout)
 
@@ -393,6 +439,28 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)
 
   lBBintegral(1:Nvol) = lBBintegral(1:Nvol) * half
 
+  voltotal = 0.0
+
+  do vvol = 1, Mvol
+    WCALL(dforce, IsMyVolume, (vvol))
+    if( IsMyVolumeValue .EQ. 0 ) then
+      cycle
+    else if( IsMyVolumeValue .EQ. -1) then
+        FATAL(dforce, .true., Unassociated volume)
+    endif
+
+    vvolume(vvol) = vvolume(vvol)
+
+    press = adiabatic(vvol) * pscale / vvolume(vvol)**gamma
+
+    betavol(vvol) = press * vvolume(vvol) / lBBintegral(vvol)
+    betavol(vvol) = betavol(vvol) * vvolume(vvol)
+    voltotal = voltotal+vvolume(vvol)
+    !write(*,*) "Calc beta: ", betavol(vvol), vvolume(vvol)
+  enddo
+
+  BetaTotal = sum(betavol(1:Nvol))/voltotal ! Calculate total beta which is obtained from individual betas
+
   Energy = sum( lBBintegral(1:Nvol) ) ! should also compute beta;
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
@@ -901,6 +969,17 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)
   !call MPI_BARRIER( MPI_COMM_WORLD, ierr )
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+
+   if(Lhessianallocated .and. Igeometry.eq.1) then
+     if( myid.eq.0 ) then ; cput = GETTIME ; write(ounit,'("hesian : ",f10.2," : LHmatrix="L2" ;")')cput-cpus, LHmatrix ;
+       write(*,*) "Writing .hessian file..."
+       open(munit, file=trim(ext)//".sp.hessian", status="unknown", form="unformatted")
+       write(munit) NGdof
+       write(munit) hessian(1:NGdof,1:NGdof)
+       close(munit)
+
+     endif
+   endif
 
   RETURN(dforce)
 

src/dfp100.f90

@@ -63,7 +63,7 @@ subroutine dfp100(Ndofgl, x, Fvec, LComputeDerivatives)
                         DDtzcc, DDtzcs, DDtzsc, DDtzss, &
                         DDzzcc, DDzzcs, DDzzsc, DDzzss, &
                         dMA, dMB, dMD, dMG, MBpsi, solution, &
-                        Nt, Nz, LILUprecond, Lsavedguvij, NOTMatrixFree, guvijsave, izbs
+                        Nt, Nz, LILUprecond, Lsavedguvij, NOTMatrixFree, guvijsave, izbs, total_pflux
 
   LOCALS
   !------
@@ -205,6 +205,19 @@ subroutine dfp100(Ndofgl, x, Fvec, LComputeDerivatives)
             Fvec(1:Mvol-1) = IPDt(1:Mvol-1) - Isurf(1:Mvol-1)
         endif
 
+	! Set the total pflux to 0
+        if(Igeometry.eq.1) then
+
+          IPDtDpf(1,Mvol) = -pi2 * Bt00(1,1,2)
+          IPdtdPf(2:Mvol,Mvol) = 0.0
+          IPDtDpf(Mvol,1:Mvol) = 1.0
+
+          ! Constraint the total pflux (pflux(Mvol))
+          ! zero for symmetric initial profiles
+          ! non-zero for asymmetric profiles
+          Fvec(Mvol) = sum(dpflux(1:Mvol)) - total_pflux 
+        endif
+
         ! Compute poloidal linking current constraint as well in case of free boundary computation
         if ( Lfreebound.eq.1 ) then
 

src/global.f90

@@ -251,11 +251,14 @@ module allglobal
   REAL                 :: ForceErr !< total force-imbalance
   REAL                 :: Energy   !< MHD energy
   REAL                 :: BnsErr   !< (in freeboundary) error in self-consistency of field on plasma boundary (Picard iteration)
+  REAL                 :: BetaTotal = 0.0   !< Beta, averaged over entire domain
 
   REAL   , allocatable :: IPDt(:), IPDtDpf(:,:)  !< Toroidal pressure-driven current
 
   INTEGER              :: Mvol
 
+  REAL :: total_pflux ! used when Lconstraint=3, Igeometry=1
+
   LOGICAL              :: YESstellsym !< internal shorthand copies of Istellsym, which is an integer input;
   LOGICAL              :: NOTstellsym !< internal shorthand copies of Istellsym, which is an integer input;
 
@@ -1675,7 +1678,7 @@ subroutine wrtend
   write(iunit,'(" adiabatic   = ",257es23.15)') adiabatic(1:Mvol)
   write(iunit,'(" mu          = ",257es23.15)') mu(1:Mvol)
   write(iunit,'(" Ivolume     = ",257es23.15)') Ivolume(1:Mvol)
-  write(iunit,'(" Isurf       = ",257es23.15)') Isurf(1:Mvol-1), 0.0
+  write(iunit,'(" Isurf       = ",257es23.15)') IPDt(1:Mvol) ! Prints the actual surf current, not the targeted one
   write(iunit,'(" Lconstraint = ",i9        )') Lconstraint
   write(iunit,'(" pl          = ",257i23    )') pl(0:Mvol)
   write(iunit,'(" ql          = ",257i23    )') ql(0:Mvol)

src/hesian.f90

@@ -22,7 +22,7 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof )
   use inputlist, only : Wmacros, Whesian, Igeometry, Nvol, pflux, helicity, mu, Lfreebound, &
                         LHevalues, LHevectors, LHmatrix, &
                         Lperturbed, dpp, dqq, &
-                        Lcheck, Lfindzero
+                        Lcheck, Lfindzero, Lconstraint
 
   use cputiming, only : Thesian
 
@@ -37,7 +37,9 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof )
                         Lhessian3Dallocated,denergydrr, denergydrz,denergydzr,denergydzz, &
 					            	LocalConstraint
 
-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
+  use sphdf5, only : write_stability
+
+  !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
   LOCALS
 
@@ -51,7 +53,7 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof )
   REAL                :: xx(0:NGdof,-2:2), ff(0:NGdof,-2:2), df(1:NGdof)!, deriv
 
   INTEGER             :: vvol, idof, ii, mi, ni, irz, issym, isymdiff, lvol, ieval(1:1), igdof, ifd
-  REAL                :: oldEnergy(-2:2), error, cpul
+  REAL                :: oldEnergy(-2:2), error, cpul 
 
   REAL                :: oldBB(1:Mvol,-2:2), oBBdRZ(1:Mvol,0:1,1:LGdof), ohessian(1:NGdof,1:NGdof)
 
@@ -90,6 +92,11 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof )
 
   BEGIN(hesian)
 
+  ! Only makes sense to compute the Hessian if helicity is constrained
+  if(Lconstraint.ne.2) then
+    return
+  endif
+
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
   lmu(1:Nvol) = mu(1:Nvol) ; lpflux(1:Nvol) = pflux(1:Nvol) ; lhelicity(1:Nvol) = helicity(1:Nvol) ! save original profile information; 20 Jun 14;
@@ -229,7 +236,7 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof )
   !if (LHmatrix) then
     Lhessian3Dallocated = .true.
   !else
-   ! Lhessianallocated = .true.
+  Lhessianallocated = .true.
   !endif
    !This step cleared.
 
@@ -585,6 +592,10 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof )
 
   endif ! end of if(                 ( LHevalues .or. LHevectors ) )
 
+  ! output hessian, eigenvalues, and eigenvectors to the .h5 file
+  if( LHmatrix .and. Lconstraint.eq.2) then
+    WCALL( hesian, write_stability, (ohessian, evalr, evali, evecr, NGdof) )
+  endif
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
   if( myid.eq.0 .and. Lperturbed.eq.1 ) then

src/newton.f90

@@ -97,7 +97,7 @@ subroutine newton( NGdof, position, ihybrd )
 
   INTEGER                :: ML, MU ! required for only Lc05ndf;
 
-  LOGICAL                :: Lexit = .true. ! perhaps this could be made user input;
+  LOGICAL                :: Lexit = .false. ! if false, hessian is printed even with initial force balance
   LOGICAL                :: LComputeAxis
 
   INTEGER                :: nprint = 1, nfev, njev

src/preset.f90

@@ -508,8 +508,14 @@ subroutine preset
   SALLOCATE( dtflux, (1:Mvol), zero )
   SALLOCATE( dpflux, (1:Mvol), zero )
 
+  ! with Igeom=1, Lcons=3, need to provide total pflux
+  if(Lconstraint.eq.3 .and. Igeometry.eq.1) then
+      total_pflux = pflux(Mvol) * phiedge / pi2
+      pflux(Mvol) = 0
+  endif
+
   select case( Igeometry )
-  case( 1   ) ; dtflux(1) = tflux(1) ; dpflux(1) = pflux(1) ! Cartesian              ; this is the "inverse" operation defined in xspech; 09 Mar 17;
+  case( 1   ) ; dtflux(1) = tflux(1) ; dpflux(1) = pflux(1) ! Cartesian ; this is the "inverse" operation defined in xspech; 09 Mar 17;
   case( 2:3 ) ; dtflux(1) = tflux(1) ; dpflux(1) = zero     ! cylindrical or toroidal;
   end select
 
@@ -807,7 +813,12 @@ subroutine preset
   if( Lfreebound.eq.1 ) then
     SALLOCATE( IPDtDpf, (1:Mvol  , 1:Mvol  ), zero)
   else
-    SALLOCATE( IPDtDpf, (1:Mvol-1, 1:Mvol-1), zero)
+    if(Igeometry.eq.1) then
+      ! add an additional constraint to make the total pflux = 0
+      SALLOCATE( IPDtDpf, (1:Mvol, 1:Mvol), zero) 
+    else
+      SALLOCATE( IPDtDpf, (1:Mvol-1, 1:Mvol-1), zero)
+    endif
   endif
 
   SALLOCATE( cheby, (0:Mrad,0:2), zero )