Advection operator 1D for templated geometry

Add an advection operator which can solve the following equation :

$\\partial_t f  + A \\partial\_{x_i} f  = 0$

The advection direction can be spatial or along the velocity dimension. We need to provide an advection field. The function \\( f \\) can be only on a full domain (with Species, Space and Velocity), on a space/velocity domain (Space and Velocity) or on a space domain (Space).

See merge request gysela-developpers/gyselalibxx!481

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

Co-authored-by: Pauline Vidal <[email protected]>

src/advection/CMakeLists.txt

@@ -17,6 +17,7 @@ target_link_libraries("advection"
         sll::SLL
         gslx::interpolation
         gslx::speciesinfo
+        gslx::timestepper
         gslx::utils
 )
 

src/advection/README.md

@@ -15,3 +15,109 @@ Here the purpose is the advection along a direction on the velocity space dimens
 The dynamics of the motion on the velocity dimension are governed by the following equation, where E is the electric field.
 
 $$ \frac{df_s}{dt}= q_s \sqrt{\frac{m_e}{m_s}} E \frac{\partial f_s}{\partial v} $$
+
+
+## 1D advection with a given advection field
+The purpose of the BslAdvection1D operator is an advection along a given direction of the phase space. The advection field is given as input. 
+The dynamics of the motion are governed by the following equation. 
+```math 
+    \partial_t f(t,x) + A\partial_{x_i}(x')f(t,x) = 0, 
+    \qquad x \in \Omega, x' \in \Omega',
+```
+
+with 
+* $`f`$ the function to advect. It is defined on $`\Omega`$ domain and the time dimension; 
+* $`A`$ the advection field. It could be defined on a subdomain $`\Omega'\subset \Omega`$; 
+* and $`x_i`$ a given direction. The advection field domain has to be defined on this dimension for the time integration method.  
+
+
+### Example of use
+Here are some examples of equation types the BslAdvection1D operator can solve:
+* Equation on a 1D domain (and time dimension): 
+```math 
+    \partial_t f(t,x) + A_x(x)\partial_{x}f(t,x) = 0,
+    \qquad x \in \Omega,
+```
+
+Here $`\Omega' = \Omega \in \mathbb{R}`$. In the code, it would correspond to 
+```cpp
+DFieldX f(x_dom); 
+DFieldX A(x_dom); 
+using IDimInterest = IDimX; 
+```
+
+
+* Equation on a 2D domain (and time dimension): 
+```math 
+    \partial_t f(t,x,y) + A_x(x,y)\partial_{x}f(t,x,y) = 0,
+```
+
+Here $`\Omega' = \Omega \in \mathbb{R}^2`$. In the code, it would correspond to 
+```cpp
+DFieldXY f(xy_dom); 
+DFieldXY A(xy_dom); 
+using IDimInterest = IDimX; 
+```
+
+
+* Equation on a 2Dx2V domain (and time dimension): 
+```math 
+    \partial_t f(t,x,y,v_x,v_y) + A_x(x,y)\partial_{x}f(t,x,y,v_x,v_y) = 0,
+```
+
+Here $`\Omega' \in \mathbb{R}^2`$ and $`\Omega \in \mathbb{R}^4`$. In the code, it would correspond to 
+```cpp
+DFieldXYVxVy f(xyvxvy_dom); 
+DFieldXY A(xy_dom); 
+using IDimInterest = IDimX; 
+```
+
+
+*  Equation on a 1Dx1V domain (and time dimension): 
+```math 
+\begin{aligned}
+    & \partial_t f(t,x,v_x) + A_x(x,v_x)\partial_{x}f(t,x,v_x) = 0, 
+    \qquad \text{ with for instance, } A_x(x,v_x) = v_x, \\
+    & \text{or } \partial_t f(t,x,v_x) + A_{v_x}(x,v_x)\partial_{v_x}f(t,x,v_x) = 0, 
+    \qquad \text{ with for instance, } A_{v_x}(x,v_x) = E(x), 
+\end{aligned}
+```
+
+Here $`\Omega' = \Omega \in \mathbb{R}^2`$. In the code, it would correspond to 
+```cpp
+DFieldXVx f(xvx_dom); 
+DFieldXVx A(xvx_dom); 
+using IDimInterest = IDimX; 
+```
+or 
+```cpp
+using IDimInterest = IDimVx; 
+```
+
+*  Equation on a 1Dx1V domain with species dimension (and time dimension): 
+```math 
+\begin{aligned}
+    & \partial_t f_s(t,x,v_x) + A_{s,x}(x)\partial_{x}f_s(t,x,v_x) = 0, 
+\end{aligned}
+```
+
+Here $`\Omega' = \Omega \in \mathbb{R}^2`$. In the code, it would correspond to 
+```cpp
+DFieldSpXVx f(sp_xvx_dom); 
+DFieldSpX A(sp_x_dom); 
+using IDimInterest = IDimX; 
+```
+
+
+### Parameters
+The operator takes as templated parameters: 
+* IDimInterest: a dimension of interest (or advection dimension) wich refers to the dimension along wich we advect the function; 
+* AdvectionDomain: an advection domain, which refers to the domain where the advection field is defined. It has to contain the dimension of interest for the interpolation of the advection field in the time integration method; 
+* FunctionDomain: the full domain where the function we want to advect is defined; 
+* AdvectionFieldBuilder: a spline builder type for the advection field. 
+* AdvectionFieldEvaluator: a spline evaluator type for the advection field. 
+* TimeStepper: a time integration method (see [timestepper](./../timestepper/README.md)) to solve the characteristic equation. It has to be defined on the advection field domain. The feet have to be a Chunk of coordinates of the dimension of interest defined on the advection field domain.
+
+**Remark/Warning:** the BslAdvection1D operator is built with builder and evaluator for the advection field and interpolator for the function we want to advect. Theses operators have to be defined on the same domain as the advection field and function. For instance, if the advection field and/or the function are defined on the species dimension, then the interpolators have to contain the species dimension in its batched dimensions (see tests in the `tests/advection/` folder). 
+
+**Remark/Warning:** The advection field need to use interpolation on B-splines. So we cannot use other type of interpolator for the advection field. However there is no constraint on the interpolator of the advected function. 

