Add Euler and Crank-Nicolson methods

src/quadrature/README.md

@@ -9,6 +9,10 @@ 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()
+-  spline_quadrature_coefficients()
 
 
 Additionally the function quadrature_coeffs_nd() helps define multi-dimensional quadrature methods from 1D methods.
+
+In the compute_norms.hpp file, the functions compute_L1_norm() and compute_L2_norm() return the norms of a function with a given quadrature method. 
+The function compute_coeffs_on_mapping() add the Jacobian determinant of the mapping as factor of the quadrature coefficients. 

src/timestepper/README.md

@@ -8,6 +8,8 @@ y(t_0) = y_0
 $$
 
 The implemented methods are:
+- Explicit Euler (first order) (Euler)
+- Crank-Nicolson (second order) (CrankNicolson)
 - Second order Runge Kutta (RK2)
 - Third order Runge Kutta (RK3)
 - Fourth order Runge Kutta (RK4)

src/timestepper/crank_nicolson.hpp

@@ -0,0 +1,202 @@
+#pragma once
+#include <ddc_helper.hpp>
+#include <vector_field_common.hpp>
+
+#include "itimestepper.hpp"
+#include "utils_tools.hpp"
+
+
+/**
+ * @brief A class which provides an implementation of a Crank-Nicolson method.
+ *
+ * A class which provides an implementation of a Crank-Nicolson method in
+ * order to evolve values over time. The values may be either scalars or vectors. In the
+ * case of vectors the appropriate dimensions must be passed as template parameters.
+ * The values which evolve are defined on a domain.
+ *
+ * For the following ODE :
+ * @f$\partial_t y(t) = f(t, y(t)) @f$,
+ *
+ * the Crank-Nicolson method is given by :
+ * @f$ y^{k} =  y^{n} + \frac{dt}{2} \left(f(t^{n}, y^{n}) + f(t^{k}, y^{k}) \right)@f$.
+ *
+ * The method is an implicit method.
+ * If @f$ |y^{k+1} -  y^{k}| < \varepsilon @f$, then we set @f$ y^{n+1} = y^{k+1} @f$.
+ *
+ * The method is order 2.
+ *
+ */
+template <class ValChunk, class DerivChunk = ValChunk>
+class CrankNicolson : public ITimeStepper
+{
+private:
+    static_assert(ddc::is_chunk_v<ValChunk> or is_field_v<ValChunk>);
+    static_assert(ddc::is_chunk_v<DerivChunk> or is_field_v<DerivChunk>);
+
+    static_assert(
+            std::is_same_v<typename ValChunk::mdomain_type, typename DerivChunk::mdomain_type>);
+
+    using Domain = typename ValChunk::mdomain_type;
+
+    using Index = typename Domain::discrete_element_type;
+
+    using ValSpan = typename ValChunk::span_type;
+    using ValView = typename ValChunk::view_type;
+
+    using DerivSpan = typename DerivChunk::span_type;
+    using DerivView = typename DerivChunk::view_type;
+
+    ValChunk m_y_init;
+    ValChunk m_y_old;
+    DerivChunk m_k1;
+    DerivChunk m_k_new;
+    DerivChunk m_k_total;
+
+    int const m_max_counter;
+    double const m_epsilon;
+
+public:
+    /**
+     * @brief Create a CrankNicolson object.
+     * @param[in] dom
+     *      The domain on which the points which evolve over time are defined.
+     * @param[in] counter
+     *      The maximal number of loops for the implicit method.
+     * @param[in] epsilon
+     *      The @f$ \varepsilon @f$ upperbound of the difference of two steps
+     *      in the implicit method: @f$ |y^{k+1} -  y^{k}| < \varepsilon @f$.
+     */
+    CrankNicolson(Domain dom, int const counter = int(20), double const epsilon = 1e-12)
+        : m_max_counter(counter)
+        , m_epsilon(epsilon)
+    {
+        m_k1 = DerivChunk(dom);
+        m_k_new = DerivChunk(dom);
+        m_k_total = DerivChunk(dom);
+        m_y_init = ValChunk(dom);
+        m_y_old = ValChunk(dom);
+    }
+
+    /**
+     * @brief Carry out one step of the Crank-Nicolson scheme.
+     *
+     * This function is a wrapper around the update function below. The values of the function are
+     * updated using the trivial method $f += df * dt$. This is the standard method however some
+     * cases may need a more complex update function which is why the more explicit method is
+     * also provided.
+     *
+     * @param[inout] y
+     *     The value(s) which should be evolved over time defined on each of the dimensions at each point
+     *     of the domain.
+     * @param[in] dt
+     *     The time step over which the values should be evolved.
+     * @param[in] dy
+     *     The function describing how the derivative of the evolve function is calculated.
+     */
+    void update(ValSpan y, double dt, std::function<void(DerivSpan, ValView)> dy)
+    {
+        static_assert(ddc::is_chunk_v<ValChunk>);
+        update(y, dt, dy, [&](ValSpan y, DerivView dy, double dt) {
+            ddc::for_each(y.domain(), [&](Index const idx) { y(idx) = y(idx) + dy(idx) * dt; });
+        });
+    }
+
+    /**
+     * @brief Carry out one step of the Crank-Nicolson scheme.
+     *
+     * @param[inout] y
+     *     The value(s) which should be evolved over time defined on each of the dimensions at each point
+     *     of the domain.
+     * @param[in] dt
+     *     The time step over which the values should be evolved.
+     * @param[in] dy
+     *     The function describing how the derivative of the evolve function is calculated.
+     * @param[in] y_update
+     *     The function describing how the value(s) are updated using the derivative.
+     */
+    void update(
+            ValSpan y,
+            double dt,
+            std::function<void(DerivSpan, ValView)> dy,
+            std::function<void(ValSpan, DerivView, double)> y_update)
+    {
+        copy(m_y_init, y);
+
+        // --------- Calculate k1 ------------
+        // Calculate k1 = f(y_n)
+        dy(m_k1, y);
+
+        // -------- Calculate k_new ----------
+        bool not_converged = true;
+        int counter = 0;
+        do {
+            counter++;
+
+            // Calculate k_new = f(y_new)
+            dy(m_k_new, y);
+
+            // Calculation of step
+            ddc::for_each(m_k_total.domain(), [&](Index const i) {
+                // k_total = k1 + k_new
+                if constexpr (is_field_v<DerivChunk>) {
+                    fill_k_total(i, m_k1(i) + m_k_new(i));
+                } else {
+                    m_k_total(i) = m_k1(i) + m_k_new(i);
+                }
+            });
+
+            // Save the old characteristic feet
+            copy(m_y_old, y);
+
+            // Re-initiliase the characteristic feet
+            copy(y, m_y_init);
+
+            // Calculate y_new := y_n + h/2*(k_1 + k_new)
+            y_update(y, m_k_total, 0.5 * dt);
+
+
+            // Check convergence
+            not_converged = not have_converged(m_y_old, y);
+
+
+        } while (not_converged and (counter < m_max_counter));
+    };
+
+private:
+    void copy(ValSpan copy_to, ValView copy_from)
+    {
+        if constexpr (is_field_v<ValSpan>) {
+            ddcHelper::deepcopy(copy_to, copy_from);
+        } else {
+            ddc::deepcopy(copy_to, copy_from);
+        }
+    }
+
+    template <class... DDims>
+    void fill_k_total(Index i, ddc::Coordinate<DDims...> new_val)
+    {
+        ((ddcHelper::get<DDims>(m_k_total)(i) = ddc::get<DDims>(new_val)), ...);
+    }
+
+
+    // Check if the relative difference of the function between
+    // two times steps is below epsilon.
+    bool have_converged(ValView y_old, ValView y_new)
+    {
+        auto const dom = y_old.domain();
+
+        double norm_old = 0.;
+        ddc::for_each(dom, [&](Index const idx) {
+            const double abs_val = norm_inf(y_old(idx));
+            norm_old = norm_old > abs_val ? norm_old : abs_val;
+        });
+
+        double max_diff = 0.;
+        ddc::for_each(dom, [&](Index const idx) {
+            const double abs_diff = norm_inf(y_old(idx) - y_new(idx));
+            max_diff = max_diff > abs_diff ? max_diff : abs_diff;
+        });
+
+        return (max_diff / norm_old) < m_epsilon;
+    }
+};

