GPEC - NEW FEATURE - Emit optional perturbed Boozer and Hamada field strength

.github/workflows/make_test.yml

@@ -221,6 +221,14 @@ jobs:
           ulimit -s unlimited
           ../../../bin/dcon && ../../../bin/gpec && ../../../bin/stride && ../../../bin/pentrc
 
+      - name: Check perturbed-coordinate gauge operator
+        run: |
+          export LD_LIBRARY_PATH=$PREFIX/lib:$LD_LIBRARY_PATH
+          gfortran -O0 -fopenmp -Igpec -Iequil regression/test/pmodb_target_gauge.f90 gpec/pmodb_gauge.o gpec/pmodb_gauge_kernel.o gpec/ismath.o gpec/gpglobal.o -Ldeps/lib -lequil -llapack -lblas -o /tmp/pmodb_target_gauge
+          /tmp/pmodb_target_gauge
+          if /tmp/pmodb_target_gauge resonant > /tmp/resonant.log 2>&1; then exit 1; fi
+          python3 -c "assert 'Nonzero resonant pmodb gauge source' in open('/tmp/resonant.log').read()"
+
       - name: Run Soloviev kinetic example
         run: |
           export GPECHOME=$(pwd)

docs/namelists.rst

@@ -45,3 +45,5 @@ This namelist sets the inputs, controls the behavior, and determines the outputs
 .. include:: ../input/pentrc.in
    :literal:
 
+
+See :doc:`pmodb_target_gauge` for the optional ``pmodb_gauge_flag`` output.

docs/outputs.rst

@@ -290,3 +290,5 @@ solutions.bin, solutions<n>.bin
     **Flag** sol_num (match.in, MATCH)
 
     **Info** Ideal Euler-Lagrange solution built from DCON eigenfunction sol_num (1 is the most unstable). sol_num=1 writes solutions.bin; sol_num=n>1 writes solutions<n>.bin so several eigenfunctions can be kept side by side.
+
+See :doc:`pmodb_target_gauge` for ``b_lag_perp`` and optional perturbed Boozer/Hamada field strength.

docs/pmodb_target_gauge.rst

