From 4b8e4b21327dd0cd150eb8b7495bd69b33f45ec7 Mon Sep 17 00:00:00 2001 From: CharlesCNorton <135471798+CharlesCNorton@users.noreply.github.com> Date: Tue, 8 Sep 2026 15:14:40 -0400 Subject: [PATCH] Size the T_l recurrence passes of the analytic add-back by their growth rate --- .../singular_integrals/singular_integrals.cc | 168 ++++++------------ .../singular_integrals_test.cc | 63 ++++--- 2 files changed, 94 insertions(+), 137 deletions(-) diff --git a/src/vmecpp/cpp/vmecpp/free_boundary/singular_integrals/singular_integrals.cc b/src/vmecpp/cpp/vmecpp/free_boundary/singular_integrals/singular_integrals.cc index 1115e4e7d..d302949e0 100644 --- a/src/vmecpp/cpp/vmecpp/free_boundary/singular_integrals/singular_integrals.cc +++ b/src/vmecpp/cpp/vmecpp/free_boundary/singular_integrals/singular_integrals.cc @@ -8,10 +8,66 @@ #include #include #include +#include #include namespace vmecpp { +namespace { + +// Logarithm of the growth of the homogeneous solutions of the T_l recurrence +// over a forward pass, sqrt(B/A)^kL, above which the recurrence is run +// backward instead; a forward pass amplifies the rounding of its inputs by at +// most this factor, ten. +constexpr double kMaxForwardLogGrowth = std::numbers::ln10; + +// The zero seed of a backward pass contaminates T_l by (A/B)^{(top - l)/2} of +// T_top; the pass starts far enough above kL to bring that below 1e-17 at kL. +constexpr double kMinSeedLogDecay = 17.0 * std::numbers::ln10; + +// T_l = int_{-1}^{1} t^l / sqrt(A t^2 + 2 d t + B) dt for l = 0, ..., kL at +// tangential grid point kl from the three-term recurrence +// (l + 1) A T_{l+1} + (2 l + 1) d T_l + l B T_{l-1} = sqrtc2 - (-1)^l sqrta2 +// with T_0 given. The characteristic roots of the homogeneous recurrence are a +// complex pair of modulus sqrt(B/A), so a forward pass amplifies the rounding +// of T_0 and of the right-hand sides by sqrt(B/A)^l and is used while that +// stays below exp(kMaxForwardLogGrowth); otherwise the recurrence runs +// backward from a zero seed, whose contamination decays by sqrt(A/B) per step. +void ComputeTl(double A, double B, double d, double sqrtc2, double sqrta2, + double T0, int kL, int kl, std::vector& m_T) { + const auto rhs = [sqrtc2, sqrta2](int l) { + return sqrtc2 + (l % 2 == 0 ? -sqrta2 : sqrta2); + }; + m_T[0][kl] = T0; + const double log_ratio = std::log(B / A); + if (kL * log_ratio <= 2.0 * kMaxForwardLogGrowth) { + double T_prev = 0.0; // T_{-1} + for (int l = 0; l < kL; ++l) { + const double T_next = + (rhs(l) - (2 * l + 1) * d * m_T[l][kl] - l * B * T_prev) / + ((l + 1) * A); + T_prev = m_T[l][kl]; + m_T[l + 1][kl] = T_next; + } + return; + } + const int tail = + static_cast(std::ceil(2.0 * kMinSeedLogDecay / log_ratio)); + double T_hi = 0.0; // T_{l+1} + double T_cur = 0.0; // T_l + for (int l = kL + tail; l >= 1; --l) { + const double T_lo = + (rhs(l) - (2 * l + 1) * d * T_cur - (l + 1) * A * T_hi) / (l * B); + T_hi = T_cur; + T_cur = T_lo; + if (1 <= l - 1 && l - 1 <= kL) { + m_T[l - 1][kl] = T_lo; + } + } +} + +} // namespace + SingularIntegrals::SingularIntegrals(const Sizes* s, const FourierBasisFastToroidal* fb, const TangentialPartitioning* tp, @@ -220,117 +276,9 @@ void SingularIntegrals::prepareUpdate( (sqrtam * sqrta2[kl] - am[kl] + d[kl])) / sqrtam; - // Fill all Tlp[0..L] and Tlm[0..L] by picking the numerically stable - // direction of the three-term recurrence on a per-(+/-), per-kl basis. - // - // The characteristic roots of the homogeneous recurrence satisfy - // A*r^2 + 2*d*r + B = 0 -> |r1 r2| = B/A. - // For T^+: (A, B) = (ap, am), so |r1 r2| = am/ap. - // For T^-: (A, B) = (am, ap), so |r1 r2| = ap/am. - // If B > A (at least one |r| > 1), forward iteration is unstable and - // backward (Miller's algorithm) is used instead; otherwise forward is fine. - // - // T^{\pm}_0 is analytic (above); T^{\pm}_{-1} = 0. Forward produces - // T_{l+1} from T_l and T_{l-1}; backward produces T_{l-1} from T_l and - // T_{l+1} via the same recurrence solved in reverse. For backward, - // iteration starts from a zero seed far above the required L; the result - // is then normalized to match the analytic T^{\pm}_0. - // - // rhs(l+1) = sqrtc2 + (-1)^{l+1}*sqrta2 (same for T^+ and T^-). const int kL = mf + nf; - // The spurious solution is damped by (A/B)^kTailExtra per pass. - // For the worst realistic ratio (A/B ~ 0.5) suppression is ~0.5^50 ~ 1e-16. - const int kTailExtra = 50; - const int kLtail = kL + kTailExtra; - - // Only switch to backward when the forward spurious-mode growth - // (|r1 r2| = B/A) would actually exceed double precision over kL steps. - // Threshold: forward is considered stable as long as (B/A)^kL < 1e10, - // i.e. spurious amplitude stays within ~1e10 of the particular solution. - // Near-degenerate kl (|r1|~|r2|~1) fall in the forward branch, where - // zero-seed Miller is known to misconverge (spurious modes never damp). - // Formula: kL * ln(B/A) < ln(1e10) -> B/A < exp(ln(1e10)/kL). - constexpr double kLogGrowthThreshold = 10.0 * 2.30258509299; // ln(1e10) - const double logRatioP = - (am[kl] > ap[kl] && ap[kl] > 0.0) ? std::log(am[kl] / ap[kl]) : 0.0; - const bool useBackwardP = - static_cast(kL) * logRatioP > kLogGrowthThreshold; - const double logRatioM = - (ap[kl] > am[kl] && am[kl] > 0.0) ? std::log(ap[kl] / am[kl]) : 0.0; - const bool useBackwardM = - static_cast(kL) * logRatioM > kLogGrowthThreshold; - - // --- T^+: A = ap, B = am --- - Tlp[0][kl] = T0p; - if (useBackwardP) { - // forward unstable -> use backward recurrence. - double T_hi = 0.0; - double T_cur = 1.0e-300; - for (int l = kLtail; l >= 1; --l) { - const double rhs = sqrtc2[kl] + (l % 2 == 0 ? -1.0 : 1.0) * sqrta2[kl]; - const double T_lo = - (rhs - (2 * l + 1) * d[kl] * T_cur - (l + 1) * ap[kl] * T_hi) / - (l * am[kl]); - T_hi = T_cur; - T_cur = T_lo; - if (l - 1 <= kL) { - Tlp[l - 1][kl] = T_lo; - } - } - const double scaleP = T0p / Tlp[0][kl]; - for (int l = 0; l <= kL; ++l) { - Tlp[l][kl] *= scaleP; - } - } else { - // forward stable. - double T_prev = 0.0; // T^+_{-1} - int sgn = 1; - for (int fl = 0; fl < kL; ++fl) { - sgn = -sgn; - const double rhs = sqrtc2[kl] + sgn * sqrta2[kl]; - const double T_next = - (rhs - (2 * fl + 1) * d[kl] * Tlp[fl][kl] - fl * am[kl] * T_prev) / - (ap[kl] * (fl + 1)); - T_prev = Tlp[fl][kl]; - Tlp[fl + 1][kl] = T_next; - } - } - - // --- T^-: A = am, B = ap --- - Tlm[0][kl] = T0m; - if (useBackwardM) { - // forward unstable -> use backward recurrence. - double T_hi = 0.0; - double T_cur = 1.0e-300; - for (int l = kLtail; l >= 1; --l) { - const double rhs = sqrtc2[kl] + (l % 2 == 0 ? -1.0 : 1.0) * sqrta2[kl]; - const double T_lo = - (rhs - (2 * l + 1) * d[kl] * T_cur - (l + 1) * am[kl] * T_hi) / - (l * ap[kl]); - T_hi = T_cur; - T_cur = T_lo; - if (l - 1 <= kL) { - Tlm[l - 1][kl] = T_lo; - } - } - const double scaleM = T0m / Tlm[0][kl]; - for (int l = 0; l <= kL; ++l) { - Tlm[l][kl] *= scaleM; - } - } else { - // forward stable. - double T_prev = 0.0; // T^-_{-1} - int sgn = 1; - for (int fl = 0; fl < kL; ++fl) { - sgn = -sgn; - const double rhs = sqrtc2[kl] + sgn * sqrta2[kl]; - const double T_next = - (rhs - (2 * fl + 1) * d[kl] * Tlm[fl][kl] - fl * ap[kl] * T_prev) / - (am[kl] * (fl + 1)); - T_prev = Tlm[fl][kl]; - Tlm[fl + 1][kl] = T_next; - } - } + ComputeTl(ap[kl], am[kl], d[kl], sqrtc2[kl], sqrta2[kl], T0p, kL, kl, Tlp); + ComputeTl(am[kl], ap[kl], d[kl], sqrtc2[kl], sqrta2[kl], T0m, kL, kl, Tlm); } // kl } // prepareUpdate diff --git a/src/vmecpp/cpp/vmecpp/free_boundary/singular_integrals/singular_integrals_test.cc b/src/vmecpp/cpp/vmecpp/free_boundary/singular_integrals/singular_integrals_test.cc index 452c41081..4c4c756c8 100644 --- a/src/vmecpp/cpp/vmecpp/free_boundary/singular_integrals/singular_integrals_test.cc +++ b/src/vmecpp/cpp/vmecpp/free_boundary/singular_integrals/singular_integrals_test.cc @@ -164,24 +164,33 @@ TEST(TestSingularIntegrals, CheckConstants) { } // CheckConstants // Verify that prepareUpdate computes T^{+/-}_l accurately for ALL l in [0, L] -// at a range of Fourier resolutions. +// at a range of Fourier resolutions and metric coefficients. // -// Resolutions span the low-kL regime (forward recurrence is accurate, so the -// forward branch of prepareUpdate is exercised) all the way up to high-kL -// (Miller backward recurrence fires; forward would lose all digits). +// The cases span the low-kL regime, where the forward recurrence is accurate, +// up to kL = 45, where one direction of each recurrence loses all digits and +// the other has to run backward from a seed far above kL. // // Reference: 64-point Gauss-Legendre quadrature of the defining integral. // This is independent of both recurrence directions and reaches ~1e-13 // relative accuracy at the chosen coefficients up to l = 45. -class TlpTlmAccuracyTest - : public ::testing::TestWithParam> {}; +struct TlCase { + int mpol; + int ntor; + // metric coefficients (a, b2, c) = (guu, guv, gvv) of the tangent plane + double a; + double b2; + double c; +}; + +class TlpTlmAccuracyTest : public ::testing::TestWithParam {}; TEST_P(TlpTlmAccuracyTest, MatchesQuadrature) { // GL-64 reference noise at l near kL is ~1e-13; 1e-11 leaves a safe margin - // while staying far below the forward-recurrence error that Miller corrects. + // while staying far below the error of a recurrence run in the wrong + // direction or started too close above kL. static constexpr double kTolerance = 1.0e-11; - const auto [mpol, ntor] = GetParam(); + const auto [mpol, ntor, a_val, b2_val, c_val] = GetParam(); const bool lasym = false; const int nfp = 5; @@ -200,21 +209,14 @@ TEST_P(TlpTlmAccuracyTest, MatchesQuadrature) { const int kL = mf + nf; SingularIntegrals si(&s, &fb, &tp, &sg, nf, mf); - // Geometry coefficients chosen so am/ap ~ 4.7 and the discriminant - // d^2 - ap*am = 0.16 - 2.31 < 0 (so the integrand is smooth on [-1, 1]). - // ap = a + b2 + c = 0.7 - // am = a - b2 + c = 3.3 - // d = c - a = 0.4 - // At these coefficients the forward recurrence loses ~kL * log10(am/ap) - // ~ 0.67 * kL significant digits; for kL >= 15 it loses > 10 digits and the - // Miller activation threshold (growth > 1e10) triggers. - const double a_val = 0.8; - const double b2_val = -1.3; - const double c_val = 1.2; const double ap = a_val + b2_val + c_val; const double am = a_val - b2_val + c_val; const double d = c_val - a_val; - ASSERT_GT(am, ap) << "test setup: need am > ap to exercise Miller path"; + // The homogeneous solutions of the T^+ (T^-) recurrence grow by + // sqrt(am/ap) (sqrt(ap/am)) per forward step, so one of the two has to run + // backward once kL is large enough. + ASSERT_GT(std::max(am / ap, ap / am), 2.0) + << "test setup: need one unstable forward direction"; ASSERT_LT(d * d, ap * am) << "test setup: need smooth integrand on [-1,1]"; const int numLocal = tp.ztMax - tp.ztMin; @@ -253,15 +255,22 @@ TEST_P(TlpTlmAccuracyTest, MatchesQuadrature) { } } -// (mpol, ntor) pairs spanning the forward-stable regime (kL=15, Miller just -// barely fires) up to high-kL where forward would lose >30 digits. +// The first three cases share the coefficients ap = 0.7, am = 3.3, d = 0.4 +// (am/ap = 4.7) at kL = 15, 27 and 45; the last one takes the metric of the +// cth_like_free_bdy boundary at the point of its largest cross term +// (ap/am = 2.2) at kL = 33, where a backward pass seeded 50 steps above kL +// leaves 2e-9 of the seed in T^-_33. INSTANTIATE_TEST_SUITE_P( ResolutionSweep, TlpTlmAccuracyTest, - ::testing::Values(std::pair{6, 8}, std::pair{12, 14}, - std::pair{20, 24}), - [](const ::testing::TestParamInfo>& info) { - return "mpol" + std::to_string(info.param.first) + "_ntor" + - std::to_string(info.param.second); + ::testing::Values(TlCase{6, 8, 0.8, -1.3, 1.2}, + TlCase{12, 14, 0.8, -1.3, 1.2}, + TlCase{20, 24, 0.8, -1.3, 1.2}, + TlCase{16, 16, 1.360e-02, 2.241e-02, 4.626e-02}), + [](const ::testing::TestParamInfo& info) { + return "mpol" + std::to_string(info.param.mpol) + "_ntor" + + std::to_string(info.param.ntor) + "_ap" + + std::to_string(static_cast( + 1000.0 * (info.param.a + info.param.b2 + info.param.c))); }); } // namespace vmecpp