Implement Bernstein's breakdown for the perturbed potential energy

.github/workflows/make_test.yml

@@ -231,6 +231,14 @@ jobs:
           ulimit -s unlimited
           ../../../bin/rdcon && ../../../bin/rmatch && ../../../bin/gpec && ../../../bin/stride
 
+      - name: Run Soloviev energy decomposition example
+        run: |
+          export GPECHOME=$(pwd)
+          export LD_LIBRARY_PATH=$PREFIX/lib:$LD_LIBRARY_PATH
+          cd docs/examples/regression_solovev_recon_example
+          ulimit -s unlimited
+          ../../../bin/dcon && ../../../bin/gpec && python3 check_recon.py
+
   make_gpec_ifort:
     runs-on: ubuntu-latest
     needs: build-netcdf-ifort
@@ -310,3 +318,13 @@ jobs:
           cd docs/examples/regression_solovev_resistive_example
           ulimit -s unlimited
           ../../../bin/rdcon && ../../../bin/rmatch && ../../../bin/gpec && ../../../bin/stride
+
+      - name: Run Soloviev energy decomposition example
+        run: |
+          source /opt/intel/oneapi/setvars.sh
+          export GPECHOME=$(pwd)
+          export LD_LIBRARY_PATH=$PREFIX/lib:$LD_LIBRARY_PATH
+          export FOR_DISABLE_STACK_TRACE=TRUE
+          cd docs/examples/regression_solovev_recon_example
+          ulimit -s unlimited
+          ../../../bin/dcon && ../../../bin/gpec && python3 check_recon.py

docs/examples/regression_solovev_recon_example/README

@@ -0,0 +1,27 @@
+Soloviev Energy Decomposition Regression Example
+================================================
+
+This example checks the GPEC ideal energy decomposition (effective_field_flag) against the DCON eigenvalue it reconstructs. It reuses the Soloviev equilibrium of regression_solovev_ideal_example (sol.in, equil.in, vac.in are copies). dcon.in differs only in delta_mlow=delta_mhigh=8: the effective-field form relies on identities that hold only for untruncated fields, so the default band leaves a ~21% error, while 8 extra harmonics on each side bring it to 0.4%. (standard_form_flag matches ep(1) to 0.02% even with the default band.)
+
+gpec.in selects the least stable eigenmode (mode_flag=t, mode=1) with reg_flag=f and no external field, turns off the usual outputs, and sets effective_field_flag=t and curl_check_flag=t.
+
+Run from this directory::
+
+  ../../../bin/dcon
+  ../../../bin/gpec
+  python3 check_recon.py
+
+check_recon.py reads the last GPEC_RECON block of gpec.log (gpec.log is appended to, so delete it before re-running by hand) and passes when
+
+- the recon psi range equals the DCON psi range (psilow to psilim),
+- dW_p(recon-normalize) matches DCON ep(1) within 2% relative, and
+- the curl check reports PASS (a consistency check on C; it passes even when the truncated run misses the energy).
+
+Reference values (gfortran 15, macOS arm64)::
+
+  recon psi range       = 1.0E-04  9.9999E-01   (same as DCON psi range)
+  dW_p(recon-normalize) = -1.04353063E+001
+  DCON ep(1)            = -1.03926066E+001   (relative difference 4.1E-03)
+  C verify: PASS        (rel 2.8E-13 against tol 1E-08)
+
+The 2% tolerance leaves room for compiler and optimization differences while still catching a broken decomposition: the truncated 8-harmonic run above misses by 21%.

docs/examples/regression_solovev_recon_example/check_recon.py