@@ -0,0 +1,80 @@
+Optional perturbed-coordinate field strength
+=============================================
+
+With ``pmodb_flag=t`` and ``reg_flag=f``, NetCDF ``b_lag_perp`` supplies
+``eulb+xi_perp.grad B``. It is the perpendicular-displacement value called
+``llagb`` in ``gpdiag_pmodb``. It is omitted for regularized solutions, whose reconstructed surface
+displacement need not equal xms/chi1. Existing ASCII columns keep their format.
+This scalar alone is not the perturbed-Boozer input to NEO-RT or NEO-2.
+
+Set ``pmodb_gauge_flag=t`` in ``GPEC_OUTPUT`` for an additional NetCDF
+field ``b_lag_boozer`` or ``b_lag_hamada``. This requires ``reg_flag=f``
+and a matching ``jac_type='boozer'`` or ``'hamada'``. The gauge is selected
+by the working chart; ``jac_out`` selects only the harmonic representation.
+The field has units Tesla and the same complex harmonic layout as ``b_lag``.
+The output's ``gauge`` attribute records the physical gauge separately from
+its ``jacobian`` representation. Use Boozer for NEO-type scalar-only
+consumers and Hamada for Shaing-type formulas with that gauge contract.
+
+Verified derivation
+-------------------
+
+GPEC angles are in turns and modes are exp(2 pi i(m theta-n zeta)). Write
+``D=partial_theta+q partial_zeta``, so ``B.grad=(chi1/J)D``. On the matching
+base ``J B^a=H(psi)``, with a=2 for Boozer or a=0 for Hamada. Under the
+parallel shift ``xi_target=xi_X+lambda B``, ideal induction is unchanged,
+``b_target=b_X+lambda B.grad B`` and ``d_target=d_X+B.grad lambda``.
+The target condition for a nonaxisymmetric mode is ``d_target+a b_target/B=0``.
+The product rule gives
+
+.. math::
+
+   T=d_X+a b_X/B,\qquad
+   T_{target}=T+\frac{\mathbf B\cdot\nabla(\lambda B^a)}{B^a},\qquad
+   w=\chi_1\lambda/J,\qquad
+   Dw=-T,\qquad w_m=\frac{iT_m}{2\pi(m-nq)},\qquad
+   b_{target}=b_X+w\partial_\theta B.
+
+The implementation forms T in real space, truncates to the existing GPEC
+poloidal spectrum, inverts each nonresonant mode, and forms the scalar in
+real space before the existing output-coordinate transform. Identities are
+exact; numerical output retains GPEC's finite spectral/spline resolution.
+A nonzero resonant source is rejected because no periodic solution exists.
+At a resonant zero source the free homogeneous shift is set to zero.
+Numerically, abs(m-nq)<1e-10 is classified as resonant; a source coefficient
+below 1e-12 times the largest source coefficient is treated as zero there.
+The near-resonant result can be large: this feature applies no smoothing.
+
+The input pair is ``b_X=b_lag_perp`` and the ordinary perpendicular
+coordinate divergence, using the derivative selected by gpeq_sol. Regularized fields/derivatives are deliberately
+excluded because they need not form one ideal-induction displacement.
+
+Verification and scope
+-----------------------
+
+``gen/pmodb_gauge/prove.py`` proves Q1-Q8 and K5/P1 using FortSym's native
+zero test, with no probe fallback (verified at ae635c80). The source, complex inverse
+coefficient and scalar correction are separately proved and emitted by
+``gen/pmodb_gauge/generate.py`` to ``gpec/pmodb_gauge_kernel.f90``.
+Regenerate and check the resulting diff; ordinary builds require no CAS.
+
+``regression/test/pmodb_target_gauge.f90`` exercises the actual operator:
+parallel displacement gives zero in both gauges, a zero target source
+leaves the scalar unchanged, and a nonzero resonant source fails. It also
+writes a Hamada scalar transformed on Boozer equilibrium angles.
+``regression/test/pmodb_target_compare.py`` compares it with rmp_torque's
+c3614334/#55 kernel (external source supplied through PYTHONPATH).
+
+For the Hamada scalar already relabeled onto Boozer equilibrium angles,
+T=2 b_H/B; both paths use the same inverse-field-line operator. Equilibrium
+relabeling must precede this conversion. PEST output is not supported as a
+target; its J/R^2 constraint also involves radial geometric displacement.
+This feature does not infer a new equilibrium chart from a scalar spectrum.
+
+Source contracts: ``gpec/gpeq.f:131`` selects the derivative;
+``gpec/gpdiag.f:1061`` defines the existing perpendicular scalar;
+``gpec/gpout.f`` computes it before output transformation. The scalar-only
+consumer formulas are those derived in rmp_torque c3614334,
+``docs/boozer_gauge_ntv_derivation.md`` Eqs. 4.3-4.8 and 4.13-4.15.
+Park 2011 PoP 18 110702 Eqs. 10-11 explain why a general kinetic action
+with an explicit Jacobian term has a different input contract.

docs/site.rst

@@ -16,3 +16,4 @@ GPEC documentation website site map
    source_documentation.rst
    installation.rst
 
+   pmodb_target_gauge.rst

gen/pmodb_gauge/generate.py

