Add ND quadrature.

Generalise `Quadrature` class to accept N dimensions. Add a helper function `quadrature_coeffs_nd` to return N-D coefficients from 1D functions defining the quadrature along each of the dimensions.

.github/workflows/tests.yml

@@ -118,7 +118,7 @@ jobs:
           git add .
           git config --global user.name "GitHub CI Documentation builder"
           git config --global user.email "[email protected]"
-          git commit -m "Update docs"
+          git commit -m "Update docs" || true
           git push
         shell: bash
         env:

docs/Adding_docs.md

@@ -52,7 +52,7 @@ If other parts of the code are relevant to your description or may be relevant i
 
 ### Example
 
-```cpp
+```
 /**
  * @brief A class which provides an interpolating function.
  *
@@ -100,6 +100,19 @@ If the folder doesn't contain a `.private` file (indicating that the folder is n
 - [my_new_folder](./my_new_folder/README.md) : Short description of contents.
 ```
 
+## Documenting functions
+
+By default Doxygen only documents classes and class methods. In order to also document the functions in a file an additional descriptor must be added to the top of that file.
+The following is an example:
+```
+/**
+ * @file my_file.hpp
+ * Description of file.
+ */
+```
+
+This can be seen in action in the files in the folder `src/quadrature/`.
+
 ## Mathematical notation in documentation
 
 Mathematical notation can be used in Doxygen output.

src/README.md

@@ -8,5 +8,5 @@ The `src/` folder contains all the code necessary to build a gyrokinetic semi-La
 - [geometryXYVxVy](./geometryXYVxVy/README.md) - Code describing methods which are specific to a simulation with 2 spatial dimensions and 2 velocity dimension.
 - [interpolation](./interpolation/README.md) - Code describing interpolation methods.
 <!-- - [paraconfpp](./paraconfpp/README.md) - Paraconf utility functions. -->
-<!-- - [quadrature](./quadrature/README.md) - Code describing different quadrature methods. -->
+- [quadrature](./quadrature/README.md) - Code describing different quadrature methods.
 <!-- - [speciesinfo](./speciesinfo/README.md) - Code used to describe the different species. -->

src/quadrature/README.md

@@ -0,0 +1,14 @@
+# Quadrature Methods
+
+Quadrature methods are methods for calculating the integral of a function from its value at a given set of points $p_1, ..., p_n$.
+These methods can be written as a scalar product of the value of the function at the points $f(p_1), ..., f(p_n)$, with the quadrature coefficients $q_1, ..., q_n$.
+
+This folder provides the class Quadrature which integrates a function from its values at the points defined by the discrete domain(s) on which the class is defined.
+The class should be initialised with the quadrature coefficients.
+
+Helper functions provide the quadrature coefficients obtained using different quadrature methods.
+The methods currently implemented are:
+-  trapezoid_quadrature_coefficients()
+
+
+Additionally the function quadrature_coeffs_nd() helps define multi-dimensional quadrature methods from 1D methods.

src/quadrature/quadrature.hpp

@@ -6,30 +6,49 @@
 
 #include <ddc/ddc.hpp>
 
-template <class IDim>
+/**
+ * @brief A class providing an operator for integrating functions defined on a discrete domain.
+ */
+template <class... IDim>
 class Quadrature
 {
 private:
-    ddc::Chunk<double, ddc::DiscreteDomain<IDim>> m_coefficients;
+    ddc::Chunk<double, ddc::DiscreteDomain<IDim...>> m_coefficients;
 
 public:
-    explicit Quadrature(ddc::Chunk<double, ddc::DiscreteDomain<IDim>>&& coeffs)
+    /**
+     * @brief Create a Quadrature object.
+     * @param coeffs
+     * 	      The coefficients of the quadrature.
+     */
+    explicit Quadrature(ddc::Chunk<double, ddc::DiscreteDomain<IDim...>>&& coeffs)
         : m_coefficients(std::move(coeffs))
     {
     }
 
+    /**
+     * @brief Create a Quadrature object by copy.
+     * @param rhs The object being copied.
+     */
     Quadrature(Quadrature&& rhs) = default;
 
     ~Quadrature() = default;
 
-    double operator()(ddc::ChunkSpan<const double, ddc::DiscreteDomain<IDim>> const values) const
+    /**
+     * @brief An operator for calculating the integral of a function defined on a discrete domain.
+     * @param[in] values
+     *        The values of the function on the points of the discrete domain.
+     *
+     * @returns The integral of the function over the domain.
+     */
+    double operator()(ddc::ChunkSpan<const double, ddc::DiscreteDomain<IDim...>> const values) const
     {
-        assert(ddc::get_domain<IDim>(values) == ddc::get_domain<IDim>(m_coefficients));
+        assert(ddc::get_domain<IDim...>(values) == ddc::get_domain<IDim...>(m_coefficients));
         return ddc::transform_reduce(
                 values.domain(),
                 0.0,
                 ddc::reducer::sum<double>(),
-                [&](ddc::DiscreteElement<IDim> const ix) {
+                [&](ddc::DiscreteElement<IDim...> const ix) {
                     return m_coefficients(ix) * values(ix);
                 });
     }

src/quadrature/quadrature_coeffs_nd.hpp

@@ -0,0 +1,72 @@
+// SPDX-License-Identifier: MIT
+/**
+ * @file quadrature_coeffs_nd.hpp
+ * File providing helper functions for defining multi-dimensional quadrature methods.
+ */
+#pragma once
+
+
+/**
+ * @brief Helper function which creates ND dimensions from N 1D quadrature coefficient functions.
+ *
+ * @param domain
+ * 	The domain on which the coefficients will be defined.
+ * @param func
+ * 	The function which defines quadrature coefficients in the first dimension.
+ * @param o_funcs
+ * 	The functions which define quadrature coefficients in the subsequent dimensions.
+ *
+ * @returns The coefficients which define the quadrature method in ND.
+ */
+template <class IDim, class... ODims>
+ddc::Chunk<double, ddc::DiscreteDomain<IDim, ODims...>> quadrature_coeffs_nd(
+        ddc::DiscreteDomain<IDim, ODims...> const& domain,
+        std::function<ddc::Chunk<double, ddc::DiscreteDomain<IDim>>(ddc::DiscreteDomain<IDim>)>
+                func,
+        std::function<ddc::Chunk<double, ddc::DiscreteDomain<ODims>>(
+                ddc::DiscreteDomain<ODims>)>... o_funcs)
+{
+    // Get domains
+    ddc::DiscreteDomain<IDim> const current_dim_domain = ddc::select<IDim>(domain);
+    ddc::DiscreteDomain<ODims...> const other_dims_domain = ddc::select<ODims...>(domain);
+
+    // Get coefficients in the first dimension
+    ddc::Chunk<double, ddc::DiscreteDomain<IDim>> current_dim_coeffs = func(current_dim_domain);
+    // Get coefficients over all other dimensions
+    ddc::Chunk<double, ddc::DiscreteDomain<ODims...>> other_dim_coeffs
+            = quadrature_coeffs_nd(other_dims_domain, o_funcs...);
+
+    using DElemC = ddc::DiscreteElement<IDim>;
+    using DElemO = ddc::DiscreteElement<ODims...>;
+
+    // Calculate the combined coefficients over all dimensions
+    ddc::Chunk<double, ddc::DiscreteDomain<IDim, ODims...>> coefficients(domain);
+
+    ddc::for_each(current_dim_domain, [&](DElemC const idim) {
+        ddc::for_each(other_dims_domain, [&](DElemO const odim) {
+            coefficients(idim, odim) = current_dim_coeffs(idim) * other_dim_coeffs(odim);
+        });
+    });
+
+    return coefficients;
+}
+
+/**
+ * @brief Specialised 1D version of quadrature_coeffs_nd<IDim, ODims...>.
+ *
+ * @param domain
+ * 	The domain on which the coefficients will be defined.
+ * @param func
+ * 	The function which defines quadrature coefficients.
+ *
+ * @returns The coefficients which define the quadrature method in 1D.
+ */
+#include <ddc/ddc.hpp>
+template <class IDim>
+ddc::Chunk<double, ddc::DiscreteDomain<IDim>> quadrature_coeffs_nd(
+        ddc::DiscreteDomain<IDim> const& domain,
+        std::function<ddc::Chunk<double, ddc::DiscreteDomain<IDim>>(ddc::DiscreteDomain<IDim>)>
+                func)
+{
+    return func(domain);
+}

src/quadrature/trapezoid_quadrature.hpp

@@ -1,9 +1,23 @@
 // SPDX-License-Identifier: MIT
-
+/** @file trapezoid_quadrature.hpp
+ * File providing quadrature coefficients via the trapezoidal method.
+ */
 #pragma once
 
 #include <ddc/ddc.hpp>
 
+#include "quadrature_coeffs_nd.hpp"
+
+/**
+ * @brief Get the trapezoid coefficients in 1D.
+ *
+ * Calculate the quadrature coefficients for the trapezoid method defined on the provided domain.
+ *
+ * @param[in] domain
+ * 	The domain on which the quadrature will be carried out.
+ *
+ * @return The quadrature coefficients for the trapezoid method defined on the provided domain.
+ */
 template <class IDim>
 ddc::Chunk<double, ddc::DiscreteDomain<IDim>> trapezoid_quadrature_coefficients(
         ddc::DiscreteDomain<IDim> const& domain)
@@ -25,3 +39,20 @@ ddc::Chunk<double, ddc::DiscreteDomain<IDim>> trapezoid_quadrature_coefficients(
 
     return coefficients;
 }
+
+/**
+ * @brief Get the trapezoid coefficients in ND.
+ *
+ * Calculate the quadrature coefficients for the trapezoid method defined on the provided domain.
+ *
+ * @param[in] domain
+ * 	The domain on which the quadrature will be carried out.
+ *
+ * @return The quadrature coefficients for the trapezoid method defined on the provided domain.
+ */
+template <class... ODims>
+ddc::Chunk<double, ddc::DiscreteDomain<ODims...>> trapezoid_quadrature_coefficients(
+        ddc::DiscreteDomain<ODims...> const& domain)
+{
+    return quadrature_coeffs_nd(domain, (trapezoid_quadrature_coefficients<ODims>)...);
+}

tests/geometryXYVxVy/CMakeLists.txt

@@ -2,5 +2,25 @@
 
 cmake_minimum_required(VERSION 3.15)
 
+add_executable(unit_tests_xy_vxvy
+    ../main.cpp
+    quadrature.cpp
+)
+target_compile_features(unit_tests_xy_vxvy PUBLIC cxx_std_17)
+target_link_libraries(unit_tests_xy_vxvy
+    PUBLIC
+        DDC::PDI_Wrapper
+        GTest::gtest
+        GTest::gmock
+        paraconf::paraconf
+        vcx::geometry_xyvxvy
+        sll::splines
+        vcx::advection
+        vcx::quadrature
+)
+
+gtest_discover_tests(unit_tests_xy_vxvy
+    TEST_SUFFIX "_xy_vxvy")
+
 
 add_subdirectory(landau)

tests/geometryXYVxVy/quadrature.cpp

@@ -0,0 +1,119 @@
+// SPDX-License-Identifier: MIT
+
+#include <sll/bsplines_non_uniform.hpp>
+
+#include <gmock/gmock.h>
+#include <gtest/gtest.h>
+
+#include "geometry.hpp"
+#include "quadrature.hpp"
+#include "trapezoid_quadrature.hpp"
+
+TEST(QuadratureTest, ExactForConstantFunc)
+{
+    CoordX const x_min(0.0);
+    CoordX const x_max(M_PI);
+    IVectX const x_size(10);
+
+    CoordY const y_min(0.0);
+    CoordY const y_max(20.0);
+    IVectY const y_size(10);
+
+    // Creating mesh & supports
+    ddc::init_discrete_space<BSplinesX>(x_min, x_max, x_size);
+    ddc::init_discrete_space<BSplinesY>(y_min, y_max, y_size);
+
+    ddc::init_discrete_space<IDimX>(SplineInterpPointsX::get_sampling());
+    ddc::init_discrete_space<IDimY>(SplineInterpPointsY::get_sampling());
+
+    IDomainX const gridx(SplineInterpPointsX::get_domain());
+    IDomainY const gridy(SplineInterpPointsY::get_domain());
+
+    IDomainXY const gridxy(gridx, gridy);
+
+    Quadrature<IDimX, IDimY> const integrate(trapezoid_quadrature_coefficients(gridxy));
+
+    DFieldXY values(gridxy);
+
+    ddc::for_each(gridxy, [&](ddc::DiscreteElement<IDimX, IDimY> const idx) { values(idx) = 1.0; });
+    double integral = integrate(values);
+    double expected_val = (x_max - x_min) * (y_max - y_min);
+    EXPECT_LE(abs(integral - expected_val), 1e-9);
+}
+
+template <std::size_t N>
+struct X
+{
+    static bool constexpr PERIODIC = false;
+};
+
+template <std::size_t N>
+struct Y
+{
+    static bool constexpr PERIODIC = false;
+};
+
+template <std::size_t N>
+double compute_error(int n_elems)
+{
+    using DimX = X<N>;
+    using DimY = Y<N>;
+    using BSplinesX = UniformBSplines<DimX, 3>;
+    using BSplinesY = UniformBSplines<DimY, 3>;
+    using GrevillePointsX
+            = GrevilleInterpolationPoints<BSplinesX, BoundCond::GREVILLE, BoundCond::GREVILLE>;
+    using GrevillePointsY
+            = GrevilleInterpolationPoints<BSplinesY, BoundCond::GREVILLE, BoundCond::GREVILLE>;
+    using IDimX = typename GrevillePointsX::interpolation_mesh_type;
+    using IDimY = typename GrevillePointsY::interpolation_mesh_type;
+    using IDomainX = ddc::DiscreteDomain<IDimX>;
+    using IDomainY = ddc::DiscreteDomain<IDimY>;
+    using IDomainXY = ddc::DiscreteDomain<IDimX, IDimY>;
+    using DFieldXY = ddc::Chunk<double, IDomainXY>;
+
+    ddc::Coordinate<DimX> const x_min(0.0);
+    ddc::Coordinate<DimX> const x_max(M_PI);
+
+    ddc::Coordinate<DimY> const y_min(0.0);
+    ddc::Coordinate<DimY> const y_max(M_PI);
+
+    ddc::init_discrete_space<BSplinesX>(x_min, x_max, n_elems);
+    ddc::init_discrete_space<BSplinesY>(y_min, y_max, n_elems);
+
+    ddc::init_discrete_space<IDimX>(GrevillePointsX::get_sampling());
+    ddc::init_discrete_space<IDimY>(GrevillePointsY::get_sampling());
+    IDomainX const gridx(GrevillePointsX::get_domain());
+    IDomainY const gridy(GrevillePointsY::get_domain());
+    IDomainXY const gridxy(gridx, gridy);
+
+    Quadrature<IDimX, IDimY> const integrate(trapezoid_quadrature_coefficients(gridxy));
+
+    DFieldXY values(gridxy);
+
+    ddc::for_each(gridxy, [&](ddc::DiscreteElement<IDimX, IDimY> const idx) {
+        double const y_cos = cos(ddc::get<DimY>(ddc::coordinate(idx)));
+        values(idx) = sin(ddc::get<DimX>(ddc::coordinate(idx))) * y_cos * y_cos;
+    });
+    double integral = integrate(values);
+    return std::abs(integral - M_PI);
+}
+
+template <std::size_t... Is>
+std::array<double, sizeof...(Is)> compute_errors(std::index_sequence<Is...>, int n_elems)
+{
+    return std::array<double, sizeof...(Is)> {compute_error<Is>(n_elems *= 2)...};
+}
+
+TEST(QuadratureTest, UniformConverge)
+{
+    constexpr int NTESTS(4);
+
+    std::array<double, NTESTS> error = compute_errors(std::make_index_sequence<NTESTS>(), 50);
+
+    for (int i(1); i < NTESTS; ++i) {
+        EXPECT_LE(error[i], error[i - 1]);
+        double order = log(error[i - 1] / error[i]) / log(2.0);
+        double order_error = abs(2 - order);
+        EXPECT_LE(order_error, 1e-1);
+    }
+}