From d374101c561828bcac87d66eed497df15d3f03ed Mon Sep 17 00:00:00 2001 From: "Franz R. Sattler" Date: Sat, 18 Jan 2025 00:00:28 +0100 Subject: [PATCH] feat: Added Cartesian integrators to mathematica package feat: Added 1D Cartesian integrator Signed-off-by: Franz R. Sattler --- .../integrator_1D_cartesian_cpu.hh | 96 ++ .../integrator_1D_cartesian_gpu.hh | 181 +++ .../tests/physics/integration/CMakeLists.txt | 7 + .../integrator_1D_cartesian_cpu.cc | 79 + .../integrator_1D_cartesian_gpu.cu | 79 + Mathematica/DiFfRG/CodeTools.m | 179 +- Mathematica/DiFfRG/CodeTools.nb | 1443 +++++++++++------ 7 files changed, 1496 insertions(+), 568 deletions(-) create mode 100644 DiFfRG/include/DiFfRG/physics/integration/integrator_1D_cartesian_cpu.hh create mode 100644 DiFfRG/include/DiFfRG/physics/integration/integrator_1D_cartesian_gpu.hh create mode 100644 DiFfRG/tests/physics/integration/integrator_1D_cartesian_cpu.cc create mode 100644 DiFfRG/tests/physics/integration/integrator_1D_cartesian_gpu.cu diff --git a/DiFfRG/include/DiFfRG/physics/integration/integrator_1D_cartesian_cpu.hh b/DiFfRG/include/DiFfRG/physics/integration/integrator_1D_cartesian_cpu.hh new file mode 100644 index 0000000..c23c5b1 --- /dev/null +++ b/DiFfRG/include/DiFfRG/physics/integration/integrator_1D_cartesian_cpu.hh @@ -0,0 +1,96 @@ +#pragma once + +// standard library +#include + +// external libraries +#include + +// DiFfRG +#include + +namespace DiFfRG +{ + template class Integrator1DCartesianTBB + { + public: + using ctype = typename get_type::ctype; + Integrator1DCartesianTBB(QuadratureProvider &quadrature_provider, const std::array 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 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 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(grid_size)), + x_quadrature_w(quadrature_provider.get_weights(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 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(0, grid_size), NT(0), + [&](const tbb::blocked_range &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 std::future 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 &x_quadrature_p; + const std::vector &x_quadrature_w; + }; +} // namespace DiFfRG \ No newline at end of file diff --git a/DiFfRG/include/DiFfRG/physics/integration/integrator_1D_cartesian_gpu.hh b/DiFfRG/include/DiFfRG/physics/integration/integrator_1D_cartesian_gpu.hh new file mode 100644 index 0000000..2b0c385 --- /dev/null +++ b/DiFfRG/include/DiFfRG/physics/integration/integrator_1D_cartesian_gpu.hh @@ -0,0 +1,181 @@ +#pragma once + +#ifdef __CUDACC__ + +// standard library +#include + +// external libraries +#include +#include +#include +#include +#include + +// DiFfRG +#include +#include + +namespace DiFfRG +{ + template class Integrator1DCartesianGPU + { + public: + using ctype = typename get_type::ctype; + template 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; + }; + + Integrator1DCartesianGPU(QuadratureProvider &quadrature_provider, const std::array 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 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 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(grid_size); + ptr_x_quadrature_w = quadrature_provider.get_device_weights(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 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(0), + thrust::make_counting_iterator(grid_size), + functor(ptr_x_quadrature_p, ptr_x_quadrature_w, qx_min, qx_extent, k, t...), + NT(0), thrust::plus()); + } + + template std::future 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 class Integrator1DCartesianGPU; +} + +#else + +#include + +namespace DiFfRG +{ + template class Integrator1DCartesianGPU : public IntegratorTBB + { + public: + using ctype = typename get_type::ctype; + + Integrator1DCartesianGPU(QuadratureProvider &quadrature_provider, const std::array 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(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(quadrature_provider, grid_size, x_extent, x_extent, qx_min, qx_max) + { + (void)max_block_size; + } + + Integrator1DCartesianGPU(QuadratureProvider &quadrature_provider, const std::array grid_size, + const ctype x_extent, const JSONValue &json) + : Integrator1DCartesianTBB(quadrature_provider, grid_size, x_extent, json) + { + } + }; +} // namespace DiFfRG + +#endif + +#endif \ No newline at end of file diff --git a/DiFfRG/tests/physics/integration/CMakeLists.txt b/DiFfRG/tests/physics/integration/CMakeLists.txt index 632b6dc..e927712 100644 --- a/DiFfRG/tests/physics/integration/CMakeLists.txt +++ b/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) diff --git a/DiFfRG/tests/physics/integration/integrator_1D_cartesian_cpu.cc b/DiFfRG/tests/physics/integration/integrator_1D_cartesian_cpu.cc new file mode 100644 index 0000000..5824ad7 --- /dev/null +++ b/DiFfRG/tests/physics/integration/integrator_1D_cartesian_cpu.cc @@ -0,0 +1,79 @@ +#define CATCH_CONFIG_MAIN +#include + +#include +#include +#include +#include + +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 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)); + } +} \ No newline at end of file diff --git a/DiFfRG/tests/physics/integration/integrator_1D_cartesian_gpu.cu b/DiFfRG/tests/physics/integration/integrator_1D_cartesian_gpu.cu new file mode 100644 index 0000000..b4db323 --- /dev/null +++ b/DiFfRG/tests/physics/integration/integrator_1D_cartesian_gpu.cu @@ -0,0 +1,79 @@ +#define CATCH_CONFIG_MAIN +#include + +#include +#include +#include +#include + +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 gpu 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; + Integrator1DCartesianGPU 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)); + } +} \ No newline at end of file diff --git a/Mathematica/DiFfRG/CodeTools.m b/Mathematica/DiFfRG/CodeTools.m index a7ed6c0..f6cdac4 100644 --- a/Mathematica/DiFfRG/CodeTools.m +++ b/Mathematica/DiFfRG/CodeTools.m @@ -359,12 +359,23 @@ This Function creates an integrator that evaluates (constantFlow + \[Integral]in ] -GridSelector=<| -0->"grid_size_int", -1->"grid_sizes_angle_int", -2->"grid_sizes_3D_int", -3->"grid_sizes_4D_int" -|>; +GridSelector[kernel_]:=Module[{}, +Return@If[kernel["Type"]=="CartesianQuadrature", +Switch[kernel["d"], +1,"grid_size_int", +2,"grid_sizes_2D_cartesian_int", +3,"grid_sizes_3D_cartesian_int", +_,Print["CartesianQuadrature has wrong dimension!"];Abort[]; +], +Switch[kernel["Angles"], +0,"grid_size_int", +1,"grid_sizes_angle_int", +2,"grid_sizes_3D_int", +3,"grid_sizes_4D_int" +] +]; +]; + GridSelectorFiniteT=<| 0->"grid_sizes_int_fT", 1->"grid_sizes_angle_int", @@ -420,7 +431,7 @@ This Function creates an integrator that evaluates (constantFlow + \[Integral]in ]; -knownTypes={"Quadrature","QMC","Quadrature","Quadratureq0","Quadraturex0"}; +knownTypes={"Quadrature","QMC","Quadrature","Quadratureq0","Quadraturex0","CartesianQuadrature","CartesianQuadratureq0"}; TypeTest[type_]:=Module[{}, If[MemberQ[knownTypes,type],Return[True],Print["Unkown kernel Type: ",type, "\nKnown types: ",knownTypes];Return[False]] ]; @@ -786,6 +797,9 @@ This Function creates an integrator that evaluates (constantFlow + \[Integral]in const std::array grid_sizes_3D_int; const std::array grid_sizes_4D_int; + const std::array grid_sizes_2D_cartesian_int; + const std::array grid_sizes_3D_cartesian_int; + public: ::DiFfRG::QuadratureProvider quadrature_provider;"<> StringJoin[Map["\n ::DiFfRG::Flows::"<>#["Name"]<>"_integrator "<>#["Name"]<>"_integrator;"&,integratorList]]<>" @@ -798,8 +812,10 @@ This Function creates an integrator that evaluates (constantFlow + \[Integral]in grid_size_int{{x_quadrature_order}}, grid_sizes_angle_int{{x_quadrature_order, 2 * angle_quadrature_order}}, grid_sizes_3D_int{{x_quadrature_order, angle_quadrature_order, angle_quadrature_order}}, - grid_sizes_4D_int{{x_quadrature_order, angle_quadrature_order, angle_quadrature_order, angle_quadrature_order}}"<> -StringJoin[Map[",\n "<>#["Name"]<>"_integrator(quadrature_provider, "<>GridSelector[#["Angles"]]<>", x_extent, json)"&,integratorList]]<>" + grid_sizes_4D_int{{x_quadrature_order, angle_quadrature_order, angle_quadrature_order, angle_quadrature_order}}, + grid_sizes_2D_cartesian_int{{x_quadrature_order, x_quadrature_order}}, + grid_sizes_3D_cartesian_int{{x_quadrature_order, x_quadrature_order, x_quadrature_order}}"<> +StringJoin[Map[",\n "<>#["Name"]<>"_integrator(quadrature_provider, "<>GridSelector[#]<>", x_extent, json)"&,integratorList]]<>" { } "; @@ -824,9 +840,9 @@ This Function creates an integrator that evaluates (constantFlow + \[Integral]in If[name=="",Print["DiFfRG::CodeTools::MakeFlowClassFiniteT: Please provide a valid name to MakeFlowClass."];Abort[]]; -integratorList=Select[kernels,#["Type"]=="Quadrature"||#["Type"]=="Function1D"||#["Type"]=="Constant"&]; -integratorq0List=Select[kernels,#["Type"]=="Quadratureq0"||#["Type"]=="Function1Dq0"&]; -integratorx0List=Select[kernels,#["Type"]=="Quadraturex0"||#["Type"]=="Function1Dx0"&]; +integratorList=Select[kernels,#["Type"]=="Quadrature"&]; +integratorq0List=Select[kernels,#["Type"]=="Quadratureq0"&]; +integratorx0List=Select[kernels,#["Type"]=="Quadraturex0"&]; includeList=StringJoin[Map["#include \""<>#["Path"]<>"/"<>#["Name"]<>".hh\"\n"&,kernels]]; @@ -865,6 +881,9 @@ This Function creates an integrator that evaluates (constantFlow + \[Integral]in const std::array grid_sizes_angle_int; const std::array grid_sizes_4D_int; + const std::array grid_sizes_2D_cartesian_int; + const std::array grid_sizes_3D_cartesian_int; + public: QuadratureProvider quadrature_provider;"<> StringJoin[Map["\n Flows::"<>#["Name"]<>"_integrator "<>#["Name"]<>"_integrator;"&,Join[integratorList,integratorx0List,integratorq0List]]]<>" @@ -880,7 +899,9 @@ This Function creates an integrator that evaluates (constantFlow + \[Integral]in grid_size_int{{x_quadrature_order}}, grid_sizes_int_fT{{x_quadrature_order, x0_quadrature_order}}, grid_sizes_angle_int{{x_quadrature_order, 2 * angle_quadrature_order, x0_quadrature_order}}, - grid_sizes_4D_int{{x_quadrature_order, angle_quadrature_order, angle_quadrature_order, x0_quadrature_order}}"<> + grid_sizes_4D_int{{x_quadrature_order, angle_quadrature_order, angle_quadrature_order, x0_quadrature_order}}, + grid_sizes_2D_cartesian_int{{x_quadrature_order, x_quadrature_order}}, + grid_sizes_3D_cartesian_int{{x_quadrature_order, x_quadrature_order, x_quadrature_order}}"<> StringJoin[Map[",\n "<>#["Name"]<>"_integrator(quadrature_provider, grid_size_int, x_extent, json)"&,integratorList]]<> StringJoin[Map[",\n "<>#["Name"]<>"_integrator(quadrature_provider, "<>GridSelectorFiniteT[#["Angles"]]<>", x_extent, x0_extent, x0_summands, json)"&,integratorx0List]]<> StringJoin[Map[",\n "<>#["Name"]<>"_integrator(quadrature_provider, "<>GridSelectorFiniteT[#["Angles"]]<>", x_extent, q0_extent, x0_summands, json)"&,integratorq0List]]<>" @@ -914,6 +935,63 @@ This Function creates an integrator that evaluates (constantFlow + \[Integral]in ]; +MakeKernelCartesian1DMethod[flow_,parameterList_List,definitions_String,computeType_String]:=Module[{templateList,head, body}, +If[Not@IsValidParameterList[parameterList],Print["Invalid parameter List!"];Abort[]]; + +templateList=StringDrop[StringJoin[Table["typename T"<>ToString[i]<>", ",{i,1,Length[parameterList]}]],-2]; +templateList=" template <"<>templateList<>">"; +head=" static __forceinline__ __host__ __device__ auto + kernel(const "<>computeType<>" q, const "<>computeType<>" k, "; +head= +head<>StringDrop[StringJoin[Table["const T"<>ToString[i]<>ArgType[parameterList[[i]]["Type"]]<>" "<>ToString[parameterList[[i]]["Name"]]<>", ",{i,1,Length[parameterList]}]],-2]<>")"; +body=""; +If[definitions=!="", +body=body<>definitions<>"\n" +]; +body=body<>GetOptimizedKernelCode[flow,parameterList,computeType]; +body=IndentCode[body,2]; +templateList<>"\n"<>head<>"\n {\n"<>body<>"\n }" +]; + + +MakeKernelCartesian2DMethod[flow_,parameterList_List,definitions_String,computeType_String]:=Module[{templateList,head, body}, +If[Not@IsValidParameterList[parameterList],Print["Invalid parameter List!"];Abort[]]; + +templateList=StringDrop[StringJoin[Table["typename T"<>ToString[i]<>", ",{i,1,Length[parameterList]}]],-2]; +templateList=" template <"<>templateList<>">"; +head=" static __forceinline__ __host__ __device__ auto + kernel(const "<>computeType<>" qx, const "<>computeType<>" qy, const "<>computeType<>" k, "; +head= +head<>StringDrop[StringJoin[Table["const T"<>ToString[i]<>ArgType[parameterList[[i]]["Type"]]<>" "<>ToString[parameterList[[i]]["Name"]]<>", ",{i,1,Length[parameterList]}]],-2]<>")"; +body=""; +If[definitions=!="", +body=body<>definitions<>"\n" +]; +body=body<>GetOptimizedKernelCode[flow,parameterList,computeType]; +body=IndentCode[body,2]; +templateList<>"\n"<>head<>"\n {\n"<>body<>"\n }" +]; + + +MakeKernelCartesian3DMethod[flow_,parameterList_List,definitions_String,computeType_String]:=Module[{templateList,head, body}, +If[Not@IsValidParameterList[parameterList],Print["Invalid parameter List!"];Abort[]]; + +templateList=StringDrop[StringJoin[Table["typename T"<>ToString[i]<>", ",{i,1,Length[parameterList]}]],-2]; +templateList=" template <"<>templateList<>">"; +head=" static __forceinline__ __host__ __device__ auto + kernel(const "<>computeType<>" qx, const "<>computeType<>" qy, const "<>computeType<>" qz, const "<>computeType<>" k, "; +head= +head<>StringDrop[StringJoin[Table["const T"<>ToString[i]<>ArgType[parameterList[[i]]["Type"]]<>" "<>ToString[parameterList[[i]]["Name"]]<>", ",{i,1,Length[parameterList]}]],-2]<>")"; +body=""; +If[definitions=!="", +body=body<>definitions<>"\n" +]; +body=body<>GetOptimizedKernelCode[flow,parameterList,computeType]; +body=IndentCode[body,2]; +templateList<>"\n"<>head<>"\n {\n"<>body<>"\n }" +]; + + (* ::Input::Initialization:: *) MakeKernel0AngleMethod[flow_,parameterList_List,definitions_String,computeType_String]:=Module[{templateList,head, body}, If[Not@IsValidParameterList[parameterList],Print["Invalid parameter List!"];Abort[]]; @@ -1034,8 +1112,15 @@ This Function creates an integrator that evaluates (constantFlow + \[Integral]in If[kernel["Angles"]<0,Print["DiFfRG::CodeTools::MakeKernelClass: T = 0 kernel cannot have < 0 angles!"];Abort[]]; If[kernel["Type"]==="Constant", -body=MakeKernelConstMethod[integrandFlow,parameterList,integrandDefinitions,computeType] -, +body=MakeKernelConstMethod[integrandFlow,parameterList,integrandDefinitions,computeType], +If[kernel["Type"]==="CartesianQuadrature", +If[kernel["Angles"]!=0,Print["Cartesian integrator cannot have angles!"];Abort[];]; +Switch[kernel["d"], +1,body=MakeKernelCartesian1DMethod[integrandFlow,parameterList,integrandDefinitions,computeType], +2,body=MakeKernelCartesian2DMethod[integrandFlow,parameterList,integrandDefinitions,computeType], +3,body=MakeKernelCartesian3DMethod[integrandFlow,parameterList,integrandDefinitions,computeType], +_,Print["Cartesian integrator for dimension "<>ToString[kernel["d"]]<>" not implemented!"];Abort[]; +], If[kernel["Angles"]==0,body=MakeKernel0AngleMethod[integrandFlow,parameterList,integrandDefinitions,computeType]]; If[kernel["Angles"]==1,body=MakeKernel1AngleMethod[integrandFlow,parameterList,integrandDefinitions,computeType]]; If[kernel["Angles"]==2&&kernel["d"]==4,body=MakeKernel2AngleMethodAlt[integrandFlow,parameterList,integrandDefinitions,computeType]]; @@ -1043,6 +1128,7 @@ This Function creates an integrator that evaluates (constantFlow + \[Integral]in If[kernel["Angles"]==3,body=MakeKernel3AngleMethod[integrandFlow,parameterList,integrandDefinitions,computeType]]; If[kernel["Angles"]>3,Print["DiFfRG::CodeTools::MakeKernelClass: T = 0 kernel cannot have > 3 angles!"];Abort[]]; ]; +]; body=body<>"\n\n"<>MakeConstantMethod[constantFlow,parameterList,constantDefinitions,computeType]; @@ -1059,27 +1145,28 @@ This Function creates an integrator that evaluates (constantFlow + \[Integral]in kernelName=ToString[kernel["Name"]]<>"_kernel"; -suffix=If[kernel["Type"]==="QMC", -"QMC", -If[kernel["Type"]==="Constant", -If[kernel["Angles"]=!=0,Print["Constant integrator cannot have angles!"];Abort[]]; -"Constant", -If[kernel["Type"]==="Quadrature", -DeviceChoice[kernel["Device"]], -Print["Unknown integration type "<>kernel["Type"]<>"!"];Abort[] -] -] +suffix=Switch[kernel["Type"], +"QMC","QMC", +"Constant",If[kernel["Angles"]=!=0,Print["Constant integrator cannot have angles!"];Abort[]];"Constant", +"Quadrature",DeviceChoice[kernel["Device"]], +"CartesianQuadrature",If[kernel["Angles"]=!=0,Print["Cartesian integrator cannot have angles!"];Abort[]];"Cartesian"<>DeviceChoice[kernel["Device"]], +_,Print["Unknown integration type "<>kernel["Type"]<>"!"];Abort[] ]; -If[kernel["Angles"]==0, +If[kernel["Type"]=="CartesianQuadrature", +integrator="DiFfRG::Integrator"<>ToString[kernel["d"]]<>"D"<>suffix<>"<"<>computeType<>", "<>kernelName<>"<__REGULATOR__>>" ; +integratorAD="DiFfRG::Integrator"<>ToString[kernel["d"]]<>"D"<>suffix<>"<"<>"autodiff::real, "<>kernelName<>"<__REGULATOR__>>" , + +Switch[kernel["Angles"], +0, integrator="DiFfRG::Integrator"<>suffix<>"<"<>ToString[kernel["d"]]<>", "<>computeType<>", "<>kernelName<>"<__REGULATOR__>>" ; -integratorAD="DiFfRG::Integrator"<>suffix<>"<"<>ToString[kernel["d"]]<>", autodiff::real, "<>kernelName<>"<__REGULATOR__>>" -]; -If[kernel["Angles"]==1, +integratorAD="DiFfRG::Integrator"<>suffix<>"<"<>ToString[kernel["d"]]<>", autodiff::real, "<>kernelName<>"<__REGULATOR__>>" , + +1, integrator="DiFfRG::IntegratorAngle"<>suffix<>"<"<>ToString[kernel["d"]]<>", "<>computeType<>", "<>kernelName<>"<__REGULATOR__>>" ; -integratorAD="DiFfRG::IntegratorAngle"<>suffix<>"<"<>ToString[kernel["d"]]<>", autodiff::real, "<>kernelName<>"<__REGULATOR__>>" -]; -If[kernel["Angles"]==2, +integratorAD="DiFfRG::IntegratorAngle"<>suffix<>"<"<>ToString[kernel["d"]]<>", autodiff::real, "<>kernelName<>"<__REGULATOR__>>" , + +2, If[kernel["d"]!=3&&kernel["d"]!=4,Print["Inconsistent dimensions!"];Abort[]]; If[kernel["d"]==3, integrator="DiFfRG::Integrator3D"<>suffix<>"<"<>computeType<>", "<>kernelName<>"<__REGULATOR__>>" ; @@ -1087,21 +1174,22 @@ This Function creates an integrator that evaluates (constantFlow + \[Integral]in , integrator="DiFfRG::Integrator4D2Ang"<>suffix<>"<"<>computeType<>", "<>kernelName<>"<__REGULATOR__>>" ; integratorAD="DiFfRG::Integrator4D2Ang"<>suffix<>"kernelName<>"<__REGULATOR__>>" ; -]; -]; -If[kernel["Angles"]==3, +], + +3, If[kernel["d"]!=4,Print["Inconsistent dimensions!"];Abort[]]; integrator="DiFfRG::Integrator4D"<>suffix<>"<"<>computeType<>", "<>kernelName<>"<__REGULATOR__>>" ; integratorAD="DiFfRG::Integrator4D"<>suffix<>"kernelName<>"<__REGULATOR__>>" ]; +]; {integrator,integratorAD} -] +]; (* ::Input::Initialization:: *) MakeKernelIntegrator[kernel_Association,parameterList_List]:= -Module[{computeType,kernelName,className,integrator,integratorAD,argList,paramList,paramListAD,kernelFile,hhFile,ccFile,cuFile,hh,cc,cu}, +Module[{computeType,kernelName,className,integrator,integratorAD,gridSizes,argList,paramList,paramListAD,kernelFile,hhFile,ccFile,cuFile,hh,cc,cu}, If[Not@IsValidKernelSpec[kernel],Print["Invalid kernel!"];Abort[]]; If[Not@IsValidParameterList[parameterList],Print["Invalid parameter List!"];Abort[]]; @@ -1112,6 +1200,11 @@ This Function creates an integrator that evaluates (constantFlow + \[Integral]in {integrator,integratorAD}=MakeIntegratorTypes[kernel]; +gridSizes=If[kernel["Type"]=!="CartesianQuadrature", +kernel["Angles"]+1, +kernel["d"] +]; + argList=StringDrop[StringJoin[Table[""<>ToString[parameterList[[i]]["Name"]]<>", ",{i,1,Length[parameterList]}]],-2]; paramList=StringDrop[StringJoin[Table["const "<>CppType[computeType][parameterList[[i]]["Type"]]<>ArgType[parameterList[[i]]["Type"]]<>" "<>ToString[parameterList[[i]]["Name"]]<>", ",{i,1,Length[parameterList]}]],-2]; paramListAD=StringDrop[StringJoin[Table[ @@ -1142,7 +1235,7 @@ This Function creates an integrator that evaluates (constantFlow + \[Integral]in class "<>className<>" { public: - "<>className<>"(QuadratureProvider &quadrature_provider, std::arrayToString[kernel["Angles"]+1]<>"> grid_sizes, const "<>computeType<>" x_extent, const JSONValue& json); + "<>className<>"(QuadratureProvider &quadrature_provider, std::arrayToString[gridSizes]<>"> grid_sizes, const "<>computeType<>" x_extent, const JSONValue& json); "<>className<>"(const "<>className<>"& other); ~"<>className<>"(); @@ -1186,8 +1279,8 @@ NT get(T&&... t) ]<>" QuadratureProvider& quadrature_provider; - const std::arrayToString[kernel["Angles"]+1]<>"> grid_sizes; - std::arrayToString[kernel["Angles"]+1]<>"> jac_grid_sizes; + const std::arrayToString[gridSizes]<>"> grid_sizes; + std::arrayToString[gridSizes]<>"> jac_grid_sizes; const "<>computeType<>" x_extent; const "<>computeType<>" jacobian_quadrature_factor; const JSONValue json; @@ -1211,12 +1304,12 @@ NT get(T&&... t) { namespace Flows { - "<>className<>"::"<>className<>"(QuadratureProvider &quadrature_provider, std::arrayToString[kernel["Angles"]+1]<>"> grid_sizes, const "<>computeType<>" x_extent, const JSONValue& json) + "<>className<>"::"<>className<>"(QuadratureProvider &quadrature_provider, std::arrayToString[gridSizes]<>"> grid_sizes, const "<>computeType<>" x_extent, const JSONValue& json) : quadrature_provider(quadrature_provider), grid_sizes(grid_sizes), x_extent(x_extent), jacobian_quadrature_factor(json.get_double(\"/integration/jacobian_quadrature_factor\")), json(json) { integrator = std::make_unique<"<>integrator<>">(quadrature_provider, grid_sizes, x_extent, json);"<> If[kernel["AD"]," - for(uint i = 0; i < "<>ToString[kernel["Angles"]+1]<>"; ++i) + for(uint i = 0; i < "<>ToString[gridSizes]<>"; ++i) jac_grid_sizes[i] = uint(jacobian_quadrature_factor * grid_sizes[i]); integrator_AD = std::make_unique<"<>integratorAD<>">(quadrature_provider, jac_grid_sizes, x_extent, json);", "" @@ -1773,7 +1866,7 @@ NT get(T&&... t) ]; -QuadTypes={"QMC","Constant","Quadrature"} +QuadTypes={"QMC","Constant","Quadrature","CartesianQuadrature"} MakeKernel[kernel_Association, parameterList_List,integrandFlow_, constantFlow_:0., integrandDefinitions_String:"", constantDefinitions_String:""] := Module[{}, If[Not@IsValidKernelSpecList[kernels],Print["DiFfRG::CodeTools::MakeKernel: Invalid kernels list!"];Abort[]]; diff --git a/Mathematica/DiFfRG/CodeTools.nb b/Mathematica/DiFfRG/CodeTools.nb index f5ae0ea..3e322a7 100644 --- a/Mathematica/DiFfRG/CodeTools.nb +++ b/Mathematica/DiFfRG/CodeTools.nb @@ -10,10 +10,10 @@ NotebookFileLineBreakTest NotebookFileLineBreakTest NotebookDataPosition[ 158, 7] -NotebookDataLength[ 335539, 7585] -NotebookOptionsPosition[ 323993, 7373] -NotebookOutlinePosition[ 324781, 7399] -CellTagsIndexPosition[ 324738, 7396] +NotebookDataLength[ 354547, 7978] +NotebookOptionsPosition[ 342703, 7763] +NotebookOutlinePosition[ 343491, 7789] +CellTagsIndexPosition[ 343448, 7786] WindowFrame->Normal*) (* Beginning of Notebook Content *) @@ -1361,16 +1361,43 @@ FermionicCoordinates1DFiniteT>&\>\""}]}], "\n", "|>"}]}], Cell[BoxData[{ RowBox[{ - RowBox[{"GridSelector", "=", - RowBox[{"<|", "\[IndentingNewLine]", - RowBox[{ - RowBox[{"0", "->", "\"\\""}], ",", "\[IndentingNewLine]", - RowBox[{"1", "->", "\"\\""}], ",", - "\[IndentingNewLine]", - RowBox[{"2", "->", "\"\\""}], ",", - "\[IndentingNewLine]", - RowBox[{"3", "->", "\"\\""}]}], - "\[IndentingNewLine]", "|>"}]}], ";"}], "\[IndentingNewLine]", + RowBox[{ + RowBox[{ + RowBox[{"GridSelector", "[", "kernel_", "]"}], ":=", + RowBox[{"Module", "[", + RowBox[{ + RowBox[{"{", "}"}], ",", "\[IndentingNewLine]", + RowBox[{ + RowBox[{"Return", "@", + RowBox[{"If", "[", + RowBox[{ + RowBox[{ + RowBox[{"kernel", "[", "\"\\"", "]"}], "==", + "\"\\""}], ",", "\[IndentingNewLine]", + RowBox[{"Switch", "[", + RowBox[{ + RowBox[{"kernel", "[", "\"\\"", "]"}], ",", + "\[IndentingNewLine]", "1", ",", "\"\\"", ",", + "\[IndentingNewLine]", "2", ",", + "\"\\"", ",", "\[IndentingNewLine]", + "3", ",", "\"\\"", ",", + "\[IndentingNewLine]", "_", ",", + RowBox[{ + RowBox[{ + "Print", "[", "\"\\"", + "]"}], ";", + RowBox[{"Abort", "[", "]"}], ";"}]}], "\[IndentingNewLine]", + "]"}], ",", "\[IndentingNewLine]", + RowBox[{"Switch", "[", + RowBox[{ + RowBox[{"kernel", "[", "\"\\"", "]"}], ",", + "\[IndentingNewLine]", "0", ",", "\"\\"", ",", + "\[IndentingNewLine]", "1", ",", "\"\\"", + ",", "\[IndentingNewLine]", "2", ",", "\"\\"", + ",", "\[IndentingNewLine]", "3", ",", + "\"\\""}], "\[IndentingNewLine]", "]"}]}], + "\[IndentingNewLine]", "]"}]}], ";"}]}], "\[IndentingNewLine]", + "]"}]}], ";"}], "\[IndentingNewLine]"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"GridSelectorFiniteT", "=", @@ -1395,8 +1422,11 @@ Cell[BoxData[{ RowBox[{ RowBox[{"DeviceChoice", "[", "\"\\"", "]"}], "=", "\"\\""}], ";"}]}], "Input", - CellChangeTimes->{{3.946110633791431*^9, - 3.946110634591062*^9}},ExpressionUUID->"f45f2f7c-00b0-422c-b938-\ + CellChangeTimes->{{3.946110633791431*^9, 3.946110634591062*^9}, { + 3.9461412308609247`*^9, 3.946141353385625*^9}, {3.946141452429474*^9, + 3.946141506361285*^9}, {3.9461415370135098`*^9, 3.946141537993731*^9}, { + 3.946141591118038*^9, 3.9461416099933*^9}, {3.946141680869418*^9, + 3.9461417268180323`*^9}},ExpressionUUID->"f45f2f7c-00b0-422c-b938-\ 377d2da7520e"] }, Open ]], @@ -1448,7 +1478,7 @@ Cell[BoxData["\<\"Flow output directory: \ "During evaluation of \ In[195]:=",ExpressionUUID->"65d59a2b-50e4-43c5-9ad1-e9e34093a585"] }, Open ]] -}, Open ]], +}, Closed]], Cell[CellGroupData[{ @@ -1591,8 +1621,9 @@ Cell[BoxData[{ RowBox[{"{", RowBox[{ "\"\\"", ",", "\"\\"", ",", "\"\\"", ",", - "\"\\"", ",", "\"\\""}], "}"}]}], - ";"}], "\[IndentingNewLine]", + "\"\\"", ",", "\"\\"", ",", + "\"\\"", ",", "\"\\""}], + "}"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"TypeTest", "[", "type_", "]"}], ":=", @@ -1677,8 +1708,9 @@ Cell[BoxData[{ "\[IndentingNewLine]", "]"}]}], ";"}]}], "Input", CellChangeTimes->{{3.9387869406479607`*^9, 3.938786984053474*^9}, { 3.938787018053985*^9, 3.9387871811947002`*^9}, {3.938787211495647*^9, - 3.9387873319640017`*^9}, {3.9387888587572727`*^9, - 3.938788868407303*^9}},ExpressionUUID->"2ea8b7a0-b679-4c15-8bdc-\ + 3.9387873319640017`*^9}, {3.9387888587572727`*^9, 3.938788868407303*^9}, { + 3.9461399223375998`*^9, + 3.946139929329183*^9}},ExpressionUUID->"2ea8b7a0-b679-4c15-8bdc-\ 98e6abdfdf6f"], Cell[BoxData[{ @@ -1845,7 +1877,7 @@ Cell[BoxData[{ 3.9185701047583847`*^9, 3.9185701094640303`*^9}, 3.9387875479475117`*^9},ExpressionUUID->"abcfc1d3-2e59-4f0a-89d2-\ 0541b475bac8"] -}, Closed]], +}, Open ]], Cell[CellGroupData[{ @@ -3584,9 +3616,10 @@ Cell[BoxData[{ "\"\ grid_size_int;\n const std::array \ grid_sizes_angle_int;\n const std::array grid_sizes_3D_int;\n \ -const std::array grid_sizes_4D_int;\n\npublic:\n \ -::DiFfRG::QuadratureProvider quadrature_provider;\>\"", "<>", - "\[IndentingNewLine]", +const std::array grid_sizes_4D_int;\n\n const std::array \ +grid_sizes_2D_cartesian_int;\n const std::array \ +grid_sizes_3D_cartesian_int;\n\npublic:\n ::DiFfRG::QuadratureProvider \ +quadrature_provider;\>\"", "<>", "\[IndentingNewLine]", RowBox[{"StringJoin", "[", RowBox[{"Map", "[", RowBox[{ @@ -3610,8 +3643,10 @@ grid_sizes_angle_int{{x_quadrature_order, 2 * angle_quadrature_order}},\n \ grid_sizes_3D_int{{x_quadrature_order, \ angle_quadrature_order, angle_quadrature_order}},\n \ grid_sizes_4D_int{{x_quadrature_order, angle_quadrature_order, \ -angle_quadrature_order, angle_quadrature_order}}\>\"", "<>", - "\[IndentingNewLine]", +angle_quadrature_order, angle_quadrature_order}},\n \ +grid_sizes_2D_cartesian_int{{x_quadrature_order, x_quadrature_order}},\n \ + grid_sizes_3D_cartesian_int{{x_quadrature_order, \ +x_quadrature_order, x_quadrature_order}}\>\"", "<>", "\[IndentingNewLine]", RowBox[{"StringJoin", "[", RowBox[{"Map", "[", RowBox[{ @@ -3619,8 +3654,7 @@ angle_quadrature_order, angle_quadrature_order}}\>\"", "<>", RowBox[{"\"\<,\\n \>\"", "<>", RowBox[{"#", "[", "\"\\"", "]"}], "<>", "\"\<_integrator(quadrature_provider, \>\"", "<>", - RowBox[{"GridSelector", "[", - RowBox[{"#", "[", "\"\\"", "]"}], "]"}], "<>", + RowBox[{"GridSelector", "[", "#", "]"}], "<>", "\"\<, x_extent, json)\>\""}], "&"}], ",", "integratorList"}], "]"}], "]"}], "<>", "\"\<\n{\n}\n\>\""}]}], ";", "\[IndentingNewLine]", "\[IndentingNewLine]", @@ -3689,10 +3723,13 @@ angle_quadrature_order, angle_quadrature_order}}\>\"", "<>", 3.933304602209537*^9, 3.9333046277690268`*^9}, {3.93330466336119*^9, 3.933304789864666*^9}, {3.93330484332067*^9, 3.933304853136348*^9}, 3.933304964720351*^9, {3.933305064176303*^9, 3.933305066319943*^9}, { - 3.937407106337494*^9, - 3.9374071214025993`*^9}},ExpressionUUID->"343ea549-3ec9-4b86-9bc0-\ + 3.937407106337494*^9, 3.9374071214025993`*^9}, {3.946141076862392*^9, + 3.94614108756663*^9}, {3.946141119447261*^9, 3.946141172258122*^9}, { + 3.946141203314022*^9, 3.946141207650319*^9}, 3.946141369226582*^9, { + 3.946141833634643*^9, + 3.946141834758081*^9}},ExpressionUUID->"343ea549-3ec9-4b86-9bc0-\ 73c18f46ff9e"] -}, Closed]], +}, Open ]], Cell[CellGroupData[{ @@ -3776,39 +3813,24 @@ valid name to MakeFlowClass.\>\"", "]"}], ";", RowBox[{"kernels", ",", RowBox[{ RowBox[{ - RowBox[{ - RowBox[{"#", "[", "\"\\"", "]"}], "==", - "\"\\""}], "||", - RowBox[{ - RowBox[{"#", "[", "\"\\"", "]"}], "==", - "\"\\""}], "||", - RowBox[{ - RowBox[{"#", "[", "\"\\"", "]"}], "==", - "\"\\""}]}], "&"}]}], "]"}]}], ";", + RowBox[{"#", "[", "\"\\"", "]"}], "==", + "\"\\""}], "&"}]}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"integratorq0List", "=", RowBox[{"Select", "[", RowBox[{"kernels", ",", RowBox[{ RowBox[{ - RowBox[{ - RowBox[{"#", "[", "\"\\"", "]"}], "==", - "\"\\""}], "||", - RowBox[{ - RowBox[{"#", "[", "\"\\"", "]"}], "==", - "\"\\""}]}], "&"}]}], "]"}]}], ";", + RowBox[{"#", "[", "\"\\"", "]"}], "==", + "\"\\""}], "&"}]}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"integratorx0List", "=", RowBox[{"Select", "[", RowBox[{"kernels", ",", RowBox[{ RowBox[{ - RowBox[{ - RowBox[{"#", "[", "\"\\"", "]"}], "==", - "\"\\""}], "||", - RowBox[{ - RowBox[{"#", "[", "\"\\"", "]"}], "==", - "\"\\""}]}], "&"}]}], "]"}]}], ";", + RowBox[{"#", "[", "\"\\"", "]"}], "==", + "\"\\""}], "&"}]}], "]"}]}], ";", "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"includeList", "=", RowBox[{"StringJoin", "[", @@ -3849,7 +3871,9 @@ valid name to MakeFlowClass.\>\"", "]"}], ";", "\"\ grid_size_int;\n const std::array \ grid_sizes_int_fT;\n const std::array grid_sizes_angle_int;\n \ -const std::array grid_sizes_4D_int;\n\npublic:\n QuadratureProvider \ +const std::array grid_sizes_4D_int;\n\n const std::array \ +grid_sizes_2D_cartesian_int;\n const std::array \ +grid_sizes_3D_cartesian_int;\n\npublic:\n QuadratureProvider \ quadrature_provider;\>\"", "<>", "\[IndentingNewLine]", RowBox[{"StringJoin", "[", RowBox[{"Map", "[", @@ -3881,15 +3905,16 @@ grid_sizes_int_fT{{x_quadrature_order, x0_quadrature_order}},\n \ grid_sizes_angle_int{{x_quadrature_order, 2 * angle_quadrature_order, \ x0_quadrature_order}},\n \ grid_sizes_4D_int{{x_quadrature_order, angle_quadrature_order, \ -angle_quadrature_order, x0_quadrature_order}}\>\"", "<>", - "\[IndentingNewLine]", +angle_quadrature_order, x0_quadrature_order}},\n \ +grid_sizes_2D_cartesian_int{{x_quadrature_order, x_quadrature_order}},\n \ + grid_sizes_3D_cartesian_int{{x_quadrature_order, \ +x_quadrature_order, x_quadrature_order}}\>\"", "<>", "\[IndentingNewLine]", RowBox[{"StringJoin", "[", RowBox[{"Map", "[", RowBox[{ RowBox[{ RowBox[{"\"\<,\\n \>\"", "<>", RowBox[{"#", "[", "\"\\"", "]"}], "<>", - "\"\<_integrator(quadrature_provider, grid_size_int, x_extent, \ json)\>\""}], "&"}], ",", "integratorList"}], "]"}], "]"}], "<>", "\[IndentingNewLine]", @@ -3999,11 +4024,13 @@ json)\>\""}], "&"}], ",", "integratorList"}], "]"}], "]"}], "<>", 3.933304593801551*^9}, {3.933304772097981*^9, 3.933304857320385*^9}, 3.93330497050459*^9, 3.933305016784762*^9, {3.933305056513434*^9, 3.933305081688542*^9}, {3.9366886751026497`*^9, 3.936688679955563*^9}, { - 3.9374071275789423`*^9, - 3.9374071419234467`*^9}},ExpressionUUID->"1ef8a59a-d834-40ef-86a8-\ + 3.9374071275789423`*^9, 3.9374071419234467`*^9}, {3.9461384391740026`*^9, + 3.946138445169738*^9}, {3.946141167871611*^9, 3.946141181999564*^9}, { + 3.946141827090726*^9, + 3.946141827750345*^9}},ExpressionUUID->"1ef8a59a-d834-40ef-86a8-\ 99dc4a1d5a67"] -}, Closed]] -}, Open ]], +}, Open ]] +}, Closed]], Cell[CellGroupData[{ @@ -4121,7 +4148,7 @@ Cell[BoxData[ Cell[BoxData[ RowBox[{ RowBox[{ - RowBox[{"MakeKernel0AngleMethod", "[", + RowBox[{"MakeKernelCartesian1DMethod", "[", RowBox[{ "flow_", ",", "parameterList_List", ",", "definitions_String", ",", "computeType_String"}], "]"}], ":=", @@ -4207,32 +4234,14 @@ Cell[BoxData[ "templateList", "<>", "\"\<\\n\>\"", "<>", "head", "<>", "\"\<\\n {\\n\>\"", "<>", "body", "<>", "\"\<\\n }\>\""}]}]}], "\n", "]"}]}], ";"}]], "Input", - InitializationCell->True, - CellChangeTimes->{{3.9165597608463287`*^9, 3.916559972877629*^9}, { - 3.91656021195391*^9, 3.916560229952181*^9}, {3.917167721293209*^9, - 3.917167743815441*^9}, {3.918381371576629*^9, 3.918381403712631*^9}, { - 3.918381456584001*^9, 3.918381462016645*^9}, {3.91856961273317*^9, - 3.918569639974669*^9}, {3.918622775852169*^9, 3.918622806300225*^9}, { - 3.919102418098164*^9, 3.919102663047504*^9}, {3.919927050145592*^9, - 3.919927155870949*^9}, {3.919927198619097*^9, 3.919927212232974*^9}, { - 3.9199273964262667`*^9, 3.919927401569687*^9}, {3.919995491191346*^9, - 3.919995518261923*^9}, {3.919995771654359*^9, 3.919995785682316*^9}, { - 3.922147469460735*^9, 3.922147487609593*^9}, {3.925272931670052*^9, - 3.925272959908493*^9}, 3.925273143560604*^9, {3.925708447994582*^9, - 3.925708561339304*^9}, {3.9257085956106243`*^9, 3.925708672230053*^9}, { - 3.9257093172950478`*^9, 3.925709336609571*^9}, {3.925709545845115*^9, - 3.925709588340704*^9}, {3.927459263128355*^9, 3.927459273627915*^9}, { - 3.927459341017798*^9, 3.927459355297276*^9}, {3.9367184101902227`*^9, - 3.9367184132994833`*^9}, {3.936718450395412*^9, 3.93671852914922*^9}, { - 3.936726723136002*^9, 3.9367267298239594`*^9}, {3.938783633680307*^9, - 3.938783666400959*^9}, {3.938783704349848*^9, 3.938783705657412*^9}, - 3.938785926852797*^9},ExpressionUUID->"edd3ec15-a535-47c7-8995-\ -12f62cc99a8a"], + CellChangeTimes->{{3.946140642393536*^9, + 3.946140645933696*^9}},ExpressionUUID->"6056740c-a58d-4f7a-86cc-\ +8745d1fcdd03"], Cell[BoxData[ RowBox[{ RowBox[{ - RowBox[{"MakeKernel1AngleMethod", "[", + RowBox[{"MakeKernelCartesian2DMethod", "[", RowBox[{ "flow_", ",", "parameterList_List", ",", "definitions_String", ",", "computeType_String"}], "]"}], ":=", @@ -4270,11 +4279,11 @@ Cell[BoxData[ ";", "\[IndentingNewLine]", RowBox[{"head", "=", RowBox[{ - "\"\< static __forceinline__ __host__ __device__ auto kernel(const \>\ -\"", "<>", "computeType", "<>", "\"\< q, const \>\"", "<>", "computeType", - "<>", "\"\< cos1, const \>\"", "<>", "computeType", "<>", + "\"\< static __forceinline__ __host__ __device__ auto\n kernel(const \ +\>\"", "<>", "computeType", "<>", "\"\< qx, const \>\"", "<>", "computeType", + "<>", "\"\< qy, const \>\"", "<>", "computeType", "<>", "\"\< k, \>\""}]}], ";", "\[IndentingNewLine]", - RowBox[{"head", "=", + RowBox[{"head", "=", "\[IndentingNewLine]", RowBox[{"head", "<>", RowBox[{"StringDrop", "[", RowBox[{ @@ -4319,14 +4328,14 @@ Cell[BoxData[ "templateList", "<>", "\"\<\\n\>\"", "<>", "head", "<>", "\"\<\\n {\\n\>\"", "<>", "body", "<>", "\"\<\\n }\>\""}]}]}], "\n", "]"}]}], ";"}]], "Input", - CellChangeTimes->{{3.938783719587192*^9, 3.9387837206671677`*^9}, - 3.9387859291810207`*^9},ExpressionUUID->"814b8edf-c0ef-4ddf-a457-\ -f04ed08a9fd7"], + CellChangeTimes->{{3.946140642393536*^9, + 3.946140673961491*^9}},ExpressionUUID->"4a990ef2-746c-4f20-981a-\ +4357ab6c2e63"], Cell[BoxData[ RowBox[{ RowBox[{ - RowBox[{"MakeKernel2AngleMethod", "[", + RowBox[{"MakeKernelCartesian3DMethod", "[", RowBox[{ "flow_", ",", "parameterList_List", ",", "definitions_String", ",", "computeType_String"}], "]"}], ":=", @@ -4364,12 +4373,12 @@ Cell[BoxData[ ";", "\[IndentingNewLine]", RowBox[{"head", "=", RowBox[{ - "\"\< static __forceinline__ __host__ __device__ auto kernel(const \>\ -\"", "<>", "computeType", "<>", "\"\< q, const \>\"", "<>", "computeType", - "<>", "\"\< cos1, const \>\"", "<>", "computeType", "<>", - "\"\< phi, const \>\"", "<>", "computeType", "<>", "\"\< k, \>\""}]}], - ";", "\[IndentingNewLine]", - RowBox[{"head", "=", + "\"\< static __forceinline__ __host__ __device__ auto\n kernel(const \ +\>\"", "<>", "computeType", "<>", "\"\< qx, const \>\"", "<>", "computeType", + "<>", "\"\< qy, const \>\"", "<>", "computeType", "<>", + "\"\< qz, const \>\"", "<>", "computeType", "<>", "\"\< k, \>\""}]}], + ";", "\[IndentingNewLine]", + RowBox[{"head", "=", "\[IndentingNewLine]", RowBox[{"head", "<>", RowBox[{"StringDrop", "[", RowBox[{ @@ -4387,7 +4396,7 @@ Cell[BoxData[ RowBox[{ RowBox[{"parameterList", "[", RowBox[{"[", "i", "]"}], "]"}], "[", "\"\\"", "]"}], - "]"}], "<>", "\"\<, \>\""}], ",", "\[IndentingNewLine]", + "]"}], "<>", "\"\<, \>\""}], ",", RowBox[{"{", RowBox[{"i", ",", "1", ",", RowBox[{"Length", "[", "parameterList", "]"}]}], "}"}]}], @@ -4414,14 +4423,14 @@ Cell[BoxData[ "templateList", "<>", "\"\<\\n\>\"", "<>", "head", "<>", "\"\<\\n {\\n\>\"", "<>", "body", "<>", "\"\<\\n }\>\""}]}]}], "\n", "]"}]}], ";"}]], "Input", - CellChangeTimes->{{3.938783722033516*^9, 3.938783722662479*^9}, - 3.9387859316149683`*^9},ExpressionUUID->"e2db6879-74a3-42e6-be87-\ -6f113e5ec7af"], + CellChangeTimes->{{3.946140642393536*^9, + 3.946140697241465*^9}},ExpressionUUID->"9b93d1c4-999f-46c4-89c4-\ +72268839be50"], Cell[BoxData[ RowBox[{ RowBox[{ - RowBox[{"MakeKernel2AngleMethodAlt", "[", + RowBox[{"MakeKernel0AngleMethod", "[", RowBox[{ "flow_", ",", "parameterList_List", ",", "definitions_String", ",", "computeType_String"}], "]"}], ":=", @@ -4459,12 +4468,10 @@ Cell[BoxData[ ";", "\[IndentingNewLine]", RowBox[{"head", "=", RowBox[{ - "\"\< static __forceinline__ __host__ __device__ auto kernel(const \>\ -\"", "<>", "computeType", "<>", "\"\< q, const \>\"", "<>", "computeType", - "<>", "\"\< cos1, const \>\"", "<>", "computeType", "<>", - "\"\< cos2, const \>\"", "<>", "computeType", "<>", - "\"\< k, \>\""}]}], ";", "\[IndentingNewLine]", - RowBox[{"head", "=", + "\"\< static __forceinline__ __host__ __device__ auto\n kernel(const \ +\>\"", "<>", "computeType", "<>", "\"\< q, const \>\"", "<>", "computeType", + "<>", "\"\< k, \>\""}]}], ";", "\[IndentingNewLine]", + RowBox[{"head", "=", "\[IndentingNewLine]", RowBox[{"head", "<>", RowBox[{"StringDrop", "[", RowBox[{ @@ -4482,7 +4489,7 @@ Cell[BoxData[ RowBox[{ RowBox[{"parameterList", "[", RowBox[{"[", "i", "]"}], "]"}], "[", "\"\\"", "]"}], - "]"}], "<>", "\"\<, \>\""}], ",", "\[IndentingNewLine]", + "]"}], "<>", "\"\<, \>\""}], ",", RowBox[{"{", RowBox[{"i", ",", "1", ",", RowBox[{"Length", "[", "parameterList", "]"}]}], "}"}]}], @@ -4509,14 +4516,32 @@ Cell[BoxData[ "templateList", "<>", "\"\<\\n\>\"", "<>", "head", "<>", "\"\<\\n {\\n\>\"", "<>", "body", "<>", "\"\<\\n }\>\""}]}]}], "\n", "]"}]}], ";"}]], "Input", - CellChangeTimes->{{3.938783722934884*^9, 3.938783723692741*^9}, - 3.93878593381638*^9},ExpressionUUID->"f22375e1-43dd-4190-a425-\ -5a521600f912"], + InitializationCell->True, + CellChangeTimes->{{3.9165597608463287`*^9, 3.916559972877629*^9}, { + 3.91656021195391*^9, 3.916560229952181*^9}, {3.917167721293209*^9, + 3.917167743815441*^9}, {3.918381371576629*^9, 3.918381403712631*^9}, { + 3.918381456584001*^9, 3.918381462016645*^9}, {3.91856961273317*^9, + 3.918569639974669*^9}, {3.918622775852169*^9, 3.918622806300225*^9}, { + 3.919102418098164*^9, 3.919102663047504*^9}, {3.919927050145592*^9, + 3.919927155870949*^9}, {3.919927198619097*^9, 3.919927212232974*^9}, { + 3.9199273964262667`*^9, 3.919927401569687*^9}, {3.919995491191346*^9, + 3.919995518261923*^9}, {3.919995771654359*^9, 3.919995785682316*^9}, { + 3.922147469460735*^9, 3.922147487609593*^9}, {3.925272931670052*^9, + 3.925272959908493*^9}, 3.925273143560604*^9, {3.925708447994582*^9, + 3.925708561339304*^9}, {3.9257085956106243`*^9, 3.925708672230053*^9}, { + 3.9257093172950478`*^9, 3.925709336609571*^9}, {3.925709545845115*^9, + 3.925709588340704*^9}, {3.927459263128355*^9, 3.927459273627915*^9}, { + 3.927459341017798*^9, 3.927459355297276*^9}, {3.9367184101902227`*^9, + 3.9367184132994833`*^9}, {3.936718450395412*^9, 3.93671852914922*^9}, { + 3.936726723136002*^9, 3.9367267298239594`*^9}, {3.938783633680307*^9, + 3.938783666400959*^9}, {3.938783704349848*^9, 3.938783705657412*^9}, + 3.938785926852797*^9},ExpressionUUID->"edd3ec15-a535-47c7-8995-\ +12f62cc99a8a"], Cell[BoxData[ RowBox[{ RowBox[{ - RowBox[{"MakeKernel3AngleMethod", "[", + RowBox[{"MakeKernel1AngleMethod", "[", RowBox[{ "flow_", ",", "parameterList_List", ",", "definitions_String", ",", "computeType_String"}], "]"}], ":=", @@ -4557,9 +4582,7 @@ Cell[BoxData[ "\"\< static __forceinline__ __host__ __device__ auto kernel(const \>\ \"", "<>", "computeType", "<>", "\"\< q, const \>\"", "<>", "computeType", "<>", "\"\< cos1, const \>\"", "<>", "computeType", "<>", - "\"\< cos2, const \>\"", "<>", "computeType", "<>", - "\"\< phi, const \>\"", "<>", "computeType", "<>", "\"\< k, \>\""}]}], - ";", "\[IndentingNewLine]", + "\"\< k, \>\""}]}], ";", "\[IndentingNewLine]", RowBox[{"head", "=", RowBox[{"head", "<>", RowBox[{"StringDrop", "[", @@ -4605,14 +4628,14 @@ Cell[BoxData[ "templateList", "<>", "\"\<\\n\>\"", "<>", "head", "<>", "\"\<\\n {\\n\>\"", "<>", "body", "<>", "\"\<\\n }\>\""}]}]}], "\n", "]"}]}], ";"}]], "Input", - CellChangeTimes->{{3.938783723889003*^9, 3.938783724743881*^9}, - 3.938785935812346*^9},ExpressionUUID->"844d1ae0-39fa-48c5-a378-\ -1a543ed454ba"], + CellChangeTimes->{{3.938783719587192*^9, 3.9387837206671677`*^9}, + 3.9387859291810207`*^9},ExpressionUUID->"814b8edf-c0ef-4ddf-a457-\ +f04ed08a9fd7"], Cell[BoxData[ RowBox[{ RowBox[{ - RowBox[{"MakeConstantMethod", "[", + RowBox[{"MakeKernel2AngleMethod", "[", RowBox[{ "flow_", ",", "parameterList_List", ",", "definitions_String", ",", "computeType_String"}], "]"}], ":=", @@ -4650,9 +4673,11 @@ Cell[BoxData[ ";", "\[IndentingNewLine]", RowBox[{"head", "=", RowBox[{ - "\"\< static __forceinline__ __host__ __device__ auto constant(const \ -\>\"", "<>", "computeType", "<>", "\"\< k, \>\""}]}], ";", - "\[IndentingNewLine]", + "\"\< static __forceinline__ __host__ __device__ auto kernel(const \>\ +\"", "<>", "computeType", "<>", "\"\< q, const \>\"", "<>", "computeType", + "<>", "\"\< cos1, const \>\"", "<>", "computeType", "<>", + "\"\< phi, const \>\"", "<>", "computeType", "<>", "\"\< k, \>\""}]}], + ";", "\[IndentingNewLine]", RowBox[{"head", "=", RowBox[{"head", "<>", RowBox[{"StringDrop", "[", @@ -4671,7 +4696,7 @@ Cell[BoxData[ RowBox[{ RowBox[{"parameterList", "[", RowBox[{"[", "i", "]"}], "]"}], "[", "\"\\"", "]"}], - "]"}], "<>", "\"\<, \>\""}], ",", + "]"}], "<>", "\"\<, \>\""}], ",", "\[IndentingNewLine]", RowBox[{"{", RowBox[{"i", ",", "1", ",", RowBox[{"Length", "[", "parameterList", "]"}]}], "}"}]}], @@ -4698,154 +4723,488 @@ Cell[BoxData[ "templateList", "<>", "\"\<\\n\>\"", "<>", "head", "<>", "\"\<\\n {\\n\>\"", "<>", "body", "<>", "\"\<\\n }\>\""}]}]}], "\n", "]"}]}], ";"}]], "Input", - CellChangeTimes->{{3.9387837249584923`*^9, 3.938783729652721*^9}, - 3.938785937577661*^9},ExpressionUUID->"c8c30020-0969-4ec0-90cc-\ -af86941facd1"] -}, Closed]], - -Cell[CellGroupData[{ - -Cell["Kernel class creation", "Subsection", - CellChangeTimes->{{3.916559983910474*^9, 3.916559989385697*^9}, { - 3.917169856535161*^9, - 3.917169861537467*^9}},ExpressionUUID->"2a678e8d-3220-4dc3-80ad-\ -773d02937bfb"], + CellChangeTimes->{{3.938783722033516*^9, 3.938783722662479*^9}, + 3.9387859316149683`*^9},ExpressionUUID->"e2db6879-74a3-42e6-be87-\ +6f113e5ec7af"], Cell[BoxData[ RowBox[{ RowBox[{ - RowBox[{"MakeKernelClass", "[", + RowBox[{"MakeKernel2AngleMethodAlt", "[", RowBox[{ - "kernel_Association", ",", "parameterList_List", ",", "integrandFlow_", - ",", "constantFlow_", ",", "integrandDefinitions_String", ",", - "constantDefinitions_String"}], "]"}], ":=", + "flow_", ",", "parameterList_List", ",", "definitions_String", ",", + "computeType_String"}], "]"}], ":=", RowBox[{"Module", "[", RowBox[{ RowBox[{"{", - RowBox[{"computeType", ",", "head", ",", " ", "body", ",", "tail"}], - "}"}], ",", "\[IndentingNewLine]", + RowBox[{"templateList", ",", "head", ",", " ", "body"}], "}"}], ",", + "\[IndentingNewLine]", RowBox[{ - RowBox[{"If", "[", - RowBox[{ - RowBox[{"Not", "@", - RowBox[{"IsValidKernelSpec", "[", "kernel", "]"}]}], ",", - RowBox[{ - RowBox[{ - "Print", "[", - "\"\\"", - "]"}], ";", - RowBox[{"Abort", "[", "]"}]}]}], "]"}], ";", "\[IndentingNewLine]", RowBox[{"If", "[", RowBox[{ RowBox[{"Not", "@", RowBox[{"IsValidParameterList", "[", "parameterList", "]"}]}], ",", RowBox[{ - RowBox[{ - "Print", "[", - "\"\\"", "]"}], ";", + RowBox[{"Print", "[", "\"\\"", "]"}], ";", + RowBox[{"Abort", "[", "]"}]}]}], "]"}], ";", "\[IndentingNewLine]", "\[IndentingNewLine]", - RowBox[{"computeType", "=", - RowBox[{"kernel", "[", "\"\\"", "]"}]}], ";", - "\[IndentingNewLine]", "\[IndentingNewLine]", + RowBox[{"templateList", "=", + RowBox[{"StringDrop", "[", + RowBox[{ + RowBox[{"StringJoin", "[", + RowBox[{"Table", "[", + RowBox[{ + RowBox[{"\"\\"", "<>", + RowBox[{"ToString", "[", "i", "]"}], "<>", "\"\<, \>\""}], ",", + RowBox[{"{", + RowBox[{"i", ",", "1", ",", + RowBox[{"Length", "[", "parameterList", "]"}]}], "}"}]}], "]"}], + "]"}], ",", + RowBox[{"-", "2"}]}], "]"}]}], ";", "\[IndentingNewLine]", + RowBox[{"templateList", "=", + RowBox[{ + "\"\< template <\>\"", "<>", "templateList", "<>", "\"\<>\>\""}]}], + ";", "\[IndentingNewLine]", RowBox[{"head", "=", - RowBox[{"\"\