feat: Added Cartesian integrators to mathematica package

feat: Added 1D Cartesian integrator

Signed-off-by: Franz R. Sattler <[email protected]>

DiFfRG/include/DiFfRG/physics/integration/integrator_1D_cartesian_cpu.hh

@@ -0,0 +1,96 @@
+#pragma once
+
+// standard library
+#include <future>
+
+// external libraries
+#include <tbb/tbb.h>
+
+// DiFfRG
+#include <DiFfRG/physics/integration/quadrature_provider.hh>
+
+namespace DiFfRG
+{
+  template <typename NT, typename KERNEL> class Integrator1DCartesianTBB
+  {
+  public:
+    using ctype = typename get_type::ctype<NT>;
+    Integrator1DCartesianTBB(QuadratureProvider &quadrature_provider, const std::array<uint, 1> grid_size,
+                             const ctype x_extent = 0., const uint max_block_size = 0, const ctype qx_min = -M_PI,
+                             const ctype qx_max = M_PI)
+        : Integrator1DCartesianTBB(quadrature_provider, grid_size[0], x_extent, qx_min, qx_max)
+    {
+    }
+
+    Integrator1DCartesianTBB(QuadratureProvider &quadrature_provider, const std::array<uint, 1> grid_size,
+                             const JSONValue &json)
+        : Integrator1DCartesianTBB(quadrature_provider, grid_size[0], 0.,
+                                   json.get_double("/discretization/integration/qx_min", -M_PI),
+                                   json.get_double("/discretization/integration/qx_max", M_PI))
+    {
+    }
+
+    Integrator1DCartesianTBB(QuadratureProvider &quadrature_provider, const std::array<uint, 1> grid_size,
+                             const ctype x_extent, const JSONValue &json)
+        : Integrator1DCartesianTBB(quadrature_provider, grid_size[0], x_extent,
+                                   json.get_double("/discretization/integration/qx_min", -M_PI),
+                                   json.get_double("/discretization/integration/qx_max", M_PI))
+    {
+    }
+
+    Integrator1DCartesianTBB(QuadratureProvider &quadrature_provider, const uint grid_size, const ctype x_extent = 0.,
+                             const ctype qx_min = -M_PI, const ctype qx_max = M_PI)
+        : grid_size(grid_size), x_quadrature_p(quadrature_provider.get_points<ctype>(grid_size)),
+          x_quadrature_w(quadrature_provider.get_weights<ctype>(grid_size))
+    {
+      this->qx_min = qx_min;
+      this->qx_extent = qx_max - qx_min;
+    }
+
+    /**
+     * @brief Set the minimum value of the qx integration range.
+     */
+    void set_qx_min(const ctype qx_min)
+    {
+      this->qx_extent = this->qx_extent - qx_min + this->qx_min;
+      this->qx_min = qx_min;
+    }
+
+    /**
+     * @brief Set the maximum value of the qx integration range.
+     */
+    void set_qx_max(const ctype qx_max) { this->qx_extent = qx_max - qx_min; }
+
+    template <typename... T> NT get(const ctype k, const T &...t) const
+    {
+      constexpr int d = 1;
+      constexpr ctype int_element = powr<-d>(2 * (ctype)M_PI); // fourier factor
+
+      return KERNEL::constant(k, t...) + tbb::parallel_reduce(
+                                             tbb::blocked_range<uint>(0, grid_size), NT(0),
+                                             [&](const tbb::blocked_range<uint> &r, NT value) -> NT {
+                                               for (uint idx_x = r.begin(); idx_x != r.end(); ++idx_x) {
+                                                 const ctype q = qx_min + qx_extent * x_quadrature_p[idx_x];
+                                                 const ctype weight = qx_extent * x_quadrature_w[idx_x];
+                                                 value += int_element * weight * KERNEL::kernel(q, k, t...);
+                                               }
+                                               return value;
+                                             },
+                                             [&](NT x, NT y) -> NT { return x + y; });
+    }
+
+    template <typename... T> std::future<NT> request(const ctype k, const T &...t) const
+    {
+      return std::async(std::launch::deferred, [=, this]() { return get(k, t...); });
+    }
+
+  private:
+    const uint grid_size;
+
+    ctype qx_min = -M_PI;
+    ctype qx_extent = 2. * M_PI;
+
+    const std::vector<ctype> &x_quadrature_p;
+    const std::vector<ctype> &x_quadrature_w;
+  };
+} // namespace DiFfRG

