PartialDerivative operator

Rk: Development performed in GitHub

See merge request gysela-developpers/gyselalibxx!821

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

Co-authored-by: Emily Bourne <[email protected]>
Co-authored-by: Virginie Grandgirard <[email protected]>

.github/actions/build_code/action.yml

@@ -11,6 +11,6 @@ runs:
     - name: Compile code
       run: |
         cmake -DCMAKE_TOOLCHAIN=${{ inputs.toolchain }} -B build -S .
-        cmake --build build
+        cmake --build build --parallel 4
       shell: bash
 

.github/actions/run_tests/action.yml

@@ -5,10 +5,12 @@ runs:
   steps:
     - name: "Run tests"
       run: |
-        ctest -j 2 --timeout 5 --output-junit tests.xml
+        ctest -j 2 --timeout 5 --output-junit tests.xml --output-on-failure
       shell: bash
-      working-directory:
-        ./build
+      working-directory: ./build
+      env:
+        OMPI_ALLOW_RUN_AS_ROOT: 1
+        OMPI_ALLOW_RUN_AS_ROOT_CONFIRM: 1
     - name: Publish Test Report
       uses: mikepenz/action-junit-report@v4
       if: success() || failure() # always run even if the previous step fails

src/math_tools/README.md

@@ -2,6 +2,12 @@
 
 This folder contains mathematical classes and functions.
 
+## Derivative tools
+
+Functions for calculating derivatives with different methods:
+
+- `partial_derivatives.hpp` : calculate derivatives using spline interpolation
+
 ## Utility tools
 
 The l\_norm\_tools.hpp file contains functions computing the infinity norm. For now, it computes the infinity norm of 

src/math_tools/partial_derivatives.hpp

@@ -0,0 +1,81 @@
+
+// SPDX-License-Identifier: MIT
+/**
+ * @file partial_derivatives.hpp
+ * File containing functions to compute the partial derivatives
+ */
+
+#pragma once
+
+#include <ddc/ddc.hpp>
+
+#include "ddc_alias_inline_functions.hpp"
+#include "ddc_aliases.hpp"
+
+
+/**
+ * @brief A class which implementes a partial derivative operator
+ */
+template <class FieldXiBuilderBatched, class FieldXiEvaluatorBatched>
+class PartialDerivative
+{
+    static_assert(std::is_same_v<
+                  typename FieldXiBuilderBatched::batched_spline_domain_type,
+                  typename FieldXiEvaluatorBatched::batched_spline_domain_type>);
+    static_assert(std::is_same_v<
+                  typename FieldXiBuilderBatched::batched_interpolation_domain_type,
+                  typename FieldXiEvaluatorBatched::batched_evaluation_domain_type>);
+
+public:
+    /// The dimension Xi on which the partial derivative is calculated.
+    using DerivativeDirection = typename FieldXiBuilderBatched::continuous_dimension_type;
+
+    /// The index range on which this operator acts.
+    using IdxRangeFieldVal = typename FieldXiBuilderBatched::batched_interpolation_domain_type;
+
+    /// The type of the object that will be differentiated.
+    using DFieldVal = DField<IdxRangeFieldVal>;
+
+    /// The type of the calculated derivative.
+    using DConstFieldVal = DConstField<IdxRangeFieldVal>;
+
+private:
+    // Type for spline representation of the field
+    using IdxRangeBSFieldXi = typename FieldXiBuilderBatched::batched_spline_domain_type;
+    using FieldXiSplineMem = DFieldMem<IdxRangeBSFieldXi>;
+    using FieldXiSplineCoeffs = DField<IdxRangeBSFieldXi>;
+
+    FieldXiBuilderBatched const& m_fieldxi_builder;
+    FieldXiEvaluatorBatched const& m_fieldxi_evaluator;
+
+public:
+    /**
+    * @brief Construct an instance of the class PartialDerivative.
+    *
+    * @param fieldxi_builder Builder for intermediate interpolation representation.
+    * @param fieldxi_evaluator Evaluator for intermediate interpolation representation.
+    */
+    explicit PartialDerivative(
+            FieldXiBuilderBatched const& fieldxi_builder,
+            FieldXiEvaluatorBatched const& fieldxi_evaluator)
+        : m_fieldxi_builder(fieldxi_builder)
+        , m_fieldxi_evaluator(fieldxi_evaluator)
+    {
+    }
+
+    /**
+    * @brief Compute the partial derivative of @f$ F(X1,..,Xn)@f$ in Xi direction.
+    *
+    * @param[out] dfieldval_dxi Partial derivatives in Xi direction.
+    * @param[in] fieldval Values of the field @f$ F(X1,..,Xn)@f$.
+    */
+    void operator()(DFieldVal dfieldval_dxi, DConstFieldVal fieldval)
+    {
+        // Build spline representation of the field ....................................
+        FieldXiSplineMem fieldxi_coefs_alloc(m_fieldxi_builder.batched_spline_domain());
+        FieldXiSplineCoeffs fieldxi_coefs = get_field(fieldxi_coefs_alloc);
+
+        m_fieldxi_builder(fieldxi_coefs, get_const_field(fieldval));
+        m_fieldxi_evaluator.deriv(dfieldval_dxi, get_const_field(fieldxi_coefs));
+    }
+};

