ForceFreeStates - REFACTOR! - Run deck matching through a MatchProblem clause in main

docs/src/galerkin.md

@@ -17,7 +17,9 @@ whichever formalism produced it.  These are the outer-region inputs to resistive
 Select it with `integrator = "galerkin"` in `[ForceFreeStates]`. It replaces the radial ODE
 integration rather than supplementing it: the run computes its own vacuum response at the
 control surface (when `vac_flag`) and produces no free-boundary energies or ODE trace.
-Setting `gal_match_flag` additionally matches the inner layer, giving a driven ξ solution.
+Setting `gal_match_flag` additionally matches the inner layer, giving a driven ξ solution. The
+composite inner-region profiles are a scripting diagnostic: they need `Galerkin(; cut_solution=true)`
+followed by `solve(MatchProblem(...))`, and a deck run does not write them.
 
 The implementation lives in `src/ForceFreeStates/Galerkin/`:
 

src/ForceFreeStates/CoreTypes.jl

@@ -201,9 +201,8 @@ gpec.toml.
     gal_sing_order_ceiling::Bool = true  # auto-raise order by ceil(2·Re(α)) per surface (high Mercier index)
     gal_rpec_flag::Bool = false     # append mpert coil-response columns to the Δ′ solve (RDCON rpec_flag): unit boundary sources whose plasma response is recorded; needed for the driven (resistive perturbed-equilibrium) Δ_gw
     gal_edge_onesided::Bool = false # pack the two end intervals one-sided toward their single rational end (vs the Fortran symmetric "both" pack); avoids the fine edge cell that inflates cond(A). Default false = faithful to gal.f.
-    gal_cut_solution::Bool = false  # also keep the cut solution for a later match's composite inner profiles; implied by gal_match_flag
     # --- DRIVEN (RPEC) outer↔inner asymptotic matching (rmatch match_rpec port) ---
-    gal_match_flag::Bool = false    # enable the RPEC inner-layer matching: solve the coil-driven matched ξ(ψ) from the gal Δ′ + the inner-layer Δ(Q). Requires gal_rpec_flag=true.
+    gal_match_flag::Bool = false    # match the inner layer on the finished solve (coil-driven ξ from Δ′ + inner Δ(Q)); needs gal_rpec_flag=true, or riccati with vac_flag=true
     gal_ideal_flag::Bool = false    # within the match, build the IDEAL solution: skip the inner-layer Δ, use bare coil columns (cout=0). Mirrors Fortran rmatch coil%ideal_flag (the EL reference). eta/rho/rotation ignored.
     gal_inner_solver::String = "ray" # inner-layer Δ backend for the match: "ray" (rotated contour + inward march; recommended for every Q) or "galerkin" (real-axis Hermite FEM; accuracy depends on the surface and arg Q)
     # --- Inner-layer "galerkin" backend overrides (gal_inner_solver = "galerkin"; unset keeps the backend default) ---

src/ForceFreeStates/Galerkin/GalerkinSolve.jl

@@ -46,18 +46,19 @@ end
 
 """
     galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines,
-                   intr::ForceFreeStatesInternal; wv=nothing) -> (GalerkinResult, Union{Nothing,DeltaPrimeData})
+                   intr::ForceFreeStatesInternal; wv=nothing, cut_solution=false) -> (GalerkinResult, Union{Nothing,DeltaPrimeData})
 
 Compute the outer-region Δ′ matching matrix by the singular Galerkin method. Port of `gal_solve`
 (gal.f). Single toroidal mode only (`intr.npert == 1`). Returns the solver internals as a
 `GalerkinResult` alongside the Δ′ payload as the shared [`DeltaPrimeData`](@ref); if there are no
 resonant surfaces in the domain it returns an empty result and `nothing`.
 
 `wv` is the vacuum energy matrix from `free_run`, which supplies the free-boundary edge term
-`wv_edge = wv · psio²`; pass `nothing` for a fixed-boundary edge.
+`wv_edge = wv · psio²`; pass `nothing` for a fixed-boundary edge. `cut_solution` also keeps the
+cut solution, for a later `MatchProblem`'s composite inner profiles.
 """
 function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines,
-    intr::ForceFreeStatesInternal; wv=nothing)
+    intr::ForceFreeStatesInternal; wv=nothing, cut_solution::Bool=false)
 
     intr.npert == 1 || error("galerkin_solve: only single-n (npert == 1) is supported")
     ctrl.gal_solver in ("LU", "cholesky") || error("galerkin_solve: gal_solver must be \"LU\" or \"cholesky\"")
