diff --git a/.github/workflows/tests-ubuntu.yaml b/.github/workflows/tests-ubuntu.yaml index 685f0e23d..43fc734ab 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 diff --git a/CMakeLists.txt b/CMakeLists.txt index 345e89e5d..9583340f9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -282,8 +282,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 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..10f39392b --- /dev/null +++ b/include/ddc/kernels/splines/spline_builder_3d.hpp @@ -0,0 +1,708 @@ +// 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" + +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; + + /// @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; + + /// @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; + + /** + * @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_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 SplineBuilder3D 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_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 SplineBuilder3D 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 SplineBuilder3D. + */ + 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 3D interpolation mesh used by this class. + * + * This is 3D because it is defined along the dimensions of interest. + * + * @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_domain_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 3D 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 3D domain on which spline coefficients are defined. + * + * The 3D spline domain corresponding to the dimensions of interest. + * + * @return The 3D 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 3D spline approximation of a function. + * + * Use the values of a function (defined on + * SplineBuilder3D::batched_interpolation_domain) and the derivatives of the + * function at the boundaries (in the case of BoundCond::HERMITE only) + * to calculate a 3D 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] 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()( + ddc::ChunkSpan< + double, + batched_spline_domain_type, + Layout, + memory_space> spline, + ddc::ChunkSpan vals, + std::optional, + Layout, + memory_space>> derivs_min1 + = std::nullopt, + std::optional, + Layout, + memory_space>> derivs_max1 + = std::nullopt, + std::optional, + Layout, + memory_space>> derivs_min2 + = std::nullopt, + std::optional, + Layout, + memory_space>> derivs_max2 + = std::nullopt, + std::optional, + Layout, + memory_space>> derivs_min3 + = std::nullopt, + std::optional, + Layout, + memory_space>> derivs_max3 + = 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, + [[maybe_unused]] std::optional, + Layout, + memory_space>> derivs_min1, + [[maybe_unused]] std::optional, + Layout, + memory_space>> derivs_max1, + [[maybe_unused]] std::optional, + Layout, + memory_space>> derivs_min2, + [[maybe_unused]] std::optional, + Layout, + memory_space>> derivs_max2, + [[maybe_unused]] std::optional, + Layout, + memory_space>> derivs_min3, + [[maybe_unused]] std::optional, + Layout, + memory_space>> derivs_max3, + [[maybe_unused]] std::optional, + Layout, + memory_space>> mixed_derivs_min1_min2_min3, + [[maybe_unused]] std::optional, + Layout, + memory_space>> mixed_derivs_max1_min2_min3, + [[maybe_unused]] std::optional, + Layout, + memory_space>> mixed_derivs_min1_max2_min3, + [[maybe_unused]] std::optional, + Layout, + memory_space>> mixed_derivs_max1_max2_min3, + [[maybe_unused]] std::optional, + Layout, + memory_space>> mixed_derivs_min1_min2_max3, + [[maybe_unused]] std::optional, + Layout, + memory_space>> mixed_derivs_max1_min2_max3, + [[maybe_unused]] std::optional, + Layout, + memory_space>> mixed_derivs_min1_max2_max3, + [[maybe_unused]] 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..83006ee0f --- /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_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 + */ + 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..ca68b1a96 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..d84fa4677 --- /dev/null +++ b/tests/splines/batched_3d_spline_builder.cpp @@ -0,0 +1,849 @@ +// 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) + // 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, 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_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 + 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: 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 + 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: This is the 2D version, modify when adding the Hermite boundary conditions + 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), + 5e-10 * 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..5fdac7213 --- /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(k), 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..b2eff23ab 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, 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); + } + + 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, 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 - 1, 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, 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); + } + + 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));