Tearing.IslandFlattening - FEATURE - Flatten the pressure across islands of prescribed width and report the stored energy they remove

docs/src/inner_layer.md

@@ -12,6 +12,7 @@ Two layer models are available, selected by `[Tearing] inner_model`:
 
 For each surface the `Dispersion` search scans the complex frequency plane, extracts the roots of the dispersion relation and classifies them, and `[Tearing.search]` sets the scan.
 With `[Tearing.CriticalResonantField]` the stage also computes, for the `SLAYER` layer, the resonant field at which a static error field overcomes the layer's shielding: the error-field penetration threshold.
+`[Tearing.IslandFlattening]` flattens the equilibrium pressure across islands of given width and reports the stored energy and pressure they remove; it needs no layer solve and runs alone when it is the only part of `[Tearing]`.
 
 The same GGJ layer serves the `[Match]` stage of the [Workflow](workflow.md), which does the complementary job: instead of finding the free modes it imposes the layer response on the driven force-free states and writes the matched solution to `SingularSurfaces/Match/`.
 
@@ -48,6 +49,24 @@ Each surface uses its own stability parameter ``\alpha = S^{-1/3}(\Delta_\mathrm
 Penetration is the first resolved saddle-node bifurcation on the path from ``Q_0`` toward the electron fluid at ``Q_e``, above which the shielded branch no longer exists; small numerical extrema and the unshielded zone next to ``Q_e`` are excluded, and a surface without a resolved bifurcation reports NaN.
 `br_crit` is normalized by the toroidal field, and `br_crit_tesla = br_crit·|B0|` is the quantity to compare with `PerturbedEquilibrium/SingularCoupling/resonant_area_weighted_field`, the total resonant field of the ideal response; the comparison is approximate, since the theory's ``b_r`` is a cylindrical field at the layer.
 
+### Island flattening
+
+Inside a magnetic island the pressure is flattened along the field, so ``dp/d\psi = 0`` there; outside it the gradient is unchanged.
+For a set of islands of full width ``w`` (ψ_N) centred on their rational surfaces, the flattened intervals are merged where islands overlap, so an overlap is counted once, and the pressure is anchored at its value on the last profile point and integrated inward, as in Fortran GPEC's saturated-island profile.
+For merged intervals ``[a_k, b_k]`` the integral is exact:
+
+```math
+p_\mathrm{mod}(\psi) = p(\psi) + \sum_{k:\, b_k > \psi} \left[p(b_k) - p(\max(\psi, a_k))\right],
+\qquad
+\Delta W = \tfrac{3}{2}\int (p - p_\mathrm{mod})\, dV .
+```
+
+Every point inward of an interval drops by ``p(a_k) - p(b_k)``, the pedestal-top pressure change when the interval is at the edge.
+Given a bootstrap profile (the `jbs_B` keyword of `island_flattening`), the result also carries the part the islands remove: to leading order ``j_\mathrm{bs} \propto p'``, so it vanishes on the flattened intervals and is unchanged elsewhere, and `jbs_fraction_lost` is the fraction of its volume integral on them.
+For ``\langle j_\mathrm{bs}\cdot B\rangle`` that is a fraction of ``\int\langle j_\mathrm{bs}\cdot B\rangle\,dV``; for the toroidal bootstrap current pass ``j_\phi\langle R^{-1}\rangle``, since ``dA = \langle R^{-1}\rangle\,dV/2\pi`` on axisymmetric surfaces.
+The pressure is anchored at the last point of the equilibrium profile grid, so for an island reaching the edge the result depends on where that grid ends.
+The flattening is static: the equilibrium is not re-solved with the modified pressure, and island widths are inputs rather than evolved.
+
 ## Running it
 
 From `examples/DIIID-like_SLAYER_example/gpec.toml`, after a `[ForceFreeStates]` solve with `integrator = "Riccati"` and `vac_flag = true` (or `"Galerkin"` with `rpec_flag = true`), which provides the Δ′ matrix, and with a kinetic file on the equilibrium:
@@ -78,6 +97,15 @@ D_c_type = "none"                 # critical-Δ offset from parallel transport (
 P_perp_model = "D_perp"           # how the perpendicular Prandtl number is built (below)
 ```
 
+The flattening takes full island widths in ψ_N keyed by surface; every q = m/n surface of the equilibrium gets its key's width:
+
+```toml
+[Tearing.IslandFlattening]
+widths = { "2/1" = 0.05, "3/2" = 0.03 }
+```
+
+In a script, `island_flattening(ffs.equil; widths=Dict("2/1" => 0.05))` gives the same result.
+
 In a script, `solve(TearingProblem(ffs; critical_resonant_field=CriticalResonantFieldControl()), SLAYERModel(); config=TearingConfig(...))` runs the same stage, and `solve_inner(GGJModel(), p, γ)` evaluates one GGJ layer at a growth rate ``\gamma``.
 GGJ's η and ρ default from the kinetic data (η by the resistivity closure SLAYER uses, ρ = μ·m_p·n_e); pass `eta` and `rho` to `TearingProblem` to override them.
 
@@ -103,6 +131,9 @@ The DIII-D-like kinetic file of the examples carries analytic χ profiles of the
 | `Tearing/Diagnostics/ValidRoots`, `Poles`, `FilteredRoots` | the root classification of the scan |
 | `Tearing/Scan/Surface_<k>/` | the scanned Δ(Q) with `store_scan = true` |
 | `Tearing/CriticalResonantField/br_crit`, `br_crit_tesla`, `Q_peak`, `Q0` | the penetration threshold as ``b_r/B_\phi`` and in tesla, the Q of the balance maximum and the natural rotation |
+| `Tearing/IslandFlattening/W`, `Delta_W`, `Delta_W_over_W`, `Delta_p0` | the stored energy, the part the islands remove and its fraction, and the on-axis pressure drop |
+| `Tearing/IslandFlattening/Islands/`, `Intervals/` | each island (`rational_m`, `rational_n`, `rational_psi`, its clipped interval) and the merged flattened intervals with the pressure drop `Delta_p` across each |
+| `Tearing/IslandFlattening/Profiles/psi`, `p`, `p_mod` | the equilibrium and flattened pressure (Pa), with `jbs_B` and `jbs_B_mod` when a bootstrap profile was given |
 
 ## Numerical notes
 
@@ -175,3 +206,9 @@ Modules = [GeneralizedPerturbedEquilibrium.Tearing.Dispersion]
 ```@autodocs
 Modules = [GeneralizedPerturbedEquilibrium.Tearing.CriticalResonantField]
 ```
+
+### IslandFlattening
+
+```@autodocs
+Modules = [GeneralizedPerturbedEquilibrium.Tearing.IslandFlattening]
+```

docs/src/workflow.md

@@ -161,7 +161,7 @@ The keys each section accepts are the fields of its control struct, listed in th
 | `[ForcingTerms]`, `[[ForcingTerms.coil_set]]` | the external field: a file or coil sets | [Forcing Terms](forcing_terms.md) |
 | `[PerturbedEquilibrium]` | the plasma response stage and its singular-coupling options | [Perturbed Equilibrium](perturbed_equilibrium.md) |
 | `[KineticForces]` | the NTV torque stage and the kinetic solve of `ForceFreeStates` | [Kinetic Forces](kinetic_forces.md) |
-| `[Tearing]`, `[Tearing.search]`, `[Tearing.CriticalResonantField]`, `[SLAYER]` | the tearing stage, its root search, the critical field and the SLAYER layer model | [Tearing](inner_layer.md) |
+| `[Tearing]`, `[Tearing.search]`, `[Tearing.CriticalResonantField]`, `[Tearing.IslandFlattening]`, `[SLAYER]` | the tearing stage, its root search, the critical field, the island pressure flattening and the SLAYER layer model | [Tearing](inner_layer.md) |
 | `[ErrorFields]`, `[ErrorFields.MonteCarlo]`, `[ErrorFields.scenario]`, `[ErrorFields.Risk]`, `[ErrorFields.NTV]` | the error-field tolerance stage | [Error Fields](error_fields.md) |
 
 ## Running it
@@ -216,7 +216,7 @@ All results are written to one HDF5 file (`[Run] output_file`, default `gpec.h5`
 | `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/` |
 | `KineticForces/` | `<method>/` torque and energy profiles; `PerSpecies/<species>/<method>/` for multi-ion runs |
 | `ErrorFields/` | `CoilSensitivities/`, `MonteCarlo/`, `Risk/`, `NTV/` |
-| `Tearing/` | `PerSurface/`, `Roots/`, `LayerWidths/`, `Diagnostics/`, `CriticalResonantField/` |
+| `Tearing/` | `PerSurface/`, `Roots/`, `LayerWidths/`, `Diagnostics/`, `CriticalResonantField/`, `IslandFlattening/` |
 | `SurfaceGeometries/` | Plasma and wall surface point clouds for visualization |
 
 ### Replaying a run

regression-harness/cases/diiid_n1_island_flattening.toml

@@ -0,0 +1,68 @@
+# Regression case: DIII-D-like n=1 with static pressure flattening across prescribed islands.
+# Same deck as diiid_n1 plus [Tearing.IslandFlattening] on the 3/2, 2/1 and 3/1 surfaces, so the
+# merged intervals, the stored energy they remove and the flattened pressure are tracked.
+# Each [quantities.*] block names an HDF5 path in the run output, how to extract it,
+# and the noise floor below which a difference is treated as zero.
+[case]
+name = "diiid_n1_island_flattening"
+description = "DIII-D-like equilibrium, n=1, pressure flattened across 3/2, 2/1 and 3/1 islands"
+example_dir = "examples/DIIID-like_ideal_example"
+
+[overrides]
+"Tearing.IslandFlattening.widths" = { "3/2" = 0.04, "2/1" = 0.06, "3/1" = 0.03 }
+
+[quantities.W]
+h5path = "Tearing/IslandFlattening/W"
+type = "real_scalar"
+extract = "value"
+label = "stored energy W"
+noise_threshold = 1e-6
+order = 10
+
+[quantities.delta_W]
+h5path = "Tearing/IslandFlattening/Delta_W"
+type = "real_scalar"
+extract = "value"
+label = "stored energy removed"
+noise_threshold = 1e-6
+order = 11
+
+[quantities.delta_W_over_W]
+h5path = "Tearing/IslandFlattening/Delta_W_over_W"
+type = "real_scalar"
+extract = "value"
+label = "fraction of W removed"
+noise_threshold = 1e-12
+order = 12
+
+[quantities.delta_p0]
+h5path = "Tearing/IslandFlattening/Delta_p0"
+type = "real_scalar"
+extract = "value"
+label = "on-axis pressure drop"
+noise_threshold = 1e-6
+order = 13
+
+[quantities.island_psi]
+h5path = "Tearing/IslandFlattening/Islands/rational_psi"
+type = "real_vector"
+extract = "all_real"
+label = "island rational surfaces"
+noise_threshold = 1e-12
+order = 20
+
+[quantities.interval_delta_p]
+h5path = "Tearing/IslandFlattening/Intervals/Delta_p"
+type = "real_vector"
+extract = "all_real"
+label = "pressure drop per merged interval"
+noise_threshold = 1e-6
+order = 21
+
+[quantities.p_mod]
+h5path = "Tearing/IslandFlattening/Profiles/p_mod"
+type = "real_vector"
+extract = "checksum"
+label = "flattened pressure (checksum)"
+noise_threshold = 0
+order = 30

src/Deck.jl

@@ -360,7 +360,9 @@ end
 """
 function tearing_from_deck(inputs::AbstractDict)
     haskey(inputs, "Tearing") || return nothing
-    table = inputs["Tearing"]
+    # A [Tearing] holding only [Tearing.IslandFlattening] runs the flattening without a layer solve.
+    table = _without(inputs["Tearing"], "IslandFlattening")
+    isempty(table) && haskey(inputs["Tearing"], "IslandFlattening") && return nothing
     get(table, "omega_E_kHz", Dict{String,Any}()) isa AbstractDict ||
         throw(ArgumentError("[Tearing] omega_E_kHz must be a table keyed by m/n, e.g. omega_E_kHz = { \"2/1\" = 3.0 }"))
     fields = _struct_fields(TearingProblem, (:ffs, :critical_resonant_field))
@@ -377,6 +379,21 @@ function tearing_from_deck(inputs::AbstractDict)
     return (; problem, model=inner_model == "SLAYER" ? slayer_model(inputs) : ggj_model(inputs), config)
 end
 
+"""
+    island_flattening_from_deck(inputs) -> Union{Nothing,Dict{String,Float64}}
+
+The full island widths (ψ_N) keyed by `"m/n"` of `[Tearing.IslandFlattening] widths`; `nothing` without that section.
+"""
+function island_flattening_from_deck(inputs::AbstractDict)
+    table = get(get(inputs, "Tearing", Dict{String,Any}()), "IslandFlattening", nothing)
+    table === nothing && return nothing
+    kwargs = section_keywords(table, "Tearing.IslandFlattening", Dict{Symbol,Any}(:widths => AbstractDict))
+    widths = get(kwargs, :widths, nothing)
+    (widths isa AbstractDict && !isempty(widths) && all(v -> v isa Real, values(widths))) ||
+        throw(ArgumentError("[Tearing.IslandFlattening] widths must be a table of full widths in ψ_N keyed by m/n, e.g. widths = { \"2/1\" = 0.05 }"))
+    return Dict{String,Float64}(String(k) => Float64(v) for (k, v) in widths)
+end
+
 """
     perturbed_equilibrium_from_deck(inputs) -> Union{Nothing,NamedTuple}
 
@@ -471,6 +488,7 @@ function parse_deck(inputs::AbstractDict, path::AbstractString)
             error("[Match] needs Δ′: use integrator = \"Galerkin\" with [Galerkin] rpec_flag = true, or \"Riccati\" with vac_flag = true")
     end
     tearing = tearing_from_deck(inputs)
+    island_flattening = island_flattening_from_deck(inputs)
     perturbed_equilibrium = perturbed_equilibrium_from_deck(inputs)
     forcing = forcing_terms_control(inputs, path)
     if haskey(inputs, "ErrorFields")
@@ -483,7 +501,7 @@ function parse_deck(inputs::AbstractDict, path::AbstractString)
         error("[GGJ] is set but neither [Match] nor [Tearing] inner_model = \"GGJ\" uses it")
     haskey(inputs, "SLAYER") && tearing_model != "SLAYER" && !(has_ffs && get(inputs["ForceFreeStates"], "psilim_from_layer_overlap", false) === true) &&
         error("[SLAYER] is set but neither [Tearing] inner_model = \"SLAYER\" nor [ForceFreeStates] psilim_from_layer_overlap = true uses it")
-    return (; run=run_settings(inputs, path), force_balance, local_stability=local_stability_from_deck(inputs), match, tearing,
+    return (; run=run_settings(inputs, path), force_balance, local_stability=local_stability_from_deck(inputs), match, tearing, island_flattening,
         kinetic_forces=haskey(inputs, "KineticForces") ? kinetic_forces_control(inputs) : nothing, perturbed_equilibrium, forcing,
         error_fields=error_fields_from_deck(inputs, path))
 end
@@ -578,6 +596,10 @@ function deck_inputs(deck::NamedTuple, config::Equilibrium.EquilibriumConfig, an
         merge!(inputs, _tearing_sections(filter(p -> p.first != :critical_resonant_field, problem), get(problem, :critical_resonant_field, nothing),
             model, config))
     end
+    if deck.island_flattening !== nothing
+        tearing = get!(inputs, "Tearing", Dict{String,Any}())
+        tearing["IslandFlattening"] = Dict{String,Any}("widths" => Dict{String,Any}(deck.island_flattening))
+    end
     deck.kinetic_forces === nothing || merge!(inputs, deck_sections(deck.kinetic_forces))
     pe = deck.perturbed_equilibrium
     pe === nothing || (inputs["PerturbedEquilibrium"] = _deck_table(merge(pe.response, pe.drive)))
@@ -626,6 +648,8 @@ deck_sections(m::ForceFreeStates.MatchResult) =
     merge(deck_sections(m.model), Dict{String,Any}("Match" => merge(_deck_table(m.problem, (:ffs,)), Dict{String,Any}("inner_model" => "GGJ"))))
 deck_sections(kf::KineticForcesResult) = deck_sections(kf.control)
 deck_sections(t::TearingResult) = t.model === nothing ? Dict{String,Any}() : _tearing_sections(t.problem, t.problem.critical_resonant_field, t.model, t.config)
+deck_sections(r::IslandFlatteningResult) =
+    Dict{String,Any}("Tearing" => Dict{String,Any}("IslandFlattening" => Dict{String,Any}("widths" => Dict{String,Any}(r.widths))))
 deck_sections(ef::ErrorFieldsResult) = _error_fields_sections(ef.settings, ef.tolerances === nothing ? "" : _RECORDED_TOLERANCES)
 
 function deck_sections(ffs::ForceFreeStatesResult)

src/GeneralizedPerturbedEquilibrium.jl

@@ -74,7 +74,7 @@ include("Rerun.jl")
 using .ForceFreeStates: AbstractIntegrator, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, ForceFreeStatesResult
 using .ForceFreeStates: MatchProblem
 using .InnerLayer: GGJModel, SLAYERModel
-using .Tearing: TearingProblem, TearingConfig, TearingResult, CriticalResonantFieldControl
+using .Tearing: TearingProblem, TearingConfig, TearingResult, CriticalResonantFieldControl, IslandFlatteningResult, island_flattening
 using .LocalStability: LocalStabilityProblem, Ballooning, LocalStabilityResult
 using .Equilibrium: PlasmaEquilibrium, KineticProfileData, IonSpecies, attach_kinetic_profiles!
 using .ForcingTerms: RMPField
@@ -140,7 +140,7 @@ equilibrium and writes nothing. The raw forcing file contents are captured in
 `Input/RawInputs/ForcingTerms/`, and `preloaded_coil_sets` lets the rerun path inject the coil
 geometry stored in `Input/RawInputs/Coils/`.
 
-Returns `(; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing)`, each
+Returns `(; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing, island_flattening)`, each
 `nothing` when its section is absent. `equil` is the equilibrium the stages used (the solve's, on
 the two-pass grid), `ffs` the result the later stages read (the matched one when the deck has
 `[Match]`), and `match` its `MatchResult`.
@@ -158,7 +158,7 @@ function main_from_inputs(
     runtimes = Vector{Pair{String,Float64}}()
     deck = parse_deck(inputs, path)
     (; output_file, verbose) = deck.run
-    local_stability = ffs = match = resp = pe = kinetic = error_fields = tearing = nothing
+    local_stability = ffs = match = resp = pe = kinetic = error_fields = tearing = island_flat = nothing
 
     equil = _timed_stage(runtimes, "equilibrium", "Equilibrium") do
         Equilibrium.PlasmaEquilibrium(eq_config, additional_input)
@@ -224,12 +224,20 @@ function main_from_inputs(
         save_to_hdf5(output_file, tearing)
     end
 
+    # After the tearing save, which rewrites Tearing/.
+    if deck.island_flattening !== nothing
+        island_flat = _timed_stage(runtimes, "island_flattening", "Island Flattening") do
+            island_flattening(equil; widths=deck.island_flattening)
+        end
+        save_to_hdf5(output_file, island_flat)
+    end
+
     total_dt = time() - total_start
     push!(runtimes, "total" => total_dt)
     ffs === nothing || _write_runtimes!(output_file, runtimes)
     @info "\n$_BANNER\n  GPEC completed successfully in $(@sprintf("%.3f", total_dt)) s\n$_BANNER"
 
-    return (; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing)
+    return (; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing, island_flattening=island_flat)
 end
 
 """
@@ -306,6 +314,7 @@ export solve, plasma_response, perturbed_equilibrium, kinetic_torque, error_fiel
 export KineticProfileData, IonSpecies, attach_kinetic_profiles!
 export PlasmaEquilibrium, ForceBalanceProblem, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, kinetic_matrices, ForceFreeStatesResult, RMPField, PlasmaResponse
 export MatchProblem, TearingProblem, TearingConfig, TearingResult, GGJModel, SLAYERModel
+export IslandFlatteningResult, island_flattening
 export LocalStabilityProblem, Ballooning, LocalStabilityResult
 export KineticForcesControl, KineticForcesResult, ErrorFieldsResult