@@ -183,7 +184,7 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines
     # solution only when a later match needs the composite inner profiles.
     ctrl.verbose && @info "Reconstructing outer-region ξ and analytic ξ′ on the gal grid"
     solution = gal_output_solution(ws, asymps, sings, intr, equil.profiles, psihigh;
-        delta=((ctrl.gal_cut_solution || ctrl.gal_match_flag) ? delta : nothing))
+        delta=(cut_solution ? delta : nothing))
 
     sing_psi = [s.psifac for s in sings]
     sing_q = [s.q for s in sings]

src/ForceFreeStates/Integrators.jl

@@ -122,6 +122,7 @@ end
 function _apply_alg!(kwargs::Dict{Symbol,Any}, alg::Galerkin)
     _set_ctrl!(kwargs, :integrator, String(_integrator_symbol(alg)), alg)
     for name in fieldnames(Galerkin)
+        name === :cut_solution && continue    # a solve option, not a control field
         _set_ctrl!(kwargs, Symbol(:gal_, name), getfield(alg, name), alg)
     end
     return kwargs

src/GeneralizedPerturbedEquilibrium.jl

@@ -257,7 +257,7 @@ function main_from_inputs(
     metric, mats, layer_overlap = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles;
         species=kf_species,
         overlap_inputs=_overlap_slayer_inputs(inputs, intr.dir_path))
