Feature: add SLAYER and GGJ tearing growth rates

Project.toml

@@ -7,9 +7,11 @@ version = "0.1.0"
 [deps]
 AdaptiveArrayPools = "4f381ef7-9af0-4cbe-99d4-cf36d7b0f233"
 Contour = "d38c429a-6771-53c6-b99e-75d170b6e991"
+DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FastInterpolations = "9ea80cae-fc13-4c00-8066-6eaedb12f34b"
@@ -23,6 +25,7 @@ PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
@@ -34,9 +37,11 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 [compat]
 AdaptiveArrayPools = "0.3.5"
 Contour = "0.6.3"
+DelaunayTriangulation = "1.6.6"
 DelimitedFiles = "1.9.1"
 DiffEqCallbacks = "4.9.0"
 Documenter = "1.14.1"
+DoubleFloats = "1.6.2"
 FFTW = "1.9.0"
 FastGaussQuadrature = "1.1.0"
 FastInterpolations = "0.4"
@@ -50,6 +55,7 @@ PlotlyJS = "0.18.17"
 Plots = "1.40.15"
 Printf = "1"
 QuadGK = "2.11.3"
+Random = "1"
 Roots = "2.2.13"
 SparseArrays = "1"
 SpecialFunctions = "2.5.1"

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_integration_paths.jl

@@ -0,0 +1,148 @@
+#!/usr/bin/env julia
+"""
+Benchmark the three integration paths (standard, riccati, parallel) on Solovev and DIIID examples.
+Runs in a single Julia process to avoid measuring compilation overhead.
+Produces accuracy and performance tables similar to PR #178.
+
+Usage:
+    julia --project=. -t4 benchmarks/benchmark_integration_paths.jl
+"""
+
+using GeneralizedPerturbedEquilibrium
+using HDF5, Printf, TOML
+
+const PROJECT_ROOT = abspath(joinpath(@__DIR__, ".."))
+
+struct BenchResult
+    example::String
+    path::String
+    et1::Float64
+    nsteps::Int
+    runtime::Float64
+end
+
+function run_one(example_dir::String, path_name::String; num_warm::Int=2)
+    abs_dir = abspath(example_dir)
+    gpec_toml = joinpath(abs_dir, "gpec.toml")
+
+    # Read and modify config
+    config = TOML.parsefile(gpec_toml)
+    ffs = get(config, "ForceFreeStates", Dict{String,Any}())
+    if path_name == "standard"
+        ffs["use_riccati"] = false
+        ffs["use_parallel"] = false
+    elseif path_name == "riccati"
+        ffs["use_riccati"] = true
+        ffs["use_parallel"] = false
+    elseif path_name == "parallel"
+        ffs["use_riccati"] = false
+        ffs["use_parallel"] = true
+    end
+    config["ForceFreeStates"] = ffs
+
+    # Write modified config in-place, restore after
+    original_toml = read(gpec_toml, String)
+
+    try
+        open(gpec_toml, "w") do f
+            TOML.print(f, config)
+        end
+
+        # JIT warmup
+        println("  [$path_name] JIT warmup...")
+        GeneralizedPerturbedEquilibrium.main([abs_dir])
+
+        # Timed runs
+        runtimes = Float64[]
+        for i in 1:num_warm
+            println("  [$path_name] Warm run $i/$num_warm...")
+            t0 = time()
+            GeneralizedPerturbedEquilibrium.main([abs_dir])
+            push!(runtimes, time() - t0)
+            @printf("    %.2f s\n", runtimes[end])
+        end
+
+        # Read results
+        gpec_h5 = joinpath(abs_dir, "gpec.h5")
+        et1, nsteps = h5open(gpec_h5, "r") do h5
+            et = read(h5["vacuum/et"])
+            ns = read(h5["integration/nstep"])
+            (real(et[1]), ns)
+        end
+
+        avg_t = sum(runtimes) / length(runtimes)
+        return BenchResult(basename(example_dir), path_name, et1, nsteps, avg_t)
+    finally
+        write(gpec_toml, original_toml)
+    end
+end
+
+function main()
+    examples = [
+        joinpath(PROJECT_ROOT, "examples", "Solovev_ideal_example"),
+        joinpath(PROJECT_ROOT, "examples", "DIIID-like_ideal_example"),
+    ]
+    paths = ["standard", "riccati", "parallel"]
+
+    results = BenchResult[]
+    for ex in examples
+        println("\n" * "="^60)
+        println("Example: $(basename(ex))")
+        println("="^60)
+        for p in paths
+            r = run_one(ex, p)
+            push!(results, r)
+            @printf("  → et[1]=%.5f  steps=%d  time=%.2fs\n", r.et1, r.nsteps, r.runtime)
+        end
+    end
+
+    # Print Accuracy table
+    println("\n\n## Accuracy\n")
+    println("| Example | Path | et[1] | Error vs std |")
+    println("|---------|------|-------|--------------|")
+    for ex in unique(r.example for r in results)
+        group = filter(r -> r.example == ex, results)
+        std_et1 = group[1].et1
+        N = 0
+        toml_path = joinpath(PROJECT_ROOT, "examples", ex, "gpec.toml")
+        if isfile(toml_path)
+            cfg = TOML.parsefile(toml_path)
+            ffs_cfg = get(cfg, "ForceFreeStates", Dict())
+            mlow = get(ffs_cfg, "delta_mlow", 8)
+            mhigh = get(ffs_cfg, "delta_mhigh", 8)
+            N = mlow + mhigh
+        end
+        for r in group
+            err_str = r.path == "standard" ? "—" : @sprintf("%.3f%%", 100*abs(r.et1 - std_et1)/abs(std_et1))
+            short_ex = startswith(r.example, "Solovev") ? "Solovev N=$N" : "DIIID N=$N"
+            @printf("| %s | %s | %.5f | %s |\n", short_ex, r.path, r.et1, err_str)
+        end
+    end
+
+    # Print Performance table
+    nthreads = Threads.nthreads()
+    println("\n## Performance ($nthreads threads)\n")
+    println("| Example | Path | Time | Speedup |")
+    println("|---------|------|------|---------|")
+    for ex in unique(r.example for r in results)
+        group = filter(r -> r.example == ex, results)
+        std_time = group[1].runtime
+        N = 0
+        toml_path = joinpath(PROJECT_ROOT, "examples", ex, "gpec.toml")
+        if isfile(toml_path)
+            cfg = TOML.parsefile(toml_path)
+            ffs_cfg = get(cfg, "ForceFreeStates", Dict())
+            mlow = get(ffs_cfg, "delta_mlow", 8)
+            mhigh = get(ffs_cfg, "delta_mhigh", 8)
+            N = mlow + mhigh
+        end
+        for r in group
+            speedup = std_time / r.runtime
+            short_ex = startswith(r.example, "Solovev") ? "Solovev N=$N" : "DIIID N=$N"
+            speedup_str = r.path == "standard" ? "1.00×" : @sprintf("**%.2f×**", speedup)
+            @printf("| %s | %s | %.2fs | %s |\n", short_ex, r.path, r.runtime, speedup_str)
+        end
+    end
+end
+
+main()