From 583e5ca9af0c92a73fc933efc40c51a0c598f938 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. Drop it from the dispatch and the order tables; next_algebraic_order now rounds order 39 up to 41. The following commit restores the size once the cause is known. 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 Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- .../quadratures/s2/ahrens_beylkin.hpp | 11 ++-- .../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, 105 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 248710d..58bc366 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 #include @@ -65,8 +66,6 @@ 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); else if (npts == 612) detail::copy_grid>(points, weights); else if (npts == 672) @@ -163,6 +162,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); } @@ -183,8 +184,6 @@ inline static int64_t npts_by_algebraic_order(int64_t order) { return 432; case 37: return 492; - case 39: - return 552; case 41: return 612; case 44: @@ -302,8 +301,6 @@ inline static int64_t algebraic_order_by_npts(int64_t npts) { return 35; case 492: return 37; - case 552: - return 39; case 612: return 41; case 672: @@ -421,8 +418,6 @@ inline static int64_t next_algebraic_order(int64_t order) { return 35; else if (order <= 37) return 37; - else if (order <= 39) - return 39; 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 f440904..80aa430 100644 --- a/include/integratorxx/quadratures/s2/delley.hpp +++ b/include/integratorxx/quadratures/s2/delley.hpp @@ -2,6 +2,7 @@ #include #include +#include #include #include #include @@ -86,6 +87,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 290917e..6d118db 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 #include @@ -106,8 +107,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++) @@ -150,6 +155,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; } @@ -189,6 +195,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; } @@ -228,6 +235,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 b0a6315..5aae6da 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 @@ -185,6 +186,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( @@ -320,6 +323,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,