From bcab41cf311875f882c4e0bf74a1683f2ea55811 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tr=C3=A9vis=20Morvany?= Date: Fri, 26 Sep 2025 09:55:31 +0200 Subject: [PATCH 1/9] Add a 3D spline builder for the periodic case --- include/ddc/kernels/splines.hpp | 2 + .../ddc/kernels/splines/spline_builder_3d.hpp | 739 +++++ .../kernels/splines/spline_evaluator_3d.hpp | 2425 +++++++++++++++++ tests/splines/CMakeLists.txt | 20 + tests/splines/batched_3d_spline_builder.cpp | 852 ++++++ tests/splines/evaluator_3d.hpp | 109 + tests/splines/spline_error_bounds.hpp | 128 + 7 files changed, 4275 insertions(+) create mode 100644 include/ddc/kernels/splines/spline_builder_3d.hpp create mode 100644 include/ddc/kernels/splines/spline_evaluator_3d.hpp create mode 100644 tests/splines/batched_3d_spline_builder.cpp create mode 100644 tests/splines/evaluator_3d.hpp diff --git a/include/ddc/kernels/splines.hpp b/include/ddc/kernels/splines.hpp index ca4777461..4148b9059 100644 --- a/include/ddc/kernels/splines.hpp +++ b/include/ddc/kernels/splines.hpp @@ -18,8 +18,10 @@ #include "splines/spline_boundary_conditions.hpp" #include "splines/spline_builder.hpp" #include "splines/spline_builder_2d.hpp" +#include "splines/spline_builder_3d.hpp" #include "splines/spline_evaluator.hpp" #include "splines/spline_evaluator_2d.hpp" +#include "splines/spline_evaluator_3d.hpp" #include "splines/spline_traits.hpp" #include "splines/splines_linear_problem.hpp" #include "splines/splines_linear_problem_2x2_blocks.hpp" diff --git a/include/ddc/kernels/splines/spline_builder_3d.hpp b/include/ddc/kernels/splines/spline_builder_3d.hpp new file mode 100644 index 000000000..a3f88ea1c --- /dev/null +++ b/include/ddc/kernels/splines/spline_builder_3d.hpp @@ -0,0 +1,739 @@ +// Copyright (C) The DDC development team, see COPYRIGHT.md file +// +// SPDX-License-Identifier: MIT + +#pragma once + +#include +#include +#include +#include + +#include + +#include "spline_builder.hpp" +#include "spline_builder_2d.hpp" + +namespace ddc { + +/** + * @brief A class for creating a 3D spline approximation of a function. + * + * A class which contains an operator () which can be used to build a 3D spline approximation + * of a function. A 3D spline approximation uses a cross-product between three 1D SplineBuilder. + * + * @see SplineBuilder + */ +template < + class ExecSpace, + class MemorySpace, + class BSpline1, + class BSpline2, + class BSpline3, + class DDimI1, + class DDimI2, + class DDimI3, + ddc::BoundCond BcLower1, + ddc::BoundCond BcUpper1, + ddc::BoundCond BcLower2, + ddc::BoundCond BcUpper2, + ddc::BoundCond BcLower3, + ddc::BoundCond BcUpper3, + ddc::SplineSolver Solver> +class SplineBuilder3D +{ +public: + /// @brief The type of the Kokkos execution space used by this class. + using exec_space = ExecSpace; + + /// @brief The type of the Kokkos memory space used by this class. + using memory_space = MemorySpace; + + /// @brief The type of the SplineBuilder used by this class to spline-approximate along first dimension. + using builder_type1 = ddc:: + SplineBuilder; + + /// @brief The type of the SplineBuilder used by this class to spline-approximate along second dimension. + using builder_type2 = ddc:: + SplineBuilder; + + /// @brief The type of the SplineBuilder used by this class to spline-approximate along third dimension. + using builder_type3 = ddc:: + SplineBuilder; + + // FIXME: change this when implementing the derivatives + /// @brief The type of the SplineBuilder used by this class to spline-approximate the second-dimension-derivatives along first dimension. + using builder_deriv_type1 = ddc::SplineBuilder2D< + ExecSpace, + MemorySpace, + BSpline1, + BSpline2, + DDimI1, + DDimI2, + BcLower1, + BcUpper1, + BcLower2, + BcLower2, + Solver>; + + /// @brief The type of the first interpolation continuous dimension. + using continuous_dimension_type1 = typename builder_type1::continuous_dimension_type; + + /// @brief The type of the second interpolation continuous dimension. + using continuous_dimension_type2 = typename builder_type2::continuous_dimension_type; + + /// @brief The type of the third interpolation continuous dimension. + using continuous_dimension_type3 = typename builder_type3::continuous_dimension_type; + + /// @brief The type of the first interpolation discrete dimension. + using interpolation_discrete_dimension_type1 = + typename builder_type1::interpolation_discrete_dimension_type; + + /// @brief The type of the second interpolation discrete dimension. + using interpolation_discrete_dimension_type2 = + typename builder_type2::interpolation_discrete_dimension_type; + + /// @brief The type of the third interpolation discrete dimension. + using interpolation_discrete_dimension_type3 = + typename builder_type3::interpolation_discrete_dimension_type; + + /// @brief The type of the B-splines in the first dimension. + using bsplines_type1 = typename builder_type1::bsplines_type; + + /// @brief The type of the B-splines in the second dimension. + using bsplines_type2 = typename builder_type2::bsplines_type; + + /// @brief The type of the B-splines in the third dimension. + using bsplines_type3 = typename builder_type3::bsplines_type; + + /// @brief The type of the Deriv domain on boundaries in the first dimension. + using deriv_type1 = typename builder_type1::deriv_type; + + /// @brief The type of the Deriv domain on boundaries in the second dimension. + using deriv_type2 = typename builder_type2::deriv_type; + + /// @brief The type of the Deriv domain on boundaries in the third dimension. + using deriv_type3 = typename builder_type3::deriv_type; + + /// @brief The type of the domain for the interpolation mesh in the first dimension. + using interpolation_domain_type1 = + typename builder_type1::interpolation_discrete_dimension_type; + + /// @brief The type of the domain for the interpolation mesh in the second dimension. + using interpolation_domain_type2 = + typename builder_type2::interpolation_discrete_dimension_type; + + /// @brief The type of the domain for the interpolation mesh in the third dimension. + using interpolation_domain_type3 = + typename builder_type3::interpolation_discrete_dimension_type; + + /// @brief The type of the domain for the interpolation mesh in the 3D dimension. + using interpolation_domain_type = ddc::DiscreteDomain< + interpolation_discrete_dimension_type1, + interpolation_discrete_dimension_type2, + interpolation_discrete_dimension_type3>; + + /** + * @brief The type of the whole domain representing interpolation points. + * + * @tparam The batched discrete domain on which the interpolation points are defined. + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using batched_interpolation_domain_type = BatchedInterpolationDDom; + + /** + * @brief The type of the batch domain (obtained by removing the dimensions of interest + * from the whole domain). + * + * @tparam The batched discrete domain on which the interpolation points are defined. + * + * Example: For batched_interpolation_domain_type = DiscreteDomain and dimensions of interest X, Y and Z, + * this is DiscreteDomain. + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using batch_domain_type = ddc::remove_dims_of_t< + BatchedInterpolationDDom, + interpolation_discrete_dimension_type1, + interpolation_discrete_dimension_type2, + interpolation_discrete_dimension_type3>; + + /** + * @brief The type of the whole spline domain (cartesian product of 3D spline domain + * and batch domain) preserving the order of dimensions. + * + * @tparam The batched discrete domain on which the interpolation points are defined. + * + * Example: For batched_interpolation_domain_type = DiscreteDomain and dimensions of interest X, Y and Z + * (associated to B-splines tags BSplinesX, BSplinesY and BSplinesZ), this is DiscreteDomain + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using batched_spline_domain_type + = ddc::detail::convert_type_seq_to_discrete_domain_t, + ddc::detail::TypeSeq< + interpolation_discrete_dimension_type1, + interpolation_discrete_dimension_type2, + interpolation_discrete_dimension_type3>, + ddc::detail::TypeSeq>>; + + /** + * @brief The type of the whole Derivs domain (cartesian product of the 1D Deriv domain + * and the associated batch domain) in the first dimension, preserving the order of dimensions. + * + * @tparam The batched discrete domain on which the interpolation points are defined. + * + * Example: For batched_interpolation_domain_type = DiscreteDomain and dimensions of interest X, Y and Z, + * this is DiscreteDomain, Y, Z, T>. + */ + // template < + // class BatchedInterpolationDDom, + // class = std::enable_if_t>> + // using batched_derivs_domain_type1 = + // typename builder_type1::template batched_derivs_domain_type; + + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using batched_derivs_domain_type_1_2 = ddc::replace_dim_of_t< + ddc::replace_dim_of_t< + BatchedInterpolationDDom, + interpolation_discrete_dimension_type1, + deriv_type1>, + interpolation_discrete_dimension_type2, + deriv_type2>; + + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using batched_derivs_domain_type_2_3 = ddc::replace_dim_of_t< + ddc::replace_dim_of_t< + BatchedInterpolationDDom, + interpolation_discrete_dimension_type2, + deriv_type2>, + interpolation_discrete_dimension_type3, + deriv_type3>; + + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using batched_derivs_domain_type_1_3 = ddc::replace_dim_of_t< + ddc::replace_dim_of_t< + BatchedInterpolationDDom, + interpolation_discrete_dimension_type3, + deriv_type3>, + interpolation_discrete_dimension_type3, + deriv_type3>; + + /** + * @brief The type of the whole Derivs domain (cartesian product of the 1D Deriv domain + * and the associated batch domain) in the second dimension, preserving the order of dimensions. + * + * @tparam The batched discrete domain on which the interpolation points are defined. + * + * Example: For batched_interpolation_domain_type = DiscreteDomain and dimensions of interest X, Y, and Z + * this is DiscreteDomain, Z, T>. + */ + // template < + // class BatchedInterpolationDDom, + // class = std::enable_if_t>> + // using batched_derivs_domain_type2 = ddc::replace_dim_of_t< + // BatchedInterpolationDDom, + // interpolation_discrete_dimension_type2, + // deriv_type2>; + + /** + * @brief The type of the whole Derivs domain (cartesian product of the 1D Deriv domain + * and the associated batch domain) in the third dimension, preserving the order of dimensions. + * + * @tparam The batched discrete domain on which the interpolation points are defined. + * + * Example: For batched_interpolation_domain_type = DiscreteDomain and dimensions of interest X, Y, and Z + * this is DiscreteDomain, T>. + */ + // template < + // class BatchedInterpolationDDom, + // class = std::enable_if_t>> + // using batched_derivs_domain_type3 = ddc::replace_dim_of_t< + // BatchedInterpolationDDom, + // interpolation_discrete_dimension_type3, + // deriv_type3>; + + /** + * @brief The type of the whole Derivs domain (cartesian product of the 3D Deriv domain + * and the batch domain), preserving the order of dimensions. + * + * @tparam The batched discrete domain on which the interpolation points are defined. + * + * Example: For batched_interpolation_domain_type = DiscreteDomain and dimensions of interest X, Y and Z, + * this is DiscreteDomain, Deriv, Deriv, T>. + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using batched_derivs_domain_type + = ddc::detail::convert_type_seq_to_discrete_domain_t, + ddc::detail::TypeSeq< + interpolation_discrete_dimension_type1, + interpolation_discrete_dimension_type2, + interpolation_discrete_dimension_type3>, + ddc::detail::TypeSeq>>; + +private: + builder_type1 m_spline_builder1; + builder_deriv_type1 m_spline_builder_deriv1; + builder_type2 m_spline_builder2; + builder_type3 m_spline_builder3; + +public: + /** + * @brief Build a SplineBuilder2D acting on interpolation_domain. + * + * @param interpolation_domain The domain on which the interpolation points are defined, without the batch dimensions. + * + * @param cols_per_chunk A parameter used by the slicer (internal to the solver) to define the size + * of a chunk of right-hand-sides of the linear problem to be computed in parallel (chunks are treated + * by the linear solver one-after-the-other). + * This value is optional. If no value is provided then the default value is chosen by the requested solver. + * + * @param preconditioner_max_block_size A parameter used by the slicer (internal to the solver) to + * define the size of a block used by the Block-Jacobi preconditioner. + * This value is optional. If no value is provided then the default value is chosen by the requested solver. + * + * @see SplinesLinearProblemSparse + */ + explicit SplineBuilder3D( + interpolation_domain_type const& interpolation_domain, + std::optional cols_per_chunk = std::nullopt, + std::optional preconditioner_max_block_size = std::nullopt) + : m_spline_builder1(interpolation_domain, cols_per_chunk, preconditioner_max_block_size) + , m_spline_builder_deriv1(interpolation_domain) + , m_spline_builder2(interpolation_domain, cols_per_chunk, preconditioner_max_block_size) + , m_spline_builder3(interpolation_domain, cols_per_chunk, preconditioner_max_block_size) + { + } + + /** + * @brief Build a SplineBuilder2D acting on the interpolation domain contained in batched_interpolation_domain. + * + * @param batched_interpolation_domain The domain on which the interpolation points are defined. + * + * @param cols_per_chunk A parameter used by the slicer (internal to the solver) to define the size + * of a chunk of right-hand-sides of the linear problem to be computed in parallel (chunks are treated + * by the linear solver one-after-the-other). + * This value is optional. If no value is provided then the default value is chosen by the requested solver. + * + * @param preconditioner_max_block_size A parameter used by the slicer (internal to the solver) to + * define the size of a block used by the Block-Jacobi preconditioner. + * This value is optional. If no value is provided then the default value is chosen by the requested solver. + * + * @see SplinesLinearProblemSparse + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + explicit SplineBuilder3D( + BatchedInterpolationDDom const& batched_interpolation_domain, + std::optional cols_per_chunk = std::nullopt, + std::optional preconditioner_max_block_size = std::nullopt) + : SplineBuilder3D( + interpolation_domain_type(batched_interpolation_domain), + cols_per_chunk, + preconditioner_max_block_size) + { + } + + /// @brief Copy-constructor is deleted. + SplineBuilder3D(SplineBuilder3D const& x) = delete; + + /** + * @brief Move-constructs. + * + * @param x An rvalue to another SplineBuilder2D. + */ + SplineBuilder3D(SplineBuilder3D&& x) = default; + + /// @brief Destructs. + ~SplineBuilder3D() = default; + + /// @brief Copy-assignment is deleted. + SplineBuilder3D& operator=(SplineBuilder3D const& x) = delete; + + /** @brief Move-assigns. + * + * @param x An rvalue to another SplineBuilder. + * @return A reference to this object. + */ + SplineBuilder3D& operator=(SplineBuilder3D&& x) = default; + + /** + * @brief Get the domain for the 2D interpolation mesh used by this class. + * + * This is 2D because it is defined along the dimensions of interest. + * + * @return The 2D domain for the interpolation mesh. + */ + interpolation_domain_type interpolation_domain() const noexcept + { + return ddc::DiscreteDomain< + interpolation_domain_type1, + interpolation_domain_type2, + interpolation_discrete_dimension_type3>( + m_spline_builder1.interpolation_domain(), + m_spline_builder2.interpolation_domain(), + m_spline_builder3.interpolation_domain()); + } + + /** + * @brief Get the whole domain representing interpolation points. + * + * Values of the function must be provided on this domain in order + * to build a spline representation of the function (cartesian product of 2D interpolation_domain and batch_domain). + * + * @param batched_interpolation_domain The whole domain on which the interpolation points are defined. + * + * @return The domain for the interpolation mesh. + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + BatchedInterpolationDDom batched_interpolation_domain( + BatchedInterpolationDDom const& batched_interpolation_domain) const noexcept + { + assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain)); + return batched_interpolation_domain; + } + + /** + * @brief Get the batch domain. + * + * Obtained by removing the dimensions of interest from the whole interpolation domain. + * + * @param batched_interpolation_domain The whole domain on which the interpolation points are defined. + * + * @return The batch domain. + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + batch_domain_type batch_domain( + BatchedInterpolationDDom const& batched_interpolation_domain) const noexcept + { + assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain)); + return ddc::remove_dims_of(batched_interpolation_domain, interpolation_domain()); + } + + /** + * @brief Get the 2D domain on which spline coefficients are defined. + * + * The 2D spline domain corresponding to the dimensions of interest. + * + * @return The 2D domain for the spline coefficients. + */ + ddc::DiscreteDomain spline_domain() + const noexcept + { + return ddc::DiscreteDomain( + ddc::discrete_space().full_domain(), + ddc::discrete_space().full_domain(), + ddc::discrete_space().full_domain()); + } + + /** + * @brief Get the whole domain on which spline coefficients are defined. + * + * Spline approximations (spline-transformed functions) are computed on this domain. + * + * @param batched_interpolation_domain The whole domain on which the interpolation points are defined. + * + * @return The domain for the spline coefficients. + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + batched_spline_domain_type batched_spline_domain( + BatchedInterpolationDDom const& batched_interpolation_domain) const noexcept + { + assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain)); + return ddc::replace_dim_of( + ddc::replace_dim_of( + ddc::replace_dim_of< + interpolation_discrete_dimension_type3, + bsplines_type3>(batched_interpolation_domain, spline_domain()), + spline_domain()), + spline_domain()); + } + + /** + * @brief Compute a 2D spline approximation of a function. + * + * Use the values of a function (defined on + * SplineBuilder2D::batched_interpolation_domain) and the derivatives of the + * function at the boundaries (in the case of BoundCond::HERMITE only) + * to calculate a 2D spline approximation of this function. + * + * The spline approximation is stored as a ChunkSpan of coefficients + * associated with B-splines. + * + * @param[out] spline + * The coefficients of the spline computed by this SplineBuilder. + * @param[in] vals + * The values of the function at the interpolation mesh. + * @param[in] derivs_min1 + * The values of the derivatives at the lower boundary in the first dimension. + * @param[in] derivs_max1 + * The values of the derivatives at the upper boundary in the first dimension. + * @param[in] derivs_min2 + * The values of the derivatives at the lower boundary in the second dimension. + * @param[in] derivs_max2 + * The values of the derivatives at the upper boundary in the second dimension. + * @param[in] mixed_derivs_min1_min2 + * The values of the the cross-derivatives at the lower boundary in the first dimension + * and the lower boundary in the second dimension. + * @param[in] mixed_derivs_max1_min2 + * The values of the the cross-derivatives at the upper boundary in the first dimension + * and the lower boundary in the second dimension. + * @param[in] mixed_derivs_min1_max2 + * The values of the the cross-derivatives at the lower boundary in the first dimension + * and the upper boundary in the second dimension. + * @param[in] mixed_derivs_max1_max2 + * The values of the the cross-derivatives at the upper boundary in the first dimension + * and the upper boundary in the second dimension. + */ + template + void operator()( + ddc::ChunkSpan< + double, + batched_spline_domain_type, + Layout, + memory_space> spline, + ddc::ChunkSpan vals, + std::optional, + Layout, + memory_space>> derivs_min_1_2 + = std::nullopt, + std::optional, + Layout, + memory_space>> derivs_max_1_2 + = std::nullopt, + std::optional, + Layout, + memory_space>> derivs_min_2_3 + = std::nullopt, + std::optional, + Layout, + memory_space>> derivs_max_2_3 + = std::nullopt, + std::optional, + Layout, + memory_space>> derivs_min_1_3 + = std::nullopt, + std::optional, + Layout, + memory_space>> derivs_max_1_3 + = std::nullopt, + std::optional, + Layout, + memory_space>> mixed_derivs_min1_min2_min3 + = std::nullopt, + std::optional, + Layout, + memory_space>> mixed_derivs_max1_min2_min3 + = std::nullopt, + std::optional, + Layout, + memory_space>> mixed_derivs_min1_max2_min3 + = std::nullopt, + std::optional, + Layout, + memory_space>> mixed_derivs_max1_max2_min3 + = std::nullopt, + std::optional, + Layout, + memory_space>> mixed_derivs_min1_min2_max3 + = std::nullopt, + std::optional, + Layout, + memory_space>> mixed_derivs_max1_min2_max3 + = std::nullopt, + std::optional, + Layout, + memory_space>> mixed_derivs_min1_max2_max3 + = std::nullopt, + std::optional, + Layout, + memory_space>> mixed_derivs_max1_max2_max3 + = std::nullopt) const; +}; + + +template < + class ExecSpace, + class MemorySpace, + class BSpline1, + class BSpline2, + class BSpline3, + class DDimI1, + class DDimI2, + class DDimI3, + ddc::BoundCond BcLower1, + ddc::BoundCond BcUpper1, + ddc::BoundCond BcLower2, + ddc::BoundCond BcUpper2, + ddc::BoundCond BcLower3, + ddc::BoundCond BcUpper3, + ddc::SplineSolver Solver> +template +void SplineBuilder3D< + ExecSpace, + MemorySpace, + BSpline1, + BSpline2, + BSpline3, + DDimI1, + DDimI2, + DDimI3, + BcLower1, + BcUpper1, + BcLower2, + BcUpper2, + BcLower3, + BcUpper3, + Solver>:: +operator()( + ddc::ChunkSpan< + double, + batched_spline_domain_type, + Layout, + memory_space> spline, + ddc::ChunkSpan vals, + std::optional, + Layout, + memory_space>> derivs_min_1_2, + std::optional, + Layout, + memory_space>> derivs_max_1_2, + std::optional, + Layout, + memory_space>> derivs_min_2_3, + std::optional, + Layout, + memory_space>> derivs_max_2_3, + std::optional, + Layout, + memory_space>> derivs_min_1_3, + std::optional, + Layout, + memory_space>> derivs_max_1_3, + std::optional, + Layout, + memory_space>> mixed_derivs_min1_min2_min3, + std::optional, + Layout, + memory_space>> mixed_derivs_max1_min2_min3, + std::optional, + Layout, + memory_space>> mixed_derivs_min1_max2_min3, + std::optional, + Layout, + memory_space>> mixed_derivs_max1_max2_min3, + std::optional, + Layout, + memory_space>> mixed_derivs_min1_min2_max3, + std::optional, + Layout, + memory_space>> mixed_derivs_max1_min2_max3, + std::optional, + Layout, + memory_space>> mixed_derivs_min1_max2_max3, + std::optional, + Layout, + memory_space>> mixed_derivs_max1_max2_max3) const +{ + auto const batched_interpolation_domain = vals.domain(); + + assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain)); + + // Spline1-approximate vals (to spline1) + ddc::Chunk spline1_alloc( + m_spline_builder1.batched_spline_domain(batched_interpolation_domain), + ddc::KokkosAllocator()); + ddc::ChunkSpan const spline1 = spline1_alloc.span_view(); + + m_spline_builder1(spline1, vals); + + // Spline2-approximate spline1 (to spline2) + ddc::Chunk spline2_alloc( + m_spline_builder2.batched_spline_domain(spline1.domain()), + ddc::KokkosAllocator()); + ddc::ChunkSpan const spline2 = spline2_alloc.span_view(); + + m_spline_builder2(spline2, spline1.span_cview()); + + // Spline3-approximate spline2 + m_spline_builder3(spline, spline2.span_cview()); +} + +} // namespace ddc diff --git a/include/ddc/kernels/splines/spline_evaluator_3d.hpp b/include/ddc/kernels/splines/spline_evaluator_3d.hpp new file mode 100644 index 000000000..de0692439 --- /dev/null +++ b/include/ddc/kernels/splines/spline_evaluator_3d.hpp @@ -0,0 +1,2425 @@ +// Copyright (C) The DDC development team, see COPYRIGHT.md file +// +// SPDX-License-Identifier: MIT + +#pragma once + +#include +#include +#include + +#include + +#include + +#include "integrals.hpp" +#include "periodic_extrapolation_rule.hpp" + +namespace ddc { + +/** + * @brief A class to evaluate, differentiate or integrate a 3D spline function. + * + * A class which contains an operator () which can be used to evaluate, differentiate or integrate a 3D spline function. + * + * @tparam ExecSpace The Kokkos execution space on which the spline evaluation is performed. + * @tparam MemorySpace The Kokkos memory space on which the data (spline coefficients and evaluation) is stored. + * @tparam BSplines1 The discrete dimension representing the B-splines along the first dimension of interest. + * @tparam BSplines2 The discrete dimension representing the B-splines along the second dimension of interest. + * @tparam BSplines3 The discrete dimension representing the B-splines along the third dimension of interest. + * @tparam EvaluationDDim1 The first discrete dimension on which evaluation points are defined. + * @tparam EvaluationDDim2 The second discrete dimension on which evaluation points are defined. + * @tparam EvaluationDDim3 The third discrete dimension on which evaluation points are defined. + * @tparam LowerExtrapolationRule1 The lower extrapolation rule type along first dimension of interest. + * @tparam UpperExtrapolationRule1 The upper extrapolation rule type along first dimension of interest. + * @tparam LowerExtrapolationRule2 The lower extrapolation rule type along second dimension of interest. + * @tparam UpperExtrapolationRule2 The upper extrapolation rule type along second dimension of interest. + * @tparam LowerExtrapolationRule3 The lower extrapolation rule type along third dimension of interest. + * @tparam UpperExtrapolationRule3 The upper extrapolation rule type along third dimension of interest. + */ +template < + class ExecSpace, + class MemorySpace, + class BSplines1, + class BSplines2, + class BSplines3, + class EvaluationDDim1, + class EvaluationDDim2, + class EvaluationDDim3, + class LowerExtrapolationRule1, + class UpperExtrapolationRule1, + class LowerExtrapolationRule2, + class UpperExtrapolationRule2, + class LowerExtrapolationRule3, + class UpperExtrapolationRule3> +class SplineEvaluator3D +{ +private: + /** + * @brief Tag to indicate that the value of the spline should be evaluated. + */ + struct eval_type + { + }; + + /** + * @brief Tag to indicate that derivative of the spline should be evaluated. + */ + struct eval_deriv_type + { + }; + +public: + /// @brief The type of the first evaluation continuous dimension used by this class. + using continuous_dimension_type1 = typename BSplines1::continuous_dimension_type; + + /// @brief The type of the second evaluation continuous dimension used by this class. + using continuous_dimension_type2 = typename BSplines2::continuous_dimension_type; + + /// @brief The type of the third evaluation continuous dimension used by this class. + using continuous_dimension_type3 = typename BSplines3::continuous_dimension_type; + + /// @brief The type of the Kokkos execution space used by this class. + using exec_space = ExecSpace; + + /// @brief The type of the Kokkos memory space used by this class. + using memory_space = MemorySpace; + + /// @brief The type of the first discrete dimension of interest used by this class. + using evaluation_discrete_dimension_type1 = EvaluationDDim1; + + /// @brief The type of the second discrete dimension of interest used by this class. + using evaluation_discrete_dimension_type2 = EvaluationDDim2; + + /// @brief The type of the third discrete dimension of interest used by this class. + using evaluation_discrete_dimension_type3 = EvaluationDDim3; + + /// @brief The discrete dimension representing the B-splines along first dimension. + using bsplines_type1 = BSplines1; + + /// @brief The discrete dimension representing the B-splines along second dimension. + using bsplines_type2 = BSplines2; + + /// @brief The discrete dimension representing the B-splines along third dimension. + using bsplines_type3 = BSplines3; + + /// @brief The type of the domain for the 1D evaluation mesh along first dimension used by this class. + using evaluation_domain_type1 = ddc::DiscreteDomain; + + /// @brief The type of the domain for the 1D evaluation mesh along second dimension used by this class. + using evaluation_domain_type2 = ddc::DiscreteDomain; + + /// @brief The type of the domain for the 1D evaluation mesh along third dimension used by this class. + using evaluation_domain_type3 = ddc::DiscreteDomain; + + /// @brief The type of the domain for the 3D evaluation mesh used by this class. + using evaluation_domain_type = ddc::DiscreteDomain< + evaluation_discrete_dimension_type1, + evaluation_discrete_dimension_type2, + evaluation_discrete_dimension_type3>; + + /** + * @brief The type of the whole domain representing evaluation points. + * + * @tparam The batched discrete domain on which the interpolation points are defined. + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using batched_evaluation_domain_type = BatchedInterpolationDDom; + + /// @brief The type of the 1D spline domain corresponding to the first dimension of interest. + using spline_domain_type1 = ddc::DiscreteDomain; + + /// @brief The type of the 1D spline domain corresponding to the second dimension of interest. + using spline_domain_type2 = ddc::DiscreteDomain; + + /// @brief The type of the 1D spline domain corresponding to the third dimension of interest. + using spline_domain_type3 = ddc::DiscreteDomain; + + /// @brief The type of the 3D spline domain corresponding to the dimensions of interest. + using spline_domain_type = ddc::DiscreteDomain; + + /** + * @brief The type of the batch domain (obtained by removing the dimensions of interest + * from the whole domain). + * + * @tparam The batched discrete domain on which the interpolation points are defined. + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using batch_domain_type = typename ddc::remove_dims_of_t< + BatchedInterpolationDDom, + evaluation_discrete_dimension_type1, + evaluation_discrete_dimension_type2, + evaluation_discrete_dimension_type3>; + + /** + * @brief The type of the whole spline domain (cartesian product of 3D spline domain + * and batch domain) preserving the underlying memory layout (order of dimensions). + * + * @tparam The batched discrete domain on which the interpolation points are defined. + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using batched_spline_domain_type = + typename ddc::detail::convert_type_seq_to_discrete_domain_t, + ddc::detail::TypeSeq< + evaluation_discrete_dimension_type1, + evaluation_discrete_dimension_type2, + evaluation_discrete_dimension_type3>, + ddc::detail::TypeSeq>>; + + /// @brief The type of the extrapolation rule at the lower boundary along the first dimension. + using lower_extrapolation_rule_1_type = LowerExtrapolationRule1; + + /// @brief The type of the extrapolation rule at the upper boundary along the first dimension. + using upper_extrapolation_rule_1_type = UpperExtrapolationRule1; + + /// @brief The type of the extrapolation rule at the lower boundary along the second dimension. + using lower_extrapolation_rule_2_type = LowerExtrapolationRule2; + + /// @brief The type of the extrapolation rule at the upper boundary along the second dimension. + using upper_extrapolation_rule_2_type = UpperExtrapolationRule2; + + /// @brief The type of the extrapolation rule at the lower boundary along the third dimension. + using lower_extrapolation_rule_3_type = LowerExtrapolationRule3; + + /// @brief The type of the extrapolation rule at the upper boundary along the third dimension. + using upper_extrapolation_rule_3_type = UpperExtrapolationRule3; + +private: + LowerExtrapolationRule1 m_lower_extrap_rule_1; + + UpperExtrapolationRule1 m_upper_extrap_rule_1; + + LowerExtrapolationRule2 m_lower_extrap_rule_2; + + UpperExtrapolationRule2 m_upper_extrap_rule_2; + + LowerExtrapolationRule3 m_lower_extrap_rule_3; + + UpperExtrapolationRule3 m_upper_extrap_rule_3; + +public: + static_assert( + std::is_same_v> + == bsplines_type1::is_periodic() + && std::is_same_v< + UpperExtrapolationRule1, + typename ddc::PeriodicExtrapolationRule> + == bsplines_type1::is_periodic() + && std::is_same_v< + LowerExtrapolationRule2, + typename ddc::PeriodicExtrapolationRule> + == bsplines_type2::is_periodic() + && std::is_same_v< + UpperExtrapolationRule2, + typename ddc::PeriodicExtrapolationRule> + == bsplines_type2::is_periodic() + && std::is_same_v< + LowerExtrapolationRule3, + typename ddc::PeriodicExtrapolationRule> + == bsplines_type3::is_periodic() + && std::is_same_v< + UpperExtrapolationRule3, + typename ddc::PeriodicExtrapolationRule> + == bsplines_type3::is_periodic(), + "PeriodicExtrapolationRule has to be used if and only if dimension is periodic"); + static_assert( + std::is_invocable_r_v< + double, + LowerExtrapolationRule1, + ddc::Coordinate, + ddc::ChunkSpan< + double const, + spline_domain_type, + Kokkos::layout_right, + memory_space>>, + "LowerExtrapolationRule1::operator() has to be callable " + "with usual arguments."); + static_assert( + std::is_invocable_r_v< + double, + UpperExtrapolationRule1, + ddc::Coordinate, + ddc::ChunkSpan< + double const, + spline_domain_type, + Kokkos::layout_right, + memory_space>>, + "UpperExtrapolationRule1::operator() has to be callable " + "with usual arguments."); + static_assert( + std::is_invocable_r_v< + double, + LowerExtrapolationRule2, + ddc::Coordinate, + ddc::ChunkSpan< + double const, + spline_domain_type, + Kokkos::layout_right, + memory_space>>, + "LowerExtrapolationRule2::operator() has to be callable " + "with usual arguments."); + static_assert( + std::is_invocable_r_v< + double, + UpperExtrapolationRule2, + ddc::Coordinate, + ddc::ChunkSpan< + double const, + spline_domain_type, + Kokkos::layout_right, + memory_space>>, + "UpperExtrapolationRule2::operator() has to be callable " + "with usual arguments."); + static_assert( + std::is_invocable_r_v< + double, + LowerExtrapolationRule3, + ddc::Coordinate, + ddc::ChunkSpan< + double const, + spline_domain_type, + Kokkos::layout_right, + memory_space>>, + "LowerExtrapolationRule3::operator() has to be callable " + "with usual arguments."); + static_assert( + std::is_invocable_r_v< + double, + UpperExtrapolationRule3, + ddc::Coordinate, + ddc::ChunkSpan< + double const, + spline_domain_type, + Kokkos::layout_right, + memory_space>>, + "UpperExtrapolationRule3::operator() has to be callable " + "with usual arguments."); + + /** + * @brief Build a SplineEvaluator3D acting on batched_spline_domain. + * + * @param lower_extrap_rule1 The extrapolation rule at the lower boundary along the first dimension. + * @param upper_extrap_rule1 The extrapolation rule at the upper boundary along the first dimension. + * @param lower_extrap_rule2 The extrapolation rule at the lower boundary along the second dimension. + * @param upper_extrap_rule2 The extrapolation rule at the upper boundary along the second dimension. + * @param lower_extrap_rule2 The extrapolation rule at the lower boundary along the third dimension. + * @param upper_extrap_rule2 The extrapolation rule at the upper boundary along the third dimension. + * + * @see NullExtrapolationRule ConstantExtrapolationRule PeriodicExtrapolationRule + */ + explicit SplineEvaluator3D( + LowerExtrapolationRule1 const& lower_extrap_rule1, + UpperExtrapolationRule1 const& upper_extrap_rule1, + LowerExtrapolationRule2 const& lower_extrap_rule2, + UpperExtrapolationRule2 const& upper_extrap_rule2, + LowerExtrapolationRule3 const& lower_extrap_rule3, + UpperExtrapolationRule3 const& upper_extrap_rule3) + : m_lower_extrap_rule_1(lower_extrap_rule1) + , m_upper_extrap_rule_1(upper_extrap_rule1) + , m_lower_extrap_rule_2(lower_extrap_rule2) + , m_upper_extrap_rule_2(upper_extrap_rule2) + , m_lower_extrap_rule_3(lower_extrap_rule3) + , m_upper_extrap_rule_3(upper_extrap_rule3) + { + } + + /** + * @brief Copy-constructs. + * + * @param x A reference to another SplineEvaluator. + */ + SplineEvaluator3D(SplineEvaluator3D const& x) = default; + + /** + * @brief Move-constructs. + * + * @param x An rvalue to another SplineEvaluator. + */ + SplineEvaluator3D(SplineEvaluator3D&& x) = default; + + /// @brief Destructs. + ~SplineEvaluator3D() = default; + + /** + * @brief Copy-assigns. + * + * @param x A reference to another SplineEvaluator. + * @return A reference to this object. + */ + SplineEvaluator3D& operator=(SplineEvaluator3D const& x) = default; + + /** + * @brief Move-assigns. + * + * @param x An rvalue to another SplineEvaluator. + * @return A reference to this object. + */ + SplineEvaluator3D& operator=(SplineEvaluator3D&& x) = default; + + /** + * @brief Get the lower extrapolation rule along the first dimension. + * + * Extrapolation rules are functors used to define the behavior of the SplineEvaluator out of the domain where the break points of the B-splines are defined. + * + * @return The lower extrapolation rule along the first dimension. + * + * @see NullExtrapolationRule ConstantExtrapolationRule PeriodicExtrapolationRule + */ + lower_extrapolation_rule_1_type lower_extrapolation_rule_dim_1() const + { + return m_lower_extrap_rule_1; + } + + /** + * @brief Get the upper extrapolation rule along the first dimension. + * + * Extrapolation rules are functors used to define the behavior of the SplineEvaluator out of the domain where the break points of the B-splines are defined. + * + * @return The upper extrapolation rule along the first dimension. + * + * @see NullExtrapolationRule ConstantExtrapolationRule PeriodicExtrapolationRule + */ + upper_extrapolation_rule_1_type upper_extrapolation_rule_dim_1() const + { + return m_upper_extrap_rule_1; + } + + /** + * @brief Get the lower extrapolation rule along the second dimension. + * + * Extrapolation rules are functors used to define the behavior of the SplineEvaluator out of the domain where the break points of the B-splines are defined. + * + * @return The lower extrapolation rule along the second dimension. + * + * @see NullExtrapolationRule ConstantExtrapolationRule PeriodicExtrapolationRule + */ + lower_extrapolation_rule_2_type lower_extrapolation_rule_dim_2() const + { + return m_lower_extrap_rule_2; + } + + /** + * @brief Get the upper extrapolation rule along the second dimension. + * + * Extrapolation rules are functors used to define the behavior of the SplineEvaluator out of the domain where the break points of the B-splines are defined. + * + * @return The upper extrapolation rule along the second dimension. + * + * @see NullExtrapolationRule ConstantExtrapolationRule PeriodicExtrapolationRule + */ + upper_extrapolation_rule_2_type upper_extrapolation_rule_dim_2() const + { + return m_upper_extrap_rule_2; + } + + /** + * @brief Get the lower extrapolation rule along the third dimension. + * + * Extrapolation rules are functors used to define the behavior of the SplineEvaluator out of the domain where the break points of the B-splines are defined. + * + * @return The lower extrapolation rule along the third dimension. + * + * @see NullExtrapolationRule ConstantExtrapolationRule PeriodicExtrapolationRule + */ + lower_extrapolation_rule_3_type lower_extrapolation_rule_dim_3() const + { + return m_lower_extrap_rule_3; + } + + /** + * @brief Get the upper extrapolation rule along the third dimension. + * + * Extrapolation rules are functors used to define the behavior of the SplineEvaluator out of the domain where the break points of the B-splines are defined. + * + * @return The upper extrapolation rule along the third dimension. + * + * @see NullExtrapolationRule ConstantExtrapolationRule PeriodicExtrapolationRule + */ + upper_extrapolation_rule_3_type upper_extrapolation_rule_dim_3() const + { + return m_upper_extrap_rule_3; + } + + /** + * @brief Evaluate 3D spline function (described by its spline coefficients) at a given coordinate. + * + * The spline coefficients represent a 3D spline function defined on a B-splines (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * Remark: calling SplineBuilder3D then SplineEvaluator3D corresponds to a 3D spline interpolation. + * + * @param coord_eval The coordinate where the spline is evaluated. Note that only the components along the dimensions of interest are used. + * @param spline_coef A ChunkSpan storing the 3D spline coefficients. + * + * @return The value of the spline function at the desired coordinate. + */ + template + KOKKOS_FUNCTION double operator()( + ddc::Coordinate const& coord_eval, + ddc::ChunkSpan const + spline_coef) const + { + return eval(coord_eval, spline_coef); + } + + /** + * @brief Evaluate 3D spline function (described by its spline coefficients) on a mesh. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a nD evaluation. This is a batched 3D evaluation. This means that for each slice of coordinates + * identified by a batch_domain_type::discrete_element_type, the evaluation is performed with the 3D set of + * spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * Remark: calling SplineBuilder3D then SplineEvaluator3D corresponds to a 3D spline interpolation. + * + * @param[out] spline_eval The values of the 3D spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. + * @param[in] coords_eval The coordinates where the spline is evaluated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 3D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template < + class Layout1, + class Layout2, + class Layout3, + class BatchedInterpolationDDom, + class... CoordsDims> + void operator()( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + ddc::Coordinate const, + BatchedInterpolationDDom, + Layout2, + memory_space> const coords_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const + { + batch_domain_type const batch_domain(coords_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); + evaluation_domain_type3 const evaluation_domain3(spline_eval.domain()); + ddc::parallel_for_each( + "ddc_splines_evaluate_3d", + exec_space(), + batch_domain, + KOKKOS_CLASS_LAMBDA( + typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { + auto const spline_eval_3D = spline_eval[j]; + auto const coords_eval_3D = coords_eval[j]; + auto const spline_coef_3D = spline_coef[j]; + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { + for (auto const i3 : evaluation_domain3) { + spline_eval_3D(i1, i2, i3) + = eval(coords_eval_3D(i1, i2, i3), spline_coef_3D); + } + } + } + }); + } + + /** + * @brief Evaluate 3D spline function (described by its spline coefficients) on a mesh. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a multidimensional evaluation. This is a batched 3D evaluation. + * This means that for each slice of spline_eval the evaluation is performed with + * the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * Remark: calling SplineBuilder3D then SplineEvaluator3D corresponds to a 3D spline interpolation. + * + * @param[out] spline_eval The values of the 3D spline function at their coordinates. + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template + void operator()( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const + { + batch_domain_type const batch_domain(spline_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); + evaluation_domain_type3 const evaluation_domain3(spline_eval.domain()); + ddc::parallel_for_each( + "ddc_splines_evaluate_3d", + exec_space(), + batch_domain, + KOKKOS_CLASS_LAMBDA( + typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { + auto const spline_eval_3D = spline_eval[j]; + auto const spline_coef_3D = spline_coef[j]; + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { + for (auto const i3 : evaluation_domain3) { + ddc::Coordinate< + continuous_dimension_type1, + continuous_dimension_type2, + continuous_dimension_type3> + coord_eval_3D( + ddc::coordinate(i1), + ddc::coordinate(i2), + ddc::coordinate(i3)); + spline_eval_3D(i1, i2, i3) + = eval(coord_eval_3D(i1, i2, i3), spline_coef_3D); + } + } + } + }); + } + + /** + * @brief Differentiate 3D spline function (described by its spline coefficients) at a given coordinate along first dimension of interest. + * + * The spline coefficients represent a 3D spline function defined on a B-splines (basis splines). They can be + * obtained via various methods, such as using a SplineBuilder3D. + * + * @param coord_eval The coordinate where the spline is differentiated. Note that only the components along the dimensions of interest are used. + * @param spline_coef A ChunkSpan storing the 3D spline coefficients. + * + * @return The derivative of the spline function at the desired coordinate. + */ + template + KOKKOS_FUNCTION double deriv_dim_1( + ddc::Coordinate const& coord_eval, + ddc::ChunkSpan const + spline_coef) const + { + return eval_no_bc(coord_eval, spline_coef); + } + + /** + * @brief Differentiate 3D spline function (described by its spline coefficients) at a given coordinate along second dimension of interest. + * + * The spline coefficients represent a 3D spline function defined on a B-splines (basis splines). They can be + * obtained via various methods, such as using a SplineBuilder3D. + * + * @param coord_eval The coordinate where the spline is differentiated. Note that only the components along the dimensions of interest are used. + * @param spline_coef A ChunkSpan storing the 3D spline coefficients. + * + * @return The derivative of the spline function at the desired coordinate. + */ + template + KOKKOS_FUNCTION double deriv_dim_2( + ddc::Coordinate const& coord_eval, + ddc::ChunkSpan const + spline_coef) const + { + return eval_no_bc(coord_eval, spline_coef); + } + + /** + * @brief Differentiate 3D spline function (described by its spline coefficients) at a given coordinate along third dimension of interest. + * + * The spline coefficients represent a 3D spline function defined on a B-splines (basis splines). They can be + * obtained via various methods, such as using a SplineBuilder3D. + * + * @param coord_eval The coordinate where the spline is differentiated. Note that only the components along the dimensions of interest are used. + * @param spline_coef A ChunkSpan storing the 3D spline coefficients. + * + * @return The derivative of the spline function at the desired coordinate. + */ + template + KOKKOS_FUNCTION double deriv_dim_3( + ddc::Coordinate const& coord_eval, + ddc::ChunkSpan const + spline_coef) const + { + return eval_no_bc(coord_eval, spline_coef); + } + + /** + * @brief Cross-differentiate 3D spline function (described by its spline coefficients) at a given coordinate along the first and second dimensions. + * + * The spline coefficients represent a 3D spline function defined on a B-splines (basis splines). They can be + * obtained via various methods, such as using a SplineBuilder3D. + * + * @param coord_eval The coordinate where the spline is differentiated. Note that only the components along the dimensions of interest are used. + * @param spline_coef A ChunkSpan storing the 3D spline coefficients. + * + * @return The derivative of the spline function at the desired coordinate. + */ + template + KOKKOS_FUNCTION double deriv_1_and_2( + ddc::Coordinate const& coord_eval, + ddc::ChunkSpan const + spline_coef) const + { + return eval_no_bc(coord_eval, spline_coef); + } + + /** + * @brief Cross-differentiate 3D spline function (described by its spline coefficients) at a given coordinate along the second and third dimensions. + * + * The spline coefficients represent a 3D spline function defined on a B-splines (basis splines). They can be + * obtained via various methods, such as using a SplineBuilder3D. + * + * @param coord_eval The coordinate where the spline is differentiated. Note that only the components along the dimensions of interest are used. + * @param spline_coef A ChunkSpan storing the 3D spline coefficients. + * + * @return The derivative of the spline function at the desired coordinate. + */ + template + KOKKOS_FUNCTION double deriv_2_and_3( + ddc::Coordinate const& coord_eval, + ddc::ChunkSpan const + spline_coef) const + { + return eval_no_bc(coord_eval, spline_coef); + } + + /** + * @brief Cross-differentiate 3D spline function (described by its spline coefficients) at a given coordinate along the first and third dimensions. + * + * The spline coefficients represent a 3D spline function defined on a B-splines (basis splines). They can be + * obtained via various methods, such as using a SplineBuilder3D. + * + * @param coord_eval The coordinate where the spline is differentiated. Note that only the components along the dimensions of interest are used. + * @param spline_coef A ChunkSpan storing the 3D spline coefficients. + * + * @return The derivative of the spline function at the desired coordinate. + */ + template + KOKKOS_FUNCTION double deriv_1_and_3( + ddc::Coordinate const& coord_eval, + ddc::ChunkSpan const + spline_coef) const + { + return eval_no_bc(coord_eval, spline_coef); + } + + /** + * @brief Cross-differentiate 3D spline function (described by its spline coefficients) at a given coordinate along the dimensions of interest. + * + * The spline coefficients represent a 3D spline function defined on a B-splines (basis splines). They can be + * obtained via various methods, such as using a SplineBuilder3D. + * + * @param coord_eval The coordinate where the spline is differentiated. Note that only the components along the dimensions of interest are used. + * @param spline_coef A ChunkSpan storing the 3D spline coefficients. + * + * @return The derivative of the spline function at the desired coordinate. + */ + template + KOKKOS_FUNCTION double deriv_1_2_3( + ddc::Coordinate const& coord_eval, + ddc::ChunkSpan const + spline_coef) const + { + return eval_no_bc< + eval_deriv_type, + eval_deriv_type, + eval_deriv_type>(coord_eval, spline_coef); + } + + /** + * @brief Differentiate 3D spline function (described by its spline coefficients) at a given coordinate along a specified dimension of interest. + * + * The spline coefficients represent a 3D spline function defined on a B-splines (basis splines). They can be + * obtained via various methods, such as using a SplineBuilder3D. + * + * @tparam InterestDim Dimension along which differentiation is performed. + * + * @param coord_eval The coordinate where the spline is differentiated. Note that only the components along the dimensions of interest are used. + * @param spline_coef A ChunkSpan storing the 3D spline coefficients. + * + * @return The derivative of the spline function at the desired coordinate. + */ + template + KOKKOS_FUNCTION double deriv( + ddc::Coordinate const& coord_eval, + ddc::ChunkSpan const + spline_coef) const + { + static_assert( + std::is_same_v + || std::is_same_v + || std::is_same_v); + if constexpr (std::is_same_v< + InterestDim, + typename evaluation_discrete_dimension_type1:: + continuous_dimension_type>) { + return deriv_dim_1(coord_eval, spline_coef); + } else if constexpr (std::is_same_v< + InterestDim, + typename evaluation_discrete_dimension_type2:: + continuous_dimension_type>) { + return deriv_dim_2(coord_eval, spline_coef); + } else if constexpr (std::is_same_v< + InterestDim, + typename evaluation_discrete_dimension_type3:: + continuous_dimension_type>) { + return deriv_dim_3(coord_eval, spline_coef); + } + } + + /** + * @brief Double-differentiate 3D spline function (described by its spline coefficients) at a given coordinate along specified dimensions of interest. + * + * The spline coefficients represent a 3D spline function defined on a B-splines (basis splines). They can be + * obtained via various methods, such as using a SplineBuilder3D. + * + * Note: double-differentiation other than cross-differentiation is not supported atm. See #440 + * + * @tparam InterestDim1 First dimension along which differentiation is performed. + * @tparam InterestDim2 Second dimension along which differentiation is performed. + * + * @param coord_eval The coordinate where the spline is double-differentiated. Note that only the components along the dimensions of interest are used. + * @param spline_coef A ChunkSpan storing the 3D spline coefficients. + * + * @return The derivative of the spline function at the desired coordinate. + */ + template + KOKKOS_FUNCTION double deriv2( + ddc::Coordinate const& coord_eval, + ddc::ChunkSpan const + spline_coef) const + { + static_assert( + (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type2::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type2::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v)); + + if constexpr ( + (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v)) { + return deriv_1_and_2(coord_eval, spline_coef); + } else if constexpr ( + (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type2::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type2::continuous_dimension_type> + && std::is_same_v)) { + return deriv_2_and_3(coord_eval, spline_coef); + } else if constexpr ( + (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v)) { + return deriv_1_and_3(coord_eval, spline_coef); + } + } + + /** + * @brief Triple-differentiate 3D spline function (described by its spline coefficients) at a given coordinate along specified dimensions of interest. + * + * The spline coefficients represent a 3D spline function defined on a B-splines (basis splines). They can be + * obtained via various methods, such as using a SplineBuilder3D. + * + * Note: triple-differentiation other than cross-differentiation is not supported atm. See #440 + * + * @tparam InterestDim1 First dimension along which differentiation is performed. + * @tparam InterestDim2 Second dimension along which differentiation is performed. + * @tparam InterestDim3 Third dimension along which differentiation is performed. + * + * @param coord_eval The coordinate where the spline is triple-differentiated. Note that only the components along the dimensions of interest are used. + * @param spline_coef A ChunkSpan storing the 3D spline coefficients. + * + * @return The derivative of the spline function at the desired coordinate. + */ + template < + class InterestDim1, + class InterestDim2, + class InterestDim3, + class Layout, + class... CoordsDims> + KOKKOS_FUNCTION double deriv3( + ddc::Coordinate const& coord_eval, + ddc::ChunkSpan const + spline_coef) const + { + static_assert( + (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v + && std::is_same_v) + || (std::is_same_v< + InterestDim3, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v + && std::is_same_v)); + + return deriv_1_2_3(coord_eval, spline_coef); + } + + /** + * @brief Differentiate 3D spline function (described by its spline coefficients) on a mesh along first dimension of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a nD evaluation. This is a batched 3D differentiation. + * This means that for each slice of coordinates identified by a batch_domain_type::discrete_element_type, + * the differentiation is performed with the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @param[out] spline_eval The derivatives of the 3D spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 3D spline coefficients retained to perform the evaluation). + * @param[in] coords_eval The coordinates where the spline is differentiated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 3D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template < + class Layout1, + class Layout2, + class Layout3, + class BatchedInterpolationDDom, + class... CoordsDims> + void deriv_dim_1( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + ddc::Coordinate const, + BatchedInterpolationDDom, + Layout2, + memory_space> const coords_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const + { + batch_domain_type const batch_domain(coords_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); + evaluation_domain_type3 const evaluation_domain3(spline_eval.domain()); + ddc::parallel_for_each( + "ddc_splines_differentiate_3d_dim_1", + exec_space(), + batch_domain, + KOKKOS_CLASS_LAMBDA( + typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { + auto const spline_eval_3D = spline_eval[j]; + auto const coords_eval_3D = coords_eval[j]; + auto const spline_coef_3D = spline_coef[j]; + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { + for (auto const i3 : evaluation_domain3) { + spline_eval_3D(i1, i2, i3) = eval_no_bc< + eval_deriv_type, + eval_type, + eval_type>(coords_eval_3D(i1, i2, i3), spline_coef_3D); + } + } + } + }); + } + + /** + * @brief Differentiate 3D spline function (described by its spline coefficients) on a mesh along first dimension of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a multidimensional evaluation. This is a batched 3D evaluation. + * This means that for each slice of spline_eval the evaluation is performed with + * the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @param[out] spline_eval The derivatives of the 3D spline function at the desired coordinates. + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template + void deriv_dim_1( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const + { + batch_domain_type const batch_domain(spline_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); + evaluation_domain_type3 const evaluation_domain3(spline_eval.domain()); + ddc::parallel_for_each( + "ddc_splines_differentiate_3d_dim_1", + exec_space(), + batch_domain, + KOKKOS_CLASS_LAMBDA( + typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { + auto const spline_eval_3D = spline_eval[j]; + auto const spline_coef_3D = spline_coef[j]; + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { + for (auto const i3 : evaluation_domain3) { + ddc::Coordinate< + continuous_dimension_type1, + continuous_dimension_type2, + continuous_dimension_type3> + coord_eval_3D( + ddc::coordinate(i1), + ddc::coordinate(i2), + ddc::coordinate(i3)); + spline_eval_3D(i1, i2, i3) = eval_no_bc< + eval_deriv_type, + eval_type, + eval_type>(coord_eval_3D, spline_coef_3D); + } + } + } + }); + } + + /** + * @brief Differentiate 3D spline function (described by its spline coefficients) on a mesh along second dimension of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a nD differentiation. This is a batched 3D differentiation. + * This means that for each slice of coordinates identified by a batch_domain_type::discrete_element_type, + * the differentiation is performed with the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @param[out] spline_eval The derivatives of the 3D spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. + * @param[in] coords_eval The coordinates where the spline is differentiated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 3D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template < + class Layout1, + class Layout2, + class Layout3, + class BatchedInterpolationDDom, + class... CoordsDims> + void deriv_dim_2( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + ddc::Coordinate const, + BatchedInterpolationDDom, + Layout2, + memory_space> const coords_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const + { + batch_domain_type const batch_domain(coords_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); + evaluation_domain_type3 const evaluation_domain3(spline_eval.domain()); + ddc::parallel_for_each( + "ddc_splines_differentiate_3d_dim_2", + exec_space(), + batch_domain, + KOKKOS_CLASS_LAMBDA( + typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { + auto const spline_eval_3D = spline_eval[j]; + auto const coords_eval_3D = coords_eval[j]; + auto const spline_coef_3D = spline_coef[j]; + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { + for (auto const i3 : evaluation_domain3) { + spline_eval_3D(i1, i2, i3) = eval_no_bc< + eval_type, + eval_deriv_type, + eval_type>(coords_eval_3D(i1, i2, i3), spline_coef_3D); + } + } + } + }); + } + + /** + * @brief Differentiate 3D spline function (described by its spline coefficients) on a mesh along second dimension of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a multidimensional evaluation. This is a batched 3D evaluation. + * This means that for each slice of spline_eval the evaluation is performed with + * the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @param[out] spline_eval The derivatives of the 3D spline function at the desired coordinates. + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template + void deriv_dim_2( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const + { + batch_domain_type const batch_domain(spline_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); + evaluation_domain_type3 const evaluation_domain3(spline_eval.domain()); + ddc::parallel_for_each( + "ddc_splines_differentiate_3d_dim_2", + exec_space(), + batch_domain, + KOKKOS_CLASS_LAMBDA( + typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { + auto const spline_eval_3D = spline_eval[j]; + auto const spline_coef_3D = spline_coef[j]; + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { + for (auto const i3 : evaluation_domain3) { + ddc::Coordinate< + continuous_dimension_type1, + continuous_dimension_type2, + continuous_dimension_type3> + coord_eval_3D( + ddc::coordinate(i1), + ddc::coordinate(i2), + ddc::coordinate(i3)); + spline_eval_3D(i1, i2, i3) = eval_no_bc< + eval_type, + eval_deriv_type, + eval_type>(coord_eval_3D, spline_coef_3D); + } + } + } + }); + } + + /** + * @brief Differentiate 3D spline function (described by its spline coefficients) on a mesh along third dimension of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a nD differentiation. This is a batched 3D differentiation. + * This means that for each slice of coordinates identified by a batch_domain_type::discrete_element_type, + * the differentiation is performed with the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @param[out] spline_eval The derivatives of the 3D spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. + * @param[in] coords_eval The coordinates where the spline is differentiated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 3D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template < + class Layout1, + class Layout2, + class Layout3, + class BatchedInterpolationDDom, + class... CoordsDims> + void deriv_dim_3( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + ddc::Coordinate const, + BatchedInterpolationDDom, + Layout2, + memory_space> const coords_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const + { + batch_domain_type const batch_domain(coords_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); + evaluation_domain_type3 const evaluation_domain3(spline_eval.domain()); + ddc::parallel_for_each( + "ddc_splines_differentiate_3d_dim_3", + exec_space(), + batch_domain, + KOKKOS_CLASS_LAMBDA( + typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { + auto const spline_eval_3D = spline_eval[j]; + auto const coords_eval_3D = coords_eval[j]; + auto const spline_coef_3D = spline_coef[j]; + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { + for (auto const i3 : evaluation_domain3) { + spline_eval_3D(i1, i2, i3) + = eval_no_bc( + coords_eval_3D(i1, i2, i3), + spline_coef_3D); + } + } + } + }); + } + + /** + * @brief Differentiate 3D spline function (described by its spline coefficients) on a mesh along third dimension of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a multidimensional evaluation. This is a batched 3D evaluation. + * This means that for each slice of spline_eval the evaluation is performed with + * the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @param[out] spline_eval The derivatives of the 3D spline function at the desired coordinates. + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template + void deriv_dim_3( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const + { + batch_domain_type const batch_domain(spline_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); + evaluation_domain_type3 const evaluation_domain3(spline_eval.domain()); + ddc::parallel_for_each( + "ddc_splines_differentiate_3d_dim_3", + exec_space(), + batch_domain, + KOKKOS_CLASS_LAMBDA( + typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { + auto const spline_eval_3D = spline_eval[j]; + auto const spline_coef_3D = spline_coef[j]; + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { + for (auto const i3 : evaluation_domain3) { + ddc::Coordinate< + continuous_dimension_type1, + continuous_dimension_type2, + continuous_dimension_type3> + coord_eval_3D( + ddc::coordinate(i1), + ddc::coordinate(i2), + ddc::coordinate(i3)); + spline_eval_3D(i1, i2, i3) = eval_no_bc< + eval_type, + eval_type, + eval_deriv_type>(coord_eval_3D, spline_coef_3D); + } + } + } + }); + } + + /** + * @brief Cross-differentiate 3D spline function (described by its spline coefficients) on a mesh along the first and second dimensions of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a nD cross-differentiation. This is a batched 3D cross-differentiation. + * This means that for each slice of coordinates identified by a batch_domain_type::discrete_element_type, + * the cross-differentiation is performed with the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @param[out] spline_eval The cross-derivatives of the 3D spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. + * @param[in] coords_eval The coordinates where the spline is differentiated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 3D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template < + class Layout1, + class Layout2, + class Layout3, + class BatchedInterpolationDDom, + class... CoordsDims> + void deriv_1_and_2( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + ddc::Coordinate const, + BatchedInterpolationDDom, + Layout2, + memory_space> const coords_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const + { + batch_domain_type const batch_domain(coords_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); + evaluation_domain_type3 const evaluation_domain3(spline_eval.domain()); + ddc::parallel_for_each( + "ddc_splines_cross_differentiate_3d_dim1_2", + exec_space(), + batch_domain, + KOKKOS_CLASS_LAMBDA( + typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { + auto const spline_eval_3D = spline_eval[j]; + auto const coords_eval_3D = coords_eval[j]; + auto const spline_coef_3D = spline_coef[j]; + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { + for (auto const i3 : evaluation_domain3) { + spline_eval_3D(i1, i2, i3) = eval_no_bc< + eval_deriv_type, + eval_deriv_type, + eval_type>(coords_eval_3D(i1, i2, i3), spline_coef_3D); + } + } + } + }); + } + + /** + * @brief Cross-differentiate 3D spline function (described by its spline coefficients) on a mesh along the first and second dimensions of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a multidimensional evaluation. This is a batched 3D evaluation. + * This means that for each slice of spline_eval the evaluation is performed with + * the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @param[out] spline_eval The cross-derivatives of the 3D spline function at the desired coordinates. + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template + void deriv_1_and_2( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const + { + batch_domain_type const batch_domain(spline_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); + evaluation_domain_type3 const evaluation_domain3(spline_eval.domain()); + ddc::parallel_for_each( + "ddc_splines_cross_differentiate_3d_dim1_2", + exec_space(), + batch_domain, + KOKKOS_CLASS_LAMBDA( + typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { + auto const spline_eval_3D = spline_eval[j]; + auto const spline_coef_3D = spline_coef[j]; + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { + for (auto const i3 : evaluation_domain3) { + ddc::Coordinate< + continuous_dimension_type1, + continuous_dimension_type2, + continuous_dimension_type3> + coord_eval_3D( + ddc::coordinate(i1), + ddc::coordinate(i2), + ddc::coordinate(i3)); + spline_eval_3D(i1, i2, i3) = eval_no_bc< + eval_deriv_type, + eval_deriv_type, + eval_type>(coord_eval_3D, spline_coef_3D); + } + } + } + }); + } + + /** + * @brief Cross-differentiate 3D spline function (described by its spline coefficients) on a mesh along the second and third dimensions of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a nD cross-differentiation. This is a batched 3D cross-differentiation. + * This means that for each slice of coordinates identified by a batch_domain_type::discrete_element_type, + * the cross-differentiation is performed with the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @param[out] spline_eval The cross-derivatives of the 3D spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. + * @param[in] coords_eval The coordinates where the spline is differentiated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 3D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template < + class Layout1, + class Layout2, + class Layout3, + class BatchedInterpolationDDom, + class... CoordsDims> + void deriv_2_and_3( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + ddc::Coordinate const, + BatchedInterpolationDDom, + Layout2, + memory_space> const coords_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const + { + batch_domain_type const batch_domain(coords_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); + evaluation_domain_type3 const evaluation_domain3(spline_eval.domain()); + ddc::parallel_for_each( + "ddc_splines_cross_differentiate_3d_dim2_3", + exec_space(), + batch_domain, + KOKKOS_CLASS_LAMBDA( + typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { + auto const spline_eval_3D = spline_eval[j]; + auto const coords_eval_3D = coords_eval[j]; + auto const spline_coef_3D = spline_coef[j]; + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { + for (auto const i3 : evaluation_domain3) { + spline_eval_3D(i1, i2, i3) + = eval_no_bc( + coords_eval_3D(i1, i2, i3), + spline_coef_3D); + } + } + } + }); + } + + /** + * @brief Cross-differentiate 3D spline function (described by its spline coefficients) on a mesh along the second and third dimensions of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a multidimensional evaluation. This is a batched 3D evaluation. + * This means that for each slice of spline_eval the evaluation is performed with + * the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @param[out] spline_eval The cross-derivatives of the 3D spline function at the desired coordinates. + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template + void deriv_2_and_3( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const + { + batch_domain_type const batch_domain(spline_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); + evaluation_domain_type3 const evaluation_domain3(spline_eval.domain()); + ddc::parallel_for_each( + "ddc_splines_cross_differentiate_3d_dim2_3", + exec_space(), + batch_domain, + KOKKOS_CLASS_LAMBDA( + typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { + auto const spline_eval_3D = spline_eval[j]; + auto const spline_coef_3D = spline_coef[j]; + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { + for (auto const i3 : evaluation_domain3) { + ddc::Coordinate< + continuous_dimension_type1, + continuous_dimension_type2, + continuous_dimension_type3> + coord_eval_3D( + ddc::coordinate(i1), + ddc::coordinate(i2), + ddc::coordinate(i3)); + spline_eval_3D(i1, i2, i3) = eval_no_bc< + eval_type, + eval_deriv_type, + eval_deriv_type>(coord_eval_3D, spline_coef_3D); + } + } + } + }); + } + + /** + * @brief Cross-differentiate 3D spline function (described by its spline coefficients) on a mesh along the first and third dimensions of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a nD cross-differentiation. This is a batched 3D cross-differentiation. + * This means that for each slice of coordinates identified by a batch_domain_type::discrete_element_type, + * the cross-differentiation is performed with the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @param[out] spline_eval The cross-derivatives of the 3D spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. + * @param[in] coords_eval The coordinates where the spline is differentiated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 3D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template < + class Layout1, + class Layout2, + class Layout3, + class BatchedInterpolationDDom, + class... CoordsDims> + void deriv_1_and_3( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + ddc::Coordinate const, + BatchedInterpolationDDom, + Layout2, + memory_space> const coords_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const + { + batch_domain_type const batch_domain(coords_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); + evaluation_domain_type3 const evaluation_domain3(spline_eval.domain()); + ddc::parallel_for_each( + "ddc_splines_cross_differentiate_3d_dim1_3", + exec_space(), + batch_domain, + KOKKOS_CLASS_LAMBDA( + typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { + auto const spline_eval_3D = spline_eval[j]; + auto const coords_eval_3D = coords_eval[j]; + auto const spline_coef_3D = spline_coef[j]; + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { + for (auto const i3 : evaluation_domain3) { + spline_eval_3D(i1, i2, i3) + = eval_no_bc( + coords_eval_3D(i1, i2, i3), + spline_coef_3D); + } + } + } + }); + } + + /** + * @brief Cross-differentiate 3D spline function (described by its spline coefficients) on a mesh along the first and third dimensions of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a multidimensional evaluation. This is a batched 3D evaluation. + * This means that for each slice of spline_eval the evaluation is performed with + * the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @param[out] spline_eval The cross-derivatives of the 3D spline function at the desired coordinates. + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template + void deriv_1_and_3( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const + { + batch_domain_type const batch_domain(spline_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); + evaluation_domain_type3 const evaluation_domain3(spline_eval.domain()); + ddc::parallel_for_each( + "ddc_splines_cross_differentiate_3d_dim1_3", + exec_space(), + batch_domain, + KOKKOS_CLASS_LAMBDA( + typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { + auto const spline_eval_3D = spline_eval[j]; + auto const spline_coef_3D = spline_coef[j]; + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { + for (auto const i3 : evaluation_domain3) { + ddc::Coordinate< + continuous_dimension_type1, + continuous_dimension_type2, + continuous_dimension_type3> + coord_eval_3D( + ddc::coordinate(i1), + ddc::coordinate(i2), + ddc::coordinate(i3)); + spline_eval_3D(i1, i2, i3) = eval_no_bc< + eval_deriv_type, + eval_type, + eval_deriv_type>(coord_eval_3D, spline_coef_3D); + } + } + } + }); + } + + /** + * @brief Cross-differentiate 3D spline function (described by its spline coefficients) on a mesh along the dimensions of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a nD cross-differentiation. This is a batched 3D cross-differentiation. + * This means that for each slice of coordinates identified by a batch_domain_type::discrete_element_type, + * the cross-differentiation is performed with the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @param[out] spline_eval The cross-derivatives of the 3D spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. + * @param[in] coords_eval The coordinates where the spline is differentiated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 3D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template < + class Layout1, + class Layout2, + class Layout3, + class BatchedInterpolationDDom, + class... CoordsDims> + void deriv_1_2_3( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + ddc::Coordinate const, + BatchedInterpolationDDom, + Layout2, + memory_space> const coords_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const + { + batch_domain_type const batch_domain(coords_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); + evaluation_domain_type3 const evaluation_domain3(spline_eval.domain()); + ddc::parallel_for_each( + "ddc_splines_cross_differentiate_3d_dim1_2_3", + exec_space(), + batch_domain, + KOKKOS_CLASS_LAMBDA( + typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { + auto const spline_eval_3D = spline_eval[j]; + auto const coords_eval_3D = coords_eval[j]; + auto const spline_coef_3D = spline_coef[j]; + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { + for (auto const i3 : evaluation_domain3) { + spline_eval_3D(i1, i2, i3) = eval_no_bc< + eval_deriv_type, + eval_deriv_type, + eval_deriv_type>( + coords_eval_3D(i1, i2, i3), + spline_coef_3D); + } + } + } + }); + } + + /** + * @brief Cross-differentiate 3D spline function (described by its spline coefficients) on a mesh along the dimensions of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a multidimensional evaluation. This is a batched 3D evaluation. + * This means that for each slice of spline_eval the evaluation is performed with + * the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @param[out] spline_eval The cross-derivatives of the 3D spline function at the desired coordinates. + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template + void deriv_1_2_3( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const + { + batch_domain_type const batch_domain(spline_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); + evaluation_domain_type3 const evaluation_domain3(spline_eval.domain()); + ddc::parallel_for_each( + "ddc_splines_cross_differentiate_3d_dim1_2_3", + exec_space(), + batch_domain, + KOKKOS_CLASS_LAMBDA( + typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { + auto const spline_eval_3D = spline_eval[j]; + auto const spline_coef_3D = spline_coef[j]; + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { + for (auto const i3 : evaluation_domain3) { + ddc::Coordinate< + continuous_dimension_type1, + continuous_dimension_type2, + continuous_dimension_type3> + coord_eval_3D( + ddc::coordinate(i1), + ddc::coordinate(i2), + ddc::coordinate(i3)); + spline_eval_3D(i1, i2, i3) = eval_no_bc< + eval_deriv_type, + eval_deriv_type, + eval_deriv_type>(coord_eval_3D, spline_coef_3D); + } + } + } + }); + } + + /** + * @brief Differentiate spline function (described by its spline coefficients) on a mesh along a specified dimension of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a nD evaluation. This is a batched 3D differentiation. + * This means that for each slice of coordinates identified by a batch_domain_type::discrete_element_type, + * the differentiation is performed with the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @tparam InterestDim Dimension along which differentiation is performed. + * @param[out] spline_eval The derivatives of the 3D spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. + * @param[in] coords_eval The coordinates where the spline is differentiated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 3D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template < + class InterestDim, + class Layout1, + class Layout2, + class Layout3, + class BatchedInterpolationDDom, + class... CoordsDims> + void deriv( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + ddc::Coordinate const, + BatchedInterpolationDDom, + Layout2, + memory_space> const coords_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const + { + static_assert( + std::is_same_v< + InterestDim, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + || std::is_same_v + || std::is_same_v); + if constexpr (std::is_same_v< + InterestDim, + typename evaluation_discrete_dimension_type1:: + continuous_dimension_type>) { + return deriv_dim_1(spline_eval, coords_eval, spline_coef); + } else if constexpr (std::is_same_v< + InterestDim, + typename evaluation_discrete_dimension_type2:: + continuous_dimension_type>) { + return deriv_dim_2(spline_eval, coords_eval, spline_coef); + } else if constexpr (std::is_same_v< + InterestDim, + typename evaluation_discrete_dimension_type3:: + continuous_dimension_type>) { + return deriv_dim_3(spline_eval, coords_eval, spline_coef); + } + } + + /** + * @brief Differentiate spline function (described by its spline coefficients) on a mesh along a specified dimension of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a multidimensional evaluation. This is a batched 3D evaluation. + * This means that for each slice of spline_eval the evaluation is performed with + * the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @tparam InterestDim Dimension along which differentiation is performed. + * @param[out] spline_eval The derivatives of the 3D spline function at the desired coordinates. + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template + void deriv( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const + { + static_assert( + std::is_same_v< + InterestDim, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + || std::is_same_v + || std::is_same_v); + if constexpr (std::is_same_v< + InterestDim, + typename evaluation_discrete_dimension_type1:: + continuous_dimension_type>) { + return deriv_dim_1(spline_eval, spline_coef); + } else if constexpr (std::is_same_v< + InterestDim, + typename evaluation_discrete_dimension_type2:: + continuous_dimension_type>) { + return deriv_dim_2(spline_eval, spline_coef); + } else if constexpr (std::is_same_v< + InterestDim, + typename evaluation_discrete_dimension_type3:: + continuous_dimension_type>) { + return deriv_dim_3(spline_eval, spline_coef); + } + } + + /** + * @brief Double-differentiate 3D spline function (described by its spline coefficients) on a mesh along specified dimensions of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a nD evaluation. This is a batched 3D differentiation. + * This means that for each slice of coordinates identified by a batch_domain_type::discrete_element_type, + * the differentiation is performed with the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * Note: double-differentiation other than cross-differentiation is not supported atm. See #440 + * + * @tparam InterestDim1 First dimension along which differentiation is performed. + * @tparam InterestDim2 Second dimension along which differentiation is performed. + * + * @param[out] spline_eval The derivatives of the 3D spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. + * @param[in] coords_eval The coordinates where the spline is differentiated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 3D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template < + class InterestDim1, + class InterestDim2, + class Layout1, + class Layout2, + class Layout3, + class BatchedInterpolationDDom, + class... CoordsDims> + void deriv2( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + ddc::Coordinate const, + BatchedInterpolationDDom, + Layout2, + memory_space> const coords_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const + { + static_assert( + (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type2::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type2::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v)); + + if constexpr ( + (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v)) { + return deriv_1_and_2(spline_eval, coords_eval, spline_coef); + } else if constexpr ( + (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type2::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type2::continuous_dimension_type> + && std::is_same_v)) { + return deriv_2_and_3(spline_eval, coords_eval, spline_coef); + } else if constexpr ( + (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v)) { + return deriv_1_and_3(spline_eval, coords_eval, spline_coef); + } + } + + /** + * @brief Double-differentiate 3D spline function (described by its spline coefficients) on a mesh along specified dimensions of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a multidimensional evaluation. This is a batched 3D evaluation. + * This means that for each slice of spline_eval the evaluation is performed with + * the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * Note: double-differentiation other than cross-differentiation is not supported atm. See #440 + * + * @tparam InterestDim1 First dimension along which differentiation is performed. + * @tparam InterestDim2 Second dimension along which differentiation is performed. + * + * @param[out] spline_eval The derivatives of the 3D spline function at the desired coordinates. + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template < + class InterestDim1, + class InterestDim2, + class Layout1, + class Layout2, + class BatchedInterpolationDDom> + void deriv2( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const + { + static_assert( + (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type2::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type2::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v)); + + if constexpr ( + (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v)) { + return deriv_1_and_2(spline_eval, spline_coef); + } else if constexpr ( + (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type2::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type2::continuous_dimension_type> + && std::is_same_v)) { + return deriv_2_and_3(spline_eval, spline_coef); + } else if constexpr ( + (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v)) { + return deriv_1_and_3(spline_eval, spline_coef); + } + } + + /** + * @brief Differentiate spline function (described by its spline coefficients) on a mesh along a specified dimensions of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a nD evaluation. This is a batched 3D differentiation. + * This means that for each slice of coordinates identified by a batch_domain_type::discrete_element_type, + * the differentiation is performed with the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * Note: triple-differentiation other than cross-differentiation is not supported atm. See #440 + * + * @tparam InterestDim1 First dimension along which differentiation is performed. + * @tparam InterestDim2 Second dimension along which differentiation is performed. + * @tparam InterestDim3 Third dimension along which differentiation is performed. + * + * @param[out] spline_eval The derivatives of the 3D spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. + * @param[in] coords_eval The coordinates where the spline is differentiated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 3D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template < + class InterestDim1, + class InterestDim2, + class InterestDim3, + class Layout1, + class Layout2, + class Layout3, + class BatchedInterpolationDDom, + class... CoordsDims> + void deriv3( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + ddc::Coordinate const, + BatchedInterpolationDDom, + Layout2, + memory_space> const coords_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const + { + static_assert( + (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v + && std::is_same_v) + || (std::is_same_v< + InterestDim3, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v + && std::is_same_v)); + + return deriv_1_2_3(spline_eval, coords_eval, spline_coef); + } + + /** + * @brief Differentiate spline function (described by its spline coefficients) on a mesh along specified dimensions of interest. + * + * The spline coefficients represent a 3D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a multidimensional evaluation. This is a batched 3D evaluation. + * This means that for each slice of spline_eval the evaluation is performed with + * the 3D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * Note: triple-differentiation other than cross-differentiation is not supported atm. See #440 + * + * @tparam InterestDim1 First dimension along which differentiation is performed. + * @tparam InterestDim2 Second dimension along which differentiation is performed. + * @tparam InterestDim3 Third dimension along which differentiation is performed. + * + * @param[out] spline_eval The derivatives of the 3D spline function at the desired coordinates. + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template < + class InterestDim1, + class InterestDim2, + class InterestDim3, + class Layout1, + class Layout2, + class BatchedInterpolationDDom> + void deriv3( + ddc::ChunkSpan const + spline_eval, + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const + { + static_assert( + (std::is_same_v< + InterestDim1, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v + && std::is_same_v) + || (std::is_same_v< + InterestDim3, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v + && std::is_same_v) + || (std::is_same_v< + InterestDim2, + typename evaluation_discrete_dimension_type1::continuous_dimension_type> + && std::is_same_v + && std::is_same_v)); + + return deriv_1_2_3(spline_eval, spline_coef); + } + + /** @brief Perform batched 3D integrations of a spline function (described by its spline coefficients) along the dimensions of interest and store results on a subdomain of batch_domain. + * + * The spline coefficients represent a 3D spline function defined on a B-splines (basis splines). They can be obtained via various methods, such as using a SplineBuilder3D. + * + * This is not a nD integration. This is a batched 3D integration. + * This means that for each element of integrals, the integration is performed with the 3D set of + * spline coefficients identified by the same DiscreteElement. + * + * @param[out] integrals The integrals of the 3D spline function on the subdomain of batch_domain. For practical reasons those are + * stored in a ChunkSpan defined on a batch_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant. + * @param[in] spline_coef A ChunkSpan storing the 3D spline coefficients. + */ + template + void integrate( + ddc::ChunkSpan const integrals, + ddc::ChunkSpan const + spline_coef) const + { + static_assert( + ddc::type_seq_contains_v< + ddc::detail::TypeSeq, + to_type_seq_t>, + "The spline coefficients domain must contain the bsplines dimensions"); + using batch_domain_type = ddc:: + remove_dims_of_t; + static_assert( + std::is_same_v, + "The integrals domain must only contain the batch dimensions"); + + batch_domain_type batch_domain(integrals.domain()); + ddc::Chunk values1_alloc( + ddc::DiscreteDomain(spline_coef.domain()), + ddc::KokkosAllocator()); + ddc::ChunkSpan values1 = values1_alloc.span_view(); + ddc::integrals(exec_space(), values1); + ddc::Chunk values2_alloc( + ddc::DiscreteDomain(spline_coef.domain()), + ddc::KokkosAllocator()); + ddc::ChunkSpan values2 = values2_alloc.span_view(); + ddc::integrals(exec_space(), values2); + ddc::Chunk values3_alloc( + ddc::DiscreteDomain(spline_coef.domain()), + ddc::KokkosAllocator()); + ddc::ChunkSpan values3 = values3_alloc.span_view(); + ddc::integrals(exec_space(), values3); + + ddc::parallel_for_each( + "ddc_splines_integrate_bsplines", + exec_space(), + batch_domain, + KOKKOS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { + integrals(j) = 0; + for (typename spline_domain_type1::discrete_element_type const i1 : + values1.domain()) { + for (typename spline_domain_type2::discrete_element_type const i2 : + values2.domain()) { + for (typename spline_domain_type3::discrete_element_type const i3 : + values3.domain()) { + integrals(j) += spline_coef(i1, i2, i3, j) * values1(i1) + * values2(i2) * values3(i3); + } + } + } + }); + } + +private: + /** + * @brief Evaluate the function on B-splines at the coordinate given. + * + * This function firstly deals with the boundary conditions and calls the SplineEvaluator3D::eval_no_bc function + * to evaluate. + * + * @param[in] coord_eval The 3D coordinate where we want to evaluate. + * @param[in] spline_coef The B-splines coefficients of the function we want to evaluate. + * @param[out] vals1 A ChunkSpan with the not-null values of each function of the spline in the first dimension. + * @param[out] vals2 A ChunkSpan with the not-null values of each function of the spline in the second dimension. + * + * @return A double with the value of the function at the coordinate given. + * + * @see SplineBoundaryValue + */ + template + KOKKOS_INLINE_FUNCTION double eval( + ddc::Coordinate coord_eval, + ddc::ChunkSpan const + spline_coef) const + { + using Dim1 = continuous_dimension_type1; + using Dim2 = continuous_dimension_type2; + using Dim3 = continuous_dimension_type3; + if constexpr (bsplines_type1::is_periodic()) { + if (ddc::get(coord_eval) < ddc::discrete_space().rmin() + || ddc::get(coord_eval) > ddc::discrete_space().rmax()) { + ddc::get(coord_eval) + -= Kokkos::floor( + (ddc::get(coord_eval) + - ddc::discrete_space().rmin()) + / ddc::discrete_space().length()) + * ddc::discrete_space().length(); + } + } + if constexpr (bsplines_type2::is_periodic()) { + if (ddc::get(coord_eval) < ddc::discrete_space().rmin() + || ddc::get(coord_eval) > ddc::discrete_space().rmax()) { + ddc::get(coord_eval) + -= Kokkos::floor( + (ddc::get(coord_eval) + - ddc::discrete_space().rmin()) + / ddc::discrete_space().length()) + * ddc::discrete_space().length(); + } + } + if constexpr (bsplines_type3::is_periodic()) { + if (ddc::get(coord_eval) < ddc::discrete_space().rmin() + || ddc::get(coord_eval) > ddc::discrete_space().rmax()) { + ddc::get(coord_eval) + -= Kokkos::floor( + (ddc::get(coord_eval) + - ddc::discrete_space().rmin()) + / ddc::discrete_space().length()) + * ddc::discrete_space().length(); + } + } + if constexpr (!bsplines_type1::is_periodic()) { + if (ddc::get(coord_eval) < ddc::discrete_space().rmin()) { + return m_lower_extrap_rule_1(coord_eval, spline_coef); + } + if (ddc::get(coord_eval) > ddc::discrete_space().rmax()) { + return m_upper_extrap_rule_1(coord_eval, spline_coef); + } + } + if constexpr (!bsplines_type2::is_periodic()) { + if (ddc::get(coord_eval) < ddc::discrete_space().rmin()) { + return m_lower_extrap_rule_2(coord_eval, spline_coef); + } + if (ddc::get(coord_eval) > ddc::discrete_space().rmax()) { + return m_upper_extrap_rule_2(coord_eval, spline_coef); + } + } + if constexpr (!bsplines_type3::is_periodic()) { + if (ddc::get(coord_eval) < ddc::discrete_space().rmin()) { + return m_lower_extrap_rule_3(coord_eval, spline_coef); + } + if (ddc::get(coord_eval) > ddc::discrete_space().rmax()) { + return m_upper_extrap_rule_3(coord_eval, spline_coef); + } + } + return eval_no_bc( + ddc::Coordinate< + continuous_dimension_type1, + continuous_dimension_type2, + continuous_dimension_type3>( + ddc::get(coord_eval), + ddc::get(coord_eval), + ddc::get(coord_eval)), + spline_coef); + } + + /** + * @brief Evaluate the function or its derivative at the coordinate given. + * + * @param[in] coord_eval The coordinate where we want to evaluate. + * @param[in] splne_coef The B-splines coefficients of the function we want to evaluate. + * @tparam EvalType1 A flag indicating if we evaluate the function or its derivative in the first dimension. The type of this object is either `eval_type` or `eval_deriv_type`. + * @tparam EvalType2 A flag indicating if we evaluate the function or its derivative in the second dimension. The type of this object is either `eval_type` or `eval_deriv_type`. + */ + template + KOKKOS_INLINE_FUNCTION double eval_no_bc( + ddc::Coordinate const& coord_eval, + ddc::ChunkSpan const + spline_coef) const + { + static_assert( + std::is_same_v || std::is_same_v); + static_assert( + std::is_same_v || std::is_same_v); + static_assert( + std::is_same_v || std::is_same_v); + ddc::DiscreteElement jmin1; + ddc::DiscreteElement jmin2; + ddc::DiscreteElement jmin3; + + std::array vals1_ptr; + Kokkos::mdspan> const + vals1(vals1_ptr.data()); + std::array vals2_ptr; + Kokkos::mdspan> const + vals2(vals2_ptr.data()); + std::array vals3_ptr; + Kokkos::mdspan> const + vals3(vals3_ptr.data()); + ddc::Coordinate const coord_eval_interest1(coord_eval); + ddc::Coordinate const coord_eval_interest2(coord_eval); + ddc::Coordinate const coord_eval_interest3(coord_eval); + + if constexpr (std::is_same_v) { + jmin1 = ddc::discrete_space().eval_basis(vals1, coord_eval_interest1); + } else if constexpr (std::is_same_v) { + jmin1 = ddc::discrete_space().eval_deriv(vals1, coord_eval_interest1); + } + if constexpr (std::is_same_v) { + jmin2 = ddc::discrete_space().eval_basis(vals2, coord_eval_interest2); + } else if constexpr (std::is_same_v) { + jmin2 = ddc::discrete_space().eval_deriv(vals2, coord_eval_interest2); + } + if constexpr (std::is_same_v) { + jmin3 = ddc::discrete_space().eval_basis(vals3, coord_eval_interest3); + } else if constexpr (std::is_same_v) { + jmin3 = ddc::discrete_space().eval_deriv(vals3, coord_eval_interest3); + } + + double y = 0.0; + for (std::size_t i = 0; i < bsplines_type1::degree() + 1; ++i) { + for (std::size_t j = 0; j < bsplines_type2::degree() + 1; ++j) { + for (std::size_t k = 0; k < bsplines_type3::degree() + 1; ++k) { + y += spline_coef( + ddc::DiscreteElement< + bsplines_type1, + bsplines_type2, + bsplines_type3>(jmin1 + i, jmin2 + j, jmin3 + k)) + * vals1[i] * vals2[j] * vals3[k]; + } + } + } + return y; + } +}; + +} // namespace ddc diff --git a/tests/splines/CMakeLists.txt b/tests/splines/CMakeLists.txt index 2a82e763b..a328c1251 100644 --- a/tests/splines/CMakeLists.txt +++ b/tests/splines/CMakeLists.txt @@ -144,6 +144,26 @@ foreach(BC "BC_PERIODIC" "BC_GREVILLE" "BC_HERMITE") endforeach() endforeach() +foreach(BC "BC_PERIODIC") #"BC_GREVILLE" "BC_HERMITE") + foreach(EVALUATOR "EVALUATOR_POLYNOMIAL") + foreach(DEGREE RANGE "${DDC_SPLINES_TEST_DEGREE_MIN}" "${DDC_SPLINES_TEST_DEGREE_MAX}") + foreach(BSPLINES_TYPE "BSPLINES_TYPE_UNIFORM" "BSPLINES_TYPE_NON_UNIFORM") + set(test_name + "3d_splines_tests_BATCHED_DEGREE_${DEGREE}_${BSPLINES_TYPE}_${EVALUATOR}_${BC}" + ) + add_executable("${test_name}" ../main.cpp batched_3d_spline_builder.cpp) + target_compile_definitions( + "${test_name}" + PUBLIC -DDEGREE=${DEGREE} -D${BSPLINES_TYPE} -D${EVALUATOR} -D${BC} + ) + target_compile_features("${test_name}" PUBLIC cxx_std_17) + target_link_libraries("${test_name}" PUBLIC DDC::core DDC::splines GTest::gtest) + gtest_discover_tests("${test_name}" DISCOVERY_MODE PRE_TEST) + endforeach() + endforeach() + endforeach() +endforeach() + foreach(SOLVER "GINKGO" "LAPACK") foreach(BC "BC_PERIODIC" "BC_GREVILLE" "BC_HERMITE") foreach(DEGREE_X RANGE "${DDC_SPLINES_TEST_DEGREE_MIN}" "${DDC_SPLINES_TEST_DEGREE_MAX}") diff --git a/tests/splines/batched_3d_spline_builder.cpp b/tests/splines/batched_3d_spline_builder.cpp new file mode 100644 index 000000000..9f3c95b93 --- /dev/null +++ b/tests/splines/batched_3d_spline_builder.cpp @@ -0,0 +1,852 @@ +// Copyright (C) The DDC development team, see COPYRIGHT.md file +// +// SPDX-License-Identifier: MIT + +#include +#include +#if defined(BC_HERMITE) +# include +#endif +#if defined(BSPLINES_TYPE_UNIFORM) +# include +#endif +#include + +#include +#include + +#include + +#include + +#include "evaluator_3d.hpp" +#if defined(BC_PERIODIC) +# include "cosine_evaluator.hpp" +#else +# include "polynomial_evaluator.hpp" +#endif +#include "spline_error_bounds.hpp" + +inline namespace anonymous_namespace_workaround_batched_3d_spline_builder_cpp { + +#if defined(BC_PERIODIC) +struct DimX +{ + static constexpr bool PERIODIC = true; +}; + +struct DimY +{ + static constexpr bool PERIODIC = true; +}; + +struct DimZ +{ + static constexpr bool PERIODIC = true; +}; +#else + +struct DimX +{ + static constexpr bool PERIODIC = false; +}; + +struct DimY +{ + static constexpr bool PERIODIC = false; +}; + +struct DimZ +{ + static constexpr bool PERIODIC = false; +}; +#endif + +struct DDimBatch +{ +}; + +constexpr std::size_t s_degree = DEGREE; + +#if defined(BC_PERIODIC) +constexpr ddc::BoundCond s_bcl = ddc::BoundCond::PERIODIC; +constexpr ddc::BoundCond s_bcr = ddc::BoundCond::PERIODIC; +#elif defined(BC_GREVILLE) +constexpr ddc::BoundCond s_bcl = ddc::BoundCond::GREVILLE; +constexpr ddc::BoundCond s_bcr = ddc::BoundCond::GREVILLE; +#elif defined(BC_HERMITE) +constexpr ddc::BoundCond s_bcl = ddc::BoundCond::HERMITE; +constexpr ddc::BoundCond s_bcr = ddc::BoundCond::HERMITE; +#endif + +template +using GrevillePoints = ddc::GrevilleInterpolationPoints; + +#if defined(BSPLINES_TYPE_UNIFORM) +template +struct BSplines : ddc::UniformBSplines +{ +}; +#elif defined(BSPLINES_TYPE_NON_UNIFORM) +template +struct BSplines : ddc::NonUniformBSplines +{ +}; +#endif + +// In the dimensions of interest, the discrete dimension is deduced from the Greville points type. +template +struct DDimGPS : GrevillePoints>::interpolation_discrete_dimension_type +{ +}; + +#if defined(BC_PERIODIC) +template +using evaluator_type = Evaluator3D::Evaluator< + CosineEvaluator::Evaluator, + CosineEvaluator::Evaluator, + CosineEvaluator::Evaluator>; +#else +template +using evaluator_type = Evaluator3D::Evaluator< + PolynomialEvaluator::Evaluator, + PolynomialEvaluator::Evaluator, + PolynomialEvaluator::Evaluator>; +#endif + +template +using DElem = ddc::DiscreteElement; +template +using DVect = ddc::DiscreteVector; +template +using Coord = ddc::Coordinate; + +// Templated function giving first coordinate of the mesh in given dimension. +template +KOKKOS_FUNCTION Coord x0() +{ + return Coord(0.); +} + +// Templated function giving last coordinate of the mesh in given dimension. +template +KOKKOS_FUNCTION Coord xN() +{ + return Coord(1.); +} + +// Templated function giving step of the mesh in given dimension. +template +double dx(std::size_t ncells) +{ + return (xN() - x0()) / ncells; +} + +// Templated function giving break points of mesh in given dimension for non-uniform case. +template +std::vector> breaks(std::size_t ncells) +{ + std::vector> out(ncells + 1); + for (std::size_t i(0); i < ncells + 1; ++i) { + out[i] = x0() + i * dx(ncells); + } + return out; +} + +template +void InterestDimInitializer(std::size_t const ncells) +{ + using CDim = typename DDim::continuous_dimension_type; +#if defined(BSPLINES_TYPE_UNIFORM) + ddc::init_discrete_space>(x0(), xN(), ncells); +#elif defined(BSPLINES_TYPE_NON_UNIFORM) + ddc::init_discrete_space>(breaks(ncells)); +#endif + ddc::init_discrete_space(GrevillePoints>::template get_sampling()); +} + +// Checks that when evaluating the spline at interpolation points one +// recovers values that were used to build the spline +template < + typename ExecSpace, + typename MemorySpace, + typename DDimI1, + typename DDimI2, + typename DDimI3, + typename... DDims> +void Batched3dSplineTest() +{ + using I1 = typename DDimI1::continuous_dimension_type; + using I2 = typename DDimI2::continuous_dimension_type; + using I3 = typename DDimI3::continuous_dimension_type; + + // Instantiate execution spaces and initialize spaces + ExecSpace const exec_space; + std::size_t const ncells = 10; + InterestDimInitializer(ncells); + InterestDimInitializer(ncells); + InterestDimInitializer(ncells); + + // Create the values domain (mesh) + ddc::DiscreteDomain const interpolation_domain1 + = GrevillePoints>::template get_domain(); + ddc::DiscreteDomain const interpolation_domain2 + = GrevillePoints>::template get_domain(); + ddc::DiscreteDomain const interpolation_domain3 + = GrevillePoints>::template get_domain(); + ddc::DiscreteDomain const interpolation_domain( + interpolation_domain1, + interpolation_domain2, + interpolation_domain3); + // The following line creates a discrete domain over all dimensions (DDims...) except DDimI1, DDimI2 and DDimI3. + auto const dom_vals_tmp + = ddc::remove_dims_of_t, DDimI1, DDimI2, DDimI3>( + ddc::DiscreteDomain(DElem(0), DVect(ncells))...); + ddc::DiscreteDomain const dom_vals( + dom_vals_tmp, + interpolation_domain1, + interpolation_domain2, + interpolation_domain3); + +#if defined(BC_HERMITE) + // Create the derivs domain + ddc::DiscreteDomain> const + derivs_domain1(DElem>(1), DVect>(s_degree / 2)); + ddc::DiscreteDomain> const + derivs_domain2(DElem>(1), DVect>(s_degree / 2)); + ddc::DiscreteDomain> const + derivs_domain3(DElem>(1), DVect>(s_degree / 2)); + ddc::DiscreteDomain, ddc::Deriv, ddc::Deriv> const + derivs_domain(derivs_domain1, derivs_domain2, derivs_domain3); + + auto const dom_derivs_1d + = ddc::replace_dim_of>(dom_vals, derivs_domain1); + auto const dom_derivs2 = ddc::replace_dim_of>(dom_vals, derivs_domain2); + auto const dom_derivs3 = ddc::replace_dim_of>(dom_vals, derivs_domain3); + auto const dom_derivs = ddc::replace_dim_of>( + ddc::replace_dim_of>(dom_derivs_1d, derivs_domain2), + derivs_domain3); + // FIXME: check later the correctness of this +#endif + + // Create a SplineBuilder over BSplines and batched along other dimensions using some boundary conditions + ddc::SplineBuilder3D< + ExecSpace, + MemorySpace, + BSplines, + BSplines, + BSplines, + DDimI1, + DDimI2, + DDimI3, + s_bcl, + s_bcr, + s_bcl, + s_bcr, + s_bcl, + s_bcr, + ddc::SplineSolver::GINKGO> const spline_builder(interpolation_domain); + + // Compute useful domains (dom_interpolation, dom_batch, dom_bsplines and dom_spline) + ddc::DiscreteDomain const dom_interpolation + = spline_builder.interpolation_domain(); + auto const dom_spline = spline_builder.batched_spline_domain(dom_vals); + + // Allocate and fill a chunk containing values to be passed as input to spline_builder. Those are values of cosine along interest dimension duplicated along batch dimensions + ddc::Chunk vals_1d_host_alloc(dom_interpolation, ddc::HostAllocator()); + ddc::ChunkSpan const vals_1d_host = vals_1d_host_alloc.span_view(); + evaluator_type const evaluator(dom_interpolation); + evaluator(vals_1d_host); + auto vals_1d_alloc = ddc::create_mirror_view_and_copy(exec_space, vals_1d_host); + ddc::ChunkSpan const vals_1d = vals_1d_alloc.span_view(); + + ddc::Chunk vals_alloc(dom_vals, ddc::KokkosAllocator()); + ddc::ChunkSpan const vals = vals_alloc.span_view(); + ddc::parallel_for_each( + exec_space, + vals.domain(), + KOKKOS_LAMBDA(DElem const e) { + vals(e) = vals_1d(DElem(e)); + }); + +#if defined(BC_HERMITE) + // FIXME: handle this later + + // Allocate and fill a chunk containing derivs to be passed as input to spline_builder. + int const shift = s_degree % 2; // shift = 0 for even order, 1 for odd order + ddc::Chunk derivs_1d_lhs_alloc(dom_derivs_1d, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs_1d_lhs = derivs_1d_lhs_alloc.span_view(); + if (s_bcl == ddc::BoundCond::HERMITE) { + ddc::Chunk derivs_1d_lhs1_host_alloc( + ddc::DiscreteDomain, DDimI2>(derivs_domain1, interpolation_domain2), + ddc::HostAllocator()); + ddc::ChunkSpan const derivs_1d_lhs1_host = derivs_1d_lhs1_host_alloc.span_view(); + ddc::for_each( + derivs_1d_lhs1_host.domain(), + KOKKOS_LAMBDA(ddc::DiscreteElement, DDimI2> const e) { + auto deriv_idx = ddc::DiscreteElement>(e).uid(); + auto x2 = ddc::coordinate(ddc::DiscreteElement(e)); + derivs_1d_lhs1_host(e) + = evaluator.deriv(x0(), x2, deriv_idx + shift - 1, 0); + }); + auto derivs_1d_lhs1_alloc + = ddc::create_mirror_view_and_copy(exec_space, derivs_1d_lhs1_host); + ddc::ChunkSpan const derivs_1d_lhs1 = derivs_1d_lhs1_alloc.span_view(); + + ddc::parallel_for_each( + exec_space, + derivs_1d_lhs.domain(), + KOKKOS_LAMBDA( + typename decltype(derivs_1d_lhs.domain())::discrete_element_type const e) { + derivs_1d_lhs(e) = derivs_1d_lhs1(DElem, DDimI2>(e)); + }); + } + + ddc::Chunk derivs_1d_rhs_alloc(dom_derivs_1d, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs_1d_rhs = derivs_1d_rhs_alloc.span_view(); + if (s_bcl == ddc::BoundCond::HERMITE) { + ddc::Chunk derivs_1d_rhs1_host_alloc( + ddc::DiscreteDomain, DDimI2>(derivs_domain1, interpolation_domain2), + ddc::HostAllocator()); + ddc::ChunkSpan const derivs_1d_rhs1_host = derivs_1d_rhs1_host_alloc.span_view(); + ddc::for_each( + derivs_1d_rhs1_host.domain(), + KOKKOS_LAMBDA(ddc::DiscreteElement, DDimI2> const e) { + auto deriv_idx = ddc::DiscreteElement>(e).uid(); + auto x2 = ddc::coordinate(ddc::DiscreteElement(e)); + derivs_1d_rhs1_host(e) + = evaluator.deriv(xN(), x2, deriv_idx + shift - 1, 0); + }); + auto derivs_1d_rhs1_alloc + = ddc::create_mirror_view_and_copy(exec_space, derivs_1d_rhs1_host); + ddc::ChunkSpan const derivs_1d_rhs1 = derivs_1d_rhs1_alloc.span_view(); + + ddc::parallel_for_each( + exec_space, + derivs_1d_rhs.domain(), + KOKKOS_LAMBDA( + typename decltype(derivs_1d_rhs.domain())::discrete_element_type const e) { + derivs_1d_rhs(e) = derivs_1d_rhs1(DElem, DDimI2>(e)); + }); + } + + ddc::Chunk derivs2_lhs_alloc(dom_derivs2, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs2_lhs = derivs2_lhs_alloc.span_view(); + if (s_bcl == ddc::BoundCond::HERMITE) { + ddc::Chunk derivs2_lhs1_host_alloc( + ddc::DiscreteDomain>(interpolation_domain1, derivs_domain2), + ddc::HostAllocator()); + ddc::ChunkSpan const derivs2_lhs1_host = derivs2_lhs1_host_alloc.span_view(); + ddc::for_each( + derivs2_lhs1_host.domain(), + KOKKOS_LAMBDA(ddc::DiscreteElement> const e) { + auto x1 = ddc::coordinate(ddc::DiscreteElement(e)); + auto deriv_idx = ddc::DiscreteElement>(e).uid(); + derivs2_lhs1_host(e) = evaluator.deriv(x1, x0(), 0, deriv_idx + shift - 1); + }); + + auto derivs2_lhs1_alloc = ddc::create_mirror_view_and_copy(exec_space, derivs2_lhs1_host); + ddc::ChunkSpan const derivs2_lhs1 = derivs2_lhs1_alloc.span_view(); + + ddc::parallel_for_each( + exec_space, + derivs2_lhs.domain(), + KOKKOS_LAMBDA( + typename decltype(derivs2_lhs.domain())::discrete_element_type const e) { + derivs2_lhs(e) = derivs2_lhs1(DElem>(e)); + }); + } + + ddc::Chunk derivs2_rhs_alloc(dom_derivs2, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs2_rhs = derivs2_rhs_alloc.span_view(); + if (s_bcl == ddc::BoundCond::HERMITE) { + ddc::Chunk derivs2_rhs1_host_alloc( + ddc::DiscreteDomain>(interpolation_domain1, derivs_domain2), + ddc::HostAllocator()); + ddc::ChunkSpan const derivs2_rhs1_host = derivs2_rhs1_host_alloc.span_view(); + ddc::for_each( + derivs2_rhs1_host.domain(), + KOKKOS_LAMBDA(ddc::DiscreteElement> const e) { + auto x1 = ddc::coordinate(ddc::DiscreteElement(e)); + auto deriv_idx = ddc::DiscreteElement>(e).uid(); + derivs2_rhs1_host(e) = evaluator.deriv(x1, xN(), 0, deriv_idx + shift - 1); + }); + + auto derivs2_rhs1_alloc = ddc::create_mirror_view_and_copy(exec_space, derivs2_rhs1_host); + ddc::ChunkSpan const derivs2_rhs1 = derivs2_rhs1_alloc.span_view(); + + ddc::parallel_for_each( + exec_space, + derivs2_rhs.domain(), + KOKKOS_LAMBDA( + typename decltype(derivs2_rhs.domain())::discrete_element_type const e) { + derivs2_rhs(e) = derivs2_rhs1(DElem>(e)); + }); + } + + ddc::Chunk derivs_mixed_lhs_lhs_alloc(dom_derivs, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs_mixed_lhs_lhs = derivs_mixed_lhs_lhs_alloc.span_view(); + ddc::Chunk derivs_mixed_rhs_lhs_alloc(dom_derivs, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs_mixed_rhs_lhs = derivs_mixed_rhs_lhs_alloc.span_view(); + ddc::Chunk derivs_mixed_lhs_rhs_alloc(dom_derivs, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs_mixed_lhs_rhs = derivs_mixed_lhs_rhs_alloc.span_view(); + ddc::Chunk derivs_mixed_rhs_rhs_alloc(dom_derivs, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs_mixed_rhs_rhs = derivs_mixed_rhs_rhs_alloc.span_view(); + + if (s_bcl == ddc::BoundCond::HERMITE) { + ddc::Chunk derivs_mixed_lhs_lhs1_host_alloc(derivs_domain, ddc::HostAllocator()); + ddc::ChunkSpan const derivs_mixed_lhs_lhs1_host + = derivs_mixed_lhs_lhs1_host_alloc.span_view(); + ddc::Chunk derivs_mixed_rhs_lhs1_host_alloc(derivs_domain, ddc::HostAllocator()); + ddc::ChunkSpan const derivs_mixed_rhs_lhs1_host + = derivs_mixed_rhs_lhs1_host_alloc.span_view(); + ddc::Chunk derivs_mixed_lhs_rhs1_host_alloc(derivs_domain, ddc::HostAllocator()); + ddc::ChunkSpan const derivs_mixed_lhs_rhs1_host + = derivs_mixed_lhs_rhs1_host_alloc.span_view(); + ddc::Chunk derivs_mixed_rhs_rhs1_host_alloc(derivs_domain, ddc::HostAllocator()); + ddc::ChunkSpan const derivs_mixed_rhs_rhs1_host + = derivs_mixed_rhs_rhs1_host_alloc.span_view(); + + for (std::size_t ii = 1; + ii < static_cast(derivs_domain.template extent>()) + 1; + ++ii) { + for (std::size_t jj = 1; + jj < static_cast(derivs_domain.template extent>()) + 1; + ++jj) { + derivs_mixed_lhs_lhs1_host( + typename decltype(derivs_domain)::discrete_element_type(ii, jj)) + = evaluator.deriv(x0(), x0(), ii + shift - 1, jj + shift - 1); + derivs_mixed_rhs_lhs1_host( + typename decltype(derivs_domain)::discrete_element_type(ii, jj)) + = evaluator.deriv(xN(), x0(), ii + shift - 1, jj + shift - 1); + derivs_mixed_lhs_rhs1_host( + typename decltype(derivs_domain)::discrete_element_type(ii, jj)) + = evaluator.deriv(x0(), xN(), ii + shift - 1, jj + shift - 1); + derivs_mixed_rhs_rhs1_host( + typename decltype(derivs_domain)::discrete_element_type(ii, jj)) + = evaluator.deriv(xN(), xN(), ii + shift - 1, jj + shift - 1); + } + } + auto derivs_mixed_lhs_lhs1_alloc + = ddc::create_mirror_view_and_copy(exec_space, derivs_mixed_lhs_lhs1_host); + ddc::ChunkSpan const derivs_mixed_lhs_lhs1 = derivs_mixed_lhs_lhs1_alloc.span_view(); + auto derivs_mixed_rhs_lhs1_alloc + = ddc::create_mirror_view_and_copy(exec_space, derivs_mixed_rhs_lhs1_host); + ddc::ChunkSpan const derivs_mixed_rhs_lhs1 = derivs_mixed_rhs_lhs1_alloc.span_view(); + auto derivs_mixed_lhs_rhs1_alloc + = ddc::create_mirror_view_and_copy(exec_space, derivs_mixed_lhs_rhs1_host); + ddc::ChunkSpan const derivs_mixed_lhs_rhs1 = derivs_mixed_lhs_rhs1_alloc.span_view(); + auto derivs_mixed_rhs_rhs1_alloc + = ddc::create_mirror_view_and_copy(exec_space, derivs_mixed_rhs_rhs1_host); + ddc::ChunkSpan const derivs_mixed_rhs_rhs1 = derivs_mixed_rhs_rhs1_alloc.span_view(); + + ddc::parallel_for_each( + exec_space, + dom_derivs, + KOKKOS_LAMBDA(typename decltype(dom_derivs)::discrete_element_type const e) { + derivs_mixed_lhs_lhs(e) + = derivs_mixed_lhs_lhs1(DElem, ddc::Deriv>(e)); + derivs_mixed_rhs_lhs(e) + = derivs_mixed_rhs_lhs1(DElem, ddc::Deriv>(e)); + derivs_mixed_lhs_rhs(e) + = derivs_mixed_lhs_rhs1(DElem, ddc::Deriv>(e)); + derivs_mixed_rhs_rhs(e) + = derivs_mixed_rhs_rhs1(DElem, ddc::Deriv>(e)); + }); + } +#endif + + // Instantiate chunk of spline coefs to receive output of spline_builder + ddc::Chunk coef_alloc(dom_spline, ddc::KokkosAllocator()); + ddc::ChunkSpan const coef = coef_alloc.span_view(); + + // Finally compute the spline by filling `coef` +#if defined(BC_HERMITE) + // FIXME: handle this later + spline_builder( + coef, + vals.span_cview(), + std::optional(derivs_1d_lhs.span_cview()), + std::optional(derivs_1d_rhs.span_cview()), + std::optional(derivs2_lhs.span_cview()), + std::optional(derivs2_rhs.span_cview()), + std::optional(derivs_mixed_lhs_lhs.span_cview()), + std::optional(derivs_mixed_rhs_lhs.span_cview()), + std::optional(derivs_mixed_lhs_rhs.span_cview()), + std::optional(derivs_mixed_rhs_rhs.span_cview())); +#else + spline_builder(coef, vals.span_cview()); +#endif + + // Instantiate a SplineEvaluator over interest dimension and batched along other dimensions +#if defined(BC_PERIODIC) + using extrapolation_rule_1_type = ddc::PeriodicExtrapolationRule; + using extrapolation_rule_2_type = ddc::PeriodicExtrapolationRule; + using extrapolation_rule_3_type = ddc::PeriodicExtrapolationRule; +#else + using extrapolation_rule_1_type = ddc::NullExtrapolationRule; + using extrapolation_rule_2_type = ddc::NullExtrapolationRule; + using extrapolation_rule_3_type = ddc::NullExtrapolationRule; +#endif + extrapolation_rule_1_type const extrapolation_rule_1; + extrapolation_rule_2_type const extrapolation_rule_2; + extrapolation_rule_3_type const extrapolation_rule_3; + + ddc::SplineEvaluator3D< + ExecSpace, + MemorySpace, + BSplines, + BSplines, + BSplines, + DDimI1, + DDimI2, + DDimI3, + extrapolation_rule_1_type, + extrapolation_rule_1_type, + extrapolation_rule_2_type, + extrapolation_rule_2_type, + extrapolation_rule_3_type, + extrapolation_rule_3_type> const + spline_evaluator( + extrapolation_rule_1, + extrapolation_rule_1, + extrapolation_rule_2, + extrapolation_rule_2, + extrapolation_rule_3, + extrapolation_rule_3); + + // Instantiate chunk of coordinates of dom_interpolation + ddc::Chunk coords_eval_alloc(dom_vals, ddc::KokkosAllocator, MemorySpace>()); + ddc::ChunkSpan const coords_eval = coords_eval_alloc.span_view(); + ddc::parallel_for_each( + exec_space, + coords_eval.domain(), + KOKKOS_LAMBDA(DElem const e) { + coords_eval(e) = ddc::coordinate(DElem(e)); + }); + + + // Instantiate chunks to receive outputs of spline_evaluator + ddc::Chunk spline_eval_alloc(dom_vals, ddc::KokkosAllocator()); + ddc::ChunkSpan const spline_eval = spline_eval_alloc.span_view(); + ddc::Chunk spline_eval_deriv1_alloc(dom_vals, ddc::KokkosAllocator()); + ddc::ChunkSpan const spline_eval_deriv1 = spline_eval_deriv1_alloc.span_view(); + ddc::Chunk spline_eval_deriv2_alloc(dom_vals, ddc::KokkosAllocator()); + ddc::ChunkSpan const spline_eval_deriv2 = spline_eval_deriv2_alloc.span_view(); + ddc::Chunk spline_eval_deriv3_alloc(dom_vals, ddc::KokkosAllocator()); + ddc::ChunkSpan const spline_eval_deriv3 = spline_eval_deriv3_alloc.span_view(); + ddc::Chunk spline_eval_deriv12_alloc(dom_vals, ddc::KokkosAllocator()); + ddc::ChunkSpan const spline_eval_deriv12 = spline_eval_deriv12_alloc.span_view(); + ddc::Chunk spline_eval_deriv23_alloc(dom_vals, ddc::KokkosAllocator()); + ddc::ChunkSpan const spline_eval_deriv23 = spline_eval_deriv23_alloc.span_view(); + ddc::Chunk spline_eval_deriv13_alloc(dom_vals, ddc::KokkosAllocator()); + ddc::ChunkSpan const spline_eval_deriv13 = spline_eval_deriv13_alloc.span_view(); + ddc::Chunk spline_eval_deriv123_alloc(dom_vals, ddc::KokkosAllocator()); + ddc::ChunkSpan const spline_eval_deriv123 = spline_eval_deriv123_alloc.span_view(); + + // Call spline_evaluator on the same mesh we started with + spline_evaluator(spline_eval, coords_eval.span_cview(), coef.span_cview()); + spline_evaluator + .template deriv(spline_eval_deriv1, coords_eval.span_cview(), coef.span_cview()); + spline_evaluator + .template deriv(spline_eval_deriv2, coords_eval.span_cview(), coef.span_cview()); + spline_evaluator + .template deriv(spline_eval_deriv3, coords_eval.span_cview(), coef.span_cview()); + spline_evaluator.template deriv2< + I1, + I2>(spline_eval_deriv12, coords_eval.span_cview(), coef.span_cview()); + spline_evaluator.template deriv2< + I2, + I3>(spline_eval_deriv23, coords_eval.span_cview(), coef.span_cview()); + spline_evaluator.template deriv2< + I1, + I3>(spline_eval_deriv13, coords_eval.span_cview(), coef.span_cview()); + spline_evaluator.template deriv3< + I1, + I2, + I3>(spline_eval_deriv123, coords_eval.span_cview(), coef.span_cview()); + + // Checking errors (we recover the initial values) + double const max_norm_error = ddc::parallel_transform_reduce( + exec_space, + spline_eval.domain(), + 0., + ddc::reducer::max(), + KOKKOS_LAMBDA(DElem const e) { + return Kokkos::abs(spline_eval(e) - vals(e)); + }); + double const max_norm_error_diff1 = ddc::parallel_transform_reduce( + exec_space, + spline_eval_deriv1.domain(), + 0., + ddc::reducer::max(), + KOKKOS_LAMBDA(DElem const e) { + Coord const x = ddc::coordinate(DElem(e)); + Coord const y = ddc::coordinate(DElem(e)); + Coord const z = ddc::coordinate(DElem(e)); + return Kokkos::abs(spline_eval_deriv1(e) - evaluator.deriv(x, y, z, 1, 0, 0)); + }); + double const max_norm_error_diff2 = ddc::parallel_transform_reduce( + exec_space, + spline_eval_deriv2.domain(), + 0., + ddc::reducer::max(), + KOKKOS_LAMBDA(DElem const e) { + Coord const x = ddc::coordinate(DElem(e)); + Coord const y = ddc::coordinate(DElem(e)); + Coord const z = ddc::coordinate(DElem(e)); + return Kokkos::abs(spline_eval_deriv2(e) - evaluator.deriv(x, y, z, 0, 1, 0)); + }); + double const max_norm_error_diff3 = ddc::parallel_transform_reduce( + exec_space, + spline_eval_deriv3.domain(), + 0., + ddc::reducer::max(), + KOKKOS_LAMBDA(DElem const e) { + Coord const x = ddc::coordinate(DElem(e)); + Coord const y = ddc::coordinate(DElem(e)); + Coord const z = ddc::coordinate(DElem(e)); + return Kokkos::abs(spline_eval_deriv3(e) - evaluator.deriv(x, y, z, 0, 0, 1)); + }); + double const max_norm_error_diff12 = ddc::parallel_transform_reduce( + exec_space, + spline_eval_deriv12.domain(), + 0., + ddc::reducer::max(), + KOKKOS_LAMBDA(DElem const e) { + Coord const x = ddc::coordinate(DElem(e)); + Coord const y = ddc::coordinate(DElem(e)); + Coord const z = ddc::coordinate(DElem(e)); + return Kokkos::abs(spline_eval_deriv12(e) - evaluator.deriv(x, y, z, 1, 1, 0)); + }); + double const max_norm_error_diff23 = ddc::parallel_transform_reduce( + exec_space, + spline_eval_deriv23.domain(), + 0., + ddc::reducer::max(), + KOKKOS_LAMBDA(DElem const e) { + Coord const x = ddc::coordinate(DElem(e)); + Coord const y = ddc::coordinate(DElem(e)); + Coord const z = ddc::coordinate(DElem(e)); + return Kokkos::abs(spline_eval_deriv23(e) - evaluator.deriv(x, y, z, 0, 1, 1)); + }); + double const max_norm_error_diff13 = ddc::parallel_transform_reduce( + exec_space, + spline_eval_deriv13.domain(), + 0., + ddc::reducer::max(), + KOKKOS_LAMBDA(DElem const e) { + Coord const x = ddc::coordinate(DElem(e)); + Coord const y = ddc::coordinate(DElem(e)); + Coord const z = ddc::coordinate(DElem(e)); + return Kokkos::abs(spline_eval_deriv13(e) - evaluator.deriv(x, y, z, 1, 0, 1)); + }); + double const max_norm_error_diff123 = ddc::parallel_transform_reduce( + exec_space, + spline_eval_deriv123.domain(), + 0., + ddc::reducer::max(), + KOKKOS_LAMBDA(DElem const e) { + Coord const x = ddc::coordinate(DElem(e)); + Coord const y = ddc::coordinate(DElem(e)); + Coord const z = ddc::coordinate(DElem(e)); + return Kokkos::abs(spline_eval_deriv123(e) - evaluator.deriv(x, y, z, 1, 1, 1)); + }); + + double const max_norm = evaluator.max_norm(); + double const max_norm_diff1 = evaluator.max_norm(1, 0, 0); + double const max_norm_diff2 = evaluator.max_norm(0, 1, 0); + double const max_norm_diff3 = evaluator.max_norm(0, 0, 1); + double const max_norm_diff12 = evaluator.max_norm(1, 1, 0); + double const max_norm_diff23 = evaluator.max_norm(0, 1, 1); + double const max_norm_diff13 = evaluator.max_norm(1, 0, 1); + double const max_norm_diff123 = evaluator.max_norm(1, 1, 1); + + SplineErrorBounds> const error_bounds(evaluator); + EXPECT_LE( + max_norm_error, + std:: + max(error_bounds.error_bound( + dx(ncells), + dx(ncells), + dx(ncells), + s_degree, + s_degree, + s_degree), + 1.0e-14 * max_norm)); + EXPECT_LE( + max_norm_error_diff1, + std:: + max(error_bounds.error_bound_on_deriv_1( + dx(ncells), + dx(ncells), + dx(ncells), + s_degree, + s_degree, + s_degree), + 1e-12 * max_norm_diff1)); + EXPECT_LE( + max_norm_error_diff2, + std:: + max(error_bounds.error_bound_on_deriv_2( + dx(ncells), + dx(ncells), + dx(ncells), + s_degree, + s_degree, + s_degree), + 1e-12 * max_norm_diff2)); + EXPECT_LE( + max_norm_error_diff3, + std:: + max(error_bounds.error_bound_on_deriv_3( + dx(ncells), + dx(ncells), + dx(ncells), + s_degree, + s_degree, + s_degree), + 1e-12 * max_norm_diff3)); + EXPECT_LE( + max_norm_error_diff12, + std:: + max(error_bounds.error_bound_on_deriv_12( + dx(ncells), + dx(ncells), + dx(ncells), + s_degree, + s_degree, + s_degree), + 1e-11 * max_norm_diff12)); + EXPECT_LE( + max_norm_error_diff23, + std:: + max(error_bounds.error_bound_on_deriv_23( + dx(ncells), + dx(ncells), + dx(ncells), + s_degree, + s_degree, + s_degree), + 1e-11 * max_norm_diff23)); + EXPECT_LE( + max_norm_error_diff13, + std:: + max(error_bounds.error_bound_on_deriv_13( + dx(ncells), + dx(ncells), + dx(ncells), + s_degree, + s_degree, + s_degree), + 1e-11 * max_norm_diff13)); + EXPECT_LE( + max_norm_error_diff123, + std:: + max(error_bounds.error_bound_on_deriv_123( + dx(ncells), + dx(ncells), + dx(ncells), + s_degree, + s_degree, + s_degree), + 1e-11 * max_norm_diff123)); +} + +} // namespace anonymous_namespace_workaround_batched_3d_spline_builder_cpp + +#if defined(BC_PERIODIC) && defined(BSPLINES_TYPE_UNIFORM) +# define SUFFIX(name) name##Periodic##Uniform +#elif defined(BC_PERIODIC) && defined(BSPLINES_TYPE_NON_UNIFORM) +# define SUFFIX(name) name##Periodic##NonUniform +#elif defined(BC_GREVILLE) && defined(BSPLINES_TYPE_UNIFORM) +# define SUFFIX(name) name##Greville##Uniform +#elif defined(BC_GREVILLE) && defined(BSPLINES_TYPE_NON_UNIFORM) +# define SUFFIX(name) name##Greville##NonUniform +#elif defined(BC_HERMITE) && defined(BSPLINES_TYPE_UNIFORM) +# define SUFFIX(name) name##Hermite##Uniform +#elif defined(BC_HERMITE) && defined(BSPLINES_TYPE_NON_UNIFORM) +# define SUFFIX(name) name##Hermite##NonUniform +#endif + +TEST(SUFFIX(Batched3dSplineHost), 3DXYZ) +{ + Batched3dSplineTest< + Kokkos::DefaultHostExecutionSpace, + Kokkos::DefaultHostExecutionSpace::memory_space, + DDimGPS, + DDimGPS, + DDimGPS, + DDimGPS, + DDimGPS, + DDimGPS>(); +} + +TEST(SUFFIX(Batched3dSplineDevice), 3DXYZ) +{ + Batched3dSplineTest< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + DDimGPS, + DDimGPS, + DDimGPS, + DDimGPS, + DDimGPS, + DDimGPS>(); +} + +TEST(SUFFIX(Batched3dSplineHost), 4DXYZB) +{ + Batched3dSplineTest< + Kokkos::DefaultHostExecutionSpace, + Kokkos::DefaultHostExecutionSpace::memory_space, + DDimGPS, + DDimGPS, + DDimGPS, + DDimGPS, + DDimGPS, + DDimGPS, + DDimBatch>(); +} + +TEST(SUFFIX(Batched3dSplineHost), 4DXBYZ) +{ + Batched3dSplineTest< + Kokkos::DefaultHostExecutionSpace, + Kokkos::DefaultHostExecutionSpace::memory_space, + DDimGPS, + DDimGPS, + DDimGPS, + DDimGPS, + DDimBatch, + DDimGPS, + DDimGPS>(); +} + +TEST(SUFFIX(Batched3dSplineHost), 4DXYBZ) +{ + Batched3dSplineTest< + Kokkos::DefaultHostExecutionSpace, + Kokkos::DefaultHostExecutionSpace::memory_space, + DDimGPS, + DDimGPS, + DDimGPS, + DDimGPS, + DDimGPS, + DDimBatch, + DDimGPS>(); +} + +TEST(SUFFIX(Batched3dSplineHost), 4DBXYZ) +{ + Batched3dSplineTest< + Kokkos::DefaultHostExecutionSpace, + Kokkos::DefaultHostExecutionSpace::memory_space, + DDimGPS, + DDimGPS, + DDimGPS, + DDimBatch, + DDimGPS, + DDimGPS, + DDimGPS>(); +} diff --git a/tests/splines/evaluator_3d.hpp b/tests/splines/evaluator_3d.hpp new file mode 100644 index 000000000..80c16db05 --- /dev/null +++ b/tests/splines/evaluator_3d.hpp @@ -0,0 +1,109 @@ +// Copyright (C) The DDC development team, see COPYRIGHT.md file +// +// SPDX-License-Identifier: MIT + +#pragma once + +#include + +struct Evaluator3D +{ + template + class Evaluator + { + private: + Eval1 eval_func1; + Eval2 eval_func2; + Eval3 eval_func3; + + public: + template + explicit Evaluator(Domain domain) + : eval_func1(ddc::DiscreteDomain(domain)) + , eval_func2(ddc::DiscreteDomain(domain)) + , eval_func3(ddc::DiscreteDomain(domain)) + { + } + + KOKKOS_FUNCTION double operator()(double const x, double const y, double const z) + const noexcept + { + return eval_func1(x) * eval_func2(y) * eval_func3(z); + } + + template + KOKKOS_FUNCTION double operator()( + ddc::Coordinate const x) const noexcept + { + return eval_func1(ddc::get(x)) * eval_func2(ddc::get(x)) + * eval_func3(ddc::get(x)); + } + + template + void operator()( + ddc::ChunkSpan> chunk) const + { + ddc::DiscreteDomain const domain = chunk.domain(); + + for (ddc::DiscreteElement const i : ddc::DiscreteDomain(domain)) { + for (ddc::DiscreteElement const j : ddc::DiscreteDomain(domain)) { + for (ddc::DiscreteElement const k : ddc::DiscreteDomain(domain)) { + chunk(i, j, k) = eval_func1(ddc::coordinate(i)) + * eval_func2(ddc::coordinate(j)) + * eval_func3(ddc::coordinate(k)); + } + } + } + } + + KOKKOS_FUNCTION double deriv( + double const x, + double const y, + double const z, + int const derivative_x, + int const derivative_y, + int const derivative_z) const noexcept + { + return eval_func1.deriv(x, derivative_x) * eval_func2.deriv(y, derivative_y) + * eval_func3.deriv(z, derivative_z); + } + + template + KOKKOS_FUNCTION double deriv( + ddc::Coordinate const x, + int const derivative_x, + int const derivative_y, + int const derivative_z) const noexcept + { + return eval_func1.deriv(ddc::get(x), derivative_x) + * eval_func2.deriv(ddc::get(x), derivative_y) + * eval_func3.deriv(ddc::get(x), derivative_z); + } + + template + void deriv( + ddc::ChunkSpan> chunk, + int const derivative_x, + int const derivative_y, + int const derivative_z) const + { + ddc::DiscreteDomain const domain = chunk.domain(); + + for (ddc::DiscreteElement const i : ddc::DiscreteDomain(domain)) { + for (ddc::DiscreteElement const j : ddc::DiscreteDomain(domain)) { + for (ddc::DiscreteElement const k : ddc::DiscreteDomain(domain)) { + chunk(i, j, k) = eval_func1.deriv(ddc::coordinate(i), derivative_x) + * eval_func2.deriv(ddc::coordinate(j), derivative_y) + * eval_func3.deriv(ddc::coordinate(j), derivative_z); + } + } + } + } + + KOKKOS_FUNCTION double max_norm(int diff1 = 0, int diff2 = 0, int diff3 = 0) const + { + return eval_func1.max_norm(diff1) * eval_func2.max_norm(diff2) + * eval_func3.max_norm(diff3); + } + }; +}; diff --git a/tests/splines/spline_error_bounds.hpp b/tests/splines/spline_error_bounds.hpp index d957fd4c6..0ee41f2ce 100644 --- a/tests/splines/spline_error_bounds.hpp +++ b/tests/splines/spline_error_bounds.hpp @@ -64,6 +64,22 @@ class SplineErrorBounds + tihomirov_error_bound(cell_width2, degree2, norm2); } + double error_bound( + double cell_width1, + double cell_width2, + double cell_width3, + int degree1, + int degree2, + int degree3) const + { + double const norm1 = m_evaluator.max_norm(degree1 + 1, 0, 0); + double const norm2 = m_evaluator.max_norm(0, degree2 + 1, 0); + double const norm3 = m_evaluator.max_norm(0, 0, degree3 + 1); + return tihomirov_error_bound(cell_width1, degree1, norm1) + + tihomirov_error_bound(cell_width2, degree2, norm2) + + tihomirov_error_bound(cell_width3, degree3, norm3); + } + double error_bound_on_deriv(double cell_width, int degree) const { return tihomirov_error_bound(cell_width, degree - 1, m_evaluator.max_norm(degree + 1)); @@ -78,6 +94,22 @@ class SplineErrorBounds + tihomirov_error_bound(cell_width2, degree2, norm2); } + double error_bound_on_deriv_1( + double cell_width1, + double cell_width2, + double cell_width3, + int degree1, + int degree2, + int degree3) const + { + double const norm1 = m_evaluator.max_norm(degree1 + 1, 0, 0); + double const norm2 = m_evaluator.max_norm(0, degree2 + 1, 0); + double const norm3 = m_evaluator.max_norm(0, 0, degree3 + 1); + return tihomirov_error_bound(cell_width1, degree1 - 1, norm1) + + tihomirov_error_bound(cell_width2, degree2, norm2) + + tihomirov_error_bound(cell_width3, degree3, norm3); + } + double error_bound_on_deriv_2(double cell_width1, double cell_width2, int degree1, int degree2) const { @@ -87,6 +119,38 @@ class SplineErrorBounds + tihomirov_error_bound(cell_width2, degree2 - 1, norm2); } + double error_bound_on_deriv_2( + double cell_width1, + double cell_width2, + double cell_width3, + int degree1, + int degree2, + int degree3) const + { + double const norm1 = m_evaluator.max_norm(degree1 + 1, 0, 0); + double const norm2 = m_evaluator.max_norm(0, degree2 + 1, 0); + double const norm3 = m_evaluator.max_norm(0, 0, degree3 + 1); + return tihomirov_error_bound(cell_width1, degree1, norm1) + + tihomirov_error_bound(cell_width2, degree2 - 1, norm2) + + tihomirov_error_bound(cell_width3, degree3, norm3); + } + + double error_bound_on_deriv_3( + double cell_width1, + double cell_width2, + double cell_width3, + int degree1, + int degree2, + int degree3) const + { + double const norm1 = m_evaluator.max_norm(degree1 + 1, 0, 0); + double const norm2 = m_evaluator.max_norm(0, degree2 + 1, 0); + double const norm3 = m_evaluator.max_norm(0, 0, degree3 + 1); + return tihomirov_error_bound(cell_width1, degree1, norm1) + + tihomirov_error_bound(cell_width2, degree2, norm2) + + tihomirov_error_bound(cell_width3, degree3 - 1, norm3); + } + /******************************************************************************* * NOTE: The following estimates have no theoretical justification but capture * the correct asympthotic rate of convergence. @@ -101,6 +165,70 @@ class SplineErrorBounds + tihomirov_error_bound(cell_width2, degree2 - 1, norm2); } + double error_bound_on_deriv_12( + double cell_width1, + double cell_width2, + double cell_width3, + int degree1, + int degree2, + int degree3) const + { + double const norm1 = m_evaluator.max_norm(degree1 + 1, 1, 1); + double const norm2 = m_evaluator.max_norm(1, degree2 + 1, 1); + double const norm3 = m_evaluator.max_norm(1, 1, degree3 + 1); + return tihomirov_error_bound(cell_width1, degree1 - 1, norm1) + + tihomirov_error_bound(cell_width2, degree2 - 1, norm2) + + tihomirov_error_bound(cell_width3, degree3, norm3); + } + + double error_bound_on_deriv_23( + double cell_width1, + double cell_width2, + double cell_width3, + int degree1, + int degree2, + int degree3) const + { + double const norm1 = m_evaluator.max_norm(degree1 + 1, 1, 1); + double const norm2 = m_evaluator.max_norm(1, degree2 + 1, 1); + double const norm3 = m_evaluator.max_norm(1, 1, degree3 + 1); + return tihomirov_error_bound(cell_width1, degree1, norm1) + + tihomirov_error_bound(cell_width2, degree2 - 1, norm2) + + tihomirov_error_bound(cell_width3, degree3 - 3, norm3); + } + + double error_bound_on_deriv_13( + double cell_width1, + double cell_width2, + double cell_width3, + int degree1, + int degree2, + int degree3) const + { + double const norm1 = m_evaluator.max_norm(degree1 + 1, 1, 1); + double const norm2 = m_evaluator.max_norm(1, degree2 + 1, 1); + double const norm3 = m_evaluator.max_norm(1, 1, degree3 + 1); + return tihomirov_error_bound(cell_width1, degree1 - 1, norm1) + + tihomirov_error_bound(cell_width2, degree2, norm2) + + tihomirov_error_bound(cell_width3, degree3 - 1, norm3); + } + + double error_bound_on_deriv_123( + double cell_width1, + double cell_width2, + double cell_width3, + int degree1, + int degree2, + int degree3) const + { + double const norm1 = m_evaluator.max_norm(degree1 + 1, 1, 1); + double const norm2 = m_evaluator.max_norm(1, degree2 + 1, 1); + double const norm3 = m_evaluator.max_norm(1, 1, degree3 + 1); + return tihomirov_error_bound(cell_width1, degree1 - 1, norm1) + + tihomirov_error_bound(cell_width2, degree2 - 1, norm2) + + tihomirov_error_bound(cell_width3, degree3 - 1, norm3); + } + double error_bound_on_int(double cell_width, int degree) const { return tihomirov_error_bound(cell_width, degree + 1, m_evaluator.max_norm(degree + 1)); From 86c827ef215b1052acc1336f2b7fc0bfe55a55e3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tr=C3=A9vis=20Morvany?= Date: Mon, 29 Sep 2025 22:27:38 +0200 Subject: [PATCH 2/9] Changes to the derivs domain - rename some of the type aliases - use only one deriv dimension in the derivs_domain type aliases - enable the greville tests - lower the norm threshold for the triple derivative to 1e-10 --- .../ddc/kernels/splines/spline_builder_3d.hpp | 199 ++++++++---------- tests/splines/CMakeLists.txt | 2 +- tests/splines/batched_3d_spline_builder.cpp | 2 +- 3 files changed, 86 insertions(+), 117 deletions(-) diff --git a/include/ddc/kernels/splines/spline_builder_3d.hpp b/include/ddc/kernels/splines/spline_builder_3d.hpp index a3f88ea1c..6c2ba818e 100644 --- a/include/ddc/kernels/splines/spline_builder_3d.hpp +++ b/include/ddc/kernels/splines/spline_builder_3d.hpp @@ -61,20 +61,12 @@ class SplineBuilder3D using builder_type3 = ddc:: SplineBuilder; - // FIXME: change this when implementing the derivatives - /// @brief The type of the SplineBuilder used by this class to spline-approximate the second-dimension-derivatives along first dimension. - using builder_deriv_type1 = ddc::SplineBuilder2D< - ExecSpace, - MemorySpace, - BSpline1, - BSpline2, - DDimI1, - DDimI2, - BcLower1, - BcUpper1, - BcLower2, - BcLower2, - Solver>; + // TODO: add documentation comment + using builder_deriv_type1 = ddc:: + SplineBuilder; + + using builder_deriv_type2 = ddc:: + SplineBuilder; /// @brief The type of the first interpolation continuous dimension. using continuous_dimension_type1 = typename builder_type1::continuous_dimension_type; @@ -191,44 +183,11 @@ class SplineBuilder3D * Example: For batched_interpolation_domain_type = DiscreteDomain and dimensions of interest X, Y and Z, * this is DiscreteDomain, Y, Z, T>. */ - // template < - // class BatchedInterpolationDDom, - // class = std::enable_if_t>> - // using batched_derivs_domain_type1 = - // typename builder_type1::template batched_derivs_domain_type; - template < class BatchedInterpolationDDom, class = std::enable_if_t>> - using batched_derivs_domain_type_1_2 = ddc::replace_dim_of_t< - ddc::replace_dim_of_t< - BatchedInterpolationDDom, - interpolation_discrete_dimension_type1, - deriv_type1>, - interpolation_discrete_dimension_type2, - deriv_type2>; - - template < - class BatchedInterpolationDDom, - class = std::enable_if_t>> - using batched_derivs_domain_type_2_3 = ddc::replace_dim_of_t< - ddc::replace_dim_of_t< - BatchedInterpolationDDom, - interpolation_discrete_dimension_type2, - deriv_type2>, - interpolation_discrete_dimension_type3, - deriv_type3>; - - template < - class BatchedInterpolationDDom, - class = std::enable_if_t>> - using batched_derivs_domain_type_1_3 = ddc::replace_dim_of_t< - ddc::replace_dim_of_t< - BatchedInterpolationDDom, - interpolation_discrete_dimension_type3, - deriv_type3>, - interpolation_discrete_dimension_type3, - deriv_type3>; + using batched_derivs_domain_type1 = + typename builder_type1::template batched_derivs_domain_type; /** * @brief The type of the whole Derivs domain (cartesian product of the 1D Deriv domain @@ -239,13 +198,13 @@ class SplineBuilder3D * Example: For batched_interpolation_domain_type = DiscreteDomain and dimensions of interest X, Y, and Z * this is DiscreteDomain, Z, T>. */ - // template < - // class BatchedInterpolationDDom, - // class = std::enable_if_t>> - // using batched_derivs_domain_type2 = ddc::replace_dim_of_t< - // BatchedInterpolationDDom, - // interpolation_discrete_dimension_type2, - // deriv_type2>; + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using batched_derivs_domain_type2 = ddc::replace_dim_of_t< + BatchedInterpolationDDom, + interpolation_discrete_dimension_type2, + deriv_type2>; /** * @brief The type of the whole Derivs domain (cartesian product of the 1D Deriv domain @@ -256,13 +215,13 @@ class SplineBuilder3D * Example: For batched_interpolation_domain_type = DiscreteDomain and dimensions of interest X, Y, and Z * this is DiscreteDomain, T>. */ - // template < - // class BatchedInterpolationDDom, - // class = std::enable_if_t>> - // using batched_derivs_domain_type3 = ddc::replace_dim_of_t< - // BatchedInterpolationDDom, - // interpolation_discrete_dimension_type3, - // deriv_type3>; + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using batched_derivs_domain_type3 = ddc::replace_dim_of_t< + BatchedInterpolationDDom, + interpolation_discrete_dimension_type3, + deriv_type3>; /** * @brief The type of the whole Derivs domain (cartesian product of the 3D Deriv domain @@ -287,13 +246,14 @@ class SplineBuilder3D private: builder_type1 m_spline_builder1; - builder_deriv_type1 m_spline_builder_deriv1; builder_type2 m_spline_builder2; builder_type3 m_spline_builder3; + builder_deriv_type1 m_spline_builder_deriv1; + builder_deriv_type2 m_spline_builder_deriv2; public: /** - * @brief Build a SplineBuilder2D acting on interpolation_domain. + * @brief Build a SplineBuilder3D acting on interpolation_domain. * * @param interpolation_domain The domain on which the interpolation points are defined, without the batch dimensions. * @@ -313,14 +273,15 @@ class SplineBuilder3D std::optional cols_per_chunk = std::nullopt, std::optional preconditioner_max_block_size = std::nullopt) : m_spline_builder1(interpolation_domain, cols_per_chunk, preconditioner_max_block_size) - , m_spline_builder_deriv1(interpolation_domain) , m_spline_builder2(interpolation_domain, cols_per_chunk, preconditioner_max_block_size) , m_spline_builder3(interpolation_domain, cols_per_chunk, preconditioner_max_block_size) + , m_spline_builder_deriv1(interpolation_domain) + , m_spline_builder_deriv2(interpolation_domain) { } /** - * @brief Build a SplineBuilder2D acting on the interpolation domain contained in batched_interpolation_domain. + * @brief Build a SplineBuilder3D acting on the interpolation domain contained in batched_interpolation_domain. * * @param batched_interpolation_domain The domain on which the interpolation points are defined. * @@ -355,7 +316,7 @@ class SplineBuilder3D /** * @brief Move-constructs. * - * @param x An rvalue to another SplineBuilder2D. + * @param x An rvalue to another SplineBuilder3D. */ SplineBuilder3D(SplineBuilder3D&& x) = default; @@ -373,18 +334,18 @@ class SplineBuilder3D SplineBuilder3D& operator=(SplineBuilder3D&& x) = default; /** - * @brief Get the domain for the 2D interpolation mesh used by this class. + * @brief Get the domain for the 3D interpolation mesh used by this class. * - * This is 2D because it is defined along the dimensions of interest. + * This is 3D because it is defined along the dimensions of interest. * - * @return The 2D domain for the interpolation mesh. + * @return The 3D domain for the interpolation mesh. */ interpolation_domain_type interpolation_domain() const noexcept { return ddc::DiscreteDomain< interpolation_domain_type1, interpolation_domain_type2, - interpolation_discrete_dimension_type3>( + interpolation_domain_type3>( m_spline_builder1.interpolation_domain(), m_spline_builder2.interpolation_domain(), m_spline_builder3.interpolation_domain()); @@ -394,7 +355,7 @@ class SplineBuilder3D * @brief Get the whole domain representing interpolation points. * * Values of the function must be provided on this domain in order - * to build a spline representation of the function (cartesian product of 2D interpolation_domain and batch_domain). + * to build a spline representation of the function (cartesian product of 3D interpolation_domain and batch_domain). * * @param batched_interpolation_domain The whole domain on which the interpolation points are defined. * @@ -430,11 +391,11 @@ class SplineBuilder3D } /** - * @brief Get the 2D domain on which spline coefficients are defined. + * @brief Get the 3D domain on which spline coefficients are defined. * - * The 2D spline domain corresponding to the dimensions of interest. + * The 3D spline domain corresponding to the dimensions of interest. * - * @return The 2D domain for the spline coefficients. + * @return The 3D domain for the spline coefficients. */ ddc::DiscreteDomain spline_domain() const noexcept @@ -471,12 +432,12 @@ class SplineBuilder3D } /** - * @brief Compute a 2D spline approximation of a function. + * @brief Compute a 3D spline approximation of a function. * * Use the values of a function (defined on - * SplineBuilder2D::batched_interpolation_domain) and the derivatives of the + * SplineBuilder3D::batched_interpolation_domain) and the derivatives of the * function at the boundaries (in the case of BoundCond::HERMITE only) - * to calculate a 2D spline approximation of this function. + * to calculate a 3D spline approximation of this function. * * The spline approximation is stored as a ChunkSpan of coefficients * associated with B-splines. @@ -493,18 +454,26 @@ class SplineBuilder3D * The values of the derivatives at the lower boundary in the second dimension. * @param[in] derivs_max2 * The values of the derivatives at the upper boundary in the second dimension. - * @param[in] mixed_derivs_min1_min2 - * The values of the the cross-derivatives at the lower boundary in the first dimension - * and the lower boundary in the second dimension. - * @param[in] mixed_derivs_max1_min2 - * The values of the the cross-derivatives at the upper boundary in the first dimension - * and the lower boundary in the second dimension. - * @param[in] mixed_derivs_min1_max2 - * The values of the the cross-derivatives at the lower boundary in the first dimension - * and the upper boundary in the second dimension. - * @param[in] mixed_derivs_max1_max2 - * The values of the the cross-derivatives at the upper boundary in the first dimension - * and the upper boundary in the second dimension. + * @param[in] derivs_min3 + * The values of the derivatives at the lower boundary in the third dimension. + * @param[in] derivs_max3 + * The values of the derivatives at the upper boundary in the third dimension. + * @param[in] mixed_derivs_min1_min2_min3 + * The values of the the cross-derivatives at the lower boundary in the first dimension, the lower boundary in the second dimension and the lower boundary in the third dimension. + * @param[in] mixed_derivs_max1_min2_min3 + * The values of the the cross-derivatives at the upper boundary in the first dimension, the lower boundary in the second dimension and the lower boundary in the third dimension. + * @param[in] mixed_derivs_min1_max2_min3 + * The values of the the cross-derivatives at the lower boundary in the first dimension, the upper boundary in the second dimension and the lower boundary in the third dimension. + * @param[in] mixed_derivs_max1_max2_min3 + * The values of the the cross-derivatives at the upper boundary in the first dimension, the upper boundary in the second dimension and the lower boundary in the third dimension. + * @param[in] mixed_derivs_min1_min2_max3 + * The values of the the cross-derivatives at the lower boundary in the first dimension, the lower boundary in the second dimension and the upper boundary in the third dimension. + * @param[in] mixed_derivs_max1_min2_max3 + * The values of the the cross-derivatives at the upper boundary in the first dimension, the lower boundary in the second dimension and the upper boundary in the third dimension. + * @param[in] mixed_derivs_min1_max2_max3 + * The values of the the cross-derivatives at the lower boundary in the first dimension, the upper boundary in the second dimension and the upper boundary in the third dimension. + * @param[in] mixed_derivs_max1_max2_max3 + * The values of the the cross-derivatives at the upper boundary in the first dimension, the upper boundary in the second dimension and the upper boundary in the third dimension. */ template void operator()( @@ -516,39 +485,39 @@ class SplineBuilder3D ddc::ChunkSpan vals, std::optional, + batched_derivs_domain_type1, Layout, - memory_space>> derivs_min_1_2 + memory_space>> derivs_min1 = std::nullopt, std::optional, + batched_derivs_domain_type1, Layout, - memory_space>> derivs_max_1_2 + memory_space>> derivs_max1 = std::nullopt, std::optional, + batched_derivs_domain_type2, Layout, - memory_space>> derivs_min_2_3 + memory_space>> derivs_min2 = std::nullopt, std::optional, + batched_derivs_domain_type2, Layout, - memory_space>> derivs_max_2_3 + memory_space>> derivs_max2 = std::nullopt, std::optional, + batched_derivs_domain_type3, Layout, - memory_space>> derivs_min_1_3 + memory_space>> derivs_min3 = std::nullopt, std::optional, + batched_derivs_domain_type3, Layout, - memory_space>> derivs_max_1_3 + memory_space>> derivs_max3 = std::nullopt, std::optional vals, std::optional, + batched_derivs_domain_type1, Layout, - memory_space>> derivs_min_1_2, + memory_space>> derivs_min1, std::optional, + batched_derivs_domain_type1, Layout, - memory_space>> derivs_max_1_2, + memory_space>> derivs_max1, std::optional, + batched_derivs_domain_type2, Layout, - memory_space>> derivs_min_2_3, + memory_space>> derivs_min2, std::optional, + batched_derivs_domain_type2, Layout, - memory_space>> derivs_max_2_3, + memory_space>> derivs_max2, std::optional, + batched_derivs_domain_type3, Layout, - memory_space>> derivs_min_1_3, + memory_space>> derivs_min3, std::optional, + batched_derivs_domain_type3, Layout, - memory_space>> derivs_max_1_3, + memory_space>> derivs_max3, std::optional, diff --git a/tests/splines/CMakeLists.txt b/tests/splines/CMakeLists.txt index a328c1251..abd6c1fe6 100644 --- a/tests/splines/CMakeLists.txt +++ b/tests/splines/CMakeLists.txt @@ -144,7 +144,7 @@ foreach(BC "BC_PERIODIC" "BC_GREVILLE" "BC_HERMITE") endforeach() endforeach() -foreach(BC "BC_PERIODIC") #"BC_GREVILLE" "BC_HERMITE") +foreach(BC "BC_PERIODIC" "BC_GREVILLE")# "BC_HERMITE") foreach(EVALUATOR "EVALUATOR_POLYNOMIAL") foreach(DEGREE RANGE "${DDC_SPLINES_TEST_DEGREE_MIN}" "${DDC_SPLINES_TEST_DEGREE_MAX}") foreach(BSPLINES_TYPE "BSPLINES_TYPE_UNIFORM" "BSPLINES_TYPE_NON_UNIFORM") diff --git a/tests/splines/batched_3d_spline_builder.cpp b/tests/splines/batched_3d_spline_builder.cpp index 9f3c95b93..1bb3e7129 100644 --- a/tests/splines/batched_3d_spline_builder.cpp +++ b/tests/splines/batched_3d_spline_builder.cpp @@ -750,7 +750,7 @@ void Batched3dSplineTest() s_degree, s_degree, s_degree), - 1e-11 * max_norm_diff123)); + 1e-10 * max_norm_diff123)); } } // namespace anonymous_namespace_workaround_batched_3d_spline_builder_cpp From 1aae704f6f32a1d8fb531ef177ef280d0eb20c30 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tr=C3=A9vis=20Morvany?= Date: Tue, 30 Sep 2025 10:33:02 +0200 Subject: [PATCH 3/9] Reword the FIXMEs in the 3D test --- .../ddc/kernels/splines/spline_builder_3d.hpp | 1 - tests/splines/batched_3d_spline_builder.cpp | 19 ++++++++----------- 2 files changed, 8 insertions(+), 12 deletions(-) diff --git a/include/ddc/kernels/splines/spline_builder_3d.hpp b/include/ddc/kernels/splines/spline_builder_3d.hpp index 6c2ba818e..72b9c3510 100644 --- a/include/ddc/kernels/splines/spline_builder_3d.hpp +++ b/include/ddc/kernels/splines/spline_builder_3d.hpp @@ -12,7 +12,6 @@ #include #include "spline_builder.hpp" -#include "spline_builder_2d.hpp" namespace ddc { diff --git a/tests/splines/batched_3d_spline_builder.cpp b/tests/splines/batched_3d_spline_builder.cpp index 1bb3e7129..c47ad36a3 100644 --- a/tests/splines/batched_3d_spline_builder.cpp +++ b/tests/splines/batched_3d_spline_builder.cpp @@ -209,24 +209,21 @@ void Batched3dSplineTest() interpolation_domain3); #if defined(BC_HERMITE) + // FIXME: This is the 2D version, modify when adding the Hermite boundary conditions + // Create the derivs domain ddc::DiscreteDomain> const derivs_domain1(DElem>(1), DVect>(s_degree / 2)); ddc::DiscreteDomain> const derivs_domain2(DElem>(1), DVect>(s_degree / 2)); - ddc::DiscreteDomain> const - derivs_domain3(DElem>(1), DVect>(s_degree / 2)); - ddc::DiscreteDomain, ddc::Deriv, ddc::Deriv> const - derivs_domain(derivs_domain1, derivs_domain2, derivs_domain3); + ddc::DiscreteDomain, ddc::Deriv> const + derivs_domain(derivs_domain1, derivs_domain2); auto const dom_derivs_1d = ddc::replace_dim_of>(dom_vals, derivs_domain1); auto const dom_derivs2 = ddc::replace_dim_of>(dom_vals, derivs_domain2); - auto const dom_derivs3 = ddc::replace_dim_of>(dom_vals, derivs_domain3); - auto const dom_derivs = ddc::replace_dim_of>( - ddc::replace_dim_of>(dom_derivs_1d, derivs_domain2), - derivs_domain3); - // FIXME: check later the correctness of this + auto const dom_derivs + = ddc::replace_dim_of>(dom_derivs_1d, derivs_domain2); #endif // Create a SplineBuilder over BSplines and batched along other dimensions using some boundary conditions @@ -270,7 +267,7 @@ void Batched3dSplineTest() }); #if defined(BC_HERMITE) - // FIXME: handle this later + // FIXME: This is the 2D version, modify when adding the Hermite boundary conditions // Allocate and fill a chunk containing derivs to be passed as input to spline_builder. int const shift = s_degree % 2; // shift = 0 for even order, 1 for odd order @@ -462,7 +459,7 @@ void Batched3dSplineTest() // Finally compute the spline by filling `coef` #if defined(BC_HERMITE) - // FIXME: handle this later + // FIXME: This is the 2D version, modify when adding the Hermite boundary conditions spline_builder( coef, vals.span_cview(), From bf46583a86ee404cafe3fb0637af629ff20ea2ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tr=C3=A9vis=20Morvany?= Date: Tue, 30 Sep 2025 10:36:44 +0200 Subject: [PATCH 4/9] Fix a typo --- include/ddc/kernels/splines/spline_evaluator_3d.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/ddc/kernels/splines/spline_evaluator_3d.hpp b/include/ddc/kernels/splines/spline_evaluator_3d.hpp index de0692439..83006ee0f 100644 --- a/include/ddc/kernels/splines/spline_evaluator_3d.hpp +++ b/include/ddc/kernels/splines/spline_evaluator_3d.hpp @@ -310,8 +310,8 @@ class SplineEvaluator3D * @param upper_extrap_rule1 The extrapolation rule at the upper boundary along the first dimension. * @param lower_extrap_rule2 The extrapolation rule at the lower boundary along the second dimension. * @param upper_extrap_rule2 The extrapolation rule at the upper boundary along the second dimension. - * @param lower_extrap_rule2 The extrapolation rule at the lower boundary along the third dimension. - * @param upper_extrap_rule2 The extrapolation rule at the upper boundary along the third dimension. + * @param lower_extrap_rule3 The extrapolation rule at the lower boundary along the third dimension. + * @param upper_extrap_rule3 The extrapolation rule at the upper boundary along the third dimension. * * @see NullExtrapolationRule ConstantExtrapolationRule PeriodicExtrapolationRule */ From 423b8f7e502ccc641e83fb3ae1acd9f5453f0a07 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tr=C3=A9vis=20Morvany?= Date: Tue, 30 Sep 2025 10:44:04 +0200 Subject: [PATCH 5/9] Format the splines CMakeLists.txt --- tests/splines/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/splines/CMakeLists.txt b/tests/splines/CMakeLists.txt index abd6c1fe6..ca68b1a96 100644 --- a/tests/splines/CMakeLists.txt +++ b/tests/splines/CMakeLists.txt @@ -144,7 +144,7 @@ foreach(BC "BC_PERIODIC" "BC_GREVILLE" "BC_HERMITE") endforeach() endforeach() -foreach(BC "BC_PERIODIC" "BC_GREVILLE")# "BC_HERMITE") +foreach(BC "BC_PERIODIC" "BC_GREVILLE") # "BC_HERMITE") foreach(EVALUATOR "EVALUATOR_POLYNOMIAL") foreach(DEGREE RANGE "${DDC_SPLINES_TEST_DEGREE_MIN}" "${DDC_SPLINES_TEST_DEGREE_MAX}") foreach(BSPLINES_TYPE "BSPLINES_TYPE_UNIFORM" "BSPLINES_TYPE_NON_UNIFORM") From 01bc97b37ce79244738742830c10ae9bcd639388 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tr=C3=A9vis=20Morvany?= Date: Tue, 30 Sep 2025 11:21:12 +0200 Subject: [PATCH 6/9] Add missing doc comments and add unused annotations to the derivs parameters --- .../ddc/kernels/splines/spline_builder_3d.hpp | 31 ++++++++++--------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/include/ddc/kernels/splines/spline_builder_3d.hpp b/include/ddc/kernels/splines/spline_builder_3d.hpp index 72b9c3510..10f39392b 100644 --- a/include/ddc/kernels/splines/spline_builder_3d.hpp +++ b/include/ddc/kernels/splines/spline_builder_3d.hpp @@ -60,10 +60,11 @@ class SplineBuilder3D using builder_type3 = ddc:: SplineBuilder; - // TODO: add documentation comment + /// @brief The type of the SplineBuilder used by this class to spline-approximate the second-dimension-derivatives along first dimension. using builder_deriv_type1 = ddc:: SplineBuilder; + /// @brief The type of the SplineBuilder used by this class to spline-approximate the third-dimension-derivatives along second dimension. using builder_deriv_type2 = ddc:: SplineBuilder; @@ -609,72 +610,72 @@ operator()( Layout, memory_space> spline, ddc::ChunkSpan vals, - std::optional, Layout, memory_space>> derivs_min1, - std::optional, Layout, memory_space>> derivs_max1, - std::optional, Layout, memory_space>> derivs_min2, - std::optional, Layout, memory_space>> derivs_max2, - std::optional, Layout, memory_space>> derivs_min3, - std::optional, Layout, memory_space>> derivs_max3, - std::optional, Layout, memory_space>> mixed_derivs_min1_min2_min3, - std::optional, Layout, memory_space>> mixed_derivs_max1_min2_min3, - std::optional, Layout, memory_space>> mixed_derivs_min1_max2_min3, - std::optional, Layout, memory_space>> mixed_derivs_max1_max2_min3, - std::optional, Layout, memory_space>> mixed_derivs_min1_min2_max3, - std::optional, Layout, memory_space>> mixed_derivs_max1_min2_max3, - std::optional, Layout, memory_space>> mixed_derivs_min1_max2_max3, - std::optional, Layout, From e69153b79a21e1e54ad5e0d5df77c4572761d6ee Mon Sep 17 00:00:00 2001 From: Thomas Padioleau Date: Tue, 30 Sep 2025 11:49:33 +0200 Subject: [PATCH 7/9] Declare new header files to CMake --- CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index ab426fc06..361c68ee7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -279,8 +279,10 @@ if("${DDC_BUILD_KERNELS_SPLINES}") include/ddc/kernels/splines/spline_boundary_conditions.hpp include/ddc/kernels/splines/spline_builder.hpp include/ddc/kernels/splines/spline_builder_2d.hpp + include/ddc/kernels/splines/spline_builder_3d.hpp include/ddc/kernels/splines/spline_evaluator.hpp include/ddc/kernels/splines/spline_evaluator_2d.hpp + include/ddc/kernels/splines/spline_evaluator_3d.hpp include/ddc/kernels/splines/spline_traits.hpp include/ddc/kernels/splines/splines_linear_problem.hpp include/ddc/kernels/splines/splines_linear_problem_2x2_blocks.hpp From e4cf534d8920b2399567408d37ae84d5fa4e9b65 Mon Sep 17 00:00:00 2001 From: Thomas Padioleau Date: Wed, 1 Oct 2025 19:58:16 +0200 Subject: [PATCH 8/9] Increase timeout --- .github/workflows/tests-ubuntu.yaml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/tests-ubuntu.yaml b/.github/workflows/tests-ubuntu.yaml index ad3b72fd4..1506bba92 100644 --- a/.github/workflows/tests-ubuntu.yaml +++ b/.github/workflows/tests-ubuntu.yaml @@ -345,7 +345,7 @@ jobs: cmake --build build if [ 'xcpu' = 'x${{matrix.backend.name}}' ] then - ctest --test-dir build --output-on-failure --timeout 10 --output-junit tests.xml + ctest --test-dir build --output-on-failure --timeout 20 --output-junit tests.xml cp build/tests.xml . ./build/examples/characteristics_advection ./build/examples/game_of_life @@ -454,7 +454,7 @@ jobs: . /src/sanitizer_env.sh fi - ctest --test-dir build --output-on-failure --timeout 10 --output-junit tests.xml + ctest --test-dir build --output-on-failure --timeout 20 --output-junit tests.xml cp build/tests.xml . ./build/examples/characteristics_advection ./build/examples/game_of_life From 2a0154110b6ceafd36b9fbf5f3df68990db8d6ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tr=C3=A9vis=20Morvany?= Date: Thu, 2 Oct 2025 10:08:12 +0200 Subject: [PATCH 9/9] Modify 3D error bound computation --- tests/splines/evaluator_3d.hpp | 2 +- tests/splines/spline_error_bounds.hpp | 20 ++++++++++---------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/tests/splines/evaluator_3d.hpp b/tests/splines/evaluator_3d.hpp index 80c16db05..5fdac7213 100644 --- a/tests/splines/evaluator_3d.hpp +++ b/tests/splines/evaluator_3d.hpp @@ -94,7 +94,7 @@ struct Evaluator3D for (ddc::DiscreteElement const k : ddc::DiscreteDomain(domain)) { chunk(i, j, k) = eval_func1.deriv(ddc::coordinate(i), derivative_x) * eval_func2.deriv(ddc::coordinate(j), derivative_y) - * eval_func3.deriv(ddc::coordinate(j), derivative_z); + * eval_func3.deriv(ddc::coordinate(k), derivative_z); } } } diff --git a/tests/splines/spline_error_bounds.hpp b/tests/splines/spline_error_bounds.hpp index 0ee41f2ce..b2eff23ab 100644 --- a/tests/splines/spline_error_bounds.hpp +++ b/tests/splines/spline_error_bounds.hpp @@ -173,9 +173,9 @@ class SplineErrorBounds int degree2, int degree3) const { - double const norm1 = m_evaluator.max_norm(degree1 + 1, 1, 1); - double const norm2 = m_evaluator.max_norm(1, degree2 + 1, 1); - double const norm3 = m_evaluator.max_norm(1, 1, degree3 + 1); + double const norm1 = m_evaluator.max_norm(degree1 + 1, 1, 0); + double const norm2 = m_evaluator.max_norm(1, degree2 + 1, 0); + double const norm3 = m_evaluator.max_norm(0, 0, degree3 + 1); return tihomirov_error_bound(cell_width1, degree1 - 1, norm1) + tihomirov_error_bound(cell_width2, degree2 - 1, norm2) + tihomirov_error_bound(cell_width3, degree3, norm3); @@ -189,12 +189,12 @@ class SplineErrorBounds int degree2, int degree3) const { - double const norm1 = m_evaluator.max_norm(degree1 + 1, 1, 1); - double const norm2 = m_evaluator.max_norm(1, degree2 + 1, 1); - double const norm3 = m_evaluator.max_norm(1, 1, degree3 + 1); + double const norm1 = m_evaluator.max_norm(degree1 + 1, 0, 0); + double const norm2 = m_evaluator.max_norm(0, degree2 + 1, 1); + double const norm3 = m_evaluator.max_norm(0, 1, degree3 + 1); return tihomirov_error_bound(cell_width1, degree1, norm1) + tihomirov_error_bound(cell_width2, degree2 - 1, norm2) - + tihomirov_error_bound(cell_width3, degree3 - 3, norm3); + + tihomirov_error_bound(cell_width3, degree3 - 1, norm3); } double error_bound_on_deriv_13( @@ -205,9 +205,9 @@ class SplineErrorBounds int degree2, int degree3) const { - double const norm1 = m_evaluator.max_norm(degree1 + 1, 1, 1); - double const norm2 = m_evaluator.max_norm(1, degree2 + 1, 1); - double const norm3 = m_evaluator.max_norm(1, 1, degree3 + 1); + double const norm1 = m_evaluator.max_norm(degree1 + 1, 0, 1); + double const norm2 = m_evaluator.max_norm(0, degree2 + 1, 0); + double const norm3 = m_evaluator.max_norm(1, 0, degree3 + 1); return tihomirov_error_bound(cell_width1, degree1 - 1, norm1) + tihomirov_error_bound(cell_width2, degree2, norm2) + tihomirov_error_bound(cell_width3, degree3 - 1, norm3);