@@ -0,0 +1,21 @@
+"""Prove and emit the real leaves; Fourier indexing stays in pmodb_gauge_mod."""
+from pathlib import Path
+import fortsym
+import fortsym.sympy as fs
+s, Tr, Ti = fs.symbols('s Tr Ti')
+weight = 1/(2*fs.pi*s)
+wr, wi = -weight*Ti, weight*Tr
+assert (2*fs.pi*s*weight-1).is_zero is True
+assert (2*fs.pi*fs.I*s*(wr+fs.I*wi)+Tr+fs.I*Ti).is_zero is True
+a, b, bt, dr, di, vr, vi, ur, ui = fs.symbols('a b bt dr di vr vi ur ui')
+source = {'tr': vr+a*dr/b, 'ti': vi+a*di/b}
+target = {'out_r': dr+ur*bt, 'out_i': di+ui*bt}
+assert (source['tr']+fs.I*source['ti']-(vr+fs.I*vi)-a*(dr+fs.I*di)/b).is_zero is True
+assert (target['out_r']+fs.I*target['out_i']-(dr+fs.I*di)-(ur+fs.I*ui)*bt).is_zero is True
+kernels = [
+    fortsym.emit_kernel('pmodb_gauge_inverse', ['s','Tr','Ti'], {'wr':wr,'wi':wi}, generator='gen/pmodb_gauge/generate.py'),
+    fortsym.emit_kernel('pmodb_gauge_source', ['a','b','dr','di','vr','vi'], source, generator='gen/pmodb_gauge/generate.py'),
+    fortsym.emit_kernel('pmodb_gauge_apply', ['dr','di','ur','ui','bt'], target, generator='gen/pmodb_gauge/generate.py')]
+root = Path(__file__).resolve().parents[2]
+(root/'gpec/pmodb_gauge_kernel.f90').write_text('module pmodb_gauge_kernel\ncontains\n'+'\n'.join(kernels)+'\nend module pmodb_gauge_kernel\n')
+print('PROVED and emitted: complex source, inverse field-line coefficient, scalar correction')

gen/pmodb_gauge/prove.py

@@ -0,0 +1,65 @@
+import fortsym.sympy as sp
+
+def gauge_coordinates(sp):
+    """Local Jacobian and field-strength gauge equations in angles in turns.
+
+    B.grad = chi1/J (d_theta + q d_zeta). For a target chart J B**a=H(psi),
+    a=0 (Hamada) or 2 (Boozer), w=chi1 lambda B**a/H obeys D w=-T,
+    T=div xi + a delta B_L/B, for n!=0 (zero flux-function perturbation).
+    """
+    B, Bt, J, chi, q, lam, Lt, Lz, d, b, c = sp.symbols(
+        "B Bt J chi q lam Lt Lz d b c")
+    m,n=sp.symbols("m n")
+    out=[]
+    for a,name in ((0,"hamada"),(2,"boozer")):
+        # First jet of lambda B**a; axisymmetric equilibrium, div B=0.
+        Ftheta = Lt*B**a + a*lam*B**(a-1)*Bt
+        Fzeta = Lz*B**a
+        direct = d + chi/J*(Lt+q*Lz) + a*(b+lam*chi/J*Bt)/B
+        product = d+a*b/B+chi/(J*B**a)*(Ftheta+q*Fzeta)
+        out.append(("Q1."+name,"T_a(xi+lambda B) = T_a(xi)+B.grad(lambda B**a)/B**a",direct-product))
+        # Solving the Fourier equation, with modes exp(2 pi i(m theta-n zeta)).
+        T,H=sp.symbols("T H")
+        w=-T/(2*sp.pi*sp.I*(m-n*q))
+        out.append(("Q2."+name,"D w = -T; w_m = -T_m/[2 pi i(m-nq)]",2*sp.pi*sp.I*(m-n*q)*w+T))
+        # On a target base J B**a=H, correction to field strength is w d_theta B.
+        lambda_w = w*J/chi
+        out.append(("Q3."+name,"delta B_target = delta B_X + w d_theta B (J B**a flux function)",
+                    lambda_w*chi/J*Bt-w*Bt))
+    # Transform the base angles: theta_B=theta_H+iota h, phi_B=phi_H+h.
+    h,Dh,iota,Hh,Hb=sp.symbols("h Dh iota Hh Hb")
+    ht,hp=sp.symbols("ht hp")
+    determinant=(1+iota*ht)*(1+hp)-ht*iota*hp
+    out.append(("Q4","det d(theta_H+iota h,phi_H+h)/d(theta_H,phi_H)=1+D h",
+                determinant-1-iota*ht-hp))
+    # Same physical scalar and relabeling, in the convention of rmp_torque E4.14.
+    bX,jX,u,v=sp.symbols("bX jX u v")
+    dbB=bX+u*Bt; dbH=bX+v*Bt
+    out.append(("Q5","Hamada->Boozer equals direct X->Boozer when w_HB=w_XB-w_XH",
+                dbH+(u-v)*Bt-dbB))
+    # Perpendicular reconstruction is unchanged by adding a field-aligned vector.
+    out.append(("Q6","field-aligned input: delta B_E=0, xi_perp=0 => target delta B=0",
+                lam*chi/J*Bt + (-lam*chi/J)*Bt))
+    out.append(("Q7","xi_perp limit: zero target source implies no gauge correction", (bX+u*Bt-bX).subs(u,0)))
+    # PEST J/R**2 is a flux function. Its linearized constraint has radial motion.
+    R,xiR,divxi,db,dr=sp.symbols("R xiR divxi db dr")
+    eps=sp.Symbol("eps")
+    ratio=J*(1+eps*divxi)/(R+eps*xiR)**2
+    out.append(("Q8","delta(J/R**2)/(J/R**2)=div xi-2 xi_R/R",
+                ratio.diff(eps).subs(eps,0)*R**2/J-divxi+2*xiR/R))
+    # Complex CGL quadratic form: real and imaginary components independently.
+    Dr,Di,Kr,Ki=sp.symbols("Dr Di Kr Ki")
+    pressure=0; code=0
+    for D,K in ((Dr,Kr),(Di,Ki)):
+        b=-(D+K); pp=-D+b; pa=-3*D-2*b
+        pressure-=((pa-pp)*b+pa*D)/2
+        code+=(sp.Rational(5,3)*D**2+sp.Rational(1,3)*(D+3*K)**2)/2
+    out.append(("K5","Complex CGL pressure work equals [5/3 |D|**2+1/3 |D+3K|**2]/2",pressure-code))
+    # Pressure term in curvature is mu0 p' xi^psi/B, not mu0 p' xi^psi.
+    P,X=sp.symbols("P X")
+    out.append(("P1","(Bk - (llagb-eulb))/B = mu0 p' xi^psi/B**2", (P*X/B)/B-P*X/B**2))
+    return out
+
+for label,equation,residual in gauge_coordinates(sp):
+    assert sp.sympify(residual).is_zero is True, label
+    print(label,"PROVED")

