From e75531408f0a3d0f50d9034f90977b1df3e5da41 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Wed, 29 Nov 2023 11:30:28 -0800 Subject: [PATCH 1/6] Refactoring generator headers to include radial/s2 generators --- .../generators/radial_generator.hpp | 40 +++++++++++++++++ .../integratorxx/generators/s2_generator.hpp | 26 +++++++++++ .../generators/spherical_factory.hpp | 44 +------------------ .../quadratures/radial/radial_transform.hpp | 15 ++++--- 4 files changed, 76 insertions(+), 49 deletions(-) create mode 100644 include/integratorxx/generators/radial_generator.hpp create mode 100644 include/integratorxx/generators/s2_generator.hpp diff --git a/include/integratorxx/generators/radial_generator.hpp b/include/integratorxx/generators/radial_generator.hpp new file mode 100644 index 0000000..054aa6e --- /dev/null +++ b/include/integratorxx/generators/radial_generator.hpp @@ -0,0 +1,40 @@ +#pragma once +#include + +namespace IntegratorXX { + +/// High-level specification of radial quadratures +enum class RadialQuad : uint32_t { + Becke = 0x0010, + MurrayHandyLaming = 0x0020, + MuraKnowles = 0x0030, + TreutlerAhlrichs = 0x0040 +}; + +template +RadialQuad radial_from_type() { + if constexpr (detail::is_becke_v) return RadialQuad::Becke; + if constexpr (detail::is_mk_v ) return RadialQuad::MuraKnowles; + if constexpr (detail::is_mhl_v) return RadialQuad::MurrayHandyLaming; + if constexpr (detail::is_ta_v) return RadialQuad::TreutlerAhlrichs; + + throw std::runtime_error("Unrecognized Radial Quadrature"); +}; + +RadialQuad radial_from_string(std::string name); + +#if 0 +struct RadialFactory { + + using radial_grid_ptr = std::shared_ptr< + QuadratureBase< + std::vector, + std::vector + >; + + //radial_grid_ptr generate(RadialTraits traits); + +}; +#endif + +} diff --git a/include/integratorxx/generators/s2_generator.hpp b/include/integratorxx/generators/s2_generator.hpp new file mode 100644 index 0000000..978dd04 --- /dev/null +++ b/include/integratorxx/generators/s2_generator.hpp @@ -0,0 +1,26 @@ +#pragma once +#include + +namespace IntegratorXX { + +/// High-level specification of angular quadratures +enum class AngularQuad : uint32_t { + AhrensBeylkin = 0x0100, + Delley = 0x0200, + LebedevLaikov = 0x0300, + Womersley = 0x0400 +}; + +template +AngularQuad angular_from_type() { + if constexpr (detail::is_ahrens_beyklin_v) return AngularQuad::AhrensBeylkin; + if constexpr (detail::is_delley_v ) return AngularQuad::Delley; + if constexpr (detail::is_lebedev_laikov_v) return AngularQuad::LebedevLaikov; + if constexpr (detail::is_womersley_v) return AngularQuad::Womersley; + + throw std::runtime_error("Unrecognized Angular Quadrature"); +}; + +AngularQuad angular_from_string(std::string name); + +} diff --git a/include/integratorxx/generators/spherical_factory.hpp b/include/integratorxx/generators/spherical_factory.hpp index f4722ed..6d08748 100644 --- a/include/integratorxx/generators/spherical_factory.hpp +++ b/include/integratorxx/generators/spherical_factory.hpp @@ -2,53 +2,13 @@ #include #include #include -#include -#include +#include +#include #include #include namespace IntegratorXX { -/// High-level specification of radial quadratures -enum class RadialQuad : uint32_t { - Becke = 0x0010, - MurrayHandyLaming = 0x0020, - MuraKnowles = 0x0030, - TreutlerAhlrichs = 0x0040 -}; - -template -RadialQuad radial_from_type() { - if constexpr (detail::is_becke_v) return RadialQuad::Becke; - if constexpr (detail::is_mk_v ) return RadialQuad::MuraKnowles; - if constexpr (detail::is_mhl_v) return RadialQuad::MurrayHandyLaming; - if constexpr (detail::is_ta_v) return RadialQuad::TreutlerAhlrichs; - - throw std::runtime_error("Unrecognized Radial Quadrature"); -}; - -RadialQuad radial_from_string(std::string name); - -/// High-level specification of angular quadratures -enum class AngularQuad : uint32_t { - AhrensBeylkin = 0x0100, - Delley = 0x0200, - LebedevLaikov = 0x0300, - Womersley = 0x0400 -}; - -template -AngularQuad angular_from_type() { - if constexpr (detail::is_ahrens_beyklin_v) return AngularQuad::AhrensBeylkin; - if constexpr (detail::is_delley_v ) return AngularQuad::Delley; - if constexpr (detail::is_lebedev_laikov_v) return AngularQuad::LebedevLaikov; - if constexpr (detail::is_womersley_v) return AngularQuad::Womersley; - - throw std::runtime_error("Unrecognized Angular Quadrature"); -}; - -AngularQuad angular_from_string(std::string name); - /// High-level specification of pruning schemes for spherical quadratures enum class PruningScheme { Unpruned, /// Unpruned quadrature diff --git a/include/integratorxx/quadratures/radial/radial_transform.hpp b/include/integratorxx/quadratures/radial/radial_transform.hpp index 93f8096..aac2129 100644 --- a/include/integratorxx/quadratures/radial/radial_transform.hpp +++ b/include/integratorxx/quadratures/radial/radial_transform.hpp @@ -4,11 +4,11 @@ namespace IntegratorXX { -template +template struct RadialTransformQuadrature : - public Quadrature> { + public Quadrature> { - using base_type = Quadrature>; + using base_type = Quadrature>; public: @@ -17,7 +17,7 @@ struct RadialTransformQuadrature : using point_container = typename base_type::point_container; using weight_container = typename base_type::weight_container; - RadialTransformQuadrature(size_t npts, const RadialTraits& traits = RadialTraits()) : + RadialTransformQuadrature(size_t npts, const RadialTraitsType& traits = RadialTraitsType()) : base_type( npts, traits ) { } RadialTransformQuadrature( const RadialTransformQuadrature& ) = default; @@ -25,16 +25,17 @@ struct RadialTransformQuadrature : }; -template -struct quadrature_traits> { +template +struct quadrature_traits> { using point_type = typename BaseQuad::point_type; using weight_type = typename BaseQuad::weight_type; + using traits_type = RadialTraitsType; 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 ) { + generate( size_t npts, const RadialTraitsType& traits ) { using base_quad_traits = quadrature_traits; From 7c74827c142616eeb2a563580138425f8788bde0 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Wed, 29 Nov 2023 12:56:26 -0800 Subject: [PATCH 2/6] Move grid order into RadialTraits - add base type to allow for polymorphic instantiation in runtime generators --- .../integratorxx/quadratures/radial/becke.hpp | 10 ++++--- .../integratorxx/quadratures/radial/mhl.hpp | 8 ++++-- .../quadratures/radial/muraknowles.hpp | 7 +++-- .../quadratures/radial/radial_transform.hpp | 26 ++++++++++++++++--- .../quadratures/radial/treutlerahlrichs.hpp | 11 +++++--- test/muraknowles.cxx | 4 +-- 6 files changed, 50 insertions(+), 16 deletions(-) diff --git a/include/integratorxx/quadratures/radial/becke.hpp b/include/integratorxx/quadratures/radial/becke.hpp index c0606b1..8cc38fe 100644 --- a/include/integratorxx/quadratures/radial/becke.hpp +++ b/include/integratorxx/quadratures/radial/becke.hpp @@ -13,8 +13,10 @@ namespace IntegratorXX { * J. Chem. Phys. 88, 2547 (1988) * DOI: https://doi.org/10.1063/1.454033 */ -class BeckeRadialTraits { - double R_; +class BeckeRadialTraits : public RadialTraits { + + size_t npts_; ///< Number of grid points + double R_; ///< Radial scaling factor public: /** @@ -24,7 +26,9 @@ class BeckeRadialTraits { * * @param[in] R Radial scaling factor */ - BeckeRadialTraits(double R = 1.0) : R_(R) {} + BeckeRadialTraits(size_t npts, double R = 1.0) : npts_(npts), R_(R) {} + + size_t npts() const noexcept { return npts_; } /** * @brief Transformation rule for the Becke radial quadratures diff --git a/include/integratorxx/quadratures/radial/mhl.hpp b/include/integratorxx/quadratures/radial/mhl.hpp index 3347de3..cc57b36 100644 --- a/include/integratorxx/quadratures/radial/mhl.hpp +++ b/include/integratorxx/quadratures/radial/mhl.hpp @@ -17,13 +17,17 @@ namespace IntegratorXX { * Typically taken to be 2. */ template -class MurrayHandyLamingRadialTraits { +class MurrayHandyLamingRadialTraits : public RadialTraits { + size_t npts_; ///< Number of grid points double R_; ///< Radial scaling factor + public: - MurrayHandyLamingRadialTraits(double R = 1.0) : R_(R) {} + MurrayHandyLamingRadialTraits(size_t npts, double R = 1.0) : npts_(npts), R_(R) {} + + size_t npts() const noexcept { return npts_; } /** * @brief Transformation rule for the MHL radial quadrature diff --git a/include/integratorxx/quadratures/radial/muraknowles.hpp b/include/integratorxx/quadratures/radial/muraknowles.hpp index ae9d890..1582b16 100644 --- a/include/integratorxx/quadratures/radial/muraknowles.hpp +++ b/include/integratorxx/quadratures/radial/muraknowles.hpp @@ -136,13 +136,16 @@ struct quadrature_traits< #else -class MuraKnowlesRadialTraits { +class MuraKnowlesRadialTraits : public RadialTraits { + size_t npts_; ///< Number of grid points double R_; ///< Radial scaling factor public: - MuraKnowlesRadialTraits(double R = 1.0) : R_(R) { } + MuraKnowlesRadialTraits(size_t npts, double R = 1.0) : npts_(npts), R_(R) { } + + size_t npts() const noexcept { return npts_; } template inline auto radial_transform(PointType x) const noexcept { diff --git a/include/integratorxx/quadratures/radial/radial_transform.hpp b/include/integratorxx/quadratures/radial/radial_transform.hpp index aac2129..a8daf2e 100644 --- a/include/integratorxx/quadratures/radial/radial_transform.hpp +++ b/include/integratorxx/quadratures/radial/radial_transform.hpp @@ -4,6 +4,16 @@ namespace IntegratorXX { +// Base type for all radial traits +struct RadialTraits { + virtual ~RadialTraits() noexcept = default; +}; + +template +const RadialTraitsType& radial_traits_cast( const RadialTraits& traits ) { + return dynamic_cast(traits); +} + template struct RadialTransformQuadrature : public Quadrature> { @@ -17,8 +27,17 @@ struct RadialTransformQuadrature : using point_container = typename base_type::point_container; using weight_container = typename base_type::weight_container; - RadialTransformQuadrature(size_t npts, const RadialTraitsType& traits = RadialTraitsType()) : - base_type( npts, traits ) { } + RadialTransformQuadrature(const RadialTraitsType& traits = RadialTraitsType()) : + base_type(traits) { } + RadialTransformQuadrature(const RadialTraits& traits) : + RadialTransformQuadrature(radial_traits_cast(traits)) { } + + template 1)> + > + RadialTransformQuadrature(Args&&... args) : + RadialTransformQuadrature(RadialTraitsType(std::forward(args)...)) { } + RadialTransformQuadrature( const RadialTransformQuadrature& ) = default; RadialTransformQuadrature( RadialTransformQuadrature&& ) noexcept = default; @@ -35,10 +54,11 @@ struct quadrature_traits> using weight_container = typename BaseQuad::weight_container; inline static std::tuple - generate( size_t npts, const RadialTraitsType& traits ) { + generate( const RadialTraitsType& traits ) { using base_quad_traits = quadrature_traits; + const auto npts = traits.npts(); point_container points( npts ); weight_container weights( npts ); diff --git a/include/integratorxx/quadratures/radial/treutlerahlrichs.hpp b/include/integratorxx/quadratures/radial/treutlerahlrichs.hpp index e2788f2..c10cd99 100644 --- a/include/integratorxx/quadratures/radial/treutlerahlrichs.hpp +++ b/include/integratorxx/quadratures/radial/treutlerahlrichs.hpp @@ -13,9 +13,10 @@ namespace IntegratorXX { * J. Chem. Phys. 102, 346 (1995) * DOI: https://doi.org/10.1063/1.469408 */ -class TreutlerAhlrichsRadialTraits { +class TreutlerAhlrichsRadialTraits : public RadialTraits { - double R_; + size_t npts_; ///< Number of grid points + double R_; ///< Radial scaling factor double alpha_; public: @@ -34,8 +35,10 @@ class TreutlerAhlrichsRadialTraits { * @param[in] R Radial scaling factor * @param[in] alpha TA exponential factor */ - TreutlerAhlrichsRadialTraits(double R = 1.0, double alpha = 0.6) : - R_(R), alpha_(alpha) { } + TreutlerAhlrichsRadialTraits(size_t npts, double R = 1.0, double alpha = 0.6) : + npts_(npts), R_(R), alpha_(alpha) { } + + size_t npts() const noexcept { return npts_; } /** * @brief Transformation rule for the TA M3+M4 radial quadratures diff --git a/test/muraknowles.cxx b/test/muraknowles.cxx index 0e178f5..52f378c 100644 --- a/test/muraknowles.cxx +++ b/test/muraknowles.cxx @@ -50,8 +50,8 @@ TEST_CASE("48 point MuraKnowles", "[1d-quad]") { 63.288684869795212, 101.72526373304323, 174.33040888875871, 330.69364765060254, 753.72105349957747, 2659.6767864598341}; - IntegratorXX::MuraKnowlesRadialTraits traits{7.0}; - IntegratorXX::MuraKnowles quad(48, traits); + IntegratorXX::MuraKnowlesRadialTraits traits{48, 7.0}; + IntegratorXX::MuraKnowles quad(traits); const auto& pts = quad.points(); const auto& wgt = quad.weights(); From b0ec520835314e3f0a14edfe7fffc21f97ec1b3a Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Thu, 30 Nov 2023 10:16:07 -0800 Subject: [PATCH 3/6] Move Spherical Generator over to RadialTraits --- .../generators/radial_generator.hpp | 39 ++++++++ .../generators/spherical_factory.hpp | 54 +++++++--- .../integratorxx/quadratures/radial/becke.hpp | 13 +++ .../integratorxx/quadratures/radial/mhl.hpp | 13 +++ .../quadratures/radial/muraknowles.hpp | 13 +++ .../quadratures/radial/radial_transform.hpp | 5 + .../quadratures/radial/treutlerahlrichs.hpp | 13 +++ src/spherical_factory.cxx | 70 +++++++------ test/spherical_generator.cxx | 98 ++++++++++--------- 9 files changed, 229 insertions(+), 89 deletions(-) diff --git a/include/integratorxx/generators/radial_generator.hpp b/include/integratorxx/generators/radial_generator.hpp index 054aa6e..a60c71e 100644 --- a/include/integratorxx/generators/radial_generator.hpp +++ b/include/integratorxx/generators/radial_generator.hpp @@ -23,6 +23,45 @@ RadialQuad radial_from_type() { RadialQuad radial_from_string(std::string name); +namespace detail { + +template +std::unique_ptr make_radial_traits(Args&&... args) { + using traits_type = RadialTraitsType; + if constexpr (std::is_constructible_v) + return std::make_unique(std::forward(args)...); + else return nullptr; +} + +} + +template +std::unique_ptr make_radial_traits(RadialQuad rq, Args&&... args) { + std::unique_ptr ptr; + switch(rq) { + case RadialQuad::Becke: + ptr = + detail::make_radial_traits(std::forward(args)...); + break; + case RadialQuad::MurrayHandyLaming: + ptr = + detail::make_radial_traits>(std::forward(args)...); + break; + case RadialQuad::MuraKnowles: + ptr = + detail::make_radial_traits(std::forward(args)...); + break; + case RadialQuad::TreutlerAhlrichs: + ptr = + detail::make_radial_traits(std::forward(args)...); + break; + } + + if(!ptr) throw std::runtime_error("RadialTraits Construction Failed"); + return ptr; +} + + #if 0 struct RadialFactory { diff --git a/include/integratorxx/generators/spherical_factory.hpp b/include/integratorxx/generators/spherical_factory.hpp index 6d08748..1b6f8f8 100644 --- a/include/integratorxx/generators/spherical_factory.hpp +++ b/include/integratorxx/generators/spherical_factory.hpp @@ -17,22 +17,29 @@ enum class PruningScheme { }; // TODO: Make these strong (non-convertible) types -using RadialScale = double; -using RadialSize = size_t; +//using RadialScale = double; +//using RadialSize = size_t; using AngularSize = size_t; - +using radial_traits_ptr = std::unique_ptr; /// Generic specification of an unpruned spherical quadrature struct UnprunedSphericalGridSpecification { - RadialQuad radial_quad; ///< Radial quadrature specification - RadialSize radial_size; ///< Number of radial quadrature points - RadialScale radial_scale; ///< Radial scaling factor + RadialQuad radial_quad; ///< Radial quadrature specification + radial_traits_ptr radial_traits; ///< Radial traits (order, scaling factors, etc) AngularQuad angular_quad; /// Angular quadrature specification AngularSize angular_size; /// Number of angular quadrature points + + UnprunedSphericalGridSpecification(RadialQuad, const RadialTraits&, AngularQuad, + AngularSize); + + UnprunedSphericalGridSpecification(const UnprunedSphericalGridSpecification& other) : + radial_quad(other.radial_quad), radial_traits(other.radial_traits ? other.radial_traits->clone() : nullptr), + angular_quad(other.angular_quad), angular_size(other.angular_size) {}; }; + /// Specification of a pruned region of an spherical quadrature struct PruningRegion { size_t idx_st; ///< Starting radial index for pruned region @@ -51,15 +58,34 @@ struct PruningRegion { struct PrunedSphericalGridSpecification { RadialQuad radial_quad; ///< Radial quadrature specification - RadialSize radial_size; ///< Number of radial quadrature points - RadialScale radial_scale; ///< Radial scaling factor + radial_traits_ptr radial_traits; ///< Radial traits (order, scaling factors, etc) std::vector pruning_regions; ///< List of pruning regions over the radial quadrature + + PrunedSphericalGridSpecification() = default; + + template + PrunedSphericalGridSpecification(RadialQuad rq, radial_traits_ptr&& traits, Arg&&... arg) : + radial_quad(rq), radial_traits(std::move(traits)), pruning_regions(std::forward(arg)...) { } + template + PrunedSphericalGridSpecification(RadialQuad rq, const RadialTraits& traits, Arg&&... arg) : + PrunedSphericalGridSpecification(rq, traits.clone(), std::forward(arg)...) { } + + PrunedSphericalGridSpecification(const PrunedSphericalGridSpecification& other) : + PrunedSphericalGridSpecification(other.radial_quad, + other.radial_traits ? other.radial_traits->clone() : nullptr, + other.pruning_regions) { } + + PrunedSphericalGridSpecification& operator=(const PrunedSphericalGridSpecification& other) { + radial_quad = other.radial_quad; + radial_traits = other.radial_traits ? other.radial_traits->clone() : nullptr; + pruning_regions = other.pruning_regions; + return *this; + } inline bool operator==(const PrunedSphericalGridSpecification& other) const noexcept { return radial_quad == other.radial_quad and - radial_size == other.radial_size and - radial_scale == other.radial_scale and + (radial_traits ? (other.radial_traits and radial_traits->compare(*other.radial_traits)) : !other.radial_traits) and pruning_regions == other.pruning_regions; } }; @@ -136,10 +162,10 @@ struct SphericalGridFactory { } - static spherical_grid_ptr generate_unpruned_grid( RadialQuad, RadialSize, - RadialScale, AngularQuad, AngularSize ); - static spherical_grid_ptr generate_pruned_grid( RadialQuad, RadialSize, - RadialScale, const std::vector&); + static spherical_grid_ptr generate_unpruned_grid( RadialQuad, const RadialTraits&, + AngularQuad, AngularSize ); + static spherical_grid_ptr generate_pruned_grid( RadialQuad, const RadialTraits&, + const std::vector&); static spherical_grid_ptr generate_grid(UnprunedSphericalGridSpecification gs); diff --git a/include/integratorxx/quadratures/radial/becke.hpp b/include/integratorxx/quadratures/radial/becke.hpp index 8cc38fe..98fb2f9 100644 --- a/include/integratorxx/quadratures/radial/becke.hpp +++ b/include/integratorxx/quadratures/radial/becke.hpp @@ -30,6 +30,19 @@ class BeckeRadialTraits : public RadialTraits { size_t npts() const noexcept { return npts_; } + std::unique_ptr clone() const { + return std::make_unique(*this); + } + + bool compare(const RadialTraits& other) const noexcept { + auto ptr = dynamic_cast(&other); + return ptr ? *this == *ptr : false; + } + + bool operator==(const BeckeRadialTraits& other) const noexcept { + return npts_ == other.npts_ and R_ == other.R_; + } + /** * @brief Transformation rule for the Becke radial quadratures * diff --git a/include/integratorxx/quadratures/radial/mhl.hpp b/include/integratorxx/quadratures/radial/mhl.hpp index cc57b36..7065ff1 100644 --- a/include/integratorxx/quadratures/radial/mhl.hpp +++ b/include/integratorxx/quadratures/radial/mhl.hpp @@ -29,6 +29,19 @@ class MurrayHandyLamingRadialTraits : public RadialTraits { size_t npts() const noexcept { return npts_; } + std::unique_ptr clone() const { + return std::make_unique(*this); + } + + bool compare(const RadialTraits& other) const noexcept { + auto ptr = dynamic_cast(&other); + return ptr ? *this == *ptr : false; + } + + bool operator==(const MurrayHandyLamingRadialTraits& other) const noexcept { + return npts_ == other.npts_ and R_ == other.R_; + } + /** * @brief Transformation rule for the MHL radial quadrature * diff --git a/include/integratorxx/quadratures/radial/muraknowles.hpp b/include/integratorxx/quadratures/radial/muraknowles.hpp index 1582b16..9d242e4 100644 --- a/include/integratorxx/quadratures/radial/muraknowles.hpp +++ b/include/integratorxx/quadratures/radial/muraknowles.hpp @@ -147,6 +147,19 @@ class MuraKnowlesRadialTraits : public RadialTraits { size_t npts() const noexcept { return npts_; } + std::unique_ptr clone() const { + return std::make_unique(*this); + } + + bool compare(const RadialTraits& other) const noexcept { + auto ptr = dynamic_cast(&other); + return ptr ? *this == *ptr : false; + } + + bool operator==(const MuraKnowlesRadialTraits& other) const noexcept { + return npts_ == other.npts_ and R_ == other.R_; + } + template inline auto radial_transform(PointType x) const noexcept { return -R_ * std::log(1.0 - x*x*x); diff --git a/include/integratorxx/quadratures/radial/radial_transform.hpp b/include/integratorxx/quadratures/radial/radial_transform.hpp index a8daf2e..a319a49 100644 --- a/include/integratorxx/quadratures/radial/radial_transform.hpp +++ b/include/integratorxx/quadratures/radial/radial_transform.hpp @@ -1,5 +1,7 @@ #pragma once +#include + #include namespace IntegratorXX { @@ -7,6 +9,9 @@ namespace IntegratorXX { // Base type for all radial traits struct RadialTraits { virtual ~RadialTraits() noexcept = default; + virtual std::unique_ptr clone() const = 0; + virtual size_t npts() const = 0; + virtual bool compare( const RadialTraits& ) const noexcept = 0; }; template diff --git a/include/integratorxx/quadratures/radial/treutlerahlrichs.hpp b/include/integratorxx/quadratures/radial/treutlerahlrichs.hpp index c10cd99..c22dab0 100644 --- a/include/integratorxx/quadratures/radial/treutlerahlrichs.hpp +++ b/include/integratorxx/quadratures/radial/treutlerahlrichs.hpp @@ -40,6 +40,19 @@ class TreutlerAhlrichsRadialTraits : public RadialTraits { size_t npts() const noexcept { return npts_; } + std::unique_ptr clone() const { + return std::make_unique(*this); + } + + bool compare(const RadialTraits& other) const noexcept { + auto ptr = dynamic_cast(&other); + return ptr ? *this == *ptr : false; + } + + bool operator==(const TreutlerAhlrichsRadialTraits& other) const noexcept { + return npts_ == other.npts_ and R_ == other.R_ and alpha_ == other.alpha_; + } + /** * @brief Transformation rule for the TA M3+M4 radial quadratures * diff --git a/src/spherical_factory.cxx b/src/spherical_factory.cxx index c270866..34007a4 100644 --- a/src/spherical_factory.cxx +++ b/src/spherical_factory.cxx @@ -7,6 +7,11 @@ namespace IntegratorXX { +UnprunedSphericalGridSpecification::UnprunedSphericalGridSpecification( + RadialQuad rq, const RadialTraits& traits, AngularQuad aq, AngularSize as) : + radial_quad(rq), radial_traits(traits.clone()), angular_quad(aq), angular_size(as) { } + + RadialQuad radial_from_string(std::string name) { std::transform(name.begin(), name.end(), name.begin(), ::toupper); if(name == "BECKE") return RadialQuad::Becke; @@ -50,22 +55,22 @@ using wo_type = Womersley; /**** Unpruned Grids ****/ /************************/ template -spherical_grid_ptr generate_unpruned_grid_impl(RadialQuad rq, RadialSize nrad, - RadialScale rscal, AngularQuadType&& ang_quad) { +spherical_grid_ptr generate_unpruned_grid_impl(RadialQuad rq, const RadialTraits& traits, + AngularQuadType&& ang_quad) { switch( rq ) { case RadialQuad::Becke: - return SphericalGridFactory::generate_unpruned_grid( bk_type(nrad, rscal), std::forward(ang_quad) ); + return SphericalGridFactory::generate_unpruned_grid( bk_type(traits), std::forward(ang_quad) ); case RadialQuad::MuraKnowles: - return SphericalGridFactory::generate_unpruned_grid( mk_type(nrad, rscal), std::forward(ang_quad) ); + return SphericalGridFactory::generate_unpruned_grid( mk_type(traits), std::forward(ang_quad) ); case RadialQuad::MurrayHandyLaming: - return SphericalGridFactory::generate_unpruned_grid( mhl_type(nrad, rscal), std::forward(ang_quad) ); + return SphericalGridFactory::generate_unpruned_grid( mhl_type(traits), std::forward(ang_quad) ); case RadialQuad::TreutlerAhlrichs: - return SphericalGridFactory::generate_unpruned_grid( ta_type(nrad, rscal), std::forward(ang_quad) ); + return SphericalGridFactory::generate_unpruned_grid( ta_type(traits), std::forward(ang_quad) ); default: throw std::runtime_error("Unsupported Radial Quadrature"); @@ -77,17 +82,17 @@ spherical_grid_ptr generate_unpruned_grid_impl(RadialQuad rq, RadialSize nrad, } spherical_grid_ptr SphericalGridFactory::generate_unpruned_grid( RadialQuad rq, - RadialSize nrad, RadialScale rscal, AngularQuad aq, AngularSize nang) { + const RadialTraits& traits, AngularQuad aq, AngularSize nang) { switch(aq) { case AngularQuad::AhrensBeylkin: - return generate_unpruned_grid_impl(rq, nrad, rscal, ah_type(nang)); + return generate_unpruned_grid_impl(rq, traits, ah_type(nang)); case AngularQuad::Delley: - return generate_unpruned_grid_impl(rq, nrad, rscal, de_type(nang)); + return generate_unpruned_grid_impl(rq, traits, de_type(nang)); case AngularQuad::LebedevLaikov: - return generate_unpruned_grid_impl(rq, nrad, rscal, ll_type(nang)); + return generate_unpruned_grid_impl(rq, traits, ll_type(nang)); case AngularQuad::Womersley: - return generate_unpruned_grid_impl(rq, nrad, rscal, wo_type(nang)); + return generate_unpruned_grid_impl(rq, traits, wo_type(nang)); default: throw std::runtime_error("Unsupported Angular Quadrature"); abort(); @@ -96,7 +101,8 @@ spherical_grid_ptr SphericalGridFactory::generate_unpruned_grid( RadialQuad rq, spherical_grid_ptr SphericalGridFactory::generate_grid( UnprunedSphericalGridSpecification gs ) { - return generate_unpruned_grid(gs.radial_quad, gs.radial_size, gs.radial_scale, + if(!gs.radial_traits) throw std::runtime_error("RadialTraits Not Set"); + return generate_unpruned_grid(gs.radial_quad, *gs.radial_traits, gs.angular_quad, gs.angular_size); } @@ -152,19 +158,19 @@ auto make_pruned_grid(const RadialQuadType& rq, } spherical_grid_ptr SphericalGridFactory::generate_pruned_grid( - RadialQuad rq, RadialSize nrad, RadialScale rscal, + RadialQuad rq, const RadialTraits& traits, const std::vector& pruning_regions) { switch( rq ) { case RadialQuad::Becke: - return make_pruned_grid( bk_type(nrad, rscal), pruning_regions ); + return make_pruned_grid( bk_type(traits), pruning_regions ); case RadialQuad::MuraKnowles: - return make_pruned_grid( mk_type(nrad, rscal), pruning_regions ); + return make_pruned_grid( mk_type(traits), pruning_regions ); case RadialQuad::MurrayHandyLaming: - return make_pruned_grid( mhl_type(nrad, rscal), pruning_regions ); + return make_pruned_grid( mhl_type(traits), pruning_regions ); case RadialQuad::TreutlerAhlrichs: - return make_pruned_grid( ta_type(nrad, rscal), pruning_regions ); + return make_pruned_grid( ta_type(traits), pruning_regions ); default: throw std::runtime_error("Unsupported Radial Quadrature"); @@ -176,8 +182,9 @@ spherical_grid_ptr SphericalGridFactory::generate_pruned_grid( spherical_grid_ptr SphericalGridFactory::generate_grid( PrunedSphericalGridSpecification gs ) { - return generate_pruned_grid( gs.radial_quad, gs.radial_size, - gs.radial_scale, gs.pruning_regions ); + if(!gs.radial_traits) throw std::runtime_error("RadialTraits Not Set"); + return generate_pruned_grid( gs.radial_quad, *gs.radial_traits, + gs.pruning_regions ); } @@ -212,7 +219,7 @@ PrunedSphericalGridSpecification robust_psi4_pruning_scheme_impl( size_t low_sz, size_t med_sz, AngularQuad angular_quad, UnprunedSphericalGridSpecification unp ) { - const size_t rsz = unp.radial_size; + const size_t rsz = unp.radial_traits->npts(); const size_t r_div_4 = rsz / 4ul + 1ul; const size_t r_div_2 = rsz / 2ul + 1ul; std::vector pruning_regions = { @@ -221,9 +228,9 @@ PrunedSphericalGridSpecification robust_psi4_pruning_scheme_impl( {r_div_2, rsz, angular_quad, unp.angular_size} }; - return PrunedSphericalGridSpecification{ - unp.radial_quad, unp.radial_size, unp.radial_scale, pruning_regions - }; + return PrunedSphericalGridSpecification( + unp.radial_quad, unp.radial_traits->clone(), pruning_regions + ); } @@ -272,7 +279,7 @@ PrunedSphericalGridSpecification treutler_pruning_scheme_impl( size_t low_sz, size_t med_sz, AngularQuad angular_quad, UnprunedSphericalGridSpecification unp ) { - const size_t rsz = unp.radial_size; + const size_t rsz = unp.radial_traits->npts(); const size_t r_div_3 = rsz / 3ul + 1ul; const size_t r_div_2 = rsz / 2ul + 1ul; std::vector pruning_regions = { @@ -281,9 +288,9 @@ PrunedSphericalGridSpecification treutler_pruning_scheme_impl( {r_div_2, rsz, angular_quad, unp.angular_size} }; - return PrunedSphericalGridSpecification{ - unp.radial_quad, unp.radial_size, unp.radial_scale, pruning_regions - }; + return PrunedSphericalGridSpecification( + unp.radial_quad, unp.radial_traits->clone(), pruning_regions + ); } @@ -321,6 +328,7 @@ PrunedSphericalGridSpecification create_pruned_spec( PruningScheme scheme, UnprunedSphericalGridSpecification unp ) { + if(!unp.radial_traits) throw std::runtime_error("RadialTraits Not Set"); switch(scheme) { case PruningScheme::Robust: return robust_psi4_pruning_scheme(unp); @@ -331,11 +339,11 @@ PrunedSphericalGridSpecification create_pruned_spec( case PruningScheme::Unpruned: default: std::vector pruning_regions = { - {0ul, (size_t)unp.radial_size, unp.angular_quad, unp.angular_size} - }; - return PrunedSphericalGridSpecification{ - unp.radial_quad, unp.radial_size, unp.radial_scale, pruning_regions + {0ul, unp.radial_traits->npts(), unp.angular_quad, unp.angular_size} }; + return PrunedSphericalGridSpecification( + unp.radial_quad, unp.radial_traits->clone(), pruning_regions + ); } } diff --git a/test/spherical_generator.cxx b/test/spherical_generator.cxx index 650a766..acf6290 100644 --- a/test/spherical_generator.cxx +++ b/test/spherical_generator.cxx @@ -80,37 +80,43 @@ TEST_CASE( "Pruning Schemes", "[sph-gen]" ) { using namespace IntegratorXX; + auto rad_traits = make_radial_traits(RadialQuad::MuraKnowles, 99, 1.0); + SECTION("Ahrens-Beylkin") { - UnprunedSphericalGridSpecification unp{ - RadialQuad::MuraKnowles, 99, 1.0, + //UnprunedSphericalGridSpecification unp = make_unpruned_spec( + // AngularQuad::AhrensBeylkin, 372, + // RadialQuad::MuraKnowles, 99, 1.0 + //); + UnprunedSphericalGridSpecification unp( + RadialQuad::MuraKnowles, *rad_traits, AngularQuad::AhrensBeylkin, 372 - }; + ); PrunedSphericalGridSpecification pruning_spec, ref_pruning_spec; SECTION("Robust") { pruning_spec = robust_psi4_pruning_scheme(unp); - ref_pruning_spec = PrunedSphericalGridSpecification{ - RadialQuad::MuraKnowles, 99, 1.0, - { + ref_pruning_spec = PrunedSphericalGridSpecification( + RadialQuad::MuraKnowles, *rad_traits, + std::vector{ {0 , 25, AngularQuad::AhrensBeylkin, 72}, {25, 50, AngularQuad::AhrensBeylkin, 312}, {50, 99, AngularQuad::AhrensBeylkin, 372}, } - }; + ); } SECTION("Treutler") { pruning_spec = treutler_pruning_scheme(unp); - ref_pruning_spec = PrunedSphericalGridSpecification{ - RadialQuad::MuraKnowles, 99, 1.0, - { + ref_pruning_spec = PrunedSphericalGridSpecification( + RadialQuad::MuraKnowles, *rad_traits, + std::vector{ {0 , 34, AngularQuad::AhrensBeylkin, 72}, {34, 50, AngularQuad::AhrensBeylkin, 72}, {50, 99, AngularQuad::AhrensBeylkin, 372}, } - }; + ); } REQUIRE(pruning_spec == ref_pruning_spec); @@ -119,7 +125,7 @@ TEST_CASE( "Pruning Schemes", "[sph-gen]" ) { SECTION("Delley") { UnprunedSphericalGridSpecification unp{ - RadialQuad::MuraKnowles, 99, 1.0, + RadialQuad::MuraKnowles, *rad_traits, AngularQuad::Delley, 302 }; @@ -127,26 +133,26 @@ TEST_CASE( "Pruning Schemes", "[sph-gen]" ) { SECTION("Robust") { pruning_spec = robust_psi4_pruning_scheme(unp); - ref_pruning_spec = PrunedSphericalGridSpecification{ - RadialQuad::MuraKnowles, 99, 1.0, - { + ref_pruning_spec = PrunedSphericalGridSpecification( + RadialQuad::MuraKnowles, *rad_traits, + std::vector{ {0 , 25, AngularQuad::Delley, 26}, {25, 50, AngularQuad::Delley, 194}, {50, 99, AngularQuad::Delley, 302}, } - }; + ); } SECTION("Treutler") { pruning_spec = treutler_pruning_scheme(unp); - ref_pruning_spec = PrunedSphericalGridSpecification{ - RadialQuad::MuraKnowles, 99, 1.0, - { + ref_pruning_spec = PrunedSphericalGridSpecification( + RadialQuad::MuraKnowles, *rad_traits, + std::vector{ {0 , 34, AngularQuad::Delley, 26}, {34, 50, AngularQuad::Delley, 50}, {50, 99, AngularQuad::Delley, 302}, } - }; + ); } REQUIRE(pruning_spec == ref_pruning_spec); @@ -155,7 +161,7 @@ TEST_CASE( "Pruning Schemes", "[sph-gen]" ) { SECTION("Lebedev-Laikov") { UnprunedSphericalGridSpecification unp{ - RadialQuad::MuraKnowles, 99, 1.0, + RadialQuad::MuraKnowles, *rad_traits, AngularQuad::LebedevLaikov, 302 }; @@ -163,26 +169,26 @@ TEST_CASE( "Pruning Schemes", "[sph-gen]" ) { SECTION("Robust") { pruning_spec = robust_psi4_pruning_scheme(unp); - ref_pruning_spec = PrunedSphericalGridSpecification{ - RadialQuad::MuraKnowles, 99, 1.0, - { + ref_pruning_spec = PrunedSphericalGridSpecification( + RadialQuad::MuraKnowles, *rad_traits, + std::vector{ {0 , 25, AngularQuad::LebedevLaikov, 26}, {25, 50, AngularQuad::LebedevLaikov, 194}, {50, 99, AngularQuad::LebedevLaikov, 302}, } - }; + ); } SECTION("Treutler") { pruning_spec = treutler_pruning_scheme(unp); - ref_pruning_spec = PrunedSphericalGridSpecification{ - RadialQuad::MuraKnowles, 99, 1.0, - { + ref_pruning_spec = PrunedSphericalGridSpecification( + RadialQuad::MuraKnowles, *rad_traits, + std::vector{ {0 , 34, AngularQuad::LebedevLaikov, 26}, {34, 50, AngularQuad::LebedevLaikov, 50}, {50, 99, AngularQuad::LebedevLaikov, 302}, } - }; + ); } REQUIRE(pruning_spec == ref_pruning_spec); @@ -191,7 +197,7 @@ TEST_CASE( "Pruning Schemes", "[sph-gen]" ) { SECTION("Womersley") { UnprunedSphericalGridSpecification unp{ - RadialQuad::MuraKnowles, 99, 1.0, + RadialQuad::MuraKnowles, *rad_traits, AngularQuad::Womersley, 339 }; @@ -199,26 +205,26 @@ TEST_CASE( "Pruning Schemes", "[sph-gen]" ) { SECTION("Robust") { pruning_spec = robust_psi4_pruning_scheme(unp); - ref_pruning_spec = PrunedSphericalGridSpecification{ - RadialQuad::MuraKnowles, 99, 1.0, - { + ref_pruning_spec = PrunedSphericalGridSpecification( + RadialQuad::MuraKnowles, *rad_traits, + std::vector{ {0 , 25, AngularQuad::Womersley, 32}, {25, 50, AngularQuad::Womersley, 201}, {50, 99, AngularQuad::Womersley, 339}, } - }; + ); } SECTION("Treutler") { pruning_spec = treutler_pruning_scheme(unp); - ref_pruning_spec = PrunedSphericalGridSpecification{ - RadialQuad::MuraKnowles, 99, 1.0, - { + ref_pruning_spec = PrunedSphericalGridSpecification( + RadialQuad::MuraKnowles, *rad_traits, + std::vector{ {0 , 34, AngularQuad::Womersley, 32}, {34, 50, AngularQuad::Womersley, 72}, {50, 99, AngularQuad::Womersley, 339}, } - }; + ); } REQUIRE(pruning_spec == ref_pruning_spec); @@ -269,10 +275,12 @@ TEMPLATE_LIST_TEST_CASE("Unpruned", "[sph-gen]", sph_test_types) { // Generate via runtime API - UnprunedSphericalGridSpecification unp{ - radial_from_type(), nrad, 1.0, + auto rad_spec = radial_from_type(); + auto rad_traits = make_radial_traits(rad_spec, nrad, 1.0); + UnprunedSphericalGridSpecification unp( + rad_spec, *rad_traits, angular_from_type(), nang - }; + ); auto sph = SphericalGridFactory::generate_grid(unp); @@ -307,10 +315,12 @@ TEMPLATE_LIST_TEST_CASE("Pruned", "[sph-gen]", sph_test_types) { angular_traits::next_algebraic_order(29)); // Smallest possible angular grid // Generate pruning scheme - UnprunedSphericalGridSpecification unp{ - radial_from_type(), nrad, 1.0, + auto rad_spec = radial_from_type(); + auto rad_traits = make_radial_traits(rad_spec, nrad, 1.0); + UnprunedSphericalGridSpecification unp( + rad_spec, *rad_traits, angular_from_type(), nang - }; + ); auto pruning_spec = robust_psi4_pruning_scheme(unp); From 82d297123288d058720ddefc5ba89d20ba49a871 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Thu, 30 Nov 2023 10:19:02 -0800 Subject: [PATCH 4/6] Add MuraKnowles test to CTest --- test/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 67f2e00..7adc743 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -53,5 +53,6 @@ add_test( NAME QUADRATURES_COMPOSITE COMMAND composite_quadratures ) add_test( NAME QUADRATURES_SPH_GEN COMMAND spherical_generator ) add_test( NAME QUADRATURES_LEGENDRE COMMAND gausslegendre ) add_test( NAME QUADRATURES_LOBATTO COMMAND gausslobatto ) +add_test( NAME QUADRATURES_MURAKNOWLES COMMAND muraknowles ) add_test( NAME QUADRATURES_CHEBYSHEV1 COMMAND gausschebyshev1 ) add_test( NAME QUADRATURES_CHEBYSHEV2 COMMAND gausschebyshev2 ) From bbf67a8a5e00098bd489d850273fecef35590605 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Thu, 30 Nov 2023 11:10:21 -0800 Subject: [PATCH 5/6] Add RadialFactory + UTs --- ...adial_generator.hpp => radial_factory.hpp} | 7 ++-- .../{s2_generator.hpp => s2_factory.hpp} | 0 .../generators/spherical_factory.hpp | 4 +- src/CMakeLists.txt | 5 ++- src/radial_factory.cxx | 41 +++++++++++++++++++ src/radial_types.hpp | 12 ++++++ src/spherical_factory.cxx | 20 +-------- test/spherical_generator.cxx | 27 ++++++++++++ 8 files changed, 91 insertions(+), 25 deletions(-) rename include/integratorxx/generators/{radial_generator.hpp => radial_factory.hpp} (96%) rename include/integratorxx/generators/{s2_generator.hpp => s2_factory.hpp} (100%) create mode 100644 src/radial_factory.cxx create mode 100644 src/radial_types.hpp diff --git a/include/integratorxx/generators/radial_generator.hpp b/include/integratorxx/generators/radial_factory.hpp similarity index 96% rename from include/integratorxx/generators/radial_generator.hpp rename to include/integratorxx/generators/radial_factory.hpp index a60c71e..2d4f0ab 100644 --- a/include/integratorxx/generators/radial_generator.hpp +++ b/include/integratorxx/generators/radial_factory.hpp @@ -62,18 +62,17 @@ std::unique_ptr make_radial_traits(RadialQuad rq, Args&&... args) } -#if 0 struct RadialFactory { using radial_grid_ptr = std::shared_ptr< QuadratureBase< std::vector, std::vector - >; + > + >; - //radial_grid_ptr generate(RadialTraits traits); + static radial_grid_ptr generate(RadialQuad rq, const RadialTraits& traits); }; -#endif } diff --git a/include/integratorxx/generators/s2_generator.hpp b/include/integratorxx/generators/s2_factory.hpp similarity index 100% rename from include/integratorxx/generators/s2_generator.hpp rename to include/integratorxx/generators/s2_factory.hpp diff --git a/include/integratorxx/generators/spherical_factory.hpp b/include/integratorxx/generators/spherical_factory.hpp index 1b6f8f8..91f391d 100644 --- a/include/integratorxx/generators/spherical_factory.hpp +++ b/include/integratorxx/generators/spherical_factory.hpp @@ -2,8 +2,8 @@ #include #include #include -#include -#include +#include +#include #include #include diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 6806f6b..52501da 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,4 +1,7 @@ -add_library(integratorxx spherical_factory.cxx) +add_library(integratorxx + spherical_factory.cxx + radial_factory.cxx +) # Target features target_compile_features( integratorxx PUBLIC cxx_std_17 ) diff --git a/src/radial_factory.cxx b/src/radial_factory.cxx new file mode 100644 index 0000000..aa0352b --- /dev/null +++ b/src/radial_factory.cxx @@ -0,0 +1,41 @@ +#include +#include "radial_types.hpp" + +#include + +namespace IntegratorXX { + +RadialQuad radial_from_string(std::string name) { + std::transform(name.begin(), name.end(), name.begin(), ::toupper); + if(name == "BECKE") return RadialQuad::Becke; + if(name == "MURAKNOWLES") return RadialQuad::MuraKnowles; + if(name == "MK") return RadialQuad::MuraKnowles; + if(name == "MURRAYHANDYLAMING") return RadialQuad::MurrayHandyLaming; + if(name == "MHL") return RadialQuad::MurrayHandyLaming; + if(name == "TREUTLERAHLRICHS") return RadialQuad::TreutlerAhlrichs; + if(name == "TA") return RadialQuad::TreutlerAhlrichs; + + throw std::runtime_error("Unrecognized Radial Quadrature"); +} + +using radial_grid_ptr = RadialFactory::radial_grid_ptr; + +radial_grid_ptr RadialFactory::generate(RadialQuad rq, const RadialTraits& traits) { + + switch(rq) { + case RadialQuad::Becke: + return std::make_unique(traits); + case RadialQuad::MuraKnowles: + return std::make_unique(traits); + case RadialQuad::MurrayHandyLaming: + return std::make_unique(traits); + case RadialQuad::TreutlerAhlrichs: + return std::make_unique(traits); + default: + throw std::runtime_error("Unsupported Radial Quadrature"); + abort(); + } +} + + +} diff --git a/src/radial_types.hpp b/src/radial_types.hpp new file mode 100644 index 0000000..4edb1da --- /dev/null +++ b/src/radial_types.hpp @@ -0,0 +1,12 @@ +#pragma once +#include + +namespace IntegratorXX { + +using bk_type = Becke; +using mk_type = MuraKnowles; +using mhl_type = MurrayHandyLaming; +using ta_type = TreutlerAhlrichs; + +} + diff --git a/src/spherical_factory.cxx b/src/spherical_factory.cxx index 34007a4..1ecad8f 100644 --- a/src/spherical_factory.cxx +++ b/src/spherical_factory.cxx @@ -1,7 +1,8 @@ #include -#include #include +#include "radial_types.hpp" + #include namespace IntegratorXX { @@ -12,19 +13,6 @@ UnprunedSphericalGridSpecification::UnprunedSphericalGridSpecification( radial_quad(rq), radial_traits(traits.clone()), angular_quad(aq), angular_size(as) { } -RadialQuad radial_from_string(std::string name) { - std::transform(name.begin(), name.end(), name.begin(), ::toupper); - if(name == "BECKE") return RadialQuad::Becke; - if(name == "MURAKNOWLES") return RadialQuad::MuraKnowles; - if(name == "MK") return RadialQuad::MuraKnowles; - if(name == "MURRAYHANDYLAMING") return RadialQuad::MurrayHandyLaming; - if(name == "MHL") return RadialQuad::MurrayHandyLaming; - if(name == "TREUTLERAHLRICHS") return RadialQuad::TreutlerAhlrichs; - if(name == "TA") return RadialQuad::TreutlerAhlrichs; - - throw std::runtime_error("Unrecognized Radial Quadrature"); -} - AngularQuad angular_from_string(std::string name) { std::transform(name.begin(), name.end(), name.begin(), ::toupper); if(name == "AHRENSBEYLKIN") return AngularQuad::AhrensBeylkin; @@ -41,10 +29,6 @@ AngularQuad angular_from_string(std::string name) { using spherical_grid_ptr = SphericalGridFactory::spherical_grid_ptr; -using bk_type = Becke; -using mk_type = MuraKnowles; -using mhl_type = MurrayHandyLaming; -using ta_type = TreutlerAhlrichs; using ah_type = AhrensBeylkin; using de_type = Delley; diff --git a/test/spherical_generator.cxx b/test/spherical_generator.cxx index acf6290..893dd3e 100644 --- a/test/spherical_generator.cxx +++ b/test/spherical_generator.cxx @@ -6,6 +6,7 @@ #include #include #include +#include using bk_type = IntegratorXX::Becke; using mk_type = IntegratorXX::MuraKnowles; @@ -232,6 +233,32 @@ TEST_CASE( "Pruning Schemes", "[sph-gen]" ) { } +using radial_test_types = std::tuple< + bk_type, mk_type, mhl_type, ta_type +>; + +TEMPLATE_LIST_TEST_CASE("Radial Generator", "[sph-gen]", radial_test_types) { + using namespace IntegratorXX; + using radial_type = TestType; + + size_t npts = 10; + radial_type rq(npts, 1.0); + auto rad_spec = radial_from_type(); + auto rad_traits = make_radial_traits(rad_spec, npts, 1.0); + auto rad_grid = RadialFactory::generate(rad_spec, *rad_traits); + + REQUIRE(rad_grid->npts() == npts); + for(auto i = 0; i < npts; ++i) { + auto pt = rad_grid->points()[i]; + auto pt_ref = rq.points()[i]; + REQUIRE_THAT(pt, Catch::Matchers::WithinAbs(pt_ref,1e-15)); + + auto w = rad_grid->weights()[i]; + auto w_ref = rq.weights()[i]; + REQUIRE_THAT(w, Catch::Matchers::WithinAbs(w_ref,1e-15)); + } +} + using sph_test_types = std::tuple< std::tuple, std::tuple, From a71ec70fe57b8a3e39143e6f0c3ec468a8d5aa73 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Thu, 30 Nov 2023 11:24:15 -0800 Subject: [PATCH 6/6] Added S2Factory + UTs --- .../integratorxx/generators/s2_factory.hpp | 15 +++++++ src/CMakeLists.txt | 1 + src/s2_factory.cxx | 41 +++++++++++++++++++ src/s2_types.hpp | 11 +++++ src/spherical_factory.cxx | 20 +-------- test/spherical_generator.cxx | 31 ++++++++++++++ 6 files changed, 100 insertions(+), 19 deletions(-) create mode 100644 src/s2_factory.cxx create mode 100644 src/s2_types.hpp diff --git a/include/integratorxx/generators/s2_factory.hpp b/include/integratorxx/generators/s2_factory.hpp index 978dd04..1e550e4 100644 --- a/include/integratorxx/generators/s2_factory.hpp +++ b/include/integratorxx/generators/s2_factory.hpp @@ -1,6 +1,8 @@ #pragma once #include +#include + namespace IntegratorXX { /// High-level specification of angular quadratures @@ -23,4 +25,17 @@ AngularQuad angular_from_type() { AngularQuad angular_from_string(std::string name); +struct S2Factory { + + using s2_grid_ptr = std::shared_ptr< + QuadratureBase< + std::vector>, + std::vector + > + >; + + static s2_grid_ptr generate(AngularQuad aq, size_t npts); + +}; + } diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 52501da..ea07e5c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,6 +1,7 @@ add_library(integratorxx spherical_factory.cxx radial_factory.cxx + s2_factory.cxx ) # Target features diff --git a/src/s2_factory.cxx b/src/s2_factory.cxx new file mode 100644 index 0000000..e3a1d93 --- /dev/null +++ b/src/s2_factory.cxx @@ -0,0 +1,41 @@ +#include + +#include "s2_types.hpp" + +#include + +namespace IntegratorXX { + +AngularQuad angular_from_string(std::string name) { + std::transform(name.begin(), name.end(), name.begin(), ::toupper); + if(name == "AHRENSBEYLKIN") return AngularQuad::AhrensBeylkin; + if(name == "AB") return AngularQuad::AhrensBeylkin; + if(name == "DELLEY") return AngularQuad::Delley; + if(name == "LEBEDEVLAIKOV") return AngularQuad::LebedevLaikov; + if(name == "LEBEDEV") return AngularQuad::LebedevLaikov; + if(name == "LL") return AngularQuad::LebedevLaikov; + if(name == "WOMERSLEY") return AngularQuad::Womersley; + + throw std::runtime_error("Unrecognized Angular Quadrature"); +} + +using s2_grid_ptr = S2Factory::s2_grid_ptr; + +s2_grid_ptr S2Factory::generate(AngularQuad aq, size_t npts) { + + switch(aq) { + case AngularQuad::AhrensBeylkin: + return std::make_unique(npts); + case AngularQuad::Delley: + return std::make_unique(npts); + case AngularQuad::LebedevLaikov: + return std::make_unique(npts); + case AngularQuad::Womersley: + return std::make_unique(npts); + default: + throw std::runtime_error("Unsupported Angular Quadrature"); + abort(); + } +} + +} diff --git a/src/s2_types.hpp b/src/s2_types.hpp new file mode 100644 index 0000000..801d617 --- /dev/null +++ b/src/s2_types.hpp @@ -0,0 +1,11 @@ +#include + + +namespace IntegratorXX { + +using ah_type = AhrensBeylkin; +using de_type = Delley; +using ll_type = LebedevLaikov; +using wo_type = Womersley; + +} diff --git a/src/spherical_factory.cxx b/src/spherical_factory.cxx index 1ecad8f..513f0e9 100644 --- a/src/spherical_factory.cxx +++ b/src/spherical_factory.cxx @@ -1,9 +1,7 @@ #include -#include #include "radial_types.hpp" - -#include +#include "s2_types.hpp" namespace IntegratorXX { @@ -13,27 +11,11 @@ UnprunedSphericalGridSpecification::UnprunedSphericalGridSpecification( radial_quad(rq), radial_traits(traits.clone()), angular_quad(aq), angular_size(as) { } -AngularQuad angular_from_string(std::string name) { - std::transform(name.begin(), name.end(), name.begin(), ::toupper); - if(name == "AHRENSBEYLKIN") return AngularQuad::AhrensBeylkin; - if(name == "AB") return AngularQuad::AhrensBeylkin; - if(name == "DELLEY") return AngularQuad::Delley; - if(name == "LEBEDEVLAIKOV") return AngularQuad::LebedevLaikov; - if(name == "LEBEDEV") return AngularQuad::LebedevLaikov; - if(name == "LL") return AngularQuad::LebedevLaikov; - if(name == "WOMERSLEY") return AngularQuad::Womersley; - - throw std::runtime_error("Unrecognized Angular Quadrature"); -} using spherical_grid_ptr = SphericalGridFactory::spherical_grid_ptr; -using ah_type = AhrensBeylkin; -using de_type = Delley; -using ll_type = LebedevLaikov; -using wo_type = Womersley; /************************/ /**** Unpruned Grids ****/ diff --git a/test/spherical_generator.cxx b/test/spherical_generator.cxx index 893dd3e..8fbe482 100644 --- a/test/spherical_generator.cxx +++ b/test/spherical_generator.cxx @@ -259,6 +259,37 @@ TEMPLATE_LIST_TEST_CASE("Radial Generator", "[sph-gen]", radial_test_types) { } } +using s2_test_types = std::tuple< + ah_type, de_type, ll_type, wo_type +>; + +TEMPLATE_LIST_TEST_CASE("S2 Generator", "[sph-gen]", s2_test_types) { + using namespace IntegratorXX; + using angular_type = TestType; + using angular_traits = quadrature_traits; + + size_t npts = angular_traits::npts_by_algebraic_order( + angular_traits::next_algebraic_order(1)); // Smallest possible angular grid + + angular_type aq(npts); + + auto ang_spec = angular_from_type(); + auto ang_grid = S2Factory::generate(ang_spec, npts); + REQUIRE(ang_grid->npts() == npts); + for(auto i = 0; i < npts; ++i) { + auto pt = ang_grid->points()[i]; + auto pt_ref = aq.points()[i]; + REQUIRE_THAT(pt[0], Catch::Matchers::WithinAbs(pt_ref[0],1e-15)); + REQUIRE_THAT(pt[1], Catch::Matchers::WithinAbs(pt_ref[1],1e-15)); + REQUIRE_THAT(pt[2], Catch::Matchers::WithinAbs(pt_ref[2],1e-15)); + + auto w = ang_grid->weights()[i]; + auto w_ref = aq.weights()[i]; + REQUIRE_THAT(w, Catch::Matchers::WithinAbs(w_ref,1e-15)); + } + +} + using sph_test_types = std::tuple< std::tuple, std::tuple,