tests/math_tools/CMakeLists.txt

@@ -5,6 +5,7 @@ include(GoogleTest)
 
 add_executable(unit_tests_math_tools
     test_lnorm_tools.cpp
+    test_partial_derivatives.cpp
     ../main.cpp
 )
 target_link_libraries(unit_tests_math_tools
@@ -13,7 +14,7 @@ target_link_libraries(unit_tests_math_tools
         GTest::gtest
         GTest::gmock
         paraconf::paraconf
-	gslx::math_tools
+        gslx::math_tools
         gslx::utils
 )
 

tests/math_tools/test_partial_derivatives.cpp

@@ -0,0 +1,142 @@
+// SPDX-License-Identifier: MIT
+#include <ddc/ddc.hpp>
+
+#include <gtest/gtest.h>
+
+#include "ddc_aliases.hpp"
+#include "partial_derivatives.hpp"
+
+namespace {
+
+struct X
+{
+    static bool constexpr PERIODIC = false;
+};
+using CoordX = Coord<X>;
+struct BSplinesX : ddc::UniformBSplines<X, 3>
+{
+};
+auto constexpr SplineXBoundary = ddc::BoundCond::GREVILLE;
+using SplineInterpPointsX
+        = ddc::GrevilleInterpolationPoints<BSplinesX, SplineXBoundary, SplineXBoundary>;
+
+struct GridX : SplineInterpPointsX::interpolation_discrete_dimension_type
+{
+};
+using IdxX = Idx<GridX>;
+using IdxStepX = IdxStep<GridX>;
+using IdxRangeX = IdxRange<GridX>;
+
+struct Y
+{
+    static bool constexpr PERIODIC = false;
+};
+using CoordY = Coord<Y>;
+struct BSplinesY : ddc::UniformBSplines<Y, 3>
+{
+};
+auto constexpr SplineYBoundary = ddc::BoundCond::GREVILLE;
+using SplineInterpPointsY
+        = ddc::GrevilleInterpolationPoints<BSplinesY, SplineYBoundary, SplineYBoundary>;
+struct GridY : SplineInterpPointsY::interpolation_discrete_dimension_type
+{
+};
+using IdxY = Idx<GridY>;
+using IdxStepY = IdxStep<GridY>;
+using IdxRangeY = IdxRange<GridY>;
+
+using IdxXY = Idx<GridX, GridY>;
+using IdxRangeXY = IdxRange<GridX, GridY>;
+using DFieldMemXY = DFieldMem<IdxRangeXY>;
+using DFieldXY = DField<IdxRangeXY>;
+using DConstFieldXY = DConstField<IdxRangeXY>;
+
+// --- Operators ---
+using SplineXBuilder = ddc::SplineBuilder<
+        Kokkos::DefaultExecutionSpace,
+        Kokkos::DefaultExecutionSpace::memory_space,
+        BSplinesX,
+        GridX,
+        SplineXBoundary,
+        SplineXBoundary,
+        ddc::SplineSolver::LAPACK,
+        GridX,
+        GridY>;
+using SplineXEvaluator = ddc::SplineEvaluator<
+        Kokkos::DefaultExecutionSpace,
+        Kokkos::DefaultExecutionSpace::memory_space,
+        BSplinesX,
+        GridX,
+        ddc::ConstantExtrapolationRule<X>,
+        ddc::ConstantExtrapolationRule<X>,
+        GridX,
+        GridY>;
+
+
+static void TestPartialDerivativeDx()
+{
+    int n_elems_x(10);
+    int n_elems_y(20);
+
+    Coord<X> const x_min(0.0);
+    Coord<X> const x_max(1.0);
+    IdxStepX x_ncells(n_elems_x);
+
+    Coord<Y> const y_min(0.0);
+    Coord<Y> const y_max(2.0);
+    IdxStepY y_ncells(n_elems_y);
+
+    ddc::init_discrete_space<BSplinesX>(x_min, x_max, x_ncells);
+    ddc::init_discrete_space<GridX>(SplineInterpPointsX::get_sampling<GridX>());
+    IdxRangeX idxrange_x(SplineInterpPointsX::get_domain<GridX>());
+
+    ddc::init_discrete_space<BSplinesY>(y_min, y_max, y_ncells);
+    ddc::init_discrete_space<GridY>(SplineInterpPointsY::get_sampling<GridY>());
+    IdxRangeY idxrange_y(SplineInterpPointsY::get_domain<GridY>());
+
+    IdxRangeXY idxrange_xy(idxrange_x, idxrange_y);
+
+    SplineXBuilder const builder_x(idxrange_xy);
+    ddc::ConstantExtrapolationRule<X> bv_x_min(x_min);
+    ddc::ConstantExtrapolationRule<X> bv_x_max(x_max);
+    SplineXEvaluator const spline_evaluator_x(bv_x_min, bv_x_max);
+
+    DFieldMemXY field_xy_alloc(idxrange_xy);
+    DFieldXY field_xy = get_field(field_xy_alloc);
+
+    ddc::parallel_for_each(
+            Kokkos::DefaultExecutionSpace(),
+            idxrange_xy,
+            KOKKOS_LAMBDA(IdxXY const idx_xy) {
+                IdxX idx_x(idx_xy);
+                IdxY idx_y(idx_xy);
+                field_xy(idx_xy)
+                        = ddc::coordinate(idx_x) * ddc::coordinate(idx_x) * ddc::coordinate(idx_y);
+            });
+
+    PartialDerivative<SplineXBuilder, SplineXEvaluator> partial_dx(builder_x, spline_evaluator_x);
+    DFieldMemXY dfield_dx_xy_alloc(idxrange_xy);
+    DFieldXY dfield_dx_xy = get_field(dfield_dx_xy_alloc);
+    partial_dx(dfield_dx_xy, get_const_field(field_xy));
+
+    double max_error = ddc::parallel_transform_reduce(
+            Kokkos::DefaultExecutionSpace(),
+            idxrange_xy,
+            0.,
+            ddc::reducer::max<double>(),
+            KOKKOS_LAMBDA(IdxXY const idx_xy) {
+                IdxX idx_x(idx_xy);
+                IdxY idx_y(idx_xy);
+                double const dfield_dx_anal = 2. * ddc::coordinate(idx_x) * ddc::coordinate(idx_y);
+                return Kokkos::abs(dfield_dx_xy(idx_xy) - dfield_dx_anal);
+            });
+    EXPECT_LE(max_error, 1e-12);
+}
+
+
+TEST(PartialDerivative, PartialDerivativeDx)
+{
+    TestPartialDerivativeDx();
+}
+
+} // namespace

tests/timestepper/CMakeLists.txt

@@ -20,7 +20,7 @@ target_link_libraries(unit_tests_timestepper
         GTest::gtest
         GTest::gmock
         paraconf::paraconf
-	gslx::math_tools
+        gslx::math_tools
         gslx::multipatch_data_types
         gslx::multipatch_geometries
         gslx::timestepper