diff --git a/include/integratorxx/quadratures/gausschebyshev1.hpp b/include/integratorxx/quadratures/gausschebyshev1.hpp index 41377f3..d489976 100644 --- a/include/integratorxx/quadratures/gausschebyshev1.hpp +++ b/include/integratorxx/quadratures/gausschebyshev1.hpp @@ -39,8 +39,7 @@ class GaussChebyshev1 using point_container = typename base_type::point_container; using weight_container = typename base_type::weight_container; - GaussChebyshev1(size_t npts, point_type lo, point_type up) - : base_type(npts, lo, up) {} + GaussChebyshev1(size_t npts) : base_type(npts) {} GaussChebyshev1(const GaussChebyshev1 &) = default; GaussChebyshev1(GaussChebyshev1 &&) noexcept = default; @@ -54,8 +53,7 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; - inline static std::tuple generate( - size_t npts, point_type lo, point_type up) { + inline static std::tuple generate(size_t npts) { point_container points(npts); weight_container weights(npts); diff --git a/include/integratorxx/quadratures/gausschebyshev2.hpp b/include/integratorxx/quadratures/gausschebyshev2.hpp index ca06e7c..9b3271f 100644 --- a/include/integratorxx/quadratures/gausschebyshev2.hpp +++ b/include/integratorxx/quadratures/gausschebyshev2.hpp @@ -42,8 +42,7 @@ class GaussChebyshev2 using point_container = typename base_type::point_container; using weight_container = typename base_type::weight_container; - GaussChebyshev2(size_t npts, point_type lo, point_type up) - : base_type(npts, lo, up) {} + GaussChebyshev2(size_t npts) : base_type(npts) {} GaussChebyshev2(const GaussChebyshev2 &) = default; GaussChebyshev2(GaussChebyshev2 &&) noexcept = default; @@ -56,9 +55,10 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; + + inline static constexpr bool bound_inclusive = false; - inline static std::tuple generate( - size_t npts, point_type lo, point_type up) { + inline static std::tuple generate(size_t npts) { const weight_type pi_ov_npts_p_1 = M_PI / (npts + 1); weight_container weights(npts); @@ -96,6 +96,7 @@ struct quadrature_traits> { return std::make_tuple(points, weights); } + }; } // namespace IntegratorXX diff --git a/include/integratorxx/quadratures/gausschebyshev2modified.hpp b/include/integratorxx/quadratures/gausschebyshev2modified.hpp index 1d8a91d..5213597 100644 --- a/include/integratorxx/quadratures/gausschebyshev2modified.hpp +++ b/include/integratorxx/quadratures/gausschebyshev2modified.hpp @@ -43,8 +43,7 @@ class GaussChebyshev2Modified using point_container = typename base_type::point_container; using weight_container = typename base_type::weight_container; - GaussChebyshev2Modified(size_t npts, point_type lo, point_type up) - : base_type(npts, lo, up) {} + GaussChebyshev2Modified(size_t npts) : base_type(npts) {} GaussChebyshev2Modified(const GaussChebyshev2Modified &) = default; GaussChebyshev2Modified(GaussChebyshev2Modified &&) noexcept = default; @@ -58,8 +57,7 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; - inline static std::tuple generate( - size_t npts, point_type lo, point_type up) { + inline static std::tuple generate(size_t npts) { const weight_type oonpp = 1.0 / (npts + 1); point_container points(npts); diff --git a/include/integratorxx/quadratures/gausschebyshev3.hpp b/include/integratorxx/quadratures/gausschebyshev3.hpp index 9b8f799..2e4684d 100644 --- a/include/integratorxx/quadratures/gausschebyshev3.hpp +++ b/include/integratorxx/quadratures/gausschebyshev3.hpp @@ -39,8 +39,7 @@ class GaussChebyshev3 using point_container = typename base_type::point_container; using weight_container = typename base_type::weight_container; - GaussChebyshev3(size_t npts, point_type lo, point_type up) - : base_type(npts, lo, up) {} + GaussChebyshev3(size_t npts) : base_type(npts) {} GaussChebyshev3(const GaussChebyshev3 &) = default; GaussChebyshev3(GaussChebyshev3 &&) noexcept = default; @@ -54,8 +53,7 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; - inline static std::tuple generate( - size_t npts, point_type lo, point_type up) { + inline static std::tuple generate(size_t npts) { const weight_type pi_ov_2n_p_1 = M_PI / (2 * npts + 1); weight_container weights(npts); @@ -77,6 +75,7 @@ struct quadrature_traits> { wi *= std::sqrt((1.0 - xi) / xi); // Finally, convert from [0,1] to [-1,1] + // TODO: This should be done externally points[idx] = 2 * xi - 1.0; weights[idx] = 2.0 * wi; } diff --git a/include/integratorxx/quadratures/mhl.hpp b/include/integratorxx/quadratures/mhl.hpp index cec02dc..1359eb5 100644 --- a/include/integratorxx/quadratures/mhl.hpp +++ b/include/integratorxx/quadratures/mhl.hpp @@ -2,142 +2,78 @@ #include #include +#include namespace IntegratorXX { /** - * @brief Implementation of the Murray-Handy-Laming radial quadrature. - * - * Generates a quadrature on the bounds (0, inf). Suitable for integrands - * which tend to zero as their argument tends to 0 and inf. Tailored for - * radial integrands, i.e. r^2 * f(r), with lim_{r->inf} f(r) = 0. + * @brief Implementation of the Murray-Handy-Laming radial quadrature + * transformation rules. * * Reference: * Molecular Physics, 78:4, 997-1014, * DOI: https://doi.org/10.1080/00268979300100651 * - * @tparam PointType Type describing the quadrature points - * @tparam WeightType Type describing the quadrature weights + * @tparam M Integer to modulate the MHL transformation. + * Typically taken to be 2. */ -template -class MurrayHandyLaming : - public Quadrature> { +template +class MurrayHandyLamingRadialTraits { - using base_type = Quadrature>; + double R_; ///< Radial scaling factor public: - using point_type = typename base_type::point_type; - using weight_type = typename base_type::weight_type; - using point_container = typename base_type::point_container; - using weight_container = typename base_type::weight_container; - + MurrayHandyLamingRadialTraits(double R = 1.0) : R_(R) {} + /** - * @brief Construct the Murray-Handy-Laming radial quadrature + * @brief Transformation rule for the MHL radial quadrature + * + * @param[in] x Point in (0,1) + * @return r = R * (x / (1-x))^M + */ + template + inline auto radial_transform(PointType x) const noexcept { + return R_ * std::pow( x / (1.0 - x), M ); + } + + /** + * @brief Jacobian of the MHL radial transformation * - * @param[in] npts Number of quadrature points to generate - * @param[in] R Radial scaling factor. Suggested to be set to - * correspond to the Bragg-Slater radius of the - * appropriate nucleus - * @param[in] m Exponent factor for the quadrature. 2 was the - * default suggested in the reference. + * @param[in] x Point in (0,1) + * @returns dr/dx (see `radial_transform`) */ - MurrayHandyLaming(size_t npts, weight_type R = 1., int m = 2): - base_type( npts, R, m ) { } + template + inline auto radial_jacobian(PointType x) const noexcept { + return R_ * M * std::pow(x, M-1) / std::pow(1.0 - x, M+1); + } - MurrayHandyLaming( const MurrayHandyLaming& ) = default; - MurrayHandyLaming( MurrayHandyLaming&& ) noexcept = default; }; - - - - /** - * @brief Quadrature traits for the Murray-Handy-Laming quadrature + * @brief Implementation of the Murray-Handy-Laming radial quadrature. + * + * Taken as the convolution of the Uniform (Trapezoid) quadrature with + * the MHL radial transformation. See MurrayHandyLamingRadialTraits for + * details. + * + * Suitable for integrands which tend to zero as their argument tends + * to 0 and inf. Tailored for radial integrands, i.e. r^2 * f(r), with + * lim_{r->inf} f(r) = 0. + * + * Reference: + * Molecular Physics, 78:4, 997-1014, + * DOI: https://doi.org/10.1080/00268979300100651 * * @tparam PointType Type describing the quadrature points * @tparam WeightType Type describing the quadrature weights */ template -struct quadrature_traits< - MurrayHandyLaming, - std::enable_if_t< - std::is_floating_point_v and - std::is_floating_point_v - > -> { - - using point_type = PointType; - using weight_type = WeightType; - - using point_container = std::vector< point_type >; - using weight_container = std::vector< weight_type >; - - - /** - * @brief Generator for the Murray-Handy-Laming quadrature - * - * @param[in] npts Number of quadrature points to generate - * @param[in] R Radial scaling factor. Suggested to be set to - * correspond to the Bragg-Slater radius of the - * appropriate nucleus - * @param[in] m Exponent fector for the quadrature. 2 was the - * default suggested in the original reference. - * - * @returns Tuple of quadrature points and weights - */ - inline static std::tuple - generate( size_t npts, weight_type R, int m ) { - - point_container points( npts ); - weight_container weights( npts ); - - - using base_quad_traits = - quadrature_traits>; - - /* - * Generate uniform trapezoid points on [0,1] - * ux(j) = j/(m-1) - * uw(j) = 1/(m-1) - * m = npts + 2, j in [0, npts+2) - */ - auto [ux, uw] = base_quad_traits::generate( npts+2, 0., 1. ); - - /* - * Perform Murray, Handy and Laming transformation - * - * Original reference: - * Molecular Physics, 78:4, 997-1014, - * DOI: https://doi.org/10.1080/00268979300100651 - * - * Closed form for points / weights obtained from: - * Journal of Computational Chemistry, 24: 732–740, 2003 - * DOI: https://doi.org/10.1002/jcc.10211 - * - * x(i) = ux(i+1) - * r(i) = [x(i) / (1 - x(i))]^m - * w(i) = uw(i+1) * m * x(i)^(m-1) * [1 - x(i)]^(-m-1) - * i in [0, npts) - * - * XXX: i+1 offset on trapezoid points ignores enpoints of trapezoid - * quadrature (f(r) = 0 with r in {0,inf}) - */ - for( size_t i = 0; i < npts; ++i ) { - const auto xi = ux[i+1]; - const auto one_m_xi = 1. - xi; - points[i] = R * std::pow( xi / one_m_xi, m ); - weights[i] = R * uw[i+1] * m * std::pow( xi, m-1 ) / std::pow( one_m_xi, m+1); - } - - - - return std::make_tuple( points, weights ); - - } - -}; +using MurrayHandyLaming = RadialTransformQuadrature< + UniformTrapezoid, + MurrayHandyLamingRadialTraits<2> +>; } + diff --git a/include/integratorxx/quadratures/muraknowles.hpp b/include/integratorxx/quadratures/muraknowles.hpp index 5847ce0..534070a 100644 --- a/include/integratorxx/quadratures/muraknowles.hpp +++ b/include/integratorxx/quadratures/muraknowles.hpp @@ -2,9 +2,11 @@ #include #include +#include namespace IntegratorXX { +#if 0 /** * @brief Implementation of the Mura-Knowles radial quadrature. * @@ -132,4 +134,35 @@ struct quadrature_traits< }; +#else + +class MuraKnowlesRadialTraits { + + double R_; ///< Radial scaling factor + +public: + + MuraKnowlesRadialTraits(double R = 1.0) : R_(R) { } + + template + inline auto radial_transform(PointType x) const noexcept { + return -R_ * std::log(1.0 - x*x*x); + } + + template + inline auto radial_jacobian(PointType x) const noexcept { + const auto x2 = x*x; + return R_ * 3.0 * x2 / (1.0 - x2 * x); + } + +}; + +template +using MuraKnowles = RadialTransformQuadrature< + UniformTrapezoid, + MuraKnowlesRadialTraits +>; + +#endif + } diff --git a/include/integratorxx/quadratures/radial_transform.hpp b/include/integratorxx/quadratures/radial_transform.hpp new file mode 100644 index 0000000..93f8096 --- /dev/null +++ b/include/integratorxx/quadratures/radial_transform.hpp @@ -0,0 +1,63 @@ +#pragma once + +#include + +namespace IntegratorXX { + +template +struct RadialTransformQuadrature : + public Quadrature> { + + using base_type = Quadrature>; + +public: + + using point_type = typename base_type::point_type; + using weight_type = typename base_type::weight_type; + using point_container = typename base_type::point_container; + using weight_container = typename base_type::weight_container; + + RadialTransformQuadrature(size_t npts, const RadialTraits& traits = RadialTraits()) : + base_type( npts, traits ) { } + + RadialTransformQuadrature( const RadialTransformQuadrature& ) = default; + RadialTransformQuadrature( RadialTransformQuadrature&& ) noexcept = default; + +}; + +template +struct quadrature_traits> { + + using point_type = typename BaseQuad::point_type; + using weight_type = typename BaseQuad::weight_type; + using point_container = typename BaseQuad::point_container; + using weight_container = typename BaseQuad::weight_container; + + inline static std::tuple + generate( size_t npts, const RadialTraits& traits ) { + + using base_quad_traits = quadrature_traits; + + point_container points( npts ); + weight_container weights( npts ); + + + const auto npts_base = base_quad_traits::bound_inclusive ? npts+2 : npts; + auto [base_x, base_w] = base_quad_traits::generate(npts_base); + + const auto ipts_offset = !!base_quad_traits::bound_inclusive; + for(size_t i = 0; i < npts; ++i) { + const auto xi = base_x[i + ipts_offset]; + const auto wi = base_w[i + ipts_offset]; + points[i] = traits.radial_transform(xi); + weights[i] = wi * traits.radial_jacobian(xi); + } + + return std::make_tuple(points, weights); + } + +}; + + + +} diff --git a/include/integratorxx/quadratures/treutlerahlrichs.hpp b/include/integratorxx/quadratures/treutlerahlrichs.hpp index 48ceb91..5759d36 100644 --- a/include/integratorxx/quadratures/treutlerahlrichs.hpp +++ b/include/integratorxx/quadratures/treutlerahlrichs.hpp @@ -1,126 +1,95 @@ #pragma once #include +#include +#include namespace IntegratorXX { - /** - * @brief Implementation of the Treutler-Ahlrichs radial quadrature. - * - * Generates a quadrature on the bounds (0, inf). Suitable for integrands - * which tend to zero as their argument tends to 0 and inf. Tailored for - * radial integrands, i.e. r^2 * f(r), with lim_{r->inf} f(r) = 0. + * @brief Implementation of the Treutler-Ahlrichs M3+M4 radial + * quadrature transformation rules. * * Reference: * J. Chem. Phys. 102, 346 (1995) * DOI: https://doi.org/10.1063/1.469408 - * - * @tparam PointType Type describing the quadrature points - * @tparam WeightType Type describing the quadrature weights */ -template -class TreutlerAhlrichs : - public Quadrature> { +class TreutlerAhlrichsRadialTraits { - using base_type = Quadrature>; + double R_; + double alpha_; public: - using point_type = typename base_type::point_type; - using weight_type = typename base_type::weight_type; - using point_container = typename base_type::point_container; - using weight_container = typename base_type::weight_container; + inline static constexpr double a = 1.0; + inline static constexpr double ln_2 = 0.693147180559945309417232; /** - * @brief Construct the Treutler-Ahlrichs radial quadrature + * Specify Treutler-Ahlrichs quadrature parameters + * + * M3: Equation (18) of J. Chem. Phys. 102, 346 (1995) + * M4: Equation (19) of J. Chem. Phys. 102, 346 (1995) + * + * Default to M4 (alpha = 0.6). M3 resolved with alpha = 0.0 * - * @param[in] npts Number of quadrature points to generate - * @param[in] R Radial scaling factor. Table for suggested - * values is given in the original reference. - * @param[in] alpha Exponent factor for the quadrature. 0.6 was the - * default suggested in the reference. + * @param[in] R Radial scaling factor + * @param[in] alpha TA exponential factor */ - TreutlerAhlrichs(size_t npts, weight_type R = 1., weight_type alpha = 0.6): - base_type( npts, R, alpha ) { } - - TreutlerAhlrichs( const TreutlerAhlrichs& ) = default; - TreutlerAhlrichs( TreutlerAhlrichs&& ) noexcept = default; -}; + TreutlerAhlrichsRadialTraits(double R = 1.0, double alpha = 0.6) : + R_(R), alpha_(alpha) { } + /** + * @brief Transformation rule for the TA M3+M4 radial quadratures + * + * @param[in] x Point in (-1,1) + * @return r = (a+x)^alpha * log((a+1)/(1-x)) / ln(2) + */ + template + inline auto radial_transform(PointType x) const noexcept { + const auto pow_term = std::pow(a + x, alpha_); + const auto log_term = std::log((a + 1.0) / (1.0 - x)); + return R_ * pow_term * log_term / ln_2; + }; + /** + * @brief Jacobian of the TA M3+M4 radial transformations + * + * @param[in] x Point in (-1,1) + * @returns dr/dx (see `radial_transform`) + */ + template + inline auto radial_jacobian(PointType x) const noexcept { + const auto pow_term = std::pow(a + x, alpha_); + const auto log_term = std::log((a + 1.0) / (1.0 - x)); + return R_ * pow_term / ln_2 * ( alpha_ * log_term / (a+x) + (1./(1.-x)) ); + } +}; /** - * @brief Quadrature traits for the Treutler-Ahlrichs quadrature + * @brief Implementation of the Treutler-Ahlrichs M4 radial quadrature. + * + * Taken as the convolution of the Gauss-Chebyshev (second kind) quadrature + * with the TA M4 radial transformation. See TreutlerAhlrichsRadialTraits for + * details. + * + * Suitable for integrands which tend to zero as their argument tends + f(r), with + * lim_{r->inf} f(r) = 0. + * + * Reference: + * J. Chem. Phys. 102, 346 (1995) + * DOI: https://doi.org/10.1063/1.469408 * - * @tparam PointType Type describing the quadrature points - * @tparam WeightType Type describing the quadrature weights + * @tparam PointType Type describing the quadrature points + * @tparam WeightType Type describing the quadrature weights */ - template -struct quadrature_traits< - TreutlerAhlrichs -> { - - using point_type = PointType; - using weight_type = WeightType; - - using point_container = std::vector< point_type >; - using weight_container = std::vector< weight_type >; - - /** - * @brief Generator for the Treutler-Ahlrichs quadrature - * - * @param[in] npts Number of quadrature points to generate - * @param[in] R Radial scaling factor. Table for suggested - * values is given in the original reference. - * @param[in] alpha Exponent factor for the quadrature. 0.6 was the - * - * @returns Tuple of quadrature points and weights - */ - inline static std::tuple - generate( size_t npts, weight_type R, weight_type alpha ) { - - - const point_type ln_2 = std::log(2.); - const point_type pi_ov_npts_p1 = M_PI / (npts + 1); - - point_container points( npts ); - weight_container weights( npts ); - - /* - * Treutler-Ahlrichs quadrature - * - * Original reference: - * J. Chem. Phys. 102, 346 (1995) - * DOI: https://doi.org/10.1063/1.469408 - * - * Closed form for points / weights obtained from: - * Journal of Computational Chemistry, 24: 732–740, 2003 - * DOI: https://doi.org/10.1002/jcc.10211 - * - */ - for( size_t i = 0; i < npts; ++i ) { - const auto xi = std::cos( (i+1) * pi_ov_npts_p1 ); - - const auto pow_term = std::pow( 1. + xi, alpha ) / ln_2; - const auto log_term = std::log( (1. - xi)/2. ); - const auto sqrt_term = std::sqrt( (1.+xi)/(1.-xi) ); - - points[i] = -R * pow_term * log_term; - weights[i] = R * pi_ov_npts_p1 * pow_term * - ( sqrt_term - alpha * log_term / sqrt_term ); - } - - - - return std::make_tuple( points, weights ); - - } - -}; +using TreutlerAhlrichs = RadialTransformQuadrature< + GaussChebyshev2, + TreutlerAhlrichsRadialTraits +>; } diff --git a/include/integratorxx/quadratures/uniform.hpp b/include/integratorxx/quadratures/uniform.hpp index e173fe2..27720f1 100644 --- a/include/integratorxx/quadratures/uniform.hpp +++ b/include/integratorxx/quadratures/uniform.hpp @@ -42,6 +42,8 @@ struct quadrature_traits< using point_container = std::vector< point_type >; using weight_container = std::vector< weight_type >; + inline static constexpr bool bound_inclusive = true; + inline static std::tuple generate( size_t npts, point_type lo, point_type up ) { @@ -60,6 +62,8 @@ struct quadrature_traits< } + inline static std::tuple + generate( size_t npts ) { return generate(npts, 0.0, 1.0); } }; } diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index c6c4f47..9374166 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -198,32 +198,32 @@ TEST_CASE( "Gauss-Chebyshev Quadratures", "[1d-quad]") { }; SECTION("First Kind") { - IntegratorXX::GaussChebyshev1 quad_even(200, -1., 1.); + IntegratorXX::GaussChebyshev1 quad_even(200); integrate(quad_even); - IntegratorXX::GaussChebyshev1 quad_odd(201, -1., 1.); + IntegratorXX::GaussChebyshev1 quad_odd(201); integrate(quad_odd); } SECTION("Second Kind") { - IntegratorXX::GaussChebyshev2 quad_even(200, -1., 1.); + IntegratorXX::GaussChebyshev2 quad_even(200); integrate(quad_even); - IntegratorXX::GaussChebyshev2 quad_odd(201, -1., 1.); + IntegratorXX::GaussChebyshev2 quad_odd(201); integrate(quad_odd); } SECTION("Second Kind (Modified)") { - IntegratorXX::GaussChebyshev2Modified quad_even(200, -1., 1.); + IntegratorXX::GaussChebyshev2Modified quad_even(200); integrate(quad_even); - IntegratorXX::GaussChebyshev2Modified quad_odd(201, -1., 1.); + IntegratorXX::GaussChebyshev2Modified quad_odd(201); integrate(quad_odd); } SECTION("Third Kind") { - IntegratorXX::GaussChebyshev3 quad(200, -1., 1.); + IntegratorXX::GaussChebyshev3 quad(200); integrate(quad); } } TEST_CASE( "Euler-Maclaurin Quadratures", "[1d-quad]" ) { - IntegratorXX::MurrayHandyLaming quad( 150 ); + IntegratorXX::MurrayHandyLaming quad(150); const auto& pts = quad.points(); const auto& wgt = quad.weights(); @@ -238,9 +238,9 @@ TEST_CASE( "Euler-Maclaurin Quadratures", "[1d-quad]" ) { } -TEST_CASE( "Ahlrichs Quadratures", "[1d-quad]" ) { +TEST_CASE( "Treutler-Ahlrichs Quadratures", "[1d-quad]" ) { - IntegratorXX::TreutlerAhlrichs quad( 150 ); + IntegratorXX::TreutlerAhlrichs quad(150); const auto& pts = quad.points(); const auto& wgt = quad.weights(); @@ -258,7 +258,7 @@ TEST_CASE( "Ahlrichs Quadratures", "[1d-quad]" ) { TEST_CASE( "Knowles Quadratures", "[1d-quad]" ) { - IntegratorXX::MuraKnowles quad( 150 ); + IntegratorXX::MuraKnowles quad(150); const auto& pts = quad.points(); const auto& wgt = quad.weights(); @@ -266,8 +266,9 @@ TEST_CASE( "Knowles Quadratures", "[1d-quad]" ) { auto f = [=]( double x ){ return gaussian(x); }; double res = 0.; - for( auto i = 0; i < quad.npts(); ++i ) + for( auto i = 0; i < quad.npts(); ++i ) res += wgt[i] * f(pts[i]); + CHECK( res == Catch::Approx(ref_gaussian_int(0.,inf)) ); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 543719a..fdb1566 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -30,6 +30,8 @@ target_link_libraries( gausslegendre PUBLIC Catch2::Catch2WithMain integratorxx add_executable( gausslobatto gausslobatto.cxx ) target_link_libraries( gausslobatto PUBLIC Catch2::Catch2WithMain integratorxx ) +add_executable( muraknowles muraknowles.cxx ) +target_link_libraries( muraknowles PUBLIC Catch2::Catch2WithMain integratorxx ) add_executable( composite_quadratures composite_quadratures.cxx ) target_link_libraries( composite_quadratures PUBLIC Catch2::Catch2WithMain integratorxx ) diff --git a/test/muraknowles.cxx b/test/muraknowles.cxx new file mode 100644 index 0000000..198eb4a --- /dev/null +++ b/test/muraknowles.cxx @@ -0,0 +1,65 @@ +// clang-format off +#include "catch2/catch_all.hpp" +#include +#include +#include +// clang-format on + +const double x_tolerance = std::numeric_limits::epsilon(); +const double w_tolerance = std::numeric_limits::epsilon(); + +// Reference data provided by Peter Knowles +// https://github.com/wavefunction91/IntegratorXX/pull/27#discussion_r1260345076 +// ALPHA = 7.0 +// M = 3 +// NPTS = 48 +TEST_CASE("48 point MuraKnowles", "[1d-quad]") { + std::array ref_pts = { + 5.9499271133939913E-005, 4.7600833032953915E-004, 1.6066578611874361E-003, + 3.8089732866334377E-003, 7.4413311256047504E-003, 1.2863600152286464E-002, + 2.0437970723343253E-002, 3.0529977777593031E-002, 4.3509725746882791E-002, + 5.9753327005997357E-002, 7.9644569620570760E-002, 0.10357683514824530, + 0.13195529327830460, 0.16519940738611444, 0.20374579393428979, + 0.24805148947543415, 0.29859769233640943, 0.35589406260236178, + 0.42048368472012887, 0.49294882318131050, 0.57391763505319449, + 0.66407204595981006, 0.76415705172765447, 0.87499178083047879, + 0.99748274937600856, 1.1326398697743505, 1.2815959495630309, + 1.4456306574120197, 1.6262002677869758, 1.8249749675453268, + 2.0438861838399740, 2.2851873784253796, 2.5515332181111825, + 2.8460842535535971, 3.1726476906379575, 3.5358703430273275, + 3.9415088946346315, 4.3968179569776069, 4.9111235220893743, + 5.4966995004301378, 6.1701626127089160, 6.9548035249514397, + 7.8847270471719106, 9.0128012057250508, 10.427589422901853, + 12.295111781246881, 14.988083164786680, 19.695799083317340}; + + std::array ref_wgt = { + 6.3191408757371128E-013, 1.6178925855460011E-010, 4.1478214302682156E-009, + 4.1457509371301411E-008, 2.4736257856667049E-007, 1.0652626596586621E-006, + 3.6641284849845252E-006, 1.0694442455644551E-005, 2.7541538616163128E-005, + 6.4278015117145217E-005, 1.3857027486190222E-004, 2.7986264953446042E-004, + 5.3525184281221503E-004, 9.7758247883188328E-004, 1.7164520965350302E-003, + 2.9130257653514745E-003, 4.7998471771280651E-003, 7.7072235182964691E-003, + 1.2098296752658244E-002, 1.8615653176199027E-002, 2.8143351896444287E-002, + 4.1889696901313278E-002, 6.1498123098954999E-002, 8.9196494268593482E-002, + 0.12799934660813353, 0.18198381348671958, 0.25666916599665018, + 0.35954374153109814, 0.50080416940484207, 0.69440463378479145, + 0.95956585458830923, 1.3229773146815327, 1.8220646710507786, + 2.5099284939617785, 3.4629679021480624, 4.7929337979063913, + 6.6665154642672082, 9.3381941641618038, 13.207427859666060, + 18.922636048645305, 27.580468647622130, 41.132941288449700, + 63.288684869795212, 101.72526373304323, 174.33040888875871, + 330.69364765060254, 753.72105349957747, 2659.6767864598341}; + + IntegratorXX::MuraKnowlesRadialTraits traits{7.0}; + IntegratorXX::MuraKnowles quad(48, traits); + + const auto& pts = quad.points(); + const auto& wgt = quad.weights(); + for(auto i = 0ul; i < 2; i++) { + REQUIRE_THAT(pts[i], Catch::Matchers::WithinAbs(ref_pts[i], x_tolerance)); + // Reference data contains the spherical Jacobian (r^2) + const auto rsq = pts[i] * pts[i]; + REQUIRE_THAT(rsq * wgt[i], + Catch::Matchers::WithinAbs(ref_wgt[i], w_tolerance)); + } +}