ForceFreeStates: Dual Riccati + Parallel FM integration + Δ' output

Project.toml

@@ -20,6 +20,7 @@ Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -43,6 +44,7 @@ Pkg = "1.11.0"
 Plots = "1.40.15"
 Printf = "1.11.0"
 SparseArrays = "1.11.0"
+Random = "1.11.0"
 Roots = "2.2.13"
 SpecialFunctions = "2.5.1"
 StaticArrays = "1.9.15"

benchmarks/benchmark_delta_prime_methods.jl

@@ -0,0 +1,95 @@
+# Sanity check: compute_delta_prime_from_ca! vs inline Δ' from riccati_cross_ideal_singular_surf!
+#
+# riccati_cross_ideal_singular_surf! computes Δ' inline at each singular surface crossing
+# using the diagonal formula (no Gaussian reduction permutation):
+#   Δ'[s] = (ca_r[ipert_res, ipert_res, 2, s] - ca_l[ipert_res, ipert_res, 2, s]) / (4π²·ψ₀)
+#
+# compute_delta_prime_from_ca! applies the identical formula post-hoc from the stored
+# ca_l/ca_r arrays. Since both operate on the same data with the same formula, results
+# should match to floating-point precision (not just approximately — exactly).
+#
+# This verifies that compute_delta_prime_from_ca! is a correct standalone implementation
+# of the Δ' formula that can be used for testing or alternative integration drivers.
+#
+# Usage (from JPEC_main root):
+#   julia --project=. benchmarks/benchmark_delta_prime_methods.jl
+
+using LinearAlgebra, Printf, TOML
+using GeneralizedPerturbedEquilibrium
+
+const FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates
+
+function setup_and_run_solovev()
+    ex = joinpath(@__DIR__, "..", "test", "test_data", "regression_solovev_ideal_example")
+    inputs = TOML.parsefile(joinpath(ex, "gpec.toml"))
+    inputs["ForceFreeStates"]["verbose"] = false
+    inputs["ForceFreeStates"]["use_riccati"] = true
+    intr = FFS.ForceFreeStatesInternal(; dir_path=ex)
+    ctrl = FFS.ForceFreeStatesControl(;
+        (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
+    eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex)
+    equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config)
+    intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
+        (Symbol(k) => v for (k, v) in inputs["Wall"])...)
+    FFS.sing_lim!(intr, ctrl, equil)
+    intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; intr.npert = 1
+    FFS.sing_find!(intr, equil)
+    intr.mlow  = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow
+    intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh
+    intr.mpert = intr.mhigh - intr.mlow + 1
+    intr.mband = intr.mpert - 1
+    intr.numpert_total = intr.mpert * intr.npert
+    metric = FFS.make_metric(equil; mband=intr.mband, fft_flag=ctrl.fft_flag)
+    ffit = FFS.make_matrix(equil, intr, metric)
+    odet = FFS.riccati_eulerlagrange_integration(ctrl, equil, ffit, intr)
+    return ctrl, equil, ffit, intr, odet
+end
+
+println("\n=== compute_delta_prime_from_ca! consistency check ===")
+println("Verifies the standalone Δ' formula matches the inline Riccati crossing computation.")
+println("Expected error: exactly zero (same formula, same data).\n")
+
+ctrl, equil, ffit, intr, odet = setup_and_run_solovev()
+msing = intr.msing
+
+# Capture Δ' values set inline by riccati_cross_ideal_singular_surf! during integration
+delta_prime_inline = [copy(intr.sing[s].delta_prime) for s in 1:msing]
+
+# Now call compute_delta_prime_from_ca! — it reads the same ca_l/ca_r arrays and
+# overwrites intr.sing[s].delta_prime using the identical diagonal formula
+FFS.compute_delta_prime_from_ca!(odet, intr, equil)
+
+println("  N=$(intr.numpert_total) modes, $msing singular surfaces\n")
+@printf("  %6s  %4s  %4s  %22s  %22s  %12s\n",
+        "Surf", "m", "n", "Δ' (inline)", "Δ' (from_ca)", "abs diff")
+println("  " * "-"^76)
+
+max_absdiff = let max_absdiff = 0.0
+    for s in 1:msing
+        sing = intr.sing[s]
+        dp_from_ca = intr.sing[s].delta_prime
+        for i in eachindex(delta_prime_inline[s])
+            dp_il  = delta_prime_inline[s][i]
+            dp_fc  = dp_from_ca[i]
+            absdiff = abs(dp_fc - dp_il)
+            max_absdiff = max(max_absdiff, absdiff)
+            @printf("  %6d  %4d  %4d  %22.6f%+.6fi  %22.6f%+.6fi  %12.4e\n",
+                    s, sing.m[i], sing.n[i],
+                    real(dp_il), imag(dp_il),
+                    real(dp_fc), imag(dp_fc),
+                    absdiff)
+        end
+    end
+    max_absdiff
+end
+
+println()
+if max_absdiff == 0.0
+    println("PASSED — Δ' values are bit-for-bit identical (max abs diff = 0.0)")
+elseif max_absdiff < 1e-14
+    @printf("PASSED — max abs diff = %.2e (floating-point rounding only)\n", max_absdiff)
+else
+    @printf("FAILED — max abs diff = %.2e (expected exact agreement)\n", max_absdiff)
+    exit(1)
+end
+println()

