Bugfix in slab geometry when using Lconstraint=3, and output of stability variables

.github/workflows/build.yml

@@ -33,10 +33,16 @@ jobs:
         make BUILD_ENV=gfortran_ubuntu dspec
     - name: run_fast_cartesian
       run: |
-        cd ${SPEC_PATH}/ci/G1V03L2Fi
+        cd ${SPEC_PATH}/ci/G1V03L0Fi
         export OMP_NUM_THREADS=1
-        mpiexec -n 2 --allow-run-as-root ${SPEC_PATH}/xspec G1V03L2Fi.001.sp
-        python3 -m py_spec.ci.test compare.h5 G1V03L2Fi.001.sp.h5 
+        mpiexec -n 2 --allow-run-as-root ${SPEC_PATH}/xspec G1V03L0Fi.sp
+        python3 -m py_spec.ci.test compare.h5 G1V03L0Fi.sp.h5 
+    - name: run_fast_cartesian_lcons3
+      run: |
+        cd ${SPEC_PATH}/ci/G1V03L3Fi
+        export OMP_NUM_THREADS=1
+        mpiexec -n 2 --allow-run-as-root ${SPEC_PATH}/xspec G1V03L3Fi.sp
+        python3 -m py_spec.ci.test compare.h5 G1V03L3Fi.sp.h5 
     - name: run_fast_cylinder
       run: |
         cd ${SPEC_PATH}/ci/G2V32L1Fi

.github/workflows/build_cmake.yml

@@ -39,10 +39,16 @@ jobs:
     #    cmake --install .
     - name: run_fast_cartesian
       run: |
-        cd ${SPEC_PATH}/ci/G1V03L2Fi
+        cd ${SPEC_PATH}/ci/G1V03L0Fi
         export OMP_NUM_THREADS=1
-        mpiexec -n 2 --allow-run-as-root $SPEC_PATH/install/bin/xspec G1V03L2Fi.001.sp
-        python3 -m py_spec.ci.test compare.h5 G1V03L2Fi.001.sp.h5 
+        mpiexec -n 2 --allow-run-as-root $SPEC_PATH/install/bin/xspec G1V03L0Fi.sp
+        python3 -m py_spec.ci.test compare.h5 G1V03L0Fi.sp.h5 
+    - name: run_fast_cartesian_lcons3
+      run: |
+        cd ${SPEC_PATH}/ci/G1V03L3Fi
+        export OMP_NUM_THREADS=1
+        mpiexec -n 2 --allow-run-as-root $SPEC_PATH/install/bin/xspec G1V03L3Fi.sp
+        python3 -m py_spec.ci.test compare.h5 G1V03L3Fi.sp.h5 
     - name: run_fast_cylinder
       run: |
         cd ${SPEC_PATH}/ci/G2V32L1Fi

.gitlab-ci.yml

@@ -31,19 +31,38 @@ stable_force_free_sheet:
   image: csmiet/testspec
   script:
   - source source.sh
-  - cd ci/G1V03L2Fi
-  - xspec G1V03L2Fi.001.sp
-  - python3 -m py_spec.ci.test compare.h5 G1V03L2Fi.001.sp.h5 
+  - cd ci/G1V03L0Fi
+  - xspec G1V03L0Fi.sp
+  - python3 -m py_spec.ci.test compare.h5 G1V03L0Fi.sp.h5 
   #comparison python3 -m py_spec.ci.testspec f1 f2
   artifacts:
     when: always
     paths:
-    - ci/G1V03L2Fi/G1V03L2Fi.001.sp.h5
+    - ci/G1V03L0Fi/G1V03L0Fi.sp.h5
     expire_in: 1 week
   dependencies:
   - build:centos7
 # reproduce:
-# docker run -v `pwd`:`pwd` -w `pwd` csmiet/testspec sh -c 'source ./source.sh;cd ci/G1V03L2Fi; xspec G1V03L2Fi.001.sp; python3 -m py_spec.ci.test compare.h5 G1V03L2Fi.001.sp.h5'
+# docker run -v `pwd`:`pwd` -w `pwd` csmiet/testspec sh -c 'source ./source.sh;cd ci/G1V03L0Fi; xspec G1V03L0Fi.sp; python3 -m py_spec.ci.test compare.h5 G1V03L0Fi.sp.h5'
+
+stable_force_free_sheet_lcons3:
+  stage: test
+  image: csmiet/testspec
+  script:
+  - source source.sh
+  - cd ci/G1V03L3Fi
+  - xspec G1V03L3Fi.sp
+  - python3 -m py_spec.ci.test compare.h5 G1V03L3Fi.sp.h5 
+  #comparison python3 -m py_spec.ci.testspec f1 f2
+  artifacts:
+    when: always
+    paths:
+    - ci/G1V03L3Fi/G1V03L3Fi.sp.h5
+    expire_in: 1 week
+  dependencies:
+  - build:centos7
+# reproduce:
+# docker run -v `pwd`:`pwd` -w `pwd` csmiet/testspec sh -c 'source ./source.sh;cd ci/G1V03L3Fi; xspec G1V03L3Fi.sp; python3 -m py_spec.ci.test compare.h5 G1V03L3Fi.sp.h5'
 
 cylindrical_many_volumes:
   stage: test

ci/G1V03L3Fi/G1V03L3Fi.sp

