Equilibrium - FEATURE - 🌱 Add the Cerfon-Freidberg diverted analytic equilibrium

docs/development/references.md

@@ -101,6 +101,19 @@ The KineticForces module (formerly PENTRC) implements neoclassical toroidal visc
   - Published: PhD Thesis, Princeton University (2015)
   - Describes: Complete NTV theory and implementation. **Chapter 7** details the hybrid drift-kinetic MHD eigenfunction calculation: 6 kinetic matrices Ak,Bk,Ck,Dk,Ek,Hk (Eqs 7.30-7.35) as energy-space integrals of perturbed action operators WX,WY,WZ; hybrid Euler-Lagrange equations; resonance splitting/suppression where Fh=(Q-P†)F̄(Q-P)+... shifts singularities away from rational surfaces (Eq 7.46); convergence to ideal limit. **Appendix C** derives the DCON matrix form of the perturbed action (Eqs C.1-C.11) used to compute the kinetic coefficient matrices. **Appendix D** details numerical treatment of integrable singularities in bounce averages.
 
+## Equilibrium Module
+
+- **Cerfon and Freidberg (2010)**: "One size fits all analytic solutions to the Grad-Shafranov equation"
+  - Published: Physics of Plasmas **17**, 032502 (2010)
+  - Link: https://doi.org/10.1063/1.3328818
+  - Describes: The Solov'ev choice of constant `p'` and `FF'` linearizes Grad-Shafranov to
+    `Δ*ψ̂ = (1−A)x² + A` in normalized coordinates, solved by a particular solution plus a
+    twelve-term homogeneous basis (seven up-down symmetric, five antisymmetric). Point-wise
+    conditions on the `ψ̂ = 0` surface set the plasma shape and place magnetic nulls on the
+    boundary. Implemented in `src/Equilibrium/CerfonEquilibrium.jl` as `eq_type = "cerfon"`,
+    the diverted generalization of the `"sol"` model (which is the `A = 0`, no-null member of
+    the same family).
+
 ## Additional References
 
 - **Park et al. (2009)**: "Nonambipolar Transport by Trapped Particles in Tokamaks"

docs/src/equilibrium.md

@@ -134,9 +134,13 @@ Modules = [GeneralizedPerturbedEquilibrium.Equilibrium]
 - `PlasmaEquilibrium` — the runtime structure containing spline fields,
 	geometry, profiles, and computed diagnostics (q-profile, separatrix,
 	etc.).
-- `LargeAspectRatioConfig`, `SolovevConfig` — convenience structs to
-	construct analytic/model equilibria when using `eq_type = "lar"` or
-	`eq_type = "sol"`.
+- `LargeAspectRatioConfig`, `SolovevConfig`, `CerfonConfig` — convenience
+	structs to construct analytic/model equilibria when using
+	`eq_type = "lar"`, `"sol"` or `"cerfon"`. `CerfonConfig` is the diverted
+	Cerfon-Freidberg generalization of `SolovevConfig`: a magnetic null sits on
+	the plasma boundary, so `q` diverges logarithmically as ψ_N → 1 and no
+	choice of `psihigh` converges — which makes it the fixture for truncation
+	sensitivity studies.
 
 ## Key functions
 

examples/Cerfon_diverted_example/gpec.toml