@@ -0,0 +1,66 @@
+#!/usr/bin/env python3
+"""
+Regression check for the GPEC energy decomposition (effective_field_flag).
+
+Run after ``dcon`` and ``gpec`` in this directory. Reads the last
+GPEC_RECON block of gpec.log and checks that
+
+  * the recon psi range covers the DCON psi range (psilow to psilim),
+  * dW_p(recon-normalize) matches DCON ep(mode) within REL_TOL, and
+  * the curl check reports PASS.
+
+Exits 0 on success and 1 on failure, so it can be chained in CI.
+"""
+
+import sys
+
+REL_TOL = 0.02  # reference run agrees to 0.4%; see README
+RANGE_TOL = 1e-6
+
+
+def last_recon_block(path):
+    with open(path) as f:
+        lines = f.read().splitlines()
+    starts = [i for i, l in enumerate(lines)
+              if l.startswith("GPEC_RECON final results:")]
+    if not starts:
+        sys.exit("FAIL: no GPEC_RECON block in " + path)
+    block = []
+    for l in lines[starts[-1] + 1:]:
+        if l.startswith("GPEC_RECON"):
+            break
+        block.append(l)
+    return block
+
+
+def values(block, key):
+    for l in block:
+        if key in l:
+            return [float(v) for v in l.split("=")[-1].split()]
+    sys.exit("FAIL: '{}' not found in GPEC_RECON block".format(key))
+
+
+def main():
+    block = last_recon_block("gpec.log")
+    dw = values(block, "dW_p(recon-normalize)")[0]
+    ep = values(block, "DCON ep(")[0]
+    rlo, rhi = values(block, "recon psi range")
+    dlo, dhi = values(block, "DCON psi range")
+    rel = abs(dw - ep) / abs(ep)
+    verdict = [l.split(":")[-1].strip() for l in block if "C verify" in l]
+    range_ok = abs(rlo - dlo) <= RANGE_TOL and abs(rhi - dhi) <= RANGE_TOL
+    print("recon psi range       = {:.6e} {:.6e}".format(rlo, rhi))
+    print("DCON psi range        = {:.6e} {:.6e}".format(dlo, dhi))
+    print("dW_p(recon-normalize) = {: .8e}".format(dw))
+    print("DCON ep               = {: .8e}".format(ep))
+    print("relative difference   = {:.3e} (tol {:.1e})".format(rel, REL_TOL))
+    print("curl check            = {}".format(verdict[-1] if verdict
+                                              else "not run"))
+    ok = (range_ok and rel <= REL_TOL and bool(verdict)
+          and verdict[-1] == "PASS")
+    print("PASS" if ok else "FAIL")
+    return 0 if ok else 1
+
+
+if __name__ == "__main__":
+    sys.exit(main())

docs/examples/regression_solovev_recon_example/dcon.in

@@ -0,0 +1,61 @@
+&DCON_CONTROL
+    bal_flag=f           ! Ideal MHD ballooning criterion for short wavelengths
+    mat_flag=t           ! Construct coefficient matrices for diagnostic purposes
+    ode_flag=t           ! Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC)
+    vac_flag=t           ! Compute plasma, vacuum, and total energies for free-boundary modes
+    mer_flag=t           ! Evaluate the Mercier criterian
+
+    sas_flag=f           ! Safety factor (q) limit determined as q_ir+dmlim where q_ir is the equil outermost rational
+    dmlim=0.2            ! See sas_flag
+    qlow=1.02            ! Integration initiated at q determined by minimum of qlow and q0 from equil
+    qhigh=1e3            ! Integration terminated at q limit determined by minimum of qhigh and qa from equil
+    sing_start=0         ! Start integration at the sing_start'th rational from the axis (psilow)
+
+    nn=1                 ! Toroidal mode number
+    delta_mlow=8         ! Expands lower bound of Fourier harmonics
+    delta_mhigh=8        ! Expands upper bound of  Fourier harmonics
+    delta_mband=0        ! Integration keeps only this wide a band of solutions along the diagonal in m,m'
+    mthvac=64            ! Number of points used in splines over poloidal angle at plasma-vacuum interface. Overrides vac.in mth.
+    thmax0=1             ! Linear multiplier on the automatic choice of theta integration bounds for high-n ideal ballooning stability computation (strictly, -inf to -inf)
+
+    kin_flag = f         ! Kinetic EL equation (default: false)
+    con_flag = f         ! Continue integration through layers (default: false)
+    kinfac1 = 1.0        ! Scale factor for energy contribution (default : 1.0)
+    kinfac2 = 1.0        ! Scale factor for torque contribution (default : 1.0)
+    kingridtype = 0      ! Regular grid method (default : 0)
+    passing_flag = f     ! Includes passing particle effects (default: false)
+    ktanh_flag = t       ! Ignore kinetic effects in the core smoothly (default: false)
+    ktc = 0.1            ! Parameter activated by ktanh_flag: roughly corresponds to core width ignored (default: 0.1)
+    ktw = 50.0           ! Parameter activated by ktanh_flag: width of hyper-tangential functions (default: 50.0)
+    ion_flag = t         ! Include ion dW_k when kin_flag is true (summed with electron contribution if electron_flag true)
+    electron_flag = f    ! Include electron dW_k when kin_flag is true (summed with ion contribution if ion_flag true)
+
+    tol_nr=1e-6          ! Relative tolerance of dynamic integration steps away from rationals
+    tol_r=1e-7           ! Relative tolerance of dynamic integration steps near rationals
+    crossover=1e-2       ! Fractional distance from rational q at which tolerance is switched to tol_r
+    singfac_min=1e-4     ! Fractional distance from rational q at which ideal jump condition is enforced
+    ucrit=1e3            ! Maximum fraction of solutions allowed before re-normalized
+
+    wv_farwall_flag = f  ! Force nowall gpec calculations while calculating mutual inductance with the wall, when set true.
+    termbycross_flag = f ! Terminate ODE solver in the event of a zero crossing
+
+    use_classic_splines = f ! Use a classical cubic spline instead of tri-diagonal solution for splines with extrapolation boundary conditions
+/
+&DCON_OUTPUT
+    crit_break=t         ! Color of the crit curve changes when crossing a singular surface
+
+    ahb_flag=f           ! Output normal magnetic field eigenvalues and eigenfunctons at plasma-vacuum interface (must be false for GPEC)
+    msol_ahb=1           ! Number of eigenfunctions output by ahb_flag=t ?
+    mthsurf0=1           ! Linear multiplier on number of boundary points used to display surface eigenfunctions for ahgb_flag=t
+
+    bin_euler=t          ! Output M psi-by-M euler-lagrange solutions to binary file euler.bin
+    bin_vac=f            ! Output deprecated binary vacuum.bin containing vacuum energy matrices
+    euler_stride=1       ! Output only every euler_stride'th psi step to binary file
+
+    out_bal1=f           ! Ascii output for bal_flag poloidal functions
+    bin_bal1=f           ! Binary output for bal_flag poloidal functions
+    out_bal2=f           ! Ascii output for bal_flag functions
+    bin_bal2=f           ! Binary output for bal_flag functions
+
+    netcdf_out=t         ! Replicate ascii dcon.out information in a netcdf file
+/

