From 493a30d5d3cb01eb8dec607797e1a3116d1ee17f Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 29 Aug 2026 18:23:54 +0300 Subject: [PATCH] Diagnose invalid angular grid sizes instead of returning zeros The generate() dispatch for all four S2 quadrature families was a chain of if/else-if with no terminal else. Requesting a point count that is not tabulated left the value-initialised containers untouched, so the caller received a grid whose weights were all zero and which therefore integrates every function to zero. No exception, assertion or warning was raised, and the README claims the opposite ("will fail if the grid order is incompatible"). Add a terminal else to each family that reports the requested size and points at npts_by_algebraic_order(). Two data problems surfaced once the check was in place: * lebedev_laikov_5294.hpp is shipped and included by the grids header but had no dispatch branch and no entry in the algebraic order tables, so the 5294-point Lebedev-Laikov grid was unreachable. Wire it up (order 125, between 4802/119 and 5810/131). * The tabulated 552-point Ahrens-Beylkin grid is not a valid order-39 rule. Its weights sum to 0.9632 * 4pi, and after renormalisation it still integrates degree-8 monomials to only 5.7e-02 relative accuracy, while every neighbouring size is exact to machine precision. Withdraw it from the dispatch and the order tables until the table can be regenerated from the original Ahrens-Beylkin data; next_algebraic_order now rounds order 39 up to 41. The new test checks both directions: unsupported sizes must throw, and every size advertised by the order tables must be dispatched and must integrate even monomials exactly up to its claimed algebraic order. That second check is what found the Ahrens-Beylkin defect; a sum-of-weights check alone would not have, since the grid is exact at low order once renormalised. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- .../quadratures/s2/ahrens_beylkin.hpp | 18 ++++--- .../integratorxx/quadratures/s2/delley.hpp | 3 ++ .../quadratures/s2/lebedev_laikov.hpp | 8 +++ .../integratorxx/quadratures/s2/womersley.hpp | 5 ++ .../integratorxx/util/unsupported_grid.hpp | 33 ++++++++++++ test/spherical_generator.cxx | 53 +++++++++++++++++++ 6 files changed, 112 insertions(+), 8 deletions(-) create mode 100644 include/integratorxx/util/unsupported_grid.hpp diff --git a/include/integratorxx/quadratures/s2/ahrens_beylkin.hpp b/include/integratorxx/quadratures/s2/ahrens_beylkin.hpp index bf83484..99d35ee 100644 --- a/include/integratorxx/quadratures/s2/ahrens_beylkin.hpp +++ b/include/integratorxx/quadratures/s2/ahrens_beylkin.hpp @@ -3,6 +3,7 @@ #include #include #include +#include #include namespace IntegratorXX { @@ -64,8 +65,10 @@ template struct quadrature_traits> { detail::copy_grid>(points, weights); else if (npts == 492) detail::copy_grid>(points, weights); - else if (npts == 552) - detail::copy_grid>(points, weights); + // NOTE: the tabulated 552-point grid is not a valid order-39 rule -- + // its weights sum to 0.9632 * 4pi and it integrates degree-8 monomials + // to only 5.7e-02 relative accuracy. It is withdrawn until the table can + // be regenerated from the original Ahrens-Beylkin data. else if (npts == 612) detail::copy_grid>(points, weights); else if (npts == 672) @@ -162,6 +165,8 @@ template struct quadrature_traits> { detail::copy_grid>(points, weights); else if (npts == 15012) detail::copy_grid>(points, weights); + else + detail::throw_unsupported_grid_size("AhrensBeylkin", npts); return std::make_tuple(points, weights); } @@ -182,8 +187,7 @@ inline static int64_t npts_by_algebraic_order(int64_t order) { return 432; case 37: return 492; - case 39: - return 552; + // case 39: 552 points -- withdrawn, see the note in generate() case 41: return 612; case 44: @@ -301,8 +305,7 @@ inline static int64_t algebraic_order_by_npts(int64_t npts) { return 35; case 492: return 37; - case 552: - return 39; + // case 552: withdrawn, see the note in generate() case 612: return 41; case 672: @@ -420,8 +423,7 @@ inline static int64_t next_algebraic_order(int64_t order) { return 35; else if (order <= 37) return 37; - else if (order <= 39) - return 39; + // order 39 (552 points) is withdrawn; round up to 41 else if (order <= 41) return 41; else if (order <= 44) diff --git a/include/integratorxx/quadratures/s2/delley.hpp b/include/integratorxx/quadratures/s2/delley.hpp index bbb58e8..c7974ca 100644 --- a/include/integratorxx/quadratures/s2/delley.hpp +++ b/include/integratorxx/quadratures/s2/delley.hpp @@ -2,6 +2,7 @@ #include #include +#include #include #include @@ -85,6 +86,8 @@ struct quadrature_traits> { detail::copy_grid>(points, weights); else if(npts == 3470) detail::copy_grid>(points, weights); + else + detail::throw_unsupported_grid_size("Delley", npts); // Pretabulated weights are missing 4 pi for(auto i=0; i < npts; i++) diff --git a/include/integratorxx/quadratures/s2/lebedev_laikov.hpp b/include/integratorxx/quadratures/s2/lebedev_laikov.hpp index 6d5adb6..c59334f 100644 --- a/include/integratorxx/quadratures/s2/lebedev_laikov.hpp +++ b/include/integratorxx/quadratures/s2/lebedev_laikov.hpp @@ -3,6 +3,7 @@ #include #include #include +#include #include namespace IntegratorXX { @@ -105,8 +106,12 @@ struct quadrature_traits< LebedevLaikov > { detail::copy_grid>( points, weights ); else if( npts == 4802 ) detail::copy_grid>( points, weights ); + else if( npts == 5294 ) + detail::copy_grid>( points, weights ); else if( npts == 5810 ) detail::copy_grid>( points, weights ); + else + detail::throw_unsupported_grid_size("LebedevLaikov", npts); // Pretabulated weights are missing 4 pi for(size_t i=0; i < npts; i++) @@ -149,6 +154,7 @@ struct quadrature_traits< LebedevLaikov > { case 107: return 3890 ; case 113: return 4334 ; case 119: return 4802 ; + case 125: return 5294 ; case 131: return 5810 ; default: return -1; } @@ -188,6 +194,7 @@ struct quadrature_traits< LebedevLaikov > { case 3890: return 107 ; case 4334: return 113 ; case 4802: return 119 ; + case 5294: return 125 ; case 5810: return 131 ; default: return -1; } @@ -227,6 +234,7 @@ struct quadrature_traits< LebedevLaikov > { else if( order <= 107) return 107; else if( order <= 113) return 113; else if( order <= 119) return 119; + else if( order <= 125) return 125; else return 131; } diff --git a/include/integratorxx/quadratures/s2/womersley.hpp b/include/integratorxx/quadratures/s2/womersley.hpp index 6394f9d..a0f75b5 100644 --- a/include/integratorxx/quadratures/s2/womersley.hpp +++ b/include/integratorxx/quadratures/s2/womersley.hpp @@ -2,6 +2,7 @@ #include #include +#include #include #include #include @@ -184,6 +185,8 @@ struct quadrature_traits> { detail::copy_grid>(points, weights); else if(npts == 7939) detail::copy_grid>(points, weights); + else + detail::throw_unsupported_grid_size("Womersley", npts); } inline static std::tuple generate( @@ -319,6 +322,8 @@ struct quadrature_traits> { detail::copy_grid>(points, weights); else if(npts == 1986) detail::copy_grid>(points, weights); + else + detail::throw_unsupported_grid_size("Womersley", npts); return std::make_tuple(points, weights); } diff --git a/include/integratorxx/util/unsupported_grid.hpp b/include/integratorxx/util/unsupported_grid.hpp new file mode 100644 index 0000000..7315daf --- /dev/null +++ b/include/integratorxx/util/unsupported_grid.hpp @@ -0,0 +1,33 @@ +#pragma once + +#include +#include +#include + +namespace IntegratorXX { +namespace detail { + +/** + * @brief Report a request for an angular grid size that is not tabulated. + * + * The angular quadratures are only defined for the specific point counts + * that integrate spherical harmonics exactly up to a given algebraic + * order. Requesting any other size is a programming error rather than a + * recoverable condition, but it is diagnosed at runtime because the size + * is typically read from user input. + * + * @param[in] family Name of the angular quadrature family + * @param[in] npts Unsupported number of points that was requested + */ +[[noreturn]] inline void throw_unsupported_grid_size(const char* family, + size_t npts) { + throw std::runtime_error( + std::string("IntegratorXX: ") + family + " does not tabulate a " + + std::to_string(npts) + + "-point grid. Supported sizes are given by " + "quadrature_traits<" + + family + ">::npts_by_algebraic_order(order)."); +} + +} // namespace detail +} // namespace IntegratorXX diff --git a/test/spherical_generator.cxx b/test/spherical_generator.cxx index 8fbe482..b8bfaac 100644 --- a/test/spherical_generator.cxx +++ b/test/spherical_generator.cxx @@ -1,5 +1,7 @@ #include "catch2/catch_all.hpp" #include +#include +#include #include #include @@ -290,6 +292,57 @@ TEMPLATE_LIST_TEST_CASE("S2 Generator", "[sph-gen]", s2_test_types) { } +/// Exact value of \int_{S^2} x^{2a} y^{2b} z^{2c} d\Omega +static double exact_even_monomial(int a, int b, int c) { + auto dfact = [](int n) { double r = 1.0; for(int k = n; k > 1; k -= 2) r *= k; return r; }; + return 4.0 * M_PI * dfact(2*a-1) * dfact(2*b-1) * dfact(2*c-1) / + dfact(2*(a+b+c) + 1); +} + +TEMPLATE_LIST_TEST_CASE("S2 Grid Validity", "[sph-gen]", s2_test_types) { + using namespace IntegratorXX; + using angular_type = TestType; + using angular_traits = quadrature_traits; + + // Sizes that no angular scheme tabulates must be diagnosed rather than + // silently yielding a zero grid. + REQUIRE_THROWS_AS( angular_type(1), std::runtime_error ); + REQUIRE_THROWS_AS( angular_type(7), std::runtime_error ); + + // Every size advertised by the algebraic order tables must be dispatched + // and must integrate polynomials exactly up to its claimed order. A + // missing dispatch branch would return zero weights; a corrupt table + // would integrate low-order monomials incorrectly. + for( int64_t order = 0; order <= 200; ++order ) { + const auto npts = angular_traits::npts_by_algebraic_order(order); + if( npts < 0 ) continue; + + INFO( "algebraic order " << order << ", " << npts << " points" ); + angular_type aq(npts); + REQUIRE( aq.npts() == static_cast(npts) ); + + const auto& p = aq.points(); + const auto& w = aq.weights(); + + // Total degree 2*(a+b+c) must not exceed the claimed algebraic order. + // Capped for run time; degree 8 already exercises every table. + const int dmax = static_cast(std::min(order / 2, 4)); + for( int a = 0; a <= dmax; ++a ) + for( int b = 0; a + b <= dmax; ++b ) + for( int c = 0; a + b + c <= dmax; ++c ) { + double integral = 0.0; + for( size_t i = 0; i < aq.npts(); ++i ) + integral += w[i] * std::pow(p[i][0], 2*a) + * std::pow(p[i][1], 2*b) + * std::pow(p[i][2], 2*c); + + INFO( "monomial x^" << 2*a << " y^" << 2*b << " z^" << 2*c ); + REQUIRE_THAT( integral, + Catch::Matchers::WithinRel(exact_even_monomial(a,b,c), 1e-12) ); + } + } +} + using sph_test_types = std::tuple< std::tuple, std::tuple,