From 1dc6138aad217c3e5e71cc9f071924607303c831 Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Mon, 12 Oct 2026 16:46:23 -0400 Subject: [PATCH 1/3] =?UTF-8?q?ForceFreeStates=20-=20BUGFIX!=20-=20Trigger?= =?UTF-8?q?=20the=20Gaussian=20reduction=20on=20the=20conditioning=20of=20?= =?UTF-8?q?U=E2=82=81=20and=20factor=20U=E2=82=81=20in=20the=20plasma=20re?= =?UTF-8?q?sponse=20matrix?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The forward integrator reduced only when the U₁ column-norm spread exceeded ucrit, which misses columns that grow alike and turn near-parallel. It now reduces when the 1-norm condition estimate of the column-normalized U₁ exceeds ucrit (ucrit keeps its entry-size meaning in the Riccati path). The reduction storage grows on demand instead of overflowing. wp = U₂ / U₁ factored U₁ᵀ, so partial pivoting ran across U₁'s columns, whose norms span ~1e110 right after a reduction; it now factors U₁ itself. Co-Authored-By: Claude Opus 5.5 --- examples/DIIID-like_SLAYER_example/gpec.toml | 2 +- .../DIIID-like_error_field_example/gpec.toml | 2 +- .../gpec.toml | 2 +- .../gpec.toml | 2 +- examples/DIIID-like_ideal_example/gpec.toml | 2 +- .../gpec.toml | 2 +- .../gpec.toml | 2 +- examples/LAR_beta_scan/gpec.toml | 2 +- examples/LAR_epsilon_scan/gpec.toml | 2 +- examples/LAR_ideal_match_test/gpec.toml | 2 +- examples/LAR_resistive_match_test/gpec.toml | 2 +- examples/Solovev_ideal_example/gpec.toml | 2 +- examples/Solovev_ideal_example_3D/gpec.toml | 2 +- .../Solovev_ideal_example_multi_n/gpec.toml | 2 +- .../single_n_1/gpec.toml | 2 +- .../single_n_2/gpec.toml | 2 +- .../Solovev_kinetic_NTV_example/gpec.toml | 2 +- .../gpec.toml | 2 +- .../gpec.toml | 2 +- examples/a10_kinetic_example/gpec.toml | 2 +- src/ForceFreeStates/CoreTypes.jl | 4 +- src/ForceFreeStates/EulerLagrange.jl | 47 +++++++++---------- src/ForceFreeStates/Free.jl | 8 ++-- src/ForceFreeStates/Result.jl | 2 +- src/ForceFreeStates/Utils.jl | 15 ++++++ test/runtests_eulerlagrange.jl | 40 ++++++++++++---- test/runtests_fullruns.jl | 4 +- .../gpec.toml | 2 +- .../gpec.toml | 2 +- .../gpec.toml | 2 +- .../gpec.toml | 2 +- .../gpec.toml | 2 +- .../gpec.toml | 2 +- 33 files changed, 103 insertions(+), 69 deletions(-) diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml index a8a99f4c7..8081fa30c 100644 --- a/examples/DIIID-like_SLAYER_example/gpec.toml +++ b/examples/DIIID-like_SLAYER_example/gpec.toml @@ -49,7 +49,7 @@ kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 eulerlagrange_tolerance = 1e-10# Relative tolerance for ODE integration of Euler-Lagrange equations save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) integrator = "riccati" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) diff --git a/examples/DIIID-like_error_field_example/gpec.toml b/examples/DIIID-like_error_field_example/gpec.toml index 9a35ba9d0..dace63f5e 100644 --- a/examples/DIIID-like_error_field_example/gpec.toml +++ b/examples/DIIID-like_error_field_example/gpec.toml @@ -52,7 +52,7 @@ kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 eulerlagrange_tolerance = 1e-10 # Relative tolerance for ODE integration of Euler-Lagrange equations save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) diff --git a/examples/DIIID-like_gal_resistive_example/gpec.toml b/examples/DIIID-like_gal_resistive_example/gpec.toml index 71d16afcb..ef1b9d65b 100644 --- a/examples/DIIID-like_gal_resistive_example/gpec.toml +++ b/examples/DIIID-like_gal_resistive_example/gpec.toml @@ -48,7 +48,7 @@ kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 eulerlagrange_tolerance = 1e-10# Relative tolerance for ODE integration of Euler-Lagrange equations save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) integrator = "galerkin" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) diff --git a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml index 71e584905..7f5aa37dc 100644 --- a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml +++ b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml @@ -48,7 +48,7 @@ kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 eulerlagrange_tolerance = 1e-10# Relative tolerance for ODE integration of Euler-Lagrange equations save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) integrator = "galerkin" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml index 2aca73d8a..394e0f4da 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/examples/DIIID-like_ideal_example/gpec.toml @@ -47,7 +47,7 @@ kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 eulerlagrange_tolerance = 1e-10 # Relative tolerance for ODE integration of Euler-Lagrange equations save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) diff --git a/examples/DIIID-like_ideal_multi_n_example/gpec.toml b/examples/DIIID-like_ideal_multi_n_example/gpec.toml index 002596c21..d5c53a986 100644 --- a/examples/DIIID-like_ideal_multi_n_example/gpec.toml +++ b/examples/DIIID-like_ideal_multi_n_example/gpec.toml @@ -49,7 +49,7 @@ kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 eulerlagrange_tolerance = 1e-10 # Relative tolerance for ODE integration of Euler-Lagrange equations save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) diff --git a/examples/DIIID-like_riccati_deltaprime_example/gpec.toml b/examples/DIIID-like_riccati_deltaprime_example/gpec.toml index 99a6f80c9..1e7bb9917 100644 --- a/examples/DIIID-like_riccati_deltaprime_example/gpec.toml +++ b/examples/DIIID-like_riccati_deltaprime_example/gpec.toml @@ -49,7 +49,7 @@ kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 eulerlagrange_tolerance = 1e-10 # Relative tolerance for ODE integration of Euler-Lagrange equations save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) integrator = "riccati" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) diff --git a/examples/LAR_beta_scan/gpec.toml b/examples/LAR_beta_scan/gpec.toml index 88770631b..a40afbd75 100644 --- a/examples/LAR_beta_scan/gpec.toml +++ b/examples/LAR_beta_scan/gpec.toml @@ -49,7 +49,7 @@ mthvac = 960 # Number of points used in splines over poloidal eulerlagrange_tolerance = 1e-12 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization sing_order = 6 # Order of the singular-surface (Frobenius) series expansion integrator = "riccati" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) diff --git a/examples/LAR_epsilon_scan/gpec.toml b/examples/LAR_epsilon_scan/gpec.toml index 393d093b6..c19e4afb0 100644 --- a/examples/LAR_epsilon_scan/gpec.toml +++ b/examples/LAR_epsilon_scan/gpec.toml @@ -50,7 +50,7 @@ mthvac = 960 # Number of points used in splines over poloidal eulerlagrange_tolerance = 1e-12 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization sing_order = 6 # Order of the singular-surface (Frobenius) series expansion integrator = "riccati" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) diff --git a/examples/LAR_ideal_match_test/gpec.toml b/examples/LAR_ideal_match_test/gpec.toml index 0706863b7..0daae5d79 100644 --- a/examples/LAR_ideal_match_test/gpec.toml +++ b/examples/LAR_ideal_match_test/gpec.toml @@ -46,7 +46,7 @@ mthvac = 960 # Number of points used in splines over poloidal eulerlagrange_tolerance = 1e-12 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization sing_order = 6 # Power-series order for the singular-surface asymptotics save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. diff --git a/examples/LAR_resistive_match_test/gpec.toml b/examples/LAR_resistive_match_test/gpec.toml index 0571301c7..c9a595c0d 100644 --- a/examples/LAR_resistive_match_test/gpec.toml +++ b/examples/LAR_resistive_match_test/gpec.toml @@ -47,7 +47,7 @@ mthvac = 960 # Number of points used in splines over poloidal eulerlagrange_tolerance = 1e-12 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization sing_order = 6 # Power-series order for the singular-surface asymptotics save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. diff --git a/examples/Solovev_ideal_example/gpec.toml b/examples/Solovev_ideal_example/gpec.toml index f74d26fa6..44d104f2b 100644 --- a/examples/Solovev_ideal_example/gpec.toml +++ b/examples/Solovev_ideal_example/gpec.toml @@ -68,7 +68,7 @@ kinetic_source = "fixed" # Kinetic matrix source: "fixed" test matrices, or kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) diff --git a/examples/Solovev_ideal_example_3D/gpec.toml b/examples/Solovev_ideal_example_3D/gpec.toml index a5c3e1943..53ae020d7 100644 --- a/examples/Solovev_ideal_example_3D/gpec.toml +++ b/examples/Solovev_ideal_example_3D/gpec.toml @@ -34,7 +34,7 @@ kinetic_source = "fixed" # Kinetic matrix source: "fixed" test matrices, or kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) diff --git a/examples/Solovev_ideal_example_multi_n/gpec.toml b/examples/Solovev_ideal_example_multi_n/gpec.toml index 7670bd5cf..f51fbd8a3 100644 --- a/examples/Solovev_ideal_example_multi_n/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/gpec.toml @@ -45,7 +45,7 @@ kinetic_source = "fixed" # Kinetic matrix source: "fixed" test matrices, or kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details). # The multi-n Δ' matrix has open issues and is skipped with a warning; the energies stay valid. diff --git a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml index 59760b231..58de244cf 100644 --- a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml @@ -44,7 +44,7 @@ kinetic_source = "fixed" # Kinetic matrix source: "fixed" test matrices, or kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization write_outputs_to_HDF5 = true # Write outputs to the HDF5 file HDF5_filename = "euler_n1.h5" # Output HDF5 filename # Solovev analytic equilibrium parameters (eq_type = "sol"); see SolovevConfig in src/Equilibrium. diff --git a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml index a80ac7741..ae19cb3a7 100644 --- a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml @@ -44,7 +44,7 @@ kinetic_source = "fixed" # Kinetic matrix source: "fixed" test matrices, or kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization write_outputs_to_HDF5 = true # Write outputs to the HDF5 file HDF5_filename = "euler_n2.h5" # Output HDF5 filename # Solovev analytic equilibrium parameters (eq_type = "sol"); see SolovevConfig in src/Equilibrium. diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml index 833a5ba22..5264cd928 100644 --- a/examples/Solovev_kinetic_NTV_example/gpec.toml +++ b/examples/Solovev_kinetic_NTV_example/gpec.toml @@ -65,7 +65,7 @@ kinetic_source = "fixed" # Kinetic matrix source: "fixed" test matrices, or kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) diff --git a/examples/Solovev_kinetic_calculated_example/gpec.toml b/examples/Solovev_kinetic_calculated_example/gpec.toml index 24824aa7c..873ae4a03 100644 --- a/examples/Solovev_kinetic_calculated_example/gpec.toml +++ b/examples/Solovev_kinetic_calculated_example/gpec.toml @@ -61,7 +61,7 @@ kinetic_factor = 1.0 # Full-strength kinetic matrices (the "calculated" integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization write_outputs_to_HDF5 = true # Write outputs to the HDF5 file verbose = false # Enable verbose logging diff --git a/examples/Solovev_kinetic_multiion_example/gpec.toml b/examples/Solovev_kinetic_multiion_example/gpec.toml index b90e628ee..2999441f1 100644 --- a/examples/Solovev_kinetic_multiion_example/gpec.toml +++ b/examples/Solovev_kinetic_multiion_example/gpec.toml @@ -66,7 +66,7 @@ kinetic_source = "fixed" # Kinetic matrix source: "fixed" test matrices, or kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. # Integration-range truncation (see ForceFreeStatesControl docstring for details) diff --git a/examples/a10_kinetic_example/gpec.toml b/examples/a10_kinetic_example/gpec.toml index eb32a2b1e..bde1ad3b3 100644 --- a/examples/a10_kinetic_example/gpec.toml +++ b/examples/a10_kinetic_example/gpec.toml @@ -42,7 +42,7 @@ kinetic_factor = 1.0 # Scaling of kinetic matrices (0 = ideal path; >0 integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization write_outputs_to_HDF5 = true # Write outputs to the HDF5 file force_termination = true # Stop after force-free states — the kinetic stability eigenvalues are the result of interest diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index 5544b0e60..f3edec6f8 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/src/ForceFreeStates/CoreTypes.jl @@ -131,9 +131,9 @@ gpec.toml. - `nstep::Int` - Maximum number of integration steps (not yet implemented) - `ksing::Int` - Singular surface handling parameter - `eulerlagrange_tolerance::Float64` - Relative tolerance for ODE integration of Euler-Lagrange equations - - `ucrit::Float64` - Critical value of unorm ratio to trigger solution normalization. In the standard path it triggers Gaussian reduction; in the Riccati path it triggers `renormalize_riccati_inplace!`. Default `1e3`. The trigger compares column-norm growth, not conditioning, so near-parallel columns can leave an error in the free-boundary energies that does not shrink with `eulerlagrange_tolerance`. + - `ucrit::Float64` - Normalization threshold, default `1e3`. Forward path: the condition number of the column-normalized U₁ above which a Gaussian reduction is applied; the reduction pivots on U₁, so its rounding error and that of the reconstructed `u_store` scale as eps·`ucrit` (values ≥ 1e10 lose digits). Riccati path: the largest entry of U₁ or U₂ above which `renormalize_riccati_inplace!` runs. - `numsteps_init::Int` - Initial array size for ODE data storage - - `numunorms_init::Int` - Initial array size for solution normalization data + - `numunorms_init::Int` - Initial capacity of the Gaussian-reduction storage; it grows as needed - `singfac_min::Float64` - Fractional distance from rational q at which ideal jump condition is enforced - `set_psilim_via_dmlim::Bool` - Truncate the integration domain at `(last_rational_q + dmlim) / n` rather than at `qhigh` / `psihigh`. Fortran STRIDE found that truncating ~20 % above the outermost rational (`dmlim = 0.2`) avoids a numerical kink instability in δW that appears when the integration ends too close to or just below a rational surface. **For diverted equilibria where q → ∞ at the separatrix** (e.g. DIII-D geqdsks, the bulk of production use) this costs negligible physical domain because rationals get arbitrarily dense near the LCFS — `set_psilim_via_dmlim = true` is the safe and recommended default. **For limited circular / analytical equilibria with finite q at the edge** (Solovev, LAR scans), rationals are sparse and 20 % above the last rational chops off too much edge, so set `set_psilim_via_dmlim = false` and let `qhigh` / `psihigh` control the truncation. Multi-`n` runs are not supported by this truncation (the "outermost rational + dmlim / n" depends on which `n`); when `set_psilim_via_dmlim = true` with `nn_low != nn_high`, `sing_lim!` warns and falls back to `qhigh` / `psihigh`. Default `true`. - `dmlim::Float64` - Distance beyond last rational surface (normalised ∈ [0,1) in units of 1/n). Only used when `set_psilim_via_dmlim` is true. Fortran STRIDE convention is 0.2 (truncate 20 % of one rational-surface spacing above the last surface), retained here. diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 9a6cb836a..12ab6dfb6 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -44,7 +44,7 @@ and a small set of temporary matrices and factors used to compute singular-layer - `numpert_total::Int` - Total number of Fourier mode combinations (m × n) used in the calculation. - - `numunorms_init::Int` - Initial allocation size for the number of normalization operations recorded. + - `numunorms_init::Int` - Initial capacity for recorded Gaussian reductions. - `msing::Int` - Number of singular surfaces in the equilibrium (used to size asymptotic coefficient arrays). @@ -376,7 +376,7 @@ end forward_eulerlagrange_integration(ctrl, equil, mats, intr; verbose=ctrl.verbose) -> (odet, nothing, nothing, nothing) Forward branch of [`eulerlagrange_integration`](@ref): integrates chunk by chunk from the axis, -applying Gaussian reduction whenever a solution norm ratio exceeds `ctrl.ucrit` and undoing it +applying Gaussian reduction whenever the column-normalized U₁ is worse conditioned than `ctrl.ucrit` and undoing it via `transform_u!` at the end, so `odet.u_store` comes back dense in the axis basis. Call directly to force this branch regardless of `ctrl.integrator`; `verbose` overrides `ctrl.verbose` for progress logging. @@ -987,13 +987,11 @@ end """ compute_solution_norms!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal, sing_flag::Bool; resonant_rows=Int[]) -Computes norms of the solution vectors of the array `u` and normalizes them -if this is not the first call after a fixup. Formerly `ode_unorm!`. -Throws an error if any vector norm is zero. It then compares the variation in norms -relative to initial values after a fixup, and applies the Gaussian reduction via -`apply_gaussian_reduction!` if the variation exceeds `ctrl.ucrit` or if `sing_flag` is true. -Performs the same function as `ode_unorm` in the Fortran code, with minor differences in indexing -and array handling. +Computes the U₁ column norms of `u`, relative to the last reduction unless this is the first call after one +(they set the reduction's pivot order). Formerly `ode_unorm!`. Throws an error if any norm is zero. Applies +`apply_gaussian_reduction!` at a singular surface or when the condition number of the column-normalized U₁ +exceeds `ctrl.ucrit`. Fortran `ode_unorm` triggers on the column-norm spread instead, which misses columns that +grow alike and turn near-parallel. Note that we pass `u` as an argument to reduce the number of copies that are necessary in the integration callback, since otherwise we'd need to copy from @@ -1005,10 +1003,6 @@ operate on `u` directly without extra copies. - u: Current solution vector array, updated in-place if fixfac is called - sing_flag: Indicates if normalization is occuring at a singular surface or not - resonant_rows: Resonant rows at an ideal crossing; when given, the reduction always runs and pivots on them - -### TODOs - -Add resizing logic for unorm arrays when ifix exceeds allocated size """ function compute_solution_norms!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal, sing_flag::Bool; resonant_rows=Int[]) @@ -1032,36 +1026,37 @@ function compute_solution_norms!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::F if !odet.new odet.unorm ./= odet.unorm0 end - uratio = maximum(odet.unorm) / minimum(odet.unorm) - # Perform Gaussian reduction if the ucrit ratio is reached or at a singular surface - if uratio > ctrl.ucrit || sing_flag - # TODO: add resizing logic here as well - if odet.ifix < ctrl.numunorms_init - odet.ifix += 1 - else - @warn "unorm storage reached, no longer saving fixfac data. Stability outputs and unorming will be correct, but cannot reconstruct `u`. \n - Increase `numunorms_init` if needed. Automatic resizing will be added in a future version." - end + # Reduce at a singular surface or once the U₁ columns lose independence + if sing_flag || column_condition(@view u[:, :, 1]) > ctrl.ucrit + odet.ifix == size(odet.fixfac, 3) && grow_reduction_storage!(odet) + odet.ifix += 1 apply_gaussian_reduction!(u, odet, intr, sing_flag; resonant_rows) odet.new = true end end +""" + column_condition(A::AbstractMatrix) -> Float64 + +1-norm condition-number estimate of `A` with its columns scaled to unit norm, from the `R` factor of its QR. +""" +column_condition(A::AbstractMatrix) = cond(UpperTriangular(qr!(A ./ norm.(eachcol(A))').R), 1) + """ apply_gaussian_reduction!(u::Array{ComplexF64,3}, odet::OdeState, intr::ForceFreeStatesInternal, sing_flag::Bool; resonant_rows=Int[]) Applies Gaussian reduction to orthogonalize solution vectors in `u`. Formerly `ode_fixup!`. Performs the same function as `ode_fixup` in the Fortran code, except now relevant `fixfac` data are stored in memory instead of dumped to `euler.bin`. -Used when the spread in norms exceeds a threshold or when a rational surface is reached. +Used when U₁ becomes ill-conditioned or when a rational surface is reached. This will update both `u` and relevant fields in `odet` in-place. See the description of `compute_solution_norms!` for more details on the benefits of in-place `u` updates. At an ideal crossing, the first pivots are the `resonant_rows`, each on the remaining column largest -there. Growth order cannot pick the resonant solution if a `ucrit` reduction fired just before the -crossing, since every column has then grown alike. +there. Growth order cannot pick the resonant solution if a reduction fired just before the crossing, +since every column has then grown alike. """ function apply_gaussian_reduction!(u::Array{ComplexF64,3}, odet::OdeState, intr::ForceFreeStatesInternal, sing_flag::Bool; resonant_rows=Int[]) diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index e3bb3fdd8..fc5e12596 100644 --- a/src/ForceFreeStates/Free.jl +++ b/src/ForceFreeStates/Free.jl @@ -130,9 +130,9 @@ calculations and data dumping. # Evaluate dV/dpsi at the plasma edge dV_dpsi = equil.profiles.dVdpsi_spline(psilim) - # Compute plasma response matrix W = U₂ * U₁⁻¹ + # Compute plasma response matrix W = U₂ * U₁⁻¹; factor U₁ itself, as its column norms can span many decades after a reduction wp = zeros(ComplexF64, numpert_total, numpert_total) - @views wp .= (odet.u[:, :, 2] / odet.u[:, :, 1]) ./ equil.psio^2 + @views wp .= (odet.u[:, :, 2] / lu(odet.u[:, :, 1])) ./ equil.psio^2 wv, vac = compute_scaled_wv(ctrl, equil, intr) @@ -275,8 +275,8 @@ wv matrix spline to `free_compute_wv_spline` and pass it in `odet.edge_scan.wvma dV_dpsi = equil.profiles.dVdpsi_spline(odet.psifac) - # Compute plasma response matrix - @views wp .= (odet.u[:, :, 2] / odet.u[:, :, 1]) ./ equil.psio^2 + # Compute plasma response matrix (lu(U₁): see free_run) + @views wp .= (odet.u[:, :, 2] / lu(odet.u[:, :, 1])) ./ equil.psio^2 # Retrieve raw vacuum matrix from spline and apply singfac analytically at the local q. # Singfac is not pre-applied in the spline (see free_compute_wv_spline) to avoid interpolating diff --git a/src/ForceFreeStates/Result.jl b/src/ForceFreeStates/Result.jl index da8fe3648..2c1fad194 100644 --- a/src/ForceFreeStates/Result.jl +++ b/src/ForceFreeStates/Result.jl @@ -240,7 +240,7 @@ function build_result( wp = if free_energies !== nothing free_energies.wp elseif odet !== nothing - (odet.u[:, :, 2] / odet.u[:, :, 1]) ./ equil.psio^2 + (odet.u[:, :, 2] / lu(odet.u[:, :, 1])) ./ equil.psio^2 else nothing end diff --git a/src/ForceFreeStates/Utils.jl b/src/ForceFreeStates/Utils.jl index 0d46b94bd..b2f19a002 100644 --- a/src/ForceFreeStates/Utils.jl +++ b/src/ForceFreeStates/Utils.jl @@ -41,6 +41,21 @@ function trim_storage!(odet::OdeState) odet.u_store = odet.u_store[:, :, :, 1:odet.step] end +""" + grow_reduction_storage!(odet::OdeState) + +Doubles the capacity of the per-reduction arrays (`fixfac`, `index`, `sing_flag`, `fixstep`, `zeroed_idx`). +""" +function grow_reduction_storage!(odet::OdeState) + nfix = size(odet.fixfac, 3) + odet.fixfac = cat(odet.fixfac, zeros(ComplexF64, size(odet.fixfac, 1), size(odet.fixfac, 2), nfix); dims=3) + odet.index = hcat(odet.index, zeros(Int, size(odet.index, 1), nfix)) + append!(odet.sing_flag, falses(nfix)) + append!(odet.fixstep, zeros(Int64, nfix)) + append!(odet.zeroed_idx, [Int[] for _ in 1:nfix]) + return odet +end + """ truncate_integration!(odet, intr, equil, step) diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index fa492d368..0f39279fb 100644 --- a/test/runtests_eulerlagrange.jl +++ b/test/runtests_eulerlagrange.jl @@ -330,7 +330,7 @@ end mpert = 2 odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(mpert, 10, 10, 10) intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal(; mpert=mpert) - ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; ucrit=10.0) + ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl() # Case 1: Basic norm computation odet.u = zeros(ComplexF64, 2, 2, 2) @@ -356,12 +356,18 @@ end GeneralizedPerturbedEquilibrium.ForceFreeStates.compute_solution_norms!(odet.u, odet, ctrl, intr, false) @test odet.unorm[1:intr.mpert] ≈ [1, 1] - # Case 4: Trigger fixup via ucrit + # Case 4: A large norm spread with independent columns does not reduce odet.unorm0 = ones(intr.mpert) - odet.u[:, 1, 1] .= [1000, 0] # large norm - odet.u[:, 2, 1] .= [1, 0] # small norm + odet.u[:, 1, 1] .= [1e8, 0] + odet.u[:, 2, 1] .= [0, 1] GeneralizedPerturbedEquilibrium.ForceFreeStates.compute_solution_norms!(odet.u, odet, ctrl, intr, false) - @test odet.new == true # implies fixup ran + @test odet.new == false && odet.ifix == 0 + + # Case 4b: Near-parallel columns reduce even though they grew alike + odet.u[:, 1, 1] .= [1, 0] + odet.u[:, 2, 1] .= [1, 1e-8] + GeneralizedPerturbedEquilibrium.ForceFreeStates.compute_solution_norms!(odet.u, odet, ctrl, intr, false) + @test odet.new == true && odet.ifix == 1 # Case 5: Trigger fixup via sing_flag odet.new = false @@ -371,6 +377,24 @@ end @test odet.new == true # fixup triggered end + @testset "column_condition and reduction storage growth" begin + FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates + @test FFS.column_condition(ComplexF64[1 0; 0 1e6]) ≈ 1 # column scaling does not count + @test FFS.column_condition(ComplexF64[1 1; 0 1e-8]) > 1e7 + intr = FFS.ForceFreeStatesInternal(; numpert_total=2) + ctrl = FFS.ForceFreeStatesControl() + odet = FFS.OdeState(2, 10, 1, 1) + for k in 1:3 + odet.u[:, :, 1] .= ComplexF64[1 0.5; 0 1] + odet.new = false + odet.unorm0 .= 1 + FFS.compute_solution_norms!(odet.u, odet, ctrl, intr, true) + end + @test odet.ifix == 3 && size(odet.fixfac, 3) == 4 + @test length(odet.sing_flag) == length(odet.fixstep) == length(odet.zeroed_idx) == size(odet.index, 2) == 4 + @test all(odet.sing_flag[1:3]) && all(odet.fixfac[i, i, k] == 1 for i in 1:2, k in 1:3) + end + @testset "ideal-crossing reduction isolates the resonant harmonic" begin FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates # Column 1 has the largest norm but no resonant (row 3) component; column 2 carries the resonance. @@ -381,8 +405,8 @@ end u0[:, 3, 1] .= [0.2, 2.0, 1.5] u0[:, :, 2] .= 1.0 intr = FFS.ForceFreeStatesInternal(; numpert_total=3) - ctrl = FFS.ForceFreeStatesControl(; ucrit=1e3) - fresh(u) = (odet = FFS.OdeState(3, 10, 10, 1); odet.u .= u; odet) + ctrl = FFS.ForceFreeStatesControl() + fresh(u) = (odet=FFS.OdeState(3, 10, 10, 1); odet.u.=u; odet) odet = fresh(u0) odet.new = false @@ -393,7 +417,7 @@ end @test odet.u[3, 1, 1] == 0 && odet.u[3, 3, 1] == 0 # only the lead column keeps the resonant harmonic @test odet.sing_flag[1] - # A crossing right after a ucrit reduction (new = true) still reduces on the resonant row. + # A crossing right after another reduction (new = true) still reduces on the resonant row. odet = fresh(u0) FFS.compute_solution_norms!(odet.u, odet, ctrl, intr, true; resonant_rows=[3]) @test odet.ifix == 1 && odet.index[1, 1] == 2 && odet.new diff --git a/test/runtests_fullruns.jl b/test/runtests_fullruns.jl index 305a23a64..d163e5368 100644 --- a/test/runtests_fullruns.jl +++ b/test/runtests_fullruns.jl @@ -10,7 +10,7 @@ using TOML true end - # A ucrit reduction landing just before an ideal crossing must not change which solution the crossing removes. + # A reduction landing just before an ideal crossing must not change which solution the crossing removes. @info "Running Solovev ideal example with frequent Gaussian reductions" @testset "et[1] is independent of the reduction frequency" begin deck = TOML.parsefile(joinpath(ex1, "gpec.toml")) @@ -24,7 +24,7 @@ using TOML end end et_ref = total_energy(1e3) - for ucrit in (2.0, 3.0, 5.0, 10.0) + for ucrit in (3.0, 10.0, 1e2, 1e6) @test total_energy(ucrit) ≈ et_ref rtol = 1e-6 end end diff --git a/test/test_data/regression_solovev_ideal_example/gpec.toml b/test/test_data/regression_solovev_ideal_example/gpec.toml index 9b3d346ac..0c9674e00 100644 --- a/test/test_data/regression_solovev_ideal_example/gpec.toml +++ b/test/test_data/regression_solovev_ideal_example/gpec.toml @@ -43,7 +43,7 @@ kinetic_source = "fixed" # Kinetic matrix source: "fixed" test matrices, or kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) diff --git a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml index 58c4a6cd1..1d6e0fecd 100644 --- a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml +++ b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml @@ -43,7 +43,7 @@ kinetic_source = "fixed" # Kinetic matrix source: "fixed" test matrices, or kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) diff --git a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml index 1a6bb6f13..229fc6539 100644 --- a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml @@ -45,7 +45,7 @@ kinetic_factor = 1.0 # Scaling of kinetic matrices (0 = ideal path; >0 integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization write_outputs_to_HDF5 = true # Write outputs to the HDF5 file verbose = false # Enable verbose logging # Solovev analytic equilibrium parameters (eq_type = "sol"); see SolovevConfig in src/Equilibrium. diff --git a/test/test_data/regression_solovev_kinetic_example/gpec.toml b/test/test_data/regression_solovev_kinetic_example/gpec.toml index 67692a315..0c0e9a2c8 100644 --- a/test/test_data/regression_solovev_kinetic_example/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_example/gpec.toml @@ -43,7 +43,7 @@ kinetic_source = "fixed" # Kinetic matrix source: "fixed" test matrices, or kinetic_factor = 1e-9 # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode; tiny value just exercises the path) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) diff --git a/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml b/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml index b2bdfb2b1..f301b03d0 100644 --- a/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml @@ -43,7 +43,7 @@ kinetic_source = "fixed" # Kinetic matrix source: "fixed" test matrices, or kinetic_factor = 1e-9 # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode; tiny value just exercises the path) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) diff --git a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml index 7debb9e8c..cdcbc3e14 100644 --- a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml @@ -45,7 +45,7 @@ kinetic_factor = 1.0 # Full-strength kinetic matrices (the "calculated" integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced -ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +ucrit = 1e3 # Forward: U₁ condition number that triggers Gaussian reduction; Riccati: entry size that triggers renormalization write_outputs_to_HDF5 = true # Write outputs to the HDF5 file verbose = false # Enable verbose logging From c1d081afb8abb39a6591e8530f0b05d37a547ba4 Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Mon, 12 Oct 2026 17:14:33 -0400 Subject: [PATCH 2/3] =?UTF-8?q?ForceFreeStates=20-=20REFACTOR=20-=20Inline?= =?UTF-8?q?=20the=20U=E2=82=81=20condition=20estimate=20at=20the=20reducti?= =?UTF-8?q?on=20check?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Drop the column_condition helper and its linear-algebra unit tests. Co-Authored-By: Claude Opus 5.5 --- src/ForceFreeStates/EulerLagrange.jl | 11 ++--------- test/runtests_eulerlagrange.jl | 4 +--- 2 files changed, 3 insertions(+), 12 deletions(-) diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 12ab6dfb6..0a722282f 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -1027,8 +1027,8 @@ function compute_solution_norms!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::F odet.unorm ./= odet.unorm0 end - # Reduce at a singular surface or once the U₁ columns lose independence - if sing_flag || column_condition(@view u[:, :, 1]) > ctrl.ucrit + # Reduce at a singular surface or once the U₁ columns lose independence (1-norm condition estimate of the column-normalized U₁) + @views if sing_flag || cond(UpperTriangular(qr!(u[:, :, 1] ./ norm.(eachcol(u[:, :, 1]))').R), 1) > ctrl.ucrit odet.ifix == size(odet.fixfac, 3) && grow_reduction_storage!(odet) odet.ifix += 1 apply_gaussian_reduction!(u, odet, intr, sing_flag; resonant_rows) @@ -1036,13 +1036,6 @@ function compute_solution_norms!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::F end end -""" - column_condition(A::AbstractMatrix) -> Float64 - -1-norm condition-number estimate of `A` with its columns scaled to unit norm, from the `R` factor of its QR. -""" -column_condition(A::AbstractMatrix) = cond(UpperTriangular(qr!(A ./ norm.(eachcol(A))').R), 1) - """ apply_gaussian_reduction!(u::Array{ComplexF64,3}, odet::OdeState, intr::ForceFreeStatesInternal, sing_flag::Bool; resonant_rows=Int[]) diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index 0f39279fb..89642dca1 100644 --- a/test/runtests_eulerlagrange.jl +++ b/test/runtests_eulerlagrange.jl @@ -377,10 +377,8 @@ end @test odet.new == true # fixup triggered end - @testset "column_condition and reduction storage growth" begin + @testset "reduction storage growth" begin FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates - @test FFS.column_condition(ComplexF64[1 0; 0 1e6]) ≈ 1 # column scaling does not count - @test FFS.column_condition(ComplexF64[1 1; 0 1e-8]) > 1e7 intr = FFS.ForceFreeStatesInternal(; numpert_total=2) ctrl = FFS.ForceFreeStatesControl() odet = FFS.OdeState(2, 10, 1, 1) From 746c6d2578f9ec0af2fd46a501331fb995ed815f Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Mon, 12 Oct 2026 17:14:33 -0400 Subject: [PATCH 3/3] =?UTF-8?q?ForceFreeStates=20-=20REFACTOR=20-=20Inline?= =?UTF-8?q?=20the=20U=E2=82=81=20condition=20estimate=20at=20the=20reducti?= =?UTF-8?q?on=20check?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Drop the column_condition helper and its linear-algebra unit tests. Co-Authored-By: Claude Opus 5.5 --- src/ForceFreeStates/EulerLagrange.jl | 11 ++--------- src/ForceFreeStates/Free.jl | 4 ++-- test/runtests_eulerlagrange.jl | 4 +--- 3 files changed, 5 insertions(+), 14 deletions(-) diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 12ab6dfb6..0a722282f 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -1027,8 +1027,8 @@ function compute_solution_norms!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::F odet.unorm ./= odet.unorm0 end - # Reduce at a singular surface or once the U₁ columns lose independence - if sing_flag || column_condition(@view u[:, :, 1]) > ctrl.ucrit + # Reduce at a singular surface or once the U₁ columns lose independence (1-norm condition estimate of the column-normalized U₁) + @views if sing_flag || cond(UpperTriangular(qr!(u[:, :, 1] ./ norm.(eachcol(u[:, :, 1]))').R), 1) > ctrl.ucrit odet.ifix == size(odet.fixfac, 3) && grow_reduction_storage!(odet) odet.ifix += 1 apply_gaussian_reduction!(u, odet, intr, sing_flag; resonant_rows) @@ -1036,13 +1036,6 @@ function compute_solution_norms!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::F end end -""" - column_condition(A::AbstractMatrix) -> Float64 - -1-norm condition-number estimate of `A` with its columns scaled to unit norm, from the `R` factor of its QR. -""" -column_condition(A::AbstractMatrix) = cond(UpperTriangular(qr!(A ./ norm.(eachcol(A))').R), 1) - """ apply_gaussian_reduction!(u::Array{ComplexF64,3}, odet::OdeState, intr::ForceFreeStatesInternal, sing_flag::Bool; resonant_rows=Int[]) diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index fc5e12596..28a6073f4 100644 --- a/src/ForceFreeStates/Free.jl +++ b/src/ForceFreeStates/Free.jl @@ -130,7 +130,7 @@ calculations and data dumping. # Evaluate dV/dpsi at the plasma edge dV_dpsi = equil.profiles.dVdpsi_spline(psilim) - # Compute plasma response matrix W = U₂ * U₁⁻¹; factor U₁ itself, as its column norms can span many decades after a reduction + # Compute plasma response matrix W = U₂ * U₁⁻¹ wp = zeros(ComplexF64, numpert_total, numpert_total) @views wp .= (odet.u[:, :, 2] / lu(odet.u[:, :, 1])) ./ equil.psio^2 @@ -275,7 +275,7 @@ wv matrix spline to `free_compute_wv_spline` and pass it in `odet.edge_scan.wvma dV_dpsi = equil.profiles.dVdpsi_spline(odet.psifac) - # Compute plasma response matrix (lu(U₁): see free_run) + # Compute plasma response matrix @views wp .= (odet.u[:, :, 2] / lu(odet.u[:, :, 1])) ./ equil.psio^2 # Retrieve raw vacuum matrix from spline and apply singfac analytically at the local q. diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index 0f39279fb..89642dca1 100644 --- a/test/runtests_eulerlagrange.jl +++ b/test/runtests_eulerlagrange.jl @@ -377,10 +377,8 @@ end @test odet.new == true # fixup triggered end - @testset "column_condition and reduction storage growth" begin + @testset "reduction storage growth" begin FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates - @test FFS.column_condition(ComplexF64[1 0; 0 1e6]) ≈ 1 # column scaling does not count - @test FFS.column_condition(ComplexF64[1 1; 0 1e-8]) > 1e7 intr = FFS.ForceFreeStatesInternal(; numpert_total=2) ctrl = FFS.ForceFreeStatesControl() odet = FFS.OdeState(2, 10, 1, 1)