DiFfRG/include/DiFfRG/physics/integration/integrator_1D_cartesian_gpu.hh

@@ -0,0 +1,181 @@
+#pragma once
+
+#ifdef __CUDACC__
+
+// standard library
+#include <future>
+
+// external libraries
+#include <rmm/cuda_stream_pool.hpp>
+#include <rmm/device_uvector.hpp>
+#include <rmm/mr/device/pool_memory_resource.hpp>
+#include <thrust/reduce.h>
+#include <thrust/transform_reduce.h>
+
+// DiFfRG
+#include <DiFfRG/common/cuda_prefix.hh>
+#include <DiFfRG/physics/integration/quadrature_provider.hh>
+
+namespace DiFfRG
+{
+  template <typename NT, typename KERNEL> class Integrator1DCartesianGPU
+  {
+  public:
+    using ctype = typename get_type::ctype<NT>;
+    template <typename... T> struct functor {
+    public:
+      functor(const ctype *x_quadrature_p, const ctype *x_quadrature_w, const ctype qx_min, const ctype qx_extent,
+              const ctype k, T... t)
+          : x_quadrature_p(x_quadrature_p), x_quadrature_w(x_quadrature_w), qx_min(qx_min), qx_extent(qx_extent), k(k),
+            t(t...)
+      {
+      }
+
+      __device__ NT operator()(const uint idx) const
+      {
+        constexpr int d = 1;
+        constexpr ctype int_element = powr<-d>(2 * (ctype)M_PI); // fourier factor
+
+        const ctype q = qx_min + qx_extent * x_quadrature_p[idx];
+        const ctype weight = qx_extent * x_quadrature_w[idx];
+
+        const NT res = std::apply([&](auto &&...args) { return KERNEL::kernel(q, k, args...); }, t);
+        return int_element * res * weight;
+      }
+
+    private:
+      const ctype *x_quadrature_p;
+      const ctype *x_quadrature_w;
+      const ctype qx_min, qx_extent;
+      const ctype k;
+      const std::tuple<T...> t;
+    };
+
+    Integrator1DCartesianGPU(QuadratureProvider &quadrature_provider, const std::array<uint, 1> grid_size,
+                             const ctype x_extent = 0., const uint max_block_size = 0, const ctype qx_min = -M_PI,
+                             const ctype qx_max = M_PI)
+        : Integrator1DCartesianGPU(quadrature_provider, grid_size[0], x_extent, qx_min, qx_max)
+    {
+    }
+
+    Integrator1DCartesianGPU(QuadratureProvider &quadrature_provider, const std::array<uint, 1> grid_size,
+                             const JSONValue &json)
+        : Integrator1DCartesianGPU(quadrature_provider, grid_size[0], 0.,
+                                   json.get_double("/discretization/integration/qx_min", -M_PI),
+                                   json.get_double("/discretization/integration/qx_max", M_PI))
+    {
+    }
+
+    Integrator1DCartesianGPU(QuadratureProvider &quadrature_provider, const std::array<uint, 1> grid_size,
+                             const ctype x_extent, const JSONValue &json)
+        : Integrator1DCartesianGPU(quadrature_provider, grid_size[0], x_extent,
+                                   json.get_double("/discretization/integration/qx_min", -M_PI),
+                                   json.get_double("/discretization/integration/qx_max", M_PI))
+    {
+    }
+
+    Integrator1DCartesianGPU(QuadratureProvider &quadrature_provider, const uint grid_size, const ctype x_extent = 0.,
+                             const ctype qx_min = -M_PI, const ctype qx_max = M_PI)
+        : grid_size(grid_size)
+    {
+      ptr_x_quadrature_p = quadrature_provider.get_device_points<ctype>(grid_size);
+      ptr_x_quadrature_w = quadrature_provider.get_device_weights<ctype>(grid_size);
+
+      this->qx_min = qx_min;
+      this->qx_extent = qx_max - qx_min;
+    }
+
+    Integrator1DCartesianGPU(const Integrator1DCartesianGPU &other)
+        : grid_size(other.grid_size), ptr_x_quadrature_p(other.ptr_x_quadrature_p),
+          ptr_x_quadrature_w(other.ptr_x_quadrature_w)
+    {
+      qx_min = other.qx_min;
+      qx_extent = other.qx_extent;
+    }
+
+    /**
+     * @brief Set the minimum value of the qx integration range.
+     */
+    void set_qx_min(const ctype qx_min)
+    {
+      this->qx_extent = this->qx_extent - qx_min + this->qx_min;
+      this->qx_min = qx_min;
+    }
+
+    /**
+     * @brief Set the maximum value of the qx integration range.
+     */
+    void set_qx_max(const ctype qx_max) { this->qx_extent = qx_max - qx_min; }
+
+    template <typename... T> NT get(const ctype k, const T &...t) const
+    {
+      return KERNEL::constant(k, t...) +
+             thrust::transform_reduce(thrust::cuda::par.on(rmm::cuda_stream_per_thread.value()),
+                                      thrust::make_counting_iterator<uint>(0),
+                                      thrust::make_counting_iterator<uint>(grid_size),
+                                      functor<T...>(ptr_x_quadrature_p, ptr_x_quadrature_w, qx_min, qx_extent, k, t...),
+                                      NT(0), thrust::plus<NT>());
+    }
+
+    template <typename... T> std::future<NT> request(const ctype k, const T &...t) const
+    {
+      return std::async(std::launch::deferred, [=, this]() { return get(k, t...); });
+    }
+
+  private:
+    const uint grid_size;
+
+    ctype qx_min = -M_PI;
+    ctype qx_extent = 2. * M_PI;
+
+    const ctype *ptr_x_quadrature_p;
+    const ctype *ptr_x_quadrature_w;
+  };
+} // namespace DiFfRG
+
+#else
+
+#ifdef USE_CUDA
+
+namespace DiFfRG
+{
+  template <typename NT, typename KERNEL> class Integrator1DCartesianGPU;
+}
+
+#else
+
+#include <DiFfRG/physics/integration/integrator_1D_cartesian_cpu.hh>
+
+namespace DiFfRG
+{
+  template <typename NT, typename KERNEL> class Integrator1DCartesianGPU : public IntegratorTBB<d, NT, KERNEL>
+  {
+  public:
+    using ctype = typename get_type::ctype<NT>;
+
+    Integrator1DCartesianGPU(QuadratureProvider &quadrature_provider, const std::array<uint, 1> grid_size,
+                             const ctype x_extent, const uint max_block_size = 256, const ctype qx_min = -M_PI,
+                             const ctype qx_max = M_PI)
+        : Integrator1DCartesianTBB<NT, KERNEL>(quadrature_provider, grid_size[0], x_extent, qx_min, qx_max)
+    {
+      (void)max_block_size;
+    }
+
+    Integrator1DCartesianGPU(QuadratureProvider &quadrature_provider, const uint grid_size, const ctype x_extent,
+                             const uint max_block_size = 256, const ctype qx_min = -M_PI, const ctype qx_max = M_PI)
+        : Integrator1DCartesianTBB<NT, KERNEL>(quadrature_provider, grid_size, x_extent, x_extent, qx_min, qx_max)
+    {
+      (void)max_block_size;
+    }
+
+    Integrator1DCartesianGPU(QuadratureProvider &quadrature_provider, const std::array<uint, 1> grid_size,
+                             const ctype x_extent, const JSONValue &json)
+        : Integrator1DCartesianTBB<NT, KERNEL>(quadrature_provider, grid_size, x_extent, json)
+    {
+    }
+  };
+} // namespace DiFfRG
+
+#endif
+
+#endif

