Add Simpson quadrature with trapzoidal extension

Add an additional Simpson quadrature which uses trapzoid quadrature over one cell when the number of grid points is not compatible with the  Simpson quadrature. Use this quadrature for the VparMu collision test.

See merge request gysela-developpers/gyselalibxx!760

--------------------------------------------

simulations/geometryVparMu/collisions/test_collSpVparMu.cpp

@@ -22,6 +22,7 @@
 #include "paraconfpp.hpp"
 #include "params.yaml.hpp"
 #include "pdi_out.yml.hpp"
+#include "simpson_quadrature.hpp"
 #include "species_info.hpp"
 #include "species_init.hpp"
 
@@ -90,8 +91,8 @@ int main(int argc, char** argv)
     DFieldMemVpar const coeff_intdvpar(neumann_spline_quadrature_coefficients<
                                        Kokkos::DefaultExecutionSpace>(idxrange_vpar, builder_vpar));
     SplineMuBuilder const builder_mu(idxrange_mu);
-    DFieldMemMu const coeff_intdmu(neumann_spline_quadrature_coefficients<
-                                   Kokkos::DefaultExecutionSpace>(idxrange_mu, builder_mu));
+    DFieldMemMu const coeff_intdmu(simpson_trapezoid_quadrature_coefficients_1d<
+                                   Kokkos::DefaultExecutionSpace>(idxrange_mu, Extremity::BACK));
     CollisionInfo const collision_info(conf_collision);
     CollisionSpVparMu<CollisionInfo, IdxRangeSpVparMu, GridVpar, GridMu, double> collision_operator(
             collision_info,

src/collisions/CollisionSpVparMu.hpp

@@ -172,11 +172,11 @@ class CollisionSpVparMu /* : public IRightHandSide */
         : m_operator_handle {}
         , m_hat_As {"m_hat_As", collisions_dimensions::get_idx_range<Species>(fdistrib_idx_range)}
         , m_hat_Zs {"m_hat_Zs", collisions_dimensions::get_idx_range<Species>(fdistrib_idx_range)}
+        , m_coeff_AD {"m_coeff_AD", collisions_dimensions::get_idx_range<GridR>(fdistrib_idx_range)}
         , m_mask_buffer_r {"m_mask_buffer_r", collisions_dimensions::get_idx_range<GridR>(fdistrib_idx_range)}
         , m_mask_LIM {"m_mask_LIM", IdxRangeThetaR {collisions_dimensions::get_idx_range<GridTheta, GridR>(fdistrib_idx_range)}}
         , m_B_norm {"m_B_norm", IdxRangeThetaR {collisions_dimensions::get_idx_range<GridTheta, GridR>(fdistrib_idx_range)}}
         , m_Bstar_s {"m_Bstar_s", IdxRangeSpThetaRVpar {collisions_dimensions::get_idx_range<Species, GridTheta, GridR, GridVpar>(fdistrib_idx_range)}}
-        , m_coeff_AD {"m_coeff_AD", collisions_dimensions::get_idx_range<GridR>(fdistrib_idx_range)}
         , m_mug {"m_mug", ddc::select<GridMu>(fdistrib_idx_range)}
         , m_vparg {"m_vparg", ddc::select<GridVpar>(fdistrib_idx_range)}
     {

src/multipatch/connectivity/edge.hpp

@@ -1,12 +1,9 @@
 // SPDX-License-Identifier: MIT
-
 #pragma once
 #include <ddc/ddc.hpp>
 
 #include "ddc_aliases.hpp"
-
-
-enum Extremity { FRONT, BACK };
+#include "geometry_descriptors.hpp"
 
 
 /**

src/quadrature/simpson_quadrature.hpp

@@ -4,7 +4,9 @@
 #include <ddc/ddc.hpp>
 
 #include "ddc_aliases.hpp"
+#include "geometry_descriptors.hpp"
 #include "quadrature_coeffs_nd.hpp"
+#include "trapezoid_quadrature.hpp"
 
 /**
  * @brief Get the Simpson coefficients in 1D.
@@ -14,15 +16,14 @@
  *
  * (x_3-x_1)(2-(x_3-x_2)/(x_2-x_1)) / 6 f(x_1) + (x_3-x_1)^3 / 6 / (x_3-x_2) / (x_2-x_1) f(x_2) + (x_3-x_1)(2-(x_2-x_1)/(x_3-x_2)) / 6 f(x_3)
  *
- * @param[in] idx_range
- * 	The index range on which the quadrature will be carried out.
- *
- * @return The quadrature coefficients for the Simpson method defined on the provided index range.
+ * @param[in] coefficients
+ * 	The field where the quadrature coefficients should be saved.
  */
 template <class ExecSpace, class Grid1D>
-DFieldMem<IdxRange<Grid1D>, typename ExecSpace::memory_space> simpson_quadrature_coefficients_1d(
-        IdxRange<Grid1D> const& idx_range)
+void fill_simpson_quadrature_coefficients_1d(
+        DField<IdxRange<Grid1D>, typename ExecSpace::memory_space> coefficients)
 {
+    IdxRange<Grid1D> idx_range = get_idx_range(coefficients);
     if constexpr (Grid1D::continuous_dimension_type::PERIODIC) {
         if (idx_range.size() % 2 != 0) {
             throw std::runtime_error(
@@ -34,9 +35,6 @@ DFieldMem<IdxRange<Grid1D>, typename ExecSpace::memory_space> simpson_quadrature
                                      "non-periodic direction.");
         }
     }
-    DFieldMem<IdxRange<Grid1D>, typename ExecSpace::memory_space> coefficients_alloc(idx_range);
-    DField<IdxRange<Grid1D>, typename ExecSpace::memory_space> const coefficients
-            = get_field(coefficients_alloc);
 
     Kokkos::parallel_for(
             "centre_left",
@@ -92,10 +90,88 @@ DFieldMem<IdxRange<Grid1D>, typename ExecSpace::memory_space> simpson_quadrature
                     coefficients(idx_range.front()) += dx_sum * (2. * dx_r - dx_l) / (6. * dx_r);
                 });
     }
+}
 