gpec/gpec.f

@@ -33,6 +33,7 @@ PROGRAM gpec_main
       CHARACTER(8) :: filter_types
       CHARACTER(128) :: infile
       LOGICAL :: singcoup_flag,singfld_flag,vsingfld_flag,pmodb_flag,
+     $     pmodb_gauge_flag,
      $     xbcontra_flag,xbnormal_flag,vbnormal_flag,xbnobo_flag,
      $     d3_flag,xbst_flag,rzphibx_flag,jacfac_flag,dw_flag,
      $     radvar_flag,eigen_flag,magpot_flag,xbtangent_flag,
@@ -64,7 +65,8 @@ PROGRAM gpec_main
      $     reg_spot,chebyshev_flag,nche,nchr,nchz,use_classic_splines,
      $     use_res_spot
       NAMELIST/gpec_output/resp_flag,singcoup_flag,nrzeq_flag,nr,nz,
-     $     singfld_flag,pmodb_flag,xbnormal_flag,rstep,jsurf_out,
+     $     singfld_flag,pmodb_flag,pmodb_gauge_flag,xbnormal_flag,
+     $     rstep,jsurf_out,
      $     jac_out,power_bout,power_rout,power_bpout,power_rcout,
      $     tmag_out,mlim_out,eqbrzphi_flag,brzphi_flag,xrzphi_flag,
      $     vbrzphi_flag,vvbrzphi_flag,divzero_flag,dw_flag,opsi1,opsi2,
@@ -144,6 +146,7 @@ PROGRAM gpec_main
       singfld_flag=.TRUE.
       vsingfld_flag=.FALSE.
       pmodb_flag=.FALSE.
+      pmodb_gauge_flag=.FALSE.
       xbnormal_flag=.TRUE.
       vbnormal_flag=.FALSE.
       xbtangent_flag=.FALSE.
@@ -622,7 +625,7 @@ PROGRAM gpec_main
          CALL gpout_dw_matrix(coil_flag)
       ENDIF
       IF (pmodb_flag) THEN
-         CALL gpout_pmodb(mode,xspmn)
+         CALL gpout_pmodb(mode,xspmn,pmodb_gauge_flag)
       ENDIF
       IF (xbnormal_flag) THEN
          CALL gpout_xbnormal(mode,xspmn,sing_spots,sing_interpspot,