DiFfRG/tests/physics/integration/CMakeLists.txt

@@ -41,6 +41,9 @@ setup_test(test_integrator_4D_2ang_cpu)
 add_executable(test_integrator_4D_qmc integrator_4D_qmc.cc)
 setup_test(test_integrator_4D_qmc)
 
+add_executable(test_integrator_1D_cartesian_cpu integrator_1D_cartesian_cpu.cc)
+setup_test(test_integrator_1D_cartesian_cpu)
+
 add_executable(test_integrator_2D_cartesian_cpu integrator_2D_cartesian_cpu.cc)
 setup_test(test_integrator_2D_cartesian_cpu)
 
@@ -88,6 +91,10 @@ if(USE_CUDA)
   add_executable(test_integrator_4D_qmc_gpu integrator_4D_qmc.cu)
   setup_test(test_integrator_4D_qmc_gpu)
 
+  add_executable(test_integrator_1D_cartesian_gpu
+                 integrator_1D_cartesian_gpu.cu)
+  setup_test(test_integrator_1D_cartesian_gpu)
+
   add_executable(test_integrator_2D_cartesian_gpu
                  integrator_2D_cartesian_gpu.cu)
   setup_test(test_integrator_2D_cartesian_gpu)

DiFfRG/tests/physics/integration/integrator_1D_cartesian_cpu.cc

