From de6d76781f3ab7f86cb057e2efbcf0ddd6e2e023 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Fri, 13 Nov 2026 22:14:41 -0500 Subject: [PATCH 1/3] Tearing.IslandFlattening - FEATURE - Flatten the pressure across islands of prescribed width and report the stored energy they remove New Tearing.IslandFlattening: dp/dpsi = 0 across islands of given full width at every q = m/n surface, the union flattened once and the pressure anchored at the boundary, in closed form. Reports W, delta_W, the pressure drop across each merged interval and on axis, and, given a profile, the bootstrap fraction removed. Deck section [Tearing.IslandFlattening] widths runs without a layer solve; output Tearing/IslandFlattening/. record_inputs! keeps the rest of a section when a result states only its sub-tables. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae --- docs/src/inner_layer.md | 34 ++++ docs/src/workflow.md | 4 +- .../cases/diiid_n1_island_flattening.toml | 68 +++++++ src/Deck.jl | 28 ++- src/GeneralizedPerturbedEquilibrium.jl | 17 +- src/Tearing/HDF5Output.jl | 67 ++++++- src/Tearing/IslandFlattening/Flattening.jl | 166 ++++++++++++++++++ .../IslandFlattening/IslandFlattening.jl | 15 ++ src/Tearing/Tearing.jl | 8 +- src/Utilities/Output.jl | 13 +- test/runtests_island_flattening.jl | 91 ++++++++++ test/runtests_solve_api.jl | 21 +++ test/suites.jl | 1 + 13 files changed, 519 insertions(+), 14 deletions(-) create mode 100644 regression-harness/cases/diiid_n1_island_flattening.toml create mode 100644 src/Tearing/IslandFlattening/Flattening.jl create mode 100644 src/Tearing/IslandFlattening/IslandFlattening.jl create mode 100644 test/runtests_island_flattening.jl diff --git a/docs/src/inner_layer.md b/docs/src/inner_layer.md index e67b984ffd..09b9aa4f9f 100644 --- a/docs/src/inner_layer.md +++ b/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,22 @@ 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 ``\langle j_\mathrm{bs}\cdot B\rangle(\psi)`` (the `jbs_B` keyword of `island_flattening`), the result also carries the bootstrap current removed: to leading order ``j_\mathrm{bs} \propto p'``, so it vanishes on the flattened intervals and is unchanged elsewhere. +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 +95,14 @@ 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 +128,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_/` | 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 (m, n, ψ_N, w, its clipped interval) and the merged flattened intervals with the pressure drop 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 +203,9 @@ Modules = [GeneralizedPerturbedEquilibrium.Tearing.Dispersion] ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.Tearing.CriticalResonantField] ``` + +### IslandFlattening + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Tearing.IslandFlattening] +``` diff --git a/docs/src/workflow.md b/docs/src/workflow.md index 38b4cdd8ba..fe9c911324 100644 --- a/docs/src/workflow.md +++ b/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/` | `/` torque and energy profiles; `PerSpecies///` 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 diff --git a/regression-harness/cases/diiid_n1_island_flattening.toml b/regression-harness/cases/diiid_n1_island_flattening.toml new file mode 100644 index 0000000000..e87d10be10 --- /dev/null +++ b/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/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 diff --git a/src/Deck.jl b/src/Deck.jl index ea0f4b17ce..ca5a0b48e3 100644 --- a/src/Deck.jl +++ b/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) && 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) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index bda2b7c8d8..1a8dc16d53 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/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 replaces 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 diff --git a/src/Tearing/HDF5Output.jl b/src/Tearing/HDF5Output.jl index dbf1ed5b9f..a69c273c07 100644 --- a/src/Tearing/HDF5Output.jl +++ b/src/Tearing/HDF5Output.jl @@ -1,9 +1,9 @@ import ..Utilities: save_to_hdf5 -import ..Utilities.HDF5Annotations: replace_group!, write_ragged! +import ..Utilities.HDF5Annotations: replace_group!, ensure_group, write_ragged! # HDF5Output.jl # -# Writer for a `TearingResult`. Layout under the caller's parent group (inputs live only under Input/): +# Writers for a `TearingResult` and an `IslandFlatteningResult`. Layout under the caller's parent group (inputs live only under Input/): # # Tearing/ # ├── PerSurface/ -- layer parameters (SLAYER or GGJ set), Δ' matrix, omega_E, q_shift @@ -12,7 +12,8 @@ import ..Utilities.HDF5Annotations: replace_group!, write_ragged! # ├── Diagnostics/ -- ValidRoots, Poles, FilteredRoots (+ *Winding), flat-plus-offsets ragged; # │ uncoupled entries are in each surface's E×B frame # ├── Scan/ -- optional: full Q/Δ scan data -# └── CriticalResonantField/ -- optional: torque-balance critical b_r/B_φ per surface +# ├── CriticalResonantField/ -- optional: torque-balance critical b_r/B_φ per surface +# └── IslandFlattening/ -- its own result; a TearingResult save replaces Tearing/, so it is written after using HDF5 @@ -366,3 +367,63 @@ function _write_single_scan!(g, data::AMRResult) g["truncated"] = Int(data.truncated) return nothing end + +# ---------- island flattening ---------- +const ISLAND_FLATTENING_H5_ANNOTATIONS = [ + "Islands/m" => (; long_name="poloidal mode number m of each island", dims=("island",)), + "Islands/n" => (; long_name="toroidal mode number n of each island", dims=("island",)), + "Islands/psi" => (; long_name="rational surface ψ_N of each island", dims=("island",)), + "Islands/w" => (; long_name="full island width in ψ_N", dims=("island",)), + "Islands/psi_in" => (; long_name="inner end of the flattened interval (ψ_N, clipped to the profile grid)", dims=("island",)), + "Islands/psi_out" => (; long_name="outer end of the flattened interval (ψ_N, clipped to the profile grid)", dims=("island",)), + "Islands/clipped" => (; long_name="flag: the island reaches past an end of the profile grid", dims=("island",)), + "Intervals/psi_in" => (; long_name="inner end of each merged flattened interval (ψ_N)", dims=("interval",)), + "Intervals/psi_out" => (; long_name="outer end of each merged flattened interval (ψ_N)", dims=("interval",)), + "Intervals/delta_p" => (; long_name="pressure drop p(ψ_in) − p(ψ_out) across each merged interval", units="Pa", dims=("interval",)), + "Profiles/psi" => (; long_name="normalized poloidal flux ψ_N (profile grid with the interval ends added)", dims=("psi",)), + "Profiles/p" => (; long_name="equilibrium pressure", units="Pa", dims=("psi",)), + "Profiles/p_mod" => (; long_name="pressure flattened across the islands, anchored at the boundary", units="Pa", dims=("psi",)), + "Profiles/jbs_B" => (; long_name="given bootstrap ⟨j_bs·B⟩", units="A*T/m^2", dims=("psi",)), + "Profiles/jbs_B_mod" => (; long_name="bootstrap ⟨j_bs·B⟩ with the flattened intervals removed (j_bs ∝ p′ to leading order)", units="A*T/m^2", dims=("psi",)), + "W" => (; long_name="stored energy (3/2)∫p dV", units="J"), + "delta_W" => (; long_name="stored energy removed by the flattening (3/2)∫(p − p_mod) dV", units="J"), + "delta_W_over_W" => (; long_name="fraction of the stored energy removed"), + "delta_p0" => (; long_name="on-axis pressure drop", units="Pa"), + "jbs_fraction_lost" => (; long_name="volume-integrated fraction of ⟨j_bs·B⟩ removed (NaN without a bootstrap profile)") +] + +""" + save_to_hdf5(parent, result::IslandFlatteningResult) + +Write `result` into `Tearing/IslandFlattening/` of `parent`, replacing any earlier one and +keeping the rest of `Tearing/`. +""" +function save_to_hdf5(parent::Union{HDF5.File,HDF5.Group}, r::IslandFlatteningResult) + g = replace_group!(ensure_group(parent, "Tearing"), "IslandFlattening") + isl = r.islands + g["Islands/m"] = [i.m for i in isl] + g["Islands/n"] = [i.n for i in isl] + g["Islands/psi"] = [i.psi for i in isl] + g["Islands/w"] = [i.w for i in isl] + g["Islands/psi_in"] = [i.psi_in for i in isl] + g["Islands/psi_out"] = [i.psi_out for i in isl] + g["Islands/clipped"] = Int[i.clipped for i in isl] + g["Intervals/psi_in"] = first.(r.intervals) + g["Intervals/psi_out"] = last.(r.intervals) + g["Intervals/delta_p"] = r.delta_p + g["Profiles/psi"] = r.psi + g["Profiles/p"] = r.p + g["Profiles/p_mod"] = r.p_mod + if !isempty(r.jbs_B) + g["Profiles/jbs_B"] = r.jbs_B + g["Profiles/jbs_B_mod"] = r.jbs_B_mod + end + g["W"] = r.W + g["delta_W"] = r.delta_W + g["delta_W_over_W"] = r.delta_W / r.W + g["delta_p0"] = r.delta_p0 + g["jbs_fraction_lost"] = r.jbs_fraction_lost + Utilities.HDF5Annotations.annotate!(g, ISLAND_FLATTENING_H5_ANNOTATIONS) + Utilities.record_inputs!(parent, r) + return g +end diff --git a/src/Tearing/IslandFlattening/Flattening.jl b/src/Tearing/IslandFlattening/Flattening.jl new file mode 100644 index 0000000000..38ae221eb4 --- /dev/null +++ b/src/Tearing/IslandFlattening/Flattening.jl @@ -0,0 +1,166 @@ +# Flattening.jl +# +# The pressure gradient vanishes inside each island and is unchanged outside it; the +# pressure keeps its boundary value at the last profile point and is integrated inward, +# as in the Fortran GPEC saturated-island profile. Islands that overlap flatten their +# union once. With the flattened intervals [a_k, b_k] merged, the integral is exact: +# +# p_mod(ψ) = p(ψ) + Σ_{k: b_k > ψ} [p(b_k) − p(max(ψ, a_k))], +# +# so every point inward of an interval drops by p(a_k) − p(b_k), and points inside it sit +# at p(b_k). The stored energy removed is ΔW = (3/2) ∫ (p − p_mod) dV. + +""" + FlattenedIsland + +One island of the set: mode numbers `m`, `n`, rational surface `psi` (ψ_N), full width `w` +(ψ_N) and the flattened interval `psi_in ≤ ψ ≤ psi_out`, clipped to the profile grid. +`clipped` is true when the island reaches past either end of the grid. +""" +struct FlattenedIsland + m::Int + n::Int + psi::Float64 + w::Float64 + psi_in::Float64 + psi_out::Float64 + clipped::Bool +end + +""" + IslandFlatteningResult + +Static flattening of the equilibrium pressure across `islands`. + + - `widths`: the requested full widths keyed by `"m/n"` (ψ_N); every rational surface of each + `m/n` gets the width + - `islands`: one `FlattenedIsland` per rational surface found + - `intervals`: the flattened intervals, merged where islands overlap, inward to outward + - `delta_p`: pressure drop across each merged interval, p(a_k) − p(b_k) (Pa); for an edge + interval it is the pedestal-top pressure change + - `psi`, `p`, `p_mod`: the profile grid with the interval ends added, the equilibrium and + flattened pressures on it (Pa) + - `W`, `delta_W`: stored energy (3/2)∫p dV and the part the islands remove (J) + - `delta_p0`: on-axis pressure drop, the sum of `delta_p` (Pa) + - `jbs_B`, `jbs_B_mod`, `jbs_fraction_lost`: the given bootstrap ⟨j_bs·B⟩ on `psi`, the + same with the flattened intervals removed, and the volume-integrated fraction removed; + empty and NaN when no profile was given. To leading order j_bs ∝ p′, so it vanishes + where p′ does and is unchanged elsewhere. +""" +struct IslandFlatteningResult + widths::Dict{String,Float64} + islands::Vector{FlattenedIsland} + intervals::Vector{Tuple{Float64,Float64}} + delta_p::Vector{Float64} + psi::Vector{Float64} + p::Vector{Float64} + p_mod::Vector{Float64} + W::Float64 + delta_W::Float64 + delta_p0::Float64 + jbs_B::Vector{Float64} + jbs_B_mod::Vector{Float64} + jbs_fraction_lost::Float64 +end + +""" + merged_intervals(islands) -> Vector{Tuple{Float64,Float64}} + +The union of the islands' flattened intervals as disjoint intervals, inward to outward. +Intervals that touch or overlap merge. +""" +function merged_intervals(islands::AbstractVector{FlattenedIsland}) + spans = sort!([(i.psi_in, i.psi_out) for i in islands if i.psi_out > i.psi_in]) + merged = Tuple{Float64,Float64}[] + for (a, b) in spans + if !isempty(merged) && a <= merged[end][2] + merged[end] = (merged[end][1], max(merged[end][2], b)) + else + push!(merged, (a, b)) + end + end + return merged +end + +""" + flattened_pressure(p, intervals, psi) -> Float64 + +The flattened pressure at `psi` for a pressure profile `p(ψ)` and disjoint flattened +`intervals`, anchored at the boundary value outside every interval. +""" +function flattened_pressure(p, intervals, psi::Real) + pm = p(psi) + for (a, b) in intervals + b > psi && (pm += p(b) - p(max(psi, a))) + end + return pm +end + +""" + island_flattening(equil; widths, jbs_B=nothing) -> IslandFlatteningResult + +Flatten the pressure of `equil` across an island of full width `widths["m/n"]` (ψ_N) at +every rational surface q = m/n. `jbs_B`, a function of ψ_N giving ⟨j_bs·B⟩, adds the +bootstrap current the flattening removes. +""" +function island_flattening(equil::PlasmaEquilibrium; widths::AbstractDict, jbs_B=nothing) + prof = equil.profiles + lo, hi = prof.xs[1], prof.xs[end] + p(psi) = prof.P_spline(psi) / mu0 + dV(psi) = prof.dVdpsi_spline(psi) + + requested = Dict{String,Float64}() + islands = FlattenedIsland[] + for (label, w) in widths + m, n = _mode_numbers(label) + w >= 0 || throw(ArgumentError("island width for $label must be non-negative, got $w")) + requested[label] = w + found = [s.psifac for s in _find_rational_surfaces(equil, n, n) if s.m == m] + isempty(found) && @warn "No q = $label surface in the equilibrium; island skipped" + for psi in found + a, b = psi - w / 2, psi + w / 2 + push!(islands, FlattenedIsland(m, n, psi, w, max(a, lo), min(b, hi), a < lo || b > hi)) + end + end + sort!(islands; by=i -> i.psi) + intervals = merged_intervals(islands) + delta_p = [p(a) - p(b) for (a, b) in intervals] + + grid = sort!(unique!(vcat(prof.xs, (x for iv in intervals for x in iv)...))) + p_eq = p.(grid) + p_mod = [flattened_pressure(p, intervals, x) for x in grid] + + # Exact for the piecewise-cubic integrands: nodes on every knot and interval end. + x, wq = gausslegendre(4) + segment(f, a, b) = (b - a) / 2 * sum(wq[j] * f((b - a) / 2 * x[j] + (a + b) / 2) for j in eachindex(x)) + integrate(f) = sum(segment(f, grid[i], grid[i+1]) for i in 1:(length(grid)-1)) + W = 1.5 * integrate(psi -> p(psi) * dV(psi)) + delta_W = 1.5 * integrate(psi -> (p(psi) - flattened_pressure(p, intervals, psi)) * dV(psi)) + + inside(psi) = any(a < psi < b for (a, b) in intervals) + if jbs_B === nothing + jbs, jbs_mod, lost = Float64[], Float64[], NaN + else + jbs = Float64.(jbs_B.(grid)) + jbs_mod = [inside(xi) ? 0.0 : j for (xi, j) in zip(grid, jbs)] + total = integrate(psi -> jbs_B(psi) * dV(psi)) + lost = sum(integrate_between(psi -> jbs_B(psi) * dV(psi), grid, a, b, segment) for (a, b) in intervals; init=0.0) / total + end + return IslandFlatteningResult(requested, islands, intervals, delta_p, grid, p_eq, p_mod, W, delta_W, sum(delta_p; init=0.0), + jbs, jbs_mod, lost) +end + +# Integral of f over [a, b], whose ends are nodes of `grid`. +function integrate_between(f, grid, a, b, segment) + i0, i1 = searchsortedfirst(grid, a), searchsortedfirst(grid, b) + return sum(segment(f, grid[i], grid[i+1]) for i in i0:(i1-1); init=0.0) +end + +function _mode_numbers(label::AbstractString) + parts = split(label, "/") + ok = length(parts) == 2 && all(s -> tryparse(Int, strip(s)) !== nothing, parts) + ok || throw(ArgumentError("island widths are keyed by \"m/n\", e.g. \"2/1\"; got \"$label\"")) + m, n = parse.(Int, strip.(parts)) + (m > 0 && n > 0) || throw(ArgumentError("island \"$label\" needs positive m and n")) + return m, n +end diff --git a/src/Tearing/IslandFlattening/IslandFlattening.jl b/src/Tearing/IslandFlattening/IslandFlattening.jl new file mode 100644 index 0000000000..9e9e3a2949 --- /dev/null +++ b/src/Tearing/IslandFlattening/IslandFlattening.jl @@ -0,0 +1,15 @@ +# Static pressure flattening across a set of magnetic islands of given full width, with the +# stored energy and pressure each set of islands removes. No time evolution. + +module IslandFlattening + +using FastGaussQuadrature: gausslegendre + +using ...Equilibrium: PlasmaEquilibrium, mu0 +using ...ForceFreeStates: _find_rational_surfaces + +include("Flattening.jl") + +export FlattenedIsland, IslandFlatteningResult, island_flattening, merged_intervals, flattened_pressure + +end # module IslandFlattening diff --git a/src/Tearing/Tearing.jl b/src/Tearing/Tearing.jl index 0055cb2739..d597b7ee8e 100644 --- a/src/Tearing/Tearing.jl +++ b/src/Tearing/Tearing.jl @@ -6,6 +6,8 @@ # extraction (consumes any InnerLayerModel) # CriticalResonantField -- real-Q torque-balance scan for the critical # resonant field of error-field penetration +# IslandFlattening -- static pressure flattening across islands of given +# width, and the stored energy it removes # TearingProblem / TearingConfig / TearingResult -- the public solve, # `solve(TearingProblem(ffs; ...), model; config)`, with one # method per inner-layer model (SLAYERModel, GGJModel) @@ -39,9 +41,11 @@ using .InnerLayer.SLAYER: _slayer_profiles, _layer_zeff, _layer_ion_mass, _layer include("LayerOverlap.jl") include("Dispersion/Dispersion.jl") include("CriticalResonantField/CriticalResonantField.jl") +include("IslandFlattening/IslandFlattening.jl") import .Dispersion as Dispersion import .CriticalResonantField as CriticalResonantField +import .IslandFlattening as IslandFlattening using .Dispersion: SurfaceCoupling, surface_coupling, MultiSurfaceCoupling, multi_surface_coupling, ScanResult, brute_force_scan, @@ -49,6 +53,7 @@ using .Dispersion: SurfaceCoupling, surface_coupling, MultiBoxAMRResult, multi_box_amr_scan, as_amr_result, GrowthRateResult, find_growth_rates using .CriticalResonantField: TorqueBalance, torque_balance_scan +using .IslandFlattening: IslandFlatteningResult, island_flattening include("Config.jl") include("Result.jl") @@ -57,12 +62,13 @@ include("SLAYERTearing.jl") include("GGJTearing.jl") include("HDF5Output.jl") -export InnerLayer, Dispersion, CriticalResonantField +export InnerLayer, Dispersion, CriticalResonantField, IslandFlattening export resistive_layer_overlap, layer_overlap_scan, LayerOverlapScan export TearingProblem, TearingConfig, TearingResult export run_slayer_from_inputs export save_to_hdf5 export CriticalResonantFieldControl, CriticalResonantFieldScan, CriticalResonantFieldResult export run_critical_resonant_field +export IslandFlatteningResult, island_flattening end # module Tearing diff --git a/src/Utilities/Output.jl b/src/Utilities/Output.jl index bf75c28f85..ce30956f19 100644 --- a/src/Utilities/Output.jl +++ b/src/Utilities/Output.jl @@ -29,13 +29,22 @@ function deck_sections end record_inputs!(target, result) Merge the deck sections of `result` into the run record `Input/gpec_toml_raw` of `target`, replacing any section of -the same name, so every saved result adds the deck it was computed from. +the same name, so every saved result adds the deck it was computed from. A result that states only sub-tables of a +section (such as `[Tearing.IslandFlattening]`) replaces those sub-tables and keeps the rest of the section. """ function record_inputs!(target::Union{HDF5.File,HDF5.Group}, result) path = "Input/gpec_toml_raw" deck = haskey(target, path) ? TOML.parse(read(target[path])) : Dict{String,Any}() - merge!(deck, deck_sections(result)) + for (name, section) in deck_sections(result) + if _only_subtables(section) && get(deck, name, nothing) isa AbstractDict + merge!(deck[name], section) + else + deck[name] = section + end + end haskey(target, path) && HDF5.delete_object(target, path) target[path] = sprint(TOML.print, deck) return nothing end + +_only_subtables(section) = section isa AbstractDict && !isempty(section) && all(v -> v isa AbstractDict, values(section)) diff --git a/test/runtests_island_flattening.jl b/test/runtests_island_flattening.jl new file mode 100644 index 0000000000..42c2cda8f4 --- /dev/null +++ b/test/runtests_island_flattening.jl @@ -0,0 +1,91 @@ +using TOML +using HDF5 + +# Static pressure flattening across islands of given width, on the coarse TJ fixture equilibrium. +@testset "island flattening" begin + GPEC = GeneralizedPerturbedEquilibrium + IF = GPEC.Tearing.IslandFlattening + template = joinpath(@__DIR__, "test_data", "regression_tj_ideal_example") + deck = TOML.parsefile(joinpath(template, "gpec.toml")) + equil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(deck, template), GPEC.analytic_config(deck)) + prof = equil.profiles + p(psi) = prof.P_spline(psi) / GPEC.Equilibrium.mu0 + dV(psi) = prof.dVdpsi_spline(psi) + one = island_flattening(equil; widths=Dict("2/1" => 0.05)) + psi_s = only(one.islands).psi + # Volume inside ψ_N = x, by the midpoint rule. + volume_inside(x; n=4000) = sum(dV(prof.xs[1] + (x - prof.xs[1]) * (k - 0.5) / n) for k in 1:n) * (x - prof.xs[1]) / n + + @testset "the Fortran trapezoid rule converges to the closed form" begin + # Fortran GPEC: zero dp/dψ on the island, anchor at the edge, integrate inward on a grid. + grid = range(prof.xs[1], prof.xs[end]; length=200_001) + dp = [abs(x - psi_s) <= 0.025 ? 0.0 : prof.P_deriv(x) / GPEC.Equilibrium.mu0 for x in grid] + pm = similar(dp) + pm[end] = p(grid[end]) + for i in (length(grid)-1):-1:1 + pm[i] = pm[i+1] - 0.5 * (dp[i+1] + dp[i]) * step(grid) + end + @test pm[1] ≈ IF.flattened_pressure(p, one.intervals, grid[1]) rtol = 1e-6 + dW = 1.5 * sum((p(x) - y) * dV(x) for (x, y) in zip(grid, pm)) * step(grid) + @test dW ≈ one.delta_W rtol = 1e-3 + end + + @testset "narrow islands: ΔW → 0 and ΔW ∝ w" begin + @test island_flattening(equil; widths=Dict("2/1" => 0.0)).delta_W == 0 + small = [island_flattening(equil; widths=Dict("2/1" => w)).delta_W for w in (1e-3, 2e-3)] + @test small[2] / small[1] ≈ 2 rtol = 1e-2 + # Leading order: the drop −p′w across the island lowers everything inward of it. + @test small[1] ≈ 1.5 * (-prof.P_deriv(psi_s) / GPEC.Equilibrium.mu0) * 1e-3 * volume_inside(psi_s) rtol = 2e-2 + end + + @testset "overlapping islands are flattened once" begin + # The 2/1 (ψ_N ≈ 0.48) and 5/2 (≈ 0.79) islands overlap, so they flatten one interval. + widths = Dict("2/1" => 0.4, "5/2" => 0.3) + two = island_flattening(equil; widths) + (a, b) = only(two.intervals) + @test (a, b) == (two.islands[1].psi_in, two.islands[2].psi_out) + alone = [island_flattening(equil; widths=Dict(k => v)).delta_W for (k, v) in widths] + @test two.delta_W < sum(alone) + n = 200_000 + grid = range(prof.xs[1], prof.xs[end]; length=n + 1) + mids = (grid[1:end-1] .+ grid[2:end]) ./ 2 + @test two.delta_W ≈ 1.5 * sum((p(x) - IF.flattened_pressure(p, [(a, b)], x)) * dV(x) for x in mids) * step(grid) rtol = 1e-6 + @test two.delta_p0 ≈ only(two.delta_p) ≈ p(a) - p(b) + # Inside the island the pressure sits at its outer-edge value. + @test IF.flattened_pressure(p, two.intervals, (a + b) / 2) ≈ p(b) + end + + @testset "edges, bootstrap fraction and bad input" begin + wide = island_flattening(equil; widths=Dict("2/1" => 4.0)) + @test only(wide.islands).clipped + # The whole profile flattens to its boundary value. + @test wide.delta_W ≈ wide.W - 1.5 * p(prof.xs[end]) * volume_inside(prof.xs[end]) rtol = 1e-8 + # A uniform ⟨j_bs·B⟩ loses the volume fraction of the flattened interval. + frac = island_flattening(equil; widths=Dict("2/1" => 0.05), jbs_B=_ -> 1.0) + (a, b) = only(frac.intervals) + @test frac.jbs_fraction_lost ≈ (volume_inside(b) - volume_inside(a)) / volume_inside(prof.xs[end]) rtol = 1e-3 + @test all(iszero, frac.jbs_B_mod[a. 0.05)) + @test_throws ArgumentError island_flattening(equil; widths=Dict("2/1" => -0.1)) + end + + @testset "HDF5 output keeps the rest of Tearing/" begin + mktempdir() do dir + path = joinpath(dir, "gpec.h5") + h5open(path, "w") do f + f["Tearing/enabled"] = 1 + f["Input/gpec_toml_raw"] = sprint(TOML.print, Dict("Tearing" => Dict("inner_model" => "SLAYER"))) + end + save_to_hdf5(path, one) + h5open(path, "r") do f + @test read(f["Tearing/enabled"]) == 1 + @test read(f["Tearing/IslandFlattening/delta_W"]) == one.delta_W + @test read(f["Tearing/IslandFlattening/Profiles/p_mod"]) == one.p_mod + recorded = TOML.parse(read(f["Input/gpec_toml_raw"]))["Tearing"] + @test recorded["inner_model"] == "SLAYER" + @test recorded["IslandFlattening"]["widths"] == Dict("2/1" => 0.05) + end + end + end +end diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 4f657bbd23..2851760f7e 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -645,5 +645,26 @@ include("display_check.jl") @test minimum(raised.kinetic.total.method_results["fgar"].psi_grid) >= raised.pe.psi_grid[1] end end + + @testset "island flattening" begin + # [Tearing.IslandFlattening] alone runs no layer solve; it is island_flattening on the solve's equilibrium. + widths = Dict{String,Any}("2/1" => 0.05, "3/1" => 0.02) + extra = Dict{String,Any}("Tearing" => Dict{String,Any}("IslandFlattening" => Dict{String,Any}("widths" => widths))) + mktempdir() do dir + reference = GPEC.main([_stage_deck(dir; extra_sections=extra)]) + @test reference.tearing === nothing + api = island_flattening(reference.ffs.equil; widths) + @test api.delta_W == reference.island_flattening.delta_W + @test api.p_mod == reference.island_flattening.p_mod + h5open(joinpath(dir, "gpec.h5")) do f + @test read(f["Tearing/IslandFlattening/delta_W"]) == api.delta_W + replay = GPEC.parse_deck(TOML.parse(read(f["Input/gpec_toml_raw"])), dir) + @test replay.island_flattening == Dict{String,Float64}(widths) + @test replay.tearing === nothing + end + @test_throws ArgumentError GPEC.island_flattening_from_deck(Dict("Tearing" => Dict("IslandFlattening" => Dict("widths" => 0.05)))) + @test_throws ArgumentError GPEC.island_flattening_from_deck(Dict("Tearing" => Dict("IslandFlattening" => Dict("width" => widths)))) + end + end end end diff --git a/test/suites.jl b/test/suites.jl index dfe79f9e4f..781be84d26 100644 --- a/test/suites.jl +++ b/test/suites.jl @@ -30,6 +30,7 @@ const TEST_SUITES = [ "runtests_tj_analytic.jl" => "core", "runtests_kinetic_profiles.jl" => "core", "runtests_resist_eval.jl" => "core", + "runtests_island_flattening.jl" => "api", "runtests_localstability.jl" => "core", "runtests_analysis_plots.jl" => "api", "runtests_slayer_params.jl" => "errorfields", From 5d0483f38402b4d64b2aa600718ed8113acb6959 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Fri, 13 Nov 2026 23:53:31 -0500 Subject: [PATCH 2/3] Tearing.IslandFlattening - MINOR - Name the outputs by the HDF5 conventions and keep an empty [Tearing] an error Delta_W, Delta_W_over_W, Delta_p0, Intervals/Delta_p and Islands/rational_{m,n,psi}; the input widths are no longer echoed into Islands/. A [Tearing] with nothing in it stops on inner_model again; only one holding [Tearing.IslandFlattening] skips the layer solve. The bootstrap fraction is stated as a fraction of the given profile's volume integral, and the docs say how to get the toroidal current fraction. Tests for a peaked bootstrap profile and inner-edge clipping. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae --- docs/src/inner_layer.md | 9 +++-- .../cases/diiid_n1_island_flattening.toml | 10 +++--- src/Deck.jl | 2 +- src/Tearing/HDF5Output.jl | 36 +++++++++---------- src/Tearing/IslandFlattening/Flattening.jl | 15 ++++---- test/runtests_island_flattening.jl | 13 +++++-- test/runtests_result_struct.jl | 4 +-- test/runtests_solve_api.jl | 2 +- 8 files changed, 52 insertions(+), 39 deletions(-) diff --git a/docs/src/inner_layer.md b/docs/src/inner_layer.md index 09b9aa4f9f..f41ed8fec4 100644 --- a/docs/src/inner_layer.md +++ b/docs/src/inner_layer.md @@ -62,7 +62,9 @@ p_\mathrm{mod}(\psi) = p(\psi) + \sum_{k:\, b_k > \psi} \left[p(b_k) - p(\max(\p ``` 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 ``\langle j_\mathrm{bs}\cdot B\rangle(\psi)`` (the `jbs_B` keyword of `island_flattening`), the result also carries the bootstrap current removed: to leading order ``j_\mathrm{bs} \propto p'``, so it vanishes on the flattened intervals and is unchanged elsewhere. +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 @@ -103,6 +105,7 @@ 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. @@ -128,8 +131,8 @@ 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_/` | 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 (m, n, ψ_N, w, its clipped interval) and the merged flattened intervals with the pressure drop across each | +| `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 diff --git a/regression-harness/cases/diiid_n1_island_flattening.toml b/regression-harness/cases/diiid_n1_island_flattening.toml index e87d10be10..31ac6fdc61 100644 --- a/regression-harness/cases/diiid_n1_island_flattening.toml +++ b/regression-harness/cases/diiid_n1_island_flattening.toml @@ -20,7 +20,7 @@ noise_threshold = 1e-6 order = 10 [quantities.delta_W] -h5path = "Tearing/IslandFlattening/delta_W" +h5path = "Tearing/IslandFlattening/Delta_W" type = "real_scalar" extract = "value" label = "stored energy removed" @@ -28,7 +28,7 @@ noise_threshold = 1e-6 order = 11 [quantities.delta_W_over_W] -h5path = "Tearing/IslandFlattening/delta_W_over_W" +h5path = "Tearing/IslandFlattening/Delta_W_over_W" type = "real_scalar" extract = "value" label = "fraction of W removed" @@ -36,7 +36,7 @@ noise_threshold = 1e-12 order = 12 [quantities.delta_p0] -h5path = "Tearing/IslandFlattening/delta_p0" +h5path = "Tearing/IslandFlattening/Delta_p0" type = "real_scalar" extract = "value" label = "on-axis pressure drop" @@ -44,7 +44,7 @@ noise_threshold = 1e-6 order = 13 [quantities.island_psi] -h5path = "Tearing/IslandFlattening/Islands/psi" +h5path = "Tearing/IslandFlattening/Islands/rational_psi" type = "real_vector" extract = "all_real" label = "island rational surfaces" @@ -52,7 +52,7 @@ noise_threshold = 1e-12 order = 20 [quantities.interval_delta_p] -h5path = "Tearing/IslandFlattening/Intervals/delta_p" +h5path = "Tearing/IslandFlattening/Intervals/Delta_p" type = "real_vector" extract = "all_real" label = "pressure drop per merged interval" diff --git a/src/Deck.jl b/src/Deck.jl index ca5a0b48e3..2f99ec7951 100644 --- a/src/Deck.jl +++ b/src/Deck.jl @@ -362,7 +362,7 @@ function tearing_from_deck(inputs::AbstractDict) haskey(inputs, "Tearing") || return nothing # A [Tearing] holding only [Tearing.IslandFlattening] runs the flattening without a layer solve. table = _without(inputs["Tearing"], "IslandFlattening") - isempty(table) && return nothing + 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)) diff --git a/src/Tearing/HDF5Output.jl b/src/Tearing/HDF5Output.jl index a69c273c07..66c6a8e4e3 100644 --- a/src/Tearing/HDF5Output.jl +++ b/src/Tearing/HDF5Output.jl @@ -370,26 +370,25 @@ end # ---------- island flattening ---------- const ISLAND_FLATTENING_H5_ANNOTATIONS = [ - "Islands/m" => (; long_name="poloidal mode number m of each island", dims=("island",)), - "Islands/n" => (; long_name="toroidal mode number n of each island", dims=("island",)), - "Islands/psi" => (; long_name="rational surface ψ_N of each island", dims=("island",)), - "Islands/w" => (; long_name="full island width in ψ_N", dims=("island",)), + "Islands/rational_m" => (; long_name="poloidal mode number m of each island", dims=("island",)), + "Islands/rational_n" => (; long_name="toroidal mode number n of each island", dims=("island",)), + "Islands/rational_psi" => (; long_name="rational surface ψ_N of each island", dims=("island",)), "Islands/psi_in" => (; long_name="inner end of the flattened interval (ψ_N, clipped to the profile grid)", dims=("island",)), "Islands/psi_out" => (; long_name="outer end of the flattened interval (ψ_N, clipped to the profile grid)", dims=("island",)), "Islands/clipped" => (; long_name="flag: the island reaches past an end of the profile grid", dims=("island",)), "Intervals/psi_in" => (; long_name="inner end of each merged flattened interval (ψ_N)", dims=("interval",)), "Intervals/psi_out" => (; long_name="outer end of each merged flattened interval (ψ_N)", dims=("interval",)), - "Intervals/delta_p" => (; long_name="pressure drop p(ψ_in) − p(ψ_out) across each merged interval", units="Pa", dims=("interval",)), + "Intervals/Delta_p" => (; long_name="pressure drop p(ψ_in) − p(ψ_out) across each merged interval", units="Pa", dims=("interval",)), "Profiles/psi" => (; long_name="normalized poloidal flux ψ_N (profile grid with the interval ends added)", dims=("psi",)), "Profiles/p" => (; long_name="equilibrium pressure", units="Pa", dims=("psi",)), "Profiles/p_mod" => (; long_name="pressure flattened across the islands, anchored at the boundary", units="Pa", dims=("psi",)), - "Profiles/jbs_B" => (; long_name="given bootstrap ⟨j_bs·B⟩", units="A*T/m^2", dims=("psi",)), - "Profiles/jbs_B_mod" => (; long_name="bootstrap ⟨j_bs·B⟩ with the flattened intervals removed (j_bs ∝ p′ to leading order)", units="A*T/m^2", dims=("psi",)), + "Profiles/jbs_B" => (; long_name="given bootstrap profile (⟨j_bs·B⟩, or a current density for a current fraction)", dims=("psi",)), + "Profiles/jbs_B_mod" => (; long_name="given bootstrap profile with the flattened intervals removed (j_bs ∝ p′ to leading order)", dims=("psi",)), "W" => (; long_name="stored energy (3/2)∫p dV", units="J"), - "delta_W" => (; long_name="stored energy removed by the flattening (3/2)∫(p − p_mod) dV", units="J"), - "delta_W_over_W" => (; long_name="fraction of the stored energy removed"), - "delta_p0" => (; long_name="on-axis pressure drop", units="Pa"), - "jbs_fraction_lost" => (; long_name="volume-integrated fraction of ⟨j_bs·B⟩ removed (NaN without a bootstrap profile)") + "Delta_W" => (; long_name="stored energy removed by the flattening (3/2)∫(p − p_mod) dV", units="J"), + "Delta_W_over_W" => (; long_name="fraction of the stored energy removed"), + "Delta_p0" => (; long_name="on-axis pressure drop", units="Pa"), + "jbs_fraction_lost" => (; long_name="fraction of the volume integral of the given bootstrap profile removed (NaN without one)") ] """ @@ -401,16 +400,15 @@ keeping the rest of `Tearing/`. function save_to_hdf5(parent::Union{HDF5.File,HDF5.Group}, r::IslandFlatteningResult) g = replace_group!(ensure_group(parent, "Tearing"), "IslandFlattening") isl = r.islands - g["Islands/m"] = [i.m for i in isl] - g["Islands/n"] = [i.n for i in isl] - g["Islands/psi"] = [i.psi for i in isl] - g["Islands/w"] = [i.w for i in isl] + g["Islands/rational_m"] = [i.m for i in isl] + g["Islands/rational_n"] = [i.n for i in isl] + g["Islands/rational_psi"] = [i.psi for i in isl] g["Islands/psi_in"] = [i.psi_in for i in isl] g["Islands/psi_out"] = [i.psi_out for i in isl] g["Islands/clipped"] = Int[i.clipped for i in isl] g["Intervals/psi_in"] = first.(r.intervals) g["Intervals/psi_out"] = last.(r.intervals) - g["Intervals/delta_p"] = r.delta_p + g["Intervals/Delta_p"] = r.delta_p g["Profiles/psi"] = r.psi g["Profiles/p"] = r.p g["Profiles/p_mod"] = r.p_mod @@ -419,9 +417,9 @@ function save_to_hdf5(parent::Union{HDF5.File,HDF5.Group}, r::IslandFlatteningRe g["Profiles/jbs_B_mod"] = r.jbs_B_mod end g["W"] = r.W - g["delta_W"] = r.delta_W - g["delta_W_over_W"] = r.delta_W / r.W - g["delta_p0"] = r.delta_p0 + g["Delta_W"] = r.delta_W + g["Delta_W_over_W"] = r.delta_W / r.W + g["Delta_p0"] = r.delta_p0 g["jbs_fraction_lost"] = r.jbs_fraction_lost Utilities.HDF5Annotations.annotate!(g, ISLAND_FLATTENING_H5_ANNOTATIONS) Utilities.record_inputs!(parent, r) diff --git a/src/Tearing/IslandFlattening/Flattening.jl b/src/Tearing/IslandFlattening/Flattening.jl index 38ae221eb4..9e1ca90ec5 100644 --- a/src/Tearing/IslandFlattening/Flattening.jl +++ b/src/Tearing/IslandFlattening/Flattening.jl @@ -42,10 +42,11 @@ Static flattening of the equilibrium pressure across `islands`. flattened pressures on it (Pa) - `W`, `delta_W`: stored energy (3/2)∫p dV and the part the islands remove (J) - `delta_p0`: on-axis pressure drop, the sum of `delta_p` (Pa) - - `jbs_B`, `jbs_B_mod`, `jbs_fraction_lost`: the given bootstrap ⟨j_bs·B⟩ on `psi`, the - same with the flattened intervals removed, and the volume-integrated fraction removed; - empty and NaN when no profile was given. To leading order j_bs ∝ p′, so it vanishes - where p′ does and is unchanged elsewhere. + - `jbs_B`, `jbs_B_mod`, `jbs_fraction_lost`: the given bootstrap profile on `psi`, the same + with the flattened intervals removed, and the fraction of its volume integral on them; + empty and NaN when no profile was given. To leading order j_bs ∝ p′, so it vanishes where + p′ does and is unchanged elsewhere. For ⟨j_bs·B⟩ the fraction is of ∫⟨j_bs·B⟩dV; for the + toroidal bootstrap current pass j_φ⟨1/R⟩, since dA = ⟨1/R⟩dV/2π. """ struct IslandFlatteningResult widths::Dict{String,Float64} @@ -100,8 +101,10 @@ end island_flattening(equil; widths, jbs_B=nothing) -> IslandFlatteningResult Flatten the pressure of `equil` across an island of full width `widths["m/n"]` (ψ_N) at -every rational surface q = m/n. `jbs_B`, a function of ψ_N giving ⟨j_bs·B⟩, adds the -bootstrap current the flattening removes. +every rational surface q = m/n (m, n > 0). `jbs_B`, a function of ψ_N giving a bootstrap +profile, adds the part the flattening removes. 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. """ function island_flattening(equil::PlasmaEquilibrium; widths::AbstractDict, jbs_B=nothing) prof = equil.profiles diff --git a/test/runtests_island_flattening.jl b/test/runtests_island_flattening.jl index 42c2cda8f4..ae13052986 100644 --- a/test/runtests_island_flattening.jl +++ b/test/runtests_island_flattening.jl @@ -25,7 +25,7 @@ using HDF5 for i in (length(grid)-1):-1:1 pm[i] = pm[i+1] - 0.5 * (dp[i+1] + dp[i]) * step(grid) end - @test pm[1] ≈ IF.flattened_pressure(p, one.intervals, grid[1]) rtol = 1e-6 + @test pm[1] ≈ IF.flattened_pressure(p, one.intervals, grid[1]) rtol = 1e-5 dW = 1.5 * sum((p(x) - y) * dV(x) for (x, y) in zip(grid, pm)) * step(grid) @test dW ≈ one.delta_W rtol = 1e-3 end @@ -65,6 +65,15 @@ using HDF5 (a, b) = only(frac.intervals) @test frac.jbs_fraction_lost ≈ (volume_inside(b) - volume_inside(a)) / volume_inside(prof.xs[end]) rtol = 1e-3 @test all(iszero, frac.jbs_B_mod[a. 0.05), jbs_B=g) + mid(f, lo, hi; n=20_000) = sum(f(lo + (hi - lo) * (k - 0.5) / n) for k in 1:n) * (hi - lo) / n + @test peaked.jbs_fraction_lost ≈ mid(x -> g(x) * dV(x), a, b) / mid(x -> g(x) * dV(x), prof.xs[1], prof.xs[end]) rtol = 1e-6 + @test peaked.jbs_fraction_lost < frac.jbs_fraction_lost + # An island reaching inside the grid's first point is clipped there. + core = only(island_flattening(equil; widths=Dict("2/1" => 2 * only(one.islands).psi + 0.1)).islands) + @test core.clipped && core.psi_in == prof.xs[1] @test isnan(one.jbs_fraction_lost) && isempty(one.jbs_B) @test_throws ArgumentError island_flattening(equil; widths=Dict("2-1" => 0.05)) @test_throws ArgumentError island_flattening(equil; widths=Dict("2/1" => -0.1)) @@ -80,7 +89,7 @@ using HDF5 save_to_hdf5(path, one) h5open(path, "r") do f @test read(f["Tearing/enabled"]) == 1 - @test read(f["Tearing/IslandFlattening/delta_W"]) == one.delta_W + @test read(f["Tearing/IslandFlattening/Delta_W"]) == one.delta_W @test read(f["Tearing/IslandFlattening/Profiles/p_mod"]) == one.p_mod recorded = TOML.parse(read(f["Input/gpec_toml_raw"]))["Tearing"] @test recorded["inner_model"] == "SLAYER" diff --git a/test/runtests_result_struct.jl b/test/runtests_result_struct.jl index f38ada7aa3..1baa914476 100644 --- a/test/runtests_result_struct.jl +++ b/test/runtests_result_struct.jl @@ -45,8 +45,8 @@ include("display_check.jl") ret = GeneralizedPerturbedEquilibrium.main([dir]) # One field per stage, nothing for the sections the deck does not hold. - @test keys(ret) == (:equil, :local_stability, :ffs, :match, :resp, :pe, :kinetic, :error_fields, :tearing) - @test ret.equil === ret.ffs.equil && all(isnothing, (ret.match, ret.resp, ret.pe, ret.kinetic, ret.error_fields, ret.tearing)) + @test keys(ret) == (:equil, :local_stability, :ffs, :match, :resp, :pe, :kinetic, :error_fields, :tearing, :island_flattening) + @test ret.equil === ret.ffs.equil && all(isnothing, (ret.match, ret.resp, ret.pe, ret.kinetic, ret.error_fields, ret.tearing, ret.island_flattening)) ffs = ret.ffs forward_result[] = ffs @test ffs isa FFS.ForceFreeStatesResult diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 2851760f7e..ba9e419471 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -657,7 +657,7 @@ include("display_check.jl") @test api.delta_W == reference.island_flattening.delta_W @test api.p_mod == reference.island_flattening.p_mod h5open(joinpath(dir, "gpec.h5")) do f - @test read(f["Tearing/IslandFlattening/delta_W"]) == api.delta_W + @test read(f["Tearing/IslandFlattening/Delta_W"]) == api.delta_W replay = GPEC.parse_deck(TOML.parse(read(f["Input/gpec_toml_raw"])), dir) @test replay.island_flattening == Dict{String,Float64}(widths) @test replay.tearing === nothing From f6aa34b8a9578ccf3aaf33c865088771b0a6f721 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 14 Nov 2026 00:24:16 -0500 Subject: [PATCH 3/3] Utilities - MINOR - Keep recorded sub-tables a result does not state, by the rule the MRE branch uses record_inputs! now matches the modified Rutherford branch line for line: a result owns its section's scalars and the sub-tables it states, and recorded sub-tables it does not state are kept. This replaces the sub-table heuristic of the first commit, so the two branches merge in either order. Two comments reworded. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae --- src/GeneralizedPerturbedEquilibrium.jl | 2 +- src/Tearing/HDF5Output.jl | 2 +- src/Utilities/Output.jl | 20 +++++++++----------- 3 files changed, 11 insertions(+), 13 deletions(-) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 1a8dc16d53..f451c77bd4 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -224,7 +224,7 @@ function main_from_inputs( save_to_hdf5(output_file, tearing) end - # After the tearing save, which replaces Tearing/. + # 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) diff --git a/src/Tearing/HDF5Output.jl b/src/Tearing/HDF5Output.jl index 66c6a8e4e3..6418cc7702 100644 --- a/src/Tearing/HDF5Output.jl +++ b/src/Tearing/HDF5Output.jl @@ -13,7 +13,7 @@ import ..Utilities.HDF5Annotations: replace_group!, ensure_group, write_ragged! # │ uncoupled entries are in each surface's E×B frame # ├── Scan/ -- optional: full Q/Δ scan data # ├── CriticalResonantField/ -- optional: torque-balance critical b_r/B_φ per surface -# └── IslandFlattening/ -- its own result; a TearingResult save replaces Tearing/, so it is written after +# └── IslandFlattening/ -- its own result, written after the TearingResult in a deck run using HDF5 diff --git a/src/Utilities/Output.jl b/src/Utilities/Output.jl index ce30956f19..f255e717d1 100644 --- a/src/Utilities/Output.jl +++ b/src/Utilities/Output.jl @@ -28,23 +28,21 @@ function deck_sections end """ record_inputs!(target, result) -Merge the deck sections of `result` into the run record `Input/gpec_toml_raw` of `target`, replacing any section of -the same name, so every saved result adds the deck it was computed from. A result that states only sub-tables of a -section (such as `[Tearing.IslandFlattening]`) replaces those sub-tables and keeps the rest of the section. +Merge the deck sections of `result` into the run record `Input/gpec_toml_raw` of `target`, so every saved result adds +the deck it was computed from. A result replaces its section's scalars and the sub-tables it states, and keeps recorded +sub-tables it does not state (a section stated only as sub-tables keeps the recorded scalars too). """ function record_inputs!(target::Union{HDF5.File,HDF5.Group}, result) path = "Input/gpec_toml_raw" deck = haskey(target, path) ? TOML.parse(read(target[path])) : Dict{String,Any}() - for (name, section) in deck_sections(result) - if _only_subtables(section) && get(deck, name, nothing) isa AbstractDict - merge!(deck[name], section) - else - deck[name] = section - end + # A result owns its section's scalars and the sub-tables it states; recorded sub-tables it does not state stay. + for (k, v) in deck_sections(result) + old = get(deck, k, nothing) + deck[k] = + !(v isa AbstractDict && old isa AbstractDict) ? v : + all(x -> x isa AbstractDict, values(v)) ? merge(old, v) : merge(filter(p -> p.second isa AbstractDict, old), v) end haskey(target, path) && HDF5.delete_object(target, path) target[path] = sprint(TOML.print, deck) return nothing end - -_only_subtables(section) = section isa AbstractDict && !isempty(section) && all(v -> v isa AbstractDict, values(section))