diff --git a/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md b/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md index 9d17a0bdd4..33c9ba6f92 100644 --- a/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md +++ b/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md @@ -1,7 +1,7 @@ # Agent Memory Index - [Fortran-Julia Correspondence Map](fortran_correspondence_map.md) — Maps Fortran GPEC source files to Julia JPEC modules -- [Known Physics Issues in PerturbedEquilibrium](known_issues_perteq.md) — xss_mn not implemented, xi_modes placeholders +- [Known Physics Issues in PerturbedEquilibrium](known_issues_perteq.md) — xi_modes placeholders, the xi_n Jacobian shortcut - [KineticForces (NTV) Audit Checklist](kinetic_ntv_map.md) — pentrc->KineticForces map, Logan 2015 matrices, what to verify - [InnerLayer (Resistive) Audit Checklist](resistive_layer_map.md) — rmatch->InnerLayer map, GGJ Wasow basis / Δ′, what to verify - [Galerkin Δ′ Assembly Map](galerkin_assembly_map.md) — gal.f<->GalerkinAssembly.jl; resonant sign chain verified; PASS diff --git a/.claude/agent-memory/fortran-physics-reviewer/known_issues_perteq.md b/.claude/agent-memory/fortran-physics-reviewer/known_issues_perteq.md index b8d5b3b947..f403d97026 100644 --- a/.claude/agent-memory/fortran-physics-reviewer/known_issues_perteq.md +++ b/.claude/agent-memory/fortran-physics-reviewer/known_issues_perteq.md @@ -4,12 +4,6 @@ description: Documented physics issues found in perturbed equilibrium module dur type: project --- -## xss_mn (toroidal displacement) not implemented -`FieldReconstruction.jl:316` sets `xss = 0`. This makes b_theta_modes and b_zeta_modes incorrect. -b_psi_modes is unaffected. The Fortran computes xss_mn from idcon_matrix A/B/C matrices. -**Why:** Not yet ported from Fortran gpeq_sol (lines 86-96). -**How to apply:** Any downstream use of b_theta or b_zeta (e.g., tangential displacement, parallel field) requires implementing xss_mn first. - ## xi_modes.theta and xi_modes.zeta are zero placeholders `FieldReconstruction.jl:136-137`. Fortran computes these via gpeq_contra Jacobian convolution. diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index d82d998a2b..bfa26d1d7b 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -77,5 +77,5 @@ repos: - id: toml-no-deprecated-keys name: 'TOML conventions: no deprecated config keys' language: pygrep - entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|reform_eq_with_psilim|use_riccati|use_parallel|parallel_threads|populate_dense_xi|gal_flag|nstep|diagnose_ca|diagnose|ksing|power_bp|power_b|power_r|power_rc|use_galgrid|moment|eq_out|theta_out|xlmda_out|eqpsi_out|save_records|gpec_file|fnml_flag|ellip_flag|diag_flag|force_xialpha|fixed_boundary|output_eigenmodes|filter_modes|singular_point_method|compute_response|local_stability_flag|profile_file|profile_group)\s*=' + entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|reform_eq_with_psilim|use_riccati|use_parallel|parallel_threads|populate_dense_xi|gal_flag|nstep|diagnose_ca|diagnose|ksing|power_bp|power_b|power_r|power_rc|use_galgrid|moment|eq_out|theta_out|xlmda_out|eqpsi_out|save_records|gpec_file|fnml_flag|ellip_flag|diag_flag|force_xialpha|fixed_boundary|output_eigenmodes|filter_modes|singular_point_method|compute_response|local_stability_flag|gal_solver|gal_sing_order|gap_kHz_threshold|force_termination|sing_start|HDF5_filename|write_outputs_to_HDF5|output_filename|kinetic_source|nn_low|nn_high|gal_[A-Za-z0-9_]+|output_benchmark_data|ion_charge|ion_mass|chi_perp_e|chi_perp_i|chi_tor|profile_file|profile_group)\s*=|^\[(DEBUG|SOL_INPUT|LAR_INPUT|TJ_ANALYTIC_INPUT|SLAYER\.scan_grid|SLAYER\.CriticalResonantField)\]|^eq_type\s*=\s*"(sol|lar|tj_analytic|tj_analytic_direct)"' files: ^(examples/.*\.toml|test/test_data/.*\.toml)$ diff --git a/CLAUDE.md b/CLAUDE.md index beedcf59e2..f5c1ad8846 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -104,6 +104,7 @@ Roster, the recommended review pipeline (physics fidelity → readability → pe - **Regenerate the CI dependency pins when `Project.toml` changes** - CI does not resolve dependencies freshly; it copies a pinned resolve from `ci/manifests/` so runs are reproducible and the restored depot cache stays valid. After editing `[deps]` or `[compat]`, run `julia +1.11 ci/manifests/update.jl` and `julia +1.12 ci/manifests/update.jl` and commit the result. See [`ci/manifests/README.md`](ci/manifests/README.md). - **Indexing**: The codebase uses 0-based indexing in many places to match Fortran conventions, then converts to 1-based Julia indexing - **No step numbering in code comments** - Avoid annotations like "Step 1: do this" followed by "Step 2: do that". These get out of sync as code changes. Just describe the action without numbering. +- **Deck and API parity** - Every capability must be reachable from both a `gpec.toml` deck (via `main`) and the scripting API (`solve` and the problem/result types). Every new feature, deck key or API keyword adds a case to the `"deck and API parity"` testset in `test/runtests_solve_api.jl`, which runs the workflow from a deck and from the API and checks that the results agree. Mark a known drift `@test_broken`; don't drop it. See `docs/src/developer_notes.md`. - **Documentation coverage** - When adding a new module or submodule with public docstrings, add a corresponding `@autodocs` block in `docs/src/`. Documenter CI will fail with a `missing_docs` error if any exported docstring is not covered. The analysis submodule docs live in `docs/src/analysis.md`. - **Docstrings are rendered as Markdown** - Documenter parses docstrings as CommonMark, so `[text](...)` patterns become hyperlinks and will fail CI with `invalid local link/image` if the target doesn't exist. Common pitfall: unit annotations like `[degrees] (or [m] if ...)` parse as `[degrees](or [m] if ...)` — a broken link. Use plain words (`in degrees`) or backticks (`` `[m]` ``) for unit labels inside docstrings. Same rule for bracketed array-shape hints (`[ncoil]`) followed by parentheses. When in doubt, preview with `julia --project=docs docs/make.jl` before pushing. - **Keep code comments concise** - A comment should be one line where possible. Do not write multi-line block comments explaining the current session's investigation, what was tried, what was wrong before, or why a specific file/path behaves differently. State what the code does and why at a general level. Example of too much detail: a 6-line block explaining that efit_by_inversion uses psilow>0 while CHEASE starts at 0, that the old code was removed, and that spline spikes result. Preferred: `# Replicate Fortran inverse.f: overwrite deta at axis (r²=0) by extrapolating from innermost surfaces.` diff --git a/benchmarks/benchmark_against_fortran_run.jl b/benchmarks/benchmark_against_fortran_run.jl index 2b51564bb5..138fc91fb1 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, "flux_surface_rtol = 1e-7") println(io) - println(io, "[Wall]") - println(io, "shape = \"nowall\"") - println(io) println(io, "[LocalStability]") println(io) println(io, "[ForceFreeStates]") @@ -418,13 +415,14 @@ function write_gpec_toml_coil( # (= outermost rational q + dmlim). psiedge from dcon.in (=1.0 → no edge dW scan). @printf(io, "qhigh = %.4f\n", qhigh) @printf(io, "dW_edge_scan_start = %.4f\n", p.psiedge) - println(io, "nn_low = $(p.nn)") - println(io, "nn_high = $(p.nn)") + println(io, "nn = $(p.nn)") println(io, "delta_mlow = $(p.delta_mlow)") println(io, "delta_mhigh = $(p.delta_mhigh)") println(io, "mtheta_vacuum = 512") - println(io, "eulerlagrange_tolerance = 1e-7") println(io, "singfac_min = 1e-4") + println(io) + println(io, "[Riccati]") + println(io, "eulerlagrange_tolerance = 1e-7") println(io, "save_interval = $(p.euler_stride)") println(io) println(io, "[ForcingTerms]") @@ -441,8 +439,6 @@ function write_gpec_toml_coil( println(io) println(io, "[PerturbedEquilibrium]") println(io, "compute_singular_coupling = true") - println(io, "verbose = true") - println(io, "write_outputs_to_HDF5 = true") end end @@ -468,22 +464,20 @@ function write_gpec_toml_file( println(io, "mtheta = $(p.mtheta_equil)") println(io, "flux_surface_rtol = 1e-7") println(io) - println(io, "[Wall]") - println(io, "shape = \"nowall\"") - println(io) println(io, "[LocalStability]") println(io) println(io, "[ForceFreeStates]") println(io, "vac_flag = true") @printf(io, "qlow = %.4f\n", p.qlow) @printf(io, "dW_edge_scan_start = %.4f\n", p.psiedge) - println(io, "nn_low = $(p.nn)") - println(io, "nn_high = $(p.nn)") + println(io, "nn = $(p.nn)") println(io, "delta_mlow = $(p.delta_mlow)") println(io, "delta_mhigh = $(p.delta_mhigh)") println(io, "mtheta_vacuum = 512") - println(io, "eulerlagrange_tolerance = 1e-7") println(io, "singfac_min = 1e-4") + println(io) + println(io, "[Riccati]") + println(io, "eulerlagrange_tolerance = 1e-7") println(io, "save_interval = $(p.euler_stride)") println(io) println(io, "[ForcingTerms]") @@ -492,8 +486,6 @@ function write_gpec_toml_file( println(io) println(io, "[PerturbedEquilibrium]") println(io, "compute_singular_coupling = true") - println(io, "verbose = true") - println(io, "write_outputs_to_HDF5 = true") end end @@ -511,7 +503,7 @@ function setup_equil(p::FortranRunParams, fortran_dir::String) "grid_type" => p.grid_type, "flux_surface_rtol" => 1e-7 ) - eq_config = Equilibrium.EquilibriumConfig(eq_dict, fortran_dir) + eq_config = GeneralizedPerturbedEquilibrium.equilibrium_config(Dict{String,Any}("Equilibrium" => eq_dict), fortran_dir) return Equilibrium.PlasmaEquilibrium(eq_config) end diff --git a/benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl b/benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl index 8da5616830..ceea29e317 100644 --- a/benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl +++ b/benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl @@ -206,7 +206,7 @@ t_equil = @elapsed begin "psi_accuracy" => 0.001, "flux_surface_rtol" => 1e-7, ) - eq_config = Equilibrium.EquilibriumConfig(eq_dict, run_dir) + eq_config = GeneralizedPerturbedEquilibrium.equilibrium_config(Dict{String,Any}("Equilibrium" => eq_dict), run_dir) equil = Equilibrium.PlasmaEquilibrium(eq_config) end @printf " Done in %.1f s (bt_sign=%d)\n" t_equil equil.params.bt_sign diff --git a/benchmarks/benchmark_delta_prime_methods.jl b/benchmarks/benchmark_delta_prime_methods.jl deleted file mode 100644 index 94c7dc0eee..0000000000 --- a/benchmarks/benchmark_delta_prime_methods.jl +++ /dev/null @@ -1,96 +0,0 @@ -# Sanity check: compute_delta_prime_from_ca! vs inline Δ' from riccati_cross_ideal_singular_surf! -# -# riccati_cross_ideal_singular_surf! computes Δ' inline at each singular surface crossing -# using the diagonal formula (no Gaussian reduction permutation): -# Δ'[s] = (ca_r[ipert_res, ipert_res, 2, s] - ca_l[ipert_res, ipert_res, 2, s]) / (4π²·ψ₀) -# -# compute_delta_prime_from_ca! applies the identical formula post-hoc from the stored -# ca_l/ca_r arrays. Since both operate on the same data with the same formula, results -# should match to floating-point precision (not just approximately — exactly). -# -# This verifies that compute_delta_prime_from_ca! is a correct standalone implementation -# of the Δ' formula that can be used for testing or alternative integration drivers. -# -# Usage (from JPEC_main root): -# julia --project=. benchmarks/benchmark_delta_prime_methods.jl - -using LinearAlgebra, Printf, TOML -using GeneralizedPerturbedEquilibrium - -const FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates - -function setup_and_run_solovev() - ex = joinpath(@__DIR__, "..", "test", "test_data", "regression_solovev_ideal_example") - inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) - inputs["ForceFreeStates"]["verbose"] = false - inputs["ForceFreeStates"]["integrator"] = "riccati" - intr = FFS.ForceFreeStatesInternal(; dir_path=ex) - ctrl = FFS.ForceFreeStatesControl(; - (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) - eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex) - equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config) - intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(; - (Symbol(k) => v for (k, v) in inputs["Wall"])...) - FFS.sing_lim!(intr, ctrl, equil) - intr.nlow = ctrl.nn_low; - intr.nhigh = ctrl.nn_high; - intr.npert = 1 - FFS.sing_find!(intr, equil) - intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow - intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh - intr.mpert = intr.mhigh - intr.mlow + 1 - intr.numpert_total = intr.mpert * intr.npert - metric = FFS.make_metric(equil, intr.mpert) - mats = FFS.build_matrix_splines(equil, intr, metric) - odet, _, _, _ = FFS.riccati_eulerlagrange_integration(ctrl, equil, mats, intr) - return ctrl, equil, mats, intr, odet -end - -println("\n=== compute_delta_prime_from_ca! consistency check ===") -println("Verifies the standalone Δ' formula matches the inline Riccati crossing computation.") -println("Expected error: exactly zero (same formula, same data).\n") - -ctrl, equil, mats, intr, odet = setup_and_run_solovev() -msing = intr.msing - -# Capture Δ' values set inline by riccati_cross_ideal_singular_surf! during integration -delta_prime_inline = [copy(intr.sing[s].delta_prime) for s in 1:msing] - -# Now call compute_delta_prime_from_ca! — it reads the same ca_l/ca_r arrays and -# overwrites intr.sing[s].delta_prime using the identical diagonal formula -FFS.compute_delta_prime_from_ca!(odet, intr, equil) - -println(" N=$(intr.numpert_total) modes, $msing singular surfaces\n") -@printf(" %6s %4s %4s %22s %22s %12s\n", - "Surf", "m", "n", "Δ' (inline)", "Δ' (from_ca)", "abs diff") -println(" " * "-"^76) - -max_absdiff = let max_absdiff = 0.0 - for s in 1:msing - sing = intr.sing[s] - dp_from_ca = intr.sing[s].delta_prime - for i in eachindex(delta_prime_inline[s]) - dp_il = delta_prime_inline[s][i] - dp_fc = dp_from_ca[i] - absdiff = abs(dp_fc - dp_il) - max_absdiff = max(max_absdiff, absdiff) - @printf(" %6d %4d %4d %22.6f%+.6fi %22.6f%+.6fi %12.4e\n", - s, sing.m[i], sing.n[i], - real(dp_il), imag(dp_il), - real(dp_fc), imag(dp_fc), - absdiff) - end - end - max_absdiff -end - -println() -if max_absdiff == 0.0 - println("PASSED — Δ' values are bit-for-bit identical (max abs diff = 0.0)") -elseif max_absdiff < 1e-14 - @printf("PASSED — max abs diff = %.2e (floating-point rounding only)\n", max_absdiff) -else - @printf("FAILED — max abs diff = %.2e (expected exact agreement)\n", max_absdiff) - exit(1) -end -println() diff --git a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl index 8cad18ab2d..5b14c41b52 100644 --- a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl +++ b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl @@ -125,31 +125,24 @@ mtheta = 256 newq0 = 0 flux_surface_rtol = 1e-7 -[Wall] -shape = "nowall" - [LocalStability] [ForceFreeStates] vac_flag = true -force_termination = false dW_edge_scan_start = 1.0 # No edge-scan truncation (rational_offset_fraction truncation off) qlow = 1.02 qhigh = 1e3 -sing_start = 0 -nn_low = 1 -nn_high = 1 +nn = 1 delta_mlow = 8 delta_mhigh = 8 mtheta_vacuum = 512 +singfac_min = 1e-4 -kinetic_source = "fixed" -kinetic_factor = 0.0 +[Riccati] eulerlagrange_tolerance = 1e-7 save_interval = 3 -singfac_min = 1e-4 renorm_threshold = 1e4 """ ) @@ -280,10 +273,9 @@ function run_benchmark(fortran_dir::String=default_fortran_dir()) # for the final full-validation run. tgar_enabled = get(ENV, "BENCHMARK_TGAR", "0") == "1" kf_ctrl = KF.KineticForcesControl(; - fgar_flag=true, tgar_flag=tgar_enabled, nn=1, bounce_harmonic_max=4, - ion_mass=2, ion_charge=1, nutype="harmonic", f0type="maxwellian", - collisionality_factor=1, psilims=[0.0, 1.0], - verbose=false) + fgar_flag=true, tgar_flag=tgar_enabled, bounce_harmonic_max=4, + nutype="harmonic", f0type="maxwellian", + collisionality_factor=1, psilims=[0.0, 1.0]) # Initialize KF internal state kf_intr = KF.KineticForcesInternal(equil; verbose=false) @@ -309,9 +301,9 @@ function run_benchmark(fortran_dir::String=default_fortran_dir()) _p("\n--- Torque integration ---") t2 = time() - kinetic_profiles = Eq.kinetic_profile_splines( - Eq.read_kinetic_file(inputs.kin_file); - zi=kf_ctrl.ion_charge, zimp=kf_ctrl.impurity_charge, mi=kf_ctrl.ion_mass, mimp=kf_ctrl.impurity_mass) + # The default species (one deuterium main ion, carbon impurity) match pentrc.in. + Eq.attach_kinetic_profiles!(equil, inputs.kin_file) + kinetic_profiles = KF.ntv_profiles(kf_ctrl, equil) kf_state = KF.KineticForcesState() KF.compute_torque_all_methods!(kf_state, kf_intr, kf_ctrl, equil, kinetic_profiles) diff --git a/benchmarks/benchmark_diiid_kinetic_stability.jl b/benchmarks/benchmark_diiid_kinetic_stability.jl index f576994b1f..d9e77e1e55 100644 --- a/benchmarks/benchmark_diiid_kinetic_stability.jl +++ b/benchmarks/benchmark_diiid_kinetic_stability.jl @@ -5,7 +5,7 @@ Benchmark the self-consistent kinetic-DCON path (KF → FFS → PE) against Fortran GPEC's kinetic DCON reference. -Runs `GPE.main()` with `kinetic_source="calculated"` and `kinetic_factor=1.0` +Runs `GPE.main()` with `kinetic = "KineticForces"` against the EFIT g-file and `.kin` profile taken from a Fortran GPEC kinetic example directory, then compares the least-stable total-energy eigenvalue `ForceFreeStates/FreeBoundaryStability/eigenmode_energies[1]` against `W_t_eigenvalue[:, 0]` in the Fortran @@ -123,8 +123,6 @@ function build_benchmark_tomldir(eq_file::String, kin_file::String) [Equilibrium] eq_filename = "$eq_name" kinetic_file = "$kin_file" -eq_type = "efit" -jac_type = "hamada" grid_type = "pow1" # Matches Fortran DIIID_kinetic_example (equil/grid.f90 powspace) psilow = 0.01 psihigh = 0.993 @@ -133,39 +131,30 @@ mtheta = 256 newq0 = 0 flux_surface_rtol = 1e-7 -[Wall] -shape = "nowall" - [LocalStability] [ForceFreeStates] vac_flag = true -force_termination = true # Skip PE+KF post-processing — we only need FFS eigenvalues dW_edge_scan_start = 1.0 # No edge-scan truncation (rational_offset_fraction truncation off) qlow = 1.02 qhigh = 1e3 -sing_start = 0 -nn_low = 1 -nn_high = 1 +nn = 1 delta_mlow = 8 delta_mhigh = 8 mtheta_vacuum = 512 -thmax0 = 1 -kinetic_source = "calculated" -kinetic_factor = 1.0 -eulerlagrange_tolerance = 1e-7 +integrator = "Forward" +kinetic = "KineticForces" singfac_min = 1e-4 + +[Forward] +eulerlagrange_tolerance = 1e-7 renorm_threshold = 1e4 -write_outputs_to_HDF5 = true [KineticForces] -nn = 1 -nl = 4 -zi = 1 -mi = 2 +bounce_harmonic_max = 4 nutype = "harmonic" f0type = "maxwellian" """ @@ -192,7 +181,7 @@ function run_benchmark(fortran_dir::String=default_fortran_dir()) tomldir = build_benchmark_tomldir(inputs.eq_file, inputs.kin_file) - _p("\n--- Running GPE.main with kinetic_source=\"calculated\" ---") + _p("\n--- Running GPE.main with kinetic = \"KineticForces\" ---") t0 = time() GPE.main([tomldir]) wall = time() - t0 diff --git a/benchmarks/benchmark_q_vs_iota_edge.jl b/benchmarks/benchmark_q_vs_iota_edge.jl index ff1957c16f..c4a5b7d0ed 100644 --- a/benchmarks/benchmark_q_vs_iota_edge.jl +++ b/benchmarks/benchmark_q_vs_iota_edge.jl @@ -22,7 +22,7 @@ const EXAMPLE_DIR = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_examp function reference_q() inputs, _, additional_input = GPE.build_inputs_from_toml(EXAMPLE_DIR) equil_dict = merge(inputs["Equilibrium"], Dict{String,Any}("grid_type" => "core_edge_packed", "mpsi" => 1024)) - eq_config = GPE.Equilibrium.EquilibriumConfig(equil_dict, EXAMPLE_DIR) + eq_config = GPE.equilibrium_config(Dict{String,Any}("Equilibrium" => equil_dict), EXAMPLE_DIR) equil = GPE.Equilibrium.PlasmaEquilibrium(eq_config, additional_input) # The dict rides along so the mpsi scan can rebuild configs; EquilibriumConfig is immutable. return equil, equil_dict @@ -65,7 +65,7 @@ function main() for N in (16, 32, 64, 128, 256) # Same log_asymptotic knot layout the auto grid uses for fixed mpsi - cfg = GPE.Equilibrium.EquilibriumConfig(merge(equil_dict, Dict{String,Any}("mpsi" => N)), EXAMPLE_DIR) + cfg = GPE.equilibrium_config(Dict{String,Any}("Equilibrium" => merge(equil_dict, Dict{String,Any}("mpsi" => N))), EXAMPLE_DIR) knots = GPE.Equilibrium._build_psi_grid(cfg, psilow, psihigh) q_nodes = [q_ref_itp(p) for p in knots] diff --git a/benchmarks/benchmark_riccati_der.jl b/benchmarks/benchmark_riccati_der.jl index eaaa64f895..9a1d6411b5 100644 --- a/benchmarks/benchmark_riccati_der.jl +++ b/benchmarks/benchmark_riccati_der.jl @@ -13,7 +13,7 @@ # the two expressions differ — so this script compares riccati_der! against the # explicit formula rather than against sing_der!. # -# Usage (from JPEC_main root): +# Usage (from the repository root): # julia --project=. benchmarks/benchmark_riccati_der.jl using LinearAlgebra, Random, Printf, TOML @@ -22,27 +22,12 @@ using GeneralizedPerturbedEquilibrium const FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates function setup_solovev() + GPEC = GeneralizedPerturbedEquilibrium ex = joinpath(@__DIR__, "..", "test", "test_data", "regression_solovev_ideal_example") inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) - inputs["ForceFreeStates"]["verbose"] = false - intr = FFS.ForceFreeStatesInternal(; dir_path=ex) - ctrl = FFS.ForceFreeStatesControl(; - (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) - eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex) - equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config) - intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(; - (Symbol(k) => v for (k, v) in inputs["Wall"])...) - FFS.sing_lim!(intr, ctrl, equil) - intr.nlow = ctrl.nn_low; - intr.nhigh = ctrl.nn_high; - intr.npert = 1 - FFS.sing_find!(intr, equil) - intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow - intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh - intr.mpert = intr.mhigh - intr.mlow + 1 - intr.numpert_total = intr.mpert * intr.npert - metric = FFS.make_metric(equil, intr.mpert) - mats = FFS.build_matrix_splines(equil, intr, metric) + equil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(inputs, ex), GPEC.analytic_config(inputs)) + prob = GPEC.ForceBalanceProblem(equil; GPEC.force_balance_keywords(inputs, ex)...) + ctrl, intr, equil, _, mats = GPEC._prepare_solve(prob, GPEC.force_balance_algorithm(inputs)) return ctrl, equil, mats, intr end @@ -80,8 +65,8 @@ println("Test state: Hermitian S (physical constraint). Expected error: ~machine ctrl, equil, mats, intr = setup_solovev() N = intr.numpert_total -odet = FFS.OdeState(N, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing) -FFS.initialize_el_at_axis!(odet, ctrl, equil.profiles, intr) +odet = FFS.OdeState(N, ctrl.algorithm.numsteps_init, ctrl.algorithm.numunorms_init, intr.msing) +FFS.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) chunks = FFS.chunk_el_integration_bounds(odet, ctrl, intr) # 30% into each chunk: well inside the interval, away from singularities at psi_end diff --git a/benchmarks/benchmark_solovev_kinetic_stability.jl b/benchmarks/benchmark_solovev_kinetic_stability.jl index 955f5f2cb2..9d0fa56291 100644 --- a/benchmarks/benchmark_solovev_kinetic_stability.jl +++ b/benchmarks/benchmark_solovev_kinetic_stability.jl @@ -2,7 +2,7 @@ """ benchmark_solovev_kinetic_stability.jl -Benchmark the self-consistent kinetic-DCON path (`kinetic_source="calculated"`) +Benchmark the self-consistent kinetic-DCON path (`kinetic = "KineticForces"`) for the analytic **Solovev** equilibrium against Fortran GPEC's kinetic DCON. Unlike the DIIID benchmark (which reads a pre-existing `dcon_output_n1.nc`), @@ -103,7 +103,7 @@ function build_matched_fortran_deck(workdir::String) # dcon.in — kinetic flags already match; align mode band, edge truncation, n. dconin = joinpath(workdir, "dcon.in") - replace_namelist_value!(dconin, "nn", string(ffs["nn_low"])) + replace_namelist_value!(dconin, "nn", string(first(ffs["nn"]))) replace_namelist_value!(dconin, "delta_mlow", string(get(ffs, "delta_mlow", 0))) replace_namelist_value!(dconin, "delta_mhigh", string(get(ffs, "delta_mhigh", 0))) replace_namelist_value!(dconin, "qlow", string(ffs["qlow"])) @@ -171,7 +171,7 @@ function run_benchmark() fort_re, fort_im = run_fortran_reference(workdir) _p(@sprintf(" Fortran W_t[1]: Re = %.6f Im = %.6f", fort_re, fort_im)) - _p("\n--- Running GPE.main with kinetic_source=\"calculated\" ---") + _p("\n--- Running GPE.main with kinetic = \"KineticForces\" ---") jul_re, jul_im, wall = run_julia_reference() _p(@sprintf(" GPE.main completed in %.1f s", wall)) diff --git a/benchmarks/benchmark_threads.jl b/benchmarks/benchmark_threads.jl index 8a9a05dc8f..a757a85467 100644 --- a/benchmarks/benchmark_threads.jl +++ b/benchmarks/benchmark_threads.jl @@ -1,39 +1,20 @@ # Thread-scaling benchmark for the Riccati chunked propagator integration. -# Runs the Solovev (N=8) and DIIID-like (N=26) examples with integrator="riccati" +# Runs the Solovev (N=8) and DIIID-like (N=26) examples with the Riccati integrator # across 1, 2, 4, 8 threads and compares against the forward path. The chunk list is # thread-independent, so et[1] must not move with the thread count — only wall-clock. # -# Usage (from JPEC_main root): +# Usage (from the repository root): # for t in 1 2 4 8; do julia -t $t --project=. benchmarks/benchmark_threads.jl; done using GeneralizedPerturbedEquilibrium, TOML, Printf, Statistics +# One solve of the example deck with `integrator` (Forward or Riccati), sharing the deck's ODE numerics. function run_ffs(ex; integrator) - inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) - inputs["ForceFreeStates"]["verbose"] = false - inputs["ForceFreeStates"]["integrator"] = integrator - inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false - intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal(; dir_path=ex) - ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; - (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) - eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex) - equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config) - intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(; - (Symbol(k) => v for (k, v) in inputs["Wall"])...) - GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil) - intr.nlow = ctrl.nn_low; - intr.nhigh = ctrl.nn_high; - intr.npert = 1 - GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_find!(intr, equil) - intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow - intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh - intr.mpert = intr.mhigh - intr.mlow + 1 - intr.numpert_total = intr.mpert * intr.npert - metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert) - mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.build_matrix_splines(equil, intr, metric) - odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, mats, intr) - vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, mats, intr) - return real(vac.et[1]), intr.numpert_total + GPE = GeneralizedPerturbedEquilibrium + inputs, config, analytic = GPE.build_inputs_from_toml(ex) + alg = integrator(; GPE.section_keywords(inputs["Forward"], "Forward", GPE._struct_fields(Forward))...) + ffs = solve(PlasmaEquilibrium(config, analytic), alg; GPE.force_balance_keywords(inputs, ex)...) + return real(ffs.free_boundary.et[1]), ffs.numpert_total end function timed_run(ex; integrator, nwarm=1, nrep=2) @@ -60,8 +41,8 @@ diiid_ex = joinpath(root, "examples", "DIIID-like_ideal_example") println("\n=== Thread-scaling benchmark ($(nthreads) thread(s)) ===\n") for (label, ex) in [("Solovev", sol_ex), ("DIIID-like", diiid_ex)] - t_fwd, et_fwd, N = timed_run(ex; integrator="forward") - t_ric, et_ric, _ = timed_run(ex; integrator="riccati") + t_fwd, et_fwd, N = timed_run(ex; integrator=Forward) + t_ric, et_ric, _ = timed_run(ex; integrator=Riccati) err_ric = abs(et_ric - et_fwd) / abs(et_fwd) * 100 diff --git a/benchmarks/compare_gal_vs_el.jl b/benchmarks/compare_gal_vs_el.jl deleted file mode 100644 index 98c8e82efe..0000000000 --- a/benchmarks/compare_gal_vs_el.jl +++ /dev/null @@ -1,78 +0,0 @@ -# Overlay the IDEAL gal matched ξ(ψ) against the EL total-energy eigenmode ξ(ψ), for the same edge -# eigenvector w. The ideal gal solution should reproduce the EL (DCON) ideal solution. -# -# EL : ξ_EL(ψ) = U_EL(ψ) · (U_EL_edge \ w) (fundamental matrix, ForceFreeStates/Solutions/ForwardIntegration/xi_psi) -# gal: ξ_gal(ψ) = U_gal(ψ) · w (identity-at-edge ⇒ coefficient is w itself) -# w = eigenvector of the total energy operator W = W_plasma + W_vacuum (FreeBoundaryStability). -# -# Needs ONE gpec.h5 from a run with integrator="forward" (EL dense u_store), gal_match_flag=true, -# gal_ideal_flag=true (gal ideal matched), vac_flag=true (energy operator). -# Usage: julia --project=. benchmarks/compare_gal_vs_el.jl [gpec.h5] [out.png] [mode] -# mode = "highest" (default, most stable), "lowest" (most unstable), or an integer eigenmode index. - -using HDF5, LinearAlgebra, Plots, Printf - -h5path = length(ARGS) >= 1 ? ARGS[1] : "/tmp/gal_ideal_test/gpec.h5" -outpng = length(ARGS) >= 2 ? ARGS[2] : joinpath(@__DIR__, "compare_gal_vs_el.png") -ksel = length(ARGS) >= 3 ? ARGS[3] : "highest" - -to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a) - -et, wt, u1, psiE, gxi, psiG, mlow, sing_psi = h5open(h5path) do f - (to_c(read(f["ForceFreeStates/FreeBoundaryStability/eigenmode_energies"])), - to_c(read(f["ForceFreeStates/FreeBoundaryStability/W_freeboundary_eigenmodes"])), - to_c(read(f["ForceFreeStates/Solutions/ForwardIntegration/xi_psi"])), read(f["ForceFreeStates/Solutions/ForwardIntegration/psi"]), - to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/xi_psi"])), read(f["ForceFreeStates/Solutions/GalerkinIntegration/psi"]), - read(f["Info/mlow"]), - read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"])) -end - -mpert = size(u1, 1) -k = ksel == "highest" ? argmax(real.(et)) : - ksel == "lowest" ? argmin(real.(et)) : - parse(Int, ksel) -w = wt[:, k] -@printf("comparing eigenmode k=%d δW = %.4e %+.4ei\n", k, real(et[k]), imag(et[k])) - -# EL profile -cEL = u1[:, :, end] \ w -xiE = reduce(hcat, (u1[:, :, ip] * cEL for ip in 1:size(u1, 3))) # (mpert, nE) - -# gal ideal profile (identity-at-edge ⇒ coefficient = w); the closed-profile grid -# already excludes on-surface points -psiGk = psiG -xiG = reduce(hcat, (gxi[:, :, ip] * w for ip in eachindex(psiG))) # (mpert, nGk) - -ms = mlow .+ (0:mpert-1) -peak = [maximum(abs, @view xiE[i, :]) for i in 1:mpert] -order = sortperm(peak; rev=true) -ndom = min(5, mpert) - -# quantitative bulk agreement on the dominant harmonic (linear-interp EL onto gal grid, off rationals/edge) -lininterp(xq, x, y) = (j = clamp(searchsortedlast(x, xq), 1, length(x) - 1); -t = (xq - x[j]) / (x[j+1] - x[j]); y[j] * (1 - t) + y[j+1] * t) -idom = order[1] -inbulk(p) = 0.1 <= p <= 0.92 && all(abs(p - ps) > 8e-3 for ps in sing_psi) -sel = [ip for ip in eachindex(psiGk) if inbulk(psiGk[ip])] -num = sum(abs2(xiG[idom, ip] - lininterp(psiGk[ip], psiE, xiE[idom, :])) for ip in sel) -den = sum(abs2(xiG[idom, ip]) for ip in sel) -@printf("dominant m=%d: bulk relative L2 ‖gal−EL‖/‖gal‖ = %.3e (n=%d, off-layer/off-edge)\n", - ms[idom], sqrt(num / den), length(sel)) - -plt = plot(; size=(1000, 650), xlabel="ψ_N", ylabel="|ξ^ψ_m(ψ)|", - title=@sprintf("Ideal gal vs EL total-energy eigenmode (δW=%.2e) — solid: EL, dashed: gal", real(et[k])), - legend=:topleft, left_margin=12Plots.mm, bottom_margin=5Plots.mm, right_margin=4Plots.mm) -palette = theme_palette(:auto) -for (c, i) in enumerate(order[1:ndom]) - col = palette[mod1(c, length(palette))] - plot!(plt, psiE, abs.(@view xiE[i, :]); color=col, lw=2.2, label="EL m=$(ms[i])") - plot!(plt, psiGk, abs.(@view xiG[i, :]); color=col, lw=1.6, ls=:dash, label="gal m=$(ms[i])") -end -for (j, ps) in enumerate(sing_psi) - vline!(plt, [ps]; color=:black, ls=:dot, lw=1, label=(j == 1 ? "rational q" : "")) -end - -display(plt) -savefig(plt, outpng) -println("dominant harmonics (m): ", [ms[i] for i in order[1:ndom]]) -println("saved: ", abspath(outpng)) diff --git a/benchmarks/compare_jbgradpsi_m2.jl b/benchmarks/compare_jbgradpsi_m2.jl index 35816d6667..e840afe3ed 100644 --- a/benchmarks/compare_jbgradpsi_m2.jl +++ b/benchmarks/compare_jbgradpsi_m2.jl @@ -1,7 +1,7 @@ # Compare the area-normalized b^ψ (PerturbedEquilibrium/Response/psi_area = b^ψ/⟨J·|∇ψ|⟩_θ) for one # poloidal harmonic between two GPEC runs that are identical except for which ξ feeds PerturbedEquilibrium: -# (1) IDEAL galerkin matched ξ (gal_match_flag=true, gal_ideal_flag=true) -# (2) FORWARD ξ (gal_match_flag=false) +# (1) IDEAL galerkin matched ξ (integrator = "Galerkin" with [Match] ideal = true) +# (2) FORWARD ξ (integrator = "Forward", no [Match]) # # PE writes no ψ grid, so it's reconstructed: gal-ideal → ForceFreeStates/Solutions/GalerkinIntegration/psi (already excludes on-surface points); # forward → ForceFreeStates/Solutions/ForwardIntegration/psi. diff --git a/benchmarks/equil_method_comparison.jl b/benchmarks/equil_method_comparison.jl index 27fe93f95e..5f83146009 100644 --- a/benchmarks/equil_method_comparison.jl +++ b/benchmarks/equil_method_comparison.jl @@ -38,7 +38,7 @@ function make_config(path::String, eq_type::String; psihigh=nothing) raw["Equilibrium"]["eq_type"] = eq_type psihigh !== nothing && (raw["Equilibrium"]["psihigh"] = psihigh) base = dirname(path) - return Equilibrium.EquilibriumConfig(raw["Equilibrium"], base) + return GeneralizedPerturbedEquilibrium.equilibrium_config(raw, base) end # ─── Run each method (3 warm runs, report average of last 2) ───────────────── diff --git a/benchmarks/equil_numerical_params.jl b/benchmarks/equil_numerical_params.jl index e6d8e86168..42c2e45fd3 100644 --- a/benchmarks/equil_numerical_params.jl +++ b/benchmarks/equil_numerical_params.jl @@ -7,7 +7,7 @@ Sweeps: - mpsi ∈ {64, 128, 256} (all methods) - mtheta ∈ {128, 256, 512} (all methods) - ODE reltol ∈ {1e-5, 1e-7, 1e-9} (efit, efit_arclength) - - Contour grid refinement ∈ {2, 4, 8} (efit_by_inversion, via REFINE env var override) + - Contour grid resolution_factor ∈ {2, 4, 8} (efit_by_inversion) Accuracy metric: max roundtrip error over a sparse (ψ, θ) sample grid. Performance metric: wall-clock runtime. @@ -28,7 +28,7 @@ function make_config(path, eq_type, mpsi, mtheta) raw["Equilibrium"]["eq_type"] = eq_type raw["Equilibrium"]["mpsi"] = mpsi raw["Equilibrium"]["mtheta"] = mtheta - return Equilibrium.EquilibriumConfig(raw["Equilibrium"], dirname(path)) + return GeneralizedPerturbedEquilibrium.equilibrium_config(raw, dirname(path)) end function roundtrip_error(pe, raw_profile) @@ -121,13 +121,13 @@ println("ODE reltol scan — NOTE: not yet exposed via config; using default rel println("(Add reltol to EquilibriumConfig and threading to enable this scan)") println("=" ^ 65) -# ── Contour grid refinement scan (efit_by_inversion only) ───────────────────── -# The refine factor in equilibrium_solver_by_inversion is a keyword arg. +# ── Contour grid resolution scan (efit_by_inversion only) ───────────────────── +# resolution_factor scales the adaptive contour grid of equilibrium_solver_by_inversion. # We call it directly here, bypassing the PlasmaEquilibrium constructor. println("\n" * "=" ^ 65) -println("Contour grid refinement scan (efit_by_inversion only, mpsi=128, mtheta=256)") +println("Contour grid resolution scan (efit_by_inversion only, mpsi=128, mtheta=256)") println("=" ^ 65) -for refine in [2, 4, 8] +for resolution_factor in [2.0, 4.0, 8.0] cfg = make_config(config_path, "efit_by_inversion", 128, 256) raw_input = Equilibrium.read_efit(cfg) local pe = nothing @@ -136,7 +136,7 @@ for refine in [2, 4, 8] local rt_err = NaN try t = @elapsed begin - pe_raw = Equilibrium.equilibrium_solver_by_inversion(raw_input; refine=refine) + pe_raw = Equilibrium.equilibrium_solver_by_inversion(raw_input; resolution_factor) pe = pe_raw Equilibrium.equilibrium_global_parameters!(pe) Equilibrium.equilibrium_qfind!(pe) @@ -146,10 +146,10 @@ for refine in [2, 4, 8] rt_err = roundtrip_error(pe, raw_for_check) success = true catch e - @warn " refine=$refine FAILED: $e" + @warn " resolution_factor=$resolution_factor FAILED: $e" end - @printf(" efit_by_inversion refine=%2d t=%.2fs rt=%.2e\n", refine, runtime, rt_err) - push!(rows, (sweep="refine", method="efit_by_inversion", param=Float64(refine), param2=128.0, + @printf(" efit_by_inversion resolution_factor=%3.1f t=%.2fs rt=%.2e\n", resolution_factor, runtime, rt_err) + push!(rows, (sweep="resolution_factor", method="efit_by_inversion", param=resolution_factor, param2=128.0, success=success, runtime=runtime, roundtrip_max=rt_err)) end diff --git a/benchmarks/equil_psihigh_scan.jl b/benchmarks/equil_psihigh_scan.jl index 9bb89f6f24..9a8c9dc78a 100644 --- a/benchmarks/equil_psihigh_scan.jl +++ b/benchmarks/equil_psihigh_scan.jl @@ -39,7 +39,7 @@ function make_config(path::String, eq_type::String, psihigh::Float64) raw["Equilibrium"]["eq_type"] = eq_type raw["Equilibrium"]["psihigh"] = psihigh base = dirname(path) - return Equilibrium.EquilibriumConfig(raw["Equilibrium"], base) + return GeneralizedPerturbedEquilibrium.equilibrium_config(raw, base) end """ @@ -52,9 +52,8 @@ function run_ffs_et1(config_path::String, eq_type::String, psihigh::Float64)::Fl raw = TOML.parsefile(config_path) raw["Equilibrium"]["eq_type"] = eq_type raw["Equilibrium"]["psihigh"] = psihigh - raw["ForceFreeStates"]["force_termination"] = true # write gpec.h5 after FFS, skip PE raw["ForceFreeStates"]["vac_flag"] = true - raw["ForceFreeStates"]["verbose"] = false + foreach(section -> delete!(raw, section), ("PerturbedEquilibrium", "ErrorFields")) # stop after FFS mktempdir() do tmpdir # Write modified config diff --git a/benchmarks/equil_psilow_scan.jl b/benchmarks/equil_psilow_scan.jl index 882f87ed32..c7288949ce 100644 --- a/benchmarks/equil_psilow_scan.jl +++ b/benchmarks/equil_psilow_scan.jl @@ -29,7 +29,7 @@ function make_config(path::String, eq_type::String, psilow::Float64) raw["Equilibrium"]["eq_type"] = eq_type raw["Equilibrium"]["psilow"] = psilow base = dirname(path) - return Equilibrium.EquilibriumConfig(raw["Equilibrium"], base) + return GeneralizedPerturbedEquilibrium.equilibrium_config(raw, base) end # Read EFIT q-axis value from header (simag line) diff --git a/benchmarks/equil_spline_comparison.jl b/benchmarks/equil_spline_comparison.jl index 22ffacbc0b..8dbbb521b5 100644 --- a/benchmarks/equil_spline_comparison.jl +++ b/benchmarks/equil_spline_comparison.jl @@ -48,7 +48,7 @@ function make_config(path::String, eq_type::String, psihigh::Float64) raw = TOML.parsefile(path) raw["Equilibrium"]["eq_type"] = eq_type raw["Equilibrium"]["psihigh"] = psihigh - return Equilibrium.EquilibriumConfig(raw["Equilibrium"], dirname(path)) + return GeneralizedPerturbedEquilibrium.equilibrium_config(raw, dirname(path)) end # ─── Load all three equilibria ───────────────────────────────────────────────── diff --git a/benchmarks/plot_xi_eigenmode.jl b/benchmarks/plot_xi_eigenmode.jl deleted file mode 100644 index 78e26ee823..0000000000 --- a/benchmarks/plot_xi_eigenmode.jl +++ /dev/null @@ -1,61 +0,0 @@ -# Plot the ξ(ψ) profile of the eigenmode with the highest total-energy eigenvalue. -# -# The total energy operator W = W_plasma + W_vacuum (free_run, Free.jl); its eigenvectors are the -# free-boundary edge displacement patterns and the eigenvalues are δW. This picks the eigenmode with the -# largest Re(eigenvalue) and reconstructs its radial profile by projecting the EL fundamental matrix -# (ForceFreeStates/Solutions/ForwardIntegration/xi_psi) onto that edge eigenvector: c = U_edge \ w, ξ(ψ) = U(ψ)·c. -# -# Requires a run with integrator="forward" so ForceFreeStates/Solutions/ForwardIntegration/xi_psi is the dense axis-basis fundamental -# matrix (not the Riccati S-matrices). -# Usage: julia --project=. benchmarks/plot_xi_eigenmode.jl [path/to/gpec.h5] [out.png] - -using HDF5, LinearAlgebra, Plots, Printf - -h5path = length(ARGS) >= 1 ? ARGS[1] : "examples/DIIID-like_ideal_example/gpec.h5" -outpng = length(ARGS) >= 2 ? ARGS[2] : joinpath(@__DIR__, "xi_eigenmode.png") - -to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a) - -et, wt, u1, psi, mlow, sing_psi = h5open(h5path) do f - (to_c(read(f["ForceFreeStates/FreeBoundaryStability/eigenmode_energies"])), - to_c(read(f["ForceFreeStates/FreeBoundaryStability/W_freeboundary_eigenmodes"])), - to_c(read(f["ForceFreeStates/Solutions/ForwardIntegration/xi_psi"])), - read(f["ForceFreeStates/Solutions/ForwardIntegration/psi"]), - read(f["Info/mlow"]), - haskey(f, "ForceFreeStates/Solutions/GalerkinIntegration/rational_psi") ? read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"]) : Float64[]) -end - -mpert, _, nstep = size(u1) -k = argmax(real.(et)) # eigenmode with the highest total-energy eigenvalue -@printf("highest-eigenvalue mode: k=%d δW = %.4e %+.4ei (of %d modes)\n", - k, real(et[k]), imag(et[k]), length(et)) - -# Reconstruct ξ(ψ) for the eigenmode whose edge displacement is the eigenvector w = wt[:,k]. -w = wt[:, k] -U_edge = u1[:, :, nstep] # edge fundamental matrix (psi inner→edge, last index = edge) -c = U_edge \ w # U_edge·c = w -xi = Array{ComplexF64}(undef, mpert, nstep) -for ip in 1:nstep - @views mul!(xi[:, ip], u1[:, :, ip], c) -end - -ms = mlow .+ (0:mpert-1) -peak = [maximum(abs, @view xi[i, :]) for i in 1:mpert] -order = sortperm(peak; rev=true) # dominant harmonics first - -plt = plot(; size=(950, 620), xlabel="ψ_N", ylabel="|ξ^ψ_m(ψ)|", - title=@sprintf("Highest total-energy eigenmode δW = %.3e", real(et[k])), - legend=:topleft, left_margin=12Plots.mm, bottom_margin=5Plots.mm, right_margin=4Plots.mm) -for i in 1:mpert # all harmonics, faint - plot!(plt, psi, abs.(@view xi[i, :]); color=:gray85, lw=0.6, label="") -end -for i in order[1:min(8, mpert)] # dominant harmonics, bold + labeled - plot!(plt, psi, abs.(@view xi[i, :]); lw=2, label="m=$(ms[i])") -end -for (j, ps) in enumerate(sing_psi) # rational surfaces - vline!(plt, [ps]; color=:black, ls=:dash, lw=1, label=(j == 1 ? "rational q" : "")) -end - -savefig(plt, outpng) -println("dominant harmonics (m): ", [ms[i] for i in order[1:min(8, mpert)]]) -println("saved: ", abspath(outpng)) diff --git a/benchmarks/verify_gal_ideal.jl b/benchmarks/verify_gal_ideal.jl deleted file mode 100644 index 4eb981a07b..0000000000 --- a/benchmarks/verify_gal_ideal.jl +++ /dev/null @@ -1,33 +0,0 @@ -# Verify the gal_ideal_flag path: the matched solution is the bare ideal coil column (cout=0). -# 1. cout / cin / deltar all zero (no resistive plasma combination, no inner layer) -# 2. matched ξ_j == the gal coil column sols(:,:,2·msing+j) (and ξ′ likewise) -# Usage: julia --project=. benchmarks/verify_gal_ideal.jl [path/to/gpec.h5] -# Requires [DEBUG] gal_basis_output = true in the deck (compares against the raw-basis dump). -using HDF5, Printf, LinearAlgebra - -h5path = length(ARGS) >= 1 ? ARGS[1] : "/tmp/gal_ideal_test/gpec.h5" -to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a) - -cout, deltar, mxi, mdxi, sols, sols_d, iss, sing_psi = h5open(h5path) do f - (to_c(read(f["SingularSurfaces/Match/cout"])), to_c(read(f["SingularSurfaces/Match/Delta_r"])), - to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/xi_psi"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/dxi_psidpsi"])), - to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/xi_psi"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/dxi_psidpsi"])), - Bool.(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/is_rational"])), - read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"])) -end -msing = length(sing_psi) -mpert, mcoil, ngrid = size(mxi) -keep = .!iss - -@printf("[1] ‖cout‖ = %.2e, ‖deltar‖ = %.2e %s\n", norm(cout), norm(deltar), - (norm(cout) == 0 && norm(deltar) == 0) ? "✓ no resistive combination (ideal)" : "✗") - -# matched ξ_j should equal the coil column sols(:,:,2msing+j) -err = maximum(1:mcoil) do j - csol = 2msing + j - max(norm(mxi[:, j, :] - sols[:, csol, keep]), norm(mdxi[:, j, :] - sols_d[:, csol, keep])) -end -@printf("[2] max‖ξ_matched − coil column‖ = %.2e %s\n", err, - err < 1e-12 ? "✓ matched ξ == bare ideal coil column" : "✗") - -println((norm(cout) == 0 && err < 1e-12) ? "\n✓ gal_ideal_flag path verified" : "\n✗ investigate") diff --git a/benchmarks/verify_gal_match.jl b/benchmarks/verify_gal_match.jl deleted file mode 100644 index 507f5bfaaa..0000000000 --- a/benchmarks/verify_gal_match.jl +++ /dev/null @@ -1,50 +0,0 @@ -# Piece 2 verification: RPEC outer↔inner matched solution (closed profiles under -# GalerkinIntegration/xi_psi + matching diagnostics under SingularSurfaces/Match/*). -# 1. linear-solve residual ‖mat·cof − rmat‖/‖rmat‖ -# 2. matched ξ / ξ′ finiteness -# 3. edge column == identity basis: each coil drive j must give ξ_edge = e_j (the j-th harmonic), -# since plasma columns vanish at the edge and coil column j carries the unit edge source -# 4. inner-layer Δ(Q) per surface (finite, nonzero); forced eigenvalues γ_s = 2πi·n·f -# Usage: julia --project=. benchmarks/verify_gal_match.jl [path/to/gpec.h5] -using HDF5, Printf, LinearAlgebra - -h5path = length(ARGS) >= 1 ? ARGS[1] : "examples/DIIID-like_gal_resistive_example/gpec.h5" -@info "Reading $h5path" - -xi, dxi, cout, cin, deltar, eig, resid, sing_psi = h5open(h5path) do f - (read(f["ForceFreeStates/Solutions/GalerkinIntegration/xi_psi"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/dxi_psidpsi"]), - read(f["SingularSurfaces/Match/cout"]), read(f["SingularSurfaces/Match/cin"]), - read(f["SingularSurfaces/Match/Delta_r"]), read(f["SingularSurfaces/Match/rpec_eig"]), - read(f["SingularSurfaces/Match/residual"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"])) -end -# HDF5 stores ComplexF64 as a compound (re,im); convert if needed -to_c(a) = eltype(a) <: Complex ? a : map(x -> ComplexF64(x.re, x.im), a) -xi = to_c(xi); dxi = to_c(dxi); cout = to_c(cout); cin = to_c(cin); deltar = to_c(deltar); eig = to_c(eig) - -mpert, mcoil, ngrid = size(xi) -msing = size(deltar, 1) -@printf("matched solution: mpert=%d ngrid=%d mcoil=%d msing=%d\n", mpert, ngrid, mcoil, msing) - -@printf("\n[1] linear-solve residual ‖mat·cof − rmat‖/‖rmat‖ = %.3e %s\n", - resid, resid < 1e-8 ? "✓" : "✗ (large)") - -nbad = count(!isfinite, xi) + count(!isfinite, dxi) -@printf("[2] finiteness: non-finite entries = %d |ξ|max=%.3e |ξ′|max=%.3e %s\n", - nbad, maximum(abs, xi), maximum(abs, dxi), nbad == 0 ? "✓" : "✗") - -# [3] edge column == identity (last grid point = psihigh edge) -edge = xi[:, :, ngrid] # (mpert mode, mcoil drive) -id_err = norm(edge - Matrix{ComplexF64}(I, mpert, mcoil)) / sqrt(mpert) -@printf("[3] edge basis: ‖ξ(edge) − I‖/√mpert = %.3e %s\n", - id_err, id_err < 1e-8 ? "✓ identity-at-edge" : "✗") - -# [4] inner-layer Δ and forced eigenvalues -@printf("\n[4] per-surface inner-layer Δ(Q) and forced eigenvalue γ = 2πi·n·f:\n") -@printf(" %-4s %-10s %-26s %-26s\n", "s", "ψ_s", "Δ₁", "Δ₂") -for s in 1:msing - @printf(" %-4d %-10.5f %+.4e%+.4ei %+.4e%+.4ei γ=%+.3f%+.3fi\n", - s, sing_psi[s], real(deltar[s, 1]), imag(deltar[s, 1]), - real(deltar[s, 2]), imag(deltar[s, 2]), real(eig[s]), imag(eig[s])) -end -allgood = resid < 1e-8 && nbad == 0 && id_err < 1e-8 && all(isfinite, deltar) && all(!iszero, deltar) -println(allgood ? "\n✓ Piece 2 matching checks PASS" : "\n✗ some checks failed — investigate") diff --git a/benchmarks/verify_gal_solution.jl b/benchmarks/verify_gal_solution.jl deleted file mode 100644 index 6404b65959..0000000000 --- a/benchmarks/verify_gal_solution.jl +++ /dev/null @@ -1,78 +0,0 @@ -# Piece 1 verification: reconstructed gal ξ(ψ) and analytic ξ′(ψ). -# 1. shapes / finiteness sanity of the raw-basis dump (GalerkinIntegration/Basis; requires -# [DEBUG] gal_basis_output = true in the deck) -# 2. analytic ξ′ vs centered finite-difference of ξ — agree in the smooth interior, diverge at the -# packed edge (the spline-endpoint-derivative artifact we deliberately avoid) -# Usage: julia --project=. verify_gal_solution.jl [path/to/gpec.h5] -using HDF5, Printf, Statistics - -h5path = length(ARGS) >= 1 ? ARGS[1] : "examples/DIIID-like_gal_resistive_example/gpec.h5" -@info "Reading $h5path" - -psi, issing, xi, dxi, sing_psi = h5open(h5path) do f - (read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/psi"]), - read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/is_rational"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/xi_psi"]), - read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/dxi_psidpsi"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"])) -end -issing = Bool.(issing) -mpert, nsol, ngrid = size(xi) -@printf("grid: mpert=%d ngrid=%d nsol=%d psi∈[%.4f, %.4f]\n", mpert, ngrid, nsol, psi[1], psi[end]) -@printf("singular surfaces at psi = %s\n", join((@sprintf("%.4f", p) for p in sing_psi), ", ")) - -# --- finiteness (skip on-surface points, which are intentionally left zero) --- -good = .!issing -nbad = count(!isfinite, xi[:, :, good]) + count(!isfinite, dxi[:, :, good]) -@printf("non-finite entries (off-surface): %d |xi|max=%.3e |dxi|max=%.3e\n", - nbad, maximum(abs, xi[:, :, good]), maximum(abs, dxi[:, :, good])) - -# --- analytic ξ′ vs centered finite difference of ξ --- -# For each column, centered diff at interior grid points (using off-surface neighbours), compared to the -# analytic xi_deriv. Report median relative error binned by distance to the nearest singular surface and -# to the packed edge. -dist_sing(p) = isempty(sing_psi) ? Inf : minimum(abs.(p .- sing_psi)) -relerr = Float64[] -dsurf = Float64[] -dedge = Float64[] -for isol in 1:nsol - for ip in 2:(ngrid-1) - (issing[ip] || issing[ip-1] || issing[ip+1]) && continue - h1 = psi[ip] - psi[ip-1] - h2 = psi[ip+1] - psi[ip] - (h1 <= 0 || h2 <= 0) && continue - for m in 1:mpert - # nonuniform centered difference - fd = (xi[m, isol, ip+1] - xi[m, isol, ip-1]) / (h2 + h1) - an = dxi[m, isol, ip] - scale = max(abs(an), abs(fd), 1e-30) - push!(relerr, abs(fd - an) / scale) - push!(dsurf, dist_sing(psi[ip])) - push!(dedge, psi[end] - psi[ip]) - end - end -end - -@printf("\nFinite-diff(ξ) vs analytic ξ′ (%d samples)\n", length(relerr)) -println(" binned by distance to nearest singular surface:") -for (lo, hi) in [(0.0, 1e-3), (1e-3, 1e-2), (1e-2, 1e-1), (1e-1, Inf)] - sel = (dsurf .>= lo) .& (dsurf .< hi) - n = count(sel) - n == 0 && continue - @printf(" |ψ-ψ_s|∈[%-7.0e,%-7.0e) n=%-7d median relerr=%.2e p90=%.2e\n", - lo, hi, n, median(relerr[sel]), quantile(relerr[sel], 0.9)) -end -println(" binned by distance to packed edge ψ=psihigh:") -for (lo, hi) in [(0.0, 1e-4), (1e-4, 1e-3), (1e-3, 1e-2), (1e-2, Inf)] - sel = (dedge .>= lo) .& (dedge .< hi) - n = count(sel) - n == 0 && continue - @printf(" edge dist∈[%-7.0e,%-7.0e) n=%-7d median relerr=%.2e p90=%.2e\n", - lo, hi, n, median(relerr[sel]), quantile(relerr[sel], 0.9)) -end - -# Interior (away from both layer and edge) should agree well; near edge/layer the finite diff degrades. -interior = (dsurf .> 1e-2) .& (dedge .> 1e-2) -@printf("\nINTERIOR (>1e-2 from layer AND edge): median relerr = %.2e (n=%d)\n", - median(relerr[interior]), count(interior)) -println(median(relerr[interior]) < 1e-3 ? - " ✓ analytic ξ′ matches finite-diff in the smooth interior" : - " ✗ interior mismatch — investigate") diff --git a/ci/conventions/__pycache__/check_subject.cpython-310.pyc b/ci/conventions/__pycache__/check_subject.cpython-310.pyc new file mode 100644 index 0000000000..6743b0f19c Binary files /dev/null and b/ci/conventions/__pycache__/check_subject.cpython-310.pyc differ diff --git a/ci/conventions/check_subject.py b/ci/conventions/check_subject.py index 543844ccc9..ecb8b4d8b4 100644 --- a/ci/conventions/check_subject.py +++ b/ci/conventions/check_subject.py @@ -48,6 +48,7 @@ "CI", "Docs", "Examples", + "HDF5Schema", "Regression", "Repo", "Test", @@ -137,7 +138,7 @@ def module_areas(root: Path) -> list[str]: if any(sub.glob("*.jl")): areas.append(f"{module.name}.{sub.name}") - # Top-level sources such as Rerun.jl and HDF5Schema.jl are components in + # Top-level sources such as Rerun.jl are components in # their own right and are committed against directly. for source in sorted(src.glob("*.jl")): if source.name != PACKAGE_ENTRY: diff --git a/docs/development/architecture.md b/docs/development/architecture.md index 2f83974d3f..781f0b1559 100644 --- a/docs/development/architecture.md +++ b/docs/development/architecture.md @@ -30,8 +30,8 @@ Splines are provided by the external `FastInterpolations` package rather than by - Supports multiple equilibrium types: - `efit` - EFIT g-file format - `chease`, `chease2` - CHEASE equilibrium code formats - - `lar` - Large Aspect Ratio analytical model - - `sol` - Solovev analytical equilibrium + - analytic models, selected by their deck section rather than `eq_type`: `[Solovev]`, + `[LargeAspectRatio]`, `[TJAnalytic]` - Key files, by pipeline stage: - `Config.jl` - `EquilibriumConfig`, `EquilibriumParameters` - `PlasmaEquilibrium.jl` - The result type and its profile and geometry splines @@ -116,7 +116,8 @@ Splines are provided by the external `FastInterpolations` package rather than by 9. **Tearing** (`src/Tearing/`) - Tearing mode dispersion and drivers - `Dispersion/` - Dispersion relation solvers - - `Runner/` - Orchestration across surfaces and toroidal mode numbers + - `TearingProblem.jl`, `GGJTearing.jl`, `SLAYERTearing.jl` - The tearing problem and its solve with each inner-layer model + - `CriticalResonantField/` - Critical resonant field (field penetration threshold) - Re-binds `InnerLayer` and exposes it alongside its own submodules 10. **KineticForces** (`src/KineticForces/`) - Kinetic contributions to the force balance @@ -146,9 +147,9 @@ All GPEC modules are configured via a single TOML file with the following sectio - `[PerturbedEquilibrium]` - Perturbed equilibrium settings - `[ForcingTerms]` - External field specification -Key parameters: -- `force_termination` - Set to `true` to exit after equilibrium/stability (skip perturbed equilibrium) -- `output_file` - Output filename (default: `gpec.h5`) +The `[Run]` section holds the run-wide settings: `output_file` (default `gpec.h5`), the one file every +stage saves into, and `verbose`. A stage runs when its section is present, so a deck stops early by +omitting the later sections. Example configuration files are provided in: - `examples/Solovev_ideal_example/gpec.toml` @@ -204,7 +205,7 @@ The complete GPEC analysis pipeline: ### Vacuum - `VacuumInput` - Input parameters for vacuum calculations -- `WallShapeSettings` - Wall geometry configuration +- `WallSettings` - Wall geometry configuration ### Stability - `SingType` - Singular surface data including: diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md index 857172b5e5..99d1887460 100644 --- a/docs/development/hdf5-conventions.md +++ b/docs/development/hdf5-conventions.md @@ -20,7 +20,7 @@ These rules govern `gpec.h5` (and any future GPEC-produced HDF5 output); harness - **Data-driven tokens are stored verbatim**: coil-set names under `Input/RawInputs/Coils/`, KineticForces method tokens (`fgar`, …), NTV species labels under `KineticForces/PerSpecies/` (`ion_z1_m2`, `impurity_z6_m12`, `electron` — charge and mass identify the species, with a numeric discriminator appended only if a run repeats a `(z, m)` pair), scan indices (`Surface_`, `psi_`). - **Word-valued names and boolean flags**: multi-word dataset names are snake_case English (`resonance_psi`, `trajectory_offsets`, `layer_widths`), never CamelCase — CamelCase is reserved for groups. A boolean flag is named for the state it asserts when true, with an `is_` prefix only where the bare word would read as a noun or collide with a data family: `is_rational` (bare `rational` would clash with the `rational_*` coordinate family) versus `enabled`, `truncated`, `no_root`, which already read as predicates. - **Literature capitalization for physics symbols**: names match the standard literature — `D_I`, `D_R`, `Delta_prime`, `tau_R`, `tau_A` (not `di`, `dr`, `delta_prime`, `taur`); lowercase stays where the literature is lowercase (`alpha`, `q`, `beta*`, `delta_s`). -- **Scalar equilibrium parameters spell out the physics**: `R_axis`, `Z_axis`, `B_T_axis`, `a_mean`, `aspect_ratio`, `I_p`, `q_edge`, `beta_N`, `delta_upper`/`delta_lower`. The qualifier is a trailing subscript (`_axis`, `_edge`, `_wall`, `_min`, `_max`, `_upper`, `_lower`, `_extremum`), and a numbered literature definition keeps its number as the last subscript (`beta_p_1`, `l_i_2`). Fortran-era contractions (`bt0`, `amean`, `crnt`, `qa`, `betan`, `li1`) survive only as `Equilibrium.EquilibriumParameters` struct fields: `EQUIL_H5_NAMES` in `src/HDF5Schema.jl` maps each field to its dataset name, and `EQUIL_H5_SKIP` drops fields that duplicate another dataset or echo a control flag. +- **Scalar equilibrium parameters spell out the physics**: `R_axis`, `Z_axis`, `B_T_axis`, `a_mean`, `aspect_ratio`, `I_p`, `q_edge`, `beta_N`, `delta_upper`/`delta_lower`. The qualifier is a trailing subscript (`_axis`, `_edge`, `_wall`, `_min`, `_max`, `_upper`, `_lower`, `_extremum`), and a numbered literature definition keeps its number as the last subscript (`beta_p_1`, `l_i_2`). Fortran-era contractions (`bt0`, `amean`, `crnt`, `qa`, `betan`, `li1`) survive only as `Equilibrium.EquilibriumParameters` struct fields: `EQUIL_H5_NAMES` in `src/Utilities/HDF5Schema.jl` maps each field to its dataset name, and `EQUIL_H5_SKIP` drops fields that duplicate another dataset or echo a control flag. - **Derivatives are `ddpsi`** (`dTdpsi`, `dVdpsi`, `dqdpsi`, `dxidpsi`) — never Fortran `1` suffixes or `_deriv`. - **"rational" over "singular"** in dataset names (`rational_psi`, `rational_q`, `rational_m`, `rational_n`, `rational_index`, `rational_count`) — kinetic/resistive runs are not singular at the rationals. Specifier order is standardized specifier-first (`rational_psi`, never `psi_rational`). - **Vector components** follow `[d][_]_[dpsi]` — the variable always comes first and the coordinate is always the trailing subscript. A bare coordinate suffix is the **contravariant** component (`xi_psi` = ξ^ψ), `_cov_` marks the **covariant** one (`b_cov_theta` = b_θ), a leading `J` marks a **Jacobian-weighted** component (`Jxi_theta` = J·ξ^θ), and other representations sit in the same slot (`xi_clebsch_psi`, `dxi_clebsch_psidpsi`). Never drop the variable: `clebsch_psi` is wrong because it does not say *what* is being represented. There is no HDF5/netCDF standard for super- vs subscripts — flat `_` names are universal — so the typeset form always appears in the dataset's `long_name`. @@ -29,7 +29,7 @@ These rules govern `gpec.h5` (and any future GPEC-produced HDF5 output); harness ## Inputs live only under `Input/` -`Input/gpec_toml_raw` stores the full merged TOML, and `Input/RawInputs/` stores the raw equilibrium/forcing/coil data — together they make `gpec.h5` a self-contained rerun snapshot (`Rerun.jl` reconstructs every control struct from them; the writer and rerun reader carry cross-reference comments marking the mirrored path pair). **Never echo TOML flags or control-struct values into any other group** — every group outside `Input/` is derived output. (The former `kinetic/` and `slayer/settings/` echoes were removed under this rule.) +`Input/gpec_toml_raw` stores the deck that reproduces the run: each result's save adds the sections stating the problem, algorithm and controls it was computed with (`Utilities.record_inputs!`), so a deck run and a script run of the same workflow record the same deck. `Input/RawInputs/` stores the raw equilibrium/forcing/coil data and tolerance file — together they make `gpec.h5` a self-contained rerun snapshot (`Rerun.jl` reconstructs every control struct from them; the writer and rerun reader carry cross-reference comments marking the mirrored path pair). **Never echo TOML flags or control-struct values into any other group** — every group outside `Input/` is derived output. (The former `kinetic/` and `slayer/settings/` echoes were removed under this rule.) ## Schema @@ -62,7 +62,7 @@ Every dataset outside `Input/` (raw snapshot) must answer "what is this, in what Root-level file attributes: `schema_version` (currently `"2.0"`; bump on breaking schema changes — readers dispatch on it), `Conventions = "GPEC-HDF5-2.0"`, `references`, `title` (run description), `date_created` (ISO 8601 UTC). The code version stays in `Info/git_version`. -Mechanism: writers stay table-driven — each writer keeps a `path => (; long_name, units, dims, scale, attach)` table next to it (`scale` marks a coordinate dataset as a dimension scale; `attach` binds axes to scales — see the `Utilities.HDF5Annotations.annotate!` docstring) (`src/HDF5Schema.jl` for the main writer; alongside `write_galerkin!`, the PerturbedEquilibrium writer, `KineticForces/Output.jl`, and `Tearing/Runner/HDF5Output.jl` for the rest) and applies it post-write via `Utilities.HDF5Annotations.annotate!`. Entries for conditionally-written datasets are simply skipped when absent. When adding a dataset, add its table entry in the same commit — the schema test fails otherwise. +Mechanism: writers stay table-driven — each writer keeps a `path => (; long_name, units, dims, scale, attach)` table next to it (`scale` marks a coordinate dataset as a dimension scale; `attach` binds axes to scales — see the `Utilities.HDF5Annotations.annotate!` docstring) (`src/Utilities/HDF5Schema.jl` for the main writer; alongside `write_galerkin!`, the PerturbedEquilibrium writer, `KineticForces/Output.jl`, and `Tearing/HDF5Output.jl` for the rest) and applies it post-write via `Utilities.HDF5Annotations.annotate!`. Entries for conditionally-written datasets are simply skipped when absent. When adding a dataset, add its table entry in the same commit — the schema test fails otherwise. ## File-wide conventions diff --git a/docs/development/naming.md b/docs/development/naming.md index 16763ac6d2..30800cbb0a 100644 --- a/docs/development/naming.md +++ b/docs/development/naming.md @@ -56,11 +56,12 @@ Modules, with the abbreviation to use in prose: | Area | Prose | Source | |---|---|---| | `Analysis` | — | `src/Analysis/` | +| `Deck` | — | `src/Deck.jl` | | `Equilibrium` | EQUIL | `src/Equilibrium/` | | `ErrorFields` | EF | `src/ErrorFields/` | | `ForceFreeStates` | FFS | `src/ForceFreeStates/` | | `ForcingTerms` | FT | `src/ForcingTerms/` | -| `HDF5Schema` | — | `src/HDF5Schema.jl` | +| `HDF5Schema` | — | `src/Utilities/HDF5Schema.jl` | | `InnerLayer` | IL | `src/InnerLayer/` | | `KineticForces` | KF | `src/KineticForces/` | | `LocalStability` | LS | `src/LocalStability/` | @@ -73,7 +74,7 @@ Modules, with the abbreviation to use in prose: `EQUIL`, `VAC`, `FFS`, and `PE` were already in common use; the rest are new. Short module names take no abbreviation. -A submodule may be named where it helps triage: `ForceFreeStates.Galerkin`, `InnerLayer.GGJ`, `InnerLayer.SLAYER`, `Tearing.Dispersion`, `Tearing.Runner`. Only directories holding Julia sources qualify, so data directories are not Areas. +A submodule may be named where it helps triage: `ForceFreeStates.Galerkin`, `InnerLayer.GGJ`, `InnerLayer.SLAYER`, `Tearing.Dispersion`. Only directories holding Julia sources qualify, so data directories are not Areas. Areas outside `src/`: `Benchmarks`, `Build`, `CI`, `Docs`, `Examples`, `Regression`, `Repo`, `Test`. diff --git a/docs/development/toml-conventions.md b/docs/development/toml-conventions.md index 43bd03a08b..dfab01507f 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`, `WallSettings`, `ForceFreeStatesControl`, `ForcingTermsControl`, the `plasma_response`/`perturbed_equilibrium` docstrings, `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 @@ -48,7 +48,7 @@ fails when its pattern matches a violation) lint-check the covered TOMLs on comm - `toml-inline-annotations` — every `key = value` line carries a trailing `# description` (`examples/*` and `test/test_data/*` only, per the rule-2 exception above). - `toml-no-deprecated-keys` — no deprecated config keys; its key list mirrors the - `_DEPRECATED_*_KEYS` lists in `src/GeneralizedPerturbedEquilibrium.jl`, which stop any run that + `_DEPRECATED_*_KEYS` and `_MOVED_KEYS` lists in `src/Deck.jl`, which stop any run that sets one, so extend the hook's pattern whenever a key is added there. Run them manually with: diff --git a/docs/src/api.md b/docs/src/api.md index 854a2e274c..2c63384863 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -16,16 +16,36 @@ pe = perturbed_equilibrium(resp, rmp) ``` `solve(eq, alg; kwargs...)` is sugar for the canonical problem form: a -[`EulerLagrangeProblem`](@ref) names WHAT is being solved (the perturbed-plasma -Euler-Lagrange system posed on this equilibrium, with this mode range, wall and closure) +[`ForceBalanceProblem`](@ref) names WHAT is being solved (the perturbed-plasma +force-balance equation posed on this equilibrium, with this mode range, wall and closure) and the integrator names HOW. A `PlasmaEquilibrium` hosts many possible problems; the problem type keeps `solve` unambiguous as other problem classes appear. ```julia -prob = EulerLagrangeProblem(eq; nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true) +prob = ForceBalanceProblem(eq; nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true) ffs = solve(prob, Riccati()) +save_to_hdf5("gpec.h5", ffs) ``` +`solve` returns its result without writing a file. Call `save_to_hdf5(path, ffs)` +or pass an open HDF5 file or group to save it explicitly. + +The same `save_to_hdf5(target, result)` function saves an equilibrium, local-stability +result, perturbed-equilibrium state, kinetic-force state, tearing result, or error-field +result. Each module owns the datasets for its result. Paths open or create an HDF5 +file; an open file or group remains owned by the caller. Stage computations return +results without writing anything. Each save also records, in `Input/gpec_toml_raw`, the deck +sections stating the problem, algorithm and controls the result was computed with, so a script +run's file replays like a deck run's: `main(["gpec.h5"])`. A forcing no deck can state (a weighted +field, a spectrum or a sum) is replayed from the control-surface modes the save records. + +A deck runs through these same calls. `main` translates each section of `gpec.toml` into the +problem, algorithm or control it names, calls the function above for every section present, +and saves each result into the one `[Run] output_file`, so a deck can never take a path a script +cannot. It returns `(; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing)`, +each `nothing` when its section is absent; `ffs` is the result the later stages read (matched +when the deck has `[Match]`) and `match` its `MatchResult`. + The objects map one-to-one onto the pipeline stages: - [`PlasmaEquilibrium`](@ref) reads and processes the equilibrium — the `[Equilibrium]` @@ -41,11 +61,17 @@ The objects map one-to-one onto the pipeline stages: refined two-pass grid the deck uses; `Equilibrium.reform(eq, psi_nodes)` re-forms it on any grid. - `solve` runs the force-free-states stage — the `[ForceFreeStates]` section — and returns a `ForceFreeStatesResult`, exactly the object the TOML driver publishes. - - [`RMPField`](@ref) describes the external field — the `[ForcingTerms]` section — without + - [`RMPField`](@ref GeneralizedPerturbedEquilibrium.ForcingTerms.RMPField) describes the external field — the `[ForcingTerms]` section — without reading anything from disk until it is applied. - [`plasma_response`](@ref) builds the plasma-response and singular-coupling operators of a `ForceFreeStatesResult` once, and `perturbed_equilibrium(resp, rmp)` drives them with one field — together the `[PerturbedEquilibrium]` section. + - [`kinetic_torque`](@ref) evaluates the kinetic (NTV) torque of a perturbed equilibrium with the + model settings of a `KineticForcesControl` — the `[KineticForces]` section — and returns a + `KineticForcesResult`: `kf = kinetic_torque(pe, KineticForcesControl())`. + - [`error_fields`](@ref) assesses the error field of a set of coils on a plasma response — the + `[ErrorFields]` section — and returns an `ErrorFieldsResult`: + `ef = error_fields(resp, pe.coil_sets, ForcingTerms.CoilConfig(ft_ctrl); tolerances)`. ## Kinetic profiles @@ -88,11 +114,12 @@ with a failed conductor is `nominal - failed_coil`, not `0.9 * nominal`. ## Choosing an integrator The second argument of `solve` picks the formalism, and its fields are that formalism's -tunables. Everything else is a `ForceFreeStatesControl` keyword, so the TOML keys and the -`solve` keywords are the same knobs: +numerical settings. Physical inputs belong on `ForceBalanceProblem`; ODE tolerances and +storage settings belong on `Forward` or `Riccati`, and `extended_precision_bvp` belongs only +on `Riccati`. The current deck reader translates its existing section layout into those objects: ```julia -ffs = solve(eq, Forward(); nn=1, vac_flag=true) # dense ξ profiles for PerturbedEquilibrium +ffs = solve(eq, Forward(; eulerlagrange_tolerance=1e-8); nn=1, vac_flag=true) # dense ξ profiles for PerturbedEquilibrium ffs = solve(eq, Riccati(; nchunks=40); nn=1, vac_flag=true) # chunked propagators, Δ′ matrix ffs = solve(eq, Galerkin(; nx=512); nn=1) # RDCON outer-region Galerkin Δ′ ``` @@ -101,6 +128,11 @@ Which products each formalism can supply differs; a result carries `nothing` in its integrator does not produce and consumers warn and skip rather than erroring. See the [Stability Analysis](stability.md) page for the result struct and its capability gates. +`singfac_min`, `frobenius_order`, `frobenius_psi_max`, `truncate_at_dW_peak`, and +`dW_edge_scan_start` remain problem settings. `Galerkin` rejects non-default values of the +last three. Its per-surface order adjustment is controlled by `frobenius_order_ceiling`, +and its factorization is selected automatically. + ## Inner-layer matching Inner-layer matching is its own problem, posed on a finished solve: a `MatchProblem` holds @@ -127,38 +159,58 @@ Riccati (`vac_flag=true`) results; a Riccati-fed match fills `bpen` and the reso accepted — `GGJModel()` today; `SLAYERModel()` is slab-only and drives the free-eigenvalue tearing solve instead. +A deck poses the match with a `[Match]` section (`eta`, `rho`, `rotation`, `gamma`, `ideal`, +`inner_model = "GGJ"`), and `[GGJ]` holds the `GGJModel` keywords (`solver` and its options). + ## Tearing stability -The free-eigenvalue tearing solve is the second flavor of inner-layer matching: instead of -prescribing the layer rotation, a `TearingProblem` holds the outer Δ′ fixed and root-finds -the growth rate where the inner-layer response matches it. The same model slot applies — -`SLAYERModel()` is the slab layer that exists for exactly this problem. GGJ growth-rate -extraction is not implemented yet, so `GGJModel()` errors here. +A `TearingProblem` holds the outer solution and the surface coupling. The plasma is the +equilibrium's kinetic data (profiles, species, transport defaults, resistivity closure), the +inner-layer model holds its own closures, and `TearingConfig` holds the scan and root-search +settings. ```julia -ffs = solve(eq, Riccati(); nn=1, vac_flag=true) # Δ′ matrix for the dispersion -tear = solve(TearingProblem(ffs), SLAYERModel()) # coupling_mode=:coupled returns one global root instead -tear.gamma_Hz, tear.rational_q # root-found rate at each rational surface +tear = solve(TearingProblem(ffs), SLAYERModel(; P_perp_model=:chi_perp_e); + config=TearingConfig(; Q_re_range=(-2.0, 2.0), Q_im_range=(-0.5, 3.0))) +ggj = solve(TearingProblem(ffs; coupling_mode=:coupled), GGJModel(; solver=:ray); + config=TearingConfig(; gamma_range=(1.0, 1e4), omega_range=(-1e4, 1e4))) +save_to_hdf5("gpec.h5", tear) ``` -Keyword arguments of `TearingProblem` are the `[SLAYER]` deck section's procedure knobs -(scan mode and Q-domain, coupling mode, critical-Δ convention, extraction filters); -kinetic profiles come from the equilibrium the solve was posed on (see Kinetic profiles above). - -Kinetic runs (`kinetic_factor > 0`) take `Forward()`. The `"calculated"` source reads the -equilibrium's kinetic profiles, with the default `KineticForcesControl` species and collision -settings. +GGJ needs the raw Δ′ matrix and validates the surface ordering. Its η and ρ default from the +kinetic data: η by the resistivity closure SLAYER uses, ρ = μ·m_p·n_e as SLAYER takes it; pass `eta` +and `rho` to `TearingProblem` to override them. Both models report lab-frame frequencies with +per-surface E×B shifts: a coupled solve evaluates each layer in its own E×B frame (GGJ's s_k = s + i·n·Ω_E, +Johnson & Greene 1967 Appendix B), and an uncoupled root is moved from its layer's frame by n·Ω_E. Ω_E comes +from the kinetic data's `omega_E`, or from `omega_E_kHz` for the m/n it lists. Both expose +`omega` in rad/s and `gamma` in 1/s, positive for instability, and the HDF5 roots group records +`frequency_frame`. SLAYER searches a normalized Q box, GGJ a physical (γ, ω) box. + +A deck runs the solve with a `[Tearing]` section: the `TearingProblem` keywords and +`inner_model = "SLAYER"` or `"GGJ"`, with `[Tearing.search]` for the `TearingConfig` keys and +`[Tearing.CriticalResonantField]` for the penetration threshold. The model's keywords sit in +`[SLAYER]` or `[GGJ]`; `[SLAYER]` also gives the layer-overlap cap its closures. + +Kinetic runs take `Forward()`. The problem's `kinetic` keyword selects the source: +`kinetic=KineticForcesControl(; ...)` evaluates the KineticForces drift-kinetic +matrices with those settings on the equilibrium's kinetic profiles and species, and +`kinetic=TestKineticMatrices(; scale)` uses synthetic test matrices. `kinetic_factor` (default 1) +scales either. `km = kinetic_matrices(prob)` evaluates the source once, and `kinetic=km` reuses those +[`KineticMatrices`](@ref GeneralizedPerturbedEquilibrium.ForceFreeStates.KineticMatrices) in later solves, for example a +`kinetic_factor` scan. A deck selects the KineticForces model with `[ForceFreeStates] kinetic = "KineticForces"`, using +its `[KineticForces]` section. ## Local stability The Mercier, resistive interchange and ballooning profiles are their own stage, posed on the -equilibrium alone. A [`LocalStabilityProblem`](@ref) holds the equilibrium, the -[`Ballooning`](@ref) algorithm holds the scan knobs, and `solve` returns a -[`LocalStabilityResult`](@ref) with the profiles as named fields. A deck runs the same solve -through its `[LocalStability]` section, whose keys are the `Ballooning` keywords. +equilibrium. A [`LocalStabilityProblem`](@ref) holds the equilibrium, ballooning angle +`theta_k`, and `alpha_boundary_scan` choice. [`Ballooning`](@ref) holds numerical tolerances +and scan settings. `solve` returns a [`LocalStabilityResult`](@ref); its `verbose` keyword +defaults to `false`. A deck states the problem in `[LocalStability]` and the numerics in `[Ballooning]`. ```julia -ls = solve(LocalStabilityProblem(eq), Ballooning(; n_scan=64)) +problem = LocalStabilityProblem(eq; theta_k=0.0, alpha_boundary_scan=true) +ls = solve(problem, Ballooning(; n_scan=64); verbose=false) ls.psi, ls.D_I, ls.D_R, ls.ballooning_Delta_prime # profiles on the equilibrium ψ grid ls.ballooning_psi, ls.alpha, ls.alpha_critical # first ballooning boundary, NaN outside its scan window ``` @@ -166,8 +218,9 @@ ls.ballooning_psi, ls.alpha, ls.alpha_critical # first ballooning boundar ## Entry points ```@docs -GeneralizedPerturbedEquilibrium.EulerLagrangeProblem +GeneralizedPerturbedEquilibrium.ForceBalanceProblem GeneralizedPerturbedEquilibrium.solve +GeneralizedPerturbedEquilibrium.kinetic_matrices GeneralizedPerturbedEquilibrium.perturbed_equilibrium ``` diff --git a/docs/src/conventions.md b/docs/src/conventions.md index d0e910679c..f0a97b2777 100644 --- a/docs/src/conventions.md +++ b/docs/src/conventions.md @@ -135,7 +135,7 @@ q = \frac{F}{2\pi}\oint \frac{J}{B_p}\,dl_p . ### Toroidal Mode Number ``n`` -- Set as a **positive integer** via `nn_low`/`nn_high` in `gpec.toml`. ``n < 1`` is not supported +- Set as a **positive integer** via `[ForceFreeStates] nn` in `gpec.toml`. ``n < 1`` is not supported and is clamped to 1 with a warning (`src/GeneralizedPerturbedEquilibrium.jl`). Results for multiple ``n`` can be superposed. diff --git a/docs/src/developer_notes.md b/docs/src/developer_notes.md index 416cf86a25..f25c828e1a 100644 --- a/docs/src/developer_notes.md +++ b/docs/src/developer_notes.md @@ -31,6 +31,16 @@ regenerate-only-when-`depends`-change rule — is in The `Ballooning and Mercier Local Stability` and `Inner Layer Module` pages are the reference exemplars. +## Deck and API parity + +Every capability must be reachable both from a `gpec.toml` deck (through `main`) and from the +scripting API (`solve` and the problem and result types). A feature that exists on only one path +is not finished. A new feature, deck key, or API keyword also adds a case to the `"deck and API parity"` +testset in `test/runtests_solve_api.jl`. That case runs the workflow once from a deck and once +from the same inputs through the API, and checks that the two results agree. A known drift +between the paths is marked `@test_broken`, not dropped, so the change that fixes it shows up +as an unexpected pass. + ## Regression Testing The regression harness **must be run on every pull request before merging into `develop`**. It is the project's primary safeguard for tracking how numerical results evolve across changes, so it is only useful if every PR exercises it. When you open a PR, paste the regression report into the PR thread so reviewers can see what moved (and what did not). If your change touches a quantity that is not yet tracked, add a new regression case — or extend an existing one — in the same PR. diff --git a/docs/src/equilibrium.md b/docs/src/equilibrium.md index 99f9584258..a7215f1de1 100644 --- a/docs/src/equilibrium.md +++ b/docs/src/equilibrium.md @@ -57,7 +57,7 @@ tuning is needed. ## Radial grid packing -With `grid_type = "auto"` and `mpsi = 0` (the defaults; `"log_asymptotic"` is a legacy alias), the radial grid is +With `grid_type = "auto"` and `mpsi = 0` (the defaults), the radial grid is built by a **two-pass measured-curvature refinement** driven by the single accuracy knob `psi_accuracy` (τ): @@ -131,9 +131,10 @@ Modules = [GeneralizedPerturbedEquilibrium.Equilibrium] - `PlasmaEquilibrium` — the runtime structure containing spline fields, geometry, profiles, and computed diagnostics (q-profile, separatrix, etc.). -- `LargeAspectRatioConfig`, `SolovevConfig` — convenience structs to - construct analytic/model equilibria when using `eq_type = "lar"` or - `eq_type = "sol"`. +- `SolovevConfig`, `LargeAspectRatioConfig`, `TJAnalyticConfig` — the parameters of the + analytic equilibria. A deck states one as a `[Solovev]`, `[LargeAspectRatio]` or + `[TJAnalytic]` section, which selects the equilibrium, so it sets no `eq_type`. + `[TJAnalytic] solver = "inverse" | "direct"` picks the TJ pipeline. ## Key functions @@ -174,7 +175,7 @@ Analytic / testing example: construct a large-aspect-ratio model using GeneralizedPerturbedEquilibrium larcfg = GeneralizedPerturbedEquilibrium.Equilibrium.LargeAspectRatioConfig(r0=10.0, a=1.0, beta0=1e-3) -pe = PlasmaEquilibrium(larcfg) # eq_type = "lar" is implied by the config type +pe = PlasmaEquilibrium(larcfg) # the config type implies the equilibrium kind println("Built LAR equilibrium with a = ", larcfg.a) ``` diff --git a/docs/src/error_fields.md b/docs/src/error_fields.md index 66862d8c01..993b4a7471 100644 --- a/docs/src/error_fields.md +++ b/docs/src/error_fields.md @@ -46,8 +46,6 @@ coupling: fd_step_shift_m = 1e-3 # Central-difference step for the rigid shifts [m] fd_step_tilt_deg = 0.1 # Central-difference step for the rigid tilts [degrees] rotation_center = "conductor" # Tilt pivot: each conductor's own centre ("conductor") or the whole set's ("set") -write_outputs_to_HDF5 = true # Write ErrorFields/CoilSensitivities/ to the output file -verbose = false # Log per-coil-set progress and linearity diagnostics exclude_coils = [] # Coil sets that are not error-field sources (the NTV efc_coils are excluded automatically) ``` @@ -58,6 +56,18 @@ curvature diagnostic per tap, the current pattern the spectra were evaluated at, its shift and tilt sensitivities, their direction-averaged in-plane magnitudes, and the in-plane shift and tilt that would cancel `delta_as_designed`). +From a script, the same stage is one call on the plasma response and the coil sets of a coil-forced +perturbed equilibrium; it returns an `ErrorFieldsResult`, and products not asked for are `nothing`: + +```julia +resp = plasma_response(ffs) +pe = perturbed_equilibrium(resp, RMPField(ft_ctrl)) +ef = error_fields(resp, pe.coil_sets, ForcingTerms.CoilConfig(ft_ctrl); + tolerances=ErrorFields.read_tolerance_toml("tolerances.toml"), + scenario=ErrorFields.ScenarioParameters(ffs.equil; n_e=5.0)) +save_to_hdf5("gpec.h5", ef) +``` + `examples/DIIID-like_error_field_example/` is a complete case: the DIII-D-like equilibrium with the C-coil as the nominal n = 1 source and the eighteen DIII-D F coils added as single-filament hoops at the centroids of their winding packs (from OpenFUSIONToolkit's TokaMaker @@ -76,7 +86,7 @@ either way. Manufacturing tolerances are engineering data, reviewed and versioned independently of any run, so they live in their own TOML file named by `tolerance_file` in `[ErrorFields]` (a path -relative to the run directory). The run validates the file, echoes its text into +relative to the deck directory). The run validates the file, echoes its text into `Input/RawInputs/ErrorFields/tolerance_toml_raw` for replay, and leaves the device numbers where they belong: outside the repository. The sampling and Monte Carlo stages that consume these tolerances follow in later releases; this release fixes the format. @@ -192,7 +202,7 @@ error. `without_coils(table, names)` is the same operation on a table in memory, ```toml [ErrorFields] -tolerance_file = "tolerances.toml" # Manufacturing-tolerance TOML, relative to the run directory +tolerance_file = "tolerances.toml" # Manufacturing-tolerance TOML, relative to the deck directory [ErrorFields.MonteCarlo] nsample = 1000000 # Samples per batch @@ -441,7 +451,7 @@ table.tilt_sensitivity_per_deg[1, :] # ∂δ/∂θx per degree, pe # In memory, from the run's returned state (no I/O) rc = PerturbedEquilibrium.ResonantCoupling(run.pe.response) dom = PerturbedEquilibrium.dominant_coupling(rc; psi_low=0.5) -table = EF.sensitivity_table(run.coil_sensitivities, dom; mode=1) +table = EF.sensitivity_table(run.error_fields.sensitivities, dom; mode=1) ``` ### How far the linear model holds diff --git a/docs/src/galerkin.md b/docs/src/galerkin.md index e488db3793..05ecabf956 100644 --- a/docs/src/galerkin.md +++ b/docs/src/galerkin.md @@ -19,20 +19,19 @@ The solve produces the inter-surface Δ′ matrix and the PEST-3 matching blocks whichever formalism produced it. These are the outer-region inputs to resistive matched-asymptotic stability analysis [Glasser 2016, Phys. Plasmas **23**, 072505]. -Select it with `integrator = "galerkin"` in `[ForceFreeStates]`. It replaces the radial ODE +Select it with `integrator = "Galerkin"` in `[ForceFreeStates]`; its settings go in `[Galerkin]`. It replaces the radial ODE integration rather than supplementing it: the run computes its own vacuum response at the control surface (when `vac_flag`) and produces no free-boundary energies or ODE trace. -Setting `gal_match_flag` additionally matches the inner layer, giving a driven ξ solution. The -composite inner-region profiles are a scripting diagnostic: they need `Galerkin(; cut_solution=true)` -followed by `solve(MatchProblem(...))`, and a deck run does not write them. +A `[Match]` section additionally matches the inner layer, giving a driven ξ solution. The +composite inner-region profiles need the cut solution, `[Galerkin] cut_solution = true` (or +`Galerkin(; cut_solution=true)` followed by `solve(MatchProblem(...))`). -`gal_solver` picks the banded solver. `"auto"` (the default) tries a Hermitian Cholesky +The banded solve automatically tries a Hermitian Cholesky factorization of the lower band (LAPACK `zpbtrf`/`zpbtrs`, with the Dirichlet rows' columns moved to the right-hand side). The free block's off-diagonal entries must be Hermitian within `1e-10 * sqrt(abs(A_ii)) * sqrt(abs(A_jj))`; the diagonal's imaginary part is deliberately discarded. Unequal bandwidths, off-diagonal asymmetry, or failure of positive definiteness -trigger a reported fallback to Jacobi-scaled LU. `"LU"` skips the Cholesky attempt; -`"cholesky"` stops with an error instead of falling back. The factorization actually used is +trigger a reported fallback to Jacobi-scaled LU. The factorization actually used is recorded as `ForceFreeStates/Solutions/GalerkinIntegration/solver` in `gpec.h5`. The implementation lives in `src/ForceFreeStates/Galerkin/`: diff --git a/docs/src/inner_layer.md b/docs/src/inner_layer.md index fa459f4a02..dc8d93fa90 100644 --- a/docs/src/inner_layer.md +++ b/docs/src/inner_layer.md @@ -4,12 +4,11 @@ The `Tearing` module groups the resistive tearing-mode analysis stack: `InnerLayer` (per-surface inner-layer matching data Δ(Q) for the GGJ and SLAYER models), `Dispersion` (physics-agnostic complex-plane scan and contour-intersection root extraction), `CriticalResonantField` (torque-balance -error-field penetration threshold), and `Runner` (user-facing TOML -configuration, profile loading, and HDF5 output). +error-field penetration threshold), and the tearing solve itself, +`solve(TearingProblem(ffs; ...), model; config=TearingConfig(...))`, with its +`TearingResult` and HDF5 output. -## Layer Inputs - -Equilibrium/ForceFreeStates glue that assembles per-surface inner-layer inputs. +## Tearing solve and layer inputs ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.Tearing] @@ -204,7 +203,7 @@ res.Δ, res.breaks, res.bc_cond # res.Δ = (Δ₊, Δ₋) = (interchange, t ## Critical resonant field -With `[SLAYER.CriticalResonantField] enabled = true`, each SLAYER surface also gets the +With a `[Tearing.CriticalResonantField]` section (`TearingProblem(ffs; critical_resonant_field=CriticalResonantFieldControl())`), each SLAYER surface also gets the critical resonant field for error-field penetration, from the steady-state torque balance of Cole and Fitzpatrick, Phys. Plasmas **13**, 032503 (2006). The electromagnetic torque on the layer balances the viscous restoring torque when (Cole Eq. 61) @@ -235,8 +234,8 @@ and balance is lost above (Cole Eq. 62) modes (Fortran fixes α = 10⁻²). Cole and Fitzpatrick have no Δ_crit; including it treats the curvature offset as independent of frequency. It is zero unless `D_c_type` is set. Cole's threshold needs a tearing-stable surface, so a surface with α ≤ 0 is skipped with a warning. -- **Q axis.** Q is on Cole's axis, with the diamagnetic poles at Q_e and Q_i; `solve_inner` - mirrors the real axis, so the scan uses conj(Δ(−Q)). +- **Q axis.** Q is on Cole's axis, with the diamagnetic poles at Q_e and Q_i; it runs opposite to the dispersion + scans' Re Q, so the scan uses the layer response at −Q. - **Threshold.** Penetration is associated with a saddle-node bifurcation (Fitzpatrick, Nucl. Fusion 33, 1049 (1993), Sec. 7.4): as the field grows, Q moves from Q0 toward the electron fluid at Q_e, where the torque vanishes, and the threshold is the first resolved saddle-node bifurcation on that path, above @@ -300,7 +299,7 @@ Modules = [GeneralizedPerturbedEquilibrium.InnerLayer.SLAYER] ## Dispersion ```@autodocs -Modules = [GeneralizedPerturbedEquilibrium.Dispersion] +Modules = [GeneralizedPerturbedEquilibrium.Tearing.Dispersion] ``` ## CriticalResonantField @@ -308,9 +307,3 @@ Modules = [GeneralizedPerturbedEquilibrium.Dispersion] ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.Tearing.CriticalResonantField] ``` - -## Runner - -```@autodocs -Modules = [GeneralizedPerturbedEquilibrium.Runner] -``` diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md index 0ef98bdb28..091492b543 100644 --- a/docs/src/kinetic_forces.md +++ b/docs/src/kinetic_forces.md @@ -29,10 +29,24 @@ extension: The NTV calculation consumes `n_i, n_e, T_i, T_e, omega_E`; `chi_e`/`chi_i`/`chi_phi`/`zeff`, when present, are carried for the resistive-layer (SLAYER) analysis and ignored here. That analysis takes ``Z_\mathrm{eff}(\psi)`` from the `zeff` dataset when the file -has one, otherwise from `n_i/n_e` with the `[SLAYER]` `impurity_charge`, otherwise -from the scalar `[SLAYER] zeff`. Supply `zeff` whenever the impurity content is not a +has one, otherwise from `n_i/n_e` with the kinetic data's `impurity_charge`, otherwise +from its `zeff_default`. Supply `zeff` whenever the impurity content is not a single species fixed by `n_i/n_e`. +The plasma's species composition (main ions, impurity charge and mass), the defaults that stand +in for absent χ and Z_eff profiles, the confinement time and the resistivity closure belong to the +equilibrium's kinetic data, not to any one physics stage, so KineticForces and the tearing layers +see one plasma. In a script, set them where the profiles are attached: + +```julia +attach_kinetic_profiles!(equil, "kinetic.h5"; ion_species=[IonSpecies(; z=1, m=2, fraction=0.5), + IonSpecies(; z=1, m=3, fraction=0.5)], impurity_charge=6, impurity_mass=12) +``` + +or copy loaded data with new settings, `KineticProfileData(equil.kinetic; impurity_charge=6)`. +The default is one deuterium main ion carrying all of `n_i`, with a carbon impurity. In a deck +they are `[Equilibrium]` keys: `[[Equilibrium.ion_species]]`, `impurity_charge` and `impurity_mass`. + ## Multi-ion runs A plasma may declare an arbitrary list of main-ion species; the NTV is computed @@ -40,34 +54,38 @@ per species under one shared full-composition ``Z_\mathrm{eff}`` and summed (``\tau = \sum_s \tau_s``) over the main ions, the quasineutrality-closing impurity, and (with `electron = true`) the electrons. This applies to both NTV paths: the post-PE ψ-quadrature diagnostic and the self-consistent -`kinetic_source = "calculated"` matrices. +kinetic matrices of a `kinetic = "KineticForces"` force-balance solve. ```toml [Equilibrium] kinetic_file = "kinetic.h5" # n_i column/dataset = TOTAL main-ion density - -[KineticForces] -electron = true # add electron NTV in addition to the ion species impurity_charge = 6 # impurity charge (closes quasineutrality) impurity_mass = 12 # impurity mass -[[KineticForces.ion_species]] +[[Equilibrium.ion_species]] z = 1 m = 2 fraction = 0.5 # this species' share of the total n_i -[[KineticForces.ion_species]] +[[Equilibrium.ion_species]] z = 1 m = 3 density = "n_T" # or: explicit n_* dataset from the HDF5 kinetic file + +[KineticForces] +electron = true # add electron NTV in addition to the ion species ``` +The composition describes the plasma, so it sits on the equilibrium's kinetic data beside +`kinetic_file`, where every kinetic consumer (NTV, the tearing layers) reads it. + Each species sets exactly one of `fraction` (share of the file's total `n_i`) or `density` (a named `n_*` dataset). An all-fraction list must sum to 1; in a mixed list the impurity content is set by the file's `n_i`/`n_e` deficit, not by a fraction shortfall, and fractions may sum below (never above) 1. Every -main-ion charge must satisfy `z < impurity_charge`. An empty `ion_species` list runs the -single main ion from `ion_charge`/`ion_mass` (with `electron = true` still adding the -electron species — the electron flag always means *in addition to* the ions). +main-ion charge must satisfy `z < impurity_charge`. A single main ion carrying all of `n_i` +runs the single-species NTV (with `electron = true` still adding the electron species — the +electron flag, a `KineticForcesControl` setting, always means *in addition to* the ions). Without +an `ion_species` list the plasma is one deuterium main ion. The summed total is written to `KineticForces//` exactly as in a single-ion run, and each species' own contribution to `KineticForces/PerSpecies/