docs/examples/regression_solovev_recon_example/equil.in

@@ -0,0 +1,30 @@
+&EQUIL_CONTROL
+    eq_type="soloviev"       ! Type of the input 2D equilibrium file. Accepts efit, chease, fluxgrid, transp, jsolver, lar, sol, etc.
+    eq_filename = "sol.in"   ! Path to input file
+
+    jac_type="pest"      ! Working coordinate system for all DCON and GPEC calculations. Overrides individual powers. Accepts hamada, pest, boozer, equal_arc
+    power_bp=0           ! del.B ~ B_p**power_bp * B**power_b / R**power_r
+    power_b=0            ! del.B ~ B_p**power_bp * B**power_b / R**power_r
+    power_r=0            ! del.B ~ B_p**power_bp * B**power_b / R**power_r
+
+    grid_type="ldp"      ! Radial grid packing of equilibrium quantities. Accepts rho, ldp, pow1, pow2, or original. ldp packs points near the core and edge. pow* packs near the edge.
+    psilow=1e-4          ! Minimum value of psi, normalized from 0 to 1
+    psihigh=0.99999      ! Maximum value of psi, normalized from 0 to 1
+    mpsi=16             ! Number of radial grid intervals for equilibrium quantities
+    mtheta=256           ! Number of equally spaced poloidal grid intervals for all splines
+    newq0=0              ! Grad-Shafranov solution invariant adjustment of the q profile to give the specified value of q at the axis. Default 0 uses input file value.
+    use_classic_splines = f ! Use a classical cubic spline instead of tri-diagonal solution for splines with extrapolation boundary conditions
+
+    input_only=f         ! Generate information about the input and then quit with no further calculation
+/
+&EQUIL_OUTPUT
+    gse_flag=f           ! Produces diagnostic output for accuracy of solution to Grad-Shafranov equation
+    out_eq_1d=f          ! Ascii output of 1D equilibrium file data
+    bin_eq_1d=f          ! Binary output of 1D equilibrium file data
+    out_eq_2d=f          ! Ascii output of 2D equilibrium file data
+    bin_eq_2d=t          ! Binary output of 2D equilibrium file data (set true for GPEC)
+    out_2d=f             ! Ascii output of processed 2D data
+    bin_2d=f             ! Binary output of processed 2D data
+    dump_flag=f          ! Binary dump of basic equilibrium data and 2D rzphi spline
+    out_ahg2msc=f        ! Output deprecated ahg2msc.out files. This is used to communicate with vacuum, but is now done through memory.
+/