src/advection/bsl_advection_1d.hpp

@@ -0,0 +1,255 @@
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+#include <ddc/ddc.hpp>
+#include <ddc/kernels/splines/deriv.hpp>
+
+#include "ddc_helper.hpp"
+#include "euler.hpp"
+#include "iinterpolator.hpp"
+
+
+/**
+ * @brief A class which computes the advection along the dimension of interest IDimInterest. 
+ * 
+ * This operator solves the following equation type
+ * 
+ * @f$ \partial_t f_s(t,x) + A_{s, x_i} (x') \cdot \partial_{x_i} f_s (t, x) = 0, \qquad x\in \Omega, x'\in\Omega'@f$
+ * 
+ * with 
+ *  * @f$ f @f$, a function defined on a domain @f$ \Omega @f$; 
+ *  * @f$ A @f$, an advection field defined on subdomain @f$ \Omega'\subset \Omega @f$; 
+ *  * @f$ x_i @f$, an advection dimension.
+ * 
+ * 
+ * The characteristic equation is solved on the advection domain @f$ \Omega'@f$. 
+ * Then the feet on @f$ \Omega@f$ are computed from the characteristic feet on @f$ \Omega'@f$ and the function 
+ * @f$ f @f$ is interpolated at the feet in @f$ \Omega @f$. 
+ * 
+ * The characteristic equation is solved using a time integration method (ITimeStepper). 
+ * 
+ * 
+ * @tparam IDimInterest
+ *          The dimension along which the advection is computed. 
+ *          It refers to the dimension of @f$ x_i @f$ in the equation. 
+ * @tparam AdvectionDomain
+ *          The domain @f$ \Omega' @f$ where the characteristic equation is solved. 
+ *          It also refers to the domain of the advection field. 
+ *          It had to also be defined on the IDimInterest for the time integration method. 
+ * @tparam FunctionDomain
+ *          The domain @f$ \Omega @f$ where allfdistribu is defined. 
+ * @tparam AdvectionFieldBuilder
+ *          The type of the spline builder for the advection field (see SplineBuilder). 
+ * @tparam AdvectionFieldEvaluator
+ *          The type of the spline evalutor for the advection field (see SplineEvaluator).  
+ * @tparam TimeStepper 
+ *          The time integration method applied to solve the characteristic equation. 
+ *          The method is picked among the child classes of ITimeStepper. 
+ */
+template <
+        class IDimInterest,
+        class AdvectionDomain,
+        class FunctionDomain,
+        class AdvectionFieldBuilder,
+        class AdvectionFieldEvaluator,
+        class TimeStepper
+        = Euler<device_t<ddc::Chunk<
+                        ddc::Coordinate<typename IDimInterest::continuous_dimension_type>,
+                        AdvectionDomain>>,
+                device_t<ddc::Chunk<double, AdvectionDomain>>>>
+class BslAdvection1D
+{
+private:
+    // Advection domain element:
+    using AdvectionIndex = typename AdvectionDomain::discrete_element_type;
+
+    // Full domain element:
+    using FunctionIndex = typename FunctionDomain::discrete_element_type;
+
+    // Advection dimension (or Interest dimension):
+    using RDimInterest = typename IDimInterest::continuous_dimension_type;
+    using CoordInterest = ddc::Coordinate<RDimInterest>;
+    using DomainInterest = ddc::DiscreteDomain<IDimInterest>;
+    using IndexInterest = typename DomainInterest::discrete_element_type;
+
+    // Type for the feet and advection field:
+    using FeetChunk = device_t<ddc::Chunk<CoordInterest, AdvectionDomain>>;
+    using FeetSpan = typename FeetChunk::span_type;
+    using FeetView = typename FeetChunk::view_type;
+
+    using AdvecFieldChunk = device_t<ddc::Chunk<double, AdvectionDomain>>;
+    using AdvecFieldSpan = typename AdvecFieldChunk::span_type;
+
+    using FunctionSpan = device_t<ddc::ChunkSpan<double, FunctionDomain>>;
+
+    // Type for spline representation of the advection field
+    using BSAdvectionDomain = typename AdvectionFieldBuilder::spline_domain_type;
+    using AdvecFieldSplineChunk = device_t<ddc::Chunk<double, BSAdvectionDomain>>;
+    using AdvecFieldSplineSpan = device_t<ddc::ChunkSpan<double, BSAdvectionDomain>>;
+
+    // Type for the derivatives of the advection field
+    using DerivDim = ddc::Deriv<RDimInterest>;
+    using DomainInterestDeriv = ddc::DiscreteDomain<DerivDim>;
+    using AdvecFieldDerivDomain = decltype(ddc::replace_dim_of<IDimInterest, DerivDim>(
+            std::declval<AdvectionDomain>(),
+            std::declval<DomainInterestDeriv>()));
+    using AdvecFieldDerivView = device_t<ddc::ChunkSpan<const double, AdvecFieldDerivDomain>>;
+
+
+    // Interpolators:
+    using FunctionPreallocatableInterpolatorType
+            = interpolator_on_domain_t<IPreallocatableInterpolator, IDimInterest, FunctionDomain>;
+    using FunctionInterpolatorType
+            = interpolator_on_domain_t<IInterpolator, IDimInterest, FunctionDomain>;
+
+
+    FunctionPreallocatableInterpolatorType const& m_function_interpolator;
+
+    AdvectionFieldBuilder const& m_adv_field_builder;
+    AdvectionFieldEvaluator const& m_adv_field_evaluator;
+
+    TimeStepper const& m_time_stepper;
+
+public:
+    /**
+     * @brief Constructor when the advection domain and the function domain are different. 
+     * 
+     * When AdvectionDomain and FunctionDomain are different, we need one interpolator for 
+     * each domain. 
+     * 
+     * We can also use it when we want two differents interpolators but defined on the same 
+     * domain (e.g. different boundary conditions for the evaluators).
+     * 
+     * @param[in] function_interpolator interpolator along the IDimInterest direction to interpolate 
+     *          the advected function (allfdistribu) on the domain of the function.
+     * @param[in] adv_field_builder builder along the IDimInterest direction to build a spline representation
+     *          of the advection field on the function domain. 
+     * @param[in] adv_field_evaluator evaluator along the IDimInterest direction to evaluate 
+     *          the advection field spline representation on the function domain.  
+     * @param[in] time_stepper time integration method for the characteristic equation. 
+     */
+    explicit BslAdvection1D(
+            FunctionPreallocatableInterpolatorType const& function_interpolator,
+            AdvectionFieldBuilder const& adv_field_builder,
+            AdvectionFieldEvaluator const& adv_field_evaluator,
+            TimeStepper const& time_stepper)
+        : m_function_interpolator(function_interpolator)
+        , m_adv_field_builder(adv_field_builder)
+        , m_adv_field_evaluator(adv_field_evaluator)
+        , m_time_stepper(time_stepper)
+    {
+    }
+
+    ~BslAdvection1D() = default;
+
+    /**
+     * @brief Advects allfdistribu along the advection dimension IDimInterest for a duration dt.
+     * 
+     * @param[in, out] allfdistribu Reference to the advected function, allocated on the device 
+     * @param[in] advection_field Reference to the advection field, allocated on the device.
+     * @param[in] advection_field_derivatives_min Reference to the advection field 
+     *              derivatives at the left side of the interest dimension, allocated on the device.
+     * @param[in] advection_field_derivatives_max Reference to the advection field 
+     *              derivatives at the right side of the interest dimension, allocated on the device.
+     * @param[in] dt Time step.
+     * 
+     * @return A reference to the allfdistribu array after advection on dt. 
+     */
+    FunctionSpan operator()(
+            FunctionSpan const allfdistribu,
+            AdvecFieldSpan const advection_field,
+            double const dt,
+            std::optional<AdvecFieldDerivView> const advection_field_derivatives_min = std::nullopt,
+            std::optional<AdvecFieldDerivView> const advection_field_derivatives_max
+            = std::nullopt) const
+    {
+        Kokkos::Profiling::pushRegion("BslAdvection1D");
+
+        // Get domains and operators .............................................................
+        FunctionDomain const function_dom = allfdistribu.domain();
+        AdvectionDomain const advection_dom = advection_field.domain();
+        auto const batch_dom = ddc::remove_dims_of(function_dom, advection_dom);
+        DomainInterest const interest_dom(advection_dom);
+
+        std::unique_ptr<FunctionInterpolatorType> const function_interpolator_ptr
+                = m_function_interpolator.preallocate();
+        FunctionInterpolatorType const& function_interpolator = *function_interpolator_ptr;
+
+
+        // Build spline representation of the advection field ....................................
+        AdvecFieldSplineChunk advection_field_coefs_alloc(m_adv_field_builder.spline_domain());
+        AdvecFieldSplineSpan advection_field_coefs = advection_field_coefs_alloc.span_view();
+
+        m_adv_field_builder(
+                advection_field_coefs,
+                advection_field.span_cview(),
+                advection_field_derivatives_min,
+                advection_field_derivatives_max);
+
+
+        // Initialise the characteristics on the mesh points .....................................
+        /*
+            For the time integration solver, the function we advect (here the characteristics)
+            need to be defined on the same domain as the advection field. We then work on space
+            slices of the characteristic feet.  
+        */
+        FeetChunk slice_feet_alloc(advection_dom);
+        FeetSpan slice_feet = slice_feet_alloc.span_view();
+        ddc::parallel_for_each(
+                Kokkos::DefaultExecutionSpace(),
+                advection_dom,
+                KOKKOS_LAMBDA(AdvectionIndex const idx) {
+                    slice_feet(idx) = ddc::coordinate(IndexInterest(idx));
+                });
+
+
+        // Compute the characteristic feet .......................................................
+        /*
+            We use a time stepper method to solve the charateristic equation. 
+            A TimeStepper needs a function to compute the updated advection field and a function to 
+            compute the updated feet. 
+                * update_adv_field: evaluate the advection field spline at the updated feet. 
+        */
+        // The function describing how the derivative of the evolve function is calculated.
+        std::function<void(AdvecFieldSpan, FeetView)> update_adv_field
+                = [&](AdvecFieldSpan updated_advection_field, FeetView slice_feet) {
+                      m_adv_field_evaluator(
+                              updated_advection_field,
+                              slice_feet,
+                              advection_field_coefs.span_cview());
+                  };
+
+
+        // Solve the characteristic equation with a time integration method
+        m_time_stepper
+                .update(Kokkos::DefaultExecutionSpace(),
+                        slice_feet.span_view(),
+                        -dt,
+                        update_adv_field);
+
+
+        // Interpolate the function ..............................................................
+        /*
+            To interpolate the function we want to advect, we build for the feet a Chunk defined 
+            on the domain where the function is defined. 
+        */
+        device_t<ddc::Chunk<CoordInterest, FunctionDomain>> feet_alloc(function_dom);
+        ddc::ChunkSpan feet = feet_alloc.span_view();
+        ddc::parallel_for_each(
+                Kokkos::DefaultExecutionSpace(),
+                function_dom,
+                KOKKOS_LAMBDA(FunctionIndex const idx) {
+                    AdvectionIndex slice_foot_index(idx);
+                    feet(idx) = slice_feet(slice_foot_index);
+                });
+
+
+        // Interpolate the function at the characteristic feet
+        function_interpolator(allfdistribu, feet);
+
+
+        Kokkos::Profiling::popRegion();
+        return allfdistribu;
+    }
+};

