diff --git a/docs/development/architecture.md b/docs/development/architecture.md index 781f0b155..c54e02c8d 100644 --- a/docs/development/architecture.md +++ b/docs/development/architecture.md @@ -210,7 +210,7 @@ The complete GPEC analysis pipeline: ### Stability - `SingType` - Singular surface data including: - Rational surface location (ψ, ρ, q = m/n, dq/dψ) - - Δ' (tearing stability parameter) — **stub**; the valid Δ' is `ForceFreeStatesResult.delta_prime.matrix` + - Δ' (tearing stability parameter) — **stub**; the valid Δ' is `ForceFreeStatesResult.delta_prime` (`DeltaPrimeData`) - Asymptotic solution bases at the inner-layer boundaries - `ForceFreeStatesResult` - Published product of a solve: mode space, metric/matrix fits, singular surfaces, and the per-integrator products (ξ solution and its basis, free-boundary energies, diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md index 403fb7ffc..7ec9aa01c 100644 --- a/docs/development/hdf5-conventions.md +++ b/docs/development/hdf5-conventions.md @@ -6,7 +6,7 @@ Conventions for the structure and naming of the `gpec.h5` output file. Every wri **The schema must be intuitive to a plasma physicist who is not a developer of this code.** The top level largely mirrors the TOML sections / major `src/` modules — which are themselves organized along physics lines — but **physics intuition wins whenever the two diverge**. Worked examples of that rule: -- All per-rational-surface stability results consolidate under `SingularSurfaces/`, regardless of which algorithm produced them: the GGJ coefficients sit beside the Δ′/PEST-3 matching results, and the Riccati BVP and the Galerkin outer-region solve write the *same* `Delta_prime_matrix`/`Delta_prime_raw`/`Delta_coil`/`pest3_*` paths rather than each owning a producer-named subgroup. +- All per-rational-surface stability results consolidate under `SingularSurfaces/`, regardless of which algorithm produced them: the GGJ coefficients sit beside the Δ′/PEST-3 matching results, and the Riccati BVP and the Galerkin outer-region solve write the *same* `Delta_prime_*`/`Delta_EF`/`pest3_*` paths rather than each owning a producer-named subgroup. - `EulerLagrangeMatrices` over a bare `Matrices` — group names must say what the data *is*, not which array it came from. Same reasoning renamed `matrices_/` → `KineticMatrices/`. Physics-topic groups elevated to top level (rather than nested under their producer): `Info/`, `Input/`, `SingularSurfaces/`, `LocalStability/`, `SurfaceGeometries/`. @@ -42,7 +42,7 @@ Top level (11 groups): | `Equilibrium/` | Scalars (`beta_N`, `q_axis`, `q_95`, `I_p`, …) plus `Profiles/` (1-D on `psi`: 2piF, mu0p, dVdpsi, q) and `Geometry/` (2-D on `psi`×`theta`: rcoords, offset, nu, jac) | | `ForceFreeStates/` | `Solutions/ForwardIntegration/` (u-solutions), `Solutions/GalerkinIntegration/` (closed ξ profiles in the shared layout, the gal surface list, `gal_basis_output`-gated `Basis/`), `write_el_matrices`-gated `EulerLagrangeMatrices/{Ideal,Kinetic}`, `FreeBoundaryStability/`, `EdgeScan/` | | `LocalStability/` | Mercier `D_I`, resistive interchange `D_R`, `ballooning_Delta_prime` on `psi`; the ballooning α boundary on `ballooning_psi` | -| `SingularSurfaces/` | Per-rational-surface data: `rational_psi`/`rational_q`/`rational_m`/`rational_n`, GGJ coefficients, `Delta_prime_matrix`/`Delta_prime_raw`/`Delta_coil`/`pest3_A`/`pest3_B`/`pest3_Gamma` (Riccati or Galerkin alike), the inner-layer `Match/` results, `Kinetic/` | +| `SingularSurfaces/` | Per-rational-surface data: `rational_psi`/`rational_q`/`rational_m`/`rational_n`, GGJ coefficients, `Delta_prime_raw_freeboundary`/`Delta_prime_raw_fixedboundary`/`Delta_EF`/`edge_response`, the free-boundary `Delta_prime_matrix`/`pest3_A`/`pest3_B`/`pest3_Gamma` and their `_fixedboundary` counterparts (Riccati or Galerkin alike), the inner-layer `Match/` results, `Kinetic/` | | `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/`, `TorqueResponse/` (kinetic forward solves: ψ-resolved `T_xe`, `T_source`), control-surface spectra | | `KineticForces/` | `/` (torque/energy profiles, `dTdpsi_ell` resolved by bounce harmonic on the `ell` scale, `KineticMatrices/`); multi-ion runs add `PerSpecies///` with the same per-method layout, summing to the top-level total | | `ErrorFields/` | `CoilSensitivities/` (per-coil-set control-surface spectra and their rigid shift/tilt derivatives, `DominantMode/` full-window projection); `MonteCarlo/` (intrinsic and corrected `\|δ\|` histograms over the sampled tolerances, per batch and averaged); `Risk/` (threshold density, `P(lock\|δ)`, locking probabilities, `ToleranceScan/`); `NTV/` (correction-coil overlap and NTV torque couplings per kAt) | diff --git a/docs/src/api.md b/docs/src/api.md index 4cb8938dd..ab0df25b9 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -129,7 +129,7 @@ ffs = solve(eq, Riccati(; nchunks=40); nn=1, vac_flag=true) ffs = solve(eq, Galerkin(; nx=512); nn=1) # finite-element outer-region Δ′ ``` -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 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, or a `Galerkin(; error_field_response=true)` solve with `vac_flag=true`, with `wall=nothing` (it errors on a wall, any other Galerkin result or a matched one, 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 @@ -137,13 +137,13 @@ A kinetic solve takes `Forward()` with the problem's `kinetic` keyword: `kinetic 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) +ffs = solve(eq, Galerkin(; error_field_response=true, cut_solution=true); nn=1) matched = solve(MatchProblem(ffs; eta=1e-6, rho=1e-7, rotation=0.0), GGJModel()) 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 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`. +The problem needs a Δ′ payload with error-field (EF) columns, so it accepts Galerkin and Riccati results with `error_field_response=true` (the default; Riccati also needs `vac_flag=true`); a Riccati-fed match fills `bpen`, the resonant data and the matched `wp` but keeps `solution === nothing`. `GGJModel()` is the closure-capable model; `SLAYERModel()` is slab-only and drives the tearing solve instead. ## Tearing stability diff --git a/docs/src/galerkin.md b/docs/src/galerkin.md index b3e59a745..4f559e791 100644 --- a/docs/src/galerkin.md +++ b/docs/src/galerkin.md @@ -23,24 +23,24 @@ 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] +error_field_response = true # error-field (EF) columns (Delta_EF): the ideal ξ solution and W_plasma, needed by [Match] ``` `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). +It replaces the radial integration rather than supplementing it: the run computes its own vacuum response at the control surface (with `vac_flag`), produces no ODE trace and no `crit` scan, and does not support the kinetic solve. With `error_field_response = true` and `vac_flag = true` it also produces the free-boundary energies of its ideal solution, so it can drive `[PerturbedEquilibrium]` like a `"Forward"` solve; a `[Match]`ed result cannot yet. +With `error_field_response = true` the error-field (EF) columns, one per edge mode with unit edge displacement and no inner-layer drive, are solved with the edge held fixed; their small-solution coefficients are `Delta_EF`, they are the ideal (perfectly shielded) ``\xi`` solution, and the force they exert on the edge is the fixed-boundary plasma energy matrix ``W_P`` of Eq. 66 of [Glasser (2016)](citations.md#ForceFreeStates). +A `[Match]` section matches the inner layer to its Δ′ and replaces them with the driven ``\xi`` solution, whose edge response ``W_P`` is complex: its anti-Hermitian part is the torque absorbed in the layers. 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'`` | +| `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; free-boundary with `vac_flag = true` (with `error_field_response = true` a second factorization re-solves the big-solution columns with the vacuum edge, so the EF columns never enter it; `free_boundary_delta_prime = false` skips it for matching-only scans, leaving no tearing Δ′) and fixed-boundary otherwise | +| `SingularSurfaces/Delta_prime_raw_freeboundary`, `Delta_prime_raw_fixedboundary`, `Delta_EF`, `edge_response` | the side-major matching matrices with the vacuum and the fixed edge, the error-field (EF) response (`error_field_response = true`), and the edge momentum of each fixed-edge big-solution column | +| `SingularSurfaces/pest3_A`, `pest3_B`, `pest3_Gamma` (and `_fixedboundary`) | 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/`. +| `ForceFreeStates/FreeBoundaryStability/` | `W_plasma`, the plasma energy matrix ``W_P`` of the closed solution (`error_field_response = true`), and with `vac_flag = true` the vacuum and free-boundary energies and eigenmodes as for a `"Forward"` solve | ## Numerical notes diff --git a/docs/src/index.md b/docs/src/index.md index aa6f53553..65c06a44c 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -24,7 +24,7 @@ The [Workflow](workflow.md) page maps every stage to its physics, deck sections | `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/` | +| `PerturbedEquilibrium` — the plasma response to the field of `[ForcingTerms]`: resonant fields, islands, overlap; needs a `"Forward"` (or `error_field_response` `"Galerkin"`) 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/` | diff --git a/docs/src/inner_layer.md b/docs/src/inner_layer.md index e67b984ff..fd2fa7374 100644 --- a/docs/src/inner_layer.md +++ b/docs/src/inner_layer.md @@ -50,7 +50,7 @@ Penetration is the first resolved saddle-node bifurcation on the path from ``Q_0 ## 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: +From `examples/DIIID-like_SLAYER_example/gpec.toml`, after a `[ForceFreeStates]` solve with `integrator = "Riccati"` or `"Galerkin"` and `vac_flag = true`, which provides the Δ′ matrix, and with a kinetic file on the equilibrium: ```toml [Tearing] diff --git a/docs/src/perturbed_equilibrium.md b/docs/src/perturbed_equilibrium.md index b5b9d27a9..db15d4a46 100644 --- a/docs/src/perturbed_equilibrium.md +++ b/docs/src/perturbed_equilibrium.md @@ -92,7 +92,7 @@ The coupling matrices `C_*` map the applied ``\tilde b`` to each of these quanti ## Running it The response needs a force-free-states solve with `integrator = "Forward"` and no `[Wall]`. -`"Forward"` is the only integrator that retains the dense ``\xi`` solution the operators are built from: a `"Riccati"` solve (the deck default) makes `plasma_response` warn that the response is skipped and return empty operators, so the failure surfaces only when a field is applied, and a `"Galerkin"` or `[Match]`ed solve, which publishes no free-boundary δW, is an error. +The operators are built from a dense ``\xi`` solution with its free-boundary δW, which a `"Forward"` solve and a `"Galerkin"` solve with `error_field_response = true` and `vac_flag = true` both provide: a `"Riccati"` solve (the deck default) makes `plasma_response` warn that the response is skipped and return empty operators, so the failure surfaces only when a field is applied, while a `"Galerkin"` solve without those flags or a `[Match]`ed solve, whose singular coupling is not defined yet, is an error. The ψ-resolved torque response (`compute_torque_response`) needs a `"Forward"` kinetic solve. A solve with a wall is an error too, because the surface inductance is the free-space one and ``\Lambda`` from a wall run would mix two vacuum models; the wall case of a stability study is therefore a separate deck. From `examples/DIIID-like_ideal_example/gpec.toml`, after the `[ForceFreeStates]` solve and a `[ForcingTerms]` drive: diff --git a/docs/src/stability.md b/docs/src/stability.md index aa86df2ec..87f9adcc2 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -59,7 +59,7 @@ qlow = 1.02 # inner edge of the integration domain qhigh = 1e3 # outer bound of qlim (see Integration domain) delta_mlow = 8 # extra poloidal modes below the resonant range delta_mhigh = 8 # extra poloidal modes above it -integrator = "Forward" # "Riccati" is the default; "Forward" (and no [Wall]) is required for [PerturbedEquilibrium] +integrator = "Forward" # "Riccati" is the default; "Forward" or an error_field_response "Galerkin" (and no [Wall]) is required for [PerturbedEquilibrium] truncate_at_rational_offset = true rational_offset_fraction = 0.2 @@ -81,7 +81,7 @@ ffs = solve(eq, Riccati(; eulerlagrange_tolerance=1e-10); nn=1, delta_mlow=8, de real(ffs.free_boundary.et[1]) # least free-boundary eigenmode energy: negative is unstable for (s, sing) in enumerate(ffs.surfaces) - println("q = $(sing.m[1])/$(sing.n[1]) at ψ_N = $(sing.psifac): Δ′ = $(real(ffs.delta_prime.matrix[s, s]))") + println("q = $(sing.m[1])/$(sing.n[1]) at ψ_N = $(sing.psifac): Δ′ = $(real(ffs.delta_prime.freeboundary.Delta[s, s]))") end ``` @@ -103,7 +103,7 @@ The number of Julia threads is set at startup (`-t N` or `JULIA_NUM_THREADS`); t | `SingularSurfaces/rational_psi`, `rational_q`, `rational_m`, `rational_n`, `dqdpsi` | the rational surfaces of the domain, core to edge | | `SingularSurfaces/D_I`, `E`, `F`, `G`, `H`, `K`, `M` | the Mercier index and the Glasser–Greene–Johnson layer coefficients at each surface | | `SingularSurfaces/ca_left`, `ca_right` | the asymptotic (large, small) coefficients on either side of each surface, from ideal crossings | -| `SingularSurfaces/Delta_prime_matrix` | the Δ′ matrix in the PEST-3 convention, from `"Riccati"` or `"Galerkin"`; `Delta_prime_raw` is the side-major (left, right) matrix it is folded from, `Delta_coil` the edge coil-response columns | +| `SingularSurfaces/Delta_prime_matrix` | the free-boundary Δ′ matrix in the PEST-3 convention, from `"Riccati"` or `"Galerkin"`, with `pest3_A`, `pest3_B`, `pest3_Gamma`; `Delta_prime_raw_freeboundary` is the side-major (left, right) matrix they are folded from. The same set with a `_fixedboundary` suffix is the Δ′ with a conformal wall on the plasma, `Delta_EF` the error-field response, and `edge_response` the edge momentum of each fixed-boundary big-solution column (what matching needs for the matched `W_plasma`) | | `SingularSurfaces/Kinetic/` | the surfaces of a kinetic solve | The energy matrices are power-normalized: the eigenmode energies are the generalized eigenvalues of the pencil ``W v = \lambda N v``, with ``N`` the surface-norm matrix ``\xi^\dagger N \xi = \langle|\xi|^2\rangle`` built from the Jacobian Fourier coefficients, so each is an energy per unit surface-averaged squared displacement rather than joules, and the spectrum is invariant to the working coordinate because ``W`` and ``N`` transform by the same congruence. diff --git a/docs/src/workflow.md b/docs/src/workflow.md index 38b4cdd8b..8a491300b 100644 --- a/docs/src/workflow.md +++ b/docs/src/workflow.md @@ -29,7 +29,7 @@ Every stage after `Equilibrium` reads the force-free states, and `KineticForces` | What are the tearing growth rate and critical resonant field? | `[Tearing]` (+ `[SLAYER]` or `[GGJ]`) | `Tearing/Roots/gamma`, `Tearing/Roots/omega`, `Tearing/CriticalResonantField/` | | What does a resistive inner layer do to the response? | `[Match]` (+ `[GGJ]`) | `SingularSurfaces/Match/` | -`[PerturbedEquilibrium]` and the stages built on it need a `[ForceFreeStates]` solve with `integrator = "Forward"` and no `[Wall]`, while the Δ′ matrix, `[Tearing]` and wall stability need `"Riccati"` or `"Galerkin"` or a `[Wall]`, so one deck cannot do both: run the response and the Δ′ or wall cases as separate decks (`DIIID-like_ideal_example` and `DIIID-like_riccati_deltaprime_example` are such a pair). +`[PerturbedEquilibrium]` and the stages built on it need a `[ForceFreeStates]` solve with `integrator = "Forward"`, or `"Galerkin"` with `error_field_response = true` and `vac_flag = true`, and no `[Wall]`, while the Δ′ matrix, `[Tearing]` and wall stability need `"Riccati"` or `"Galerkin"` or a `[Wall]`. Only a `"Galerkin"` deck without a `[Wall]` does both; otherwise run the response and the Δ′ or wall cases as separate decks (`DIIID-like_ideal_example` and `DIIID-like_riccati_deltaprime_example` are such a pair). ### Headline outputs @@ -68,7 +68,7 @@ Page: [Local Stability](ballooning.md). ### Match -Matches the outer-region Δ′ to a resistive inner layer (`inner_model = "GGJ"`, the only model, with its options in an optional `[GGJ]` section) and rebuilds the force-free states with the resulting layer response, writing the matched solution and the penetrated resonant fields to `SingularSurfaces/Match/`; `ideal = true` keeps the shielded limit for comparison. +Matches the outer-region Δ′ to a resistive inner layer (`inner_model = "GGJ"`, the only model, with its options in an optional `[GGJ]` section) and rebuilds the force-free states with the resulting layer response, writing the matched solution and the penetrated resonant fields to `SingularSurfaces/Match/`; the Galerkin solve without `[Match]` is the shielded limit for comparison. Outputs: `SingularSurfaces/Match/`. Page: [Tearing](inner_layer.md). @@ -76,7 +76,7 @@ Page: [Tearing](inner_layer.md). Computes the plasma response to an external non-axisymmetric field: `plasma_response` builds the permeability, inductance and reluctance operators on the control surface from the force-free states and the vacuum matrices, and `perturbed_equilibrium` applies them to the drive that `ForcingTerms` supplies, either a normal-field spectrum from a file or the Biot–Savart field of `[[ForcingTerms.coil_set]]` coils. The resonant fields, island half-widths and Chirikov parameter at each rational surface follow from the singular coupling (`compute_singular_coupling`), and `regularization_width` controls the regularization of the ideal response there. -The response is built only from a `"Forward"` solve without a wall: a `"Riccati"` result (the default) makes `plasma_response` skip the operators with a warning, and a `"Galerkin"` or matched result or a wall is an error. +The response is built from a `"Forward"` solve, or a `"Galerkin"` solve with `error_field_response = true` and `vac_flag = true`, without a wall: a `"Riccati"` result (the default) makes `plasma_response` skip the operators with a warning, and any other `"Galerkin"` result, a matched result or a wall is an error. Outputs: `PerturbedEquilibrium/ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/`. Key references: [Park et al. (2007a)](citations.md#PerturbedEquilibrium), [Park et al. (2009)](citations.md#PerturbedEquilibrium). Pages: [Forcing Terms](forcing_terms.md), [Perturbed Equilibrium](perturbed_equilibrium.md). @@ -100,7 +100,7 @@ Page: [Error Fields](error_fields.md). Solves the tearing dispersion relation of each rational surface by matching the outer-region Δ′ to an inner-layer model: the two-fluid slab `SLAYER` (`[SLAYER]`) or the toroidal resistive-MHD `GGJ` (`[GGJ]`), selected by `[Tearing] inner_model`. `[Tearing.search]` sets the root search in the complex frequency plane, and `[Tearing.CriticalResonantField]` adds the resonant field at which the layer's shielding is overcome. -Needs the Δ′ matrix of `"Riccati"` with `vac_flag = true` or `"Galerkin"` with `[Galerkin] rpec_flag = true`, and the kinetic profiles of `[Equilibrium] kinetic_file` (`SLAYER` reads its densities, temperatures, rotation and diffusivities there). +Needs the free-boundary Δ′ matrix of `"Riccati"` or `"Galerkin"` with `vac_flag = true` (or the fixed-boundary one with `fixed_boundary = true`), and the kinetic profiles of `[Equilibrium] kinetic_file` (`SLAYER` reads its densities, temperatures, rotation and diffusivities there). Outputs: `Tearing/PerSurface/`, `Roots/`, `LayerWidths/`, `Diagnostics/`, `CriticalResonantField/`. Key references: [Glasser, Wang and Park (2016)](citations.md#Tearing), [Park (2022)](citations.md#Tearing). Page: [Tearing](inner_layer.md). @@ -127,7 +127,7 @@ alpha_boundary_scan = true [ForceFreeStates] nn = 1 vac_flag = true -integrator = "Forward" # "Forward": [PerturbedEquilibrium] needs the dense ξ solution and no [Wall] ("Riccati" is the default) +integrator = "Forward" # "Forward": [PerturbedEquilibrium] needs a dense ξ solution and no [Wall] ("Riccati" is the default) [ForcingTerms] forcing_data_format = "coil" @@ -191,7 +191,7 @@ From Julia, the same pipeline is `GeneralizedPerturbedEquilibrium.main(["example | `Solovev_ideal_example_3D` | the same plasma with the vacuum response on a toroidal grid | | `Solovev_ideal_example_multi_n` | the same plasma swept over n = 1 and 2 in one run, with single-n cross-checks | | `LAR_reversed_shear_example` | a `[LargeAspectRatio]` equilibrium with reversed shear, q = 2 crossed twice | -| `LAR_ideal_match_test` | a `[TJAnalytic]` equilibrium with `[Match] ideal = true`, the shielded reference for the resistive case | +| `LAR_ideal_match_test` | a `[TJAnalytic]` equilibrium solved by Galerkin without `[Match]`, the shielded reference for the resistive case | | `LAR_resistive_match_test` | the same equilibrium with resistive GGJ matching | | `LAR_beta_scan`, `LAR_epsilon_scan` | `run_scan.jl` sweeps of the `[TJAnalytic]` pressure factor and inverse aspect ratio | | `TJ_ideal_example` | a `[TJAnalytic]` equilibrium: n = 1 stability and a file-driven perturbed equilibrium | diff --git a/examples/DIIID-like_gal_resistive_example/gpec.toml b/examples/DIIID-like_gal_resistive_example/gpec.toml index 0abc47dc1..909d79cc4 100644 --- a/examples/DIIID-like_gal_resistive_example/gpec.toml +++ b/examples/DIIID-like_gal_resistive_example/gpec.toml @@ -1,6 +1,6 @@ # DIII-D-like H-mode resistive-MHD example: outer-region singular Galerkin Δ′ solve. # n=1 tearing stability on the bundled geqdsk: the gal solver computes the Δ′ -# matching matrix (PEST-3 convention) plus rpec coil-response columns (delta_coil). +# matching matrix (PEST-3 convention) plus error-field (EF) columns (Delta_EF). # The core_edge_packed radial grid is densest at the axis and edge; psihigh = 0.993 # stays clear of the separatrix. No forcing and no PerturbedEquilibrium stage. [Run] @@ -53,8 +53,8 @@ rational_offset_fraction = 0.2 # Truncate integration at (last_rational_q + rati frobenius_order = 6 # Order of singular layer (Frobenius) expansion at rational surfaces -# Outer-region singular Galerkin Δ′ solver with rpec coil columns. -# This case exercises the gal solve AND the rpec coil/edge path (delta_coil) for regression tracking. +# Outer-region singular Galerkin Δ′ solver with error-field (EF) columns. +# This case exercises the gal solve AND the EF/edge path (Delta_EF) for regression tracking. [Galerkin] nx = 256 # Hermite-cubic elements per interval between singular surfaces nq = 6 # Gauss-Lobatto quadrature order per element @@ -67,7 +67,7 @@ tol = 1e-10 # Resonant-cell adaptive-quadrature (QuadGK) tole gnstep = 5000 # Maximum resonant-integration steps dx1dx2_flag = true # Enable the special dx1/dx2 element sizing for resonant/extension cells frobenius_order_ceiling = true # Auto-raise the series order by ceil(2·Re(α)) on high-Mercier-index surfaces -rpec_flag = true # Append mpert coil-response columns (unit edge sources) to the solve → delta_coil block +error_field_response = true # Add the error-field (EF) columns, a unit edge displacement per mode: gives Delta_EF, the ideal solution and wp # Driven inner-layer matching on the finished solve. Per-surface resistive inputs (core→edge) for a DIII-D-like case: # η=8e-8, ρ=3.3e-7 kg/m³ (ni≈1e20 m^-3), rotation=1 Hz, Γ=5/3. 4 surfaces. diff --git a/examples/LAR_ideal_match_test/gpec.toml b/examples/LAR_ideal_match_test/gpec.toml index d266d6af9..5717eaa9e 100644 --- a/examples/LAR_ideal_match_test/gpec.toml +++ b/examples/LAR_ideal_match_test/gpec.toml @@ -1,9 +1,8 @@ # LAR ideal-limit reference: large-aspect-ratio analytic TJ equilibrium, epsilon = a/R0 = 0.2, n = 1. # -# The ideal-limit companion to LAR_resistive_match_test, identical except that [Match] ideal = true skips -# the resistive inner layer entirely and builds the bare (perfectly shielded) coil columns. Running -# the two side by side isolates what the inner layer contributes on this geometry. The per-surface -# resistivity, density and rotation inputs are omitted here because the ideal limit never reads them. +# The ideal-limit companion to LAR_resistive_match_test, identical except that it has no [Match] section: +# the Galerkin solve's EF columns are already the perfectly shielded solution. Running the two side by +# side isolates what the inner layer contributes on this geometry. # This run terminates after the stability stage. [Run] @@ -63,7 +62,7 @@ integrator = "Galerkin" # Integration formalism, whose numerical settings are it truncate_at_rational_offset = false # Keep psilim at psihigh (do not truncate at last_rational_q + rational_offset_fraction) rational_offset_fraction = 0.2 # Truncate integration at (last_rational_q + rational_offset_fraction) / n (only used when truncate_at_rational_offset = true) -# Outer-region singular Galerkin Delta-prime solver with rpec coil columns. +# Outer-region singular Galerkin Delta-prime solver with error-field (EF) columns. [Galerkin] nx = 256 # Hermite-cubic elements per interval between singular surfaces nq = 6 # Gauss-Lobatto quadrature order per element @@ -76,8 +75,4 @@ tol = 1e-10 # Resonant-cell adaptive-quadrature (QuadGK) tole gnstep = 5000 # Maximum resonant-integration steps dx1dx2_flag = true # Enable the special dx1/dx2 element sizing for resonant/extension cells frobenius_order_ceiling = true # Auto-raise the series order by ceil(2*Re(alpha)) on high-Mercier-index surfaces -rpec_flag = true # Append mpert coil-response columns (unit edge sources) to the solve, giving the delta_coil block - -# Driven inner-layer matching on the finished solve in the ideal limit: the inner layer is skipped entirely. -[Match] -ideal = true # true = build the IDEAL solution (skip the inner layer, bare coil columns) +error_field_response = true # Add the error-field (EF) columns, a unit edge displacement per mode: gives Delta_EF, the ideal solution and wp diff --git a/examples/LAR_resistive_match_test/gpec.toml b/examples/LAR_resistive_match_test/gpec.toml index 248f1cd14..4fb3f0da0 100644 --- a/examples/LAR_resistive_match_test/gpec.toml +++ b/examples/LAR_resistive_match_test/gpec.toml @@ -64,7 +64,7 @@ integrator = "Galerkin" # Integration formalism, whose numerical settings are it truncate_at_rational_offset = false # Keep psilim at psihigh (do not truncate at last_rational_q + rational_offset_fraction) rational_offset_fraction = 0.2 # Truncate integration at (last_rational_q + rational_offset_fraction) / n (only used when truncate_at_rational_offset = true) -# Outer-region singular Galerkin Delta-prime solver with rpec coil columns. +# Outer-region singular Galerkin Delta-prime solver with error-field (EF) columns. [Galerkin] nx = 256 # Hermite-cubic elements per interval between singular surfaces nq = 6 # Gauss-Lobatto quadrature order per element @@ -77,7 +77,7 @@ tol = 1e-10 # Resonant-cell adaptive-quadrature (QuadGK) tole gnstep = 5000 # Maximum resonant-integration steps dx1dx2_flag = true # Enable the special dx1/dx2 element sizing for resonant/extension cells frobenius_order_ceiling = true # Auto-raise the series order by ceil(2*Re(alpha)) on high-Mercier-index surfaces -rpec_flag = true # Append mpert coil-response columns (unit edge sources) to the solve, giving the delta_coil block +error_field_response = true # Add the error-field (EF) columns, a unit edge displacement per mode: gives Delta_EF, the ideal solution and wp # Driven inner-layer matching on the finished solve. Per-surface resistive inputs (core to edge) for the two # rational surfaces q = 2 and q = 3: eta=8e-8, rho=3.3e-7 kg/m^3, rotation=1 Hz, Gamma=5/3. diff --git a/ext/GPECMakieExt/forcefreestates.jl b/ext/GPECMakieExt/forcefreestates.jl index cb7d33a75..e306f9159 100644 --- a/ext/GPECMakieExt/forcefreestates.jl +++ b/ext/GPECMakieExt/forcefreestates.jl @@ -214,16 +214,16 @@ function FFS._plot_delta_prime_raw(dp; part=:real) end end -function FFS._plot_delta_coil(dp; part=:real) +function FFS._plot_delta_ef(dp; part=:real) return themed() do ns = length(dp.labels) - me = dp.mlow .+ (0:(size(dp.coil, 1)-1)) + me = dp.mlow .+ (0:(size(dp.ef, 1)-1)) fig = Figure(; size=(900, 560)) ax = _surface_side_axes!(fig[1, 1], dp.labels; xgroups=false, xlabel="applied edge harmonic m", - title=rich(_part_label(part, sub("Δ", "coil")), ": surface response to each edge harmonic"), + title=rich(_part_label(part, sub("Δ", "EF")), ": surface response to each edge harmonic"), limits=((me[1] - 0.5, me[end] + 0.5), (0.5, 2ns + 0.5))) - hm = _complex_heatmap!(ax, me, 1:(2ns), _part(dp.coil, part), part) - Colorbar(fig[1, 2], hm; label=_part_label(part, sub("Δ", "coil"))) + hm = _complex_heatmap!(ax, me, 1:(2ns), _part(dp.ef, part), part) + Colorbar(fig[1, 2], hm; label=_part_label(part, sub("Δ", "EF"))) fig end end diff --git a/regression-harness/cases/diiid_n1_riccati.toml b/regression-harness/cases/diiid_n1_riccati.toml index 4c2642145..8718787ce 100644 --- a/regression-harness/cases/diiid_n1_riccati.toml +++ b/regression-harness/cases/diiid_n1_riccati.toml @@ -21,18 +21,27 @@ noise_threshold = 1e-8 order = 10 [quantities.delta_prime_raw] -h5path = "SingularSurfaces/Delta_prime_raw" +h5path = "SingularSurfaces/Delta_prime_raw_freeboundary" type = "complex_matrix" extract = "all_complex" label = "delta prime (raw side-major)" noise_threshold = 1e-8 order = 11 -[quantities.delta_coil] -h5path = "SingularSurfaces/Delta_coil" +# Fixed-boundary (conformal wall on the plasma) PEST-3 Δ′ diagonal +[quantities.delta_prime_fixedboundary_diag] +h5path = "SingularSurfaces/Delta_prime_matrix_fixedboundary" +type = "complex_matrix" +extract = "diagonal_complex" +label = "delta prime fixed boundary (BVP diagonal)" +noise_threshold = 1e-8 +order = 13 + +[quantities.delta_ef] +h5path = "SingularSurfaces/Delta_EF" type = "complex_matrix" extract = "all_complex" -label = "edge coil response delta_coil" +label = "edge coil response delta_ef" noise_threshold = 1e-8 order = 12 diff --git a/regression-harness/cases/gal_resistive_diiid.toml b/regression-harness/cases/gal_resistive_diiid.toml index 691b9c949..8b8dceecd 100644 --- a/regression-harness/cases/gal_resistive_diiid.toml +++ b/regression-harness/cases/gal_resistive_diiid.toml @@ -1,12 +1,12 @@ # Regression case: DIII-D-like n=1 outer-region singular Galerkin Δ′ solve. # Tracks the resonant-surface inventory, the PEST-3 Δ matching-matrix diagonal # (the physics-meaningful tearing Δ′), norms of the full outer Δ′ matrix and the -# rpec coil-response (Δ_coil) block, and the inner-layer matching outputs. +# error-field (Δ_EF) block, and the inner-layer matching outputs. # Each [quantities.*] block names an HDF5 path in the run output, how to extract it, # and the noise floor below which a difference is treated as zero. [case] name = "gal_resistive_diiid" -description = "DIII-D-like, n=1, RDCON outer-region Galerkin Δ′ with rpec coil columns (delta_coil)" +description = "DIII-D-like, n=1, RDCON outer-region Galerkin Δ′ with error-field (EF) columns (Delta_EF)" example_dir = "examples/DIIID-like_gal_resistive_example" # Number of resonant surfaces in the gal domain @@ -37,19 +37,28 @@ order = 20 # Raw outer Δ′ matrix (2·msing × 2·msing) — Frobenius norm catches any element drift [quantities.gal_delta_norm] -h5path = "SingularSurfaces/Delta_prime_raw" +h5path = "SingularSurfaces/Delta_prime_raw_freeboundary" type = "complex_matrix" extract = "norm" label = "||gal Δ′ matrix||" noise_threshold = 1e-6 order = 21 -# rpec coil-response block (edge mode × surface-side) — Frobenius norm -[quantities.gal_delta_coil_norm] -h5path = "SingularSurfaces/Delta_coil" +# Fixed-boundary (conformal wall on the plasma) PEST-3 Δ′ diagonal +[quantities.gal_pest3_delta_fixedboundary_diag] +h5path = "SingularSurfaces/Delta_prime_matrix_fixedboundary" +type = "complex_matrix" +extract = "diagonal_complex" +label = "gal PEST3 Δ diagonal, fixed boundary" +noise_threshold = 1e-6 +order = 24 + +# error-field (Δ_EF) block (edge mode × surface-side) — Frobenius norm +[quantities.gal_delta_ef_norm] +h5path = "SingularSurfaces/Delta_EF" type = "complex_matrix" extract = "norm" -label = "||gal Δ_coil block||" +label = "||gal Δ_EF block||" noise_threshold = 1e-6 order = 22 @@ -85,7 +94,7 @@ label = "||gal inner-layer Δ||" noise_threshold = 1e-6 order = 41 -# Outer-region matched coefficients cout (2·msing × mcoil) — the 4·msing matching assembly + solve. +# Outer-region matched coefficients cout (2·msing × nef) — the 4·msing matching assembly + solve. [quantities.gal_match_cout_norm] h5path = "SingularSurfaces/Match/cout" type = "complex_matrix" @@ -94,7 +103,7 @@ label = "||gal match cout||" noise_threshold = 1e-6 order = 42 -# Penetrated resonant field at each layer center (msing × mcoil) — the headline matching output. +# Penetrated resonant field at each layer center (msing × nef) — the headline matching output. [quantities.gal_match_bpen_norm] h5path = "SingularSurfaces/Match/bpen" type = "complex_matrix" diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl index cb38650ac..b2bd09f3d 100644 --- a/src/Analysis/ForceFreeStates.jl +++ b/src/Analysis/ForceFreeStates.jl @@ -72,23 +72,27 @@ function _edge_scan_data(h5path::AbstractString) end # The Δ′ family with "m/n" labels for the surfaces it spans: Galerkin solves across its in-band subset of the rational surfaces. -# `coil` is (edge harmonic × surface side), empty when the run computed no coil response. +# The free-boundary Δ′ and raw D′, or the fixed-boundary ones when the run has no vacuum. `ef` is Δ_EF as +# (edge harmonic × surface side), empty when the run computed none. function _delta_prime_data(r::ForceFreeStatesResult) dp = r.delta_prime - k = size(dp.matrix, 1) + blocks, raw = dp.freeboundary === nothing ? (dp.fixedboundary, dp.raw_fixedboundary) : (dp.freeboundary, dp.raw_freeboundary) + k = size(blocks.Delta, 1) g = r.galerkin m, n = g !== nothing && length(g.sing_m) == k ? (g.sing_m, g.sing_n) : ([s.m[1] for s in r.surfaces[1:k]], [s.n[1] for s in r.surfaces[1:k]]) - return (; matrix=dp.matrix, raw=dp.raw, coil=permutedims(dp.coil), labels=["$mi/$ni" for (mi, ni) in zip(m, n)], mlow=r.mlow) + return (; matrix=blocks.Delta, raw, ef=permutedims(dp.ef), labels=["$mi/$ni" for (mi, ni) in zip(m, n)], mlow=r.mlow) end function _delta_prime_data(h5path::AbstractString) h5open(h5path, "r") do fid - matrix = read(fid["SingularSurfaces/Delta_prime_matrix"]) + suffix = haskey(fid, "SingularSurfaces/Delta_prime_matrix") ? "" : "_fixedboundary" + matrix = read(fid["SingularSurfaces/Delta_prime_matrix$suffix"]) k = size(matrix, 1) group = haskey(fid, "$_GAL/rational_m") && length(fid["$_GAL/rational_m"]) == k ? _GAL : "SingularSurfaces" m, n = vec(read(fid["$group/rational_m"]))[1:k], vec(read(fid["$group/rational_n"]))[1:k] - coil = haskey(fid, "SingularSurfaces/Delta_coil") ? read(fid["SingularSurfaces/Delta_coil"]) : zeros(ComplexF64, 0, 2k) - (; matrix, raw=read(fid["SingularSurfaces/Delta_prime_raw"]), coil, labels=["$mi/$ni" for (mi, ni) in zip(m, n)], mlow=read(fid["Info/mlow"])) + ef = haskey(fid, "SingularSurfaces/Delta_EF") ? read(fid["SingularSurfaces/Delta_EF"]) : zeros(ComplexF64, 0, 2k) + raw = read(fid["SingularSurfaces/Delta_prime_raw_$(isempty(suffix) ? "freeboundary" : "fixedboundary")"]) + (; matrix, raw, ef, labels=["$mi/$ni" for (mi, ni) in zip(m, n)], mlow=read(fid["Info/mlow"])) end end @@ -117,7 +121,7 @@ function _plot_eigenvalues end function _plot_edge_stability_scan end function _plot_delta_prime end function _plot_delta_prime_raw end -function _plot_delta_coil end +function _plot_delta_ef end function _plot_cond_fbar end function _plot_ballooning_alpha_boundary end function _plot_ballooning_alpha_boundaries end @@ -298,14 +302,14 @@ function plot_delta_prime_raw(src::Source; part=:real, save_path=nothing) end """ - plot_delta_coil(src; part=:real, save_path=nothing) + plot_delta_ef(src; part=:real, save_path=nothing) -Plot Δ_coil, the response at each side of each rational surface to a unit source on each edge poloidal harmonic, -as a heatmap. +Plot Δ_EF, the small-solution coefficient at each side of each rational surface for a unit edge displacement in each +edge poloidal harmonic, as a heatmap. ### Arguments - - `src`: A Riccati `ForceFreeStatesResult` or a Galerkin one solved with `rpec_flag = true`, or the path of a `gpec.h5` it was saved to + - `src`: A Riccati or Galerkin `ForceFreeStatesResult` solved with `error_field_response = true`, or the path of a `gpec.h5` it was saved to ### Keyword arguments @@ -316,10 +320,10 @@ as a heatmap. A Makie `Figure`. """ -function plot_delta_coil(src::Source; part=:real, save_path=nothing) +function plot_delta_ef(src::Source; part=:real, save_path=nothing) _check_part(part) - _require(plot_delta_coil, src) - return _draw(_plot_delta_coil, _delta_prime_data(src); part, save_path) + _require(plot_delta_ef, src) + return _draw(_plot_delta_ef, _delta_prime_data(src); part, save_path) end """ @@ -436,7 +440,7 @@ register!(plot_energy_eigenvectors, (ForceFreeStatesResult, AbstractString), (:f register!(plot_edge_stability_scan, (ForceFreeStatesResult, AbstractString), (:edge_scan,), "edge energies against the truncation q") register!(plot_delta_prime, (ForceFreeStatesResult, AbstractString), (:delta_prime,), "Δ′ matrix between rational surfaces") register!(plot_delta_prime_raw, (ForceFreeStatesResult, AbstractString), (:delta_prime,), "raw outer-region matching matrix D′") -register!(plot_delta_coil, (ForceFreeStatesResult, AbstractString), (:delta_coil,), "Δ_coil: surface response to each edge harmonic") +register!(plot_delta_ef, (ForceFreeStatesResult, AbstractString), (:delta_ef,), "Δ_EF: surface response to each edge harmonic") register!(plot_cond_fbar, (ForceFreeStatesResult, AbstractString), (:kinetic_scan,), "cond(F̄) scan for kinetic singular surfaces") register!(plot_ballooning_alpha_boundary, (LocalStabilityResult, AbstractString), (:ballooning_scan,), "ballooning α against its first stability boundary") diff --git a/src/Analysis/Registry.jl b/src/Analysis/Registry.jl index df1b6a852..aeb60ce80 100644 --- a/src/Analysis/Registry.jl +++ b/src/Analysis/Registry.jl @@ -43,12 +43,12 @@ What each capability needs, worded as what to change in the run to get it. const REASONS = Dict{Symbol,String}( :equilibrium => "needs an equilibrium", :singular_surfaces => "needs rational surfaces in the integration domain", - :solution_profiles => "needs stored ξ profiles: use the Forward integrator or a matched Galerkin solve (Riccati and unmatched Galerkin store none)", + :solution_profiles => "needs stored ξ profiles: use the Forward integrator, or Galerkin with error_field_response = true (Riccati stores none)", :fixed_boundary_criterion => "needs the critical determinant scan: use the Forward or Riccati integrator (Galerkin computes none)", - :free_boundary => "needs the free-boundary energies: run with vac_flag = true and the Forward or Riccati integrator", + :free_boundary => "needs the free-boundary energies: run with vac_flag = true (Galerkin also needs error_field_response = true)", :edge_scan => "needs the edge stability scan: set dW_edge_scan_start below psilim", :delta_prime => "needs the Δ′ matrix: use the Riccati or Galerkin integrator (Forward computes none)", - :delta_coil => "needs Δ_coil: use the Riccati integrator, or Galerkin with rpec_flag = true", + :delta_ef => "needs Δ_EF: use Galerkin, or Riccati with vac_flag = true, with error_field_response = true", :kinetic_scan => "needs the kinetic singular-surface scan: run with a kinetic source and singfac_min > 0", :ballooning_scan => "needs the ballooning scan: add a [LocalStability] section", :singular_coupling => "needs the singular coupling: run PerturbedEquilibrium with compute_singular_coupling = true", @@ -87,8 +87,8 @@ function capabilities(r::ForceFreeStatesResult) diag === nothing || isempty(diag.crit_store) || push!(caps, :fixed_boundary_criterion) diag === nothing || isempty(diag.edge_scan.psi) || push!(caps, :edge_scan) r.free_boundary === nothing || isempty(r.free_boundary.et) || push!(caps, :free_boundary) - r.delta_prime === nothing || isempty(r.delta_prime.matrix) || push!(caps, :delta_prime) - r.delta_prime === nothing || isempty(r.delta_prime.coil) || push!(caps, :delta_coil) + r.delta_prime === nothing || push!(caps, :delta_prime) + r.delta_prime === nothing || isempty(r.delta_prime.ef) || push!(caps, :delta_ef) isempty(r.kinetic.scan_psi) || push!(caps, :kinetic_scan) return caps end @@ -131,8 +131,8 @@ const _H5_CAPABILITIES = ( :fixed_boundary_criterion => "ForceFreeStates/Solutions/ForwardIntegration/crit", :edge_scan => "ForceFreeStates/EdgeScan/psi", :free_boundary => "ForceFreeStates/FreeBoundaryStability/eigenmode_energies", - :delta_prime => "SingularSurfaces/Delta_prime_matrix", - :delta_coil => "SingularSurfaces/Delta_coil", + :delta_prime => ("SingularSurfaces/Delta_prime_matrix", "SingularSurfaces/Delta_prime_matrix_fixedboundary"), + :delta_ef => "SingularSurfaces/Delta_EF", :kinetic_scan => "SingularSurfaces/Kinetic/scan_psi", :ballooning_scan => "LocalStability/alpha", :singular_coupling => "PerturbedEquilibrium/SingularCoupling/resonant_area_weighted_field", diff --git a/src/Deck.jl b/src/Deck.jl index ea0f4b17c..1a91b0a00 100644 --- a/src/Deck.jl +++ b/src/Deck.jl @@ -32,13 +32,14 @@ const _MOVED_KEYS = Dict{String,Dict{String,String}}( "extended_precision_bvp" => "moved to [Riccati]", "nchunks" => "moved to [Riccati]", "gal_match_flag" => "add a [Match] section to match the inner layer", - "gal_ideal_flag" => "moved to [Match] as ideal", + "gal_ideal_flag" => "removed: a Galerkin solve with [Galerkin] error_field_response = true is the ideal solution", "gal_inner_solver" => "moved to [GGJ] as solver", (k => "moved to the integrator's section, [Forward] or [Riccati]" for k in _ODE_DECK_KEYS)..., ("gal_$k" => "moved to [Match] as $k" for k in ("eta", "rho", "rotation", "gamma"))..., ("gal_inner_$k" => "moved to [GGJ] as $k" for k in ("xfac", "nx", "nq", "cutoff", "kmax"))..., ("gal_$k" => "moved to [Galerkin] as $k" for k in ("nx", "nq", "pfac", "dx0", "dx1", "dx2", "cutoff", "tol", "gnstep", - "dx1dx2_flag", "rpec_flag", "edge_onesided"))..., + "dx1dx2_flag", "edge_onesided"))..., + "gal_rpec_flag" => "moved to [Galerkin] as error_field_response", "gal_sing_order_ceiling" => "moved to [Galerkin] as frobenius_order_ceiling"), "LocalStability" => Dict{String,String}( k => "moved to [Ballooning]" for k in ("ode_tol", "n_scan", "max_alpha_scale", "alpha_tol", "psi_window")), @@ -344,7 +345,7 @@ end """ function match_from_deck(inputs::AbstractDict) haskey(inputs, "Match") || return nothing - fields = Dict{Symbol,Any}(:eta => Any, :rho => Any, :rotation => Any, :gamma => Any, :ideal => Bool, :inner_model => String) + fields = Dict{Symbol,Any}(:eta => Any, :rho => Any, :rotation => Any, :gamma => Any, :inner_model => String) problem = section_keywords(inputs["Match"], "Match", fields) _selector(inputs["Match"], "Match", "inner_model", ("GGJ",), "GGJ") delete!(problem, :inner_model) @@ -467,8 +468,8 @@ function parse_deck(inputs::AbstractDict, path::AbstractString) if match !== nothing (; keywords, algorithm) = force_balance first(keywords[:nn]) == last(keywords[:nn]) || error("[Match] is single-n: set [ForceFreeStates] nn to one toroidal mode number") - (algorithm isa Galerkin && algorithm.rpec_flag) || (algorithm isa Riccati && get(keywords, :vac_flag, false)) || - error("[Match] needs Δ′: use integrator = \"Galerkin\" with [Galerkin] rpec_flag = true, or \"Riccati\" with vac_flag = true") + hasproperty(algorithm, :error_field_response) && algorithm.error_field_response || + error("[Match] needs Δ_EF: use the \"Galerkin\" or \"Riccati\" integrator with error_field_response = true") end tearing = tearing_from_deck(inputs) perturbed_equilibrium = perturbed_equilibrium_from_deck(inputs) diff --git a/src/Display.jl b/src/Display.jl index 662cbbd47..d8aa68fbb 100644 --- a/src/Display.jl +++ b/src/Display.jl @@ -134,8 +134,9 @@ function Base.show(io::IO, fb::FreeBoundaryResult) end function Base.show(io::IO, d::DeltaPrimeData) - source = d.A === nothing ? "Riccati" : "Galerkin" - return print(io, "DeltaPrimeData(", source, ", ", _shape(d.matrix), ")") + boundaries = join([b for (b, blocks) in (("free", d.freeboundary), ("fixed", d.fixedboundary)) if blocks !== nothing], " + ") + msing = size(something(d.freeboundary, d.fixedboundary).Delta, 1) + return print(io, "DeltaPrimeData(", boundaries, " boundary, ", _plural(msing, "surface"), isempty(d.ef) ? "" : ", Δ_EF", ")") end function Base.show(io::IO, g::GalerkinResult) @@ -204,7 +205,7 @@ end # m/n of each row of the Δ′ matrix: Galerkin solves across its in-band subset of the surfaces. function _delta_prime_labels(r::ForceFreeStates.ForceFreeStatesResult) - k = minimum(size(r.delta_prime.matrix)) + k = size(something(r.delta_prime.freeboundary, r.delta_prime.fixedboundary).Delta, 1) g = r.galerkin g !== nothing && length(g.sing_m) == k && return ["$(m)/$(n)" for (m, n) in zip(g.sing_m, g.sing_n)] length(r.surfaces) == k && return _surface_label.(r.surfaces) @@ -221,8 +222,9 @@ function Base.show(io::IO, ::MIME"text/plain", r::ForceFreeStates.ForceFreeState _row(io, "Fixed boundary", _fixed_boundary_note(r)) _row(io, "Free boundary", _free_boundary_note(r)) if r.delta_prime !== nothing - dp = r.delta_prime.matrix - _row(io, "Δ′_ii (ψ_N)", _list(["$l: $(_fmt(real(dp[i, i])))" for (i, l) in enumerate(_delta_prime_labels(r))]; max=4)) + free = r.delta_prime.freeboundary !== nothing + dp = something(r.delta_prime.freeboundary, r.delta_prime.fixedboundary).Delta + _row(io, free ? "Δ′_ii (ψ_N)" : "Δ′_ii fixed boundary (ψ_N)", _list(["$l: $(_fmt(real(dp[i, i])))" for (i, l) in enumerate(_delta_prime_labels(r))]; max=4)) end _row(io, "Solution", _solution_note(r)) r.galerkin === nothing || _row(io, "Galerkin", "$(_plural(r.galerkin.msing, "surface")), $(r.galerkin.solver) factorization") @@ -241,7 +243,7 @@ function Base.show(io::IO, ::MIME"text/plain", m::MatchResult) print(io, "MatchResult ", InnerLayer.model_name(m.model)) _row(io, "Surfaces", "$(size(m.deltar, 1)) matched, $(size(m.cout, 2)) edge drives") _row(io, "Residual", _fmt(m.residual)) - _row(io, "Forced γ [rad/s]", isempty(m.rpec_eig) ? "none" : _list(_fmt.(m.rpec_eig); max=3)) + _row(io, "Forced γ [rad/s]", isempty(m.forced_eigenvalue) ? "none" : _list(_fmt.(m.forced_eigenvalue); max=3)) outer = isempty(m.xi) ? "not computed" : "$(_shape(m.xi)) (one column per edge drive)" _row(io, "Stored", "matched outer ξ $outer, inner profiles $(isempty(m.inner_psi) ? "not computed" : "on $(length(m.inner_psi)) surfaces")") return nothing diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index 1cb794f65..433cb4ec8 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/src/ForceFreeStates/CoreTypes.jl @@ -60,25 +60,14 @@ Solve-time scratch state for a force-balance solve: the resolved mode space, sin wall_settings::Union{Nothing,Vacuum.WallSettings} = nothing verbose::Bool = false """ - Inter-surface Δ' matrix, `msing × msing`, PEST3 convention: the odd-parity projection of - `delta_prime_raw`, set by `compute_delta_prime_matrix!` (Riccati path only). + Riccati Δ' BVP results (`compute_delta_prime_matrix!`): raw side-major `2msing × 2msing` D' with the vacuum edge + and with ξ(edge) = 0, Δ_EF (`2msing × numpert_total`), and the edge momentum of the ξ(edge) = 0 solutions + (`numpert_total × 2msing`). Empty when not computed. """ - delta_prime_matrix::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0) - - """ - Edge coil-response matrix, `2msing × numpert_total`: column k is the resonant small-solution - response at each surface side to a unit source on edge mode k (Eq. 37 rpec edge condition, - `_solve_bvp_edge_coil`). Empty unless the S-axis vacuum-edge BVP was assembled. - """ - delta_coil::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0) - - """ - Raw `2msing × 2msing` outer-region matching matrix D' (Pletzer–Dewar 1991, before parity - rotation), side-major ordering `[L_s1, R_s1, …, L_sm, R_sm]`. `pest3_decompose` recovers the - A', B', Γ', Δ' blocks. No ½ prefactor (matches Fortran rdcon; the paper applies ½). - Empty unless the Riccati integrator was used. - """ - delta_prime_raw::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0) + delta_prime_raw_freeboundary::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0) + delta_prime_raw_fixedboundary::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0) + delta_ef::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0) + delta_prime_edge_response::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0) end """ @@ -124,7 +113,7 @@ Physical and workflow controls for a force-balance solve, together with its sele """ @kwdef struct ForceFreeStatesControl{A<:AbstractIntegrator} algorithm::A = Riccati() - vac_flag::Bool = false + vac_flag::Bool = true mtheta_vacuum::Int = 480 nzeta_vacuum::Int = 1 nn_low::Int = 0 diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index d40700059..e4fb18083 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -135,7 +135,7 @@ OdeState(numpert_total::Int, numsteps_init::Int, numunorms_init::Int, msing::Int **Stub, not physically valid.** Per-surface Δ' estimate `(ca_r[i,i,2,s] - ca_l[i,i,2,s]) / (4π²·psio)` from the asymptotic coefficients. It ignores inter-surface coupling and the vacuum, so it does not -agree with the physical `ForceFreeStatesInternal.delta_prime_matrix`. Used only by tests. +agree with the physical the BVP Δ' (`DeltaPrimeData`). Used only by tests. """ function compute_delta_prime_from_ca!(odet::OdeState, intr::ForceFreeStatesInternal, equil::Equilibrium.PlasmaEquilibrium) denom = (2π)^2 * equil.psio # = twopi * chi1 in SingularCoupling.jl diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index 0d971fa07..762538971 100644 --- a/src/ForceFreeStates/Free.jl +++ b/src/ForceFreeStates/Free.jl @@ -87,7 +87,7 @@ downstream code consumes the stored ξ profiles. end """ - compute_scaled_wv(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal) -> (wv, vac) + compute_scaled_wv(ctrl, equil, modes) -> (wv, vac) Vacuum response matrix at the control surface `intr.psilim`, scaled by the singular factors `(m - n·q)(m' - n'·q)` [Chance Phys. Plasmas 1997 2161 eq. 126]. Handles 2D single-n, 2D @@ -98,7 +98,7 @@ share it. Returns the scaled `wv` alongside the full vacuum response `vac`, whos point clouds the free-boundary result carries to HDF5. `wv` aliases `vac.wv`, which is scaled in place. """ -function compute_scaled_wv(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal) +function compute_scaled_wv(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ModeSpace) (; mlow, mhigh, nlow, nhigh, psilim, qlim, wall_settings) = intr vac_inputs = Vacuum.VacuumInput(equil, psilim, ctrl.mtheta_vacuum, ctrl.nzeta_vacuum, mlow:mhigh, nlow:nhigh) vac = Vacuum.compute_vacuum_response(vac_inputs, wall_settings) @@ -109,16 +109,18 @@ function compute_scaled_wv(ctrl::ForceFreeStatesControl, equil::Equilibrium.Plas end """ - free_run(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal) -> FreeBoundaryResult + free_run(wp, ctrl, equil, mats, modes) -> FreeBoundaryResult -Free-boundary energies at `intr.psilim`: the plasma response from the ODE edge state plus the -singfac-scaled vacuum matrix, diagonalized as the pencil (W, N). Returns the eigenmodes and their -plasma/vacuum/total energies as a [`FreeBoundaryResult`](@ref). +Free-boundary energies and eigenmodes of `W = wp + wv` at `modes.psilim`, for a plasma energy matrix `wp` +(Hermitian or not) and the singfac-scaled vacuum matrix `wv` of `ctrl`'s wall. """ -@with_pool pool function free_run(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal) +@with_pool pool function free_run(wp::Matrix{ComplexF64}, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, + mats::MatrixSplines, modes::ModeSpace; verbose::Bool=false) + + wv, vac = compute_scaled_wv(ctrl, equil, modes) # Initializations and allocations - (; mpert, numpert_total, psilim, npert) = intr + (; mpert, numpert_total, psilim, npert) = modes wpt = zeros!(pool, ComplexF64, numpert_total, numpert_total) wvt = zeros!(pool, ComplexF64, numpert_total, numpert_total) tmp_mat = zeros!(pool, ComplexF64, numpert_total, numpert_total) @@ -126,12 +128,6 @@ plasma/vacuum/total energies as a [`FreeBoundaryResult`](@ref). # Evaluate dV/dpsi at the plasma edge dV_dpsi = equil.profiles.dVdpsi_spline(psilim) - # Compute plasma response matrix W = U₂ * U₁⁻¹ - wp = zeros(ComplexF64, numpert_total, numpert_total) - @views wp .= (odet.u[:, :, 2] / odet.u[:, :, 1]) ./ equil.psio^2 - - wv, vac = compute_scaled_wv(ctrl, equil, intr) - # Power-normalization matrix N at the plasma edge: ξ†·N·ξ = ⟨|ξ|²⟩ (see power_norm_matrix!). # The Jacobian band is evaluated at psilim (same surface as W), not at the last grid surface. Nmat = zeros!(pool, ComplexF64, numpert_total, numpert_total) @@ -187,7 +183,7 @@ plasma/vacuum/total energies as a [`FreeBoundaryResult`](@ref). ev = diag(wvt) # Write energies to screen - if intr.verbose + if verbose @info "Least Stable Eigenmode Energies:\n" * " Plasma = $((@sprintf "%+.3e %+.3ei" real(ep[1]) imag(ep[1])))\n" * " Vacuum = $((@sprintf "%+.3e %+.3ei" real(ev[1]) imag(ev[1])))\n" * diff --git a/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl b/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl index 55d093194..5c453a1c8 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl @@ -409,21 +409,20 @@ function gal_assemble_rhs!(ws::GalWorkspace, mpert::Int) end """ - gal_set_boundary!(ws, mpert, wv_edge) + gal_set_boundary!(ws, mpert, wv_edge; nef=ws.nsol - 2msing) Apply the axis (`u(0)=0` identity) and edge boundary conditions to the boundary cells' `mat`, plus the -rpec coil-source RHS injection, before global assembly. Port of `gal_set_boundary` (gal.f). The edge +error-field (EF) source RHS injection, before global assembly. Port of `gal_set_boundary` (gal.f). The edge branch is one of: - - rpec (`ncoil = ws.nsol - 2*msing > 0`): identity edge and a unit source at the edge value DOF of coil + - EF columns (`nef > 0`, `error_field_response`): identity edge and a unit source at the edge value DOF of EF column `2*msing+ipert` for each poloidal mode `ipert`. `gal_assemble_rhs!` only fills columns `1:2*msing`, so the injection is not clobbered. - free boundary (`wv_edge !== nothing`): add the vacuum block `wvac·psio²` at the edge value DOF. - fixed boundary (else): identity edge. """ -function gal_set_boundary!(ws::GalWorkspace, mpert::Int, wv_edge::Union{Nothing,Matrix{ComplexF64}}) +function gal_set_boundary!(ws::GalWorkspace, mpert::Int, wv_edge::Union{Nothing,Matrix{ComplexF64}}; nef::Int=ws.nsol - 2 * (length(ws.intvl) - 1)) msing = length(ws.intvl) - 1 - ncoil = ws.nsol - 2 * msing # axis u(0): first cell of interval 0, Hermite value DOF ip=0 (Julia index 1) cell = ws.intvl[1].cells[1] @@ -434,8 +433,8 @@ function gal_set_boundary!(ws::GalWorkspace, mpert::Int, wv_edge::Union{Nothing, # edge u(1): last cell of last interval, edge-value DOF ip=3 (Julia index 4, see gauss_quad swap) cell = ws.intvl[msing+1].cells[ws.nx] - if ncoil > 0 - # rpec: identity edge + unit coil sources (gal.f) + if nef > 0 + # EF columns: identity edge + unit EF sources (gal.f) cell.mat[:, :, 4, :] .= 0 for idx in 1:mpert cell.mat[idx, idx, 4, 4] = 1 diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl index 5d4f6d2df..6b7eea41d 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl @@ -1,20 +1,21 @@ # GalerkinSolve.jl # # Top-level driver for the RDCON outer-region singular Galerkin Δ′ solve, plus the per-cell assembly -# orchestration, the banded solve, Δ′ extraction, PEST-3 blocks, and HDF5 output. -# Ports gal_make_arrays (gal.f), gal_solve (gal.f), and gal_write_pest3_data +# orchestration, the banded solve, Δ′ extraction, and HDF5 output. Ports gal_make_arrays and gal_solve # (gal.f). Inner-layer matching is applied afterwards by MatchProblem. """ - gal_make_arrays!(ws, ctrl, equil, mats, intr, asymps, sings, nn, wv_edge) + gal_make_arrays!(ws, ctrl, equil, mats, intr, asymps, sings, nn) -Assemble the global banded matrix and RHS. For each cell: Gauss-Lobatto Hermite stiffness -(`gal_gauss_quad!`), then the resonant (`gal_resonant!`) or extension (`gal_extension!`) contributions; -then the boundary conditions and the scatter into `ws.mat`/`ws.rhs`. Port of `gal_make_arrays`. +Assemble the global banded matrix and RHS with the fixed edge (plus the EF sources): per-cell stiffness and +resonant/extension contributions, the boundary conditions, and the scatter into `ws.mat`/`ws.rhs`. Port of +`gal_make_arrays`. + +Returns the edge cell's stiffness before the boundary conditions, `(mpert, mpert, np+1, np+1)`. """ function gal_make_arrays!(ws::GalWorkspace, ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines, intr::ForceFreeStatesInternal, asymps::Vector{GalSingAsymp}, sings::Vector{SingType}, - nn::Int, wv_edge::Union{Nothing,Matrix{ComplexF64}}) + nn::Int) mpert = intr.numpert_total msing = length(ws.intvl) - 1 @@ -33,17 +34,61 @@ function gal_make_arrays!(ws::GalWorkspace, ctrl::ForceFreeStatesControl, equil, end end - gal_set_boundary!(ws, mpert, wv_edge) + edge_mat = copy(ws.intvl[msing+1].cells[ws.nx].mat) + gal_set_boundary!(ws, mpert, nothing) gal_assemble_mat!(ws, mpert) gal_assemble_rhs!(ws, mpert) - return ws + return edge_mat +end + +""" + gal_vacuum_edge_solve(ws, edge_mat, wv_edge, mpert) -> Matrix{ComplexF64} + +Solve the 2·msing big-solution columns of an assembled workspace again with the vacuum edge `wv_edge` +instead of the fixed edge, restoring the edge cell from `edge_mat`. Returns the `(ndim, 2msing)` +solution, which gives the free-boundary Δ′. +""" +function gal_vacuum_edge_solve(ws::GalWorkspace, edge_mat::Array{ComplexF64,4}, wv_edge::Matrix{ComplexF64}, mpert::Int) + msing = length(ws.intvl) - 1 + ws.intvl[msing+1].cells[ws.nx].mat .= edge_mat + gal_set_boundary!(ws, mpert, wv_edge; nef=0) + fill!(ws.mat, 0) + gal_assemble_mat!(ws, mpert) + sol = ws.rhs[:, 1:(2msing)] + gal_banded_solve!(ws.mat, sol, ws.kl, ws.ku, "auto") + return sol +end + +""" + gal_small_coefficients(ws, sol) -> Matrix{ComplexF64} + +Small-solution coefficients of the columns of `sol`, `(ncol, 2msing)`, surface sides in side-major order. +""" +function gal_small_coefficients(ws::GalWorkspace, sol::AbstractMatrix{ComplexF64}) + rows = [cell.emap for intvl in ws.intvl for cell in intvl.cells if cell.etype == GCT_RES] + return permutedims(sol[rows, :]) +end + +""" + gal_edge_response(edge_rows, edge_map, sol) -> Matrix{ComplexF64} + +Edge momentum FΞ′ + KΞ, `(mpert, nsol)`, of each column of `sol`, from the unconstrained edge rows +`edge_rows` of the edge cell mapped by `edge_map`. On the columns with `ξ(edge) = eⱼ` it is `W_P·ψ₀²` +[Glasser 2016 PoP 23 072505, Eq. 66]. +""" +function gal_edge_response(edge_rows::Array{ComplexF64,3}, edge_map::Matrix{Int}, sol::Matrix{ComplexF64}) + force = zeros(ComplexF64, size(edge_rows, 1), size(sol, 2)) + for jp in axes(edge_rows, 3) + mul!(force, @view(edge_rows[:, :, jp]), @view(sol[edge_map[:, jp], :]), true, true) + end + return force end """ Empty `GalerkinResult` for a domain with no resonant surfaces. """ function empty_galerkin_result() - return GalerkinResult(0, Float64[], Float64[], Int[], Int[], Float64[], ComplexF64[], nothing, "none") + return GalerkinResult(0, Float64[], Float64[], Int[], Int[], Float64[], ComplexF64[], nothing, "none", Matrix{ComplexF64}(undef, 0, 0)) end """ @@ -55,9 +100,9 @@ Compute the outer-region Δ′ matching matrix by the singular Galerkin method. `GalerkinResult` alongside the Δ′ payload as the shared [`DeltaPrimeData`](@ref); if there are no resonant surfaces in the domain it returns an empty result and `nothing`. -`wv` is the vacuum energy matrix from `free_run`, which supplies the free-boundary edge term -`wv_edge = wv · psio²`; pass `nothing` for a fixed-boundary edge. `cut_solution` also keeps the -cut solution, for a later `MatchProblem`'s composite inner profiles. +The fixed-boundary Δ′ is always computed; `wv`, the singfac-scaled vacuum matrix, adds the free-boundary +Δ′ unless `free_boundary_delta_prime` is off. `cut_solution` also keeps the cut solution, for a later +`MatchProblem`'s composite inner profiles. """ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines, intr::ForceFreeStatesInternal; wv=nothing, cut_solution::Bool=false) @@ -99,10 +144,10 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines kl = mpert * (np + 1) ku = kl ldab = 2kl + ku + 1 - # rpec_flag (RDCON, gal.f): append mpert coil-response columns. Each is a unit source at - # the edge value DOF in one poloidal mode; the recorded plasma response is the coil block of Δ_gw. - ncoil = ctrl.algorithm.rpec_flag ? mpert : 0 - nsol = 2 * msing + ncoil + # error_field_response (RDCON rpec_flag, gal.f): append mpert error-field (EF) columns, each a unit edge displacement in one + # poloidal mode with no big-solution drive; their small coefficients are the Δ_EF block. + nef = ctrl.algorithm.error_field_response ? mpert : 0 + nsol = 2 * msing + nef intvl = [GalInterval(zeros(Float64, nx + 1), zeros(Float64, nx + 1), [GalCell(mpert) for _ in 1:nx]) for _ in 0:msing] ws = GalWorkspace("auto", nx, ctrl.algorithm.nq, np, 0, kl, ku, ldab, nsol, @@ -117,15 +162,7 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines ws.rhs = zeros(ComplexF64, ws.ndim, nsol) ws.sol = zeros(ComplexF64, ws.ndim, nsol) - # Free-boundary edge term wvac·psio² (wv is already singfac-scaled at qlim in free_run). - # rpec passes wv_edge=nothing; gal_set_boundary! then applies the identity edge AND injects the coil - # unit sources (its three-way branch). The vacuum block is only built for the non-rpec free case. - wv_edge = nothing - if ctrl.vac_flag && wv !== nothing && ncoil == 0 - wv_edge = Matrix{ComplexF64}(wv .* equil.psio^2) - end - - gal_make_arrays!(ws, ctrl, equil, mats, intr, asymps, sings, nn, wv_edge) + edge_mat = gal_make_arrays!(ws, ctrl, equil, mats, intr, asymps, sings, nn) if intr.verbose offdbg = ws.kl + ws.ku + 1 @@ -143,43 +180,37 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines ws.sol .= ws.rhs ws.solver = gal_banded_solve!(ws.mat, ws.sol, ws.kl, ws.ku, "auto") intr.verbose && @info "Galerkin banded solve: $(ws.solver)" + edge_response = gal_edge_response(edge_mat[:, :, 4, :], ws.intvl[msing+1].cells[ws.nx].map, ws.sol) # --- extract Δ′ from the small resonant coefficients (gal.f) --- - delta = zeros(ComplexF64, nsol, 2 * msing) - jsol = 0 - for ising in 0:msing - for cell in ws.intvl[ising+1].cells - if cell.etype == GCT_RES - jsol += 1 - @views delta[:, jsol] .= ws.sol[cell.emap, :] - end - end - end - - Ap, Bp, Gammap, Deltap = gal_pest3_blocks(delta, msing) + delta = gal_small_coefficients(ws, ws.sol) # Mercier index D_I = -Re(α²). α is the half-separation of the resonant exponents. di = [real(-asymps[i].right.alpha[1]^2) for i in 1:msing] alpha = [asymps[i].right.alpha[1] for i in 1:msing] - # Pack the Δ′ payload in the shared layout: the raw D′ is the leading 2msing×2msing side-major - # block, and the coil-response rows 2*msing+1 : 2*msing+mpert (rpec_flag) transpose into the - # (surface-side × edge mode) orientation the Riccati BVP also produces. - dp_raw = delta[1:(2*msing), 1:(2*msing)] - dp_coil = ncoil > 0 ? permutedims(delta[(2*msing+1):(2*msing+ncoil), :]) : Matrix{ComplexF64}(undef, 0, 0) - dp = DeltaPrimeData(Deltap, dp_raw, dp_coil, Ap, Bp, Gammap) - # Reconstruct ξ(ψ) AND analytic ξ′(ψ) on the gal-native grid (gal_output_solution); the cut # solution only when a later match needs the composite inner profiles. intr.verbose && @info "Reconstructing outer-region ξ and analytic ξ′ on the gal grid" solution = gal_output_solution(ws, asymps, sings, intr, equil.profiles, psihigh; delta=(cut_solution ? delta : nothing)) + # The fixed-edge solve gives the fixed-boundary D′ and, in the (surface-side × edge mode) orientation, Δ_EF. + # The free-boundary D′ comes from re-solving the big-solution columns with the vacuum edge. + raw_fixed = delta[1:(2msing), :] + ef = nef > 0 ? permutedims(delta[(2msing+1):end, :]) : Matrix{ComplexF64}(undef, 0, 0) + raw_free = Matrix{ComplexF64}(undef, 0, 0) + if ctrl.vac_flag && wv !== nothing && ctrl.algorithm.free_boundary_delta_prime + # wv is singfac-scaled at qlim; the Galerkin edge term is wv·ψ₀². + raw_free = gal_small_coefficients(ws, gal_vacuum_edge_solve(ws, edge_mat, Matrix{ComplexF64}(wv .* equil.psio^2), mpert)) + end + dp = DeltaPrimeData(raw_free, raw_fixed, ef, edge_response[:, 1:(2msing)]) + sing_psi = [s.psifac for s in sings] sing_q = [s.q for s in sings] sing_m = [s.m[1] for s in sings] sing_n = [s.n[1] for s in sings] - return GalerkinResult(msing, sing_psi, sing_q, sing_m, sing_n, di, alpha, solution, ws.solver), dp + return GalerkinResult(msing, sing_psi, sing_q, sing_m, sing_n, di, alpha, solution, ws.solver, edge_response), dp end """ @@ -204,30 +235,6 @@ function gal_scaled_lu_solve!(mat::Matrix{ComplexF64}, sol::AbstractVecOrMat{Com return sol end -""" - gal_pest3_blocks(delta, msing) -> (Ap, Bp, Gammap, Deltap) - -PEST-3 matching blocks (each `msing×msing`) as ± combinations of the left/right Δ′ entries. Port of -`gal_write_pest3_data`. Row index `2i-1`/`2i` = left/right of surface `i`. -""" -function gal_pest3_blocks(delta::Matrix{ComplexF64}, msing::Int) - Ap = zeros(ComplexF64, msing, msing) - Bp = zeros(ComplexF64, msing, msing) - Gammap = zeros(ComplexF64, msing, msing) - Deltap = zeros(ComplexF64, msing, msing) - for i in 1:msing, j in 1:msing - d11 = delta[2i-1, 2j-1] - d12 = delta[2i-1, 2j] - d21 = delta[2i, 2j-1] - d22 = delta[2i, 2j] - Ap[i, j] = d22 + d21 + d12 + d11 - Bp[i, j] = d22 - d21 + d12 - d11 - Gammap[i, j] = d22 + d21 - d12 - d11 - Deltap[i, j] = d22 - d21 - d12 + d11 - end - return Ap, Bp, Gammap, Deltap -end - """ write_galerkin!(out_h5, result::GalerkinResult; basis_output=false) @@ -350,7 +357,7 @@ end gal_cholesky_solve!(mat, sol, kl, ku) -> Bool Banded Cholesky solve of the full-band system `mat`, leaving `mat` untouched. Rows that are -identity rows (Dirichlet DOFs: axis, fixed or coil edge) have their columns moved to the +identity rows (Dirichlet DOFs: axis, fixed or EF edge) have their columns moved to the right-hand side. The free block's off-diagonal entries must agree with their conjugate partners within `1e-10 * sqrt(abs(A_ii)) * sqrt(abs(A_jj))`. Cholesky uses the lower band and deliberately projects the diagonal onto its real part. Returns `false`, with `sol` diff --git a/src/ForceFreeStates/Galerkin/GalerkinStructs.jl b/src/ForceFreeStates/Galerkin/GalerkinStructs.jl index e0cb88afc..db0c31fc8 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinStructs.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinStructs.jl @@ -150,6 +150,8 @@ Solver internals and FEM diagnostics of the outer-region Galerkin solve. The Δ - `solution::Union{Nothing,GalerkinSolution}` — reconstructed radial ξ(ψ) and analytic ξ′(ψ) on the gal-native grid; `nothing` if no resonant surfaces. - `solver::String` — factorization actually used: `"cholesky"`, `"LU"`, or `"none"` when no resonant surfaces were present. + - `edge_response::Matrix{ComplexF64}` — `(mpert, nsol)` edge momentum FΞ′ + KΞ of every solution column + (`gal_edge_response`); the EF columns give `W_P·ψ₀²`. Empty if no resonant surfaces. """ struct GalerkinResult msing::Int @@ -161,4 +163,5 @@ struct GalerkinResult alpha::Vector{ComplexF64} solution::Union{Nothing,GalerkinSolution} solver::String + edge_response::Matrix{ComplexF64} end diff --git a/src/ForceFreeStates/Integrators.jl b/src/ForceFreeStates/Integrators.jl index 8210f4f8c..e6cd24f3b 100644 --- a/src/ForceFreeStates/Integrators.jl +++ b/src/ForceFreeStates/Integrators.jl @@ -41,6 +41,8 @@ Chunked Riccati integration and the singular-surface boundary-value solve. - `save_interval::Int` - Save every Nth ODE step; always retain steps near rational surfaces. - `nchunks::Int` - Target chunk count; zero derives the count from the singular surfaces. - `extended_precision_bvp::Bool` - Use Complex{Double64} for the Delta-prime BVP and PEST3 combination to reduce cancellation error. + - `error_field_response::Bool` - Compute `Delta_EF`, the response to a unit edge displacement per mode (needs `vac_flag`). + Required by `MatchProblem`. """ @kwdef struct Riccati <: AbstractIntegrator eulerlagrange_tolerance::Float64 = 1e-8 @@ -50,6 +52,7 @@ Chunked Riccati integration and the singular-surface boundary-value solve. save_interval::Int = 1 nchunks::Int = 0 extended_precision_bvp::Bool = true + error_field_response::Bool = true end """ @@ -73,7 +76,10 @@ corresponding discretization setting. Factorization is selected automatically. - `gnstep::Int` - Maximum resonant-quadrature evaluations. - `dx1dx2_flag::Bool` - Enable the special dx1/dx2 treatment of resonant and extension elements. - `frobenius_order_ceiling::Bool` - Auto-raise the order per surface for a high Mercier index. - - `rpec_flag::Bool` - Append the mpert coil-response columns to the Δ′ solve. Required for a later `MatchProblem` solve on the result. + - `error_field_response::Bool` - Add the error-field (EF) columns, a unit edge displacement per mode: gives `Delta_EF`, the ideal solution + and `wp`. Required by `MatchProblem`. + - `free_boundary_delta_prime::Bool` - With `error_field_response` and `vac_flag`, compute the free-boundary (tearing) Δ′ by a second solve + (about a quarter of the solve time). Off: no tearing Δ′. - `edge_onesided::Bool` - Pack the two end intervals one-sided toward their single rational end instead of the Fortran symmetric pack. - `basis_output::Bool` - Save the raw Galerkin basis for development diagnostics. - `cut_solution::Bool` - Also keep the cut solution, for a later `MatchProblem`'s composite inner profiles. @@ -90,7 +96,8 @@ corresponding discretization setting. Factorization is selected automatically. gnstep::Int = 20000 dx1dx2_flag::Bool = true frobenius_order_ceiling::Bool = true - rpec_flag::Bool = false + error_field_response::Bool = true + free_boundary_delta_prime::Bool = true edge_onesided::Bool = false cut_solution::Bool = false basis_output::Bool = false diff --git a/src/ForceFreeStates/Matching/DeltaPrime.jl b/src/ForceFreeStates/Matching/DeltaPrime.jl index 7358dceee..69628a772 100644 --- a/src/ForceFreeStates/Matching/DeltaPrime.jl +++ b/src/ForceFreeStates/Matching/DeltaPrime.jl @@ -1,35 +1,67 @@ """ - DeltaPrimeData + PEST3Blocks -The solve's Δ′/outer-region matching payload, in the PEST-3 convention shared by the Riccati -boundary-value problem and the Galerkin solve, so consumers never branch on which integrator ran. +PEST-3 parity blocks of a raw D′ [Glasser-Kolemen 2018 PoP 25, 032501, Eq. 31], each (msing × msing). -Matrices are indexed by the surfaces the formalism solved across, core→edge: `result.surfaces` for -Riccati, its in-domain in-band subset for Galerkin (`gal_resonant_surfaces`). Side-major orderings -run `[L_s1, R_s1, L_s2, R_s2, …]`. +## Fields + + - `A` - Interchange↔interchange block. + - `B` - Interchange↔tearing block. + - `Gamma` - Tearing↔interchange block. + - `Delta` - Tearing↔tearing block, the tearing Δ′. +""" +struct PEST3Blocks + A::Matrix{ComplexF64} + B::Matrix{ComplexF64} + Gamma::Matrix{ComplexF64} + Delta::Matrix{ComplexF64} +end + +""" + pest3_decompose(raw) -> PEST3Blocks + +PEST-3 parity blocks of the side-major raw D′ `raw` (2msing × 2msing). No Pletzer–Dewar ½ prefactor, as in +Fortran RDCON. +""" +function pest3_decompose(raw::AbstractMatrix{ComplexF64}) + m = size(raw, 1) ÷ 2 + size(raw) == (2m, 2m) || throw(ArgumentError("pest3_decompose: raw must be 2m × 2m, got $(size(raw))")) + LL = raw[1:2:end, 1:2:end] + LR = raw[1:2:end, 2:2:end] + RL = raw[2:2:end, 1:2:end] + RR = raw[2:2:end, 2:2:end] + return PEST3Blocks(RR + RL + LR + LL, RR - RL + LR - LL, RR + RL - LR - LL, RR - RL - LR + LL) +end + +""" + DeltaPrimeData(raw_freeboundary, raw_fixedboundary, ef, edge_response) + +The solve's Δ′/outer-region matching payload. Matrices are indexed by the surfaces the formalism solved +across, core→edge; side-major orderings run `[L_s1, R_s1, L_s2, R_s2, …]`. An empty matrix (or `nothing`) +means the solve did not compute it. The PEST-3 blocks are built from the raw matrices. ## Fields - - `matrix::Matrix{ComplexF64}` - Inter-surface Δ′ of shape (msing × msing) in PEST3 - convention, the tearing↔tearing parity projection of `raw`. Same as `Delta` of the - PEST-3 block set. Both formalisms. - - `raw::Matrix{ComplexF64}` - Raw outer-region matching matrix D′ of shape - (2msing × 2msing), side-major on both axes. Both formalisms. - - `coil::Matrix{ComplexF64}` - Edge coil-response matrix of shape - (2msing × numpert_total); column k is the resonant small-solution response at each - surface side to a unit source on edge poloidal mode k. Riccati fills it from the - vacuum-edge BVP, Galerkin from its `rpec_flag` columns. Empty when neither ran. - - `A::Union{Nothing,Matrix{ComplexF64}}` - PEST-3 interchange↔interchange block - (msing × msing). Galerkin only; `nothing` for Riccati, which persists only `raw` and - recovers the blocks on demand via `pest3_decompose`. - - `B::Union{Nothing,Matrix{ComplexF64}}` - PEST-3 interchange↔tearing block; see `A`. - - `Gamma::Union{Nothing,Matrix{ComplexF64}}` - PEST-3 tearing↔interchange block; see `A`. + - `raw_freeboundary::Matrix{ComplexF64}` - Raw D′ (2msing × 2msing) with the vacuum (or wall) at the edge. + - `raw_fixedboundary::Matrix{ComplexF64}` - Raw D′ (2msing × 2msing) with ξ(edge) = 0, a conformal wall on the plasma. + - `ef::Matrix{ComplexF64}` - Error-field response Δ_EF (2msing × numpert_total): column k is the small-solution + coefficient at each surface side for a unit edge displacement in mode k with the edge otherwise fixed. It pairs + with `raw_fixedboundary` in matching. + - `edge_response::Matrix{ComplexF64}` - Edge momentum Fξ′ + Kξ (numpert_total × 2msing) of the `raw_fixedboundary` + solutions: column k has unit big-solution coefficient on surface side k and ξ(edge) = 0. + - `freeboundary::Union{Nothing,PEST3Blocks}` - PEST-3 blocks of `raw_freeboundary`. + - `fixedboundary::Union{Nothing,PEST3Blocks}` - PEST-3 blocks of `raw_fixedboundary`. """ struct DeltaPrimeData - matrix::Matrix{ComplexF64} - raw::Matrix{ComplexF64} - coil::Matrix{ComplexF64} - A::Union{Nothing,Matrix{ComplexF64}} - B::Union{Nothing,Matrix{ComplexF64}} - Gamma::Union{Nothing,Matrix{ComplexF64}} + raw_freeboundary::Matrix{ComplexF64} + raw_fixedboundary::Matrix{ComplexF64} + ef::Matrix{ComplexF64} + edge_response::Matrix{ComplexF64} + freeboundary::Union{Nothing,PEST3Blocks} + fixedboundary::Union{Nothing,PEST3Blocks} end + +DeltaPrimeData(raw_freeboundary::AbstractMatrix, raw_fixedboundary::AbstractMatrix, ef::AbstractMatrix, edge_response::AbstractMatrix) = + DeltaPrimeData(raw_freeboundary, raw_fixedboundary, ef, edge_response, + isempty(raw_freeboundary) ? nothing : pest3_decompose(raw_freeboundary), + isempty(raw_fixedboundary) ? nothing : pest3_decompose(raw_fixedboundary)) diff --git a/src/ForceFreeStates/Matching/MatchProblem.jl b/src/ForceFreeStates/Matching/MatchProblem.jl index f80f52eb4..eb783aacf 100644 --- a/src/ForceFreeStates/Matching/MatchProblem.jl +++ b/src/ForceFreeStates/Matching/MatchProblem.jl @@ -3,13 +3,13 @@ # Driven inner-layer matching on a finished solve (rmatch match_rpec + match_output_solution, match.f). """ - MatchProblem(ffs; eta=nothing, rho=nothing, rotation=nothing, gamma=5/3, ideal=false) + MatchProblem(ffs; eta=nothing, rho=nothing, rotation=nothing, gamma=5/3) Driven inner-layer matching on a finished force-free-states solve `ffs`: each rational surface is matched at the forced eigenvalue γ = 2πi·n·f, with f the plasma rotation frequency there. `solve(prob, GGJModel())` returns a copy of `ffs` with the matched solution and the penetrated -field `bpen`. `ffs` needs the Δ′ coil block, from `Galerkin(; rpec_flag=true)` or `Riccati()` -with `vac_flag=true`. +field `bpen`. `ffs` needs `Delta_EF`, from `Galerkin()` or `Riccati()` with `error_field_response` (Riccati +also needs `vac_flag`). # Keywords @@ -17,17 +17,16 @@ with `vac_flag=true`. in Hz, a scalar for every surface or one value per `ffs.surfaces` entry (core to edge). Each one not given comes from the equilibrium's kinetic data, as the tearing problem's do: η from [`InnerLayer.layer_eta`](@ref), ρ from [`InnerLayer.layer_rho`](@ref), and the - rotation f = Ω_E/2π from the E×B rotation of [`InnerLayer.layer_omega_E`](@ref). Unused - with `ideal=true`. + rotation f = Ω_E/2π from the E×B rotation of [`InnerLayer.layer_omega_E`](@ref). - `gamma`: ratio of specific heats. - - `ideal`: skip the inner layer and return the perfectly shielding solution. + +The perfectly shielding (ideal) solution needs no match: it is the Galerkin solve's own result. ## Fields - `ffs` - The outer solve being matched; held by reference, not copied. - `eta`, `rho`, `rotation` - Per-surface η in Ω·m, ρ in kg/m³ and rotation frequency in Hz, following `ffs.surfaces`. - `gamma` - Ratio of specific heats. - - `ideal` - Whether the inner layer is skipped. Ref: Fortran rmatch `match_rpec` (match.f). """ @@ -37,7 +36,6 @@ struct MatchProblem{R<:ForceFreeStatesResult} rho::Vector{Float64} rotation::Vector{Float64} gamma::Float64 - ideal::Bool end function MatchProblem( @@ -45,16 +43,15 @@ function MatchProblem( eta::Union{Nothing,Real,AbstractVector{<:Real}}=nothing, rho::Union{Nothing,Real,AbstractVector{<:Real}}=nothing, rotation::Union{Nothing,Real,AbstractVector{<:Real}}=nothing, - gamma::Real=5 / 3, - ideal::Bool=false + gamma::Real=5 / 3 ) ffs.npert == 1 || error("MatchProblem: matching is single-n; the result spans n = $(ffs.nlow):$(ffs.nhigh)") ffs.closure === :matched && error("MatchProblem: the result is already matched; pose the match on the outer solve") dp = ffs.delta_prime dp === nothing && - error("MatchProblem: the $(ffs.integrator) result has no Δ′; solve with Galerkin(; rpec_flag=true), or Riccati() with vac_flag=true") - isempty(dp.coil) && - error("MatchProblem: the Δ′ coil block is empty; solve with Galerkin(; rpec_flag=true), or Riccati() with vac_flag=true") + error("MatchProblem: the $(ffs.integrator) result has no Δ′; solve with Galerkin() or Riccati() (Riccati needs vac_flag=true)") + (isempty(dp.ef) || isempty(dp.raw_fixedboundary)) && + error("MatchProblem: needs Δ_EF and the fixed-boundary Δ′; solve with error_field_response=true (Riccati also needs vac_flag=true)") msing = length(ffs.surfaces) if ffs.galerkin !== nothing @@ -62,10 +59,9 @@ function MatchProblem( isempty(dropped) || error("MatchProblem: the Galerkin solve dropped the surfaces at ψ = $dropped (on the domain edge or " * "outside mlow:mhigh); matching needs every surface") end - size(dp.raw, 1) == 2msing || - error("MatchProblem: Δ′ raw block is $(size(dp.raw, 1))×$(size(dp.raw, 2)) but the result carries $msing surfaces") + size(dp.raw_fixedboundary, 1) == 2msing || + error("MatchProblem: Δ′ raw block is $(size(dp.raw_fixedboundary, 1))×$(size(dp.raw_fixedboundary, 2)) but the result carries $msing surfaces") - ideal && return MatchProblem(ffs, zeros(msing), zeros(msing), zeros(msing), Float64(gamma), true) InnerLayer._require_layer_inputs("MatchProblem", ffs.equil; eta, rho, rotation) eta === nothing && (eta = InnerLayer.layer_eta(ffs.equil, ffs.surfaces)) rho === nothing && (rho = InnerLayer.layer_rho(ffs.equil, ffs.surfaces)) @@ -78,7 +74,7 @@ function MatchProblem( v isa Real ? fill(Float64(v), msing) : length(v) == msing ? collect(Float64, v) : error("MatchProblem: $name has length $(length(v)), expected $msing (one per surface, core to edge)") - return MatchProblem(ffs, per_surface("eta", eta), per_surface("rho", rho), per_surface("rotation", rotation), Float64(gamma), false) + return MatchProblem(ffs, per_surface("eta", eta), per_surface("rho", rho), per_surface("rotation", rotation), Float64(gamma)) end """ @@ -94,16 +90,16 @@ closure_capable(::InnerLayer.SLAYERModel) = false """ solve(prob::MatchProblem, model) -> ForceFreeStatesResult -Match `prob.ffs` to the inner layer of `model`. The result has `closure = :matched` (`:ideal` -with `ideal=true`), `bpen`, the [`MatchResult`](@ref) in `result.match` and, for a Galerkin -solve, the matched ξ. A `:matched` result drops the ideal δW (`wp`, `free_boundary`), so -PerturbedEquilibrium skips its response rather than mixing ideal and matched physics. +Match `prob.ffs` to the inner layer of `model`. The result has `closure = :matched`, `bpen`, the +[`MatchResult`](@ref) in `result.match`, the matched edge response `wp` (see [`ForceFreeStatesResult`](@ref)) +and, when the solve had `vac_flag`, the free-boundary energies of that response. A Galerkin solve also gives +the matched ξ. """ function CommonSolve.solve(prob::MatchProblem, model::InnerLayer.InnerLayerModel) closure_capable(model) || error("$(nameof(typeof(model))) cannot solve a MatchProblem (use GGJModel(); SLAYERModel is for TearingProblem)") match = _compute_match(prob, model) - return _matched_result(prob.ffs, match, prob.ideal) + return _matched_result(prob.ffs, match) end # Inner-layer Δ per surface, the matching system, then the matched outer and inner profiles. @@ -112,88 +108,74 @@ function _compute_match(prob::MatchProblem, model::InnerLayer.GGJModel) dp = ffs.delta_prime sings = ffs.surfaces msing = length(sings) - mcoil = size(dp.coil, 2) + nef = size(dp.ef, 2) nn = ffs.nlow equil = ffs.equil chi1 = 2π * equil.psio gal_sol = ffs.galerkin === nothing ? nothing : ffs.galerkin.solution - if prob.ideal - # Ideal limit (rmatch coil%ideal_flag): no inner layer, cout = 0. - cout = zeros(ComplexF64, 2msing, mcoil) - cin = zeros(ComplexF64, 2msing, mcoil) - deltar = zeros(ComplexF64, msing, 2) - bpen = zeros(ComplexF64, msing, mcoil) - inner_psi = Vector{Float64}[] - inner_xi = Matrix{ComplexF64}[] - inner_b = Matrix{ComplexF64}[] - inner_params = InnerLayer.GGJParameters[] - rpec_eig = zeros(ComplexF64, msing) - residual = 0.0 - else - # Inner-layer Δ and profiles per surface (match.f deltac_run, intotsol). - deltar = zeros(ComplexF64, msing, 2) - rpec_eig = zeros(ComplexF64, msing) - # Layer-center field weight per parity (match.f intotsol_b). - pen = zeros(ComplexF64, msing, 2) - inner_psi = Vector{Vector{Float64}}(undef, msing) # ψ_s ± X·dψdx, ascending - inner_odd = Vector{Vector{ComplexF64}}(undef, msing) # Ξ₋, odd about ψ_s - inner_even = Vector{Vector{ComplexF64}}(undef, msing) # Ξ₊, even - inner_bodd = Vector{Vector{ComplexF64}}(undef, msing) # b^ψ from Ψ₋, even - inner_beven = Vector{Vector{ComplexF64}}(undef, msing) # b^ψ from Ψ₊, odd - inner_params = Vector{InnerLayer.GGJParameters}(undef, msing) - profiles = true - for i in 1:msing - params = ggj_parameters(sings[i], equil; eta=prob.eta[i], rho=prob.rho[i], - gamma=prob.gamma, ising=i) - inner_params[i] = params - γ = 2π * im * nn * prob.rotation[i] # forced eigenvalue, f in Hz - rpec_eig[i] = γ - inner = InnerLayer.solve_inner_profile(model, params, γ) # parity columns (Δ₊, Δ₋) - deltar[i, 1] = inner.Δ[1] - deltar[i, 2] = inner.Δ[2] - profΨ, profΞ, xg = inner.Ψ, inner.Ξ, inner.x - resc = inner.rescale # inner big-branch amplitude → outer δψ normalization - # b_m = 2πi·χ₁·(m−nq)·ξ_m with m−nq = −n·q′·dψdx·X and Ψ = XΞ at leading order (GWP2016 Eqs. 16, 27; Ψ normalization Eq. A17) ⇒ b_m = scale·resc·Ψ. - scale = -2π * chi1 * im * nn * sings[i].q1 * inner.dψdx - pen[i, 1] = scale * profΨ[1, 1] * resc # X = 0, even parity (Ψ(0) = 0) - pen[i, 2] = scale * profΨ[1, 2] * resc # odd parity - # A backend returning only layer-center values (real-ψ profiles not requested) skips the inner profiles. - if length(xg) == 1 - profiles = false - continue - end - xvar = xg .* inner.dψdx # X → ψ distance - inner_psi[i] = vcat(reverse(sings[i].psifac .- xvar), sings[i].psifac .+ xvar) - inner_odd[i] = resc .* vcat(reverse(.-profΞ[:, 2]), profΞ[:, 2]) - inner_even[i] = resc .* vcat(reverse(profΞ[:, 1]), profΞ[:, 1]) - inner_bodd[i] = (scale * resc) .* vcat(reverse(profΨ[:, 2]), profΨ[:, 2]) - inner_beven[i] = (scale * resc) .* vcat(reverse(.-profΨ[:, 1]), profΨ[:, 1]) + # Inner-layer Δ and profiles per surface (match.f deltac_run, intotsol). + deltar = zeros(ComplexF64, msing, 2) + forced_eigenvalue = zeros(ComplexF64, msing) + # Layer-center field weight per parity (match.f intotsol_b). + pen = zeros(ComplexF64, msing, 2) + inner_psi = Vector{Vector{Float64}}(undef, msing) # ψ_s ± X·dψdx, ascending + inner_odd = Vector{Vector{ComplexF64}}(undef, msing) # Ξ₋, odd about ψ_s + inner_even = Vector{Vector{ComplexF64}}(undef, msing) # Ξ₊, even + inner_bodd = Vector{Vector{ComplexF64}}(undef, msing) # b^ψ from Ψ₋, even + inner_beven = Vector{Vector{ComplexF64}}(undef, msing) # b^ψ from Ψ₊, odd + inner_params = Vector{InnerLayer.GGJParameters}(undef, msing) + profiles = true + for i in 1:msing + params = ggj_parameters(sings[i], equil; eta=prob.eta[i], rho=prob.rho[i], + gamma=prob.gamma, ising=i) + inner_params[i] = params + γ = 2π * im * nn * prob.rotation[i] # forced eigenvalue, f in Hz + forced_eigenvalue[i] = γ + inner = InnerLayer.solve_inner_profile(model, params, γ) # parity columns (Δ₊, Δ₋) + deltar[i, 1] = inner.Δ[1] + deltar[i, 2] = inner.Δ[2] + profΨ, profΞ, xg = inner.Ψ, inner.Ξ, inner.x + resc = inner.rescale # inner big-branch amplitude → outer δψ normalization + # b_m = 2πi·χ₁·(m−nq)·ξ_m with m−nq = −n·q′·dψdx·X and Ψ = XΞ at leading order (GWP2016 Eqs. 16, 27; Ψ normalization Eq. A17) ⇒ b_m = scale·resc·Ψ. + scale = -2π * chi1 * im * nn * sings[i].q1 * inner.dψdx + pen[i, 1] = scale * profΨ[1, 1] * resc # X = 0, even parity (Ψ(0) = 0) + pen[i, 2] = scale * profΨ[1, 2] * resc # odd parity + # A backend returning only layer-center values (real-ψ profiles not requested) skips the inner profiles. + if length(xg) == 1 + profiles = false + continue end + xvar = xg .* inner.dψdx # X → ψ distance + inner_psi[i] = vcat(reverse(sings[i].psifac .- xvar), sings[i].psifac .+ xvar) + inner_odd[i] = resc .* vcat(reverse(.-profΞ[:, 2]), profΞ[:, 2]) + inner_even[i] = resc .* vcat(reverse(profΞ[:, 1]), profΞ[:, 1]) + inner_bodd[i] = (scale * resc) .* vcat(reverse(profΨ[:, 2]), profΨ[:, 2]) + inner_beven[i] = (scale * resc) .* vcat(reverse(.-profΨ[:, 1]), profΨ[:, 1]) + end - cout, cin, residual = _match_system(dp.raw, dp.coil, deltar) + cout, cin, residual = _match_system(dp.raw_fixedboundary, dp.ef, deltar) - # Penetrated field at each layer center (match.f intotsol_b); cin rows per surface are (d₊, d₋). - bpen = zeros(ComplexF64, msing, mcoil) - for i in 1:msing, j in 1:mcoil - bpen[i, j] = pen[i, 1] * cin[2i-1, j] + pen[i, 2] * cin[2i, j] - end + # Penetrated field at each layer center (match.f intotsol_b); cin rows per surface are (d₊, d₋). + bpen = zeros(ComplexF64, msing, nef) + for i in 1:msing, j in 1:nef + bpen[i, j] = pen[i, 1] * cin[2i-1, j] + pen[i, 2] * cin[2i, j] + end - # Matched inner ξ_ψ and b^ψ per coil: odd parity × cin[2i], even × cin[2i-1] (match.f intotsol). - if profiles - inner_xi = [inner_odd[i] * transpose(cin[2i, :]) .+ inner_even[i] * transpose(cin[2i-1, :]) for i in 1:msing] - inner_b = [inner_bodd[i] * transpose(cin[2i, :]) .+ inner_beven[i] * transpose(cin[2i-1, :]) for i in 1:msing] - else - inner_psi = Vector{Float64}[] - inner_xi = Matrix{ComplexF64}[] - inner_b = Matrix{ComplexF64}[] - end + # Matched inner ξ_ψ and b^ψ per EF column: odd parity × cin[2i], even × cin[2i-1] (match.f intotsol). + if profiles + inner_xi = [inner_odd[i] * transpose(cin[2i, :]) .+ inner_even[i] * transpose(cin[2i-1, :]) for i in 1:msing] + inner_b = [inner_bodd[i] * transpose(cin[2i, :]) .+ inner_beven[i] * transpose(cin[2i-1, :]) for i in 1:msing] + else + inner_psi = Vector{Float64}[] + inner_xi = Matrix{ComplexF64}[] + inner_b = Matrix{ComplexF64}[] end - reconnected_flux = Matrix{ComplexF64}(dp.coil .+ transpose(dp.raw) * cout) + reconnected_flux = Matrix{ComplexF64}(dp.ef .+ transpose(dp.raw_fixedboundary) * cout) - # Matched outer ξ per coil: coil column + Σ cout·plasma columns (needs the Galerkin basis). + # Matched outer ξ per EF column: EF column + Σ cout·plasma columns (needs the Galerkin basis). if gal_sol === nothing xi = Array{ComplexF64,3}(undef, 0, 0, 0) xi_deriv = Array{ComplexF64,3}(undef, 0, 0, 0) @@ -205,9 +187,9 @@ function _compute_match(prob::MatchProblem, model::InnerLayer.GGJModel) sols_deriv = gal_sol.xi_deriv mpert = size(sols, 1) ngrid = size(sols, 2) - xi = zeros(ComplexF64, mpert, ngrid, mcoil) - xi_deriv = zeros(ComplexF64, mpert, ngrid, mcoil) - for j in 1:mcoil + xi = zeros(ComplexF64, mpert, ngrid, nef) + xi_deriv = zeros(ComplexF64, mpert, ngrid, nef) + for j in 1:nef csol = 2msing + j @views xi[:, :, j] .= sols[:, :, csol] @views xi_deriv[:, :, j] .= sols_deriv[:, :, csol] @@ -247,16 +229,16 @@ function _compute_match(prob::MatchProblem, model::InnerLayer.GGJModel) inner_xi[i] = inner_xi[i][inside, :] inner_b[i] = inner_b[i][inside, :] - # Cut outer background of the resonant harmonic, per coil. - cutmn = Matrix{ComplexF64}(undef, length(psi_keep), mcoil) - for j in 1:mcoil + # Cut outer background of the resonant harmonic, per EF column. + cutmn = Matrix{ComplexF64}(undef, length(psi_keep), nef) + for j in 1:nef @views cutmn[:, j] .= sols_cut[ires, keep, 2msing+j] for isol in 1:2msing @views cutmn[:, j] .+= cout[isol, j] .* sols_cut[ires, keep, isol] end end itp = cubic_interp(psi_keep, Series(cutmn); extrap=ExtendExtrap()) - buf = Vector{ComplexF64}(undef, mcoil) + buf = Vector{ComplexF64}(undef, nef) hint = Ref(1) for (ip, psi_p) in enumerate(inner_psi[i]) itp(buf, psi_p; hint=hint) @@ -269,19 +251,28 @@ function _compute_match(prob::MatchProblem, model::InnerLayer.GGJModel) end end - return MatchResult(cout, cin, deltar, rpec_eig, residual, bpen, reconnected_flux, + return MatchResult(cout, cin, deltar, forced_eigenvalue, residual, bpen, reconnected_flux, xi, xi_deriv, inner_psi, inner_xi, inner_b, inner_params, model, prob) end -# Copy of ffs with the match applied; the ideal reference keeps closure :ideal, its zero bpen and δW. -function _matched_result(ffs::ForceFreeStatesResult, match::MatchResult, ideal::Bool) +# Copy of ffs with the match applied. A Galerkin match publishes the matched ξ, its edge response (the matched +# columns are EF + plasma·cout, so their edge response is too) and the free-boundary energies of that response with +# the solve's own vacuum. For Riccati match these remain TODO. +function _matched_result(ffs::ForceFreeStatesResult, match::MatchResult) gal = ffs.galerkin - solution = (gal !== nothing && gal.solution !== nothing && !isempty(match.xi)) ? - _matched_gal_profiles(gal, match, ffs.mats, ffs) : ffs.solution - ideal && return _with(ffs; match=match, solution=solution) - # The ideal δW does not describe the matched plasma; PE gates on its absence. - return _with(ffs; closure=:matched, bpen=match.bpen, match=match, solution=solution, - wp=nothing, free_boundary=nothing) + solution = nothing + if gal !== nothing && gal.solution !== nothing && !isempty(match.xi) + solution = _gal_profiles(gal.solution, match.xi, match.xi_deriv, ffs.mats, ffs) + end + # Matched edge momentum: the ideal EF response plus the big-solution columns weighted by cout. + wp = nothing + free_boundary = nothing + dp = ffs.delta_prime + if ffs.wp !== nothing && !isempty(dp.edge_response) + wp = ffs.wp + dp.edge_response * match.cout ./ ffs.equil.psio^2 + ffs.free_boundary === nothing || (free_boundary = free_run(wp, ffs.control, ffs.equil, ffs.mats, ffs)) + end + return _with(ffs; closure=:matched, bpen=match.bpen, match=match, solution=solution, wp=wp, free_boundary=free_boundary) end # A copy of the immutable `x` with the named fields replaced. diff --git a/src/ForceFreeStates/Matching/ResonantMatch.jl b/src/ForceFreeStates/Matching/ResonantMatch.jl index 3c542f332..80e8a74e5 100644 --- a/src/ForceFreeStates/Matching/ResonantMatch.jl +++ b/src/ForceFreeStates/Matching/ResonantMatch.jl @@ -6,17 +6,17 @@ MatchResult Matching coefficients and matched profiles, found in `result.match`. The profile -fields are empty when the solve kept no Galerkin basis; the inner fields are empty with `ideal`. +fields are empty when the solve kept no Galerkin basis. ## Fields - - `cout`, `cin`: `(2msing, ncoil)` outer and inner coefficients. + - `cout`, `cin`: `(2msing, nef)` outer and inner coefficients. - `deltar`: `(msing, 2)` inner-layer Δ per surface, columns (Δ₊, Δ₋) = (interchange, tearing). - - `rpec_eig`: forced eigenvalue γ = 2πi·n·f per surface. + - `forced_eigenvalue`: forced eigenvalue γ = 2πi·n·f per surface. - `residual`: relative residual of the linear solve. - - `bpen`: `(msing, ncoil)` penetrated resonant field at each layer center. - - `reconnected_flux`: `(2msing, ncoil)` `Δ_coil + Δ_outᵀ·cout`. - - `xi`, `xi_deriv`: `(mpert, ngrid, ncoil)` matched outer ξ and ξ′, one column per coil. + - `bpen`: `(msing, nef)` penetrated resonant field at each layer center. + - `reconnected_flux`: `(2msing, nef)` `Δ_EF + Δ_outᵀ·cout`. + - `xi`, `xi_deriv`: `(mpert, ngrid, nef)` matched outer ξ and ξ′, one column per EF column. - `inner_psi`, `inner_xi`, `inner_b`: per surface, the inner ψ grid and the composite ξ_ψ and b^ψ on it; empty unless the inner model returns real-ψ profiles (`GGJModel(; solver=:galerkin)` or `GGJModel(; real_profiles=true)`). - `inner_params`: per-surface `GGJParameters` used. @@ -27,7 +27,7 @@ struct MatchResult{P} cout::Matrix{ComplexF64} cin::Matrix{ComplexF64} deltar::Matrix{ComplexF64} - rpec_eig::Vector{ComplexF64} + forced_eigenvalue::Vector{ComplexF64} residual::Float64 bpen::Matrix{ComplexF64} reconnected_flux::Matrix{ComplexF64} @@ -42,13 +42,13 @@ struct MatchResult{P} end # Solve the 4·msing system mat·[cout; cin] = rmat (rmatch match_rpec). -function _match_system(dp_raw::AbstractMatrix, dp_coil::AbstractMatrix, deltar::AbstractMatrix) +function _match_system(dp_raw::AbstractMatrix, dp_ef::AbstractMatrix, deltar::AbstractMatrix) msing = size(dp_raw, 1) ÷ 2 - ncoil = size(dp_coil, 2) + nef = size(dp_ef, 2) mat = zeros(ComplexF64, 4msing, 4msing) - rmat = zeros(ComplexF64, 4msing, ncoil) + rmat = zeros(ComplexF64, 4msing, nef) @views mat[(2msing+1):4msing, 1:2msing] .= transpose(dp_raw) # Δ_out - @views rmat[(2msing+1):4msing, :] .= .-dp_coil # −Δ_coil + @views rmat[(2msing+1):4msing, :] .= .-dp_ef # −Δ_EF for ising in 1:msing idx1 = 2ising - 1 idx2 = 2ising @@ -82,7 +82,7 @@ function save_to_hdf5(out_h5::Union{HDF5.File,HDF5.Group}, m::MatchResult) out_h5["$g/Delta_r"] = m.deltar out_h5["$g/bpen"] = m.bpen out_h5["$g/reconnected_flux"] = m.reconnected_flux - out_h5["$g/rpec_eig"] = m.rpec_eig + out_h5["$g/forced_eigenvalue"] = m.forced_eigenvalue # Per-surface inner-layer ξ_ψ(ψ) (match.f intotsol); ragged grids → one dataset pair per surface. for i in eachindex(m.inner_psi) out_h5["$g/Inner/psi_$i"] = m.inner_psi[i] @@ -107,20 +107,20 @@ function save_to_hdf5(out_h5::Union{HDF5.File,HDF5.Group}, m::MatchResult) return nothing end -# Matching quantities are per unit coil-response column (the Δ′ coil block's normalization). +# Matching quantities are per unit EF column (the Δ′ EF block's normalization). const _MATCH_SURF = (; dims=("surface",)) const MATCH_H5_ANNOTATIONS = [ "SingularSurfaces/Match/cout" => - (; long_name="outer-region matching coefficients per surface parity, per unit coil-response column", dims=("surface_parity", "coil")), + (; long_name="outer-region matching coefficients per surface parity, per unit EF column", dims=("surface_parity", "ef")), "SingularSurfaces/Match/cin" => - (; long_name="inner-region matching coefficients per surface parity, per unit coil-response column", dims=("surface_parity", "coil")), + (; long_name="inner-region matching coefficients per surface parity, per unit EF column", dims=("surface_parity", "ef")), "SingularSurfaces/Match/Delta_r" => (; long_name="inner-layer Δ per surface and parity, (Δ₊ interchange, Δ₋ tearing)", dims=("surface", "parity")), "SingularSurfaces/Match/bpen" => - (; long_name="penetrated resonant field b^ψ at each layer center, per unit coil-response column", dims=("surface", "coil")), + (; long_name="penetrated resonant field b^ψ at each layer center, per unit EF column", dims=("surface", "ef")), "SingularSurfaces/Match/reconnected_flux" => - (; long_name="reconnected flux Δ_coil + Δ_outᵀ·cout per surface parity, per unit coil-response column", dims=("surface_parity", "coil")), - "SingularSurfaces/Match/rpec_eig" => (; long_name="forced eigenvalue γ = 2πi·n·f at each surface", units="rad/s", _MATCH_SURF...), + (; long_name="reconnected flux Δ_EF + Δ_outᵀ·cout per surface parity, per unit EF column", dims=("surface_parity", "ef")), + "SingularSurfaces/Match/forced_eigenvalue" => (; long_name="forced eigenvalue γ = 2πi·n·f at each surface", units="rad/s", _MATCH_SURF...), "SingularSurfaces/Match/residual" => (; long_name="relative residual of the matching linear solve"), "SingularSurfaces/Match/InnerParams/E" => (; long_name="Glasser-Greene-Johnson coefficient E per surface", _MATCH_SURF...), "SingularSurfaces/Match/InnerParams/F" => (; long_name="Glasser-Greene-Johnson coefficient F per surface", _MATCH_SURF...), @@ -134,6 +134,6 @@ const MATCH_H5_ANNOTATIONS = [ ] const _MATCH_INNER_H5_ANNOTATIONS = [ "psi" => (; long_name="normalized poloidal flux ψ_N grid across the inner layer of this surface", dims=("inner_psi",)), - "xi" => (; long_name="composite inner-layer ξ^ψ per unit coil-response column", dims=("inner_psi", "coil")), - "b" => (; long_name="composite inner-layer b^ψ per unit coil-response column", dims=("inner_psi", "coil")) + "xi" => (; long_name="composite inner-layer ξ^ψ per unit EF column", dims=("inner_psi", "ef")), + "b" => (; long_name="composite inner-layer b^ψ per unit EF column", dims=("inner_psi", "ef")) ] diff --git a/src/ForceFreeStates/Output.jl b/src/ForceFreeStates/Output.jl index ba8480148..0beea4dea 100644 --- a/src/ForceFreeStates/Output.jl +++ b/src/ForceFreeStates/Output.jl @@ -180,22 +180,19 @@ function save_to_hdf5( # Galerkin is the in-domain in-band subset of `SingularSurfaces/`. dp = result.delta_prime if msing > 0 && dp !== nothing - # Inter-surface Δ' matrix, shape [msing × msing] — PEST3-convention deltap. - out_h5["SingularSurfaces/Delta_prime_matrix"] = dp.matrix - - # Edge coil-response matrix, stored (numpert_total × 2msing) = (edge mode, surface-side) - # so H5Web heatmaps read x = edge mode, y = surface-side. The carried matrix stays - # (2msing × numpert_total); transpose only at write. - isempty(dp.coil) || (out_h5["SingularSurfaces/Delta_coil"] = permutedims(dp.coil)) - - # Raw 2msing×2msing outer-region D' matrix, side-major [L_s1, R_s1, L_s2, R_s2, …]; pest3_decompose recovers (A', B', Γ', Δ'). - isempty(dp.raw) || (out_h5["SingularSurfaces/Delta_prime_raw"] = dp.raw) - - # Remaining PEST-3 parity blocks, when the formalism persisted them (Galerkin); - # Riccati recovers them from Delta_prime_raw via pest3_decompose. - dp.A === nothing || (out_h5["SingularSurfaces/pest3_A"] = dp.A) - dp.B === nothing || (out_h5["SingularSurfaces/pest3_B"] = dp.B) - dp.Gamma === nothing || (out_h5["SingularSurfaces/pest3_Gamma"] = dp.Gamma) + ss = "SingularSurfaces" + isempty(dp.raw_freeboundary) || (out_h5["$ss/Delta_prime_raw_freeboundary"] = dp.raw_freeboundary) + isempty(dp.raw_fixedboundary) || (out_h5["$ss/Delta_prime_raw_fixedboundary"] = dp.raw_fixedboundary) + # Δ_EF stored (edge mode × surface-side) so H5Web heatmaps read x = edge mode. + isempty(dp.ef) || (out_h5["$ss/Delta_EF"] = permutedims(dp.ef)) + isempty(dp.edge_response) || (out_h5["$ss/edge_response"] = dp.edge_response) + for (blocks, suffix) in ((dp.freeboundary, ""), (dp.fixedboundary, "_fixedboundary")) + blocks === nothing && continue + out_h5["$ss/Delta_prime_matrix$suffix"] = blocks.Delta + out_h5["$ss/pest3_A$suffix"] = blocks.A + out_h5["$ss/pest3_B$suffix"] = blocks.B + out_h5["$ss/pest3_Gamma$suffix"] = blocks.Gamma + end end result.match === nothing || save_to_hdf5(out_h5, result.match) @@ -213,7 +210,7 @@ function save_to_hdf5( # Jacobian-invariant); eigenmode columns are sorted most-unstable first, unit power norm, largest entry real-positive. fbs = "ForceFreeStates/FreeBoundaryStability" out_h5["$fbs/W_freeboundary"] = free_energies !== nothing ? free_energies.wt0 : ComplexF64[] - out_h5["$fbs/W_plasma"] = free_energies !== nothing ? free_energies.wp : ComplexF64[] + out_h5["$fbs/W_plasma"] = result.wp !== nothing ? result.wp : ComplexF64[] out_h5["$fbs/W_vacuum"] = free_energies !== nothing ? free_energies.wv : ComplexF64[] out_h5["$fbs/W_freeboundary_eigenmodes"] = free_energies !== nothing ? free_energies.wt : ComplexF64[] out_h5["$fbs/eigenmode_energies"] = free_energies !== nothing ? free_energies.et : ComplexF64[] diff --git a/src/ForceFreeStates/Result.jl b/src/ForceFreeStates/Result.jl index 2d2affa45..539032ec8 100644 --- a/src/ForceFreeStates/Result.jl +++ b/src/ForceFreeStates/Result.jl @@ -8,7 +8,7 @@ formalism produced the solution. ## Fields - `basis::Symbol` - Which formalism's basis the profiles are in: `:el_axis` (forward - integrator, axis Euler-Lagrange basis) or `:gal_native` (inner-layer-matched Galerkin + integrator, axis Euler-Lagrange basis) or `:gal_native` (ideal or inner-layer-matched Galerkin solution, identity-at-edge basis). - `step::Int` - Number of stored radial nodes. - `psi_store::Vector{Float64}` - ψ at each node. @@ -40,9 +40,11 @@ Optional fields are `Union{Nothing,T}` and follow a presence-equals-capability r consumer that needs one gates on [`require`](@ref) (or [`require_solution`](@ref) for the ξ profiles) and warns-and-skips rather than erroring, so a result from an integrator that cannot supply a given product still flows through the pipeline. PerturbedEquilibrium refuses -Galerkin and `:matched` results instead, since they publish no `wp` or `free_boundary` yet. +`:matched` results, whose singular coupling is not defined yet, and Galerkin results without `free_boundary`. Every final result records how its basis crosses the rational surfaces, so bpen and closure are always present. +A Forward or Riccati solve applies the ideal jump condition at the rationals; a Galerkin solve +is ideal-closed by its EF columns (`error_field_response`), and is unclosed without them. ## Fields @@ -58,16 +60,19 @@ Every final result records how its basis crosses the rational surfaces, so bpen - `kinetic::NamedTuple` - Kinetic singular-surface scan (`kmsing`, `kinsing`, `scan_psi`, `scan_cond`, `scan_threshold`); empty unless the finder ran. - `closure::Symbol` - How the basis crosses the rationals: `:ideal` (the ideal jump condition was imposed), `:regular` (kinetic terms regularize the force balance equations so no closure is needed) or - `:matched` (an inner-layer solution was matched in). A `:matched` result carries no ideal δW: `wp` and `free_boundary` are `nothing`. + `:matched` (an inner-layer solution was matched in). A `:matched` result carries the matched edge response as `wp`. - `bpen::Matrix{ComplexF64}` - `(msing × numpert_total)` penetrated resonant field from the inner layer, per surface (rows follow `surfaces`) and driving mode; all zeros under `:ideal` and `:regular` closure. - `match::Union{Nothing,MatchResult}` - The [`MatchResult`](@ref) of a `MatchProblem` solve; `nothing` otherwise. - `solution::Union{Nothing,SolutionProfiles}` - THE solve's ξ solution; `nothing` when the - formalism produced none (Riccati, and Galerkin without a match). + formalism produced none (Riccati, and Galerkin without `error_field_response`). - `diagnostics::Union{Nothing,OdeState}` - The integrator's raw ODE state (ψ trace, `crit`, edge scan, asymptotic coefficients). Written to HDF5; not a consumable solution. - - `wp::Union{Nothing,Matrix{ComplexF64}}` - Fixed-boundary plasma energy matrix `W_p = (U₂·U₁⁻¹)/ψ₀²` at - `psilim`. Present for any Euler-Lagrange sweep regardless of `vac_flag`; equals `free_boundary.wp` when that ran. + - `wp::Union{Nothing,Matrix{ComplexF64}}` - Fixed-boundary plasma energy matrix at `psilim`, the edge + displacement → edge momentum map over ψ₀² [Glasser 2016 PoP 23 072505, Eq. 66]: `(U₂·U₁⁻¹)/ψ₀²` for an + Euler-Lagrange sweep (regardless of `vac_flag`; equals `free_boundary.wp` when that ran), the EF-column + edge response over ψ₀² for Galerkin. Under `:matched` closure it is the complex matched edge response, whose + anti-Hermitian part is the torque the inner layers absorb [Logan 2013, eq. 19]. - `free_boundary::Union{Nothing,FreeBoundaryResult}` - Free-boundary energies and eigenmodes; `nothing` when the vacuum step was skipped or the formalism computes none. - `delta_prime::Union{Nothing,DeltaPrimeData}` - Δ′/outer-region matching payload from whichever formalism ran (Riccati BVP or Galerkin); `nothing` when neither produced one. @@ -143,17 +148,22 @@ function require_solution(result::ForceFreeStatesResult, calc::AbstractString) return false end -# Pack the matched Galerkin solution into `SolutionProfiles`: on-surface `issing` points are dropped, Ξ′ is the -# analytic Galerkin derivative, and Ξ_s = −A⁻¹(B·Ξ′ + C·Ξ) as in `sing_der!`. Grid runs inner→edge. -function _matched_gal_profiles(gal_result::GalerkinResult, m::MatchResult, mats::MatrixSplines, intr::ModeSpace) - sol = gal_result.solution +# Columns of the Galerkin error-field (EF) drives (error_field_response). +gal_ef_columns(gal::GalerkinResult) = (2gal.msing+1):size(gal.edge_response, 2) + +# Plasma energy matrix W_P of the Galerkin EF columns, in the units of `(U₂·U₁⁻¹)/ψ₀²`. +gal_wp(gal::GalerkinResult, psio::Float64) = gal.edge_response[:, gal_ef_columns(gal)] ./ psio^2 + +# `SolutionProfiles` of closed Galerkin columns `xi`, `xi_deriv` (npert, ngrid, ncol) on the grid of `sol`, without on-surface points. +function _gal_profiles(sol::GalerkinSolution, xi::AbstractArray{ComplexF64,3}, xi_deriv::AbstractArray{ComplexF64,3}, mats::MatrixSplines, + intr::ModeSpace) npert = intr.numpert_total keep = .!sol.issing psi_f = sol.psi[keep] q_f = sol.q[keep] - xi_f = m.xi[:, keep, :] # (npert, ngrid_f, mcoil) - dxi_f = m.xi_deriv[:, keep, :] + xi_f = xi[:, keep, :] # (npert, ngrid_f, ncol) + dxi_f = xi_deriv[:, keep, :] ngrid_f = length(psi_f) # u_store[:, :, 2] stays zero: the conjugate momentum has no consumer on this path and no @@ -179,8 +189,9 @@ end -> ForceFreeStatesResult Assemble the published result once the solve is finished, with the ideal or regular closure (a -`MatchProblem` solve transforms it afterwards). Materializes the forward path's derivative stores and forms `W_p` -when the free-boundary stage did not run; everything else is copied or aliased. +`MatchProblem` solve transforms it afterwards). Materializes the forward path's derivative stores, packs the +Galerkin EF columns as the solution, and forms `W_p` when the free-boundary stage did not run; everything else +is copied or aliased. `odet` is the ODE state (`nothing` for Galerkin); `gal_data`/`gal_dp` are the Galerkin internals and its Δ′ payload (`nothing` otherwise). @@ -197,22 +208,27 @@ function build_result( gal_data::Union{Nothing,GalerkinResult}, gal_dp::Union{Nothing,DeltaPrimeData} ) + # The Galerkin EF columns (ξ = eⱼ at the edge, no big-solution drive) are its ideal closure. + gal_ef = gal_data === nothing || gal_data.solution === nothing ? (1:0) : gal_ef_columns(gal_data) + # The forward sweep is the only formalism whose stores need materializing; doing it here # keeps `SolutionProfiles.du_store`/`xi_s_store` populated by construction. solution = if integrator === :forward && odet !== nothing materialize_derivative_stores!(odet, equil, mats, intr) SolutionProfiles(:el_axis, odet.step, odet.psi_store, odet.q_store, odet.u_store, odet.du_store, odet.xi_s_store) + elseif !isempty(gal_ef) + gs = gal_data.solution + _gal_profiles(gs, @view(gs.xi[:, :, gal_ef]), @view(gs.xi_deriv[:, :, gal_ef]), mats, intr) else nothing end - # One Δ′ payload whichever formalism produced it; Riccati leaves the PEST-3 blocks unfilled - # because it persists only the raw D′ and recovers them on demand via `pest3_decompose`. + # One Δ′ payload whichever formalism produced it. delta_prime = if gal_dp !== nothing gal_dp - elseif !isempty(intr.delta_prime_matrix) - DeltaPrimeData(intr.delta_prime_matrix, intr.delta_prime_raw, intr.delta_coil, nothing, nothing, nothing) + elseif !isempty(intr.delta_prime_raw_freeboundary) || !isempty(intr.delta_prime_raw_fixedboundary) + DeltaPrimeData(intr.delta_prime_raw_freeboundary, intr.delta_prime_raw_fixedboundary, intr.delta_ef, intr.delta_prime_edge_response) else nothing end @@ -228,6 +244,8 @@ function build_result( free_energies.wp elseif odet !== nothing (odet.u[:, :, 2] / odet.u[:, :, 1]) ./ equil.psio^2 + elseif !isempty(gal_ef) + gal_wp(gal_data, equil.psio) else nothing end diff --git a/src/ForceFreeStates/Riccati/Crossings.jl b/src/ForceFreeStates/Riccati/Crossings.jl index d0b157d4a..7d9644e92 100644 --- a/src/ForceFreeStates/Riccati/Crossings.jl +++ b/src/ForceFreeStates/Riccati/Crossings.jl @@ -106,7 +106,7 @@ function _capture_right_crossing_data!(odet::OdeState, singp::SingType, sing_asy end # STUB: per-surface ca-based Δ' (not physically valid; see SingType.delta_prime docstring). -# The canonical Δ' is intr.delta_prime_matrix from compute_delta_prime_matrix!. +# The canonical Δ' is the BVP Δ' (`DeltaPrimeData`) from compute_delta_prime_matrix!. function _stash_per_surface_delta_prime_stub!(odet::OdeState, intr::ForceFreeStatesInternal, ising::Int, ipert_res, sing_asymp_right, equil::Equilibrium.PlasmaEquilibrium, diff --git a/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl b/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl index fcac7b115..fe74aaac2 100644 --- a/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl +++ b/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl @@ -26,22 +26,17 @@ which requires forward crossing chunks. ## Output: PEST3-convention Δ' (deltap) -The raw BVP solution is a 2·msing × 2·msing matrix `dp` with left/right -sub-indices at each surface, stored in `intr.delta_prime_raw`. The PEST3-convention Δ' is -[Chance, PPPL-2527]: - -``` - deltap(i,j) = dp(2i,2j) - dp(2i,2j-1) - dp(2i-1,2j) + dp(2i-1,2j-1) -``` - -stored in `intr.delta_prime_matrix` (msing × msing). +The raw BVP solution is a 2·msing × 2·msing matrix with left/right sub-indices at each surface: +`intr.delta_prime_raw_freeboundary` with the vacuum edge, `intr.delta_prime_raw_fixedboundary` with +ξ(edge) = 0 (from the vacuum-edge matrix with a Dirichlet edge, which also gives `intr.delta_ef`, or +directly when `wv = nothing`). `pest3_decompose` gives the PEST-3 blocks [Chance, PPPL-2527]. ## Limitations Assumes one resonant mode per surface: if any surface carries several (a multi-`n` run where one q satisfies two `(m, n)`, e.g. q = 2 with (2, 1) and (4, 2)), it warns and skips the BVP. -`intr.delta_prime_matrix` is the only physically valid Δ'; the per-surface ca-based +The BVP Δ' is the only physically valid one; the per-surface ca-based `intr.sing[*].delta_prime` is a diagnostic placeholder and need not agree with it. """ function compute_delta_prime_matrix!( @@ -106,18 +101,21 @@ function compute_delta_prime_matrix!( @info "Δ' BVP: nMat=$nMat, rank(M)=$(rank(M)), cond(M)=$(@sprintf("%.2e", cond(M)))" end - # rpec coil-response block: needs the S-axis row layout that `_solve_bvp_edge_coil` assumes - # for its edge rows, and a vacuum edge in the assembled matrix (col_edge in junc_rows). - if use_S_axis && wv !== nothing - intr.delta_coil = _solve_bvp_edge_coil(M, col_edge, msing, N, ipert_all) + # With a conducting-wall edge the edge unknowns are U₂(edge) of each big-solution drive. + raw, edge = _solve_bvp_raw(M, msing, N, nMat, use_S_axis, ipert_all, col_edge, ctrl, debug) + if wv === nothing + intr.delta_prime_raw_fixedboundary = raw + intr.delta_prime_edge_response = edge + else + intr.delta_prime_raw_freeboundary = raw + # The fixed edge and the EF block need the S-axis row layout of the edge rows. + if use_S_axis + Mc = _fixed_edge_bvp(M, msing, N, col_edge) + intr.delta_prime_raw_fixedboundary, intr.delta_prime_edge_response = + _solve_bvp_raw(Mc, msing, N, nMat, use_S_axis, ipert_all, col_edge, ctrl, debug) + (ctrl === nothing || ctrl.algorithm.error_field_response) && (intr.delta_ef = _solve_bvp_edge_ef(Mc, msing, N, ipert_all)) + end end - - deltap, dp_raw_persisted = _solve_bvp_and_combine_pest3( - M, msing, N, nMat, use_S_axis, ipert_all, col_edge, ctrl, debug) - - # The raw side-major D' (Fortran gal_write_delta layout) feeds `pest3_decompose`. - intr.delta_prime_matrix = deltap - intr.delta_prime_raw = dp_raw_persisted end # Column index helpers for the BVP matrix. j is the 1-based singular-surface index, @@ -393,26 +391,34 @@ function _assemble_bvp_S_axis(uShootR::Vector{Matrix{ComplexF64}}, return M, nMat, col_edge end -# Coil-response block for the Eq. (37) edge [Glasser-Kolemen 2018 PoP 25, 032501]: impose the rpec edge BC -# (RDCON `gal_set_boundary`, unit source per poloidal mode) and read the small-solution coefficient at each -# surface side, solving all N modes with one factorization. Returns delta_coil (2·msing × N). -function _solve_bvp_edge_coil(M::Matrix{ComplexF64}, col_edge, msing::Int, N::Int, ipert_all::Vector{Int}) +# The S-axis BVP matrix `M` with its Eq. (38) vacuum-edge rows [Glasser-Kolemen 2018 PoP 25, 032501] replaced by the +# conducting-wall edge of the `wv = nothing` assembly: U₁(edge) equals the right-hand side of the first N edge rows +# (0 for the fixed-boundary D′) and the edge unknowns hold U₂(edge). +function _fixed_edge_bvp(M::Matrix{ComplexF64}, msing::Int, N::Int, col_edge) nMat = size(M, 1) - bot = (nMat-2msing-N+1):(nMat-2msing) # Eq. (38) bottom rows, carrying the W_V block + edge_rows = (nMat-2msing-2N+1):(nMat-2msing) Mc = copy(M) - Mc[bot, :] .= 0 - Mc[bot, col_edge] .= I(N) # Dirichlet edge: the edge coefficients equal the source + Mc[edge_rows, col_edge] .= 0 + Mc[edge_rows[(N+1):end], col_edge] .= -I(N) + return Mc +end + +# Δ_EF (2·msing × N) from the fixed-edge BVP matrix `Mc`: the small-solution coefficient at each surface side for a +# unit edge displacement in each poloidal mode. +function _solve_bvp_edge_ef(Mc::Matrix{ComplexF64}, msing::Int, N::Int, ipert_all::Vector{Int}) + nMat = size(Mc, 1) + u1_rows = (nMat-2msing-2N+1):(nMat-2msing-N) B = zeros(ComplexF64, nMat, N) - B[bot, :] .= I(N) # unit drive, one column per edge poloidal mode + B[u1_rows, :] .= I(N) # unit edge displacement, one column per edge poloidal mode X = lu(Mc) \ B - delta_coil = zeros(ComplexF64, 2msing, N) + delta_ef = zeros(ComplexF64, 2msing, N) for j in 1:msing row_left = _col_left(j, N)[ipert_all[j]+N] row_right = _col_right(j, N)[ipert_all[j]+N] - @views delta_coil[2j-1, :] .= X[row_left, :] - @views delta_coil[2j, :] .= X[row_right, :] + @views delta_ef[2j-1, :] .= X[row_left, :] + @views delta_ef[2j, :] .= X[row_right, :] end - return delta_coil + return delta_ef end # Fallback BVP assembly with FM-based axis BC (used when no Riccati S matrices are available). @@ -468,8 +474,8 @@ function _assemble_bvp_FM_axis(Phi_L_mats::Vector{Matrix{ComplexF64}}, return M, nMat, col_edge end -# Solve the BVP and apply the PEST3 four-term combination. -function _solve_bvp_and_combine_pest3(M::Matrix{ComplexF64}, msing::Int, N::Int, nMat::Int, +# Solve the BVP for the raw side-major D′ (2·msing × 2·msing) and the edge unknowns of each drive (N × 2·msing). +function _solve_bvp_raw(M::Matrix{ComplexF64}, msing::Int, N::Int, nMat::Int, use_S_axis::Bool, ipert_all::Vector{Int}, col_edge, ctrl, debug::Bool) s2 = 2 * msing @@ -483,6 +489,7 @@ function _solve_bvp_and_combine_pest3(M::Matrix{ComplexF64}, msing::Int, N::Int, end dp_raw = zeros(Tc, s2, s2) + edge = zeros(Tc, N, s2) b = zeros(Tc, nMat) for jsing in 1:msing, side in 1:2 dRow = 2jsing - (2 - side) @@ -499,17 +506,11 @@ function _solve_bvp_and_combine_pest3(M::Matrix{ComplexF64}, msing::Int, N::Int, dp_raw[dRow, 2ksing-1] = x[_col_left(ksing, N)[ipert_k+N]] dp_raw[dRow, 2ksing] = x[_col_right(ksing, N)[ipert_k+N]] end + edge[:, dRow] .= x[col_edge] end - # PEST3 four-term combination [Chance PPPL-2527; Glasser-Kolemen 2018 PoP 25, 032501 Eq. 31]. - # Δ'[i,j] = (NW − NE − SW + SE) on each 2×2 block of dp_raw. - deltap = zeros(Tc, msing, msing) - for i in 1:msing, j in 1:msing - deltap[i, j] = dp_raw[2i, 2j] - dp_raw[2i, 2j-1] - dp_raw[2i-1, 2j] + dp_raw[2i-1, 2j-1] - end - - debug && _log_bvp_pest3(dp_raw, deltap, s2, msing) - return ComplexF64.(deltap), ComplexF64.(dp_raw) + debug && _log_bvp_pest3(dp_raw, pest3_decompose(ComplexF64.(dp_raw)).Delta, s2, msing) + return ComplexF64.(dp_raw), ComplexF64.(edge) end # Logging helpers for `compute_delta_prime_matrix!`. Called only when debug=true. @@ -617,60 +618,3 @@ function _log_bvp_pest3(dp_raw, deltap, s2, msing) @info "Δ' BVP: Raw dp diagonal = $([@sprintf("%.4f%+.4fi", Float64(real(dp_raw[i,i])), Float64(imag(dp_raw[i,i]))) for i in 1:s2])" @info "Δ' BVP: deltap diagonal = $([@sprintf("%.4f%+.4fi", real(deltap[i,i]), imag(deltap[i,i])) for i in 1:msing])" end - -""" - pest3_decompose(dp_raw::AbstractMatrix) -> NamedTuple{(:A, :B, :Γ, :Δ)} - -Rotate the raw 2m×2m outer-region matching matrix `dp_raw` (side-major -ordering `[L_s1, R_s1, L_s2, R_s2, …]`) into the Pletzer–Dewar 1991 parity -blocks. Given rows and columns paired by surface (odd index = left, even -index = right), the Fortran RDCON parity combination is - -``` -A'(i,j) = RR + RL + LR + LL (even-i, even-j) — interchange↔interchange -B'(i,j) = RR − RL + LR − LL (even-i, odd-j) — interchange↔tearing -Γ'(i,j) = RR + RL − LR − LL (odd-i, even-j) — tearing↔interchange -Δ'(i,j) = RR − RL − LR + LL (odd-i, odd-j) — tearing↔tearing -``` - -where `RR = dp_raw[2i, 2j]`, `RL = dp_raw[2i, 2j−1]`, -`LR = dp_raw[2i−1, 2j]`, `LL = dp_raw[2i−1, 2j−1]`. Each block is m×m. - -No Pletzer–Dewar ½ prefactor, matching Fortran RDCON; it cancels in `det(D' − D(γ)) = 0`. - -The Δ' block returned here equals `intr.delta_prime_matrix` (the m×m PEST3 -tearing projection computed inside `compute_delta_prime_matrix!`). - -# Arguments - - - `dp_raw` — 2m×2m complex matrix (typically `intr.delta_prime_raw`). - -# Returns - -Named tuple `(A, B, Γ, Δ)` of m×m matrices, the blocks of `D' = [[A' B'] [Γ' Δ']]` in the -`det(D' − D(γ)) = 0` problem (interchange upper-left, tearing lower-right). -""" -function pest3_decompose(dp_raw::AbstractMatrix) - s2 = size(dp_raw, 1) - size(dp_raw, 2) == s2 || - throw(ArgumentError("pest3_decompose: dp_raw must be square, got $(size(dp_raw))")) - iseven(s2) || - throw(ArgumentError("pest3_decompose: dp_raw side must be 2m for integer m, got $s2")) - m = s2 ÷ 2 - Tc = eltype(dp_raw) - Ap = zeros(Tc, m, m) - Bp = zeros(Tc, m, m) - Gp = zeros(Tc, m, m) - Dp = zeros(Tc, m, m) - for i in 1:m, j in 1:m - LL = dp_raw[2i-1, 2j-1] - LR = dp_raw[2i-1, 2j] - RL = dp_raw[2i, 2j-1] - RR = dp_raw[2i, 2j] - Ap[i, j] = RR + RL + LR + LL - Bp[i, j] = RR - RL + LR - LL - Gp[i, j] = RR + RL - LR - LL - Dp[i, j] = RR - RL - LR + LL - end - return (A=Ap, B=Bp, Γ=Gp, Δ=Dp) -end diff --git a/src/ForceFreeStates/Run.jl b/src/ForceFreeStates/Run.jl index 7e440b1f5..f004b72d5 100644 --- a/src/ForceFreeStates/Run.jl +++ b/src/ForceFreeStates/Run.jl @@ -56,8 +56,8 @@ function run_force_free_states( nstring = _mode_range_label(intr) # The three formalisms are exclusive. Galerkin solves the same Euler-Lagrange system - # variationally rather than by radial ODE integration, so it replaces both the integration - # and the free-boundary energies, and supplies its own vacuum response at psilim. + # variationally rather than by radial ODE integration, and supplies its own vacuum response at + # psilim; its EF columns (error_field_response) give the plasma response the free-boundary energies need. odet = nothing free_energies = nothing gal_data = nothing @@ -70,6 +70,9 @@ function run_force_free_states( gal_data, gal_dp = galerkin_solve(ctrl, equil, mats, intr; wv=wv, cut_solution=cut_solution) gal_dt = time() - gal_start @info "Galerkin solve completed in $(@sprintf("%.3f", gal_dt)) s" + if ctrl.vac_flag && gal_data.solution !== nothing && !isempty(gal_ef_columns(gal_data)) + free_energies = free_run(gal_wp(gal_data, equil.psio), ctrl, equil, mats, intr; verbose=intr.verbose) + end else if intr.verbose @info "Integrating Euler-Lagrange equation" @@ -84,7 +87,7 @@ function run_force_free_states( wall_desc = intr.wall_settings === nothing ? "no wall" : nameof(typeof(intr.wall_settings)) @info "Computing free boundary energies ($wall_desc)" end - free_energies = free_run(odet, ctrl, equil, mats, intr) + free_energies = free_run((odet.u[:, :, 2] / odet.u[:, :, 1]) ./ equil.psio^2, ctrl, equil, mats, intr; verbose=intr.verbose) normalize_eigenfunctions!(odet, free_energies.wt, equil.psio) if real(free_energies.et[1]) < 0 if intr.verbose diff --git a/src/ForceFreeStates/Surfaces/Types.jl b/src/ForceFreeStates/Surfaces/Types.jl index 0099590ab..447cf5a65 100644 --- a/src/ForceFreeStates/Surfaces/Types.jl +++ b/src/ForceFreeStates/Surfaces/Types.jl @@ -14,7 +14,7 @@ A mutable struct holding data related to the singular surfaces in the equilibriu - `q::Float64` - Safety factor (= m/n) - `q1::Float64` - Derivative of safety factor with respect to ψ - `delta_prime::Vector{ComplexF64}` - **Debug stub, not physically valid.** Per-surface ca-based Δ' estimate; - use `ForceFreeStatesInternal.delta_prime_matrix` (STRIDE global BVP, Glasser 2018 PoP 25, 032501) for tearing stability. + use the BVP Δ' (`DeltaPrimeData`) (STRIDE global BVP, Glasser 2018 PoP 25, 032501) for tearing stability. - `delta_prime_col::Matrix{ComplexF64}` - **Debug stub, not physically valid.** Per-surface ca-based Δ' columns, `(numpert_total, n_res_modes)`, `delta_prime_col[j, i] = (ca_r[j,ipert_res_i,2] - ca_l[j,ipert_res_i,2]) / (4π²·psio)`; diagonal matches `delta_prime`. Riccati/parallel paths only. """ diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl index 2beaa5916..259564589 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl @@ -51,9 +51,10 @@ their dominant modes (Park, Phys. Plasmas 16, 056115 (2009)), and with `compute_ with `perturbed_equilibrium(resp, rmp)` as many times as needed. Products the producing integrator could not supply (free-boundary energies, a ξ solution) leave the -corresponding operators empty: the step warns and is skipped. Galerkin and `:matched` results publish -no free-boundary δW yet, so building a response from one is an error. So is a solve with a wall: the -surface inductance is free-space, and the response needs Λ from the same no-wall vacuum. +corresponding operators empty: the step warns and is skipped. A Galerkin result carries a free-boundary +δW only with `error_field_response` and `vac_flag`, so building a response from one without it is an error. So is a +`:matched` result, whose singular coupling is not defined yet, and a solve with a wall: the surface +inductance is free-space, and the response needs Λ from the same no-wall vacuum. """ function plasma_response(ffs::ForceFreeStatesResult; compute_singular_coupling::Bool=true, compute_torque_response::Bool=false, verbose::Bool=false)::PlasmaResponse @@ -62,11 +63,14 @@ function plasma_response(ffs::ForceFreeStatesResult; compute_singular_coupling:: "plasma_response needs a force-free-states solve without a wall: its surface inductance L is free-space, " * "so Λ from a wall run would mix two vacuum models. Solve with wall=nothing (no [Wall] section)." ) - # Galerkin and matched results publish no plasma edge response wp, so no operator could be built. - if ffs.free_boundary === nothing && (ffs.integrator === :galerkin || ffs.closure === :matched) - error("plasma_response needs the free-boundary δW, which a $(ffs.integrator) result with $(ffs.closure) closure " * - "does not publish yet; use a forward solve") - end + # The singular coupling evaluates ideal-closure resonant metrics off each surface, where a matched solution's + # big solution is singular. + ffs.closure === :matched && + error("plasma_response does not support a matched result yet: its singular-coupling metrics assume the ideal closure; " * + "use the unmatched solve") + ffs.free_boundary === nothing && ffs.integrator === :galerkin && + error("plasma_response needs the free-boundary δW, which this Galerkin result does not carry: solve with " * + "Galerkin(; error_field_response=true) and vac_flag=true") m_modes, n_modes = mode_arrays(ffs) # A run has at most one ξ solution, already closed at the rationals and carrying populated diff --git a/src/Tearing/Dispersion/CoupledFullMatch.jl b/src/Tearing/Dispersion/CoupledFullMatch.jl index 19e7e52ba..91e3c877b 100644 --- a/src/Tearing/Dispersion/CoupledFullMatch.jl +++ b/src/Tearing/Dispersion/CoupledFullMatch.jl @@ -11,7 +11,7 @@ # Eq. (36)-(40); original matching construction Pletzer & Dewar 1991, # J. Plasma Phys. **45**, 427). # -# The outer D' (`intr.delta_prime_raw`) is in the side-major L/R basis while the inner layer returns +# The outer D' (`DeltaPrimeData` raw D') is in the side-major L/R basis while the inner layer returns # Δ in the parity ± basis, so both amplitude sets are kept as unknowns and joined by the matching # identity rows, giving a 4m × 4m system instead of a basis-mixing 2m × 2m determinant. # @@ -64,7 +64,7 @@ end inner_kwargs=NamedTuple()) -> MultiSurfaceCouplingFull Construct the 4m × 4m dispersion matrix driver. `dp_raw` must be the -2m × 2m matrix in side-major ordering (the `intr.delta_prime_raw` +2m × 2m matrix in side-major ordering (the `DeltaPrimeData` raw D' field populated by `ForceFreeStates.compute_delta_prime_matrix!` on the Riccati path). Surface k's inner layer is evaluated at `Q·ratio_k` Doppler-shifted by its E×B offset `q_shift_k`, exactly as in the reduced `multi_surface_coupling`: added to Q for a diff --git a/src/Tearing/GGJTearing.jl b/src/Tearing/GGJTearing.jl index c6780aac9..b075c8398 100644 --- a/src/Tearing/GGJTearing.jl +++ b/src/Tearing/GGJTearing.jl @@ -4,18 +4,25 @@ # determinant (Dispersion.ggj_growth_rates), each surface on its own when uncoupled, # the leading `msing_max` surfaces together when coupled. -# The raw side-major Δ′ the determinant needs, after checking it is aligned with ordered surfaces. -function _ggj_delta_prime_raw(ffs) +# The raw side-major Δ′ of the problem's boundary that the determinant needs, aligned with ordered surfaces. +function _ggj_delta_prime_raw(prob::TearingProblem) + ffs = prob.ffs ffs.npert == 1 || throw(ArgumentError("GGJ tearing is single-n; the result spans n = $(ffs.nlow):$(ffs.nhigh)")) surfaces = ffs.surfaces psi = Float64[s.psifac for s in surfaces] issorted(psi; lt=<=) || throw(ArgumentError("GGJ tearing needs distinct surfaces ordered core to edge; got ψ = $psi")) dp = ffs.delta_prime - dp === nothing && throw(ArgumentError("GGJ tearing needs the raw Δ′ matrix; solve with Riccati() and vac_flag=true, or Galerkin(; rpec_flag=true)")) - size(dp.raw) == (2length(surfaces), 2length(surfaces)) || - throw(ArgumentError("GGJ tearing: Δ′ raw block is $(size(dp.raw, 1))×$(size(dp.raw, 2)) but the result carries $(length(surfaces)) surfaces")) - all(isfinite, dp.raw) || throw(ArgumentError("GGJ tearing: the raw Δ′ matrix has non-finite entries")) - return dp.raw + dp === nothing && throw(ArgumentError("GGJ tearing needs the raw Δ′ matrix; solve with Riccati() or Galerkin()")) + raw = prob.fixed_boundary ? dp.raw_fixedboundary : dp.raw_freeboundary + size(raw) == (2length(surfaces), 2length(surfaces)) || + throw( + ArgumentError( + "GGJ tearing: the $(prob.fixed_boundary ? "fixed" : "free")-boundary Δ′ raw block is " * + "$(size(raw, 1))×$(size(raw, 2)) but the result carries $(length(surfaces)) surfaces" + ) + ) + all(isfinite, raw) || throw(ArgumentError("GGJ tearing: the raw Δ′ matrix has non-finite entries")) + return raw end # η and ρ given on the problem, as MatchProblem takes them: a scalar or one value per surface, core to edge. @@ -70,7 +77,7 @@ function CommonSolve.solve(prob::TearingProblem, model::GGJModel; config::Tearin ffs = prob.ffs surfaces = ffs.surfaces isempty(surfaces) && return empty_tearing_result(prob, config, model) - raw = _ggj_delta_prime_raw(ffs) + raw = _ggj_delta_prime_raw(prob) count = length(surfaces) eta, rho = _ggj_layer_inputs(prob) params = [ggj_parameters(s, ffs.equil; eta=eta[k], rho=rho[k], gamma=prob.gamma, ising=k) for (k, s) in enumerate(surfaces)] @@ -118,7 +125,8 @@ function CommonSolve.solve(prob::TearingProblem, model::GGJModel; config::Tearin end end return TearingResult(prob, true, config, model, params, - Float64[s.psifac for s in surfaces], Float64[s.q for s in surfaces], Matrix{ComplexF64}(ffs.delta_prime.matrix), + Float64[s.psifac for s in surfaces], Float64[s.q for s in surfaces], + Matrix{ComplexF64}((prob.fixed_boundary ? ffs.delta_prime.fixedboundary : ffs.delta_prime.freeboundary).Delta), Q_root, omega, Float64[e.gamma for e in extractions], coupled ? GrowthRateResult[] : extractions, coupled ? only(extractions) : nothing, LayerWidths[], Ω_E, zeros(count), scans, empty_critical_resonant_field_result()) diff --git a/src/Tearing/HDF5Output.jl b/src/Tearing/HDF5Output.jl index dbf1ed5b9..e64081e82 100644 --- a/src/Tearing/HDF5Output.jl +++ b/src/Tearing/HDF5Output.jl @@ -38,7 +38,8 @@ function save_to_hdf5(parent::Union{HDF5.File,HDF5.Group}, _write_per_surface!(g, result.params, result.dp_matrix) if result.model isa GGJModel && result.problem !== nothing - g["PerSurface/Delta_prime_raw"] = result.problem.ffs.delta_prime.raw + dp = result.problem.ffs.delta_prime + g["PerSurface/Delta_prime_raw"] = result.problem.fixed_boundary ? dp.raw_fixedboundary : dp.raw_freeboundary end # E×B rotation resolved per surface, and the Doppler offset it applies to the inner-layer Q. g["PerSurface/omega_E"] = result.omega_E diff --git a/src/Tearing/TearingProblem.jl b/src/Tearing/TearingProblem.jl index 1d1b5a75d..e62cbc6c4 100644 --- a/src/Tearing/TearingProblem.jl +++ b/src/Tearing/TearingProblem.jl @@ -11,7 +11,7 @@ GGJ layer inputs; the SLAYER layer closures belong to [`SLAYERModel`](@ref) and Z_eff and transport defaults, confinement time and resistivity closure) is the equilibrium's kinetic data, `ffs.equil.kinetic` (see [`KineticProfileData`](@ref)), so every layer model sees one plasma. Solve with `solve(problem, GGJModel(); config)` or `solve(problem, SLAYERModel(); config)`. -GGJ needs the raw Δ′ matrix (Riccati with `vac_flag=true`, or Galerkin with `rpec_flag=true`). +Both read the free-boundary Δ′ of `ffs.delta_prime`, or the fixed-boundary one with `fixed_boundary`. Both models report lab-frame roots from each surface's E×B rotation. ## Fields @@ -24,6 +24,7 @@ Both models report lab-frame roots from each surface's E×B rotation. (core to edge), as in [`MatchProblem`](@ref GeneralizedPerturbedEquilibrium.ForceFreeStates.MatchProblem). `nothing` (default) derives them from the kinetic data: η from its resistivity closure, the one SLAYER uses, and ρ = μ·m_p·n_e as SLAYER takes it. Unused by `SLAYERModel`. - `gamma` - GGJ ratio of specific heats. + - `fixed_boundary` - Use the fixed-boundary Δ′ (a conformal wall on the plasma) instead of the free-boundary one. """ @kwdef struct TearingProblem{R,E,D} ffs::R @@ -34,6 +35,7 @@ Both models report lab-frame roots from each surface's E×B rotation. msing_max::Int = 3 omega_E_kHz::Dict{String,Float64} = Dict{String,Float64}() critical_resonant_field::Union{Nothing,CriticalResonantFieldControl} = nothing + fixed_boundary::Bool = false end TearingProblem(ffs; eta=nothing, rho=nothing, kwargs...) = @@ -90,6 +92,10 @@ function CommonSolve.solve(prob::TearingProblem, model::SLAYERModel; config::Tea validate(prob, config) isempty(surfaces) && return empty_tearing_result(prob, config, model) + dp = ffs.delta_prime + blocks = dp === nothing ? nothing : prob.fixed_boundary ? dp.fixedboundary : dp.freeboundary + dp !== nothing && blocks === nothing && + error("TearingProblem: the solve did not compute the $(prob.fixed_boundary ? "fixed" : "free")-boundary Δ′") equil.kinetic === nothing && error("TearingProblem: the equilibrium carries no kinetic profiles; set kinetic_file in [Equilibrium] or attach_kinetic_profiles!") @@ -99,7 +105,7 @@ function CommonSolve.solve(prob::TearingProblem, model::SLAYERModel; config::Tea params = _tearing_params(model, equil, surfaces, loaded) # Prefer the full inter-surface Δ' matrix; fall back to a diagonal of each surface's scalar delta_prime. - dpm = ffs.delta_prime === nothing ? Matrix{ComplexF64}(undef, 0, 0) : ffs.delta_prime.matrix + dpm = blocks === nothing ? Matrix{ComplexF64}(undef, 0, 0) : blocks.Delta dp = if !isempty(dpm) && size(dpm) == (length(params), length(params)) Matrix{ComplexF64}(dpm) else diff --git a/src/Utilities/HDF5Schema.jl b/src/Utilities/HDF5Schema.jl index f17b53511..e14f2a9fe 100644 --- a/src/Utilities/HDF5Schema.jl +++ b/src/Utilities/HDF5Schema.jl @@ -244,14 +244,25 @@ const MAIN_H5_ANNOTATIONS = [ (; long_name="μ0 × dp/dψ_N at each surface", units="T^2", dims=("surface",), attach=(1 => "SingularSurfaces/rational_psi", 1 => "SingularSurfaces/rational_q")), "SingularSurfaces/dVdpsi" => (; long_name="dV/dψ_N at each surface", units="m^3", dims=("surface",), attach=(1 => "SingularSurfaces/rational_psi", 1 => "SingularSurfaces/rational_q")), - "SingularSurfaces/Delta_prime_matrix" => - (; long_name="inter-surface Δ' matrix (PEST-3 tearing↔tearing block; STRIDE BVP or Galerkin outer region)", dims=("surface_row", "surface_col")), - "SingularSurfaces/Delta_prime_raw" => - (; long_name="raw 2msing×2msing outer-region D' matrix, side-major ordering [L_s1, R_s1, ...]", dims=("surface_side_row", "surface_side_col")), - "SingularSurfaces/Delta_coil" => (; long_name="edge coil-response matrix (edge mode × surface-side)", dims=("mode", "surface_side")), - "SingularSurfaces/pest3_A" => (; long_name="PEST-3 matching block A' (interchange↔interchange)", dims=("surface_row", "surface_col")), - "SingularSurfaces/pest3_B" => (; long_name="PEST-3 matching block B' (interchange↔tearing)", dims=("surface_row", "surface_col")), - "SingularSurfaces/pest3_Gamma" => (; long_name="PEST-3 matching block Γ' (tearing↔interchange)", dims=("surface_row", "surface_col")), + "SingularSurfaces/Delta_prime_raw_freeboundary" => + (; long_name="raw outer-region D' with the vacuum edge, side-major [L_s1, R_s1, ...]", dims=("surface_side_row", "surface_side_col")), + "SingularSurfaces/Delta_prime_raw_fixedboundary" => + (; long_name="raw outer-region D' with a fixed edge (conformal wall on the plasma), side-major [L_s1, R_s1, ...]", + dims=("surface_side_row", "surface_side_col")), + "SingularSurfaces/Delta_EF" => + (; long_name="error-field response Δ_EF: small-solution coefficient per unit edge displacement (edge mode × surface-side)", dims=("mode", "surface_side")), + "SingularSurfaces/edge_response" => + (; long_name="edge momentum Fξ′+Kξ of the fixed-edge solution with a unit big-solution coefficient on each surface side (edge mode × surface-side)", + dims=("mode", "surface_side")), + "SingularSurfaces/Delta_prime_matrix" => (; long_name="free-boundary Δ' matrix (PEST-3 tearing↔tearing block)", dims=("surface_row", "surface_col")), + "SingularSurfaces/pest3_A" => (; long_name="free-boundary PEST-3 block A' (interchange↔interchange)", dims=("surface_row", "surface_col")), + "SingularSurfaces/pest3_B" => (; long_name="free-boundary PEST-3 block B' (interchange↔tearing)", dims=("surface_row", "surface_col")), + "SingularSurfaces/pest3_Gamma" => (; long_name="free-boundary PEST-3 block Γ' (tearing↔interchange)", dims=("surface_row", "surface_col")), + "SingularSurfaces/Delta_prime_matrix_fixedboundary" => + (; long_name="fixed-boundary Δ' matrix (PEST-3 tearing↔tearing block)", dims=("surface_row", "surface_col")), + "SingularSurfaces/pest3_A_fixedboundary" => (; long_name="fixed-boundary PEST-3 block A' (interchange↔interchange)", dims=("surface_row", "surface_col")), + "SingularSurfaces/pest3_B_fixedboundary" => (; long_name="fixed-boundary PEST-3 block B' (interchange↔tearing)", dims=("surface_row", "surface_col")), + "SingularSurfaces/pest3_Gamma_fixedboundary" => (; long_name="fixed-boundary PEST-3 block Γ' (tearing↔interchange)", dims=("surface_row", "surface_col")), # --- SingularSurfaces/Kinetic/ --- "SingularSurfaces/Kinetic/rational_count" => (; long_name="number of kinetic singular surfaces (det(F̄) near-zeros)"), "SingularSurfaces/Kinetic/rational_psi" => (; long_name="normalized poloidal flux ψ_N of kinetic singular surfaces", scale="psi_kinetic_rational"), diff --git a/test/runtests_analysis_plots.jl b/test/runtests_analysis_plots.jl index 77caa1211..794bb9622 100644 --- a/test/runtests_analysis_plots.jl +++ b/test/runtests_analysis_plots.jl @@ -146,7 +146,7 @@ using TOML @test A.ForceFreeStates.plot_ffs_summary(ret.ffs) isa Figure # The Δ′ family draws any part of its complex entries and refuses others. @test A.ForceFreeStates.plot_delta_prime(ret.ffs; part=:abs) isa Figure - @test A.ForceFreeStates.plot_delta_coil(h5; part=:imag) isa Figure + @test A.ForceFreeStates.plot_delta_ef(h5; part=:imag) isa Figure @test_throws ArgumentError A.ForceFreeStates.plot_delta_prime(ret.ffs; part=:phase) end end diff --git a/test/runtests_galerkin_cholesky.jl b/test/runtests_galerkin_cholesky.jl index 11a38764f..e1fff1e8f 100644 --- a/test/runtests_galerkin_cholesky.jl +++ b/test/runtests_galerkin_cholesky.jl @@ -2,7 +2,7 @@ using Test using LinearAlgebra # gal_banded_solve! on a full-band (gbtrf layout) system built like the Galerkin matrix: Hermitian, with -# Dirichlet identity rows (row zeroed, column kept) and nonzero prescribed values, as in the rpec coil columns. +# Dirichlet identity rows (row zeroed, column kept) and nonzero prescribed values, as in the EF columns. @testset "Galerkin banded solve: Cholesky first, LU fallback" begin FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates @test FFS.empty_galerkin_result().solver == "none" diff --git a/test/runtests_parallel_integration.jl b/test/runtests_parallel_integration.jl index 7ddba18d9..f7cba83d9 100644 --- a/test/runtests_parallel_integration.jl +++ b/test/runtests_parallel_integration.jl @@ -29,12 +29,11 @@ end x = factor \ rhs raw[row, :] = x[[3, 5]] end - delta = reshape([raw[2, 2] - raw[2, 1] - raw[1, 2] + raw[1, 1]], 1, 1) - return ComplexF64.(delta), ComplexF64.(raw) + return ComplexF64.(raw) end - default = FFS._solve_bvp_and_combine_pest3(M, 1, 1, 6, false, [1], nothing, nothing, false) - extended = FFS._solve_bvp_and_combine_pest3(M, 1, 1, 6, false, [1], nothing, FFS.ForceFreeStatesControl(), false) - ordinary = FFS._solve_bvp_and_combine_pest3(M, 1, 1, 6, false, [1], nothing, FFS.ForceFreeStatesControl(; algorithm=Riccati(; extended_precision_bvp=false)), false) + default = FFS._solve_bvp_raw(M, 1, 1, 6, false, [1], 6:6, nothing, false)[1] + extended = FFS._solve_bvp_raw(M, 1, 1, 6, false, [1], 6:6, FFS.ForceFreeStatesControl(), false)[1] + ordinary = FFS._solve_bvp_raw(M, 1, 1, 6, false, [1], 6:6, FFS.ForceFreeStatesControl(; algorithm=Riccati(; extended_precision_bvp=false)), false)[1] @test default == reference @test extended == reference @test ordinary != reference @@ -254,7 +253,7 @@ end function run_solovev(integrator) ctrl, intr, equil, metric, mats = _prepared_solovev(integrator) odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, mats, intr) - vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, mats, intr) + vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run((odet.u[:, :, 2] / odet.u[:, :, 1]) ./ equil.psio^2, ctrl, equil, mats, intr) return real(vac.et[1]), intr end diff --git a/test/runtests_result_struct.jl b/test/runtests_result_struct.jl index f38ada7aa..0f3a895d3 100644 --- a/test/runtests_result_struct.jl +++ b/test/runtests_result_struct.jl @@ -110,13 +110,16 @@ include("display_check.jl") @test iszero(ffs.bpen) @test ffs.delta_prime !== nothing - @test size(ffs.delta_prime.matrix) == (length(ffs.surfaces), length(ffs.surfaces)) - @test size(ffs.delta_prime.raw) == (2 * length(ffs.surfaces), 2 * length(ffs.surfaces)) - # Riccati persists only the raw D′ and recovers the parity blocks on demand. - @test ffs.delta_prime.A === nothing - @test ffs.delta_prime.B === nothing - @test ffs.delta_prime.Gamma === nothing - @test ffs.delta_prime.matrix ≈ FFS.pest3_decompose(ffs.delta_prime.raw).Δ + dp = ffs.delta_prime + @test size(dp.raw_freeboundary) == (2 * length(ffs.surfaces), 2 * length(ffs.surfaces)) + @test size(dp.raw_fixedboundary) == size(dp.raw_freeboundary) && dp.raw_fixedboundary != dp.raw_freeboundary + # The PEST-3 blocks of both boundaries are stored, whichever formalism ran. + for (blocks, raw) in ((dp.freeboundary, dp.raw_freeboundary), (dp.fixedboundary, dp.raw_fixedboundary)) + ref = FFS.pest3_decompose(raw) + @test blocks.A ≈ ref.A && blocks.B ≈ ref.B && blocks.Gamma ≈ ref.Gamma && blocks.Delta ≈ ref.Delta + end + @test size(dp.ef) == (2 * length(ffs.surfaces), ffs.numpert_total) + @test size(dp.edge_response) == (ffs.numpert_total, 2 * length(ffs.surfaces)) && !iszero(dp.edge_response) @test ffs.free_boundary !== nothing foreach(_check_display, (ffs, ffs.delta_prime)) @@ -179,9 +182,9 @@ include("display_check.jl") end end - # Standalone Galerkin with the RPEC inner-layer match: the matched outer solution becomes - # the result's ξ solution, and the closure records whether an inner layer was matched in. - @testset "matched Galerkin publishes a gal-native solution" begin + # Standalone Galerkin publishes its ideal-closed EF columns as the ξ solution; with the RPEC + # inner-layer match the matched outer solution replaces them, and the closure records which. + @testset "Galerkin publishes a gal-native solution" begin for (deck, ideal) in (("LAR_ideal_match_test", true), ("LAR_resistive_match_test", false)) @testset "$deck" begin mktempdir() do dir @@ -191,7 +194,7 @@ include("display_check.jl") @test ffs.integrator === :galerkin @test ffs.galerkin !== nothing - @test ffs.match !== nothing + @test (ffs.match === nothing) == ideal sol = ffs.solution @test sol isa FFS.SolutionProfiles @@ -200,44 +203,41 @@ include("display_check.jl") @test !iszero(sol.du_store) # analytic galerkin Ξ′ @test !iszero(sol.xi_s_store) - # Galerkin runs no Euler-Lagrange sweep and computes no free-boundary energies. + # Galerkin runs no Euler-Lagrange sweep; the edge response of its closed columns gives the + # free-boundary energies (vac_flag = true in both decks). @test ffs.diagnostics === nothing - @test ffs.free_boundary === nothing - @test ffs.wp === nothing + @test size(ffs.wp) == (ffs.numpert_total, ffs.numpert_total) + @test ffs.free_boundary !== nothing && ffs.free_boundary.wp == ffs.wp + @test ffs.free_boundary.wt0 ≈ ffs.wp + ffs.free_boundary.wv rtol = 1e-12 # Galerkin publishes its Δ′ through the same unified field the Riccati BVP # uses, in the same PEST-3 convention and the same shared layout. dp = ffs.delta_prime @test dp !== nothing msing = ffs.galerkin.msing - @test size(dp.matrix) == (msing, msing) - @test size(dp.raw) == (2msing, 2msing) - @test dp.matrix ≈ FFS.pest3_decompose(dp.raw).Δ - - # Unlike Riccati, the Galerkin solve persists the remaining parity blocks. - for (blk, ref) in ((dp.A, :A), (dp.B, :B), (dp.Gamma, :Γ)) - @test blk !== nothing - @test size(blk) == (msing, msing) - @test blk ≈ getfield(FFS.pest3_decompose(dp.raw), ref) + for (blocks, raw) in ((dp.freeboundary, dp.raw_freeboundary), (dp.fixedboundary, dp.raw_fixedboundary)) + @test size(raw) == (2msing, 2msing) + ref = FFS.pest3_decompose(raw) + @test blocks.A ≈ ref.A && blocks.B ≈ ref.B && blocks.Gamma ≈ ref.Gamma && blocks.Delta ≈ ref.Delta end - # The rpec coil block is normalized to the Riccati orientation at pack time. - @test size(dp.coil) == (2msing, ffs.numpert_total) - @test !iszero(dp.coil) + # The EF block is normalized to the Riccati orientation at pack time. + @test size(dp.ef) == (2msing, ffs.numpert_total) + @test !iszero(dp.ef) + @test dp.edge_response == ffs.galerkin.edge_response[:, 1:(2msing)] _check_main_display(ret) foreach(_check_display, (ffs, ffs.galerkin, ffs.solution)) - # The ideal-flag match skips the inner layer, so its basis is ideal-closed. + # Without a match the EF columns are the ideal closure, with a Hermitian δW. if ideal @test ffs.closure === :ideal @test iszero(ffs.bpen) + @test norm(ffs.wp - ffs.wp') <= 1e-8 * norm(ffs.wp) else @test ffs.closure === :matched @test ffs.bpen == ffs.match.bpen @test !iszero(ffs.bpen) - # The ideal δW does not describe the matched plasma. - @test ffs.wp === nothing && ffs.free_boundary === nothing end # The closed profiles land in the shared Solutions layout: same names and @@ -251,16 +251,19 @@ include("display_check.jl") @test read(f["$gal/xi_psi"]) == sol.u_store[:, :, 1, :] @test read(f["$gal/dxi_psidpsi"]) == sol.du_store @test read(f["$gal/xi_s"]) == sol.xi_s_store + @test read(f["ForceFreeStates/FreeBoundaryStability/W_plasma"]) == ffs.wp # Matching diagnostics sit with the other per-surface results, and the # raw outer basis is debug-gated off by default. - @test haskey(f, "SingularSurfaces/Match/cout") + @test haskey(f, "SingularSurfaces/Match/cout") != ideal @test !haskey(f, "$gal/Match") - # Match/ is self-describing: rows follow SingularSurfaces/rational_psi, and it names the model that ran. - mg = f["SingularSurfaces/Match"] - @test size(read(mg["bpen"]), 1) == length(read(f["SingularSurfaces/rational_psi"])) - @test read(mg["reconnected_flux"]) == ffs.match.reconnected_flux - @test attrs(mg)["inner_model"] == GeneralizedPerturbedEquilibrium.InnerLayer.model_name(ffs.match.model) - @test isempty(_collect_metadata_violations(mg)) + if !ideal + # Match/ is self-describing: rows follow SingularSurfaces/rational_psi, and it names the model that ran. + mg = f["SingularSurfaces/Match"] + @test size(read(mg["bpen"]), 1) == length(read(f["SingularSurfaces/rational_psi"])) + @test read(mg["reconnected_flux"]) == ffs.match.reconnected_flux + @test attrs(mg)["inner_model"] == GeneralizedPerturbedEquilibrium.InnerLayer.model_name(ffs.match.model) + @test isempty(_collect_metadata_violations(mg)) + end @test !haskey(f, "$gal/Basis") @test !haskey(f, "$gal/Solution") end diff --git a/test/runtests_riccati.jl b/test/runtests_riccati.jl index b3904b970..d044eb3dc 100644 --- a/test/runtests_riccati.jl +++ b/test/runtests_riccati.jl @@ -102,13 +102,13 @@ end # (which overwrites intr_ric.sing[s].delta_prime) delta_prime_inline = [copy(intr_ric.sing[s].delta_prime) for s in 1:intr_ric.msing] - vac_ric = FFS.free_run(odet_ric, ctrl, equil, mats, intr_ric) + vac_ric = FFS.free_run((odet_ric.u[:, :, 2] / odet_ric.u[:, :, 1]) ./ equil.psio^2, ctrl, equil, mats, intr_ric) et_ric = real(vac_ric.et[1]) # Forward integration (needed only for energy comparison). _, intr_fwd, _, _, _ = prepared_solovev() odet_fwd, _, _, _ = FFS.forward_eulerlagrange_integration(ctrl, equil, mats, intr_fwd) - vac_fwd = FFS.free_run(odet_fwd, ctrl, equil, mats, intr_fwd) + vac_fwd = FFS.free_run((odet_fwd.u[:, :, 2] / odet_fwd.u[:, :, 1]) ./ equil.psio^2, ctrl, equil, mats, intr_fwd) et_fwd = real(vac_fwd.et[1]) # ───────────────────────────────────────────────────────────────────────── diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index 434047642..467e145ba 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -139,7 +139,7 @@ end # Δ′ is unified across formalisms, so a Galerkin run feeds SLAYER exactly as a Riccati one - # does: `result.delta_prime.matrix` is populated and already sized to the surface list, which + # does: `result.delta_prime.freeboundary` is populated and already sized to the surface list, which # is the predicate the tearing solve uses to accept it over the per-surface diagonal stub. @testset "Galerkin-fed SLAYER: gal Δ' drives the coupled solve" begin mktempdir() do dir @@ -158,7 +158,7 @@ ffs = GeneralizedPerturbedEquilibrium.main([dir]).ffs @test ffs.integrator === :galerkin - dpm = ffs.delta_prime.matrix + dpm = ffs.delta_prime.freeboundary.Delta @test size(dpm) == (length(ffs.surfaces), length(ffs.surfaces)) @test size(dpm, 1) == 2 diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 2e5901a52..8dadabaf3 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -76,18 +76,18 @@ include("display_check.jl") @test api.solution === nothing @test api.delta_prime !== nothing msing = length(api.surfaces) - @test size(api.delta_prime.matrix) == (msing, msing) - @test api.delta_prime.matrix ≈ FFS.pest3_decompose(api.delta_prime.raw).Δ + @test size(api.delta_prime.freeboundary.Delta) == (msing, msing) + @test api.delta_prime.freeboundary.Delta ≈ FFS.pest3_decompose(api.delta_prime.raw_freeboundary).Delta # Chunking is a decomposition of the same problem, so Δ′ tracks the TOML run. - @test size(reference.delta_prime.matrix) == size(api.delta_prime.matrix) - @test api.delta_prime.matrix ≈ reference.delta_prime.matrix rtol = 1e-6 + @test size(reference.delta_prime.freeboundary.Delta) == size(api.delta_prime.freeboundary.Delta) + @test api.delta_prime.freeboundary.Delta ≈ reference.delta_prime.freeboundary.Delta rtol = 1e-6 # With the deck's own (auto) chunking the two runs are the same computation. same = solve(prob, riccati) - @test same.delta_prime.matrix ≈ reference.delta_prime.matrix rtol = 1e-12 - @test same.delta_prime.raw ≈ reference.delta_prime.raw rtol = 1e-12 - @test same.delta_prime.coil ≈ reference.delta_prime.coil rtol = 1e-12 + @test same.delta_prime.freeboundary.Delta ≈ reference.delta_prime.freeboundary.Delta rtol = 1e-12 + @test same.delta_prime.raw_freeboundary ≈ reference.delta_prime.raw_freeboundary rtol = 1e-12 + @test same.delta_prime.ef ≈ reference.delta_prime.ef rtol = 1e-12 end end @@ -193,13 +193,13 @@ include("display_check.jl") @test api.integrator === :galerkin @test api.control.algorithm.nx == 32 @test api.galerkin !== nothing - @test api.free_boundary === nothing # galerkin computes no free-boundary energies + # By default the EF columns and the vacuum give the ideal solution and its free-boundary energies. + @test api.control.algorithm.error_field_response && api.free_boundary !== nothing dp = api.delta_prime @test dp !== nothing msing = api.galerkin.msing - @test size(dp.matrix) == (msing, msing) - @test dp.matrix ≈ FFS.pest3_decompose(dp.raw).Δ - @test dp.A !== nothing # the parity blocks the galerkin solve persists + @test size(dp.fixedboundary.Delta) == (msing, msing) + @test dp.fixedboundary.Delta ≈ FFS.pest3_decompose(dp.raw_fixedboundary).Delta reference = GPEC.main([ @@ -207,10 +207,13 @@ include("display_check.jl") extra_sections=Dict{String,Any}("Galerkin" => Dict{String,Any}("nx" => 32))) ]).ffs @test reference.galerkin.msing == msing - @test dp.matrix ≈ reference.delta_prime.matrix rtol = 1e-12 - @test dp.raw ≈ reference.delta_prime.raw rtol = 1e-12 - # With no free-boundary δW, no plasma response can be built from it. - @test_throws ErrorException plasma_response(api) + @test dp.raw_fixedboundary ≈ reference.delta_prime.raw_fixedboundary rtol = 1e-12 + @test dp.raw_freeboundary ≈ reference.delta_prime.raw_freeboundary rtol = 1e-12 + # Without the EF columns there is no free-boundary δW to build a plasma response from. + noef = solve(equil, Galerkin(; nx=32, error_field_response=false); nn=1, wall=wall, ffs_kwargs...) + @test noef.free_boundary === nothing && noef.solution === nothing && isempty(noef.delta_prime.ef) + @test noef.delta_prime.raw_freeboundary == dp.raw_freeboundary + @test_throws ErrorException plasma_response(noef) end end @@ -385,7 +388,7 @@ include("display_check.jl") @testset "solves leave their inputs untouched" begin snap = deepcopy(equil) mktempdir() do dir - gal = solve(equil, Galerkin(; nx=32, rpec_flag=true, cut_solution=true); nn=1, ffs_kwargs...) + gal = solve(equil, Galerkin(; nx=32, error_field_response=true, cut_solution=true); nn=1, ffs_kwargs...) @test first_diff(equil, snap) === nothing # Matching reuses the outer solve: repeated match solves must not alter it. gal_snap = deepcopy(gal) @@ -400,7 +403,7 @@ include("display_check.jl") end @testset "MatchProblem defaults its layer inputs from the kinetic data" begin - gal = solve(equil, Galerkin(; nx=32, rpec_flag=true, cut_solution=true); nn=1, ffs_kwargs...) + gal = solve(equil, Galerkin(; nx=32, error_field_response=true, cut_solution=true); nn=1, ffs_kwargs...) # Without kinetic data the error names exactly the inputs still to give. @test_throws "MatchProblem: give eta, rho and rotation, or load kinetic profiles" MatchProblem(gal) @test_throws "MatchProblem: give rho, or load kinetic profiles" MatchProblem(gal; eta=1e-7, rotation=1.0) @@ -435,7 +438,7 @@ include("display_check.jl") @testset "matched inner profiles: opt-in, and parity columns agree across inner models" begin mktempdir() do dir - gal = solve(equil, Galerkin(; nx=32, rpec_flag=true, cut_solution=true); nn=1, ffs_kwargs...) + gal = solve(equil, Galerkin(; nx=32, error_field_response=true, cut_solution=true); nn=1, ffs_kwargs...) prob = MatchProblem(gal; eta=1e-7, rho=1e-7, rotation=1.0) center = solve(prob, GGJModel()).match ray = solve(prob, GGJModel(; real_profiles=true)).match @@ -456,19 +459,23 @@ include("display_check.jl") end end - @testset "a resistive match keeps its coefficients and drops the ideal δW" begin + @testset "a Riccati match keeps its coefficients and publishes its matched edge response" begin mktempdir() do dir ric = solve(equil, riccati; nn=1, ffs_kwargs...) - @test ric.free_boundary !== nothing && !isempty(ric.delta_prime.coil) + @test ric.free_boundary !== nothing && !isempty(ric.delta_prime.ef) matched = solve(MatchProblem(ric; eta=1e-7, rho=1e-7, rotation=1.0), GGJModel()) @test matched.closure === :matched @test matched.match !== nothing && matched.bpen == matched.match.bpen - @test matched.wp === nothing && matched.free_boundary === nothing @test matched.solution === nothing # no Galerkin basis to build a matched ξ from - # It carries no ideal δW, so no plasma response can be built from it. - @test_throws ErrorException plasma_response(matched) + # The edge response is W_P plus the big-solution edge momenta weighted by cout; the layers make it non-Hermitian. + ψ₀² = ric.equil.psio^2 + @test matched.wp ≈ ric.wp + ric.delta_prime.edge_response * matched.match.cout / ψ₀² rtol = 1e-12 + @test norm(matched.wp - matched.wp') > 1e-8 * norm(matched.wp) + @test matched.free_boundary.wt0 ≈ matched.wp + ric.free_boundary.wv rtol = 1e-12 + # The singular coupling of a matched solution is not defined yet. + @test_throws "matched" plasma_response(matched) # A matched result is a new solution, not an outer solve to match again. - @test_throws ErrorException MatchProblem(matched; ideal=true) + @test_throws ErrorException MatchProblem(matched; eta=1e-7, rho=1e-7, rotation=1.0) end end @@ -476,24 +483,62 @@ include("display_check.jl") mktempdir() do dir # A Forward result carries no Δ′ payload, so the problem is unconstructible. fwd = solve(equil, forward; nn=1, ffs_kwargs...) - @test_throws ErrorException MatchProblem(fwd; ideal=true) + @test_throws ErrorException MatchProblem(fwd; eta=1e-7, rho=1e-7, rotation=1.0) # A slab model can never close a matched solution. - gal = solve(equil, Galerkin(; nx=32, rpec_flag=true); nn=1, ffs_kwargs...) - prob = MatchProblem(gal; ideal=true) + gal = solve(equil, Galerkin(; nx=32, error_field_response=true); nn=1, ffs_kwargs...) + prob = MatchProblem(gal; eta=1e-7, rho=1e-7, rotation=1.0) @test_throws ErrorException solve(prob, SLAYERModel()) # A resistive match needs every layer input. @test_throws "give rho and rotation" MatchProblem(gal; eta=1e-7) @test_throws ErrorException MatchProblem(gal; eta=fill(1e-7, length(gal.surfaces) + 1), rho=1e-7, rotation=1.0) # The tearing model is the solve argument, never a TearingProblem keyword. @test_throws MethodError TearingProblem(gal; inner_model=:ggj_ray) - # The ideal reference match keeps the ideal closure and replaces the solution with - # the bare coil columns in the identity-at-edge basis. + # The Galerkin solve is itself ideal-closed: its EF columns, identity at the edge, with their δW. + @test gal.closure === :ideal + @test gal.solution !== nothing && gal.solution.basis === :gal_native + N = gal.numpert_total + @test gal.solution.u_store[:, :, 1, end] ≈ Matrix{ComplexF64}(I, N, N) atol = 1e-10 + @test size(gal.wp) == (N, N) && norm(gal.wp - gal.wp') <= 1e-8 * norm(gal.wp) + # A matched Galerkin result keeps the identity-at-edge basis and publishes its complex edge response. matched = solve(prob, GGJModel()) - @test matched.closure === :ideal - @test matched.match !== nothing - # The ideal reference keeps the ideal δW (only a :matched closure drops it). - @test matched.wp === gal.wp && matched.free_boundary === gal.free_boundary - @test matched.solution !== nothing && matched.solution.basis === :gal_native + @test matched.closure === :matched + @test matched.solution.u_store[:, :, 1, end] ≈ Matrix{ComplexF64}(I, N, N) atol = 1e-10 + ψ₀² = gal.equil.psio^2 + F = gal.galerkin.edge_response + @test matched.wp ≈ gal.wp + F[:, 1:(2length(gal.surfaces))] * matched.match.cout / ψ₀² rtol = 1e-12 + # Its free-boundary energies use the matched response with the solve's own vacuum. + @test matched.free_boundary.wt0 ≈ matched.wp + gal.free_boundary.wv rtol = 1e-12 + # The singular coupling of a matched solution is not defined yet. + @test_throws "matched" plasma_response(matched) + # The EF columns do not change either Δ′: both equal those of a solve without them. + plain = solve(equil, Galerkin(; nx=32, error_field_response=false); nn=1, ffs_kwargs...) + @test gal.delta_prime.raw_freeboundary == plain.delta_prime.raw_freeboundary + @test gal.delta_prime.raw_fixedboundary == plain.delta_prime.raw_fixedboundary + @test gal.delta_prime.raw_fixedboundary != gal.delta_prime.raw_freeboundary + # Skipping the free-boundary solve leaves only the fixed-boundary Δ′: matching is unchanged, and tearing + # runs on the fixed boundary only when asked to. + nofree = solve(equil, Galerkin(; nx=32, error_field_response=true, free_boundary_delta_prime=false); nn=1, ffs_kwargs...) + @test isempty(nofree.delta_prime.raw_freeboundary) && nofree.delta_prime.freeboundary === nothing + @test nofree.delta_prime.raw_fixedboundary == gal.delta_prime.raw_fixedboundary && nofree.delta_prime.ef == gal.delta_prime.ef + @test solve(MatchProblem(nofree; eta=1e-7, rho=1e-7, rotation=1.0), GGJModel()).match.cout == matched.match.cout + @test_throws "free-boundary Δ′" solve(TearingProblem(nofree), SLAYERModel()) + end + end + + @testset "a Galerkin solve drives the plasma response like a forward solve" begin + mktempdir() do dir + fwd = solve(equil, forward; nn=1, ffs_kwargs...) + gal = solve(equil, Galerkin(; nx=32, error_field_response=true); nn=1, ffs_kwargs...) + @test gal.psilim == fwd.psilim + @test gal.wp ≈ fwd.wp rtol = 1e-5 + @test gal.free_boundary.wt0 ≈ gal.wp + gal.free_boundary.wv rtol = 1e-12 + @test gal.free_boundary.et[1] ≈ fwd.free_boundary.et[1] rtol = 1e-5 + # The edge permeability P = Λ·L⁻¹ does not depend on the solution basis. + @test plasma_response(gal).permeability_flux ≈ plasma_response(fwd).permeability_flux rtol = 5e-4 + # Without the vacuum there is no free-boundary δW to build Λ from. + novac = solve(equil, Galerkin(; nx=32, error_field_response=true); nn=1, merge(ffs_kwargs, Dict{Symbol,Any}(:vac_flag => false))...) + @test novac.wp !== nothing && novac.free_boundary === nothing + @test_throws "vac_flag=true" plasma_response(novac) end end @@ -546,11 +591,11 @@ include("display_check.jl") "Match" => Dict{String,Any}("eta" => fill(layer.eta, n), "rho" => fill(layer.rho, n), "rotation" => fill(layer.rotation, n))) mktempdir() do dir - gal = solve(equil, Galerkin(; nx=32, rpec_flag=true); nn=1, wall=wall, ffs_kwargs...) + gal = solve(equil, Galerkin(; nx=32, error_field_response=true); nn=1, wall=wall, ffs_kwargs...) # The galerkin inner model returns real-ψ profiles, so the inner profiles are compared too. api = solve(MatchProblem(gal; layer...), GGJModel(; solver=:galerkin)) n = length(gal.surfaces) - sections = merge(Dict{String,Any}("Galerkin" => Dict{String,Any}("nx" => 32, "rpec_flag" => true)), match_sections(n, "galerkin")) + sections = merge(Dict{String,Any}("Galerkin" => Dict{String,Any}("nx" => 32, "error_field_response" => true)), match_sections(n, "galerkin")) reference = GPEC.main([_stage_deck(dir; ffs_overrides=Dict{String,Any}("integrator" => "Galerkin"), extra_sections=sections)]).ffs @test reference.closure === api.closure === :matched @@ -559,7 +604,7 @@ include("display_check.jl") @test api.match.xi ≈ reference.match.xi rtol = 1e-12 # The composite inner profiles need the cut solution, which Galerkin keeps only on request. @test isempty(reference.match.inner_xi) && isempty(api.match.inner_xi) - cut = solve(equil, Galerkin(; nx=32, rpec_flag=true, cut_solution=true); nn=1, wall=wall, ffs_kwargs...) + cut = solve(equil, Galerkin(; nx=32, error_field_response=true, cut_solution=true); nn=1, wall=wall, ffs_kwargs...) @test length(solve(MatchProblem(cut; layer...), GGJModel(; solver=:galerkin)).match.inner_xi) == n end @@ -576,7 +621,7 @@ include("display_check.jl") # A forward solve has no Δ′ to match against, so the deck stops at parse time instead of after the solve. mktempdir() do dir - deck_dir = _stage_deck(dir; extra_sections=Dict{String,Any}("Match" => Dict{String,Any}("ideal" => true))) + deck_dir = _stage_deck(dir; extra_sections=Dict{String,Any}("Match" => Dict{String,Any}("eta" => 1e-7, "rho" => 1e-7, "rotation" => 1.0))) @test_throws "[Match] needs Δ′" GPEC.parse_deck(TOML.parsefile(joinpath(deck_dir, "gpec.toml")), deck_dir) end