docs/examples/regression_solovev_recon_example/gpec.in

@@ -0,0 +1,92 @@
+&GPEC_INPUT
+   dcon_dir   = './'                 ! Path to DCON outputs. Individual file paths are appended to this
+   idconfile  = 'euler.bin'          ! Eigenfunctions output by DCON run
+   ieqfile    = 'psi_in.bin'         ! Flux surfaces output by DCON run
+   rdconfile  = 'globalsol.bin'      ! Galerkin eigenfunctions output by resistive DCON run
+
+   gal_flag = f             ! Use the Galerkin eigenfunctions from resistive DCON instead of the ideal DCON shooting method solutions
+
+   jac_in="pest"            ! Coordinate system of external perturbation (affects data_flag and harmonic_flag). Options include: pest,boozer,hamada
+   jsurf_in=0               ! True(1) if external perturbation is area weighted (i.e. flux)
+   tmag_in =0               ! True(1) if external perturbation has toroidal angle defined by jac_in. False(0) if machine angle.
+   mthsurf =0               ! Number of poloidal angle gridpoints. Minimum of twice the Nyquist condition from DCON m limits is enforced.
+
+   coil_flag=f              ! Calculate external perturbation from coils in coil.in
+
+   data_flag=f              ! Apply perturbation from thrid party code output
+   data_type="surfmn"       ! Name of said code. Choose: surfmn,??
+   infile="path/to/file"    ! Path to said output
+   nmin=0                   ! Extra parameters needed to read surfmn files.
+   nmax=64                  ! Extra parameters needed to read surfmn files.
+   mmin=-64                 ! Extra parameters needed to read surfmn files.
+   mmax=64                  ! Extra parameters needed to read surfmn files.
+
+   harmonic_flag=f          ! Apply self-defined perturbations
+   cosmn(3)=1e-4            ! Real part of m=3 component (n defined by DCON run).
+   sinmn(6)=1e-4            ! Imaginary part of m=6 component. etc.
+
+   displacement_flag=f      ! Perturbations are displacement (m). Default is field (T).
+   fixed_boundary_flag=f    ! Total perturbation = external perturbation
+
+   mode_flag = t            ! Return only a single DCON eigenfunction defined by mode (ignores external pert.)
+   mode  = 1                ! DCON eigen function index (1 is least stable)
+   filter_types = 'w'       ! Isolate the filter_modes of energy (w), reluctance (r), permiability (p) or singular coupling (s) eigenvectors/singular-vectors
+   filter_modes = 0         ! Number of modes isolated, negative values isolate backwards from the last mode
+/
+&GPEC_CONTROL
+   resp_index=0             ! Integer 0-4, determining method of inductance matrix formation. 0 uses energy identity, and is recommended.
+   sing_spot=5e-4           ! Minimum distance to singular surface in (m-nq) space
+   sing_resspot=1.0         ! If use_res_spot=t, use inverse cube root of lundquist number scaled by sing_resspot in psi space as spot size
+   use_res_spot=f           ! Use resistive singular layer width instead of fixed sing_spot
+   sing_interpspot=1.0      ! Scaling of the interpolation width, in units of spot size, for smoothed outputs of singular surface quantities
+   sing_npsi=100            ! Number of points between rational surfaces used for b-normal splines across sin_spot
+   reg_flag=f               ! Regularize the solution near resonant surfaces
+   reg_spot=5e-2            ! Regularize by (m-nq)^2/((m-nq)^2+reg_spot^2)
+   chebyshev_flag=f         ! Fit Chebyshev polynomials to the displacement.
+   nche=20                  ! Number of Chebyshev polynomials used
+   use_classic_splines = f  ! Use a classical cubic spline instead of tri-diagonal solution for splines with extrapolation boundary conditions
+/
+&GPEC_OUTPUT
+   effective_field_flag=t   ! Bernstein effective-field reconstruction of dW_p (needs mode_flag=t, reg_flag=f, coil_flag=f)
+   curl_check_flag=t        ! Check (curl C).grad(psi)=0 and report PASS/WARN
+   recon_out=t              ! Write recon files and the gpec.log summary read by check_recon.py
+   jac_out="pest"           ! Coordinate system of outputs
+   jsurf_out=0              ! True(1) for area weighted outputs (only appends standard gpec_control output)
+   tmag_out=1               ! True(1) for magnetic toroidal angle outputs (affects singcoup and singfld, appends gpec_control)
+   mlim_out=64              ! If positive, Fourier decomposed outputs will span -mlim_out<=m<=mlim_out. Default (<=0) outputs on internal m range.
+
+   resp_flag=f              ! Output energy, reluctance, inductance, and permeability eigenvalues, eigenvectors, and matrices
+   filter_flag=f            ! Outputs energy, reluctance, permeability, and (optionally) the singular-coupling eigenmodes on the control surface
+   singcoup_flag=f          ! Calculate coupling of each m to resonant surfaces
+   singfld_flag=f           ! Output of resonant surface quantities (flux,current,etc.)
+   vsingfld_flag=f          ! Output of vacuum resonant surface quantities
+
+   pmodb_flag=f             ! Outputs delta-B_Lagrangian
+   xbnormal_flag=f          ! Outputs normal displacement and field in plasma
+   vbnormal_flag=f          ! Outputs vacuum normal displacement and field in plasma
+   xclebsch_flag=f          ! Outputs clebsch coordinate displacement for PENT
+   dw_flag=f                ! Outputs self-consistent energy and torque profiles (when kin_flag=t from DCON)
+
+   eqbrzphi_flag=f          ! Outputs equilibrium field on regular (r,z) grid
+   brzphi_flag=f            ! Outputs perturbed field on regular (r,z) grid
+   xrzphi_flag=f            ! Outputs displacement on regular (r,z) grid
+   vbrzphi_flag=f           ! Outputs field on regular (r,z) grid due to total boundary surface current (not a real field)
+   vvbrzphi_flag=f          ! Outputs field on regular (r,z) grid due to external perturbation surface current (not a real field)
+   vsbrzphi_flag=f          ! Outputs the field on regular (r,z) grid due to singular surface currents
+   ss_flag(3)=f             ! Outputs the field on regular (r,z) grid due to the 3rd singular surface current, etc.
+   xbrzphifun_flag=f        ! Outputs r,z, and perturbation vectors on psi,theta grid
+   arzphifun_flag=f         ! Outputs r,z, and vector potential on psi,theta grid
+
+   nr = 3                  ! Number of radial grid points in rzphi outputs
+   nz = 3                 ! Number of vertical grid points in rzphi outputs
+   fun_flag=f               ! Write outputs to (r,z) in addition to default (psi,m)
+   flux_flag=f              ! Write outputs to (psi,theta) in addition to default (psi,m)
+   bin_flag=f               ! Write binary outputs for use with xdraw
+   bin_2d_flag=f            ! Write 2D binary outputs for use with xdraw
+   max_linesout=0           ! Maximum length of ascii data tables enforced by increasing radial step size (iff .gt. 0).
+   verbose=t                ! Print run log to terminal
+/
+&GPEC_DIAGNOSE
+   timeit=t                 ! Print timer splits for major subroutines
+   radvar_flag=f            ! Map various radial variables (rho,psi_tor) on psi_n grid
+/

docs/examples/regression_solovev_recon_example/sol.in

@@ -0,0 +1,14 @@
+! Construct Soloviev equilibrium
+&SOL_INPUT
+   mr = 128      ! number of radial grid zones
+   mz = 128      ! number of axial grid zones
+   ma = 128      ! number of flux grid zones
+   e = 1.6       ! elongation
+   a = 0.33      ! minor radius
+   r0 = 1.0      ! major radius
+   q0 = 1.9      ! safety factor at the o-point
+   p0fac=1       ! scale on-axis pressure (P-> P+P0*p0fac. beta changes. Phi,q constant)
+   b0fac=1       ! scale toroidal field at constant beta (s*Phi,s*f,s^2*P. bt changes. Shape,beta constant)
+   f0fac=1       ! scale toroidal field at constant pressure (s*f. beta,q changes. Phi,p,bp constant)
+/
+