diff --git a/.claude/agents/regression-guardian.md b/.claude/agents/regression-guardian.md index b74ceb16b..e98ecac11 100644 --- a/.claude/agents/regression-guardian.md +++ b/.claude/agents/regression-guardian.md @@ -36,8 +36,8 @@ Useful invocations: cases are added over time, so never assume the list. (At time of writing: `diiid_n1`, `diiid_n1_riccati`, `diiid_slayer_n1`, `efit_fixedbdy_separatrix`, `gal_resistive_diiid`, `gal_resistive_pe`, `ggj_ray_q500i`, `ggj_reference`, `solovev_n1`, `solovev_multi_n`, - `solovev_kinetic_calculated`, `solovev_kinetic_multiion`, `solovev_kinetic_ntv`, - `solovev_kinetic_nuzero`.) + `solovev_kinetic_nuzero`, `tj_n1`, `tj_kinetic_calculated`, `tj_kinetic_multiion`, + `tj_kinetic_ntv`.) - `--cases --refs ` — compare two refs (branches/commits), e.g. `--refs develop,my-feature-branch`. `local` names the uncommitted working tree; prefer a branch ref (see "Which refs to compare" below). @@ -49,7 +49,7 @@ Useful invocations: Pick the case(s) by the module touched: - Equilibrium / ideal stability / perturbed equilibrium → `diiid_n1` (broadest, ~31 quantities). - Analytical equilibrium + ideal stability → `solovev_n1`, `solovev_multi_n`. -- KineticForces (NTV) changes → `solovev_kinetic_calculated`. +- KineticForces (NTV) changes → `tj_kinetic_calculated`. - InnerLayer (GGJ resistive layer) changes → `ggj_reference`. When in doubt, run `diiid_n1` plus whichever case targets the module you changed. diff --git a/examples/Solovev_ideal_example/analyze_example.jl b/examples/Solovev_ideal_example/analyze_example.jl index 4c5a8eb7c..9f08e1469 100644 --- a/examples/Solovev_ideal_example/analyze_example.jl +++ b/examples/Solovev_ideal_example/analyze_example.jl @@ -9,7 +9,6 @@ h5path = joinpath(@__DIR__, "gpec.h5") # Summary plots p_eq = Analysis.Equilibrium.plot_equilibrium_summary(h5path) p_ffs = Analysis.ForceFreeStates.plot_ffs_summary(h5path) -p_pe = Analysis.PerturbedEquilibrium.plot_perturbed_equilibrium_summary(h5path) display(p_eq); Plots.savefig(p_eq, joinpath(@__DIR__, "equilibrium_summary.png")) display(p_ffs); Plots.savefig(p_ffs, joinpath(@__DIR__, "ffs_summary.png")) diff --git a/examples/Solovev_ideal_example/gpec.toml b/examples/Solovev_ideal_example/gpec.toml index 330f1bba6..7cc873c46 100644 --- a/examples/Solovev_ideal_example/gpec.toml +++ b/examples/Solovev_ideal_example/gpec.toml @@ -1,4 +1,4 @@ -# Solovev analytical equilibrium — n=1 ideal stability + perturbed equilibrium. +# Solovev analytical equilibrium — n=1 ideal stability. # Near-marginal n=1 external kink stabilized by a close conformal wall (see [Wall] note). # Integration-pipeline smoke-test and et[1] regression fixture; the equilibrium is # generated analytically from the embedded [Solovev] section. @@ -27,14 +27,6 @@ shape = "conformal" # Wall shape (conformal, elliptical, dee standoff_fraction = 0.2415 # Wall standoff in plasma minor radii equal_arc_wall = true # Equal arc length distribution of nodes on wall -[ForcingTerms] -forcing_data_file = "forcing.dat" # Manual mode table (n, m, amplitude_real, amplitude_imag) -forcing_data_format = "ascii" # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires) - -[PerturbedEquilibrium] -compute_singular_coupling = true # Compute singular layer coupling metrics -regularization_width = 0.0 # Regularization width for singular surfaces (0 = disabled) - # Note: SLAYER tearing analysis is not run on the Solovev analytic equilibrium # — its Δ' / inner-layer dispersion does not yield meaningful tearing roots. # The SLAYER example and regression case use the DIII-D-like geqdsk instead diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml deleted file mode 100644 index 6c363ffb5..000000000 --- a/examples/Solovev_kinetic_NTV_example/gpec.toml +++ /dev/null @@ -1,94 +0,0 @@ -# Solovev analytical equilibrium — n=1 ideal stability + perturbed equilibrium + NTV torque. -# Extends the Solovev ideal example with a [KineticForces] section, so after the plasma -# response is computed the neoclassical toroidal viscosity torque is integrated over ψ as a -# post-processing diagnostic (adaptive quadrature paneled at the rational surfaces). -# The kinetic profiles (kinetic.dat) are constructed to be consistent with this equilibrium's -# pressure; the equilibrium is generated analytically from the embedded [Solovev] section. - -[Run] -verbose = true # Log each stage's progress - -[Equilibrium] -kinetic_file = "kinetic.dat" # Kinetic-profile file (HDF5 GPEC kinetic schema or 6-column ASCII), read once into the equilibrium -jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, other) -grid_type = "core_edge_packed" # Radial grid packing type -psilow = 1e-4 # Lower limit of normalized poloidal flux -psihigh = 0.9995 # Upper limit of normalized poloidal flux -mpsi = 128 # Number of radial grid points (0 = auto-compute from psi_accuracy) -mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) -flux_surface_rtol = 1e-7 # Relative tolerance of the ODE solves that trace each flux surface - -[Wall] -# Close conformal wall is required to stabilize this Solovev fixture's n=1 external kink: -# with nowall, et[1] = -6.8 (strongly unstable); with this wall, et[1] = +0.24 (barely stable). -# The plasma is near marginal stability, so the BVP Δ' matrix values are pathological -# (dpm magnitudes ~ 10¹¹, |Im/Re| ≫ 1). This fixture's role is integration-pipeline -# smoke testing + et[1] regression, NOT BVP Δ' regression — DIIID-like is the canonical -# Δ'-matrix fixture (stable et[1] = +1.6, clean BVP Δ'). -shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file) -standoff_fraction = 0.2415 # Wall standoff in plasma minor radii -equal_arc_wall = true # Equal arc length distribution of nodes on wall - -[ForcingTerms] -forcing_data_file = "forcing.dat" # Manual mode table (n, m, amplitude_real, amplitude_imag) -forcing_data_format = "ascii" # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires) - -[PerturbedEquilibrium] -compute_singular_coupling = true # Compute singular layer coupling metrics -regularization_width = 0.05 # Regularization width for singular surfaces (0 = disabled) - -# Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the LocalStabilityProblem settings, numerics in [Ballooning] -[LocalStability] -theta_k = 0.0 # Ballooning angle θ_k, used by both the profile scan and the α-boundary scan -alpha_boundary_scan = true # Run the first ballooning α-boundary scan (a root find per ψ surface) - -# Numerical settings of the ballooning scans -[Ballooning] -ode_tol = 1e-8 # Relative and absolute tolerance of the ballooning ODE integration behind Δ' -n_scan = 24 # Resolution of the α scan for the first boundary -max_alpha_scale = 8.0 # Upper end of the α search, as a multiple of the experimental α -alpha_tol = 1e-3 # Bisection tolerance of the critical α, in the same multiples -psi_window = [0.1, 0.99] # ψ_N range of the α scan; surfaces outside it are NaN - -[ForceFreeStates] -vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes - -dW_edge_scan_start = 0.99 # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable -qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) -qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) - -nn = 1 # Toroidal mode number, or consecutive numbers as a list -delta_mlow = 8 # Expands lower bound of Fourier harmonics -delta_mhigh = 8 # Expands upper bound of Fourier harmonics -mtheta_vacuum = 960 # Number of points used in splines over poloidal angle at the plasma-vacuum interface - -singfac_min = 1e-6 # Fractional distance from rational q at which ideal jump enforced - -# Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) -integrator = "Forward" # Integration formalism, whose numerical settings are its own section: "Forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "Riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "Galerkin" (outer-region singular Galerkin Δ′) -truncate_at_rational_offset = false # FALSE for limited/analytical equilibria — rationals sparse, rational_offset_fraction would chop too much edge -rational_offset_fraction = 0.2 # Truncate integration at (last_rational_q + rational_offset_fraction)/n (used when truncate_at_rational_offset = true) - -# Numerical settings of the forward integrator -[Forward] -eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations -renorm_threshold = 1e3 # Column-norm threshold that triggers solution renormalization -save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. - -# Solovev analytic equilibrium parameters; this section selects the equilibrium. See SolovevConfig in src/Equilibrium. -[Solovev] -mr = 128 # Number of radial grid zones -mz = 128 # Number of axial grid zones -ma = 128 # Number of flux grid zones -elongation = 1.6 # Elongation -a = 0.33 # Minor radius -r0 = 1.0 # Major radius -q0 = 1.9 # Safety factor at the magnetic axis (O-point) -p0fac = 1 # Scale on-axis pressure (β changes; Φ, q fixed) -b0fac = 1 # Scale toroidal field at constant β (Bt changes; shape, β fixed) -f0fac = 1 # Scale toroidal field at constant pressure (β, q change; Φ, p, Bp fixed) - -[KineticForces] -fcgl_flag = true # Fluid Chew–Goldberger–Low energy (l = 0 fluid limit) alongside the default fgar method -axis_validity_suppression = true # Suppress kinetic terms inside the near-axis orbit-width validity limit diff --git a/examples/Solovev_kinetic_NTV_example/kinetic.dat b/examples/Solovev_kinetic_NTV_example/kinetic.dat deleted file mode 100644 index 96fd9fa3e..000000000 --- a/examples/Solovev_kinetic_NTV_example/kinetic.dat +++ /dev/null @@ -1,23 +0,0 @@ -# Kinetic profiles consistent with this Solovev equilibrium's pressure. -# Constraint: P(psi_n) = n_e*T_e + n_i*T_i must equal the equilibrium pressure -# P(psi_n) = P0*(1 - psi_n), P0 = 4.273e4 Pa (linear, zero at the edge). -# Construction (quasineutral, n_e=n_i, T_e=T_i): -# T(psi_n) = 100 + 2600*(1 - psi_n) eV (core 2700 eV, edge 100 eV) -# n(psi_n) = P0*(1 - psi_n) / (2*e*T) (core ~4.9e19; product n0*T0 = P0/2e fixed) -# Density carries the edge falloff to a small floor at psi_n=1 so T stays finite -# (avoids T->0 in the collisionality / drift-frequency formulas); the residual edge -# pressure from the floor is ~1e-6 of P0. omega_E is a free (non-pressure) profile. -psi_n n_i_m3 n_e_m3 T_i_eV T_e_eV omega_E_rads -0.00 4.939e19 4.939e19 2700.0 2700.0 1.0e4 -0.10 4.919e19 4.919e19 2440.0 2440.0 9.0e3 -0.20 4.894e19 4.894e19 2180.0 2180.0 8.0e3 -0.30 4.862e19 4.862e19 1920.0 1920.0 7.0e3 -0.40 4.820e19 4.820e19 1660.0 1660.0 6.0e3 -0.50 4.763e19 4.763e19 1400.0 1400.0 5.0e3 -0.60 4.679e19 4.679e19 1140.0 1140.0 4.0e3 -0.70 4.546e19 4.546e19 880.0 880.0 3.0e3 -0.80 4.302e19 4.302e19 620.0 620.0 2.0e3 -0.90 3.705e19 3.705e19 360.0 360.0 1.0e3 -0.95 2.899e19 2.899e19 230.0 230.0 5.0e2 -0.99 1.058e19 1.058e19 126.0 126.0 1.0e2 -1.00 2.000e18 2.000e18 100.0 100.0 0.0 diff --git a/examples/Solovev_kinetic_calculated_example/gpec.toml b/examples/Solovev_kinetic_calculated_example/gpec.toml deleted file mode 100644 index 416c8ee9a..000000000 --- a/examples/Solovev_kinetic_calculated_example/gpec.toml +++ /dev/null @@ -1,94 +0,0 @@ -# Solovev analytical equilibrium — n=1 ideal stability with calculated kinetic matrices. -# Exercises the real KineticForces (NTV) physics path at full strength, with the main ion -# declared as a 50/50 D-T mix so the self-consistent kinetic matrices are accumulated over -# the multi-ion species loop (pure plasma, n_i = n_e, so no impurity species arises). The -# equilibrium is generated analytically from the embedded [Solovev] section, with n·T tied -# to the Solovev P₀≈4.27e4 Pa via the local kinetic.dat. -# Also the only deck that runs [PerturbedEquilibrium] on a kinetic solve, so it builds the -# ψ-resolved torque response matrices (T_xe on the applied field, T_source for the forcing file). - -[Run] -verbose = true # Log each stage's progress - -[Equilibrium] -kinetic_file = "kinetic.dat" # Kinetic-profile file (HDF5 GPEC kinetic schema or 6-column ASCII), read once into the equilibrium -jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) -grid_type = "core_edge_packed" # Radial grid packing type -psilow = 1e-4 # Lower limit of normalized flux coordinate -psihigh = 0.9995 # Upper limit of normalized flux coordinate -mpsi = 16 # Number of radial grid points -mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) -flux_surface_rtol = 1e-7 # Relative tolerance of the ODE solves that trace each flux surface - -[[Equilibrium.ion_species]] -z = 1 # Species charge (fundamental units) -m = 2 # Species mass (proton masses) — Deuterium -fraction = 0.5 # This species' share of the total main-ion density n_i - -[[Equilibrium.ion_species]] -z = 1 # Species charge (fundamental units) -m = 3 # Species mass (proton masses) — Tritium -fraction = 0.5 # This species' share of the total main-ion density n_i - -[Wall] -shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file) -standoff_fraction = 0.2415 # Wall standoff in plasma minor radii -equal_arc_wall = true # Equal arc length distribution of nodes on wall - -[ForcingTerms] -forcing_data_file = "forcing.dat" # Manual mode table (n, m, amplitude_real, amplitude_imag) -forcing_data_format = "ascii" # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires) - -[PerturbedEquilibrium] -compute_singular_coupling = true # Compute singular layer coupling metrics -compute_torque_response = true # Build the ψ-resolved torque response matrices T_xe and T_source - -# Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the LocalStabilityProblem settings, numerics in [Ballooning] -[LocalStability] -theta_k = 0.0 # Ballooning angle θ_k, used by both the profile scan and the α-boundary scan -alpha_boundary_scan = true # Run the first ballooning α-boundary scan (a root find per ψ surface) - -# Numerical settings of the ballooning scans -[Ballooning] -ode_tol = 1e-8 # Relative and absolute tolerance of the ballooning ODE integration behind Δ' -n_scan = 24 # Resolution of the α scan for the first boundary -max_alpha_scale = 8.0 # Upper end of the α search, as a multiple of the experimental α -alpha_tol = 1e-3 # Bisection tolerance of the critical α, in the same multiples -psi_window = [0.1, 0.99] # ψ_N range of the α scan; surfaces outside it are NaN - -[ForceFreeStates] -vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes - -dW_edge_scan_start = 0.99 # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable -qlow = 1.02 # Integration initiated at q determined by min(q0, qlow)... -qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh)... - -nn = 1 # Toroidal mode number, or consecutive numbers as a list -delta_mlow = 8 # Expands lower bound of Fourier harmonics -delta_mhigh = 8 # Expands upper bound of Fourier harmonics -mtheta_vacuum = 64 # Number of points used in splines over poloidal angle at plasma-vacuum interface. - -kinetic = "KineticForces" # Kinetic source: the drift-kinetic matrices of the [KineticForces] model (omit for an ideal run) -integrator = "Forward" # Integration formalism, whose numerical settings are its own section: "Forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "Riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "Galerkin" (outer-region singular Galerkin Δ′) -singfac_min = 1e-6 # Fractional distance from rational q at which ideal jump enforced - -# Numerical settings of the forward integrator -[Forward] -eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations -renorm_threshold = 1e3 # Column-norm threshold that triggers solution renormalization - -[KineticForces] - -# Solovev analytic equilibrium parameters; this section selects the equilibrium. See SolovevConfig in src/Equilibrium. -[Solovev] -mr = 128 # Number of radial grid zones -mz = 128 # Number of axial grid zones -ma = 128 # Number of flux grid zones -elongation = 1.6 # Elongation -a = 0.33 # Minor radius -r0 = 1.0 # Major radius -q0 = 1.9 # Safety factor at the magnetic axis (O-point) -p0fac = 1 # Scale on-axis pressure (β changes; Φ, q fixed) -b0fac = 1 # Scale toroidal field at constant β (Bt changes; shape, β fixed) -f0fac = 1 # Scale toroidal field at constant pressure (β, q change; Φ, p, Bp fixed) diff --git a/examples/Solovev_kinetic_calculated_example/kinetic.dat b/examples/Solovev_kinetic_calculated_example/kinetic.dat deleted file mode 100644 index 96fd9fa3e..000000000 --- a/examples/Solovev_kinetic_calculated_example/kinetic.dat +++ /dev/null @@ -1,23 +0,0 @@ -# Kinetic profiles consistent with this Solovev equilibrium's pressure. -# Constraint: P(psi_n) = n_e*T_e + n_i*T_i must equal the equilibrium pressure -# P(psi_n) = P0*(1 - psi_n), P0 = 4.273e4 Pa (linear, zero at the edge). -# Construction (quasineutral, n_e=n_i, T_e=T_i): -# T(psi_n) = 100 + 2600*(1 - psi_n) eV (core 2700 eV, edge 100 eV) -# n(psi_n) = P0*(1 - psi_n) / (2*e*T) (core ~4.9e19; product n0*T0 = P0/2e fixed) -# Density carries the edge falloff to a small floor at psi_n=1 so T stays finite -# (avoids T->0 in the collisionality / drift-frequency formulas); the residual edge -# pressure from the floor is ~1e-6 of P0. omega_E is a free (non-pressure) profile. -psi_n n_i_m3 n_e_m3 T_i_eV T_e_eV omega_E_rads -0.00 4.939e19 4.939e19 2700.0 2700.0 1.0e4 -0.10 4.919e19 4.919e19 2440.0 2440.0 9.0e3 -0.20 4.894e19 4.894e19 2180.0 2180.0 8.0e3 -0.30 4.862e19 4.862e19 1920.0 1920.0 7.0e3 -0.40 4.820e19 4.820e19 1660.0 1660.0 6.0e3 -0.50 4.763e19 4.763e19 1400.0 1400.0 5.0e3 -0.60 4.679e19 4.679e19 1140.0 1140.0 4.0e3 -0.70 4.546e19 4.546e19 880.0 880.0 3.0e3 -0.80 4.302e19 4.302e19 620.0 620.0 2.0e3 -0.90 3.705e19 3.705e19 360.0 360.0 1.0e3 -0.95 2.899e19 2.899e19 230.0 230.0 5.0e2 -0.99 1.058e19 1.058e19 126.0 126.0 1.0e2 -1.00 2.000e18 2.000e18 100.0 100.0 0.0 diff --git a/examples/Solovev_kinetic_multiion_example/gpec.toml b/examples/Solovev_kinetic_multiion_example/gpec.toml deleted file mode 100644 index fc96aa120..000000000 --- a/examples/Solovev_kinetic_multiion_example/gpec.toml +++ /dev/null @@ -1,107 +0,0 @@ -# Solovev analytical equilibrium — n=1 ideal stability + perturbed equilibrium + NTV torque. -# Extends the Solovev ideal example with a [KineticForces] section, so after the plasma -# response is computed the neoclassical toroidal viscosity torque is integrated over ψ as a -# post-processing diagnostic (adaptive quadrature paneled at the rational surfaces). -# The kinetic profiles (kinetic.dat) are constructed to be consistent with this equilibrium's -# pressure; the equilibrium is generated analytically from the embedded [Solovev] section. - -[Run] -verbose = true # Log each stage's progress - -[Equilibrium] -kinetic_file = "kinetic.dat" # Kinetic-profile file (HDF5 GPEC kinetic schema or 6-column ASCII), read once into the equilibrium -jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, other) -grid_type = "core_edge_packed" # Radial grid packing type -psilow = 1e-4 # Lower limit of normalized poloidal flux -psihigh = 0.9995 # Upper limit of normalized poloidal flux -mpsi = 128 # Number of radial grid points (0 = auto-compute from psi_accuracy) -mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) -flux_surface_rtol = 1e-7 # Relative tolerance of the ODE solves that trace each flux surface -impurity_charge = 6 # Impurity charge (fundamental units); closes quasineutrality -impurity_mass = 12 # Impurity mass (proton masses) - -[[Equilibrium.ion_species]] -z = 1 # Species charge (fundamental units) -m = 2 # Species mass (proton masses) — Deuterium -fraction = 0.5 # This species' share of the total main-ion density n_i - -[[Equilibrium.ion_species]] -z = 1 # Species charge (fundamental units) -m = 3 # Species mass (proton masses) — Tritium -fraction = 0.5 # This species' share of the total main-ion density n_i - -[Wall] -# Close conformal wall is required to stabilize this Solovev fixture's n=1 external kink: -# with nowall, et[1] = -6.8 (strongly unstable); with this wall, et[1] = +0.24 (barely stable). -# The plasma is near marginal stability, so the BVP Δ' matrix values are pathological -# (dpm magnitudes ~ 10¹¹, |Im/Re| ≫ 1). This fixture's role is integration-pipeline -# smoke testing + et[1] regression, NOT BVP Δ' regression — DIIID-like is the canonical -# Δ'-matrix fixture (stable et[1] = +1.6, clean BVP Δ'). -shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file) -standoff_fraction = 0.2415 # Wall standoff in plasma minor radii -equal_arc_wall = true # Equal arc length distribution of nodes on wall - -[ForcingTerms] -forcing_data_file = "forcing.dat" # Manual mode table (n, m, amplitude_real, amplitude_imag) -forcing_data_format = "ascii" # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires) - -[PerturbedEquilibrium] -compute_singular_coupling = true # Compute singular layer coupling metrics -regularization_width = 0.05 # Regularization width for singular surfaces (0 = disabled) - -# Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the LocalStabilityProblem settings, numerics in [Ballooning] -[LocalStability] -theta_k = 0.0 # Ballooning angle θ_k, used by both the profile scan and the α-boundary scan -alpha_boundary_scan = true # Run the first ballooning α-boundary scan (a root find per ψ surface) - -# Numerical settings of the ballooning scans -[Ballooning] -ode_tol = 1e-8 # Relative and absolute tolerance of the ballooning ODE integration behind Δ' -n_scan = 24 # Resolution of the α scan for the first boundary -max_alpha_scale = 8.0 # Upper end of the α search, as a multiple of the experimental α -alpha_tol = 1e-3 # Bisection tolerance of the critical α, in the same multiples -psi_window = [0.1, 0.99] # ψ_N range of the α scan; surfaces outside it are NaN - -[ForceFreeStates] -integrator = "Forward" # Integration formalism, whose numerical settings are its own section: "Forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "Riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "Galerkin" (outer-region singular Galerkin Δ′) -vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes - -dW_edge_scan_start = 0.99 # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable -qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) -qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) - -nn = 1 # Toroidal mode number, or consecutive numbers as a list -delta_mlow = 8 # Expands lower bound of Fourier harmonics -delta_mhigh = 8 # Expands upper bound of Fourier harmonics -mtheta_vacuum = 960 # Number of points used in splines over poloidal angle at the plasma-vacuum interface - -singfac_min = 1e-6 # Fractional distance from rational q at which ideal jump enforced - -# Integration-range truncation (see ForceFreeStatesControl docstring for details) -truncate_at_rational_offset = false # FALSE for limited/analytical equilibria — rationals sparse, rational_offset_fraction would chop too much edge -rational_offset_fraction = 0.2 # Truncate integration at (last_rational_q + rational_offset_fraction)/n (used when truncate_at_rational_offset = true) - -# Numerical settings of the forward integrator -[Forward] -eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations -renorm_threshold = 1e3 # Column-norm threshold that triggers solution renormalization -save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. - -# Solovev analytic equilibrium parameters; this section selects the equilibrium. See SolovevConfig in src/Equilibrium. -[Solovev] -mr = 128 # Number of radial grid zones -mz = 128 # Number of axial grid zones -ma = 128 # Number of flux grid zones -elongation = 1.6 # Elongation -a = 0.33 # Minor radius -r0 = 1.0 # Major radius -q0 = 1.9 # Safety factor at the magnetic axis (O-point) -p0fac = 1 # Scale on-axis pressure (β changes; Φ, q fixed) -b0fac = 1 # Scale toroidal field at constant β (Bt changes; shape, β fixed) -f0fac = 1 # Scale toroidal field at constant pressure (β, q change; Φ, p, Bp fixed) - -[KineticForces] -# Multi-ion NTV: the analytic Solovev main ion is modelled as a 50/50 D-T mix (n_i = total), -# so the summed NTV runs over D + T + the neutrality-closing impurity + electrons. -electron = true # Add electron NTV in addition to the ion species diff --git a/examples/Solovev_kinetic_multiion_example/kinetic.dat b/examples/Solovev_kinetic_multiion_example/kinetic.dat deleted file mode 100644 index 96fd9fa3e..000000000 --- a/examples/Solovev_kinetic_multiion_example/kinetic.dat +++ /dev/null @@ -1,23 +0,0 @@ -# Kinetic profiles consistent with this Solovev equilibrium's pressure. -# Constraint: P(psi_n) = n_e*T_e + n_i*T_i must equal the equilibrium pressure -# P(psi_n) = P0*(1 - psi_n), P0 = 4.273e4 Pa (linear, zero at the edge). -# Construction (quasineutral, n_e=n_i, T_e=T_i): -# T(psi_n) = 100 + 2600*(1 - psi_n) eV (core 2700 eV, edge 100 eV) -# n(psi_n) = P0*(1 - psi_n) / (2*e*T) (core ~4.9e19; product n0*T0 = P0/2e fixed) -# Density carries the edge falloff to a small floor at psi_n=1 so T stays finite -# (avoids T->0 in the collisionality / drift-frequency formulas); the residual edge -# pressure from the floor is ~1e-6 of P0. omega_E is a free (non-pressure) profile. -psi_n n_i_m3 n_e_m3 T_i_eV T_e_eV omega_E_rads -0.00 4.939e19 4.939e19 2700.0 2700.0 1.0e4 -0.10 4.919e19 4.919e19 2440.0 2440.0 9.0e3 -0.20 4.894e19 4.894e19 2180.0 2180.0 8.0e3 -0.30 4.862e19 4.862e19 1920.0 1920.0 7.0e3 -0.40 4.820e19 4.820e19 1660.0 1660.0 6.0e3 -0.50 4.763e19 4.763e19 1400.0 1400.0 5.0e3 -0.60 4.679e19 4.679e19 1140.0 1140.0 4.0e3 -0.70 4.546e19 4.546e19 880.0 880.0 3.0e3 -0.80 4.302e19 4.302e19 620.0 620.0 2.0e3 -0.90 3.705e19 3.705e19 360.0 360.0 1.0e3 -0.95 2.899e19 2.899e19 230.0 230.0 5.0e2 -0.99 1.058e19 1.058e19 126.0 126.0 1.0e2 -1.00 2.000e18 2.000e18 100.0 100.0 0.0 diff --git a/examples/Solovev_ideal_example/forcing.dat b/examples/TJ_ideal_example/forcing.dat similarity index 64% rename from examples/Solovev_ideal_example/forcing.dat rename to examples/TJ_ideal_example/forcing.dat index 42eb94703..ce7aee865 100644 --- a/examples/Solovev_ideal_example/forcing.dat +++ b/examples/TJ_ideal_example/forcing.dat @@ -1,5 +1,6 @@ # Forcing data for perturbed equilibrium calculations # normalization: normal_field_T # Format: n m amplitude_real amplitude_imag -# Single mode test case: n=1, m=2, amplitude=1e-4 T +# n=1, m=2 and m=3, amplitude=1e-4 T on the two resonant surfaces 1 2 1e-4 0.0 +1 3 1e-4 0.0 diff --git a/examples/TJ_ideal_example/gpec.toml b/examples/TJ_ideal_example/gpec.toml new file mode 100644 index 000000000..7516d63c3 --- /dev/null +++ b/examples/TJ_ideal_example/gpec.toml @@ -0,0 +1,65 @@ +# TJ-analytic large-aspect-ratio equilibrium (q₀ = 1.6, q_a = 3.2, β_N = 0.42) with the 2/1 and 3/1 surfaces, +# stable at n = 1 with no wall (et[1] = +0.435), so a perturbed equilibrium can run without one. +# Integration-pipeline smoke test and PE regression fixture for an analytic equilibrium; the +# equilibrium is generated analytically from the embedded [TJAnalytic] section. + +[Run] +verbose = true # Log each stage's progress + +[Equilibrium] +jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) +grid_type = "core_edge_packed" # Radial grid packing type +psilow = 0.01 # Lower limit of normalized poloidal flux +psihigh = 0.99 # Upper limit of normalized poloidal flux; et[1] depends on it (+0.502 at 0.995), so it is pinned +mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) +mtheta = 256 # Number of poloidal grid points + +[ForcingTerms] +forcing_data_file = "forcing.dat" # Manual mode table (n, m, amplitude_real, amplitude_imag) +forcing_data_format = "ascii" # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires) + +[PerturbedEquilibrium] +compute_singular_coupling = true # Compute singular layer coupling metrics +regularization_width = 0.0 # Regularization width for singular surfaces (0 = disabled) + +# Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the LocalStabilityProblem settings, numerics in [Ballooning] +[LocalStability] +theta_k = 0.0 # Ballooning angle θ_k, used by both the profile scan and the α-boundary scan +alpha_boundary_scan = true # Run the first ballooning α-boundary scan (a root find per ψ surface) + +# Numerical settings of the ballooning scans +[Ballooning] +ode_tol = 1e-8 # Relative and absolute tolerance of the ballooning ODE integration behind Δ' +n_scan = 24 # Resolution of the α scan for the first boundary +max_alpha_scale = 8.0 # Upper end of the α search, as a multiple of the experimental α +alpha_tol = 1e-3 # Bisection tolerance of the critical α, in the same multiples +psi_window = [0.1, 0.99] # ψ_N range of the α scan; surfaces outside it are NaN + +[ForceFreeStates] +vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes +qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) +qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) +nn = 1 # Toroidal mode number, or consecutive numbers as a list +delta_mlow = 8 # Expands lower bound of Fourier harmonics +delta_mhigh = 8 # Expands upper bound of Fourier harmonics +mtheta_vacuum = 960 # Number of points used in splines over poloidal angle at the plasma-vacuum interface +singfac_min = 1e-6 # Fractional distance from rational q at which ideal jump enforced +integrator = "Forward" # Integration formalism, whose numerical settings are its own section: "Forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "Riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "Galerkin" (outer-region singular Galerkin Δ′) +truncate_at_rational_offset = false # Truncate integration at (last_rational_q + rational_offset_fraction)/n; false for analytic equilibria + +# Numerical settings of the forward integrator +[Forward] +eulerlagrange_tolerance = 1e-9 # Relative tolerance for ODE integration of Euler-Lagrange equations + +# TJ-analytic equilibrium parameters (cf. R. Fitzpatrick, TJ code, https://github.com/rfitzp/TJ); this section selects the +# equilibrium. Keep p0 at or below 0.2: et[1] stays near +0.42 up to there and is not monotonic above it. +[TJAnalytic] +r0 = 10.0 # Major radius R₀ [m] (centerline radius of the magnetic axis) +a = 1.0 # Minor radius a [m] (plasma half-width at the midplane; here ε = 0.1) +q0 = 1.6 # On-axis safety factor q(r=0) +qa = 3.2 # Edge safety factor q(r=a) +p0 = 0.1 # Normalized on-axis pressure +pressure_peaking = 2.0 # Pressure-profile peaking exponent in p₂(r) = p0·(1−r²)^μ +B0 = 1.0 # On-axis toroidal field [T] +ma = 128 # TJ-analytic internal radial grid resolution (shape-ODE nodes) +mtheta = 128 # TJ-analytic internal poloidal grid resolution (θ-spline nodes) diff --git a/examples/Solovev_kinetic_NTV_example/run_example.jl b/examples/TJ_ideal_example/run_example.jl similarity index 100% rename from examples/Solovev_kinetic_NTV_example/run_example.jl rename to examples/TJ_ideal_example/run_example.jl diff --git a/examples/Solovev_kinetic_NTV_example/forcing.dat b/examples/TJ_kinetic_NTV_example/forcing.dat similarity index 64% rename from examples/Solovev_kinetic_NTV_example/forcing.dat rename to examples/TJ_kinetic_NTV_example/forcing.dat index 42eb94703..ce7aee865 100644 --- a/examples/Solovev_kinetic_NTV_example/forcing.dat +++ b/examples/TJ_kinetic_NTV_example/forcing.dat @@ -1,5 +1,6 @@ # Forcing data for perturbed equilibrium calculations # normalization: normal_field_T # Format: n m amplitude_real amplitude_imag -# Single mode test case: n=1, m=2, amplitude=1e-4 T +# n=1, m=2 and m=3, amplitude=1e-4 T on the two resonant surfaces 1 2 1e-4 0.0 +1 3 1e-4 0.0 diff --git a/examples/TJ_kinetic_NTV_example/gpec.toml b/examples/TJ_kinetic_NTV_example/gpec.toml new file mode 100644 index 000000000..b9d3e6659 --- /dev/null +++ b/examples/TJ_kinetic_NTV_example/gpec.toml @@ -0,0 +1,58 @@ +# TJ-analytic large-aspect-ratio equilibrium (q₀ = 1.6, q_a = 3.2, β_N = 0.42) with the 2/1 and 3/1 surfaces, +# stable at n = 1 with no wall (et[1] = +0.435), so a perturbed equilibrium can run without one. +# Adds a [KineticForces] section, so after the plasma response is computed the neoclassical toroidal +# viscosity torque is integrated over ψ as a post-processing diagnostic (adaptive quadrature paneled at +# the rational surfaces). The kinetic profiles (kinetic.dat) are consistent with this equilibrium's pressure. + +[Run] +verbose = true # Log each stage's progress + +[Equilibrium] +kinetic_file = "kinetic.dat" # Kinetic-profile file (HDF5 GPEC kinetic schema or 6-column ASCII), read once into the equilibrium +jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) +grid_type = "core_edge_packed" # Radial grid packing type +psilow = 0.01 # Lower limit of normalized poloidal flux +psihigh = 0.99 # Upper limit of normalized poloidal flux; et[1] depends on it (+0.502 at 0.995), so it is pinned +mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) +mtheta = 256 # Number of poloidal grid points + +[ForcingTerms] +forcing_data_file = "forcing.dat" # Manual mode table (n, m, amplitude_real, amplitude_imag) +forcing_data_format = "ascii" # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires) + +[PerturbedEquilibrium] +compute_singular_coupling = true # Compute singular layer coupling metrics +regularization_width = 0.05 # Regularization width for singular surfaces (0 = disabled) + +[ForceFreeStates] +vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes +qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) +qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) +nn = 1 # Toroidal mode number, or consecutive numbers as a list +delta_mlow = 8 # Expands lower bound of Fourier harmonics +delta_mhigh = 8 # Expands upper bound of Fourier harmonics +mtheta_vacuum = 960 # Number of points used in splines over poloidal angle at the plasma-vacuum interface +singfac_min = 1e-6 # Fractional distance from rational q at which ideal jump enforced +integrator = "Forward" # Integration formalism, whose numerical settings are its own section: "Forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "Riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "Galerkin" (outer-region singular Galerkin Δ′) +truncate_at_rational_offset = false # Truncate integration at (last_rational_q + rational_offset_fraction)/n; false for analytic equilibria + +# Numerical settings of the forward integrator +[Forward] +eulerlagrange_tolerance = 1e-9 # Relative tolerance for ODE integration of Euler-Lagrange equations + +# TJ-analytic equilibrium parameters (cf. R. Fitzpatrick, TJ code, https://github.com/rfitzp/TJ); this section selects the +# equilibrium. Keep p0 at or below 0.2: et[1] stays near +0.42 up to there and is not monotonic above it. +[TJAnalytic] +r0 = 10.0 # Major radius R₀ [m] (centerline radius of the magnetic axis) +a = 1.0 # Minor radius a [m] (plasma half-width at the midplane; here ε = 0.1) +q0 = 1.6 # On-axis safety factor q(r=0) +qa = 3.2 # Edge safety factor q(r=a) +p0 = 0.1 # Normalized on-axis pressure +pressure_peaking = 2.0 # Pressure-profile peaking exponent in p₂(r) = p0·(1−r²)^μ +B0 = 1.0 # On-axis toroidal field [T] +ma = 128 # TJ-analytic internal radial grid resolution (shape-ODE nodes) +mtheta = 128 # TJ-analytic internal poloidal grid resolution (θ-spline nodes) + +[KineticForces] +fcgl_flag = true # Fluid Chew–Goldberger–Low energy (l = 0 fluid limit) alongside the default fgar method +axis_validity_suppression = true # Suppress kinetic terms inside the near-axis orbit-width validity limit diff --git a/examples/TJ_kinetic_NTV_example/kinetic.dat b/examples/TJ_kinetic_NTV_example/kinetic.dat new file mode 100644 index 000000000..e3a393077 --- /dev/null +++ b/examples/TJ_kinetic_NTV_example/kinetic.dat @@ -0,0 +1,45 @@ +# Kinetic profiles consistent with the TJ-analytic equilibrium's pressure (P₀ = 784 Pa). +# Constraint: P(psi_n) = n_e*T_e + n_i*T_i, with n_e = n_i and T_e = T_i = 20 + 180*(1 - psi_n) eV. +# Omega_E = 2e4*(1 - psi_n) rad/s. +# psi_n n_i n_e T_i[eV] T_e[eV] omega_E +0.000000 1.223086e+19 1.223086e+19 200.0000 200.0000 20000.0000 +0.025000 1.222817e+19 1.222817e+19 195.5000 195.5000 19500.0000 +0.050000 1.203522e+19 1.203522e+19 191.0000 191.0000 19000.0000 +0.075000 1.183797e+19 1.183797e+19 186.5000 186.5000 18500.0000 +0.100000 1.163625e+19 1.163625e+19 182.0000 182.0000 18000.0000 +0.125000 1.142988e+19 1.142988e+19 177.5000 177.5000 17500.0000 +0.150000 1.121868e+19 1.121868e+19 173.0000 173.0000 17000.0000 +0.175000 1.100245e+19 1.100245e+19 168.5000 168.5000 16500.0000 +0.200000 1.078097e+19 1.078097e+19 164.0000 164.0000 16000.0000 +0.225000 1.055403e+19 1.055403e+19 159.5000 159.5000 15500.0000 +0.250000 1.032142e+19 1.032142e+19 155.0000 155.0000 15000.0000 +0.275000 1.008285e+19 1.008285e+19 150.5000 150.5000 14500.0000 +0.300000 9.838102e+18 9.838102e+18 146.0000 146.0000 14000.0000 +0.325000 9.586866e+18 9.586866e+18 141.5000 141.5000 13500.0000 +0.350000 9.328821e+18 9.328821e+18 137.0000 137.0000 13000.0000 +0.375000 9.063714e+18 9.063714e+18 132.5000 132.5000 12500.0000 +0.400000 8.791260e+18 8.791260e+18 128.0000 128.0000 12000.0000 +0.425000 8.511117e+18 8.511117e+18 123.5000 123.5000 11500.0000 +0.450000 8.222884e+18 8.222884e+18 119.0000 119.0000 11000.0000 +0.475000 7.926125e+18 7.926125e+18 114.5000 114.5000 10500.0000 +0.500000 7.620476e+18 7.620476e+18 110.0000 110.0000 10000.0000 +0.525000 7.305553e+18 7.305553e+18 105.5000 105.5000 9500.0000 +0.550000 6.980924e+18 6.980924e+18 101.0000 101.0000 9000.0000 +0.575000 6.646101e+18 6.646101e+18 96.5000 96.5000 8500.0000 +0.600000 6.300621e+18 6.300621e+18 92.0000 92.0000 8000.0000 +0.625000 5.944099e+18 5.944099e+18 87.5000 87.5000 7500.0000 +0.650000 5.576149e+18 5.576149e+18 83.0000 83.0000 7000.0000 +0.675000 5.196380e+18 5.196380e+18 78.5000 78.5000 6500.0000 +0.700000 4.804420e+18 4.804420e+18 74.0000 74.0000 6000.0000 +0.725000 4.399955e+18 4.399955e+18 69.5000 69.5000 5500.0000 +0.750000 3.982956e+18 3.982956e+18 65.0000 65.0000 5000.0000 +0.775000 3.553750e+18 3.553750e+18 60.5000 60.5000 4500.0000 +0.800000 3.113088e+18 3.113088e+18 56.0000 56.0000 4000.0000 +0.825000 2.662383e+18 2.662383e+18 51.5000 51.5000 3500.0000 +0.850000 2.204114e+18 2.204114e+18 47.0000 47.0000 3000.0000 +0.875000 1.742718e+18 1.742718e+18 42.5000 42.5000 2500.0000 +0.900000 1.285949e+18 1.285949e+18 38.0000 38.0000 2000.0000 +0.925000 8.470999e+17 8.470999e+17 33.5000 33.5000 1500.0000 +0.950000 4.497124e+17 4.497124e+17 29.0000 29.0000 1000.0000 +0.975000 1.378556e+17 1.378556e+17 24.5000 24.5000 500.0000 +1.000000 2.762506e+16 2.762506e+16 20.0000 20.0000 0.0000 diff --git a/examples/Solovev_kinetic_multiion_example/run_example.jl b/examples/TJ_kinetic_NTV_example/run_example.jl similarity index 100% rename from examples/Solovev_kinetic_multiion_example/run_example.jl rename to examples/TJ_kinetic_NTV_example/run_example.jl diff --git a/examples/Solovev_kinetic_calculated_example/forcing.dat b/examples/TJ_kinetic_calculated_example/forcing.dat similarity index 64% rename from examples/Solovev_kinetic_calculated_example/forcing.dat rename to examples/TJ_kinetic_calculated_example/forcing.dat index 42eb94703..ce7aee865 100644 --- a/examples/Solovev_kinetic_calculated_example/forcing.dat +++ b/examples/TJ_kinetic_calculated_example/forcing.dat @@ -1,5 +1,6 @@ # Forcing data for perturbed equilibrium calculations # normalization: normal_field_T # Format: n m amplitude_real amplitude_imag -# Single mode test case: n=1, m=2, amplitude=1e-4 T +# n=1, m=2 and m=3, amplitude=1e-4 T on the two resonant surfaces 1 2 1e-4 0.0 +1 3 1e-4 0.0 diff --git a/examples/TJ_kinetic_calculated_example/gpec.toml b/examples/TJ_kinetic_calculated_example/gpec.toml new file mode 100644 index 000000000..3b8a45e27 --- /dev/null +++ b/examples/TJ_kinetic_calculated_example/gpec.toml @@ -0,0 +1,70 @@ +# TJ-analytic large-aspect-ratio equilibrium (q₀ = 1.6, q_a = 3.2, β_N = 0.42) with the 2/1 and 3/1 surfaces, +# stable at n = 1 with no wall (et[1] = +0.435), so a perturbed equilibrium can run without one. +# Exercises the real KineticForces (NTV) physics path at full strength, with the main ion declared as a +# 50/50 D-T mix so the self-consistent kinetic matrices are accumulated over the multi-ion species loop +# (pure plasma, n_i = n_e, so no impurity species arises). The kinetic profiles (kinetic.dat) are consistent +# with this equilibrium's pressure. Also the only deck that runs [PerturbedEquilibrium] on a kinetic solve, +# so it builds the ψ-resolved torque response matrices (T_xe on the applied field, T_source for the forcing file). + +[Run] +verbose = true # Log each stage's progress + +[Equilibrium] +kinetic_file = "kinetic.dat" # Kinetic-profile file (HDF5 GPEC kinetic schema or 6-column ASCII), read once into the equilibrium +jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) +grid_type = "core_edge_packed" # Radial grid packing type +psilow = 0.01 # Lower limit of normalized poloidal flux +psihigh = 0.99 # Upper limit of normalized poloidal flux; et[1] depends on it (+0.502 at 0.995), so it is pinned +mpsi = 32 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) +mtheta = 256 # Number of poloidal grid points + +[[Equilibrium.ion_species]] +z = 1 # Species charge (fundamental units) +m = 2 # Species mass (proton masses) — Deuterium +fraction = 0.5 # This species' share of the total main-ion density n_i + +[[Equilibrium.ion_species]] +z = 1 # Species charge (fundamental units) +m = 3 # Species mass (proton masses) — Tritium +fraction = 0.5 # This species' share of the total main-ion density n_i + +[ForcingTerms] +forcing_data_file = "forcing.dat" # Manual mode table (n, m, amplitude_real, amplitude_imag) +forcing_data_format = "ascii" # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires) + +[PerturbedEquilibrium] +compute_singular_coupling = true # Compute singular layer coupling metrics +compute_torque_response = true # Build the ψ-resolved torque response matrices T_xe and T_source + +[ForceFreeStates] +vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes +qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) +qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) +nn = 1 # Toroidal mode number, or consecutive numbers as a list +delta_mlow = 8 # Expands lower bound of Fourier harmonics +delta_mhigh = 8 # Expands upper bound of Fourier harmonics +mtheta_vacuum = 64 # Number of points used in splines over poloidal angle at the plasma-vacuum interface +singfac_min = 1e-6 # Fractional distance from rational q at which ideal jump enforced +kinetic = "KineticForces" # Kinetic source: the drift-kinetic matrices of the [KineticForces] model (omit for an ideal run) +integrator = "Forward" # Integration formalism, whose numerical settings are its own section: "Forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "Riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "Galerkin" (outer-region singular Galerkin Δ′) +truncate_at_rational_offset = false # Truncate integration at (last_rational_q + rational_offset_fraction)/n; false for analytic equilibria + +# Numerical settings of the forward integrator +[Forward] +eulerlagrange_tolerance = 1e-9 # Relative tolerance for ODE integration of Euler-Lagrange equations + +# TJ-analytic equilibrium parameters (cf. R. Fitzpatrick, TJ code, https://github.com/rfitzp/TJ); this section selects the +# equilibrium. Keep p0 at or below 0.2: et[1] stays near +0.42 up to there and is not monotonic above it. +[TJAnalytic] +r0 = 10.0 # Major radius R₀ [m] (centerline radius of the magnetic axis) +a = 1.0 # Minor radius a [m] (plasma half-width at the midplane; here ε = 0.1) +q0 = 1.6 # On-axis safety factor q(r=0) +qa = 3.2 # Edge safety factor q(r=a) +p0 = 0.1 # Normalized on-axis pressure +pressure_peaking = 2.0 # Pressure-profile peaking exponent in p₂(r) = p0·(1−r²)^μ +B0 = 1.0 # On-axis toroidal field [T] +ma = 128 # TJ-analytic internal radial grid resolution (shape-ODE nodes) +mtheta = 128 # TJ-analytic internal poloidal grid resolution (θ-spline nodes) + +[KineticForces] +axis_validity_suppression = true # Suppress kinetic terms inside the near-axis orbit-width validity limit diff --git a/examples/TJ_kinetic_calculated_example/kinetic.dat b/examples/TJ_kinetic_calculated_example/kinetic.dat new file mode 100644 index 000000000..e3a393077 --- /dev/null +++ b/examples/TJ_kinetic_calculated_example/kinetic.dat @@ -0,0 +1,45 @@ +# Kinetic profiles consistent with the TJ-analytic equilibrium's pressure (P₀ = 784 Pa). +# Constraint: P(psi_n) = n_e*T_e + n_i*T_i, with n_e = n_i and T_e = T_i = 20 + 180*(1 - psi_n) eV. +# Omega_E = 2e4*(1 - psi_n) rad/s. +# psi_n n_i n_e T_i[eV] T_e[eV] omega_E +0.000000 1.223086e+19 1.223086e+19 200.0000 200.0000 20000.0000 +0.025000 1.222817e+19 1.222817e+19 195.5000 195.5000 19500.0000 +0.050000 1.203522e+19 1.203522e+19 191.0000 191.0000 19000.0000 +0.075000 1.183797e+19 1.183797e+19 186.5000 186.5000 18500.0000 +0.100000 1.163625e+19 1.163625e+19 182.0000 182.0000 18000.0000 +0.125000 1.142988e+19 1.142988e+19 177.5000 177.5000 17500.0000 +0.150000 1.121868e+19 1.121868e+19 173.0000 173.0000 17000.0000 +0.175000 1.100245e+19 1.100245e+19 168.5000 168.5000 16500.0000 +0.200000 1.078097e+19 1.078097e+19 164.0000 164.0000 16000.0000 +0.225000 1.055403e+19 1.055403e+19 159.5000 159.5000 15500.0000 +0.250000 1.032142e+19 1.032142e+19 155.0000 155.0000 15000.0000 +0.275000 1.008285e+19 1.008285e+19 150.5000 150.5000 14500.0000 +0.300000 9.838102e+18 9.838102e+18 146.0000 146.0000 14000.0000 +0.325000 9.586866e+18 9.586866e+18 141.5000 141.5000 13500.0000 +0.350000 9.328821e+18 9.328821e+18 137.0000 137.0000 13000.0000 +0.375000 9.063714e+18 9.063714e+18 132.5000 132.5000 12500.0000 +0.400000 8.791260e+18 8.791260e+18 128.0000 128.0000 12000.0000 +0.425000 8.511117e+18 8.511117e+18 123.5000 123.5000 11500.0000 +0.450000 8.222884e+18 8.222884e+18 119.0000 119.0000 11000.0000 +0.475000 7.926125e+18 7.926125e+18 114.5000 114.5000 10500.0000 +0.500000 7.620476e+18 7.620476e+18 110.0000 110.0000 10000.0000 +0.525000 7.305553e+18 7.305553e+18 105.5000 105.5000 9500.0000 +0.550000 6.980924e+18 6.980924e+18 101.0000 101.0000 9000.0000 +0.575000 6.646101e+18 6.646101e+18 96.5000 96.5000 8500.0000 +0.600000 6.300621e+18 6.300621e+18 92.0000 92.0000 8000.0000 +0.625000 5.944099e+18 5.944099e+18 87.5000 87.5000 7500.0000 +0.650000 5.576149e+18 5.576149e+18 83.0000 83.0000 7000.0000 +0.675000 5.196380e+18 5.196380e+18 78.5000 78.5000 6500.0000 +0.700000 4.804420e+18 4.804420e+18 74.0000 74.0000 6000.0000 +0.725000 4.399955e+18 4.399955e+18 69.5000 69.5000 5500.0000 +0.750000 3.982956e+18 3.982956e+18 65.0000 65.0000 5000.0000 +0.775000 3.553750e+18 3.553750e+18 60.5000 60.5000 4500.0000 +0.800000 3.113088e+18 3.113088e+18 56.0000 56.0000 4000.0000 +0.825000 2.662383e+18 2.662383e+18 51.5000 51.5000 3500.0000 +0.850000 2.204114e+18 2.204114e+18 47.0000 47.0000 3000.0000 +0.875000 1.742718e+18 1.742718e+18 42.5000 42.5000 2500.0000 +0.900000 1.285949e+18 1.285949e+18 38.0000 38.0000 2000.0000 +0.925000 8.470999e+17 8.470999e+17 33.5000 33.5000 1500.0000 +0.950000 4.497124e+17 4.497124e+17 29.0000 29.0000 1000.0000 +0.975000 1.378556e+17 1.378556e+17 24.5000 24.5000 500.0000 +1.000000 2.762506e+16 2.762506e+16 20.0000 20.0000 0.0000 diff --git a/examples/Solovev_kinetic_multiion_example/forcing.dat b/examples/TJ_kinetic_multiion_example/forcing.dat similarity index 64% rename from examples/Solovev_kinetic_multiion_example/forcing.dat rename to examples/TJ_kinetic_multiion_example/forcing.dat index 42eb94703..ce7aee865 100644 --- a/examples/Solovev_kinetic_multiion_example/forcing.dat +++ b/examples/TJ_kinetic_multiion_example/forcing.dat @@ -1,5 +1,6 @@ # Forcing data for perturbed equilibrium calculations # normalization: normal_field_T # Format: n m amplitude_real amplitude_imag -# Single mode test case: n=1, m=2, amplitude=1e-4 T +# n=1, m=2 and m=3, amplitude=1e-4 T on the two resonant surfaces 1 2 1e-4 0.0 +1 3 1e-4 0.0 diff --git a/examples/TJ_kinetic_multiion_example/gpec.toml b/examples/TJ_kinetic_multiion_example/gpec.toml new file mode 100644 index 000000000..5b32e39f7 --- /dev/null +++ b/examples/TJ_kinetic_multiion_example/gpec.toml @@ -0,0 +1,69 @@ +# TJ-analytic large-aspect-ratio equilibrium (q₀ = 1.6, q_a = 3.2, β_N = 0.42) with the 2/1 and 3/1 surfaces, +# stable at n = 1 with no wall (et[1] = +0.435), so a perturbed equilibrium can run without one. +# Multi-ion NTV: the main ion is a 50/50 D-T mix (n_i = total), so the summed NTV runs over D + T + +# the neutrality-closing impurity + electrons. The kinetic profiles (kinetic.dat) are consistent with +# this equilibrium's pressure. + +[Run] +verbose = true # Log each stage's progress + +[Equilibrium] +kinetic_file = "kinetic.dat" # Kinetic-profile file (HDF5 GPEC kinetic schema or 6-column ASCII), read once into the equilibrium +jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) +grid_type = "core_edge_packed" # Radial grid packing type +psilow = 0.01 # Lower limit of normalized poloidal flux +psihigh = 0.99 # Upper limit of normalized poloidal flux; et[1] depends on it (+0.502 at 0.995), so it is pinned +mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) +mtheta = 256 # Number of poloidal grid points +impurity_charge = 6 # Impurity charge (fundamental units); closes quasineutrality +impurity_mass = 12 # Impurity mass (proton masses) + +[[Equilibrium.ion_species]] +z = 1 # Species charge (fundamental units) +m = 2 # Species mass (proton masses) — Deuterium +fraction = 0.5 # This species' share of the total main-ion density n_i + +[[Equilibrium.ion_species]] +z = 1 # Species charge (fundamental units) +m = 3 # Species mass (proton masses) — Tritium +fraction = 0.5 # This species' share of the total main-ion density n_i + +[ForcingTerms] +forcing_data_file = "forcing.dat" # Manual mode table (n, m, amplitude_real, amplitude_imag) +forcing_data_format = "ascii" # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires) + +[PerturbedEquilibrium] +compute_singular_coupling = true # Compute singular layer coupling metrics +regularization_width = 0.05 # Regularization width for singular surfaces (0 = disabled) + +[ForceFreeStates] +vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes +qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) +qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) +nn = 1 # Toroidal mode number, or consecutive numbers as a list +delta_mlow = 8 # Expands lower bound of Fourier harmonics +delta_mhigh = 8 # Expands upper bound of Fourier harmonics +mtheta_vacuum = 960 # Number of points used in splines over poloidal angle at the plasma-vacuum interface +singfac_min = 1e-6 # Fractional distance from rational q at which ideal jump enforced +integrator = "Forward" # Integration formalism, whose numerical settings are its own section: "Forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "Riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "Galerkin" (outer-region singular Galerkin Δ′) +truncate_at_rational_offset = false # Truncate integration at (last_rational_q + rational_offset_fraction)/n; false for analytic equilibria + +# Numerical settings of the forward integrator +[Forward] +eulerlagrange_tolerance = 1e-9 # Relative tolerance for ODE integration of Euler-Lagrange equations + +# TJ-analytic equilibrium parameters (cf. R. Fitzpatrick, TJ code, https://github.com/rfitzp/TJ); this section selects the +# equilibrium. Keep p0 at or below 0.2: et[1] stays near +0.42 up to there and is not monotonic above it. +[TJAnalytic] +r0 = 10.0 # Major radius R₀ [m] (centerline radius of the magnetic axis) +a = 1.0 # Minor radius a [m] (plasma half-width at the midplane; here ε = 0.1) +q0 = 1.6 # On-axis safety factor q(r=0) +qa = 3.2 # Edge safety factor q(r=a) +p0 = 0.1 # Normalized on-axis pressure +pressure_peaking = 2.0 # Pressure-profile peaking exponent in p₂(r) = p0·(1−r²)^μ +B0 = 1.0 # On-axis toroidal field [T] +ma = 128 # TJ-analytic internal radial grid resolution (shape-ODE nodes) +mtheta = 128 # TJ-analytic internal poloidal grid resolution (θ-spline nodes) + +[KineticForces] +electron = true # Add electron NTV in addition to the ion species diff --git a/examples/TJ_kinetic_multiion_example/kinetic.dat b/examples/TJ_kinetic_multiion_example/kinetic.dat new file mode 100644 index 000000000..e3a393077 --- /dev/null +++ b/examples/TJ_kinetic_multiion_example/kinetic.dat @@ -0,0 +1,45 @@ +# Kinetic profiles consistent with the TJ-analytic equilibrium's pressure (P₀ = 784 Pa). +# Constraint: P(psi_n) = n_e*T_e + n_i*T_i, with n_e = n_i and T_e = T_i = 20 + 180*(1 - psi_n) eV. +# Omega_E = 2e4*(1 - psi_n) rad/s. +# psi_n n_i n_e T_i[eV] T_e[eV] omega_E +0.000000 1.223086e+19 1.223086e+19 200.0000 200.0000 20000.0000 +0.025000 1.222817e+19 1.222817e+19 195.5000 195.5000 19500.0000 +0.050000 1.203522e+19 1.203522e+19 191.0000 191.0000 19000.0000 +0.075000 1.183797e+19 1.183797e+19 186.5000 186.5000 18500.0000 +0.100000 1.163625e+19 1.163625e+19 182.0000 182.0000 18000.0000 +0.125000 1.142988e+19 1.142988e+19 177.5000 177.5000 17500.0000 +0.150000 1.121868e+19 1.121868e+19 173.0000 173.0000 17000.0000 +0.175000 1.100245e+19 1.100245e+19 168.5000 168.5000 16500.0000 +0.200000 1.078097e+19 1.078097e+19 164.0000 164.0000 16000.0000 +0.225000 1.055403e+19 1.055403e+19 159.5000 159.5000 15500.0000 +0.250000 1.032142e+19 1.032142e+19 155.0000 155.0000 15000.0000 +0.275000 1.008285e+19 1.008285e+19 150.5000 150.5000 14500.0000 +0.300000 9.838102e+18 9.838102e+18 146.0000 146.0000 14000.0000 +0.325000 9.586866e+18 9.586866e+18 141.5000 141.5000 13500.0000 +0.350000 9.328821e+18 9.328821e+18 137.0000 137.0000 13000.0000 +0.375000 9.063714e+18 9.063714e+18 132.5000 132.5000 12500.0000 +0.400000 8.791260e+18 8.791260e+18 128.0000 128.0000 12000.0000 +0.425000 8.511117e+18 8.511117e+18 123.5000 123.5000 11500.0000 +0.450000 8.222884e+18 8.222884e+18 119.0000 119.0000 11000.0000 +0.475000 7.926125e+18 7.926125e+18 114.5000 114.5000 10500.0000 +0.500000 7.620476e+18 7.620476e+18 110.0000 110.0000 10000.0000 +0.525000 7.305553e+18 7.305553e+18 105.5000 105.5000 9500.0000 +0.550000 6.980924e+18 6.980924e+18 101.0000 101.0000 9000.0000 +0.575000 6.646101e+18 6.646101e+18 96.5000 96.5000 8500.0000 +0.600000 6.300621e+18 6.300621e+18 92.0000 92.0000 8000.0000 +0.625000 5.944099e+18 5.944099e+18 87.5000 87.5000 7500.0000 +0.650000 5.576149e+18 5.576149e+18 83.0000 83.0000 7000.0000 +0.675000 5.196380e+18 5.196380e+18 78.5000 78.5000 6500.0000 +0.700000 4.804420e+18 4.804420e+18 74.0000 74.0000 6000.0000 +0.725000 4.399955e+18 4.399955e+18 69.5000 69.5000 5500.0000 +0.750000 3.982956e+18 3.982956e+18 65.0000 65.0000 5000.0000 +0.775000 3.553750e+18 3.553750e+18 60.5000 60.5000 4500.0000 +0.800000 3.113088e+18 3.113088e+18 56.0000 56.0000 4000.0000 +0.825000 2.662383e+18 2.662383e+18 51.5000 51.5000 3500.0000 +0.850000 2.204114e+18 2.204114e+18 47.0000 47.0000 3000.0000 +0.875000 1.742718e+18 1.742718e+18 42.5000 42.5000 2500.0000 +0.900000 1.285949e+18 1.285949e+18 38.0000 38.0000 2000.0000 +0.925000 8.470999e+17 8.470999e+17 33.5000 33.5000 1500.0000 +0.950000 4.497124e+17 4.497124e+17 29.0000 29.0000 1000.0000 +0.975000 1.378556e+17 1.378556e+17 24.5000 24.5000 500.0000 +1.000000 2.762506e+16 2.762506e+16 20.0000 20.0000 0.0000 diff --git a/examples/TJ_kinetic_multiion_example/run_example.jl b/examples/TJ_kinetic_multiion_example/run_example.jl new file mode 100644 index 000000000..fa0a9e53d --- /dev/null +++ b/examples/TJ_kinetic_multiion_example/run_example.jl @@ -0,0 +1,4 @@ +using Pkg; +Pkg.activate(joinpath(@__DIR__, "../..")) +using GeneralizedPerturbedEquilibrium +GeneralizedPerturbedEquilibrium.main([dirname(@__FILE__)]) diff --git a/regression-harness/cases/solovev_kinetic_nuzero.toml b/regression-harness/cases/solovev_kinetic_nuzero.toml index 2a261ab33..b7bb8615d 100644 --- a/regression-harness/cases/solovev_kinetic_nuzero.toml +++ b/regression-harness/cases/solovev_kinetic_nuzero.toml @@ -1,17 +1,20 @@ # Regression case: Solovev n=1 collisionless kinetic run (nutype = "zero"). -# Reuses the Solovev_kinetic_calculated_example deck (50/50 D-T multi-ion) via -# [overrides] at a mode-converged band; locks the collisionless real-x-space -# energy-integral eigenvalue plus equilibrium and singular-surface anchors. +# Reuses the nuzero test fixture via [overrides] at a mode-converged band with 50/50 D-T +# multi-ion matrices; locks the collisionless real-x-space energy-integral eigenvalue +# plus equilibrium and singular-surface anchors. # 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 = "solovev_kinetic_nuzero" -description = "Solovev, n=1, kinetic = 'KineticForces', nutype=zero (collisionless energy integral), 50/50 D-T multi-ion kinetic matrices, mode-converged band (delta_m=8, mpert=32). Reuses the calculated example deck via override; pressure-consistent kinetic.dat." -example_dir = "examples/Solovev_kinetic_calculated_example" +description = "Solovev, n=1, kinetic = 'KineticForces', nutype=zero (collisionless energy integral), 50/50 D-T multi-ion kinetic matrices, mode-converged band (delta_m=8, mpert=32). Reuses the nuzero test fixture via override; pressure-consistent kinetic.dat." +example_dir = "test/test_data/regression_solovev_kinetic_nuzero" -# Run the shared calculated deck collisionless (no separate nuzero example directory). +# The fixture is coarse for test speed; widen it to the converged band and the D-T ion mix. [overrides] -"KineticForces.nutype" = "zero" +"Equilibrium.ion_species" = [{ z = 1, m = 2, fraction = 0.5 }, { z = 1, m = 3, fraction = 0.5 }] +"ForceFreeStates.delta_mlow" = 8 +"ForceFreeStates.delta_mhigh" = 8 +"ForceFreeStates.singfac_min" = 1e-6 # Energies — leading eigenvalues [quantities.et_real] diff --git a/regression-harness/cases/solovev_kinetic_calculated.toml b/regression-harness/cases/tj_kinetic_calculated.toml similarity index 92% rename from regression-harness/cases/solovev_kinetic_calculated.toml rename to regression-harness/cases/tj_kinetic_calculated.toml index 973002241..af99a5f6d 100644 --- a/regression-harness/cases/solovev_kinetic_calculated.toml +++ b/regression-harness/cases/tj_kinetic_calculated.toml @@ -1,13 +1,13 @@ -# Regression case: Solovev n=1 with the "calculated" kinetic path at full strength. +# Regression case: TJ-analytic n=1 with the "calculated" kinetic path at full strength. # Guards the real KineticForces (PENTRC-style) physics path against drift, with the main # ion declared as a 50/50 D-T mix so the self-consistent kinetic-matrix species loop # (multi-ion accumulation) is pinned end-to-end. 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 = "solovev_kinetic_calculated" -description = "Solovev, n=1, kinetic = 'KineticForces', 50/50 D-T multi-ion kinetic matrices, mode-converged band (delta_m=8, mpert=32); pressure-consistent kinetic.dat (n·T tied to the Solovev P₀≈4.27e4 Pa)." -example_dir = "examples/Solovev_kinetic_calculated_example" +name = "tj_kinetic_calculated" +description = "TJ-analytic, n=1, no wall, kinetic = 'KineticForces', 50/50 D-T multi-ion kinetic matrices; pressure-consistent kinetic.dat (P₀≈784 Pa)." +example_dir = "examples/TJ_kinetic_calculated_example" # Energies — leading eigenvalues [quantities.et_real] diff --git a/regression-harness/cases/solovev_kinetic_multiion.toml b/regression-harness/cases/tj_kinetic_multiion.toml similarity index 86% rename from regression-harness/cases/solovev_kinetic_multiion.toml rename to regression-harness/cases/tj_kinetic_multiion.toml index faf23acfc..c6238b821 100644 --- a/regression-harness/cases/solovev_kinetic_multiion.toml +++ b/regression-harness/cases/tj_kinetic_multiion.toml @@ -1,12 +1,12 @@ -# Regression case: Solovev n=1 + PE + MULTI-ION NTV torque. The Solovev main ion is modelled as a +# Regression case: TJ-analytic n=1 + PE + MULTI-ION NTV torque. The main ion is modelled as a # 50/50 D-T mix, so the KineticForces quadrature sums over D + T + electrons (this fixture is a pure # plasma, n_e = n_i, so no impurity species is present). Tracks the summed total AND each per-species # contribution, pinning the multi-species loop (resolve_ntv_species + combine_species_states) and the # full-composition Zeff. The impurity-species path is covered by the unit tests + the TC-24 case. [case] -name = "solovev_kinetic_multiion" -description = "Solovev n=1, ideal + PE + multi-ion NTV (D + T + electron summed)" -example_dir = "examples/Solovev_kinetic_multiion_example" +name = "tj_kinetic_multiion" +description = "TJ-analytic n=1, no wall, ideal + PE + multi-ion NTV (D + T + electron summed)" +example_dir = "examples/TJ_kinetic_multiion_example" # Summed total NTV over all species — the headline multi-ion quantity. [quantities.ntv_total] diff --git a/regression-harness/cases/solovev_kinetic_ntv.toml b/regression-harness/cases/tj_kinetic_ntv.toml similarity index 89% rename from regression-harness/cases/solovev_kinetic_ntv.toml rename to regression-harness/cases/tj_kinetic_ntv.toml index d0a95c529..e4310bc8d 100644 --- a/regression-harness/cases/solovev_kinetic_ntv.toml +++ b/regression-harness/cases/tj_kinetic_ntv.toml @@ -1,12 +1,12 @@ -# Regression case: Solovev n=1 ideal stability + perturbed equilibrium + post-PE NTV torque. +# Regression case: TJ-analytic n=1 ideal stability + perturbed equilibrium + post-PE NTV torque. # The only case exercising the KineticForces outer ψ torque quadrature (fgar method): # rational-surface paneling, maxevals_psi cap, and the atol_psi/rtol_psi defaults. # 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 = "solovev_kinetic_ntv" -description = "Solovev analytical equilibrium, n=1, ideal + PE + NTV torque quadrature" -example_dir = "examples/Solovev_kinetic_NTV_example" +name = "tj_kinetic_ntv" +description = "TJ-analytic equilibrium, n=1, no wall, ideal + PE + NTV torque quadrature" +example_dir = "examples/TJ_kinetic_NTV_example" # NTV torque — the primary tracked quantity of this case. [Re, Im] pair; the ψ quadrature # converges to rtol_psi=1e-2, but the result is deterministic for fixed code + inputs, so diff --git a/regression-harness/cases/tj_n1.toml b/regression-harness/cases/tj_n1.toml new file mode 100644 index 000000000..ad9d79974 --- /dev/null +++ b/regression-harness/cases/tj_n1.toml @@ -0,0 +1,225 @@ +# Regression case: TJ-analytic n=1 ideal stability and perturbed equilibrium, no wall. +# Tracks the leading free-boundary energies, the stability quantities and the singular coupling of +# a stable analytic equilibrium. 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 = "tj_n1" +description = "TJ-analytic equilibrium, n=1, no wall, ideal stability + perturbed equilibrium" +example_dir = "examples/TJ_ideal_example" + +# Energies — leading generalized (W,N) pencil eigenvalues at the final truncation (psilim): power-normalized (⟨|ξ|²⟩ = 1 metric) and invariant to the working-coordinate Jacobian. +[quantities.et_real] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies" +type = "complex_vector" +extract = "real_first" +label = "total energy Re(et[1])" +noise_threshold = 1e-10 +order = 10 + +[quantities.et_imag] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies" +type = "complex_vector" +extract = "imag_first" +label = "total energy Im(et[1])" +noise_threshold = 1e-10 +order = 11 + +[quantities.ep_real] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_plasma_energies" +type = "complex_vector" +extract = "real_first" +label = "plasma energy Re(ep[1])" +noise_threshold = 1e-10 +order = 12 + +[quantities.ev_real] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_vacuum_energies" +type = "complex_vector" +extract = "real_first" +label = "vacuum energy Re(ev[1])" +noise_threshold = 1e-10 +order = 13 + +[quantities.vacuum_eigenvalue] +h5path = "ForceFreeStates/FreeBoundaryStability/vacuum_eigenvalue" +type = "real_scalar" +extract = "value" +label = "vacuum matrix min eigenvalue" +noise_threshold = 1e-10 +order = 14 + +# Energies — full eigenvalue arrays +[quantities.ep_all] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_plasma_energies" +type = "complex_vector" +extract = "all_complex" +label = "plasma energy (all)" +noise_threshold = 1e-10 +order = 20 + +[quantities.ev_all] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_vacuum_energies" +type = "complex_vector" +extract = "all_complex" +label = "vacuum energy (all)" +noise_threshold = 1e-10 +order = 21 + +[quantities.et_all] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies" +type = "complex_vector" +extract = "all_complex" +label = "total energy (all)" +noise_threshold = 1e-10 +order = 22 + +# Integration +[quantities.nstep] +h5path = "ForceFreeStates/Solutions/ForwardIntegration/nstep" +type = "int_scalar" +extract = "value" +label = "ODE steps (saved)" +noise_threshold = 0 +order = 30 + +[quantities.nstep_total] +h5path = "ForceFreeStates/Solutions/ForwardIntegration/nstep_total" +type = "int_scalar" +extract = "value" +label = "ODE steps (total)" +noise_threshold = 0 +order = 31 + +# Equilibrium +[quantities.q0] +h5path = "Equilibrium/q_axis" +type = "real_scalar" +extract = "value" +label = "q0" +noise_threshold = 1e-12 +order = 40 + +[quantities.q95] +h5path = "Equilibrium/q_95" +type = "real_scalar" +extract = "value" +label = "q95" +noise_threshold = 1e-12 +order = 41 + +[quantities.betat] +h5path = "Equilibrium/beta_t" +type = "real_scalar" +extract = "value" +label = "beta_t" +noise_threshold = 1e-12 +order = 42 + +[quantities.betan] +h5path = "Equilibrium/beta_N" +type = "real_scalar" +extract = "value" +label = "beta_n" +noise_threshold = 1e-12 +order = 43 + +# Singular surfaces +[quantities.msing] +h5path = "SingularSurfaces/rational_count" +type = "int_scalar" +extract = "value" +label = "# singular surfaces" +noise_threshold = 0 +order = 50 + +[quantities.sing_psi] +h5path = "SingularSurfaces/rational_psi" +type = "real_vector" +extract = "all_real" +label = "singular psi locations" +noise_threshold = 1e-8 +order = 51 + +[quantities.sing_q] +h5path = "SingularSurfaces/rational_q" +type = "real_vector" +extract = "all_real" +label = "singular q values" +noise_threshold = 1e-8 +order = 52 + +# Mode numbers +[quantities.mpert] +h5path = "Info/mpert" +type = "int_scalar" +extract = "value" +label = "mpert" +noise_threshold = 0 +order = 60 + +[quantities.npert] +h5path = "Info/npert" +type = "int_scalar" +extract = "value" +label = "npert" +noise_threshold = 0 +order = 61 + +# Profiles (checksums for fast comparison of large arrays) +[quantities.q_profile] +h5path = "Equilibrium/Profiles/q" +type = "real_vector" +extract = "checksum" +label = "q profile (checksum)" +noise_threshold = 0 +order = 70 + +[quantities.pressure_profile] +h5path = "Equilibrium/Profiles/mu0p" +type = "real_vector" +extract = "checksum" +label = "pressure profile (checksum)" +noise_threshold = 0 +order = 71 + +# Perturbed equilibrium at the resonant surfaces, with no wall +[quantities.resonant_field] +h5path = "PerturbedEquilibrium/SingularCoupling/resonant_area_weighted_field" +type = "complex_vector" +extract = "all_complex" +label = "resonant area-weighted field" +noise_threshold = 1e-10 +order = 70 + +[quantities.island_half_width] +h5path = "PerturbedEquilibrium/SingularCoupling/island_half_width" +type = "real_vector" +extract = "all_real" +label = "island half-width" +noise_threshold = 1e-10 +order = 71 + +[quantities.chirikov] +h5path = "PerturbedEquilibrium/SingularCoupling/chirikov_parameter" +type = "real_vector" +extract = "all_real" +label = "Chirikov parameter" +noise_threshold = 1e-10 +order = 72 + +[quantities.toroidal_torque] +h5path = "PerturbedEquilibrium/Energies/toroidal_torque" +type = "real_scalar" +extract = "value" +label = "toroidal torque [N·m]" +noise_threshold = 1e-10 +order = 73 + +# Runtime (special: not from H5) +[quantities.runtime] +h5path = "" +type = "runtime" +extract = "value" +label = "Runtime (s)" +noise_threshold = 0.0 +order = 999 diff --git a/src/Deck.jl b/src/Deck.jl index cc68e19c9..af593faa0 100644 --- a/src/Deck.jl +++ b/src/Deck.jl @@ -446,6 +446,8 @@ function parse_deck(inputs::AbstractDict, path::AbstractString) end force_balance = has_ffs ? (; keywords=force_balance_keywords(inputs, path), algorithm=force_balance_algorithm(inputs)) : nothing kinetic_source = has_ffs && get(inputs["ForceFreeStates"], "kinetic", nothing) == "KineticForces" + haskey(inputs, "PerturbedEquilibrium") && haskey(inputs, "Wall") && + error("[PerturbedEquilibrium] needs a run without a wall (its surface inductance is free-space): remove the [Wall] section") haskey(inputs, "ForcingTerms") && !haskey(inputs, "PerturbedEquilibrium") && error("[ForcingTerms] is set but there is no [PerturbedEquilibrium] section to apply it") haskey(inputs, "KineticForces") && !haskey(inputs, "PerturbedEquilibrium") && !kinetic_source && diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl index d8493440a..b930fabcf 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl @@ -52,10 +52,14 @@ 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. +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. """ function plasma_response(ffs::ForceFreeStatesResult; compute_singular_coupling::Bool=true, compute_torque_response::Bool=false, verbose::Bool=false)::PlasmaResponse + ffs.wall_settings === nothing || + error("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 " * diff --git a/test/runtests_analysis_plots.jl b/test/runtests_analysis_plots.jl index 4a237bbb4..d17239ec6 100644 --- a/test/runtests_analysis_plots.jl +++ b/test/runtests_analysis_plots.jl @@ -10,7 +10,7 @@ using TOML A = GPEC.Analysis # Copy a coarse fixture deck into `dir` with HDF5 output on and the given extra sections; return its gpec.h5. - function _run_fixture(dir, fixture; sections=Dict{String,Any}()) + function _run_fixture(dir, fixture; sections=Dict{String,Any}(), ffs=Dict{String,Any}()) template = joinpath(@__DIR__, "test_data", fixture) for name in readdir(template) cp(joinpath(template, name), joinpath(dir, name)) @@ -18,6 +18,7 @@ using TOML toml_path = joinpath(dir, "gpec.toml") inputs = TOML.parsefile(toml_path) merge!(inputs, sections) + merge!(inputs["ForceFreeStates"], ffs) open(io -> TOML.print(io, inputs), toml_path, "w") GPEC.main([dir]) return joinpath(dir, "gpec.h5") @@ -27,9 +28,9 @@ using TOML @testset "ideal run with coil forcing and the perturbed-equilibrium response" begin mktempdir() do dir - h5 = _run_fixture(dir, "regression_solovev_ideal_example"; sections=Dict{String,Any}( + h5 = _run_fixture(dir, "regression_tj_ideal_example"; ffs=Dict{String,Any}("dW_edge_scan_start" => 0.95), sections=Dict{String,Any}( "ForcingTerms" => Dict{String,Any}("forcing_data_format" => "coil", "mtheta_coil" => 240, "nzeta_coil" => 32, - "coil_set" => [Dict{String,Any}("name" => "hoop", "source" => "pf_hoop", "radius" => 1.5, "height" => 0.4, + "coil_set" => [Dict{String,Any}("name" => "hoop", "source" => "pf_hoop", "radius" => 12.0, "height" => 1.2, "currents" => [2.0e3], "tiltx" => [3.0])]), "PerturbedEquilibrium" => Dict{String,Any}("compute_singular_coupling" => true))) diff --git a/test/runtests_deck_roundtrip.jl b/test/runtests_deck_roundtrip.jl index f84abc3c8..d0df7e00a 100644 --- a/test/runtests_deck_roundtrip.jl +++ b/test/runtests_deck_roundtrip.jl @@ -43,7 +43,7 @@ end recorded(path) = h5open(f -> TOML.parse(read(f["Input/gpec_toml_raw"])), path) stages(deck, dir) = Base.structdiff(GPEC.parse_deck(deck, dir), (; run=nothing)) mktempdir() do dir - example = joinpath(@__DIR__, "..", "examples", "Solovev_ideal_example") + example = joinpath(@__DIR__, "..", "examples", "TJ_ideal_example") foreach(name -> cp(joinpath(example, name), joinpath(dir, name)), readdir(example)) inputs = TOML.parsefile(joinpath(dir, "gpec.toml")) inputs["Run"]["verbose"] = false diff --git a/test/runtests_dominant_coupling.jl b/test/runtests_dominant_coupling.jl index db95ea471..9a4b0fad1 100644 --- a/test/runtests_dominant_coupling.jl +++ b/test/runtests_dominant_coupling.jl @@ -61,8 +61,8 @@ include("h5_metadata_check.jl") end @testset "perturbed-equilibrium run: summary, ResonantCoupling, and normalization" begin - template = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") - forcing = joinpath(@__DIR__, "..", "examples", "Solovev_ideal_example", "forcing.dat") + template = joinpath(@__DIR__, "test_data", "regression_tj_ideal_example") + forcing = joinpath(@__DIR__, "..", "examples", "TJ_ideal_example", "forcing.dat") mktempdir() do dir for name in readdir(template) @@ -88,7 +88,7 @@ include("h5_metadata_check.jl") @test !isempty(resp.dominant_singular_values) rank = length(resp.dominant_singular_values) rows = resp.dominant_rational_index - @test rows == 1:length(resp.rational_psi) + @test rows == findall(ψ -> PE.CORE_PSI_LOW <= ψ <= PE.CORE_PSI_HIGH, resp.rational_psi) @test issorted(resp.dominant_singular_values; rev=true) V, σ, U = resp.dominant_right_singular_vectors, resp.dominant_singular_values, resp.dominant_left_singular_vectors @test size(V) == (length(pe.forcing_b_rootarea), rank) @@ -120,7 +120,7 @@ include("h5_metadata_check.jl") # Windowing is a post-hoc choice on the same object. if length(rc.rational_psi) >= 2 lo = rc.rational_psi[2] - dom_w = PE.dominant_coupling(rc; psi_low=lo) + dom_w = PE.dominant_coupling(rc; psi_low=lo, psi_high=1.0) @test dom_w.rational_index == findall(>=(lo), rc.rational_psi) @test dom_w.singular_values ≈ svd(rc.C[dom_w.rational_index, :]).S end diff --git a/test/runtests_error_fields.jl b/test/runtests_error_fields.jl index c3d37e02e..8f2f1769e 100644 --- a/test/runtests_error_fields.jl +++ b/test/runtests_error_fields.jl @@ -60,8 +60,8 @@ include("h5_metadata_check.jl") @test lower.locking_probability_percent > straddling.locking_probability_percent end - @testset "coil-forced Solovev run: identities, output, and post-hoc entry point" begin - template = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") + @testset "coil-forced TJ run: identities, output, and post-hoc entry point" begin + template = joinpath(@__DIR__, "test_data", "regression_tj_ideal_example") mktempdir() do dir for name in readdir(template) @@ -69,22 +69,22 @@ include("h5_metadata_check.jl") end toml_path = joinpath(dir, "gpec.toml") inputs = TOML.parsefile(toml_path) - # Three PF hoops outside the plasma (Solovev R ≈ 0.67–1.33 m): one tilted so the run has + # Three PF hoops outside the plasma (R ≈ 9–11 m): one tilted so the run has # a nominal n=1 forcing, one axisymmetric so its rigid-motion response is known exactly, # and one tilted the other way as the correction array, which is swept but is not an # error-field source. inputs["ForcingTerms"] = Dict{String,Any}( "forcing_data_format" => "coil", "mtheta_coil" => 240, "nzeta_coil" => 32, "coil_set" => [ - Dict{String,Any}("name" => "hoop_tilted", "source" => "pf_hoop", "radius" => 1.5, "height" => 0.4, "currents" => [2.0e3], "tiltx" => [3.0]), - Dict{String,Any}("name" => "hoop_axi", "source" => "pf_hoop", "radius" => 1.5, "height" => -0.4, "currents" => [2.0e3]), - Dict{String,Any}("name" => "hoop_efc", "source" => "pf_hoop", "radius" => 1.5, "height" => 0.0, "currents" => [2.0e3], "tilty" => [2.0]) + Dict{String,Any}("name" => "hoop_tilted", "source" => "pf_hoop", "radius" => 12.0, "height" => 1.2, "currents" => [2.0e3], "tiltx" => [3.0]), + Dict{String,Any}("name" => "hoop_axi", "source" => "pf_hoop", "radius" => 12.0, "height" => -1.2, "currents" => [2.0e3]), + Dict{String,Any}("name" => "hoop_efc", "source" => "pf_hoop", "radius" => 12.0, "height" => 0.0, "currents" => [2.0e3], "tilty" => [2.0]) ]) inputs["PerturbedEquilibrium"] = Dict{String,Any}( "compute_singular_coupling" => true) cp(joinpath(@__DIR__, "test_data", "ErrorFields", "tolerances_two_hoops.toml"), joinpath(dir, "tolerances.toml")) - # Kinetic profiles consistent with this Solovev pressure, for the NTV torque couplings. - cp(joinpath(@__DIR__, "..", "examples", "Solovev_kinetic_NTV_example", "kinetic.dat"), joinpath(dir, "kinetic.dat")) + # Kinetic profiles consistent with this pressure, for the NTV torque couplings. + cp(joinpath(@__DIR__, "..", "examples", "TJ_kinetic_NTV_example", "kinetic.dat"), joinpath(dir, "kinetic.dat")) inputs["Equilibrium"]["kinetic_file"] = "kinetic.dat" # Exercise the kinetic pipeline at smoke-test precision, retaining all bounce harmonics. inputs["KineticForces"] = Dict{String,Any}( @@ -163,7 +163,7 @@ include("h5_metadata_check.jl") @test_throws ArgumentError EF.sensitivity_table(sens, shifted) # The bare-matrix construction carries no basis, so it stays permitted. @test EF.sensitivity_table(sens, PE.dominant_coupling(rc.C, rc.rational_psi)) isa EF.SensitivityTable - windowed = EF.sensitivity_table(sens, PE.dominant_coupling(rc; psi_low=rc.rational_psi[end])) + windowed = EF.sensitivity_table(sens, PE.dominant_coupling(rc; psi_low=rc.rational_psi[1])) @test length(windowed.delta_as_designed) == 3 # The correction array is swept with the others but is not an error-field source: the run # excludes the NTV efc_coils from the error field, and the table helper drops them. @@ -204,7 +204,7 @@ include("h5_metadata_check.jl") @test mc_file.abs_delta_pdf == mc.abs_delta_pdf && mc_file.abs_delta_bin_edges == mc.abs_delta_bin_edges && mc_file.nsample == 512 rerun = EF.run_monte_carlo(h5path; nsample=512, nbatch=2, seed=5, nbins=16) @test rerun.abs_delta_pdf == mc.abs_delta_pdf - windowed_mc = EF.run_monte_carlo(h5path; psi_low=rc.rational_psi[end], nsample=256, nbatch=1, seed=5, nbins=16) + windowed_mc = EF.run_monte_carlo(h5path; psi_low=rc.rational_psi[1], nsample=256, nbatch=1, seed=5, nbins=16) @test windowed_mc.abs_delta_total_as_designed ≈ abs(sum(EF.without_coils(windowed, ["hoop_efc"]).delta_as_designed)) # Locking risk and tolerance scan: written, bounded, and reproducible from the file. @@ -640,7 +640,7 @@ include("h5_metadata_check.jl") end @testset "[Run] output_file holds every stage and the runtimes" begin - template = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") + template = joinpath(@__DIR__, "test_data", "regression_tj_ideal_example") mktempdir() do dir for name in readdir(template) cp(joinpath(template, name), joinpath(dir, name)) diff --git a/test/runtests_frozen_equilibrium.jl b/test/runtests_frozen_equilibrium.jl index 96aa9dc14..e6814cb16 100644 --- a/test/runtests_frozen_equilibrium.jl +++ b/test/runtests_frozen_equilibrium.jl @@ -117,7 +117,7 @@ using GeneralizedPerturbedEquilibrium @test analytic_snapshot.profiles.q_spline.y == profile_values[4] # The stored deck's kinetic_file resolves against the h5's directory and is re-read; a missing file warns and leaves none. kinetic_path = joinpath(dir, "kinetic.dat") - cp(joinpath(@__DIR__, "..", "examples", "Solovev_kinetic_NTV_example", "kinetic.dat"), kinetic_path) + cp(joinpath(@__DIR__, "..", "examples", "TJ_kinetic_NTV_example", "kinetic.dat"), kinetic_path) analytic_inputs["Equilibrium"]["kinetic_file"] = "kinetic.dat" h5open(path, "r+") do f delete_object(f, "Input/gpec_toml_raw") diff --git a/test/runtests_fullruns.jl b/test/runtests_fullruns.jl index 67c90396e..245b7cdf7 100644 --- a/test/runtests_fullruns.jl +++ b/test/runtests_fullruns.jl @@ -45,7 +45,7 @@ using TOML h5open(test_kinetic_run(dir, "regression_solovev_ideal_example"), "r") do h5 # Smoke test: this nerfed deck (mpsi=16, delta_m=0) is for "does it run / no NaN", # not numeric regression — the mode-converged physical value is pinned in the - # regression harness (examples/Solovev_kinetic_calculated_example). Assert only + # regression harness (examples/TJ_kinetic_calculated_example). Assert only # nerfed-grid-robust facts: finite and positive (Solovev is stable). et = read(h5["ForceFreeStates/FreeBoundaryStability/eigenmode_energies"]) @test isfinite(real(et[1])) @@ -84,7 +84,7 @@ using TOML # KF→FFS kinetic-MHD path end-to-end. NO numeric value is pinned here — the prior # imag(et[1]) ≈ -0.711 rtol=0.08 pin was platform-fragile (failed on macOS aarch64 # at -0.856; issue #273). Physical regression is tracked in the harness on the - # mode-converged deck (examples/Solovev_kinetic_calculated_example, mpert=32). + # mode-converged deck (examples/TJ_kinetic_calculated_example, mpert=32). # Assert only nerfed-grid-robust physics: finite, positive total energy, # negative kinetic damping (imag sign is robust across configs; only its # magnitude is FP-sensitive). diff --git a/test/runtests_rerun_from_h5.jl b/test/runtests_rerun_from_h5.jl index e09df1e84..952a01005 100644 --- a/test/runtests_rerun_from_h5.jl +++ b/test/runtests_rerun_from_h5.jl @@ -2,7 +2,7 @@ using HDF5 using TOML # Tests for the gpec.h5 snapshot + replay feature, kept deliberately light: one full-pipeline -# source run + replay (the analytic Solovev path), plus pipeline-free ingest round-trip checks +# source run + replay (the analytic TJ path), plus pipeline-free ingest round-trip checks # for the direct (EFIT) and inverse (CHEASE) equilibrium kinds and the override flags. # Collect every leaf dataset path under an open HDF5 file, skipping the groups/paths that @@ -50,19 +50,19 @@ function _rerun_dataset_mismatches(source_h5, rerun_h5) return n_compared, n_mismatched end -# Analytic (Solovev) equilibrium with ASCII file forcing: the single end-to-end source + replay. +# Analytic (TJ) equilibrium with ASCII file forcing: the single end-to-end source + replay. # The forcing file is deleted before replay and the replay runs from a temp directory, proving # the snapshot is self-contained (no g-file or forcing data needed), and reproduces the source # bit-for-bit across every output dataset, the normalized perturbed equilibrium included. -@testset "Rerun from gpec.h5 (analytic Solovev)" begin - template_dir = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") +@testset "Rerun from gpec.h5 (analytic TJ)" begin + template_dir = joinpath(@__DIR__, "test_data", "regression_tj_ideal_example") mktempdir() do snapshot_dir for name in readdir(template_dir) cp(joinpath(template_dir, name), joinpath(snapshot_dir, name)) end forcing_dat = joinpath(snapshot_dir, "forcing.dat") - cp(joinpath(@__DIR__, "..", "examples", "Solovev_ideal_example", "forcing.dat"), forcing_dat) + cp(joinpath(@__DIR__, "..", "examples", "TJ_ideal_example", "forcing.dat"), forcing_dat) toml_path = joinpath(snapshot_dir, "gpec.toml") inputs = TOML.parsefile(toml_path) inputs["ForcingTerms"] = Dict{String,Any}("forcing_data_file" => "forcing.dat", "forcing_data_format" => "ascii") @@ -136,7 +136,7 @@ end end @testset "--coil-source coils errors when source has no coil snapshot" begin - # This Solovev source run used ASCII forcing, so there is no + # This source run used ASCII forcing, so there is no # Input/RawInputs/Coils group to replay from. mktempdir() do replay_dir @test_throws ErrorException GeneralizedPerturbedEquilibrium.main([ @@ -153,7 +153,7 @@ end # the source run and the replay: the default (recompute-from-TOML) path then # fails, while `--coil-source coils` succeeds from the stored geometry. @testset "Rerun from gpec.h5: coil geometry snapshot" begin - template_dir = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") + template_dir = joinpath(@__DIR__, "test_data", "regression_tj_ideal_example") src_dat = joinpath(@__DIR__, "..", "src", "ForcingTerms", "coil_geometries", "d3d_il.dat") read_resfield(path) = diff --git a/test/runtests_result_struct.jl b/test/runtests_result_struct.jl index 97fe0724b..9c1597e56 100644 --- a/test/runtests_result_struct.jl +++ b/test/runtests_result_struct.jl @@ -29,7 +29,7 @@ include("h5_metadata_check.jl") merge!(inputs["Equilibrium"], Dict("mpsi" => 32, "mtheta" => 64)) merge!(inputs["TJAnalytic"], Dict("ma" => 32, "mtheta" => 32)) merge!(inputs["ForceFreeStates"], Dict("delta_mlow" => 2, "delta_mhigh" => 2, "mtheta_vacuum" => 64)) - inputs["Galerkin"]["nx"] = 32 + haskey(inputs, "Galerkin") && (inputs["Galerkin"]["nx"] = 32) end open(io -> TOML.print(io, inputs), toml_path, "w") return dir @@ -127,7 +127,7 @@ include("h5_metadata_check.jl") @testset "PerturbedEquilibrium warn-skips on a riccati result" begin mktempdir() do dir - _stage_deck(dir, template; + _stage_deck(dir, joinpath(@__DIR__, "test_data", "regression_tj_ideal_example"); ffs_overrides=Dict{String,Any}("integrator" => "Riccati"), extra_sections=Dict{String,Any}( "ForcingTerms" => Dict{String,Any}( @@ -135,7 +135,7 @@ include("h5_metadata_check.jl") "forcing_data_format" => "ascii"), "PerturbedEquilibrium" => Dict{String,Any}( "compute_singular_coupling" => true))) - cp(joinpath(@__DIR__, "..", "examples", "Solovev_ideal_example", "forcing.dat"), + cp(joinpath(@__DIR__, "..", "examples", "TJ_ideal_example", "forcing.dat"), joinpath(dir, "forcing.dat")) local ret diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 7b56441b1..6a1b63015 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -2,12 +2,12 @@ using TOML using HDF5 # The scripting API: `PlasmaEquilibrium` / `solve(eq, alg)` / `perturbed_equilibrium` must run the -# same pipeline the TOML driver does. Every assertion is anchored on the coarse Solovev fixture +# same pipeline the TOML driver does. Every assertion is anchored on the coarse TJ fixture # deck (mpsi=16, mtheta_vacuum=64, delta_m=0), so an API run is compared against `main` on the same deck. @testset "solve API" begin GPEC = GeneralizedPerturbedEquilibrium FFS = GPEC.ForceFreeStates - template = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") + template = joinpath(@__DIR__, "test_data", "regression_tj_ideal_example") deck = TOML.parsefile(joinpath(template, "gpec.toml")) # The fixture equilibrium and wall, built exactly as `build_inputs_from_toml` builds them. @@ -132,6 +132,8 @@ using HDF5 @test api.free_boundary.et[1] ≈ reference.free_boundary.et[1] rtol = 1e-10 @test api.diagnostics.nzero == reference.diagnostics.nzero @test api.kinetic.kmsing == reference.kinetic.kmsing + # This fixture solves with a wall, and the free-space surface inductance cannot pair with its Λ. + @test_throws "without a wall" plasma_response(api) # The KineticForces matrices evaluated once solve identically when reused, and only on the problem they were made for. alg = GPEC.force_balance_algorithm(kdeck) @@ -192,7 +194,7 @@ using HDF5 end @testset "perturbed_equilibrium round-trips a forward result" begin - forcing = joinpath(@__DIR__, "..", "examples", "Solovev_ideal_example", "forcing.dat") + forcing = joinpath(@__DIR__, "..", "examples", "TJ_ideal_example", "forcing.dat") pe_section = Dict{String,Any}( "compute_singular_coupling" => true) pe_kw = (; verbose=false) @@ -448,7 +450,7 @@ using HDF5 pe_kwargs = (; compute_singular_coupling=true) mktempdir() do dir modes = GPEC.ForcingTerms.ForcingMode[] - GPEC.ForcingTerms.load_forcing_data!(modes, joinpath(@__DIR__, "..", "examples", "Solovev_ideal_example", "forcing.dat"), "ascii", false) + GPEC.ForcingTerms.load_forcing_data!(modes, joinpath(@__DIR__, "..", "examples", "TJ_ideal_example", "forcing.dat"), "ascii", false) h5open(f -> GPEC.ForcingTerms.save_forcing_to_h5(modes, f), joinpath(dir, "forcing.h5"), "w") _stage_deck(dir; extra_sections=Dict{String,Any}( @@ -523,6 +525,14 @@ using HDF5 deck_dir = _stage_deck(dir; extra_sections=Dict{String,Any}("Match" => Dict{String,Any}("ideal" => true))) @test_throws "[Match] needs Δ′" GPEC.parse_deck(TOML.parsefile(joinpath(deck_dir, "gpec.toml")), deck_dir) end + + # A perturbed equilibrium with a wall stops at parse time. + mktempdir() do dir + walled = Dict{String,Any}("Wall" => Dict{String,Any}("shape" => "conformal", "standoff_fraction" => 0.2), + "PerturbedEquilibrium" => Dict{String,Any}("compute_singular_coupling" => true)) + deck_dir = _stage_deck(dir; extra_sections=walled) + @test_throws "without a wall" GPEC.parse_deck(TOML.parsefile(joinpath(deck_dir, "gpec.toml")), deck_dir) + end end @testset "deck kinetic source names" begin @@ -537,7 +547,7 @@ using HDF5 @testset "kinetic torque" begin # The deck's [KineticForces] stage is kinetic_torque on its perturbed equilibrium, saved by its save method. - ntv_dir = joinpath(@__DIR__, "..", "examples", "Solovev_kinetic_NTV_example") + ntv_dir = joinpath(@__DIR__, "..", "examples", "TJ_kinetic_NTV_example") mktempdir() do dir for name in readdir(ntv_dir) cp(joinpath(ntv_dir, name), joinpath(dir, name)) diff --git a/test/runtests_torque_response.jl b/test/runtests_torque_response.jl index fd3f3c14b..8ceb87036 100644 --- a/test/runtests_torque_response.jl +++ b/test/runtests_torque_response.jl @@ -2,9 +2,9 @@ using Test using HDF5 using LinearAlgebra -# ψ-resolved torque response matrices of a kinetic forward solve: the Solovev calculated-kinetic +# ψ-resolved torque response matrices of a kinetic forward solve: the TJ calculated-kinetic # fixture driven by two analytic window-pane coil sets, so T_xe, T_source and the run's own -# forcing are all exercised on one (cheap, mpsi=16) solve. +# forcing are all exercised on one (cheap, mpsi=32) solve. const _GPE_TR = GeneralizedPerturbedEquilibrium const _PE_TR = _GPE_TR.PerturbedEquilibrium @@ -22,7 +22,7 @@ nzeta_coil = 32 # Toroidal grid (0 = auto = 32·n) name = "pane_upper" # Analytic window-pane array above the midplane source = "window_pane" # Rectangular picture-frame coils placed by explicit corners n_coils = 6 # Number of frames around the torus -rz_corners = [[1.5, 0.05], [1.5, 0.45]] # Frame corners (R, Z) [m] +rz_corners = [[12.0, 0.1], [12.0, 1.5]] # Frame corners (R, Z) [m] gap_fraction = 0.15 # Toroidal gap between neighbouring frames as a fraction of the pitch currents = [1000.0, 500.0, -500.0, -1000.0, -500.0, 500.0] # n=1 cosine phasing, 1 kA peak [A] @@ -30,7 +30,7 @@ currents = [1000.0, 500.0, -500.0, -1000.0, -500.0, 500.0] # n=1 cosine phasing name = "pane_lower" # Analytic window-pane array below the midplane source = "window_pane" # Rectangular picture-frame coils placed by explicit corners n_coils = 6 # Number of frames around the torus -rz_corners = [[1.5, -0.45], [1.5, -0.05]] # Frame corners (R, Z) [m] +rz_corners = [[12.0, -1.5], [12.0, -0.1]] # Frame corners (R, Z) [m] gap_fraction = 0.15 # Toroidal gap between neighbouring frames as a fraction of the pitch currents = [1000.0, 0.0, -1000.0, -1000.0, 0.0, 1000.0] # n=1 phasing shifted by 60 degrees, 1 kA peak [A] @@ -40,7 +40,7 @@ compute_torque_response = true # Build the ψ-resolved torque response """ @testset "Torque response matrices" begin - template = joinpath(@__DIR__, "test_data", "regression_solovev_kinetic_calculated") + template = joinpath(@__DIR__, "test_data", "regression_tj_kinetic_calculated") mktempdir() do dir for name in readdir(template) diff --git a/test/test_data/regression_tj_ideal_example/gpec.toml b/test/test_data/regression_tj_ideal_example/gpec.toml new file mode 100644 index 000000000..4d75cae03 --- /dev/null +++ b/test/test_data/regression_tj_ideal_example/gpec.toml @@ -0,0 +1,52 @@ +# Cheap TJ-analytic equilibrium for the test suite: stable at n = 1 with no wall, so the tests can add +# [PerturbedEquilibrium] and [ForcingTerms] sections. A coarse copy of examples/TJ_ideal_example. + +[Equilibrium] +jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) +grid_type = "core_edge_packed" # Radial grid packing type +psilow = 0.01 # Lower limit of normalized poloidal flux +psihigh = 0.99 # Upper limit of normalized poloidal flux; et[1] depends on it (+0.502 at 0.995), so it is pinned +mpsi = 32 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) +mtheta = 256 # Number of poloidal grid points + +# Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the LocalStabilityProblem settings, numerics in [Ballooning] +[LocalStability] +theta_k = 0.0 # Ballooning angle θ_k, used by both the profile scan and the α-boundary scan +alpha_boundary_scan = true # Run the first ballooning α-boundary scan (a root find per ψ surface) + +# Numerical settings of the ballooning scans +[Ballooning] +ode_tol = 1e-8 # Relative and absolute tolerance of the ballooning ODE integration behind Δ' +n_scan = 24 # Resolution of the α scan for the first boundary +max_alpha_scale = 8.0 # Upper end of the α search, as a multiple of the experimental α +alpha_tol = 1e-3 # Bisection tolerance of the critical α, in the same multiples +psi_window = [0.1, 0.99] # ψ_N range of the α scan; surfaces outside it are NaN + +[ForceFreeStates] +vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes +qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) +qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) +nn = 1 # Toroidal mode number, or consecutive numbers as a list +delta_mlow = 4 # Expands lower bound of Fourier harmonics +delta_mhigh = 4 # Expands upper bound of Fourier harmonics +mtheta_vacuum = 64 # Number of points used in splines over poloidal angle at the plasma-vacuum interface +singfac_min = 1e-6 # Fractional distance from rational q at which ideal jump enforced +integrator = "Forward" # Integration formalism, whose numerical settings are its own section: "Forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "Riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "Galerkin" (outer-region singular Galerkin Δ′) +truncate_at_rational_offset = false # Truncate integration at (last_rational_q + rational_offset_fraction)/n; false for analytic equilibria + +# Numerical settings of the forward integrator +[Forward] +eulerlagrange_tolerance = 1e-9 # Relative tolerance for ODE integration of Euler-Lagrange equations + +# TJ-analytic equilibrium parameters (cf. R. Fitzpatrick, TJ code, https://github.com/rfitzp/TJ); this section selects the +# equilibrium. Keep p0 at or below 0.2: et[1] stays near +0.42 up to there and is not monotonic above it. +[TJAnalytic] +r0 = 10.0 # Major radius R₀ [m] (centerline radius of the magnetic axis) +a = 1.0 # Minor radius a [m] (plasma half-width at the midplane; here ε = 0.1) +q0 = 1.6 # On-axis safety factor q(r=0) +qa = 3.2 # Edge safety factor q(r=a) +p0 = 0.1 # Normalized on-axis pressure +pressure_peaking = 2.0 # Pressure-profile peaking exponent in p₂(r) = p0·(1−r²)^μ +B0 = 1.0 # On-axis toroidal field [T] +ma = 128 # TJ-analytic internal radial grid resolution (shape-ODE nodes) +mtheta = 128 # TJ-analytic internal poloidal grid resolution (θ-spline nodes) diff --git a/test/test_data/regression_tj_kinetic_calculated/gpec.toml b/test/test_data/regression_tj_kinetic_calculated/gpec.toml new file mode 100644 index 000000000..3fe8fdfad --- /dev/null +++ b/test/test_data/regression_tj_kinetic_calculated/gpec.toml @@ -0,0 +1,57 @@ +# Cheap TJ-analytic equilibrium with calculated kinetic matrices for the test suite: stable at n = 1 +# with no wall. A coarse copy of examples/TJ_kinetic_calculated_example. + +[Equilibrium] +kinetic_file = "kinetic.dat" # Kinetic-profile file (HDF5 GPEC kinetic schema or 6-column ASCII), read once into the equilibrium +jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) +grid_type = "core_edge_packed" # Radial grid packing type +psilow = 0.01 # Lower limit of normalized poloidal flux +psihigh = 0.99 # Upper limit of normalized poloidal flux; et[1] depends on it (+0.502 at 0.995), so it is pinned +mpsi = 32 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) +mtheta = 256 # Number of poloidal grid points + +# Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the LocalStabilityProblem settings, numerics in [Ballooning] +[LocalStability] +theta_k = 0.0 # Ballooning angle θ_k, used by both the profile scan and the α-boundary scan +alpha_boundary_scan = true # Run the first ballooning α-boundary scan (a root find per ψ surface) + +# Numerical settings of the ballooning scans +[Ballooning] +ode_tol = 1e-8 # Relative and absolute tolerance of the ballooning ODE integration behind Δ' +n_scan = 24 # Resolution of the α scan for the first boundary +max_alpha_scale = 8.0 # Upper end of the α search, as a multiple of the experimental α +alpha_tol = 1e-3 # Bisection tolerance of the critical α, in the same multiples +psi_window = [0.1, 0.99] # ψ_N range of the α scan; surfaces outside it are NaN + +[ForceFreeStates] +vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes +qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) +qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) +nn = 1 # Toroidal mode number, or consecutive numbers as a list +delta_mlow = 4 # Expands lower bound of Fourier harmonics +delta_mhigh = 4 # Expands upper bound of Fourier harmonics +mtheta_vacuum = 64 # Number of points used in splines over poloidal angle at the plasma-vacuum interface +singfac_min = 1e-6 # Fractional distance from rational q at which ideal jump enforced +kinetic = "KineticForces" # Kinetic source: the drift-kinetic matrices of the [KineticForces] model (omit for an ideal run) +integrator = "Forward" # Integration formalism, whose numerical settings are its own section: "Forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "Riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "Galerkin" (outer-region singular Galerkin Δ′) +truncate_at_rational_offset = false # Truncate integration at (last_rational_q + rational_offset_fraction)/n; false for analytic equilibria + +# Numerical settings of the forward integrator +[Forward] +eulerlagrange_tolerance = 1e-9 # Relative tolerance for ODE integration of Euler-Lagrange equations + +# TJ-analytic equilibrium parameters (cf. R. Fitzpatrick, TJ code, https://github.com/rfitzp/TJ); this section selects the +# equilibrium. Keep p0 at or below 0.2: et[1] stays near +0.42 up to there and is not monotonic above it. +[TJAnalytic] +r0 = 10.0 # Major radius R₀ [m] (centerline radius of the magnetic axis) +a = 1.0 # Minor radius a [m] (plasma half-width at the midplane; here ε = 0.1) +q0 = 1.6 # On-axis safety factor q(r=0) +qa = 3.2 # Edge safety factor q(r=a) +p0 = 0.1 # Normalized on-axis pressure +pressure_peaking = 2.0 # Pressure-profile peaking exponent in p₂(r) = p0·(1−r²)^μ +B0 = 1.0 # On-axis toroidal field [T] +ma = 128 # TJ-analytic internal radial grid resolution (shape-ODE nodes) +mtheta = 128 # TJ-analytic internal poloidal grid resolution (θ-spline nodes) + +[KineticForces] +axis_validity_suppression = true # Suppress kinetic terms inside the near-axis orbit-width validity limit diff --git a/test/test_data/regression_tj_kinetic_calculated/kinetic.dat b/test/test_data/regression_tj_kinetic_calculated/kinetic.dat new file mode 100644 index 000000000..e3a393077 --- /dev/null +++ b/test/test_data/regression_tj_kinetic_calculated/kinetic.dat @@ -0,0 +1,45 @@ +# Kinetic profiles consistent with the TJ-analytic equilibrium's pressure (P₀ = 784 Pa). +# Constraint: P(psi_n) = n_e*T_e + n_i*T_i, with n_e = n_i and T_e = T_i = 20 + 180*(1 - psi_n) eV. +# Omega_E = 2e4*(1 - psi_n) rad/s. +# psi_n n_i n_e T_i[eV] T_e[eV] omega_E +0.000000 1.223086e+19 1.223086e+19 200.0000 200.0000 20000.0000 +0.025000 1.222817e+19 1.222817e+19 195.5000 195.5000 19500.0000 +0.050000 1.203522e+19 1.203522e+19 191.0000 191.0000 19000.0000 +0.075000 1.183797e+19 1.183797e+19 186.5000 186.5000 18500.0000 +0.100000 1.163625e+19 1.163625e+19 182.0000 182.0000 18000.0000 +0.125000 1.142988e+19 1.142988e+19 177.5000 177.5000 17500.0000 +0.150000 1.121868e+19 1.121868e+19 173.0000 173.0000 17000.0000 +0.175000 1.100245e+19 1.100245e+19 168.5000 168.5000 16500.0000 +0.200000 1.078097e+19 1.078097e+19 164.0000 164.0000 16000.0000 +0.225000 1.055403e+19 1.055403e+19 159.5000 159.5000 15500.0000 +0.250000 1.032142e+19 1.032142e+19 155.0000 155.0000 15000.0000 +0.275000 1.008285e+19 1.008285e+19 150.5000 150.5000 14500.0000 +0.300000 9.838102e+18 9.838102e+18 146.0000 146.0000 14000.0000 +0.325000 9.586866e+18 9.586866e+18 141.5000 141.5000 13500.0000 +0.350000 9.328821e+18 9.328821e+18 137.0000 137.0000 13000.0000 +0.375000 9.063714e+18 9.063714e+18 132.5000 132.5000 12500.0000 +0.400000 8.791260e+18 8.791260e+18 128.0000 128.0000 12000.0000 +0.425000 8.511117e+18 8.511117e+18 123.5000 123.5000 11500.0000 +0.450000 8.222884e+18 8.222884e+18 119.0000 119.0000 11000.0000 +0.475000 7.926125e+18 7.926125e+18 114.5000 114.5000 10500.0000 +0.500000 7.620476e+18 7.620476e+18 110.0000 110.0000 10000.0000 +0.525000 7.305553e+18 7.305553e+18 105.5000 105.5000 9500.0000 +0.550000 6.980924e+18 6.980924e+18 101.0000 101.0000 9000.0000 +0.575000 6.646101e+18 6.646101e+18 96.5000 96.5000 8500.0000 +0.600000 6.300621e+18 6.300621e+18 92.0000 92.0000 8000.0000 +0.625000 5.944099e+18 5.944099e+18 87.5000 87.5000 7500.0000 +0.650000 5.576149e+18 5.576149e+18 83.0000 83.0000 7000.0000 +0.675000 5.196380e+18 5.196380e+18 78.5000 78.5000 6500.0000 +0.700000 4.804420e+18 4.804420e+18 74.0000 74.0000 6000.0000 +0.725000 4.399955e+18 4.399955e+18 69.5000 69.5000 5500.0000 +0.750000 3.982956e+18 3.982956e+18 65.0000 65.0000 5000.0000 +0.775000 3.553750e+18 3.553750e+18 60.5000 60.5000 4500.0000 +0.800000 3.113088e+18 3.113088e+18 56.0000 56.0000 4000.0000 +0.825000 2.662383e+18 2.662383e+18 51.5000 51.5000 3500.0000 +0.850000 2.204114e+18 2.204114e+18 47.0000 47.0000 3000.0000 +0.875000 1.742718e+18 1.742718e+18 42.5000 42.5000 2500.0000 +0.900000 1.285949e+18 1.285949e+18 38.0000 38.0000 2000.0000 +0.925000 8.470999e+17 8.470999e+17 33.5000 33.5000 1500.0000 +0.950000 4.497124e+17 4.497124e+17 29.0000 29.0000 1000.0000 +0.975000 1.378556e+17 1.378556e+17 24.5000 24.5000 500.0000 +1.000000 2.762506e+16 2.762506e+16 20.0000 20.0000 0.0000