@@ -0,0 +1,75 @@
+# Cerfon-Freidberg diverted analytic equilibrium — n=1 ideal stability, lower single null.
+# The plasma boundary carries a magnetic null, so q diverges logarithmically as psi_N -> 1
+# and no choice of psihigh converges. This is the truncation-sensitivity fixture: raise
+# psihigh by a decade and q at the TRUNCATION SURFACE climbs by a fixed increment forever
+# (about 0.31 per decade here, and twice that with null = "dn"). q at psi_N = 1 is infinite
+# in every configuration — it is never reached, psihigh always cuts the profile short.
+# q0 = 1.5 is a wall-stabilized operating point: n=1 is unstable with no wall and stable
+# with the conformal wall out to a gap of about 0.46, so a = 0.25 sits at roughly half the
+# critical gap. The equilibrium is generated analytically from the [CERFON_INPUT] section.
+
+[Equilibrium]
+eq_type = "cerfon"                      # Type of the input 2D equilibrium file
+jac_type = "pest"                       # Coordinate system (hamada, pest, boozer, equal_arc, park, custom)
+grid_type = "ldp"                       # Radial grid packing type
+psilow = 1e-4                           # Lower limit of normalized poloidal flux
+psihigh = 0.999                         # Upper limit of normalized poloidal flux
+mpsi = 256                              # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
+mtheta = 256                            # Number of poloidal grid points
+newq0 = 0                               # Override for on-axis safety factor (0 = use input value)
+etol = 1e-8                             # Error tolerance for equilibrium solver
+force_termination = false               # Terminate after equilibrium setup (skip stability calculations)
+
+[Wall]
+shape = "conformal"                     # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath)
+a = 0.25                                # Distance from plasma (conformal) or shape parameter
+aw = 0.05                               # Half-thickness parameter for Dee-shaped walls
+bw = 1.5                                # Elongation parameter for wall shapes
+cw = 0                                  # Offset of wall center from major radius
+dw = 0.5                                # Triangularity parameter for wall shapes
+tw = 0.05                               # Sharpness of wall corners (try 0.05 as initial value)
+equal_arc_wall = true                   # Equal arc length distribution of nodes on wall
+
+[ForceFreeStates]
+local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
+vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
+
+psiedge = 0.99                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable
+qlow = 1.02                   # Integration initiated at q determined by min(q0, qlow)
+qhigh = 1e3                   # Integration terminated at q limit determined by min(qa, qhigh)
+sing_start = 0                # Start integration at the sing_start'th rational from the axis (psilow)
+
+nn_low = 1                    # Smallest toroidal mode number to include
+nn_high = 1                   # Largest toroidal mode number to include
+delta_mlow = 8                # Expands lower bound of Fourier harmonics
+delta_mhigh = 8               # Expands upper bound of Fourier harmonics
+mthvac = 960                  # Number of points used in splines over poloidal angle at the plasma-vacuum interface
+
+kinetic_source = "fixed"      # Kinetic matrix source: "fixed" test matrices, or "calculated" from the kinetic NTV model
+kinetic_factor = 0.0          # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode)
+eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations
+singfac_min = 1e-4            # Fractional distance from rational q at which ideal jump enforced
+ucrit = 1e3                   # Column-norm threshold that triggers solution renormalization
+save_interval = 3             # Save every Nth ODE step (1=all). Always saves near rational surfaces.
+
+use_parallel = true           # Run parallel FM-propagator BVP path (unlocks singular/delta_prime_matrix)
+parallel_threads = 2          # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates)
+populate_dense_xi = false     # Append serial-EL pass so dense ξ is stored (required with a [PerturbedEquilibrium] section)
+set_psilim_via_dmlim = false  # FALSE here: psihigh is the truncation knob under study, so nothing else may move it
+dmlim = 0.2                   # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true)
+
+# Cerfon-Freidberg analytic equilibrium parameters (eq_type = "cerfon"); see CerfonConfig in src/Equilibrium.
+[CERFON_INPUT]
+mr = 384             # Number of radial grid zones in the tabulated psi(R,Z) map
+mz = 512             # Number of axial grid zones in the tabulated psi(R,Z) map
+ma = 256             # Number of flux grid zones for the F and mu0*p profile spline
+r0 = 1.0             # Major radius R0 in m
+b0 = 2.0             # Vacuum toroidal field at r0 in T
+epsilon = 0.32       # Inverse aspect ratio eps = a/R0
+kappa = 1.7          # Elongation kappa
+delta = 0.33         # Triangularity delta
+A = -0.155           # Solov'ev parameter A setting the p'/FF' split (0 = pressure-driven only)
+q0 = 1.5             # Target safety factor on axis; fixes the flux scale
+null = "lsn"         # Null topology: "lsn" (lower single null) or "dn" (double null)
+xsep = 1.1           # Places the null at (1 - xsep*delta*eps, -xsep*kappa*eps)
+box_margin = 0.12    # psi(R,Z) box padding beyond the boundary, in units of eps

regression-harness/cases/cerfon_lsn_n1.toml