src/timestepper/euler.hpp

@@ -0,0 +1,109 @@
+#pragma once
+#include <array>
+
+#include <ddc_helper.hpp>
+#include <vector_field_common.hpp>
+
+#include "itimestepper.hpp"
+
+/**
+ * @brief A class which provides an implementation of an explicit Euler method.
+ *
+ * A class which provides an implementation of an explicit Euler method in
+ * order to evolve values over time. The values may be either scalars or vectors. In the
+ * case of vectors the appropriate dimensions must be passed as template parameters.
+ * The values which evolve are defined on a domain.
+ *
+ * For the following ODE :
+ * @f$\partial_t y(t) = f(t, y(t)) @f$,
+ *
+ * the explicit Euler method is given by :
+ * @f$ y^{n+1} =  y^{n} + dt f(t^{n}, y^{n})@f$.
+ *
+ * The method is order 1.
+ *
+ */
+template <class ValChunk, class DerivChunk = ValChunk>
+class Euler : public ITimeStepper
+{
+private:
+    static_assert(ddc::is_chunk_v<ValChunk> or is_field_v<ValChunk>);
+    static_assert(ddc::is_chunk_v<DerivChunk> or is_field_v<DerivChunk>);
+
+    static_assert(
+            std::is_same_v<typename ValChunk::mdomain_type, typename DerivChunk::mdomain_type>);
+
+    using Domain = typename ValChunk::mdomain_type;
+
+    using Index = typename Domain::discrete_element_type;
+
+    using ValSpan = typename ValChunk::span_type;
+    using ValView = typename ValChunk::view_type;
+
+    using DerivSpan = typename DerivChunk::span_type;
+    using DerivView = typename DerivChunk::view_type;
+
+    DerivChunk m_k1;
+
+public:
+    /**
+     * @brief Create a Euler object.
+     * @param[in] dom The domain on which the points which evolve over time are defined.
+     */
+    Euler(Domain dom)
+    {
+        m_k1 = DerivChunk(dom);
+    }
+
+    /**
+     * @brief Carry out one step of the explicit Euler scheme.
+     *
+     * This function is a wrapper around the update function below. The values of the function are
+     * updated using the trivial method $f += df * dt$. This is the standard method however some
+     * cases may need a more complex update function which is why the more explicit method is
+     * also provided.
+     *
+     * @param[inout] y
+     *     The value(s) which should be evolved over time defined on each of the dimensions at each point
+     *     of the domain.
+     * @param[in] dt
+     *     The time step over which the values should be evolved.
+     * @param[in] dy
+     *     The function describing how the derivative of the evolve function is calculated.
+     */
+    void update(ValSpan y, double dt, std::function<void(DerivSpan, ValView)> dy)
+    {
+        static_assert(ddc::is_chunk_v<ValChunk>);
+        update(y, dt, dy, [&](ValSpan y, DerivView dy, double dt) {
+            ddc::for_each(y.domain(), [&](Index const idx) { y(idx) = y(idx) + dy(idx) * dt; });
+        });
+    }
+
+    /**
+     * @brief Carry out one step of the explicit Euler scheme.
+     *
+     * @param[inout] y
+     *     The value(s) which should be evolved over time defined on each of the dimensions at each point
+     *     of the domain.
+     * @param[in] dt
+     *     The time step over which the values should be evolved.
+     * @param[in] dy
+     *     The function describing how the derivative of the evolve function is calculated.
+     * @param[in] y_update
+     *     The function describing how the value(s) are updated using the derivative.
+     */
+    void update(
+            ValSpan y,
+            double dt,
+            std::function<void(DerivSpan, ValView)> dy,
+            std::function<void(ValSpan, DerivView, double)> y_update)
+    {
+        // --------- Calculate k1 ------------
+        // Calculate k1 = f(y_n)
+        dy(m_k1, y);
+
+        // ----------- Update y --------------
+        // Calculate y_new := y_n + h*k_1
+        y_update(y, m_k1, dt);
+    };
+};