-    ffs_result = run_force_free_states(ctrl, equil, mats, intr, metric)
+    ffs_result = run_deck_match(run_force_free_states(ctrl, equil, mats, intr, metric), ctrl)
 
     if ctrl.write_outputs_to_HDF5
         write_outputs_to_HDF5(
@@ -747,18 +747,19 @@ function prepare_force_free_states!(
 end
 
 """
-    run_force_free_states(ctrl, equil, mats, intr, metric) -> ForceFreeStatesResult
+    run_force_free_states(ctrl, equil, mats, intr, metric; cut_solution=false) -> ForceFreeStatesResult
 
 Run the formalism selected by `ctrl.integrator` — the standalone Galerkin solve, or the
 Euler-Lagrange sweep with its free-boundary energies and Δ′ BVP — and publish its products as a
-`ForceFreeStatesResult`.
+`ForceFreeStatesResult`. `cut_solution` is the scripting-only `Galerkin(; cut_solution)` option.
 """
 function run_force_free_states(
     ctrl::ForceFreeStatesControl,
     equil::Equilibrium.PlasmaEquilibrium,
     mats,
     intr::ForceFreeStatesInternal,
-    metric
+    metric;
+    cut_solution::Bool=false
 )
     nstring = _mode_range_label(intr)
 
@@ -774,7 +775,7 @@ function run_force_free_states(
             error("integrator = \"galerkin\" does not support kinetic runs (kinetic_factor > 0); use integrator = \"forward\".")
         gal_start = time()
         wv = ctrl.vac_flag ? first(ForceFreeStates.compute_scaled_wv(ctrl, equil, intr)) : nothing
-        gal_data, gal_dp = galerkin_solve(ctrl, equil, mats, intr; wv=wv)
+        gal_data, gal_dp = galerkin_solve(ctrl, equil, mats, intr; wv=wv, cut_solution=cut_solution)
         gal_dt = time() - gal_start
         @info "Galerkin solve completed in $(@sprintf("%.3f", gal_dt)) s"
     else
@@ -822,35 +823,40 @@ function run_force_free_states(
     # Publish the solve: from here on the downstream stages read the result, never `intr`.
     result = build_result(Symbol(ctrl.integrator), ctrl, equil, intr, metric, mats, odet, free_energies, gal_data, gal_dp)
 
-    # Deck-driven matching, through the same MatchProblem as the API.
-    if ctrl.integrator == "galerkin" && ctrl.gal_match_flag
-        ctrl.gal_rpec_flag || error("gal_match_flag=true requires gal_rpec_flag=true")
-        ctrl.gal_ideal_flag || !(isempty(ctrl.gal_eta) || isempty(ctrl.gal_rho) || isempty(ctrl.gal_rotation)) ||
-            error("gal_match_flag needs gal_eta, gal_rho and gal_rotation (one value per matched surface, core to edge)")
-        ctrl.gal_inner_solver in ("ray", "galerkin") ||
-            error("gal_inner_solver = \"$(ctrl.gal_inner_solver)\" (expected \"ray\" or \"galerkin\")")
-        ctrl.verbose && @info(
-            ctrl.gal_ideal_flag ?
-            "RPEC matching: IDEAL solution (inner layer skipped, bare coil columns)" :
-            "RPEC matching: inner-layer Δ(Q) + outer↔inner solve for the coil-driven ξ"
-        )
-        # The gal_inner_* knobs override the "galerkin" backend's defaults; unset keys keep them.
-        model = if ctrl.gal_inner_solver == "galerkin"
-            knobs = (xfac=ctrl.gal_inner_xfac, nx=ctrl.gal_inner_nx, nq=ctrl.gal_inner_nq,
-                cutoff=ctrl.gal_inner_cutoff, kmax=ctrl.gal_inner_kmax)
-            InnerLayer.GGJModel(; solver=:galerkin, (k => v for (k, v) in pairs(knobs) if v !== nothing)...)
-        else
-            InnerLayer.GGJModel(; solver=:ray)
-        end
-        prob = ForceFreeStates.MatchProblem(result; eta=ctrl.gal_eta, rho=ctrl.gal_rho, rotation=ctrl.gal_rotation,
-            gamma=ctrl.gal_gamma, ideal=ctrl.gal_ideal_flag)
-        result = solve(prob, model)
-        ctrl.gal_ideal_flag || (ctrl.verbose && @info "RPEC matching: linear-solve residual = $(result.match.residual)")
-    end
-
     return result
 end
 
+"""
+    run_deck_match(result, ctrl) -> ForceFreeStatesResult
+
+Solve the deck's inner-layer match (`gal_match_flag`) on a finished solve, posing the same
+`MatchProblem` a script would from the `gal_*` keys. Returns `result` unchanged when the deck
+requests no match.
+"""
+function run_deck_match(result::ForceFreeStatesResult, ctrl::ForceFreeStatesControl)
+    ctrl.gal_match_flag || return result
+    ctrl.gal_inner_solver in ("ray", "galerkin") ||
+        error("gal_inner_solver = \"$(ctrl.gal_inner_solver)\" (expected \"ray\" or \"galerkin\")")
+    ctrl.verbose && @info(
+        ctrl.gal_ideal_flag ?
+        "RPEC matching: IDEAL solution (inner layer skipped, bare coil columns)" :
+        "RPEC matching: inner-layer Δ(Q) + outer↔inner solve for the coil-driven ξ"
+    )
+    # The gal_inner_* knobs override the "galerkin" backend's defaults; unset keys keep them.
+    model = if ctrl.gal_inner_solver == "galerkin"
+        knobs = (xfac=ctrl.gal_inner_xfac, nx=ctrl.gal_inner_nx, nq=ctrl.gal_inner_nq,
+            cutoff=ctrl.gal_inner_cutoff, kmax=ctrl.gal_inner_kmax)
+        InnerLayer.GGJModel(; solver=:galerkin, (k => v for (k, v) in pairs(knobs) if v !== nothing)...)
+    else
+        InnerLayer.GGJModel(; solver=:ray)
+    end
+    prob = ForceFreeStates.MatchProblem(result; eta=ctrl.gal_eta, rho=ctrl.gal_rho, rotation=ctrl.gal_rotation,
+        gamma=ctrl.gal_gamma, ideal=ctrl.gal_ideal_flag)
+    matched = solve(prob, model)
+    ctrl.gal_ideal_flag || (ctrl.verbose && @info "RPEC matching: linear-solve residual = $(matched.match.residual)")
+    return matched
+end
+
 """
     EulerLagrangeProblem(equil; nn, wall=Vacuum.WallShapeSettings(),
                          dir_path=".", debug=DebugSettings(), kwargs...)
@@ -939,6 +945,8 @@ function solve(prob::EulerLagrangeProblem, alg::ForceFreeStates.AbstractIntegrat
     intr.debug_settings = prob.debug
 
     resolve_mode_space!(intr, ctrl)
+    ctrl.gal_match_flag &&
+        error("gal_match_flag is the deck's match request; match the result instead: solve(MatchProblem(ffs; eta, rho, rotation), GGJModel())")
 
     # The API path never reads kinetic profiles, so the KineticForces control is only the
     # placeholder `prepare_force_free_states!` threads into its (unused) callback.
@@ -953,7 +961,7 @@ function solve(prob::EulerLagrangeProblem, alg::ForceFreeStates.AbstractIntegrat
 
     locstab, ballooning_boundary = run_local_stability(ctrl, equil)
     metric, mats, _ = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, nothing)
-    result = run_force_free_states(ctrl, equil, mats, intr, metric)
+    result = run_force_free_states(ctrl, equil, mats, intr, metric; cut_solution=alg isa ForceFreeStates.Galerkin && alg.cut_solution)
 
     if ctrl.write_outputs_to_HDF5
         write_outputs_to_HDF5(result; locstab=locstab, ballooning_boundary=ballooning_boundary)

test/runtests_solve_api.jl

@@ -226,6 +226,7 @@ using HDF5
         @test_throws ErrorException solve(equil, Forward(); nn=1, dir_path=".", ffs_kwargs..., integrator="riccati")
         @test_throws ErrorException solve(equil, Riccati(); nn=1, dir_path=".", ffs_kwargs..., nchunks=8)
         @test_throws ErrorException solve(equil, Forward(); nn=1, dir_path=".", ffs_kwargs..., nn_low=2)
+        @test_throws ErrorException solve(equil, Galerkin(; nx=32, rpec_flag=true); nn=1, dir_path=".", ffs_kwargs..., gal_match_flag=true)
     end
 
     # Path of the first field where `a` and `b` differ (recursing into structs, arrays and
@@ -390,27 +391,42 @@ using HDF5
 
         @testset "deck matching under gal_match_flag" begin
             layer = (eta=1e-7, rho=1e-7, rotation=1.0)
-            # The galerkin inner model returns real-ψ profiles, so the composite inner profiles are compared too.
+            # The deck's match keys for `n` surfaces on top of the integrator's own.
+            match_keys(n, solver) = Dict{String,Any}("gal_match_flag" => true, "gal_inner_solver" => solver,
+                "gal_eta" => fill(layer.eta, n), "gal_rho" => fill(layer.rho, n), "gal_rotation" => fill(layer.rotation, n))
+
             mktempdir() do dir
-                gal = solve(equil, Galerkin(; nx=32, rpec_flag=true, cut_solution=true); nn=1, wall=wall, dir_path=dir, ffs_kwargs...)
+                gal = solve(equil, Galerkin(; nx=32, rpec_flag=true); nn=1, wall=wall, dir_path=dir, ffs_kwargs...)
+                # The galerkin inner model returns real-ψ profiles, so the inner profiles are compared too.
                 api = solve(MatchProblem(gal; layer...), GGJModel(; solver=:galerkin))
                 n = length(gal.surfaces)
-                reference = GPEC.main([_stage_deck(dir;
-                    ffs_overrides=Dict{String,Any}(
-                        "integrator" => "galerkin", "gal_nx" => 32, "gal_rpec_flag" => true, "gal_match_flag" => true,
-                        "gal_inner_solver" => "galerkin", "gal_eta" => fill(layer.eta, n), "gal_rho" => fill(layer.rho, n),
-                        "gal_rotation" => fill(layer.rotation, n)))]).ffs
+                overrides = merge(Dict{String,Any}("integrator" => "galerkin", "gal_nx" => 32, "gal_rpec_flag" => true), match_keys(n, "galerkin"))
+                reference = GPEC.main([_stage_deck(dir; ffs_overrides=overrides)]).ffs
 
                 @test reference.closure === api.closure === :matched
                 @test api.bpen ≈ reference.bpen rtol = 1e-12
                 @test api.match.xi ≈ reference.match.xi rtol = 1e-12
-                @test length(api.match.inner_xi) == length(reference.match.inner_xi) == n
-                @test all(isapprox.(api.match.inner_xi, reference.match.inner_xi; rtol=1e-12))
+                # The composite inner profiles are API-only: a deck never keeps the cut solution.
+                @test isempty(reference.match.inner_xi) && isempty(api.match.inner_xi)
+                cut = solve(equil, Galerkin(; nx=32, rpec_flag=true, cut_solution=true); nn=1, wall=wall, dir_path=dir, ffs_kwargs...)
+                @test length(solve(MatchProblem(cut; layer...), GGJModel(; solver=:galerkin)).match.inner_xi) == n
+            end
+
+            mktempdir() do dir
+                ric = solve(equil, Riccati(); nn=1, wall=wall, dir_path=dir, ffs_kwargs...)
+                api = solve(MatchProblem(ric; layer...), GGJModel(; solver=:ray))
+                overrides = merge(Dict{String,Any}("integrator" => "riccati"), match_keys(length(ric.surfaces), "ray"))
+                reference = GPEC.main([_stage_deck(dir; ffs_overrides=overrides)]).ffs
 
-                uncut = solve(equil, Galerkin(; nx=32, rpec_flag=true); nn=1, wall=wall, dir_path=dir, ffs_kwargs...)
-                matched = solve(MatchProblem(uncut; layer...), GGJModel(; solver=:galerkin))
-                # Known drift: gal_match_flag implies the cut solution on the deck; the API keeps it only when asked.
-                @test_broken length(matched.match.inner_xi) == length(reference.match.inner_xi)
+                @test reference.closure === :matched
+                @test api.bpen ≈ reference.bpen rtol = 1e-12
+                @test api.match.cout ≈ reference.match.cout rtol = 1e-12
+            end
+
+            # A forward solve has no Δ′ to match against, so the deck errors instead of ignoring the flag.
+            mktempdir() do dir
+                @test_throws ErrorException GPEC.main([_stage_deck(dir;
+                    ffs_overrides=Dict{String,Any}("gal_match_flag" => true, "gal_ideal_flag" => true))])
             end
         end
     end