@@ -0,0 +1,201 @@
+# Regression case: Cerfon-Freidberg diverted analytic equilibrium, n=1 ideal stability.
+# Lower single null with the magnetic null on the plasma boundary, so q is unbounded as
+# psi_N -> 1 and the run is truncated by psihigh rather than by a finite edge q. That makes
+# this the fixture that notices changes in edge/truncation handling: anything touching the
+# radial grid near psihigh, the flux-surface tracer, or the singular-surface search moves
+# q at the truncation surface and the outermost rational here long before it shows up in the
+# limited fixtures. The deck sits at a wall-stabilized operating point (n=1 stable with the
+# conformal wall, unstable without), so a sign change in the energies is a real signal, not
+# a marginal case wobbling about zero.
+[case]
+name = "cerfon_lsn_n1"
+description = "Cerfon-Freidberg diverted (LSN) equilibrium, n=1, ideal stability"
+example_dir = "examples/Cerfon_diverted_example"
+
+# Energies — leading generalized (W,N) pencil eigenvalues at the final truncation (psilim): power-normalized (⟨|ξ|²⟩ = 1 metric) and invariant to the working-coordinate Jacobian.
+[quantities.et_real]
+h5path = "FreeBoundaryStability/eigenmode_energies"
+type = "complex_vector"
+extract = "real_first"
+label = "total energy Re(et[1])"
+noise_threshold = 1e-10
+order = 10
+class = "physics_converged"
+
+[quantities.et_imag]
+h5path = "FreeBoundaryStability/eigenmode_energies"
+type = "complex_vector"
+extract = "imag_first"
+label = "total energy Im(et[1])"
+noise_threshold = 1e-10
+order = 11
+class = "physics_converged"
+
+[quantities.ep_real]
+h5path = "FreeBoundaryStability/eigenmode_plasma_energies"
+type = "complex_vector"
+extract = "real_first"
+label = "plasma energy Re(ep[1])"
+noise_threshold = 1e-10
+order = 12
+class = "physics_converged"
+
+[quantities.ev_real]
+h5path = "FreeBoundaryStability/eigenmode_vacuum_energies"
+type = "complex_vector"
+extract = "real_first"
+label = "vacuum energy Re(ev[1])"
+noise_threshold = 1e-10
+order = 13
+class = "physics_converged"
+
+[quantities.vacuum_eigenvalue]
+h5path = "FreeBoundaryStability/vacuum_eigenvalue"
+type = "real_scalar"
+extract = "value"
+label = "vacuum matrix min eigenvalue"
+noise_threshold = 1e-10
+order = 14
+class = "physics_converged"
+
+[quantities.et_all]
+h5path = "FreeBoundaryStability/eigenmode_energies"
+type = "complex_vector"
+extract = "all_complex"
+label = "total energy (all)"
+noise_threshold = 1e-10
+order = 20
+class = "physics_converged"
+
+# Integration
+[quantities.nstep]
+h5path = "integration/nstep"
+type = "int_scalar"
+extract = "value"
+label = "ODE steps (saved)"
+noise_threshold = 0
+order = 30
+class = "diagnostic"
+
+[quantities.nstep_total]
+h5path = "integration/nstep_total"
+type = "int_scalar"
+extract = "value"
+label = "ODE steps (total)"
+noise_threshold = 0
+order = 31
+class = "diagnostic"
+
+# Equilibrium — q95 and q0 pin the flux scale, which is solved from the requested on-axis q.
+[quantities.q0]
+h5path = "equil/q0"
+type = "real_scalar"
+extract = "value"
+label = "q0"
+noise_threshold = 1e-12
+order = 40
+class = "equilibrium_scalar"
+
+[quantities.q95]
+h5path = "equil/q95"
+type = "real_scalar"
+extract = "value"
+label = "q95"
+noise_threshold = 1e-12
+order = 41
+class = "equilibrium_scalar"
+
+[quantities.betat]
+h5path = "equil/betat"
+type = "real_scalar"
+extract = "value"
+label = "beta_t"
+noise_threshold = 1e-12
+order = 42
+class = "equilibrium_scalar"
+
+[quantities.betan]
+h5path = "equil/betan"
+type = "real_scalar"
+extract = "value"
+label = "beta_n"
+noise_threshold = 1e-12
+order = 43
+class = "equilibrium_scalar"
+
+# Singular surfaces — the outermost rational is set by where psihigh cuts the diverging q,
+# so these are the quantities that expose a change in truncation behaviour.
+[quantities.msing]
+h5path = "singular/msing"
+type = "int_scalar"
+extract = "value"
+label = "# singular surfaces"
+noise_threshold = 0
+order = 50
+class = "topological"
+
+[quantities.sing_psi]
+h5path = "singular/psi"
+type = "real_vector"
+extract = "all_real"
+label = "singular psi locations"
+noise_threshold = 1e-8
+order = 51
+class = "physics_converged"
+
+[quantities.sing_q]
+h5path = "singular/q"
+type = "real_vector"
+extract = "all_real"
+label = "singular q values"
+noise_threshold = 1e-8
+order = 52
+class = "physics_converged"
+
+# Mode numbers
+[quantities.mpert]
+h5path = "info/mpert"
+type = "int_scalar"
+extract = "value"
+label = "mpert"
+noise_threshold = 0
+order = 60
+class = "topological"
+
+[quantities.npert]
+h5path = "info/npert"
+type = "int_scalar"
+extract = "value"
+label = "npert"
+noise_threshold = 0
+order = 61
+class = "topological"
+
+# Profiles (checksums for fast comparison of large arrays)
+[quantities.q_profile]
+h5path = "splines/profiles/q"
+type = "real_vector"
+extract = "checksum"
+label = "q profile (checksum)"
+noise_threshold = 0
+order = 70
+class = "physics_converged"
+
+[quantities.pressure_profile]
+h5path = "splines/profiles/mu0p"
+type = "real_vector"
+extract = "checksum"
+label = "pressure profile (checksum)"
+noise_threshold = 0
+order = 71
+class = "physics_converged"
+
+# Runtime (special: not from H5)
+[quantities.runtime]
+h5path = ""
+type = "runtime"
+extract = "value"
+label = "Runtime (s)"
+noise_threshold = 0.0
+order = 999
+class = "diagnostic"