@@ -0,0 +1,109 @@
+&physicslist
+ Igeometry   =         1
+ Istellsym   =         1
+ Lfreebound  =         0
+ phiedge     =   1.000000000000000E+00
+ curtor      =   0.000000000000000E+00
+ curpol      =   0.000000000000000E+00
+ gamma       =   0.000000000000000E+00
+ Nfp         =         1
+ Nvol        =         3
+ Mpol        =         1
+ Ntor        =         0
+ Lrad        =                       12                      12                      12
+ tflux       =   4.000000000000000E-01  6.000000000000001E-01  1.000000000000000E+00
+ pflux       =  -4.000000000000001E-02 -4.000000000000001E-02  0.000000000000000E+00
+ helicity    =  -8.673617379884035E-19  1.331565149512061E-03  1.734723475976807E-18
+ pscale      =   0.000000000000000E+00
+ Ladiabatic  =         0
+ pressure    =   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
+ adiabatic   =   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
+ mu          =   0.000000000000000E+00  1.000000000000000E-01  0.000000000000000E+00
+ Lconstraint =         3
+ pl          =                       0                      0                      0                      0
+ ql          =                       0                      0                      0                      0
+ pr          =                       0                      0                      0                      0
+ qr          =                       0                      0                      0                      0
+ iota        =  -1.618000000000000E+00 -1.016882954640000E+00 -6.188295464000000E-01  6.188295464000000E-01
+ lp          =                      81                     81                     81                     81
+ lq          =                       0                      0                      0                      0
+ rp          =                       0                      0                      0                      0
+ rq          =                       0                      0                      0                      0
+ oita        =  -1.618000000000000E+00 -1.016882954640000E+00 -6.188295464000000E-01  6.188295464000000E-01
+ mupftol     =   1.000000000000000E-16
+ mupfits     =         8
+ Rac         =   0.000000000000000E+00
+ Zas         =   0.000000000000000E+00
+ Ras         =   0.000000000000000E+00
+ Zac         =   0.000000000000000E+00
+Rbc(0,0)    =  1.000000000000000E+01 Zbs(0,0)    =  0.000000000000000E+00 Rbs(0,0)    =  0.000000000000000E+00 Zbc(0,0)    =  0.000000000000000E+00
+Rbc(0,1)    =  0.000000000000000E+00 Zbs(0,1)    =  0.000000000000000E+00 Rbs(0,1)    =  0.000000000000000E+00 Zbc(0,1)    =  0.000000000000000E+00
+Rbc(0,2)    =  0.000000000000000E+00 Zbs(0,2)    =  0.000000000000000E+00 Rbs(0,2)    =  0.000000000000000E+00 Zbc(0,2)    =  0.000000000000000E+00
+Rwc(0,0)    =  0.000000000000000E+00 Zws(0,0)    =  0.000000000000000E+00 Rws(0,0)    =  0.000000000000000E+00 Zwc(0,0)    =  0.000000000000000E+00
+Rwc(0,1)    =  0.000000000000000E+00 Zws(0,1)    =  0.000000000000000E+00 Rws(0,1)    =  0.000000000000000E+00 Zwc(0,1)    =  0.000000000000000E+00
+Rwc(0,2)    =  0.000000000000000E+00 Zws(0,2)    =  0.000000000000000E+00 Rws(0,2)    =  0.000000000000000E+00 Zwc(0,2)    =  0.000000000000000E+00
+Vns(0,0)    =  0.000000000000000E+00 Bns(0,0)    =  0.000000000000000E+00 Vnc(0,0)    =  0.000000000000000E+00 Bnc(0,0)    =  0.000000000000000E+00
+Vns(0,1)    =  0.000000000000000E+00 Bns(0,1)    =  0.000000000000000E+00 Vnc(0,1)    =  0.000000000000000E+00 Bnc(0,1)    =  0.000000000000000E+00
+Vns(0,2)    =  0.000000000000000E+00 Bns(0,2)    =  0.000000000000000E+00 Vnc(0,2)    =  0.000000000000000E+00 Bnc(0,2)    =  0.000000000000000E+00
+/
+&numericlist
+ Linitialize =         1
+ Lzerovac    =         0
+ Ndiscrete   =         2
+ Nquad       =        -1
+ iMpol       =        -4
+ iNtor       =        -4
+ Lsparse     =         0
+ Lsvdiota    =         0
+ imethod     =         3
+ iorder      =         2
+ iprecon     =         1
+ iotatol     =  -1.000000000000000E+00
+ Lextrap     =         0
+ Mregular    =        -1
+/
+&locallist
+ LBeltrami   =         4
+ Linitgues   =         1
+/
+&globallist
+ Lfindzero   =         2
+ escale      =   0.000000000000000E+00
+ opsilon     =   1.000000000000000E+00
+ pcondense   =   2.000000000000000E+00
+ epsilon     =   0.000000000000000E+00
+ wpoloidal   =   1.000000000000000E+00
+ upsilon     =   1.000000000000000E+00
+ forcetol    =   1.000000000000000E-12
+ c05xmax     =   1.000000000000000E-06
+ c05xtol     =   1.000000000000000E-10
+ c05factor   =   1.000000000000000E-02
+ LreadGF     =         F
+ mfreeits    =         0
+ gBntol      =   1.000000000000000E-06
+ gBnbld      =   6.660000000000000E-01
+ vcasingeps  =   1.000000000000000E-12
+ vcasingtol  =   1.000000000000000E-08
+ vcasingits  =         8
+ vcasingper  =         1
+/
+&diagnosticslist
+ odetol      =   1.000000000000000E-07
+ nPpts       =       500
+ nPtrj       =      5    10     5
+ LHevalues   =         F
+ LHevectors  =         F
+ LHmatrix    =         T
+ Lperturbed  =         0
+ dpp         =        -1
+ dqq         =        -1
+ Lcheck      =         1
+ Ltiming     =         F
+/
+&screenlist
+ Wpp00aa = T
+ Whesian = T
+/
+     0     0  4.002650778313936E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  5.997349221686633E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  1.000000000000000E+01  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
+     1     0 -3.027260137998658E-24  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00 -1.425242065616607E-28  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
+     2     0 -2.762051886527595E-25  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  2.843088394703559E-29  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00

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;
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
@@ -902,6 +970,17 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
+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 !< total number of volumes (including the vacuum region in the case of free-boundary calculations)
 
+  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;
 
@@ -1678,7 +1681,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)