diff --git a/benchmarks/benchmark_against_fortran_run.jl b/benchmarks/benchmark_against_fortran_run.jl index 6116a8058..188f9e1c9 100644 --- a/benchmarks/benchmark_against_fortran_run.jl +++ b/benchmarks/benchmark_against_fortran_run.jl @@ -406,9 +406,6 @@ function write_gpec_toml_coil( println(io, "mtheta = $(p.mtheta_equil)") println(io, "etol = 1e-7") println(io) - println(io, "[Wall]") - println(io, "shape = \"nowall\"") - println(io) println(io, "[ForceFreeStates]") println(io, "local_stability_flag = true") println(io, "vac_flag = true") @@ -469,9 +466,6 @@ function write_gpec_toml_file( println(io, "mtheta = $(p.mtheta_equil)") println(io, "etol = 1e-7") println(io) - println(io, "[Wall]") - println(io, "shape = \"nowall\"") - println(io) println(io, "[ForceFreeStates]") println(io, "local_stability_flag = true") println(io, "vac_flag = true") diff --git a/benchmarks/benchmark_delta_prime_methods.jl b/benchmarks/benchmark_delta_prime_methods.jl index 9390c1eb1..48913d36b 100644 --- a/benchmarks/benchmark_delta_prime_methods.jl +++ b/benchmarks/benchmark_delta_prime_methods.jl @@ -29,8 +29,7 @@ function setup_and_run_solovev() (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex) equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config) - intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(; - (Symbol(k) => v for (k, v) in inputs["Wall"])...) + intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.parse_wall(get(inputs, "Wall", nothing)) FFS.sing_lim!(intr, ctrl, equil) intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; diff --git a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl index 7da279b91..021abee9b 100644 --- a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl +++ b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl @@ -35,7 +35,9 @@ const KF = GPE.KineticForces const Eq = GPE.Equilibrium const PE = GPE.PerturbedEquilibrium -"Fortran GPEC DIII-D ideal example run directory, from the environment (no on-disk default)." +""" +Fortran GPEC DIII-D ideal example run directory, from the environment (no on-disk default). +""" default_fortran_dir() = get(() -> error("Set GPEC_FORTRAN_DIIID, or pass the Fortran run directory as the first argument"), ENV, "GPEC_FORTRAN_DIIID") @@ -125,9 +127,6 @@ mtheta = 256 newq0 = 0 etol = 1e-7 -[Wall] -shape = "nowall" - [ForceFreeStates] local_stability_flag = true vac_flag = true diff --git a/benchmarks/benchmark_diiid_kinetic_stability.jl b/benchmarks/benchmark_diiid_kinetic_stability.jl index 9ae37b6b3..8f29510d6 100644 --- a/benchmarks/benchmark_diiid_kinetic_stability.jl +++ b/benchmarks/benchmark_diiid_kinetic_stability.jl @@ -36,7 +36,9 @@ using GeneralizedPerturbedEquilibrium const GPE = GeneralizedPerturbedEquilibrium const AnalysisFFS = GPE.Analysis.ForceFreeStates -"Fortran GPEC DIII-D kinetic example run directory, from the environment (no on-disk default)." +""" +Fortran GPEC DIII-D kinetic example run directory, from the environment (no on-disk default). +""" default_fortran_dir() = get(() -> error("Set GPEC_FORTRAN_DIIID_DCON, or pass the Fortran run directory as the first argument"), ENV, "GPEC_FORTRAN_DIIID_DCON") @@ -132,9 +134,6 @@ mtheta = 256 newq0 = 0 etol = 1e-7 -[Wall] -shape = "nowall" - [ForceFreeStates] local_stability_flag = true vac_flag = true diff --git a/benchmarks/benchmark_riccati_der.jl b/benchmarks/benchmark_riccati_der.jl index e9d88ef82..bcaaab2e3 100644 --- a/benchmarks/benchmark_riccati_der.jl +++ b/benchmarks/benchmark_riccati_der.jl @@ -30,8 +30,7 @@ function setup_solovev() (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex) equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config) - intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(; - (Symbol(k) => v for (k, v) in inputs["Wall"])...) + intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.parse_wall(get(inputs, "Wall", nothing)) FFS.sing_lim!(intr, ctrl, equil) intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; diff --git a/benchmarks/benchmark_threads.jl b/benchmarks/benchmark_threads.jl index 1b2912a70..11834a641 100644 --- a/benchmarks/benchmark_threads.jl +++ b/benchmarks/benchmark_threads.jl @@ -18,8 +18,7 @@ function run_ffs(ex; integrator) (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex) equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config) - intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(; - (Symbol(k) => v for (k, v) in inputs["Wall"])...) + intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.parse_wall(get(inputs, "Wall", nothing)) GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil) intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; diff --git a/docs/development/architecture.md b/docs/development/architecture.md index 37aad79b3..095f2bf73 100644 --- a/docs/development/architecture.md +++ b/docs/development/architecture.md @@ -136,7 +136,7 @@ Splines are provided by the external `FastInterpolations` package rather than by All GPEC modules are configured via a single TOML file with the following sections: - `[Equilibrium]` - Equilibrium solver settings -- `[Wall]` - Wall geometry and vacuum region +- `[Wall]` - Wall geometry (omitted for no wall) - `[ForceFreeStates]` - Stability analysis parameters - `[PerturbedEquilibrium]` - Perturbed equilibrium settings - `[ForcingTerms]` - External field specification @@ -199,7 +199,7 @@ The complete GPEC analysis pipeline: ### Vacuum - `VacuumInput` - Input parameters for vacuum calculations -- `WallShapeSettings` - Wall geometry configuration +- `WallSettings` - Wall shape configuration (`ConformalWall`, `DeeWall`, ...); no wall is `nothing` ### Stability - `SingType` - Singular surface data including: diff --git a/docs/development/toml-conventions.md b/docs/development/toml-conventions.md index 5ee2c3d5d..89f7ee144 100644 --- a/docs/development/toml-conventions.md +++ b/docs/development/toml-conventions.md @@ -10,7 +10,7 @@ modify them. **Do not invent a new convention** — match what the existing exam n-range, ideal vs kinetic). No ASCII-art separators or decorative dividers. 2. **Inline annotation on every variable line.** Every `key = value` gets a trailing `# description`. Source the wording from the matching config struct's docstring - (`EquilibriumConfig`, `WallShapeSettings`, `ForceFreeStatesControl`, `ForcingTermsControl`, + (`EquilibriumConfig`, the `WallSettings` types, `ForceFreeStatesControl`, `ForcingTermsControl`, `PerturbedEquilibriumControl`, `SolovevConfig`, `TJAnalyticConfig`, `KineticForcesControl`) and keep it to one terse line. **Use the same description for the same variable across all files** — `examples/DIIID-like_ideal_example/gpec.toml` is the canonical reference for the @@ -20,7 +20,8 @@ modify them. **Do not invent a new convention** — match what the existing exam the header plus an informative block comment per quantity group instead. 3. **Section comments only when informative.** No comment is required above a section. Keep a block comment only when it carries real information (e.g. the Solovev `[Wall]` note on why a - conformal wall is needed). Never add decorative section dividers. + conformal wall is needed). Never add decorative section dividers. No wall is spelled by omitting `[Wall]` + entirely; a `[Wall]` section always names a `shape` and only that shape's parameters. 4. **No Fortran references — annotations must be self-contained and Julia-centric.** The Julia code is the production workhorse; a new user must not need to know the legacy Fortran GPEC code to read an example. Headers explain the *scenario and its unique aspect*; inline comments diff --git a/docs/src/set_up.md b/docs/src/set_up.md index fcebbc0a3..f07cd9ac1 100644 --- a/docs/src/set_up.md +++ b/docs/src/set_up.md @@ -353,7 +353,7 @@ println("Beta-N: ", equil.params.betan) GPEC uses TOML configuration files (`gpec.toml`) with the following main sections: - **`[Equilibrium]`**: Equilibrium solver settings (input file, grid resolution, coordinate system, etc.) -- **`[Wall]`**: Wall geometry for vacuum calculations (shape, size, position) +- **`[Wall]`**: Wall geometry for vacuum calculations (shape, size, position); omit for no wall - **`[ForceFreeStates]`**: Stability analysis settings (mode numbers, tolerances, flags) - **`[PerturbedEquilibrium]`**: Plasma response settings (forcing data, output options) diff --git a/docs/src/stability.md b/docs/src/stability.md index c61ca716e..cfc2e2a1a 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -313,8 +313,7 @@ equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium( GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex)) intr = FFS.ForceFreeStatesInternal(; dir_path=ex) -intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(; - (Symbol(k) => v for (k, v) in inputs["Wall"])...) +intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.parse_wall(get(inputs, "Wall", nothing)) intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; intr.npert = 1 FFS.sing_lim!(intr, ctrl, equil) FFS.sing_find!(intr, equil) diff --git a/docs/src/vacuum.md b/docs/src/vacuum.md index 6fdff539e..69109b1ca 100644 --- a/docs/src/vacuum.md +++ b/docs/src/vacuum.md @@ -22,11 +22,11 @@ Contains plasma boundary data and calculation parameters including: - Grid resolution (mtheta) for vacuum calculations - Optional kernel sign and symmetry flags -### WallShapeSettings -Specifies wall geometry configuration with options for: -- No wall, conformal, elliptical, dee-shaped, or custom walls -- Geometric parameters (a, aw, bw, cw, dw, tw) -- Equal arc length spacing option +### WallSettings +Abstract type of the wall shapes, one concrete type each with only its own parameters: +- `ConformalWall`, `EllipticalWall`, `DeeWall`, `ModDeeWall`, `FileWall` +- Equal arc length spacing option (`equal_arc_wall`) +- No wall is `nothing` rather than a shape: the `[Wall]` section is omitted from the deck, and `compute_vacuum_response(inputs, nothing)` leaves the vacuum region unbounded ## API Reference @@ -61,10 +61,9 @@ inputs = GeneralizedPerturbedEquilibrium.Vacuum.VacuumInput( kernelsign = 1.0 # Kernel sign (+1 or -1) ) -# Define wall settings -wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings( - shape = "conformal", # Wall shape type - a = 0.3, # Wall distance parameter +# Define wall settings (use `nothing` for no wall) +wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.ConformalWall( + a = 0.3, # Wall distance in units of the minor radius equal_arc_wall = true # Use equal arc length spacing ) @@ -77,5 +76,5 @@ vac = GeneralizedPerturbedEquilibrium.Vacuum.compute_vacuum_response(inputs, wal - The Julia implementation uses Bulirsch's algorithm for elliptic integrals, providing improved accuracy over polynomial approximations - For large mode numbers (nρ̂ ≥ 0.1), 32-point Gaussian quadrature is used for Legendre function evaluation - For n=0 modes with closed walls, automatic regularization is applied -- Wall shapes support: nowall, conformal, elliptical, dee, mod_dee, or custom from file +- Wall shapes: conformal, elliptical, dee, mod_dee, or custom from file; no wall is `nothing` - The vacuum response matrix wv is scaled by the singular factor (m - nq)(m' - nq) per [Chance Phys. Plasmas 1997] diff --git a/docs/src/workflow.md b/docs/src/workflow.md index 86a7ae0df..0264b7c05 100644 --- a/docs/src/workflow.md +++ b/docs/src/workflow.md @@ -51,11 +51,11 @@ The single `gpec.toml` file supplies user-selected options to every module. The **Purpose**: Compute the magnetostatic vacuum response matrices that couple the plasma boundary to external fields. These matrices are used by ForceFreeStates and PerturbedEquilibrium to apply boundary conditions and decompose fields in poloidal mode space. -**Configuration**: `[Wall]` section of `gpec.toml` +**Configuration**: `[Wall]` section of `gpec.toml`; omit it for no wall **Inputs**: - `PlasmaEquilibrium` from Stage 1 (plasma boundary coordinates, toroidal mode number n) -- Wall geometry: shape (conformal, elliptical, dee-shaped, or custom) and dimensions +- Wall geometry (optional): shape (conformal, elliptical, dee-shaped, or custom) and dimensions; with no `[Wall]` section the vacuum region is unbounded **Outputs**: - `wv` — Vacuum response matrix (scaled by the singular factor (m - nq)(m' - nq), see Chance 1997) diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml index a8a99f4c7..ca88ff03c 100644 --- a/examples/DIIID-like_SLAYER_example/gpec.toml +++ b/examples/DIIID-like_SLAYER_example/gpec.toml @@ -18,16 +18,6 @@ newq0 = 0 # Target on-axis safety factor q(0) (0 = use inpu etol = 1e-10 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) -[Wall] -shape = "nowall" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) -equal_arc_wall = true # Equal arc length distribution of nodes on wall - [ForceFreeStates] force_termination = true # Run FFS + SLAYER, skip PerturbedEquilibrium local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile diff --git a/examples/DIIID-like_error_field_example/gpec.toml b/examples/DIIID-like_error_field_example/gpec.toml index 9a35ba9d0..c66a30ed9 100644 --- a/examples/DIIID-like_error_field_example/gpec.toml +++ b/examples/DIIID-like_error_field_example/gpec.toml @@ -22,16 +22,6 @@ newq0 = 0 # Target on-axis safety factor q(0) (0 = use inpu etol = 1e-10 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) -[Wall] -shape = "nowall" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) -equal_arc_wall = true # Equal arc length distribution of nodes on wall - [ForceFreeStates] local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes diff --git a/examples/DIIID-like_gal_resistive_example/gpec.toml b/examples/DIIID-like_gal_resistive_example/gpec.toml index 71d16afcb..998a74315 100644 --- a/examples/DIIID-like_gal_resistive_example/gpec.toml +++ b/examples/DIIID-like_gal_resistive_example/gpec.toml @@ -18,16 +18,6 @@ newq0 = 0 # Target on-axis safety factor q(0) (0 = use inpu etol = 1e-10 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) -[Wall] -shape = "nowall" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) -equal_arc_wall = true # Equal arc length distribution of nodes on wall - [ForceFreeStates] local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes diff --git a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml index 71e584905..e001feb53 100644 --- a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml +++ b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml @@ -18,16 +18,6 @@ newq0 = 0 # Target on-axis safety factor q(0) (0 = use inpu etol = 1e-10 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) -[Wall] -shape = "nowall" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) -equal_arc_wall = true # Equal arc length distribution of nodes on wall - [ForceFreeStates] local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml index 2aca73d8a..04ef55818 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/examples/DIIID-like_ideal_example/gpec.toml @@ -17,16 +17,6 @@ newq0 = 0 # Target on-axis safety factor q(0) (0 = use inpu etol = 1e-10 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) -[Wall] -shape = "nowall" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) -equal_arc_wall = true # Equal arc length distribution of nodes on wall - [ForceFreeStates] local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes diff --git a/examples/DIIID-like_ideal_example_IMAS/gpec.toml b/examples/DIIID-like_ideal_example_IMAS/gpec.toml index 0893460a2..77690328f 100644 --- a/examples/DIIID-like_ideal_example_IMAS/gpec.toml +++ b/examples/DIIID-like_ideal_example_IMAS/gpec.toml @@ -17,9 +17,6 @@ mpsi = 128 # Number of radial grid intervals (0 = two-pass a mtheta = 256 # Number of poloidal grid points imas_cocos = 11 # Expected COCOS convention of the input IMAS dd (11 = IMAS standard) -[Wall] -shape = "nowall" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) - [ForceFreeStates] nn_low = 1 # Smallest toroidal mode number to include nn_high = 1 # Largest toroidal mode number to include diff --git a/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl b/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl index c4fe1637f..16590b557 100644 --- a/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl +++ b/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl @@ -3,7 +3,7 @@ IMAS Integration Benchmark — DIII-D-like H-mode equilibrium Demonstrates that GPEC produces the same stability eigenvalue when reading the equilibrium from an IMAS data dictionary (IMASdd.dd) as when reading directly -from the g-file. Both runs use identical ForceFreeStates and Wall settings so +from the g-file. Both runs use identical ForceFreeStates settings and no wall so the only difference is how the equilibrium is loaded. EFIT.jl is NOT a GPEC dependency. Install it in your own environment before running: @@ -22,7 +22,7 @@ using EFIT const GEQDSK = joinpath(@__DIR__, "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk") @assert isfile(GEQDSK) "G-file not found at $GEQDSK" -# Shared settings: both runs use identical ForceFreeStates and Wall config. +# Shared settings: both runs use identical ForceFreeStates config and no wall. # Only the [Equilibrium] section differs between IMAS and g-file paths. const SHARED_CONFIG = Dict( "ForceFreeStates" => Dict( @@ -34,8 +34,7 @@ const SHARED_CONFIG = Dict( "verbose" => false, "delta_mhigh" => 16, "qlow" => 1.02 - ), - "Wall" => Dict("shape" => "nowall") + ) ) const EQ_COMMON = Dict( diff --git a/examples/DIIID-like_ideal_multi_n_example/gpec.toml b/examples/DIIID-like_ideal_multi_n_example/gpec.toml index 002596c21..09ed15472 100644 --- a/examples/DIIID-like_ideal_multi_n_example/gpec.toml +++ b/examples/DIIID-like_ideal_multi_n_example/gpec.toml @@ -18,16 +18,6 @@ newq0 = 0 # Target on-axis safety factor q(0) (0 = use inpu etol = 1e-10 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) -[Wall] -shape = "nowall" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) -equal_arc_wall = true # Equal arc length distribution of nodes on wall - [ForceFreeStates] local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes diff --git a/examples/DIIID-like_riccati_deltaprime_example/gpec.toml b/examples/DIIID-like_riccati_deltaprime_example/gpec.toml index 99a6f80c9..6612d2c9f 100644 --- a/examples/DIIID-like_riccati_deltaprime_example/gpec.toml +++ b/examples/DIIID-like_riccati_deltaprime_example/gpec.toml @@ -18,16 +18,6 @@ newq0 = 0 # Override for on-axis safety factor (0 = use inp etol = 1e-10 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) -[Wall] -shape = "nowall" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) -equal_arc_wall = true # Equal arc length distribution of nodes on wall - [ForceFreeStates] local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes diff --git a/examples/LAR_beta_scan/gpec.toml b/examples/LAR_beta_scan/gpec.toml index 88770631b..3a53547c1 100644 --- a/examples/LAR_beta_scan/gpec.toml +++ b/examples/LAR_beta_scan/gpec.toml @@ -30,8 +30,8 @@ ma = 128 # TJ-analytic internal radial grid resolution (sh mtau = 128 # TJ-analytic internal poloidal grid resolution (θ-spline nodes) [Wall] -shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 20 # Distance from plasma (conformal) or shape parameter (20× → effectively no wall) +shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file) +a = 20 # Distance from plasma in units of the minor radius (20× → effectively no wall) [ForceFreeStates] local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile diff --git a/examples/LAR_epsilon_scan/gpec.toml b/examples/LAR_epsilon_scan/gpec.toml index 393d093b6..3371948e5 100644 --- a/examples/LAR_epsilon_scan/gpec.toml +++ b/examples/LAR_epsilon_scan/gpec.toml @@ -31,8 +31,8 @@ ma = 128 # TJ-analytic internal radial grid resolution (sh mtau = 128 # TJ-analytic internal poloidal grid resolution (θ-spline nodes) [Wall] -shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 20 # Distance from plasma (conformal) or shape parameter (20× → effectively no wall) +shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file) +a = 20 # Distance from plasma in units of the minor radius (20× → effectively no wall) [ForceFreeStates] local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile diff --git a/examples/LAR_ideal_match_test/gpec.toml b/examples/LAR_ideal_match_test/gpec.toml index 0706863b7..9a23f46f0 100644 --- a/examples/LAR_ideal_match_test/gpec.toml +++ b/examples/LAR_ideal_match_test/gpec.toml @@ -27,8 +27,8 @@ ma = 128 # Radial grid points mtau = 128 # Poloidal grid points [Wall] -shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 20 # Distance from plasma (conformal) or shape parameter +shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file) +a = 20 # Distance from plasma in units of the minor radius [ForceFreeStates] HDF5_filename = "gpec_lar_ideal.h5" # Output filename for the stability-stage results diff --git a/examples/LAR_resistive_match_test/gpec.toml b/examples/LAR_resistive_match_test/gpec.toml index 0571301c7..0a81661c1 100644 --- a/examples/LAR_resistive_match_test/gpec.toml +++ b/examples/LAR_resistive_match_test/gpec.toml @@ -28,8 +28,8 @@ ma = 128 # Radial grid points mtau = 128 # Poloidal grid points [Wall] -shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 20 # Distance from plasma (conformal) or shape parameter +shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file) +a = 20 # Distance from plasma in units of the minor radius [ForceFreeStates] HDF5_filename = "gpec_lar_resistive.h5" # Output filename for the stability-stage results diff --git a/examples/Solovev_ideal_example/gpec.toml b/examples/Solovev_ideal_example/gpec.toml index f74d26fa6..a9e21b502 100644 --- a/examples/Solovev_ideal_example/gpec.toml +++ b/examples/Solovev_ideal_example/gpec.toml @@ -22,13 +22,8 @@ force_termination = false # Terminate after equilibrium setup (ski # (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 (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) +shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file) +a = 0.2415 # Distance from plasma in units of the minor radius equal_arc_wall = true # Equal arc length distribution of nodes on wall [ForcingTerms] diff --git a/examples/Solovev_ideal_example_3D/gpec.toml b/examples/Solovev_ideal_example_3D/gpec.toml index a5c3e1943..a44950743 100644 --- a/examples/Solovev_ideal_example_3D/gpec.toml +++ b/examples/Solovev_ideal_example_3D/gpec.toml @@ -44,13 +44,8 @@ dmlim = 0.2 # Truncate integration at (last_rational_q + dmli [Wall] # See examples/Solovev_ideal_example/gpec.toml for the rationale behind the close conformal wall. -shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) +shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file) +a = 0.2415 # Distance from plasma in units of the minor radius equal_arc_wall = false # Equal arc length distribution of nodes on wall # Uncomment to enable perturbed equilibrium calculations: diff --git a/examples/Solovev_ideal_example_multi_n/gpec.toml b/examples/Solovev_ideal_example_multi_n/gpec.toml index 7670bd5cf..5cdc338d7 100644 --- a/examples/Solovev_ideal_example_multi_n/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/gpec.toml @@ -17,13 +17,8 @@ force_termination = false # Terminate after equilibrium setup (ski [Wall] # See examples/Solovev_ideal_example/gpec.toml for the rationale behind the close conformal wall. -shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) +shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file) +a = 0.2415 # Distance from plasma in units of the minor radius equal_arc_wall = true # Equal arc length distribution of nodes on wall [ForceFreeStates] diff --git a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml index 59760b231..120ac994e 100644 --- a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml @@ -16,13 +16,8 @@ force_termination = false # Terminate after equilibrium setup (ski [Wall] # See examples/Solovev_ideal_example/gpec.toml for the rationale behind the close conformal wall. -shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) +shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file) +a = 0.2415 # Distance from plasma in units of the minor radius equal_arc_wall = true # Equal arc length distribution of nodes on wall [ForceFreeStates] diff --git a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml index a80ac7741..efefd817f 100644 --- a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml @@ -16,13 +16,8 @@ force_termination = false # Terminate after equilibrium setup (ski [Wall] # See examples/Solovev_ideal_example/gpec.toml for the rationale behind the close conformal wall. -shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) +shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file) +a = 0.2415 # Distance from plasma in units of the minor radius equal_arc_wall = true # Equal arc length distribution of nodes on wall [ForceFreeStates] diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml index 833a5ba22..200ec3ef4 100644 --- a/examples/Solovev_kinetic_NTV_example/gpec.toml +++ b/examples/Solovev_kinetic_NTV_example/gpec.toml @@ -24,13 +24,8 @@ force_termination = false # Terminate after equilibrium setup (ski # (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 (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) +shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file) +a = 0.2415 # Distance from plasma in units of the minor radius equal_arc_wall = true # Equal arc length distribution of nodes on wall [ForcingTerms] diff --git a/examples/Solovev_kinetic_calculated_example/gpec.toml b/examples/Solovev_kinetic_calculated_example/gpec.toml index 24824aa7c..aaf9150f0 100644 --- a/examples/Solovev_kinetic_calculated_example/gpec.toml +++ b/examples/Solovev_kinetic_calculated_example/gpec.toml @@ -19,13 +19,8 @@ etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) [Wall] -shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) +shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file) +a = 0.2415 # Distance from plasma in units of the minor radius equal_arc_wall = true # Equal arc length distribution of nodes on wall [ForcingTerms] diff --git a/examples/Solovev_kinetic_multiion_example/gpec.toml b/examples/Solovev_kinetic_multiion_example/gpec.toml index b90e628ee..09b1ca52f 100644 --- a/examples/Solovev_kinetic_multiion_example/gpec.toml +++ b/examples/Solovev_kinetic_multiion_example/gpec.toml @@ -24,13 +24,8 @@ force_termination = false # Terminate after equilibrium setup (ski # (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 (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) +shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file) +a = 0.2415 # Distance from plasma in units of the minor radius equal_arc_wall = true # Equal arc length distribution of nodes on wall [ForcingTerms] diff --git a/examples/a10_kinetic_example/gpec.toml b/examples/a10_kinetic_example/gpec.toml index eb32a2b1e..1dad070a0 100644 --- a/examples/a10_kinetic_example/gpec.toml +++ b/examples/a10_kinetic_example/gpec.toml @@ -18,9 +18,6 @@ mtheta = 256 # Number of poloidal grid points newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-7 # Error tolerance for equilibrium solver -[Wall] -shape = "nowall" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) - [ForceFreeStates] local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index 5544b0e60..77b77b231 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/src/ForceFreeStates/CoreTypes.jl @@ -51,7 +51,7 @@ A mutable struct holding internal state variables for stability calculations. - `psilim::Float64` - Flux limit for integration - `qlim::Float64` - Safety factor at psilim - `q1lim::Float64` - Safety factor derivative at psilim - - `wall_settings::Vacuum.WallShapeSettings` - Wall shape settings for vacuum calculations + - `wall_settings::Union{Nothing,Vacuum.WallSettings}` - Wall shape for vacuum calculations; `nothing` for no wall """ @kwdef mutable struct ForceFreeStatesInternal <: ModeSpace dir_path::String = "" @@ -74,7 +74,7 @@ A mutable struct holding internal state variables for stability calculations. qlim::Float64 = 0.0 q1lim::Float64 = 0.0 debug_settings::DebugSettings = DebugSettings() - wall_settings::Vacuum.WallShapeSettings = Vacuum.WallShapeSettings() + wall_settings::Union{Nothing,Vacuum.WallSettings} = nothing """ Inter-surface Δ' matrix of shape (msing × msing) in PEST3 convention. Computed by `compute_delta_prime_matrix!` (parallel FM path only) using the STRIDE diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index 32509e77a..1c5a24dbd 100644 --- a/src/ForceFreeStates/Free.jl +++ b/src/ForceFreeStates/Free.jl @@ -16,7 +16,7 @@ power-normalized and invariant to the working (Jacobian) coordinate. - `et::Vector{ComplexF64}` - Total energy eigenvalues of the pencil (W, N); et = ep + ev per mode - `n_tor_idx::Vector{Int}` - 0-based toroidal mode number index of each sorted eigenvalue - `vacuum_eigenvalue::Float64` - Least stable (minimum) eigenvalue of the pencil (wv, N), clamped to zero - - `plasma_pts`, `wall_pts::Matrix{Float64}` - Cartesian (x, y, z) surface coordinates, `numpoints × 3`, retained for HDF5 output + - `plasma_pts`, `wall_pts::Matrix{Float64}` - Cartesian (x, y, z) surface coordinates, `numpoints × 3`, retained for HDF5 output; `wall_pts` has no rows when there is no wall """ struct FreeBoundaryResult wt::Matrix{ComplexF64} @@ -202,7 +202,7 @@ calculations and data dumping. " Total = $((@sprintf "%+.3e %+.3ei" real(et[1]) imag(et[1])))" end - return FreeBoundaryResult(wt, wt0, wp, wv, ep, ev, et, n_tor_idx, vacuum_eigenvalue, vac.plasma_pts, vac.wall_pts) + return FreeBoundaryResult(wt, wt0, wp, wv, ep, ev, et, n_tor_idx, vacuum_eigenvalue, vac.plasma_pts, isnothing(intr.wall_settings) ? zeros(0, 3) : vac.wall_pts) end """ diff --git a/src/ForceFreeStates/Result.jl b/src/ForceFreeStates/Result.jl index da8fe3648..2f99d55f1 100644 --- a/src/ForceFreeStates/Result.jl +++ b/src/ForceFreeStates/Result.jl @@ -52,7 +52,7 @@ so bpen and closure are always present. - `mlow`, `mhigh`, `mpert`, `nlow`, `nhigh`, `npert`, `numpert_total` - Resolved (m, n) mode space; see `ModeSpace`. - `psilow`, `psilim`, `qlim`, `q1lim` - Integration bounds, and q with its ψ-derivative at `psilim`. - `dir_path::String` - Working directory of the run. - - `wall_settings::Vacuum.WallShapeSettings` - Wall shape used by the vacuum calculation. + - `wall_settings::Union{Nothing,Vacuum.WallSettings}` - Wall shape used by the vacuum calculation; `nothing` for no wall. - `debug_settings::DebugSettings` - Diagnostic dump settings (the `[DEBUG]` deck section / `debug=` API keyword). - `metric::MetricData`, `mats::MatrixSplines` - Metric data and Euler-Lagrange matrix interpolants. - `surfaces::Vector{SingType}` - Ideal singular surfaces in the integration domain, with asymptotic bases and GGJ coefficients. @@ -94,7 +94,7 @@ struct ForceFreeStatesResult{E<:Equilibrium.PlasmaEquilibrium,F<:MatrixSplines} qlim::Float64 q1lim::Float64 dir_path::String - wall_settings::Vacuum.WallShapeSettings + wall_settings::Union{Nothing,Vacuum.WallSettings} debug_settings::DebugSettings # Assembly products, always present. diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 23de76dad..0a30c9bcf 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -232,11 +232,7 @@ function main_from_inputs( return end - if "Wall" in keys(inputs) - intr.wall_settings = Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in inputs["Wall"])...) - else - intr.wall_settings = Vacuum.WallShapeSettings() - end + intr.wall_settings = Vacuum.parse_wall(get(inputs, "Wall", nothing)) if "DEBUG" in keys(inputs) intr.debug_settings = DebugSettings(; (Symbol(k) => v for (k, v) in inputs["DEBUG"])...) @@ -693,7 +689,7 @@ function run_force_free_states( # Compute free boundary energies. if ctrl.vac_flag && !(ctrl.ksing > 0 && ctrl.ksing <= intr.msing + 1) if ctrl.verbose - wall_desc = intr.wall_settings.shape == "nowall" ? "no wall" : intr.wall_settings.shape + wall_desc = isnothing(intr.wall_settings) ? "no wall" : nameof(typeof(intr.wall_settings)) @info "Computing free boundary energies ($wall_desc)" end free_energies = free_run(odet, ctrl, equil, mats, intr) @@ -727,7 +723,7 @@ function run_force_free_states( end """ - EulerLagrangeProblem(equil; nn, wall=Vacuum.WallShapeSettings(), match=nothing, + EulerLagrangeProblem(equil; nn, wall=nothing, match=nothing, dir_path=".", debug=DebugSettings(), kwargs...) The perturbed-plasma Euler-Lagrange problem posed on an equilibrium: the extremization of @@ -736,7 +732,7 @@ kinetic) perturbed states. This is the WHAT of a stability solve; the integrator [`solve`](@ref) is the HOW. A `PlasmaEquilibrium` hosts many possible problems — this type names this one, so `solve` stays unambiguous as other problem classes appear. -`nn` is the toroidal mode number or range. `wall` is the vacuum wall shape, `match` an +`nn` is the toroidal mode number or range. `wall` is the vacuum `Vacuum.WallSettings`, or `nothing` for no wall, `match` an optional [`ResistiveMatch`](@ref) closing the basis with an inner-layer solution instead of the ideal jump, `dir_path` the working directory outputs are written to, and `debug` the diagnostic dump settings of the DEBUG deck section. Any remaining keyword is a @@ -746,7 +742,7 @@ knobs. `nn_low`/`nn_high` are rejected — they come from `nn`. ## Fields - `equil::Equilibrium.PlasmaEquilibrium` - The equilibrium the problem is posed on. - - `wall::Vacuum.WallShapeSettings` - Vacuum wall shape for the free-boundary energies. + - `wall::Union{Nothing,Vacuum.WallSettings}` - Vacuum wall shape for the free-boundary energies; `nothing` for no wall. - `match::Union{Nothing,ForceFreeStates.ResistiveMatch}` - Optional inner-layer closure. - `dir_path::String` - Working directory for outputs. - `debug::DebugSettings` - Diagnostic dump settings. @@ -754,7 +750,7 @@ knobs. `nn_low`/`nn_high` are rejected — they come from `nn`. """ struct EulerLagrangeProblem equil::Equilibrium.PlasmaEquilibrium - wall::Vacuum.WallShapeSettings + wall::Union{Nothing,Vacuum.WallSettings} match::Union{Nothing,ForceFreeStates.ResistiveMatch} dir_path::String debug::DebugSettings @@ -764,7 +760,7 @@ end function EulerLagrangeProblem( equil::Equilibrium.PlasmaEquilibrium; nn::Union{Int,AbstractUnitRange{Int}}, - wall::Vacuum.WallShapeSettings=Vacuum.WallShapeSettings(), + wall::Union{Nothing,Vacuum.WallSettings}=nothing, match::Union{Nothing,ForceFreeStates.ResistiveMatch}=nothing, dir_path::AbstractString=".", debug::DebugSettings=DebugSettings(), @@ -1217,7 +1213,9 @@ function efc_couplings( ψk = kp.xs ω_E = [kp.omegaE_spline(ψ) for ψ in ψk] ω_star_T = [-2π * kp.Ti_deriv(ψ) / (chrg * chi1) for ψ in ψk] - ω_star_n = [(n = kp.ni_spline(ψ); n > 0 ? -2π * kp.Ti_spline(ψ) * kp.ni_deriv(ψ) / (chrg * chi1 * n) : 0.0) for ψ in ψk] + ω_star_n = [ + (n=kp.ni_spline(ψ); n > 0 ? -2π * kp.Ti_spline(ψ) * kp.ni_deriv(ψ) / (chrg * chi1 * n) : 0.0) for ψ in ψk + ] ω_φ = ω_E .+ ω_star_n .+ ω_star_T # Momentum weight n_i·dV/dψ with ⟨R²⟩ taken as R₀², inside the plasma only. w = [ψ <= ffs.psilim ? max(kp.ni_spline(ψ), 0.0) * equil.profiles.dVdpsi_spline(clamp(ψ, equil.profiles.xs[1], equil.profiles.xs[end])) : 0.0 for ψ in ψk] diff --git a/src/PerturbedEquilibrium/ResponseMatrices.jl b/src/PerturbedEquilibrium/ResponseMatrices.jl index 405b3bebf..07fc7a3d6 100644 --- a/src/PerturbedEquilibrium/ResponseMatrices.jl +++ b/src/PerturbedEquilibrium/ResponseMatrices.jl @@ -299,8 +299,7 @@ function calc_surface_inductance( n_modes::Union{Integer,AbstractVector{<:Integer}} )::Matrix{ComplexF64} vac_input = Vacuum.VacuumInput(equil, ψ, mtheta, 1, m_modes, vcat(n_modes)) - wall_settings = Vacuum.WallShapeSettings(; shape="nowall") - I_v = Vacuum.compute_vacuum_response(vac_input, wall_settings; compute_Iv=true).I_v + I_v = Vacuum.compute_vacuum_response(vac_input, nothing; compute_Iv=true).I_v μ₀ = 4π * 1e-7 return inv(I_v) .* (μ₀ * (2π)^2) end diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl index a817e9a27..4f7a58897 100644 --- a/src/Vacuum/DataTypes.jl +++ b/src/Vacuum/DataTypes.jl @@ -158,7 +158,8 @@ boundary-integral solve produces along the way. `μ₀`/`4π²` normalization: the physical surface inductance is `μ₀(2π)²·I_v⁻¹` (see `PerturbedEquilibrium.calc_surface_inductance`). In 3D it differences two solves, so it needs a finer toroidal grid than `wv` to converge. - - `plasma_pts`, `wall_pts::Matrix{Float64}`: Cartesian surface coordinates (`num_points × 3`) + - `plasma_pts`, `wall_pts::Matrix{Float64}`: Cartesian surface coordinates (`num_points × 3`); `wall_pts` is zero + when there is no wall """ struct VacuumResponse wv::Matrix{ComplexF64} @@ -184,63 +185,146 @@ function VacuumResponse(inputs::VacuumInput) end """ - WallShapeSettings + WallSettings -Struct containing input settings for vacuum wall geometry. +Abstract supertype of the vacuum wall shapes: [`ConformalWall`](@ref), [`EllipticalWall`](@ref), +[`DeeWall`](@ref), [`ModDeeWall`](@ref) and [`FileWall`](@ref). A run with no wall has no +`WallSettings` at all; the wall is `nothing`, and the deck omits its `[Wall]` section. -# Fields +Every shape carries `equal_arc_wall::Bool` (default `true`), which redistributes the wall nodes +to equal arc length (recommended unless the wall is very close to the plasma). Re-parameterizing +the contour breaks the plasma/wall grid index correspondence that the singular quadrature relies +on, so it is ignored outright for `nzeta_in > 1` walls. + +A non-axisymmetric boundary (`nzeta_in > 1`) supports only [`ConformalWall`](@ref); the other +shapes are poloidal contours that get revolved and so need `nzeta_in == 1`. +""" +abstract type WallSettings end - - `shape::String`: String selecting wall shape. Options are: +""" + ConformalWall(; a=0.3, equal_arc_wall=true) + +Wall offset from the plasma surface along its normal. + +## Fields - + `"nowall"`: No wall - + `"conformal"`: Wall conformal to plasma surface at distance `a` - + `"elliptical"`: Elliptical wall - + `"dee"`: Dee-shaped wall - + `"mod_dee"`: Modified Dee-shaped wall - + `"filepath"`: Custom wall shape from the file you specify + - `a::Float64`: Distance of the wall from the plasma in units of the minor radius `0.5(max R - min R)`. + On a non-axisymmetric boundary the extrema are taken over the whole torus, so `a` scales with half + the global R-extent rather than a cross-section minor radius. + - `equal_arc_wall::Bool`: Equal arc length node distribution, see [`WallSettings`](@ref) +""" +@kwdef struct ConformalWall <: WallSettings + a::Float64 = 0.3 + equal_arc_wall::Bool = true +end - A non-axisymmetric boundary (`nzeta_in > 1`) supports `"nowall"` and `"conformal"`; the others - are poloidal contours that get revolved and so need `nzeta_in == 1`. +""" + EllipticalWall(; a=0.3, equal_arc_wall=true) - - `a::Float64`: Distance of wall from plasma in units of the minor radius `0.5(max R - min R)` - (conformal), or shape parameter (others). On a non-axisymmetric boundary the extrema are taken - over the whole torus, so `a` scales with half the global R-extent rather than a cross-section - minor radius; it reduces to the axisymmetric definition when the boundary is axisymmetric. +Ellipse confocal with the plasma's bounding ellipse. - - `aw::Float64`: Half-thickness parameter for Dee-shaped walls +## Fields - - `bw::Float64`: Elongation parameter for wall shapes + - `a::Float64`: Horizontal semi-axis in meters; the vertical semi-axis follows from the confocality + - `equal_arc_wall::Bool`: Equal arc length node distribution, see [`WallSettings`](@ref) +""" +@kwdef struct EllipticalWall <: WallSettings + a::Float64 = 0.3 + equal_arc_wall::Bool = true +end - - `cw::Float64`: Offset of the center of the wall from the major radius +""" + DeeWall(; a=0.3, aw=0.05, bw=1.5, cw=0.0, dw=0.5, tw=0.05, equal_arc_wall=true) - - `dw::Float64`: Triangularity parameter for wall shapes +Dee-shaped wall scaled to the plasma: minor radius `r_minor (1 + a - cw)`, centered at +`r_major + cw r_minor`. - - `tw::Float64`: Sharpness of the corners of the wall (try 0.05 as initial value) +## Fields - # Core shape selection + - `a::Float64`: Distance of the wall from the plasma in units of the minor radius `0.5(max R - min R)` + - `aw::Float64`: Half-thickness parameter + - `bw::Float64`: Elongation parameter + - `cw::Float64`: Offset of the wall center from the major radius, in units of the minor radius + - `dw::Float64`: Triangularity parameter + - `tw::Float64`: Sharpness of the corners (try 0.05 as initial value) + - `equal_arc_wall::Bool`: Equal arc length node distribution, see [`WallSettings`](@ref) +""" +@kwdef struct DeeWall <: WallSettings + a::Float64 = 0.3 + aw::Float64 = 0.05 + bw::Float64 = 1.5 + cw::Float64 = 0.0 + dw::Float64 = 0.5 + tw::Float64 = 0.05 + equal_arc_wall::Bool = true +end - - `equal_arc_wall::Bool`: Flag to enforce equal arc length distribution of nodes on the wall - (recommended unless wall is very close to plasma). Re-parameterizing the contour breaks the - plasma/wall grid index correspondence that the singular quadrature relies on; it is ignored - outright for `nzeta_in > 1` walls. """ -@kwdef struct WallShapeSettings + ModDeeWall(; a=0.3, aw=0.05, bw=1.5, cw=0.0, dw=0.5, tw=0.05, equal_arc_wall=true) - # Core shape selection - shape::String = "nowall" +Dee-shaped wall in absolute coordinates, independent of the plasma. Same shape parameters as +[`DeeWall`](@ref), except that `a` and `cw` are lengths in meters. - # Standard geometric parameters for Dee/Mod-Dee +## Fields + + - `a::Float64`: Minor radius of the wall in meters + - `aw::Float64`: Half-thickness parameter + - `bw::Float64`: Elongation parameter + - `cw::Float64`: Major radius of the wall center in meters + - `dw::Float64`: Triangularity parameter + - `tw::Float64`: Sharpness of the corners (try 0.05 as initial value) + - `equal_arc_wall::Bool`: Equal arc length node distribution, see [`WallSettings`](@ref) +""" +@kwdef struct ModDeeWall <: WallSettings a::Float64 = 0.3 aw::Float64 = 0.05 bw::Float64 = 1.5 cw::Float64 = 0.0 dw::Float64 = 0.5 tw::Float64 = 0.05 + equal_arc_wall::Bool = true +end + +""" + FileWall(; file, equal_arc_wall=true) - # Algorithmic options +Wall contour read from a text file: the point count, a reference major radius and a header line, +then one `index R Z` row per point, with at least `mtheta` points. + +## Fields + + - `file::String`: Path to the wall file + - `equal_arc_wall::Bool`: Equal arc length node distribution, see [`WallSettings`](@ref) +""" +@kwdef struct FileWall <: WallSettings + file::String equal_arc_wall::Bool = true end +const WALL_SHAPES = Dict("conformal" => ConformalWall, "elliptical" => EllipticalWall, "dee" => DeeWall, "mod_dee" => ModDeeWall, "file" => FileWall) + +""" + parse_wall(section) -> Union{Nothing,WallSettings} + +Build the wall from the `[Wall]` section of a deck, or return `nothing` (no wall) when there is no +section. The section's `shape` picks the [`WallSettings`](@ref) type. Keys that type does not use are +ignored with a warning, and a parameter it requires (such as `file` for `FileWall`) is an error if missing. +""" +parse_wall(::Nothing) = nothing +function parse_wall(section::AbstractDict) + params = Dict{String,Any}(section) + shape = pop!(params, "shape", nothing) + haskey(WALL_SHAPES, shape) || error("[Wall] needs a `shape` of $(join(sort!(collect(keys(WALL_SHAPES))), ", ")), got $(repr(shape)). For no wall, omit the [Wall] section.") + T = WALL_SHAPES[shape] + unused = setdiff(keys(params), string.(fieldnames(T))) + isempty(unused) || @warn "[Wall] shape \"$shape\" does not use $(join(sort!(collect(unused)), ", ")); ignoring it. Its parameters are $(join(fieldnames(T), ", "))." + try + return T(; (Symbol(k) => v for (k, v) in params if !(k in unused))...) + catch err + err isa UndefKeywordError ? error("[Wall] shape \"$shape\" requires `$(err.var)`.") : rethrow() + end +end + """ PlasmaGeometry @@ -267,12 +351,10 @@ Struct holding wall geometry data for vacuum calculations. Arrays are of length # Fields - - `nowall::Bool`: Boolean flag indicating if there is no wall - `x::Vector{Float64}`: Wall R-coordinates - `z::Vector{Float64}`: Wall Z-coordinates """ struct WallGeometry - nowall::Bool x::Vector{Float64} z::Vector{Float64} end @@ -429,144 +511,142 @@ function PlasmaGeometry3D(inputs::VacuumInput) end """ - WallGeometry(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_settings::WallShapeSettings) -> WallGeometry - -Constructor to initialize the wall geometry based on the provided vacuum inputs and wall shape settings. - -This performs functionality similar to portions of the `arrays` function in the original -Fortran VACUUM code. It returns a `WallGeometry` struct containing the necessary wall -surface data for vacuum calculations. - -# Arguments - - - `inputs::VacuumInput`: Struct containing vacuum calculation parameters - - `plasma_surf::PlasmaGeometry`: Struct with plasma surface geometry (used for reference) - - `wall_settings::WallShapeSettings`: Struct specifying wall shape and parameters + wall_contour(wall::WallSettings, plasma_surf::PlasmaGeometry) -> (x_wall, z_wall) -# Returns +Poloidal wall contour on the plasma's θ grid, before the equal-arc redistribution and validity +checks that [`WallGeometry`](@ref) applies to every shape. +""" +function wall_contour end - - `WallGeometry`: Struct containing wall surface coordinates and derivatives +# Minor and major radius from the plasma's R-extent +function minor_major_radius(plasma_surf::PlasmaGeometry) + xmin, xmax = extrema(plasma_surf.x) + return 0.5 * (xmax - xmin), 0.5 * (xmax + xmin) +end -# Notes +poloidal_angles(mtheta::Int) = (0:(mtheta-1)) .* (2π / mtheta) - - Supports multiple wall shapes: nowall, conformal, elliptical, dee, mod_dee, from_file - - Optionally redistributes wall points to equal arc length spacing if `equal_arc_wall=true` -""" -function WallGeometry(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_settings::WallShapeSettings) +function wall_contour(wall::ConformalWall, plasma_surf::PlasmaGeometry) + (; a) = wall + x_plasma = plasma_surf.x + z_plasma = plasma_surf.z + mtheta = length(x_plasma) + r_minor, _ = minor_major_radius(plasma_surf) + dx = a * r_minor + @info "Calculating conformal wall shape $((@sprintf "%.2e" dx)) m from plasma surface." - # Output wall coordinate arrays - mtheta = inputs.mtheta x_wall = zeros(mtheta) z_wall = zeros(mtheta) + centerstack_min = min(0.1, 0.1 * minimum(x_plasma)) # Avoid wall crossing R=0 axis + for i in 1:mtheta + j = mod1(i - 1, mtheta) + k = mod1(i + 1, mtheta) + # Normal vector calculation + alph = atan(x_plasma[j] - x_plasma[k], z_plasma[k] - z_plasma[j]) + x_wall[i] = max(centerstack_min, x_plasma[i] + a * r_minor * cos(alph)) + z_wall[i] = z_plasma[i] + a * r_minor * sin(alph) + end - if wall_settings.shape == "nowall" - return WallGeometry( - true, - x_wall, - z_wall - ) + if any(x_wall .<= centerstack_min + eps(Float64)) + @warn "Conformal wall with a=$a would cross R=0 axis; forcing minimum wall R to $(@sprintf "%.2e" centerstack_min) m to avoid unphysical geometry." end + return x_wall, z_wall +end - # Compute plasma surface quantities - x_plasma = plasma_surf.x - z_plasma = plasma_surf.z +function wall_contour(wall::EllipticalWall, plasma_surf::PlasmaGeometry) + (; a) = wall + @info "Calculating elliptical wall shape with a = $((@sprintf "%.2e" a)) m." + r_minor, r_major = minor_major_radius(plasma_surf) + zmin, zmax = extrema(plasma_surf.z) + + zrad = 0.5 * (zmax - zmin) + zh = sqrt(abs(zrad^2 - r_minor^2)) + zmuw = log((a/zh) + sqrt((a/zh)^2 + 1)) + bw_eff = (zh * cosh(zmuw)) / a + + the = poloidal_angles(length(plasma_surf.x)) + x_wall = @. r_major + a * cos(the) + z_wall = @. bw_eff * a * sin(the) + return x_wall, z_wall +end - # Common geometric parameters - xmin = minimum(x_plasma) - xmax = maximum(x_plasma) - zmin = minimum(z_plasma) - zmax = maximum(z_plasma) +function wall_contour(wall::DeeWall, plasma_surf::PlasmaGeometry) + (; a, aw, bw, cw, dw, tw) = wall + r_minor, r_major = minor_major_radius(plasma_surf) + wcentr = r_major + cw * r_minor + @info "Calculating dee-shaped wall with R = $((@sprintf "%.2e" wcentr)) + $((@sprintf "%.2e" r_minor)) * (1.0 + $((@sprintf "%.2e" a)) - $((@sprintf "%.2e" cw))) * cos(θ + $((@sprintf "%.2e" dw)) * sin(θ)), Z = $((@sprintf "%.2e" bw)) * $((@sprintf "%.2e" r_minor)) * (1.0 + $((@sprintf "%.2e" a)) - $((@sprintf "%.2e" cw))) * sin(θ + $((@sprintf "%.2e" tw)) * sin(2θ)) + $((@sprintf "%.2e" aw)) * $((@sprintf "%.2e" r_minor)) * sin(2θ)." - r_minor = 0.5 * (xmax - xmin) - r_major = 0.5 * (xmax + xmin) + the = poloidal_angles(length(plasma_surf.x)) + x_wall = @. wcentr + r_minor * (1.0 + a - cw) * cos(the + dw * sin(the)) + z_wall = @. bw * r_minor * (1.0 + a - cw) * sin(the + tw * sin(2.0*the)) + aw * r_minor * sin(2.0*the) + return x_wall, z_wall +end - # Destructuring settings for readability - (; aw, bw, cw, dw, tw, a) = wall_settings - wcentr = 0.0 # Initialize +function wall_contour(wall::ModDeeWall, plasma_surf::PlasmaGeometry) + (; a, aw, bw, cw, dw, tw) = wall + @info "Calculating modified dee-shaped wall with R = $((@sprintf "%.2e" cw)) + $((@sprintf "%.2e" a)) * cos(θ + $((@sprintf "%.2e" dw)) * sin(θ)), Z = $((@sprintf "%.2e" bw)) * $((@sprintf "%.2e" a)) * sin(θ + $((@sprintf "%.2e" tw)) * sin(2θ)) + $((@sprintf "%.2e" aw)) * sin(2θ)." - if wall_settings.shape == "conformal" - dx = a * r_minor - @info "Calculating conformal wall shape $((@sprintf "%.2e" dx)) m from plasma surface." - wcentr = r_major - centerstack_min = min(0.1, 0.1 * minimum(x_plasma)) # Avoid wall crossing R=0 axis - for i in 1:mtheta - j = mod1(i - 1, mtheta) - k = mod1(i + 1, mtheta) - # Normal vector calculation - alph = atan(x_plasma[j] - x_plasma[k], z_plasma[k] - z_plasma[j]) - x_wall[i] = max(centerstack_min, x_plasma[i] + a * r_minor * cos(alph)) - z_wall[i] = z_plasma[i] + a * r_minor * sin(alph) - end + the = poloidal_angles(length(plasma_surf.x)) + x_wall = @. cw + a * cos(the + dw * sin(the)) + z_wall = @. bw * a * sin(the + tw * sin(2.0*the)) + aw * sin(2.0*the) + return x_wall, z_wall +end - if any(x_wall .<= centerstack_min + eps(Float64)) - @warn "Conformal wall with a=$a would cross R=0 axis; forcing minimum wall R to $(@sprintf "%.2e" centerstack_min) m to avoid unphysical geometry." - end +function wall_contour(wall::FileWall, plasma_surf::PlasmaGeometry) + (; file) = wall + isfile(file) || error("Wall geometry file $file does not exist.") - elseif wall_settings.shape == "elliptical" - @info "Calculating elliptical wall shape with a = $((@sprintf "%.2e" a)) m." - wcentr = r_major + mtheta = length(plasma_surf.x) + x_wall = zeros(mtheta) + z_wall = zeros(mtheta) + open(file, "r") do io + npts = parse(Int, readline(io)) # Number of points in file + readline(io) # Reference major radius, unused + readline(io) # Skip header/comment line - zrad = 0.5 * (zmax - zmin) - zh = sqrt(abs(zrad^2 - r_minor^2)) - zmuw = log((a/zh) + sqrt((a/zh)^2 + 1)) - bw_eff = (zh * cosh(zmuw)) / a + npts < mtheta && error("Wall geometry file $file contains fewer points ($npts) than mtheta ($mtheta).") for i in 1:mtheta - the = (i - 1) * (2π / mtheta) - x_wall[i] = r_major + a * cos(the) - z_wall[i] = bw_eff * a * sin(the) + line = split(readline(io)) + # Assumes file format: [index R_coord Z_coord] + x_wall[i] = parse(Float64, line[2]) + z_wall[i] = parse(Float64, line[3]) end + end + # Wall files written for the old clockwise convention are turned round, keeping their first point + if signed_area(x_wall, z_wall) < 0 + @warn "Wall file $file runs clockwise; reversing its point order to GPEC's counter-clockwise θ." + reverse!(@view x_wall[2:end]) + reverse!(@view z_wall[2:end]) + end + return x_wall, z_wall +end - elseif wall_settings.shape == "dee" - wcentr = r_major + cw * r_minor - @info "Calculating dee-shaped wall with R = $((@sprintf "%.2e" wcentr)) + $((@sprintf "%.2e" r_minor)) * (1.0 + $((@sprintf "%.2e" a)) - $((@sprintf "%.2e" cw))) * cos(θ + $((@sprintf "%.2e" dw)) * sin(θ)), Z = $((@sprintf "%.2e" bw)) * $((@sprintf "%.2e" r_minor)) * (1.0 + $((@sprintf "%.2e" a)) - $((@sprintf "%.2e" cw))) * sin(θ + $((@sprintf "%.2e" tw)) * sin(2θ)) + $((@sprintf "%.2e" aw)) * $((@sprintf "%.2e" r_minor)) * sin(2θ)." - for i in 1:mtheta - the = (i - 1) * (2π / mtheta) - x_wall[i] = wcentr + r_minor * (1.0 + a - cw) * cos(the + dw * sin(the)) - z_wall[i] = bw * r_minor * (1.0 + a - cw) * sin(the + tw * sin(2.0*the)) + aw * r_minor * sin(2.0*the) - end +""" + WallGeometry(plasma_surf::PlasmaGeometry, wall_settings::WallSettings) -> WallGeometry - elseif wall_settings.shape == "mod_dee" - @info "Calculating modified dee-shaped wall with R = $((@sprintf "%.2e" cw)) + $((@sprintf "%.2e" a)) * cos(θ + $((@sprintf "%.2e" dw)) * sin(θ)), Z = $((@sprintf "%.2e" bw)) * $((@sprintf "%.2e" a)) * sin(θ + $((@sprintf "%.2e" tw)) * sin(2θ)) + $((@sprintf "%.2e" aw)) * sin(2θ)." - wcentr = cw - for i in 1:mtheta - the = (i - 1) * (2π / mtheta) - x_wall[i] = cw + a * cos(the + dw * sin(the)) - z_wall[i] = bw * a * sin(the + tw * sin(2.0*the)) + aw * sin(2.0*the) - end +Constructor to initialize the wall geometry from the plasma surface and the wall shape. - else - filepath = wall_settings.shape - if !isfile(filepath) - error( - "Wall geometry file $filepath does not exist. " * - "Please set the wall shape parameter to a valid file path or a built-in shape " * - "(nowall, conformal, elliptical, dee, mod_dee)." - ) - end +This performs functionality similar to portions of the `arrays` function in the original +Fortran VACUUM code. It returns a `WallGeometry` struct containing the necessary wall +surface data for vacuum calculations. - wcentr = 0.0 - open(wall_settings.shape, "r") do io - npots0 = parse(Int, readline(io)) # Number of points in file - wcentr = parse(Float64, readline(io)) - readline(io) # Skip header/comment line +# Arguments - npots0 < mtheta && error("Wall geometry file $filepath contains fewer points ($npots0) than mtheta ($mtheta).") + - `plasma_surf::PlasmaGeometry`: Struct with plasma surface geometry (used for reference) + - `wall_settings::WallSettings`: Struct specifying wall shape and parameters - for i in 1:mtheta - line = split(readline(io)) - # Assumes file format: [index R_coord Z_coord] - x_wall[i] = parse(Float64, line[2]) - z_wall[i] = parse(Float64, line[3]) - end - end - # Wall files written for the old clockwise convention are turned round, keeping their first point - if signed_area(x_wall, z_wall) < 0 - @warn "Wall file $filepath runs clockwise; reversing its point order to GPEC's counter-clockwise θ." - reverse!(@view x_wall[2:end]) - reverse!(@view z_wall[2:end]) - end - end +# Returns + + - `WallGeometry`: Struct containing wall surface coordinates and derivatives + +# Notes + + - The contour comes from [`wall_contour`](@ref) for the shape in `wall_settings` + - Optionally redistributes wall points to equal arc length spacing if `equal_arc_wall=true` +""" +function WallGeometry(plasma_surf::PlasmaGeometry, wall_settings::WallSettings) + x_wall, z_wall = wall_contour(wall_settings, plasma_surf) # Optional: Re-parameterization for equal arc length spacing of wall points if wall_settings.equal_arc_wall @@ -578,7 +658,7 @@ function WallGeometry(inputs::VacuumInput, plasma_surf::PlasmaGeometry, wall_set any(x_wall .<= 0.0) && error("Wall R-coordinates contain non-physical values (R <= 0). Check wall geometry.") assert_counterclockwise(x_wall, z_wall, "Wall") - return WallGeometry(false, x_wall, z_wall) + return WallGeometry(x_wall, z_wall) end """ @@ -589,7 +669,6 @@ ordered with the poloidal index fastest. # Fields - - `nowall::Bool`: Boolean flag indicating if there is no wall - `is_closed_toroidal::Bool`: Boolean flag indicating if the wall is a closed toroidal surface - `mtheta::Int`: Number of poloidal grid points - `nzeta::Int`: Number of toroidal grid points @@ -599,7 +678,6 @@ ordered with the poloidal index fastest. - `normal::Matrix{Float64}`: Outward normal vectors at wall """ struct WallGeometry3D - nowall::Bool is_closed_toroidal::Bool mtheta::Int nzeta::Int @@ -611,21 +689,21 @@ struct WallGeometry3D end """ - WallGeometry3D(inputs::VacuumInput, plasma_surf::PlasmaGeometry3D, wall_settings::WallShapeSettings) -> WallGeometry3D + WallGeometry3D(inputs::VacuumInput, plasma_surf::PlasmaGeometry3D, wall_settings::WallSettings) -> WallGeometry3D Constructor to initialize the 3D wall geometry based on the provided vacuum inputs, plasma surface and wall shape settings. This is the 3D counterpart of [`WallGeometry`](@ref) and selects the shape the same way: an axisymmetric boundary builds the 2D poloidal contour and revolves it, so every shape [`WallGeometry`](@ref) offers is available; a non-axisymmetric boundary offsets the plasma surface. -Expects a full-torus boundary on the same `mtheta`/`nzeta` grid as `plasma_surf`, and `nfp == 1`; +Expects a full-torus boundary on the same `mtheta`/`nzeta` grid as `plasma_surf`; [`expand_field_periods`](@ref) produces both from a per-period non-axisymmetric boundary. # Arguments - `inputs::VacuumInput`: Struct containing vacuum calculation parameters - `plasma_surf::PlasmaGeometry3D`: Plasma surface the wall is built around (used by the conformal wall) - - `wall_settings::WallShapeSettings`: Struct specifying wall shape and parameters + - `wall_settings::WallSettings`: Struct specifying wall shape and parameters # Returns @@ -633,8 +711,8 @@ Expects a full-torus boundary on the same `mtheta`/`nzeta` grid as `plasma_surf` # Notes - - Axisymmetric boundaries (`nzeta_in == 1`) support nowall, conformal, elliptical, dee, mod_dee and - from_file; non-axisymmetric boundaries support nowall and conformal + - Axisymmetric boundaries (`nzeta_in == 1`) support every [`WallSettings`](@ref) shape; + non-axisymmetric boundaries support only [`ConformalWall`](@ref) - The non-axisymmetric conformal wall displaces each plasma point along its own normal, so wall grid index `(i, j)` is the closest wall point to plasma index `(i, j)` — the correspondence the near-field patch of `compute_3D_kernel_matrices!` assumes. The revolved branch does not share it: the wall sits @@ -645,10 +723,7 @@ Expects a full-torus boundary on the same `mtheta`/`nzeta` grid as `plasma_surf` Both fold tests are local, so a global self-intersection between distant lobes — where the offset surface meets itself with aligned normals and a positive area element — passes them """ -function WallGeometry3D(inputs::VacuumInput, plasma_surf::PlasmaGeometry3D, wall_settings::WallShapeSettings) - - # Basic wall flags - nowall = wall_settings.shape == "nowall" +function WallGeometry3D(inputs::VacuumInput, plasma_surf::PlasmaGeometry3D, wall_settings::WallSettings) is_closed_toroidal = true (; mtheta, nzeta) = inputs @@ -659,14 +734,6 @@ function WallGeometry3D(inputs::VacuumInput, plasma_surf::PlasmaGeometry3D, wall "but this wall is $(mtheta)×$(nzeta). Build both from the same VacuumInput." ) - # expand_field_periods tiles a non-axisymmetric boundary to the full torus and returns nfp = 1, so - # nfp > 1 reaches here only from an axisymmetric boundary, which has no field periodicity to declare. - inputs.nfp > 1 && error( - "3D wall geometry requires nfp = 1 on an axisymmetric (nzeta_in = 1) boundary: an axisymmetric " * - "surface has no field periodicity, and nfp > 1 leaves the solve reading $(mtheta)×$(nzeta) points " * - "as a full torus of $(mtheta)×$(nzeta * inputs.nfp)." - ) - dθ = 2π / mtheta dζ = 2π / nzeta θ_grid = range(; start=0, length=mtheta, step=dθ) @@ -680,28 +747,13 @@ function WallGeometry3D(inputs::VacuumInput, plasma_surf::PlasmaGeometry3D, wall dr_dζ = zeros(num_points, 3) normal_orient = 1 - if nowall - @info "Using no wall" - return WallGeometry3D( - nowall, - is_closed_toroidal, - mtheta, - nzeta, - r, - dr_dθ, - dr_dζ, - normal, - normal_orient - ) - end - # Plasma-wall separation at each grid point, measured the way each branch makes exact gap = zeros(num_points) if inputs.nzeta_in == 1 # Axisymmetric boundary: build the 2D poloidal contour and revolve it toroidally surf_2D = PlasmaGeometry(inputs) - wall_2D = WallGeometry(inputs, surf_2D, wall_settings) + wall_2D = WallGeometry(surf_2D, wall_settings) for i in 1:mtheta, (j, ϕ) in enumerate(ϕ_grid) r[i+mtheta*(j-1), :] .= [wall_2D.x[i] * cos(ϕ), wall_2D.x[i] * sin(ϕ), wall_2D.z[i]] end @@ -713,12 +765,12 @@ function WallGeometry3D(inputs::VacuumInput, plasma_surf::PlasmaGeometry3D, wall Z_p = plasma_surf.r[idx, 3] gap[idx] = minimum(hypot(R_p - wall_2D.x[i], Z_p - wall_2D.z[i]) for i in 1:mtheta) end - elseif wall_settings.shape == "conformal" + elseif wall_settings isa ConformalWall # Displace every plasma point outward along its own normal wall_settings.equal_arc_wall && @info "equal_arc_wall is ignored for non-axisymmetric (nzeta_in > 1) walls: it re-parameterizes a 2D contour and would break the plasma/wall index alignment the near-field patch relies on." - # Same logic as 2D, with the extrema taken over the whole torus (see WallShapeSettings.a) + # Same logic as 2D, with the extrema taken over the whole torus (see ConformalWall.a) R_plasma = [hypot(plasma_surf.r[idx, 1], plasma_surf.r[idx, 2]) for idx in axes(plasma_surf.r, 1)] offset_gap = wall_settings.a * 0.5 * (maximum(R_plasma) - minimum(R_plasma)) @info "Calculating conformal wall shape $((@sprintf "%.2e" offset_gap)) m from plasma surface." @@ -734,7 +786,7 @@ function WallGeometry3D(inputs::VacuumInput, plasma_surf::PlasmaGeometry3D, wall # Every point moves exactly offset_gap along its own normal, so the separation needs no search fill!(gap, offset_gap) else - error("Wall shape $(wall_settings.shape) is not available for a non-axisymmetric boundary (nzeta_in > 1).") + error("$(nameof(typeof(wall_settings))) is not available for a non-axisymmetric boundary (nzeta_in > 1); only ConformalWall is.") end # Compute tangent vectors and normal vectors via periodic bicubic splines @@ -755,7 +807,7 @@ function WallGeometry3D(inputs::VacuumInput, plasma_surf::PlasmaGeometry3D, wall @views normal .*= normal_orient # Fold check needs a pointwise normal offset (same-index pair) - equal_arc_wall re-parameterizes and breaks that pairing - if wall_settings.shape == "conformal" && (inputs.nzeta_in > 1 || !wall_settings.equal_arc_wall) + if wall_settings isa ConformalWall && (inputs.nzeta_in > 1 || !wall_settings.equal_arc_wall) # The offset stays regular while the normals stay aligned and the area element has not collapsed; both fail past the local concave radius of curvature min_align = Inf # min n̂_wall · n̂_plasma_outward; +1 healthy, ≤0 folded min_area_ratio = Inf # min ||n_wall||/||n_plasma||; →0 at a caustic @@ -795,7 +847,6 @@ function WallGeometry3D(inputs::VacuumInput, plasma_surf::PlasmaGeometry3D, wall end return WallGeometry3D( - nowall, is_closed_toroidal, mtheta, nzeta, diff --git a/src/Vacuum/Symmetry3D.jl b/src/Vacuum/Symmetry3D.jl index ffdd6a5e2..27687d7aa 100644 --- a/src/Vacuum/Symmetry3D.jl +++ b/src/Vacuum/Symmetry3D.jl @@ -86,11 +86,11 @@ end """ stell_sym_map(plasma, wall, nfp) -> Union{Vector{Int},Nothing} -Within-period map of `(θ, ζ) → (−θ, −ζ)` if both surfaces are stellarator symmetric to 1e-10 of the largest +Within-period map of `(θ, ζ) → (−θ, −ζ)` if the plasma and the wall (when there is one) are stellarator symmetric to 1e-10 of the largest coordinate, `nothing` otherwise. A nearly symmetric boundary (e.g. one extracted from an equilibrium, off by ~1e-3) is rejected and takes the general path; this test is the only guard on the symmetric path. """ -function stell_sym_map(plasma::PlasmaGeometry3D, wall::WallGeometry3D, nfp::Int) +function stell_sym_map(plasma::PlasmaGeometry3D, wall::Union{Nothing,WallGeometry3D}, nfp::Int) # This must be called after expand_field_periods, so nzeta_full is the full-torus nzeta mtheta, nzeta_full = plasma.mtheta, plasma.nzeta @assert nzeta_full % nfp == 0 "full-torus nzeta ($nzeta_full) must be divisible by nfp ($nfp)" @@ -105,7 +105,7 @@ function stell_sym_map(plasma::PlasmaGeometry3D, wall::WallGeometry3D, nfp::Int) # Check if the boundary(s) are stellarator symmetric # return the mapping if so, otherwise return nothing - if is_symmetric(plasma.r) && (wall.nowall || is_symmetric(wall.r)) + if is_symmetric(plasma.r) && (isnothing(wall) || is_symmetric(wall.r)) nzeta = nzeta_full ÷ nfp return [reflect_index(p, mtheta, nzeta) for p in 1:(mtheta*nzeta)] else diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index b0e008575..5d7ff0797 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -18,7 +18,7 @@ include("PnQuadCache.jl") include("Kernel2D.jl") include("Kernel3D.jl") -export VacuumInput, VacuumResponse, WallShapeSettings +export VacuumInput, VacuumResponse, WallSettings export compute_vacuum_response, compute_vacuum_response! export extract_plasma_surface_at_psi export PlasmaGeometry @@ -50,7 +50,7 @@ function _warn_and_symmetrize!(mat::AbstractMatrix, name::String) end """ - _compute_vacuum_response_2d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) + _compute_vacuum_response_2d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::Union{Nothing,WallSettings}; compute_Iv=false) 2D (axisymmetric) vacuum response calculation. @@ -59,7 +59,7 @@ building the double-/single-layer operators, solving the exterior system for `wv optionally the interior system to build `I_v` when `compute_Iv=true`. Green's functions are internal scratch only. """ -@with_pool pool function _compute_vacuum_response_2d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) +@with_pool pool function _compute_vacuum_response_2d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::Union{Nothing,WallSettings}; compute_Iv::Bool=false) mpert = length(inputs.m_modes) mlow = inputs.m_modes[1] @@ -70,7 +70,7 @@ Green's functions are internal scratch only. # Form the plasma and wall geometries plasma_surf = PlasmaGeometry(inputs) - wall = WallGeometry(inputs, plasma_surf, wall_settings) + wall = isnothing(wall_settings) ? nothing : WallGeometry(plasma_surf, wall_settings) # Loop over all decoupled toroidal modes for (idx_n, n) in enumerate(inputs.n_modes) @@ -82,7 +82,7 @@ Green's functions are internal scratch only. wv_block = @view vac_data.wv[block_idx, block_idx] # Active rows for computation (plasma only if no wall, plasma+wall if wall present) - num_points_total = wall.nowall ? num_points_surf : 2 * num_points_surf + num_points_total = isnothing(wall) ? num_points_surf : 2 * num_points_surf # Local work matrices grad_green = zeros!(pool, num_points_total, num_points_total) @@ -95,7 +95,7 @@ Green's functions are internal scratch only. # Project plasma observer onto source basis exp(i*(mθ - nζ)) mul!(view(grre, 1:num_points_surf, :), green_temp, ft.basis') - if !wall.nowall + if !isnothing(wall) # Plasma–Wall block compute_2D_kernel_matrices!(grad_green, green_temp, plasma_surf, wall, n) # Wall–Wall block @@ -148,9 +148,9 @@ Green's functions are internal scratch only. vac_data.plasma_pts[:, 1] .= plasma_surf.x vac_data.plasma_pts[:, 2] .= 0.0 vac_data.plasma_pts[:, 3] .= plasma_surf.z - vac_data.wall_pts[:, 1] .= wall.x + vac_data.wall_pts[:, 1] .= isnothing(wall) ? 0.0 : wall.x vac_data.wall_pts[:, 2] .= 0.0 - vac_data.wall_pts[:, 3] .= wall.z + vac_data.wall_pts[:, 3] .= isnothing(wall) ? 0.0 : wall.z end end @@ -226,7 +226,7 @@ Real `D_d`, `S_d` give `D̂₋ₖ = conj(D̂ₖ)`, so class `mod(nfp - k, nfp)` a conjugated mode basis. A stellarator-symmetric boundary gets a [`StellSymBasis`](@ref), making the class operator real and splitting a self-conjugate class into two half-size blocks. """ -@with_pool pool function _compute_vacuum_response_3d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) +@with_pool pool function _compute_vacuum_response_3d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::Union{Nothing,WallSettings}; compute_Iv::Bool=false) (; mtheta, nzeta, nfp, m_modes, n_modes) = inputs fill!(vac_data.wv, 0) @@ -234,12 +234,18 @@ class operator real and splitting a self-conjugate class into two half-size bloc # Full-torus geometry for sources; observers are one field period full = expand_field_periods(inputs) + # Expansion leaves nfp > 1 only on an axisymmetric boundary, which has no field periodicity + full.nfp > 1 && error( + "3D vacuum requires nfp = 1 on an axisymmetric (nzeta_in = 1) boundary: an axisymmetric " * + "surface has no field periodicity, and nfp > 1 leaves the solve reading $(mtheta)×$(nzeta) points " * + "as a full torus of $(mtheta)×$(nzeta * nfp)." + ) plasma_surf = PlasmaGeometry3D(full) - wall = WallGeometry3D(full, plasma_surf, wall_settings) + wall = isnothing(wall_settings) ? nothing : WallGeometry3D(full, plasma_surf, wall_settings) num_points_per_fp = mtheta * nzeta mpert = length(m_modes) - nb = wall.nowall ? 1 : 2 # plasma, or [plasma; wall] + nb = isnothing(wall) ? 1 : 2 # plasma, or [plasma; wall] n_obs = nb * num_points_per_fp num_modes = mpert * length(n_modes) @@ -298,7 +304,7 @@ class operator real and splitting a self-conjugate class into two half-size bloc # Plasma–Plasma block compute_3D_kernel_matrices!(grad_blocks, green_blocks, plasma_surf, plasma_surf, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases, sym_basis) - if !wall.nowall + if !isnothing(wall) # Plasma–Wall block compute_3D_kernel_matrices!(grad_blocks, green_blocks, plasma_surf, wall, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases, sym_basis) # Wall–Wall block @@ -388,27 +394,28 @@ class operator real and splitting a self-conjugate class into two half-size bloc # Populate coordinate arrays vac_data.plasma_pts .= plasma_surf.r - vac_data.wall_pts .= wall.r + vac_data.wall_pts .= isnothing(wall) ? 0.0 : wall.r end """ - compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) -> VacuumResponse + compute_vacuum_response(inputs::VacuumInput, wall_settings::Union{Nothing,WallSettings}; compute_Iv=false) -> VacuumResponse Allocating version of [`compute_vacuum_response!`](@ref). """ -function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) +function compute_vacuum_response(inputs::VacuumInput, wall_settings::Union{Nothing,WallSettings}; compute_Iv::Bool=false) vac = VacuumResponse(inputs) compute_vacuum_response!(vac, inputs, wall_settings; compute_Iv) return vac end """ - compute_vacuum_response!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) + compute_vacuum_response!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::Union{Nothing,WallSettings}; compute_Iv=false) Compute the vacuum response into `vac_data`. 2D (`inputs.nzeta == 1`) and 3D otherwise. Pass -`compute_Iv=true` to also fill the surface-current matrix `I_v`. +`compute_Iv=true` to also fill the surface-current matrix `I_v`. `wall_settings` is a +[`WallSettings`](@ref), or `nothing` for no wall (the vacuum region extends to infinity). """ -function compute_vacuum_response!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) +function compute_vacuum_response!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::Union{Nothing,WallSettings}; compute_Iv::Bool=false) if inputs.nzeta == 1 _compute_vacuum_response_2d!(vac_data, inputs, wall_settings; compute_Iv) else diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index 98fca1a63..c5da1a8fb 100644 --- a/test/runtests_eulerlagrange.jl +++ b/test/runtests_eulerlagrange.jl @@ -562,7 +562,7 @@ end eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], example_dir) sol_cfg = haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, sol_cfg) - intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in inputs["Wall"])...) + intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.parse_wall(get(inputs, "Wall", nothing)) FFS.sing_lim!(intr, ctrl, equil) intr.nlow = ctrl.nn_low intr.nhigh = ctrl.nn_high diff --git a/test/runtests_parallel_integration.jl b/test/runtests_parallel_integration.jl index b427c3795..a9e6defad 100644 --- a/test/runtests_parallel_integration.jl +++ b/test/runtests_parallel_integration.jl @@ -114,7 +114,10 @@ using TOML ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex) - equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing) + equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium( + eq_config, + haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing + ) intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; intr.npert = 1 @@ -201,7 +204,10 @@ using TOML ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex) - equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing) + equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium( + eq_config, + haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing + ) intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; intr.npert = 1 @@ -262,9 +268,11 @@ using TOML ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex) - equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing) - intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(; - (Symbol(k) => v for (k, v) in inputs["Wall"])...) + equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium( + eq_config, + haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing + ) + intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.parse_wall(get(inputs, "Wall", nothing)) intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; intr.npert = 1 @@ -311,7 +319,10 @@ using TOML ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex) - equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing) + equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium( + eq_config, + haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing + ) # Apply the two-pass auto grid exactly as the main driver does (the example ships # grid_type="auto", mpsi=0): measured-curvature refinement with rational surfaces # pinned as mandatory knots, re-formed from the captured ingest. The pinned values @@ -322,8 +333,7 @@ using TOML rerun_input = GeneralizedPerturbedEquilibrium.Equilibrium.build_direct_from_ingest(eq_config, equil.ingest) equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, rerun_input; override_psi_nodes=psi_nodes) end - intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(; - (Symbol(k) => v for (k, v) in inputs["Wall"])...) + intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.parse_wall(get(inputs, "Wall", nothing)) intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; intr.npert = 1 @@ -372,7 +382,10 @@ using TOML ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex) - equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing) + equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium( + eq_config, + haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing + ) intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; intr.npert = 1 @@ -423,9 +436,11 @@ using TOML ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], example_dir) - equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing) - intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(; - (Symbol(k) => v for (k, v) in inputs["Wall"])...) + equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium( + eq_config, + haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing + ) + intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.parse_wall(get(inputs, "Wall", nothing)) intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; intr.npert = 1 @@ -476,7 +491,10 @@ using TOML ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex) - equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing) + equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium( + eq_config, + haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing + ) # Apply the two-pass auto grid (measured-curvature refinement, rational surfaces pinned # as mandatory knots) exactly as the main driver does (see the FM testset above); the # pinned values below are for this grid, the production default. @@ -486,8 +504,7 @@ using TOML rerun_input = GeneralizedPerturbedEquilibrium.Equilibrium.build_direct_from_ingest(eq_config, equil.ingest) equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, rerun_input; override_psi_nodes=psi_nodes) end - intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(; - (Symbol(k) => v for (k, v) in inputs["Wall"])...) + intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.parse_wall(get(inputs, "Wall", nothing)) intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; intr.npert = 1 diff --git a/test/runtests_riccati.jl b/test/runtests_riccati.jl index e5989aed8..bebe399e3 100644 --- a/test/runtests_riccati.jl +++ b/test/runtests_riccati.jl @@ -7,8 +7,7 @@ const FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates # setup because intr is mutated by each integration (sing[s].delta_prime etc.). function make_solovev_intr(inputs, ctrl, equil, ex) intr = FFS.ForceFreeStatesInternal(; dir_path=ex) - intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(; - (Symbol(k) => v for (k, v) in inputs["Wall"])...) + intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.parse_wall(get(inputs, "Wall", nothing)) intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; intr.npert = 1 diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index d04eb9065..003f5cc8a 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -14,12 +14,12 @@ using TOML GPEC.Equilibrium.EquilibriumConfig(deck["Equilibrium"], template), GPEC.Equilibrium.SolovevConfig(deck["SOL_INPUT"]) ) - wall = GPEC.Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in deck["Wall"])...) + wall = GPEC.Vacuum.parse_wall(get(deck, "Wall", nothing)) # The deck's `[ForceFreeStates]` block as `solve` keywords: everything except the knobs the # integrator object and the `nn` keyword own. ffs_kwargs = Dict(Symbol(k) => v for (k, v) in deck["ForceFreeStates"] - if !(k in ("integrator", "nn_low", "nn_high"))) + if !(k in ("integrator", "nn_low", "nn_high"))) # Copy the fixture deck into `dir` and apply `[ForceFreeStates]` overrides / extra sections. function _stage_deck(dir; ffs_overrides=Dict{String,Any}(), extra_sections=Dict{String,Any}()) @@ -160,7 +160,7 @@ using TOML @test d.terms[2].scale == 1.0 + 0.0im @test (im * a).scale == im @test (a * 3).scale == 3.0 + 0.0im - @test (a - b).terms[2].scale == -0.5 + 0.0im + @test (a-b).terms[2].scale == -0.5 + 0.0im @test (-a).scale == -1.0 + 0.0im end diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl index 68b9728d4..1f79950e9 100644 --- a/test/runtests_vacuum.jl +++ b/test/runtests_vacuum.jl @@ -1,3 +1,5 @@ +import GeneralizedPerturbedEquilibrium.Vacuum: ConformalWall, EllipticalWall, DeeWall, ModDeeWall, FileWall, parse_wall + @testset "Vacuum.jl Unit Tests" begin @testset "Vacuum.jl (2D)" begin @@ -28,20 +30,33 @@ end end - @testset "WallShapeSettings" begin + @testset "WallSettings" begin @testset "default constructor" begin - w = WallShapeSettings() - @test w.shape == "nowall" + w = ConformalWall() @test w.a == 0.3 @test w.equal_arc_wall == true end @testset "keyword constructor" begin - w = WallShapeSettings(; shape="conformal", a=0.2, equal_arc_wall=false) - @test w.shape == "conformal" + w = ConformalWall(a=0.2, equal_arc_wall=false) @test w.a == 0.2 @test w.equal_arc_wall == false end + + @testset "parse_wall" begin + # No [Wall] section means no wall + @test parse_wall(nothing) === nothing + @test parse_wall(Dict("shape" => "conformal", "a" => 0.2)) == ConformalWall(a=0.2) + @test parse_wall(Dict("shape" => "dee", "a" => 1, "cw" => 0)) == DeeWall(a=1.0, cw=0.0) + @test parse_wall(Dict("shape" => "file", "file" => "wall.dat")) == FileWall(file="wall.dat") + # A section needs a shape, and "nowall" is spelled by omitting the section + @test_throws "needs a `shape`" parse_wall(Dict("a" => 0.2)) + @test_throws "omit the [Wall] section" parse_wall(Dict("shape" => "nowall")) + # Keys the shape does not use are ignored with a warning; a required one that is missing is an error + w = @test_logs (:warn, r"does not use aw") parse_wall(Dict("shape" => "conformal", "a" => 0.2, "aw" => 0.05)) + @test w == ConformalWall(a=0.2) + @test_throws "requires `file`" parse_wall(Dict("shape" => "file")) + end end @testset "PlasmaGeometry" begin @@ -97,22 +112,11 @@ nzeta=1 ) - @testset "nowall" begin - inputs = _circle_inputs(16) - plasma_surf = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry(inputs) - wall_settings = WallShapeSettings(; shape="nowall") - wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry(inputs, plasma_surf, wall_settings) - @test wall.nowall == true - @test length(wall.x) == 16 - @test length(wall.z) == 16 - end - @testset "conformal" begin inputs = _circle_inputs(16) plasma_surf = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry(inputs) - wall_settings = WallShapeSettings(; shape="conformal", a=0.2) - wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry(inputs, plasma_surf, wall_settings) - @test wall.nowall == false + wall_settings = ConformalWall(a=0.2) + wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry(plasma_surf, wall_settings) @test length(wall.x) == 16 @test length(wall.z) == 16 @test all(isfinite, wall.x) @@ -125,9 +129,8 @@ @testset "elliptical" begin inputs = _circle_inputs(16) plasma_surf = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry(inputs) - wall_settings = WallShapeSettings(; shape="elliptical", a=0.5) - wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry(inputs, plasma_surf, wall_settings) - @test wall.nowall == false + wall_settings = EllipticalWall(a=0.5) + wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry(plasma_surf, wall_settings) @test length(wall.x) == 16 @test all(isfinite, wall.x) @test all(isfinite, wall.z) @@ -136,9 +139,8 @@ @testset "dee" begin inputs = _circle_inputs(16) plasma_surf = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry(inputs) - wall_settings = WallShapeSettings(; shape="dee", a=0.1, cw=0.0) - wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry(inputs, plasma_surf, wall_settings) - @test wall.nowall == false + wall_settings = DeeWall(a=0.1, cw=0.0) + wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry(plasma_surf, wall_settings) @test length(wall.x) == 16 @test all(wall.x .> 0) end @@ -159,8 +161,8 @@ # Use a "dee" wall shape with parameters that will produce R < 0 # Setting cw (offset) to a large negative value will shift the wall left past R=0 - wall_settings = WallShapeSettings(; shape="dee", cw=-1.5, a=0.1) - @test_throws ErrorException GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry(inputs, plasma_surf_near_zero, wall_settings) + wall_settings = DeeWall(cw=-1.5, a=0.1) + @test_throws ErrorException GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry(plasma_surf_near_zero, wall_settings) end # Test that conformal wall R-coordinates are clamped by centerstack_min @@ -179,8 +181,8 @@ nzeta=1 ) ) - wall_settings = WallShapeSettings(; shape="conformal", a=10.0, equal_arc_wall=false) - wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry(inputs, plasma_surf_near_zero, wall_settings) + wall_settings = ConformalWall(a=10.0, equal_arc_wall=false) + wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry(plasma_surf_near_zero, wall_settings) expected_min = min(0.1, 0.1 * minimum(plasma_surf_near_zero.x)) @test all(wall.x .>= expected_min) @test any(wall.x .<= expected_min + 1e-10) @@ -360,9 +362,9 @@ mtheta=mtheta ) - @testset "nowall" begin + @testset "no wall" begin inputs = _make_inputs() - wall_settings = WallShapeSettings(; shape="nowall") + wall_settings = nothing (; wv, I_v, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta @@ -376,14 +378,15 @@ @test size(plasma_pts) == (numpoints, 3) @test all(isfinite, plasma_pts) @test size(wall_pts) == (numpoints, 3) + @test all(iszero, wall_pts) # no wall to report # wv is always Hermitian-symmetrized @test isapprox(wv, wv', rtol=1e-12) end - @testset "nowall compute_Iv=true" begin + @testset "no wall compute_Iv=true" begin inputs = _make_inputs() - wall_settings = WallShapeSettings(; shape="nowall") + wall_settings = nothing (; wv, I_v, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings; compute_Iv=true) num_modes = length(inputs.m_modes) * length(inputs.n_modes) @@ -397,7 +400,7 @@ @testset "conformal wall" begin inputs = _make_inputs() - wall_settings = WallShapeSettings(; shape="conformal", a=0.5) + wall_settings = ConformalWall(a=0.5) (; wv, I_v, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings) num_modes = length(inputs.m_modes) * length(inputs.n_modes) @@ -413,7 +416,7 @@ @testset "edge: single poloidal mode mpert=1" begin inputs = _make_inputs(; m_modes=[1], n_modes=[1]) - wall_settings = WallShapeSettings(; shape="nowall") + wall_settings = nothing (; wv, I_v, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings) @test size(wv) == (1, 1) @test all(isfinite, wv) @@ -423,7 +426,7 @@ @testset "edge: small mtheta" begin # Keep mtheta_eq=17 so boundary has enough points for periodic spline inputs = _make_inputs(; mtheta=16, mtheta_eq=17) - wall_settings = WallShapeSettings(; shape="nowall") + wall_settings = nothing (; wv, I_v, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings) @test size(wv) == (2, 2) @test size(I_v) == (2, 2) @@ -433,7 +436,7 @@ @testset "in-place compute_vacuum_response! matches wrapper" begin # The allocating wrapper is a thin caller of the in-place routine; verify the # in-place entry populates caller-owned storage identically. - for wall_settings in (WallShapeSettings(; shape="nowall"), WallShapeSettings(; shape="conformal", a=0.5)) + for wall_settings in (nothing, ConformalWall(a=0.5)) inputs = _make_inputs() ref = compute_vacuum_response(inputs, wall_settings; compute_Iv=true) @@ -453,11 +456,11 @@ inputs = _make_inputs() vac = VacuumResponse(inputs) - compute_vacuum_response!(vac, inputs, WallShapeSettings(; shape="conformal", a=0.5); compute_Iv=true) + compute_vacuum_response!(vac, inputs, ConformalWall(a=0.5); compute_Iv=true) @test !all(iszero, vac.I_v) - compute_vacuum_response!(vac, inputs, WallShapeSettings(; shape="nowall")) - fresh = compute_vacuum_response(inputs, WallShapeSettings(; shape="nowall")) + compute_vacuum_response!(vac, inputs, nothing) + fresh = compute_vacuum_response(inputs, nothing) @test all(iszero, vac.I_v) @test vac.wv ≈ fresh.wv @@ -607,24 +610,13 @@ @test maximum(abs, surf.r[:, 2]) > 0 end - @testset "WallGeometry3D nowall" begin - inputs = _make_3d_inputs(; mtheta=32, nzeta=32, mtheta_eq=17) - wall = - GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(inputs, GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(inputs), WallShapeSettings(; shape="nowall")) - @test wall.nowall == true - @test wall.mtheta == 32 - @test wall.nzeta == 32 - @test size(wall.r) == (32 * 32, 3) - end - @testset "WallGeometry3D conformal" begin inputs = _make_3d_inputs(; mtheta=32, nzeta=32, mtheta_eq=17) wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D( inputs, GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(inputs), - WallShapeSettings(; shape="conformal", a=0.2) + ConformalWall(a=0.2) ) - @test wall.nowall == false @test wall.mtheta == 32 @test wall.nzeta == 32 num_points = 32 * 32 @@ -635,7 +627,7 @@ end # Corrugated torus on an endpoint-excluded grid. `_make_3d_nonaxis_inputs` repeats the seam point, which is - # harmless for a nowall response but leaves a kink that folds the offset surface, so wall tests need their own. + # harmless for a no-wall response but leaves a kink that folds the offset surface, so wall tests need their own. _make_3d_periodic_inputs(; mtheta=24, nzeta=24, mtheta_in=16, nzeta_in=16, eps=0.05) = begin θ_in = range(; start=0, length=mtheta_in, step=2π / mtheta_in) ζ_in = range(; start=0, length=nzeta_in, step=2π / nzeta_in) @@ -654,12 +646,11 @@ @testset "WallGeometry3D conformal, non-axisymmetric input" begin inputs = _make_3d_periodic_inputs(; mtheta=24, nzeta=24) - settings = WallShapeSettings(; shape="conformal", a=0.2, equal_arc_wall=false) + settings = ConformalWall(a=0.2, equal_arc_wall=false) plasma = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(inputs) wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(inputs, plasma, settings) num_points = 24 * 24 - @test wall.nowall == false @test size(wall.r) == (num_points, 3) @test all(isfinite, wall.r) && all(isfinite, wall.normal) @@ -690,7 +681,7 @@ inputs = _make_3d_periodic_inputs(; mtheta=24, nzeta=24, mtheta_in=32, nzeta_in=16, eps=0.0) plasma = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(inputs) for a in (0.2, 0.5) - wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(inputs, plasma, WallShapeSettings(; shape="conformal", a=a, equal_arc_wall=false)) + wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(inputs, plasma, ConformalWall(a=a, equal_arc_wall=false)) radii = [hypot(hypot(wall.r[i, 1], wall.r[i, 2]) - R0, wall.r[i, 3]) for i in axes(wall.r, 1)] # Residual is the 32 -> 24 poloidal resampling of the boundary spline, not the offset @test all(isapprox(rad, r_minor * (1 + a); atol=5e-4) for rad in radii) @@ -700,35 +691,36 @@ @testset "WallGeometry3D non-axisymmetric error paths" begin inputs = _make_3d_periodic_inputs(; mtheta=24, nzeta=24) plasma = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(inputs) - for shape in ("elliptical", "dee", "mod_dee", "some_wall_file.dat") - @test_throws ErrorException GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(inputs, plasma, WallShapeSettings(; shape=shape)) + for wall in (EllipticalWall(), DeeWall(), ModDeeWall(), FileWall(file="some_wall_file.dat")) + @test_throws ErrorException GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(inputs, plasma, wall) end # equal_arc_wall has no meaning without a 2D contour and is ignored with a notice. It defaults to # true, so this fires on runs that never set it: @info, not @warn. @test_logs (:info, r"equal_arc_wall is ignored") match_mode = :any GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D( inputs, plasma, - WallShapeSettings(; shape="conformal", a=0.2, equal_arc_wall=true) + ConformalWall(a=0.2, equal_arc_wall=true) ) - # nfp > 1 can only reach the constructor from an axisymmetric boundary, and is rejected there - per_period = VacuumInput(; x=inputs.x, y=inputs.y, z=inputs.z, mtheta_in=16, nzeta_in=16, - m_modes=inputs.m_modes, n_modes=inputs.n_modes, mtheta=24, nzeta=24, nfp=3) - @test_throws ErrorException GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(per_period, plasma, WallShapeSettings(; shape="conformal", a=0.2)) - # nowall does not skip that guard - @test_throws ErrorException GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(per_period, plasma, WallShapeSettings(; shape="nowall")) + # nfp > 1 on an axisymmetric boundary has no field periodicity to declare, with or without a wall + axi = _make_3d_inputs(; mtheta=32, nzeta=32, mtheta_eq=17) + axi_nfp = VacuumInput(; x=axi.x, z=axi.z, ν=axi.ν, mtheta_in=axi.mtheta_in, nzeta_in=1, + m_modes=axi.m_modes, n_modes=axi.n_modes, mtheta=axi.mtheta, nzeta=axi.nzeta, nfp=3) + for wall in (nothing, ConformalWall(a=0.2)) + @test_throws "requires nfp = 1" compute_vacuum_response(axi_nfp, wall) + end # a plasma surface from a different grid is rejected rather than indexed off the end other_grid = _make_3d_periodic_inputs(; mtheta=16, nzeta=16) @test_throws ErrorException GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D( inputs, GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(other_grid), - WallShapeSettings(; shape="conformal", a=0.2) + ConformalWall(a=0.2) ) # an offset that folds the surface is rejected rather than silently returned folded = _make_3d_nonaxis_inputs(; mtheta=24, nzeta=24, mtheta_in=12, nzeta_in=12) @test_throws ErrorException GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D( folded, GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(folded), - WallShapeSettings(; shape="conformal", a=0.2, equal_arc_wall=false) + ConformalWall(a=0.2, equal_arc_wall=false) ) end @@ -739,9 +731,9 @@ axi = VacuumInput(; mtheta_in=32, nzeta_in=1, x=collect(1.7 .+ 0.3 .* cos.(θ_eq)), z=collect(0.45 .* sin.(θ_eq)), ν=zeros(32), m_modes=[1, 2], n_modes=[0, 1], mtheta=32, nzeta=32) - for (shape, settings) in (("elliptical", WallShapeSettings(; shape="elliptical", a=0.5)), - ("dee", WallShapeSettings(; shape="dee", a=0.3)), - ("mod_dee", WallShapeSettings(; shape="mod_dee", a=0.5, cw=1.7))) + for (shape, settings) in (("elliptical", EllipticalWall(a=0.5)), + ("dee", DeeWall(a=0.3)), + ("mod_dee", ModDeeWall(a=0.5, cw=1.7))) wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(axi, GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(axi), settings) @test size(wall.r) == (32 * 32, 3) @test all(isfinite, wall.normal) @@ -752,19 +744,19 @@ # Gap warning uses the measured same-index separation, not a conformal-only offset that # is 0 for these shapes. a = 1 m sits well outside one coarser cell, so no sub-cell warning. @test_logs min_level = Base.CoreLogging.Warn GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D( - axi, GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(axi), WallShapeSettings(; shape="elliptical", a=1.0) + axi, GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(axi), EllipticalWall(a=1.0) ) # an unreadable wall file is still reported by the 2D reader @test_throws ErrorException GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D( axi, GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(axi), - WallShapeSettings(; shape="no_such_wall.dat") + FileWall(file="no_such_wall.dat") ) end - @testset "compute_vacuum_response 3D nowall" begin + @testset "compute_vacuum_response 3D no wall" begin inputs = _make_3d_inputs(; mtheta=32, nzeta=32, mtheta_eq=17) - wall_settings = WallShapeSettings(; shape="nowall") + wall_settings = nothing (; wv, I_v, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta @@ -777,6 +769,7 @@ @test size(plasma_pts) == (numpoints, 3) @test all(isfinite, plasma_pts) @test size(wall_pts) == (numpoints, 3) + @test all(iszero, wall_pts) # no wall to report # 3D plasma_pts are (X,Y,Z) Cartesian @test isapprox(plasma_pts[1, 1]^2 + plasma_pts[1, 2]^2, (1.7 + 0.3)^2, rtol=0.1) @test isapprox(plasma_pts[1, 3], 0.0, atol=0.1) @@ -785,7 +778,7 @@ @testset "compute_vacuum_response 3D nonaxisymmetric boundary" begin inputs = _make_3d_nonaxis_inputs(; mtheta=24, nzeta=24, mtheta_in=12, nzeta_in=12, mpert=2, nlow=0, npert=2) - wall_settings = WallShapeSettings(; shape="nowall") + wall_settings = nothing (; wv, I_v, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta @@ -804,7 +797,7 @@ @testset "compute_vacuum_response 3D conformal wall" begin inputs = _make_3d_inputs(; mtheta=32, nzeta=32, mtheta_eq=17) - wall_settings = WallShapeSettings(; shape="conformal", a=0.3) + wall_settings = ConformalWall(a=0.3) (; wv, I_v, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings) num_modes = length(inputs.m_modes) * length(inputs.n_modes) @@ -820,7 +813,7 @@ @testset "compute_vacuum_response 3D compute_Iv=true" begin inputs = _make_3d_inputs(; mtheta=32, nzeta=32, mtheta_eq=17) # Both wall branches: I_v was left zeroed in 3D before, and the wall adds a second source block - for wall_settings in (WallShapeSettings(; shape="nowall"), WallShapeSettings(; shape="conformal", a=0.3)) + for wall_settings in (nothing, ConformalWall(a=0.3)) (; I_v) = compute_vacuum_response(inputs, wall_settings; compute_Iv=true) @test all(isfinite, I_v) @test !all(iszero, I_v) @@ -828,9 +821,9 @@ # A reused buffer must not keep a stale I_v from an earlier compute_Iv=true solve vac = GeneralizedPerturbedEquilibrium.Vacuum.VacuumResponse(inputs) - compute_vacuum_response!(vac, inputs, WallShapeSettings(; shape="nowall"); compute_Iv=true) + compute_vacuum_response!(vac, inputs, nothing; compute_Iv=true) @test !all(iszero, vac.I_v) - compute_vacuum_response!(vac, inputs, WallShapeSettings(; shape="nowall")) + compute_vacuum_response!(vac, inputs, nothing) @test all(iszero, vac.I_v) end @@ -845,7 +838,7 @@ Z = 0.3 .* sin.(θ) .+ 0.08 .* sin.(2θ) .+ 0.08 .* cos.(θ) make(nzeta) = VacuumInput(; x=collect(R), z=collect(Z), ν=zeros(mtheta), mtheta_in=mtheta, nzeta_in=1, m_modes=[-2, -1, 0, 1, 2], n_modes=[1], mtheta=mtheta, nzeta=nzeta) - nowall = WallShapeSettings(; shape="nowall") + nowall = nothing r2d = compute_vacuum_response(make(1), nowall; compute_Iv=true) r3d = compute_vacuum_response(make(96), nowall; compute_Iv=true) @@ -866,7 +859,7 @@ ν = 0.25 .* sin.(θ) .+ 0.1 .* cos.(2θ) .+ 0.05 make(nzeta, ν) = VacuumInput(x=collect(R), z=collect(Z), ν=collect(ν), mtheta_in=mtheta, nzeta_in=1, m_modes=collect(-3:3), n_modes=[1], mtheta=mtheta, nzeta=nzeta) - nowall = WallShapeSettings(shape="nowall") + nowall = nothing r3d = compute_vacuum_response(make(48, ν), nowall) @test norm(compute_vacuum_response(make(1, ν), nowall).wv - r3d.wv) / norm(r3d.wv) < 1e-2 @@ -882,7 +875,7 @@ Z = 0.8 .* sin.(θ) .+ 0.08 .* cos.(θ) .+ 0.1 .* cos.(2θ) m_modes = collect(-4:6) D = Diagonal(cis.(m_modes .* (2π * k / mtheta))) - nowall = WallShapeSettings(shape="nowall") + nowall = nothing # 2D only: "3D I_v matches the 2D path" ties the 3D frame to this one make(x, z) = VacuumInput(; x, z, ν=zeros(mtheta), mtheta_in=mtheta, nzeta_in=1, m_modes, n_modes=[1], mtheta, nzeta=1) ref = compute_vacuum_response(make(collect(R), collect(Z)), nowall; compute_Iv=true) @@ -898,7 +891,7 @@ θ = range(; start=0, length=mtheta, step=2π/mtheta) R, Z = 1.7 .+ 0.3 .* cos.(θ), 0.4 .* sin.(θ) make(x, z) = VacuumInput(; x, z, ν=zeros(mtheta), mtheta_in=mtheta, nzeta_in=1, m_modes=collect(-2:2), n_modes=[1], mtheta, nzeta=1) - @test_throws ArgumentError compute_vacuum_response(make(reverse(R), reverse(Z)), WallShapeSettings(shape="nowall")) + @test_throws ArgumentError compute_vacuum_response(make(reverse(R), reverse(Z)), nothing) # A wall file written clockwise is turned round and gives the same operator as its counter-clockwise twin Rw, Zw = 1.7 .+ 0.6 .* cos.(θ), 0.8 .* sin.(θ) @@ -908,7 +901,7 @@ println(io, mtheta, "\n", 1.7, "\nindex R Z") foreach(i -> println(io, i, " ", x[i], " ", z[i]), 1:mtheta) end - compute_vacuum_response(make(collect(R), collect(Z)), WallShapeSettings(shape=path, equal_arc_wall=false)) + compute_vacuum_response(make(collect(R), collect(Z)), FileWall(file=path, equal_arc_wall=false)) end ccw = wall_response(Rw, Zw) cw = @test_logs (:warn, r"runs clockwise") match_mode = :any wall_response(circshift(reverse(Rw), 1), circshift(reverse(Zw), 1)) @@ -922,15 +915,15 @@ # the near-field patch keeps its full PATCH_RAD. On a 24 x 24 grid those two demands cannot be met # by a close wall: h_ζ ~ 2π R0 / nzeta is 0.45 m here. inputs = _make_3d_periodic_inputs(; mtheta=24, nzeta=24) - settings = WallShapeSettings(; shape="conformal", a=2.5, equal_arc_wall=false) + settings = ConformalWall(a=2.5, equal_arc_wall=false) wv = compute_vacuum_response(inputs, settings).wv @test all(isfinite, wv) @test isapprox(wv, wv', rtol=1e-12) - # A conducting wall is stabilizing, so Wᵛ must exceed the nowall value. Catches a flipped normal + # A conducting wall is stabilizing, so Wᵛ must exceed the no-wall value. Catches a flipped normal # or an offset built on the wrong side, which finiteness and Hermiticity both pass. - @test real(wv[1, 1]) > real(compute_vacuum_response(inputs, WallShapeSettings(; shape="nowall")).wv[1, 1]) + @test real(wv[1, 1]) > real(compute_vacuum_response(inputs, nothing).wv[1, 1]) # This configuration clears the one-cell floor, so it draws no sub-cell or near-fold warning @test_logs min_level = Base.CoreLogging.Warn GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D( @@ -975,7 +968,7 @@ mtheta=mtheta, nzeta=nzeta_p, nfp=nfp ) - wall_settings = WallShapeSettings(; shape="nowall") + wall_settings = nothing # Reduced (block-circulant) path vac_red = compute_vacuum_response(inputs_red, wall_settings) @@ -1021,7 +1014,7 @@ # A wall adds a second source block to the operator, so repeat the check with one present: # the field-period fold has to land the wall columns in the right block for both to agree. - walled = WallShapeSettings(; shape="conformal", a=0.2, equal_arc_wall=false) + walled = ConformalWall(a=0.2, equal_arc_wall=false) vac_red_wall = compute_vacuum_response(inputs_red, walled; compute_Iv=true) vac_full_wall = compute_vacuum_response(inputs_full, walled; compute_Iv=true) @test isapprox(vac_red_wall.wv, vac_full_wall.wv; rtol=1e-6, atol=1e-7) @@ -1059,7 +1052,7 @@ end mtheta, nzeta_p = 24, 8 - walled = WallShapeSettings(; shape="conformal", a=0.2, equal_arc_wall=false) + walled = ConformalWall(a=0.2, equal_arc_wall=false) # The whole item rests on the operator inheriting the reflection symmetry, so assert it directly. inputs = _stell_inputs(; mtheta=mtheta, nzeta_p=nzeta_p, nfp=3, n_modes=[1]) @@ -1101,8 +1094,7 @@ inp = _stell_inputs(; mtheta=mtheta, nzeta_p=nzeta_p, nfp=nfp, n_modes=[k]) geom = GeneralizedPerturbedEquilibrium.Vacuum.expand_field_periods(inp) plas = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(geom) - wal = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(geom, plas, WallShapeSettings(; shape="nowall")) - mirror = GeneralizedPerturbedEquilibrium.Vacuum.stell_sym_map(plas, wal, nfp) + mirror = GeneralizedPerturbedEquilibrium.Vacuum.stell_sym_map(plas, nothing, nfp) @test mirror !== nothing nfp_pts = plas.mtheta * plas.nzeta ÷ nfp ω = ComplexF64[cis(-2π * (k * d) / nfp) for d in 0:(nfp-1)] @@ -1135,7 +1127,7 @@ # Every reduced-vs-full check above uses a symmetric boundary; this takes the general complex path, with pairing asym_red = _stell_inputs(; mtheta=mtheta, nzeta_p=nzeta_p, nfp=3, n_modes=[1, 2, 3, 4], odd=0.07) asym_ref = GeneralizedPerturbedEquilibrium.Vacuum.expand_field_periods(asym_red) - for ws in (WallShapeSettings(; shape="nowall"), walled) + for ws in (nothing, walled) red = compute_vacuum_response(asym_red, ws; compute_Iv=true) ref = compute_vacuum_response(asym_ref, ws; compute_Iv=true) @test isapprox(red.wv, ref.wv; rtol=1e-10) diff --git a/test/test_data/regression_solovev_ideal_example/gpec.toml b/test/test_data/regression_solovev_ideal_example/gpec.toml index 9b3d346ac..b4d251d2c 100644 --- a/test/test_data/regression_solovev_ideal_example/gpec.toml +++ b/test/test_data/regression_solovev_ideal_example/gpec.toml @@ -15,13 +15,8 @@ etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) [Wall] -shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath); close conformal wall stabilizes the Solovev n=1 kink (see examples/Solovev_ideal_example/gpec.toml) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) +shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file); close conformal wall stabilizes the Solovev n=1 kink (see examples/Solovev_ideal_example/gpec.toml) +a = 0.2415 # Distance from plasma in units of the minor radius equal_arc_wall = true # Equal arc length distribution of nodes on wall [ForceFreeStates] diff --git a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml index 58c4a6cd1..572c66845 100644 --- a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml +++ b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml @@ -15,13 +15,8 @@ etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) [Wall] -shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath); close conformal wall stabilizes the Solovev n=1 kink (see examples/Solovev_ideal_example/gpec.toml) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) +shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file); close conformal wall stabilizes the Solovev n=1 kink (see examples/Solovev_ideal_example/gpec.toml) +a = 0.2415 # Distance from plasma in units of the minor radius equal_arc_wall = true # Equal arc length distribution of nodes on wall [ForceFreeStates] diff --git a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml index 1a6bb6f13..7e511579d 100644 --- a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml @@ -16,13 +16,8 @@ force_termination = false # Terminate after equilibrium setup (ski [Wall] # A close conformal wall stabilizes the Solovev n=1 kink for this fixture. -shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) +shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file) +a = 0.2415 # Distance from plasma in units of the minor radius equal_arc_wall = true # Equal arc length distribution of nodes on wall [ForceFreeStates] diff --git a/test/test_data/regression_solovev_kinetic_example/gpec.toml b/test/test_data/regression_solovev_kinetic_example/gpec.toml index 67692a315..77f7b00b8 100644 --- a/test/test_data/regression_solovev_kinetic_example/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_example/gpec.toml @@ -15,13 +15,8 @@ etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) [Wall] -shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath); close conformal wall stabilizes the Solovev n=1 kink (see examples/Solovev_ideal_example/gpec.toml) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) +shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file); close conformal wall stabilizes the Solovev n=1 kink (see examples/Solovev_ideal_example/gpec.toml) +a = 0.2415 # Distance from plasma in units of the minor radius equal_arc_wall = true # Equal arc length distribution of nodes on wall [ForceFreeStates] diff --git a/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml b/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml index b2bdfb2b1..319e9c05b 100644 --- a/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml @@ -15,13 +15,8 @@ etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) [Wall] -shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath); close conformal wall stabilizes the Solovev n=1 kink (see examples/Solovev_ideal_example/gpec.toml) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) +shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file); close conformal wall stabilizes the Solovev n=1 kink (see examples/Solovev_ideal_example/gpec.toml) +a = 0.2415 # Distance from plasma in units of the minor radius equal_arc_wall = true # Equal arc length distribution of nodes on wall [ForceFreeStates] diff --git a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml index 7debb9e8c..c2046d995 100644 --- a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml @@ -16,13 +16,8 @@ etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) [Wall] -shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) -a = 0.2415 # Distance from plasma (conformal) or shape parameter -aw = 0.05 # Half-thickness parameter for Dee-shaped walls -bw = 1.5 # Elongation parameter for wall shapes -cw = 0 # Offset of wall center from major radius -dw = 0.5 # Triangularity parameter for wall shapes -tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) +shape = "conformal" # Wall shape (conformal, elliptical, dee, mod_dee, file) +a = 0.2415 # Distance from plasma in units of the minor radius equal_arc_wall = true # Equal arc length distribution of nodes on wall [ForceFreeStates]