src/geometryRTheta/advection/README.md

@@ -73,8 +73,8 @@ $X^{n+1} = X^{k+1}$ once converged.
 $`X^n = X^{n+1} - \frac{dt}{6}  \left( k_1 + 4 k_2 + k_3 \right)`$
 	- with 
 		- $`k_1 =  A^{n+1}(X^{n+1})`$, 
-		- $`k_2 =  A(t^{n+1/2}, X^{n+1} - \frac{dt}{2} k_1)`$, 
-		- $`k_3 =  A(t^{n+1/2}, X^{n+1} - dt( 2k_2 - k_1))`$.
+		- $`k_2 =  A{n+1/2} (X^{n+1} - \frac{dt}{2} k_1)`$, 
+		- $`k_3 =  A{n+1/2} (X^{n+1} - dt( 2k_2 - k_1))`$.
 
 - Convergence order : 3.
 
@@ -86,9 +86,9 @@ $`X^n = X^{n+1} - \frac{dt}{6}  \left( k_1 + 4 k_2 + k_3 \right)`$
 $`X^n = X^{n+1} - \frac{dt}{6}  \left( k_1 + 2 k_2 + 2 k_3  + k_4\right)`$
 	- with 
 		- $`k_1 =  A^{n+1}(X^{n+1})`$, 
-		- $`k_2 =  A(t^{n+1/2}, X^{n+1} - \frac{dt}{2} k_1)`$, 
-		- $`k_3 =  A(t^{n+1/2}, X^{n+1}   \frac{dt}{2} k_2)`$, 
-		- $`k_4 =  A(t^{n}, X^{n+1} - dt k_3)`$.
+		- $`k_2 =  A^{n+1/2} (X^{n+1} - \frac{dt}{2} k_1)`$, 
+		- $`k_3 =  A^{n+1/2} (X^{n+1} - \frac{dt}{2} k_2)`$, 
+		- $`k_4 =  A^{n} (X^{n+1} - dt k_3)`$.
 
 - Convergence order : 4.