From 68de6992c08db6095782f4d39d705813738b57b5 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Mon, 19 Oct 2026 23:28:21 +0200 Subject: [PATCH 1/8] DOCS - BUG FIX - State the gauge of the pmodb lagb column and how the pmodb columns relate lagb (gpout_pmodb) is the Lagrangian field strength of xi_perp + lambda B with xi^zeta = 0 in the working (jac_type) chart. It is neither the xi_perp-only value nor Park et al. (2009) Eq. 7 with GPEC's divergence-free xi, and lagb = -(Bdivxprp + Bkxprp) is not an identity: lagb + Bdivxprp + Bkxprp = -(xi_perp^zeta/q) d_theta B - xi_perp.F/B, F = (curl B) x B - mu0 grad p. The docs called lagb "the Lagrangian |B| necessary for NTV post processing" and gpdiag labelled it "with parallel components". docs/pmodb_gauge.rst derives every column in the GPEC chart for any jac_type, the regularisation terms of reg_flag, and what PENTRC and perturbed-Boozer/Hamada consumers need. The identities are proved with the FortSym native zero test (docs/derivations/pmodb_gauge, verdicts JSON); FortSym is needed only to rerun the proofs. No code change. --- docs/derivations/pmodb_gauge/derive.py | 93 ++++ .../pmodb_gauge/fortsym_verdicts.json | 275 ++++++++++ docs/derivations/pmodb_gauge/identities.py | 481 ++++++++++++++++++ docs/outputs.rst | 4 +- docs/pmodb_gauge.rst | 172 +++++++ docs/site.rst | 1 + gpec/gpdiag.f | 4 +- gpec/gpout.f | 3 + 8 files changed, 1029 insertions(+), 4 deletions(-) create mode 100644 docs/derivations/pmodb_gauge/derive.py create mode 100644 docs/derivations/pmodb_gauge/fortsym_verdicts.json create mode 100644 docs/derivations/pmodb_gauge/identities.py create mode 100644 docs/pmodb_gauge.rst diff --git a/docs/derivations/pmodb_gauge/derive.py b/docs/derivations/pmodb_gauge/derive.py new file mode 100644 index 000000000..17337b82d --- /dev/null +++ b/docs/derivations/pmodb_gauge/derive.py @@ -0,0 +1,93 @@ +"""Run the pmodb gauge identities through FortSym and record the verdicts. + +Usage (FortSym checkout with a built C ABI library):: + + PYTHONPATH=$FORTSYM/python python3 docs/derivations/pmodb_gauge/derive.py + +Writes ``fortsym_verdicts.json`` beside this file. ``PROVED`` means FortSym's +native zero decision returned true; ``UNKNOWN`` and ``NONZERO`` fail an +identity and are the required outcome for a negative control. No numerical +probe is accepted as evidence. +""" + +from __future__ import annotations + +import json +import subprocess +import sys +import time +from pathlib import Path + +HERE = Path(__file__).resolve().parent +sys.path.insert(0, str(HERE)) + +import fortsym # noqa: E402 +import fortsym.sympy as fs # noqa: E402 +import identities # noqa: E402 + + +def verdict(expr) -> str: + if not hasattr(expr, "is_zero"): + return "PROVED" if expr == 0 else "NONZERO" + value = expr.is_zero + return {True: "PROVED", False: "NONZERO"}.get(value, "UNKNOWN") + + +def counterexample(expr): + """Exact rational point at which FortSym decides the residual is non-zero.""" + symbols = sorted(expr.free_symbols, key=str) + for attempt in range(1, 6): + point, value = {}, expr + for j, sym in enumerate(symbols): + num = 7 * attempt + 3 * j + 1 + point[str(sym)] = f"{num}/{5 + j}" + value = value.subs(sym, fs.Rational(num, 5 + j)) + if value.is_zero is False: + return point + return None + + +def revision() -> str: + root = Path(fortsym.__file__).resolve().parents[2] + try: + return subprocess.run(["git", "-C", str(root), "rev-parse", "HEAD"], + capture_output=True, text=True, check=True).stdout.strip() + except (OSError, subprocess.CalledProcessError): + return "unknown" + + +def main() -> int: + records, failed = [], [] + for label, equation, residual in identities.all_identities(fs): + start = time.perf_counter() + v = verdict(residual) + records.append({"label": label, "equation": equation, "verdict": v, + "seconds": round(time.perf_counter() - start, 3)}) + if v != "PROVED": + failed.append(label) + print(f"{label:6s} {v:8s} {equation}", flush=True) + for label, equation, residual in identities.open_gaps(fs): + v = verdict(residual) + records.append({"label": label, "equation": equation, "verdict": v, + "open_gap": True}) + print(f"{label:6s} {v:8s} (open FortSym gap) {equation}", flush=True) + for label, equation, residual in identities.negative_controls(fs): + v = verdict(residual) + records.append({"label": label, "equation": equation, "verdict": v, + "negative_control": True}) + witness = counterexample(residual) if v != "NONZERO" else None + if witness is not None: + records[-1]["counterexample"] = witness + if v == "PROVED" or (v != "NONZERO" and witness is None): + failed.append(label) + print(f"{label:6s} {v:8s} (negative control" + f"{', witness ' + str(witness) if witness else ''}) {equation}", flush=True) + out = {"engine": "fortsym native zero test via C ABI (fortsym.sympy)", + "fortsym_revision": revision(), "checks": records, "failed": failed} + (HERE / "fortsym_verdicts.json").write_text(json.dumps(out, indent=2) + "\n") + print("FAILED: " + " ".join(failed) if failed else "all identities proved") + return 1 if failed else 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/docs/derivations/pmodb_gauge/fortsym_verdicts.json b/docs/derivations/pmodb_gauge/fortsym_verdicts.json new file mode 100644 index 000000000..5f0f0e2fd --- /dev/null +++ b/docs/derivations/pmodb_gauge/fortsym_verdicts.json @@ -0,0 +1,275 @@ +{ + "engine": "fortsym native zero test via C ABI (fortsym.sympy)", + "fortsym_revision": "7e8ef7abdd90c49bb5f75e2c119e7d181656bcc6", + "checks": [ + { + "label": "G1.1", + "equation": "J dB^psi = GPEC bwp (pointwise)", + "verdict": "PROVED", + "seconds": 0.01 + }, + { + "label": "G1.2", + "equation": "J dB^theta = GPEC bwt (pointwise)", + "verdict": "PROVED", + "seconds": 0.004 + }, + { + "label": "G1.3", + "equation": "J dB^zeta = GPEC bwz (pointwise)", + "verdict": "PROVED", + "seconds": 0.007 + }, + { + "label": "G1.4", + "equation": "(xi x B)_psi = chi1 (q xi^theta - xi^zeta) = xi_s", + "verdict": "PROVED", + "seconds": 0.001 + }, + { + "label": "G1.5", + "equation": "xi_s(xi + lam B) = xi_s(xi)", + "verdict": "PROVED", + "seconds": 0.001 + }, + { + "label": "G1F.1", + "equation": "bwp = chi1 (m - n q) 2 pi i xsp", + "verdict": "PROVED", + "seconds": 0.001 + }, + { + "label": "G1F.2", + "equation": "bwt = -(chi1 xsp1 + 2 pi i n xss)", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "G1F.3", + "equation": "bwz = -(chi1 (q1 xsp + q xsp1) + 2 pi i m xss)", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "G2.1", + "equation": "(d_t + q d_z)(J xi^t) = d_z(J xi_s)/chi1 - d_psi(J xi^psi) + J div xi", + "verdict": "PROVED", + "seconds": 0.001 + }, + { + "label": "G2.2", + "equation": "(d_t + q d_z)(J xi^z) = -q d_psi(J xi^psi) - d_t(J xi_s)/chi1 + q J div xi", + "verdict": "PROVED", + "seconds": 0.005 + }, + { + "label": "G2.3", + "equation": "q xwt - xwz = J xi_s / chi1", + "verdict": "PROVED", + "seconds": 0.001 + }, + { + "label": "G3", + "equation": "B.dB = chi1 (bvt + q bvz)/J (eulb times |B|)", + "verdict": "PROVED", + "seconds": 0.013 + }, + { + "label": "G4.1", + "equation": "J xi_perp^theta = xmt - (chi1/B)^2 (xvt + q xvz)/J", + "verdict": "PROVED", + "seconds": 0.27 + }, + { + "label": "G4.2", + "equation": "J xi_perp^zeta = xmz - q (chi1/B)^2 (xvt + q xvz)/J", + "verdict": "PROVED", + "seconds": 0.119 + }, + { + "label": "G4.3", + "equation": "B_i xi_perp^i = 0 for xi_perp^i = xi^i - (B_j xi^j / B_k B^k) B^i", + "verdict": "PROVED", + "seconds": 0.012 + }, + { + "label": "G4.4", + "equation": "xi_perp^theta = xms/(chi1 q) + xmz/(J q) - (chi1/(J B))^2 (xvt + q xvz) (curvfun, gpout.f)", + "verdict": "PROVED", + "seconds": 0.183 + }, + { + "label": "G5.1", + "equation": "xi_g^zeta = 0 for lam = -xi_perp^zeta / B^zeta", + "verdict": "PROVED", + "seconds": 0.232 + }, + { + "label": "G5.2", + "equation": "xi_g^theta = xi_s/(chi1 q)", + "verdict": "PROVED", + "seconds": 0.757 + }, + { + "label": "G5.2b", + "equation": "xi_g^psi = xi^psi", + "verdict": "PROVED", + "seconds": 0.008 + }, + { + "label": "G5.3", + "equation": "curl(xi_g x B) = curl(xi_perp x B) (theta component)", + "verdict": "PROVED", + "seconds": 0.076 + }, + { + "label": "G5.4", + "equation": "lagb (gpout.f) = dB_L(xi_g), the xi^zeta = 0 gauge", + "verdict": "PROVED", + "seconds": 0.026 + }, + { + "label": "G5.5", + "equation": "lagb - dB_L(xi_perp) = -(xi_perp^zeta/q) d_theta B", + "verdict": "PROVED", + "seconds": 0.017 + }, + { + "label": "G5.6", + "equation": "llagb (gpdiag.f) = dB_L(xi_perp)", + "verdict": "PROVED", + "seconds": 0.038 + }, + { + "label": "G5.7", + "equation": "dB_L(xi, div-free) = lagb + (xmz/(J q)) d_theta B", + "verdict": "PROVED", + "seconds": 0.014 + }, + { + "label": "G5.8", + "equation": "xi_perp^zeta = -(B_theta xi_s/chi1 + q B_psi xi^psi)/(B_theta + q B_zeta)", + "verdict": "PROVED", + "seconds": 0.735 + }, + { + "label": "G7.1", + "equation": "J B^2 (xi_perp.grad B + mu0 p' xi^psi/B)/B = sband xsp + tband xms/chi1, xi = e_psi", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "G7.2", + "equation": "J B^2 (xi_perp.grad B + mu0 p' xi^psi/B)/B = sband xsp + tband xms/chi1, xi = e_theta", + "verdict": "PROVED", + "seconds": 0.007 + }, + { + "label": "G7.3", + "equation": "J B^2 (xi_perp.grad B + mu0 p' xi^psi/B)/B = sband xsp + tband xms/chi1, xi = e_zeta", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "C1", + "equation": "b.curl(xi_perp x B) + xi_perp.grad B = -B(div xi_perp + kappa.xi_perp)", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "C2", + "equation": "B kappa.xi_perp = xi_perp.grad B + xi_perp.((curl B) x B)/B", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "R1", + "equation": "J div xi (reg) = (1 - R)(J R xsp1 + J' xsp)", + "verdict": "PROVED", + "seconds": 0.014 + }, + { + "label": "R2", + "equation": "J div xi (R = 1) = 0", + "verdict": "PROVED", + "seconds": 0.028 + }, + { + "label": "R3", + "equation": "J div dB (reg) = 2 pi i chi1 (m - n q)(1 - R) xsp1", + "verdict": "PROVED", + "seconds": 0.001 + }, + { + "label": "K1", + "equation": "CGL energy density = (1/2)[(5/3) D^2 + (1/3)(D + 3K)^2]", + "verdict": "PROVED", + "seconds": 0.001 + }, + { + "label": "K3", + "equation": "v [drive change under xi -> xi + lam B - d/dl (v B lam)] = 0", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "K4", + "equation": "pmodb_file-path drive change = -mu B lam dB/dl / v", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "G6", + "equation": "B.dB + xi_perp.grad B^2 + B^2 div xi_perp + xi_perp.(curl B x B) = 0", + "verdict": "UNKNOWN", + "open_gap": true + }, + { + "label": "N1", + "equation": "lagb is the xi_perp-only Lagrangian |B| (false)", + "verdict": "UNKNOWN", + "negative_control": true, + "counterexample": { + "Bm_psi": "8/5", + "Bm_theta": "11/6", + "J_0": "14/7", + "chi1": "17/8", + "g_11_0": "20/9", + "g_12_0": "23/10", + "g_13_0": "26/11", + "g_22_0": "29/12", + "g_23_0": "32/13", + "q_0": "35/14", + "xp_0": "38/15", + "xt_0": "41/16", + "xz_0": "44/17" + } + }, + { + "label": "N2", + "equation": "lagb is Park 2009 Eq. 7 with the div-free xi (false)", + "verdict": "UNKNOWN", + "negative_control": true, + "counterexample": { + "Bm_psi": "8/5", + "Bm_theta": "11/6", + "chi1": "14/7", + "q_0": "17/8", + "xp_0": "20/9", + "xt_0": "23/10", + "xz_0": "26/11" + } + }, + { + "label": "K2neg", + "equation": "fcgl with kapx = -(dbob + divx)/2 equals the CGL density", + "verdict": "UNKNOWN", + "negative_control": true, + "counterexample": { + "D": "8/5", + "K": "11/6" + } + } + ], + "failed": [] +} diff --git a/docs/derivations/pmodb_gauge/identities.py b/docs/derivations/pmodb_gauge/identities.py new file mode 100644 index 000000000..d52b3bcd0 --- /dev/null +++ b/docs/derivations/pmodb_gauge/identities.py @@ -0,0 +1,481 @@ +"""Symbolic identities behind the GPEC pmodb |B| columns. + +Every function returns ``(label, equation, residual)`` triples; a residual must +vanish identically. The module is written against a SymPy-shaped namespace +``sp`` so that the same derivation runs on FortSym (``fortsym.sympy``, the +authority used by ``derive.py``). Labels refer to ``docs/pmodb_gauge.rst``. + +GPEC/DCON coordinates (psi, theta, zeta): angles in turns, general +straight-field-line Jacobian ``J(psi, theta)`` (any ``jac_type``), + + B = chi1 grad(zeta - q theta) x grad(psi), J B^theta = chi1, + J B^zeta = chi1 q, B^psi = 0, + +and the surface displacement (Glasser 2016 Eqs. 9 and 18) + + xi_s = (xi x B)_psi = chi1 (q xi^theta - xi^zeta). + +Pointwise identities are proved on jets: every field is an affine (1-jet) +polynomial in (psi, theta, zeta) about the origin with independent symbolic +coefficients, and the residual is evaluated at the origin after +differentiation. Each identity involves at most first derivatives of the +jetted fields, so it holds for all smooth fields iff it holds for all 1-jets. +The metric is general except for ``det g = J**2``, imposed by solving for +``g_33``. +""" + +from __future__ import annotations + +COORDS = ("psi", "theta", "zeta") + + +class Dual: + """First-order jet at the origin: value ``v`` and derivatives ``d[i]``. + + ``d is None`` marks a value whose own derivatives were not carried (a + derivative of a jet); differentiating it again raises. + """ + + def __init__(self, v, d=(0, 0, 0)): + self.v, self.d = v, (None if d is None else tuple(d)) + + @staticmethod + def lift(x): + return x if isinstance(x, Dual) else Dual(x) + + def _need(self): + if self.d is None: + raise ValueError("second derivative requested from a 1-jet") + return self.d + + def __add__(self, o): + o = Dual.lift(o) + if self.d is None or o.d is None: + return Dual(self.v + o.v, None) + return Dual(self.v + o.v, [a + b for a, b in zip(self.d, o.d)]) + + __radd__ = __add__ + + def __neg__(self): + return Dual(-self.v, None if self.d is None else [-a for a in self.d]) + + def __sub__(self, o): + return self + (-Dual.lift(o)) + + def __rsub__(self, o): + return Dual.lift(o) - self + + def __mul__(self, o): + o = Dual.lift(o) + if self.d is None or o.d is None: + return Dual(self.v * o.v, None) + return Dual(self.v * o.v, [self.v * b + a * o.v for a, b in zip(self.d, o.d)]) + + __rmul__ = __mul__ + + def inv(self): + if self.d is None: + return Dual(1 / self.v, None) + return Dual(1 / self.v, [-a / self.v ** 2 for a in self.d]) + + def __truediv__(self, o): + return self * Dual.lift(o).inv() + + def __rtruediv__(self, o): + return Dual.lift(o) * self.inv() + + def __pow__(self, k): + assert isinstance(k, int) and k >= 1 + out = self + for _ in range(k - 1): + out = out * self + return out + + +class Setup: + """Jetted GPEC fields (1-jets at the origin) in a SymPy-shaped namespace.""" + + def __init__(self, sp, metric=False): + self.sp = sp + self.u = (0, 1, 2) # psi, theta, zeta + self.chi1 = Dual(sp.Symbol("chi1", positive=True)) + self.q = self.jet("q", (0,)) + self.J = self.jet("J", (0, 1)) + # Contravariant displacement components (general, incl. xi_parallel). + self.xi = [self.jet(f"x{c[0]}", self.u) for c in COORDS] + self.Bup = [Dual(0), self.chi1 / self.J, self.chi1 * self.q / self.J] + if metric: + g = {} + for (i, j) in ((0, 0), (0, 1), (0, 2), (1, 1), (1, 2)): + g[(i, j)] = self.jet(f"g_{i + 1}{j + 1}", (0, 1)) + minor = g[(0, 0)] * g[(1, 1)] - g[(0, 1)] ** 2 + rest = (-g[(0, 0)] * g[(1, 2)] ** 2 + 2 * g[(0, 1)] * g[(1, 2)] * g[(0, 2)] + - g[(1, 1)] * g[(0, 2)] ** 2) + g[(2, 2)] = (self.J ** 2 - rest) / minor + self.g = [[g[(min(i, j), max(i, j))] for j in range(3)] for i in range(3)] + + def jet(self, name, deps): + sp = self.sp + return Dual(sp.Symbol(f"{name}_0"), + [sp.Symbol(f"{name}_{COORDS[i]}") if i in deps else 0 for i in range(3)]) + + def d(self, f, i): + f = Dual.lift(f) + return Dual(f._need()[i], None) + + def at0(self, e): + return Dual.lift(e).v + + # --- tensor algebra in the (psi, theta, zeta) chart ----------------------- + @staticmethod + def eps(i, j, k): + return (i - j) * (j - k) * (k - i) // 2 + + def cross_cov(self, a, b): + """(a x b)_i = J eps_ijk a^j b^k for contravariant a, b.""" + return [self.J * sum(self.eps(i, j, k) * a[j] * b[k] + for j in range(3) for k in range(3)) for i in range(3)] + + def curl_up(self, A): + """(curl A)^i = eps^ijk d_j A_k / J for covariant A.""" + return [sum(self.eps(i, j, k) * self.d(A[k], j) + for j in range(3) for k in range(3)) / self.J for i in range(3)] + + def div(self, v): + return sum(self.d(self.J * v[i], i) for i in range(3)) / self.J + + def lower(self, v): + return [sum(self.g[i][j] * v[j] for j in range(3)) for i in range(3)] + + def dot(self, a, b): + return sum(a[i] * b[i] for i in range(3)) + + def xi_s(self, v): + return self.chi1 * (self.q * v[1] - v[2]) + + +def induction(sp): + """G1: GPEC's J dB^i (gpeq.f:107-109) is J curl(xi x B)^i for any Jacobian.""" + S = Setup(sp) + xi, B = S.xi, S.Bup + JdB = [S.J * c for c in S.curl_up(S.cross_cov(xi, B))] + xs = S.xi_s(xi) + d = S.d + gpec = [S.chi1 * (d(xi[0], 1) + S.q * d(xi[0], 2)), + d(xs, 2) - S.chi1 * d(xi[0], 0), + -d(xs, 1) - S.chi1 * d(S.q * xi[0], 0)] + out = [(f"G1.{i + 1}", f"J dB^{COORDS[i]} = GPEC bw{'ptz'[i]} (pointwise)", + S.at0(JdB[i] - gpec[i])) for i in range(3)] + # (xi x B)_psi = xi_s: the surface displacement of Glasser (2016) Eq. 9. + out.append(("G1.4", "(xi x B)_psi = chi1 (q xi^theta - xi^zeta) = xi_s", + S.at0(S.cross_cov(xi, B)[0] - xs))) + # A field-aligned displacement changes neither xi^psi nor xi_s. + lam = S.jet("lam", S.u) + shifted = [xi[i] + lam * B[i] for i in range(3)] + out.append(("G1.5", "xi_s(xi + lam B) = xi_s(xi)", + S.at0(S.xi_s(shifted) - xs))) + return out + + +def fourier_symbols(sp): + """G1F: Fourier symbols of G1 for xi ~ exp(2 pi i (m theta - n zeta)).""" + I, pi = sp.I, sp.pi + m, n, chi1 = sp.Symbol("m"), sp.Symbol("n"), sp.Symbol("chi1", positive=True) + q, q1 = sp.Symbol("q"), sp.Symbol("q1") + X, X1, Xs = sp.Symbol("xsp"), sp.Symbol("xsp1"), sp.Symbol("xss") + dth, dze = 2 * pi * I * m, -2 * pi * I * n + out = [] + out.append(("G1F.1", "bwp = chi1 (m - n q) 2 pi i xsp", + chi1 * (dth + q * dze) * X - chi1 * (m - n * q) * 2 * pi * I * X)) + out.append(("G1F.2", "bwt = -(chi1 xsp1 + 2 pi i n xss)", + (dze * Xs - chi1 * X1) + (chi1 * X1 + 2 * pi * I * n * Xs))) + out.append(("G1F.3", "bwz = -(chi1 (q1 xsp + q xsp1) + 2 pi i m xss)", + (-dth * Xs - chi1 * (q1 * X + q * X1)) + + (chi1 * (q1 * X + q * X1) + 2 * pi * I * m * Xs))) + return out + + +def incompressible_closure(sp): + """G2: GPEC's xi^theta, xi^zeta (gpeq.f gpeq_contra) follow from div xi = 0.""" + S = Setup(sp) + xi, J, q, chi1, d = S.xi, S.J, S.q, S.chi1, S.d + D = sum(d(J * xi[i], i) for i in range(3)) # J div xi + xs = S.xi_s(xi) + Jt, Jz = J * xi[1], J * xi[2] + r1 = (d(Jt, 1) + q * d(Jt, 2)) - d(J * xs, 2) / chi1 + d(J * xi[0], 0) - D + r2 = (d(Jz, 1) + q * d(Jz, 2)) + q * d(J * xi[0], 0) + d(J * xs, 1) / chi1 - q * D + r3 = q * Jt - Jz - J * xs / chi1 + return [("G2.1", "(d_t + q d_z)(J xi^t) = d_z(J xi_s)/chi1 - d_psi(J xi^psi) + J div xi", + S.at0(r1)), + ("G2.2", "(d_t + q d_z)(J xi^z) = -q d_psi(J xi^psi) - d_t(J xi_s)/chi1" + " + q J div xi", S.at0(r2)), + ("G2.3", "q xwt - xwz = J xi_s / chi1", S.at0(r3))] + + +def columns(sp): + """G3-G5: eulb, xi_perp components, the gauge of lagb, and llagb.""" + S = Setup(sp, metric=True) + xi, B, J, q, chi1 = S.xi, S.Bup, S.J, S.q, S.chi1 + Bm = S.jet("Bm", S.u[:2]) # |B|(psi, theta), axisymmetric + dB = S.curl_up(S.cross_cov(xi, B)) + dBcov = S.lower(dB) + B2 = S.dot(S.lower(B), B) + xs = S.xi_s(xi) + xv = S.lower(xi) # GPEC xvt, xvz (covariant) + xmt, xmz = J * xi[1], J * xi[2] # GPEC xmt, xmz (J xi^i) + sigma = S.dot(xv, B) / B2 + xperp = [xi[i] - sigma * B[i] for i in range(3)] + out = [] + out.append(("G3", "B.dB = chi1 (bvt + q bvz)/J (eulb times |B|)", + S.at0(S.dot(dBcov, B) - chi1 * (dBcov[1] + q * dBcov[2]) / J))) + out.append(("G4.1", "J xi_perp^theta = xmt - (chi1/B)^2 (xvt + q xvz)/J", + S.at0(J * xperp[1] - (xmt - chi1 ** 2 / B2 * (xv[1] + q * xv[2]) / J)))) + out.append(("G4.2", "J xi_perp^zeta = xmz - q (chi1/B)^2 (xvt + q xvz)/J", + S.at0(J * xperp[2] - (xmz - q * chi1 ** 2 / B2 * (xv[1] + q * xv[2]) / J)))) + # B . xi_perp = 0 as a statement about arbitrary vectors: the metric-jet + # form exceeds FortSym's polynomial-cancellation limit (issue #84), so the + # projection is proved on independent covariant B_i, contravariant B^i, xi^i. + Bl = [sp.Symbol(f"Bl_{i}") for i in range(3)] + Bu = [sp.Symbol(f"Bu_{i}") for i in range(3)] + xu = [sp.Symbol(f"xu_{i}") for i in range(3)] + sig = sum(Bl[i] * xu[i] for i in range(3)) / sum(Bl[i] * Bu[i] for i in range(3)) + out.append(("G4.3", "B_i xi_perp^i = 0 for xi_perp^i = xi^i - (B_j xi^j / B_k B^k) B^i", + sum(Bl[i] * (xu[i] - sig * Bu[i]) for i in range(3)))) + coef = xs / (chi1 * q) + xmz / (J * q) - chi1 ** 2 / (J ** 2 * B2) * (xv[1] + q * xv[2]) + out.append(("G4.4", "xi_perp^theta = xms/(chi1 q) + xmz/(J q)" + " - (chi1/(J B))^2 (xvt + q xvz) (curvfun, gpout.f)", + S.at0(xperp[1] - coef))) + # Gauge of lagb: xi_g = xi_perp + lam B with xi_g^zeta = 0. + lam = -xperp[2] / B[2] + xg = [xperp[i] + lam * B[i] for i in range(3)] + out.append(("G5.1", "xi_g^zeta = 0 for lam = -xi_perp^zeta / B^zeta", S.at0(xg[2]))) + out.append(("G5.2", "xi_g^theta = xi_s/(chi1 q)", S.at0(xg[1] - xs / (chi1 * q)))) + out.append(("G5.2b", "xi_g^psi = xi^psi", S.at0(xg[0] - xi[0]))) + dBg = S.curl_up(S.cross_cov(xg, B)) + dBp = S.curl_up(S.cross_cov(xperp, B)) + out.append(("G5.3", "curl(xi_g x B) = curl(xi_perp x B) (theta component)", + S.at0(dBg[1] - dBp[1]))) + eul = S.dot(dBcov, B) # |B| * eulb + grad = lambda v: sum(v[i] * S.d(Bm, i) for i in range(3)) + lagb = eul + Bm * (xi[0] * S.d(Bm, 0) + xs / (chi1 * q) * S.d(Bm, 1)) + lagb_g = eul + Bm * grad(xg) + llagb = eul + Bm * grad(xperp) + llagb_code = lagb + Bm * S.d(Bm, 1) * (xmz / (J * q) - chi1 ** 2 / (J ** 2 * B2) + * (xv[1] + q * xv[2])) + out.append(("G5.4", "lagb (gpout.f) = dB_L(xi_g), the xi^zeta = 0 gauge", + S.at0(lagb - lagb_g))) + out.append(("G5.5", "lagb - dB_L(xi_perp) = -(xi_perp^zeta/q) d_theta B", + S.at0(lagb - llagb + Bm * xperp[2] / q * S.d(Bm, 1)))) + out.append(("G5.6", "llagb (gpdiag.f) = dB_L(xi_perp)", S.at0(llagb_code - llagb))) + # Park (2009) Eq. 7 with the div-free xi (GPEC's own xmt, xmz). + lagb_div_free = eul + Bm * (xi[0] * S.d(Bm, 0) + xi[1] * S.d(Bm, 1)) + out.append(("G5.7", "dB_L(xi, div-free) = lagb + (xmz/(J q)) d_theta B", + S.at0(lagb_div_free - lagb - Bm * xmz / (J * q) * S.d(Bm, 1)))) + # Boozer-type closed form: B_i xi_perp^i = 0 and the xi_s definition. + Bc = S.lower(B) + out.append(("G5.8", "xi_perp^zeta = -(B_theta xi_s/chi1 + q B_psi xi^psi)" + "/(B_theta + q B_zeta)", + S.at0(xperp[2] + (Bc[1] * xs / chi1 + q * Bc[0] * xi[0]) + / (Bc[1] + q * Bc[2])))) + return out + + +def force_balance_identity(sp): + """G6: |B| times (llagb + Bdivxprp + Bkxprp) = -xi_perp . F, F = (curl B) x B - mu0 grad p. + + With pmodb's Bkxprp = xi_perp.grad B + mu0 p' xi^psi / B (force-balance form + of B kappa.xi_perp, gpout.f curvfun) and Bdivxprp = B div xi_perp, this is + B.dB + xi_perp.grad(B^2) + B^2 div xi_perp + xi_perp.((curl B) x B) = 0. + """ + S = Setup(sp, metric=True) + xi, B = S.xi, S.Bup + Bc = S.lower(B) + B2 = S.dot(Bc, B) + sigma = S.dot(S.lower(xi), B) / B2 + xperp = [xi[i] - sigma * B[i] for i in range(3)] + dB = S.curl_up(S.cross_cov(xi, B)) + curlB = S.curl_up(Bc) + jxb = S.cross_cov(curlB, B) # ((curl B) x B)_i + res = (S.dot(S.lower(dB), B) + sum(xperp[i] * S.d(B2, i) for i in range(3)) + + B2 * S.div(xperp) + S.dot(jxb, xperp)) + return [("G6", "B.dB + xi_perp.grad B^2 + B^2 div xi_perp + xi_perp.(curl B x B) = 0", + S.at0(res))] + + +def pentrc_kappa_matrices(sp): + """G7: idcon.f fmodb sband/tband give J B^2 (Bkxprp/B) from xi^psi and xi^alpha. + + PENTRC set_peq (inputs.f90) and GPEC.jl set_perturbation_data! form + jbb_kapx = smat xsp + tmat xi_alpha with xi_alpha = xms/chi1 (gpout_xclebsch). + The statement is algebraic at a point (only the values of d_psi B^2 and + d_theta B^2 enter), so the metric, Jacobian and derivatives are independent + symbols here: + sband = J (mu0 p' + d_psi B^2/2) - chi1^2 (d_theta B^2/(2 B^2)) (g12 + q g13)/J, + tband = chi1^2 (d_theta B^2/(2 B^2)) (g23 + q g33)/J. + """ + chi1, q, J, P, Dp, Dt = sp.symbols("chi1 q J P1 Dp Dt") + # B^2 = B^i B_i = (chi1/J)(B_theta + q B_zeta) needs no det g = J^2, so the + # metric components are independent here. + g11, g12, g13, g22, g23, g33 = sp.symbols("g11 g12 g13 g22 g23 g33") + g = [[g11, g12, g13], [g12, g22, g23], [g13, g23, g33]] + xi = list(sp.symbols("xp xt xz")) + B = [0, chi1 / J, chi1 * q / J] + Bl = [sum(g[i][j] * B[j] for j in range(3)) for i in range(3)] + B2 = sum(Bl[i] * B[i] for i in range(3)) + sigma = sum(Bl[i] * xi[i] for i in range(3)) / B2 + xperp = [xi[i] - sigma * B[i] for i in range(3)] + xs = chi1 * (q * xi[1] - xi[2]) + lhs = J * (xperp[0] * (P + Dp) + xperp[1] * Dt) + sband = J * (P + Dp) - chi1 ** 2 * (Dt / B2) * (g12 + q * g13) / J + tband = chi1 ** 2 * (Dt / B2) * (g23 + q * g33) / J + res = lhs - sband * xi[0] - tband * xs / chi1 + # The residual is linear in xi; proving it for each basis displacement keeps + # every zero test within FortSym's 12-variable polynomial limit (issue #84). + out = [] + for k, name in enumerate(("psi", "theta", "zeta")): + r = res + for j in range(3): + r = r.subs(xi[j], 1 if j == k else 0) + out.append((f"G7.{k + 1}", "J B^2 (xi_perp.grad B + mu0 p' xi^psi/B)/B" + f" = sband xsp + tband xms/chi1, xi = e_{name}", r)) + return out + + +def cartesian_kappa(sp): + """C1, C2: Cartesian 1-jet forms of the kinematic and curvature identities.""" + x, y, z = sp.symbols("x y z") + X = (x, y, z) + beta = sp.Symbol("beta", positive=True) + g = [[sp.Symbol(f"h_{i}{j}") for j in range(3)] for i in range(3)] + B = [[0, 0, beta][i] + sum(g[i][j] * X[j] for j in range(3)) for i in range(3)] + d = lambda f, c: f.diff(c) if hasattr(f, "diff") else 0 + dot = lambda a, b: sum(a[i] * b[i] for i in range(3)) + cross = lambda a, b: [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], + a[0] * b[1] - a[1] * b[0]] + curl = lambda v: [d(v[2], y) - d(v[1], z), d(v[0], z) - d(v[2], x), + d(v[1], x) - d(v[0], y)] + div = lambda v: sum(d(v[i], X[i]) for i in range(3)) + at0 = lambda e: e.subs(x, 0).subs(y, 0).subs(z, 0) + Bm = sp.sqrt(dot(B, B)) + b = [c / Bm for c in B] + kappa = [sum(b[j] * d(b[i], X[j]) for j in range(3)) for i in range(3)] + a = [sp.Symbol(f"a_{i}") + sum(sp.Symbol(f"a_{i}{j}") * X[j] for j in range(3)) + for i in range(3)] + xp = cross(B, a) + gradB = lambda v: sum(v[j] * d(Bm, X[j]) for j in range(3)) + dBpar = dot(b, curl(cross(xp, B))) + return [("C1", "b.curl(xi_perp x B) + xi_perp.grad B = -B(div xi_perp + kappa.xi_perp)", + at0(dBpar + gradB(xp) + Bm * (div(xp) + dot(kappa, xp)))), + ("C2", "B kappa.xi_perp = xi_perp.grad B + xi_perp.((curl B) x B)/B", + at0(Bm * dot(kappa, xp) - gradB(xp) - dot(xp, cross(curl(B), B)) / Bm))] + + +def regularisation(sp): + """R1-R3: Park-regularised components in GPEC (reg_flag) on one harmonic. + + Diagonal model: one harmonic, J independent of theta. R = s^2/(s^2 + reg^2), + s = m - n q. gpeq_sol uses xmp1 = R xsp1 in bmt/bmz and xms; gpeq_contra + builds xwt/xwz from xmp1 and multiplies by R again; gpout_pmodb uses + R xsp1 explicitly in divxprp. + """ + I, pi = sp.I, sp.pi + m, n, q, chi1 = sp.symbols("m n q chi1") + J0, J1 = sp.symbols("J_0 J_1") # J and dJ/dpsi (theta-independent) + X, X1, Xs, R = sp.symbols("xsp xsp1 xms R") + s = m - n * q + xwt = -(J0 * R * X1 + J1 * X + 2 * pi * I * n / chi1 * J0 * Xs) / (2 * pi * I * s) + xwz = -(q * J0 * R * X1 + q * J1 * X + 2 * pi * I * m / chi1 * J0 * Xs) / (2 * pi * I * s) + xmt, xmz = R * xwt, R * xwz + jdiv = J0 * R * X1 + J1 * X + 2 * pi * I * m * xmt - 2 * pi * I * n * xmz + out = [("R1", "J div xi (reg) = (1 - R)(J R xsp1 + J' xsp)", + jdiv - (1 - R) * (J0 * R * X1 + J1 * X))] + jdiv_unreg = (J0 * X1 + J1 * X + 2 * pi * I * m * xwt.subs(R, 1) + - 2 * pi * I * n * xwz.subs(R, 1)) + out.append(("R2", "J div xi (R = 1) = 0", jdiv_unreg)) + # Divergence of the regularised Eulerian field: J div(dB) = d_psi bwp + ... . + bwp1 = 2 * pi * I * chi1 * (s * X1 - n * sp.Symbol("q1") * X) + bmt = -(chi1 * R * X1 + 2 * pi * I * n * Xs) + bmz = -(chi1 * (sp.Symbol("q1") * X + q * R * X1) + 2 * pi * I * m * Xs) + out.append(("R3", "J div dB (reg) = 2 pi i chi1 (m - n q)(1 - R) xsp1", + bwp1 + 2 * pi * I * m * bmt - 2 * pi * I * n * bmz + - 2 * pi * I * chi1 * s * (1 - R) * X1)) + return out + + +def cgl(sp): + """K1-K2: the PENTRC/GPEC.jl fcgl integrand against the CGL closure. + + D = div xi_perp, K = kappa.xi_perp, dB/B = -(D + K) (C1). CGL: + dp_perp/p = -D + dB/B, dp_par/p = -3D - 2 dB/B (no xi_par). Logan (2013) + Eq. 4 energy density, -(1/2)[(dp_par - dp_perp) dB/B + dp_par D]/p. + """ + D, K = sp.symbols("D K") + b = -(D + K) + dpe, dpa = -D + b, -3 * D - 2 * b + w = -((dpa - dpe) * b + dpa * D) / 2 + code = lambda k: (sp.Rational(5, 3) * D ** 2 + sp.Rational(1, 3) * (D + 3 * k) ** 2) / 2 + return [("K1", "CGL energy density = (1/2)[(5/3) D^2 + (1/3)(D + 3K)^2]", w - code(K)), + ("K2neg", "fcgl with kapx = -(dbob + divx)/2 equals the CGL density", + w - code(-(b + D) / 2))] + + +def drive_gauge(sp): + """K3, K4: gauge (in)variance of the Logan (2013) Eq. 5 bounce drive. + + Along a field line (arc length l, M = 1): v = sqrt(2(E - mu B)), integrand + (1/v)[(2E - 3 mu B) dB/B + (2E - 2 mu B) div xi]. Under xi -> xi + lam B, + dB -> dB + lam B dB/dl and div xi -> div xi + B dlam/dl (C1, Park 2011 Eq. 4). + """ + l, E, mu = sp.Symbol("l"), sp.Symbol("E", positive=True), sp.Symbol("mu", positive=True) + b0, b1, l0, l1 = sp.symbols("b_0 b_1 lam_0 lam_1") + B = b0 + b1 * l + lam = l0 + l1 * l + v = sp.sqrt(2 * (E - mu * B)) + dBdl, dlam = B.diff(l), lam.diff(l) + change = ((2 * E - 3 * mu * B) * lam * dBdl + (2 * E - 2 * mu * B) * B * dlam) / v + total = (v * B * lam).diff(l) + # pmodb_file path: dbob = lagb/B (gauge xi_g), divx = -(lagb + Bkxprp)/B. + change_pmodb = ((2 * E - 3 * mu * B) * lam * dBdl - (2 * E - 2 * mu * B) * lam * dBdl) / v + # Multiplied by v > 0 (away from turning points): FortSym leaves the + # unmultiplied sqrt quotient UNKNOWN. + return [("K3", "v [drive change under xi -> xi + lam B - d/dl (v B lam)] = 0", + ((change - total) * v).subs(l, 0)), + ("K4", "pmodb_file-path drive change = -mu B lam dB/dl / v", + (change_pmodb + mu * B * lam * dBdl / v).subs(l, 0))] + + +def all_identities(sp): + return (induction(sp) + fourier_symbols(sp) + incompressible_closure(sp) + + columns(sp) + pentrc_kappa_matrices(sp) + + cartesian_kappa(sp) + regularisation(sp) + + [c for c in cgl(sp) if not c[0].endswith("neg")] + drive_gauge(sp)) + + +def open_gaps(sp): + """Identities FortSym does not decide; their content is covered elsewhere. + + G6 is the coordinate form of C1 + C2 (proved on Cartesian 1-jets) with + det g = J^2 imposed through g_33. The zero test exceeds the native + polynomial-term limit (lazy-fortran/fortsym#84 class); it is recorded, not + counted as proved. + """ + return force_balance_identity(sp) + + +def negative_controls(sp): + """Residuals that must NOT vanish (the claims they encode are false).""" + S = Setup(sp, metric=True) + xi, B, J, q, chi1 = S.xi, S.Bup, S.J, S.q, S.chi1 + Bm = S.jet("Bm", S.u[:2]) + B2 = S.dot(S.lower(B), B) + sigma = S.dot(S.lower(xi), B) / B2 + xperp = [xi[i] - sigma * B[i] for i in range(3)] + xs = S.xi_s(xi) + lagb_adv = xi[0] * S.d(Bm, 0) + xs / (chi1 * q) * S.d(Bm, 1) + perp_adv = sum(xperp[i] * S.d(Bm, i) for i in range(3)) + out = [("N1", "lagb is the xi_perp-only Lagrangian |B| (false)", + S.at0(lagb_adv - perp_adv))] + out.append(("N2", "lagb is Park 2009 Eq. 7 with the div-free xi (false)", + S.at0(lagb_adv - (xi[0] * S.d(Bm, 0) + xi[1] * S.d(Bm, 1))))) + out += [c for c in cgl(sp) if c[0].endswith("neg")] + return out diff --git a/docs/outputs.rst b/docs/outputs.rst index a0d80b2ad..f642bd747 100644 --- a/docs/outputs.rst +++ b/docs/outputs.rst @@ -174,11 +174,11 @@ GPEC_singfld_n#.out **Info** Additional section showing the overlap field (Tesla) and overlap percentage for each eigenmode in the singcoup_svd output. See :ref:`coupling_overlap` for definitions. -GPEC_pmod_n#.out +GPEC_pmodb_n#.out **Flag** pmodb_flag - **Info** Eulerian and Lagrangian \left|\mathbf{B}\right|(real and imaginary) for each poloidal mode number at each value of \Psi_{N} output. This output is necessary for NTV post processing. + **Info** Perturbed field strength (real and imaginary, Tesla) for each poloidal mode number at each value of \Psi_{N} output: the Eulerian \delta B (eulb), the Lagrangian \delta B of the displacement with zero \xi^\zeta in the working (jac_type) coordinates (lagb), B\nabla\cdot\xi_\perp (Bdivxprp) and the force-balance form of B\kappa\cdot\xi_\perp (Bkxprp). lagb is not the \xi_\perp-only Lagrangian \delta B; that is -(Bdivxprp+Bkxprp) up to the equilibrium force-balance residual. See :doc:`pmodb_gauge` before using these columns for NTV post processing. GPEC_xbnormal_n#.out diff --git a/docs/pmodb_gauge.rst b/docs/pmodb_gauge.rst new file mode 100644 index 000000000..577e954d3 --- /dev/null +++ b/docs/pmodb_gauge.rst @@ -0,0 +1,172 @@ +.. _pmodb_gauge: + +****************************************** +Perturbed field strength (pmodb) and gauge +****************************************** + +This page states what the columns of ``gpec_pmodb_n#.out`` are, how they are +related, and which of them a field-strength consumer needs. Every identity +quoted here is proved with the FortSym native zero test by +``docs/derivations/pmodb_gauge/derive.py``; the verdicts are stored in +``docs/derivations/pmodb_gauge/fortsym_verdicts.json`` and the labels (G1, +G5.5, ...) refer to that file. + +.. contents:: On this page + :local: + :depth: 1 + + +Coordinates and displacement +============================ + +GPEC works in the DCON chart :math:`(\psi,\theta,\zeta)` with angles in turns, +any straight-field-line Jacobian :math:`\mathcal{J}(\psi,\theta)` selected by +``jac_type``, :math:`\mathcal{J}B^\theta=\chi_1`, +:math:`\mathcal{J}B^\zeta=\chi_1 q` and Fourier factor +:math:`\exp[2\pi i(m\theta-n\zeta)]`. The Euler-Lagrange solution supplies +:math:`\xi^\psi`, its radial derivative and the surface displacement +(Glasser 2016, Eqs. 9 and 18) + +.. math:: + + \xi_s = (\boldsymbol\xi\times\mathbf B)_\psi = \chi_1(q\xi^\theta-\xi^\zeta), + +the in-surface displacement across :math:`\mathbf B` (``xms``; +``gpec_xclebsch`` writes :math:`\xi^\alpha=\xi_s/\chi_1`). A field-aligned +displacement changes neither :math:`\xi^\psi` nor :math:`\xi_s` (G1.5). The +parallel part is fixed in ``gpeq_contra`` by +:math:`\nabla\cdot\boldsymbol\xi=0` (G2), which gives ``xmt`` +:math:`=\mathcal J\xi^\theta` and ``xmz`` :math:`=\mathcal J\xi^\zeta`. + +For every ``jac_type`` (G1, G1F) the Jacobian-weighted perturbed field of +``gpeq_sol`` is :math:`\mathcal J\,\nabla\times(\boldsymbol\xi\times\mathbf B)`. + + +The pmodb columns +================= + +``eulb`` + :math:`\delta B_E=\hat{\mathbf b}\cdot\delta\mathbf B`, the Eulerian + field strength (G3). + +``lagb`` + :math:`\delta B_E+\xi^\psi\partial_\psi B+\xi_s/(\chi_1 q)\,\partial_\theta B`. + This is the Lagrangian field strength of the displacement + :math:`\boldsymbol\xi_g=\boldsymbol\xi_\perp+\lambda\mathbf B` with + :math:`\xi_g^\zeta=0` in the *working* chart (``jac_type``) (G5.1-G5.4). + It is neither the :math:`\boldsymbol\xi_\perp`-only value (negative + control N1) nor Park et al. (2009) Eq. 7 with the divergence-free + :math:`\boldsymbol\xi` that GPEC constructs (N2; G5.7 gives that value as + ``lagb`` :math:`+\,\xi^\zeta/q\,\partial_\theta B`). ``jac_out`` only + relabels the harmonics; it does not change this gauge. + +``Bdivxprp`` + :math:`B\,\nabla\cdot\boldsymbol\xi_\perp`, with + :math:`\boldsymbol\xi_\perp=\boldsymbol\xi-(\boldsymbol\xi\cdot\mathbf B)\mathbf B/B^2` + built from ``xmt``, ``xmz`` and the covariant ``xvt``, ``xvz`` (G4.1-G4.3). + +``Bkxprp`` + :math:`\boldsymbol\xi_\perp\cdot\nabla B+\mu_0p'\xi^\psi/B` (G4.4). This is + :math:`B\boldsymbol\kappa\cdot\boldsymbol\xi_\perp` only through + equilibrium force balance: exactly, + :math:`B\boldsymbol\kappa\cdot\boldsymbol\xi_\perp=\boldsymbol\xi_\perp\cdot\nabla B + +\boldsymbol\xi_\perp\cdot[(\nabla\times\mathbf B)\times\mathbf B]/B` (C2). + The PENTRC/GPEC.jl matrices ``smat``, ``tmat`` (``idcon_build``) form the + same quantity times :math:`\mathcal J B^2` from :math:`\xi^\psi` and + :math:`\xi^\alpha` (G7). + +``llagb`` (``gpec_diagnostics_pmodb``, ``pmodbmn_flag``) + :math:`\delta B_L(\boldsymbol\xi_\perp)=\delta B_E+\boldsymbol\xi_\perp\cdot\nabla B` + (G5.6). ``lllagb`` is :math:`-(`\ ``Bdivxprp``\ :math:`+`\ ``Bkxprp``\ :math:`)`. + + +Relations between the columns +============================= + +Gauge of ``lagb`` (G5.5, G5.8): + +.. math:: + + \mathrm{lagb}=\delta B_L(\boldsymbol\xi_\perp)-\frac{\xi_\perp^\zeta}{q}\partial_\theta B, + \qquad + \xi_\perp^\zeta=-\frac{B_\theta\xi_s/\chi_1+qB_\psi\xi^\psi}{B_\theta+qB_\zeta}. + +Kinematic identity (C1; any :math:`\boldsymbol\xi_\perp`, no force balance) +and its pmodb form (C2 with G4; the same statement written directly in the +GPEC chart with :math:`\det g=\mathcal J^2`, G6, exceeds the FortSym +polynomial-term limit and is recorded as an open gap, not as proved): + +.. math:: + + \delta B_L(\boldsymbol\xi_\perp)=-B(\nabla\cdot\boldsymbol\xi_\perp + +\boldsymbol\kappa\cdot\boldsymbol\xi_\perp) + =-(\mathrm{Bdivxprp}+\mathrm{Bkxprp})-\boldsymbol\xi_\perp\cdot\mathbf F/B, + \qquad \mathbf F=(\nabla\times\mathbf B)\times\mathbf B-\mu_0\nabla p. + +Hence + +.. math:: + + \mathrm{lagb}+\mathrm{Bdivxprp}+\mathrm{Bkxprp} + =-\frac{\xi_\perp^\zeta}{q}\partial_\theta B-\frac{\boldsymbol\xi_\perp\cdot\mathbf F}{B}. + +``lagb`` :math:`=-(`\ ``Bdivxprp``\ :math:`+`\ ``Bkxprp``\ :math:`)` is +therefore not an identity. The first term is a gauge term of first order in +the perturbation; it is largest where :math:`B_\theta/(qB_\zeta)`, the +Boozer :math:`B_\psi` and :math:`\partial_\theta B` are large, i.e. toward +the edge. The second term vanishes for an exact equilibrium. + +Regularisation (``reg_flag``, Park et al. 2009 Eq. 10). ``gpeq_sol`` uses +:math:`R\,\partial_\psi\xi^\psi`, :math:`R=s^2/(s^2+\sigma^2)`, +:math:`s=m-nq`, in the field and in ``xms``; ``gpeq_contra`` multiplies +``xwt``, ``xwz`` by :math:`R` again. On one harmonic with +:math:`\theta`-independent :math:`\mathcal J` (R1-R3): + +.. math:: + + \mathcal J\nabla\cdot\boldsymbol\xi=(1-R)(\mathcal J R\,\partial_\psi\xi^\psi+\mathcal J'\xi^\psi),\qquad + \mathcal J\nabla\cdot\delta\mathbf B=2\pi i\chi_1 s(1-R)\,\partial_\psi\xi^\psi. + +Inside :math:`|m-nq|\lesssim\sigma` the regularised displacement is neither +divergence free nor does it produce a divergence-free field, so +``Bdivxprp`` (which carries the explicit :math:`R\,\partial_\psi\xi^\psi`) +and ``lagb`` (which uses only ``xms``) are not related by the identities +above on the resonant harmonics. + + +What consumers need +=================== + +PENTRC, GPEC.jl ``KineticForces`` + The drive of Logan et al. (2013) Eqs. 4-5 needs :math:`\delta B` and + :math:`\nabla\cdot\boldsymbol\xi` of the *same* displacement; the + combination is then invariant under :math:`\boldsymbol\xi\to\boldsymbol\xi+\lambda\mathbf B`. + The default PENTRC path (``peq_file`` = ``gpec_xclebsch``, ``set_peq``) + uses :math:`\delta B/B=-(\nabla\cdot\boldsymbol\xi_\perp+\boldsymbol\kappa\cdot\boldsymbol\xi_\perp)` + and :math:`\nabla\cdot\boldsymbol\xi_\perp` from the ``smat``-``zmat`` + matrices, i.e. a consistent :math:`\boldsymbol\xi_\perp` pair. The + ``pmodb_file`` path (``read_pmodb``) pairs ``lagb`` with + :math:`-(\mathrm{lagb}+\mathrm{Bkxprp})/B`, which is not the divergence of + :math:`\boldsymbol\xi_g`. Under :math:`\boldsymbol\xi\to\boldsymbol\xi+\lambda\mathbf B` + the consistent drive changes by a total derivative along the field line + (K3), while this pairing changes it by + :math:`-\mu B\lambda\,\partial_l B/v_\parallel` (K4). + +Bounce-averaged codes in perturbed Boozer or Hamada coordinates (NEO-RT, NEO-2, Shaing-type) + need :math:`\delta B` in the perturbed straight-field-line coordinates. + That follows from :math:`\delta B_L(\boldsymbol\xi_\perp)` (``llagb``, or + :math:`-(\mathrm{Bdivxprp}+\mathrm{Bkxprp})` away from rational surfaces) + and :math:`\nabla\cdot\boldsymbol\xi_\perp` by a magnetic differential + equation for the parallel gauge. ``lagb`` must not be used as + :math:`\delta B_L(\boldsymbol\xi_\perp)`. + + +Reproducing the proofs +====================== + +:: + + PYTHONPATH=$FORTSYM/python python3 docs/derivations/pmodb_gauge/derive.py + +FortSym (https://github.com/lazy-fortran/fortsym) is needed only to rerun the +proofs, not to build or run GPEC. diff --git a/docs/site.rst b/docs/site.rst index 7905ced0e..a62a5e2eb 100644 --- a/docs/site.rst +++ b/docs/site.rst @@ -9,6 +9,7 @@ GPEC documentation website site map contact.rst outputs.rst sign_conventions.rst + pmodb_gauge.rst dcon_readme.rst examples.rst releases.rst diff --git a/gpec/gpdiag.f b/gpec/gpdiag.f index d5b6825ef..70d19c6b9 100644 --- a/gpec/gpdiag.f +++ b/gpec/gpdiag.f @@ -1096,8 +1096,8 @@ SUBROUTINE gpdiag_pmodb(egnum,xspmn) $ mthsurf,curv2_mn) eulbparmns(istep,:)=eulbpar_mn ! p'xsp_mn+bxcdeltamn - lagbparmns(istep,:)=lagbpar_mn ! with parallel components - llagbparmns(istep,:)=llagbpar_mn ! without parallel components + lagbparmns(istep,:)=lagbpar_mn ! xi^zeta=0 gauge (pmodb lagb) + llagbparmns(istep,:)=llagbpar_mn ! xi_perp only divxmns(istep,:)=-divx_mn ! Bdivxprp curvmns(istep,:)=-curv_mn ! Bkxprp divx2mns(istep,:)=divx2_mn ! jb-weighting diff --git a/gpec/gpout.f b/gpec/gpout.f index a02897d0d..0b61d5e0d 100644 --- a/gpec/gpout.f +++ b/gpec/gpout.f @@ -3038,6 +3038,9 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) eulbparfun(istep,itheta)= $ chi1*(bvt_fun(itheta)+sq%f(4)*bvz_fun(itheta)) $ /(rzphi%f(4)*eqfun%f(1)) +c lagb: Lagrangian dB of xi_perp + lambda B with xi^zeta=0 +c in the working coordinates, not the xi_perp-only value; see +c docs/pmodb_gauge.rst. lagbparfun(istep,itheta)= $ eulbparfun(istep,itheta)+ $ xsp_fun(itheta)*eqfun%fx(1)+ From b0e3f3dec804f84ec2449cdef3246ff1580ef7c5 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 20 Oct 2026 00:25:54 +0200 Subject: [PATCH 2/8] DOCS - BUG FIX - Complete native proofs and consumer gauge contracts --- docs/derivations/pmodb_gauge/derive.py | 6 +- .../pmodb_gauge/fortsym_verdicts.json | 166 +++++++++++++++--- docs/derivations/pmodb_gauge/identities.py | 143 ++++++++++++++- docs/pmodb_gauge.rst | 150 ++++++++++++++-- 4 files changed, 410 insertions(+), 55 deletions(-) diff --git a/docs/derivations/pmodb_gauge/derive.py b/docs/derivations/pmodb_gauge/derive.py index 17337b82d..d6f90cc3b 100644 --- a/docs/derivations/pmodb_gauge/derive.py +++ b/docs/derivations/pmodb_gauge/derive.py @@ -27,9 +27,7 @@ def verdict(expr) -> str: - if not hasattr(expr, "is_zero"): - return "PROVED" if expr == 0 else "NONZERO" - value = expr.is_zero + value = fs.sympify(expr).is_zero return {True: "PROVED", False: "NONZERO"}.get(value, "UNKNOWN") @@ -85,7 +83,7 @@ def main() -> int: out = {"engine": "fortsym native zero test via C ABI (fortsym.sympy)", "fortsym_revision": revision(), "checks": records, "failed": failed} (HERE / "fortsym_verdicts.json").write_text(json.dumps(out, indent=2) + "\n") - print("FAILED: " + " ".join(failed) if failed else "all identities proved") + print("FAILED: " + " ".join(failed) if failed else "required identities proved; see explicit open-gap records") return 1 if failed else 0 diff --git a/docs/derivations/pmodb_gauge/fortsym_verdicts.json b/docs/derivations/pmodb_gauge/fortsym_verdicts.json index 5f0f0e2fd..8de71dd18 100644 --- a/docs/derivations/pmodb_gauge/fortsym_verdicts.json +++ b/docs/derivations/pmodb_gauge/fortsym_verdicts.json @@ -1,42 +1,42 @@ { "engine": "fortsym native zero test via C ABI (fortsym.sympy)", - "fortsym_revision": "7e8ef7abdd90c49bb5f75e2c119e7d181656bcc6", + "fortsym_revision": "ae635c80b42c922f5f2028c5cec67859bb8d3929", "checks": [ { "label": "G1.1", "equation": "J dB^psi = GPEC bwp (pointwise)", "verdict": "PROVED", - "seconds": 0.01 + "seconds": 0.001 }, { "label": "G1.2", "equation": "J dB^theta = GPEC bwt (pointwise)", "verdict": "PROVED", - "seconds": 0.004 + "seconds": 0.001 }, { "label": "G1.3", "equation": "J dB^zeta = GPEC bwz (pointwise)", "verdict": "PROVED", - "seconds": 0.007 + "seconds": 0.001 }, { "label": "G1.4", "equation": "(xi x B)_psi = chi1 (q xi^theta - xi^zeta) = xi_s", "verdict": "PROVED", - "seconds": 0.001 + "seconds": 0.0 }, { "label": "G1.5", "equation": "xi_s(xi + lam B) = xi_s(xi)", "verdict": "PROVED", - "seconds": 0.001 + "seconds": 0.0 }, { "label": "G1F.1", "equation": "bwp = chi1 (m - n q) 2 pi i xsp", "verdict": "PROVED", - "seconds": 0.001 + "seconds": 0.0 }, { "label": "G1F.2", @@ -54,103 +54,103 @@ "label": "G2.1", "equation": "(d_t + q d_z)(J xi^t) = d_z(J xi_s)/chi1 - d_psi(J xi^psi) + J div xi", "verdict": "PROVED", - "seconds": 0.001 + "seconds": 0.0 }, { "label": "G2.2", "equation": "(d_t + q d_z)(J xi^z) = -q d_psi(J xi^psi) - d_t(J xi_s)/chi1 + q J div xi", "verdict": "PROVED", - "seconds": 0.005 + "seconds": 0.0 }, { "label": "G2.3", "equation": "q xwt - xwz = J xi_s / chi1", "verdict": "PROVED", - "seconds": 0.001 + "seconds": 0.0 }, { "label": "G3", "equation": "B.dB = chi1 (bvt + q bvz)/J (eulb times |B|)", "verdict": "PROVED", - "seconds": 0.013 + "seconds": 0.023 }, { "label": "G4.1", "equation": "J xi_perp^theta = xmt - (chi1/B)^2 (xvt + q xvz)/J", "verdict": "PROVED", - "seconds": 0.27 + "seconds": 0.006 }, { "label": "G4.2", "equation": "J xi_perp^zeta = xmz - q (chi1/B)^2 (xvt + q xvz)/J", "verdict": "PROVED", - "seconds": 0.119 + "seconds": 0.004 }, { "label": "G4.3", "equation": "B_i xi_perp^i = 0 for xi_perp^i = xi^i - (B_j xi^j / B_k B^k) B^i", "verdict": "PROVED", - "seconds": 0.012 + "seconds": 0.001 }, { "label": "G4.4", "equation": "xi_perp^theta = xms/(chi1 q) + xmz/(J q) - (chi1/(J B))^2 (xvt + q xvz) (curvfun, gpout.f)", "verdict": "PROVED", - "seconds": 0.183 + "seconds": 0.011 }, { "label": "G5.1", "equation": "xi_g^zeta = 0 for lam = -xi_perp^zeta / B^zeta", "verdict": "PROVED", - "seconds": 0.232 + "seconds": 0.013 }, { "label": "G5.2", "equation": "xi_g^theta = xi_s/(chi1 q)", "verdict": "PROVED", - "seconds": 0.757 + "seconds": 0.017 }, { "label": "G5.2b", "equation": "xi_g^psi = xi^psi", "verdict": "PROVED", - "seconds": 0.008 + "seconds": 0.0 }, { "label": "G5.3", "equation": "curl(xi_g x B) = curl(xi_perp x B) (theta component)", "verdict": "PROVED", - "seconds": 0.076 + "seconds": 0.049 }, { "label": "G5.4", "equation": "lagb (gpout.f) = dB_L(xi_g), the xi^zeta = 0 gauge", "verdict": "PROVED", - "seconds": 0.026 + "seconds": 0.032 }, { "label": "G5.5", "equation": "lagb - dB_L(xi_perp) = -(xi_perp^zeta/q) d_theta B", "verdict": "PROVED", - "seconds": 0.017 + "seconds": 0.027 }, { "label": "G5.6", "equation": "llagb (gpdiag.f) = dB_L(xi_perp)", "verdict": "PROVED", - "seconds": 0.038 + "seconds": 0.027 }, { "label": "G5.7", "equation": "dB_L(xi, div-free) = lagb + (xmz/(J q)) d_theta B", "verdict": "PROVED", - "seconds": 0.014 + "seconds": 0.015 }, { "label": "G5.8", "equation": "xi_perp^zeta = -(B_theta xi_s/chi1 + q B_psi xi^psi)/(B_theta + q B_zeta)", "verdict": "PROVED", - "seconds": 0.735 + "seconds": 0.011 }, { "label": "G7.1", @@ -162,7 +162,7 @@ "label": "G7.2", "equation": "J B^2 (xi_perp.grad B + mu0 p' xi^psi/B)/B = sband xsp + tband xms/chi1, xi = e_theta", "verdict": "PROVED", - "seconds": 0.007 + "seconds": 0.001 }, { "label": "G7.3", @@ -184,27 +184,135 @@ }, { "label": "R1", - "equation": "J div xi (reg) = (1 - R)(J R xsp1 + J' xsp)", + "equation": "J D_effective(reg) = (1-R)(J R xsp1 + J' xsp)", "verdict": "PROVED", - "seconds": 0.014 + "seconds": 0.001 + }, + { + "label": "R1a", + "equation": "J div xi_reconstructed = J(1-R**2)xsp1 + (1-R)J' xsp", + "verdict": "PROVED", + "seconds": 0.001 + }, + { + "label": "R4", + "equation": "chi1(q xi_reconstructed^theta-xi_reconstructed^zeta)=R xms", + "verdict": "PROVED", + "seconds": 0.0 }, { "label": "R2", "equation": "J div xi (R = 1) = 0", "verdict": "PROVED", - "seconds": 0.028 + "seconds": 0.002 }, { "label": "R3", "equation": "J div dB (reg) = 2 pi i chi1 (m - n q)(1 - R) xsp1", "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "E4.1", + "equation": "d/de J(X + e xi) = d_i(J xi^i) = J div xi", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "G8.1", + "equation": "J B**2 div xi_perp = xmat Xp + ymat X + zmat A", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "G8.2", + "equation": "idcon zband2/zband3 normalization from B**2=chi1**2*(g22+2qg23+q**2g33)/J**2", + "verdict": "PROVED", + "seconds": 0.001 + }, + { + "label": "Q1.hamada", + "equation": "T_a(xi+lambda B) = T_a(xi)+B.grad(lambda B**a)/B**a", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "Q2.hamada", + "equation": "D w = -T; w_m = -T_m/[2 pi i(m-nq)]", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "Q3.hamada", + "equation": "delta B_target = delta B_X + w d_theta B (J B**a flux function)", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "Q1.boozer", + "equation": "T_a(xi+lambda B) = T_a(xi)+B.grad(lambda B**a)/B**a", + "verdict": "PROVED", "seconds": 0.001 }, + { + "label": "Q2.boozer", + "equation": "D w = -T; w_m = -T_m/[2 pi i(m-nq)]", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "Q3.boozer", + "equation": "delta B_target = delta B_X + w d_theta B (J B**a flux function)", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "Q4", + "equation": "det d(theta_H+iota h,phi_H+h)/d(theta_H,phi_H)=1+D h", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "Q5", + "equation": "Hamada->Boozer equals direct X->Boozer when w_HB=w_XB-w_XH", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "Q6", + "equation": "field-aligned input: delta B_E=0, xi_perp=0 => target delta B=0", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "Q7", + "equation": "xi_perp limit: zero target source implies no gauge correction", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "Q8", + "equation": "delta(J/R**2)/(J/R**2)=div xi-2 xi_R/R", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "K5", + "equation": "Complex CGL pressure work equals [5/3 |D|**2+1/3 |D+3K|**2]/2", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "P1", + "equation": "(Bk - (llagb-eulb))/B = mu0 p' xi^psi/B**2", + "verdict": "PROVED", + "seconds": 0.0 + }, { "label": "K1", "equation": "CGL energy density = (1/2)[(5/3) D^2 + (1/3)(D + 3K)^2]", "verdict": "PROVED", - "seconds": 0.001 + "seconds": 0.0 }, { "label": "K3", diff --git a/docs/derivations/pmodb_gauge/identities.py b/docs/derivations/pmodb_gauge/identities.py index d52b3bcd0..f38bc384c 100644 --- a/docs/derivations/pmodb_gauge/identities.py +++ b/docs/derivations/pmodb_gauge/identities.py @@ -233,7 +233,7 @@ def columns(sp): out.append(("G4.2", "J xi_perp^zeta = xmz - q (chi1/B)^2 (xvt + q xvz)/J", S.at0(J * xperp[2] - (xmz - q * chi1 ** 2 / B2 * (xv[1] + q * xv[2]) / J)))) # B . xi_perp = 0 as a statement about arbitrary vectors: the metric-jet - # form exceeds FortSym's polynomial-cancellation limit (issue #84), so the + # form is larger than necessary, so the # projection is proved on independent covariant B_i, contravariant B^i, xi^i. Bl = [sp.Symbol(f"Bl_{i}") for i in range(3)] Bu = [sp.Symbol(f"Bu_{i}") for i in range(3)] @@ -330,7 +330,7 @@ def pentrc_kappa_matrices(sp): tband = chi1 ** 2 * (Dt / B2) * (g23 + q * g33) / J res = lhs - sband * xi[0] - tband * xs / chi1 # The residual is linear in xi; proving it for each basis displacement keeps - # every zero test within FortSym's 12-variable polynomial limit (issue #84). + # every check independently covers one displacement component. out = [] for k, name in enumerate(("psi", "theta", "zeta")): r = res @@ -386,9 +386,14 @@ def regularisation(sp): xwt = -(J0 * R * X1 + J1 * X + 2 * pi * I * n / chi1 * J0 * Xs) / (2 * pi * I * s) xwz = -(q * J0 * R * X1 + q * J1 * X + 2 * pi * I * m / chi1 * J0 * Xs) / (2 * pi * I * s) xmt, xmz = R * xwt, R * xwz - jdiv = J0 * R * X1 + J1 * X + 2 * pi * I * m * xmt - 2 * pi * I * n * xmz - out = [("R1", "J div xi (reg) = (1 - R)(J R xsp1 + J' xsp)", - jdiv - (1 - R) * (J0 * R * X1 + J1 * X))] + jdiv_effective = J0 * R * X1 + J1 * X + 2 * pi * I * m * xmt - 2 * pi * I * n * xmz + jdiv_actual = J0 * X1 + J1 * X + 2 * pi * I * m * xmt - 2 * pi * I * n * xmz + out = [("R1", "J D_effective(reg) = (1-R)(J R xsp1 + J' xsp)", + jdiv_effective - (1-R)*(J0*R*X1+J1*X)), + ("R1a", "J div xi_reconstructed = J(1-R**2)xsp1 + (1-R)J' xsp", + jdiv_actual-J0*(1-R**2)*X1-(1-R)*J1*X), + ("R4", "chi1(q xi_reconstructed^theta-xi_reconstructed^zeta)=R xms", + chi1*(q*xmt-xmz)/J0-R*Xs)] jdiv_unreg = (J0 * X1 + J1 * X + 2 * pi * I * m * xwt.subs(R, 1) - 2 * pi * I * n * xwz.subs(R, 1)) out.append(("R2", "J div xi (R = 1) = 0", jdiv_unreg)) @@ -444,19 +449,141 @@ def drive_gauge(sp): (change_pmodb + mu * B * lam * dBdl / v).subs(l, 0))] +def divergence_matrices(sp): + """G8: compare projected contravariant divergence with idcon's x/y/z bands. + + U=(g12+q*g13)/(J*B**2), V=(g23+q*g33)/(J*B**2). The perpendicular + projection gives J xi^theta=chi1**2*(-U X+V A), xi^zeta=q xi^theta-A. + U,V and their theta derivatives remain independent jets here. The field + norm supplies the only algebraic metric relation used by the z bands. + """ + J,Jp,B2,c,q,m,n,X,Xp,A,U,Ut,V,Vt=sp.symbols("J Jp B2 c q m n X Xp A U Ut V Vt") + th,ze=2*sp.pi*sp.I*m,-2*sp.pi*sp.I*n + # Differentiate J xi^theta and J xi^zeta with the product rule. + Y=c**2*(-U*X+V*A) + Yt=c**2*(-Ut*X-U*th*X+Vt*A+V*th*A) + Yz=c**2*(-U*ze*X+V*ze*A) + direct=B2*(J*Xp+Jp*X+Yt+q*Yz-J*ze*A) + bands=J*B2*Xp+(Jp*B2-c**2*B2*Ut-2*sp.pi*sp.I*(m-n*q)*c**2*B2*U)*X + bands+=(c**2*B2*Vt+2*sp.pi*sp.I*(c**2*B2*V*(m-n*q)+n*J*B2))*A + out=[("G8.1","J B**2 div xi_perp = xmat Xp + ymat X + zmat A",direct-bands)] + g22,g23,g33,g12,g13=sp.symbols("g22 g23 g33 g12 g13") + norm=c**2*(g22+2*q*g23+q**2*g33)/J**2 + zcode=c**2/J*(m*(g23+q*g33)+n*(g22+q*g23)) + zproj=(c**2*norm*((g23+q*g33)/(J*norm))*(m-n*q)+n*J*norm) + out.append(("G8.2","idcon zband2/zband3 normalization from B**2=chi1**2*(g22+2qg23+q**2g33)/J**2",zcode-zproj)) + return out + + +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 + + + +def jacobian_identity(sp): + """Eq. (4.1): d/de det(d(X + e xi)/du) = d_i(J xi^i) (Jacobi formula).""" + u = sp.symbols("u_1 u_2 u_3") + eps = sp.Symbol("epsilon") + c = [[sp.Symbol(f"c_{a}{i}") for i in range(3)] for a in range(3)] + q = [[[sp.Symbol(f"q_{a}{min(i, j)}{max(i, j)}") for j in range(3)] + for i in range(3)] for a in range(3)] + pos = [sum(c[a][i] * u[i] for i in range(3)) + + sum(q[a][i][j] * u[i] * u[j] for i in range(3) for j in range(3)) / 2 + for a in range(3)] + xi_up = [sp.Symbol(f"s_{i}") + sum(sp.Symbol(f"t_{i}{j}") * u[j] for j in range(3)) + for i in range(3)] + e_cov = [[pos[a].diff(u[i]) for a in range(3)] for i in range(3)] + xi_cart = [sum(xi_up[i] * e_cov[i][a] for i in range(3)) for a in range(3)] + moved = [pos[a] + eps * xi_cart[a] for a in range(3)] + + def det3(m): + return (m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) + - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) + + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0])) + + def at0(e): + for v in u: + e = e.subs(v, 0) + return e + + # J(X + eps xi) needs no further u-derivative: evaluate entries at u = 0 first. + jac_moved = det3([[at0(moved[a].diff(u[i])) for i in range(3)] for a in range(3)]) + jac = det3([[pos[a].diff(u[i]) for i in range(3)] for a in range(3)]) + first = jac_moved.diff(eps).subs(eps, 0) + divergence_density = at0(sum((jac * xi_up[i]).diff(u[i]) for i in range(3))) + res = first - divergence_density + return [("E4.1", "d/de J(X + e xi) = d_i(J xi^i) = J div xi", res.expand())] + def all_identities(sp): return (induction(sp) + fourier_symbols(sp) + incompressible_closure(sp) + columns(sp) + pentrc_kappa_matrices(sp) + cartesian_kappa(sp) + regularisation(sp) - + [c for c in cgl(sp) if not c[0].endswith("neg")] + drive_gauge(sp)) + + jacobian_identity(sp) + divergence_matrices(sp) + gauge_coordinates(sp) + [c for c in cgl(sp) if not c[0].endswith("neg")] + drive_gauge(sp)) def open_gaps(sp): """Identities FortSym does not decide; their content is covered elsewhere. G6 is the coordinate form of C1 + C2 (proved on Cartesian 1-jets) with - det g = J^2 imposed through g_33. The zero test exceeds the native - polynomial-term limit (lazy-fortran/fortsym#84 class); it is recorded, not + det g = J^2 imposed through g_33. On ae635c80 the zero test still returns UNKNOWN for this large rational + residual; the variable cap removed by #84 is no longer the explanation; it is recorded, not counted as proved. """ return force_balance_identity(sp) diff --git a/docs/pmodb_gauge.rst b/docs/pmodb_gauge.rst index 577e954d3..f5eedc0dd 100644 --- a/docs/pmodb_gauge.rst +++ b/docs/pmodb_gauge.rst @@ -5,8 +5,8 @@ Perturbed field strength (pmodb) and gauge ****************************************** This page states what the columns of ``gpec_pmodb_n#.out`` are, how they are -related, and which of them a field-strength consumer needs. Every identity -quoted here is proved with the FortSym native zero test by +related, and which of them a field-strength consumer needs. Verified identities +are proved with the FortSym native zero test by ``docs/derivations/pmodb_gauge/derive.py``; the verdicts are stored in ``docs/derivations/pmodb_gauge/fortsym_verdicts.json`` and the labels (G1, G5.5, ...) refer to that file. @@ -19,6 +19,9 @@ G5.5, ...) refer to that file. Coordinates and displacement ============================ +Assume smooth nested surfaces with B>0 and nonzero J, chi1 and q. +Fourier inversions below additionally exclude nonzero resonant sources. + GPEC works in the DCON chart :math:`(\psi,\theta,\zeta)` with angles in turns, any straight-field-line Jacobian :math:`\mathcal{J}(\psi,\theta)` selected by ``jac_type``, :math:`\mathcal{J}B^\theta=\chi_1`, @@ -92,9 +95,8 @@ Gauge of ``lagb`` (G5.5, G5.8): \xi_\perp^\zeta=-\frac{B_\theta\xi_s/\chi_1+qB_\psi\xi^\psi}{B_\theta+qB_\zeta}. Kinematic identity (C1; any :math:`\boldsymbol\xi_\perp`, no force balance) -and its pmodb form (C2 with G4; the same statement written directly in the -GPEC chart with :math:`\det g=\mathcal J^2`, G6, exceeds the FortSym -polynomial-term limit and is recorded as an open gap, not as proved): +and its pmodb form (C2 with G4; the optional direct contraction in the +GPEC chart with :math:`\det g=\mathcal J^2`, G6, returns UNKNOWN on ae635c80 and is recorded as an open gap, not as proved): .. math:: @@ -112,9 +114,9 @@ Hence ``lagb`` :math:`=-(`\ ``Bdivxprp``\ :math:`+`\ ``Bkxprp``\ :math:`)` is therefore not an identity. The first term is a gauge term of first order in -the perturbation; it is largest where :math:`B_\theta/(qB_\zeta)`, the -Boozer :math:`B_\psi` and :math:`\partial_\theta B` are large, i.e. toward -the edge. The second term vanishes for an exact equilibrium. +the perturbation. Its size depends on :math:`B_\theta/(qB_\zeta)`, the +Boozer :math:`B_\psi` and :math:`\partial_\theta B` and the displacement spectrum; growth toward the edge is a hypothesis +for a particular equilibrium, not a general bound. The second term vanishes for an exact equilibrium. Regularisation (``reg_flag``, Park et al. 2009 Eq. 10). ``gpeq_sol`` uses :math:`R\,\partial_\psi\xi^\psi`, :math:`R=s^2/(s^2+\sigma^2)`, @@ -124,14 +126,25 @@ Regularisation (``reg_flag``, Park et al. 2009 Eq. 10). ``gpeq_sol`` uses .. math:: - \mathcal J\nabla\cdot\boldsymbol\xi=(1-R)(\mathcal J R\,\partial_\psi\xi^\psi+\mathcal J'\xi^\psi),\qquad + \mathcal J D_{effective}=(1-R)(\mathcal J R\,\partial_\psi\xi^\psi+\mathcal J'\xi^\psi),\qquad \mathcal J\nabla\cdot\delta\mathbf B=2\pi i\chi_1 s(1-R)\,\partial_\psi\xi^\psi. -Inside :math:`|m-nq|\lesssim\sigma` the regularised displacement is neither -divergence free nor does it produce a divergence-free field, so -``Bdivxprp`` (which carries the explicit :math:`R\,\partial_\psi\xi^\psi`) -and ``lagb`` (which uses only ``xms``) are not related by the identities -above on the resonant harmonics. +Here the displayed R1 formula is the *effective* divergence diagnostic, +using :math:`R\partial_\psi\xi^\psi` for its radial derivative. The actual +divergence of the reconstructed displacement is instead (R1a) + +.. math:: + + \mathcal J\nabla\cdot\boldsymbol\xi_{rec} + =\mathcal J(1-R^2)\partial_\psi\xi^\psi+(1-R)\mathcal J'\xi^\psi. + +Also :math:`\chi_1(q\xi_{rec}^\theta-\xi_{rec}^\zeta)=R\xi_s` (R4). +These formulas are for the diagonal model only. Near rational surfaces the +regularized field, displacement and effective derivative therefore do not +form one exact ideal-induction pair. This is a limitation of the smoothing +policy, not evidence of resistive physics. The unconditional smoothing in +``gpout_pmodb`` even with ``reg_flag=f`` is a separate verified output bug. + What consumers need @@ -170,3 +183,112 @@ Reproducing the proofs FortSym (https://github.com/lazy-fortran/fortsym) is needed only to rerun the proofs, not to build or run GPEC. + + +Target gauges in GPEC units (verified Q1-Q8) +============================================ + +Let :math:`D=\partial_\theta+q\partial_\zeta`, so +:math:`\mathbf B\cdot\nabla=(\chi_1/\mathcal J)D`. +For a matching equilibrium chart :math:`\mathcal J B^a=H(\psi)`, +:math:`a=2` is Boozer and :math:`a=0` is Hamada. A displacement with +:math:`b_X=\delta B_L(\boldsymbol\xi_X)` and +:math:`d_X=\nabla\cdot\boldsymbol\xi_X` has target-gauge source + +.. math:: + + T_a=d_X+a b_X/B,\qquad + \boldsymbol\xi_{target}=\boldsymbol\xi_X+\lambda\mathbf B,\qquad + T_a(\boldsymbol\xi_{target})=T_a(\boldsymbol\xi_X) + +\frac{\mathbf B\cdot\nabla(\lambda B^a)}{B^a}. + +For nonaxisymmetric modes the target flux-function perturbation is zero. +With :math:`w=\chi_1\lambda/\mathcal J`, this gives + +.. math:: + + Dw=-T_a,\qquad w_m=\frac{i T_{a,m}}{2\pi(m-nq)},\qquad + b_{target}=b_X+w\partial_\theta B. + +The source is formed in real space before Fourier transformation. A +nonzero source at :math:`m=nq` has no periodic solution. A resonant zero +source leaves an arbitrary homogeneous shift, set to zero here. +The factor :math:`2\pi` is from angles in turns and cancels the same +factor in :math:`\partial_\theta B`. Relabeling equilibrium harmonics +with ``jac_out`` alone does not solve this equation. + +For a PEST chart :math:`\mathcal J/R^2` is a flux function; its +linear perturbation is :math:`\nabla\cdot\boldsymbol\xi-2\xi_R/R` +(Q8). PEST output is not automatically either target gauge. Transform the +equilibrium chart first, then solve the target equation. The Hamada scalar +on Boozer equilibrium angles uses :math:`T_2=2 b_H/B`, matching +rmp_torque c3614334, Eq. 4.14. Q5 proves composition of parallel shifts; +Q4 proves the angular-map determinant. The native regression compares the +finite spectral operator directly with that consumer's kernel. + + +Verified defects and numerical attribution +========================================== + +* ``gpout_pmodb`` smooths the radial derivative even with ``reg_flag=f``; + use ``xmp1_mn`` selected by ``gpeq_sol``. No change for ``reg_flag=t``. +* ``gpdiag_pmodb`` uses a stale ``jac`` in its first theta loop. The local + Jacobian must be assigned before projecting out the parallel component. +* PENTRC ``read_pmodb`` overwrites the supplied divergence with + ``-(lagb+Bkxprp)``. Use the supplied positive ``Bdivxprp`` and + ``Bkxprp`` to form ``b_perp=-(Bdivxprp+Bkxprp)`` instead. +* Both CGL torque paths recover ``kappa.xi_perp`` as + ``-(dbob+divx)/2`` although their input convention is + ``dbob=-(divx+kapx)``. Remove the extra half. K1 and K5 prove the + double-adiabatic pressure work for real and complex perturbations. + +Finding (ii) from the TC24 reconstruction is a missing division by B in +the pressure *regressor*, not a GPEC normalization error (P1). For the +retained ITER n=3 output the original fit gives 0.1677; the corrected +pressure term plus the gauge term fits coefficients (1,1,1) with relative +residual about 1e-4. This is numerical corroboration only. + +Finding (i) combines the gauge term, effective-derivative smoothing, +finite spectral/spline error, and the equilibrium force residual. The +identity above separates them; it does not prove a quantitative allocation +of the originally reported 0.5-19 percent to any one cause. The explicit +pressure term in curvature has coefficient one. Edge growth alone does +not diagnose nonideal plasma physics. Finding (iii) is consistent with the +explicit radial derivative in divergence and its smoothing near rational +surfaces. No change to the smoothing model itself is proposed here. + + +Source map and literature +========================== + +Verified code contracts (line numbers at Fortran base 3622c25c): +``gpec/idcon.f:93`` (chi1=2pi psio), ``:617`` (projection geometry), +``:675`` (pressure in Pa), ``:1026`` (curvature and divergence matrices), +``gpec/gpeq.f:114`` (regularization selection), ``:183`` (surface-displacement +closure), ``gpec/gpout.f:3035`` (eulb/lagb), ``:3047`` (divergence and +pressure curvature), ``gpec/gpdiag.f:1003`` (diagnostic projection), +``:1061`` (llagb), ``pentrc/inputs.f90:469`` (ASCII columns), +``:487`` (old gauge mixing), ``:944`` (perpendicular matrix pair), +``pentrc/torque.F90:379`` (old CGL half). + +Julia base 0304fd54c: ``src/KineticForces/KineticForcesStructs.jl:435`` +uses the correct perpendicular dbob convention; +``src/KineticForces/Torque.jl:295`` repeats the CGL half. +``src/PerturbedEquilibrium/FieldReconstruction.jl``'s xi_s is Glasser's +surface displacement, not xi^zeta. Eq. 18 describes its elimination in the +Euler-Lagrange system; Eq. 9 defines it. + +The physical definitions were checked against Park 2009, PoP 16 056115, +Eqs. 7 and 10; Park 2011, PoP 18 110702, Eqs. 4, 10 and 11; +Park and Logan 2017, PoP 24 032505, Eq. 13 and p. 6; Logan and Park 2013, +PoP 20 122507, Eqs. 4 and 5; and Boozer and Nuehrenberg 2006, +PoP 13 102501, Eqs. 42 and 43. The latter's Delta includes twice the +curvature; its Eq. 43 is Eulerian and is not the perpendicular Lagrangian +scalar. These are source contracts, separate from the native algebra proofs. + +FortSym limits: #84, #85 and #86 are present at ae635c80. The direct G6 +contraction remains UNKNOWN (33 symbols, rational metric jets); no probe +is promoted to proof. #91 (coefficient extraction after zero substitution) +and #92 (radical power reduction) are not exercised by the accepted +polynomial/rational gauge proofs. C1/C2 supply the invariant kinematic +identity independently of G6. From 97716d7e07867738963a7d0239bebd4b3b8bb8b6 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 20 Oct 2026 00:32:19 +0200 Subject: [PATCH 3/8] DOCS - BUG FIX - Distinguish regularized llagb from a perpendicular projection --- docs/derivations/pmodb_gauge/check_output.py | 120 ++++++++++++++++++ .../pmodb_gauge/fortsym_verdicts.json | 28 ++-- docs/derivations/pmodb_gauge/identities.py | 3 + docs/pmodb_gauge.rst | 15 ++- 4 files changed, 154 insertions(+), 12 deletions(-) create mode 100644 docs/derivations/pmodb_gauge/check_output.py diff --git a/docs/derivations/pmodb_gauge/check_output.py b/docs/derivations/pmodb_gauge/check_output.py new file mode 100644 index 000000000..c4c0f7843 --- /dev/null +++ b/docs/derivations/pmodb_gauge/check_output.py @@ -0,0 +1,120 @@ +"""Numeric oracle for the pmodb gauge identities on the TC24 n=3 GPEC run. + +Usage: PYTHONPATH=$RMP_TORQUE/src python3 check_output.py RUN_DIR + +Requires numpy, scipy, and the gpec_ascii readers from the referenced +rmp_torque Boozer-gauge-total branch. Numerical evidence only; see derive.py +for the native FortSym proofs. +""" +import sys +from pathlib import Path + +import numpy as np +from scipy.interpolate import CubicSpline +from rmp_torque.gpec_ascii import (read_pmodb, read_xclebsch, read_equilibrium_grid, + reconstruct_advective_spectrum) + +R = Path(sys.argv[1]) +MU0 = 4e-7 * np.pi + + +def table(path, ncol): + rows = [] + for line in open(path): + f = line.split() + if len(f) != ncol: + continue + try: + rows.append([float(x) for x in f]) + except ValueError: + pass + return np.array(rows) + + +d = table(R / "gpec_diagnostics_pmodb_n3.out", 16) +psi = np.unique(d[:, 0]); ms = np.unique(d[:, 1]).astype(int) +shape = (psi.size, ms.size) +c = lambda k: (d[:, k] + 1j * d[:, k + 1]).reshape(shape) +lagb, llagb, lll, bdiv, bk = c(2), c(4), c(6), c(8), c(12) +sp = read_pmodb(R / "gpec_pmodb_n3.out") +xc = read_xclebsch(R / "gpec_xclebsch_n3.out") +eq = read_equilibrium_grid(R / "idcon_equil.out") +q = np.interp(psi, eq.psi_poloidal, eq.q) +nres = 3 * q +far = np.abs(nres - np.round(nres)) > 0.15 +rel = lambda a, b: np.linalg.norm(a - b, axis=1) / np.linalg.norm(b, axis=1) +bands = [(0.0, 0.8), (0.8, 0.9), (0.9, 0.98)] + + +def report(name, a, b): + r = rel(a, b) + parts = [] + for lo, hi in bands: + sel = (psi >= lo) & (psi < hi) + parts.append(f"{lo:.1f}-{hi:.2f}: med {np.median(r[sel]):.4f} " + f"(far {np.median(r[sel & far]):.4f})") + print(f"{name:34s}", " | ".join(parts)) + + +print("diag grid", shape, "pmodb grid", sp.lagrangian.shape) +assert np.allclose(sp.psi_poloidal, psi) +report("pmodb lagb vs -(Bdiv+Bk) [cols]", sp.lagrangian, + -(sp.b_div_xi_perp + sp.b_kappa_xi_perp)) +report("diag lagb vs pmodb lagb", lagb, sp.lagrangian) +report("diag lagb vs lllagb=-(Bdiv+Bk)", lagb, lll) +report("diag llagb vs lllagb", llagb, lll) +report("gauge term (lagb-llagb) / llagb", lagb - llagb + llagb, llagb) + +# Finding (ii): the pressure-coefficient fit. +fl = [] +sec = False +for line in open(R / "idcon_equil.out"): + if "Flux functions:" in line: + sec = True + continue + f = line.split() + if sec and len(f) == 6: + try: + fl.append([float(x) for x in f]) + except ValueError: + pass +fl = np.array(fl) +dp = CubicSpline(fl[:, 0], fl[:, 2]).derivative()(xc.psi_poloidal) +th = eq.theta_radians +ph = np.exp(1j * np.outer(xc.modes, th)) +bs = CubicSpline(eq.psi_poloidal, eq.b_magnitude, axis=0)(xc.psi_poloidal) + + +def spec(x): + r = np.fft.fft((x @ ph) / bs, axis=1) / len(th) + return np.stack([r[:, int(m) % len(th)] for m in xc.modes], axis=1) + + +adv = reconstruct_advective_spectrum(xc, eq) # xi_g . grad B (xi^zeta = 0) +A1 = spec(adv) +A2_wrong = spec(MU0 * xc.xi_psi * dp[:, None]) # rmp_torque #56 fit: missing 1/B +_r = np.fft.fft(MU0 * dp[:, None] * (xc.xi_psi @ ph) / bs**2, axis=1) / len(th) +A2 = np.stack([_r[:, int(m) % len(th)] for m in xc.modes], axis=1) # (mu0 p' xi^psi/B)/B +A3 = spec(llagb - lagb) # (xi_perp^zeta/q) d_theta B +y = spec(sp.b_kappa_xi_perp) +sel = far[:, None] & (np.abs(xc.modes) <= 12)[None, :] +for name, cols in (("adv, p/B [#56]", (A1, A2_wrong)), ("adv, p", (A1, A2)), ("adv, p, gauge", (A1, A2, A3))): + M = np.stack([a[sel] for a in cols], axis=1) + coef, *_ = np.linalg.lstsq(M, y[sel], rcond=None) + res = np.linalg.norm(M @ coef - y[sel]) / np.linalg.norm(y[sel]) + print(f"fit Bkxprp ~ {name:14s}: coef", np.round(coef, 4), f"resid {res:.4f}") +M = A1 + A2 + A3 +print("fixed coefficients (1,1,1): resid %.4f" % ( + np.linalg.norm(M[sel] - y[sel]) / np.linalg.norm(y[sel]))) + +print("=== consistency of the regressors (far surfaces, |m|<=12)") +eul = sp.eulerian +def r2(a, b): + return np.linalg.norm(a[sel] - b[sel]) / np.linalg.norm(b[sel]) +print("adv (xclebsch) vs lagb-eulb %.4f" % r2(A1, spec(sp.lagrangian - eul))) +print("p-term vs Bk-(llagb-eulb) %.4f" % r2(A2, spec(sp.b_kappa_xi_perp - llagb + eul))) +print("Bk vs (llagb-eulb)+p-term %.4f" % r2(y, spec(llagb - eul) + A2)) + +T = spec(sp.b_kappa_xi_perp - llagb + eul) +scale = np.vdot(A2[sel], T[sel]) / np.vdot(A2[sel], A2[sel]) +print("pressure coefficient against corrected regressor:", scale) diff --git a/docs/derivations/pmodb_gauge/fortsym_verdicts.json b/docs/derivations/pmodb_gauge/fortsym_verdicts.json index 8de71dd18..1090e2204 100644 --- a/docs/derivations/pmodb_gauge/fortsym_verdicts.json +++ b/docs/derivations/pmodb_gauge/fortsym_verdicts.json @@ -72,13 +72,13 @@ "label": "G3", "equation": "B.dB = chi1 (bvt + q bvz)/J (eulb times |B|)", "verdict": "PROVED", - "seconds": 0.023 + "seconds": 0.019 }, { "label": "G4.1", "equation": "J xi_perp^theta = xmt - (chi1/B)^2 (xvt + q xvz)/J", "verdict": "PROVED", - "seconds": 0.006 + "seconds": 0.005 }, { "label": "G4.2", @@ -96,19 +96,19 @@ "label": "G4.4", "equation": "xi_perp^theta = xms/(chi1 q) + xmz/(J q) - (chi1/(J B))^2 (xvt + q xvz) (curvfun, gpout.f)", "verdict": "PROVED", - "seconds": 0.011 + "seconds": 0.008 }, { "label": "G5.1", "equation": "xi_g^zeta = 0 for lam = -xi_perp^zeta / B^zeta", "verdict": "PROVED", - "seconds": 0.013 + "seconds": 0.008 }, { "label": "G5.2", "equation": "xi_g^theta = xi_s/(chi1 q)", "verdict": "PROVED", - "seconds": 0.017 + "seconds": 0.013 }, { "label": "G5.2b", @@ -120,37 +120,37 @@ "label": "G5.3", "equation": "curl(xi_g x B) = curl(xi_perp x B) (theta component)", "verdict": "PROVED", - "seconds": 0.049 + "seconds": 0.029 }, { "label": "G5.4", "equation": "lagb (gpout.f) = dB_L(xi_g), the xi^zeta = 0 gauge", "verdict": "PROVED", - "seconds": 0.032 + "seconds": 0.024 }, { "label": "G5.5", "equation": "lagb - dB_L(xi_perp) = -(xi_perp^zeta/q) d_theta B", "verdict": "PROVED", - "seconds": 0.027 + "seconds": 0.019 }, { "label": "G5.6", "equation": "llagb (gpdiag.f) = dB_L(xi_perp)", "verdict": "PROVED", - "seconds": 0.027 + "seconds": 0.018 }, { "label": "G5.7", "equation": "dB_L(xi, div-free) = lagb + (xmz/(J q)) d_theta B", "verdict": "PROVED", - "seconds": 0.015 + "seconds": 0.013 }, { "label": "G5.8", "equation": "xi_perp^zeta = -(B_theta xi_s/chi1 + q B_psi xi^psi)/(B_theta + q B_zeta)", "verdict": "PROVED", - "seconds": 0.011 + "seconds": 0.01 }, { "label": "G7.1", @@ -200,6 +200,12 @@ "verdict": "PROVED", "seconds": 0.0 }, + { + "label": "R5", + "equation": "llagb_code - (eulb+xi_rec_perp.grad B) = (1-R) xms Bt/(chi1 q)", + "verdict": "PROVED", + "seconds": 0.001 + }, { "label": "R2", "equation": "J div xi (R = 1) = 0", diff --git a/docs/derivations/pmodb_gauge/identities.py b/docs/derivations/pmodb_gauge/identities.py index f38bc384c..05df76c13 100644 --- a/docs/derivations/pmodb_gauge/identities.py +++ b/docs/derivations/pmodb_gauge/identities.py @@ -394,6 +394,9 @@ def regularisation(sp): jdiv_actual-J0*(1-R**2)*X1-(1-R)*J1*X), ("R4", "chi1(q xi_reconstructed^theta-xi_reconstructed^zeta)=R xms", chi1*(q*xmt-xmz)/J0-R*Xs)] + Bt = sp.Symbol("Bt") + out.append(("R5", "llagb_code - (eulb+xi_rec_perp.grad B) = (1-R) xms Bt/(chi1 q)", + (Xs/(chi1*q)+xmz/(J0*q)-xmt/J0)*Bt-(1-R)*Xs*Bt/(chi1*q))) jdiv_unreg = (J0 * X1 + J1 * X + 2 * pi * I * m * xwt.subs(R, 1) - 2 * pi * I * n * xwz.subs(R, 1)) out.append(("R2", "J div xi (R = 1) = 0", jdiv_unreg)) diff --git a/docs/pmodb_gauge.rst b/docs/pmodb_gauge.rst index f5eedc0dd..7c25f69e6 100644 --- a/docs/pmodb_gauge.rst +++ b/docs/pmodb_gauge.rst @@ -48,6 +48,10 @@ For every ``jac_type`` (G1, G1F) the Jacobian-weighted perturbed field of The pmodb columns ================= +The displacement identifications in this section assume unregularized +consistent components. The literal lagb expression is valid for both modes; +the regularized closure defect is given by R4-R5 below. + ``eulb`` :math:`\delta B_E=\hat{\mathbf b}\cdot\delta\mathbf B`, the Eulerian field strength (G3). @@ -139,6 +143,14 @@ divergence of the reconstructed displacement is instead (R1a) =\mathcal J(1-R^2)\partial_\psi\xi^\psi+(1-R)\mathcal J'\xi^\psi. Also :math:`\chi_1(q\xi_{rec}^\theta-\xi_{rec}^\zeta)=R\xi_s` (R4). +The diagnostic llagb also differs from the projection of the reconstructed +contravariant displacement by (R5) + +.. math:: + + \mathrm{llagb}_{code}-(\delta B_E+\boldsymbol\xi_{rec,\perp}\cdot\nabla B) + =(1-R)\frac{\xi_s}{\chi_1 q}\partial_\theta B. + These formulas are for the diagonal model only. Near rational surfaces the regularized field, displacement and effective derivative therefore do not form one exact ideal-induction pair. This is a limitation of the smoothing @@ -246,7 +258,8 @@ Finding (ii) from the TC24 reconstruction is a missing division by B in the pressure *regressor*, not a GPEC normalization error (P1). For the retained ITER n=3 output the original fit gives 0.1677; the corrected pressure term plus the gauge term fits coefficients (1,1,1) with relative -residual about 1e-4. This is numerical corroboration only. +residual about 1e-4. This is numerical corroboration only; ``docs/derivations/pmodb_gauge/check_output.py`` +reproduces the fit from the pmodb, xclebsch, diagnostic and IDCON ASCII tables. Finding (i) combines the gauge term, effective-derivative smoothing, finite spectral/spline error, and the equilibrium force residual. The From e032054c8e5b713d1a6a843e0d01d4a4aa13edc5 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 20 Oct 2026 00:37:43 +0200 Subject: [PATCH 4/8] DOCS - MINOR - Pin source references to the audited upstream revision --- docs/derivations/pmodb_gauge/identities.py | 4 ++-- docs/pmodb_gauge.rst | 16 ++++++++-------- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/docs/derivations/pmodb_gauge/identities.py b/docs/derivations/pmodb_gauge/identities.py index 05df76c13..1dadb9506 100644 --- a/docs/derivations/pmodb_gauge/identities.py +++ b/docs/derivations/pmodb_gauge/identities.py @@ -155,7 +155,7 @@ def xi_s(self, v): def induction(sp): - """G1: GPEC's J dB^i (gpeq.f:107-109) is J curl(xi x B)^i for any Jacobian.""" + """G1: GPEC's J dB^i (gpeq.f:114-116) is J curl(xi x B)^i for any Jacobian.""" S = Setup(sp) xi, B = S.xi, S.Bup JdB = [S.J * c for c in S.curl_up(S.cross_cov(xi, B))] @@ -586,7 +586,7 @@ def open_gaps(sp): G6 is the coordinate form of C1 + C2 (proved on Cartesian 1-jets) with det g = J^2 imposed through g_33. On ae635c80 the zero test still returns UNKNOWN for this large rational - residual; the variable cap removed by #84 is no longer the explanation; it is recorded, not + residual; the removed variable cap is no longer the explanation; it is recorded, not counted as proved. """ return force_balance_identity(sp) diff --git a/docs/pmodb_gauge.rst b/docs/pmodb_gauge.rst index 7c25f69e6..bfa03b5ac 100644 --- a/docs/pmodb_gauge.rst +++ b/docs/pmodb_gauge.rst @@ -275,16 +275,16 @@ Source map and literature ========================== Verified code contracts (line numbers at Fortran base 3622c25c): -``gpec/idcon.f:93`` (chi1=2pi psio), ``:617`` (projection geometry), -``:675`` (pressure in Pa), ``:1026`` (curvature and divergence matrices), -``gpec/gpeq.f:114`` (regularization selection), ``:183`` (surface-displacement -closure), ``gpec/gpout.f:3035`` (eulb/lagb), ``:3047`` (divergence and -pressure curvature), ``gpec/gpdiag.f:1003`` (diagnostic projection), -``:1061`` (llagb), ``pentrc/inputs.f90:469`` (ASCII columns), -``:487`` (old gauge mixing), ``:944`` (perpendicular matrix pair), +``gpec/idcon.f:93`` (chi1=2pi psio), ``:616`` (projection geometry), +``:677`` (pressure in Pa), ``:1026`` (curvature and divergence matrices), +``gpec/gpeq.f:131`` (regularization selection), ``:208`` (surface-displacement +closure), ``gpec/gpout.f:3041`` (eulb/lagb), ``:3050`` (divergence and +pressure curvature), ``gpec/gpdiag.f:1012`` (diagnostic projection), +``:1061`` (llagb), ``pentrc/inputs.f90:376`` (ASCII columns), +``:396`` (old gauge mixing), ``:856`` (perpendicular matrix pair), ``pentrc/torque.F90:379`` (old CGL half). -Julia base 0304fd54c: ``src/KineticForces/KineticForcesStructs.jl:435`` +Julia base 0304fd54c: ``src/KineticForces/KineticForcesStructs.jl:437`` uses the correct perpendicular dbob convention; ``src/KineticForces/Torque.jl:295`` repeats the CGL half. ``src/PerturbedEquilibrium/FieldReconstruction.jl``'s xi_s is Glasser's From cd0dfc349c8db19ddab6973bc27f6f61f573d4ff Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 20 Oct 2026 00:45:32 +0200 Subject: [PATCH 5/8] DOCS - NEW FEATURE - Prove the general-J pmodb error budget --- .../pmodb_gauge/fortsym_verdicts.json | 56 +++++++++++-------- docs/derivations/pmodb_gauge/identities.py | 34 ++++++++++- docs/pmodb_gauge.rst | 25 +++++++++ 3 files changed, 92 insertions(+), 23 deletions(-) diff --git a/docs/derivations/pmodb_gauge/fortsym_verdicts.json b/docs/derivations/pmodb_gauge/fortsym_verdicts.json index 1090e2204..b271512dc 100644 --- a/docs/derivations/pmodb_gauge/fortsym_verdicts.json +++ b/docs/derivations/pmodb_gauge/fortsym_verdicts.json @@ -6,19 +6,19 @@ "label": "G1.1", "equation": "J dB^psi = GPEC bwp (pointwise)", "verdict": "PROVED", - "seconds": 0.001 + "seconds": 0.0 }, { "label": "G1.2", "equation": "J dB^theta = GPEC bwt (pointwise)", "verdict": "PROVED", - "seconds": 0.001 + "seconds": 0.0 }, { "label": "G1.3", "equation": "J dB^zeta = GPEC bwz (pointwise)", "verdict": "PROVED", - "seconds": 0.001 + "seconds": 0.0 }, { "label": "G1.4", @@ -72,43 +72,43 @@ "label": "G3", "equation": "B.dB = chi1 (bvt + q bvz)/J (eulb times |B|)", "verdict": "PROVED", - "seconds": 0.019 + "seconds": 0.009 }, { "label": "G4.1", "equation": "J xi_perp^theta = xmt - (chi1/B)^2 (xvt + q xvz)/J", "verdict": "PROVED", - "seconds": 0.005 + "seconds": 0.003 }, { "label": "G4.2", "equation": "J xi_perp^zeta = xmz - q (chi1/B)^2 (xvt + q xvz)/J", "verdict": "PROVED", - "seconds": 0.004 + "seconds": 0.002 }, { "label": "G4.3", "equation": "B_i xi_perp^i = 0 for xi_perp^i = xi^i - (B_j xi^j / B_k B^k) B^i", "verdict": "PROVED", - "seconds": 0.001 + "seconds": 0.0 }, { "label": "G4.4", "equation": "xi_perp^theta = xms/(chi1 q) + xmz/(J q) - (chi1/(J B))^2 (xvt + q xvz) (curvfun, gpout.f)", "verdict": "PROVED", - "seconds": 0.008 + "seconds": 0.004 }, { "label": "G5.1", "equation": "xi_g^zeta = 0 for lam = -xi_perp^zeta / B^zeta", "verdict": "PROVED", - "seconds": 0.008 + "seconds": 0.006 }, { "label": "G5.2", "equation": "xi_g^theta = xi_s/(chi1 q)", "verdict": "PROVED", - "seconds": 0.013 + "seconds": 0.007 }, { "label": "G5.2b", @@ -120,37 +120,37 @@ "label": "G5.3", "equation": "curl(xi_g x B) = curl(xi_perp x B) (theta component)", "verdict": "PROVED", - "seconds": 0.029 + "seconds": 0.017 }, { "label": "G5.4", "equation": "lagb (gpout.f) = dB_L(xi_g), the xi^zeta = 0 gauge", "verdict": "PROVED", - "seconds": 0.024 + "seconds": 0.012 }, { "label": "G5.5", "equation": "lagb - dB_L(xi_perp) = -(xi_perp^zeta/q) d_theta B", "verdict": "PROVED", - "seconds": 0.019 + "seconds": 0.008 }, { "label": "G5.6", "equation": "llagb (gpdiag.f) = dB_L(xi_perp)", "verdict": "PROVED", - "seconds": 0.018 + "seconds": 0.008 }, { "label": "G5.7", "equation": "dB_L(xi, div-free) = lagb + (xmz/(J q)) d_theta B", "verdict": "PROVED", - "seconds": 0.013 + "seconds": 0.009 }, { "label": "G5.8", "equation": "xi_perp^zeta = -(B_theta xi_s/chi1 + q B_psi xi^psi)/(B_theta + q B_zeta)", "verdict": "PROVED", - "seconds": 0.01 + "seconds": 0.005 }, { "label": "G7.1", @@ -162,7 +162,7 @@ "label": "G7.2", "equation": "J B^2 (xi_perp.grad B + mu0 p' xi^psi/B)/B = sband xsp + tband xms/chi1, xi = e_theta", "verdict": "PROVED", - "seconds": 0.001 + "seconds": 0.0 }, { "label": "G7.3", @@ -186,7 +186,7 @@ "label": "R1", "equation": "J D_effective(reg) = (1-R)(J R xsp1 + J' xsp)", "verdict": "PROVED", - "seconds": 0.001 + "seconds": 0.0 }, { "label": "R1a", @@ -204,13 +204,13 @@ "label": "R5", "equation": "llagb_code - (eulb+xi_rec_perp.grad B) = (1-R) xms Bt/(chi1 q)", "verdict": "PROVED", - "seconds": 0.001 + "seconds": 0.0 }, { "label": "R2", "equation": "J div xi (R = 1) = 0", "verdict": "PROVED", - "seconds": 0.002 + "seconds": 0.001 }, { "label": "R3", @@ -218,6 +218,18 @@ "verdict": "PROVED", "seconds": 0.0 }, + { + "label": "R6", + "equation": "H_parallel = B(X1-U)+(B_theta Delta_zeta-B_zeta Delta_theta)/(J B)", + "verdict": "PROVED", + "seconds": 0.0 + }, + { + "label": "R7", + "equation": "lagb+Bdivxprp+Bkxprp = -xi_rec_perp^zeta Bt/q + 2 Delta Bt/(chi1 q) + curl_delta + B(V-U) - xi_rec_perp.F/B", + "verdict": "PROVED", + "seconds": 0.001 + }, { "label": "E4.1", "equation": "d/de J(X + e xi) = d_i(J xi^i) = J div xi", @@ -234,7 +246,7 @@ "label": "G8.2", "equation": "idcon zband2/zband3 normalization from B**2=chi1**2*(g22+2qg23+q**2g33)/J**2", "verdict": "PROVED", - "seconds": 0.001 + "seconds": 0.0 }, { "label": "Q1.hamada", @@ -258,7 +270,7 @@ "label": "Q1.boozer", "equation": "T_a(xi+lambda B) = T_a(xi)+B.grad(lambda B**a)/B**a", "verdict": "PROVED", - "seconds": 0.001 + "seconds": 0.0 }, { "label": "Q2.boozer", diff --git a/docs/derivations/pmodb_gauge/identities.py b/docs/derivations/pmodb_gauge/identities.py index 1dadb9506..e268dee93 100644 --- a/docs/derivations/pmodb_gauge/identities.py +++ b/docs/derivations/pmodb_gauge/identities.py @@ -410,6 +410,38 @@ def regularisation(sp): return out +def regularized_column_budget(sp): + """General-J continuum budget, allowing independent field/output derivatives. + + U is the derivative selected by gpeq_sol; V is the derivative used in + gpout divergence. Delta=xi_s-chi1(q xi_rec^theta-xi_rec^zeta). + F denotes xi_rec_perp.((curl B)xB-mu0 grad p)/B, in Tesla. + Metric/field-norm consistency and C1/C2 are assumed; finite truncation + errors are not folded into the physical force residual. + """ + c,J,B,q,Bt,Bt_cov,Delta,Dt,Dz,X1,U,V,D,K,F,Z = sp.symbols( + "chi J B q Bt Btheta Delta Delta_theta Delta_zeta X1 U V D K F Z") + Bz_cov=(J*B**2/c-Bt_cov)/q + # J-weighted induction mismatch: supplied field minus curl(xi_rec x B). + Ht=c*(X1-U)+Dz + Hz=c*q*(X1-U)-Dt + Hpar=(Bt_cov*Ht+Bz_cov*Hz)/(J*B) + curl_delta=(Bt_cov*Dz-Bz_cov*Dt)/(J*B) + out=[("R6", "H_parallel = B(X1-U)+(B_theta Delta_zeta-B_zeta Delta_theta)/(J B)", + Hpar-B*(X1-U)-curl_delta)] + g=Delta*Bt/(c*q) + # C1 for the reconstructed perpendicular displacement, plus the supplied H. + LLrec=-B*(D+K)+Hpar + LLcode=LLrec+g + Lagcode=LLcode-Z*Bt/q + BDcode=B*(D+V-X1) + BKcode=B*K+g-F + budget=-Z*Bt/q+2*g+curl_delta+B*(V-U)-F + out.append(("R7", "lagb+Bdivxprp+Bkxprp = -xi_rec_perp^zeta Bt/q + 2 Delta Bt/(chi1 q) + curl_delta + B(V-U) - xi_rec_perp.F/B", + Lagcode+BDcode+BKcode-budget)) + return out + + def cgl(sp): """K1-K2: the PENTRC/GPEC.jl fcgl integrand against the CGL closure. @@ -577,7 +609,7 @@ def at0(e): def all_identities(sp): return (induction(sp) + fourier_symbols(sp) + incompressible_closure(sp) + columns(sp) + pentrc_kappa_matrices(sp) - + cartesian_kappa(sp) + regularisation(sp) + + cartesian_kappa(sp) + regularisation(sp) + regularized_column_budget(sp) + jacobian_identity(sp) + divergence_matrices(sp) + gauge_coordinates(sp) + [c for c in cgl(sp) if not c[0].endswith("neg")] + drive_gauge(sp)) diff --git a/docs/pmodb_gauge.rst b/docs/pmodb_gauge.rst index bfa03b5ac..07aef9af0 100644 --- a/docs/pmodb_gauge.rst +++ b/docs/pmodb_gauge.rst @@ -158,6 +158,31 @@ policy, not evidence of resistive physics. The unconditional smoothing in ``gpout_pmodb`` even with ``reg_flag=f`` is a separate verified output bug. +A general-J continuum error budget (verified R6-R7) is available without +assuming a diagonal regularizer. Let U be the derivative selected by +``gpeq_sol``, V the derivative used by pmodb divergence, and +:math:`\Delta_s=\xi_s-\chi_1(q\xi_{rec}^\theta-\xi_{rec}^\zeta)`. +Writing :math:`B_\theta,B_\zeta` for covariant components gives + +.. math:: + + \mathrm{lagb}+\mathrm{Bdivxprp}+\mathrm{Bkxprp} + =-\frac{\xi_{rec,\perp}^\zeta}{q}\partial_\theta B + +\frac{2\Delta_s}{\chi_1q}\partial_\theta B + +\frac{B_\theta\partial_\zeta\Delta_s-B_\zeta\partial_\theta\Delta_s} + {\mathcal J B} + +B(V-U)-\frac{\boldsymbol\xi_{rec,\perp}\cdot\mathbf F}{B}. + +This uses consistent continuum geometry and the field norm. Finite Fourier +truncation and spline errors are additional numerical errors, not part of F. +The derivative term cancels when V=U, even if both are regularized; in the +old reg_flag=f output V differs from U. With regularization enabled, +Delta_s can be nonzero and leaves the two explicit closure terms. Thus a +large derivative in Bdivxprp alone does not measure its contribution to the +net mismatch. This budget separates the candidate causes of finding (i), +without claiming a numerical allocation of the original edge residual. + + What consumers need =================== From 1abfd57fde112b543ee28e605c2f46c2361d5467 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 24 Oct 2026 00:59:27 -0400 Subject: [PATCH 6/8] DOCS - CLEANUP - Bring the pmodb gauge page in line with develop #340-#343 merged after this page was written. Describe and prove the current xi_perp construction from xi^psi and xi_s (G4.1, G4.2, G4.4, G5.6), state that both PENTRC input paths now use the xi_perp pair, and reduce the regularisation section to what the output columns use (R3, R8). Replace the list of fixed defects with a short history, and drop material that depends on files outside this repository (check_output.py, the TC24 attribution, rmp_torque and FortSym issue references). derive.py falls back to SymPy when FortSym is absent; the verdicts were regenerated with SymPy 1.9 (57 proved, G6 open, all negative controls nonzero at exact rational witnesses). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VMDGyfQKRZQMdk3NJXsKQE --- docs/derivations/pmodb_gauge/check_output.py | 120 -------- docs/derivations/pmodb_gauge/derive.py | 56 +++- .../pmodb_gauge/fortsym_verdicts.json | 160 +++++------ docs/derivations/pmodb_gauge/identities.py | 70 +++-- docs/pmodb_gauge.rst | 266 +++++++----------- 5 files changed, 242 insertions(+), 430 deletions(-) delete mode 100644 docs/derivations/pmodb_gauge/check_output.py diff --git a/docs/derivations/pmodb_gauge/check_output.py b/docs/derivations/pmodb_gauge/check_output.py deleted file mode 100644 index c4c0f7843..000000000 --- a/docs/derivations/pmodb_gauge/check_output.py +++ /dev/null @@ -1,120 +0,0 @@ -"""Numeric oracle for the pmodb gauge identities on the TC24 n=3 GPEC run. - -Usage: PYTHONPATH=$RMP_TORQUE/src python3 check_output.py RUN_DIR - -Requires numpy, scipy, and the gpec_ascii readers from the referenced -rmp_torque Boozer-gauge-total branch. Numerical evidence only; see derive.py -for the native FortSym proofs. -""" -import sys -from pathlib import Path - -import numpy as np -from scipy.interpolate import CubicSpline -from rmp_torque.gpec_ascii import (read_pmodb, read_xclebsch, read_equilibrium_grid, - reconstruct_advective_spectrum) - -R = Path(sys.argv[1]) -MU0 = 4e-7 * np.pi - - -def table(path, ncol): - rows = [] - for line in open(path): - f = line.split() - if len(f) != ncol: - continue - try: - rows.append([float(x) for x in f]) - except ValueError: - pass - return np.array(rows) - - -d = table(R / "gpec_diagnostics_pmodb_n3.out", 16) -psi = np.unique(d[:, 0]); ms = np.unique(d[:, 1]).astype(int) -shape = (psi.size, ms.size) -c = lambda k: (d[:, k] + 1j * d[:, k + 1]).reshape(shape) -lagb, llagb, lll, bdiv, bk = c(2), c(4), c(6), c(8), c(12) -sp = read_pmodb(R / "gpec_pmodb_n3.out") -xc = read_xclebsch(R / "gpec_xclebsch_n3.out") -eq = read_equilibrium_grid(R / "idcon_equil.out") -q = np.interp(psi, eq.psi_poloidal, eq.q) -nres = 3 * q -far = np.abs(nres - np.round(nres)) > 0.15 -rel = lambda a, b: np.linalg.norm(a - b, axis=1) / np.linalg.norm(b, axis=1) -bands = [(0.0, 0.8), (0.8, 0.9), (0.9, 0.98)] - - -def report(name, a, b): - r = rel(a, b) - parts = [] - for lo, hi in bands: - sel = (psi >= lo) & (psi < hi) - parts.append(f"{lo:.1f}-{hi:.2f}: med {np.median(r[sel]):.4f} " - f"(far {np.median(r[sel & far]):.4f})") - print(f"{name:34s}", " | ".join(parts)) - - -print("diag grid", shape, "pmodb grid", sp.lagrangian.shape) -assert np.allclose(sp.psi_poloidal, psi) -report("pmodb lagb vs -(Bdiv+Bk) [cols]", sp.lagrangian, - -(sp.b_div_xi_perp + sp.b_kappa_xi_perp)) -report("diag lagb vs pmodb lagb", lagb, sp.lagrangian) -report("diag lagb vs lllagb=-(Bdiv+Bk)", lagb, lll) -report("diag llagb vs lllagb", llagb, lll) -report("gauge term (lagb-llagb) / llagb", lagb - llagb + llagb, llagb) - -# Finding (ii): the pressure-coefficient fit. -fl = [] -sec = False -for line in open(R / "idcon_equil.out"): - if "Flux functions:" in line: - sec = True - continue - f = line.split() - if sec and len(f) == 6: - try: - fl.append([float(x) for x in f]) - except ValueError: - pass -fl = np.array(fl) -dp = CubicSpline(fl[:, 0], fl[:, 2]).derivative()(xc.psi_poloidal) -th = eq.theta_radians -ph = np.exp(1j * np.outer(xc.modes, th)) -bs = CubicSpline(eq.psi_poloidal, eq.b_magnitude, axis=0)(xc.psi_poloidal) - - -def spec(x): - r = np.fft.fft((x @ ph) / bs, axis=1) / len(th) - return np.stack([r[:, int(m) % len(th)] for m in xc.modes], axis=1) - - -adv = reconstruct_advective_spectrum(xc, eq) # xi_g . grad B (xi^zeta = 0) -A1 = spec(adv) -A2_wrong = spec(MU0 * xc.xi_psi * dp[:, None]) # rmp_torque #56 fit: missing 1/B -_r = np.fft.fft(MU0 * dp[:, None] * (xc.xi_psi @ ph) / bs**2, axis=1) / len(th) -A2 = np.stack([_r[:, int(m) % len(th)] for m in xc.modes], axis=1) # (mu0 p' xi^psi/B)/B -A3 = spec(llagb - lagb) # (xi_perp^zeta/q) d_theta B -y = spec(sp.b_kappa_xi_perp) -sel = far[:, None] & (np.abs(xc.modes) <= 12)[None, :] -for name, cols in (("adv, p/B [#56]", (A1, A2_wrong)), ("adv, p", (A1, A2)), ("adv, p, gauge", (A1, A2, A3))): - M = np.stack([a[sel] for a in cols], axis=1) - coef, *_ = np.linalg.lstsq(M, y[sel], rcond=None) - res = np.linalg.norm(M @ coef - y[sel]) / np.linalg.norm(y[sel]) - print(f"fit Bkxprp ~ {name:14s}: coef", np.round(coef, 4), f"resid {res:.4f}") -M = A1 + A2 + A3 -print("fixed coefficients (1,1,1): resid %.4f" % ( - np.linalg.norm(M[sel] - y[sel]) / np.linalg.norm(y[sel]))) - -print("=== consistency of the regressors (far surfaces, |m|<=12)") -eul = sp.eulerian -def r2(a, b): - return np.linalg.norm(a[sel] - b[sel]) / np.linalg.norm(b[sel]) -print("adv (xclebsch) vs lagb-eulb %.4f" % r2(A1, spec(sp.lagrangian - eul))) -print("p-term vs Bk-(llagb-eulb) %.4f" % r2(A2, spec(sp.b_kappa_xi_perp - llagb + eul))) -print("Bk vs (llagb-eulb)+p-term %.4f" % r2(y, spec(llagb - eul) + A2)) - -T = spec(sp.b_kappa_xi_perp - llagb + eul) -scale = np.vdot(A2[sel], T[sel]) / np.vdot(A2[sel], A2[sel]) -print("pressure coefficient against corrected regressor:", scale) diff --git a/docs/derivations/pmodb_gauge/derive.py b/docs/derivations/pmodb_gauge/derive.py index d6f90cc3b..abf835099 100644 --- a/docs/derivations/pmodb_gauge/derive.py +++ b/docs/derivations/pmodb_gauge/derive.py @@ -1,18 +1,21 @@ -"""Run the pmodb gauge identities through FortSym and record the verdicts. +"""Run the pmodb gauge identities through a CAS and record the verdicts. -Usage (FortSym checkout with a built C ABI library):: +Usage:: - PYTHONPATH=$FORTSYM/python python3 docs/derivations/pmodb_gauge/derive.py + python3 docs/derivations/pmodb_gauge/derive.py -Writes ``fortsym_verdicts.json`` beside this file. ``PROVED`` means FortSym's -native zero decision returned true; ``UNKNOWN`` and ``NONZERO`` fail an -identity and are the required outcome for a negative control. No numerical -probe is accepted as evidence. +FortSym (a checkout with a built C ABI library on ``PYTHONPATH``) is used +when it is importable; otherwise plain SymPy is used. Writes +``fortsym_verdicts.json`` beside this file, with the engine recorded. +``PROVED`` means the zero test returned true; ``UNKNOWN`` and ``NONZERO`` +fail an identity and are the required outcome for a negative control. No +numerical probe is accepted as evidence. """ from __future__ import annotations import json +import signal import subprocess import sys import time @@ -21,18 +24,43 @@ HERE = Path(__file__).resolve().parent sys.path.insert(0, str(HERE)) -import fortsym # noqa: E402 -import fortsym.sympy as fs # noqa: E402 +try: + import fortsym # noqa: E402 + import fortsym.sympy as fs # noqa: E402 +except ImportError: + fortsym = None + import sympy as fs # noqa: E402 import identities # noqa: E402 +TIMEOUT = 600 # seconds per SymPy zero test + + +class _Timeout(Exception): + pass + + +def _alarm(*_): + raise _Timeout + def verdict(expr) -> str: - value = fs.sympify(expr).is_zero + if fortsym is not None: + value = fs.sympify(expr).is_zero + else: + signal.signal(signal.SIGALRM, _alarm) + signal.alarm(TIMEOUT) + try: + e = fs.expand(expr) + value = True if e == 0 else fs.simplify(fs.together(e)) == 0 or None + except _Timeout: + return "UNKNOWN" + finally: + signal.alarm(0) return {True: "PROVED", False: "NONZERO"}.get(value, "UNKNOWN") def counterexample(expr): - """Exact rational point at which FortSym decides the residual is non-zero.""" + """Exact rational point at which the engine decides the residual is non-zero.""" symbols = sorted(expr.free_symbols, key=str) for attempt in range(1, 6): point, value = {}, expr @@ -46,6 +74,8 @@ def counterexample(expr): def revision() -> str: + if fortsym is None: + return "sympy " + fs.__version__ root = Path(fortsym.__file__).resolve().parents[2] try: return subprocess.run(["git", "-C", str(root), "rev-parse", "HEAD"], @@ -80,7 +110,9 @@ def main() -> int: failed.append(label) print(f"{label:6s} {v:8s} (negative control" f"{', witness ' + str(witness) if witness else ''}) {equation}", flush=True) - out = {"engine": "fortsym native zero test via C ABI (fortsym.sympy)", + engine = ("fortsym native zero test via C ABI (fortsym.sympy)" if fortsym + else "sympy expand/simplify zero test") + out = {"engine": engine, "fortsym_revision": revision(), "checks": records, "failed": failed} (HERE / "fortsym_verdicts.json").write_text(json.dumps(out, indent=2) + "\n") print("FAILED: " + " ".join(failed) if failed else "required identities proved; see explicit open-gap records") diff --git a/docs/derivations/pmodb_gauge/fortsym_verdicts.json b/docs/derivations/pmodb_gauge/fortsym_verdicts.json index b271512dc..6c74f984d 100644 --- a/docs/derivations/pmodb_gauge/fortsym_verdicts.json +++ b/docs/derivations/pmodb_gauge/fortsym_verdicts.json @@ -1,42 +1,42 @@ { - "engine": "fortsym native zero test via C ABI (fortsym.sympy)", - "fortsym_revision": "ae635c80b42c922f5f2028c5cec67859bb8d3929", + "engine": "sympy expand/simplify zero test", + "fortsym_revision": "sympy 1.9", "checks": [ { "label": "G1.1", "equation": "J dB^psi = GPEC bwp (pointwise)", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.002 }, { "label": "G1.2", "equation": "J dB^theta = GPEC bwt (pointwise)", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.003 }, { "label": "G1.3", "equation": "J dB^zeta = GPEC bwz (pointwise)", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.006 }, { "label": "G1.4", "equation": "(xi x B)_psi = chi1 (q xi^theta - xi^zeta) = xi_s", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.001 }, { "label": "G1.5", "equation": "xi_s(xi + lam B) = xi_s(xi)", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.001 }, { "label": "G1F.1", "equation": "bwp = chi1 (m - n q) 2 pi i xsp", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.002 }, { "label": "G1F.2", @@ -54,61 +54,61 @@ "label": "G2.1", "equation": "(d_t + q d_z)(J xi^t) = d_z(J xi_s)/chi1 - d_psi(J xi^psi) + J div xi", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.001 }, { "label": "G2.2", "equation": "(d_t + q d_z)(J xi^z) = -q d_psi(J xi^psi) - d_t(J xi_s)/chi1 + q J div xi", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.005 }, { "label": "G2.3", "equation": "q xwt - xwz = J xi_s / chi1", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.001 }, { "label": "G3", "equation": "B.dB = chi1 (bvt + q bvz)/J (eulb times |B|)", "verdict": "PROVED", - "seconds": 0.009 + "seconds": 0.052 }, { "label": "G4.1", - "equation": "J xi_perp^theta = xmt - (chi1/B)^2 (xvt + q xvz)/J", + "equation": "J xi_perp^theta = chi1^2 (-f2 xsp + f3 xms/chi1) (cspl 1)", "verdict": "PROVED", - "seconds": 0.003 + "seconds": 0.469 }, { "label": "G4.2", - "equation": "J xi_perp^zeta = xmz - q (chi1/B)^2 (xvt + q xvz)/J", + "equation": "J xi_perp^zeta = q J xi_perp^theta - J xms/chi1 (cspl 2)", "verdict": "PROVED", - "seconds": 0.002 + "seconds": 0.307 }, { "label": "G4.3", "equation": "B_i xi_perp^i = 0 for xi_perp^i = xi^i - (B_j xi^j / B_k B^k) B^i", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.055 }, { "label": "G4.4", - "equation": "xi_perp^theta = xms/(chi1 q) + xmz/(J q) - (chi1/(J B))^2 (xvt + q xvz) (curvfun, gpout.f)", + "equation": "xsp d_psi B + (cspl 1/J) d_theta B = xi_perp.grad B (curvfun)", "verdict": "PROVED", - "seconds": 0.004 + "seconds": 0.325 }, { "label": "G5.1", "equation": "xi_g^zeta = 0 for lam = -xi_perp^zeta / B^zeta", "verdict": "PROVED", - "seconds": 0.006 + "seconds": 0.0 }, { "label": "G5.2", "equation": "xi_g^theta = xi_s/(chi1 q)", "verdict": "PROVED", - "seconds": 0.007 + "seconds": 0.032 }, { "label": "G5.2b", @@ -120,115 +120,97 @@ "label": "G5.3", "equation": "curl(xi_g x B) = curl(xi_perp x B) (theta component)", "verdict": "PROVED", - "seconds": 0.017 + "seconds": 0.037 }, { "label": "G5.4", "equation": "lagb (gpout.f) = dB_L(xi_g), the xi^zeta = 0 gauge", "verdict": "PROVED", - "seconds": 0.012 + "seconds": 0.005 }, { "label": "G5.5", "equation": "lagb - dB_L(xi_perp) = -(xi_perp^zeta/q) d_theta B", "verdict": "PROVED", - "seconds": 0.008 + "seconds": 0.016 }, { "label": "G5.6", "equation": "llagb (gpdiag.f) = dB_L(xi_perp)", "verdict": "PROVED", - "seconds": 0.008 + "seconds": 0.308 }, { "label": "G5.7", "equation": "dB_L(xi, div-free) = lagb + (xmz/(J q)) d_theta B", "verdict": "PROVED", - "seconds": 0.009 + "seconds": 0.002 }, { "label": "G5.8", "equation": "xi_perp^zeta = -(B_theta xi_s/chi1 + q B_psi xi^psi)/(B_theta + q B_zeta)", "verdict": "PROVED", - "seconds": 0.005 + "seconds": 0.362 }, { "label": "G7.1", "equation": "J B^2 (xi_perp.grad B + mu0 p' xi^psi/B)/B = sband xsp + tband xms/chi1, xi = e_psi", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.013 }, { "label": "G7.2", "equation": "J B^2 (xi_perp.grad B + mu0 p' xi^psi/B)/B = sband xsp + tband xms/chi1, xi = e_theta", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.013 }, { "label": "G7.3", "equation": "J B^2 (xi_perp.grad B + mu0 p' xi^psi/B)/B = sband xsp + tband xms/chi1, xi = e_zeta", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.009 }, { "label": "C1", "equation": "b.curl(xi_perp x B) + xi_perp.grad B = -B(div xi_perp + kappa.xi_perp)", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.002 }, { "label": "C2", "equation": "B kappa.xi_perp = xi_perp.grad B + xi_perp.((curl B) x B)/B", "verdict": "PROVED", - "seconds": 0.0 - }, - { - "label": "R1", - "equation": "J D_effective(reg) = (1-R)(J R xsp1 + J' xsp)", - "verdict": "PROVED", - "seconds": 0.0 - }, - { - "label": "R1a", - "equation": "J div xi_reconstructed = J(1-R**2)xsp1 + (1-R)J' xsp", - "verdict": "PROVED", - "seconds": 0.001 - }, - { - "label": "R4", - "equation": "chi1(q xi_reconstructed^theta-xi_reconstructed^zeta)=R xms", - "verdict": "PROVED", - "seconds": 0.0 - }, - { - "label": "R5", - "equation": "llagb_code - (eulb+xi_rec_perp.grad B) = (1-R) xms Bt/(chi1 q)", - "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.004 }, { "label": "R2", "equation": "J div xi (R = 1) = 0", "verdict": "PROVED", - "seconds": 0.001 + "seconds": 0.028 }, { "label": "R3", "equation": "J div dB (reg) = 2 pi i chi1 (m - n q)(1 - R) xsp1", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.005 }, { "label": "R6", "equation": "H_parallel = B(X1-U)+(B_theta Delta_zeta-B_zeta Delta_theta)/(J B)", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.01 }, { "label": "R7", "equation": "lagb+Bdivxprp+Bkxprp = -xi_rec_perp^zeta Bt/q + 2 Delta Bt/(chi1 q) + curl_delta + B(V-U) - xi_rec_perp.F/B", "verdict": "PROVED", - "seconds": 0.001 + "seconds": 0.002 + }, + { + "label": "R8", + "equation": "with V = U and Delta = 0: lagb+Bdivxprp+Bkxprp = -xi_perp^zeta Bt/q - xi_perp.F/B", + "verdict": "PROVED", + "seconds": 0.003 }, { "label": "E4.1", @@ -240,13 +222,13 @@ "label": "G8.1", "equation": "J B**2 div xi_perp = xmat Xp + ymat X + zmat A", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.008 }, { "label": "G8.2", "equation": "idcon zband2/zband3 normalization from B**2=chi1**2*(g22+2qg23+q**2g33)/J**2", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.006 }, { "label": "Q1.hamada", @@ -270,7 +252,7 @@ "label": "Q1.boozer", "equation": "T_a(xi+lambda B) = T_a(xi)+B.grad(lambda B**a)/B**a", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.006 }, { "label": "Q2.boozer", @@ -288,7 +270,7 @@ "label": "Q4", "equation": "det d(theta_H+iota h,phi_H+h)/d(theta_H,phi_H)=1+D h", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.001 }, { "label": "Q5", @@ -312,13 +294,13 @@ "label": "Q8", "equation": "delta(J/R**2)/(J/R**2)=div xi-2 xi_R/R", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.002 }, { "label": "K5", "equation": "Complex CGL pressure work equals [5/3 |D|**2+1/3 |D+3K|**2]/2", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.007 }, { "label": "P1", @@ -330,19 +312,13 @@ "label": "K1", "equation": "CGL energy density = (1/2)[(5/3) D^2 + (1/3)(D + 3K)^2]", "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.004 }, { "label": "K3", "equation": "v [drive change under xi -> xi + lam B - d/dl (v B lam)] = 0", "verdict": "PROVED", - "seconds": 0.0 - }, - { - "label": "K4", - "equation": "pmodb_file-path drive change = -mu B lam dB/dl / v", - "verdict": "PROVED", - "seconds": 0.0 + "seconds": 0.008 }, { "label": "G6", @@ -356,19 +332,18 @@ "verdict": "UNKNOWN", "negative_control": true, "counterexample": { - "Bm_psi": "8/5", - "Bm_theta": "11/6", - "J_0": "14/7", - "chi1": "17/8", - "g_11_0": "20/9", - "g_12_0": "23/10", - "g_13_0": "26/11", - "g_22_0": "29/12", - "g_23_0": "32/13", - "q_0": "35/14", - "xp_0": "38/15", - "xt_0": "41/16", - "xz_0": "44/17" + "Bm_theta": "8/5", + "J_0": "11/6", + "chi1": "14/7", + "g_11_0": "17/8", + "g_12_0": "20/9", + "g_13_0": "23/10", + "g_22_0": "26/11", + "g_23_0": "29/12", + "q_0": "32/13", + "xp_0": "35/14", + "xt_0": "38/15", + "xz_0": "41/16" } }, { @@ -377,13 +352,10 @@ "verdict": "UNKNOWN", "negative_control": true, "counterexample": { - "Bm_psi": "8/5", - "Bm_theta": "11/6", - "chi1": "14/7", - "q_0": "17/8", - "xp_0": "20/9", - "xt_0": "23/10", - "xz_0": "26/11" + "Bm_theta": "8/5", + "q_0": "11/6", + "xt_0": "14/7", + "xz_0": "17/8" } }, { diff --git a/docs/derivations/pmodb_gauge/identities.py b/docs/derivations/pmodb_gauge/identities.py index e268dee93..e8a31a9b3 100644 --- a/docs/derivations/pmodb_gauge/identities.py +++ b/docs/derivations/pmodb_gauge/identities.py @@ -221,17 +221,23 @@ def columns(sp): dBcov = S.lower(dB) B2 = S.dot(S.lower(B), B) xs = S.xi_s(xi) - xv = S.lower(xi) # GPEC xvt, xvz (covariant) - xmt, xmz = J * xi[1], J * xi[2] # GPEC xmt, xmz (J xi^i) - sigma = S.dot(xv, B) / B2 + g = S.g + xmz = J * xi[2] # GPEC xmz (J xi^zeta) + sigma = S.dot(S.lower(xi), B) / B2 xperp = [xi[i] - sigma * B[i] for i in range(3)] + # gpout_pmodb/gpdiag_pmodb cspl: J xi_perp from xi^psi and xi_s with the + # eqfun coefficients f2, f3 (equil/direct.f, covariant metric). + f2 = (g[0][1] + q * g[0][2]) / (J * B2) + f3 = (g[1][2] + q * g[2][2]) / (J * B2) + jxt = chi1 ** 2 * (-f2 * xi[0] + f3 * xs / chi1) + jxz = q * jxt - J * xs / chi1 out = [] out.append(("G3", "B.dB = chi1 (bvt + q bvz)/J (eulb times |B|)", S.at0(S.dot(dBcov, B) - chi1 * (dBcov[1] + q * dBcov[2]) / J))) - out.append(("G4.1", "J xi_perp^theta = xmt - (chi1/B)^2 (xvt + q xvz)/J", - S.at0(J * xperp[1] - (xmt - chi1 ** 2 / B2 * (xv[1] + q * xv[2]) / J)))) - out.append(("G4.2", "J xi_perp^zeta = xmz - q (chi1/B)^2 (xvt + q xvz)/J", - S.at0(J * xperp[2] - (xmz - q * chi1 ** 2 / B2 * (xv[1] + q * xv[2]) / J)))) + out.append(("G4.1", "J xi_perp^theta = chi1^2 (-f2 xsp + f3 xms/chi1) (cspl 1)", + S.at0(J * xperp[1] - jxt))) + out.append(("G4.2", "J xi_perp^zeta = q J xi_perp^theta - J xms/chi1 (cspl 2)", + S.at0(J * xperp[2] - jxz))) # B . xi_perp = 0 as a statement about arbitrary vectors: the metric-jet # form is larger than necessary, so the # projection is proved on independent covariant B_i, contravariant B^i, xi^i. @@ -241,10 +247,9 @@ def columns(sp): sig = sum(Bl[i] * xu[i] for i in range(3)) / sum(Bl[i] * Bu[i] for i in range(3)) out.append(("G4.3", "B_i xi_perp^i = 0 for xi_perp^i = xi^i - (B_j xi^j / B_k B^k) B^i", sum(Bl[i] * (xu[i] - sig * Bu[i]) for i in range(3)))) - coef = xs / (chi1 * q) + xmz / (J * q) - chi1 ** 2 / (J ** 2 * B2) * (xv[1] + q * xv[2]) - out.append(("G4.4", "xi_perp^theta = xms/(chi1 q) + xmz/(J q)" - " - (chi1/(J B))^2 (xvt + q xvz) (curvfun, gpout.f)", - S.at0(xperp[1] - coef))) + out.append(("G4.4", "xsp d_psi B + (cspl 1/J) d_theta B = xi_perp.grad B (curvfun)", + S.at0(xi[0] * S.d(Bm, 0) + jxt / J * S.d(Bm, 1) + - sum(xperp[i] * S.d(Bm, i) for i in range(3))))) # Gauge of lagb: xi_g = xi_perp + lam B with xi_g^zeta = 0. lam = -xperp[2] / B[2] xg = [xperp[i] + lam * B[i] for i in range(3)] @@ -260,8 +265,7 @@ def columns(sp): lagb = eul + Bm * (xi[0] * S.d(Bm, 0) + xs / (chi1 * q) * S.d(Bm, 1)) lagb_g = eul + Bm * grad(xg) llagb = eul + Bm * grad(xperp) - llagb_code = lagb + Bm * S.d(Bm, 1) * (xmz / (J * q) - chi1 ** 2 / (J ** 2 * B2) - * (xv[1] + q * xv[2])) + llagb_code = eul + Bm * (xi[0] * S.d(Bm, 0) + jxt / J * S.d(Bm, 1)) out.append(("G5.4", "lagb (gpout.f) = dB_L(xi_g), the xi^zeta = 0 gauge", S.at0(lagb - lagb_g))) out.append(("G5.5", "lagb - dB_L(xi_perp) = -(xi_perp^zeta/q) d_theta B", @@ -371,12 +375,12 @@ def cartesian_kappa(sp): def regularisation(sp): - """R1-R3: Park-regularised components in GPEC (reg_flag) on one harmonic. + """R2-R3: Park-regularised components in GPEC (reg_flag) on one harmonic. Diagonal model: one harmonic, J independent of theta. R = s^2/(s^2 + reg^2), s = m - n q. gpeq_sol uses xmp1 = R xsp1 in bmt/bmz and xms; gpeq_contra - builds xwt/xwz from xmp1 and multiplies by R again; gpout_pmodb uses - R xsp1 explicitly in divxprp. + builds xwt/xwz from xmp1. The pmodb columns use xmp1 and xms + (gpout_pmodb), so R2-R3 are the only regularisation identities needed. """ I, pi = sp.I, sp.pi m, n, q, chi1 = sp.symbols("m n q chi1") @@ -385,18 +389,7 @@ def regularisation(sp): s = m - n * q xwt = -(J0 * R * X1 + J1 * X + 2 * pi * I * n / chi1 * J0 * Xs) / (2 * pi * I * s) xwz = -(q * J0 * R * X1 + q * J1 * X + 2 * pi * I * m / chi1 * J0 * Xs) / (2 * pi * I * s) - xmt, xmz = R * xwt, R * xwz - jdiv_effective = J0 * R * X1 + J1 * X + 2 * pi * I * m * xmt - 2 * pi * I * n * xmz - jdiv_actual = J0 * X1 + J1 * X + 2 * pi * I * m * xmt - 2 * pi * I * n * xmz - out = [("R1", "J D_effective(reg) = (1-R)(J R xsp1 + J' xsp)", - jdiv_effective - (1-R)*(J0*R*X1+J1*X)), - ("R1a", "J div xi_reconstructed = J(1-R**2)xsp1 + (1-R)J' xsp", - jdiv_actual-J0*(1-R**2)*X1-(1-R)*J1*X), - ("R4", "chi1(q xi_reconstructed^theta-xi_reconstructed^zeta)=R xms", - chi1*(q*xmt-xmz)/J0-R*Xs)] - Bt = sp.Symbol("Bt") - out.append(("R5", "llagb_code - (eulb+xi_rec_perp.grad B) = (1-R) xms Bt/(chi1 q)", - (Xs/(chi1*q)+xmz/(J0*q)-xmt/J0)*Bt-(1-R)*Xs*Bt/(chi1*q))) + out = [] jdiv_unreg = (J0 * X1 + J1 * X + 2 * pi * I * m * xwt.subs(R, 1) - 2 * pi * I * n * xwz.subs(R, 1)) out.append(("R2", "J div xi (R = 1) = 0", jdiv_unreg)) @@ -415,6 +408,10 @@ def regularized_column_budget(sp): U is the derivative selected by gpeq_sol; V is the derivative used in gpout divergence. Delta=xi_s-chi1(q xi_rec^theta-xi_rec^zeta). + This is bookkeeping under the stated model relations (BDcode, BKcode), + not a proof that the code satisfies them. In gpout_pmodb, V = U + (xmp1) and xi_perp is built from xi_s itself (G4.2), so Delta = 0 and + R8 gives the budget of the output columns. F denotes xi_rec_perp.((curl B)xB-mu0 grad p)/B, in Tesla. Metric/field-norm consistency and C1/C2 are assumed; finite truncation errors are not folded into the physical force residual. @@ -439,6 +436,8 @@ def regularized_column_budget(sp): budget=-Z*Bt/q+2*g+curl_delta+B*(V-U)-F out.append(("R7", "lagb+Bdivxprp+Bkxprp = -xi_rec_perp^zeta Bt/q + 2 Delta Bt/(chi1 q) + curl_delta + B(V-U) - xi_rec_perp.F/B", Lagcode+BDcode+BKcode-budget)) + out.append(("R8", "with V = U and Delta = 0: lagb+Bdivxprp+Bkxprp = -xi_perp^zeta Bt/q - xi_perp.F/B", + (Lagcode+BDcode+BKcode).subs({V: U, Delta: 0, Dt: 0, Dz: 0})+Z*Bt/q+F)) return out @@ -460,7 +459,7 @@ def cgl(sp): def drive_gauge(sp): - """K3, K4: gauge (in)variance of the Logan (2013) Eq. 5 bounce drive. + """K3: gauge invariance of the Logan (2013) Eq. 5 bounce drive. Along a field line (arc length l, M = 1): v = sqrt(2(E - mu B)), integrand (1/v)[(2E - 3 mu B) dB/B + (2E - 2 mu B) div xi]. Under xi -> xi + lam B, @@ -474,14 +473,10 @@ def drive_gauge(sp): dBdl, dlam = B.diff(l), lam.diff(l) change = ((2 * E - 3 * mu * B) * lam * dBdl + (2 * E - 2 * mu * B) * B * dlam) / v total = (v * B * lam).diff(l) - # pmodb_file path: dbob = lagb/B (gauge xi_g), divx = -(lagb + Bkxprp)/B. - change_pmodb = ((2 * E - 3 * mu * B) * lam * dBdl - (2 * E - 2 * mu * B) * lam * dBdl) / v # Multiplied by v > 0 (away from turning points): FortSym leaves the # unmultiplied sqrt quotient UNKNOWN. return [("K3", "v [drive change under xi -> xi + lam B - d/dl (v B lam)] = 0", - ((change - total) * v).subs(l, 0)), - ("K4", "pmodb_file-path drive change = -mu B lam dB/dl / v", - (change_pmodb + mu * B * lam * dBdl / v).subs(l, 0))] + ((change - total) * v).subs(l, 0))] def divergence_matrices(sp): @@ -542,7 +537,7 @@ def gauge_coordinates(sp): 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. + # Same physical scalar and relabeling, composed through Hamada. 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", @@ -617,9 +612,8 @@ def open_gaps(sp): """Identities FortSym does not decide; their content is covered elsewhere. G6 is the coordinate form of C1 + C2 (proved on Cartesian 1-jets) with - det g = J^2 imposed through g_33. On ae635c80 the zero test still returns UNKNOWN for this large rational - residual; the removed variable cap is no longer the explanation; it is recorded, not - counted as proved. + det g = J^2 imposed through g_33. Neither FortSym nor SymPy decides this + large rational residual; it is recorded, not counted as proved. """ return force_balance_identity(sp) diff --git a/docs/pmodb_gauge.rst b/docs/pmodb_gauge.rst index 07aef9af0..7ba8fae94 100644 --- a/docs/pmodb_gauge.rst +++ b/docs/pmodb_gauge.rst @@ -5,11 +5,10 @@ Perturbed field strength (pmodb) and gauge ****************************************** This page states what the columns of ``gpec_pmodb_n#.out`` are, how they are -related, and which of them a field-strength consumer needs. Verified identities -are proved with the FortSym native zero test by -``docs/derivations/pmodb_gauge/derive.py``; the verdicts are stored in -``docs/derivations/pmodb_gauge/fortsym_verdicts.json`` and the labels (G1, -G5.5, ...) refer to that file. +related, and which of them a field-strength consumer needs. The identities +are proved symbolically by ``docs/derivations/pmodb_gauge/derive.py``; the +verdicts are stored in ``docs/derivations/pmodb_gauge/fortsym_verdicts.json`` +and the labels (G1, G5.5, ...) refer to that file. .. contents:: On this page :local: @@ -36,10 +35,12 @@ any straight-field-line Jacobian :math:`\mathcal{J}(\psi,\theta)` selected by the in-surface displacement across :math:`\mathbf B` (``xms``; ``gpec_xclebsch`` writes :math:`\xi^\alpha=\xi_s/\chi_1`). A field-aligned -displacement changes neither :math:`\xi^\psi` nor :math:`\xi_s` (G1.5). The -parallel part is fixed in ``gpeq_contra`` by +displacement changes neither :math:`\xi^\psi` nor :math:`\xi_s` (G1.5), so +:math:`(\xi^\psi,\xi_s)` determine :math:`\boldsymbol\xi_\perp` completely. +The parallel part used elsewhere in GPEC is fixed in ``gpeq_contra`` by :math:`\nabla\cdot\boldsymbol\xi=0` (G2), which gives ``xmt`` :math:`=\mathcal J\xi^\theta` and ``xmz`` :math:`=\mathcal J\xi^\zeta`. +The pmodb columns do not use ``xmt`` or ``xmz``. For every ``jac_type`` (G1, G1F) the Jacobian-weighted perturbed field of ``gpeq_sol`` is :math:`\mathcal J\,\nabla\times(\boldsymbol\xi\times\mathbf B)`. @@ -48,10 +49,6 @@ For every ``jac_type`` (G1, G1F) the Jacobian-weighted perturbed field of The pmodb columns ================= -The displacement identifications in this section assume unregularized -consistent components. The literal lagb expression is valid for both modes; -the regularized closure defect is given by R4-R5 below. - ``eulb`` :math:`\delta B_E=\hat{\mathbf b}\cdot\delta\mathbf B`, the Eulerian field strength (G3). @@ -68,9 +65,21 @@ the regularized closure defect is given by R4-R5 below. relabels the harmonics; it does not change this gauge. ``Bdivxprp`` - :math:`B\,\nabla\cdot\boldsymbol\xi_\perp`, with - :math:`\boldsymbol\xi_\perp=\boldsymbol\xi-(\boldsymbol\xi\cdot\mathbf B)\mathbf B/B^2` - built from ``xmt``, ``xmz`` and the covariant ``xvt``, ``xvz`` (G4.1-G4.3). + :math:`B\,\nabla\cdot\boldsymbol\xi_\perp`. ``gpout_pmodb`` forms + :math:`\mathcal J\boldsymbol\xi_\perp` directly from :math:`\xi^\psi` and + :math:`\xi_s` (G4.1, G4.2), + + .. math:: + + \mathcal J\xi_\perp^\theta=\chi_1^2\left(-f_2\,\xi^\psi+f_3\,\xi_s/\chi_1\right), + \qquad + \mathcal J\xi_\perp^\zeta=q\,\mathcal J\xi_\perp^\theta-\mathcal J\xi_s/\chi_1, + + with the equilibrium coefficients + :math:`f_2=(g_{\psi\theta}+qg_{\psi\zeta})/(\mathcal JB^2)` and + :math:`f_3=(g_{\theta\zeta}+qg_{\zeta\zeta})/(\mathcal JB^2)` (``eqfun`` + columns 2 and 3, covariant metric). The same projection gives the + ``xmat``-``zmat`` matrices of ``idcon_build`` (G8). ``Bkxprp`` :math:`\boldsymbol\xi_\perp\cdot\nabla B+\mu_0p'\xi^\psi/B` (G4.4). This is @@ -78,13 +87,14 @@ the regularized closure defect is given by R4-R5 below. equilibrium force balance: exactly, :math:`B\boldsymbol\kappa\cdot\boldsymbol\xi_\perp=\boldsymbol\xi_\perp\cdot\nabla B +\boldsymbol\xi_\perp\cdot[(\nabla\times\mathbf B)\times\mathbf B]/B` (C2). - The PENTRC/GPEC.jl matrices ``smat``, ``tmat`` (``idcon_build``) form the + The PENTRC matrices ``smat``, ``tmat`` (``idcon_build``) form the same quantity times :math:`\mathcal J B^2` from :math:`\xi^\psi` and :math:`\xi^\alpha` (G7). ``llagb`` (``gpec_diagnostics_pmodb``, ``pmodbmn_flag``) :math:`\delta B_L(\boldsymbol\xi_\perp)=\delta B_E+\boldsymbol\xi_\perp\cdot\nabla B` - (G5.6). ``lllagb`` is :math:`-(`\ ``Bdivxprp``\ :math:`+`\ ``Bkxprp``\ :math:`)`. + (G5.6). ``lllagb`` is :math:`-(`\ ``Bdivxprp``\ :math:`+`\ ``Bkxprp``\ :math:`)`, + and ``Bkxprp`` :math:`=` ``llagb`` :math:`-` ``eulb`` :math:`+\mu_0p'\xi^\psi/B` (P1). Relations between the columns @@ -99,8 +109,7 @@ Gauge of ``lagb`` (G5.5, G5.8): \xi_\perp^\zeta=-\frac{B_\theta\xi_s/\chi_1+qB_\psi\xi^\psi}{B_\theta+qB_\zeta}. Kinematic identity (C1; any :math:`\boldsymbol\xi_\perp`, no force balance) -and its pmodb form (C2 with G4; the optional direct contraction in the -GPEC chart with :math:`\det g=\mathcal J^2`, G6, returns UNKNOWN on ae635c80 and is recorded as an open gap, not as proved): +and its pmodb form (C2 with G4): .. math:: @@ -109,7 +118,9 @@ GPEC chart with :math:`\det g=\mathcal J^2`, G6, returns UNKNOWN on ae635c80 and =-(\mathrm{Bdivxprp}+\mathrm{Bkxprp})-\boldsymbol\xi_\perp\cdot\mathbf F/B, \qquad \mathbf F=(\nabla\times\mathbf B)\times\mathbf B-\mu_0\nabla p. -Hence +C1 and C2 are proved on Cartesian 1-jets. Their direct contraction in the +GPEC chart with :math:`\det g=\mathcal J^2` (G6) is too large for either CAS +to decide and is recorded as an open gap, not as proved. Hence .. math:: @@ -118,112 +129,64 @@ Hence ``lagb`` :math:`=-(`\ ``Bdivxprp``\ :math:`+`\ ``Bkxprp``\ :math:`)` is therefore not an identity. The first term is a gauge term of first order in -the perturbation. Its size depends on :math:`B_\theta/(qB_\zeta)`, the -Boozer :math:`B_\psi` and :math:`\partial_\theta B` and the displacement spectrum; growth toward the edge is a hypothesis -for a particular equilibrium, not a general bound. The second term vanishes for an exact equilibrium. - -Regularisation (``reg_flag``, Park et al. 2009 Eq. 10). ``gpeq_sol`` uses -:math:`R\,\partial_\psi\xi^\psi`, :math:`R=s^2/(s^2+\sigma^2)`, -:math:`s=m-nq`, in the field and in ``xms``; ``gpeq_contra`` multiplies -``xwt``, ``xwz`` by :math:`R` again. On one harmonic with -:math:`\theta`-independent :math:`\mathcal J` (R1-R3): - -.. math:: - - \mathcal J D_{effective}=(1-R)(\mathcal J R\,\partial_\psi\xi^\psi+\mathcal J'\xi^\psi),\qquad - \mathcal J\nabla\cdot\delta\mathbf B=2\pi i\chi_1 s(1-R)\,\partial_\psi\xi^\psi. - -Here the displayed R1 formula is the *effective* divergence diagnostic, -using :math:`R\partial_\psi\xi^\psi` for its radial derivative. The actual -divergence of the reconstructed displacement is instead (R1a) +the perturbation; its size depends on :math:`B_\theta/(qB_\zeta)`, +:math:`B_\psi`, :math:`\partial_\theta B` and the displacement spectrum. +The second term vanishes for an exact equilibrium. In the output, finite +Fourier truncation and spline error add to both. -.. math:: - - \mathcal J\nabla\cdot\boldsymbol\xi_{rec} - =\mathcal J(1-R^2)\partial_\psi\xi^\psi+(1-R)\mathcal J'\xi^\psi. - -Also :math:`\chi_1(q\xi_{rec}^\theta-\xi_{rec}^\zeta)=R\xi_s` (R4). -The diagnostic llagb also differs from the projection of the reconstructed -contravariant displacement by (R5) - -.. math:: - \mathrm{llagb}_{code}-(\delta B_E+\boldsymbol\xi_{rec,\perp}\cdot\nabla B) - =(1-R)\frac{\xi_s}{\chi_1 q}\partial_\theta B. +Regularisation +============== -These formulas are for the diagonal model only. Near rational surfaces the -regularized field, displacement and effective derivative therefore do not -form one exact ideal-induction pair. This is a limitation of the smoothing -policy, not evidence of resistive physics. The unconditional smoothing in -``gpout_pmodb`` even with ``reg_flag=f`` is a separate verified output bug. - - -A general-J continuum error budget (verified R6-R7) is available without -assuming a diagonal regularizer. Let U be the derivative selected by -``gpeq_sol``, V the derivative used by pmodb divergence, and -:math:`\Delta_s=\xi_s-\chi_1(q\xi_{rec}^\theta-\xi_{rec}^\zeta)`. -Writing :math:`B_\theta,B_\zeta` for covariant components gives +With ``reg_flag`` (Park et al. 2009 Eq. 10), ``gpeq_sol`` replaces the +radial derivative by :math:`R\,\partial_\psi\xi^\psi`, +:math:`R=s^2/(s^2+\sigma^2)`, :math:`s=m-nq`, in the field and in ``xms``. +The regularised field is then not divergence free (R3): .. math:: - \mathrm{lagb}+\mathrm{Bdivxprp}+\mathrm{Bkxprp} - =-\frac{\xi_{rec,\perp}^\zeta}{q}\partial_\theta B - +\frac{2\Delta_s}{\chi_1q}\partial_\theta B - +\frac{B_\theta\partial_\zeta\Delta_s-B_\zeta\partial_\theta\Delta_s} - {\mathcal J B} - +B(V-U)-\frac{\boldsymbol\xi_{rec,\perp}\cdot\mathbf F}{B}. - -This uses consistent continuum geometry and the field norm. Finite Fourier -truncation and spline errors are additional numerical errors, not part of F. -The derivative term cancels when V=U, even if both are regularized; in the -old reg_flag=f output V differs from U. With regularization enabled, -Delta_s can be nonzero and leaves the two explicit closure terms. Thus a -large derivative in Bdivxprp alone does not measure its contribution to the -net mismatch. This budget separates the candidate causes of finding (i), -without claiming a numerical allocation of the original edge residual. + \mathcal J\nabla\cdot\delta\mathbf B=2\pi i\chi_1 s(1-R)\,\partial_\psi\xi^\psi. +``Bdivxprp`` uses the same derivative ``xmp1`` as the field, and +:math:`\boldsymbol\xi_\perp` is built from the same ``xms``, so the +regularisation cancels between the columns and the relation above holds +unchanged with the regularised quantities (R6-R8: a bookkeeping +decomposition under the stated model, not a separate proof about the code). +Near rational surfaces the regularised field and displacement do not form +an exact ideal-induction pair. This is a limitation of the smoothing +policy, not evidence of resistive physics. With ``reg_flag=f`` the +columns are unregularised and singular at rational surfaces. What consumers need =================== -PENTRC, GPEC.jl ``KineticForces`` +PENTRC The drive of Logan et al. (2013) Eqs. 4-5 needs :math:`\delta B` and :math:`\nabla\cdot\boldsymbol\xi` of the *same* displacement; the - combination is then invariant under :math:`\boldsymbol\xi\to\boldsymbol\xi+\lambda\mathbf B`. - The default PENTRC path (``peq_file`` = ``gpec_xclebsch``, ``set_peq``) - uses :math:`\delta B/B=-(\nabla\cdot\boldsymbol\xi_\perp+\boldsymbol\kappa\cdot\boldsymbol\xi_\perp)` - and :math:`\nabla\cdot\boldsymbol\xi_\perp` from the ``smat``-``zmat`` - matrices, i.e. a consistent :math:`\boldsymbol\xi_\perp` pair. The - ``pmodb_file`` path (``read_pmodb``) pairs ``lagb`` with - :math:`-(\mathrm{lagb}+\mathrm{Bkxprp})/B`, which is not the divergence of - :math:`\boldsymbol\xi_g`. Under :math:`\boldsymbol\xi\to\boldsymbol\xi+\lambda\mathbf B` - the consistent drive changes by a total derivative along the field line - (K3), while this pairing changes it by - :math:`-\mu B\lambda\,\partial_l B/v_\parallel` (K4). + combination is then invariant under + :math:`\boldsymbol\xi\to\boldsymbol\xi+\lambda\mathbf B` up to a total + derivative along the field line (K3). Both PENTRC input paths use the + :math:`\boldsymbol\xi_\perp` pair + :math:`\delta B/B=-(\nabla\cdot\boldsymbol\xi_\perp+\boldsymbol\kappa\cdot\boldsymbol\xi_\perp)` + and :math:`\nabla\cdot\boldsymbol\xi_\perp`: the default path + (``peq_file`` = ``gpec_xclebsch``, ``set_peq``) from the ``smat``-``zmat`` + matrices, and the ``pmodb_file`` path (``read_pmodb``) from + :math:`-(\mathrm{Bdivxprp}+\mathrm{Bkxprp})` and ``Bdivxprp``. The + ``lagb`` column is not read. Regenerate ``pmodb_file`` inputs with GPEC + 1.6 or later. Bounce-averaged codes in perturbed Boozer or Hamada coordinates (NEO-RT, NEO-2, Shaing-type) need :math:`\delta B` in the perturbed straight-field-line coordinates. That follows from :math:`\delta B_L(\boldsymbol\xi_\perp)` (``llagb``, or :math:`-(\mathrm{Bdivxprp}+\mathrm{Bkxprp})` away from rational surfaces) and :math:`\nabla\cdot\boldsymbol\xi_\perp` by a magnetic differential - equation for the parallel gauge. ``lagb`` must not be used as - :math:`\delta B_L(\boldsymbol\xi_\perp)`. + equation for the parallel gauge (next section). ``lagb`` must not be + used as :math:`\delta B_L(\boldsymbol\xi_\perp)`. -Reproducing the proofs -====================== - -:: - - PYTHONPATH=$FORTSYM/python python3 docs/derivations/pmodb_gauge/derive.py - -FortSym (https://github.com/lazy-fortran/fortsym) is needed only to rerun the -proofs, not to build or run GPEC. - - -Target gauges in GPEC units (verified Q1-Q8) -============================================ +Target gauges in GPEC units (Q1-Q8) +=================================== Let :math:`D=\partial_\theta+q\partial_\zeta`, so :math:`\mathbf B\cdot\nabla=(\chi_1/\mathcal J)D`. @@ -257,64 +220,42 @@ with ``jac_out`` alone does not solve this equation. For a PEST chart :math:`\mathcal J/R^2` is a flux function; its linear perturbation is :math:`\nabla\cdot\boldsymbol\xi-2\xi_R/R` (Q8). PEST output is not automatically either target gauge. Transform the -equilibrium chart first, then solve the target equation. The Hamada scalar -on Boozer equilibrium angles uses :math:`T_2=2 b_H/B`, matching -rmp_torque c3614334, Eq. 4.14. Q5 proves composition of parallel shifts; -Q4 proves the angular-map determinant. The native regression compares the -finite spectral operator directly with that consumer's kernel. - - -Verified defects and numerical attribution -========================================== - -* ``gpout_pmodb`` smooths the radial derivative even with ``reg_flag=f``; - use ``xmp1_mn`` selected by ``gpeq_sol``. No change for ``reg_flag=t``. -* ``gpdiag_pmodb`` uses a stale ``jac`` in its first theta loop. The local - Jacobian must be assigned before projecting out the parallel component. -* PENTRC ``read_pmodb`` overwrites the supplied divergence with - ``-(lagb+Bkxprp)``. Use the supplied positive ``Bdivxprp`` and - ``Bkxprp`` to form ``b_perp=-(Bdivxprp+Bkxprp)`` instead. -* Both CGL torque paths recover ``kappa.xi_perp`` as - ``-(dbob+divx)/2`` although their input convention is - ``dbob=-(divx+kapx)``. Remove the extra half. K1 and K5 prove the - double-adiabatic pressure work for real and complex perturbations. - -Finding (ii) from the TC24 reconstruction is a missing division by B in -the pressure *regressor*, not a GPEC normalization error (P1). For the -retained ITER n=3 output the original fit gives 0.1677; the corrected -pressure term plus the gauge term fits coefficients (1,1,1) with relative -residual about 1e-4. This is numerical corroboration only; ``docs/derivations/pmodb_gauge/check_output.py`` -reproduces the fit from the pmodb, xclebsch, diagnostic and IDCON ASCII tables. - -Finding (i) combines the gauge term, effective-derivative smoothing, -finite spectral/spline error, and the equilibrium force residual. The -identity above separates them; it does not prove a quantitative allocation -of the originally reported 0.5-19 percent to any one cause. The explicit -pressure term in curvature has coefficient one. Edge growth alone does -not diagnose nonideal plasma physics. Finding (iii) is consistent with the -explicit radial derivative in divergence and its smoothing near rational -surfaces. No change to the smoothing model itself is proposed here. +equilibrium chart first, then solve the target equation. Q4 proves the +angular-map determinant and Q5 the composition of parallel shifts. + + +Reproducing the proofs +====================== + +:: + + python3 docs/derivations/pmodb_gauge/derive.py + +This needs only SymPy and takes a few minutes, most of it the G6 attempt. +FortSym (https://github.com/lazy-fortran/fortsym) is used instead when it is +on ``PYTHONPATH``. Neither is needed to build or run GPEC. + + +History +======= + +Before GPEC 1.6.0, ``Bdivxprp`` and ``Bkxprp`` subtracted a +Fourier-truncated parallel displacement, ``gpout_pmodb`` regularised the +divergence even with ``reg_flag=f``, ``gpdiag_pmodb`` used a stale Jacobian, +PENTRC's ``read_pmodb`` paired ``lagb`` with an inconsistent divergence, and +both PENTRC CGL paths halved the curvature term. These were fixed in +1.6.0 (#340-#343). Source map and literature -========================== - -Verified code contracts (line numbers at Fortran base 3622c25c): -``gpec/idcon.f:93`` (chi1=2pi psio), ``:616`` (projection geometry), -``:677`` (pressure in Pa), ``:1026`` (curvature and divergence matrices), -``gpec/gpeq.f:131`` (regularization selection), ``:208`` (surface-displacement -closure), ``gpec/gpout.f:3041`` (eulb/lagb), ``:3050`` (divergence and -pressure curvature), ``gpec/gpdiag.f:1012`` (diagnostic projection), -``:1061`` (llagb), ``pentrc/inputs.f90:376`` (ASCII columns), -``:396`` (old gauge mixing), ``:856`` (perpendicular matrix pair), -``pentrc/torque.F90:379`` (old CGL half). - -Julia base 0304fd54c: ``src/KineticForces/KineticForcesStructs.jl:437`` -uses the correct perpendicular dbob convention; -``src/KineticForces/Torque.jl:295`` repeats the CGL half. -``src/PerturbedEquilibrium/FieldReconstruction.jl``'s xi_s is Glasser's -surface displacement, not xi^zeta. Eq. 18 describes its elimination in the -Euler-Lagrange system; Eq. 9 defines it. +========================= + +Code: ``gpeq_sol`` and ``gpeq_contra`` (``gpec/gpeq.f``, regularisation and +the divergence-free closure), ``gpout_pmodb`` (``gpec/gpout.f``, the pmodb +columns), ``gpdiag_pmodb`` (``gpec/gpdiag.f``, ``llagb`` and ``lllagb``), +``idcon_build`` (``gpec/idcon.f``, ``smat``-``zmat``), ``read_pmodb`` and +``set_peq`` (``pentrc/inputs.f90``), and ``direct_run`` (``equil/direct.f``, +``eqfun``). The physical definitions were checked against Park 2009, PoP 16 056115, Eqs. 7 and 10; Park 2011, PoP 18 110702, Eqs. 4, 10 and 11; @@ -322,11 +263,4 @@ Park and Logan 2017, PoP 24 032505, Eq. 13 and p. 6; Logan and Park 2013, PoP 20 122507, Eqs. 4 and 5; and Boozer and Nuehrenberg 2006, PoP 13 102501, Eqs. 42 and 43. The latter's Delta includes twice the curvature; its Eq. 43 is Eulerian and is not the perpendicular Lagrangian -scalar. These are source contracts, separate from the native algebra proofs. - -FortSym limits: #84, #85 and #86 are present at ae635c80. The direct G6 -contraction remains UNKNOWN (33 symbols, rational metric jets); no probe -is promoted to proof. #91 (coefficient extraction after zero substitution) -and #92 (radical power reduction) are not exercised by the accepted -polynomial/rational gauge proofs. C1/C2 supply the invariant kinematic -identity independently of G6. +scalar. From 09c85cddd9d5783f7c3afbe383772a29ba1322d0 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 24 Oct 2026 01:32:02 -0400 Subject: [PATCH 7/8] DOCS - CLEANUP - Correct stale and wrong statements across the docs for 1.6 Factual pass over the existing pages, each edit checked against the source; no new sections. outputs.rst: euler.bin, lowercase gpec_* file names, the real flags for vpbrzphi/vsbrzphi and bnormal_spectrum.bin, xbnormal columns, output formats and when the control/singfld outputs are written. sign_conventions.rst: the Jacobian powers give 1/J, and the jac_type table gains park and other. installation.rst and install/README: X11_HOME, NETCDF_C_DIR is derived, allf builds the same modules as all, OMPFLAG is automatic. Broken links in namelists.rst and source_documentation.rst; the DIII-D kinetic and resistive READMEs are linked from the index; example READMEs no longer claim rlar output while rlar_flag is off. Template comments: jac_type options, ahb_flag. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VMDGyfQKRZQMdk3NJXsKQE --- docs/examples.rst | 8 +++++ docs/examples/DIIID_ideal_example/README | 4 +-- docs/examples/DIIID_kinetic_example/README | 2 +- docs/index.rst | 2 ++ docs/installation.rst | 8 ++--- docs/namelists.rst | 4 +-- docs/outputs.rst | 40 +++++++++++----------- docs/sign_conventions.rst | 12 +++++-- docs/source_documentation.rst | 2 +- input/dcon.in | 2 +- input/equil.in | 2 +- install/README | 4 +-- 12 files changed, 54 insertions(+), 36 deletions(-) diff --git a/docs/examples.rst b/docs/examples.rst index e9f79fa29..f83b9830e 100644 --- a/docs/examples.rst +++ b/docs/examples.rst @@ -3,3 +3,11 @@ .. _DIIID_ideal_example: .. include:: examples/DIIID_ideal_example/README + +.. _DIIID_kinetic_example: + +.. include:: examples/DIIID_kinetic_example/README + +.. _DIIID_resistive_example: + +.. include:: examples/DIIID_resistive_example/README diff --git a/docs/examples/DIIID_ideal_example/README b/docs/examples/DIIID_ideal_example/README index bc4e65dcb..04b3fbb68 100644 --- a/docs/examples/DIIID_ideal_example/README +++ b/docs/examples/DIIID_ideal_example/README @@ -3,7 +3,7 @@ GPEC DIII-D Ideal Example This example uses a baseline DIII-D-like H-mode. The equilibrium is an artificial one made with TokaMaker to have 0D parameters matching typical DIII-D H-mode values, and the kinetic profiles are self-consistent with the equilibrium. -Note that the coil and pentrc input namelists use relative paths to the package's data files. If you would like to run the example in a different location, you must update those paths appropriately. +Note that the coil input namelist uses a relative path to the package's data files. If you would like to run the example in a different location, you must update that path appropriately. Following a typical workflow we first ``cd`` to the example's directory, then run DCON:: @@ -39,4 +39,4 @@ Finally, we run PENTRC to compute the neoclassical toroidal viscosity (NTV) torq PENTRC START => ... PENTRC STOP => Normal termination. -The full (``fgar``) and trapped-only (``tgar``) torque calculations agree reasonably well, consistent with the expectation that banana-orbit particles dominate the nonambipolar transport. The reduced large-aspect-ratio (``rlar``) result typically differs significantly from the general-aspect-ratio calculations for shaped, high-beta H-mode plasmas. +The full (``fgar``) and trapped-only (``tgar``) torque calculations agree reasonably well, consistent with the expectation that banana-orbit particles dominate the nonambipolar transport. The reduced large-aspect-ratio (``rlar``) result (enabled with ``rlar_flag=t``; off in this example) typically differs significantly from the general-aspect-ratio calculations for shaped, high-beta H-mode plasmas. diff --git a/docs/examples/DIIID_kinetic_example/README b/docs/examples/DIIID_kinetic_example/README index beb77909c..26653414e 100644 --- a/docs/examples/DIIID_kinetic_example/README +++ b/docs/examples/DIIID_kinetic_example/README @@ -30,4 +30,4 @@ Finally, we run PENTRC for the full neoclassical torque calculation:: PENTRC START => ... PENTRC STOP => Normal termination. -The full general-aspect-ratio (``fgar``) and trapped-only (``tgar``) torque results can be compared to assess the trapped-particle dominance. The ``data_dir`` in ``pentrc.in`` points to the precomputed special function matrix used by the reduced large-aspect-ratio (``rlar``) method, which serves as a cross-check. +The full general-aspect-ratio (``fgar``) and trapped-only (``tgar``) torque results can be compared to assess the trapped-particle dominance. The ``data_dir`` in ``pentrc.in`` points to the precomputed special function matrix used by the reduced large-aspect-ratio (``rlar``) method, which can be enabled with ``rlar_flag=t`` as a cross-check. diff --git a/docs/index.rst b/docs/index.rst index aa1df6a0c..09be708fc 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -67,3 +67,5 @@ The GPEC repository contains a number of examples intended to provide both consi Example console commands and output are given below. - :ref:`DIII-D Ideal Example ` +- :ref:`DIII-D Kinetic Example ` +- :ref:`DIII-D Resistive Example ` diff --git a/docs/installation.rst b/docs/installation.rst index f7b488392..43813b04e 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -50,7 +50,7 @@ Important make targets ---------------------- * all (default): build all main executables and libraries -* allf: same as all but Fortran-only modules +* allf: currently builds the same modules as all * clean: remove built objects and module files * clear: remove example outputs in docs/examples * realclean: deep clean including dependencies built into deps @@ -82,7 +82,7 @@ netCDF libraries: * NETCDF_FORTRAN_HOME: preferred when C and Fortran installs are split * NETCDFHOME or NETCDF_DIR: netCDF installation prefix * NETCDFINC: netCDF include path -* NETCDF_C_HOME or NETCDF_C_DIR: netCDF-C location if separate +* NETCDF_C_HOME: netCDF-C location if separate Dependency install location: @@ -180,8 +180,8 @@ Then build with:: X11 (xdraw) ^^^^^^^^^^^ -If xdraw fails to link on macOS, ensure XQuartz is installed and add include -and library paths in xdraw/makefile, for example:: +If xdraw fails to link on macOS, ensure XQuartz is installed and set +X11_HOME (read by install/DEFAULTS.inc), for example:: brew install libx11 export X11_HOME=/opt/homebrew diff --git a/docs/namelists.rst b/docs/namelists.rst index 4629dade7..491cd44d2 100644 --- a/docs/namelists.rst +++ b/docs/namelists.rst @@ -1,5 +1,5 @@ -Fully annotated namelist input files can be found in the $GPECHOME/input directory and any of the examples in $GPECHOME/docs/examples directories of any `public installation `. +Fully annotated namelist input files can be found in the $GPECHOME/input directory and any of the examples in $GPECHOME/docs/examples directories of any :ref:`public installation `. The major input files requiring user interaction are as follows. @@ -40,7 +40,7 @@ This namelist controls the MATCH code, including which DCON eigenfunction (sol_n The PENTRC Namelist ==================== -This namelist sets the inputs, controls the behavior, and determines the outputs of the GPEC code. +This namelist sets the inputs, controls the behavior, and determines the outputs of the PENTRC code. .. include:: ../input/pentrc.in :literal: diff --git a/docs/outputs.rst b/docs/outputs.rst index f642bd747..22d35140b 100644 --- a/docs/outputs.rst +++ b/docs/outputs.rst @@ -125,46 +125,46 @@ The binary GPEC outputs can be viewed using the commandxdraw filenamewhere filen Output Files ============ -The GPEC outputs are entirely specified by flags (bool types t or f) set in the GPEC_output section of GPEC.in. All outputs are ASCII files, and can be grouped into two major categories. +The GPEC outputs are entirely specified by flags (bool types t or f) set in the gpec_output section of gpec.in. Outputs are written as ASCII tables (ascii_flag), xdraw binary files (bin_flag, bin_2d_flag) and netcdf files (netcdf_flag), and can be grouped into the categories below. ASCII File Outputs ------------------ -A number of the flag options in GPEC.in instruct GPEC to output ASCII file data. Some of these outputs are always available. Some, however, require a input error field instead of a hard coded harmonic_flag call. Both groups are listed in detail here. +A number of the flag options in gpec.in instruct GPEC to output ASCII file data. Some of these outputs are always available. Some, however, require an external 3D field input (data_flag, coil_flag, harmonic_flag) or a single eigenmode (mode_flag). Both groups are listed in detail here. Outputs Always Available -These GPEC outputs can always be obtained from a equil.bin file output from DCON. +These GPEC outputs can always be obtained from a euler.bin file output from DCON. -GPEC_response_n#.out +gpec_response_n#.out **Flag** resp_flag **Info** Energy for DCON eigenmodes and stability indices. Eigenvalues and eigenvectors for vacuum and plasma inductance (virtual casing currents to fields), plasma permeability (external fields to total fields), and plasma reluctance. -GPEC_singcoup_matrix_n#.out +gpec_singcoup_matrix_n#.out **Flag** singcoup_flag **Info** The coupling matrix to resonant fields, coupling matrix to singular currents, and to island half-widths is given for each rational surface within the plasma (q=2, 3, etc) for each surface the real and imaginary coupling constants are given for each poloidal mode number on the control surface. -GPEC_singcoup_svd_n#.out +gpec_singcoup_svd_n#.out **Flag** singcoup_flag - **Info** The SVD singular values (s) and eigen vectors for each coupling matrix in GPEC_singcoup_matrix_n#.out. Large s corresponds to large amplification, with the largest (most important mode) listed at the top. The results should be dotted with the unweighted normal field spectrum to give physical meaning. + **Info** The SVD singular values (s) and eigen vectors for each coupling matrix in gpec_singcoup_matrix_n#.out. Large s corresponds to large amplification, with the largest (most important mode) listed at the top. The results should be dotted with the unweighted normal field spectrum to give physical meaning. Outputs Available When Error Field is Provided -These outputs are only available when an external error field file is provided as an input to GPEC. This means GPEC.in must have the data_flag turned on and a infile specified. +These outputs are only available when an external 3D field is provided as an input to GPEC. This means gpec.in must have data_flag (with an infile specified), coil_flag, harmonic_flag, or mode_flag turned on. -GPEC_control_n#.out +gpec_control_n#.out - **Flag** + **Flag** none (always written with ascii_flag) **Info** The Plasma response for an external perturbation on the control surface. This includes the vacuum energy, surface energy, plasma energy, real and imaginary vacuum input \mathbf{B}_{in} and total field on plasma boundary\mathbf{B}_{out}as a function of poloidal mode number. -GPEC_singfld_n#.out +gpec_singfld_n#.out **Flag** singfld_flag @@ -174,21 +174,21 @@ GPEC_singfld_n#.out **Info** Additional section showing the overlap field (Tesla) and overlap percentage for each eigenmode in the singcoup_svd output. See :ref:`coupling_overlap` for definitions. -GPEC_pmodb_n#.out +gpec_pmodb_n#.out **Flag** pmodb_flag **Info** Perturbed field strength (real and imaginary, Tesla) for each poloidal mode number at each value of \Psi_{N} output: the Eulerian \delta B (eulb), the Lagrangian \delta B of the displacement with zero \xi^\zeta in the working (jac_type) coordinates (lagb), B\nabla\cdot\xi_\perp (Bdivxprp) and the force-balance form of B\kappa\cdot\xi_\perp (Bkxprp). lagb is not the \xi_\perp-only Lagrangian \delta B; that is -(Bdivxprp+Bkxprp) up to the equilibrium force-balance residual. See :doc:`pmodb_gauge` before using these columns for NTV post processing. -GPEC_xbnormal_n#.out +gpec_xbnormal_n#.out **Flag** xbnormal_flag - **Info** The normal components of the displacement, magnetic field without the plasma response, and magnetic field with the plasma response included for each poloidal mode number at each value of \Psi_{N} output.?? + **Info** The normal displacement (xno), the normal perturbed magnetic field (bno), and the Jacobian-weighted contravariant psi component of the perturbed field J\mathbf{B}\cdot\nabla\psi (bwp) for each poloidal mode number at each value of \Psi_{N} output. -GPEC_*rzphi_n#.out +gpec_*rzphi_n#.out -A number of output files have a similar structure. Here the * in the file name is replaced by the appropriate leading letters of the corresponding flag. For example the xrzphi_flag for n=1 creates a GPEC_xrzphi_n1.out file. some common properties of these files are: +A number of output files have a similar structure. Here the * in the file name is replaced by the appropriate leading letters of the corresponding flag. For example the xrzphi_flag for n=1 creates a gpec_xrzphi_n1.out file. some common properties of these files are: • real and imaginary components: Output files contain two dimensional data on an \left(r,z\right) grid for a single toroidal harmonic. To translate into three dimensions, perform the transformationB\left(r,z,\phi\right)=B_{real}\left(r,z\right)\cos\left(n\phi\right)+B_{imag}\left(r,z\right)\sin\left(n\phi\right) @@ -211,7 +211,7 @@ A number of output files have a similar structure. Here the * in the file name i **Info** The false perturbed magnetic field in the vacuum region on the \left(r,z,\phi\right) grid calculated using the GPEC boundary surface current composed of both the vacuum component and the plasma response. - **Flag** vpbrzphi_flag + **Flag** pbrzphi_flag (writes gpec_vpbrzphi_n#.out) **Info** The true perturbed magnetic field in the vacuum region on the \left(r,z,\phi\right) grid due to the plasma response alone calculated from the plasma response surface condition. @@ -219,9 +219,9 @@ A number of output files have a similar structure. Here the * in the file name i **Info** The false perturbed magnetic field in the vacuum region on the \left(r,z,\phi\right) grid calculated using the GPEC boundary surface current from the external fields alone. - **Flag** ssbrzphi_flag + **Flag** vsbrzphi_flag (with ss_flag(i) selecting surfaces; requires singfld_flag) - **Info** The false perturbed magnetic field in the vacuum region on the \left(r,z,\phi\right) grid calculated using the GPEC boundary surface current from the external fields alone. + **Info** The vacuum field on the \left(r,z,\phi\right) grid due to the currents on the selected resonant surfaces. Energy Decomposition Outputs ---------------------------- @@ -281,7 +281,7 @@ pflux_re(im)_2d.bin bnormal_spectrum.bin - **Flag** bin_flag + **Flag** bin_2d_flag **Info** Surfmn type contours of the normal perturbed magnetic fields as a function of poloidal harmonic number and \Psi_{N}. diff --git a/docs/sign_conventions.rst b/docs/sign_conventions.rst index 40dfccb45..25290474a 100644 --- a/docs/sign_conventions.rst +++ b/docs/sign_conventions.rst @@ -57,8 +57,8 @@ Toroidal Coordinate :math:`\zeta` and :math:`\phi` Working Coordinate Options -------------------------- -Controlled by ``jac_type`` in ``equil.in``. The Jacobian is -:math:`J \propto B_p^{p_{bp}} \, B^{p_b} \, R^{-p_r}`: +Controlled by ``jac_type`` in ``equil.in``. The Jacobian is set by +:math:`\mathbf{B}\cdot\nabla\theta \propto 1/J \propto B_p^{p_{bp}} \, B^{p_b} \, R^{-p_r}`: .. list-table:: :header-rows: 1 @@ -84,6 +84,14 @@ Controlled by ``jac_type`` in ``equil.in``. The Jacobian is - 1 - 0 - 0 + * - Park + - 0 + - 1 + - 0 + * - Other (powers set explicitly in ``equil.in``) + - + - + - See also the coordinate discussion in :doc:`outputs`. diff --git a/docs/source_documentation.rst b/docs/source_documentation.rst index 038d69188..5c3767513 100644 --- a/docs/source_documentation.rst +++ b/docs/source_documentation.rst @@ -1,6 +1,6 @@ .. _source_documentation: -The recommended post processing for GPEC is now to use the devoted `OMFIT module `. +The recommended post processing for GPEC is now to use the devoted `OMFIT module `_. The deprecated code below has been ported there and improved, but is included here for reference. diff --git a/input/dcon.in b/input/dcon.in index 8db588137..57ec14a12 100644 --- a/input/dcon.in +++ b/input/dcon.in @@ -52,7 +52,7 @@ 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 + mthsurf0=1 ! Linear multiplier on number of boundary points used to display surface eigenfunctions for ahb_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 diff --git a/input/equil.in b/input/equil.in index 86f3b7f3c..6189d78c0 100644 --- a/input/equil.in +++ b/input/equil.in @@ -2,7 +2,7 @@ eq_type="efit" ! Type of the input 2D equilibrium file. Accepts efit, chease, fluxgrid, transp, jsolver, lar, sol, etc. eq_filename="mypath" ! Path to input file (geqdsk, etc.) Must specified (no default). - jac_type="hamada" ! Working coordinate system for all DCON and GPEC calculations. Overrides individual powers. Accepts hamada, pest, boozer, equal_arc + jac_type="hamada" ! Working coordinate system for all DCON and GPEC calculations. Overrides individual powers. Accepts hamada, pest, boozer, equal_arc, park, other 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 diff --git a/install/README b/install/README index 992eaeb01..3b3dd0942 100644 --- a/install/README +++ b/install/README @@ -45,8 +45,8 @@ Personal copies of GPEC should be cloned from the `Github repository &1 | tee make.log + make 2>&1 | tee make.log where the logging part is merely a useful suggestion. Appropriate compiler details for the major fusion clusters on which public installations are maintained is recorded in the $GPECHOME/install/DEFAULTS.inc file annotations. \ No newline at end of file From 459114934f1f17ac8c6ce220b739fa2a000e67d9 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 24 Oct 2026 05:38:29 -0400 Subject: [PATCH 8/8] DOCS - CLEANUP - Drop the pmodb derivation scripts from the release The symbolic proof scripts under docs/derivations are a new verification framework that nothing else in the repo uses, so they are out of scope for 1.6.0. The pmodb_gauge page keeps its statements and equations; only the references to the scripts and their identity labels are removed. The scripts remain in this PR's commit history. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VMDGyfQKRZQMdk3NJXsKQE --- docs/derivations/pmodb_gauge/derive.py | 123 ---- .../pmodb_gauge/fortsym_verdicts.json | 373 ---------- docs/derivations/pmodb_gauge/identities.py | 637 ------------------ docs/pmodb_gauge.rst | 77 +-- 4 files changed, 29 insertions(+), 1181 deletions(-) delete mode 100644 docs/derivations/pmodb_gauge/derive.py delete mode 100644 docs/derivations/pmodb_gauge/fortsym_verdicts.json delete mode 100644 docs/derivations/pmodb_gauge/identities.py diff --git a/docs/derivations/pmodb_gauge/derive.py b/docs/derivations/pmodb_gauge/derive.py deleted file mode 100644 index abf835099..000000000 --- a/docs/derivations/pmodb_gauge/derive.py +++ /dev/null @@ -1,123 +0,0 @@ -"""Run the pmodb gauge identities through a CAS and record the verdicts. - -Usage:: - - python3 docs/derivations/pmodb_gauge/derive.py - -FortSym (a checkout with a built C ABI library on ``PYTHONPATH``) is used -when it is importable; otherwise plain SymPy is used. Writes -``fortsym_verdicts.json`` beside this file, with the engine recorded. -``PROVED`` means the zero test returned true; ``UNKNOWN`` and ``NONZERO`` -fail an identity and are the required outcome for a negative control. No -numerical probe is accepted as evidence. -""" - -from __future__ import annotations - -import json -import signal -import subprocess -import sys -import time -from pathlib import Path - -HERE = Path(__file__).resolve().parent -sys.path.insert(0, str(HERE)) - -try: - import fortsym # noqa: E402 - import fortsym.sympy as fs # noqa: E402 -except ImportError: - fortsym = None - import sympy as fs # noqa: E402 -import identities # noqa: E402 - -TIMEOUT = 600 # seconds per SymPy zero test - - -class _Timeout(Exception): - pass - - -def _alarm(*_): - raise _Timeout - - -def verdict(expr) -> str: - if fortsym is not None: - value = fs.sympify(expr).is_zero - else: - signal.signal(signal.SIGALRM, _alarm) - signal.alarm(TIMEOUT) - try: - e = fs.expand(expr) - value = True if e == 0 else fs.simplify(fs.together(e)) == 0 or None - except _Timeout: - return "UNKNOWN" - finally: - signal.alarm(0) - return {True: "PROVED", False: "NONZERO"}.get(value, "UNKNOWN") - - -def counterexample(expr): - """Exact rational point at which the engine decides the residual is non-zero.""" - symbols = sorted(expr.free_symbols, key=str) - for attempt in range(1, 6): - point, value = {}, expr - for j, sym in enumerate(symbols): - num = 7 * attempt + 3 * j + 1 - point[str(sym)] = f"{num}/{5 + j}" - value = value.subs(sym, fs.Rational(num, 5 + j)) - if value.is_zero is False: - return point - return None - - -def revision() -> str: - if fortsym is None: - return "sympy " + fs.__version__ - root = Path(fortsym.__file__).resolve().parents[2] - try: - return subprocess.run(["git", "-C", str(root), "rev-parse", "HEAD"], - capture_output=True, text=True, check=True).stdout.strip() - except (OSError, subprocess.CalledProcessError): - return "unknown" - - -def main() -> int: - records, failed = [], [] - for label, equation, residual in identities.all_identities(fs): - start = time.perf_counter() - v = verdict(residual) - records.append({"label": label, "equation": equation, "verdict": v, - "seconds": round(time.perf_counter() - start, 3)}) - if v != "PROVED": - failed.append(label) - print(f"{label:6s} {v:8s} {equation}", flush=True) - for label, equation, residual in identities.open_gaps(fs): - v = verdict(residual) - records.append({"label": label, "equation": equation, "verdict": v, - "open_gap": True}) - print(f"{label:6s} {v:8s} (open FortSym gap) {equation}", flush=True) - for label, equation, residual in identities.negative_controls(fs): - v = verdict(residual) - records.append({"label": label, "equation": equation, "verdict": v, - "negative_control": True}) - witness = counterexample(residual) if v != "NONZERO" else None - if witness is not None: - records[-1]["counterexample"] = witness - if v == "PROVED" or (v != "NONZERO" and witness is None): - failed.append(label) - print(f"{label:6s} {v:8s} (negative control" - f"{', witness ' + str(witness) if witness else ''}) {equation}", flush=True) - engine = ("fortsym native zero test via C ABI (fortsym.sympy)" if fortsym - else "sympy expand/simplify zero test") - out = {"engine": engine, - "fortsym_revision": revision(), "checks": records, "failed": failed} - (HERE / "fortsym_verdicts.json").write_text(json.dumps(out, indent=2) + "\n") - print("FAILED: " + " ".join(failed) if failed else "required identities proved; see explicit open-gap records") - return 1 if failed else 0 - - -if __name__ == "__main__": - raise SystemExit(main()) diff --git a/docs/derivations/pmodb_gauge/fortsym_verdicts.json b/docs/derivations/pmodb_gauge/fortsym_verdicts.json deleted file mode 100644 index 6c74f984d..000000000 --- a/docs/derivations/pmodb_gauge/fortsym_verdicts.json +++ /dev/null @@ -1,373 +0,0 @@ -{ - "engine": "sympy expand/simplify zero test", - "fortsym_revision": "sympy 1.9", - "checks": [ - { - "label": "G1.1", - "equation": "J dB^psi = GPEC bwp (pointwise)", - "verdict": "PROVED", - "seconds": 0.002 - }, - { - "label": "G1.2", - "equation": "J dB^theta = GPEC bwt (pointwise)", - "verdict": "PROVED", - "seconds": 0.003 - }, - { - "label": "G1.3", - "equation": "J dB^zeta = GPEC bwz (pointwise)", - "verdict": "PROVED", - "seconds": 0.006 - }, - { - "label": "G1.4", - "equation": "(xi x B)_psi = chi1 (q xi^theta - xi^zeta) = xi_s", - "verdict": "PROVED", - "seconds": 0.001 - }, - { - "label": "G1.5", - "equation": "xi_s(xi + lam B) = xi_s(xi)", - "verdict": "PROVED", - "seconds": 0.001 - }, - { - "label": "G1F.1", - "equation": "bwp = chi1 (m - n q) 2 pi i xsp", - "verdict": "PROVED", - "seconds": 0.002 - }, - { - "label": "G1F.2", - "equation": "bwt = -(chi1 xsp1 + 2 pi i n xss)", - "verdict": "PROVED", - "seconds": 0.0 - }, - { - "label": "G1F.3", - "equation": "bwz = -(chi1 (q1 xsp + q xsp1) + 2 pi i m xss)", - "verdict": "PROVED", - "seconds": 0.0 - }, - { - "label": "G2.1", - "equation": "(d_t + q d_z)(J xi^t) = d_z(J xi_s)/chi1 - d_psi(J xi^psi) + J div xi", - "verdict": "PROVED", - "seconds": 0.001 - }, - { - "label": "G2.2", - "equation": "(d_t + q d_z)(J xi^z) = -q d_psi(J xi^psi) - d_t(J xi_s)/chi1 + q J div xi", - "verdict": "PROVED", - "seconds": 0.005 - }, - { - "label": "G2.3", - "equation": "q xwt - xwz = J xi_s / chi1", - "verdict": "PROVED", - "seconds": 0.001 - }, - { - "label": "G3", - "equation": "B.dB = chi1 (bvt + q bvz)/J (eulb times |B|)", - "verdict": "PROVED", - "seconds": 0.052 - }, - { - "label": "G4.1", - "equation": "J xi_perp^theta = chi1^2 (-f2 xsp + f3 xms/chi1) (cspl 1)", - "verdict": "PROVED", - "seconds": 0.469 - }, - { - "label": "G4.2", - "equation": "J xi_perp^zeta = q J xi_perp^theta - J xms/chi1 (cspl 2)", - "verdict": "PROVED", - "seconds": 0.307 - }, - { - "label": "G4.3", - "equation": "B_i xi_perp^i = 0 for xi_perp^i = xi^i - (B_j xi^j / B_k B^k) B^i", - "verdict": "PROVED", - "seconds": 0.055 - }, - { - "label": "G4.4", - "equation": "xsp d_psi B + (cspl 1/J) d_theta B = xi_perp.grad B (curvfun)", - "verdict": "PROVED", - "seconds": 0.325 - }, - { - "label": "G5.1", - "equation": "xi_g^zeta = 0 for lam = -xi_perp^zeta / B^zeta", - "verdict": "PROVED", - "seconds": 0.0 - }, - { - "label": "G5.2", - "equation": "xi_g^theta = xi_s/(chi1 q)", - "verdict": "PROVED", - "seconds": 0.032 - }, - { - "label": "G5.2b", - "equation": "xi_g^psi = xi^psi", - "verdict": "PROVED", - "seconds": 0.0 - }, - { - "label": "G5.3", - "equation": "curl(xi_g x B) = curl(xi_perp x B) (theta component)", - "verdict": "PROVED", - "seconds": 0.037 - }, - { - "label": "G5.4", - "equation": "lagb (gpout.f) = dB_L(xi_g), the xi^zeta = 0 gauge", - "verdict": "PROVED", - "seconds": 0.005 - }, - { - "label": "G5.5", - "equation": "lagb - dB_L(xi_perp) = -(xi_perp^zeta/q) d_theta B", - "verdict": "PROVED", - "seconds": 0.016 - }, - { - "label": "G5.6", - "equation": "llagb (gpdiag.f) = dB_L(xi_perp)", - "verdict": "PROVED", - "seconds": 0.308 - }, - { - "label": "G5.7", - "equation": "dB_L(xi, div-free) = lagb + (xmz/(J q)) d_theta B", - "verdict": "PROVED", - "seconds": 0.002 - }, - { - "label": "G5.8", - "equation": "xi_perp^zeta = -(B_theta xi_s/chi1 + q B_psi xi^psi)/(B_theta + q B_zeta)", - "verdict": "PROVED", - "seconds": 0.362 - }, - { - "label": "G7.1", - "equation": "J B^2 (xi_perp.grad B + mu0 p' xi^psi/B)/B = sband xsp + tband xms/chi1, xi = e_psi", - "verdict": "PROVED", - "seconds": 0.013 - }, - { - "label": "G7.2", - "equation": "J B^2 (xi_perp.grad B + mu0 p' xi^psi/B)/B = sband xsp + tband xms/chi1, xi = e_theta", - "verdict": "PROVED", - "seconds": 0.013 - }, - { - "label": "G7.3", - "equation": "J B^2 (xi_perp.grad B + mu0 p' xi^psi/B)/B = sband xsp + tband xms/chi1, xi = e_zeta", - "verdict": "PROVED", - "seconds": 0.009 - }, - { - "label": "C1", - "equation": "b.curl(xi_perp x B) + xi_perp.grad B = -B(div xi_perp + kappa.xi_perp)", - "verdict": "PROVED", - "seconds": 0.002 - }, - { - "label": "C2", - "equation": "B kappa.xi_perp = xi_perp.grad B + xi_perp.((curl B) x B)/B", - "verdict": "PROVED", - "seconds": 0.004 - }, - { - "label": "R2", - "equation": "J div xi (R = 1) = 0", - "verdict": "PROVED", - "seconds": 0.028 - }, - { - "label": "R3", - "equation": "J div dB (reg) = 2 pi i chi1 (m - n q)(1 - R) xsp1", - "verdict": "PROVED", - "seconds": 0.005 - }, - { - "label": "R6", - "equation": "H_parallel = B(X1-U)+(B_theta Delta_zeta-B_zeta Delta_theta)/(J B)", - "verdict": "PROVED", - "seconds": 0.01 - }, - { - "label": "R7", - "equation": "lagb+Bdivxprp+Bkxprp = -xi_rec_perp^zeta Bt/q + 2 Delta Bt/(chi1 q) + curl_delta + B(V-U) - xi_rec_perp.F/B", - "verdict": "PROVED", - "seconds": 0.002 - }, - { - "label": "R8", - "equation": "with V = U and Delta = 0: lagb+Bdivxprp+Bkxprp = -xi_perp^zeta Bt/q - xi_perp.F/B", - "verdict": "PROVED", - "seconds": 0.003 - }, - { - "label": "E4.1", - "equation": "d/de J(X + e xi) = d_i(J xi^i) = J div xi", - "verdict": "PROVED", - "seconds": 0.0 - }, - { - "label": "G8.1", - "equation": "J B**2 div xi_perp = xmat Xp + ymat X + zmat A", - "verdict": "PROVED", - "seconds": 0.008 - }, - { - "label": "G8.2", - "equation": "idcon zband2/zband3 normalization from B**2=chi1**2*(g22+2qg23+q**2g33)/J**2", - "verdict": "PROVED", - "seconds": 0.006 - }, - { - "label": "Q1.hamada", - "equation": "T_a(xi+lambda B) = T_a(xi)+B.grad(lambda B**a)/B**a", - "verdict": "PROVED", - "seconds": 0.0 - }, - { - "label": "Q2.hamada", - "equation": "D w = -T; w_m = -T_m/[2 pi i(m-nq)]", - "verdict": "PROVED", - "seconds": 0.0 - }, - { - "label": "Q3.hamada", - "equation": "delta B_target = delta B_X + w d_theta B (J B**a flux function)", - "verdict": "PROVED", - "seconds": 0.0 - }, - { - "label": "Q1.boozer", - "equation": "T_a(xi+lambda B) = T_a(xi)+B.grad(lambda B**a)/B**a", - "verdict": "PROVED", - "seconds": 0.006 - }, - { - "label": "Q2.boozer", - "equation": "D w = -T; w_m = -T_m/[2 pi i(m-nq)]", - "verdict": "PROVED", - "seconds": 0.0 - }, - { - "label": "Q3.boozer", - "equation": "delta B_target = delta B_X + w d_theta B (J B**a flux function)", - "verdict": "PROVED", - "seconds": 0.0 - }, - { - "label": "Q4", - "equation": "det d(theta_H+iota h,phi_H+h)/d(theta_H,phi_H)=1+D h", - "verdict": "PROVED", - "seconds": 0.001 - }, - { - "label": "Q5", - "equation": "Hamada->Boozer equals direct X->Boozer when w_HB=w_XB-w_XH", - "verdict": "PROVED", - "seconds": 0.0 - }, - { - "label": "Q6", - "equation": "field-aligned input: delta B_E=0, xi_perp=0 => target delta B=0", - "verdict": "PROVED", - "seconds": 0.0 - }, - { - "label": "Q7", - "equation": "xi_perp limit: zero target source implies no gauge correction", - "verdict": "PROVED", - "seconds": 0.0 - }, - { - "label": "Q8", - "equation": "delta(J/R**2)/(J/R**2)=div xi-2 xi_R/R", - "verdict": "PROVED", - "seconds": 0.002 - }, - { - "label": "K5", - "equation": "Complex CGL pressure work equals [5/3 |D|**2+1/3 |D+3K|**2]/2", - "verdict": "PROVED", - "seconds": 0.007 - }, - { - "label": "P1", - "equation": "(Bk - (llagb-eulb))/B = mu0 p' xi^psi/B**2", - "verdict": "PROVED", - "seconds": 0.0 - }, - { - "label": "K1", - "equation": "CGL energy density = (1/2)[(5/3) D^2 + (1/3)(D + 3K)^2]", - "verdict": "PROVED", - "seconds": 0.004 - }, - { - "label": "K3", - "equation": "v [drive change under xi -> xi + lam B - d/dl (v B lam)] = 0", - "verdict": "PROVED", - "seconds": 0.008 - }, - { - "label": "G6", - "equation": "B.dB + xi_perp.grad B^2 + B^2 div xi_perp + xi_perp.(curl B x B) = 0", - "verdict": "UNKNOWN", - "open_gap": true - }, - { - "label": "N1", - "equation": "lagb is the xi_perp-only Lagrangian |B| (false)", - "verdict": "UNKNOWN", - "negative_control": true, - "counterexample": { - "Bm_theta": "8/5", - "J_0": "11/6", - "chi1": "14/7", - "g_11_0": "17/8", - "g_12_0": "20/9", - "g_13_0": "23/10", - "g_22_0": "26/11", - "g_23_0": "29/12", - "q_0": "32/13", - "xp_0": "35/14", - "xt_0": "38/15", - "xz_0": "41/16" - } - }, - { - "label": "N2", - "equation": "lagb is Park 2009 Eq. 7 with the div-free xi (false)", - "verdict": "UNKNOWN", - "negative_control": true, - "counterexample": { - "Bm_theta": "8/5", - "q_0": "11/6", - "xt_0": "14/7", - "xz_0": "17/8" - } - }, - { - "label": "K2neg", - "equation": "fcgl with kapx = -(dbob + divx)/2 equals the CGL density", - "verdict": "UNKNOWN", - "negative_control": true, - "counterexample": { - "D": "8/5", - "K": "11/6" - } - } - ], - "failed": [] -} diff --git a/docs/derivations/pmodb_gauge/identities.py b/docs/derivations/pmodb_gauge/identities.py deleted file mode 100644 index e8a31a9b3..000000000 --- a/docs/derivations/pmodb_gauge/identities.py +++ /dev/null @@ -1,637 +0,0 @@ -"""Symbolic identities behind the GPEC pmodb |B| columns. - -Every function returns ``(label, equation, residual)`` triples; a residual must -vanish identically. The module is written against a SymPy-shaped namespace -``sp`` so that the same derivation runs on FortSym (``fortsym.sympy``, the -authority used by ``derive.py``). Labels refer to ``docs/pmodb_gauge.rst``. - -GPEC/DCON coordinates (psi, theta, zeta): angles in turns, general -straight-field-line Jacobian ``J(psi, theta)`` (any ``jac_type``), - - B = chi1 grad(zeta - q theta) x grad(psi), J B^theta = chi1, - J B^zeta = chi1 q, B^psi = 0, - -and the surface displacement (Glasser 2016 Eqs. 9 and 18) - - xi_s = (xi x B)_psi = chi1 (q xi^theta - xi^zeta). - -Pointwise identities are proved on jets: every field is an affine (1-jet) -polynomial in (psi, theta, zeta) about the origin with independent symbolic -coefficients, and the residual is evaluated at the origin after -differentiation. Each identity involves at most first derivatives of the -jetted fields, so it holds for all smooth fields iff it holds for all 1-jets. -The metric is general except for ``det g = J**2``, imposed by solving for -``g_33``. -""" - -from __future__ import annotations - -COORDS = ("psi", "theta", "zeta") - - -class Dual: - """First-order jet at the origin: value ``v`` and derivatives ``d[i]``. - - ``d is None`` marks a value whose own derivatives were not carried (a - derivative of a jet); differentiating it again raises. - """ - - def __init__(self, v, d=(0, 0, 0)): - self.v, self.d = v, (None if d is None else tuple(d)) - - @staticmethod - def lift(x): - return x if isinstance(x, Dual) else Dual(x) - - def _need(self): - if self.d is None: - raise ValueError("second derivative requested from a 1-jet") - return self.d - - def __add__(self, o): - o = Dual.lift(o) - if self.d is None or o.d is None: - return Dual(self.v + o.v, None) - return Dual(self.v + o.v, [a + b for a, b in zip(self.d, o.d)]) - - __radd__ = __add__ - - def __neg__(self): - return Dual(-self.v, None if self.d is None else [-a for a in self.d]) - - def __sub__(self, o): - return self + (-Dual.lift(o)) - - def __rsub__(self, o): - return Dual.lift(o) - self - - def __mul__(self, o): - o = Dual.lift(o) - if self.d is None or o.d is None: - return Dual(self.v * o.v, None) - return Dual(self.v * o.v, [self.v * b + a * o.v for a, b in zip(self.d, o.d)]) - - __rmul__ = __mul__ - - def inv(self): - if self.d is None: - return Dual(1 / self.v, None) - return Dual(1 / self.v, [-a / self.v ** 2 for a in self.d]) - - def __truediv__(self, o): - return self * Dual.lift(o).inv() - - def __rtruediv__(self, o): - return Dual.lift(o) * self.inv() - - def __pow__(self, k): - assert isinstance(k, int) and k >= 1 - out = self - for _ in range(k - 1): - out = out * self - return out - - -class Setup: - """Jetted GPEC fields (1-jets at the origin) in a SymPy-shaped namespace.""" - - def __init__(self, sp, metric=False): - self.sp = sp - self.u = (0, 1, 2) # psi, theta, zeta - self.chi1 = Dual(sp.Symbol("chi1", positive=True)) - self.q = self.jet("q", (0,)) - self.J = self.jet("J", (0, 1)) - # Contravariant displacement components (general, incl. xi_parallel). - self.xi = [self.jet(f"x{c[0]}", self.u) for c in COORDS] - self.Bup = [Dual(0), self.chi1 / self.J, self.chi1 * self.q / self.J] - if metric: - g = {} - for (i, j) in ((0, 0), (0, 1), (0, 2), (1, 1), (1, 2)): - g[(i, j)] = self.jet(f"g_{i + 1}{j + 1}", (0, 1)) - minor = g[(0, 0)] * g[(1, 1)] - g[(0, 1)] ** 2 - rest = (-g[(0, 0)] * g[(1, 2)] ** 2 + 2 * g[(0, 1)] * g[(1, 2)] * g[(0, 2)] - - g[(1, 1)] * g[(0, 2)] ** 2) - g[(2, 2)] = (self.J ** 2 - rest) / minor - self.g = [[g[(min(i, j), max(i, j))] for j in range(3)] for i in range(3)] - - def jet(self, name, deps): - sp = self.sp - return Dual(sp.Symbol(f"{name}_0"), - [sp.Symbol(f"{name}_{COORDS[i]}") if i in deps else 0 for i in range(3)]) - - def d(self, f, i): - f = Dual.lift(f) - return Dual(f._need()[i], None) - - def at0(self, e): - return Dual.lift(e).v - - # --- tensor algebra in the (psi, theta, zeta) chart ----------------------- - @staticmethod - def eps(i, j, k): - return (i - j) * (j - k) * (k - i) // 2 - - def cross_cov(self, a, b): - """(a x b)_i = J eps_ijk a^j b^k for contravariant a, b.""" - return [self.J * sum(self.eps(i, j, k) * a[j] * b[k] - for j in range(3) for k in range(3)) for i in range(3)] - - def curl_up(self, A): - """(curl A)^i = eps^ijk d_j A_k / J for covariant A.""" - return [sum(self.eps(i, j, k) * self.d(A[k], j) - for j in range(3) for k in range(3)) / self.J for i in range(3)] - - def div(self, v): - return sum(self.d(self.J * v[i], i) for i in range(3)) / self.J - - def lower(self, v): - return [sum(self.g[i][j] * v[j] for j in range(3)) for i in range(3)] - - def dot(self, a, b): - return sum(a[i] * b[i] for i in range(3)) - - def xi_s(self, v): - return self.chi1 * (self.q * v[1] - v[2]) - - -def induction(sp): - """G1: GPEC's J dB^i (gpeq.f:114-116) is J curl(xi x B)^i for any Jacobian.""" - S = Setup(sp) - xi, B = S.xi, S.Bup - JdB = [S.J * c for c in S.curl_up(S.cross_cov(xi, B))] - xs = S.xi_s(xi) - d = S.d - gpec = [S.chi1 * (d(xi[0], 1) + S.q * d(xi[0], 2)), - d(xs, 2) - S.chi1 * d(xi[0], 0), - -d(xs, 1) - S.chi1 * d(S.q * xi[0], 0)] - out = [(f"G1.{i + 1}", f"J dB^{COORDS[i]} = GPEC bw{'ptz'[i]} (pointwise)", - S.at0(JdB[i] - gpec[i])) for i in range(3)] - # (xi x B)_psi = xi_s: the surface displacement of Glasser (2016) Eq. 9. - out.append(("G1.4", "(xi x B)_psi = chi1 (q xi^theta - xi^zeta) = xi_s", - S.at0(S.cross_cov(xi, B)[0] - xs))) - # A field-aligned displacement changes neither xi^psi nor xi_s. - lam = S.jet("lam", S.u) - shifted = [xi[i] + lam * B[i] for i in range(3)] - out.append(("G1.5", "xi_s(xi + lam B) = xi_s(xi)", - S.at0(S.xi_s(shifted) - xs))) - return out - - -def fourier_symbols(sp): - """G1F: Fourier symbols of G1 for xi ~ exp(2 pi i (m theta - n zeta)).""" - I, pi = sp.I, sp.pi - m, n, chi1 = sp.Symbol("m"), sp.Symbol("n"), sp.Symbol("chi1", positive=True) - q, q1 = sp.Symbol("q"), sp.Symbol("q1") - X, X1, Xs = sp.Symbol("xsp"), sp.Symbol("xsp1"), sp.Symbol("xss") - dth, dze = 2 * pi * I * m, -2 * pi * I * n - out = [] - out.append(("G1F.1", "bwp = chi1 (m - n q) 2 pi i xsp", - chi1 * (dth + q * dze) * X - chi1 * (m - n * q) * 2 * pi * I * X)) - out.append(("G1F.2", "bwt = -(chi1 xsp1 + 2 pi i n xss)", - (dze * Xs - chi1 * X1) + (chi1 * X1 + 2 * pi * I * n * Xs))) - out.append(("G1F.3", "bwz = -(chi1 (q1 xsp + q xsp1) + 2 pi i m xss)", - (-dth * Xs - chi1 * (q1 * X + q * X1)) - + (chi1 * (q1 * X + q * X1) + 2 * pi * I * m * Xs))) - return out - - -def incompressible_closure(sp): - """G2: GPEC's xi^theta, xi^zeta (gpeq.f gpeq_contra) follow from div xi = 0.""" - S = Setup(sp) - xi, J, q, chi1, d = S.xi, S.J, S.q, S.chi1, S.d - D = sum(d(J * xi[i], i) for i in range(3)) # J div xi - xs = S.xi_s(xi) - Jt, Jz = J * xi[1], J * xi[2] - r1 = (d(Jt, 1) + q * d(Jt, 2)) - d(J * xs, 2) / chi1 + d(J * xi[0], 0) - D - r2 = (d(Jz, 1) + q * d(Jz, 2)) + q * d(J * xi[0], 0) + d(J * xs, 1) / chi1 - q * D - r3 = q * Jt - Jz - J * xs / chi1 - return [("G2.1", "(d_t + q d_z)(J xi^t) = d_z(J xi_s)/chi1 - d_psi(J xi^psi) + J div xi", - S.at0(r1)), - ("G2.2", "(d_t + q d_z)(J xi^z) = -q d_psi(J xi^psi) - d_t(J xi_s)/chi1" - " + q J div xi", S.at0(r2)), - ("G2.3", "q xwt - xwz = J xi_s / chi1", S.at0(r3))] - - -def columns(sp): - """G3-G5: eulb, xi_perp components, the gauge of lagb, and llagb.""" - S = Setup(sp, metric=True) - xi, B, J, q, chi1 = S.xi, S.Bup, S.J, S.q, S.chi1 - Bm = S.jet("Bm", S.u[:2]) # |B|(psi, theta), axisymmetric - dB = S.curl_up(S.cross_cov(xi, B)) - dBcov = S.lower(dB) - B2 = S.dot(S.lower(B), B) - xs = S.xi_s(xi) - g = S.g - xmz = J * xi[2] # GPEC xmz (J xi^zeta) - sigma = S.dot(S.lower(xi), B) / B2 - xperp = [xi[i] - sigma * B[i] for i in range(3)] - # gpout_pmodb/gpdiag_pmodb cspl: J xi_perp from xi^psi and xi_s with the - # eqfun coefficients f2, f3 (equil/direct.f, covariant metric). - f2 = (g[0][1] + q * g[0][2]) / (J * B2) - f3 = (g[1][2] + q * g[2][2]) / (J * B2) - jxt = chi1 ** 2 * (-f2 * xi[0] + f3 * xs / chi1) - jxz = q * jxt - J * xs / chi1 - out = [] - out.append(("G3", "B.dB = chi1 (bvt + q bvz)/J (eulb times |B|)", - S.at0(S.dot(dBcov, B) - chi1 * (dBcov[1] + q * dBcov[2]) / J))) - out.append(("G4.1", "J xi_perp^theta = chi1^2 (-f2 xsp + f3 xms/chi1) (cspl 1)", - S.at0(J * xperp[1] - jxt))) - out.append(("G4.2", "J xi_perp^zeta = q J xi_perp^theta - J xms/chi1 (cspl 2)", - S.at0(J * xperp[2] - jxz))) - # B . xi_perp = 0 as a statement about arbitrary vectors: the metric-jet - # form is larger than necessary, so the - # projection is proved on independent covariant B_i, contravariant B^i, xi^i. - Bl = [sp.Symbol(f"Bl_{i}") for i in range(3)] - Bu = [sp.Symbol(f"Bu_{i}") for i in range(3)] - xu = [sp.Symbol(f"xu_{i}") for i in range(3)] - sig = sum(Bl[i] * xu[i] for i in range(3)) / sum(Bl[i] * Bu[i] for i in range(3)) - out.append(("G4.3", "B_i xi_perp^i = 0 for xi_perp^i = xi^i - (B_j xi^j / B_k B^k) B^i", - sum(Bl[i] * (xu[i] - sig * Bu[i]) for i in range(3)))) - out.append(("G4.4", "xsp d_psi B + (cspl 1/J) d_theta B = xi_perp.grad B (curvfun)", - S.at0(xi[0] * S.d(Bm, 0) + jxt / J * S.d(Bm, 1) - - sum(xperp[i] * S.d(Bm, i) for i in range(3))))) - # Gauge of lagb: xi_g = xi_perp + lam B with xi_g^zeta = 0. - lam = -xperp[2] / B[2] - xg = [xperp[i] + lam * B[i] for i in range(3)] - out.append(("G5.1", "xi_g^zeta = 0 for lam = -xi_perp^zeta / B^zeta", S.at0(xg[2]))) - out.append(("G5.2", "xi_g^theta = xi_s/(chi1 q)", S.at0(xg[1] - xs / (chi1 * q)))) - out.append(("G5.2b", "xi_g^psi = xi^psi", S.at0(xg[0] - xi[0]))) - dBg = S.curl_up(S.cross_cov(xg, B)) - dBp = S.curl_up(S.cross_cov(xperp, B)) - out.append(("G5.3", "curl(xi_g x B) = curl(xi_perp x B) (theta component)", - S.at0(dBg[1] - dBp[1]))) - eul = S.dot(dBcov, B) # |B| * eulb - grad = lambda v: sum(v[i] * S.d(Bm, i) for i in range(3)) - lagb = eul + Bm * (xi[0] * S.d(Bm, 0) + xs / (chi1 * q) * S.d(Bm, 1)) - lagb_g = eul + Bm * grad(xg) - llagb = eul + Bm * grad(xperp) - llagb_code = eul + Bm * (xi[0] * S.d(Bm, 0) + jxt / J * S.d(Bm, 1)) - out.append(("G5.4", "lagb (gpout.f) = dB_L(xi_g), the xi^zeta = 0 gauge", - S.at0(lagb - lagb_g))) - out.append(("G5.5", "lagb - dB_L(xi_perp) = -(xi_perp^zeta/q) d_theta B", - S.at0(lagb - llagb + Bm * xperp[2] / q * S.d(Bm, 1)))) - out.append(("G5.6", "llagb (gpdiag.f) = dB_L(xi_perp)", S.at0(llagb_code - llagb))) - # Park (2009) Eq. 7 with the div-free xi (GPEC's own xmt, xmz). - lagb_div_free = eul + Bm * (xi[0] * S.d(Bm, 0) + xi[1] * S.d(Bm, 1)) - out.append(("G5.7", "dB_L(xi, div-free) = lagb + (xmz/(J q)) d_theta B", - S.at0(lagb_div_free - lagb - Bm * xmz / (J * q) * S.d(Bm, 1)))) - # Boozer-type closed form: B_i xi_perp^i = 0 and the xi_s definition. - Bc = S.lower(B) - out.append(("G5.8", "xi_perp^zeta = -(B_theta xi_s/chi1 + q B_psi xi^psi)" - "/(B_theta + q B_zeta)", - S.at0(xperp[2] + (Bc[1] * xs / chi1 + q * Bc[0] * xi[0]) - / (Bc[1] + q * Bc[2])))) - return out - - -def force_balance_identity(sp): - """G6: |B| times (llagb + Bdivxprp + Bkxprp) = -xi_perp . F, F = (curl B) x B - mu0 grad p. - - With pmodb's Bkxprp = xi_perp.grad B + mu0 p' xi^psi / B (force-balance form - of B kappa.xi_perp, gpout.f curvfun) and Bdivxprp = B div xi_perp, this is - B.dB + xi_perp.grad(B^2) + B^2 div xi_perp + xi_perp.((curl B) x B) = 0. - """ - S = Setup(sp, metric=True) - xi, B = S.xi, S.Bup - Bc = S.lower(B) - B2 = S.dot(Bc, B) - sigma = S.dot(S.lower(xi), B) / B2 - xperp = [xi[i] - sigma * B[i] for i in range(3)] - dB = S.curl_up(S.cross_cov(xi, B)) - curlB = S.curl_up(Bc) - jxb = S.cross_cov(curlB, B) # ((curl B) x B)_i - res = (S.dot(S.lower(dB), B) + sum(xperp[i] * S.d(B2, i) for i in range(3)) - + B2 * S.div(xperp) + S.dot(jxb, xperp)) - return [("G6", "B.dB + xi_perp.grad B^2 + B^2 div xi_perp + xi_perp.(curl B x B) = 0", - S.at0(res))] - - -def pentrc_kappa_matrices(sp): - """G7: idcon.f fmodb sband/tband give J B^2 (Bkxprp/B) from xi^psi and xi^alpha. - - PENTRC set_peq (inputs.f90) and GPEC.jl set_perturbation_data! form - jbb_kapx = smat xsp + tmat xi_alpha with xi_alpha = xms/chi1 (gpout_xclebsch). - The statement is algebraic at a point (only the values of d_psi B^2 and - d_theta B^2 enter), so the metric, Jacobian and derivatives are independent - symbols here: - sband = J (mu0 p' + d_psi B^2/2) - chi1^2 (d_theta B^2/(2 B^2)) (g12 + q g13)/J, - tband = chi1^2 (d_theta B^2/(2 B^2)) (g23 + q g33)/J. - """ - chi1, q, J, P, Dp, Dt = sp.symbols("chi1 q J P1 Dp Dt") - # B^2 = B^i B_i = (chi1/J)(B_theta + q B_zeta) needs no det g = J^2, so the - # metric components are independent here. - g11, g12, g13, g22, g23, g33 = sp.symbols("g11 g12 g13 g22 g23 g33") - g = [[g11, g12, g13], [g12, g22, g23], [g13, g23, g33]] - xi = list(sp.symbols("xp xt xz")) - B = [0, chi1 / J, chi1 * q / J] - Bl = [sum(g[i][j] * B[j] for j in range(3)) for i in range(3)] - B2 = sum(Bl[i] * B[i] for i in range(3)) - sigma = sum(Bl[i] * xi[i] for i in range(3)) / B2 - xperp = [xi[i] - sigma * B[i] for i in range(3)] - xs = chi1 * (q * xi[1] - xi[2]) - lhs = J * (xperp[0] * (P + Dp) + xperp[1] * Dt) - sband = J * (P + Dp) - chi1 ** 2 * (Dt / B2) * (g12 + q * g13) / J - tband = chi1 ** 2 * (Dt / B2) * (g23 + q * g33) / J - res = lhs - sband * xi[0] - tband * xs / chi1 - # The residual is linear in xi; proving it for each basis displacement keeps - # every check independently covers one displacement component. - out = [] - for k, name in enumerate(("psi", "theta", "zeta")): - r = res - for j in range(3): - r = r.subs(xi[j], 1 if j == k else 0) - out.append((f"G7.{k + 1}", "J B^2 (xi_perp.grad B + mu0 p' xi^psi/B)/B" - f" = sband xsp + tband xms/chi1, xi = e_{name}", r)) - return out - - -def cartesian_kappa(sp): - """C1, C2: Cartesian 1-jet forms of the kinematic and curvature identities.""" - x, y, z = sp.symbols("x y z") - X = (x, y, z) - beta = sp.Symbol("beta", positive=True) - g = [[sp.Symbol(f"h_{i}{j}") for j in range(3)] for i in range(3)] - B = [[0, 0, beta][i] + sum(g[i][j] * X[j] for j in range(3)) for i in range(3)] - d = lambda f, c: f.diff(c) if hasattr(f, "diff") else 0 - dot = lambda a, b: sum(a[i] * b[i] for i in range(3)) - cross = lambda a, b: [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], - a[0] * b[1] - a[1] * b[0]] - curl = lambda v: [d(v[2], y) - d(v[1], z), d(v[0], z) - d(v[2], x), - d(v[1], x) - d(v[0], y)] - div = lambda v: sum(d(v[i], X[i]) for i in range(3)) - at0 = lambda e: e.subs(x, 0).subs(y, 0).subs(z, 0) - Bm = sp.sqrt(dot(B, B)) - b = [c / Bm for c in B] - kappa = [sum(b[j] * d(b[i], X[j]) for j in range(3)) for i in range(3)] - a = [sp.Symbol(f"a_{i}") + sum(sp.Symbol(f"a_{i}{j}") * X[j] for j in range(3)) - for i in range(3)] - xp = cross(B, a) - gradB = lambda v: sum(v[j] * d(Bm, X[j]) for j in range(3)) - dBpar = dot(b, curl(cross(xp, B))) - return [("C1", "b.curl(xi_perp x B) + xi_perp.grad B = -B(div xi_perp + kappa.xi_perp)", - at0(dBpar + gradB(xp) + Bm * (div(xp) + dot(kappa, xp)))), - ("C2", "B kappa.xi_perp = xi_perp.grad B + xi_perp.((curl B) x B)/B", - at0(Bm * dot(kappa, xp) - gradB(xp) - dot(xp, cross(curl(B), B)) / Bm))] - - -def regularisation(sp): - """R2-R3: Park-regularised components in GPEC (reg_flag) on one harmonic. - - Diagonal model: one harmonic, J independent of theta. R = s^2/(s^2 + reg^2), - s = m - n q. gpeq_sol uses xmp1 = R xsp1 in bmt/bmz and xms; gpeq_contra - builds xwt/xwz from xmp1. The pmodb columns use xmp1 and xms - (gpout_pmodb), so R2-R3 are the only regularisation identities needed. - """ - I, pi = sp.I, sp.pi - m, n, q, chi1 = sp.symbols("m n q chi1") - J0, J1 = sp.symbols("J_0 J_1") # J and dJ/dpsi (theta-independent) - X, X1, Xs, R = sp.symbols("xsp xsp1 xms R") - s = m - n * q - xwt = -(J0 * R * X1 + J1 * X + 2 * pi * I * n / chi1 * J0 * Xs) / (2 * pi * I * s) - xwz = -(q * J0 * R * X1 + q * J1 * X + 2 * pi * I * m / chi1 * J0 * Xs) / (2 * pi * I * s) - out = [] - jdiv_unreg = (J0 * X1 + J1 * X + 2 * pi * I * m * xwt.subs(R, 1) - - 2 * pi * I * n * xwz.subs(R, 1)) - out.append(("R2", "J div xi (R = 1) = 0", jdiv_unreg)) - # Divergence of the regularised Eulerian field: J div(dB) = d_psi bwp + ... . - bwp1 = 2 * pi * I * chi1 * (s * X1 - n * sp.Symbol("q1") * X) - bmt = -(chi1 * R * X1 + 2 * pi * I * n * Xs) - bmz = -(chi1 * (sp.Symbol("q1") * X + q * R * X1) + 2 * pi * I * m * Xs) - out.append(("R3", "J div dB (reg) = 2 pi i chi1 (m - n q)(1 - R) xsp1", - bwp1 + 2 * pi * I * m * bmt - 2 * pi * I * n * bmz - - 2 * pi * I * chi1 * s * (1 - R) * X1)) - return out - - -def regularized_column_budget(sp): - """General-J continuum budget, allowing independent field/output derivatives. - - U is the derivative selected by gpeq_sol; V is the derivative used in - gpout divergence. Delta=xi_s-chi1(q xi_rec^theta-xi_rec^zeta). - This is bookkeeping under the stated model relations (BDcode, BKcode), - not a proof that the code satisfies them. In gpout_pmodb, V = U - (xmp1) and xi_perp is built from xi_s itself (G4.2), so Delta = 0 and - R8 gives the budget of the output columns. - F denotes xi_rec_perp.((curl B)xB-mu0 grad p)/B, in Tesla. - Metric/field-norm consistency and C1/C2 are assumed; finite truncation - errors are not folded into the physical force residual. - """ - c,J,B,q,Bt,Bt_cov,Delta,Dt,Dz,X1,U,V,D,K,F,Z = sp.symbols( - "chi J B q Bt Btheta Delta Delta_theta Delta_zeta X1 U V D K F Z") - Bz_cov=(J*B**2/c-Bt_cov)/q - # J-weighted induction mismatch: supplied field minus curl(xi_rec x B). - Ht=c*(X1-U)+Dz - Hz=c*q*(X1-U)-Dt - Hpar=(Bt_cov*Ht+Bz_cov*Hz)/(J*B) - curl_delta=(Bt_cov*Dz-Bz_cov*Dt)/(J*B) - out=[("R6", "H_parallel = B(X1-U)+(B_theta Delta_zeta-B_zeta Delta_theta)/(J B)", - Hpar-B*(X1-U)-curl_delta)] - g=Delta*Bt/(c*q) - # C1 for the reconstructed perpendicular displacement, plus the supplied H. - LLrec=-B*(D+K)+Hpar - LLcode=LLrec+g - Lagcode=LLcode-Z*Bt/q - BDcode=B*(D+V-X1) - BKcode=B*K+g-F - budget=-Z*Bt/q+2*g+curl_delta+B*(V-U)-F - out.append(("R7", "lagb+Bdivxprp+Bkxprp = -xi_rec_perp^zeta Bt/q + 2 Delta Bt/(chi1 q) + curl_delta + B(V-U) - xi_rec_perp.F/B", - Lagcode+BDcode+BKcode-budget)) - out.append(("R8", "with V = U and Delta = 0: lagb+Bdivxprp+Bkxprp = -xi_perp^zeta Bt/q - xi_perp.F/B", - (Lagcode+BDcode+BKcode).subs({V: U, Delta: 0, Dt: 0, Dz: 0})+Z*Bt/q+F)) - return out - - -def cgl(sp): - """K1-K2: the PENTRC/GPEC.jl fcgl integrand against the CGL closure. - - D = div xi_perp, K = kappa.xi_perp, dB/B = -(D + K) (C1). CGL: - dp_perp/p = -D + dB/B, dp_par/p = -3D - 2 dB/B (no xi_par). Logan (2013) - Eq. 4 energy density, -(1/2)[(dp_par - dp_perp) dB/B + dp_par D]/p. - """ - D, K = sp.symbols("D K") - b = -(D + K) - dpe, dpa = -D + b, -3 * D - 2 * b - w = -((dpa - dpe) * b + dpa * D) / 2 - code = lambda k: (sp.Rational(5, 3) * D ** 2 + sp.Rational(1, 3) * (D + 3 * k) ** 2) / 2 - return [("K1", "CGL energy density = (1/2)[(5/3) D^2 + (1/3)(D + 3K)^2]", w - code(K)), - ("K2neg", "fcgl with kapx = -(dbob + divx)/2 equals the CGL density", - w - code(-(b + D) / 2))] - - -def drive_gauge(sp): - """K3: gauge invariance of the Logan (2013) Eq. 5 bounce drive. - - Along a field line (arc length l, M = 1): v = sqrt(2(E - mu B)), integrand - (1/v)[(2E - 3 mu B) dB/B + (2E - 2 mu B) div xi]. Under xi -> xi + lam B, - dB -> dB + lam B dB/dl and div xi -> div xi + B dlam/dl (C1, Park 2011 Eq. 4). - """ - l, E, mu = sp.Symbol("l"), sp.Symbol("E", positive=True), sp.Symbol("mu", positive=True) - b0, b1, l0, l1 = sp.symbols("b_0 b_1 lam_0 lam_1") - B = b0 + b1 * l - lam = l0 + l1 * l - v = sp.sqrt(2 * (E - mu * B)) - dBdl, dlam = B.diff(l), lam.diff(l) - change = ((2 * E - 3 * mu * B) * lam * dBdl + (2 * E - 2 * mu * B) * B * dlam) / v - total = (v * B * lam).diff(l) - # Multiplied by v > 0 (away from turning points): FortSym leaves the - # unmultiplied sqrt quotient UNKNOWN. - return [("K3", "v [drive change under xi -> xi + lam B - d/dl (v B lam)] = 0", - ((change - total) * v).subs(l, 0))] - - -def divergence_matrices(sp): - """G8: compare projected contravariant divergence with idcon's x/y/z bands. - - U=(g12+q*g13)/(J*B**2), V=(g23+q*g33)/(J*B**2). The perpendicular - projection gives J xi^theta=chi1**2*(-U X+V A), xi^zeta=q xi^theta-A. - U,V and their theta derivatives remain independent jets here. The field - norm supplies the only algebraic metric relation used by the z bands. - """ - J,Jp,B2,c,q,m,n,X,Xp,A,U,Ut,V,Vt=sp.symbols("J Jp B2 c q m n X Xp A U Ut V Vt") - th,ze=2*sp.pi*sp.I*m,-2*sp.pi*sp.I*n - # Differentiate J xi^theta and J xi^zeta with the product rule. - Y=c**2*(-U*X+V*A) - Yt=c**2*(-Ut*X-U*th*X+Vt*A+V*th*A) - Yz=c**2*(-U*ze*X+V*ze*A) - direct=B2*(J*Xp+Jp*X+Yt+q*Yz-J*ze*A) - bands=J*B2*Xp+(Jp*B2-c**2*B2*Ut-2*sp.pi*sp.I*(m-n*q)*c**2*B2*U)*X - bands+=(c**2*B2*Vt+2*sp.pi*sp.I*(c**2*B2*V*(m-n*q)+n*J*B2))*A - out=[("G8.1","J B**2 div xi_perp = xmat Xp + ymat X + zmat A",direct-bands)] - g22,g23,g33,g12,g13=sp.symbols("g22 g23 g33 g12 g13") - norm=c**2*(g22+2*q*g23+q**2*g33)/J**2 - zcode=c**2/J*(m*(g23+q*g33)+n*(g22+q*g23)) - zproj=(c**2*norm*((g23+q*g33)/(J*norm))*(m-n*q)+n*J*norm) - out.append(("G8.2","idcon zband2/zband3 normalization from B**2=chi1**2*(g22+2qg23+q**2g33)/J**2",zcode-zproj)) - return out - - -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, composed through Hamada. - 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 - - - -def jacobian_identity(sp): - """Eq. (4.1): d/de det(d(X + e xi)/du) = d_i(J xi^i) (Jacobi formula).""" - u = sp.symbols("u_1 u_2 u_3") - eps = sp.Symbol("epsilon") - c = [[sp.Symbol(f"c_{a}{i}") for i in range(3)] for a in range(3)] - q = [[[sp.Symbol(f"q_{a}{min(i, j)}{max(i, j)}") for j in range(3)] - for i in range(3)] for a in range(3)] - pos = [sum(c[a][i] * u[i] for i in range(3)) - + sum(q[a][i][j] * u[i] * u[j] for i in range(3) for j in range(3)) / 2 - for a in range(3)] - xi_up = [sp.Symbol(f"s_{i}") + sum(sp.Symbol(f"t_{i}{j}") * u[j] for j in range(3)) - for i in range(3)] - e_cov = [[pos[a].diff(u[i]) for a in range(3)] for i in range(3)] - xi_cart = [sum(xi_up[i] * e_cov[i][a] for i in range(3)) for a in range(3)] - moved = [pos[a] + eps * xi_cart[a] for a in range(3)] - - def det3(m): - return (m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) - - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) - + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0])) - - def at0(e): - for v in u: - e = e.subs(v, 0) - return e - - # J(X + eps xi) needs no further u-derivative: evaluate entries at u = 0 first. - jac_moved = det3([[at0(moved[a].diff(u[i])) for i in range(3)] for a in range(3)]) - jac = det3([[pos[a].diff(u[i]) for i in range(3)] for a in range(3)]) - first = jac_moved.diff(eps).subs(eps, 0) - divergence_density = at0(sum((jac * xi_up[i]).diff(u[i]) for i in range(3))) - res = first - divergence_density - return [("E4.1", "d/de J(X + e xi) = d_i(J xi^i) = J div xi", res.expand())] - -def all_identities(sp): - return (induction(sp) + fourier_symbols(sp) + incompressible_closure(sp) - + columns(sp) + pentrc_kappa_matrices(sp) - + cartesian_kappa(sp) + regularisation(sp) + regularized_column_budget(sp) - + jacobian_identity(sp) + divergence_matrices(sp) + gauge_coordinates(sp) + [c for c in cgl(sp) if not c[0].endswith("neg")] + drive_gauge(sp)) - - -def open_gaps(sp): - """Identities FortSym does not decide; their content is covered elsewhere. - - G6 is the coordinate form of C1 + C2 (proved on Cartesian 1-jets) with - det g = J^2 imposed through g_33. Neither FortSym nor SymPy decides this - large rational residual; it is recorded, not counted as proved. - """ - return force_balance_identity(sp) - - -def negative_controls(sp): - """Residuals that must NOT vanish (the claims they encode are false).""" - S = Setup(sp, metric=True) - xi, B, J, q, chi1 = S.xi, S.Bup, S.J, S.q, S.chi1 - Bm = S.jet("Bm", S.u[:2]) - B2 = S.dot(S.lower(B), B) - sigma = S.dot(S.lower(xi), B) / B2 - xperp = [xi[i] - sigma * B[i] for i in range(3)] - xs = S.xi_s(xi) - lagb_adv = xi[0] * S.d(Bm, 0) + xs / (chi1 * q) * S.d(Bm, 1) - perp_adv = sum(xperp[i] * S.d(Bm, i) for i in range(3)) - out = [("N1", "lagb is the xi_perp-only Lagrangian |B| (false)", - S.at0(lagb_adv - perp_adv))] - out.append(("N2", "lagb is Park 2009 Eq. 7 with the div-free xi (false)", - S.at0(lagb_adv - (xi[0] * S.d(Bm, 0) + xi[1] * S.d(Bm, 1))))) - out += [c for c in cgl(sp) if c[0].endswith("neg")] - return out diff --git a/docs/pmodb_gauge.rst b/docs/pmodb_gauge.rst index 7ba8fae94..5da0559f9 100644 --- a/docs/pmodb_gauge.rst +++ b/docs/pmodb_gauge.rst @@ -5,10 +5,7 @@ Perturbed field strength (pmodb) and gauge ****************************************** This page states what the columns of ``gpec_pmodb_n#.out`` are, how they are -related, and which of them a field-strength consumer needs. The identities -are proved symbolically by ``docs/derivations/pmodb_gauge/derive.py``; the -verdicts are stored in ``docs/derivations/pmodb_gauge/fortsym_verdicts.json`` -and the labels (G1, G5.5, ...) refer to that file. +related, and which of them a field-strength consumer needs. .. contents:: On this page :local: @@ -35,14 +32,14 @@ any straight-field-line Jacobian :math:`\mathcal{J}(\psi,\theta)` selected by the in-surface displacement across :math:`\mathbf B` (``xms``; ``gpec_xclebsch`` writes :math:`\xi^\alpha=\xi_s/\chi_1`). A field-aligned -displacement changes neither :math:`\xi^\psi` nor :math:`\xi_s` (G1.5), so +displacement changes neither :math:`\xi^\psi` nor :math:`\xi_s`, so :math:`(\xi^\psi,\xi_s)` determine :math:`\boldsymbol\xi_\perp` completely. The parallel part used elsewhere in GPEC is fixed in ``gpeq_contra`` by -:math:`\nabla\cdot\boldsymbol\xi=0` (G2), which gives ``xmt`` +:math:`\nabla\cdot\boldsymbol\xi=0`, which gives ``xmt`` :math:`=\mathcal J\xi^\theta` and ``xmz`` :math:`=\mathcal J\xi^\zeta`. The pmodb columns do not use ``xmt`` or ``xmz``. -For every ``jac_type`` (G1, G1F) the Jacobian-weighted perturbed field of +For every ``jac_type`` the Jacobian-weighted perturbed field of ``gpeq_sol`` is :math:`\mathcal J\,\nabla\times(\boldsymbol\xi\times\mathbf B)`. @@ -51,23 +48,23 @@ The pmodb columns ``eulb`` :math:`\delta B_E=\hat{\mathbf b}\cdot\delta\mathbf B`, the Eulerian - field strength (G3). + field strength. ``lagb`` :math:`\delta B_E+\xi^\psi\partial_\psi B+\xi_s/(\chi_1 q)\,\partial_\theta B`. This is the Lagrangian field strength of the displacement :math:`\boldsymbol\xi_g=\boldsymbol\xi_\perp+\lambda\mathbf B` with - :math:`\xi_g^\zeta=0` in the *working* chart (``jac_type``) (G5.1-G5.4). - It is neither the :math:`\boldsymbol\xi_\perp`-only value (negative - control N1) nor Park et al. (2009) Eq. 7 with the divergence-free - :math:`\boldsymbol\xi` that GPEC constructs (N2; G5.7 gives that value as - ``lagb`` :math:`+\,\xi^\zeta/q\,\partial_\theta B`). ``jac_out`` only + :math:`\xi_g^\zeta=0` in the *working* chart (``jac_type``). + It is neither the :math:`\boldsymbol\xi_\perp`-only value nor Park et + al. (2009) Eq. 7 with the divergence-free :math:`\boldsymbol\xi` that + GPEC constructs (that value is ``lagb`` + :math:`+\,\xi^\zeta/q\,\partial_\theta B`). ``jac_out`` only relabels the harmonics; it does not change this gauge. ``Bdivxprp`` :math:`B\,\nabla\cdot\boldsymbol\xi_\perp`. ``gpout_pmodb`` forms :math:`\mathcal J\boldsymbol\xi_\perp` directly from :math:`\xi^\psi` and - :math:`\xi_s` (G4.1, G4.2), + :math:`\xi_s`, .. math:: @@ -79,28 +76,28 @@ The pmodb columns :math:`f_2=(g_{\psi\theta}+qg_{\psi\zeta})/(\mathcal JB^2)` and :math:`f_3=(g_{\theta\zeta}+qg_{\zeta\zeta})/(\mathcal JB^2)` (``eqfun`` columns 2 and 3, covariant metric). The same projection gives the - ``xmat``-``zmat`` matrices of ``idcon_build`` (G8). + ``xmat``-``zmat`` matrices of ``idcon_build``. ``Bkxprp`` - :math:`\boldsymbol\xi_\perp\cdot\nabla B+\mu_0p'\xi^\psi/B` (G4.4). This is + :math:`\boldsymbol\xi_\perp\cdot\nabla B+\mu_0p'\xi^\psi/B`. This is :math:`B\boldsymbol\kappa\cdot\boldsymbol\xi_\perp` only through equilibrium force balance: exactly, :math:`B\boldsymbol\kappa\cdot\boldsymbol\xi_\perp=\boldsymbol\xi_\perp\cdot\nabla B - +\boldsymbol\xi_\perp\cdot[(\nabla\times\mathbf B)\times\mathbf B]/B` (C2). + +\boldsymbol\xi_\perp\cdot[(\nabla\times\mathbf B)\times\mathbf B]/B`. The PENTRC matrices ``smat``, ``tmat`` (``idcon_build``) form the same quantity times :math:`\mathcal J B^2` from :math:`\xi^\psi` and - :math:`\xi^\alpha` (G7). + :math:`\xi^\alpha`. ``llagb`` (``gpec_diagnostics_pmodb``, ``pmodbmn_flag``) :math:`\delta B_L(\boldsymbol\xi_\perp)=\delta B_E+\boldsymbol\xi_\perp\cdot\nabla B` - (G5.6). ``lllagb`` is :math:`-(`\ ``Bdivxprp``\ :math:`+`\ ``Bkxprp``\ :math:`)`, - and ``Bkxprp`` :math:`=` ``llagb`` :math:`-` ``eulb`` :math:`+\mu_0p'\xi^\psi/B` (P1). + ``lllagb`` is :math:`-(`\ ``Bdivxprp``\ :math:`+`\ ``Bkxprp``\ :math:`)`, + and ``Bkxprp`` :math:`=` ``llagb`` :math:`-` ``eulb`` :math:`+\mu_0p'\xi^\psi/B`. Relations between the columns ============================= -Gauge of ``lagb`` (G5.5, G5.8): +Gauge of ``lagb``: .. math:: @@ -108,8 +105,8 @@ Gauge of ``lagb`` (G5.5, G5.8): \qquad \xi_\perp^\zeta=-\frac{B_\theta\xi_s/\chi_1+qB_\psi\xi^\psi}{B_\theta+qB_\zeta}. -Kinematic identity (C1; any :math:`\boldsymbol\xi_\perp`, no force balance) -and its pmodb form (C2 with G4): +Kinematic identity (any :math:`\boldsymbol\xi_\perp`, no force balance) +and its pmodb form: .. math:: @@ -118,9 +115,7 @@ and its pmodb form (C2 with G4): =-(\mathrm{Bdivxprp}+\mathrm{Bkxprp})-\boldsymbol\xi_\perp\cdot\mathbf F/B, \qquad \mathbf F=(\nabla\times\mathbf B)\times\mathbf B-\mu_0\nabla p. -C1 and C2 are proved on Cartesian 1-jets. Their direct contraction in the -GPEC chart with :math:`\det g=\mathcal J^2` (G6) is too large for either CAS -to decide and is recorded as an open gap, not as proved. Hence +Hence .. math:: @@ -141,7 +136,7 @@ Regularisation With ``reg_flag`` (Park et al. 2009 Eq. 10), ``gpeq_sol`` replaces the radial derivative by :math:`R\,\partial_\psi\xi^\psi`, :math:`R=s^2/(s^2+\sigma^2)`, :math:`s=m-nq`, in the field and in ``xms``. -The regularised field is then not divergence free (R3): +The regularised field is then not divergence free: .. math:: @@ -150,8 +145,7 @@ The regularised field is then not divergence free (R3): ``Bdivxprp`` uses the same derivative ``xmp1`` as the field, and :math:`\boldsymbol\xi_\perp` is built from the same ``xms``, so the regularisation cancels between the columns and the relation above holds -unchanged with the regularised quantities (R6-R8: a bookkeeping -decomposition under the stated model, not a separate proof about the code). +unchanged with the regularised quantities. Near rational surfaces the regularised field and displacement do not form an exact ideal-induction pair. This is a limitation of the smoothing policy, not evidence of resistive physics. With ``reg_flag=f`` the @@ -166,7 +160,7 @@ PENTRC :math:`\nabla\cdot\boldsymbol\xi` of the *same* displacement; the combination is then invariant under :math:`\boldsymbol\xi\to\boldsymbol\xi+\lambda\mathbf B` up to a total - derivative along the field line (K3). Both PENTRC input paths use the + derivative along the field line. Both PENTRC input paths use the :math:`\boldsymbol\xi_\perp` pair :math:`\delta B/B=-(\nabla\cdot\boldsymbol\xi_\perp+\boldsymbol\kappa\cdot\boldsymbol\xi_\perp)` and :math:`\nabla\cdot\boldsymbol\xi_\perp`: the default path @@ -185,8 +179,8 @@ Bounce-averaged codes in perturbed Boozer or Hamada coordinates (NEO-RT, NEO-2, used as :math:`\delta B_L(\boldsymbol\xi_\perp)`. -Target gauges in GPEC units (Q1-Q8) -=================================== +Target gauges in GPEC units +=========================== Let :math:`D=\partial_\theta+q\partial_\zeta`, so :math:`\mathbf B\cdot\nabla=(\chi_1/\mathcal J)D`. @@ -218,22 +212,9 @@ factor in :math:`\partial_\theta B`. Relabeling equilibrium harmonics with ``jac_out`` alone does not solve this equation. For a PEST chart :math:`\mathcal J/R^2` is a flux function; its -linear perturbation is :math:`\nabla\cdot\boldsymbol\xi-2\xi_R/R` -(Q8). PEST output is not automatically either target gauge. Transform the -equilibrium chart first, then solve the target equation. Q4 proves the -angular-map determinant and Q5 the composition of parallel shifts. - - -Reproducing the proofs -====================== - -:: - - python3 docs/derivations/pmodb_gauge/derive.py - -This needs only SymPy and takes a few minutes, most of it the G6 attempt. -FortSym (https://github.com/lazy-fortran/fortsym) is used instead when it is -on ``PYTHONPATH``. Neither is needed to build or run GPEC. +linear perturbation is :math:`\nabla\cdot\boldsymbol\xi-2\xi_R/R`. +PEST output is not automatically either target gauge. Transform the +equilibrium chart first, then solve the target equation. History