diff --git a/.claude/agents/regression-guardian.md b/.claude/agents/regression-guardian.md index e98ecac11f..70252c9c90 100644 --- a/.claude/agents/regression-guardian.md +++ b/.claude/agents/regression-guardian.md @@ -35,7 +35,7 @@ Useful invocations: - `--list-cases` — enumerate available cases. **Always run this first** to get the current set; cases are added over time, so never assume the list. (At time of writing: `diiid_n1`, `diiid_n1_riccati`, `diiid_slayer_n1`, `efit_fixedbdy_separatrix`, `gal_resistive_diiid`, - `gal_resistive_pe`, `ggj_ray_q500i`, `ggj_reference`, `solovev_n1`, `solovev_multi_n`, + `ggj_ray_q500i`, `ggj_reference`, `solovev_n1`, `solovev_multi_n`, `solovev_kinetic_nuzero`, `tj_n1`, `tj_kinetic_calculated`, `tj_kinetic_multiion`, `tj_kinetic_ntv`.) - `--cases --refs ` — compare two refs (branches/commits), e.g. diff --git a/README.md b/README.md index b34dd7cc58..46f7c6fc6c 100644 --- a/README.md +++ b/README.md @@ -1,12 +1,13 @@ # GPEC — Generalized Perturbed Equilibrium Code -GPEC is a Julia implementation of the [GPEC suite](https://github.com/PrincetonUniversity/GPEC) for magnetohydrodynamic (MHD) analysis of fusion plasmas. The Fortran suite (DCON, VACUUM, STRIDE, GPEC, PENTRC and SLAYER) has long been used to model the response of tokamak plasmas to small non-axisymmetric magnetic perturbations, from error-field correction and resonant magnetic perturbation optimization to kink and tearing stability. This implementation consolidates them into one code base of cooperating physics modules that share data in memory behind a unified, modern interface, which removes the hand-offs between separate codes and makes parameter scans, uncertainty quantification and integrated modeling practical. +This Julia implementation of GPEC continues the legacy of the Fortran [GPEC suite](https://github.com/PrincetonUniversity/GPEC) for magnetohydrodynamic (MHD) analysis of fusion plasmas. The Fortran suite (DCON, VACUUM, STRIDE, GPEC, PENTRC and SLAYER) has long been used to model the response of tokamak plasmas to small non-axisymmetric magnetic perturbations, from error-field correction and resonant magnetic perturbation optimization to kink and tearing stability. This implementation consolidates them into one code base of cooperating physics modules that share data in memory behind a unified, modern interface, which removes the hand-offs between separate codes and makes parameter scans, uncertainty quantification and integrated modeling practical. -It performs a three-stage analysis pipeline: +One input deck runs the stages it names, in order, into one output file: -1. **Equilibrium reconstruction** from experimental data (EFIT, CHEASE) or analytical models — solves the Grad-Shafranov equation, computes flux surfaces, safety factor q(ψ), and builds field quantity splines. -2. **Ideal MHD stability analysis** via Newcomb's criterion (DCON-style ODE integration) — identifies singular surfaces where q = m/n, computes force-free eigenfunctions ξ(ψ,θ), evaluates the tearing stability parameter Δ' at each rational surface, and checks fixed/free boundary stability. -3. **Perturbed equilibrium** — quantifies the self-consistent plasma response to external magnetic perturbations (RMPs, error fields, correction coils). Key outputs include resonant flux, island half-widths, and the Chirikov overlap parameter — the principal diagnostics for ELM suppression, error field correction, and disruption avoidance. +1. **Equilibrium** — reads an EFIT or CHEASE reconstruction, an IMAS data dictionary or an analytic model, maps it to straight-field-line flux coordinates, and finds the rational surfaces q = m/n. +2. **Stability** — ideal MHD stability by Newcomb's criterion (the DCON integration, or the STRIDE Riccati and RDCON Galerkin formulations): the force-free states ξ(ψ,θ), the free-boundary kink energies with the VACUUM response, and the tearing stability matrix Δ′. +3. **Perturbed equilibrium** — the plasma response to external non-axisymmetric fields (RMPs, error fields, correction coils): resonant fields, island half-widths and the Chirikov overlap parameter. +4. **Torque, error fields and tearing** — the neoclassical toroidal viscosity torque (PENTRC), coil-tolerance and locking-risk analysis, and tearing growth rates and critical resonant fields from inner-layer models (SLAYER, GGJ). Full documentation: [https://openfusiontoolkit.github.io/GPEC/dev/](https://openfusiontoolkit.github.io/GPEC/dev/) diff --git a/docs/DOC_STANDARD.md b/docs/DOC_STANDARD.md index a1af2043e5..bafe45902f 100644 --- a/docs/DOC_STANDARD.md +++ b/docs/DOC_STANDARD.md @@ -13,7 +13,7 @@ render it to the website. It is a contributor-facing reference, read here in the repository — which is why `developer_notes.md` links to it by GitHub URL rather than as a site-internal page. -The `Ballooning and Mercier Local Stability` and `Inner Layer Module` pages are +The `Stability` and `Local Stability` pages are the reference exemplars; new or reworked module pages should follow their shape. ## Module-page template @@ -21,24 +21,18 @@ the reference exemplars; new or reworked module pages should follow their shape. A module reference page (`docs/src/.md`) has these sections, in order. Not every section is large for every module, but the skeleton is the same. -1. **Overview** — what the module computes and where it sits in the - equilibrium → stability → perturbed-equilibrium pipeline. -2. **Governing equations** — the physics the module discretizes, written as - numbered display equations. Each equation or block cites the paper it comes - from (the PDFs in `docs/resources/`), e.g. *(GWP2016 Eq. 11)*. Citing - equation numbers is required, not optional — it is what makes the - implementation auditable against theory. -3. **Numerical method** — the discretization and algorithm: what is expanded in - what basis, how the linear/eigenvalue/BVP problem is assembled and solved, - and the *design choices* that matter (why this contour, this ordering, this - precision). This is the part that distinguishes a methods page from an API - dump. -4. **Validation & benchmarks** — evidence the method is correct: a table of - pinned values cross-checked against an independent code or analytic result, - and figures (see below) showing convergence, invariance, or agreement. -5. **Practical usage** — the configuration keys and the public entry points a - user calls, with a short worked example. -6. **API Reference** — an `@autodocs` block for the module so every exported +1. **What it computes** — the physics of the stage in at most two paragraphs, + and where it sits in the pipeline (a link back to the Workflow page). +2. **Governing equations** — the equations the module discretizes, as display + equations, each citing its paper by "Author (year)" with a link to the + citations page; equation numbers where the paper is in `docs/resources/`. +3. **Running it** — the deck section(s), as a snippet copied from a shipped + example that is named, and the matching scripting-API call. +4. **Outputs** — a table of the `gpec.h5` datasets the stage writes, by full + path, with meaning and units. +5. **Numerical notes** — the discretization and algorithm choices that matter + to a user, kept short, with the validation figures (see below). +6. **API reference** — an `@autodocs` block for the module so every exported docstring is rendered and `checkdocs=:exports` stays satisfied. ## Figures diff --git a/docs/development/regression-harness.md b/docs/development/regression-harness.md index f89d7d5ee8..32613aa13c 100644 --- a/docs/development/regression-harness.md +++ b/docs/development/regression-harness.md @@ -205,3 +205,42 @@ If the two refs in a comparison ran under different thread counts, the report fl - `0` — every run completed (and, with `--fail-on-change`, nothing changed) - `1` — a run failed, or a quantity changed under `--fail-on-change` + +## Pinned Δ′ and grid convergence + +The ideal Δ′ diagonal pinned in `test/runtests_parallel_integration.jl` and tracked by the harness is a +single-grid snapshot that catches unintended changes, not a converged value. A `psi_accuracy` scan of the +DIII-D-like SLAYER deck on the automatic grid (`implied` is the knot count `implied_knot_count` asks for +after the refined solve): + +| psi_accuracy | knots | implied | dpm[1,1] | dpm[2,2] | dpm[3,3] | gamma 2/1 (Hz) | +|---|---|---|---|---|---|---| +| 2.0e-3 | 172 | – | 6.850 | -5.783 | -16.308 | 165.3 | +| 1.0e-3 | 205 | – | 6.387 | -5.873 | -16.867 | 154.3 | +| 5.0e-4 | 252 | – | 8.003 | -6.105 | -16.413 | 192.9 | +| 2.5e-4 | 307 | 522 | 8.567 | -6.235 | -16.529 | 206.2 | +| 1.25e-4 | 372 | 666 | 8.146 | -6.285 | -16.649 | 196.2 | +| 6.25e-5 | 457 | 916 | 8.911 | -6.477 | -16.510 | 214.4 | +| 3.125e-5 | 559 | 1226 | 9.514 | -6.260 | -16.485 | 228.7 | + +The same deck on a fixed `core_edge_packed` grid: + +| mpsi | dpm[1,1] | dpm[2,2] | dpm[3,3] | gamma 2/1 (Hz) | +|---|---|---|---|---| +| 128 | 5.208 | -8.666 | -19.157 | 126.2 | +| 256 | 7.097 | -5.187 | -15.729 | 171.2 | +| 512 | 8.489 | -6.217 | -16.517 | 204.4 | +| 1024 | 8.932 | -6.298 | -16.529 | 214.9 | +| 2048 | 8.914 | -6.385 | -16.468 | 214.5 | + +- Convergence is per surface: the q = 4 `dpm[3,3]` settles near -16.5 on both grid families, while the q = 2 + `dpm[1,1]` converges to about 8.91 only on the fixed grid and keeps moving on the automatic one. +- The SLAYER growth rate is linear in Δ′ here (`gamma/dpm[1,1]` ≈ 24.1 within 0.5%), so a converged Δ′ gives a + converged growth rate. +- The automatic grid stops after its second pass, and the knot count it implies grows faster than the grid in + use as `psi_accuracy` tightens, so tightening moves it further from self-consistency. + +A plateau criterion (scan `psi_accuracy` and the truncation controls, pin the stationary value per surface) +would be a real convergence statement; it needs either the fixed grid or knot refinement iterated to a fixed +point, and is not implemented. Until then, read a pinned diagonal as a sign and order-of-magnitude signal and +re-pin it when the equilibrium grid changes. diff --git a/docs/make.jl b/docs/make.jl index 8b5fe99dbc..1bd8818aa2 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -28,15 +28,15 @@ makedocs(; "Workflow" => "workflow.md", "Scripting API" => "api.md", "Conventions Reference" => "conventions.md", - "API Reference" => [ - "Vacuum" => "vacuum.md", + "Modules" => [ "Equilibrium" => "equilibrium.md", - "Stability Analysis" => "stability.md", + "Vacuum" => "vacuum.md", + "Stability" => "stability.md", "Galerkin Solver" => "galerkin.md", - "Ballooning Local Stability" => "ballooning.md", - "KineticForces" => "kinetic_forces.md", + "Local Stability" => "ballooning.md", "Forcing Terms" => "forcing_terms.md", "Perturbed Equilibrium" => "perturbed_equilibrium.md", + "Kinetic Forces" => "kinetic_forces.md", "Error Fields" => "error_fields.md", "Tearing" => "inner_layer.md", "Analysis" => "analysis.md", diff --git a/docs/src/analysis.md b/docs/src/analysis.md index 575adf10e1..6034f2caf2 100644 --- a/docs/src/analysis.md +++ b/docs/src/analysis.md @@ -1,14 +1,18 @@ -# Analysis Module +# Analysis -The Analysis module provides post-processing and visualization utilities for GPEC simulation outputs. +`Analysis` holds the plotting and post-processing functions for GPEC results, grouped by the stage they read; each works on the in-memory result or on `gpec.h5`. ## Submodules -- `ForceFreeStates`: Plotting functions for DCON-style ideal MHD stability results -- `Equilibrium`: Plotting functions for equilibrium objects -- `CoilForcing`: Plotting functions for coil geometry and normal field spectra -- `PerturbedEquilibrium`: Plotting functions for perturbed equilibrium results -- `PerturbedEquilibriumModes`: Data helpers to convert modal GPEC output to (ψ, θ) and (ψ, θ, φ) grids +- `ForceFreeStates`: the eigenmode energies, the `crit` scan and the displacement profiles of a stability solve +- `Equilibrium`: flux surfaces, profiles and the radial grid of an equilibrium +- `CoilForcing`: coil geometry and the normal-field spectra it drives +- `PerturbedEquilibrium`: the plasma response, the resonant fields, island widths and Chirikov parameter against ψ_N or q +- `PerturbedEquilibriumModes`: helpers that convert the modal output to (ψ, θ) and (ψ, θ, φ) grids for real-space plots +- `ErrorFields`: coil sensitivities, tolerance distributions, locking risk and correction limits (see [Error Fields](error_fields.md)) + +`examples/DIIID-like_ideal_example/analyze_example.jl`, `examples/Solovev_ideal_example/analyze_example.jl` and `examples/DIIID-like_error_field_example/analyze_example.jl` are worked examples of each; the s–α maps of local stability use `salpha_reference`, `compute_ballooning_stability!` and `scan_delta_prime_map`. +Figure conventions for spectrum plots are in `docs/development/plotting.md` of the repository, and `write_imas` on the [Utilities](utilities.md) page exports results to an IMAS data dictionary. ```@docs GeneralizedPerturbedEquilibrium.Analysis diff --git a/docs/src/api.md b/docs/src/api.md index 2c63384863..c520cf44b5 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -1,98 +1,70 @@ # Scripting API -GPEC's pipeline is normally driven from a `gpec.toml` deck through -`GeneralizedPerturbedEquilibrium.main`. The scripting API exposes the same stages as -ordinary Julia functions, so a run can be built, parameterized and looped over in a script -without writing a deck. +A deck runs through `GeneralizedPerturbedEquilibrium.main`; the scripting API exposes the same stages as Julia functions, so a run can be built, parameterized and looped over without a deck, and a deck can never take a path a script cannot. +Every stage is a `solve` of a problem type with an algorithm, or a function on the previous stage's result; nothing is written until `save_to_hdf5` is called. ```julia using GeneralizedPerturbedEquilibrium -eq = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada") -ffs = solve(eq, Riccati(); nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true) -rmp = RMPField("coils.dat") -resp = plasma_response(ffs) -pe = perturbed_equilibrium(resp, rmp) +eq = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada", kinetic_file="kinetic.h5") +ffs = solve(eq, Forward(); nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true) # [ForceFreeStates]: Forward, no wall, for the response below +ls = solve(LocalStabilityProblem(eq; alpha_boundary_scan=true), Ballooning()) # [LocalStability] +resp = plasma_response(ffs; compute_singular_coupling=true) # [PerturbedEquilibrium] +pe = perturbed_equilibrium(resp, RMPField("coils.dat")) # with the [ForcingTerms] drive +kf = kinetic_torque(pe, KineticForcesControl()) # [KineticForces] +ef = error_fields(resp, pe.coil_sets, ForcingTerms.CoilConfig(ft_ctrl); # [ErrorFields]; ft_ctrl: the [ForcingTerms] control + tolerances=ErrorFields.read_tolerance_toml("tolerances.toml")) +tear = solve(TearingProblem(ffs), SLAYERModel(); config=TearingConfig()) # [Tearing] + +save_to_hdf5("gpec.h5", ffs; local_stability=ls) +foreach(r -> save_to_hdf5("gpec.h5", r), (pe, kf, ef, tear)) ``` -`solve(eq, alg; kwargs...)` is sugar for the canonical problem form: a -[`ForceBalanceProblem`](@ref) names WHAT is being solved (the perturbed-plasma -force-balance equation posed on this equilibrium, with this mode range, wall and closure) -and the integrator names HOW. A `PlasmaEquilibrium` hosts many possible problems; the -problem type keeps `solve` unambiguous as other problem classes appear. +`solve(eq, alg; kwargs...)` is sugar for the problem form `solve(ForceBalanceProblem(eq; kwargs...), alg)`: the problem names what is solved (the force-balance equation on this equilibrium, with this mode range, wall and closure) and the algorithm names how. +`main` translates each deck section into the problem, algorithm or control it names, calls these functions for every section present, saves each result into `[Run] output_file`, and returns `(; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing)`, each `nothing` when its section is absent. -```julia -prob = ForceBalanceProblem(eq; nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true) -ffs = solve(prob, Riccati()) -save_to_hdf5("gpec.h5", ffs) +| Stage | Deck section | Call | Result | +|---|---|---|---| +| Equilibrium | `[Equilibrium]` (or `[Solovev]`, `[LargeAspectRatio]`, `[TJAnalytic]`) | `PlasmaEquilibrium(file; kwargs...)`, `PlasmaEquilibrium(config)`, `PlasmaEquilibrium(dd)` | `PlasmaEquilibrium` | +| ForceFreeStates | `[ForceFreeStates]` + `[Forward]`, `[Riccati]` or `[Galerkin]` | `solve(ForceBalanceProblem(eq; ...), Forward() \| Riccati() \| Galerkin())` | `ForceFreeStatesResult` | +| LocalStability | `[LocalStability]` + `[Ballooning]` | `solve(LocalStabilityProblem(eq; ...), Ballooning(; ...))` | `LocalStabilityResult` | +| Match | `[Match]` + `[GGJ]` | `solve(MatchProblem(ffs; ...), GGJModel())` | `ForceFreeStatesResult` with `closure == :matched` | +| Drive | `[ForcingTerms]` | `RMPField(ft_ctrl)` (a `ForcingTermsControl`), `RMPField(file)`, `RMPField(modes)` | `RMPField` | +| PerturbedEquilibrium | `[PerturbedEquilibrium]` | `plasma_response(ffs; ...)`, then `perturbed_equilibrium(resp, rmp; ...)` | `PlasmaResponse`, `PerturbedEquilibriumState` | +| KineticForces | `[KineticForces]` | `kinetic_torque(pe, KineticForcesControl(; ...))` | `KineticForcesResult` | +| ErrorFields | `[ErrorFields]` | `error_fields(resp, pe.coil_sets, CoilConfig(...); ...)` | `ErrorFieldsResult` | +| Tearing | `[Tearing]` + `[SLAYER]` or `[GGJ]` | `solve(TearingProblem(ffs; ...), SLAYERModel() \| GGJModel(); config=TearingConfig(...))` | `TearingResult` | + +The [Workflow](workflow.md) page describes what each stage computes; the module pages document the keywords, which are the deck keys of the matching section. + +## Saving and rerunning + +`save_to_hdf5(target, result)` saves any stage result; `target` is a path, which opens or creates the file, or an open HDF5 file or group, which stays owned by the caller. +Each save records in `Input/gpec_toml_raw` the deck sections that state the problem, algorithm and controls the result was computed with, so a script run's file replays like a deck run's: + +```bash +./gpec gpec.h5 # main(["gpec.h5"]) in Julia +./gpec gpec.h5 --coil-source coils ``` -`solve` returns its result without writing a file. Call `save_to_hdf5(path, ffs)` -or pass an open HDF5 file or group to save it explicitly. - -The same `save_to_hdf5(target, result)` function saves an equilibrium, local-stability -result, perturbed-equilibrium state, kinetic-force state, tearing result, or error-field -result. Each module owns the datasets for its result. Paths open or create an HDF5 -file; an open file or group remains owned by the caller. Stage computations return -results without writing anything. Each save also records, in `Input/gpec_toml_raw`, the deck -sections stating the problem, algorithm and controls the result was computed with, so a script -run's file replays like a deck run's: `main(["gpec.h5"])`. A forcing no deck can state (a weighted -field, a spectrum or a sum) is replayed from the control-surface modes the save records. - -A deck runs through these same calls. `main` translates each section of `gpec.toml` into the -problem, algorithm or control it names, calls the function above for every section present, -and saves each result into the one `[Run] output_file`, so a deck can never take a path a script -cannot. It returns `(; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing)`, -each `nothing` when its section is absent; `ffs` is the result the later stages read (matched -when the deck has `[Match]`) and `match` its `MatchResult`. - -The objects map one-to-one onto the pipeline stages: - - - [`PlasmaEquilibrium`](@ref) reads and processes the equilibrium — the `[Equilibrium]` - section. Analytic equilibria take their parameters from their config object, with the - `eq_type` implied, and IMAS equilibria from a data dictionary: - - ```julia - eq = PlasmaEquilibrium(Equilibrium.SolovevConfig(; elongation=1.5); mpsi=0) - eq = PlasmaEquilibrium(dd; imas_cocos=11) - ``` - - Every equilibrium remembers its source, so with `mpsi=0` `solve` re-forms it on the same - refined two-pass grid the deck uses; `Equilibrium.reform(eq, psi_nodes)` re-forms it on any grid. - - `solve` runs the force-free-states stage — the `[ForceFreeStates]` section — and returns - a `ForceFreeStatesResult`, exactly the object the TOML driver publishes. - - [`RMPField`](@ref GeneralizedPerturbedEquilibrium.ForcingTerms.RMPField) describes the external field — the `[ForcingTerms]` section — without - reading anything from disk until it is applied. - - [`plasma_response`](@ref) builds the plasma-response and singular-coupling operators of a - `ForceFreeStatesResult` once, and `perturbed_equilibrium(resp, rmp)` drives them with one field — - together the `[PerturbedEquilibrium]` section. - - [`kinetic_torque`](@ref) evaluates the kinetic (NTV) torque of a perturbed equilibrium with the - model settings of a `KineticForcesControl` — the `[KineticForces]` section — and returns a - `KineticForcesResult`: `kf = kinetic_torque(pe, KineticForcesControl())`. - - [`error_fields`](@ref) assesses the error field of a set of coils on a plasma response — the - `[ErrorFields]` section — and returns an `ErrorFieldsResult`: - `ef = error_fields(resp, pe.coil_sets, ForcingTerms.CoilConfig(ft_ctrl); tolerances)`. - -## Kinetic profiles - -The equilibrium carries the kinetic profiles. Name the file with the `kinetic_file` (and, -for HDF5, `kinetic_group`) keyword, the same keys as the `[Equilibrium]` deck section, and it -is read once into `equil.kinetic`; every kinetic stage derives its own view from it. -[`Equilibrium.attach_kinetic_profiles!`](@ref) attaches a file or in-memory `KineticProfileData` to an -equilibrium built without one, such as an analytic equilibrium. +A coil run is replayed either from the stored deck coil configuration (`--coil-source forcing-modes`, the default, which needs the original geometry files) or from the coil geometry stored under `Input/RawInputs/Coils/` (`--coil-source coils`, which needs no original files and allows the equilibrium to be overridden). +A forcing file is replayed from its raw contents under `Input/RawInputs/ForcingTerms/`, and a forcing no deck can state (a weighted field, a spectrum or a sum) from the control-surface modes the save records. + +## Equilibrium and kinetic profiles + +Analytic equilibria take their parameters from a config object, with the `eq_type` implied, and IMAS equilibria from a data dictionary: ```julia -eq = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada", kinetic_file="kinetic.h5") +eq = PlasmaEquilibrium(Equilibrium.SolovevConfig(; elongation=1.5); mpsi=0) +eq = PlasmaEquilibrium(dd; imas_cocos=11) ``` +Every equilibrium remembers its source, so with `mpsi=0` `solve` re-forms it on the refined two-pass grid the deck uses; `Equilibrium.reform(eq, psi_nodes)` re-forms it on any grid. +The `kinetic_file` (and, for HDF5, `kinetic_group`) keyword reads the kinetic profiles once into `equil.kinetic`, from which every kinetic stage derives its view; [`Equilibrium.attach_kinetic_profiles!`](@ref) attaches a file or an in-memory `KineticProfileData` to an equilibrium built without one. + ## Combining forcing sources -`RMPField`s form a vector space: `+`, `-` and scalar `*` build lazy linear combinations that -record their terms and compute nothing until the perturbed-equilibrium stage materializes -them. Because the plasma response is linear in the forcing, driving with a combination is -exactly equivalent to combining the individually materialized fields — each term is -evaluated on the control surface and the mode amplitudes are summed. A complex scalar -phase-rotates a source. +`RMPField`s form a vector space: `+`, `-` and scalar `*` build lazy linear combinations that compute nothing until `perturbed_equilibrium` materializes them, and because the response is linear in the forcing, driving with a combination equals combining the materialized fields; a complex scalar phase-rotates a source. ```julia nominal = RMPField("nominal_efc.dat") @@ -100,119 +72,57 @@ weld_field = RMPField("weld_fields.h5") pe = perturbed_equilibrium(resp, 2.0 * nominal + 0.5 * weld_field) ``` -This is the substrate for error-field workflows: build per-unit sources independently -(a coil set per ampere, a displacement field per millimeter), then assemble physical cases -by weighting and summing — without recomputing anything per combination until the final -`perturbed_equilibrium` call, and with one `plasma_response` serving every combination. -`RMPField(modes)` wraps a control-surface spectrum computed elsewhere, already in the unit-norm -`Φ_x` convention, as a leaf that combines like any other. - -The weights are linear-combination coefficients, not physical amplitudes: sources that are -not scalar multiples of each other get their own description and the algebra. A coil set -with a failed conductor is `nominal - failed_coil`, not `0.9 * nominal`. +This is the substrate for error-field workflows: build per-unit sources (a coil set per ampere, a displacement field per millimetre), then assemble physical cases by weighting and summing, with one `plasma_response` serving every combination. +`RMPField(modes)` wraps a control-surface spectrum computed elsewhere, in the unit-norm ``\Phi_x`` convention, as a leaf. +The weights are linear-combination coefficients, not physical amplitudes: a coil set with a failed conductor is `nominal - failed_coil`, not `0.9 * nominal`. ## Choosing an integrator -The second argument of `solve` picks the formalism, and its fields are that formalism's -numerical settings. Physical inputs belong on `ForceBalanceProblem`; ODE tolerances and -storage settings belong on `Forward` or `Riccati`, and `extended_precision_bvp` belongs only -on `Riccati`. The current deck reader translates its existing section layout into those objects: +The second argument of `solve` picks the formalism and its fields are that formalism's numerical settings; physical inputs belong on `ForceBalanceProblem`: ```julia -ffs = solve(eq, Forward(; eulerlagrange_tolerance=1e-8); nn=1, vac_flag=true) # dense ξ profiles for PerturbedEquilibrium -ffs = solve(eq, Riccati(; nchunks=40); nn=1, vac_flag=true) # chunked propagators, Δ′ matrix -ffs = solve(eq, Galerkin(; nx=512); nn=1) # RDCON outer-region Galerkin Δ′ +ffs = solve(eq, Forward(; eulerlagrange_tolerance=1e-8); nn=1, vac_flag=true) # dense ξ profiles: the only input plasma_response accepts; kinetic solves +ffs = solve(eq, Riccati(; nchunks=40); nn=1, vac_flag=true) # chunked propagators, Δ′ matrix (default) +ffs = solve(eq, Galerkin(; nx=512); nn=1) # finite-element outer-region Δ′ ``` -Which products each formalism can supply differs; a result carries `nothing` in the fields -its integrator does not produce and consumers warn and skip rather than erroring. See the -[Stability Analysis](stability.md) page for the result struct and its capability gates. - -`singfac_min`, `frobenius_order`, `frobenius_psi_max`, `truncate_at_dW_peak`, and -`dW_edge_scan_start` remain problem settings. `Galerkin` rejects non-default values of the -last three. Its per-surface order adjustment is controlled by `frobenius_order_ceiling`, -and its factorization is selected automatically. +A result carries `nothing` in the fields its integrator does not produce, and consumers warn and skip rather than erroring; `plasma_response` in particular needs a `Forward()` solve with `wall=nothing` (it errors on a wall or a Galerkin result and skips with a warning on a Riccati one), see [Stability](stability.md). +A kinetic solve takes `Forward()` with the problem's `kinetic` keyword: `kinetic=KineticForcesControl(; ...)` evaluates the drift-kinetic matrices on the equilibrium's kinetic profiles, `kinetic=TestKineticMatrices(; scale)` uses synthetic matrices, `kinetic_factor` scales either, and `km = kinetic_matrices(prob)` evaluates the source once for reuse in later solves ([`KineticMatrices`](@ref GeneralizedPerturbedEquilibrium.ForceFreeStates.KineticMatrices)). ## Inner-layer matching -Inner-layer matching is its own problem, posed on a finished solve: a `MatchProblem` holds -the outer Δ′ the solve published plus the per-surface layer parameters, and the inner-layer -model passed to `solve` computes the layer response at the prescribed rotation. Solving it -returns a new result with the closure changed from `:ideal` to `:matched` and the -eigenfunctions replaced — the expensive outer solve is reused, so layer-parameter scans -cost one cheap match solve per point: +A `MatchProblem` holds the outer Δ′ a solve published plus per-surface layer parameters, and the inner-layer model passed to `solve` computes the layer response at the prescribed rotation; the result is a new `ForceFreeStatesResult` with the closure changed from `:ideal` to `:matched` and the eigenfunctions replaced, so layer-parameter scans cost one cheap match per point: ```julia ffs = solve(eq, Galerkin(; rpec_flag=true, cut_solution=true); nn=1) - matched = solve(MatchProblem(ffs; eta=1e-6, rho=1e-7, rotation=0.0), GGJModel()) -@assert matched.closure === :matched - -# A rotation scan reuses the one outer solve: bpens = [solve(MatchProblem(ffs; eta=1e-6, rho=1e-7, rotation=f), GGJModel()).bpen for f in 0.0:50.0:500.0] ``` -A scalar applies to every surface; a vector gives one value per entry of `ffs.surfaces`, -core to edge. The problem needs a Δ′ payload with coil-response columns, so it accepts Galerkin (`rpec_flag=true`) and -Riccati (`vac_flag=true`) results; a Riccati-fed match fills `bpen` and the resonant data but keeps -`solution === nothing` (no outer basis is retained). Only a closure-capable model is -accepted — `GGJModel()` today; `SLAYERModel()` is slab-only and drives the free-eigenvalue tearing -solve instead. - -A deck poses the match with a `[Match]` section (`eta`, `rho`, `rotation`, `gamma`, `ideal`, -`inner_model = "GGJ"`), and `[GGJ]` holds the `GGJModel` keywords (`solver` and its options). +A scalar applies to every surface and a vector gives one value per entry of `ffs.surfaces`, core to edge. +The problem needs a Δ′ payload with coil-response columns, so it accepts Galerkin (`rpec_flag=true`) and Riccati (`vac_flag=true`) results; a Riccati-fed match fills `bpen` and the resonant data but keeps `solution === nothing`. +`GGJModel()` is the closure-capable model; `SLAYERModel()` is slab-only and drives the tearing solve instead. ## Tearing stability -A `TearingProblem` holds the outer solution and the surface coupling. The plasma is the -equilibrium's kinetic data (profiles, species, transport defaults, resistivity closure), the -inner-layer model holds its own closures, and `TearingConfig` holds the scan and root-search -settings. +A `TearingProblem` holds the outer solution and the surface coupling, the inner-layer model its own closures, and `TearingConfig` the scan and root-search settings: ```julia tear = solve(TearingProblem(ffs), SLAYERModel(; P_perp_model=:chi_perp_e); config=TearingConfig(; Q_re_range=(-2.0, 2.0), Q_im_range=(-0.5, 3.0))) ggj = solve(TearingProblem(ffs; coupling_mode=:coupled), GGJModel(; solver=:ray); config=TearingConfig(; gamma_range=(1.0, 1e4), omega_range=(-1e4, 1e4))) -save_to_hdf5("gpec.h5", tear) ``` -GGJ needs the raw Δ′ matrix and validates the surface ordering. Its η and ρ default from the -kinetic data: η by the resistivity closure SLAYER uses, ρ = μ·m_p·n_e as SLAYER takes it; pass `eta` -and `rho` to `TearingProblem` to override them. Both models report lab-frame frequencies with -per-surface E×B shifts: a coupled solve evaluates each layer in its own E×B frame (GGJ's s_k = s + i·n·Ω_E, -Johnson & Greene 1967 Appendix B), and an uncoupled root is moved from its layer's frame by n·Ω_E. Ω_E comes -from the kinetic data's `omega_E`, or from `omega_E_kHz` for the m/n it lists. Both expose -`omega` in rad/s and `gamma` in 1/s, positive for instability, and the HDF5 roots group records -`frequency_frame`. SLAYER searches a normalized Q box, GGJ a physical (γ, ω) box. - -A deck runs the solve with a `[Tearing]` section: the `TearingProblem` keywords and -`inner_model = "SLAYER"` or `"GGJ"`, with `[Tearing.search]` for the `TearingConfig` keys and -`[Tearing.CriticalResonantField]` for the penetration threshold. The model's keywords sit in -`[SLAYER]` or `[GGJ]`; `[SLAYER]` also gives the layer-overlap cap its closures. - -Kinetic runs take `Forward()`. The problem's `kinetic` keyword selects the source: -`kinetic=KineticForcesControl(; ...)` evaluates the KineticForces drift-kinetic -matrices with those settings on the equilibrium's kinetic profiles and species, and -`kinetic=TestKineticMatrices(; scale)` uses synthetic test matrices. `kinetic_factor` (default 1) -scales either. `km = kinetic_matrices(prob)` evaluates the source once, and `kinetic=km` reuses those -[`KineticMatrices`](@ref GeneralizedPerturbedEquilibrium.ForceFreeStates.KineticMatrices) in later solves, for example a -`kinetic_factor` scan. A deck selects the KineticForces model with `[ForceFreeStates] kinetic = "KineticForces"`, using -its `[KineticForces]` section. +Both models report lab-frame frequencies (`omega` in rad/s, `gamma` in 1/s, positive for instability, with `frequency_frame` recorded on the roots group); SLAYER searches a normalized Q box and GGJ a physical (γ, ω) box. ## Local stability -The Mercier, resistive interchange and ballooning profiles are their own stage, posed on the -equilibrium. A [`LocalStabilityProblem`](@ref) holds the equilibrium, ballooning angle -`theta_k`, and `alpha_boundary_scan` choice. [`Ballooning`](@ref) holds numerical tolerances -and scan settings. `solve` returns a [`LocalStabilityResult`](@ref); its `verbose` keyword -defaults to `false`. A deck states the problem in `[LocalStability]` and the numerics in `[Ballooning]`. - ```julia problem = LocalStabilityProblem(eq; theta_k=0.0, alpha_boundary_scan=true) ls = solve(problem, Ballooning(; n_scan=64); verbose=false) ls.psi, ls.D_I, ls.D_R, ls.ballooning_Delta_prime # profiles on the equilibrium ψ grid -ls.ballooning_psi, ls.alpha, ls.alpha_critical # first ballooning boundary, NaN outside its scan window +ls.ballooning_psi, ls.alpha, ls.alpha_critical # the ballooning boundary, NaN outside its scan window ``` ## Entry points @@ -226,7 +136,7 @@ GeneralizedPerturbedEquilibrium.perturbed_equilibrium [`plasma_response`](@ref GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.plasma_response), which builds the [`PlasmaResponse`](@ref GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.PlasmaResponse) that `perturbed_equilibrium` -applies, is documented with the [perturbed equilibrium](perturbed_equilibrium.md) module. +applies, is documented with the [Perturbed Equilibrium](perturbed_equilibrium.md) module. ## Integrator selectors diff --git a/docs/src/assets/density_decomposition.png b/docs/src/assets/density_decomposition.png deleted file mode 100644 index 8ae2c7659b..0000000000 Binary files a/docs/src/assets/density_decomposition.png and /dev/null differ diff --git a/docs/src/assets/ntv_cumulative_torque.png b/docs/src/assets/ntv_cumulative_torque.png deleted file mode 100644 index 1d86252550..0000000000 Binary files a/docs/src/assets/ntv_cumulative_torque.png and /dev/null differ diff --git a/docs/src/assets/ntv_eval_clustering.png b/docs/src/assets/ntv_eval_clustering.png deleted file mode 100644 index b2276a57fd..0000000000 Binary files a/docs/src/assets/ntv_eval_clustering.png and /dev/null differ diff --git a/docs/src/assets/ntv_torque_density.png b/docs/src/assets/ntv_torque_density.png deleted file mode 100644 index 392ee06955..0000000000 Binary files a/docs/src/assets/ntv_torque_density.png and /dev/null differ diff --git a/docs/src/ballooning.md b/docs/src/ballooning.md index 8ef6cd10bc..0adcb01fe5 100644 --- a/docs/src/ballooning.md +++ b/docs/src/ballooning.md @@ -1,918 +1,77 @@ -# Ballooning and Mercier Local Stability +# Local Stability -This note describes the physical and mathematical foundations of `Bal.jl`, which -computes local Mercier and high-``n`` ballooning stability. It covers the -ballooning-mode equation, the asymptotic boundary treatment, and the local -perturbation framework used to generate the ``s``-``\alpha`` diagram. +`LocalStability` evaluates the local, high-``n`` ideal and resistive criteria on every flux surface of the equilibrium: the Mercier interchange index ``D_I``, the Glasser–Greene–Johnson resistive-interchange index ``D_R`` and the ballooning stability index, with an optional scan for the critical pressure gradient. +It runs as its own stage, after `ForceFreeStates` and on the same equilibrium, and is the local counterpart of the global analysis on the [Stability](stability.md) page. -## Ballooning-mode equation +## What it computes -### Ballooning representation +- **Mercier criterion** ``D_I``: the ideal interchange index of each surface. A surface is Mercier unstable where ``D_I > 0``; the same number sets the exponents ``-1/2 \pm \sqrt{-D_I}`` of the large and small solutions that the global integration matches at every rational surface. +- **Resistive interchange** ``D_R = D_I + (H - 1/2)^2``, with ``H`` from the flux-surface averages of [Glasser, Greene and Johnson (1975)](citations.md#LocalStability). ``D_R > 0`` indicates resistive interchange instability. +- **Ballooning stability** ``\Delta'``: the high-``n`` ballooning index, obtained by integrating the marginal ballooning equation along the field line and matching the solution to its asymptotic forms at both ends ([Connor, Hastie and Taylor (1978)](citations.md#LocalStability)). It is distinct from the tearing Δ′ of `SingularSurfaces/`; ``\Delta' = 0`` is the marginal ballooning boundary. +- **Critical pressure gradient**: with `alpha_boundary_scan = true` the pressure gradient of each surface is scaled until the ballooning boundary is reached, giving the normalized gradient ``\alpha`` and its critical value ``\alpha_\mathrm{crit}``, the s–α diagram collapsed to a profile. -In a tokamak equilibrium, physical quantities are periodic in the poloidal angle -``\theta`` and toroidal angle ``\zeta``. The perturbation function ``\varphi`` -satisfies +## Governing equations -```math -\varphi(\psi, \theta, \beta) -= \varphi(\psi, \theta, \beta + 1) -= \varphi(\psi, \theta + 1, \beta - q), -``` - -where ``q`` is the safety factor. - -Ballooning modes have long wavelengths along magnetic field lines and short -wavelengths perpendicular to the field. They are represented using a -covering-space function ``\hat{\varphi}``, defined on -``-\infty < \theta < \infty``: - -```math -\varphi(\psi, \theta, \beta) -= \sum_{l=-\infty}^{\infty} -\hat{\varphi}(\psi, \theta - l) -\exp[-2\pi i n(\beta + lq)]. -``` - -Here - -```math -\beta = \zeta - q\theta + \int q'\theta_0(\psi)\,d\psi -``` - -is the Clebsch coordinate. The wave-number vector contains a secular term from -the magnetic shear ``q'``: - -```math -|\nabla\beta|^2 -= |\nabla\zeta - q\nabla\theta|^2 -+ q'^2(\theta-\theta_0)^2|\nabla\psi|^2 -- 2q'(\theta-\theta_0)(\nabla\zeta-q\nabla\theta)\cdot\nabla\psi. -``` - -### Ideal marginal ballooning equation - -At the stability boundary, ``\gamma = 0``. Variation of the ideal-MHD energy -principle gives the ballooning ODE - -```math -\mathbf{B}\cdot\nabla -\left( - \frac{|\nabla\beta|^2}{B^2} - \mathbf{B}\cdot\nabla\hat{\varphi} -\right) -+ 2\frac{P'}{\chi'}\kappa_w\hat{\varphi} -= 0, -``` - -where ``\chi`` is the poloidal flux and ``P'`` is the pressure-gradient term used -by the local-stability equations. Written using the Jacobian ``\mathcal{J}``, - -```math -\frac{\chi'}{\mathcal{J}}\frac{\partial}{\partial\theta} -\left( - \frac{|\nabla\beta|^2}{B^2} - \frac{\chi'}{\mathcal{J}} - \frac{\partial\hat{\varphi}}{\partial\theta} -\right) -+ 2\frac{P'}{\chi'}\kappa_w\hat{\varphi} -= 0. -``` - -`Bal.jl` converts this to a first-order matrix equation, - -```math -\frac{\partial\mathbf{u}}{\partial\theta} = \mathbf{L}\mathbf{u}, -\qquad -\mathbf{u} = -\begin{pmatrix} -\hat{\varphi} \\ -\dfrac{|\nabla\beta|^2}{B^2} -\dfrac{\chi'}{\mathcal{J}} -\dfrac{\partial\hat{\varphi}}{\partial\theta} -\end{pmatrix}, -``` - -with - -```math -\mathbf{L} -= \frac{\mathcal{J}}{\chi'} -\begin{pmatrix} -0 & \dfrac{B^2}{|\nabla\beta|^2} \\ --2\dfrac{P'}{\chi'}\kappa_w & 0 -\end{pmatrix}. -``` - -### Curvature decomposition - -The instability drive is the projected curvature - -```math -\kappa_w -\equiv -\frac{\boldsymbol{\kappa}\times\mathbf{B}}{B^2}\cdot\nabla\beta. -``` - -Using -``\nabla\beta = \nabla\zeta - q\nabla\theta - q'\theta\nabla\psi``, -the curvature separates into non-secular and secular pieces: - -```math -\kappa_w = \kappa_n - q'\theta\kappa_s, -``` - -where - -```math -\kappa_s -\equiv -\frac{\boldsymbol{\kappa}\times\mathbf{B}}{B^2}\cdot\nabla\psi -= -\frac{\pi f}{\mathcal{J}} -\frac{\partial}{\partial\theta} -\left(\frac{1}{B^2}\right). -``` - -This decomposition is the main structure used to assemble the ballooning -coefficient matrix. - -## Asymptotic basis and Mercier exponent - -In the ``\theta\to\infty`` limit, the coefficient matrix diverges because -``|\nabla\beta|^2 \simeq q'^2\theta^2|\nabla\psi|^2``. `Bal.jl` applies the -canonical transformation - -```math -\mathbf{R} = \operatorname{diag}(\theta^{-1/2}, \theta^{1/2}) -``` - -and solves - -```math -\frac{\partial\mathbf{v}}{\partial\theta} = \mathbf{M}\mathbf{v}, -\qquad -\mathbf{M} -= -\begin{pmatrix} -\dfrac{1}{2\theta} -& -\dfrac{\mathcal{J}}{\chi'} -\dfrac{B^2\theta}{|\nabla\beta|^2} -\\ --\dfrac{\mathcal{J}}{\chi'} -\dfrac{2P'\kappa_w}{\chi'\theta} -& --\dfrac{1}{2\theta} -\end{pmatrix}. -``` - -To keep the lowest-order matrix finite, the code applies - -```math -\mathbf{S} -= -\begin{pmatrix} -1 & 0 \\ -\sigma & 1 -\end{pmatrix}, -\qquad -\sigma = --\frac{2\pi f P'q'}{\chi'^2 B^2}. -``` - -The transformed matrices ``\mathbf{N}_k = \mathbf{S}^{-1}\mathbf{M}_k\mathbf{S}`` -enter the recurrence relation - -```math -\frac{\partial\mathbf{w}^{(k+1)}}{\partial\theta} -= -[\mathbf{N}_0-(\alpha+k)\mathbf{I}]\mathbf{w}^{(k)} -+ \sum_{l=1}^{k}\mathbf{N}_l\mathbf{w}^{(k-l)}. -``` - -The Mercier exponent ``\alpha`` follows from the zeroth-order solvability -condition: - -```math -\det[\langle\mathbf{N}_0\rangle-\alpha\mathbf{I}] = 0, -\qquad -\alpha = \pm\sqrt{-D_I}. -``` - -### Equivalence with the surface-average Mercier formula - -The previous `Mercier.jl` evaluated the Mercier -index from flux-surface averages. The two routes are mathematically equivalent: -they are different reductions of the same large-``\theta`` ballooning -solvability condition. - -Let the flux-surface average used by the old implementation be - -```math -\langle X\rangle_{\mathrm{fs}} -\equiv -\int_0^1 \frac{\mathcal{J}}{V'} X(\theta)\,d\theta, -``` - -and define - -```math -A_1 = \left\langle\frac{B^2}{|\nabla\psi|^2}\right\rangle_{\mathrm{fs}}, -\qquad -A_2 = \left\langle\frac{1}{|\nabla\psi|^2}\right\rangle_{\mathrm{fs}}, -``` - -```math -A_3 = \left\langle\frac{1}{B^2}\right\rangle_{\mathrm{fs}}, -\qquad -A_4 = \left\langle\frac{1}{B^2|\nabla\psi|^2}\right\rangle_{\mathrm{fs}}, -\qquad -A_5 = \left\langle B^2\right\rangle_{\mathrm{fs}}. -``` - -Below, an overbar denotes the period integral used to form the `Bal.jl` -``\bar{d}_0`` matrix: - -```math -\overline{X} \equiv \int_0^1 X(\theta)\,d\theta. -``` - -With ``F = 2\pi f``, the old code evaluated - -```math -T = -\frac{F P' V'}{q'\chi'^3} A_2, -``` - -and then - -```math -D_I^{\mathrm{avg}} -= --\frac14 -+ T(1-T) -+ P'\left(\frac{V'}{q'\chi'^2}\right)^2 A_1 -\left[ - P'\left(A_3+\frac{F^2}{\chi'^2}A_4\right) - - \frac{V''}{V'} -\right]. -``` - -This is the surface-average Mercier expression that appeared in the old -`Mercier.jl` and, temporarily, in the earlier `resistive_interchange` helper. -The associated GGJ resistive-interchange quantity used - -```math -H = -\frac{F P' V'}{q'\chi'^3} -\left( - A_2 - \frac{A_1}{A_5} -\right), -\qquad -D_R = D_I + \left(H-\frac12\right)^2. -``` - -The determinant route starts from the same ballooning equation, takes the -large-``\theta`` limit, applies the ``\mathbf{R}`` and ``\mathbf{S}`` -transformations, and obtains the periodic lowest-order system - -```math -\frac{\partial\mathbf{w}^{(1)}}{\partial\theta} -= -(\mathbf{N}_0-\alpha\mathbf{I})\mathbf{w}^{(0)}. -``` - -Solvability requires the mean operator to have a null vector: - -```math -[\overline{\mathbf{N}_0}-\alpha\mathbf{I}]\mathbf{w}^{(0)}=0. -``` - -The transformed matrix is traceless by construction, - -```math -\operatorname{tr}\overline{\mathbf{N}_0} = 0, -``` - -so the characteristic equation reduces to - -```math -\alpha^2 + \det\overline{\mathbf{N}_0} = 0. -``` - -Comparing this with the Mercier exponent definition -``\alpha^2 = -D_I`` gives - -```math -D_I = \det\overline{\mathbf{N}_0}. -``` - -The reduction can be written explicitly. For the unperturbed Mercier scan, -``\mathbf{N}_0`` has the code form - -```math -\mathbf{N}_0(\theta) -= -\begin{pmatrix} -\frac12 + M_{12}\sigma & M_{12} \\ -M_{21}-\sigma-M_{12}\sigma^2 & -\frac12-M_{12}\sigma -\end{pmatrix}, -``` - -where - -```math -M_{12} -= -\frac{\mathcal{J}}{\chi'} -\frac{B^2}{q'^2|\nabla\psi|^2}, -\qquad -M_{21} -= --2\frac{\mathcal{J}}{\chi'^2}P'\kappa_n, -\qquad -\sigma -= --\frac{F P' q'}{\chi'^2 B^2}. -``` - -The averaged matrix can therefore be written as - -```math -\bar{d}_0 -\equiv -\overline{\mathbf{N}_0} -= -\begin{pmatrix} -\frac12+B & A \\ -K & -\frac12-B -\end{pmatrix}, -``` - -with - -```math -A = \overline{M_{12}}, -\qquad -B = \overline{M_{12}\sigma}, -\qquad -K = \overline{M_{21}-\sigma-M_{12}\sigma^2}. -``` - -The first two averages are direct: - -```math -A -= -\frac{V'}{q'^2\chi'}A_1, -``` - -```math -B -= -\overline{ - \frac{\mathcal{J}}{\chi'} - \frac{B^2}{q'^2|\nabla\psi|^2} - \left( - -\frac{F P' q'}{\chi'^2B^2} - \right) -} -= --\frac{F P' V'}{q'\chi'^3}A_2 -= --T. -``` - -For ``K``, separate the three pieces: - -```math -K -= -\overline{M_{21}} --\overline{\sigma} --\overline{M_{12}\sigma^2}. -``` - -The last term is also direct: - -```math -\overline{M_{12}\sigma^2} -= -\overline{ - \frac{\mathcal{J}}{\chi'} - \frac{B^2}{q'^2|\nabla\psi|^2} - \frac{F^2P'^2q'^2}{\chi'^4B^4} -} -= -\frac{F^2P'^2V'}{\chi'^5}A_4. -``` - -The remaining curvature combination reduces by the axisymmetric equilibrium -identity used in the Mercier derivation: +The stage integrates the marginal ballooning equation of [Connor, Hastie and Taylor (1978)](citations.md#LocalStability) for the perturbation ``\hat{\varphi}`` along the field line on the covering space ``-\infty < \theta < \infty``, where the shear ``q'`` enters through ``|\nabla\beta|^2`` of the Clebsch angle ``\beta`` and the projected curvature ``\kappa_w = \kappa_n - q'\theta\,\kappa_s`` splits into a non-secular part and a secular part driven by the shear (``\chi``, ``P'`` and ``\mathcal{J}`` below are the poloidal flux, the pressure gradient and the Jacobian). +The code integrates it as the first-order system ``\partial_\theta \mathbf{u} = \mathbf{L}\mathbf{u}`` in ``\mathbf{u} = (\hat{\varphi},\ |\nabla\beta|^2 \chi' \partial_\theta\hat{\varphi}/(B^2\mathcal{J}))``. -```math -\int_0^1 -\left( - 2\mathcal{J}\kappa_n - - \frac{Fq'}{B^2} -\right)d\theta -= -\frac{V'}{\chi'} -\left( - P'A_3 - \frac{V''}{V'} -\right). -``` - -Using ``M_{21}=-2\mathcal{J}P'\kappa_n/\chi'^2`` and -``\sigma=-FP'q'/(\chi'^2B^2)``, this gives - -```math -\overline{M_{21}}-\overline{\sigma} -= --\frac{P'}{\chi'^2} -\int_0^1 -\left( - 2\mathcal{J}\kappa_n - - \frac{Fq'}{B^2} -\right)d\theta -= --\frac{P'V'}{\chi'^3} -\left( - P'A_3 - \frac{V''}{V'} -\right). -``` - -Therefore - -```math -K -= --\frac{P'V'}{\chi'^3} -\left[ - P'\left(A_3+\frac{F^2}{\chi'^2}A_4\right) - - \frac{V''}{V'} -\right]. -``` - -Now take the determinant: - -```math -\det\bar{d}_0 -= -\left(\frac12+B\right)\left(-\frac12-B\right)-AK. -``` - -Substituting ``B=-T`` gives the first part of the old Mercier expression: - -```math -\left(\frac12-T\right)\left(-\frac12+T\right) -= --\frac14+T(1-T). -``` - -Substituting ``A`` and ``K`` gives the second part: - -```math --AK -= -\frac{V'}{q'^2\chi'}A_1 -\frac{P'V'}{\chi'^3} -\left[ - P'\left(A_3+\frac{F^2}{\chi'^2}A_4\right) - - \frac{V''}{V'} -\right] -``` - -```math -= -P'\left(\frac{V'}{q'\chi'^2}\right)^2 A_1 -\left[ - P'\left(A_3+\frac{F^2}{\chi'^2}A_4\right) - - \frac{V''}{V'} -\right]. -``` - -Thus the determinant route gives exactly - -```math -\det\bar{d}_0 -= --\frac14 -+ T(1-T) -+ P'\left(\frac{V'}{q'\chi'^2}\right)^2 A_1 -\left[ - P'\left(A_3+\frac{F^2}{\chi'^2}A_4\right) - - \frac{V''}{V'} -\right] -= -D_I^{\mathrm{avg}}. -``` - -So, under the same normalization and exact equilibrium identities, - -```math -D_I^{\mathrm{det}} -\equiv -\det(\bar{d}_0) -= -D_I^{\mathrm{avg}}. -``` - -The implementation difference is numerical rather than physical. The old -surface-average route evaluates the closed Mercier formula directly, while the -new `Bal.jl` route constructs the same Mercier index from the averaged -large-``\theta`` local-mode matrix. In finite-precision calculations the two can -separate near the magnetic axis or other ill-conditioned surfaces because both -forms contain factors such as ``1/q'`` and ``1/q'^2``. Away from those regions, -they should agree to the accuracy of the equilibrium derivatives and surface -quadrature. - -The divergent big solution and decaying small solution have leading forms - -```math -\mathbf{u}_{\mathrm{large}}^{(0)} -\simeq -\begin{pmatrix} -\theta^{\alpha-1/2} \\ -\theta^{\alpha+1/2}(\sigma+w_{21}) -\end{pmatrix}, -\qquad -\mathbf{u}_{\mathrm{small}}^{(0)} -\simeq -\begin{pmatrix} -\theta^{-\alpha-1/2} \\ -\theta^{-\alpha+1/2}(\sigma+w_{22}) -\end{pmatrix}. -``` - -## First-order asymptotic matching - -The zeroth-order large-``\theta`` solution alone is not accurate enough at a -finite integration boundary ``\theta_{\max}``. `Bal.jl` therefore includes the -first-order correction ``\mathbf{w}^{(1)}``: - -```math -\frac{\partial\mathbf{w}^{(1)}}{\partial\theta} -= -(\mathbf{N}_0-\alpha\mathbf{I})\mathbf{w}^{(0)} -= -\tilde{\mathbf{N}}_0(\theta)\mathbf{w}^{(0)}. -``` - -The first-order solution is written as - -```math -\mathbf{w}^{(1)}(\theta) -= -\left(\int\tilde{\mathbf{N}}_0(\theta)\,d\theta\right)\mathbf{w}^{(0)} -+ \mathbf{c}^{(1)} -\equiv -\tilde{\mathbf{w}}^{(1)}(\theta) + \mathbf{c}^{(1)}. -``` - -The constant vector is fixed by the first-order solvability condition, - -```math -[\langle\mathbf{N}_0\rangle-(\alpha+1)\mathbf{I}]\mathbf{c}^{(1)} -= -- -\left\langle - [\mathbf{N}_0-(\alpha+1)\mathbf{I}] - \tilde{\mathbf{w}}^{(1)}(\theta) - + \mathbf{N}_1\mathbf{w}^{(0)} -\right\rangle. -``` - -This gives the asymptotic boundary form - -```math -\mathbf{u}_{\mathrm{boundary}} -\simeq -\mathbf{R}\mathbf{S} -\left( - \mathbf{w}^{(0)} - + \frac{1}{\theta}\mathbf{w}^{(1)}(\theta) -\right). -``` - -## Local thin-layer ordering - -The ``s``-``\alpha`` diagram scans local instability near a flux surface -``\psi_b`` without recomputing the full equilibrium. A thin-layer parameter -``0 < \mu \ll 1`` is introduced through - -```math -s = \frac{\psi-\psi_b}{\mu}, -\qquad -\frac{\partial}{\partial\psi} -= -\frac{1}{\mu}\frac{\partial}{\partial s}. -``` - -Profiles are expanded as background plus local perturbation, - -```math -P(\psi) = P^{(0)}(\psi) + \mu P^{(1)}(s) + \cdots, -\qquad -q(\psi) = q^{(0)}(\psi) + \mu q^{(1)}(s) + \cdots. -``` - -After differentiation, the ``\mu`` amplitude cancels the ``1/\mu`` derivative -scale. Thus the derivative-level perturbations are order unity: - -```math -P_\psi = P_\psi^{(0)} + P_s^{(1)}, -\qquad -q_\psi = q_\psi^{(0)} + q_s^{(1)}. -``` - -The subscript ``s`` therefore represents the promoted local derivative -perturbation used as the practical scan variable in the ``s``-``\alpha`` diagram. -Retained first-order contributions are denoted by a superscript ``[1]``; for -example, ``I_{\mathrm{per}}^{[1]}`` is the retained local periodic variation of -the integrated shear. - -## Local ``P'`` and ``q'`` perturbations - -Using the integrated-shear function - -```math -I_{\mathrm{tot}} = q'(\theta-\theta_k) + I_{\mathrm{per}}, -``` - -the wave-number space and curvature are written - -```math -\left(\frac{|\nabla\beta|}{B}\right)^2 -= -\frac{1}{\chi'^2|\nabla\psi|^2} -+ \frac{|\nabla\psi|^2}{B^2}I_{\mathrm{tot}}^2, -``` - -```math -\kappa_w -= -\frac{\boldsymbol{\kappa}\cdot\nabla\psi} -{\chi'|\nabla\psi|^2} -- \kappa_s I_{\mathrm{tot}}. -``` - -With - -```math -I_{\mathrm{tot}} -= I_{\mathrm{tot}}^{(0)} + I_{\mathrm{tot}}^{[1]}, -\qquad -I_{\mathrm{tot}}^{[1]} -= q_s^{(1)}(\theta-\theta_k) + I_{\mathrm{per}}^{[1]}, -``` - -the retained order-unity perturbation terms are - -```math -\left(\frac{|\nabla\beta|^2}{B^2}\right)_{O(1),\mathrm{ret}} -= -\frac{|\nabla\psi|_0^2}{B_0^2} -\left[ - 2I_{\mathrm{tot}}^{(0)}I_{\mathrm{tot}}^{[1]} - + \left(I_{\mathrm{tot}}^{[1]}\right)^2 -\right], -``` - -```math -\kappa_w^{[1]} -= --\kappa_s^{(0)}I_{\mathrm{tot}}^{[1]}. -``` - -The cross term and nonlinear ``(I_{\mathrm{tot}}^{[1]})^2`` term are retained in -the matrix assembly rather than discarded by linearization. - -## ``\mathbf{N}_0`` matrix construction - -The lower-triangular element of the canonical transformation combines the local -``P'`` and ``q'`` perturbations: - -```math -\sigma -= --\frac{2\pi f}{\chi_0'^2B_0^2} -\left(P_\psi^{(0)} + P_s^{(1)}\right) -\left(q_\psi^{(0)} + q_s^{(1)}\right). -``` - -The ``(1,2)`` component of the asymptotic matrix responds to ``q_s^{(1)}``: - -```math -(\mathbf{M}_0)_{12} -= -\frac{\mathcal{J}_0}{\chi_0'} -\frac{B_0^2} -{\left(q_\psi^{(0)}+q_s^{(1)}\right)^2|\nabla\psi|_0^2}. -``` - -The ``(2,1)`` component uses the expanded non-secular curvature, - -```math -\kappa_n(\theta) -= -\kappa_n^{(0)}(\theta) -+ \kappa_s^{(0)}(\theta)q_s^{(1)}\theta_k -- \kappa_s^{(0)}(\theta)I_{\mathrm{per}}^{[1]}(\theta), -``` - -which gives - -```math -(\mathbf{M}_0)_{21} -= --\frac{2\mathcal{J}_0}{\chi_0'^2} -\left(P_\psi^{(0)} + P_s^{(1)}\right) -\left[ - \kappa_n^{(0)}(\theta) - + \kappa_s^{(0)}(\theta)q_s^{(1)}\theta_k - - \kappa_s^{(0)}(\theta)I_{\mathrm{per}}^{[1]}(\theta) -\right]. -``` - -The final matrix -``\mathbf{N}_0(\theta)=\mathbf{S}^{-1}\mathbf{M}_0(\theta)\mathbf{S}`` is -period-averaged to form ``\langle\mathbf{N}_0\rangle`` in the Mercier solvability -equation. - -## ``\mathbf{N}_1`` correction - -The first-order asymptotic boundary condition requires the -``O(\theta^{-2})`` matrix ``\mathbf{N}_1``. In the large-``\theta`` limit, only the -upper-right component of ``\mathbf{M}_1`` remains: - -```math -\mathbf{M}_1(\theta) -= -\begin{pmatrix} -0 & a_1(\theta) \\ -0 & 0 -\end{pmatrix}, -``` - -where +As ``\theta \to \infty`` the coefficients diverge like ``q'^2\theta^2|\nabla\psi|^2``. +The canonical transformation ``\mathbf{R} = \operatorname{diag}(\theta^{-1/2}, \theta^{1/2})`` followed by a shear ``\mathbf{S}`` with ``\sigma = -2\pi f P' q'/(\chi'^2 B^2)`` makes the leading-order matrix ``\mathbf{N}_0`` finite, and the period average of that matrix fixes the asymptotic exponent: ```math -a_1(\theta) -= --2\frac{\mathcal{J}_0}{\chi_0'} -\frac{B_0^2} -{|\nabla\psi|_0^2\left(q_\psi^{(0)}+q_s^{(1)}\right)^3} -\left[ - I_{\mathrm{per}}^{(0)}(\theta) - + I_{\mathrm{per}}^{[1]}(\theta) - - \left(q_\psi^{(0)}+q_s^{(1)}\right)\theta_k -\right]. +\det[\langle\mathbf{N}_0\rangle - \alpha\mathbf{I}] = 0, \qquad \alpha = \pm\sqrt{-D_I} . ``` -Applying the ``\mathbf{S}^{-1}`` transformation gives +This is the Mercier index: ``D_I = -\alpha^2`` evaluated from ``\langle\mathbf{N}_0\rangle`` equals the classical surface-average Mercier expression, and ``D_I > 0`` (complex exponents) is interchange instability. +The large and small solutions behave as ``\theta^{\alpha - 1/2}`` and ``\theta^{-\alpha - 1/2}``; the boundary condition at the finite integration limit ``\theta_\max`` includes the first-order correction of the asymptotic series, and the ballooning index is the jump in the logarithmic derivative of the solution between the two ends. -```math -\mathbf{N}_1(\theta) -= -a_1(\theta) -\begin{pmatrix} -\sigma & 1 \\ --\sigma^2 & -\sigma -\end{pmatrix}. -``` +The s–α diagram scans the local stability of one surface without recomputing the equilibrium: a thin-layer ordering promotes local perturbations of the pressure gradient ``P'_s`` and shear ``q'_s`` to order unity in the derivatives, and the integrated shear ``I_\mathrm{tot} = q'(\theta - \theta_k) + I_\mathrm{per}`` carries them into ``|\nabla\beta|^2`` and ``\kappa_w``. +The cross term and the quadratic term in the perturbed shear are kept in the matrix assembly, and the ballooning angle ``\theta_k`` is the gauge where the total accumulated shear vanishes, ``I_\mathrm{tot}(\theta_k) = 0``. -Thus ``\mathbf{N}_1`` is closed in terms of -``I_{\mathrm{per}}^{[1]}``, ``P_s^{(1)}``, and ``q_s^{(1)}``. The -``I_{\mathrm{per}}^{[1]}`` and ``q_s^{(1)}`` terms control ``a_1(\theta)``, while -``P_s^{(1)}`` enters through ``\sigma``. +## Running it -## Practical implementation in `Bal.jl` +A `[LocalStability]` section runs the stage; it needs `[ForceFreeStates]` in the same deck. +Its keys are the problem settings, and the optional `[Ballooning]` section holds the [`Ballooning`](@ref) numerical settings, as in `examples/DIIID-like_ideal_example/gpec.toml`: -### Geometric curvature formula +```toml +[LocalStability] +theta_k = 0.0 # ballooning angle +alpha_boundary_scan = true # also find the critical alpha of each surface -The non-secular curvature calculation uses the geometric expression - -```math -\kappa_n^{(0)} -= --\frac{1}{\mathcal{J}_0B_0} -\frac{\partial}{\partial\psi} -\left(\frac{\mathcal{J}_0B_0}{\chi'_0}\right) -+ -\frac{1}{\mathcal{J}_0B_0} -\frac{\partial}{\partial\theta} -\left( - \mathcal{J}_0 - \frac{\mathbf{B}_0\cdot\nabla\theta\times\nabla\zeta}{B_0} -\right) -+ -\frac{2\pi f q'_0}{B_0^2\mathcal{J}_0}. +[Ballooning] +ode_tol = 1e-8 +n_scan = 24 # alpha scan points +max_alpha_scale = 8.0 +alpha_tol = 1e-3 +psi_window = [0.1, 0.99] # surfaces scanned for the alpha boundary ``` -This avoids relying entirely on Grad-Shafranov identities, which can amplify -numerical equilibrium errors near the edge or in other regions where the input -equilibrium does not satisfy the Grad-Shafranov equation exactly. - -### Two-loop structure - -The code separates the calculation into two loops to avoid repeatedly -differentiating noisy geometric quantities. - -- **Loop 1:** Use only unperturbed background quantities, such as ``P'`` and - ``q'``, to precompute the base structural curvature ``\kappa_n^{(0)}``. -- **Loop 2:** Insert the requested local perturbations, such as - ``P_s^{(1)}``, ``q_s^{(1)}``, and ``I_{\mathrm{per}}^{[1]}``. The code adds only - the analytically derived algebraic perturbation - ``\kappa_s^{(0)}q_s^{(1)}\theta_k - - \kappa_s^{(0)}I_{\mathrm{per}}^{[1]}`` to the precomputed - ``\kappa_n^{(0)}``. - -This preserves the perturbation structure while avoiding unstable repeated -geometric differentiation. - -### ``\theta_{\max}`` and `Baloo.f`-style matching - -For large ``\theta``, the wave-number factor behaves like a quadratic polynomial, - -```math -a_2\theta^2 + a_1\theta + a_0 -= -a_2\left(\theta+\frac{a_1}{2a_2}\right)^2 -+ \left(a_0-\frac{a_1^2}{4a_2}\right). -``` +A script runs the same stage with `solve(LocalStabilityProblem(eq; theta_k=0.0, alpha_boundary_scan=true), Ballooning(); verbose=false)`; see the [Scripting API](api.md). +`examples/DIIID-like_ideal_example/analyze_example.jl` plots the ``s`` and ``\alpha`` profiles, the ``D_I`` and ballooning ``\Delta'`` profiles, and the two-dimensional s–α maps with their zero contours, using `salpha_reference`, `compute_ballooning_stability!` and `scan_delta_prime_map`. -The residual scale is ``c = a_0-a_1^2/(4a_2)``. To make the asymptotic -approximation valid, the code chooses an integration boundary where the -``\theta^2`` term dominates, using the scale ``\sqrt{|c/a_2|}``. -Because the shooting boundary is a single flux-surface quantity, `Bal.jl` -evaluates this scale on the ballooning theta grid and uses the largest finite -value on that surface. The base integration boundary is then capped as -``\theta_{\max}=\min(16.5,10\max_\theta\sqrt{|c/a_2|})``. +## Outputs -Directly extracting ``c_{a1}`` or ``\Delta'`` from the small asymptotic solution is -numerically fragile. The comparison below shows that the `Baloo.f` style, which -sets the displacement to vanish at ``\pm\theta_{\max}`` and computes ``\Delta'`` -from the resulting solution, gives the most stable continuous profile. - -![Comparison of numerical stability among methods for extracting ``\Delta'`` and ``c_{a1}``.](figures/ballooning/delta_prime_matching.png) - -`Bal.jl` therefore adopts the `Baloo.f`-style boundary replacement for the -ballooning ``\Delta'`` calculation. - -### Ballooning angle and integral-gauge correction - -The ballooning angle ``\theta_k`` is the integration constant specifying where the -total accumulated local shear is zero: - -```math -I_{\mathrm{tot}}(\theta_k) = 0. -``` - -In `Bal.jl`, - -```math -I_{\mathrm{tot}} -= -\left(q'^{(0)}+q_s^{(1)}\right)(\theta-\theta_k) -+ \left(I_{\mathrm{per}}^{(0)}(\theta) - + I_{\mathrm{per}}^{[1]}(\theta)\right). -``` - -To satisfy the gauge exactly, the periodic shear functions are shifted so that -they vanish at ``\theta_k``: - -```math -I_{\mathrm{per}}^{(0)}(\theta) -\leftarrow -I_{\mathrm{per}}^{(0)}(\theta) -- I_{\mathrm{per}}^{(0)}(\theta_k), -``` - -```math -I_{\mathrm{per}}^{[1]}(\theta) -\leftarrow -I_{\mathrm{per}}^{[1]}(\theta) -- I_{\mathrm{per}}^{[1]}(\theta_k). -``` +| Dataset | Meaning | +|---|---| +| `LocalStability/psi` | the flux surfaces, ψ_N | +| `LocalStability/D_I` | Mercier index; ``D_I > 0`` is interchange unstable | +| `LocalStability/D_R` | resistive interchange index; ``D_R > 0`` is unstable | +| `LocalStability/ballooning_Delta_prime` | ballooning stability index along `psi`; ``\Delta' = 0`` is marginal | +| `LocalStability/ballooning_psi`, `alpha`, `alpha_critical` | with `alpha_boundary_scan = true`: the scanned surfaces, their normalized pressure gradient and its critical value | -This aligns the numerical implementation with the theoretical gauge used in the -ballooning formulation. +`SingularSurfaces/D_I` repeats the Mercier index at the rational surfaces of the global solve. -### `locstab_fs` outputs +## Numerical notes -`Bal.jl` now computes both the Mercier criterion and the ballooning ``\Delta'``. -The previous standalone `Mercier.jl` path can therefore be removed from the local -stability calculation. The comparison below shows agreement between the previous -Mercier calculation and the new `Bal.jl` calculation. +- The non-secular curvature ``\kappa_n`` is evaluated from the geometric expression in ``\mathcal{J}B`` and its ``\psi`` and ``\theta`` derivatives rather than from Grad–Shafranov identities, which amplify equilibrium errors near the edge. +- Background quantities are differentiated once, in a first loop; the s–α perturbations are then added algebraically, so noisy geometric derivatives are never repeated. +- The integration limit is set where the quadratic term of the wave-number factor dominates, ``\theta_\max = \min(16.5,\ 10\max_\theta\sqrt{|c/a_2|})``, and the ballooning index is taken from a solution forced to vanish at ``\pm\theta_\max`` rather than from the coefficient of the small solution, which is numerically fragile: -![Comparison of Mercier criterion calculations from `Mercier.jl` and `Bal.jl`.](figures/ballooning/mercier_comparison.png) +![Comparison of the numerical stability of the methods for extracting the ballooning index.](figures/ballooning/delta_prime_matching.png) -The local-stability output now stores ballooning ``\Delta'`` in the fourth -`locstab_fs` entry. In the HDF5 output this is written as -`LocalStability/ballooning_Delta_prime`, distinct from the tearing ``\Delta'`` outputs -under `SingularSurfaces/` and `PerturbedEquilibrium/SingularCoupling/`. +- The periodic shear functions are shifted to vanish at ``\theta_k`` so the gauge ``I_\mathrm{tot}(\theta_k) = 0`` holds exactly. -## API +## API reference ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.LocalStability] diff --git a/docs/src/citations.md b/docs/src/citations.md index 9bf7f5a03e..014a258393 100644 --- a/docs/src/citations.md +++ b/docs/src/citations.md @@ -1,10 +1,10 @@ # Citations -The papers below provide the theoretical foundations for GPEC's algorithms. When working with specific modules, refer to the relevant papers for derivations, equation numbering, and physical context. PDFs are available locally in `docs/resources/`. +The papers below provide the theoretical foundations for GPEC's algorithms; refer to them for derivations, equation numbering and physical context. PDFs of most of them are in `docs/resources/` of the repository. --- -## Vacuum Module +## Vacuum > M. S. Chance, "Vacuum calculations in azimuthally symmetric geometry," > *Physics of Plasmas* **4**, 2161 (1997). @@ -50,7 +50,7 @@ Applies the scheme to the singular surface integrals of magnetic fusion (exterio --- -## ForceFreeStates Module +## ForceFreeStates > A. H. Glasser, "The direct criterion of Newcomb for the ideal MHD stability of an axisymmetric toroidal plasma," > *Physics of Plasmas* **23**, 072505 (2016). @@ -64,11 +64,11 @@ The primary reference for the `ForceFreeStates` module. Derives the Euler-Lagran > *Physics of Plasmas* **25**, 032507 (2018). > DOI: [10.1063/1.5007042](https://doi.org/10.1063/1.5007042) -Reformulates the DCON eigenvalue problem as a Riccati matrix ODE, enabling parallel integration across singular surfaces and faster computation. Implemented in `src/ForceFreeStates/Riccati/` and enabled via `integrator = "riccati"` in `[ForceFreeStates]`. +Reformulates the DCON eigenvalue problem as a Riccati matrix ODE, enabling parallel integration across singular surfaces and faster computation. Implemented in `src/ForceFreeStates/Riccati/` and enabled via `integrator = "Riccati"` in `[ForceFreeStates]`. --- -## PerturbedEquilibrium Module +## PerturbedEquilibrium > J.-K. Park, A. H. Boozer, and A. H. Glasser, "Computation of three-dimensional tokamak and spherical torus equilibria," > *Physics of Plasmas* **14**, 052110 (2007). @@ -110,7 +110,23 @@ Establishes the area weightings that make flux-surface Fourier spectra of pertur --- -## InnerLayer Module +## LocalStability + +> A. H. Glasser, J. M. Greene, and J. L. Johnson, "Resistive instabilities in general toroidal plasma configurations," +> *Physics of Fluids* **18**, 875 (1975). + +The resistive interchange criterion ``D_R`` and the layer coefficients E, F, G, H, K, M of each rational surface, evaluated by `LocalStability` and stored under `SingularSurfaces/`. + +--- + +> J. W. Connor, R. J. Hastie, and J. B. Taylor, "Shear, periodicity, and plasma ballooning modes," +> *Physical Review Letters* **40**, 396 (1978). + +The ballooning representation on the covering space and the marginal ballooning equation that `LocalStability` integrates along the field line. + +--- + +## Tearing > A. H. Glasser, Z. R. Wang, and J.-K. Park, "Computation of resistive instabilities by matched asymptotic expansions," > *Physics of Plasmas* **23**, 112506 (2016). @@ -148,7 +164,7 @@ Demonstrates asymptotic matching for resistive plasma response in realistic toro > *Physics of Plasmas* **29**, 072506 (2022). > DOI: [10.1063/5.0093079](https://doi.org/10.1063/5.0093079) -The basis for the presently implemented in Fortran SLAYER code which computes the inner layer response in a two-fluid slab layer model. +The SLAYER solver implemented in `InnerLayer.SLAYER`: the two-fluid slab layer equations of Cole and Fitzpatrick (2006) with a few terms dropped, solved for the inner-layer response of each rational surface. --- @@ -156,7 +172,7 @@ The basis for the presently implemented in Fortran SLAYER code which computes th > *Physics of Plasmas* **13**, 032503 (2006). > DOI: [10.1063/1.2178167](https://doi.org/10.1063/1.2178167) -The torque-balance error-field penetration threshold (Eqs. 61–62) behind the critical resonant field computed in `Tearing.CriticalResonantField`. +The two-fluid drift-MHD layer equations that SLAYER solves in reduced form, and the torque-balance error-field penetration threshold (Eqs. 61–62) behind the critical resonant field computed in `Tearing.CriticalResonantField`. --- @@ -169,10 +185,10 @@ Combines the toroidal outer-region calculation with the two-fluid slab layer mod --- > R. Fitzpatrick, "Response of a magnetically diverted tokamak plasma to a resonant magnetic perturbation," -> *Nuclear Fusion* (2025). +> *Nuclear Fusion* **66**, 046025 (2026). > DOI: [10.1088/1741-4326/ae4fdd](https://doi.org/10.1088/1741-4326/ae4fdd) (open access: arXiv:2511.07666) -Gives the diffusive-resistive layer width (Eq. 100) and, in Sect. 5.9, restricts the response calculation to the region inside the point where adjacent resonant layers first overlap. Both are implemented by the layer-overlap scan in `src/Tearing/LayerOverlap.jl`, which feeds the `psilim_from_layer_overlap` cap described in [Integration domain truncation](@ref). +Gives the diffusive-resistive layer width (Eq. 100) and, in Sect. 5.9, restricts the response calculation to the region inside the point where adjacent resonant layers first overlap. Both are implemented by the layer-overlap scan in `src/Tearing/LayerOverlap.jl`, which feeds the `psilim_from_layer_overlap` cap described in [Integration domain](@ref). --- @@ -191,7 +207,7 @@ The toroidal (Eq. 59) and large-aspect-ratio (Eq. 61) critical-Δ offsets from f --- -## ErrorFields Module +## ErrorFields > N. C. Logan et al., "Robustness of the tokamak error field correction tolerance scaling," > *Plasma Physics and Controlled Fusion* **62**, 084001 (2020). @@ -223,7 +239,23 @@ module computes (Eqs. 1–5). --- -## Kinetic Forces +> N. C. Logan et al., "SPARC tokamak error field expectations and physics-based correction coil design," +> *Nuclear Fusion* **66**, 066003 (2026). +> DOI: [10.1088/1741-4326/ae6086](https://doi.org/10.1088/1741-4326/ae6086) + +The NTV-limited error-field correction model of `ErrorFields.NTV`: the correction current is capped by the torque budget the non-resonant field of the correction array spends. + +--- + +> N. Leuthold et al., "ARC physics basis–magnetohydrodynamics," +> *Journal of Plasma Physics* **92**, E49 (2026). +> DOI: [10.1017/S0022377826101421](https://doi.org/10.1017/S0022377826101421) + +The same model applied to ARC, with the zero-rotation `torque_only` limit reported beside the residual-torque one. + +--- + +## KineticForces The following papers develop the kinetic-force and neoclassical toroidal viscosity (NTV) theory underpinning GPEC's kinetic analysis path — the energy principle with kinetic effects, the self-consistent coupling of the perturbed equilibrium to NTV, and the PENTRC (Perturbed Equilibrium Neoclassical TRansport Code) formalism. diff --git a/docs/src/conventions.md b/docs/src/conventions.md index f599fef855..6a95d70114 100644 --- a/docs/src/conventions.md +++ b/docs/src/conventions.md @@ -1,20 +1,14 @@ # Conventions Reference -This page is a comprehensive reference for the sign, coordinate, and field-amplitude conventions -used throughout GPEC. Understanding them is essential for interpreting outputs, interfacing with -other codes, and correctly setting up perturbed-equilibrium and kinetic calculations. +The sign, coordinate and field-amplitude conventions used throughout GPEC, for reading its output and interfacing it with other codes. +They are inherited from the Fortran GPEC suite and predate the COCOS standard (see [COCOS compatibility](@ref)); each section names the source that establishes the convention. -These conventions are inherited from the Fortran GPEC suite and are GPEC-native: they predate -and do not fully coincide with the COCOS standard (see [COCOS Compatibility](@ref) below). -Each section names the Julia source that establishes the convention so the documentation stays -traceable to the implementation. - -## Quick Reference +## Quick reference | Quantity | Convention | Forced? | |----------|-----------|---------| | ``\psi`` (poloidal flux) | Normalized 0 (axis) to 1 (edge), no ``2\pi`` | Yes, positive | -| ``\theta`` (poloidal angle) | Upward outboard | Fixed | +| ``\theta`` (poloidal angle) | Upward outboard, one circuit = 1 | Fixed | | ``\phi`` (toroidal angle) | CCW for LH, CW for RH | By helicity | | ``F = R B_\phi`` | Always positive | Yes, `abs` | | ``q`` (safety factor) | Computed (direct) / input (inverse), positive | Effectively positive | @@ -22,8 +16,9 @@ traceable to the implementation. | Resonant ``m`` | Always positive (since ``n>0``, ``q>0``) | By construction | | helicity | +1 RH, ``-1`` LH | Computed | | ``\omega_E`` | Positive = direction of ``\zeta`` | No | +| Field amplitudes | ``b``, ``\bar b``, ``\tilde b``, all in tesla | By area weighting | -## Coordinate System +## Coordinate system GPEC uses right-handed magnetic coordinates ``(\psi, \theta, \zeta)`` with the Fourier kernel @@ -31,48 +26,28 @@ GPEC uses right-handed magnetic coordinates ``(\psi, \theta, \zeta)`` with the F \exp\!\big(i(m\theta - n\zeta)\big), ``` -matching the Fortran ``\exp(im\theta - in\phi)`` convention. The forward analysis transform that -produces the mode coefficients ``f_m`` from a real-space profile ``f(\theta)`` uses the conjugate -kernel ``\exp(-im\theta)``; the inverse transform reconstructs with ``\exp(+im\theta)`` -(`src/Utilities/FourierTransforms.jl`). - -### Poloidal Flux ``\psi`` - -- Normalized from 0 (magnetic axis) to 1 (plasma boundary). -- ``\psi_0 = |\psi_\mathrm{bry} - \psi_\mathrm{axis}|`` is **forced positive** on read - (`read_eq_efit` in `src/Equilibrium/Sources/GFile.jl`). If the EQDSK has - ``\psi_\mathrm{bry} < \psi_\mathrm{axis}``, the 2D flux array is sign-flipped so that ``\psi`` - increases from axis to edge. -- The internal ``\psi`` carries **no** factor of ``2\pi`` (poloidal flux per radian). IMAS inputs, - which use the full poloidal flux, are divided by ``2\pi`` on read — see - [COCOS Compatibility](@ref). -- Occasionally ``\rho = \sqrt{\psi}`` is used as a radius-like variable. - -### Poloidal Angle ``\theta`` - -- "Upward outboard" convention, always: ``\theta = 0`` on the outboard midplane, increasing - upward. -- ``\theta`` is normalized to increase by 1 (not ``2\pi``) over one poloidal circuit. This is fixed - regardless of helicity or working-coordinate choice. - -### Toroidal Coordinate ``\zeta`` and ``\phi`` - -- The magnetic coordinate toroidal angle is ``\zeta = (\phi - \nu(\psi,\theta))/(2\pi)``, where ``\phi`` is GPEC's internal - toroidal angle (the physical angle is ``-\mathrm{helicity}\,\phi``, below) and ``\nu`` is a single-valued - straight-field-line offset that depends on the working coordinate. PEST coordinates have ``\nu = 0``. -- The physical toroidal angle is reconstructed as - ``\phi = -\,\mathrm{helicity}\,(2\pi\zeta + \nu)`` (`sample_boundary_grid` in - `src/ForcingTerms/CoilFourier.jl`). Thus ``\phi`` is effectively **counter-clockwise** (viewed - from above) for left-handed (LH) configurations and **clockwise** for right-handed (RH) - configurations. - -### Working Coordinate Options - -Controlled by `jac_type` in the `[Equilibrium]` section of `gpec.toml`. The Jacobian is -``J \propto B_p^{p_{bp}}\, B^{p_b}\, R^{-p_r}\, r_c^{-p_{rc}}``, where ``B_p`` is the poloidal -field magnitude, ``B`` is the total field magnitude, ``R`` is the major radius, and -``r_c = \sqrt{(R-R_0)^2 + (Z-Z_0)^2}`` is the minor radius (distance from the magnetic axis -``(R_0, Z_0)``): +matching the Fortran ``\exp(im\theta - in\phi)`` convention, with the angles in radians in the kernel (``2\pi`` times the unit-period ``\theta`` and ``\zeta`` below). +The analysis transform that produces the mode coefficients ``f_m`` of a profile ``f(\theta)`` uses the conjugate kernel ``\exp(-im\theta)``, and the synthesis reconstructs with ``\exp(+im\theta)`` (`Utilities.FourierTransforms`). + +### Poloidal flux ``\psi`` + +- Normalized from 0 (magnetic axis) to 1 (plasma boundary); ``\psi_0 = |\psi_\mathrm{bry} - \psi_\mathrm{axis}|`` is forced positive on read, and an EQDSK with ``\psi_\mathrm{bry} < \psi_\mathrm{axis}`` has its flux array sign-flipped so that ``\psi`` increases outward (`Equilibrium.Sources.GFile`). +- The internal ``\psi`` carries no factor of ``2\pi`` (poloidal flux per radian); IMAS inputs, which use the full flux, are divided by ``2\pi`` on read. +- ``\rho = \sqrt{\psi}`` serves as a radius-like variable where one is needed. + +### Poloidal angle ``\theta`` + +- Upward outboard, always: ``\theta = 0`` on the outboard midplane, increasing upward. +- Normalized to increase by 1, not ``2\pi``, over one poloidal circuit, whatever the helicity or working coordinate. + +### Toroidal coordinate ``\zeta`` and ``\phi`` + +- The magnetic toroidal angle is ``\zeta = (\phi - \nu(\psi,\theta))/(2\pi)``, with ``\phi`` the internal toroidal angle and ``\nu`` the single-valued straight-field-line offset of the working coordinate (``\nu = 0`` in PEST coordinates). +- The physical (laboratory) toroidal angle is ``\phi_\mathrm{lab} = -\,\mathrm{helicity}\,(2\pi\zeta + \nu)`` (`ForcingTerms`), a different angle from the internal ``\phi`` above, so ``\phi_\mathrm{lab}`` runs counter-clockwise viewed from above for left-handed configurations and clockwise for right-handed ones. + +### Working coordinate options + +`[Equilibrium] jac_type` selects the Jacobian ``J \propto B_p^{p_{bp}}\, B^{p_b}\, R^{-p_r}\, r_c^{-p_{rc}}``, with ``B_p`` the poloidal field, ``B`` the total field, ``R`` the major radius and ``r_c`` the distance from the magnetic axis: | Name | `power_bp` | `power_b` | `power_r` | `power_rc` | |------|-----------|-----------|-----------|------------| @@ -82,346 +57,161 @@ field magnitude, ``B`` is the total field magnitude, ``R`` is the major radius, | Equal-arc | 1 | 0 | 0 | 0 | | Park | 0 | 1 | 0 | 0 | -The powers are set automatically from `jac_type` (`EquilibriumConfig` in -`src/Equilibrium/Config.jl`) and applied in the field-line Jacobian -(`direct_fieldline_der!` in `src/Equilibrium/Axisymmetric/Direct.jl`). Setting `jac_type = "custom"` -lets you specify the four exponents manually via `jac_custom_power_bp`, `jac_custom_power_b`, -`jac_custom_power_r`, and `jac_custom_power_rc`. - -## Helicity and Handedness +`jac_type = "custom"` sets the four exponents by hand through `jac_custom_power_bp`, `jac_custom_power_b`, `jac_custom_power_r` and `jac_custom_power_rc`. -Helicity is derived from the equilibrium current and toroidal-field directions: +## Helicity and handedness -```julia -helicity = bt_sign * Int(sign(crnt)) # src/ForcingTerms/CoilFourier.jl -``` - -- `bt_sign` is the toroidal-field direction and `crnt` is the signed plasma current. -- **helicity = +1**: right-handed (RH) --- ``B_t`` and ``I_p`` in the same direction. -- **helicity = -1**: left-handed (LH) --- ``B_t`` and ``I_p`` opposed. -- Helicity sets the toroidal-grid direction ``\phi_j = -\mathrm{helicity}\times 2\pi j/n_\zeta`` and - the output sign flips (see [Spectrum Output Sign Conventions](@ref)). For DIII-D - (``B_t < 0``, ``I_p > 0`` ⇒ helicity = ``-1``) ``\phi`` increases with ``j``. +Helicity is the product of the toroidal-field and plasma-current directions, `helicity = bt_sign * sign(crnt)` (`ForcingTerms`): +1 is right-handed (``B_t`` and ``I_p`` parallel), ``-1`` left-handed (opposed). +It sets the direction of the laboratory toroidal grid, ``\phi_{\mathrm{lab},j} = -\mathrm{helicity}\times 2\pi j/n_\zeta``, and the output sign flips of [Spectrum output sign conventions](@ref); for DIII-D (``B_t < 0``, ``I_p > 0``, helicity ``-1``) ``\phi_\mathrm{lab}`` increases with ``j``. ## ``F = R B_\phi`` -The poloidal-current function ``F = R B_\phi`` from the Grad-Shafranov equation is **forced -positive** via `abs` (`read_eq_efit` in `src/Equilibrium/Sources/GFile.jl`): - -```julia -abs.(fpol_data) -``` - -The code always works with ``|F|``; the sign of ``B_t`` is carried separately (`fpol_sign`, and the -`bt_sign`/`crnt` helicity inputs). +The poloidal-current function is forced positive on read, ``F = |F|`` (`Equilibrium.Sources.GFile`); the sign of ``B_t`` is carried separately by `fpol_sign` and the helicity inputs. -## Safety Factor ``q`` +## Safety factor ``q`` -The direct solver **computes** ``q`` by field-line integration for direct (EFIT) equilibria and -overwrites the file profile (`direct_run` in `src/Equilibrium/Axisymmetric/Direct.jl`): +The direct solver computes ``q`` by field-line integration and overwrites the file profile (`Equilibrium.Axisymmetric.Direct`), ```math -q = \frac{F}{2\pi}\oint \frac{J}{B_p}\,dl_p . +q = \frac{F}{2\pi}\oint \frac{J}{B_p}\,dl_p , ``` -- Because the Jacobian ``J > 0``, the closed-loop integral is positive, and ``F = |F|`` is forced - positive, **``q`` is positive** in standard operation. The solver emits a warning if the on-axis - extrapolation yields ``q_0 \le 0``, treating it as a spline artifact rather than a physical - result. -- **Inverse** equilibria, for which only CHEASE and analytic options are currently supported, instead use the ``q`` profile supplied with the input - (`src/Equilibrium/Axisymmetric/Inverse.jl`). This follows the fortran GPEC CHEASE reading convention. +which is positive because ``J > 0`` and ``F = |F|``; a warning marks an on-axis extrapolation with ``q_0 \le 0`` as a spline artifact. +Inverse equilibria use the ``q`` profile supplied with the input (`Equilibrium.Axisymmetric.Inverse`), the convention of the Fortran GPEC CHEASE reader. -## Mode Numbers ``m`` and ``n`` +## Mode numbers ``m`` and ``n`` -### Toroidal Mode Number ``n`` +### Toroidal mode number ``n`` -- Set as a **positive integer** via `[ForceFreeStates] nn` in `gpec.toml`. ``n < 1`` is not supported - and is clamped to 1 with a warning (`src/GeneralizedPerturbedEquilibrium.jl`). Results for - multiple ``n`` can be superposed. +`[ForceFreeStates] nn` is a positive integer or a list of consecutive integers. +``n < 1`` is not supported: a run whose modes are all below 1 is an error, and a lowest mode below 1 is clamped to 1 with a warning (`ForceFreeStates.Run`). +Results for several ``n`` can be superposed. -### Poloidal Mode Range +### Poloidal mode range -The poloidal spectrum spans `mlow` to `mhigh` (`src/GeneralizedPerturbedEquilibrium.jl`): +The poloidal spectrum spans `mlow` to `mhigh` (`ForceFreeStates.ForceBalance`): ```math -m_\mathrm{low} = \min(\lfloor n\,q_\mathrm{min}\rfloor,\,0) - 4 - \Delta m_\mathrm{low}, \qquad -m_\mathrm{high} = \lfloor n\,q_\mathrm{max}\rfloor + 2\,\Delta m_\mathrm{high}. +m_\mathrm{low} = \min(\lfloor n_\mathrm{low}\,q_\mathrm{min}\rfloor,\,0) - 4 - \Delta m_\mathrm{low}, \qquad +m_\mathrm{high} = \lfloor n_\mathrm{high}\,q_\mathrm{max}\rfloor + 2\,\Delta m_\mathrm{high}. ``` -`delta_mlow` and `delta_mhigh` widen the range beyond the resonant modes. Note the asymmetric -factor of 2 on ``\Delta m_\mathrm{high}``: `delta_mhigh` is doubled before being applied (for -consistency with the Fortran DCON convention), so the user-supplied value adds twice as many -modes on the high side as `delta_mlow` does on the low side. - -### Why Positive ``m`` Is Always Resonant - -Resonant surfaces are found by `sing_find!` (`src/ForceFreeStates/Surfaces/Finding.jl`), which locates flux -surfaces where ``m = n\,q(\psi)``. Since ``n > 0`` by convention and ``q > 0`` for a standard -tokamak, the resonant ``m`` is always positive: +`delta_mlow` and `delta_mhigh` widen the range beyond the resonant modes; `delta_mhigh` is doubled before use, the Fortran DCON convention, so it adds twice as many modes on the high side as `delta_mlow` does on the low side. -```math -m_\mathrm{res} = n\,q > 0 . -``` +### Why positive ``m`` is always resonant -Negative-``m`` modes are always non-resonant. This is by design: the ``m = 2`` displacement shows -resonant behavior at ``q = 2``, etc. +Resonant surfaces are the flux surfaces where ``m = n\,q(\psi)`` (`ForceFreeStates.Surfaces`). +Since ``n > 0`` by convention and ``q > 0`` for a standard tokamak, the resonant ``m`` is always positive, ``m_\mathrm{res} = n\,q > 0``, and negative-``m`` modes are never resonant: the ``m = 2`` displacement resonates at ``q = 2``. -## Spectrum Output Sign Conventions +## Spectrum output sign conventions -For the real-space representation, GPEC takes the complex conjugate for **right-handed** -configurations. This is implemented in the mode reconstruction -(`src/Analysis/PerturbedEquilibriumModes.jl`): - -```julia -if helicity > 0 - theta_data[:, :, k] .= conj.(theta_data[:, :, k]) -end -``` - -The straight-field-line offset enters as ``\exp(i\,n\,\nu)`` before the conjugation. For the full -Fourier representation ``\exp(i(m\theta - n\zeta))`` the conjugate operation flips both the -toroidal direction and the up/down (poloidal) sense. +For the real-space representation GPEC takes the complex conjugate of the mode amplitudes for right-handed configurations (`Analysis.PerturbedEquilibriumModes`), after the straight-field-line offset has entered as ``\exp(i\,n\,\nu)``. +For the full Fourier representation ``\exp(i(m\theta - n\zeta))`` the conjugation flips both the toroidal direction and the up/down sense. ### Interfacing with SURFMN -SURFMN expands in ``\exp(-im\theta - in\phi)`` with CCW ``\phi``. To convert from GPEC: +SURFMN expands in ``\exp(-im\theta - in\phi)`` with CCW ``\phi``; to convert from GPEC, ```python m_surfmn = helicity * m_gpec b_surfmn = real(b_m) - 1j * helicity * imag(b_m) ``` -For LH configurations only the sign of ``m`` is flipped; for RH, ``m`` is unchanged but the complex -conjugate is taken. +so that for LH configurations only the sign of ``m`` is flipped and for RH the complex conjugate is taken. ### Interfacing with Vacuum -Unlike the Fortran VACUUM code, the GPEC Vacuum module works in the same angles as the rest of the code: ``\theta`` upward-outboard (counter-clockwise in the -``(R, Z)`` plane) and the Fourier basis ``\exp(-i(m\theta - 2\pi n\zeta))``. A plasma contour passed in clockwise is rejected; -a clockwise wall file is reversed with a warning. The ``\nu`` passed to Vacuum is GPEC's ``\nu = \phi - 2\pi\zeta``. -Chance 1997 eq. 97 writes ``\zeta = \phi + \nu`` for the same ``\nu`` because the Fortran VACUUM code's ``\theta`` and ``\phi`` both -run opposite to GPEC's. +The `Vacuum` module works in the same angles as the rest of GPEC, ``\theta`` upward-outboard and the analysis basis ``\exp(-i(m\theta - n\zeta))`` in radians, unlike the Fortran VACUUM code whose ``\theta`` and ``\phi`` both run the other way; a clockwise plasma contour is rejected and a clockwise wall file is reversed with a warning. +The ``\nu`` passed to it is GPEC's ``\nu = \phi - 2\pi\zeta``, which Chance (1997, Eq. 97) writes as ``\zeta = \phi + \nu``. -## Field Amplitudes and Units +## Field amplitudes and units -GPEC reports perturbed magnetic fields as Fourier spectra on flux surfaces. The spectrum of a field -on a surface depends in general on the coordinate system used. However, certain resonant components -on rational surfaces or the 2-norm of the Fourier coefficient vector on any surface can be made -coordinate-independent if the proper area weighting is used within the Fourier decomposition. The -coordinate-invariance of the weightings used in GPEC is established in Pharr (2026) (see -[Citations](citations.md)). In each case, we normalize the harmonics by scalar area factors in such -a way that the units are always Tesla. +GPEC reports perturbed fields as Fourier spectra on flux surfaces. +A raw spectrum depends on the working coordinate ([Park, Boozer and Menard (2008)](citations.md#PerturbedEquilibrium)), but the resonant component on a rational surface and the 2-norm of the spectrum on any surface become coordinate-invariant with the right area weighting in the Fourier integrand ([Pharr et al. (2026)](citations.md#PerturbedEquilibrium)). +GPEC therefore works with three decompositions of the normal field, all in tesla, distinguished only by that weight: -### Resonant Field - -The resonant field is the pitch-resonant (``m = nq``) Fourier component of the perturbed flux at a -rational surface, divided by the **scalar area** ``A^r`` of that surface: +| Pharr symbol | Integrand weight | Name | HDF5 token | +|---|---|---|---| +| ``b`` (bare) | ``1`` | normal field | `b_n` | +| ``\bar b`` (bar) | ``\mathcal{J}\,\lvert\nabla\psi\rvert`` | area-weighted field | `area` | +| ``\tilde b`` (tilde) | ``\sqrt{\mathcal{J}\,\lvert\nabla\psi\rvert}`` | root-area-weighted field | `root_area` | ```math -b^r = \frac{\Phi^r}{A^r} \qquad [\mathrm{T}]. +b_m = \oint (\mathbf{b}\cdot\hat{\mathbf n})\, e^{-i m\theta}\, d\theta, \qquad +\bar b_m = \frac{1}{A}\oint \mathcal{J}\lvert\nabla\psi\rvert\,(\mathbf{b}\cdot\hat{\mathbf n})\, e^{-i m\theta}\, d\theta, \qquad +\tilde b_m = \frac{1}{\sqrt{A}}\oint \sqrt{\mathcal{J}\lvert\nabla\psi\rvert}\,(\mathbf{b}\cdot\hat{\mathbf n})\, e^{-i m\theta}\, d\theta , ``` -Dividing the resonant flux ``\Phi^r`` by the surface area turns it into a field amplitude in Tesla -that is invariant under changes of the poloidal-angle (working) coordinate. This is the quantity -reported as `resonant_flux`. The need for this normalization is the central point of Park, Boozer, -and Menard (2008) (see [Citations](citations.md)): the raw Fourier spectrum of the perturbed field is -*spectrally asymmetric* and changes with the magnetic working coordinate, so only the resonant -harmonic scaled by the scalar surface area ``A^r`` is a coordinate-independent physical measure of -the resonant field. - -Note that there is strictly NO resonant field in ideal MHD perturbed equilibrium calculations, and -only the effective resonant field is non-zero in these cases. Effective resonant fields are computed -as above but for the resonant component of flux that is being shielded by the shielding current layer -on the rational surface. Resistive and kinetic MHD perturbed equilibria have, in general, finite -resonant fields as well as effective resonant fields. +with ``A`` the scalar surface area. +The root-area weighting is the unique one under which the spectrum transforms unitarily between working coordinates, so ``\sum_m |\tilde b_m|^2 = \langle (\mathbf{b}\cdot\hat{\mathbf n})^2 \rangle`` (Parseval) and the quantities derived from ``\tilde b``, such as the singular values of the resonant coupling matrix, are coordinate-invariant; "power-normalized field" in early GPEC literature names this property. +The area weighting is coordinate-invariant for the pitch-resonant ``m = nq`` harmonic on a rational surface: the **resonant field** ``\bar b^{\,r} = \Phi^r/A^r``, the resonant flux divided by the surface area, stored as `PerturbedEquilibrium/SingularCoupling/resonant_area_weighted_field`. +In an ideal perturbed equilibrium the surface is perfectly shielded and the stored value is the effective resonant field, the resonant flux the shielding current cancels; resistive and kinetic equilibria have finite resonant fields as well. -### Root-Area Normalized Field - -The root-area normalized field is the Fourier spectrum of the square-root-area-weighted normal field -``\sqrt{\mathcal{J}\,|\nabla\psi|}\,(\mathbf{b}\cdot\hat{\mathbf{n}})``, divided by the scalar -``\sqrt{A}`` of the surface: - -```math -\tilde{b}_m = \frac{1}{\sqrt{A}}\oint \sqrt{\mathcal{J}\,|\nabla\psi|}\; -(\mathbf{b}\cdot\hat{\mathbf{n}})\, e^{-i m\theta}\, d\theta \qquad [\mathrm{T}]. -``` - -According to Parseval's theorem, the sum of its squared mode amplitudes equals the -area-averaged squared normal field — the surface-averaged field "power", +With ``\Sigma`` the mode-space square-root-weight convolution (`sqrtamat`) and ``A`` the surface area (`jarea`), the three fields and the flux are related by (`Equilibrium.CoordinateInvariant`) ```math -\sum_m |\tilde{b}_m|^2 = \big\langle (\mathbf{b}\cdot\hat{\mathbf{n}})^2 \big\rangle. +\tilde b = \Sigma\, b, \qquad +\bar b = (\Sigma/\sqrt{A})\,\tilde b, \qquad +\Phi = A\,\bar b = \Sigma\sqrt{A}\;\tilde b , ``` -The ``\sqrt{A}`` weighting is the unique one for which the spectrum transforms unitarily between -working coordinates, so the power-normalized amplitudes — and quantities derived from them, such as -the singular values of the resonant coupling matrix — are coordinate-invariant. - -### The Three Field Amplitudes +where `rootarea_to_area_weight` ``= \Sigma/\sqrt{A}`` maps ``\tilde b \to \bar b``, `area_to_rootarea_weight` is its inverse, and ``R = \Sigma\sqrt{A}`` is the internal flux-conform operator. +Flux ``\Phi`` is never stored: it appears only when a user supplies flux-valued forcing and is recovered on demand as ``\Phi = A\,\bar b``. -Following Pharr (2026), *Coordinate-invariant flux-surface Fourier analysis in tokamaks* (see -[Citations](citations.md)), GPEC works with **three** Fourier decompositions of the normal field, **all -in field units (tesla)**, distinguished only by the **area weight applied to the Fourier integrand**. -These are the authoritative names used in prose, in code, and in the HDF5 output: +### Field representations in GPEC output -| Pharr symbol | FT integrand weight ``W`` | English name | HDF5 token | -|---|---|---|---| -| ``b`` (bare) | ``1`` | **normal field** | `b_n` | -| ``\bar b`` (bar) | ``\mathcal{J}\,\lvert\nabla\psi\rvert`` | **area-weighted field** | `area` | -| ``\tilde b`` (tilde) | ``\sqrt{\mathcal{J}\,\lvert\nabla\psi\rvert}`` | **root-area-weighted field** | `root_area` | - -All three carry units of tesla because each integral is divided by the appropriate power of the -coordinate-invariant scalar surface area ``A``: - -```math -b_m = \oint (\mathbf{b}\cdot\hat{\mathbf n})\, e^{-i m\theta}\, d\theta, \qquad -\bar b_m = \frac{1}{A}\oint \mathcal{J}\lvert\nabla\psi\rvert\,(\mathbf{b}\cdot\hat{\mathbf n})\, e^{-i m\theta}\, d\theta, \qquad -\tilde b_m = \frac{1}{\sqrt{A}}\oint \sqrt{\mathcal{J}\lvert\nabla\psi\rvert}\,(\mathbf{b}\cdot\hat{\mathbf n})\, e^{-i m\theta}\, d\theta . -``` +- `PerturbedEquilibrium/forcing_b`, `forcing_b_root_area`, `forcing_b_area` and the `response_*` triplet: the control-surface spectra in the bare, root-area-weighted and area-weighted representations. +- `PerturbedEquilibrium/Response/b_n` and `b_psi_area_weighted`: the bare normal field and the area-weighted radial field profiles. +- `PerturbedEquilibrium/SingularCoupling/resonant_area_weighted_field`, `C_resonant_area_weighted_field` and `penetrated_area_weighted_field`: the resonant area-weighted field, its coupling matrix and the penetrated field. +- `PerturbedEquilibrium/ResponseMatrices/`: the response matrices in root-area-weighted space, with `rootarea_to_area_weight_operator` ``S`` (``L_{\bar b} = S\,\tilde L\,S^\dagger``) and `surface_area` ``A`` (``\Phi = A\,\bar b``) to recover the other forms. +- `ForceFreeStates/FreeBoundaryStability/eigenmode_energies`: energies per unit surface-averaged squared displacement, a normalization described on the [Stability](stability.md) page. -Only the **square-root-area weighting** ``\tilde b`` transforms unitarily between working coordinates, so -its 2-norm (and quantities derived from it, such as the singular values of the resonant coupling matrix) -are coordinate-invariant. Note that the "root-area-weighted field" sometimes appears as "power-normalized -field" in early GPEC literature, as the latter names a Parseval *property* of ``\tilde b``. The **full-area weighting** ``\bar b`` is coordinate-invariant for the -pitch-resonant ``m = nq`` harmonic on a rational surface — the **resonant area-weighted field** -``\bar b^{\,r} = \Phi^r/A^r`` as described above. +## KineticForces rotation conventions -GPEC operates and outputs **only in these field representations — poloidal flux ``\Phi`` (weber) is never -stored.** Flux appears, briefly, only internally when a user supplies flux-valued forcing, and is -recovered on demand as the scalar product ``\Phi = A\,\bar b``. +The rotation profiles enter the `KineticForces` integrals as follows. -### Translation Operators +### ``\omega_E`` (E × B rotation) -With ``\Sigma`` ≡ `sqrtamat` (the mode-space ``\sqrt{}``weight convolution) and the scalar surface area -``A`` ≡ `jarea`, the three fields are related by (`Equilibrium/CoordinateInvariant.jl`): +``\omega_E`` is the input E × B toroidal rotation profile, positive in the direction of the toroidal coordinate ``\zeta``. +In the laboratory frame, with ``\phi`` counter-clockwise viewed from above, ``\phi = -h\,(2\pi\zeta + \nu)`` for helicity ``h`` ([Helicity and handedness](@ref)), so positive ``\omega_E`` is counter-clockwise with sign ``-h``, and its sign relative to the plasma current is ``-h\,s_I = -s_B`` with ``s_I``, ``s_B`` the signs of ``I_p`` and ``B_t`` (`ip_sign`, `bt_sign`): positive ``\omega_E`` is co-current exactly when ``B_t`` is clockwise (``s_B = -1``), whatever the helicity. -```math -\tilde b = \Sigma\, b, \qquad -\bar b = (\Sigma/\sqrt{A})\,\tilde b, \qquad -\Phi = A\,\bar b = \Sigma\sqrt{A}\;\tilde b . -``` +### Diamagnetic frequencies -- **`rootarea_to_area_weight`** ``= \Sigma/\sqrt{A}`` maps ``\tilde b \to \bar b``. -- **`area_to_rootarea_weight`** ``= \sqrt{A}\,\Sigma^{-1}`` is its inverse (``\bar b \to \tilde b``). - -The internal flux-conform operator is just ``R = \Sigma\sqrt{A} = `` `rootarea_to_area_weight` ``\cdot A``. - -### Field Representations in GPEC Output - -- **`forcing_b` / `forcing_b_root_area` / `forcing_b_area`** (and the `response_*` triplet) — the - control-surface forcing and response spectra in the bare (``b``), root-area-weighted (``\tilde b``) and - area-weighted (``\bar b``) representations, under `PerturbedEquilibrium/`. -- **`b_n`** — the bare normal field ``\mathbf{b}\cdot\hat{\mathbf n}`` (and the area-weighted radial field - `b_psi_area_weighted`), under `PerturbedEquilibrium/Response/`. -- **`resonant_area_weighted_field`** / **`C_resonant_area_weighted_field`** — the resonant area-weighted - field ``\bar b^{\,r} = \Phi^r/A^r`` and its coupling matrix, under - `PerturbedEquilibrium/SingularCoupling/`. The sibling `penetrated_area_weighted_field` follows the - same convention. -- **Root-area-weighted (``\tilde b``) space** — the control-surface response matrices (`permeability`, - `reluctance`, `plasma_inductance`, `surface_inductance`) are stored in this coordinate-invariant space - under `PerturbedEquilibrium/ResponseMatrices/`. The stored `rootarea_to_area_weight_operator` ``S`` - recovers the area-weighted field forms (``L_{\bar b} = S\,\tilde L\,S^\dagger``) and the scalar - `surface_area` ``A`` recovers flux (``\Phi = A\,\bar b``). -- **Eigenmode energies** (`ForceFreeStates/FreeBoundaryStability/eigenmode_energies`) — the generalized eigenvalues of the - pencil ``W\,v = \lambda\,N\,v``, where ``N`` is the power-normalization (surface-norm) matrix built from - the Jacobian Fourier coefficients: ``\xi^\dagger N\,\xi = \oint J\,|\xi(\theta)|^2\,d\theta\,/\,V' = - \langle|\xi|^2\rangle``. The eigenvalues are stationary values of the power quotient - ``\xi^\dagger W \xi / \xi^\dagger N \xi`` — energies per unit surface-averaged squared displacement - (a normalization convention, not raw joules). Because ``W`` and ``N`` transform by the same congruence - under a change of working (Jacobian) coordinate, the spectrum is invariant to that coordinate choice. - ``N`` is not an inertia matrix (the stability stage carries no velocities); it is the DCON power - normalization promoted from a post-hoc per-eigenvector division into the metric of the eigenproblem. - - -## Rotation Velocity Conventions (KineticForces) - -The `KineticForces` module (the Julia port of PENTRC) conventions below are used in -`src/KineticForces/Torque.jl`, `EnergyIntegration.jl`, and `PitchIntegration.jl`. - -### ``\omega_E`` (E × B Rotation) - -- ``\omega_E`` is the input E × B toroidal-rotation profile, evaluated as - `welec = omegaE_spline(psi)` (`Torque.jl`). Positive ``\omega_E`` means rotation in the direction - of the toroidal coordinate ``\zeta``. -- In the laboratory frame, where ``\phi`` increases counter-clockwise viewed from above, the toroidal - coordinate maps as ``\phi = -h\,(2\pi\zeta + \nu)`` with ``h`` the helicity - ([Helicity and Handedness](@ref)), so positive ``\omega_E`` has counter-clockwise sign ``-h``. - Relative to the plasma current its sign is ``-h\,s_I = -s_B``, where ``s_I`` and ``s_B`` are the - signs of ``I_p`` and ``B_t`` (`ip_sign`, `bt_sign`): positive ``\omega_E`` is co-current exactly when - ``B_t`` is clockwise (``s_B = -1``), whatever the helicity. For ``s_I = s_B = +1`` it is clockwise and - counter-current. - -### Diamagnetic Frequencies - -Computed at evaluation time from the kinetic-profile gradients (Logan & Park 2013, Eq. 7; -`Torque.jl`): +Computed at evaluation time from the kinetic-profile gradients (Logan et al. (2013), Eq. 7): ```math \omega_{*n} = -\frac{2\pi\, T_i}{e\, Z_i\, \chi_1\, n_i}\frac{dn_i}{d\psi_n}, \qquad \omega_{*T} = -\frac{2\pi}{e\, Z_i\, \chi_1}\frac{dT_i}{d\psi_n}, ``` -where ``\chi_1 = 2\pi\psi_0`` (`chi1 = 2π·psio`). The negative signs mean a negative -(outward-decreasing) density or temperature gradient yields a positive diamagnetic frequency. The -total toroidal rotation is the sum +with ``\chi_1 = 2\pi\psi_0``, so that an outward-decreasing density or temperature gives a positive diamagnetic frequency; the total toroidal rotation is ``\omega_\phi = \omega_E + \omega_{*n} + \omega_{*T}``. -```math -\omega_\phi = \omega_E + \omega_{*n} + \omega_{*T} . -``` +### Energy integral resonance -### Energy Integral Resonance - -In the kinetic energy integral (`integrate_energy` in `EnergyIntegration.jl`) the resonance -denominator is ``i\,\Omega(x) - \nu``, where the resonance condition is +The resonance denominator of the energy integral is ``i\,\Omega(x) - \nu`` with ```math \Omega(x) = n\,\omega_E + \ell_\mathrm{eff}\,\omega_b\sqrt{x} + n\,\omega_D\,x , +\qquad +\ell_\mathrm{eff} = \begin{cases} \ell + n q & \text{passing particles},\\ \ell & \text{trapped particles}, \end{cases} ``` -with ``x = E/T`` the normalized energy and ``\nu`` the (energy-dependent) effective collisionality. -Here: - -- ``\omega_b`` and ``\omega_D`` are the bounce and magnetic-precession frequency **coefficients at - thermal energy** (``x = 1``): ``\omega_b \propto v_\mathrm{th}`` and - ``\omega_D \propto v_\mathrm{th}^2`` (`wbbar`, `wdbar` times `bhat`, `dhat` in `BounceAveraging.jl`). Their velocity-space - dependence enters explicitly as ``\sqrt{x}`` (bounce ``\propto v``) and ``x`` (precession - ``\propto E``). -- ``\ell_\mathrm{eff}`` is the effective bounce harmonic (`PitchIntegration.jl`): - -```math -\ell_\mathrm{eff} = -\begin{cases} -\ell + n q & \text{circulating (passing) particles},\\ -\ell & \text{trapped particles}. -\end{cases} -``` - -The sign of ``\omega_E`` shifts the resonance location in velocity space. - - -## COCOS Compatibility - -The conventions above are GPEC-native and predate the COCOS standard, as defined in -[Sauter, Comp. Phys. Comm. 2013](https://doi.org/10.1016/j.cpc.2012.09.010). The internal convention: +where ``x = E/T`` is the normalized energy, ``\nu`` the energy-dependent effective collisionality, and ``\omega_b``, ``\omega_D`` the bounce and precession frequencies at thermal energy (``\omega_b \propto v_\mathrm{th}``, ``\omega_D \propto v_\mathrm{th}^2``; `BounceAveraging.jl`). +The sign of ``\omega_E`` shifts the resonance in velocity space. -- Right-handed magnetic coordinates ``(\psi, \theta, \zeta)``. -- ``\psi`` excludes the ``2\pi`` factor (poloidal flux per radian) and increases from axis to edge. -- ``q > 0`` and ``n > 0`` in standard operation; ``F = |F|``. +## COCOS compatibility -The code labels this internal convention **COCOS 2** (`read_imas` in -`src/Equilibrium/Sources/IMAS.jl`) and converts IMAS inputs accordingly: +The conventions above are GPEC-native and predate the COCOS standard of [Sauter and Medvedev (2013)](https://doi.org/10.1016/j.cpc.2012.09.010): right-handed ``(\psi, \theta, \zeta)``, ``\psi`` without the ``2\pi`` and increasing outward, ``q > 0``, ``n > 0`` and ``F = |F|``. +The code labels this internal convention COCOS 2 (`Equilibrium.Sources.IMAS`) and accepts IMAS input in two forms: ```julia imas_cocos = 11 # default: IMAS standard, ψ divided by 2π on read imas_cocos = 2 # data already in the internal convention, no conversion ``` -!!! warning "The COCOS indexing is still a beta feature" - Note that the IMAS reader converts COCOS 11 → "COCOS 2" by simply dividing ``\psi`` by ``2\pi`` - (i.e. conversion from COCOS 11 to 1) and relying on the standardly enforced sign conventions - above to eventually conform to the COCOS 2 conventions. +The conversion from COCOS 11 divides ``\psi`` by ``2\pi`` and relies on the forced sign conventions above (positive ``\psi``, ``F`` and ``q``) for the rest; no other COCOS index is converted. ## See also -- [Workflow](workflow.md) --- where each convention enters the analysis pipeline -- [Citations](citations.md) --- theoretical references for GPEC's algorithms +- [Workflow](workflow.md), where each convention enters the pipeline +- [Citations](citations.md), the references behind the conventions diff --git a/docs/src/developer_notes.md b/docs/src/developer_notes.md index f25c828e1a..e84e98d68d 100644 --- a/docs/src/developer_notes.md +++ b/docs/src/developer_notes.md @@ -1,6 +1,6 @@ # Developer Notes -## Git Workflow +## Git workflow This project uses [GitFlow](http://nvie.com/posts/a-successful-git-branching-model): @@ -8,7 +8,7 @@ This project uses [GitFlow](http://nvie.com/posts/a-successful-git-branching-mod - `main` is updated only at release-ready stages - Feature branches should come off `develop` and merge back with `--no-ff` -## Commit Message Format +## Commit message format ``` Area[.Submodule] - TAG[!] - Imperative summary @@ -21,15 +21,15 @@ Both vocabularies are closed — do not invent new ones. Pull request titles use ## Documentation standard -Module reference pages follow a shared *Journal of Computational Physics*-style -structure (governing equations → numerical method → validation figures → API), -and every documentation figure is committed together with the script that made +Module reference pages follow one skeleton (*What it computes* → *Governing equations* → +*Running it*, with a deck snippet from a shipped example and the matching script call → *Outputs*, +a table of `gpec.h5` datasets with units → *Numerical notes*, including validation figures → +*API reference*), and every documentation figure is committed together with the script that made it and a provenance stamp. The full policy — page template, figure organization under `docs/src/figures//`, provenance, and the regenerate-only-when-`depends`-change rule — is in [`docs/DOC_STANDARD.md`](https://github.com/OpenFUSIONToolkit/GPEC/blob/develop/docs/DOC_STANDARD.md). -The `Ballooning and Mercier Local Stability` and `Inner Layer Module` pages are -the reference exemplars. +The `Stability` and `Local Stability` pages are the reference exemplars. ## Deck and API parity @@ -41,103 +41,16 @@ from the same inputs through the API, and checks that the two results agree. A k between the paths is marked `@test_broken`, not dropped, so the change that fixes it shows up as an unexpected pass. -## Regression Testing +## Regression testing The regression harness **must be run on every pull request before merging into `develop`**. It is the project's primary safeguard for tracking how numerical results evolve across changes, so it is only useful if every PR exercises it. When you open a PR, paste the regression report into the PR thread so reviewers can see what moved (and what did not). If your change touches a quantity that is not yet tracked, add a new regression case — or extend an existing one — in the same PR. -### Open problem: pinning grid-sensitive Δ′ robustly - -The ideal-MHD Δ′ values pinned in `test/runtests_parallel_integration.jl` -(`delta_prime_matrix` diagonal) and tracked by the harness are **single-point -snapshots**: one value at one `psi_accuracy` on one grid. They exist to catch -unintended changes, and are explicitly *not* converged Δ′. The extraction is -intrinsically grid-sensitive — a `psi_accuracy` scan from 2e-3 to 2.5e-4 swings -the DIII-D-like `dpm[1,1]` by roughly 50% (about 6.2 to 9.9), and switching the -grid generator moves it again. Any such pin therefore encodes an arbitrary point -on a varying curve, and re-pinning is required whenever the grid changes, which -weakens it as a regression signal. - -A more defensible criterion is to pin the **plateau**: scan Δ′ across -`psi_accuracy` and the integration-truncation controls, and take the value where -the result is stationary — the mode of the resulting distribution rather than any -single sample. Where a plateau exists it is a property of the physics rather than -of the discretization, so it would survive grid changes and would be a genuine -convergence statement. - -A `psi_accuracy` scan on the DIII-D-like SLAYER deck (2e-3 down to 3.125e-5, a -64x range) shows what a plateau search is up against: - -| psi_accuracy | knots | implied | dpm[1,1] | dpm[2,2] | dpm[3,3] | gamma 2/1 (Hz) | -|---|---|---|---|---|---|---| -| 2.0e-3 | 172 | – | 6.850 | -5.783 | -16.308 | 165.3 | -| 1.0e-3 | 205 | – | 6.387 | -5.873 | -16.867 | 154.3 | -| 5.0e-4 | 252 | – | 8.003 | -6.105 | -16.413 | 192.9 | -| 2.5e-4 | 307 | 522 | 8.567 | -6.235 | -16.529 | 206.2 | -| 1.25e-4 | 372 | 666 | 8.146 | -6.285 | -16.649 | 196.2 | -| 6.25e-5 | 457 | 916 | 8.911 | -6.477 | -16.510 | 214.4 | -| 3.125e-5 | 559 | 1226 | 9.514 | -6.260 | -16.485 | 228.7 | - -Three things follow. First, convergence is per-surface: the outer surface -`dpm[3,3]` (q=4) *does* plateau under the auto grid — the last four points agree -to about 1% — while `dpm[2,2]` is marginal and `dpm[1,1]` (q=2) never settles, -still moving ~7% between the two tightest grids. A plateau detector must -therefore report per-surface rather than pass/fail for the whole diagonal. -(As the `core_edge_packed` scan below shows, q=2 is not inherently unconvergeable — it is the -auto grid that prevents it from settling.) - -Second, the growth rate is linear in Δ′: `gamma/dpm[1,1]` is 24.1 to within 0.5% -across the whole scan, so a converged Δ′ is both necessary and sufficient for a -converged SLAYER growth rate on this surface. - -Third — and this is the blocker — the `implied` column (what -`implied_knot_count` requests after the refined solve) grows monotonically -relative to the grid actually used, from 1.7x to 2.2x. The two-pass scheme stops -after pass 2, so tightening `psi_accuracy` moves it *further* from -self-consistency rather than closer, and the warning's advice to "consider -tightening psi_accuracy" is counterproductive in this regime. - -The same deck on a deterministic `core_edge_packed` grid, which skips the measure-and-re-form -step entirely, converges: - -| mpsi | dpm[1,1] | dpm[2,2] | dpm[3,3] | gamma 2/1 (Hz) | -|---|---|---|---|---| -| 128 | 5.208 | -8.666 | -19.157 | 126.2 | -| 256 | 7.097 | -5.187 | -15.729 | 171.2 | -| 512 | 8.489 | -6.217 | -16.517 | 204.4 | -| 1024 | 8.932 | -6.298 | -16.529 | 214.9 | -| 2048 | 8.914 | -6.385 | -16.468 | 214.5 | - -Successive `dpm[1,1]` increments are +1.889, +1.392, +0.443, -0.018: the last -doubling moves both Δ′ and the growth rate by 0.2%. So the q=2 Δ′ is **not** -intrinsically ill-conditioned — it converges to ≈8.91, with gamma ≈214.5 Hz, and -`dpm[3,3]` converges to ≈-16.5 on *both* grid families, which is what a physical -value should do. The non-convergence under the auto grid is an artifact of the -generator, not of the Δ′ extraction. - -Two consequences. A plateau criterion is implementable today against a fixed -`core_edge_packed` grid, without waiting on the auto-grid work. And the auto grid's answers -are biased in both directions relative to the converged value: at its default -`psi_accuracy` it gave 6.39 (28% low), at its tightest 9.51 (7% high). Anything -pinned on the auto grid should be read with that in mind. - -This is not implemented. Doing it properly needs: - - - either a fixed `core_edge_packed` grid (which already converges, see above) or knot - refinement iterated to a fixed point (repeat the measure-and-re-form - step until `implied_knot_count` stops exceeding the grid in use), since - without it the scan target keeps moving; - - a scan driver over `psi_accuracy` × truncation (`psihigh`, `rational_offset_fraction`, `qhigh`) - that records the Δ′ diagonal per configuration; - - a plateau/mode detector with an explicit stationarity tolerance, plus a - defined failure mode for surfaces where no plateau exists (the - near-separatrix surfaces are expected to fall in this class); - - a harness case that pins plateau values and their scan width, replacing the - single-point pins. - -Until then, treat the pinned diagonal as an order-of-magnitude and sign -diagnostic only, and expect to re-pin it whenever the equilibrium grid changes. - -### Regression Harness: Quick usage guide +### Pinned Δ′ values + +The Δ′ values pinned in the tests and tracked by the harness are single-grid snapshots, not converged values: the extraction is grid-sensitive, and on the DIII-D-like deck the q = 2 diagonal converges only on a fixed `core_edge_packed` grid, while the automatic grid leaves it moving by several percent between refinements. +Treat a pinned diagonal as a sign and order-of-magnitude regression signal, expect to re-pin it when the equilibrium grid changes, and see `docs/development/regression-harness.md` for the convergence scans behind this. + +### Regression harness: quick usage guide Set up an alias for convenience (optional): ```bash @@ -148,64 +61,11 @@ alias regress='julia --project=regression-harness regression-harness/regress.jl' ```bash regress --list-cases ``` -``` -Available regression cases: ----------------------------------------------------------------- - diiid_n1 DIII-D-like equilibrium, n=1, ideal + perturbed equilibrium - dir: examples/DIIID-like_ideal_example (24 quantities) - diiid_multi_n DIII-D-like equilibrium, n = 1, 2, ideal + perturbed equilibrium - dir: examples/DIIID-like_ideal_multi_n_example (21 quantities) - solovev_multi_n Solovev analytical equilibrium, multi-n, ideal stability - dir: examples/Solovev_ideal_example_multi_n (12 quantities) - solovev_n1 Solovev analytical equilibrium, n=1, ideal stability - dir: examples/Solovev_ideal_example (18 quantities) -``` - **Compare two branches/commits:** ```bash regress --cases diiid_n1 --refs develop,feature/kinetic-damping ``` -``` -================================================================ -Case: diiid_n1 — DIII-D-like equilibrium, n=1, ideal + perturbed equilibrium -================================================================ -[ Info: Cached: diiid_n1 @ 0a905a7d (2026-04-06T23:41:50+09:00) -[ Info: Cached: diiid_n1 @ 44b2494f (2026-04-08T18:30:46+09:00) - -Regression Report: diiid_n1 -================================================================================================================== -Ref 1: develop @ 0a905a7d (2026-04-06) -Ref 2: feature/kinetic-damping @ 44b2494f (2026-04-08) ------------------------------------------------------------------------------------------------------------------- -Quantity develop feature/kinetic-damping Diff Status ------------------------------------------------------------------------------------------------------------------- -beta_n -1.376214e+00 -1.376214e+00 0.0e+00 OK -beta_t 1.322850e-02 1.322850e-02 0.0e+00 OK -Chirikov parameter [4 elements] [4 elements] 0.0e+00 OK -delta prime [4 elements] [4 elements] 0.0e+00 OK -plasma energy Re(ep[1]) -8.809610e-01 -8.809610e-01 0.0e+00 OK -total energy Im(et[1]) 6.175834e-05 6.175834e-05 0.0e+00 OK -total energy Re(et[1]) 1.199597e+00 1.199597e+00 0.0e+00 OK -vacuum energy Re(ev[1]) 2.080558e+00 2.080558e+00 0.0e+00 OK -island half-widths [4 elements] [4 elements] 0.0e+00 OK -mpert 34 34 0.0e+00 OK -# singular surfaces 4 4 0.0e+00 OK -npert 1 1 0.0e+00 OK -ODE steps (saved) 740 740 0.0e+00 OK -ODE steps (total) 1348 1348 0.0e+00 OK -PE plasma energy 0.000000e+00 0.000000e+00 0.0e+00 OK -PE total energy 0.000000e+00 0.000000e+00 0.0e+00 OK -pressure profile (checksum) 657ad2329d7b... 657ad2329d7b... identical OK -q0 1.209710e+00 1.209710e+00 0.0e+00 OK -q95 4.505007e+00 4.505007e+00 0.0e+00 OK -q profile (checksum) 75912afcc351... 75912afcc351... identical OK -||resonant flux|| 4.523707e+02 4.523707e+02 0.0e+00 OK -Runtime (s) 50.9s 52.0s -- -singular psi locations [4 elements] [4 elements] 0.0e+00 OK -singular q values [4 elements] [4 elements] 0.0e+00 OK -================================================================================================================== -Summary: 23 unchanged, 3 missing/N/A -``` +The report lists every tracked quantity of the case for each ref with its difference and an `OK`/changed status. **Compare your uncommitted working tree against develop:** ```bash @@ -216,16 +76,6 @@ regress --cases solovev_n1 --refs develop,local ```bash regress --show et_real --case solovev_n1 ``` -``` -History: et_real — solovev_n1 -================================================================================ -Commit Date Value Δ from prev Status --------------------------------------------------------------------------------- -edff6e86 2026-04-02 -4.624928e-01 -- -- -0a905a7d 2026-04-06 -4.624928e-01 0.0e+00 OK -================================================================================ -``` - **Scan across a range of commits (git-bisect style):** ```bash regress --cases solovev_n1 --ref-range develop~10..develop @@ -239,3 +89,103 @@ regress --cases solovev_n1 --ref-range develop~10..develop The full flag list, the environment-pinning behaviour that keeps source code the only variable in a comparison, and the exit-status contract are documented in [`docs/development/regression-harness.md`](https://github.com/OpenFUSIONToolkit/GPEC/blob/develop/docs/development/regression-harness.md). + +## Developer tips + +### Precompile workload + +Building the package runs a small Solovev case through the solver so that a fresh Julia process +skips most first-call compilation, at the cost of a longer one-time build. Code compiled into the +package image can round differently in the last bit from code compiled on first use, and +near-marginal quantities such as Δ′ amplify that, so results are reproducible only within one +build mode. The regression harness runs with the workload disabled. To disable it yourself (for a +faster build, or to match harness numbers), set the preference by UUID from the project, without +loading GPEC (which would first build it with the workload on). `Preferences` must be installed, +e.g. in your global environment: + +```julia +using Preferences, UUIDs +set_preferences!(UUID("462872dd-e066-4d2e-b993-6468b5239634"), "precompile_workload" => false; force=true) +``` + +This writes `LocalPreferences.toml` next to the active `Project.toml` (the file is gitignored); +delete the entry to re-enable. The regression harness overrides this file for its own runs. + +The workload runs only when Julia rebuilds the package image, not at the start of each run. That +happens after any change under `src/` (a pull, a branch switch or a local edit) and after a change +of dependency versions, Julia version or this preference. Repeated runs on an unchanged checkout +reuse the cached image. Julia keeps a separate image per set of compile flags, so the first +`Pkg.test()` after such a change (which checks array bounds) runs the workload once more. If you edit `src/` and restart Julia often, disabling the workload keeps +each rebuild short; within one Revise session nothing is rebuilt. + +### Revise.jl + +When iterating on code, use [Revise.jl](https://timholy.github.io/Revise.jl/stable/) to avoid full recompilation on every change. It tracks source file modifications and recompiles only the affected code. + +Install Revise in your **global** Julia environment (not in this project's `Project.toml`): + +```julia +# In the Julia REPL (no --project flag) +using Pkg +Pkg.add("Revise") +``` + +Then load it before `GeneralizedPerturbedEquilibrium` in any session: + +```julia +using Revise +using GeneralizedPerturbedEquilibrium +``` + +For it to load automatically, add `using Revise` to your [Julia startup file](https://timholy.github.io/Revise.jl/stable/config/#Using-Revise-by-default): + +```julia +# ~/.julia/config/startup.jl +try + using Revise +catch e + @warn "Could not load Revise" exception=e +end +``` + +### Pre-commit hooks + +The repository uses pre-commit hooks to maintain code quality. They run automatically on `git commit`, and +the same hooks run in CI on every pull request (the `pre-commit` check), which fails if they would change a file. +They: +- Format Julia code according to `.JuliaFormatter.toml` settings +- Remove trailing whitespace and fix line endings +- Validate YAML/TOML syntax +- Prevent accidentally committing large files (>5MB) +- Check config-style TOML annotation conventions (header block, inline `# description` + comments, no deprecated keys — see `docs/development/toml-conventions.md`) + +**Installation:** + +```bash +# Install pre-commit (requires Python/pip) +pip install pre-commit + +# Install JuliaFormatter 1.0.62 in the default environment: the hook runs whatever version is there, +# and CI pins this one, so another version formats differently and fails the check +julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="1.0.62"))' + +# Install the git hooks in your local repository +cd /path/to/GPEC +pre-commit install +``` + +You can also run the hooks manually without committing: + +```bash +pre-commit run --all-files # run on all files +pre-commit run # run on staged files only +``` + +In rare cases where you need to bypass the hooks: + +```bash +git commit --no-verify +``` + +This is discouraged as it may introduce formatting inconsistencies. diff --git a/docs/src/equilibrium.md b/docs/src/equilibrium.md index a7215f1de1..f58f644318 100644 --- a/docs/src/equilibrium.md +++ b/docs/src/equilibrium.md @@ -1,192 +1,102 @@ -# Equilibrium Module - -The Equilibrium module provides tools to read, construct, and analyze -magnetohydrodynamic (MHD) plasma equilibria. It wraps a variety of -equilibrium file formats (EFIT, CHEASE, SOL, LAR, and others), runs the -appropriate direct or inverse solver, and returns a processed -`PlasmaEquilibrium` object ready for downstream calculations. - -## Overview - -Key responsibilities of the module: - -- Read equilibrium input files and TOML configuration (see `EquilibriumConfig`). -- Provide convenient constructors for analytic / model equilibria (Large - Aspect Ratio, Solovev). -- Build spline representations used throughout the code (1D cubic and - 2D bicubic splines). -- Run the direct or inverse equilibrium solver and post-process results - (global parameters, q-profile finding, separatrix finding, GSE checks). - -The module exposes a small public API that covers setup, configuration, -and common analyses used by other GPEC components (e.g. `ForceFreeStates`, vacuum -interfaces). - -## EFIT solver strategies - -When `eq_type` is one of the EFIT-based options, three solver strategies are available: - -| `eq_type` | Method | Best for | -|-----------|--------|----------| -| `"efit"` | Geometric-angle field-line ODE | Default; fast and robust for standard cases | -| `"efit_arclength"` | Arc-length field-line ODE | Near-separatrix surfaces where the geometric-angle ODE becomes singular | -| `"efit_by_inversion"` | Contour.jl marching-squares → inverse solver | Highest geometric accuracy; avoids ODE singularities entirely | - -**`efit`** integrates field lines using the geometric angle η (polar angle from the -magnetic axis, 0→2π) as the independent variable. Position at each step is computed as -`R = R₀ + rfac·cos(η)`, so the RHS denominator is `Bz·cos(η) − Br·sin(η)` — the -projection of B_pol onto the radial direction. Near an X-point where Bp → 0 this -denominator passes through zero, causing a coordinate singularity: the ODE steps -become extremely small or the solver fails to converge. - -**`efit_arclength`** uses arc length along the flux surface as the independent variable -instead. The tangent direction `(dR/ds, dZ/ds) = ∇ψ⊥/|∇ψ|` has unit magnitude by -construction, so the position equations have no denominator and remain well-behaved -all the way to the separatrix. The 1/Bp factors in the accumulated integrals still -diverge near X-points, but their solver tolerances are set large so they do not -restrict the step size — only the position tracking drives adaptivity. The two methods -produce identical results when both succeed; `efit_arclength` extends the reliable -range of psihigh closer to 1. - -**`efit_by_inversion`** traces flux surface level sets directly from ψ(R,Z) using marching -squares, resamples each closed curve to a uniform geometric-angle grid, and feeds the -result into the inverse equilibrium solver (the same path used by CHEASE input). -The Cartesian evaluation grid is clipped to the separatrix bounding box and its -resolution is set adaptively from a bilinear interpolation error bound, so no manual -tuning is needed. - -## Radial grid packing - -With `grid_type = "auto"` and `mpsi = 0` (the defaults), the radial grid is -built by a **two-pass measured-curvature refinement** driven by the single accuracy knob -`psi_accuracy` (τ): - -1. **Pass 1** forms the equilibrium on a coarse three-region layout (geometric in - log(ψ) at the core, uniform in the middle, geometric in log(1−ψ) at the edge). -2. A knot density is derived from **one sizing rule** — the cubic-spline derivative - error model err(f′) ≈ h³|f''''|/24 ≤ τ·|f| — applied to three sources of - curvature |f''''|: - - **Measured** (mid-radius and edge): divided differences on the pass-1 nodes of the - 1D profiles (F, P, dV/dψ, q), the 2D geometry channels (r², η-offset, ν, Jacobian) - along sampled θ-lines (the bicubic ψ-axis shares the same knots), and the kinetic - profiles (n, T, ω_E) when loaded, so steep pedestal gradients attract knots. Nodal - values come from independent field-line integrations, so the estimate is immune to - inter-knot spline ringing. Curvature is measured against ρ = √ψ, in which the - equilibrium is regular at the magnetic axis — in ψ itself the geometry channels - behave as ψ^(k/2) (R−R₀ ~ √ψ), so their ψ-curvature diverges under refinement while - their ρ-curvature converges; the ρ-spacing maps back through dψ/dρ = 2√ψ, which by - itself packs √ψ-tight toward the axis. - - **Separatrix model** (edge floor, ψ ≥ 0.9): the same rule on q ≈ −A·ln(1−ψ) gives - geometric-in-log(1−ψ) packing with uniform relative q′ error, independent of A. - Normally inactive — the pass-1 layout is already log-packed at the edge, so the - measured density dominates — it only guards a pass-1 that under-sampled the - divergence. - - **Axis model** (core, ψ ≤ 0.03): the same rule on the power-law axis form gives - geometric-in-log(ψ) packing (constant-ratio spacing ∝ ψ down to `psilow`). Here the - model *replaces* measurement: nodal data from the smallest flux surfaces is - dominated by integration and axis-extrapolation noise, which grid refinement - amplifies rather than resolves. - - **Rational surfaces**: local packing around every q = m/n in the requested n range - (pinned as mandatory knots by the main driver), at any ψ including inside the - modeled core region. This local resolution of the ψ-splined Euler-Lagrange - coefficient matrices is what converges the Δ′ boundary-value problem. -3. The density integral is equidistributed into the knot vector, mandatory - rational-surface knots are inserted with a minimum-spacing snap guard, and - **pass 2** re-forms the equilibrium on the refined grid from the in-memory input - (no file re-read). - -Every region's knot count scales as τ^(-1/3), so tightening `psi_accuracy` refines -the core, pedestal, and edge proportionally. The fixed `grid_type = "core_edge_packed"` -(sin²-spaced), `"pow1"`, `"uniform"`, and explicit `mpsi > 0` (single-pass, fixed -layout) are still supported. A `PlasmaEquilibrium` built directly with `mpsi = 0` keeps -the coarse pass-1 grid until `solve` re-forms it; to refine it yourself, pass the nodes from -`refined_psi_grid` to `Equilibrium.reform`. - -The packing on the DIII-D-like example (n=1) compared to fixed `core_edge_packed` grids — note the -coarse spacing across the smooth mid-radius, the spacing dips at each rational surface, -and the core/pedestal/edge packing (`benchmarks/plot_grid_knot_placement.jl` regenerates -this figure): - -![Radial knot packing: auto two-pass vs core_edge_packed](assets/grid_knot_placement.png) - -Decomposing the density by source on the same example shows the pedestal band -(ψ_N ≈ 0.85–0.98) is driven by *measured* curvature, not the edge floor: the pressure, -q, and dV/dψ profiles contribute comparably, and the rzphi geometry channels — the -Grad-Shafranov response to the same pedestal p′ (Shafranov-shift / angle-offset -steepening) — contribute the most at every node in the band. Because the density is -measured from the formed solution, the packing follows the pedestal wherever it sits: - -![Knot-density decomposition by source](assets/density_decomposition.png) - -## API Reference +# Equilibrium -```@autodocs -Modules = [GeneralizedPerturbedEquilibrium.Equilibrium] -``` +`Equilibrium` reads or constructs the axisymmetric MHD equilibrium, maps it to the straight-field-line coordinates every later stage works in, and returns a `PlasmaEquilibrium` holding the flux-surface splines, the profiles and the global parameters of the plasma. +It is the first stage of every run and the only one that reads physics input files: the equilibrium itself and, optionally, the kinetic profiles. -## Important types - -- `EquilibriumConfig` — configuration container, built from the `[Equilibrium]` - section of `gpec.toml` or from keywords. -- `PlasmaEquilibrium` — the runtime structure containing spline fields, - geometry, profiles, and computed diagnostics (q-profile, separatrix, - etc.). -- `SolovevConfig`, `LargeAspectRatioConfig`, `TJAnalyticConfig` — the parameters of the - analytic equilibria. A deck states one as a `[Solovev]`, `[LargeAspectRatio]` or - `[TJAnalytic]` section, which selects the equilibrium, so it sets no `eq_type`. - `[TJAnalytic] solver = "inverse" | "direct"` picks the TJ pipeline. - -## Key functions - -- `PlasmaEquilibrium(file; eq_type="efit", kwargs...)`, `PlasmaEquilibrium(analytic_config; kwargs...)` - and `PlasmaEquilibrium(dd; kwargs...)` — read or solve the equilibrium and return a - `PlasmaEquilibrium`; `kwargs` are `EquilibriumConfig` fields. -- `PlasmaEquilibrium(config::EquilibriumConfig, source=nothing)` — the form every other - constructor ends in, with `source` an analytic config, an IMAS `dd`, or `nothing` to read `eq_filename`. -- `equilibrium_separatrix_find!(pe::PlasmaEquilibrium)` — locate - separatrix and related boundary geometry in-place. -- `equilibrium_global_parameters!(pe::PlasmaEquilibrium)` — populate - common global parameters (major radius, magnetic axis, volumes, etc.). -- `equilibrium_qfind!(pe::PlasmaEquilibrium)` — compute safety factor - (q) information across the grid. -- `equilibrium_gse!(pe::PlasmaEquilibrium)` — diagnostics on the - Grad–Shafranov solution. - -## Example usage - -Basic example: read a g-file and build an equilibrium +## What it computes -```julia -using GeneralizedPerturbedEquilibrium +From the poloidal flux ``\psi(R, Z)`` and the profiles ``F = R B_\phi`` and ``p`` of a Grad–Shafranov solution, the module traces the closed flux surfaces, finds the magnetic axis and the separatrix, integrates the safety factor ``q(\psi)`` and the surface quantities ``dV/d\psi``, and builds bicubic splines of the surface geometry ``R(\psi, \theta)``, ``Z(\psi, \theta)`` and of the toroidal-angle offset ``\nu(\psi, \theta)`` in the coordinates selected by `jac_type` (`"hamada"`, `"pest"`, `"boozer"`, `"equal_arc"`, …; see the [Conventions Reference](conventions.md) page). +The normalized flux ``\psi_N \in [0, 1]`` runs from the axis to the boundary, and the rational surfaces q = m/n of the run are located on it. + +## Inputs + +`[Equilibrium] eq_type` selects the file format of `eq_filename`: + +| `eq_type` | Format | Notes | +|---|---|---| +| `"efit"` | EFIT g-file | flux surfaces traced by a field-line ODE in the geometric angle; the default | +| `"efit_arclength"` | EFIT g-file | the same traced in arc length, which stays regular near the separatrix and reaches higher `psihigh` | +| `"efit_by_inversion"` | EFIT g-file | level sets traced by marching squares and fed to the inverse solver; the most accurate geometry | +| `"chease"` (`"chease_binary"`) | CHEASE binary output | inverse solver | +| `"chease2"` (`"chease_ascii"`) | CHEASE ASCII output | inverse solver | +| `"ifile"` | inverse i-file, as TokaMaker writes | inverse solver | +| `"imas"` | an `IMASdd.dd` passed from a script | `PlasmaEquilibrium(dd; kwargs...)`; see `examples/DIIID-like_ideal_example_IMAS` | + +An analytic equilibrium is selected by its own section instead of `eq_type`: `[Solovev]` (the Solov'ev solution), `[LargeAspectRatio]` (a large-aspect-ratio model with prescribed profiles) or `[TJAnalytic]` (the profile family of Fitzpatrick's TJ code, `solver = "inverse" | "direct"`); their keys are the fields of `SolovevConfig`, `LargeAspectRatioConfig` and `TJAnalyticConfig`. + +### Kinetic profiles + +`[Equilibrium] kinetic_file` attaches the density, temperature, rotation and transport profiles that the kinetic solve, [Kinetic Forces](kinetic_forces.md) and the `SLAYER` layer of [Tearing](inner_layer.md) read, so every kinetic consumer sees one plasma. +Two formats are read, by extension: + +- **HDF5** (`.h5`, `.hdf5`), the GPEC kinetic schema: datasets at the file root, each with a `units` attribute, on the normalized flux `psi`: required `n_e`, `T_i`, `T_e`, `omega_E` (and `n_i`, defaulting to `n_e`); optional `omega_tor`, `chi_e`, `chi_i`, `chi_phi` (perpendicular heat and momentum diffusivities), `zeff`, and the supplied collision rates `lnLambda`, `nu_i` and `nu_e` (see [Kinetic Forces](kinetic_forces.md)); further `n_*` datasets are named per-species densities. Densities in m⁻³, temperatures in eV, frequencies in rad/s, collision rates in 1/s, diffusivities in m²/s. `write_kinetic_h5` writes the format. +- **ASCII** (`.gpeckf`, `.kin`, `.dat`): the six-column table `psi_n n_i n_e T_i T_e omega_E`, or with `kinetic_format = "rates-v1"` nine columns that append `lnLambda nu_i nu_e` (the Fortran PENTRC rates-v1 table); the column count must match `kinetic_format`, and header rows are skipped. -pe = PlasmaEquilibrium("examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk"; jac_type="hamada") +The ion composition is part of the same data: `[[Equilibrium.ion_species]]` lists the main ions, each with `z`, `m` and either `fraction` (its share of the file's total `n_i`) or `density` (a named `n_*` dataset), and `impurity_charge` and `impurity_mass` describe the impurity that closes quasineutrality; the default is one deuterium ion carrying all of `n_i` with a carbon impurity. +``Z_\mathrm{eff}(\psi)`` is taken from the `zeff` dataset when present, else from `n_i/n_e` with the impurity charge. +In a script the same data is attached with `attach_kinetic_profiles!(equil, "kinetic.h5"; ion_species=[IonSpecies(; z=1, m=2, fraction=0.5), IonSpecies(; z=1, m=3, fraction=0.5)], impurity_charge=6, impurity_mass=12)`. -println("Magnetic axis: ", pe.ro, ", ", pe.zo) -println("q(0) = ", pe.params.q0) +### Radial grid -# Find separatrix (in-place) and inspect results -GeneralizedPerturbedEquilibrium.Equilibrium.equilibrium_separatrix_find!(pe) -println("rsep = ", pe.params.rsep) +With `grid_type = "auto"` and `mpsi = 0` (the defaults) the radial knots are placed by a two-pass refinement driven by one accuracy knob, `psi_accuracy`: the equilibrium is formed on a coarse grid, the curvature of the profiles, of the surface geometry and of the kinetic profiles is measured on it, and the knot density that keeps the cubic-spline derivative error below the tolerance is equidistributed, with geometric packing toward the axis and the separatrix and a mandatory knot at every rational surface of the run. +Every region's knot count scales as `psi_accuracy^(-1/3)`, so tightening it refines the core, the pedestal and the edge together; the rational-surface knots are what converge the Δ′ boundary-value problem. +The fixed layouts `grid_type = "core_edge_packed"`, `"pow1"` and `"uniform"` with an explicit `mpsi` remain available. + +![Radial knot packing on the DIII-D-like example: the automatic two-pass grid against a fixed core_edge_packed grid.](assets/grid_knot_placement.png) + +## Running it + +From `examples/DIIID-like_ideal_example/gpec.toml`: + +```toml +[Equilibrium] +eq_type = "efit" +eq_filename = "TkMkr_D3Dlike_Hmode.geqdsk" +kinetic_file = "TkMkr_D3Dlike_Hmode_kinetic.h5" +jac_type = "hamada" +grid_type = "auto" +psilow = 1e-4 # innermost flux surface +psihigh = 0.995 # outermost flux surface +mpsi = 0 # 0: automatic radial grid +psi_accuracy = 0.001 # spline accuracy that sizes the automatic grid +mtheta = 256 # poloidal resolution +impurity_charge = 6 +impurity_mass = 12 + +[[Equilibrium.ion_species]] +z = 1 +m = 2 +fraction = 1.0 ``` -Analytic / testing example: construct a large-aspect-ratio model +The same from a script, with the deck keys as keywords: ```julia using GeneralizedPerturbedEquilibrium -larcfg = GeneralizedPerturbedEquilibrium.Equilibrium.LargeAspectRatioConfig(r0=10.0, a=1.0, beta0=1e-3) -pe = PlasmaEquilibrium(larcfg) # the config type implies the equilibrium kind - -println("Built LAR equilibrium with a = ", larcfg.a) +eq = PlasmaEquilibrium("examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk"; jac_type="hamada", psihigh=0.995) +eq = PlasmaEquilibrium(Equilibrium.LargeAspectRatioConfig(; r0=10.0, a=1.0, beta0=1e-3)) # an analytic equilibrium ``` -## Notes and Caveats +Relative paths in the deck are resolved against the deck directory. +A `PlasmaEquilibrium` built with `mpsi = 0` keeps the coarse first-pass grid until a `solve` re-forms it; `Equilibrium.reform` with the nodes of `refined_psi_grid` does the same by hand. -- `EquilibriumConfig` is constructed from the `[Equilibrium]` section of `gpec.toml`. - Paths that are not absolute are resolved relative to the TOML file location. -- The Equilibrium module contains readers for EFIT and CHEASE formats. Ensure the - required data files are present and paths are set correctly in `gpec.toml`. +## Outputs -## See also +The group is written once `[ForceFreeStates]` has run, on the grid of the solve. -- `docs/src/stability.md` — ideal MHD stability analysis built on top of the equilibrium +| Dataset | Meaning | Units | +|---|---|---| +| `Equilibrium/R_axis`, `Z_axis`, `B_T_axis`, `R_mean`, `a_mean`, `aspect_ratio`, `kappa`, `delta_upper`, `delta_lower` | axis position and field, mean radii and shape parameters | m, T | +| `Equilibrium/I_p`, `volume` | plasma current and volume | A, m³ | +| `Equilibrium/beta_N`, `beta_t`, `beta_p_1`, `l_i_1`, … | normalized, toroidal and poloidal beta and internal inductance | | +| `Equilibrium/q_axis`, `q_95`, `q_edge`, `q_min`, `q_max`, `q_extrema_psi` | the safety factor at the axis, the 95% surface and the edge, and its extrema | | +| `Equilibrium/psi_boundary`, `psi_total`, `I_p_sign`, `B_T_sign` | the flux normalization and the sign conventions of the file | | +| `Equilibrium/Profiles/psi`, `q`, `2piF`, `mu0p`, `dVdpsi` | the 1-D profiles on ψ_N: q, 2πF, μ₀p and dV/dψ_N | T·m (`2piF`), T² (`mu0p`), m³ (`dVdpsi`) | +| `Equilibrium/Geometry/psi`, `theta`, `rcoords`, `offset`, `nu`, `jac` | the 2-D surface geometry on (ψ_N, θ): the radial coordinate, the angle offset, the toroidal offset ν and the Jacobian | | + +## API reference + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Equilibrium] +``` diff --git a/docs/src/error_fields.md b/docs/src/error_fields.md index 993b4a7471..c108780fd6 100644 --- a/docs/src/error_fields.md +++ b/docs/src/error_fields.md @@ -1,555 +1,178 @@ -# ErrorFields Module +# Error Fields -The `ErrorFields` module quantifies how sensitive a perturbed equilibrium's resonant drive is -to the placement of each coil set. It is the foundation of an error-field tolerance -assessment: once one plasma solve exists, the question *"how much resonant field does a -misaligned coil produce?"* is linear in the coil's motion, so the module linearizes it once -and everything downstream (tolerance Monte Carlo, locking risk, allowable-tolerance scans) -becomes linear algebra on a small table. +`ErrorFields` turns a coil-driven perturbed equilibrium into an error-field tolerance assessment: how much resonant field each coil set drives when it is misaligned, how likely the assembled machine is to lock for given manufacturing tolerances, and how much correction current an error-field correction array may and must carry. +Once one plasma solve exists the resonant drive is linear in each coil's motion, so the module linearizes it once and everything downstream is linear algebra on a small table. -## What is computed +## What it computes -For every named coil set of a coil-forced run, the module evaluates the root-area-weighted -control-surface spectrum ``\tilde{b}`` of the set as built and its central-difference -derivatives with respect to the six rigid-body degrees of freedom: Cartesian shifts -``(\Delta x, \Delta y, \Delta z)`` in metres and rotations ``(\theta_x, \theta_y, \theta_z)`` -about the machine axes in degrees. The spectra are placed on the run's ``(m, n)`` ordering and -conformed with the control-surface operator, so they are exactly the vectors the resonant -coupling matrix and its singular vectors act on -(see [Dominant resonant-coupling mode](perturbed_equilibrium.md#Dominant-resonant-coupling-mode)). +1. **Coil sensitivities** (`[ErrorFields]`, `ErrorFields/CoilSensitivities/`): the resonant spectrum of each coil set as built and its derivatives with respect to rigid shifts and tilts. +2. **Tolerances** (`tolerance_file`): the manufacturing tolerances of each coil set, in their own TOML file. +3. **Monte Carlo** (`[ErrorFields.MonteCarlo]`, `ErrorFields/MonteCarlo/`): the distribution of the dominant-mode overlap ``|\delta|`` over sampled misalignments, intrinsic and with error-field correction. +4. **Locking risk** (`[ErrorFields.scenario]`, `[ErrorFields.Risk]`, `ErrorFields/Risk/`): the overlap distribution against the empirical penetration-threshold scaling, giving the locking probability and its scan against tolerance scale, from which the allowable tolerance follows. +5. **Correction limits** (`[ErrorFields.NTV]`, `ErrorFields/NTV/`): how much correction current each array may carry before the NTV torque of its non-resonant field spends the rotation budget, and how much it needs. -The overlap of a coil set with singular mode ``k`` of the coupling matrix, normalized by the -axis toroidal field, is the dimensionless +For every named coil set of a coil-forced run the stage evaluates the root-area-weighted control-surface spectrum ``\tilde b`` of the set as built and its central-difference derivatives with respect to the six rigid-body degrees of freedom, the shifts ``(\Delta x, \Delta y, \Delta z)`` in metres and the rotations ``(\theta_x, \theta_y, \theta_z)`` about the machine axes in degrees, placed on the run's (m, n) ordering so that they are exactly the vectors the resonant coupling matrix and its singular vectors act on (see [Dominant resonant-coupling mode](perturbed_equilibrium.md#Dominant-resonant-coupling-mode)). +The stored primitive is the linearization of the spectrum, not a scalar, so the resonant surfaces retained, the singular mode and the normalization are post-hoc choices, and every later step can be re-evaluated from `gpec.h5` without a new plasma solve. +In the snippets below `EF` is `GeneralizedPerturbedEquilibrium.ErrorFields` and `AEF` is `GeneralizedPerturbedEquilibrium.Analysis.ErrorFields`. -```math -\delta = \frac{V_k^{\mathrm{H}}\,\tilde{b}}{B_{T0}}, -``` - -and because the projection is linear, the derivatives of ``\tilde{b}`` project to the -derivatives of ``\delta``. A rigid shift ``(\Delta x, \Delta y)`` therefore moves the overlap -by ``S_x \Delta x + S_y \Delta y`` with complex ``S_x = \partial\delta/\partial\Delta x`` and -``S_y = \partial\delta/\partial\Delta y``, which is exact for any coil shape; for an -axisymmetric coil ``S_y = \pm i S_x`` and the response reduces to a single magnitude with a -free phase, the model the OMFIT tolerance tool used. - -The stored primitive is the linearization of the spectrum, not a scalar, so the resonant -surfaces retained, the singular mode, and the normalization are all post-hoc choices. - -## Running it +## Governing equations -Add an `[ErrorFields]` section to a deck whose `[ForcingTerms]` uses -`forcing_data_format = "coil"` and whose `[PerturbedEquilibrium]` computes the singular -coupling: +The overlap of a coil set with singular mode ``k`` of the coupling matrix, normalized by the axis toroidal field, is the dimensionless -```toml -[ErrorFields] -fd_step_shift_m = 1e-3 # Central-difference step for the rigid shifts [m] -fd_step_tilt_deg = 0.1 # Central-difference step for the rigid tilts [degrees] -rotation_center = "conductor" # Tilt pivot: each conductor's own centre ("conductor") or the whole set's ("set") -exclude_coils = [] # Coil sets that are not error-field sources (the NTV efc_coils are excluded automatically) +```math +\delta = \frac{V_k^{\mathrm{H}}\,\tilde{b}}{B_{T0}}, ``` -The stage runs after the perturbed equilibrium and writes `ErrorFields/CoilSensitivities/`: -the spectra `field_as_designed`, `shift_sensitivity_per_m`, `tilt_sensitivity_per_deg`, a finite-difference -curvature diagnostic per tap, the current pattern the spectra were evaluated at, and a -`DominantMode/` summary projected onto the run's full-window dominant mode (`delta_as_designed`, -its shift and tilt sensitivities, their direction-averaged in-plane magnitudes, and the in-plane -shift and tilt that would cancel `delta_as_designed`). +and because the projection is linear, a rigid shift ``(\Delta x, \Delta y)`` moves it by ``S_x \Delta x + S_y \Delta y`` with the complex sensitivities ``S_x = \partial\delta/\partial\Delta x``, ``S_y = \partial\delta/\partial\Delta y``, exact for any coil shape; for an axisymmetric coil ``S_y = \pm i S_x`` and the response reduces to one magnitude with a free phase. +Over sampled misalignments the overlap of the machine is the coherent sum over coil sets of ``S\,(\Delta + u) + T\,(\theta + v)`` for each coil's draw ``(\Delta, \theta)`` and placement uncertainty ``(u, v)``, plus one shared draw per coherent group and a random-direction unattributed budget; no field is recomputed, so a million samples take about a second. -From a script, the same stage is one call on the plasma response and the coil sets of a coil-forced -perturbed equilibrium; it returns an `ErrorFieldsResult`, and products not asked for are `nothing`: +An overlap distribution becomes a locking risk through the empirical ITPA penetration-threshold scalings of [Logan et al. (2020)](citations.md#ErrorFields) for n = 1 and n = 2 and [Bursch et al. (2026)](citations.md#ErrorFields) for n = 1, -```julia -resp = plasma_response(ffs) -pe = perturbed_equilibrium(resp, RMPField(ft_ctrl)) -ef = error_fields(resp, pe.coil_sets, ForcingTerms.CoilConfig(ft_ctrl); - tolerances=ErrorFields.read_tolerance_toml("tolerances.toml"), - scenario=ErrorFields.ScenarioParameters(ffs.equil; n_e=5.0)) -save_to_hdf5("gpec.h5", ef) +```math +\delta_\mathrm{thresh} = 10^{\alpha_c}\, n_e^{\alpha_n}\, B_T^{\alpha_B}\, R_0^{\alpha_R}\, (\beta_N/l_i)^{\alpha_\beta}\, I_p^{\alpha_I}, ``` -`examples/DIIID-like_error_field_example/` is a complete case: the DIII-D-like equilibrium with -the C-coil as the nominal n = 1 source and the eighteen DIII-D F coils added as single-filament -hoops at the centroids of their winding packs (from OpenFUSIONToolkit's TokaMaker -`DIIID_geom.json`). An axisymmetric hoop drives no n = 1 field as built, so each F coil's -error-field content is entirely its sensitivity to misalignment; `analyze_example.jl` ranks the -coils by error field per millimetre of shift and per tenth of a degree of tilt, over all -rational surfaces and over the edge only. - -The tilt pivot matters for multi-filament winding packs: `"conductor"` rotates each filament -about its own arc-length centre, the Fortran `coil_read` convention inherited by the OMFIT -tolerance tool, while `"set"` rotates the pack rigidly about its common centre, which is what -an engineering axis-line tolerance constrains. Single-conductor sets give identical results -either way. - -## Tolerance input - -Manufacturing tolerances are engineering data, reviewed and versioned independently of any -run, so they live in their own TOML file named by `tolerance_file` in `[ErrorFields]` (a path -relative to the deck directory). The run validates the file, echoes its text into -`Input/RawInputs/ErrorFields/tolerance_toml_raw` for replay, and leaves the device numbers -where they belong: outside the repository. The sampling and Monte Carlo stages that consume -these tolerances follow in later releases; this release fixes the format. - -```toml -# Fallbacks for keys a coil block leaves out (every key optional) -[ErrorFields.defaults] -shift_sigma_mm = 0.0 # Gaussian uncertainty added to the sampled shift [mm] -tilt_sigma = 0.0 # Gaussian uncertainty added to the sampled tilt, in tilt_units -tilt_units = "deg" # "deg", or "m" for a rim displacement converted through the nominal radius -radial_shape = "hollow" # Radial sampling density on the disk: flat, uniform_area, hollow, or ring -tolerance_model = "additive" # "additive" (independent shift and tilt) or "cylinder" (correlated axis line) -cylinder_half_height_m = 0.0 # Cylinder half-height for the cylinder model [m] -current_factor = 1.0 # Current as built relative to the run's; scales the coil's overlap and sensitivities - -# One block per coil set that moves on its own -[[ErrorFields.coil]] -name = "PF1U" # Coil set name; must match a [[ForcingTerms.coil_set]] of the run -shift_tol_mm = 0.5 # Radius of the in-plane displacement disk the coil centre may lie in [mm] -tilt_tol = 0.019 # Tilt tolerance, in tilt_units -tolerance_model = "cylinder" # Axis line confined to a cylinder; shift and tilt drawn together -cylinder_half_height_m = 1.5 # Cylinder half-height; sets the tilt reachable at the given radius [m] - -# Coil sets that move together: one shared draw per sample -[[ErrorFields.coherent_group]] -name = "upper_pf_brace" # Label of the coherently moving group -members = ["PF1U", "PF2U"] # Coil sets sharing the draw -shift_tol_mm = 2.0 # Coherent shift amplitude shared by the group [mm] -tilt_tol = 0.0 # Coherent tilt amplitude shared by the group, in tilt_units -phase_group = "braces" # Groups naming the same phase_group share the random direction -rotation_center_z_m = 0.0 # Height of the pivot of the group's rigid rotation on the machine axis [m] -tilt_lever_arm_m = 3.2 # Lever arm converting a tilt_tol in metres to the group's rotation angle [m] (optional) - -[ErrorFields.correctability] -uncorrectable_coils = ["EFCC_U"] # Coil sets the error-field correction cannot reduce -efc_factor = 2.0 # Divisor applied to correctable contributions under error-field correction - -[ErrorFields.other_field] -magnitude = 8.2e-6 # Overlap budget of sources not attributed to any coil -sigma = 0.0 # Gaussian uncertainty on that budget -radial_shape = "ring" # Radial sampling density of the budget magnitude -``` +with ``n_e`` in 10¹⁹ m⁻³, ``B_T`` in T, ``R_0`` in m and ``I_p`` in MA (the current term only in the 2026 fits). +Sampling the fitted exponents within their standard errors turns the threshold into a distribution whose cumulative distribution is the probability ``P(\mathrm{lock}\,|\,\delta)`` that an overlap locks, and the locking probability of the assembled machine is ``100 \int \mathrm{pdf}(\delta)\, P(\mathrm{lock}\,|\,\delta)\, d\delta`` over the Monte Carlo bins, per batch. -The shift tolerance is one number, the radius of the disk the coil centre may sit in; the -direction is sampled. A tilt may be given in degrees or, as legacy tolerance tables do, as a -rim displacement in metres, which `tilt_tolerance_deg` converts through the coil set's -arc-length-weighted major radius exactly as the coil loader's `tilt_in_meters` does. A coherent -group's tilt is a rigid rotation of all its members about a pivot on the machine axis at -`rotation_center_z_m`, so a member at height `z` also shifts laterally by `(z − z_pivot)·θ`. -Groups sharing a `phase_group` label draw the same random direction with independent -amplitudes. A group's tilt given in metres is converted through `tilt_lever_arm_m` when the -group has one (a reference structure's lever arm is its own, not any member coil's radius), and -otherwise through the members' common major radius. A group is correctable only if none of its -members is listed as uncorrectable. Coil sets of the run without a tolerance block contribute -their as-designed overlap only. A coil's `current_factor` scales its as-designed overlap and its -sensitivities together, since the field is linear in current: `0` switches a coil off, `0.5` runs -it at half current, a negative value reverses it; this is the knob for a risk against the current -in an uncorrectable coil. - -Every key is checked against the schema, so a misspelled key is an error rather than a silent -default. `read_tolerance_toml` parses a file into a `ToleranceSet`, `validate_tolerances` -checks its names against a run's coil sets, and `update` copies a set, a coil, a group or the -unattributed budget with named fields replaced, which is how a programmatic scan over a -quantity the file fixes (a sigma, a budget, a current factor) is written: +A correction array cancels the dominant-mode overlap at ``C_c`` per kilo-ampere-turn, but the non-resonant remainder of its field drives an NTV torque ``T I^2`` that a perfect correction does not remove. +With a torque budget ``T_0`` and the threshold taken to fall in proportion to the torque spent, the current that corrects an intrinsic overlap ``\delta_\mathrm{EF}`` solves -```julia -ts = EF.read_tolerance_snapshot("gpec.h5") -half = EF.update(ts; coils=[c.name == "F6A" ? EF.update(c; current_factor=0.5) : c for c in ts.coils]) -tight = EF.update(ts; other_field=EF.update(ts.other_field; magnitude=1e-5)) +```math +\delta_\mathrm{EF} - C_c I = s\,\delta_\mathrm{thresh}\left(1 - \frac{T_\mathrm{residual} I^2}{T_0}\right), ``` -## Sampling a tolerance +the model of [Logan et al. (2026)](citations.md#ErrorFields) for SPARC and [Leuthold et al. (2026)](citations.md#ErrorFields) for ARC. +It holds only while the plasma still rotates, so the current is capped at ``I_\max = \sqrt{T_0/T_\mathrm{residual}}``, and the largest correctable overlap is the peak of the quadratic or the overlap cancelled at ``I_\max``, whichever comes first; the two zero-rotation limits ``C_c\sqrt{T_0/T_\mathrm{residual}}`` and ``C_c\sqrt{T_0/T_\mathrm{full}}`` are reported beside it. -A tolerance is one number, the radius of the disk the coil centre may lie in; the direction is -random and the radial density is a `RadialDistribution`: `Flat` (uniform in radius, the OMFIT -tool's `flat`), `UniformArea` (uniform over the disk), `Hollow` (peaked toward the edge, the -OMFIT default), `Ring` (always on the edge), or `PowerLaw(p)`. `sample_disk` draws a point as -`Δx + iΔy`, `sample_uncertainty` adds the Gaussian uncertainty on where the coil actually sits, -and the two tolerance models combine them: `sample_additive` draws a shift (metres) and a tilt -(degrees) independently, while `sample_cylinder` confines the coil's axis line to a cylinder -of radius `R` and half-height `z_top`, drawing its two endpoints and deriving the correlated -midplane shift and lean. Tilts are the rotation angles `θx + iθy` about the machine axes in -the sense `apply_transforms` uses. Every sampler takes the random generator and optional fixed -directions, so coherent groups sharing a direction pass one phase and draw their own radii. - -```julia -rng = Random.Xoshiro(1) -Δ = EF.sample_disk(rng, 0.5e-3, EF.Hollow()) # a point in a 0.5 mm disk, m -Δ, θ = EF.sample_cylinder(rng, 0.5e-3, 1.5, EF.randpow(EF.Flat())) # correlated shift [m] and tilt [deg] -``` +## Running it -## Tolerance Monte Carlo - -With a `tolerance_file` named, the run samples every coil set's misalignment within its -tolerance and histograms the dominant-mode overlap `|δ|`: per sample and coil set the overlap -moves by `S·(Δ + u) + T·(θ + v)` for the coil's own draw (`Δ`, `θ`) and Gaussian placement -uncertainties (`u`, `v`); coherent groups add one shared draw per group, with the lateral shift -a rigid rotation about the group pivot gives each member; the unattributed budget adds a random -direction. Because the overlap is linear in the misalignments, a million samples take about a -second and no field is recomputed. Two histograms are written to `ErrorFields/MonteCarlo/`: the -intrinsic `|δ|` and the corrected one, in which every correctable term (coil sets and groups not -listed as uncorrectable, and the unattributed budget) is divided by `efc_factor`. Batches are -seeded individually, so results are bit-identical for any thread count, and their spread is the -statistical error bar of anything derived from them. - -A coil set the run energizes but that is not an error-field source — a correction array named -in `[ErrorFields.NTV] efc_coils`, or anything listed in `[ErrorFields] exclude_coils` — is swept -like the others, so its sensitivities and couplings are tabulated, but it is left out of the -as-designed error field, the Monte Carlo, the worst-case bound and the scans, in the run and in -every post-hoc entry point that rebuilds the table from the file. A tolerance that names it is an -error. `without_coils(table, names)` is the same operation on a table in memory, and -`excluded_coil_names("gpec.h5")` says which sets a run left out. +The stage needs a `[PerturbedEquilibrium]` run with `forcing_data_format = "coil"` and the singular coupling computed; `examples/DIIID-like_error_field_example` is the complete case, with the DIII-D C-coil as the n = 1 source and the eighteen F coils as single-filament hoops whose error-field content is entirely their sensitivity to misalignment. +From its deck: ```toml [ErrorFields] -tolerance_file = "tolerances.toml" # Manufacturing-tolerance TOML, relative to the deck directory +fd_step_shift_m = 1e-3 # central-difference step of the rigid shifts [m] +fd_step_tilt_deg = 0.1 # central-difference step of the rigid tilts [degrees] +rotation_center = "conductor" # tilt pivot: each conductor's own centre, or "set" for the whole set's +exclude_coils = [] # coil sets that are not error-field sources (the NTV efc_coils are excluded automatically) +tolerance_file = "tolerances.toml" # manufacturing tolerances, relative to the deck directory [ErrorFields.MonteCarlo] -nsample = 1000000 # Samples per batch -nbatch = 10 # Independent batches; their spread is the statistical error bar -seed = 1 # Base seed; batch b uses Xoshiro(hash((seed, b))) -nbins = 300 # Histogram bins, linear on [0, delta_max] -delta_max = 0.0 # Upper histogram edge; 0 = 1.5 × the worst-case alignment bound -tolerance_scale = 1.0 # Multiplies every shift and tilt tolerance (for tolerance scans) -coil_subset = [] # Coil sets whose tolerances are sampled; empty = all -scale_subset = [] # Coil sets and groups that tolerance_scale applies to; the rest hold their tolerance (empty = all) -scale_map = {} # Extra multiplier per coil set or group name, e.g. {F6A = 2.0} -``` - -`coil_subset` and `scale_subset` answer different questions. `coil_subset` samples only the -named coils and gives every other coil zero tolerance, which isolates one coil's contribution. -`scale_subset` applies `tolerance_scale` to the named coils and groups and leaves every other -name at its own tolerance, which is the product question of a design review: how far can the -tolerance of one class of coils be relaxed while the rest of the machine holds its own? Names -in either that are not coil sets of the run or groups of the tolerance file are errors. Neither -changes which random numbers are drawn, so runs with the same seed differ only through the -tolerances they scale. +nsample = 200000 # samples per batch +nbatch = 10 # independent batches; their spread is the statistical error bar +seed = 1 +nbins = 300 # histogram bins on [0, delta_max]; delta_max = 0 is 1.5 × the worst-case bound -The run's histogram is the full-window, dominant-mode summary. Any other window or mode, a -tolerance scale, or a coil subset is a post-hoc re-run of the same kernel: - -```julia -mc = EF.run_monte_carlo("gpec.h5"; psi_low=0.5, tolerance_scale=2.0, coil_subset=["PF1U", "PF2U"]) -mc.abs_delta_pdf, mc.abs_delta_bin_edges # intrinsic |δ| density -mc.abs_delta_efc_pdf # corrected -mc.abs_delta_sampled_mean, mc.abs_delta_total_as_designed -``` - -## Locking risk and allowable tolerance - -An overlap distribution becomes a locking risk through the empirical ITPA penetration-threshold -scalings ([citations](citations.md#ErrorFields-Module)): the n = 1 fits of Logan et al., -*Plasma Phys. Control. Fusion* **62**, 084001 (2020) and of Bursch et al., *Plasma Phys. Control. -Fusion* (2026), [doi:10.1088/1361-6587/aea7d6](https://doi.org/10.1088/1361-6587/aea7d6), and the n = 2 fits of Logan et al., -*Nucl. Fusion* **60**, 086010 (2020): -`δ_thresh = 10^α_c · n_e^α_n · B_T^α_B · R_0^α_R · (β_N/l_i)^α_β · I_p^α_I`, with `n_e` in -10¹⁹ m⁻³, `B_T` in T, `R_0` in m and `I_p` in MA. Each fit is chosen by `year`, `dataset` and `fit`; only the 2026 fits carry the current term: - -| Fit (`year`, `dataset`, `fit`) | α_c | α_n | α_B | α_R | α_β | α_I | -|---|---|---|---|---|---|---| -| `2026`, `"O,L"`, `"OLS"` (Eq. 7) | −4.31 ± 0.09 | 0.77 ± 0.08 | 0.19 ± 0.09 | 1.88 ± 0.16 | 0.25 ± 0.08 | −0.97 ± 0.08 | -| `2026`, `"O,L"`, `"WLS"` (Eq. 8) | −4.26 ± 0.09 | 0.56 ± 0.08 | 0.30 ± 0.10 | 1.57 ± 0.15 | 0.13 ± 0.06 | −1.01 ± 0.07 | - -The 2026 fits use only ohmic and L-mode discharges of conventional tokamaks (C-Mod, DIII-D, -EAST, JET, J-TEXT and KSTAR; no NSTX or COMPASS), so they suit conventional-aspect-ratio designs -in linear ohmic confinement or L-mode. Sampling the fitted exponents -within their standard errors turns the threshold into a distribution; its cumulative -distribution is the probability that an overlap `δ` locks, and the locking probability of the -assembled machine is `100 ∫ pdf(δ) P(lock|δ) dδ` over the Monte Carlo bins, per batch. The -operating point is an `[ErrorFields.scenario]` table: density must be given (it is not an -equilibrium output); field, major radius, β_N, l_i and the plasma current default from the -equilibrium. - -```toml [ErrorFields.scenario] -n_e = 5.0 # Electron density for the threshold scaling [1e19 m^-3] -# i_p = 1.2 # Plasma current magnitude [MA]; defaults to the equilibrium's +n_e = 5.0 # electron density of the threshold scaling [1e19 m^-3]; B_T, R_0, beta_N, l_i and I_p default from the equilibrium [ErrorFields.Risk] -year = 2020 # Publication year of the threshold fit: 2020 (n = 1 and n = 2) or 2026 (n = 1) -dataset = "O,L" # ITPA dataset of the fit: 2020 "O,L" or "O,L,H" (n = 1), "O,L", "O,L,-C", "O,L,N" (n = 2); 2026 "O,L" -fit = "WLS" # Fitting method: "OLS", "DSOLS", or "WLS" (2026: "OLS" or "WLS") -distribution = "normal" # How the fit exponents are sampled: "normal", "flat", or "normal_truncated" -nsample_threshold = 1000000 # Threshold samples -seed = 1 # Seed of the threshold sampling -scan_scales = [0.25, 0.5, 1.0, 2.0, 4.0] # Tolerance multipliers of the allowable-tolerance scan (empty: no scan) -``` - -`ErrorFields/Risk/` holds the threshold density and `P(lock|δ)` on the Monte Carlo grid, the -locking probability of the intrinsic and corrected distributions (with per-batch values), the -as-designed risk, and the sharp-threshold risk; `ErrorFields/Risk/ToleranceScan/` the risk -against tolerance scale. The scan is the stored quantity; the allowable tolerance for a target -risk is a post-hoc inversion, and every window or fit choice is re-evaluated from the file: - -```julia -scan = EF.ToleranceScan("gpec.h5") -EF.allowable_tolerance(scan, 1.0) # tolerance multiplier at 1 % locking risk -EF.allowable_tolerance(scan, 1.0; corrected=true) # with error-field correction -risk = EF.locking_risk("gpec.h5"; n_e=5.0, psi_low=0.7, risk_ctrl=EF.RiskControl(; dataset="O,L,H")) -risk26 = EF.locking_risk("gpec.h5"; n_e=5.0, risk_ctrl=EF.RiskControl(; year=2026, fit="OLS")) -scan2 = EF.tolerance_scan("gpec.h5"; n_e=5.0, scales=[0.5, 1, 2, 4], coil_subset=["F6A", "F7A"]) -scan3 = EF.tolerance_scan("gpec.h5"; n_e=5.0, scales=[0.5, 1, 2, 4], scale_subset=["F6A", "F7A"]) # scan the outboard pair, hold the rest -``` +year = 2020 # threshold fit: 2020 (n = 1 and 2) or 2026 (n = 1) +dataset = "O,L" # ITPA dataset of the fit ("O,L", "O,L,H", ...) +fit = "WLS" # "OLS", "DSOLS" or "WLS" +nsample_threshold = 1000000 +scan_scales = [0.05, 0.1, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0, 16.0] # tolerance multipliers of the allowable-tolerance scan -A threshold scaling is fitted for one toroidal mode number, so the in-run risk and every -file-path entry point refuse a run that spans several (`single_toroidal_mode`); project onto one -n's coupling first. - -A number quoted at the target risk level should be shown free of sampling and binning bias. -`risk_convergence` repeats the Monte Carlo over a list of sample counts and a list of bin -counts, holding the other at its control value, and returns the locking probability with its -batch spread for each: the spread must fall as `1/√nsample`, and the value must not move with -`nbins` by more than the spread. The spread alone cannot see a binning bias, which is why the bin -sweep is there. - -```julia -conv = EF.risk_convergence(table, ts, coil_sets, sc, scen; nsamples=[10^5, 10^6, 10^7], nbins_list=[100, 300, 3000]) -conv.locking_probability_percent_by_nsample, conv.locking_probability_spread_percent_by_nsample -conv.locking_probability_percent_by_nbins -``` - -## Plots and coil-array phasing - -`Analysis.ErrorFields` plots everything above from `gpec.h5`, and every function takes a list -of `label => path` pairs so coil-design revisions overplot on one axis: - -```julia -AEF = GeneralizedPerturbedEquilibrium.Analysis.ErrorFields -AEF.plot_coil_sensitivities(["rev A" => "revA/gpec.h5", "rev B" => "revB/gpec.h5"]; quantity=:shift) -AEF.plot_tolerance_pdf("gpec.h5"; corrected=true) -AEF.plot_locking_risk("gpec.h5"; target_percent=1.0) # marks the allowable scale -AEF.plot_threshold_scaling("gpec.h5") -AEF.plot_dominant_mode_spectrum("gpec.h5") -AEF.plot_error_field_summary("gpec.h5"; save_path="error_field_summary.png") -``` - -Four more views answer the questions the assessment's single numbers hide. `plot_overlap_phasors` -lays each coil's complex overlap head to tail, so a coil that adds to the error field is told -apart from one that cancels another's; `plot_tolerance_budget` stacks the terms of the worst-case -bound `abs_delta_worst_case` (`EF.worst_case_terms`) per coil, group and unattributed budget, so the -tolerance the budget is spent on is visible; `plot_linearity_residuals` shows the curvature the -linear sensitivity model neglects, at the finite-difference step or rescaled to each coil's own -tolerance; and `quantity=:fraction` on the sensitivity bars is the resonant share of each coil's -own spectrum, which separates a coil that is small from one that drives the wrong harmonics. - -```julia -AEF.plot_overlap_phasors("gpec.h5") # one panel per source -AEF.plot_tolerance_budget("gpec.h5"; sort=:total, top=10) # largest budget terms first -AEF.plot_linearity_residuals("gpec.h5"; at=:tolerance) # curvature over the tolerance range -AEF.plot_coil_sensitivities("gpec.h5"; quantity=:fraction) # resonant fraction, in percent -AEF.plot_coil_sensitivities("gpec.h5"; quantity=:as_designed, yscale=:log10) # stems, never bars, on a log axis -AEF.plot_tolerance_pdf("gpec.h5"; show_batches=true) # batch spread and the clamped fraction -EF.worst_case_terms("gpec.h5").total # = ErrorFields/MonteCarlo/abs_delta_worst_case -``` - -When several independently powered coil arrays share the job of correcting the error field, -the relative phases of their current patterns decide how much dominant-mode field they can -drive per ampere-turn. `phasing_map` evaluates `|Σ_k δ_k e^{iφ_k}|` per kilo-ampere-turn and -the resonant fraction of the applied field on a grid of the `N − 1` relative phases from the -stored nominal spectra, a closed form with no optimizer, and `plot_phasing_map` draws it (a -line for two arrays, a contour for three): - -```julia -pmap = EF.phasing_map("gpec.h5", ["EFCC_L", "EFCC_M", "EFCC_U"]; psi_low=0.5) -EF.extreme_phasing(pmap) # best |δ| per kAt and the phases giving it -AEF.plot_phasing_map(pmap; quantity=:resonant_fraction_percent) -``` - -## NTV limits of error-field correction - -A correction coil cancels the dominant-mode overlap at `C_c` per kilo-ampere-turn, but the -non-resonant remainder of its field drives a neoclassical toroidal viscosity (NTV) torque -`T·I²` that a perfect correction does not remove. With a torque budget `T_0` and the threshold -taken to fall in proportion to the torque spent, the current that corrects an intrinsic overlap -`δ_EF` solves `δ_EF − C_c I = s δ_thresh (1 − T_residual I²/T_0)`. This is the model used for -SPARC by Logan et al., Nucl. Fusion (2026), doi:10.1088/1741-4326/ae6086, and for ARC by -Leuthold et al., J. Plasma Phys. 92, E49 (2026), doi:10.1017/S0022377826101421. It holds only -while the plasma still rotates, so the current is capped at `I_max = √(T_0/T_residual)`, and the -largest correctable overlap is whichever comes first: the peak of the quadratic, or the overlap -cancelled at `I_max`. `max_correctable_overlap` returns that limit (`with_ntv`) together with the two -published zero-rotation limits, `C_c √(T_0/T_residual)` (SPARC, `residual_only`) and -`C_c √(T_0/T_full)` (ARC, `torque_only`). `[ErrorFields.NTV]` names the correction arrays; the run -evaluates each one's `C_c`, resonant fraction, and NTV torque per kAt² for its whole field and -for its field with the dominant mode projected out — two plasma-response evaluations of the -unit-current spectrum followed by the kinetic torque, which needs a `[KineticForces]` section — -and writes `ErrorFields/NTV/`. The spectrum is evaluated on the run's `[ForcingTerms]` boundary -grid and normalized by the magnitude of the array's ampere-turns, `|nw| × max|I|`, so the winding -sense and current pattern of the deck stay the phase reference; the torques are stored with their -sign, and the limits consume the budget with their magnitude (the sign depends on the rotation -and on conventions, so a negative torque is never read as no torque). Torque budget, threshold -and safety factor are analysis choices: - -```toml [ErrorFields.NTV] -efc_coils = ["d3d_c"] # Coil set names of the correction arrays to evaluate; never error-field sources -method = "fgar" # KineticForces torque method (must be enabled in [KineticForces]) -``` - -```julia -couplings = EF.read_efc_couplings("gpec.h5") -curve = EF.efc_current_curve(couplings[1]; delta_threshold=1.4e-4, torque_budget=4.0) -EF.max_correctable_overlap(couplings[1]; delta_threshold=1.4e-4, torque_budget=4.0) -AEF.plot_efc_ntv_limits("gpec.h5"; torque_budget=4.0) # threshold from the run's Risk/ group -# The SPARC paper's uncertainty: ±50 % on the NTV torque, T_0 = 4 ± 2 N·m -AEF.plot_efc_ntv_limits("gpec.h5"; torque_budget=4.0, torque_rtol=0.5, budget_rtol=0.5) -AEF.plot_efc_ntv_limits("gpec.h5"; torque_budget=4.0, profiles=true) # adds T(ψ) integrated from the axis +efc_coils = ["d3d_c"] # the correction arrays; never error-field sources +method = "fgar" # KineticForces torque method, enabled in [KineticForces] +rotation_scan = true # tabulate the torques against a rigid E×B rotation shift +rotation_scan_points = 5 +rotation_scan_max_points = 9 +rotation_scan_tolerance = 0.05 +rotation_span_factor = 1.0 +rotation_offset_factor = 2.0 ``` -The torque budget is the plasma's intrinsic torque. The SPARC paper uses 4 N·m for SPARC and the -ARC paper 5–20 N·m for ARC; pick a value for the machine being analysed. - -### Designing the correction - -The NTV limits say how much correction current an array *may* carry; `correction_requirement` -says how much it *needs*, and on which surfaces the correction falls short. For one array the -dominant-mode answer is the ratio `−δ_source/δ_array`; for several arrays the dominant-mode -condition alone does not fix the currents (one complex equation in several unknowns leaves a family -of solutions in which any two arrays can cancel each other, of which the result names the -member with the least total ampere-turns squared); the question worth a function is the joint one over every rational surface: the complex factors on all arrays that minimize the resonant -field `C·(b̃_source + Σ_k f_k b̃_k)` in the least-squares sense, the field left on each surface -after that and after each array's dominant-mode correction alone, and how much of the source's -spectrum each array can see. Each current comes twice: as a complex factor on the array's -current as given, and in kilo-ampere-turns through the array's `ampere_turns_kat` -(`|winding multiplier| × max|I| / 1000`, the same normalization as the EFC couplings), which is the -form to compare arrays whose patterns were given at different currents and to set beside the -NTV-limited current. Any `CoilOverlap` can be a source or an array, so an as-built assembly -(`combine_overlaps`) can be corrected with arrays that were never part of the run. Over the tolerance samples the needed -current is a distribution, and `uncorrectable_probability` is the fraction of the Monte Carlo's -overlap beyond what the NTV-limited array can correct at all: the explicit-correction counterpart -of the corrected locking probability's `efc_factor` model. The Monte Carlo leaves the -correction arrays out of its overlap, so its as-designed current and a `correction_requirement` -whose source combines the error-field coils describe the same field. +`[ErrorFields.MonteCarlo]` also takes `tolerance_scale` (multiplies every tolerance, for scans; `scale_subset` and `scale_map` restrict it) and `coil_subset` (the coil sets sampled; empty samples all), and `[ErrorFields.Risk]` takes `distribution`, how the fit exponents are sampled (`"normal"`, `"flat"`, `"normal_truncated"`). -```julia -req = EF.correction_requirement("gpec.h5", "F_coils_as_built", ["d3d_c", "iu", "il"]) -req.current_factor_least_squares, req.least_squares_rank # unique when the rank is the number of arrays -req.current_factor_dominant, req.current_factor_dominant_minimum_norm -req.current_least_squares_kat # the same in kilo-ampere-turns, comparable across arrays -abs.(req.resonant_field_least_squares_t) ./ abs.(req.resonant_field_source_t) # what is left, per surface -AEF.plot_correction_requirement(req) -mc = EF.MonteCarloResult("gpec.h5"); c = EF.read_efc_couplings("gpec.h5")[1] -EF.uncorrectable_probability(mc, c; delta_threshold=1.4e-4, torque_budget=4.0) -AEF.plot_needed_current(mc, array; coupling=c, delta_threshold=1.4e-4, torque_budget=4.0) -``` +The tilt pivot matters for multi-filament winding packs: `"conductor"` rotates each filament about its own centre, `"set"` rotates the pack rigidly about its common centre, which is what an engineering axis-line tolerance constrains. +A coil set that is not an error-field source (a correction array or anything in `exclude_coils`) is swept like the others but left out of the as-designed field, the Monte Carlo and the scans; a tolerance that names it is an error. +The NTV limits need a `[KineticForces]` section, since each array's torque per kAt² is one kinetic torque evaluation of its unit-current spectrum, with and without the dominant mode projected out. -## Analysis after the run +The tolerances are engineering data, reviewed independently of any run, so they live in their own file, echoed into `Input/RawInputs/ErrorFields/tolerance_toml_raw` for replay; every key is checked against the schema. +Its tables are `[ErrorFields.defaults]` (fallbacks for every coil block), one `[[ErrorFields.coil]]` per coil set that moves on its own, `[[ErrorFields.coherent_group]]` for coil sets that move together, `[ErrorFields.correctability]` and `[ErrorFields.other_field]`: -Window the coupling to any range of rational surfaces and project onto any singular mode -without re-running anything, from memory or from the file: +```toml +[ErrorFields.defaults] +shift_sigma_mm = 0.0 # Gaussian placement uncertainty added to the sampled shift [mm] +tilt_sigma = 0.0 # the same for the tilt, in tilt_units +tilt_units = "deg" # "deg", or "m" for a rim displacement converted through the coil's major radius +radial_shape = "hollow" # radial sampling density on the disk: flat, uniform_area, hollow or ring +tolerance_model = "additive" # "additive" (independent shift and tilt) or "cylinder" (correlated axis line) +current_factor = 1.0 # current as built relative to the run's; 0 switches a coil off -```julia -using GeneralizedPerturbedEquilibrium -EF = GeneralizedPerturbedEquilibrium.ErrorFields - -# From the file: the coupling is rebuilt, windowed to 0.5 ≤ ψ_N ≤ 1, and projected onto mode 1 -table = EF.sensitivity_table("gpec.h5"; psi_low=0.5) -table.delta_as_designed # complex overlap of each coil set -table.abs_delta_shift_per_mm # direction-averaged |∂δ/∂Δ| per millimetre of shift, per coil set -table.abs_delta_rim_per_mm # the same for tilt, as rim displacement at the coil's major radius -table.tilt_sensitivity_per_deg[1, :] # ∂δ/∂θx per degree, per coil set - -# In memory, from the run's returned state (no I/O) -rc = PerturbedEquilibrium.ResonantCoupling(run.pe.response) -dom = PerturbedEquilibrium.dominant_coupling(rc; psi_low=0.5) -table = EF.sensitivity_table(run.error_fields.sensitivities, dom; mode=1) -``` +[[ErrorFields.coil]] +name = "PF1U" # must match a [[ForcingTerms.coil_set]] of the run +shift_tol_mm = 0.5 # radius of the disk the coil centre may lie in [mm] +tilt_tol = 0.019 # tilt tolerance, in tilt_units +tolerance_model = "cylinder" +cylinder_half_height_m = 1.5 # sets the tilt reachable at the given radius -### How far the linear model holds +[[ErrorFields.coherent_group]] +name = "upper_pf_brace" +members = ["PF1U", "PF2U"] # one shared draw per sample +shift_tol_mm = 2.0 +tilt_tol = 0.0 +phase_group = "braces" # groups naming the same phase_group share the random direction +rotation_center_z_m = 0.0 # pivot of the group's rigid rotation on the machine axis [m] -Every sweep above assumes the overlap is linear in the misalignments out to the tolerance edge. -The run checks that only at the finite-difference step (the stored `*_linearity_residual` -ratios, `plot_linearity_residuals`). `linearity_check` displaces each coil set with a tolerance -block by `±scale × tolerance` along each in-plane shift and tilt axis, rotates each coherent group -as one body about its pivot, recomputes the overlap of the displaced geometry and reports the -relative error of the linear prediction per row. Include the largest scale a tolerance scan -will use, since that is the regime it relies on; the cost is one Biot–Savart pass per row. +[ErrorFields.correctability] +uncorrectable_coils = ["EFCC_U"] # coil sets the error-field correction cannot reduce +efc_factor = 2.0 # divisor applied to correctable contributions under correction -```julia -lin = EF.linearity_check("gpec.h5"; scales=(1.0, 2.0, 8.0)) -maximum(filter(!isnan, lin.relative_error)) # the worst row -AEF.plot_linearity_check(lin) # stems per coil, log axis -AEF.plot_linearity_residuals("gpec.h5"; at=:step, yscale=:log10) # the stored step residuals, readable below 1 % +[ErrorFields.other_field] +magnitude = 8.2e-6 # overlap budget of sources not attributed to any coil +sigma = 0.0 +radial_shape = "ring" ``` -A coil set whose field at the run's toroidal mode is round-off (an axisymmetric hoop at n = 1) has -no resonant fraction: `coil_overlaps`, `combine_overlaps` and `plot_coil_sensitivities(quantity=:fraction)` -report `NaN` for it rather than the ratio of two noise vectors, which would read as a fraction -between 1 % and 30 %; the floor is `1e-8` of the largest field norm among the coil sets judged -together (`resonant_fraction_percent`). - -## Comparing coil revisions +A coherent group's tilt is a rigid rotation of its members about the pivot, so a member at height ``z`` also shifts by ``(z - z_\mathrm{pivot})\,\theta``; coil sets without a tolerance block contribute their as-designed overlap only; a coil's `current_factor` scales its overlap and sensitivities together, the knob for a risk against the current in an uncorrectable coil. -Assessing a *new* coil design against an existing run costs one Biot-Savart pass per coil set and -no plasma solve: the equilibrium, the control surface and the resonant coupling all come from the -stored file. Gather them once with `ResonantDriveContext` and hand it whatever geometry you like. +From a script the stage is one call on the plasma response and the coil sets, and every product is re-evaluated post hoc from the file for any surface window, mode, fit, tolerance scale or coil subset: ```julia -ctx = EF.ResonantDriveContext("gpec.h5") # nzeta_coil=…, dat_dir=… override the deck -old = EF.coil_overlaps(ctx, ForcingTerms.load_coil_sets(old_cfg, 1)) -new = EF.coil_overlaps(ctx, ForcingTerms.load_coil_sets(new_cfg, 1)) +ef = error_fields(resp, pe.coil_sets, ForcingTerms.CoilConfig(ft_ctrl); # ft_ctrl: the [ForcingTerms] control + tolerances=EF.read_tolerance_toml("tolerances.toml"), + scenario=EF.ScenarioParameters(ffs.equil; n_e=5.0)) +save_to_hdf5("gpec.h5", ef) -new[1].delta # dimensionless overlap δ = Vᴴb̃ / B_T0 -new[1].resonant_fraction_percent # how much of this coil's own spectrum is resonant -new[1].spectrum # b̃ itself, for the diagnostics below +table = EF.sensitivity_table("gpec.h5"; psi_low=0.5) # overlaps and sensitivities per coil set, edge window +mc = EF.run_monte_carlo("gpec.h5"; psi_low=0.5, tolerance_scale=2.0) # |δ| distribution, intrinsic and corrected +risk = EF.locking_risk("gpec.h5"; n_e=5.0, risk_ctrl=EF.RiskControl(; year=2026, fit="OLS")) +EF.allowable_tolerance(EF.ToleranceScan("gpec.h5"), 1.0; corrected=true) # tolerance multiplier at 1% locking risk +EF.max_correctable_overlap(EF.read_efc_couplings("gpec.h5")[1]; delta_threshold=1.4e-4, torque_budget=4.0) +AEF.plot_error_field_summary("gpec.h5"; save_path="error_field_summary.png") ``` -Every normalization the quantity is quoted in travels on the `CoilOverlap`, so a caller -never has to work out which one a bare number was in. +`correction_requirement` gives the current each array needs, by the dominant mode or jointly over every rational surface in the least-squares sense, and `uncorrectable_probability` the fraction of the Monte Carlo overlap beyond what the NTV-limited array can correct; `phasing_map` evaluates the dominant-mode field per kilo-ampere-turn over the relative phases of several arrays; `ResonantDriveContext` and `coil_overlaps` assess a new coil design against the stored run at one Biot–Savart pass per coil set, with `combine_overlaps` applying a design current pattern afterwards. +`Analysis.ErrorFields` plots each product, and every plotting function takes `label => path` pairs so coil-design revisions overplot on one axis. -A revision can rename or split coils. Because the field is linear in the currents, a design current -pattern is applied afterwards rather than by re-running the geometry, and a split is written as the -combination it is: +## Outputs -```julia -sets = ForcingTerms.conductors(ForcingTerms.read_coil_dat("old_pair.dat")) # one file, two coils -ovs = EF.coil_overlaps(ctx, sets) # each at 1 kA, say +| Dataset | Meaning | +|---|---| +| `ErrorFields/CoilSensitivities/field_as_designed`, `shift_sensitivity_per_m`, `tilt_sensitivity_per_deg` | the control-surface spectrum of each coil set as built and its derivatives with respect to the rigid shifts and tilts, with the current pattern they were evaluated at and the finite-difference linearity residual per tap | +| `ErrorFields/CoilSensitivities/DominantMode/delta_as_designed`, … | the overlaps and sensitivities projected onto the run's full-window dominant mode, their direction-averaged in-plane magnitudes, and the shift and tilt that would cancel the as-designed overlap | +| `ErrorFields/MonteCarlo/abs_delta_pdf`, `abs_delta_efc_pdf`, `abs_delta_bin_edges` | the intrinsic and corrected ``\lvert\delta\rvert`` histograms over the sampled tolerances, per batch and averaged | +| `ErrorFields/MonteCarlo/abs_delta_total_as_designed`, `abs_delta_worst_case` | the overlap of the design as built and the worst-case bound with every term at its tolerance edge and in phase | +| `ErrorFields/Risk/locking_probability_given_delta`, threshold density | ``P(\mathrm{lock}\,|\,\delta)`` and the threshold distribution on the Monte Carlo grid | +| `ErrorFields/Risk/locking_probability_percent`, `locking_probability_efc_percent`, `*_batches_percent`, `locking_probability_as_designed_percent`, `locking_probability_fit_threshold_percent` | the locking probability of the intrinsic and corrected distributions with per-batch values, of the design as built, and at the sharp fitted threshold | +| `ErrorFields/Risk/ToleranceScan/` | the same against `tolerance_scale`; the allowable tolerance at a target risk is its inversion | +| `ErrorFields/NTV/` | each correction array's dominant-mode coupling ``C_c``, resonant fraction and NTV torque per kAt² for its full field and with the dominant mode projected out | -EF.combine_overlaps(ovs, "old_pair_1" => 20.0, "old_pair_2" => -15.0) # weights in kA -``` +Every dataset name follows one grammar, `[_][_efc][_][_]`: `delta` is the overlap and `abs_delta` its magnitude; `_as_designed`, `_sampled`, `_worst_case` name the geometry; `_efc` marks error-field correction applied; a probability stored as a percentage carries `_percent`. -Weights multiply each spectrum as it was evaluated; they are not absolute currents, since the sets -may have been swept at any current. Evaluating every conductor at 1 kA and weighting by the design -current in kA is the idiom that makes them read as currents. An unmatched name raises rather than -contributing zero, so a coil renamed between revisions cannot quietly drop out of a comparison. +## Numerical notes -For sensitivities to rigid motion as well, `compute_coil_sensitivities` takes the same context. -It costs thirteen spectrum evaluations per coil set instead of one, so reach for `coil_overlaps` -when only the overlaps are wanted. - -## Result names - -Every result struct in `ErrorFields` (`CoilSensitivities`, `SensitivityTable`, `CoilOverlap`, -`MonteCarloResult`, `RiskResult`, `ToleranceScan`, `PhasingMap`, `EFCCoupling`) names its fields by one -grammar, and every field name is also its HDF5 dataset name under `ErrorFields/`: - -``` -[_][_efc][_][_] -``` +- The sensitivities assume the overlap is linear in the misalignments out to the tolerance edge; the run checks it only at the finite-difference step, and `linearity_check` displaces each coil set by `±scale × tolerance` and recomputes its overlap to report the relative error of the linear model at the largest scale a tolerance scan will use (one Biot–Savart pass per row). +- Monte Carlo batches are seeded individually, so results are bit-identical for any thread count; `risk_convergence` repeats the Monte Carlo over sample and bin counts, and the batch spread must fall as ``1/\sqrt{n_\mathrm{sample}}`` while the value must not move with the bin count by more than the spread. +- A threshold scaling is fitted for one toroidal mode number, so the risk refuses a run that spans several; project onto one n's coupling first. +- A coil set whose field at the run's toroidal mode is round-off (an axisymmetric hoop at n = 1) has no resonant fraction and reports `NaN` rather than the ratio of two noise vectors. +- The NTV torques are stored with their sign and the limits consume the budget with their magnitude; the torque budget is the plasma's intrinsic torque and an analysis choice. -- **quantity** is the physical thing, spelled out: `delta` is the dimensionless dominant-mode overlap - `v·b̃ / B_T0`, a complex phasor per coil set, and `abs_delta` its magnitude; `threshold` is the - penetration threshold in the same units; `locking_probability` is the probability of locking (never - abbreviated, and never called a risk in a name); `resonant_fraction` is `100·|v·b̃| / ‖b̃‖`; - `field`, `spectrum`, `torque`, `current`, `shift_sensitivity` and `tilt_sensitivity` are what they say. -- **instance** says which geometry or value the quantity is evaluated for: `_as_designed`, the - unperturbed design; `_sampled`, over the Monte Carlo misalignment samples; `_worst_case`, every term - at its tolerance edge and in phase; `_fit`, a scaling-law value at the fitted exponents without the - fit's scatter; `_total`, the coherent sum over coil sets (a per-coil vector carries no qualifier and - documents its shape as `[ncoil_set]`). A result evaluated on whatever geometry the caller passed - (`CoilOverlap`, `EFCCoupling`) carries no instance qualifier. The word *nominal* is not used in names. -- **`_efc`** marks a value with error-field correction applied, as in the TOML keys `efc_factor` and - `efc_coils`. -- **statistic**: `_mean`, `_pdf`, `_batches`, `_spread` (range over the Monte Carlo batches), - `_samples`, `_bin_edges`. -- **unit** is present whenever the quantity has one: `_percent`, `_t` (tesla), `_m`, `_mm`, `_deg`, - `_per_m`, `_per_deg`, `_per_mm`, `_per_kat`, `_per_kat2`. `delta` and probabilities given as fractions - carry none; a probability stored as a percentage always carries `_percent`. - -The order is fixed, so `locking_probability_efc_batches_percent` reads as "the locking probability, -corrected, per batch, in percent", and `abs_delta_total_as_designed` as "the magnitude of the total -overlap of the design as built". - -## API Reference +## API reference ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.ErrorFields] diff --git a/docs/src/figures/ballooning/mercier_comparison.png b/docs/src/figures/ballooning/mercier_comparison.png deleted file mode 100644 index 2c8da8b6cd..0000000000 Binary files a/docs/src/figures/ballooning/mercier_comparison.png and /dev/null differ diff --git a/docs/src/figures/inner_layer/convergence_Sinvariance.png b/docs/src/figures/inner_layer/convergence_Sinvariance.png deleted file mode 100644 index cedbc59056..0000000000 Binary files a/docs/src/figures/inner_layer/convergence_Sinvariance.png and /dev/null differ diff --git a/docs/src/figures/inner_layer/make_convergence_Sinvariance.jl b/docs/src/figures/inner_layer/make_convergence_Sinvariance.jl deleted file mode 100644 index 2b4129baf5..0000000000 --- a/docs/src/figures/inner_layer/make_convergence_Sinvariance.jl +++ /dev/null @@ -1,45 +0,0 @@ -# make_convergence_Sinvariance.jl -# -# Figure: honest numerical error bar on the matching data Δ at Q = 500i. The -# `delta_convergence` battery re-solves with each numerical knob perturbed on -# an independent axis (contour angle θ, mesh resolution, series order/radius, -# march tolerance, handoff radius, purification e-folds). Δ -# is an analytic invariant of the contour, so the worst-case spread across -# these orthogonal perturbations bounds the true error. Seven solves; ~30 s. -# -# Run manually: julia --project=. docs/src/figures/inner_layer/make_convergence_Sinvariance.jl - -include(joinpath(@__DIR__, "..", "..", "..", "figure_tools.jl")) - -using GeneralizedPerturbedEquilibrium -using Printf - -const IL = GeneralizedPerturbedEquilibrium.InnerLayer - -p = IL.GGJ.q4_surface_benchmark() -Q = 500.0im - -conv = IL.delta_convergence(p, Q; verbose=false) -names = [r.name for r in conv.table] -d1 = [max(r.d1, 1e-16) for r in conv.table] # relative change of Δ₊ (interchange) -d2 = [max(r.d2, 1e-16) for r in conv.table] # relative change of Δ₋ (tearing) -n = length(names) - -plt = plot(; size=(820, 540), legend=:topright, yscale=:log10, - xticks=(1:n, names), xrotation=30, ylims=(1e-12, 1e-2), - ylabel="relative change of Δ per knob", xlabel="perturbed numerical knob", - title=@sprintf("Convergence & contour-invariance of Δ, Q = %gi (q=4)", imag(Q)), - left_margin=10Plots.mm, bottom_margin=14Plots.mm) - -scatter!(plt, 1:n, d1; marker=:circle, ms=7, color=1, label="δΔ₊ (interchange)") -scatter!(plt, 1:n, d2; marker=:diamond, ms=7, color=2, label="δΔ₋ (tearing)") - -# Worst-case spread — the reported error bar on each parity. -hline!(plt, [conv.spread[1]]; color=1, ls=:dash, lw=1.5, - label=@sprintf("worst-case Δ₊ %.1e", conv.spread[1])) -hline!(plt, [conv.spread[2]]; color=2, ls=:dash, lw=1.5, - label=@sprintf("worst-case Δ₋ %.1e", conv.spread[2])) - -save_doc_figure(plt, "inner_layer", "convergence_Sinvariance"; - script=basename(@__FILE__), - depends=["src/InnerLayer/GGJ/Ray.jl"]) diff --git a/docs/src/figures/inner_layer/make_rotated_ray_contour.jl b/docs/src/figures/inner_layer/make_rotated_ray_contour.jl deleted file mode 100644 index a23c39fba0..0000000000 --- a/docs/src/figures/inner_layer/make_rotated_ray_contour.jl +++ /dev/null @@ -1,65 +0,0 @@ -# make_rotated_ray_contour.jl -# -# Figure: the rotated integration ray in the complex layer-coordinate (x) plane. -# Shows why the :ray backend rotates the contour by θ = arg(Q)/4: the real-axis -# contour (θ = 0) runs straight through the pseudo-resonance at -# x² = −Q²(G + K F) (real and large on the imaginary-Q axis), while the rotated -# ray clears it. Geometry only — no BVP solve — so it is cheap to regenerate. -# -# Run manually: julia --project=. docs/src/figures/inner_layer/make_rotated_ray_contour.jl - -include(joinpath(@__DIR__, "..", "..", "..", "figure_tools.jl")) - -using GeneralizedPerturbedEquilibrium -using Printf - -const IL = GeneralizedPerturbedEquilibrium.InnerLayer -const GGJ = IL.GGJ - -# q = 4 rational-surface benchmark on the imaginary-Q axis (the regime that -# defeats :galerkin). Q is the scaled growth rate; θ makes the parabolic -# exponent real. -p = IL.GGJ.q4_surface_benchmark() -Q = 500.0im -θ = angle(Q) / 4 # = π/8 = 22.5° - -# Pseudo-resonance location: x² = −Q²(G + K F). On the imaginary-Q axis this is -# real and large, so it sits ON the real-x contour. -x_pr = sqrt(-Q^2 * (p.G + p.K * p.F)) -xpr = real(x_pr) # imaginary part ≈ 0 here - -# Matching radius from the series-residual criterion along the ray (annotation). -S, _, ok = IL.pick_smax(p, Q; θ=θ) - -# Draw both contours out to ~1.8× the pseudo-resonance radius so it is in frame. -smax = 1.8 * xpr -s = range(0, smax; length=400) -ray = cis(θ) .* s # x = e^{iθ} s - -plt = plot(; size=(720, 560), legend=:topleft, framestyle=:zerolines, - xlabel="Re x", ylabel="Im x", aspect_ratio=:equal, - title="Rotated integration ray, Q = 500i (q=4 surface)", - left_margin=6Plots.mm, bottom_margin=4Plots.mm) - -# Real-axis contour (θ = 0) and the rotated ray. -plot!(plt, [0, smax], [0, 0]; lw=2, ls=:dash, color=:gray, - label="real-axis contour (θ = 0)") -plot!(plt, real.(ray), imag.(ray); lw=3, color=1, - label=@sprintf("rotated ray x = e^{iθ}s, θ = %.1f°", rad2deg(θ))) - -# Pseudo-resonance on the real axis — the point the rotation steps around. -scatter!(plt, [xpr], [0.0]; marker=:xcross, ms=9, color=:red, msw=3, - label="pseudo-resonance x² = −Q²(G+KF)") - -# Clearance note above the ray near the pseudo-resonance radius, plus the true -# matching radius. Positions are tied to xpr so they stay inside the frame. -annotate!(plt, 1.02 * xpr, 0.56 * xpr, - Plots.text("ray stays clear of the\nreal-axis pseudo-resonance", 8, :left, - RGB(0.15, 0.15, 0.15))) -annotate!(plt, 0.45 * xpr, -0.14 * xpr, - Plots.text(@sprintf("matching radius S ≈ %.2g%s", S, ok ? "" : " (series tol not met)"), - 8, :gray, :left)) - -save_doc_figure(plt, "inner_layer", "rotated_ray_contour"; - script=basename(@__FILE__), - depends=["src/InnerLayer/GGJ/Ray.jl"]) diff --git a/docs/src/figures/inner_layer/make_solution_profiles.jl b/docs/src/figures/inner_layer/make_solution_profiles.jl deleted file mode 100644 index 6f89b9c66b..0000000000 --- a/docs/src/figures/inner_layer/make_solution_profiles.jl +++ /dev/null @@ -1,59 +0,0 @@ -# make_solution_profiles.jl -# -# Figure: the inner-layer fields (Ψ, Ξ, Υ) reconstructed on the rotated ray for -# the q=4 benchmark surface. The Hermite BVP solution on [0, s_m] (solid) and -# the analytic u_small + Δ·u_big asymptotic representation for s ≥ S (dashed) -# share the same power-law tail — the numeric↔asymptotic overlap the outer- -# region matching relies on. One BVP solve; a few minutes at large imaginary Q. -# -# Run manually: julia --project=. docs/src/figures/inner_layer/make_solution_profiles.jl - -include(joinpath(@__DIR__, "..", "..", "..", "figure_tools.jl")) - -using GeneralizedPerturbedEquilibrium -using Printf - -const IL = GeneralizedPerturbedEquilibrium.InnerLayer -const GGJ = IL.GGJ - -p = IL.GGJ.q4_surface_benchmark() -Q = 2.0im # imaginary axis, beyond the |Q|≪1 shooting regime, -# small enough that the BVP domain reaches -# the series radius directly (no march) — so the -# numeric and asymptotic segments join seamlessly. -isol = 2 # odd-parity (tearing, Δ₋) solution: Ψ'(0)=Ξ(0)=Υ(0)=0 - -res = IL.solve_ray(p, Q) -s_m = res.breaks[end] -@printf("solve_ray: Q=%s θ=%.3f S=%.4g s_m=%.4g cells=%d\n", - string(Q), res.θ, res.S, s_m, length(res.breaks) - 1) - -# BVP solution over [0, s_m]. -prof = IL.solution_profile(res; npc=8) - -# Analytic tail on [S, few·S] where the inps series is trusted. -srange = 10 .^ range(log10(res.S), log10(4 * res.S); length=80) -asy = IL.asymptotic_profile(p, res, srange) - -fields = ((:Ψ, prof.Ψ, asy.Ψ, 1), (:Ξ, prof.Ξ, asy.Ξ, 2), (:Υ, prof.Υ, asy.Υ, 3)) - -plt = plot(; size=(760, 560), legend=:bottomleft, xscale=:log10, yscale=:log10, - ylims=(1e-3, 1e2), - xlabel="ray parameter s (x = e^{iθ}s)", ylabel="|field| (odd parity, tearing)", - title=@sprintf("Inner-layer fields on the rotated ray, Q = %gi (q=4)", imag(Q)), - left_margin=8Plots.mm, bottom_margin=5Plots.mm) - -for (nm, col, acol, ci) in fields - plot!(plt, prof.s[2:end], abs.(col[2:end, isol]); lw=2.5, color=ci, - label="|$(nm)| BVP") - plot!(plt, asy.s, abs.(acol[:, isol]); lw=2, ls=:dash, color=ci, label="|$(nm)| asymptotic") -end - -# Matching radius: with no march the BVP edge s_m coincides with the series -# radius S, so numeric and asymptotic meet at a single point. -vline!(plt, [res.S]; color=:black, ls=:dot, lw=1, - label=@sprintf("S = s_m ≈ %.1f (match point)", res.S)) - -save_doc_figure(plt, "inner_layer", "solution_profiles"; - script=basename(@__FILE__), - depends=["src/InnerLayer/GGJ/Ray.jl", "src/InnerLayer/GGJ/InnerAsymptotics.jl"]) diff --git a/docs/src/figures/inner_layer/rotated_ray_contour.png b/docs/src/figures/inner_layer/rotated_ray_contour.png deleted file mode 100644 index 27884cf1f6..0000000000 Binary files a/docs/src/figures/inner_layer/rotated_ray_contour.png and /dev/null differ diff --git a/docs/src/figures/inner_layer/solution_profiles.png b/docs/src/figures/inner_layer/solution_profiles.png deleted file mode 100644 index de394ef0d6..0000000000 Binary files a/docs/src/figures/inner_layer/solution_profiles.png and /dev/null differ diff --git a/docs/src/figures/manifest.toml b/docs/src/figures/manifest.toml index 731585bc5d..fd45fad1e4 100644 --- a/docs/src/figures/manifest.toml +++ b/docs/src/figures/manifest.toml @@ -7,37 +7,12 @@ [ballooning.delta_prime_matching] commit = "legacy" date = "unknown" -depends = ["src/ForceFreeStates/Ballooning.jl"] +depends = ["src/LocalStability/Ballooning.jl"] note = "Predates the provenance system; compared methods for extracting the ballooning Delta_prime and c_a1. Original generator script not preserved." script = "(not preserved)" -[ballooning.mercier_comparison] -commit = "legacy" -date = "unknown" -depends = ["src/ForceFreeStates/Ballooning.jl"] -note = "Predates the provenance system; compared the since-removed standalone Mercier.jl against Bal.jl. Not reproducible from current code (Mercier.jl was removed once Bal.jl subsumed it)." -script = "(not preserved)" - [inner_layer.backend_regime_map] commit = "3a0837e3-dirty" date = "2026-07-08" depends = ["src/InnerLayer/GGJ/Ray.jl", "src/InnerLayer/GGJ/Galerkin.jl", "src/InnerLayer/GGJ/Shooting.jl"] script = "make_backend_regime_map.jl" - -[inner_layer.convergence_Sinvariance] -commit = "5d11aad3e" -date = "2026-10-21" -depends = ["src/InnerLayer/GGJ/Ray.jl"] -script = "make_convergence_Sinvariance.jl" - -[inner_layer.rotated_ray_contour] -commit = "3a0837e3-dirty" -date = "2026-07-08" -depends = ["src/InnerLayer/GGJ/Ray.jl"] -script = "make_rotated_ray_contour.jl" - -[inner_layer.solution_profiles] -commit = "5d11aad3e" -date = "2026-10-21" -depends = ["src/InnerLayer/GGJ/Ray.jl", "src/InnerLayer/GGJ/InnerAsymptotics.jl"] -script = "make_solution_profiles.jl" diff --git a/docs/src/forcing_terms.md b/docs/src/forcing_terms.md index 256fbfebc4..8ad9ae5a0f 100644 --- a/docs/src/forcing_terms.md +++ b/docs/src/forcing_terms.md @@ -1,79 +1,31 @@ -# ForcingTerms Module +# Forcing Terms -The ForcingTerms module specifies external magnetic field perturbations for the perturbed -equilibrium calculation. It supports reading from ASCII or HDF5 files, as well as computing -perturbations directly from coil Biot-Savart calculations. +`ForcingTerms` specifies the external magnetic perturbation that drives the perturbed equilibrium: the spectrum of the normal field on the control surface at `psilim`, given either as a file of (m, n) amplitudes or computed by Biot–Savart from coil geometry. +The `[ForcingTerms]` section is read when `[PerturbedEquilibrium]` runs; its result is the `RMPField` that `perturbed_equilibrium` applies to the plasma response. -## Normalization conventions +## Normalization -Two normalization tags are supported, specified inside the file (not in `gpec.toml`): +Internally every forcing is the unit-norm flux spectrum ``\Phi_x``, the Fourier amplitude of ``2\pi R\,(B_R\,\partial Z/\partial\theta - B_Z\,\partial R/\partial\theta)`` over the unit-norm angles ``\theta, \zeta \in [0, 1]``, the quantity the Fortran GPEC calls `Phi_x`. +A forcing file declares its own normalization by a tag; a coil set needs none, because the Biot–Savart pipeline produces ``\Phi_x`` directly. -| Tag | Quantity | Units | Notes | -|-----|----------|-------|-------| -| `"normal_field_T"` | Fourier modes of B·n̂ | Tesla | **Default** — most intuitive for users | -| `"sfl_flux_Wb"` | Fourier modes of R×(B_R ∂Z/∂θ − B_Z ∂R/∂θ) | T·m² | 2π-angle SFL flux (user input) | +| Tag | Quantity | Units | On load | +|---|---|---|---| +| `"normal_field_T"` (default) | Fourier modes of ``B\cdot\hat n`` over ``\theta \in [0, 2\pi]`` | T | weighted by the area element ``J\lvert\nabla\psi\rvert`` and projected onto every harmonic of the run; a poloidally varying weight mixes modes, so the file must cover its sidebands (a warning asks for a wider m range when an edge harmonic carries more than 10% of the flux) | +| `"sfl_flux_Wb"` | the straight-field-line flux ``R\,(B_R\,\partial Z/\partial\theta - B_Z\,\partial R/\partial\theta)`` over ``\theta \in [0, 2\pi]`` | T·m² | multiplied by ``(2\pi)^2`` | +| `"unit_flux"` | ``\Phi_x`` itself, as a run records its forcing for replay | T·m² | loaded as is | -### Angle convention +## File formats -Julia's ForcingTerms use the **unit-norm angle convention** (θ_norm ∈ [0, 1], ζ_norm ∈ [0, 1]) -internally, matching Fortran GPEC and the Equilibrium/ForceFreeStates modules. - -The coil Biot-Savart pipeline produces mode amplitudes directly equal to Fortran's `Phi_x`: - -``` -Julia output ≡ Fortran Phi_x -``` - -The formula is `Phi_x(θ_norm, ζ_norm) = 2π × R × (B_R ∂Z/∂θ_norm − B_Z ∂R/∂θ_norm)`, where -the `2π` comes from the toroidal Jacobian ∂r/∂ζ_norm = 2π·R·ê_φ. - -Users who provide forcing files may use the more intuitive **2π-angle Fourier amplitudes** -(e.g. the amplitude of cos(mθ − nζ) where θ ∈ [0, 2π]). These are tagged `"sfl_flux_Wb"` -and multiplied by (2π)² internally to reach unit-norm convention. - -### Automatic conversion - -Files tagged `normal_field_T` are automatically converted to unit-norm Phi_x convention -during loading by weighting with the control-surface area element J|∇ψ|. This weight varies -poloidally, so the conversion is a mode-mixing operation and the flux is projected onto -every harmonic of the run (`mlow:mhigh`), sidebands included. - -The file's spectrum must also extend far enough that the harmonics it omits are negligible, -since those would mix into the run's range too. GPEC estimates this from the file's edge -harmonics: if either one alone carries more than 10% of the converted flux, a warning asks for -a wider m range. Files with fewer than three harmonics are not checked. - -Files tagged `sfl_flux_Wb` (2π-angle SFL flux) are scaled by (2π)² to reach unit-norm. - -The coil Biot-Savart code (`forcing_data_format = "coil"`) always produces unit-norm -amplitudes directly — no conversion is applied. - -## ASCII file format +ASCII (`forcing_data_format = "ascii"`): an optional `# normalization: ` line, comment lines starting with `#`, and one `n m real(amplitude) imag(amplitude)` row per mode; all four columns are required: ``` # normalization: normal_field_T -# Format: n m real(amplitude) imag(amplitude) 1 6 1.0e-4 0.0 1 7 5.0e-5 2.0e-5 ``` -- The `# normalization: ` line is optional; default is `normal_field_T`. -- Lines starting with `#` are comments and are ignored by the data parser. -- Columns: toroidal mode n, poloidal mode m, real part of amplitude, imaginary part (optional, default 0). -- Amplitudes are in Tesla (for `normal_field_T`) or T·m² (for `sfl_flux_Wb`). - -## HDF5 file format - -Required datasets: -- `n`: integer array of toroidal mode numbers -- `m`: integer array of poloidal mode numbers -- `amplitude_real`: float array of real parts -- `amplitude_imag`: float array of imaginary parts (optional) - -Optional root attribute: -- `normalization`: string, either `"normal_field_T"` (default) or `"sfl_flux_Wb"` +HDF5 (`forcing_data_format = "hdf5"`): datasets `n`, `m`, `amplitude_real` and `amplitude_imag`, all required, with the optional root attribute `normalization`: -Example (Julia): ```julia using HDF5 h5open("forcing.h5", "w") do f @@ -85,57 +37,33 @@ h5open("forcing.h5", "w") do f end ``` -## Coil Biot-Savart format - -Set `forcing_data_format = "coil"` in `[ForcingTerms]` TOML to compute the field -directly from coil geometry. Output is always in `sfl_flux_Wb` convention. - -The coil pipeline: -1. `sample_boundary_grid`: evaluate plasma boundary geometry at (mtheta × nzeta) grid -2. `compute_biot_savart_boundary!`: compute B at all grid points via Biot-Savart law -3. `project_normal_flux!`: compute Phi_x = 2π×R×(B_R ∂Z/∂θ_norm − B_Z ∂R/∂θ_norm) -4. `fourier_decompose_bn`: Fourier decompose to get bmn in unit-norm (= Phi_x) convention - -Steps 2–4 are linear in the coils, so a `CoilForcingGrid` (step 1 plus the observation-point -layout, built once per equilibrium and `n`) can be reused to evaluate coil sets one at a time -with `coil_forcing_modes`; the per-set spectra sum to the assembly's. Use this when many coil -variants must be compared against one plasma — moved or re-wound coils, or one coil at a time — -without re-sampling the boundary: - -```julia -forcing_grid = CoilForcingGrid(equil, cfg, n; psi=ffs.psilim) -spectra = [coil_forcing_modes(cs, forcing_grid, n, m_low, m_high) for cs in coil_sets] -``` - -The toroidal angle for each observation point at SFL grid coordinate (θ_i, ζ_j) is: -``` -φ_phys(i, j) = -helicity × (2π × ζ_j + ν(ψ, θ_i)) -``` -where `helicity = sign(Bt) × sign(Ip)` and `ν(ψ, θ)` is the SFL toroidal coordinate correction -stored in `equil.rzphi_nu`. This matches Fortran's `phi = -helicity*(twopi*czeta + crzphi_f(3))`. -For DIII-D (Bt < 0, Ip > 0 → helicity = -1), ν ranges from about -3.4 to +3.4 radians at the -boundary — omitting it completely scrambles the Fourier mode decomposition. - -### Coil geometry sources +## Coils -Each `[[ForcingTerms.coil_set]]` block supplies one coil set from one of several sources, -selected by `source`. Shifts and tilts (`shiftx/y/z`, `tiltx/y/z`, `n_tilt`, ...) apply to -every source, so analytic coils can be perturbed exactly like file-based ones. +`forcing_data_format = "coil"` computes the field from coil geometry: the plasma boundary is sampled on an `mtheta_coil × nzeta_coil` grid, the Biot–Savart field of every coil set is evaluated there, projected onto the normal flux and Fourier decomposed. +The toroidal angle of each observation point includes the straight-field-line offset ``\nu(\psi, \theta)`` and the helicity ``\mathrm{sign}(B_T)\,\mathrm{sign}(I_p)`` of the equilibrium, without which the decomposition would be wrong. +Each `[[ForcingTerms.coil_set]]` block supplies one coil set, selected by `source`; shifts and tilts (`shiftx/y/z`, `tiltx/y/z`, `n_tilt`, …) apply to every source, which is what the [Error Fields](error_fields.md) stage sweeps. +From `examples/DIIID-like_ideal_example/gpec.toml` (the first block) and the other sources: ```toml -# Legacy ASCII .dat file (unchanged; resolves to {dat_dir}/{machine}_{name}.dat) +[ForcingTerms] +forcing_data_format = "coil" +machine = "d3d" +mtheta_coil = 480 +nzeta_coil = 32 + +# ASCII .dat file, resolved to {dat_dir}/{machine}_{name}.dat [[ForcingTerms.coil_set]] -name = "il" -currents = [1e4, 1e4, 1e4, 1e4, 1e4, 1e4] +name = "c" +currents = [1000.0, 500.0, -500.0, -1000.0, -500.0, 500.0] -# Modern HDF5 file; `name` selects one set from a multi-set file (omit if it holds one) +# HDF5 file; `name` selects one set of a multi-set file [[ForcingTerms.coil_set]] source = "file" h5_file = "/path/to/my_coils.h5" name = "il" currents = [1e4, 1e4, 1e4, 1e4, 1e4, 1e4] -# Analytic PF hoop — a horizontal circle; tilt breaks axisymmetry to give n ≠ 0 forcing +# Analytic PF hoop; a tilt breaks axisymmetry to give n ≠ 0 forcing [[ForcingTerms.coil_set]] source = "pf_hoop" radius = 1.8 @@ -143,7 +71,7 @@ height = 0.9 currents = [5e3] tiltx = [2.0] -# Analytic window-pane (picture-frame) array — corner mode: two opposite [R, Z] corners +# Analytic window-pane array, by two opposite [R, Z] corners ... [[ForcingTerms.coil_set]] source = "window_pane" n_coils = 8 @@ -151,10 +79,7 @@ rz_corners = [[2.2, -0.6], [2.2, 0.6]] gap_fraction = 0.15 currents = [1e3, -1e3, 1e3, -1e3, 1e3, -1e3, 1e3, -1e3] -# Analytic window-pane — standoff mode: placed by physical distance off the plasma surface. -# The frame center sits `standoff` m along the outward surface normal at cylindrical -# `poloidal_angle` (degrees); legs span `poloidal_length` m, tilted by `poloidal_tilt` -# (degrees, 0 = tangential). Requires an equilibrium (the main forcing pipeline supplies it). +# ... or by standoff off the plasma surface at a poloidal angle (degrees) [[ForcingTerms.coil_set]] source = "window_pane" n_coils = 2 @@ -165,8 +90,7 @@ poloidal_tilt = 0.0 gap_fraction = 0.15 currents = [1e4, -1e4] -# Analytic helical array on a circular cross-section torus (pitch = m/n). -# theta_lo/theta_hi (poloidal extent, degrees) default to a full 0–360 turn. +# Analytic helical array on a circular torus, pitch m_hel/n_hel [[ForcingTerms.coil_set]] source = "helical" R0 = 1.7 @@ -177,47 +101,21 @@ n_coils = 4 currents = [1e3, 1e3, 1e3, 1e3] ``` -### HDF5 coil geometry format +A coil HDF5 file stores one subgroup per coil set under a parent group (default `coils`): `//x`, `y`, `z` as `ncoil × s × nsec` Cartesian metres, an optional `currents` dataset, and the winding multiplier as the attribute `nw`; `convert_coil_dat_to_h5` and `convert_coil_h5_to_dat` translate the ASCII format, `read_coil_h5` and `write_coil_h5` read and write it. -A coil HDF5 file stores one **subgroup per coil set** under a parent group (default `coils`), -so a single file can hold many sets and the reader ignores any unrelated content elsewhere: +The field is linear in the coils, so a `CoilForcingGrid` built once per equilibrium and n evaluates coil sets one at a time, for comparing coil variants against one plasma without re-sampling the boundary: -```text -//x Float64 ncoil×s×nsec Cartesian metres -//y Float64 ncoil×s×nsec -//z Float64 ncoil×s×nsec -//currents Float64 ncoil optional -/ attr nw Float64 winding multiplier, default 1.0 +```julia +forcing_grid = CoilForcingGrid(equil, cfg, n; psi=ffs.psilim) +spectra = [coil_forcing_modes(cs, forcing_grid, n, m_low, m_high) for cs in coil_sets] ``` -Use `convert_coil_dat_to_h5` / `convert_coil_h5_to_dat` to migrate legacy files, and -`read_coil_h5` / `write_coil_h5` for direct I/O. - -### Rerun snapshot - -When a coil run writes `gpec.h5`, the coil geometry actually used (after shifts/tilts and with -currents) is captured under `Input/RawInputs/Coils/`. The run can then be replayed with -`main_from_h5` even if the original `.dat`/`.h5` files are gone: - -- `--coil-source forcing-modes` (default) recomputes the field from the stored TOML coil config - (needs the original geometry files to be present). -- `--coil-source coils` recomputes the field from the geometry stored in `gpec.h5`, so no - original files are required — and the equilibrium may be overridden on replay. - -An ASCII or HDF5 forcing file is captured raw, with its `normalization` tag, under -`Input/RawInputs/ForcingTerms/`. The replay writes it back out as `_forcing.h5` beside its -output and reads it like any forcing file, so the normalization is redone against the replay's -own equilibrium. - -## Comparison with Fortran GPEC +## Outputs -| Quantity | Julia convention | Fortran GPEC convention | -|----------|-----------------|------------------------| -| Angle range | θ_norm ∈ [0, 1], ζ_norm ∈ [0, 1] | θ_norm ∈ [0, 1], ζ_norm ∈ [0, 1] | -| Coil forcing output | unit-norm = 2π×R×(B_R ∂Z/∂θ_norm − B_Z ∂R/∂θ_norm) | `Phi_x` [Wb] | -| Ratio | Julia output ≡ Fortran Phi_x (directly comparable) | — | +The applied spectrum is written with the response: `PerturbedEquilibrium/ForcingModes/` (`m`, `n`, `amplitude`) and the control-surface spectra `PerturbedEquilibrium/forcing_b`, `forcing_b_root_area` and `forcing_b_area` (see [Perturbed Equilibrium](perturbed_equilibrium.md)). +A forcing file is captured raw with its tag under `Input/RawInputs/ForcingTerms/`, and the coil geometry as used, after shifts and tilts, under `Input/RawInputs/Coils/`, so that `main(["gpec.h5"])` can replay the run without the original files; see the [Scripting API](api.md). -## API Reference +## API reference ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.ForcingTerms] diff --git a/docs/src/galerkin.md b/docs/src/galerkin.md index 05ecabf956..b3e59a745b 100644 --- a/docs/src/galerkin.md +++ b/docs/src/galerkin.md @@ -1,57 +1,57 @@ -# Galerkin Δ′ Solver (RDCON outer region) - -The Galerkin solver is the Julia port of RDCON's `gal.f` — Dewar's singular Galerkin method for -the outer-region Δ′ computation. Instead of integrating the Euler-Lagrange ODE radially, it -discretizes the displacement on a packed radial grid of Hermite-cubic finite elements and solves -a single global banded system. Cells adjacent to each rational surface ("resonant" and -"extension" cells) carry the large/small power-law asymptotic solutions analytically, so the -singular behavior is built into the basis rather than resolved numerically. - -The ideal coefficient is reconstructed as ``F = (QL)(QL)^\dagger`` from the interpolated -Cholesky factor ``L``, with ``Q = \operatorname{diag}(m-nq)``. This uses the same coefficient -representation as the Euler–Lagrange integration and singular-surface expansions, including -between spline knots, and preserves Hermitian positive semidefiniteness of ``F``. - -The solve produces the inter-surface Δ′ matrix and the PEST-3 matching blocks -(``A'``, ``B'``, ``\Gamma'``, ``\Delta'``), written to the HDF5 output under the same -`SingularSurfaces/` paths the Riccati Δ′ uses (`Delta_prime_matrix`, `Delta_prime_raw`, -`Delta_coil`, `pest3_A`, `pest3_B`, `pest3_Gamma`) — one set of names per quantity, -whichever formalism produced it. These are the outer-region inputs to resistive matched-asymptotic stability analysis -[Glasser 2016, Phys. Plasmas **23**, 072505]. - -Select it with `integrator = "Galerkin"` in `[ForceFreeStates]`; its settings go in `[Galerkin]`. 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. -A `[Match]` section additionally matches the inner layer, giving a driven ξ solution. The -composite inner-region profiles need the cut solution, `[Galerkin] cut_solution = true` (or -`Galerkin(; cut_solution=true)` followed by `solve(MatchProblem(...))`). - -The banded solve automatically tries a Hermitian Cholesky -factorization of the lower band (LAPACK `zpbtrf`/`zpbtrs`, with the Dirichlet rows' columns moved -to the right-hand side). The free block's off-diagonal entries must be Hermitian within -`1e-10 * sqrt(abs(A_ii)) * sqrt(abs(A_jj))`; the diagonal's imaginary part is deliberately -discarded. Unequal bandwidths, off-diagonal asymmetry, or failure of positive definiteness -trigger a reported fallback to Jacobi-scaled LU. The factorization actually used is -recorded as `ForceFreeStates/Solutions/GalerkinIntegration/solver` in `gpec.h5`. - -The implementation lives in `src/ForceFreeStates/Galerkin/`: - -| File | Content | -|------|---------| -| `GalerkinStructs.jl` | Cell, interval, and workspace types mirroring the Fortran derived types | -| `GalerkinGrid.jl` | Packed grid construction and local→global DOF mapping | -| `GalerkinAssembly.jl` | Element-level assembly: Hermite basis, Gauss-Lobatto stiffness, resonant and extension cells, boundary conditions | -| `GalerkinSolution.jl` | Reconstruct ξ(ψ) and analytic ξ′(ψ) on the gal-native grid | -| `GalerkinSolve.jl` | Top-level driver `galerkin_solve`, banded solve, Δ′ extraction, PEST-3 blocks, HDF5 output | - -## API Reference +# Galerkin Solver + +The Galerkin solver is the third `[ForceFreeStates] integrator`: the singular Galerkin method of the Fortran RDCON for the outer-region Δ′ of resistive stability ([Glasser, Wang and Park (2016)](citations.md#Tearing)). + +## What it computes + +Instead of integrating the Euler–Lagrange system radially, it discretizes the displacement on a packed radial grid of Hermite-cubic finite elements and solves one global banded system. +The cells adjacent to each rational surface carry the large and small power-law asymptotic solutions analytically, so the singular behavior is built into the basis rather than resolved numerically, and the solve yields the Δ′ matrix and the PEST-3 matching blocks of every surface at once. + +## Governing equations + +The ideal coefficient is reconstructed as ``F = (QL)(QL)^\dagger`` from the interpolated Cholesky factor ``L`` with the resonance factor ``Q = \operatorname{diag}(m - nq)``, the same representation the Euler–Lagrange integration and the singular-surface expansions use, so ``F`` stays Hermitian positive semidefinite between spline knots. +The weak form of the same system as on the [Stability](stability.md) page is assembled on the Hermite basis, with the resonant and extension cells next to each rational surface carrying the asymptotic solutions ``|\psi - \psi_s|^{-1/2 \pm \sqrt{-D_I}}``. + +## Running it + +Select it in `[ForceFreeStates]` and set its options in `[Galerkin]`, as in `examples/DIIID-like_gal_resistive_example` (abridged): + +```toml +[ForceFreeStates] +integrator = "Galerkin" + +[Galerkin] +nx = 256 # Hermite-cubic elements per interval between singular surfaces +pfac = 0.03 # grid packing ratio toward the singular surfaces +rpec_flag = true # edge coil-response columns (Delta_coil), needed by [Match] and [Tearing] +``` + +`cut_solution = true` (default `false`) adds the cut solution that composite inner-region profiles need. + +It replaces the radial integration rather than supplementing it: the run computes its own vacuum response at the control surface (with `vac_flag`) but produces no free-boundary energies, no ODE trace and no `crit` scan, it does not support the kinetic solve, and its result cannot drive `[PerturbedEquilibrium]`, which needs the free-boundary δW of a `"Forward"` solve. +A `[Match]` section matches the inner layer to its Δ′ and gives a driven ``\xi`` solution; the composite inner-region profiles need `cut_solution = true` (`Galerkin(; cut_solution=true)` followed by `solve(MatchProblem(ffs))` in a script). + +## Outputs + +| Dataset | Meaning | +|---|---| +| `SingularSurfaces/Delta_prime_matrix` | the Δ′ matrix in the PEST-3 convention, on the same path as the Riccati result, so downstream stages never branch on which formalism ran | +| `SingularSurfaces/Delta_prime_raw`, `Delta_coil` | the side-major matching matrix and the coil-response columns (`rpec_flag = true`) | +| `SingularSurfaces/pest3_A`, `pest3_B`, `pest3_Gamma` | the PEST-3 interchange and mixed blocks ``A'``, ``B'``, ``\Gamma'`` | +| `ForceFreeStates/Solutions/GalerkinIntegration/` | the closed ``\xi`` profiles on the Galerkin grid, the surface list, the factorization used (`solver`) and, with `basis_output = true`, the `Basis/` group | + +Nothing is written under `ForceFreeStates/FreeBoundaryStability/`. + +## Numerical notes + +- The radial grid is packed toward the rational surfaces and the edge, and the banded system is factorized by a Hermitian Cholesky decomposition, falling back to a Jacobi-scaled LU when the band is not Hermitian positive definite; the factorization used is recorded as `solver`. +- The singular cells make the solution insensitive to the resolution at the rational surfaces, so convergence is set by the smooth part of the grid (`nx`). + +## API reference ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.ForceFreeStates] Pages = ["Galerkin/GalerkinStructs.jl", "Galerkin/GalerkinGrid.jl", "Galerkin/GalerkinAssembly.jl", "Galerkin/GalerkinSolution.jl", "Galerkin/GalerkinSolve.jl"] ``` -## See also - -- `docs/src/stability.md` — the ForceFreeStates module this solver belongs to (EL system, metric matrices, singular surfaces) -- `docs/src/inner_layer.md` — inner-layer solutions to be matched against the outer-region Δ′ +See [Stability](stability.md) for the module this solver belongs to and [Tearing](inner_layer.md) for the inner-layer models matched to its Δ′. diff --git a/docs/src/index.md b/docs/src/index.md index 287c9e46e8..aa6f53553e 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -13,18 +13,47 @@ readme = replace(readme, r"^Full documentation:.*$"m => "") Markdown.parse(readme) ``` -## Overview -- [**Introduction Slides**](https://docs.google.com/presentation/d/1Iey7bsvpDG1WDRow24OD0NOWyKYOXMs0U7QLc7Y-wj4/edit?slide=id.g2798d4bd85b_0_15#slide=id.g2798d4bd85b_0_15) - GPEC Introduction slides +## What GPEC computes +One deck, `gpec.toml`, runs the stages below in order, each when its section is present, into one output file, `gpec.h5`. +The [Workflow](workflow.md) page maps every stage to its physics, deck sections and datasets. + +| Stage | Deck section | Main output | +|---|---|---| +| `Equilibrium` — the axisymmetric equilibrium, its flux-surface coordinates and rational surfaces | `[Equilibrium]` | `Equilibrium/` | +| `ForceFreeStates` — ideal (or drift-kinetic) MHD stability by the Newcomb criterion, with the `Vacuum` response of the `[Wall]` (a wall case is a separate deck from a response case) | `[ForceFreeStates]` | `ForceFreeStates/FreeBoundaryStability/`, `SingularSurfaces/Delta_prime_matrix` | +| `LocalStability` — Mercier, resistive-interchange and ballooning criteria | `[LocalStability]` | `LocalStability/` | +| `Match` — the resistive inner layer matched to the outer-region Δ′ | `[Match]` | `SingularSurfaces/Match/` | +| `PerturbedEquilibrium` — the plasma response to the field of `[ForcingTerms]`: resonant fields, islands, overlap; needs a `"Forward"` solve without a wall | `[PerturbedEquilibrium]` | `PerturbedEquilibrium/SingularCoupling/` | +| `KineticForces` — the neoclassical toroidal viscosity torque of that response | `[KineticForces]` | `KineticForces/` | +| `ErrorFields` — coil tolerances, locking risk and correction limits | `[ErrorFields]` | `ErrorFields/Risk/` | +| `Tearing` — tearing growth rates and critical resonant fields from an inner-layer model | `[Tearing]` | `Tearing/Roots/` | + +### Coming from the Fortran codes + +| Fortran | Here | +|---|---| +| DCON | `[ForceFreeStates]` with `integrator = "Forward"` | +| STRIDE | `integrator = "Riccati"`, with Δ′ in `SingularSurfaces/Delta_prime_matrix` | +| RDCON and its `rmatch` | `integrator = "Galerkin"` and `[Match]` | +| VACUUM | the `Vacuum` module, with the wall in `[Wall]` | +| GPEC | `[ForcingTerms]` and `[PerturbedEquilibrium]` | +| PENTRC | `[KineticForces]` | +| SLAYER | `[Tearing]` with `inner_model = "SLAYER"`, options in `[SLAYER]` | +| `dcon.in`, `vac.in`, `gpec.in`, `pentrc.in`, … | one `gpec.toml`, and one `gpec.h5` for every result | + +The [Conventions Reference](conventions.md) states the coordinates, signs and normalizations, including where they differ from the Fortran outputs. ## Documentation -- [**Setup**](set_up.md) — Installation instructions for macOS and Windows (WSL) -- [**Workflow**](workflow.md) — Full pipeline description: inputs, outputs, and key physics per module -- [**API Reference**](vacuum.md) — Detailed function and type documentation for each module +- [**Setup**](set_up.md) — Installation on Linux, macOS and Windows (WSL), and a first run +- [**Workflow**](workflow.md) — The pipeline: physics of each stage, deck sections, examples and the output file +- [**Scripting API**](api.md) — The same pipeline from Julia scripts +- [**Conventions Reference**](conventions.md) — Coordinates, signs, normalizations and field amplitudes +- [**Modules**](equilibrium.md) — Governing equations, deck keys, outputs and function reference per module - [**Citations**](citations.md) — Papers underpinning GPEC's algorithms -- [**Developer Notes**](developer_notes.md) — Git workflow, commit conventions, and the regression-testing harness - +- [**Developer Notes**](developer_notes.md) — Git workflow, commit conventions and the regression-testing harness +- [**Introduction slides**](https://docs.google.com/presentation/d/1Iey7bsvpDG1WDRow24OD0NOWyKYOXMs0U7QLc7Y-wj4/edit?slide=id.g2798d4bd85b_0_15#slide=id.g2798d4bd85b_0_15) — An overview of GPEC and its physics ## Contact diff --git a/docs/src/inner_layer.md b/docs/src/inner_layer.md index dc8d93fa90..e67b984ffd 100644 --- a/docs/src/inner_layer.md +++ b/docs/src/inner_layer.md @@ -1,282 +1,150 @@ -# Tearing Module +# Tearing -The `Tearing` module groups the resistive tearing-mode analysis stack: -`InnerLayer` (per-surface inner-layer matching data Δ(Q) for the GGJ and -SLAYER models), `Dispersion` (physics-agnostic complex-plane scan and -contour-intersection root extraction), `CriticalResonantField` (torque-balance -error-field penetration threshold), and the tearing solve itself, -`solve(TearingProblem(ffs; ...), model; config=TearingConfig(...))`, with its -`TearingResult` and HDF5 output. +`Tearing` solves the resistive tearing stability of each rational surface by matching the outer-region Δ′ of the force-free states to the response of a resistive inner layer. +The outer region, ideal MHD everywhere except at the rational surfaces, gives the tearing drive Δ′; the inner layer, where resistivity, two-fluid effects and rotation act, gives the response Δ(Q) as a function of the complex mode frequency; the roots of Δ′ = Δ(Q) are the tearing modes, with their growth rates and real frequencies. -## Tearing solve and layer inputs +## What it computes -```@autodocs -Modules = [GeneralizedPerturbedEquilibrium.Tearing] -``` +Two layer models are available, selected by `[Tearing] inner_model`: -The module provides an abstract [`InnerLayerModel`](@ref InnerLayer.InnerLayerModel) interface and two -concrete models: the **GGJ (Glasser–Greene–Johnson)** layer documented below, with three -interchangeable solver backends, and the pressureless **SLAYER** Fitzpatrick Riccati layer. -Further inner-layer models (kinetic layers) plug in through the same interface. - -Two papers define the equations and the asymptotic construction; the code and this page cite -them by equation number rather than restating them: - -- **GWP2016** — A. H. Glasser, Z. R. Wang and J.-K. Park, *Computation of - resistive instabilities by matched asymptotic expansions*, Phys. Plasmas - **23**, 112506 (2016). -- **GW2020** — A. H. Glasser and Z. R. Wang, *Asymptotic solutions and - convergence studies of the resistive inner region equations*, Phys. Plasmas - **27**, 012506 (2020). - -## Equations and matching data - -The GGJ layer equations (GWP2016 Eq. 11 ≡ GW2020 Eq. 1) couple the perturbed flux ``\Psi`` and -the layer variables ``\Xi``, ``\Upsilon`` in the scaled distance ``x`` from the rational surface. -Their equilibrium coefficients ``E, F, G, H, K`` (GWP2016 Eq. A8) are held in -[`GGJParameters`](@ref InnerLayer.GGJParameters). A physical growth rate ``\gamma`` maps to the -scaled ``Q`` by [`inner_Q`](@ref InnerLayer.inner_Q) (GWP2016 Eqs. A11, A14–A16), and -[`rescale_delta`](@ref InnerLayer.rescale_delta) returns the matching data to physical units. The -model requires Mercier stability, ``D_I < 0`` ([`mercier_di`](@ref InnerLayer.mercier_di), -GWP2016 Eq. A9), so that at large ``x`` the solutions include two power laws (GW2020 Eq. 49). - -The matching datum ``\Delta`` is the ratio of the large to the small power-law coefficient. -Imposing the two parities at ``x = 0`` gives the pair ``(\Delta_+, \Delta_-)`` of GWP2016 -Eqs. (34)–(35). Parity is that of the displacement ``\Xi`` (GWP2016 Eq. 34; Pletzer, Bondeson & -Dewar, J. Comput. Phys. **115**, 530 (1994), Eqs. 32–33): ``\Delta_+`` is even, -``\Psi(0)=\Xi'(0)=\Upsilon'(0)=0``, the **interchange** channel; ``\Delta_-`` is odd, -``\Psi'(0)=\Xi(0)=\Upsilon(0)=0``, the **tearing** channel. Positional pairs, such as `deltar` in -the matching and `RaySolveResult.Δ`, are always ordered ``(\Delta_+, \Delta_-)``, like the -unknowns ``d_\pm`` of GWP2016 Eq. (38) and Pletzer, Bondeson & Dewar Eqs. 83, 85; -`InnerLayerResponse` names them `interchange` and `tearing`. - -## Numerical method - -!!! note "Benchmark surface used in the figures" - The figures use the ``q = 4`` DIII-D-like surface of - [`q4_surface_benchmark`](@ref InnerLayer.GGJ.q4_surface_benchmark) - (``S \approx 4.6\times10^{6}``, ``D_I \approx -0.31``). - -### Solver backends - -The `GGJ` model exposes three solvers through the `solver` type parameter of -[`GGJModel`](@ref InnerLayer.GGJModel); all use the same large-``x`` asymptotic basis and return -``\Delta`` in the same convention: +- `"SLAYER"`, the two-fluid drift-MHD slab layer equations of Cole and Fitzpatrick (2006), in the form Park's SLAYER solver takes them with a few terms dropped ([Park (2022)](citations.md#Tearing), [Burgess et al. (2026)](citations.md#Tearing)), with resistivity, perpendicular and toroidal transport, diamagnetic drifts and E×B rotation taken from the kinetic profiles of the equilibrium; its options are in `[SLAYER]`. +- `"GGJ"`, the toroidal resistive-MHD layer of Glasser, Greene and Johnson ([Glasser, Wang and Park (2016)](citations.md#Tearing), [Glasser and Wang (2020)](citations.md#Tearing)), whose coefficients E, F, G, H, K are the flux-surface averages stored under `SingularSurfaces/`; its solver options are in `[GGJ]`. -| backend | method | robust regime | -|:--------|:-------|:--------------| -| `:shooting` | backward stable shoot from ``X_\mathrm{max}\to 0`` | ``\lvert Q\rvert \ll 1`` | -| `:galerkin` | Hermite-cubic finite-element (real axis) | surface- and ``\arg Q``-dependent; see [Choosing a backend](@ref) | -| `:ray` (default) | rotated contour: inward-marched far-field condition + Hermite-cubic BVP | ``\lvert Q\rvert \sim 500`` on/near the imaginary axis | +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. + +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/`. -Rotation and resistivity scans live on the imaginary-``Q`` axis, where the pseudo-resonance at -``x^2 = -Q^2(G+KF)`` sits on the real-``x`` path. On the ``q = 4`` surface there `:shooting` and -`:galerkin` lose accuracy beyond ``|Q| \sim 1`` and ``|Q| \sim 4`` respectively (the poles of ``\Delta(Q)`` -themselves lie near the real-``Q`` axis): +## Governing equations -![Relative error of the :shooting and :galerkin backends against the :ray reference along the imaginary-Q axis, on the q=4 benchmark surface. Each curve's crossing of the 1% line marks that method's practical reach.](figures/inner_layer/backend_regime_map.png) +The GGJ layer equations couple the perturbed flux ``\Psi`` and the layer variables ``\Xi``, ``\Upsilon`` in the scaled distance ``x`` from the rational surface (Glasser, Wang and Park (2016), Eq. 11). +A physical growth rate ``\gamma`` maps to the scaled frequency ``Q`` by [`inner_Q`](@ref InnerLayer.inner_Q), the model requires Mercier stability ``D_I < 0`` ([`mercier_di`](@ref InnerLayer.mercier_di)), and at large ``x`` the solutions contain the two power laws of the outer region. +The matching datum ``\Delta`` is the ratio of the large to the small power-law coefficient; imposing the two parities at ``x = 0`` gives the pair ``(\Delta_+, \Delta_-)``: ``\Delta_+`` even in the displacement ``\Xi``, the **interchange** channel, ``\Delta_-`` odd, the **tearing** channel. +Positional pairs are always ordered ``(\Delta_+, \Delta_-)``, and `InnerLayerResponse` names them `interchange` and `tearing`. +The dispersion relation of a surface is then ``\Delta'_s = \Delta_-(Q)`` with ``\Delta'_s`` the surface's diagonal of `SingularSurfaces/Delta_prime_matrix` converted to the layer's normalization (`PerSurface/Delta_prime_norm`); `coupling_mode = "coupled"` solves the determinant of the full matrix over up to `msing_max` surfaces instead. -### The rotated ray +The SLAYER layer is a slab with the Lundquist number ``S``, the temperature ratio ``\tau = T_i/T_e``, the compressibility ``c_\beta``, the perpendicular and toroidal Prandtl numbers ``P_\perp``, ``P_\mathrm{tor}`` and the diamagnetic frequencies ``Q_e``, ``Q_i`` of each surface (all stored under `Tearing/PerSurface/`), solved by a Riccati integration across the layer; its layer thickness ``\delta_s`` is written to `Tearing/LayerWidths/` and feeds the layer-overlap cap of [Stability](stability.md). -`:ray` writes the layer equations as ``dv/dx = M(x)\,v`` in the plain state -``v = (\Psi, \Xi, \Upsilon, \Psi', \Xi', \Upsilon')``. ``M`` is polynomial in ``x``, so -``x = 0`` is an ordinary point and the system continues analytically onto the ray +### Critical resonant field + +With `[Tearing.CriticalResonantField]` each `SLAYER` surface gets the critical resonant field for error-field penetration, from the steady-state torque balance of Cole and Fitzpatrick (2006): the electromagnetic torque on the layer balances the viscous restoring torque when ```math -x = e^{i\theta}\, s, \qquad s \in [0, S], \qquad \theta = \tfrac14\arg Q . +\frac{\operatorname{Im}\hat\Delta(Q)}{|\alpha + \hat\Delta(Q)|^2} + = \frac{2P\,(Q_0 - Q)}{S\hat\kappa\,(b_r/B_\phi)^2}, +\qquad +\hat\kappa = \left[\frac{2}{s}\right]^2 \int_{r_s}^{a} \frac{\mu(r_s)}{\mu(r)}\,\frac{dr}{r}, ``` -This angle makes the parabolic-cylinder exponent of the outer solutions real and lifts the -pseudo-resonance off the contour; ``\theta = 0`` is a real-axis solve. +and balance is lost above -![The rotated integration ray in the complex layer-coordinate plane at Q = 500i. The real-axis contour runs through the pseudo-resonance x² = −Q²(G+KF); the ray at θ = arg(Q)/4 = 22.5° clears it.](figures/inner_layer/rotated_ray_contour.png) +```math +\left(\frac{b_r}{B_\phi}\right)^2_\mathrm{crit} + = \max_Q \frac{2P\,(Q_0 - Q)}{S\hat\kappa\,\operatorname{Im}[-1/(\alpha + \hat\Delta(Q))]} . +``` -### Hermite-Galerkin BVP +``P`` is the layer's `P_tor`, ``S`` its Lundquist number, ``Q_0 = \tau_k n \omega_E`` the natural E×B rotation of the mode (from the kinetic file's `omega_E` or the `omega_E_kHz` override, a `[Tearing]` table keyed by surface such as `omega_E_kHz = { "2/1" = 3.0 }`; a surface without one is skipped), and the viscosity integral is taken as 1/2, so ``\hat\kappa = 2/s^2``. +Each surface uses its own stability parameter ``\alpha = S^{-1/3}(\Delta_\mathrm{crit} - r_s\Delta')``, so the torque balance shares its poles with the free modes; the threshold needs a tearing-stable surface, and one with ``\alpha \le 0`` is skipped with a warning. +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. + +## 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: + +```toml +[Tearing] +inner_model = "SLAYER" # or "GGJ" +coupling_mode = "uncoupled" # one dispersion relation per surface; "coupled" uses the full Δ′ matrix + +[Tearing.search] +scan_mode = "amr" # adaptively refined scan of the complex Q plane +Q_re_range = [-2.0, 2.0] +Q_im_range = [-0.5, 3.0] +nre = 41 +nim = 31 +pole_threshold_adaptive = true +filter_above_poles = true +filter_outside_re = true +polish_roots = true +store_scan = false # true: write Tearing/Scan/Surface_/ + +[Tearing.CriticalResonantField] +n = 2000 # Q samples of the torque balance +store_scan = false + +[SLAYER] +D_c_type = "none" # critical-Δ offset from parallel transport ("none": zero) +P_perp_model = "D_perp" # how the perpendicular Prandtl number is built (below) +``` -On ``[0, s_m]`` the system is solved with the `:galerkin` backend's Hermite-cubic assembly -(GWP2016 Eq. 32) written on the ray, ``A u_{ss} + e^{i\theta}B u_s + e^{2i\theta}C u = 0``, -with the parity conditions at ``s = 0``. The far end carries the point condition +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. -```math -v(s_m) = U_s + \Delta\, U_b + c_1 E_1 + c_2 E_2 , -``` +### SLAYER inputs: perpendicular diffusivities -where ``U_s, U_b`` are the small and large power-like solutions and ``E_{1,2}`` the -forward-decaying exponential pair; ``\Delta, c_1, c_2`` replace the last node's six unknowns. -Each cell resolves the local WKB rate, with a minimum cell count for small ``|Q|``. - -### Far-field boundary and the inward march - -The boundary data come from the same `inps` Wasow basis as the other backends (GW2020 -Eqs. 3–53), evaluated at the radius ``S`` where its residual is below tolerance -([`pick_smax`](@ref InnerLayer.pick_smax)). At large ``|Q|``, ``S`` lies far beyond ``s_m``, and -the power pair is marched inward by 2-stage Radau IIA modulo the decaying pair, the subspace in -which ``\Delta`` is defined. The decaying pair grows under backward integration at a rate -``\sim\rho s``. Radau's stability function decays as ``2/z``, so steps with every growing mode at -``z > 8`` damp those directions while the slowly varying power pair stays resolved; the step count -grows only like ``\ln(S/s_m)``. Nearer ``s_m`` the march resolves the growing directions, carries -the decaying frame, and projects it out periodically. The damped zone runs in `Complex{Double64}` -because its ill-conditioned implicit solves mix the power pair at a level ``\sim\epsilon\rho S^2``. - -On the real axis with imaginary ``Q`` the switch to resolved steps lands on the -pseudo-resonance, a turning point where the decaying directions cannot be identified; the -default ``\theta = \arg(Q)/4`` keeps it off the contour. - -For matching, `solve_inner_profile` returns ``\Delta`` and the layer-center fields -``\Psi(0), \Xi(0)`` from the rotated solve; ``x = 0`` lies on every contour, so these are the -real-axis values. Fields across the layer on the real ``\psi`` axis are an opt-in diagnostic -(`GGJModel(; real_profiles=true)`, not exposed in the deck): they come from an extra ``\theta = 0`` -solve, while the ``x = 0`` values, and so the penetrated field, stay those of the rotated solve. -The ``\theta = 0`` solve's ``\Delta_\pm`` are checked against the rotated ones (`certΔ`, which -warns on disagreement, as can happen at large ``|Q|`` near the imaginary axis). - -![Inner-layer fields Ψ, Ξ, Υ on the rotated ray for the q=4 surface at Q = 2i. The BVP solution (solid) joins the asymptotic representation (dashed) seamlessly at the match point S = s_m.](figures/inner_layer/solution_profiles.png) - -## Validation and benchmarks - -**Cross-check against the Fortran `rmatch` solver.** At the Glasser & Wang -(2020) Eq. (55) operating point ``Q = 0.1234`` (real) the `:galerkin` backend -reproduces the Fortran `rmatch deltac` solver to ``\sim 10^{-8}``: - -| quantity | `:galerkin` (= Fortran to ``10^{-8}``) | `:ray` | -|:---------|:---------------------------------------|:-------| -| ``\Delta_+`` (interchange) | ``3.698368\times 10^{4}`` | ``3.69789\times 10^{4}`` | -| ``\Delta_-`` (tearing) | ``14.759721`` | ``14.759715`` | - -This point sits near a pole of ``\Delta_+(Q)``, which amplifies every solver's error; hence the -``1.3\times10^{-4}`` `:galerkin`↔`:ray` difference in ``\Delta_+`` against ``4\times10^{-7}`` in ``\Delta_-``. - -**Physical benchmark on the imaginary axis.** On the ``q=4`` rational-surface -benchmark ([`q4_surface_benchmark`](@ref InnerLayer.GGJ.q4_surface_benchmark), ``S \approx 4.58\times10^6``, -``D_I \approx -0.312``) the `:ray` backend is pinned at ``Q = 500i``, a regime -entirely beyond `:galerkin`: - -| quantity | value at ``Q = 500i`` | -|:---------|:----------------------| -| ``\Delta_+`` (interchange) | ``2.4720 + 13.3540\,i`` | -| ``\Delta_-`` (tearing) | ``0.13750 + 0.74275\,i`` | - -**Convergence and contour invariance.** ``\Delta`` does not depend on the contour angle, so -re-solving with each numerical knob perturbed in turn -([`delta_convergence`](@ref InnerLayer.delta_convergence)) gives an error bar: at ``Q = 500i`` the -worst case is ``\sim 5\times10^{-6}`` for ``\Delta_+`` and ``\sim 6\times10^{-7}`` for ``\Delta_-``. -Across machines and BLAS builds the values reproduce to ``\sim 10^{-5}`` relative, hence the five -significant figures above. - -![Relative change of Δ at Q = 500i under independent perturbations of each numerical knob (contour angle, mesh resolution, series order/radius, march tolerance, BVP edge s_m, decaying-pair separation). The worst-case spread is the reported error bar.](figures/inner_layer/convergence_Sinvariance.png) - -## Choosing a backend - -`:ray` is the **default** — `GGJModel()` constructs `GGJModel{:ray}()` — and the recommended -backend for every ``Q``. It uses the same Hermite-cubic weak form as `:galerkin`; they differ in -where and how the far-field condition is imposed. `:galerkin` imposes it on the real axis through -resonant and extension cells at the radius where the series residual reaches ``10^{-7}``, on a -fixed `nx`-cell grid. `:ray` takes the series data where the residual reaches `smax_tol` = -``10^{-9}`` ([`pick_smax`](@ref InnerLayer.pick_smax)), marches them inward to ``s_m``, where every -exponential family has decayed 45 e-folds from the origin, and resolves the local WKB rate in -every cell on ``[0, s_m]``. Both radii are computed per solve; `S`, `smax_tol` and `sm_fac` -override them, and [`delta_convergence`](@ref InnerLayer.delta_convergence) perturbs them to give -an error bar. - -`:galerkin` takes a few milliseconds per solve, 3–60 times less than `:ray`, but its accuracy -depends on the surface and on ``\arg Q`` even at small ``|Q|``. Against a tightened `:ray` -reference, on the GW2020 Eq. (55) surface it errs by ``7\times10^{-5}`` at ``Q = 0.1`` and -``3\times10^{-3}`` at ``Q = 0.3``; near ``\arg Q = 157.5^\circ`` it errs by order unity at -``|Q| = 0.1`` on that surface and by 0.4 at ``|Q| = 0.3`` on the ``q = 4`` surface. `:ray` held -``2\times10^{-6}`` on all three surfaces tested, for ``|Q|`` from 0.01 to 1. - -The backends take disjoint numerical-knob keywords: pass `GGJModel(solver=:galerkin)` explicitly -when supplying the Galerkin `nx`/`xfac` knobs; passing them to `:ray` throws a `MethodError`. +The layer widths and the critical island width of `SLAYER` take perpendicular diffusivities from the kinetic file (`chi_e`, `chi_i`, `chi_phi` profiles in m²/s). +Two formulas want different quantities: the layer width goes through `P_perp`, the diffusion of the total pressure, while the critical island width ``W_d \propto (\chi_\perp/\chi_\parallel)^{1/4}`` is set by the electron heat transport that flattens the electron temperature across the island. +`[SLAYER] P_perp_model` selects how `P_perp` is built: from the perpendicular particle diffusivity ``D_\perp`` that weights the electron and ion channels by their share of the pressure gradient (`"D_perp"`, the default), from `chi_e` alone (`"chi_perp_e"`, the single-χ convention of Fitzpatrick's TJ code), from `chi_i`, from the toroidal Prandtl number, from its ``c_\beta^2`` floor, or from a whole-plasma confinement time (`"tau_E"`, for a plasma without χ profiles). +`[SLAYER] chi_perp_Dc` sets the electron diffusivity of the critical island width alone, one value for every surface, for a run in which one χ⊥,e cannot serve both formulas; where an interpretive power balance cannot separate the species, treat it as a modeling choice and scan it. +Since the widths scale as ``\chi^{1/4}``, a factor-two change in χ moves them by about 19%, which can be enough to move the resistive-layer overlap cut of [Stability](stability.md) by one surface. +The DIII-D-like kinetic file of the examples carries analytic χ profiles of the order of power-balance values for DIII-D H-modes; they are not fits to a measured discharge, do not resolve the pedestal transport barrier, and `chi_e` is raised by a factor of two so that the default deck exercises the layer-overlap cut. -```julia -using GeneralizedPerturbedEquilibrium.InnerLayer +## Outputs -p = GGJ.q4_surface_benchmark() # GGJParameters for the q=4 benchmark surface -γ = 500im * GGJ.q0(p) # physical rate placing Q on the imaginary axis at 500i +| Dataset | Meaning | +|---|---| +| `Tearing/PerSurface/rational_psi`, `rational_q`, `m`, `n` | the analyzed surfaces | +| `Tearing/PerSurface/S`, `tau`, `c_beta`, `P_perp`, `P_tor`, `Q_e`, `Q_i`, `tau_k`, `tau_R`, `eta`, `rs`, … | the layer parameters of each surface (SLAYER); `E`, `F`, `G`, `H`, `K`, `M`, `tau_A` and `tau_R` for GGJ | +| `Tearing/PerSurface/Delta_prime_raw` and `Delta_prime_matrix_rs` (SLAYER) or `Delta_prime_matrix` (GGJ) | the outer matrix each layer was matched to | +| `Tearing/PerSurface/Delta_prime_norm`, `delta_prime_conversion` | the Δ′ normalization of the SLAYER layer | +| `Tearing/Roots/Q_root`, `omega`, `gamma`, `no_root` | the dispersion root of each surface in the lab frame (`frequency_frame` attribute), its real frequency and growth rate, and the flag of a surface with no usable root | +| `Tearing/LayerWidths/delta_s`, `delta_s_abs`, `d_beta` | the complex resistive layer thickness, its magnitude in metres, and the ion drift scale | +| `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 | -Δ = solve_inner(GGJModel(), p, γ) # Δ.interchange, Δ.tearing with the default :ray backend +## Numerical notes -# Full result (raw Δ, contour, mesh, nodal solution, boundary conditioning) via solve_ray: -res = solve_ray(p, GGJ.inner_Q(p, γ)) -res.Δ, res.breaks, res.bc_cond # res.Δ = (Δ₊, Δ₋) = (interchange, tearing) -``` +The `GGJ` model exposes three solvers through the `solver` type parameter of [`GGJModel`](@ref InnerLayer.GGJModel), all returning ``\Delta`` in the same convention: -## Critical resonant field +| backend | method | robust regime | +|:--------|:-------|:--------------| +| `:shooting` | backward stable shoot from ``X_\mathrm{max}\to 0`` | ``\lvert Q\rvert \ll 1`` | +| `:galerkin` | Hermite-cubic finite-element solve on the real axis | surface- and ``\arg Q``-dependent; see [Choosing a backend](@ref) | +| `:ray` (default) | rotated contour: inward-marched far-field condition and a Hermite-cubic BVP | ``\lvert Q\rvert \sim 500`` on or near the imaginary axis | -With a `[Tearing.CriticalResonantField]` section (`TearingProblem(ffs; critical_resonant_field=CriticalResonantFieldControl())`), each SLAYER surface also gets the -critical resonant field for error-field penetration, from the steady-state torque balance of -Cole and Fitzpatrick, Phys. Plasmas **13**, 032503 (2006). The electromagnetic torque on the -layer balances the viscous restoring torque when (Cole Eq. 61) +`:ray` continues the layer equations analytically onto the ray ``x = e^{i\theta}s`` with ``\theta = \arg(Q)/4``, which makes the exponent of the outer solutions real and lifts the pseudo-resonance at ``x^2 = -Q^2(G + KF)`` off the contour; it takes the far-field data from the asymptotic series where its residual is below tolerance ([`pick_smax`](@ref InnerLayer.pick_smax)), marches the power pair inward by an implicit Radau scheme that damps the growing exponential directions, and solves the Hermite-cubic boundary-value problem with the parity conditions at ``s = 0``. +``\Delta`` does not depend on the contour angle, so [`delta_convergence`](@ref InnerLayer.delta_convergence) re-solves with each numerical knob perturbed and reports the spread as an error bar. -```math -\frac{\operatorname{Im}\hat\Delta(Q)}{|\alpha + \hat\Delta(Q)|^2} - = \frac{2P\,(Q_0 - Q)}{S\hat\kappa\,(b_r/B_\phi)^2}, -\qquad -\hat\kappa = \left[\frac{2}{s}\right]^2 \int_{r_s}^{a} \frac{\mu(r_s)}{\mu(r)}\,\frac{dr}{r}, -``` +![Relative error of the :shooting and :galerkin backends against the :ray reference along the imaginary-Q axis, on the q = 4 DIII-D-like benchmark surface; each curve's crossing of the 1% line marks that method's practical reach.](figures/inner_layer/backend_regime_map.png) -and balance is lost above (Cole Eq. 62) +Validation: at the real operating point ``Q = 0.1234`` of Glasser and Wang (2020) the `:galerkin` backend reproduces the Fortran `rmatch` solver to ``\sim 10^{-8}`` and `:ray` agrees to ``4\times10^{-7}`` in ``\Delta_-`` (``1.3\times10^{-4}`` in ``\Delta_+``, which sits near a pole); on the ``q = 4`` benchmark surface ([`q4_surface_benchmark`](@ref InnerLayer.GGJ.q4_surface_benchmark), ``S \approx 4.6\times10^6``, ``D_I \approx -0.31``) `:ray` is pinned at ``Q = 500i`` to ``\Delta_+ = 2.4720 + 13.3540\,i``, ``\Delta_- = 0.13750 + 0.74275\,i``, with a knob-perturbation error bar of ``\sim 5\times10^{-6}`` and ``6\times10^{-7}``. +The critical-field torque balance reproduces the Fortran SLAYER's balance to all printed digits when fed the same Δ(Q) table; the thresholds differ through the layer model (Park's reduced form of the layer equations against Cole and Fitzpatrick's full ones) and the selection rule above. -```math -\left(\frac{b_r}{B_\phi}\right)^2_\mathrm{crit} - = \max_Q \frac{2P\,(Q_0 - Q)}{S\hat\kappa\,\operatorname{Im}[-1/(\alpha + \hat\Delta(Q))]} . +![Fortran and Julia inner-layer Δ(Q) and torque balance on the DIII-D-like example.](assets/crf_fortran_benchmark.png) + +### Choosing a backend + +`:ray` is the default and the recommended backend for every ``Q``: it resolves the local WKB rate in every cell and holds ``\sim 2\times10^{-6}`` for ``|Q|`` from 0.01 to 1 on every surface tested. +`:galerkin` takes a few milliseconds per solve, 3–60 times less than `:ray`, but its accuracy depends on the surface and on ``\arg Q`` even at small ``|Q|``, erring by order unity near ``\arg Q = 157.5^\circ``. +The backends take disjoint numerical keywords: pass `GGJModel(solver=:galerkin)` explicitly when supplying its `nx`/`xfac` knobs. + +```julia +using GeneralizedPerturbedEquilibrium.InnerLayer + +p = GGJ.q4_surface_benchmark() # GGJParameters for the q = 4 benchmark surface +γ = 500im * GGJ.q0(p) # physical rate placing Q on the imaginary axis at 500i +Δ = solve_inner(GGJModel(), p, γ) # Δ.interchange, Δ.tearing with the default :ray backend ``` -- **Inputs.** P is each layer's `P_tor` (so χ_φ from the kinetic file, or the scalar - `chi_tor`), S its Lundquist number, and Q0 = τ_k·n·ω_E the natural E×B rotation of the - m/n mode (Cole's ω0; the diamagnetic drifts enter through Q_e and Q_i in Δ̂). Ω_E is the - kinetic file's `omega_E` or the `omega_E_kHz` override; the analysis needs it, and a surface - without one is skipped. -- **Approximations.** As in Fortran `gslayer.f`, the viscosity integral is taken as 1/2 - (κ̂ = 2/s²). -- **Stability α.** Each surface uses its own α = S^(−1/3)(Δ_crit − r_s Δ'), so α + Δ̂ is the - dispersion relation's residual and the torque balance shares its poles with SLAYER's free - modes (Fortran fixes α = 10⁻²). Cole and Fitzpatrick have no Δ_crit; including it treats the - curvature offset as independent of frequency. It is zero unless `D_c_type` is set. Cole's - threshold needs a tearing-stable surface, so a surface with α ≤ 0 is skipped with a warning. -- **Q axis.** Q is on Cole's axis, with the diamagnetic poles at Q_e and Q_i; it runs opposite to the dispersion - scans' Re Q, so the scan uses the layer response at −Q. -- **Threshold.** Penetration is associated with a saddle-node bifurcation (Fitzpatrick, Nucl. Fusion 33, 1049 (1993), Sec. 7.4): - as the field grows, Q moves from Q0 toward the electron fluid at Q_e, where the torque - vanishes, and the threshold is the first resolved saddle-node bifurcation on that path, above - which the shielded branch followed from Q0 no longer exists. Both the rise to and fall from - an interior maximum must exceed 0.1 % of its value within a contiguous retained interval. - This excludes monotonic clipped windows and small numerical extrema; excluded samples do - not connect separate intervals. A later, higher maximum does not replace the first resolved - saddle-node bifurcation. Points where |Δ̂| ≤ α, the - unshielded zone next to Q_e, are not candidates, so the balance spike there is never taken. - With no resolved saddle-node bifurcation no discrete threshold is reported (NaN). Fortran takes the global maximum on a - window instead, which picks that spike. The search assumes Q_i ≤ 0 ≤ Q_e; a surface where - that fails is skipped with a warning. - -Outputs go to `Tearing/CriticalResonantField/` (`br_crit` = b_r/B_φ, `br_crit_tesla`, `Q_peak`, `Q0`, with P in -`PerSurface/P_tor`, and with `store_scan` the per-surface `Scan/Surface_/{Q, balance, Delta}`). - -`br_crit` is normalized by the toroidal field, as in Cole and Fitzpatrick (2006, Eqs. 61–68) and -Fitzpatrick, Phys. Plasmas **30**, 092512 (2023), Eq. 120. `br_crit_tesla` = `br_crit`·|B0|, with -the layer's B0, is for comparison with `SingularCoupling/resonant_area_weighted_field`, the total -resonant field with the ideal plasma response, which Park et al., Phys. Rev. Lett. **99**, 195003 -(2007) identify as the field that opens islands. The comparison is approximate: the theory's b_r is -a cylindrical field at the layer, and no reference maps it onto the area-weighted toroidal field. -`penetrated_area_weighted_field` already includes the layer response and is not the matching quantity. - -### Benchmark against Fortran - -On the DIII-D-like example (2/1–6/1, n = 1, P = 1, κ̂ = 2/s², Q ∈ [−5, 5] with 20001 -points), both codes were given Fortran's normalized inputs and Fortran's fixed α = 10⁻². - -- **Torque balance.** Fed Fortran's own Δ(Q) table, the Julia routine reproduces Fortran's - balance to all printed digits on every surface, so the thresholds differ only through the - selection rule above. -- **Layer model.** Fortran solves Park's drift-MHD layer and Julia Fitzpatrick's, so Δ(Q) - differs by a median 24–78 % between Q_i and Q_e on 2/1, 3/1, 4/1 and 6/1 (7 % on 5/1). - -The dots in the figure mark each code's balance maximum between Q_e and Q0 at that fixed α, -not the fold threshold above. - -![Fortran vs Julia inner-layer Δ(Q) and torque balance](assets/crf_fortran_benchmark.png) - -## API Reference +## API reference + +### Tearing + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Tearing] +``` ### InnerLayer @@ -290,19 +158,19 @@ Modules = [GeneralizedPerturbedEquilibrium.InnerLayer] Modules = [GeneralizedPerturbedEquilibrium.InnerLayer.GGJ] ``` -## SLAYER +### SLAYER ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.InnerLayer.SLAYER] ``` -## Dispersion +### Dispersion ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.Tearing.Dispersion] ``` -## CriticalResonantField +### CriticalResonantField ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.Tearing.CriticalResonantField] diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md index 5a315ff9a0..d88681a6b2 100644 --- a/docs/src/kinetic_forces.md +++ b/docs/src/kinetic_forces.md @@ -1,232 +1,106 @@ -# KineticForces Module - -Kinetic torque and energy calculations for perturbed equilibria. -Implements neoclassical toroidal viscosity (NTV) from the PENTRC formulation -([Logan & Park, 2013](citations.md); [Logan, 2015](citations.md)). - -## Kinetic profile file formats - -Kinetic profiles are read by `read_kinetic_file`, which dispatches on the file -extension: - - - **HDF5 (`.h5`/`.hdf5`)** — the GPEC kinetic schema (recommended). Datasets - at the file root: required `psi` (normalized poloidal flux), `n_e`, `T_i`, - `T_e`, `omega_E` (and `n_i`, defaulting to `n_e` if omitted); optional - `omega_tor`, `chi_e` and `chi_i` (electron and ion perpendicular heat diffusivities - ``\chi_{\perp,e}``, ``\chi_{\perp,i}``), `chi_phi` (toroidal momentum diffusivity - ``\chi_\phi``), `zeff` - (effective charge ``Z_\mathrm{eff}``, dimensionless), and the supplied collision - rates `lnLambda` (dimensionless), `nu_i` and `nu_e` (1/s). Each dataset - carries a `units` attribute; the root carries `schema_version` and - `provenance`. Densities are m⁻³, temperatures eV, frequencies rad/s, - diffusivities m²/s. Additional datasets named `n_*` (e.g. `n_D`, `n_T`) - are named per-species densities for multi-ion runs; names outside `n_*` - are reserved for future schema fields and ignored. Write files with - `write_kinetic_h5` (which round-trips the per-species densities). - - **ASCII (`.gpeckf`/`.kin`/`.dat`)** — legacy six-column whitespace table - `psi_n n_i n_e T_i[eV] T_e[eV] omega_E`, retained for backward - compatibility, or, with `kinetic_format = "rates-v1"` in `[Equilibrium]`, nine columns that append - `lnLambda nu_i[1/s] nu_e[1/s]` (the Fortran PENTRC rates-v1 table). The column count must match - `kinetic_format`. Header rows are skipped. - -Supplied collision rates replace the internal Krook collision model of the NTV calculation -(Logan & Park 2013 Eq. 6), which otherwise computes them from ``n`` and ``T``: - - - Each of `lnLambda`, `nu_i` and `nu_e` is optional and must be positive. - - They are interpolated in log space. Outside the table's ψ range the end interval is extended in log - space, so supply them over the whole ψ range. - - A supplied `nu_i` is the first main ion's rate; the other ion species and the impurity keep the - model's ``(Z_s/Z_1)^2\sqrt{m_1/m_s}`` ratio to it. - - They feed only the NTV calculation; the tearing layers' resistivity keeps its `lnLambda_form` closure. - - The density and temperature scaling factors of `[KineticForces]` do not rescale supplied rates; - `collisionality_factor` still multiplies them. - -The NTV calculation consumes `n_i, n_e, T_i, T_e, omega_E`; `chi_e`/`chi_i`/`chi_phi`/`zeff`, -when present, are carried for the resistive-layer (SLAYER) analysis and ignored -here. That analysis takes ``Z_\mathrm{eff}(\psi)`` from the `zeff` dataset when the file -has one, otherwise from `n_i/n_e` with the kinetic data's `impurity_charge`, otherwise -from its `zeff_default`. Supply `zeff` whenever the impurity content is not a -single species fixed by `n_i/n_e`. - -The plasma's species composition (main ions, impurity charge and mass), the defaults that stand -in for absent χ and Z_eff profiles, the confinement time and the resistivity closure belong to the -equilibrium's kinetic data, not to any one physics stage, so KineticForces and the tearing layers -see one plasma. In a script, set them where the profiles are attached: +# Kinetic Forces -```julia -attach_kinetic_profiles!(equil, "kinetic.h5"; ion_species=[IonSpecies(; z=1, m=2, fraction=0.5), - IonSpecies(; z=1, m=3, fraction=0.5)], impurity_charge=6, impurity_mass=12) -``` +`KineticForces` computes the neoclassical toroidal viscosity (NTV) of a perturbed equilibrium: the toroidal torque and the kinetic energy that a non-axisymmetric field exerts on the plasma through the drift-kinetic response of its particles. +It implements the PENTRC formulation of [Park et al. (2009)](citations.md#KineticForces), [Logan et al. (2013)](citations.md#KineticForces) and [Park and Logan (2017)](citations.md#KineticForces). -or copy loaded data with new settings, `KineticProfileData(equil.kinetic; impurity_charge=6)`. -The default is one deuterium main ion carrying all of `n_i`, with a carbon impurity. In a deck -they are `[Equilibrium]` keys: `[[Equilibrium.ion_species]]`, `impurity_charge` and `impurity_mass`. +## What it computes -## Multi-ion runs +A non-axisymmetric perturbation ``\delta B`` breaks the toroidal symmetry that conserves the canonical angular momentum of each particle. +The resulting non-ambipolar radial transport, dominated by trapped particles whose bounce and precession motion resonates with the E×B rotation, carries toroidal momentum out of the plasma: the NTV torque ``T_\phi``, which brakes or spins up the rotation, together with a kinetic contribution ``\delta W_k`` to the perturbed energy that stabilizes or destabilizes the response. +Both are halves of one complex quantity, ``T = T_\phi + 2in\,\delta W_k``, and the module integrates it over the plasma for one or more torque `method`s. -A plasma may declare an arbitrary list of main-ion species; the NTV is computed -per species under one shared full-composition ``Z_\mathrm{eff}`` and summed -(``\tau = \sum_s \tau_s``) over the main ions, the quasineutrality-closing -impurity, and (with `electron = true`) the electrons. This applies to both NTV -paths: the post-PE ψ-quadrature diagnostic and the self-consistent -kinetic matrices of a `kinetic = "KineticForces"` force-balance solve. +The torque enters the pipeline in two ways: -```toml -[Equilibrium] -kinetic_file = "kinetic.h5" # n_i column/dataset = TOTAL main-ion density -impurity_charge = 6 # impurity charge (closes quasineutrality) -impurity_mass = 12 # impurity mass - -[[Equilibrium.ion_species]] -z = 1 -m = 2 -fraction = 0.5 # this species' share of the total n_i -[[Equilibrium.ion_species]] -z = 1 -m = 3 -density = "n_T" # or: explicit n_* dataset from the HDF5 kinetic file +| Workflow | Deck | What it gives | +|---|---|---| +| Torque of a perturbed equilibrium | `[KineticForces]` after `[PerturbedEquilibrium]` | ``T(\psi)`` integrated over the stored response; the stage of the [Workflow](workflow.md) | +| Self-consistent kinetic solve | `[ForceFreeStates] kinetic = "KineticForces"` | the drift-kinetic energy matrices replace the ideal ones inside the stability solve, so the response itself carries the torque (see [Stability](stability.md) and `TorqueResponse/` on the [Perturbed Equilibrium](perturbed_equilibrium.md) page) | + +Both need the kinetic profiles of `[Equilibrium] kinetic_file` (n_e, T_i, T_e, ω_E, and n_i when it differs from n_e) and the ion composition of `[[Equilibrium.ion_species]]`; the file schema is on the [Equilibrium](equilibrium.md) page. + +## Governing equations + +The torque is a sum over species of three nested integrals, over flux ``\psi``, over the pitch variable ``\Lambda`` separating trapped from passing particles and over the normalized energy ``x = \varepsilon/T``, with a sum over the bounce harmonics ``\ell`` up to `bounce_harmonic_max`. +The integrand is the bounce-averaged perturbation of the field strength ``\langle \delta B_\parallel / B \rangle_\ell``, evaluated from the displacement ``\xi`` of the perturbed equilibrium, weighted by the Maxwellian and by the resonance factor ``1/(i\Omega(x) - \nu)``: the energy integral is Eq. 8 of [Logan et al. (2013)](citations.md#KineticForces) and the torque and its normalization Eq. 19, ``\operatorname{Im} T = 2n\,\delta W_k``. +The resonance frequency ``\Omega(x)``, and the diamagnetic frequencies that relate the E×B rotation ``\omega_E`` of the kinetic file to the toroidal rotation, are given with their sign conventions under [KineticForces rotation conventions](conventions.md#KineticForces-rotation-conventions). +The collision operator is selected by `nutype`: `"harmonic"` (the default), `"krook"`, `"small"` or `"zero"` (the collisionless limit), with the collision frequencies computed from the density and temperature profiles; `f0type` selects the equilibrium distribution, `"maxwellian"` (the default), `"park"` or `"cgl"`. +The kinetic file may also supply the Coulomb logarithm and collision rates `lnLambda`, `nu_i` and `nu_e`, which then replace the collision model of the NTV calculation (Logan et al. (2013), Eq. 6) that otherwise computes them from ``n`` and ``T``. +Each is optional and must be positive; they are interpolated, and extended past the table's ψ range, in log space, so supply them over the whole range. +A supplied `nu_i` is the first main ion's rate, and the other ions and the impurity keep the model's ``(Z_s/Z_1)^2\sqrt{m_1/m_s}`` ratio to it. +They feed only the NTV calculation (the tearing layers' resistivity keeps its `lnLambda_form` closure), the density and temperature scaling factors do not rescale them, and `collisionality_factor` still multiplies them. + +### Torque methods + +Each `method` is one evaluation of the same integral under a different model, switched on by its `_flag` in `[KineticForces]`; the first letter selects the particle population, `f` full, `t` trapped only, `p` passing only. +| Method | Model | +|---|---| +| `fgar` (default), `tgar`, `pgar` | general-aspect-ratio calculation of the torque and energy | +| `fcgl` | fluid Chew–Goldberger–Low limit | +| `fwmm`, `twmm`, `pwmm` | energy from the Euler–Lagrange kinetic matrix | +| `ftmm`, `ttmm`, `ptmm` | torque from the Euler–Lagrange kinetic matrix | +| `fkmm`, `tkmm`, `pkmm` | norm of the kinetic energy matrix | +| `frmm`, `trmm`, `prmm` | norm of the kinetic torque matrix | + +## Running it + +From `examples/DIIID-like_ideal_example/gpec.toml`, after the `[PerturbedEquilibrium]` section (`examples/TJ_kinetic_NTV_example` is the same stage on an analytic equilibrium): + +```toml [KineticForces] -electron = true # add electron NTV in addition to the ion species +bounce_harmonic_max = 6 # bounce harmonics summed over ℓ = -max:max +electron = false # add the electron NTV to the ion species +nutype = "harmonic" # collision operator: "zero", "small", "krook", "harmonic" +f0type = "maxwellian" # distribution function: "maxwellian", "park", "cgl" +atol_pitch = 1e-9 # absolute tolerance of the pitch integration +axis_validity_suppression = true # suppress the near-axis terms where the model is not valid ``` -The composition describes the plasma, so it sits on the equilibrium's kinetic data beside -`kinetic_file`, where every kinetic consumer (NTV, the tearing layers) reads it. - -Each species sets exactly one of `fraction` (share of the file's total `n_i`) -or `density` (a named `n_*` dataset). An all-fraction list must sum to 1; in a -mixed list the impurity content is set by the file's `n_i`/`n_e` deficit, not -by a fraction shortfall, and fractions may sum below (never above) 1. Every -main-ion charge must satisfy `z < impurity_charge`. A single main ion carrying all of `n_i` -runs the single-species NTV (with `electron = true` still adding the electron species — the -electron flag, a `KineticForcesControl` setting, always means *in addition to* the ions). Without -an `ion_species` list the plasma is one deuterium main ion. - -The summed total is written to `KineticForces//` exactly as in a single-ion -run, and each species' own contribution to `KineticForces/PerSpecies/