@@ -0,0 +1,79 @@
+#define CATCH_CONFIG_MAIN
+#include <catch2/catch_all.hpp>
+
+#include <DiFfRG/common/math.hh>
+#include <DiFfRG/common/polynomials.hh>
+#include <DiFfRG/physics/integration/integrator_1D_cartesian_cpu.hh>
+#include <DiFfRG/physics/integration/quadrature_provider.hh>
+
+using namespace DiFfRG;
+
+//--------------------------------------------
+// Quadrature integration
+
+class PolyIntegrand
+{
+public:
+  static __forceinline__ __host__ __device__ auto kernel(const double q, const double /*k*/, const double /*c*/,
+                                                         const double x0, const double x1, const double x2,
+                                                         const double x3)
+  {
+    return (x0 + x1 * powr<1>(q) + x2 * powr<2>(q) + x3 * powr<3>(q));
+  }
+
+  static __forceinline__ __host__ __device__ auto constant(const double /*k*/, const double c, const double /*x0*/,
+                                                           const double /*x1*/, const double /*x2*/,
+                                                           const double /*x3*/)
+  {
+    return c;
+  }
+};
+
+TEST_CASE("Test 2d cartesian cpu momentum integrals", "[integration][quadrature integration]")
+{
+  const double qx_min = GENERATE(take(2, random(-2., -1.)));
+  const double qx_max = GENERATE(take(2, random(1., 2.)));
+
+  QuadratureProvider quadrature_provider;
+  Integrator1DCartesianTBB<double, PolyIntegrand> integrator(quadrature_provider, {{32}}, 0., 0, qx_min, qx_max);
+
+  SECTION("Volume integral")
+  {
+    const double k = GENERATE(take(1, random(0., 1.)));
+    const double reference_integral = (qx_max - qx_min) / powr<1>(2. * M_PI);
+
+    const double integral = integrator.request(k, 0., 1., 0., 0., 0.).get();
+
+    if (!is_close(reference_integral, integral, 1e-6)) {
+      std::cerr << "reference: " << reference_integral << "| integral: " << integral
+                << "| relative error: " << std::abs(reference_integral - integral) / std::abs(reference_integral)
+                << std::endl;
+    }
+    CHECK(is_close(reference_integral, integral, 1e-6));
+  }
+
+  SECTION("Random polynomials")
+  {
+    constexpr uint take_n = 2;
+    const auto x_poly = Polynomial({
+        GENERATE(take(take_n, random(-1., 1.))), // x0
+        GENERATE(take(take_n, random(-1., 1.))), // x1
+        GENERATE(take(take_n, random(-1., 1.))), // x2
+        GENERATE(take(take_n, random(-1., 1.))), // x3
+    });
+
+    const double k = GENERATE(take(take_n, random(0., 1.)));
+    const double constant = GENERATE(take(take_n, random(-1., 1.)));
+
+    const double reference_integral = constant + x_poly.integral(qx_min, qx_max) / powr<1>(2. * M_PI);
+
+    const double integral = integrator.get(k, constant, x_poly[0], x_poly[1], x_poly[2], x_poly[3]);
+
+    if (!is_close(reference_integral, integral, 1e-6)) {
+      std::cerr << "reference: " << reference_integral << "| integral: " << integral
+                << "| relative error: " << std::abs(reference_integral - integral) / std::abs(reference_integral)
+                << std::endl;
+    }
+    CHECK(is_close(reference_integral, integral, 1e-6));
+  }
+}