benchmarks/benchmark_riccati_der.jl

@@ -0,0 +1,131 @@
+# Sanity check: riccati_der! correctly evaluates the explicit Riccati ODE.
+#
+# riccati_der! implements [Glasser 2018 Phys. Plasmas 25, 032507, Eq. 19]:
+#   dS/dψ = w†·F̄⁻¹·w - S·Ḡ·S,   w = Q - K̄·S
+#
+# where Q = diag(1/(m - n·q)), F̄ = L·L† (Cholesky), K̄ and Ḡ are the MHD
+# metric matrices evaluated at ψ.
+#
+# NOTE: The identity between this Riccati ODE and the EL chain rule
+#   dS/dψ = dU₁·U₂⁻¹ - S·dU₂·U₂⁻¹
+# holds ONLY for Hermitian S (physical states evolved from the axis, where
+# S†=S is preserved by the EL symmetry). For arbitrary non-Hermitian (U₁, U₂),
+# the two expressions differ — so this script compares riccati_der! against the
+# explicit formula rather than against sing_der!.
+#
+# Usage (from JPEC_main root):
+#   julia --project=. benchmarks/benchmark_riccati_der.jl
+
+using LinearAlgebra, Random, Printf, TOML
+using GeneralizedPerturbedEquilibrium
+
+const FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates
+
+function setup_solovev()
+    ex = joinpath(@__DIR__, "..", "test", "test_data", "regression_solovev_ideal_example")
+    inputs = TOML.parsefile(joinpath(ex, "gpec.toml"))
+    inputs["ForceFreeStates"]["verbose"] = false
+    intr = FFS.ForceFreeStatesInternal(; dir_path=ex)
+    ctrl = FFS.ForceFreeStatesControl(;
+        (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
+    eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex)
+    equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config)
+    intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
+        (Symbol(k) => v for (k, v) in inputs["Wall"])...)
+    FFS.sing_lim!(intr, ctrl, equil)
+    intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; intr.npert = 1
+    FFS.sing_find!(intr, equil)
+    intr.mlow  = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow
+    intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh
+    intr.mpert = intr.mhigh - intr.mlow + 1
+    intr.mband = intr.mpert - 1
+    intr.numpert_total = intr.mpert * intr.npert
+    metric = FFS.make_metric(equil; mband=intr.mband, fft_flag=ctrl.fft_flag)
+    ffit = FFS.make_matrix(equil, intr, metric)
+    return ctrl, equil, ffit, intr
+end
+
+# Evaluate the Riccati RHS explicitly from splines: dS = w†·F̄⁻¹·w - S·Ḡ·S
+function riccati_rhs_manual(S, psi, equil, ffit, intr)
+    N = intr.numpert_total
+    L    = zeros(ComplexF64, N, N)
+    Kmat = zeros(ComplexF64, N, N)
+    Gmat = zeros(ComplexF64, N, N)
+    ffit.fmats_lower(vec(L),    psi; hint=ffit._hint)
+    ffit.kmats(vec(Kmat), psi; hint=ffit._hint)
+    ffit.gmats(vec(Gmat), psi; hint=ffit._hint)
+
+    q = equil.profiles.q_spline(psi)
+    singfac = vec(1.0 ./ ((intr.mlow:intr.mhigh) .- q .* (intr.nlow:intr.nhigh)'))
+
+    # w = Q - K̄·S  (Q is diagonal; add only the diagonal entries)
+    w = -Kmat * S
+    for i in 1:N
+        w[i, i] += singfac[i]
+    end
+
+    # v = F̄⁻¹·w  via stored Cholesky factor L (L·L† = F̄)
+    v = copy(w)
+    ldiv!(LowerTriangular(L), v)
+    ldiv!(UpperTriangular(L'), v)
+
+    return adjoint(w) * v - S * Gmat * S
+end
+
+println("\n=== riccati_der! formula verification ===")
+println("Verifies riccati_der! output matches manual evaluation of Glasser 2018 Eq. 19.")
+println("Test state: Hermitian S (physical constraint). Expected error: ~machine epsilon.\n")
+
+ctrl, equil, ffit, intr = setup_solovev()
+N = intr.numpert_total
+
+odet = FFS.OdeState(N, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing)
+FFS.initialize_el_at_axis!(odet, ctrl, equil.profiles, intr)
+chunks = FFS.chunk_el_integration_bounds(odet, ctrl, intr)
+
+# 30% into each chunk: well inside the interval, away from singularities at psi_end
+test_psis = [c.psi_start + 0.3 * (c.psi_end - c.psi_start) for c in chunks]
+
+println("  N=$N modes, $(length(test_psis)) test ψ points (30% into each chunk)\n")
+@printf("  %8s  %14s  %14s  %12s\n", "ψ", "‖dS_manual‖", "‖dS_ric‖", "rel error")
+println("  " * "-"^54)
+
+rng = Random.MersenneTwister(42)
+threshold = 1e-10
+
+max_err = let max_err = 0.0
+    for psi in test_psis
+        # Hermitian S: physical Riccati matrix is Hermitian (preserved by EL symmetry)
+        A = randn(rng, ComplexF64, N, N)
+        S = (A + A') / 2   # Hermitian by construction
+
+        # Manual RHS
+        dS_manual = riccati_rhs_manual(S, psi, equil, ffit, intr)
+
+        # riccati_der! RHS
+        u_ric  = zeros(ComplexF64, N, N, 2)
+        du_ric = zeros(ComplexF64, N, N, 2)
+        u_ric[:, :, 1] .= S
+        u_ric[:, :, 2] .= Matrix{ComplexF64}(I, N, N)
+        dummy_chunk = FFS.IntegrationChunk(psi, psi, false, 0, 1)
+        params = (ctrl, equil, ffit, intr, odet, dummy_chunk)
+        FFS.riccati_der!(du_ric, u_ric, params, psi)
+        dS_ric = du_ric[:, :, 1]
+
+        ref = max(norm(dS_manual), 1e-10)
+        err = norm(dS_ric - dS_manual) / ref
+        max_err = max(max_err, err)
+        status = err < threshold ? "" : "  ← FAIL"
+        @printf("  %8.4f  %14.4e  %14.4e  %12.4e%s\n", psi, norm(dS_manual), norm(dS_ric), err, status)
+    end
+    max_err
+end
+
+println()
+if max_err < threshold
+    @printf("PASSED — max rel error = %.2e (threshold %.0e)\n", max_err, threshold)
+else
+    @printf("FAILED — max rel error = %.2e exceeds threshold %.0e\n", max_err, threshold)
+    exit(1)
+end
+println()

benchmarks/benchmark_threads.jl

@@ -0,0 +1,76 @@
+# Thread-scaling benchmark for the bidirectional parallel FM integration.
+# Runs the Solovev (N=8) and DIIID-like (N=26) examples with use_parallel=true
+# across 1, 2, 4, 8 threads and compares against the serial Riccati path.
+#
+# Usage (from JPEC_main root):
+#   for t in 1 2 4 8; do julia -t $t --project=. benchmarks/benchmark_threads.jl; done
+
+using GeneralizedPerturbedEquilibrium, TOML, Printf, Statistics
+
+function run_ffs(ex; use_parallel, use_riccati=false)
+    inputs = TOML.parsefile(joinpath(ex, "gpec.toml"))
+    inputs["ForceFreeStates"]["verbose"] = false
+    inputs["ForceFreeStates"]["use_parallel"] = use_parallel
+    inputs["ForceFreeStates"]["use_riccati"] = use_riccati
+    inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false
+    intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal(; dir_path=ex)
+    ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(;
+        (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
+    eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex)
+    equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config)
+    intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
+        (Symbol(k) => v for (k, v) in inputs["Wall"])...)
+    GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
+    intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; intr.npert = 1
+    GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_find!(intr, equil)
+    intr.mlow  = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow
+    intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh
+    intr.mpert = intr.mhigh - intr.mlow + 1
+    intr.mband = intr.mpert - 1
+    intr.numpert_total = intr.mpert * intr.npert
+    metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil; mband=intr.mband, fft_flag=ctrl.fft_flag)
+    ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric)
+    odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, ffit, intr)
+    vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run!(odet, ctrl, equil, ffit, intr)
+    return real(vac.et[1]), intr.numpert_total
+end
+
+function timed_run(ex; use_parallel, use_riccati=false, nwarm=1, nrep=2)
+    # Warmup
+    for _ in 1:nwarm
+        run_ffs(ex; use_parallel, use_riccati)
+    end
+    # Timed runs
+    times = Float64[]
+    local et1, N
+    for _ in 1:nrep
+        t0 = time()
+        et1, N = run_ffs(ex; use_parallel, use_riccati)
+        push!(times, time() - t0)
+    end
+    return mean(times), et1, N
+end
+
+nthreads = Threads.nthreads()
+root     = joinpath(@__DIR__, "..")
+sol_ex   = joinpath(root, "test", "test_data", "regression_solovev_ideal_example")
+diiid_ex = joinpath(root, "examples", "DIIID-like_ideal_example")
+
+println("\n=== Thread-scaling benchmark ($(nthreads) thread(s)) ===\n")
+
+for (label, ex) in [("Solovev", sol_ex), ("DIIID-like", diiid_ex)]
+    t_std,    et_std,  N = timed_run(ex; use_parallel=false, use_riccati=false)
+    t_ric,    et_ric,  _ = timed_run(ex; use_parallel=false, use_riccati=true)
+    t_par,    et_par,  _ = timed_run(ex; use_parallel=true,  use_riccati=false)
+
+    err_ric = abs(et_ric - et_std) / abs(et_std) * 100
+    err_par = abs(et_par - et_std) / abs(et_std) * 100
+
+    println("$label (N=$N, nthreads=$nthreads)")
+    @printf("  standard   et[1]=%.5f  t=%.2fs  speedup=1.00×\n", et_std, t_std)
+    @printf("  riccati    et[1]=%.5f  t=%.2fs  speedup=%.2f×  err=%.4f%%\n",
+            et_ric, t_ric, t_std/t_ric, err_ric)
+    @printf("  parallel   et[1]=%.5f  t=%.2fs  speedup=%.2f×  err=%.4f%%\n",
+            et_par, t_par, t_std/t_par, err_par)
+    println()
+end

docs/make.jl

@@ -26,6 +26,7 @@ makedocs(;
         "API Reference" => [
             "Vacuum" => "vacuum.md",
             "Equilibrium" => "equilibrium.md",
+            "Stability Analysis" => "stability.md",
             "Utilities" => "utilities.md",
             "Forcing Terms" => "forcing_terms.md",
             "Perturbed Equilibrium" => "perturbed_equilibrium.md"

docs/src/equilibrium.md

@@ -104,6 +104,7 @@ Notes:
 
 ## See also
 
+- `docs/src/stability.md` — ideal MHD stability analysis built on top of the equilibrium
 - `docs/src/splines.md` — spline helpers used by equilibrium routines
 - `docs/src/vacuum.md` — coupling between equilibrium and vacuum solvers