diff --git a/AUTHORS b/AUTHORS index 704743132..a92fb5f5c 100644 --- a/AUTHORS +++ b/AUTHORS @@ -35,3 +35,6 @@ Baptiste Legouix - CEA (baptiste.legouix@cea.fr) Thierry Antoun - CEA (thierry.antoun@cea.fr) * Work on Documentation * Work on adding new features + +Trévis Morvany - CEA (trevis.morvany@cea.fr) +* Work on splines diff --git a/benchmarks/splines.cpp b/benchmarks/splines.cpp index 9a17016f5..8d3037588 100644 --- a/benchmarks/splines.cpp +++ b/benchmarks/splines.cpp @@ -160,9 +160,7 @@ void characteristics_advection_unitary(benchmark::State& state) DDimX, ddc::BoundCond::PERIODIC, ddc::BoundCond::PERIODIC, - Backend, - DDimX, - DDimY> const spline_builder(x_mesh, cols_per_chunk, preconditioner_max_block_size); + Backend> const spline_builder(x_domain, cols_per_chunk, preconditioner_max_block_size); ddc::PeriodicExtrapolationRule const periodic_extrapolation; ddc::SplineEvaluator< ExecSpace, @@ -170,15 +168,14 @@ void characteristics_advection_unitary(benchmark::State& state) BSplinesX, DDimX, ddc::PeriodicExtrapolationRule, - ddc::PeriodicExtrapolationRule, - DDimX, - DDimY> const spline_evaluator(periodic_extrapolation, periodic_extrapolation); + ddc::PeriodicExtrapolationRule> const + spline_evaluator(periodic_extrapolation, periodic_extrapolation); ddc::Chunk coef_alloc( - spline_builder.batched_spline_domain(), + spline_builder.batched_spline_domain(x_mesh), ddc::KokkosAllocator()); ddc::ChunkSpan const coef = coef_alloc.span_view(); ddc::Chunk feet_coords_alloc( - spline_builder.batched_interpolation_domain(), + spline_builder.batched_interpolation_domain(x_mesh), ddc::KokkosAllocator, typename ExecSpace::memory_space>()); ddc::ChunkSpan const feet_coords = feet_coords_alloc.span_view(); diff --git a/examples/characteristics_advection.cpp b/examples/characteristics_advection.cpp index 3ebb25880..e97c86604 100644 --- a/examples/characteristics_advection.cpp +++ b/examples/characteristics_advection.cpp @@ -207,9 +207,7 @@ int main(int argc, char** argv) DDimX, BoundCond, BoundCond, - ddc::SplineSolver::LAPACK, - DDimX, - DDimY> const spline_builder(x_mesh); + ddc::SplineSolver::LAPACK> const spline_builder(x_domain); ExtrapolationRule const extrapolation_rule; ddc::SplineEvaluator< Kokkos::DefaultExecutionSpace, @@ -217,19 +215,19 @@ int main(int argc, char** argv) BSplinesX, DDimX, ExtrapolationRule, - ExtrapolationRule, - DDimX, - DDimY> const spline_evaluator(extrapolation_rule, extrapolation_rule); + ExtrapolationRule> const spline_evaluator(extrapolation_rule, extrapolation_rule); //! [instantiate solver] //! [instantiate intermediate chunks] // Instantiate chunk of spline coefs to receive output of spline_builder - ddc::Chunk coef_alloc(spline_builder.batched_spline_domain(), ddc::DeviceAllocator()); + ddc::Chunk coef_alloc( + spline_builder.batched_spline_domain(x_mesh), + ddc::DeviceAllocator()); ddc::ChunkSpan const coef = coef_alloc.span_view(); // Instantiate chunk to receive feet coords ddc::Chunk feet_coords_alloc( - spline_builder.batched_interpolation_domain(), + spline_builder.batched_interpolation_domain(x_mesh), ddc::DeviceAllocator>()); ddc::ChunkSpan const feet_coords = feet_coords_alloc.span_view(); //! [instantiate intermediate chunks] diff --git a/include/ddc/kernels/splines/spline_builder.hpp b/include/ddc/kernels/splines/spline_builder.hpp index 1cc9e42e0..66c9e91e6 100644 --- a/include/ddc/kernels/splines/spline_builder.hpp +++ b/include/ddc/kernels/splines/spline_builder.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include @@ -49,7 +50,6 @@ enum class SplineSolver { * @tparam BcLower The lower boundary condition. * @tparam BcUpper The upper boundary condition. * @tparam Solver The SplineSolver giving the backend used to perform the spline approximation. - * @tparam DDimX A variadic template of all the discrete dimensions forming the full space (InterpolationDDim + batched dimensions). */ template < class ExecSpace, @@ -58,8 +58,7 @@ template < class InterpolationDDim, ddc::BoundCond BcLower, ddc::BoundCond BcUpper, - SplineSolver Solver, - class... DDimX> + SplineSolver Solver> class SplineBuilder { static_assert( @@ -90,29 +89,45 @@ class SplineBuilder /// @brief The type of the domain for the 1D interpolation mesh used by this class. using interpolation_domain_type = ddc::DiscreteDomain; - /// @brief The type of the whole domain representing interpolation points. - using batched_interpolation_domain_type = ddc::DiscreteDomain; + /** + * @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 dimension 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 a dimension of interest Y, * this is DiscreteDomain */ - using batch_domain_type = ddc::remove_dims_of_t< - batched_interpolation_domain_type, - interpolation_discrete_dimension_type>; + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using batch_domain_type = ddc:: + remove_dims_of_t; /** * @brief The type of the whole spline domain (cartesian product of 1D 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. + * * Example: For batched_interpolation_domain_type = DiscreteDomain and a dimension of interest Y * (associated to a B-splines tag BSplinesY), this is DiscreteDomain. */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> using batched_spline_domain_type = ddc::replace_dim_of_t< - batched_interpolation_domain_type, + BatchedInterpolationDDom, interpolation_discrete_dimension_type, bsplines_type>; @@ -121,14 +136,19 @@ class SplineBuilder * @brief The type of the whole spline domain (cartesian product of the 1D spline domain * and the batch domain) with 1D spline dimension being the leading dimension. * + * @tparam The batched discrete domain on which the interpolation points are defined. + * * Example: For batched_interpolation_domain_type = DiscreteDomain and a dimension of interest Y * (associated to a B-splines tag BSplinesY), this is DiscreteDomain. */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> using batched_spline_tr_domain_type = typename ddc::detail::convert_type_seq_to_discrete_domain_t, ddc::type_seq_remove_t< - ddc::detail::TypeSeq, + ddc::to_type_seq_t, ddc::detail::TypeSeq>>>; public: @@ -136,11 +156,16 @@ class SplineBuilder * @brief The type of the whole Deriv domain (cartesian product of 1D Deriv domain * and batch domain) preserving the underlying memory layout (order of dimensions). * + * @tparam The batched discrete domain on which the interpolation points are defined. + * * Example: For batched_interpolation_domain_type = DiscreteDomain and a dimension of interest Y, * this is DiscreteDomain,Z> */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> using batched_derivs_domain_type = ddc::replace_dim_of_t< - batched_interpolation_domain_type, + BatchedInterpolationDDom, interpolation_discrete_dimension_type, deriv_type>; @@ -163,7 +188,7 @@ class SplineBuilder static constexpr SplineSolver s_spline_solver = Solver; private: - batched_interpolation_domain_type m_batched_interpolation_domain; + interpolation_domain_type m_interpolation_domain; int m_offset = 0; @@ -177,9 +202,9 @@ class SplineBuilder public: /** - * @brief Build a SplineBuilder acting on batched_interpolation_domain. + * @brief Build a SplineBuilder acting on interpolation_domain. * - * @param batched_interpolation_domain The domain on which the interpolation points are defined. + * @param 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 @@ -193,10 +218,10 @@ class SplineBuilder * @see MatrixSparse */ explicit SplineBuilder( - batched_interpolation_domain_type const& batched_interpolation_domain, + interpolation_domain_type const& interpolation_domain, std::optional cols_per_chunk = std::nullopt, std::optional preconditioner_max_block_size = std::nullopt) - : m_batched_interpolation_domain(batched_interpolation_domain) + : m_interpolation_domain(interpolation_domain) , m_dx((ddc::discrete_space().rmax() - ddc::discrete_space().rmin()) / ddc::discrete_space().ncells()) { @@ -205,7 +230,7 @@ class SplineBuilder "Incompatible boundary conditions"); check_valid_grid(); - compute_offset(interpolation_domain(), m_offset); + compute_offset(this->interpolation_domain(), m_offset); // Calculate block sizes int lower_block_size; @@ -224,6 +249,37 @@ class SplineBuilder preconditioner_max_block_size); } + /** + * @brief Build a SplineBuilder acting on the interpolation domain contained by batched_interpolation_domain. + * + * @param batched_interpolation_domain The whole 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 MatrixSparse + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + explicit SplineBuilder( + BatchedInterpolationDDom const& batched_interpolation_domain, + std::optional cols_per_chunk = std::nullopt, + std::optional preconditioner_max_block_size = std::nullopt) + : SplineBuilder( + interpolation_domain_type(batched_interpolation_domain), + cols_per_chunk, + preconditioner_max_block_size) + { + } + + /// @brief Copy-constructor is deleted. SplineBuilder(SplineBuilder const& x) = delete; @@ -255,7 +311,7 @@ class SplineBuilder */ interpolation_domain_type interpolation_domain() const noexcept { - return interpolation_domain_type(m_batched_interpolation_domain); + return m_interpolation_domain; } /** @@ -264,11 +320,18 @@ class SplineBuilder * Values of the function must be provided on this domain in order * to build a spline representation of the function (cartesian product of 1D 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. */ - batched_interpolation_domain_type batched_interpolation_domain() const noexcept + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + BatchedInterpolationDDom batched_interpolation_domain( + BatchedInterpolationDDom const& batched_interpolation_domain) const noexcept { - return m_batched_interpolation_domain; + assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain)); + return batched_interpolation_domain; } /** @@ -276,11 +339,16 @@ class SplineBuilder * * Obtained by removing the dimension 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. */ - batch_domain_type batch_domain() const noexcept + template + batch_domain_type batch_domain( + BatchedInterpolationDDom const& batched_interpolation_domain) const noexcept { - return ddc::remove_dims_of(batched_interpolation_domain(), interpolation_domain()); + assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain)); + return ddc::remove_dims_of(batched_interpolation_domain, interpolation_domain()); } /** @@ -300,13 +368,18 @@ class SplineBuilder * * 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. */ - batched_spline_domain_type batched_spline_domain() const noexcept + template + 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< interpolation_discrete_dimension_type, - bsplines_type>(batched_interpolation_domain(), spline_domain()); + bsplines_type>(batched_interpolation_domain, spline_domain()); } private: @@ -315,16 +388,22 @@ class SplineBuilder * * This is used internally due to solver limitation and because it may be beneficial to computation performance. For LAPACK backend and non-periodic boundary condition, we are using SplinesLinearSolver3x3Blocks which requires upper_block_size additional rows for internal operations. * + * @param batched_interpolation_domain The whole domain on which the interpolation points are defined. + * * @return The (transposed) domain for the spline coefficients. */ - batched_spline_tr_domain_type batched_spline_tr_domain() const noexcept + template + batched_spline_tr_domain_type batched_spline_tr_domain( + BatchedInterpolationDDom const& batched_interpolation_domain) const noexcept { - return batched_spline_tr_domain_type(ddc::replace_dim_of( - batched_spline_domain(), - ddc::DiscreteDomain( - ddc::DiscreteElement(0), - ddc::DiscreteVector( - matrix->required_number_of_rhs_rows())))); + assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain)); + return batched_spline_tr_domain_type( + ddc::replace_dim_of( + batched_spline_domain(batched_interpolation_domain), + ddc::DiscreteDomain( + ddc::DiscreteElement(0), + ddc::DiscreteVector( + matrix->required_number_of_rhs_rows())))); } public: @@ -333,12 +412,17 @@ class SplineBuilder * * This is only used with BoundCond::HERMITE boundary conditions. * + * @param batched_interpolation_domain The whole domain on which the interpolation points are defined. + * * @return The domain for the Derivs values. */ - batched_derivs_domain_type batched_derivs_xmin_domain() const noexcept + template + batched_derivs_domain_type batched_derivs_xmin_domain( + BatchedInterpolationDDom const& batched_interpolation_domain) const noexcept { + assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain)); return ddc::replace_dim_of( - batched_interpolation_domain(), + batched_interpolation_domain, ddc::DiscreteDomain( ddc::DiscreteElement(1), ddc::DiscreteVector(s_nbc_xmin))); @@ -349,12 +433,17 @@ class SplineBuilder * * This is only used with BoundCond::HERMITE boundary conditions. * + * @param batched_interpolation_domain The whole domain on which the interpolation points are defined. + * * @return The domain for the Derivs values. */ - batched_derivs_domain_type batched_derivs_xmax_domain() const noexcept + template + batched_derivs_domain_type batched_derivs_xmax_domain( + BatchedInterpolationDDom const& batched_interpolation_domain) const noexcept { + assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain)); return ddc::replace_dim_of( - batched_interpolation_domain(), + batched_interpolation_domain, ddc::DiscreteDomain( ddc::DiscreteElement(1), ddc::DiscreteVector(s_nbc_xmax))); @@ -379,18 +468,25 @@ class SplineBuilder * @param[in] derivs_xmax The values of the derivatives at the upper boundary * (used only with BoundCond::HERMITE upper boundary condition). */ - template + template void operator()( - ddc::ChunkSpan spline, - ddc::ChunkSpan - vals, - std::optional< - ddc::ChunkSpan> - derivs_xmin + ddc::ChunkSpan< + double, + batched_spline_domain_type, + Layout, + memory_space> spline, + ddc::ChunkSpan vals, + std::optional, + Layout, + memory_space>> derivs_xmin = std::nullopt, - std::optional< - ddc::ChunkSpan> - derivs_xmax + std::optional, + Layout, + memory_space>> derivs_xmax = std::nullopt) const; /** @@ -460,17 +556,8 @@ template < class InterpolationDDim, ddc::BoundCond BcLower, ddc::BoundCond BcUpper, - SplineSolver Solver, - class... DDimX> -void SplineBuilder< - ExecSpace, - MemorySpace, - BSplines, - InterpolationDDim, - BcLower, - BcUpper, - Solver, - DDimX...>:: + SplineSolver Solver> +void SplineBuilder:: compute_offset(interpolation_domain_type const& interpolation_domain, int& offset) { if constexpr (bsplines_type::is_periodic()) { @@ -504,17 +591,9 @@ template < class InterpolationDDim, ddc::BoundCond BcLower, ddc::BoundCond BcUpper, - SplineSolver Solver, - class... DDimX> -int SplineBuilder< - ExecSpace, - MemorySpace, - BSplines, - InterpolationDDim, - BcLower, - BcUpper, - Solver, - DDimX...>::compute_block_sizes_uniform(ddc::BoundCond const bound_cond, int const nbc) + SplineSolver Solver> +int SplineBuilder:: + compute_block_sizes_uniform(ddc::BoundCond const bound_cond, int const nbc) { if (bound_cond == ddc::BoundCond::PERIODIC) { return static_cast(bsplines_type::degree()) / 2; @@ -538,17 +617,9 @@ template < class InterpolationDDim, ddc::BoundCond BcLower, ddc::BoundCond BcUpper, - SplineSolver Solver, - class... DDimX> -int SplineBuilder< - ExecSpace, - MemorySpace, - BSplines, - InterpolationDDim, - BcLower, - BcUpper, - Solver, - DDimX...>::compute_block_sizes_non_uniform(ddc::BoundCond const bound_cond, int const nbc) + SplineSolver Solver> +int SplineBuilder:: + compute_block_sizes_non_uniform(ddc::BoundCond const bound_cond, int const nbc) { if (bound_cond == ddc::BoundCond::PERIODIC || bound_cond == ddc::BoundCond::GREVILLE) { return static_cast(bsplines_type::degree()) - 1; @@ -568,17 +639,8 @@ template < class InterpolationDDim, ddc::BoundCond BcLower, ddc::BoundCond BcUpper, - SplineSolver Solver, - class... DDimX> -void SplineBuilder< - ExecSpace, - MemorySpace, - BSplines, - InterpolationDDim, - BcLower, - BcUpper, - Solver, - DDimX...>:: + SplineSolver Solver> +void SplineBuilder:: allocate_matrix( [[maybe_unused]] int lower_block_size, [[maybe_unused]] int upper_block_size, @@ -635,17 +697,9 @@ template < class InterpolationDDim, ddc::BoundCond BcLower, ddc::BoundCond BcUpper, - SplineSolver Solver, - class... DDimX> -void SplineBuilder< - ExecSpace, - MemorySpace, - BSplines, - InterpolationDDim, - BcLower, - BcUpper, - Solver, - DDimX...>::build_matrix_system() + SplineSolver Solver> +void SplineBuilder:: + build_matrix_system() { // Hermite boundary conditions at xmin, if any if constexpr (BcLower == ddc::BoundCond::HERMITE) { @@ -738,32 +792,31 @@ template < class InterpolationDDim, ddc::BoundCond BcLower, ddc::BoundCond BcUpper, - SplineSolver Solver, - class... DDimX> -template -void SplineBuilder< - ExecSpace, - MemorySpace, - BSplines, - InterpolationDDim, - BcLower, - BcUpper, - Solver, - DDimX...>:: + SplineSolver Solver> +template +void SplineBuilder:: operator()( - ddc::ChunkSpan spline, - ddc::ChunkSpan vals, + ddc::ChunkSpan< + double, + batched_spline_domain_type, + Layout, + memory_space> spline, + ddc::ChunkSpan vals, std::optional, Layout, memory_space>> const derivs_xmin, std::optional, Layout, memory_space>> const derivs_xmax) const { + auto const batched_interpolation_domain = vals.domain(); + + assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain)); + assert(vals.template extent() == ddc::discrete_space().nbasis() - s_nbc_xmin - s_nbc_xmax); @@ -788,8 +841,10 @@ operator()( ddc::parallel_for_each( "ddc_splines_hermite_compute_lower_coefficients", exec_space(), - batch_domain(), - KOKKOS_LAMBDA(typename batch_domain_type::discrete_element_type j) { + batch_domain(batched_interpolation_domain), + KOKKOS_LAMBDA( + typename batch_domain_type::discrete_element_type + j) { for (int i = s_nbc_xmin; i > 0; --i) { spline(ddc::DiscreteElement(s_nbc_xmin - i), j) = derivs_xmin_values(ddc::DiscreteElement(i), j) @@ -831,8 +886,10 @@ operator()( ddc::parallel_for_each( "ddc_splines_hermite_compute_upper_coefficients", exec_space(), - batch_domain(), - KOKKOS_LAMBDA(typename batch_domain_type::discrete_element_type j) { + batch_domain(batched_interpolation_domain), + KOKKOS_LAMBDA( + typename batch_domain_type::discrete_element_type + j) { for (int i = 0; i < s_nbc_xmax; ++i) { spline(ddc::DiscreteElement(nbasis_proxy - s_nbc_xmax - i), j) @@ -845,14 +902,15 @@ operator()( // Allocate and fill a transposed version of spline in order to get dimension of interest as last dimension (optimal for GPU, necessary for Ginkgo). Also select only relevant rows in case of periodic boundaries auto const& offset_proxy = m_offset; ddc::Chunk spline_tr_alloc( - batched_spline_tr_domain(), + batched_spline_tr_domain(batched_interpolation_domain), ddc::KokkosAllocator()); ddc::ChunkSpan const spline_tr = spline_tr_alloc.span_view(); ddc::parallel_for_each( "ddc_splines_transpose_rhs", exec_space(), - batch_domain(), - KOKKOS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { + batch_domain(batched_interpolation_domain), + KOKKOS_LAMBDA(typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { for (std::size_t i = 0; i < nbasis_proxy; ++i) { spline_tr(ddc::DiscreteElement(i), j) = spline(ddc::DiscreteElement(i + offset_proxy), j); @@ -862,15 +920,16 @@ operator()( Kokkos::View const bcoef_section( spline_tr.data_handle(), static_cast(spline_tr.template extent()), - batch_domain().size()); + batch_domain(batched_interpolation_domain).size()); // Compute spline coef matrix->solve(bcoef_section, false); // Transpose back spline_tr into spline. ddc::parallel_for_each( "ddc_splines_transpose_back_rhs", exec_space(), - batch_domain(), - KOKKOS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { + batch_domain(batched_interpolation_domain), + KOKKOS_LAMBDA(typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { for (std::size_t i = 0; i < nbasis_proxy; ++i) { spline(ddc::DiscreteElement(i + offset_proxy), j) = spline_tr(ddc::DiscreteElement(i), j); @@ -882,8 +941,9 @@ operator()( ddc::parallel_for_each( "ddc_splines_periodic_rows_duplicate_rhs", exec_space(), - batch_domain(), - KOKKOS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { + batch_domain(batched_interpolation_domain), + KOKKOS_LAMBDA(typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { if (offset_proxy != 0) { for (int i = 0; i < offset_proxy; ++i) { spline(ddc::DiscreteElement(i), j) = spline( @@ -912,8 +972,7 @@ template < class InterpolationDDim, ddc::BoundCond BcLower, ddc::BoundCond BcUpper, - SplineSolver Solver, - class... DDimX> + SplineSolver Solver> template std::tuple< ddc::Chunk< @@ -930,15 +989,8 @@ std::tuple< ddc::DiscreteDomain< ddc::Deriv>, ddc::KokkosAllocator>> -SplineBuilder< - ExecSpace, - MemorySpace, - BSplines, - InterpolationDDim, - BcLower, - BcUpper, - Solver, - DDimX...>::quadrature_coefficients() const +SplineBuilder:: + quadrature_coefficients() const { // Compute integrals of bsplines ddc::Chunk integral_bsplines(spline_domain(), ddc::KokkosAllocator()); @@ -1055,18 +1107,9 @@ template < class InterpolationDDim, ddc::BoundCond BcLower, ddc::BoundCond BcUpper, - SplineSolver Solver, - class... DDimX> + SplineSolver Solver> template -void SplineBuilder< - ExecSpace, - MemorySpace, - BSplines, - InterpolationDDim, - BcLower, - BcUpper, - Solver, - DDimX...>:: +void SplineBuilder:: check_n_points_in_cell(int const n_points_in_cell, KnotElement const current_cell_end_idx) { if (n_points_in_cell > BSplines::degree() + 1) { @@ -1086,17 +1129,9 @@ template < class InterpolationDDim, ddc::BoundCond BcLower, ddc::BoundCond BcUpper, - SplineSolver Solver, - class... DDimX> -void SplineBuilder< - ExecSpace, - MemorySpace, - BSplines, - InterpolationDDim, - BcLower, - BcUpper, - Solver, - DDimX...>::check_valid_grid() + SplineSolver Solver> +void SplineBuilder:: + check_valid_grid() { std::size_t const n_interp_points = interpolation_domain().size(); std::size_t const expected_npoints diff --git a/include/ddc/kernels/splines/spline_builder_2d.hpp b/include/ddc/kernels/splines/spline_builder_2d.hpp index 1a2f4edc6..03dc45e10 100644 --- a/include/ddc/kernels/splines/spline_builder_2d.hpp +++ b/include/ddc/kernels/splines/spline_builder_2d.hpp @@ -4,6 +4,7 @@ #pragma once +#include #include #include #include @@ -33,8 +34,7 @@ template < ddc::BoundCond BcUpper1, ddc::BoundCond BcLower2, ddc::BoundCond BcUpper2, - ddc::SplineSolver Solver, - class... DDimX> + ddc::SplineSolver Solver> class SplineBuilder2D { public: @@ -45,40 +45,16 @@ class SplineBuilder2D 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< - ExecSpace, - MemorySpace, - BSpline1, - DDimI1, - BcLower1, - BcUpper1, - Solver, - DDimX...>; + 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< - ExecSpace, - MemorySpace, - BSpline2, - DDimI2, - BcLower2, - BcUpper2, - Solver, - std::conditional_t, BSpline1, DDimX>...>; + using builder_type2 = 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< - ExecSpace, - MemorySpace, - BSpline1, - DDimI1, - BcLower1, - BcUpper1, - Solver, - std::conditional_t< - std::is_same_v, - typename builder_type2::deriv_type, - DDimX>...>; + using builder_deriv_type1 = ddc:: + SplineBuilder; /// @brief The type of the first interpolation continuous dimension. using continuous_dimension_type1 = typename builder_type1::continuous_dimension_type; @@ -119,18 +95,30 @@ class SplineBuilder2D interpolation_discrete_dimension_type1, interpolation_discrete_dimension_type2>; - /// @brief The type of the whole domain representing interpolation points. - using batched_interpolation_domain_type = ddc::DiscreteDomain; + /** + * @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 and Y, * this is DiscreteDomain. */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> using batch_domain_type = ddc::remove_dims_of_t< - batched_interpolation_domain_type, + BatchedInterpolationDDom, interpolation_discrete_dimension_type1, interpolation_discrete_dimension_type2>; @@ -138,12 +126,17 @@ class SplineBuilder2D * @brief The type of the whole spline domain (cartesian product of 2D 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 and Y * (associated to B-splines tags BSplinesX and BSplinesY), 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::to_type_seq_t, ddc::detail::TypeSeq< interpolation_discrete_dimension_type1, interpolation_discrete_dimension_type2>, @@ -153,20 +146,31 @@ class SplineBuilder2D * @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 and Y, * this is DiscreteDomain, Y, Z>. */ - using batched_derivs_domain_type1 = typename builder_type1::batched_derivs_domain_type; + 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 and Y, * this is DiscreteDomain, Z>. */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> using batched_derivs_domain_type2 = ddc::replace_dim_of_t< - batched_interpolation_domain_type, + BatchedInterpolationDDom, interpolation_discrete_dimension_type2, deriv_type2>; @@ -174,12 +178,17 @@ class SplineBuilder2D * @brief The type of the whole Derivs domain (cartesian product of the 2D Deriv domain * and the 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 and Y, * this is DiscreteDomain, Deriv, Z>. */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> using batched_derivs_domain_type = ddc::detail::convert_type_seq_to_discrete_domain_t, + ddc::to_type_seq_t, ddc::detail::TypeSeq< interpolation_discrete_dimension_type1, interpolation_discrete_dimension_type2>, @@ -192,7 +201,33 @@ class SplineBuilder2D public: /** - * @brief Build a SplineBuilder2D acting on batched_interpolation_domain. + * @brief Build a SplineBuilder2D acting on interpolation_domain. + * + * @param interpolation_domain The domain on which the interpolation points are defined, without the batch dimensions. + * + * @param cols_per_chunk A parameter used by the slicer (internal to the solver) to define the size + * of a chunk of right-hand-sides of the linear problem to be computed in parallel (chunks are treated + * by the linear solver one-after-the-other). + * This value is optional. If no value is provided then the default value is chosen by the requested solver. + * + * @param preconditioner_max_block_size A parameter used by the slicer (internal to the solver) to + * define the size of a block used by the Block-Jacobi preconditioner. + * This value is optional. If no value is provided then the default value is chosen by the requested solver. + * + * @see SplinesLinearProblemSparse + */ + explicit SplineBuilder2D( + interpolation_domain_type const& interpolation_domain, + std::optional cols_per_chunk = std::nullopt, + std::optional preconditioner_max_block_size = std::nullopt) + : m_spline_builder1(interpolation_domain, cols_per_chunk, preconditioner_max_block_size) + , m_spline_builder_deriv1(interpolation_domain) + , m_spline_builder2(interpolation_domain, cols_per_chunk, preconditioner_max_block_size) + { + } + + /** + * @brief Build a SplineBuilder2D acting on the interpolation domain contained in batched_interpolation_domain. * * @param batched_interpolation_domain The domain on which the interpolation points are defined. * @@ -207,22 +242,15 @@ class SplineBuilder2D * * @see SplinesLinearProblemSparse */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> explicit SplineBuilder2D( - batched_interpolation_domain_type const& batched_interpolation_domain, + BatchedInterpolationDDom const& batched_interpolation_domain, std::optional cols_per_chunk = std::nullopt, std::optional preconditioner_max_block_size = std::nullopt) - : m_spline_builder1( - batched_interpolation_domain, - cols_per_chunk, - preconditioner_max_block_size) - , m_spline_builder_deriv1( - ddc::replace_dim_of( - m_spline_builder1.batched_interpolation_domain(), - ddc::DiscreteDomain( - ddc::DiscreteElement(1), - ddc::DiscreteVector(bsplines_type2::degree() / 2)))) - , m_spline_builder2( - m_spline_builder1.batched_spline_domain(), + : SplineBuilder2D( + interpolation_domain_type(batched_interpolation_domain), cols_per_chunk, preconditioner_max_block_size) { @@ -271,11 +299,18 @@ class SplineBuilder2D * Values of the function must be provided on this domain in order * to build a spline representation of the function (cartesian product of 2D interpolation_domain and batch_domain). * + * @param batched_interpolation_domain The whole domain on which the interpolation points are defined. + * * @return The domain for the interpolation mesh. */ - batched_interpolation_domain_type batched_interpolation_domain() const noexcept + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + BatchedInterpolationDDom batched_interpolation_domain( + BatchedInterpolationDDom const& batched_interpolation_domain) const noexcept { - return m_spline_builder1.batched_interpolation_domain(); + assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain)); + return batched_interpolation_domain; } /** @@ -283,11 +318,18 @@ class SplineBuilder2D * * 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. */ - batch_domain_type batch_domain() const noexcept + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + batch_domain_type batch_domain( + BatchedInterpolationDDom const& batched_interpolation_domain) const noexcept { - return ddc::remove_dims_of(batched_interpolation_domain(), interpolation_domain()); + assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain)); + return ddc::remove_dims_of(batched_interpolation_domain, interpolation_domain()); } /** @@ -309,14 +351,21 @@ class SplineBuilder2D * * 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. */ - batched_spline_domain_type batched_spline_domain() const noexcept + 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< interpolation_discrete_dimension_type2, - bsplines_type2>(batched_interpolation_domain(), spline_domain()), + bsplines_type2>(batched_interpolation_domain, spline_domain()), spline_domain()); } @@ -356,42 +405,61 @@ class SplineBuilder2D * The values of the the cross-derivatives at the upper boundary in the first dimension * and the upper boundary in the second dimension. */ - template + template void operator()( - ddc::ChunkSpan spline, - ddc::ChunkSpan - vals, - std::optional< - ddc::ChunkSpan> - derivs_min1 + 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< - ddc::ChunkSpan> - derivs_max1 + std::optional, + Layout, + memory_space>> derivs_max1 = std::nullopt, - std::optional< - ddc::ChunkSpan> - derivs_min2 + std::optional, + Layout, + memory_space>> derivs_min2 = std::nullopt, - std::optional< - ddc::ChunkSpan> - derivs_max2 + std::optional, + Layout, + memory_space>> derivs_max2 = std::nullopt, - std::optional< - ddc::ChunkSpan> - mixed_derivs_min1_min2 + std::optional, + Layout, + memory_space>> mixed_derivs_min1_min2 = std::nullopt, - std::optional< - ddc::ChunkSpan> - mixed_derivs_max1_min2 + std::optional, + Layout, + memory_space>> mixed_derivs_max1_min2 = std::nullopt, - std::optional< - ddc::ChunkSpan> - mixed_derivs_min1_max2 + std::optional, + Layout, + memory_space>> mixed_derivs_min1_max2 = std::nullopt, - std::optional< - ddc::ChunkSpan> - mixed_derivs_max1_max2 + std::optional, + Layout, + memory_space>> mixed_derivs_max1_max2 = std::nullopt) const; }; @@ -407,9 +475,8 @@ template < ddc::BoundCond BcUpper1, ddc::BoundCond BcLower2, ddc::BoundCond BcUpper2, - ddc::SplineSolver Solver, - class... DDimX> -template + ddc::SplineSolver Solver> +template void SplineBuilder2D< ExecSpace, MemorySpace, @@ -421,56 +488,71 @@ void SplineBuilder2D< BcUpper1, BcLower2, BcUpper2, - Solver, - DDimX...>:: + Solver>:: operator()( - ddc::ChunkSpan spline, - ddc::ChunkSpan vals, + ddc::ChunkSpan< + double, + batched_spline_domain_type, + Layout, + memory_space> spline, + ddc::ChunkSpan vals, std::optional, Layout, memory_space>> const derivs_min1, std::optional, Layout, memory_space>> const derivs_max1, std::optional, Layout, memory_space>> const derivs_min2, std::optional, Layout, memory_space>> const derivs_max2, std::optional, Layout, memory_space>> const mixed_derivs_min1_min2, std::optional, Layout, memory_space>> const mixed_derivs_max1_min2, std::optional, Layout, memory_space>> const mixed_derivs_min1_max2, std::optional, Layout, memory_space>> const mixed_derivs_max1_max2) const { + auto const batched_interpolation_domain = vals.domain(); + + assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain)); + // TODO: perform computations along dimension 1 on different streams ? // Spline1-approximate derivs_min2 (to spline1_deriv_min) + + auto const batched_interpolation_deriv_domain + = ddc::replace_dim_of( + batched_interpolation_domain, + ddc::DiscreteDomain( + ddc::DiscreteElement(1), + ddc::DiscreteVector(bsplines_type2::degree() / 2))); + ddc::Chunk spline1_deriv_min_alloc( - m_spline_builder_deriv1.batched_spline_domain(), + m_spline_builder_deriv1.batched_spline_domain(batched_interpolation_deriv_domain), ddc::KokkosAllocator()); auto spline1_deriv_min = spline1_deriv_min_alloc.span_view(); auto spline1_deriv_min_opt = std::optional(spline1_deriv_min.span_cview()); @@ -486,7 +568,7 @@ operator()( // Spline1-approximate vals (to spline1) ddc::Chunk spline1_alloc( - m_spline_builder1.batched_spline_domain(), + m_spline_builder1.batched_spline_domain(batched_interpolation_domain), ddc::KokkosAllocator()); ddc::ChunkSpan const spline1 = spline1_alloc.span_view(); @@ -494,7 +576,7 @@ operator()( // Spline1-approximate derivs_max2 (to spline1_deriv_max) ddc::Chunk spline1_deriv_max_alloc( - m_spline_builder_deriv1.batched_spline_domain(), + m_spline_builder_deriv1.batched_spline_domain(batched_interpolation_deriv_domain), ddc::KokkosAllocator()); auto spline1_deriv_max = spline1_deriv_max_alloc.span_view(); auto spline1_deriv_max_opt = std::optional(spline1_deriv_max.span_cview()); diff --git a/include/ddc/kernels/splines/spline_evaluator.hpp b/include/ddc/kernels/splines/spline_evaluator.hpp index 2c1c63e7b..0b60426e8 100644 --- a/include/ddc/kernels/splines/spline_evaluator.hpp +++ b/include/ddc/kernels/splines/spline_evaluator.hpp @@ -29,7 +29,6 @@ namespace ddc { * @tparam EvaluationDDim The discrete dimension on which evaluation points are defined. * @tparam LowerExtrapolationRule The lower extrapolation rule type. * @tparam UpperExtrapolationRule The upper extrapolation rule type. - * @tparam DDimX A variadic template of all the discrete dimensions forming the full space (EvaluationDDim + batched dimensions). */ template < class ExecSpace, @@ -37,8 +36,7 @@ template < class BSplines, class EvaluationDDim, class LowerExtrapolationRule, - class UpperExtrapolationRule, - class... DDimX> + class UpperExtrapolationRule> class SplineEvaluator { private: @@ -75,8 +73,15 @@ class SplineEvaluator /// @brief The type of the domain for the 1D evaluation mesh used by this class. using evaluation_domain_type = ddc::DiscreteDomain; - /// @brief The type of the whole domain representing evaluation points. - using batched_evaluation_domain_type = ddc::DiscreteDomain; + /** + * @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 dimension of interest. using spline_domain_type = ddc::DiscreteDomain; @@ -84,16 +89,26 @@ class SplineEvaluator /** * @brief The type of the batch domain (obtained by removing the dimension of interest * from the whole domain). + * + * @tparam The batched discrete domain on which the interpolation points are defined. */ - using batch_domain_type = ddc:: - remove_dims_of_t; + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using batch_domain_type + = ddc::remove_dims_of_t; /** * @brief The type of the whole spline domain (cartesian product of 1D spline domain * and batch domain) preserving the 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 = ddc::replace_dim_of_t< - batched_evaluation_domain_type, + BatchedInterpolationDDom, evaluation_discrete_dimension_type, bsplines_type>; @@ -260,26 +275,35 @@ class SplineEvaluator * the set of 1D spline coefficients retained to perform the evaluation). * @param[in] spline_coef A ChunkSpan storing the spline coefficients. */ - template + template < + class Layout1, + class Layout2, + class Layout3, + class BatchedInterpolationDDom, + class... CoordsDims> void operator()( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, ddc::ChunkSpan< ddc::Coordinate const, - batched_evaluation_domain_type, + BatchedInterpolationDDom, Layout2, memory_space> const coords_eval, - ddc::ChunkSpan const - spline_coef) const + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const { evaluation_domain_type const evaluation_domain(spline_eval.domain()); - batch_domain_type const batch_domain(spline_eval.domain()); + batch_domain_type const batch_domain(spline_eval.domain()); ddc::parallel_for_each( "ddc_splines_evaluate", exec_space(), batch_domain, - KOKKOS_CLASS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { + KOKKOS_CLASS_LAMBDA(typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { const auto spline_eval_1D = spline_eval[j]; const auto coords_eval_1D = coords_eval[j]; const auto spline_coef_1D = spline_coef[j]; @@ -307,21 +331,25 @@ class SplineEvaluator * of the mesh. * @param[in] spline_coef A ChunkSpan storing the spline coefficients. */ - template + template void operator()( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, - ddc::ChunkSpan const - spline_coef) const + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const { evaluation_domain_type const evaluation_domain(spline_eval.domain()); - batch_domain_type const batch_domain(spline_eval.domain()); + batch_domain_type const batch_domain(spline_eval.domain()); ddc::parallel_for_each( "ddc_splines_evaluate", exec_space(), batch_domain, - KOKKOS_CLASS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { + KOKKOS_CLASS_LAMBDA(typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { const auto spline_eval_1D = spline_eval[j]; const auto spline_coef_1D = spline_coef[j]; for (auto const i : evaluation_domain) { @@ -370,26 +398,35 @@ class SplineEvaluator * the set of 1D spline coefficients retained to perform the evaluation). * @param[in] spline_coef A ChunkSpan storing the spline coefficients. */ - template + template < + class Layout1, + class Layout2, + class Layout3, + class BatchedInterpolationDDom, + class... CoordsDims> void deriv( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, ddc::ChunkSpan< ddc::Coordinate const, - batched_evaluation_domain_type, + BatchedInterpolationDDom, Layout2, memory_space> const coords_eval, - ddc::ChunkSpan const - spline_coef) const + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const { evaluation_domain_type const evaluation_domain(spline_eval.domain()); - batch_domain_type const batch_domain(spline_eval.domain()); + batch_domain_type const batch_domain(spline_eval.domain()); ddc::parallel_for_each( "ddc_splines_differentiate", exec_space(), batch_domain, - KOKKOS_CLASS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { + KOKKOS_CLASS_LAMBDA(typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { const auto spline_eval_1D = spline_eval[j]; const auto coords_eval_1D = coords_eval[j]; const auto spline_coef_1D = spline_coef[j]; @@ -414,21 +451,25 @@ class SplineEvaluator * @param[out] spline_eval The derivatives of the spline function at the coordinates. * @param[in] spline_coef A ChunkSpan storing the spline coefficients. */ - template + template void deriv( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, - ddc::ChunkSpan const - spline_coef) const + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const { evaluation_domain_type const evaluation_domain(spline_eval.domain()); - batch_domain_type const batch_domain(spline_eval.domain()); + batch_domain_type const batch_domain(spline_eval.domain()); ddc::parallel_for_each( "ddc_splines_differentiate", exec_space(), batch_domain, - KOKKOS_CLASS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { + KOKKOS_CLASS_LAMBDA(typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { const auto spline_eval_1D = spline_eval[j]; const auto spline_coef_1D = spline_coef[j]; for (auto const i : evaluation_domain) { @@ -453,13 +494,21 @@ class SplineEvaluator * points represented by this domain are unused and irrelevant. * @param[in] spline_coef A ChunkSpan storing the spline coefficients. */ - template + template void integrate( - ddc::ChunkSpan const integrals, - ddc::ChunkSpan const + ddc::ChunkSpan const integrals, + ddc::ChunkSpan const spline_coef) const { - batch_domain_type const batch_domain(integrals.domain()); + static_assert( + ddc::in_tags_v>, + "The spline coefficients domain must contain the bsplines dimension"); + using batch_domain_type = ddc::remove_dims_of_t; + static_assert( + std::is_same_v, + "The integrals domain must only contain the batch dimensions"); + + BatchedDDom const batch_domain(integrals.domain()); ddc::Chunk values_alloc( ddc::DiscreteDomain(spline_coef.domain()), ddc::KokkosAllocator()); @@ -470,7 +519,7 @@ class SplineEvaluator "ddc_splines_integrate", exec_space(), batch_domain, - KOKKOS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { + KOKKOS_LAMBDA(typename BatchedDDom::discrete_element_type const j) { integrals(j) = 0; for (typename spline_domain_type::discrete_element_type const i : values.domain()) { diff --git a/include/ddc/kernels/splines/spline_evaluator_2d.hpp b/include/ddc/kernels/splines/spline_evaluator_2d.hpp index 4994929ae..681b14683 100644 --- a/include/ddc/kernels/splines/spline_evaluator_2d.hpp +++ b/include/ddc/kernels/splines/spline_evaluator_2d.hpp @@ -32,7 +32,6 @@ namespace ddc { * @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 DDimX A variadic template of all the discrete dimensions forming the full space (EvaluationDDim1 + EvaluationDDim2 + batched dimensions). */ template < class ExecSpace, @@ -44,8 +43,7 @@ template < class LowerExtrapolationRule1, class UpperExtrapolationRule1, class LowerExtrapolationRule2, - class UpperExtrapolationRule2, - class... DDimX> + class UpperExtrapolationRule2> class SplineEvaluator2D { private: @@ -99,8 +97,15 @@ class SplineEvaluator2D evaluation_discrete_dimension_type1, evaluation_discrete_dimension_type2>; - /// @brief The type of the whole domain representing evaluation points. - using batched_evaluation_domain_type = ddc::DiscreteDomain; + /** + * @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; @@ -114,19 +119,29 @@ class SplineEvaluator2D /** * @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< - batched_evaluation_domain_type, + BatchedInterpolationDDom, evaluation_discrete_dimension_type1, evaluation_discrete_dimension_type2>; /** * @brief The type of the whole spline domain (cartesian product of 2D 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::to_type_seq_t, ddc::detail::TypeSeq< evaluation_discrete_dimension_type1, evaluation_discrete_dimension_type2>, @@ -372,26 +387,35 @@ class SplineEvaluator2D * the set of 2D spline coefficients retained to perform the evaluation). * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. */ - template + template < + class Layout1, + class Layout2, + class Layout3, + class BatchedInterpolationDDom, + class... CoordsDims> void operator()( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, ddc::ChunkSpan< ddc::Coordinate const, - batched_evaluation_domain_type, + BatchedInterpolationDDom, Layout2, memory_space> const coords_eval, - ddc::ChunkSpan const - spline_coef) const + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const { - batch_domain_type const batch_domain(coords_eval.domain()); + 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()); ddc::parallel_for_each( "ddc_splines_evaluate_2d", exec_space(), batch_domain, - KOKKOS_CLASS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { + KOKKOS_CLASS_LAMBDA(typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { const auto spline_eval_2D = spline_eval[j]; const auto coords_eval_2D = coords_eval[j]; const auto spline_coef_2D = spline_coef[j]; @@ -418,21 +442,25 @@ class SplineEvaluator2D * @param[out] spline_eval The values of the 2D spline function at their coordinates. * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. */ - template + template void operator()( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, - ddc::ChunkSpan const - spline_coef) const + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const { - batch_domain_type const batch_domain(spline_eval.domain()); + 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()); ddc::parallel_for_each( "ddc_splines_evaluate_2d", exec_space(), batch_domain, - KOKKOS_CLASS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { + KOKKOS_CLASS_LAMBDA(typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { const auto spline_eval_2D = spline_eval[j]; const auto spline_coef_2D = spline_coef[j]; for (auto const i1 : evaluation_domain1) { @@ -594,26 +622,35 @@ class SplineEvaluator2D * the set of 2D spline coefficients retained to perform the evaluation). * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. */ - template + template < + class Layout1, + class Layout2, + class Layout3, + class BatchedInterpolationDDom, + class... CoordsDims> void deriv_dim_1( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, ddc::ChunkSpan< ddc::Coordinate const, - batched_evaluation_domain_type, + BatchedInterpolationDDom, Layout2, memory_space> const coords_eval, - ddc::ChunkSpan const - spline_coef) const + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const { - batch_domain_type const batch_domain(coords_eval.domain()); + 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()); ddc::parallel_for_each( "ddc_splines_differentiate_2d_dim_1", exec_space(), batch_domain, - KOKKOS_CLASS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { + KOKKOS_CLASS_LAMBDA(typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { const auto spline_eval_2D = spline_eval[j]; const auto coords_eval_2D = coords_eval[j]; const auto spline_coef_2D = spline_coef[j]; @@ -640,21 +677,25 @@ class SplineEvaluator2D * @param[out] spline_eval The derivatives of the 2D spline function at the desired coordinates. * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. */ - template + template void deriv_dim_1( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, - ddc::ChunkSpan const - spline_coef) const + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const { - batch_domain_type const batch_domain(spline_eval.domain()); + 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()); ddc::parallel_for_each( "ddc_splines_differentiate_2d_dim_1", exec_space(), batch_domain, - KOKKOS_CLASS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { + KOKKOS_CLASS_LAMBDA(typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { const auto spline_eval_2D = spline_eval[j]; const auto spline_coef_2D = spline_coef[j]; for (auto const i1 : evaluation_domain1) { @@ -687,26 +728,35 @@ class SplineEvaluator2D * the set of 2D spline coefficients retained to perform the evaluation). * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. */ - template + template < + class Layout1, + class Layout2, + class Layout3, + class BatchedInterpolationDDom, + class... CoordsDims> void deriv_dim_2( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, ddc::ChunkSpan< ddc::Coordinate const, - batched_evaluation_domain_type, + BatchedInterpolationDDom, Layout2, memory_space> const coords_eval, - ddc::ChunkSpan const - spline_coef) const + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const { - batch_domain_type const batch_domain(coords_eval.domain()); + 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()); ddc::parallel_for_each( "ddc_splines_differentiate_2d_dim_2", exec_space(), batch_domain, - KOKKOS_CLASS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { + KOKKOS_CLASS_LAMBDA(typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { const auto spline_eval_2D = spline_eval[j]; const auto coords_eval_2D = coords_eval[j]; const auto spline_coef_2D = spline_coef[j]; @@ -733,21 +783,25 @@ class SplineEvaluator2D * @param[out] spline_eval The derivatives of the 2D spline function at the desired coordinates. * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. */ - template + template void deriv_dim_2( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, - ddc::ChunkSpan const - spline_coef) const + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const { - batch_domain_type const batch_domain(spline_eval.domain()); + 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()); ddc::parallel_for_each( "ddc_splines_differentiate_2d_dim_2", exec_space(), batch_domain, - KOKKOS_CLASS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { + KOKKOS_CLASS_LAMBDA(typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { const auto spline_eval_2D = spline_eval[j]; const auto spline_coef_2D = spline_coef[j]; for (auto const i1 : evaluation_domain1) { @@ -780,26 +834,35 @@ class SplineEvaluator2D * the set of 2D spline coefficients retained to perform the evaluation). * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. */ - template + template < + class Layout1, + class Layout2, + class Layout3, + class BatchedInterpolationDDom, + class... CoordsDims> void deriv_1_and_2( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, ddc::ChunkSpan< ddc::Coordinate const, - batched_evaluation_domain_type, + BatchedInterpolationDDom, Layout2, memory_space> const coords_eval, - ddc::ChunkSpan const - spline_coef) const + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const { - batch_domain_type const batch_domain(coords_eval.domain()); + 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()); ddc::parallel_for_each( "ddc_splines_cross_differentiate", exec_space(), batch_domain, - KOKKOS_CLASS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { + KOKKOS_CLASS_LAMBDA(typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { const auto spline_eval_2D = spline_eval[j]; const auto coords_eval_2D = coords_eval[j]; const auto spline_coef_2D = spline_coef[j]; @@ -826,21 +889,25 @@ class SplineEvaluator2D * @param[out] spline_eval The cross-derivatives of the 2D spline function at the desired coordinates. * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. */ - template + template void deriv_1_and_2( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, - ddc::ChunkSpan const - spline_coef) const + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const { - batch_domain_type const batch_domain(spline_eval.domain()); + 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()); ddc::parallel_for_each( "ddc_splines_cross_differentiate", exec_space(), batch_domain, - KOKKOS_CLASS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { + KOKKOS_CLASS_LAMBDA(typename batch_domain_type< + BatchedInterpolationDDom>::discrete_element_type const j) { const auto spline_eval_2D = spline_eval[j]; const auto spline_coef_2D = spline_coef[j]; for (auto const i1 : evaluation_domain1) { @@ -874,17 +941,26 @@ class SplineEvaluator2D * the set of 2D spline coefficients retained to perform the evaluation). * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. */ - template + template < + class InterestDim, + class Layout1, + class Layout2, + class Layout3, + class BatchedInterpolationDDom, + class... CoordsDims> void deriv( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, ddc::ChunkSpan< ddc::Coordinate const, - batched_evaluation_domain_type, + BatchedInterpolationDDom, Layout2, memory_space> const coords_eval, - ddc::ChunkSpan const - spline_coef) const + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const { static_assert( std::is_same_v< @@ -918,12 +994,15 @@ class SplineEvaluator2D * @param[out] spline_eval The derivatives of the 2D spline function at the desired coordinates. * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. */ - template + template void deriv( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, - ddc::ChunkSpan const - spline_coef) const + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const { static_assert( std::is_same_v< @@ -972,17 +1051,21 @@ class SplineEvaluator2D class Layout1, class Layout2, class Layout3, + class BatchedInterpolationDDom, class... CoordsDims> void deriv2( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, ddc::ChunkSpan< ddc::Coordinate const, - batched_evaluation_domain_type, + BatchedInterpolationDDom, Layout2, memory_space> const coords_eval, - ddc::ChunkSpan const - spline_coef) const + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout3, + memory_space> const spline_coef) const { static_assert( (std::is_same_v< @@ -1014,12 +1097,20 @@ class SplineEvaluator2D * @param[out] spline_eval The derivatives of the 2D spline function at the desired coordinates. * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. */ - template + template < + class InterestDim1, + class InterestDim2, + class Layout1, + class Layout2, + class BatchedInterpolationDDom> void deriv2( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, - ddc::ChunkSpan const - spline_coef) const + ddc::ChunkSpan< + double const, + batched_spline_domain_type, + Layout2, + memory_space> const spline_coef) const { static_assert( (std::is_same_v< @@ -1046,12 +1137,23 @@ class SplineEvaluator2D * points represented by this domain are unused and irrelevant. * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. */ - template + template void integrate( - ddc::ChunkSpan const integrals, - ddc::ChunkSpan const + 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()), diff --git a/tests/splines/CMakeLists.txt b/tests/splines/CMakeLists.txt index 91208ad4c..3e1c99f83 100644 --- a/tests/splines/CMakeLists.txt +++ b/tests/splines/CMakeLists.txt @@ -142,3 +142,47 @@ foreach(BC "BC_PERIODIC" "BC_GREVILLE" "BC_HERMITE") 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}") + foreach(BSPLINES_TYPE "BSPLINES_TYPE_UNIFORM" "BSPLINES_TYPE_NON_UNIFORM") + set(test_name + "splines_tests_BATCHED_MULTIPLE_DOMAINS_DEGREE_X_${DEGREE_X}_${BSPLINES_TYPE}_${BC}_${SOLVER}" + ) + add_executable("${test_name}" ../main.cpp multiple_batch_domain_spline_builder.cpp) + target_compile_definitions( + "${test_name}" + PUBLIC -DDEGREE_X=${DEGREE_X} -D${BSPLINES_TYPE} -D${BC} -DSOLVER_${SOLVER} + ) + 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(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 + "2d_splines_tests_BATCHED_MULTIPLE_DOMAINS_DEGREE_${DEGREE}_${BSPLINES_TYPE}_${EVALUATOR}_${BC}" + ) + add_executable( + "${test_name}" + ../main.cpp + multiple_batch_domain_2d_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() diff --git a/tests/splines/batched_2d_spline_builder.cpp b/tests/splines/batched_2d_spline_builder.cpp index ce74e6938..815ade764 100644 --- a/tests/splines/batched_2d_spline_builder.cpp +++ b/tests/splines/batched_2d_spline_builder.cpp @@ -186,6 +186,8 @@ void Batched2dSplineTest() = GrevillePoints>::template get_domain(); ddc::DiscreteDomain const interpolation_domain2 = GrevillePoints>::template get_domain(); + ddc::DiscreteDomain const + interpolation_domain(interpolation_domain1, interpolation_domain2); // The following line creates a discrete domain over all dimensions (DDims...) except DDimI1 and DDimI2. auto const dom_vals_tmp = ddc::remove_dims_of_t, DDimI1, DDimI2>( ddc::DiscreteDomain(DElem(0), DVect(ncells))...); @@ -220,13 +222,12 @@ void Batched2dSplineTest() s_bcr, s_bcl, s_bcr, - ddc::SplineSolver::GINKGO, - DDims...> const spline_builder(dom_vals); + 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(); + 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()); @@ -472,8 +473,7 @@ void Batched2dSplineTest() extrapolation_rule_1_type, extrapolation_rule_1_type, extrapolation_rule_2_type, - extrapolation_rule_2_type, - DDims...> const + extrapolation_rule_2_type> const spline_evaluator( extrapolation_rule_1, extrapolation_rule_1, diff --git a/tests/splines/batched_spline_builder.cpp b/tests/splines/batched_spline_builder.cpp index f4e2ab7fa..9e3827186 100644 --- a/tests/splines/batched_spline_builder.cpp +++ b/tests/splines/batched_spline_builder.cpp @@ -204,16 +204,16 @@ void BatchedSplineTest() s_bcl, s_bcr, #if defined(SOLVER_LAPACK) - ddc::SplineSolver::LAPACK, + ddc::SplineSolver::LAPACK #elif defined(SOLVER_GINKGO) - ddc::SplineSolver::GINKGO, + ddc::SplineSolver::GINKGO #endif - DDims...> const spline_builder(dom_vals); + > 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_batch = spline_builder.batch_domain(); - auto const dom_spline = spline_builder.batched_spline_domain(); + auto const dom_batch = spline_builder.batch_domain(dom_vals); + 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()); @@ -308,8 +308,8 @@ void BatchedSplineTest() BSplines, DDimI, extrapolation_rule_type, - extrapolation_rule_type, - DDims...> const spline_evaluator_batched(extrapolation_rule, extrapolation_rule); + extrapolation_rule_type> const + spline_evaluator_batched(extrapolation_rule, extrapolation_rule); // Instantiate chunk of coordinates of dom_interpolation ddc::Chunk coords_eval_alloc(dom_vals, ddc::KokkosAllocator, MemorySpace>()); @@ -359,8 +359,8 @@ void BatchedSplineTest() spline_eval_integrals.domain(), 0., ddc::reducer::max(), - KOKKOS_LAMBDA(typename decltype(spline_builder)::batch_domain_type:: - discrete_element_type const e) { + KOKKOS_LAMBDA(typename decltype(spline_builder)::template batch_domain_type< + ddc::DiscreteDomain>::discrete_element_type const e) { return Kokkos::abs( spline_eval_integrals(e) - evaluator.deriv(xN(), -1) + evaluator.deriv(x0(), -1)); diff --git a/tests/splines/extrapolation_rule.cpp b/tests/splines/extrapolation_rule.cpp index a6574aa19..6501c57d9 100644 --- a/tests/splines/extrapolation_rule.cpp +++ b/tests/splines/extrapolation_rule.cpp @@ -191,13 +191,12 @@ void ExtrapolationRuleSplineTest() s_bcr1, s_bcl2, s_bcr2, - ddc::SplineSolver::GINKGO, - DDims...> const spline_builder(dom_vals); + 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(); + 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()); @@ -267,8 +266,7 @@ void ExtrapolationRuleSplineTest() extrapolation_rule_dim_1_type, extrapolation_rule_dim_1_type, extrapolation_rule_dim_2_type, - extrapolation_rule_dim_2_type, - DDims...> const + extrapolation_rule_dim_2_type> const spline_evaluator_batched( extrapolation_rule_left_dim_1, extrapolation_rule_right_dim_1, diff --git a/tests/splines/multiple_batch_domain_2d_spline_builder.cpp b/tests/splines/multiple_batch_domain_2d_spline_builder.cpp new file mode 100644 index 000000000..6bb4fb355 --- /dev/null +++ b/tests/splines/multiple_batch_domain_2d_spline_builder.cpp @@ -0,0 +1,757 @@ +// Copyright (C) The DDC development team, see COPYRIGHT.md file +// +// SPDX-License-Identifier: MIT + +#include +#include +#include +#if defined(BC_HERMITE) +#include +#endif +#if defined(BSPLINES_TYPE_UNIFORM) +#include +#endif +#include + +#include +#include + +#include + +#include + +#include "evaluator_2d.hpp" +#if defined(BC_PERIODIC) +#include "cosine_evaluator.hpp" +#else +#include "polynomial_evaluator.hpp" +#endif +#include "spline_error_bounds.hpp" + +namespace DDC_HIP_5_7_ANONYMOUS_NAMESPACE_WORKAROUND(BATCHED_2D_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 +{ +}; + +struct DDimBatchExtra +{ +}; + +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 = Evaluator2D:: + Evaluator, CosineEvaluator::Evaluator>; +#else +template +using evaluator_type = Evaluator2D::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()); +} + +/// Computes the evaluation error when evaluating a 2D spline on its interpolation points. +template < + typename ExecSpace, + typename MemorySpace, + typename DDimI1, + typename DDimI2, + typename Builder, + typename Evaluator, + typename... DDims> +std::tuple ComputeEvaluationError( + ExecSpace const& exec_space, + ddc::DiscreteDomain const& dom_vals, + Builder const& spline_builder, + Evaluator const& spline_evaluator, + evaluator_type const& evaluator) +{ + using I1 = typename DDimI1::continuous_dimension_type; + using I2 = typename DDimI2::continuous_dimension_type; + +#if defined(BC_HERMITE) + ddc::DiscreteDomain const interpolation_domain1(dom_vals); + ddc::DiscreteDomain const interpolation_domain2(dom_vals); + // 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 + + // 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(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) + // 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) + 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 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_deriv12_alloc(dom_vals, ddc::KokkosAllocator()); + ddc::ChunkSpan const spline_eval_deriv12 = spline_eval_deriv12_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 deriv2< + I1, + I2>(spline_eval_deriv12, 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)); + return Kokkos::abs(spline_eval_deriv1(e) - evaluator.deriv(x, y, 1, 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)); + return Kokkos::abs(spline_eval_deriv2(e) - evaluator.deriv(x, y, 0, 1)); + }); + double const max_norm_error_diff12 = 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)); + return Kokkos::abs(spline_eval_deriv12(e) - evaluator.deriv(x, y, 1, 1)); + }); + + return std::make_tuple( + max_norm_error, + max_norm_error_diff1, + max_norm_error_diff2, + max_norm_error_diff12); +} + +// 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... DDims> +void MultipleBatchDomain2dSplineTest() +{ + using I1 = typename DDimI1::continuous_dimension_type; + using I2 = typename DDimI2::continuous_dimension_type; + + // Instantiate execution spaces and initialize spaces + ExecSpace const exec_space; + std::size_t const ncells = 10; + 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_domain(interpolation_domain1, interpolation_domain2); + // The following line creates a discrete domain over all dimensions (DDims...) except DDimI1 and DDimI2. + auto const dom_vals_tmp = ddc::remove_dims_of_t, DDimI1, DDimI2>( + ddc::DiscreteDomain(DElem(0), DVect(ncells))...); + ddc::DiscreteDomain const + dom_vals(dom_vals_tmp, interpolation_domain1, interpolation_domain2); + + ddc::DiscreteDomain const + extra_batch_domain(DElem(0), DVect(ncells)); + ddc::DiscreteDomain const dom_vals_extra( + dom_vals_tmp, + interpolation_domain1, + interpolation_domain2, + extra_batch_domain); + + // Create a SplineBuilder over BSplines and BSplines and using some boundary conditions + ddc::SplineBuilder2D< + ExecSpace, + MemorySpace, + BSplines, + BSplines, + DDimI1, + DDimI2, + s_bcl, + s_bcr, + s_bcl, + s_bcr, + ddc::SplineSolver::GINKGO> const spline_builder(interpolation_domain); + + evaluator_type const evaluator(spline_builder.interpolation_domain()); + + // Instantiate a SplineEvaluator over interest dimensions +#if defined(BC_PERIODIC) + using extrapolation_rule_1_type = ddc::PeriodicExtrapolationRule; + using extrapolation_rule_2_type = ddc::PeriodicExtrapolationRule; +#else + using extrapolation_rule_1_type = ddc::NullExtrapolationRule; + using extrapolation_rule_2_type = ddc::NullExtrapolationRule; +#endif + extrapolation_rule_1_type const extrapolation_rule_1; + extrapolation_rule_2_type const extrapolation_rule_2; + + ddc::SplineEvaluator2D< + ExecSpace, + MemorySpace, + BSplines, + BSplines, + DDimI1, + DDimI2, + extrapolation_rule_1_type, + extrapolation_rule_1_type, + extrapolation_rule_2_type, + extrapolation_rule_2_type> const + spline_evaluator( + extrapolation_rule_1, + extrapolation_rule_1, + extrapolation_rule_2, + extrapolation_rule_2); + + // Check the evaluation error for the first domain + auto const [max_norm_error, max_norm_error_diff1, max_norm_error_diff2, max_norm_error_diff12] + = ComputeEvaluationError< + ExecSpace, + MemorySpace>(exec_space, dom_vals, spline_builder, spline_evaluator, evaluator); + + double const max_norm = evaluator.max_norm(); + double const max_norm_diff1 = evaluator.max_norm(1, 0); + double const max_norm_diff2 = evaluator.max_norm(0, 1); + double const max_norm_diff12 = evaluator.max_norm(1, 1); + + SplineErrorBounds> const error_bounds(evaluator); + EXPECT_LE( + max_norm_error, + std:: + max(error_bounds + .error_bound(dx(ncells), dx(ncells), 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), + 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), + s_degree, + s_degree), + 1e-12 * max_norm_diff2)); + EXPECT_LE( + max_norm_error_diff12, + std:: + max(error_bounds.error_bound_on_deriv_12( + dx(ncells), + dx(ncells), + s_degree, + s_degree), + 1e-11 * max_norm_diff12)); + + // Check the evaluation error for the domain with an additional batch dimension + auto const + [max_norm_error_extra, + max_norm_error_diff1_extra, + max_norm_error_diff2_extra, + max_norm_error_diff12_extra] + = ComputeEvaluationError( + exec_space, + dom_vals_extra, + spline_builder, + spline_evaluator, + evaluator); + + double const max_norm_extra = evaluator.max_norm(); + double const max_norm_diff1_extra = evaluator.max_norm(1, 0); + double const max_norm_diff2_extra = evaluator.max_norm(0, 1); + double const max_norm_diff12_extra = evaluator.max_norm(1, 1); + + SplineErrorBounds> const error_bounds_extra(evaluator); + EXPECT_LE( + max_norm_error_extra, + std:: + max(error_bounds_extra + .error_bound(dx(ncells), dx(ncells), s_degree, s_degree), + 1.0e-14 * max_norm_extra)); + EXPECT_LE( + max_norm_error_diff1_extra, + std:: + max(error_bounds_extra.error_bound_on_deriv_1( + dx(ncells), + dx(ncells), + s_degree, + s_degree), + 1e-12 * max_norm_diff1_extra)); + EXPECT_LE( + max_norm_error_diff2_extra, + std:: + max(error_bounds_extra.error_bound_on_deriv_2( + dx(ncells), + dx(ncells), + s_degree, + s_degree), + 1e-12 * max_norm_diff2_extra)); + EXPECT_LE( + max_norm_error_diff12_extra, + std:: + max(error_bounds_extra.error_bound_on_deriv_12( + dx(ncells), + dx(ncells), + s_degree, + s_degree), + 1e-11 * max_norm_diff12_extra)); +} + +} // namespace DDC_HIP_5_7_ANONYMOUS_NAMESPACE_WORKAROUND(BATCHED_2D_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(MultipleBatchDomain2dSpline), 2DXY) +{ + MultipleBatchDomain2dSplineTest< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + DDimGPS, + DDimGPS, + DDimGPS, + DDimGPS>(); +} + +TEST(SUFFIX(MultipleBatchDomain2dSpline), 3DXY) +{ + MultipleBatchDomain2dSplineTest< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + DDimGPS, + DDimGPS, + DDimGPS, + DDimGPS, + DDimBatch>(); +} + +TEST(SUFFIX(MultipleBatchDomain2dSpline), 3DXZ) +{ + MultipleBatchDomain2dSplineTest< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + DDimGPS, + DDimGPS, + DDimGPS, + DDimBatch, + DDimGPS>(); +} + +TEST(SUFFIX(MultipleBatchDomain2dSpline), 3DYZ) +{ + MultipleBatchDomain2dSplineTest< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + DDimGPS, + DDimGPS, + DDimBatch, + DDimGPS, + DDimGPS>(); +} diff --git a/tests/splines/multiple_batch_domain_spline_builder.cpp b/tests/splines/multiple_batch_domain_spline_builder.cpp new file mode 100644 index 000000000..1b59685dc --- /dev/null +++ b/tests/splines/multiple_batch_domain_spline_builder.cpp @@ -0,0 +1,485 @@ +// Copyright (C) The DDC development team, see COPYRIGHT.md file +// +// SPDX-License-Identifier: MIT + +#include +#include +#include +#if defined(BC_HERMITE) +#include +#endif +#if defined(BSPLINES_TYPE_UNIFORM) +#include +#endif +#include + +#include +#include + +#include + +#include + +#include "cosine_evaluator.hpp" +#include "spline_error_bounds.hpp" + +namespace DDC_HIP_5_7_ANONYMOUS_NAMESPACE_WORKAROUND(BATCHED_SPLINE_BUILDER_CPP) { + +#if defined(BC_PERIODIC) +struct DimX +{ + static constexpr bool PERIODIC = true; +}; + +struct DimY +{ + static constexpr bool PERIODIC = true; +}; +#else + +struct DimX +{ + static constexpr bool PERIODIC = false; +}; + +struct DimY +{ + static constexpr bool PERIODIC = false; +}; +#endif + +struct DDimBatch1 +{ +}; + +struct DDimBatchExtra +{ +}; + +constexpr std::size_t s_degree_x = DEGREE_X; + +#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 dimension of interest, the discrete dimension is deduced from the Greville points type. +template +struct DDimGPS : GrevillePoints>::interpolation_discrete_dimension_type +{ +}; + +template +using evaluator_type = CosineEvaluator::Evaluator; + +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()); +} + +/// Computes the evaluation error when evaluating a spline on its interpolation points. +template < + typename ExecSpace, + typename MemorySpace, + typename DDimI, + class Builder, + class Evaluator, + typename... DDims> +std::tuple ComputeEvaluationError( + ExecSpace const& exec_space, + ddc::DiscreteDomain const& dom_vals, + Builder const& spline_builder, + Evaluator const& spline_evaluator, + evaluator_type const& evaluator) +{ + using I = typename DDimI::continuous_dimension_type; + +#if defined(BC_HERMITE) + // Create the derivs domains + ddc::DiscreteDomain> const + derivs_domain(DElem>(1), DVect>(s_degree_x / 2)); + auto const dom_derivs = ddc::replace_dim_of>(dom_vals, derivs_domain); +#endif + + // Compute useful domains (dom_interpolation, dom_batch, dom_bsplines and dom_spline) + ddc::DiscreteDomain const dom_interpolation = spline_builder.interpolation_domain(); + auto const dom_batch = spline_builder.batch_domain(dom_vals); + 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(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) + // Allocate and fill a chunk containing derivs to be passed as input to spline_builder. + int const shift = s_degree_x % 2; // shift = 0 for even order, 1 for odd order + ddc::Chunk derivs_lhs_alloc(dom_derivs, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs_lhs = derivs_lhs_alloc.span_view(); + if (s_bcl == ddc::BoundCond::HERMITE) { + ddc::Chunk derivs_lhs1_host_alloc(derivs_domain, ddc::HostAllocator()); + ddc::ChunkSpan const derivs_lhs1_host = derivs_lhs1_host_alloc.span_view(); + for (int ii = 1; ii < derivs_lhs1_host.domain().template extent>() + 1; + ++ii) { + derivs_lhs1_host( + typename decltype(derivs_lhs1_host.domain())::discrete_element_type(ii)) + = evaluator.deriv(x0(), ii + shift - 1); + } + auto derivs_lhs1_alloc = ddc::create_mirror_view_and_copy(exec_space, derivs_lhs1_host); + ddc::ChunkSpan const derivs_lhs1 = derivs_lhs1_alloc.span_view(); + ddc::parallel_for_each( + exec_space, + derivs_lhs.domain(), + KOKKOS_LAMBDA( + typename decltype(derivs_lhs.domain())::discrete_element_type const e) { + derivs_lhs(e) = derivs_lhs1(DElem>(e)); + }); + } + + ddc::Chunk derivs_rhs_alloc(dom_derivs, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs_rhs = derivs_rhs_alloc.span_view(); + if (s_bcr == ddc::BoundCond::HERMITE) { + ddc::Chunk derivs_rhs1_host_alloc(derivs_domain, ddc::HostAllocator()); + ddc::ChunkSpan const derivs_rhs1_host = derivs_rhs1_host_alloc.span_view(); + for (int ii = 1; ii < derivs_rhs1_host.domain().template extent>() + 1; + ++ii) { + derivs_rhs1_host( + typename decltype(derivs_rhs1_host.domain())::discrete_element_type(ii)) + = evaluator.deriv(xN(), ii + shift - 1); + } + auto derivs_rhs1_alloc = ddc::create_mirror_view_and_copy(exec_space, derivs_rhs1_host); + ddc::ChunkSpan const derivs_rhs1 = derivs_rhs1_alloc.span_view(); + + ddc::parallel_for_each( + exec_space, + derivs_rhs.domain(), + KOKKOS_LAMBDA( + typename decltype(derivs_rhs.domain())::discrete_element_type const e) { + derivs_rhs(e) = derivs_rhs1(DElem>(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) + spline_builder( + coef, + vals.span_cview(), + std::optional(derivs_lhs.span_cview()), + std::optional(derivs_rhs.span_cview())); +#else + spline_builder(coef, vals.span_cview()); +#endif + + // 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_deriv_alloc(dom_vals, ddc::KokkosAllocator()); + ddc::ChunkSpan const spline_eval_deriv = spline_eval_deriv_alloc.span_view(); + ddc::Chunk spline_eval_integrals_alloc(dom_batch, ddc::KokkosAllocator()); + ddc::ChunkSpan const spline_eval_integrals = spline_eval_integrals_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.deriv(spline_eval_deriv, coords_eval.span_cview(), coef.span_cview()); + spline_evaluator.integrate(spline_eval_integrals, 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_diff = ddc::parallel_transform_reduce( + exec_space, + spline_eval_deriv.domain(), + 0., + ddc::reducer::max(), + KOKKOS_LAMBDA(DElem const e) { + Coord const x = ddc::coordinate(DElem(e)); + return Kokkos::abs(spline_eval_deriv(e) - evaluator.deriv(x, 1)); + }); + double const max_norm_error_integ = ddc::parallel_transform_reduce( + exec_space, + spline_eval_integrals.domain(), + 0., + ddc::reducer::max(), + KOKKOS_LAMBDA(typename Builder::template batch_domain_type< + ddc::DiscreteDomain>::discrete_element_type const e) { + return Kokkos::abs( + spline_eval_integrals(e) - evaluator.deriv(xN(), -1) + + evaluator.deriv(x0(), -1)); + }); + + return std::make_tuple(max_norm_error, max_norm_error_diff, max_norm_error_integ); +} + +// Checks that when evaluating the spline at interpolation points with +// multiple batch patterns using the same builders and evaluators, one +// recovers values that were used to build the spline +template +void MultipleBatchDomainSplineTest() +{ + using I = typename DDimI::continuous_dimension_type; + + // Instantiate execution spaces and initialize spaces + ExecSpace const exec_space; + + std::size_t const ncells = 10; + InterestDimInitializer(ncells); + + // Create the values domain (mesh) + ddc::DiscreteDomain const interpolation_domain + = GrevillePoints>::template get_domain(); + // The following line creates a discrete domain over all dimensions (DDims...) except DDimI. + auto const dom_vals_tmp = ddc::remove_dims_of_t, DDimI>( + ddc::DiscreteDomain(DElem(0), DVect(ncells))...); + ddc::DiscreteDomain const dom_vals(dom_vals_tmp, interpolation_domain); + + // Create a discrete domain with an additional batch dimension + ddc::DiscreteDomain const + extra_domain(DElem(0), DVect(ncells)); + ddc::DiscreteDomain const + dom_vals_extra(dom_vals_tmp, interpolation_domain, extra_domain); + + // Create a SplineBuilder over BSplines using some boundary conditions + ddc::SplineBuilder< + ExecSpace, + MemorySpace, + BSplines, + DDimI, + s_bcl, + s_bcr, +#if defined(SOLVER_LAPACK) + ddc::SplineSolver::LAPACK +#elif defined(SOLVER_GINKGO) + ddc::SplineSolver::GINKGO +#endif + > const spline_builder(interpolation_domain); + + // Instantiate a SplineEvaluator over interest dimension +#if defined(BC_PERIODIC) + using extrapolation_rule_type = ddc::PeriodicExtrapolationRule; +#else + using extrapolation_rule_type = ddc::NullExtrapolationRule; +#endif + extrapolation_rule_type const extrapolation_rule; + + ddc::SplineEvaluator< + ExecSpace, + MemorySpace, + BSplines, + DDimI, + extrapolation_rule_type, + extrapolation_rule_type> const + spline_evaluator_batched(extrapolation_rule, extrapolation_rule); + + evaluator_type const evaluator(spline_builder.interpolation_domain()); + + // Check the evaluation error for the first domain + auto const [max_norm_error, max_norm_error_diff, max_norm_error_integ] = ComputeEvaluationError< + ExecSpace, + MemorySpace>(exec_space, dom_vals, spline_builder, spline_evaluator_batched, evaluator); + + double const max_norm = evaluator.max_norm(); + double const max_norm_diff = evaluator.max_norm(1); + double const max_norm_int = evaluator.max_norm(-1); + + SplineErrorBounds> const error_bounds(evaluator); + EXPECT_LE( + max_norm_error, + std::max(error_bounds.error_bound(dx(ncells), s_degree_x), 1.0e-14 * max_norm)); + EXPECT_LE( + max_norm_error_diff, + std:: + max(error_bounds.error_bound_on_deriv(dx(ncells), s_degree_x), + 1e-12 * max_norm_diff)); + EXPECT_LE( + max_norm_error_integ, + std:: + max(error_bounds.error_bound_on_int(dx(ncells), s_degree_x), + 1.0e-14 * max_norm_int)); + + // Check the evaluation error for the domain with an additional batch dimension + auto const [max_norm_error_extra, max_norm_error_diff_extra, max_norm_error_integ_extra] + = ComputeEvaluationError( + exec_space, + dom_vals_extra, + spline_builder, + spline_evaluator_batched, + evaluator); + + double const max_norm_extra = evaluator.max_norm(); + double const max_norm_diff_extra = evaluator.max_norm(1); + double const max_norm_int_extra = evaluator.max_norm(-1); + + SplineErrorBounds> const error_bounds_extra(evaluator); + EXPECT_LE( + max_norm_error_extra, + std:: + max(error_bounds_extra.error_bound(dx(ncells), s_degree_x), + 1.0e-14 * max_norm_extra)); + EXPECT_LE( + max_norm_error_diff_extra, + std:: + max(error_bounds_extra.error_bound_on_deriv(dx(ncells), s_degree_x), + 1e-12 * max_norm_diff_extra)); + EXPECT_LE( + max_norm_error_integ_extra, + std:: + max(error_bounds_extra.error_bound_on_int(dx(ncells), s_degree_x), + 1.0e-14 * max_norm_int_extra)); +} + + +} // namespace DDC_HIP_5_7_ANONYMOUS_NAMESPACE_WORKAROUND(BATCHED_SPLINE_BUILDER_CPP) + +#if defined(BC_PERIODIC) && defined(BSPLINES_TYPE_UNIFORM) && defined(SOLVER_LAPACK) +#define SUFFIX(name) name##Lapack##Periodic##Uniform +#elif defined(BC_PERIODIC) && defined(BSPLINES_TYPE_NON_UNIFORM) && defined(SOLVER_LAPACK) +#define SUFFIX(name) name##Lapack##Periodic##NonUniform +#elif defined(BC_GREVILLE) && defined(BSPLINES_TYPE_UNIFORM) && defined(SOLVER_LAPACK) +#define SUFFIX(name) name##Lapack##Greville##Uniform +#elif defined(BC_GREVILLE) && defined(BSPLINES_TYPE_NON_UNIFORM) && defined(SOLVER_LAPACK) +#define SUFFIX(name) name##Lapack##Greville##NonUniform +#elif defined(BC_HERMITE) && defined(BSPLINES_TYPE_UNIFORM) && defined(SOLVER_LAPACK) +#define SUFFIX(name) name##Lapack##Hermite##Uniform +#elif defined(BC_HERMITE) && defined(BSPLINES_TYPE_NON_UNIFORM) && defined(SOLVER_LAPACK) +#define SUFFIX(name) name##Lapack##Hermite##NonUniform +#elif defined(BC_PERIODIC) && defined(BSPLINES_TYPE_UNIFORM) && defined(SOLVER_GINKGO) +#define SUFFIX(name) name##Ginkgo##Periodic##Uniform +#elif defined(BC_PERIODIC) && defined(BSPLINES_TYPE_NON_UNIFORM) && defined(SOLVER_GINKGO) +#define SUFFIX(name) name##Ginkgo##Periodic##NonUniform +#elif defined(BC_GREVILLE) && defined(BSPLINES_TYPE_UNIFORM) && defined(SOLVER_GINKGO) +#define SUFFIX(name) name##Ginkgo##Greville##Uniform +#elif defined(BC_GREVILLE) && defined(BSPLINES_TYPE_NON_UNIFORM) && defined(SOLVER_GINKGO) +#define SUFFIX(name) name##Ginkgo##Greville##NonUniform +#elif defined(BC_HERMITE) && defined(BSPLINES_TYPE_UNIFORM) && defined(SOLVER_GINKGO) +#define SUFFIX(name) name##Ginkgo##Hermite##Uniform +#elif defined(BC_HERMITE) && defined(BSPLINES_TYPE_NON_UNIFORM) && defined(SOLVER_GINKGO) +#define SUFFIX(name) name##Ginkgo##Hermite##NonUniform +#endif + +TEST(SUFFIX(MultipleBatchDomainSpline), 1DX) +{ + MultipleBatchDomainSplineTest< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + DDimGPS, + DDimGPS>(); +} + +TEST(SUFFIX(MultipleBatchDomainSpline), 2DX) +{ + MultipleBatchDomainSplineTest< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + DDimGPS, + DDimGPS, + DDimBatch1>(); +} + +TEST(SUFFIX(MultipleBatchDomainSpline), 2DY) +{ + MultipleBatchDomainSplineTest< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + DDimGPS, + DDimBatch1, + DDimGPS>(); +} diff --git a/tests/splines/non_periodic_spline_builder.cpp b/tests/splines/non_periodic_spline_builder.cpp index f8491c97a..88375a45e 100644 --- a/tests/splines/non_periodic_spline_builder.cpp +++ b/tests/splines/non_periodic_spline_builder.cpp @@ -119,8 +119,7 @@ void TestNonPeriodicSplineBuilderTestIdentity() DDimX, s_bcl, s_bcr, - ddc::SplineSolver::GINKGO, - DDimX> const spline_builder(interpolation_domain); + ddc::SplineSolver::GINKGO> const spline_builder(interpolation_domain); // 5. Allocate and fill a chunk over the interpolation domain ddc::Chunk yvals_alloc(interpolation_domain, ddc::KokkosAllocator()); @@ -177,8 +176,8 @@ void TestNonPeriodicSplineBuilderTestIdentity() BSplinesX, DDimX, ddc::NullExtrapolationRule, - ddc::NullExtrapolationRule, - DDimX> const spline_evaluator(extrapolation_rule, extrapolation_rule); + ddc::NullExtrapolationRule> const + spline_evaluator(extrapolation_rule, extrapolation_rule); ddc::Chunk coords_eval_alloc(interpolation_domain, ddc::KokkosAllocator()); @@ -200,8 +199,9 @@ void TestNonPeriodicSplineBuilderTestIdentity() spline_evaluator .deriv(spline_eval_deriv.span_view(), coords_eval.span_cview(), coef.span_cview()); - ddc::Chunk - integral(spline_builder.batch_domain(), ddc::KokkosAllocator()); + ddc::Chunk integral( + spline_builder.batch_domain(interpolation_domain), + ddc::KokkosAllocator()); spline_evaluator.integrate(integral.span_view(), coef.span_cview()); ddc::Chunk< diff --git a/tests/splines/periodic_spline_builder.cpp b/tests/splines/periodic_spline_builder.cpp index 5d61844a0..9216f91c1 100644 --- a/tests/splines/periodic_spline_builder.cpp +++ b/tests/splines/periodic_spline_builder.cpp @@ -94,8 +94,7 @@ void TestPeriodicSplineBuilderTestIdentity() DDimX, ddc::BoundCond::PERIODIC, ddc::BoundCond::PERIODIC, - ddc::SplineSolver::GINKGO, - DDimX> const spline_builder(interpolation_domain); + ddc::SplineSolver::GINKGO> const spline_builder(interpolation_domain); // 5. Allocate and fill a chunk over the interpolation domain ddc::Chunk yvals_alloc(interpolation_domain, ddc::KokkosAllocator()); @@ -117,8 +116,8 @@ void TestPeriodicSplineBuilderTestIdentity() BSplinesX, DDimX, ddc::PeriodicExtrapolationRule, - ddc::PeriodicExtrapolationRule, - DDimX> const spline_evaluator(periodic_extrapolation, periodic_extrapolation); + ddc::PeriodicExtrapolationRule> const + spline_evaluator(periodic_extrapolation, periodic_extrapolation); ddc::Chunk coords_eval_alloc(interpolation_domain, ddc::KokkosAllocator()); @@ -140,8 +139,9 @@ void TestPeriodicSplineBuilderTestIdentity() spline_evaluator .deriv(spline_eval_deriv.span_view(), coords_eval.span_cview(), coef.span_cview()); - ddc::Chunk - integral(spline_builder.batch_domain(), ddc::KokkosAllocator()); + ddc::Chunk integral( + spline_builder.batch_domain(interpolation_domain), + ddc::KokkosAllocator()); spline_evaluator.integrate(integral.span_view(), coef.span_cview()); ddc::Chunk, ddc::KokkosAllocator> diff --git a/tests/splines/periodicity_spline_builder.cpp b/tests/splines/periodicity_spline_builder.cpp index 8cf073bf3..17c2a8574 100644 --- a/tests/splines/periodicity_spline_builder.cpp +++ b/tests/splines/periodicity_spline_builder.cpp @@ -131,8 +131,7 @@ void PeriodicitySplineBuilderTest() DDim, ddc::BoundCond::PERIODIC, ddc::BoundCond::PERIODIC, - ddc::SplineSolver::GINKGO, - DDim> const spline_builder(dom_vals); + ddc::SplineSolver::GINKGO> const spline_builder(dom_vals); // Compute useful domains (dom_interpolation, dom_batch, dom_bsplines and dom_spline) ddc::DiscreteDomain> const dom_bsplines = spline_builder.spline_domain(); @@ -160,8 +159,8 @@ void PeriodicitySplineBuilderTest() BSplines, DDim, ddc::PeriodicExtrapolationRule, - ddc::PeriodicExtrapolationRule, - DDim> const spline_evaluator(extrapolation_rule, extrapolation_rule); + ddc::PeriodicExtrapolationRule> const + spline_evaluator(extrapolation_rule, extrapolation_rule); // Instantiate chunk of coordinates of dom_interpolation ddc::Chunk coords_eval_alloc(dom_vals, ddc::KokkosAllocator, MemorySpace>()); diff --git a/tests/splines/spline_builder.cpp b/tests/splines/spline_builder.cpp index cc0d114ae..3c8df451a 100644 --- a/tests/splines/spline_builder.cpp +++ b/tests/splines/spline_builder.cpp @@ -55,8 +55,7 @@ TEST(SplineBuilder, ShortInterpolationGrid) DDimX, ddc::BoundCond::PERIODIC, ddc::BoundCond::PERIODIC, - ddc::SplineSolver::GINKGO, - DDimX>(interpolation_domain)), + ddc::SplineSolver::GINKGO>(interpolation_domain)), std::runtime_error); } @@ -82,8 +81,7 @@ TEST(SplineBuilder, LongInterpolationGrid) DDimX, ddc::BoundCond::PERIODIC, ddc::BoundCond::PERIODIC, - ddc::SplineSolver::GINKGO, - DDimX>(interpolation_domain)), + ddc::SplineSolver::GINKGO>(interpolation_domain)), std::runtime_error); } @@ -109,8 +107,7 @@ TEST(SplineBuilder, BadShapeInterpolationGrid) DDimX, ddc::BoundCond::PERIODIC, ddc::BoundCond::PERIODIC, - ddc::SplineSolver::GINKGO, - DDimX>(interpolation_domain)), + ddc::SplineSolver::GINKGO>(interpolation_domain)), std::runtime_error); } @@ -134,6 +131,5 @@ TEST(SplineBuilder, CorrectInterpolationGrid) DDimX, ddc::BoundCond::PERIODIC, ddc::BoundCond::PERIODIC, - ddc::SplineSolver::GINKGO, - DDimX>(interpolation_domain))); + ddc::SplineSolver::GINKGO>(interpolation_domain))); }