+/**
+ * @brief Get the Simpson coefficients in 1D.
+ *
+ * Calculate the quadrature coefficients for the Simpson method defined on the provided index range.
+ * The non-uniform form of the Simpson quadrature is:
+ *
+ * (x_3-x_1)(2-(x_3-x_2)/(x_2-x_1)) / 6 f(x_1) + (x_3-x_1)^3 / 6 / (x_3-x_2) / (x_2-x_1) f(x_2) + (x_3-x_1)(2-(x_2-x_1)/(x_3-x_2)) / 6 f(x_3)
+ *
+ * @param[in] idx_range
+ * 	The index range on which the quadrature will be carried out.
+ *
+ * @return The quadrature coefficients for the Simpson method defined on the provided index range.
+ */
+template <class ExecSpace, class Grid1D>
+DFieldMem<IdxRange<Grid1D>, typename ExecSpace::memory_space> simpson_quadrature_coefficients_1d(
+        IdxRange<Grid1D> const& idx_range)
+{
+    DFieldMem<IdxRange<Grid1D>, typename ExecSpace::memory_space> coefficients_alloc(idx_range);
+    fill_simpson_quadrature_coefficients_1d<ExecSpace>(get_field(coefficients_alloc));
     return coefficients_alloc;
 }
 
+/**
+ * @brief Get the Simpson coefficients in 1D.
+ *
+ * If the number of grid points is not compatible with the Simpson quadrature scheme then
+ * use a trapezoid formula over one cell at the specified extremity.
+ *
+ * @param[in] idx_range
+ * 	The index range on which the quadrature will be carried out.
+ * @param[in] trapezoid_extremity
+ *  The extremity where the trapezoid quadrature may be used.
+ *
+ * @return The quadrature coefficients for the Simpson method defined on the provided index range.
+ */
+template <class ExecSpace, class Grid1D>
+DFieldMem<IdxRange<Grid1D>, typename ExecSpace::memory_space>
+simpson_trapezoid_quadrature_coefficients_1d(
+        IdxRange<Grid1D> const& idx_range,
+        Extremity trapezoid_extremity)
+{
+    static_assert(
+            !Grid1D::continuous_dimension_type::PERIODIC,
+            "The extremity is non-sensical in a Periodic dimension");
+    try {
+        return simpson_quadrature_coefficients_1d<ExecSpace>(idx_range);
+    } catch (const std::runtime_error& error) {
+        DFieldMem<IdxRange<Grid1D>, typename ExecSpace::memory_space> coefficients_alloc(idx_range);
+        DField<IdxRange<Grid1D>, typename ExecSpace::memory_space> coefficients(
+                get_field(coefficients_alloc));
+        IdxStep<Grid1D> npts_to_remove(1);
+        if (trapezoid_extremity == Extremity::FRONT) {
+            fill_simpson_quadrature_coefficients_1d<ExecSpace>(
+                    coefficients_alloc[idx_range.remove_first(npts_to_remove)]);
+            Kokkos::parallel_for(
+                    "trapezoid_region",
+                    Kokkos::RangePolicy<ExecSpace>(0, 1),
+                    KOKKOS_LAMBDA(const int i) {
+                        Idx<Grid1D> idx_r = idx_range.back();
+                        double dx_cell = distance_at_left(idx_r);
+                        coefficients(idx_r) = dx_cell;
+                        coefficients(idx_r - 1) += dx_cell;
+                    });
+        } else {
+            fill_simpson_quadrature_coefficients_1d<ExecSpace>(
+                    coefficients_alloc[idx_range.remove_last(npts_to_remove)]);
+            Kokkos::parallel_for(
+                    "trapezoid_region",
+                    Kokkos::RangePolicy<ExecSpace>(0, 1),
+                    KOKKOS_LAMBDA(const int i) {
+                        Idx<Grid1D> idx_l = idx_range.front();
+                        double dx_cell = distance_at_right(idx_l);
+                        coefficients(idx_l) = dx_cell;
+                        coefficients(idx_l + 1) += dx_cell;
+                    });
+        }
+        return coefficients_alloc;
+    }
+}
+
 /**
  * @brief Get the simpson coefficients in ND.
  *

src/utils/geometry_descriptors.hpp

@@ -0,0 +1,4 @@
+// SPDX-License-Identifier: MIT
+#pragma once
+
+enum Extremity { FRONT, BACK };