src/timestepper/rk2.hpp

@@ -13,6 +13,18 @@
  * order to evolve values over time. The values may be either scalars or vectors. In the
  * case of vectors the appropriate dimensions must be passed as template parameters.
  * The values which evolve are defined on a domain.
+ *
+ *
+ * For the following ODE :
+ * @f$\partial_t y(t) = f(t, y(t)) @f$,
+ *
+ * the Runge-Kutta 2 method is given by :
+ * @f$ y^{n+1} =  y^{n} + dt k_2 @f$,
+ *
+ * with
+ *
+ * - @f$ k_1 = f(t^{n}, y^{n}) @f$,
+ * - @f$ k_2 = f(t^{n+1/2}, y^{n} + \frac{dt}{2} k_1 ) @f$,
  */
 template <class ValChunk, class DerivChunk = ValChunk>
 class RK2 : public ITimeStepper

src/timestepper/rk3.hpp

@@ -1,5 +1,6 @@
 #pragma once
 #include <ddc_helper.hpp>
+#include <vector_field_common.hpp>
 
 #include "itimestepper.hpp"
 
@@ -11,6 +12,19 @@
  * order to evolve values over time. The values may be either scalars or vectors. In the
  * case of vectors the appropriate dimensions must be passed as template parameters.
  * The values which evolve are defined on a domain.
+ *
+ * For the following ODE :
+ * @f$\partial_t y(t) = f(t, y(t)) @f$,
+ *
+ * the Runge-Kutta 3 method is given by :
+ * @f$ y^{n+1} =  y^{n} + \frac{dt}{6} \left(k_1 + 4 k_2 + k_3 \right) @f$,
+ *
+ * with
+ *
+ * - @f$ k_1 = f(t^{n}, y^{n}) @f$,
+ * - @f$ k_2 = f(t^{n+1/2}, y^{n} + \frac{dt}{2} k_1 ) @f$,
+ * - @f$ k_3 = f(t^{n+1/2}, y^{n} + dt ( 2 k_2 - k_1) ) @f$.
+ *
  */
 template <class ValChunk, class DerivChunk = ValChunk>
 class RK3 : public ITimeStepper