diff --git a/benchmarks/splines.cpp b/benchmarks/splines.cpp index 8cee4dd1f..7d16b4386 100644 --- a/benchmarks/splines.cpp +++ b/benchmarks/splines.cpp @@ -178,6 +178,14 @@ void characteristics_advection_unitary(benchmark::State& state) spline_builder.batched_interpolation_domain(x_mesh), ddc::KokkosAllocator, typename ExecSpace::memory_space>()); ddc::ChunkSpan const feet_coords = feet_coords_alloc.span_view(); + ddc::ChunkSpan< + double const, + ddc::StridedDiscreteDomain< + DDimX, + ddc::Deriv::continuous_dimension_type>, + DDimY>, + Kokkos::layout_right, + typename ExecSpace::memory_space> const derivs {}; for (auto _ : state) { Kokkos::Profiling::pushRegion("FeetCharacteristics"); @@ -193,7 +201,7 @@ void characteristics_advection_unitary(benchmark::State& state) }); Kokkos::Profiling::popRegion(); Kokkos::Profiling::pushRegion("SplineBuilder"); - spline_builder(coef, density.span_cview()); + spline_builder(coef, density.span_cview(), derivs.span_cview()); Kokkos::Profiling::popRegion(); Kokkos::Profiling::pushRegion("SplineEvaluator"); spline_evaluator(density, feet_coords.span_cview(), coef.span_cview()); diff --git a/examples/characteristics_advection.cpp b/examples/characteristics_advection.cpp index e97c86604..eb0be4831 100644 --- a/examples/characteristics_advection.cpp +++ b/examples/characteristics_advection.cpp @@ -13,6 +13,7 @@ #include #include +#include #define PERIODIC_DOMAIN // Comment this to run non-periodic simulation @@ -230,6 +231,13 @@ int main(int argc, char** argv) spline_builder.batched_interpolation_domain(x_mesh), ddc::DeviceAllocator>()); ddc::ChunkSpan const feet_coords = feet_coords_alloc.span_view(); + + // Instantiate empty derivative chunkspan + ddc::ChunkSpan< + double const, + ddc::StridedDiscreteDomain, DDimY>, + Kokkos::layout_right, + ddc::DeviceAllocator::memory_space> const derivs {}; //! [instantiate intermediate chunks] @@ -255,7 +263,7 @@ int main(int argc, char** argv) - ddc::Coordinate(vx * ddc::step()); }); // Interpolate the values at feet on the grid - spline_builder(coef, last_density.span_cview()); + spline_builder(coef, last_density.span_cview(), derivs.span_cview()); spline_evaluator(next_density, feet_coords.span_cview(), coef.span_cview()); //! [numerical scheme] diff --git a/include/ddc/kernels/splines/spline_builder.hpp b/include/ddc/kernels/splines/spline_builder.hpp index b295e1baf..9822ef67b 100644 --- a/include/ddc/kernels/splines/spline_builder.hpp +++ b/include/ddc/kernels/splines/spline_builder.hpp @@ -17,6 +17,8 @@ #include +#include "ddc/detail/type_seq.hpp" + #include "deriv.hpp" #include "integrals.hpp" #include "math_tools.hpp" @@ -25,6 +27,145 @@ namespace ddc { +namespace detail { + +// TODO: move all these functions to another file + +/** + * @brief Create a DiscreteVector with each component initialized to the same value. + */ +template +ddc::DiscreteVector make_discrete_vector(T const& value) +{ + return ddc::DiscreteVector((ddc::DiscreteVector(value))...); +} + +/** + * @brief Convert a DiscreteDomain to a stride 1 StridedDiscreteDomain. + */ +template +ddc::StridedDiscreteDomain to_strided_ddom(ddc::DiscreteDomain const& ddom) +{ + return ddc::StridedDiscreteDomain< + DDim...>(ddom.front(), ddom.extents(), make_discrete_vector(1)); +} + +template +struct to_whole_derivs_domain; + +template +struct to_whole_derivs_domain, TypeSeq, TypeSeq> +{ + using type = convert_type_seq_to_strided_discrete_domain_t, + type_seq_replace_t, TypeSeq, TypeSeq>>>; +}; + +/** + * @brief The type of the deriv domain to be passed as argument to a spline builder. + * + * @tparam SeqDDimsI A TypeSeq containing the dimensions of interest. + * @tparam SeqDerivDims A TypeSeq containing the deriv dimensions. + * @tparam SeqDDims A TypeSeq containing all the dimensions (interpolation + batch). + */ +template +using to_whole_derivs_domain_t = + typename to_whole_derivs_domain::type; + +/** + * @brief Constructs a StridedDiscreteDomain to pass as argument to a spline builder. + * + * This function creates a StridedDiscreteDomain containing a strided dimension for the + * first and last point of each dimension of an interpolation domain, and the dimensions + * of a batched interpolation domain where the interpolation dimensions have been replaced + * by deriv dimensions. + * + * Example: For + * - interpolation_domain = DiscreteDomain(De(0,0), Dv(Nx,Nz)) + * - batched_domain = DiscreteDomain(De(0,0,0,0), Dv(Nx,Ny,Nz,Nt)) + * - DerivDims = [DX, DZ], + * this is StridedDiscreteDomain( + * De(0,0,1,0,1,0), + * Dv(2,2,nbc_x,Ny,nbc_z,Nt), + * Dv(Nx-2,Nz-2,1,1,1,1)) + * + * @tparam DerivDims... The deriv dimensions which will replace the interpolation dimensions. + * @param interpolation_dom The domain containing only the interpolation dimensions. + * @param batched_domain The domain containing both the interpolation and batch dimensions. + * @param nb_constraints The number of constraints for each of the deriv dimensions. + * + * @return A StridedDiscreteDomain with the batch dimensions, the deriv dimensions and the strided interpolation dimensions. + */ +template < + typename... DerivDims, + typename... Ints, + typename... DDimsI, + template class DDomInterp, + typename... DDims, + template class DDomBatched> +auto get_whole_derivs_domain( + DDomInterp const& interpolation_dom, + DDomBatched const& batched_domain, + Ints... nb_constraints) +{ + static_assert(sizeof...(DerivDims) == sizeof...(Ints)); + + using result_domain_t = to_whole_derivs_domain_t< + TypeSeq, + TypeSeq, + TypeSeq>; + using DElem = typename result_domain_t::discrete_element_type; + using DVect = typename result_domain_t::discrete_vector_type; + + auto const batch_domain = ddc::remove_dims_of(batched_domain); + + return result_domain_t( + DElem(interpolation_dom.front(), + batch_domain.front(), + (ddc::DiscreteElement(1))...), + DVect((ddc::DiscreteVector(2))..., + batch_domain.extents(), + (ddc::DiscreteVector(nb_constraints))...), + DVect((ddc::DiscreteVector( + interpolation_dom.template extent().value() - 1))..., + typename decltype(batch_domain)::discrete_vector_type( + (ddc::DiscreteVector(1))...), + (ddc::DiscreteVector(1))...)); +} + +/** + * @brief Constructs a StridedDiscreteDomain to pass as argument to a spline builder. + * + * This function creates a StridedDiscreteDomain containing a strided dimension for the + * first and last point of each dimension of an interpolation domain, and the derivs dimensions. + * + * Example: For + * - interpolation_domain = DiscreteDomain(De(0,0), Dv(Nx,Nz)) + * - DerivDims = [DX,DY], + * this is StridedDiscreteDomain( + * De(0,0,1,1), + * Dv(2,2,nbc_x,nbc_y), + * Dv(Nx-1,Ny-1,1,1)) + * + * @tparam DerivDims... The deriv dimensions. + * @param interpolation_dom The domain containing the interpolation dimensions. + * @param nb_constraints The number of constraints for each of the deriv dimensions. + * + * @return A StridedDiscreteDomain with the deriv dimensions and the strided interpolation dimensions. + */ +template < + typename... DerivDims, + typename... Ints, + typename... DDimsI, + template class DDomInterp> +auto get_whole_derivs_domain(DDomInterp const& interpolation_dom, Ints... nb_constraints) +{ + return get_whole_derivs_domain< + DerivDims...>(interpolation_dom, interpolation_dom, nb_constraints...); +} + +} // namespace detail + /** * @brief An enum determining the backend solver of a SplineBuilder or SplineBuilder2d. * @@ -152,6 +293,27 @@ class SplineBuilder ddc::detail::TypeSeq>>>; public: + /** + * @brief The type of the whole Deriv domain (1D dimension of interest and cartesian + * product of 1D Deriv domain and batch domain) to be passed as argument to the builder, + * 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 StridedDiscreteDomain,Z> + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using whole_derivs_domain_type + = ddc::detail::convert_type_seq_to_strided_discrete_domain_t, + ddc::to_type_seq_t>>>; + /** * @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). @@ -159,15 +321,14 @@ class SplineBuilder * @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> + * this is StridedDiscreteDomain,Z> */ template < class BatchedInterpolationDDom, class = std::enable_if_t>> - using batched_derivs_domain_type = ddc::replace_dim_of_t< - BatchedInterpolationDDom, - interpolation_discrete_dimension_type, - deriv_type>; + using batched_derivs_domain_type = ddc::remove_dims_of_t< + whole_derivs_domain_type, + interpolation_discrete_dimension_type>; /// @brief Indicates if the degree of the splines is odd or even. static constexpr bool s_odd = BSplines::degree() % 2; @@ -421,11 +582,13 @@ class SplineBuilder BatchedInterpolationDDom const& batched_interpolation_domain) const noexcept { assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain)); - return ddc::replace_dim_of( - batched_interpolation_domain, - ddc::DiscreteDomain( - ddc::DiscreteElement(1), - ddc::DiscreteVector(s_nbc_xmin))); + return ddc::StridedDiscreteDomain( + detail::to_strided_ddom( + ddc::replace_dim_of( + batched_interpolation_domain, + ddc::DiscreteDomain( + ddc::DiscreteElement(1), + ddc::DiscreteVector(s_nbc_xmin))))); } /** @@ -442,11 +605,13 @@ class SplineBuilder BatchedInterpolationDDom const& batched_interpolation_domain) const noexcept { assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain)); - return ddc::replace_dim_of( - batched_interpolation_domain, - ddc::DiscreteDomain( - ddc::DiscreteElement(1), - ddc::DiscreteVector(s_nbc_xmax))); + return ddc::StridedDiscreteDomain( + detail::to_strided_ddom( + ddc::replace_dim_of( + batched_interpolation_domain, + ddc::DiscreteDomain( + ddc::DiscreteElement(1), + ddc::DiscreteVector(s_nbc_xmax))))); } /** @@ -463,12 +628,10 @@ class SplineBuilder * * @param[out] spline The coefficients of the spline computed by this SplineBuilder. * @param[in] vals The values of the function on the interpolation mesh. - * @param[in] derivs_xmin The values of the derivatives at the lower boundary - * (used only with BoundCond::HERMITE lower boundary condition). - * @param[in] derivs_xmax The values of the derivatives at the upper boundary + * @param[in] derivs The values of the derivatives. * (used only with BoundCond::HERMITE upper boundary condition). */ - template + template void operator()( ddc::ChunkSpan< double, @@ -476,18 +639,11 @@ class SplineBuilder Layout, memory_space> spline, ddc::ChunkSpan vals, - std::optional, - Layout, - memory_space>> derivs_xmin - = std::nullopt, - std::optional, - Layout, - memory_space>> derivs_xmax - = std::nullopt) const; + whole_derivs_domain_type, + StridedLayout, + memory_space> derivs) const; /** * @brief Compute the quadrature coefficients associated to the b-splines used by this SplineBuilder. @@ -793,7 +949,7 @@ template < ddc::BoundCond BcLower, ddc::BoundCond BcUpper, SplineSolver Solver> -template +template void SplineBuilder:: operator()( ddc::ChunkSpan< @@ -802,16 +958,11 @@ operator()( Layout, memory_space> spline, ddc::ChunkSpan vals, - std::optional, - Layout, - memory_space>> const derivs_xmin, - std::optional, - Layout, - memory_space>> const derivs_xmax) const + whole_derivs_domain_type, + StridedLayout, + memory_space> const derivs) const { auto const batched_interpolation_domain = vals.domain(); @@ -820,23 +971,17 @@ operator()( assert(vals.template extent() == ddc::discrete_space().nbasis() - s_nbc_xmin - s_nbc_xmax); - assert((BcLower == ddc::BoundCond::HERMITE) - != (!derivs_xmin.has_value() || derivs_xmin->template extent() == 0)); - assert((BcUpper == ddc::BoundCond::HERMITE) - != (!derivs_xmax.has_value() || derivs_xmax->template extent() == 0)); - if constexpr (BcLower == BoundCond::HERMITE) { - assert(ddc::DiscreteElement(derivs_xmin->domain().front()).uid() == 1); - } - if constexpr (BcUpper == BoundCond::HERMITE) { - assert(ddc::DiscreteElement(derivs_xmax->domain().front()).uid() == 1); + if constexpr (BcLower == BoundCond::HERMITE || BcUpper == BoundCond::HERMITE) { + assert(derivs.template extent() > 0); + assert(ddc::DiscreteElement(derivs.domain().front()).uid() == 1); } // Hermite boundary conditions at xmin, if any // NOTE: For consistency with the linear system, the i-th derivative // provided by the user must be multiplied by dx^i if constexpr (BcLower == BoundCond::HERMITE) { - assert(derivs_xmin->template extent() == s_nbc_xmin); - auto derivs_xmin_values = *derivs_xmin; + assert(derivs.template extent() == s_nbc_xmin); + auto derivs_xmin_values = derivs[interpolation_domain().front()]; auto const dx_proxy = m_dx; ddc::parallel_for_each( "ddc_splines_hermite_compute_lower_coefficients", @@ -880,8 +1025,8 @@ operator()( // provided by the user must be multiplied by dx^i auto const& nbasis_proxy = ddc::discrete_space().nbasis(); if constexpr (BcUpper == BoundCond::HERMITE) { - assert(derivs_xmax->template extent() == s_nbc_xmax); - auto derivs_xmax_values = *derivs_xmax; + assert(derivs.template extent() == s_nbc_xmax); + auto derivs_xmax_values = derivs[interpolation_domain().back()]; auto const dx_proxy = m_dx; ddc::parallel_for_each( "ddc_splines_hermite_compute_upper_coefficients", diff --git a/include/ddc/kernels/splines/spline_builder_2d.hpp b/include/ddc/kernels/splines/spline_builder_2d.hpp index 03dc45e10..75963268d 100644 --- a/include/ddc/kernels/splines/spline_builder_2d.hpp +++ b/include/ddc/kernels/splines/spline_builder_2d.hpp @@ -15,6 +15,111 @@ namespace ddc { +namespace detail { +template +auto strided_to_discrete_domain_chunkspan( + ddc::ChunkSpan< + ElementType, + ddc::StridedDiscreteDomain, + Layout, + MemorySpace> const& cs) +{ + assert(((cs.domain().strides().template get() == 1) && ...)); + return ddc::ChunkSpan, Layout, MemorySpace>( + cs.data_handle(), + ddc::DiscreteDomain(cs.domain().front(), cs.domain().extents())); +} + +template +auto discrete_to_strided_domain_chunkspan( + ddc::ChunkSpan, Layout, MemorySpace> const& cs) +{ + return ddc::ChunkSpan, Layout, MemorySpace>( + cs.data_handle(), + ddc::StridedDiscreteDomain( + cs.domain().front(), + cs.domain().extents(), + make_discrete_vector(1))); +} + +template +struct Min +{ +}; + +template +constexpr Min dmin; + +template +struct Max +{ +}; + +template +constexpr Max dmax; + +template +auto derivs( + ddc::ChunkSpan const& derivs_span, + ddc::DiscreteElement delem) +{ + return derivs_span[delem]; +} + +template < + typename DerivDim, + typename... DerivDims, + typename... DDims, + typename T, + typename Support, + class Layout, + class MemorySpace> +auto derivs( + ddc::ChunkSpan const& derivs_span, + ddc::DiscreteElement delem, + Min, + DerivDims... deriv_dims) +{ + return derivs( + derivs_span, + ddc::DiscreteElement< + DDims..., + DerivDim>(delem, ddc::DiscreteElement(derivs_span.domain().front())), + deriv_dims...); +} + +template < + typename DerivDim, + typename... DerivDims, + typename... DDims, + typename T, + typename Support, + class Layout, + class MemorySpace> +auto derivs( + ddc::ChunkSpan const& derivs_span, + ddc::DiscreteElement delem, + Max, + DerivDims... deriv_dims) +{ + return derivs( + derivs_span, + ddc::DiscreteElement< + DDims..., + DerivDim>(delem, ddc::DiscreteElement(derivs_span.domain().back())), + deriv_dims...); +} + +template +auto derivs( + ddc::ChunkSpan const& derivs_span, + DerivDims... deriv_dims) +{ + return derivs(derivs_span, ddc::DiscreteElement<>(), deriv_dims...); +} + +} // namespace detail + /** * @brief A class for creating a 2D spline approximation of a function. * @@ -157,6 +262,22 @@ class SplineBuilder2D using batched_derivs_domain_type1 = typename builder_type1::template batched_derivs_domain_type; + /** + * @brief The type of the whole Derivs domain (1D dimension of interest and cartesian + * product of 1D Deriv domain and batch domain) in the first dimension to be passed as argument to the builder, + * 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 dimensions of interest X and Y, + * this is StridedDiscreteDomain, Y, Z> + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using whole_derivs_domain_type1 = + typename builder_type1::template whole_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. @@ -174,6 +295,22 @@ class SplineBuilder2D interpolation_discrete_dimension_type2, deriv_type2>; + /** + * @brief The type of the whole Derivs domain (1D dimension of interest and cartesian + * product of 1D Deriv domain and batch domain) in the second dimension to be passed as argument to the builder, + * 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 dimensions of interest X and Y, + * this is StridedDiscreteDomain, Z> + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using whole_derivs_domain_type2 = + typename builder_type2::template whole_derivs_domain_type; + /** * @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. @@ -194,6 +331,29 @@ class SplineBuilder2D interpolation_discrete_dimension_type2>, ddc::detail::TypeSeq>>; + /** + * @brief The type of the whole Derivs domain (all dimensions of interest and cartesian + * product of 2D Deriv domain and batch domain) to be passed as argument to the builder, + * 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 dimensions of interest X and Y, + * this is StridedDiscreteDomain, Deriv, Z> + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using whole_derivs_domain_type + = ddc::detail::convert_type_seq_to_strided_discrete_domain_t, + ddc::type_seq_replace_t< + ddc::to_type_seq_t, + ddc::detail::TypeSeq< + interpolation_discrete_dimension_type1, + interpolation_discrete_dimension_type2>, + ddc::detail::TypeSeq>>>; + private: builder_type1 m_spline_builder1; builder_deriv_type1 m_spline_builder_deriv1; @@ -384,28 +544,14 @@ class SplineBuilder2D * The coefficients of the spline computed by this SplineBuilder. * @param[in] vals * The values of the function at the interpolation mesh. - * @param[in] derivs_min1 - * The values of the derivatives at the lower boundary in the first dimension. - * @param[in] derivs_max1 - * The values of the derivatives at the upper boundary in the first dimension. - * @param[in] derivs_min2 - * The values of the derivatives at the lower boundary in the second dimension. - * @param[in] derivs_max2 - * The values of the derivatives at the upper boundary in the second dimension. - * @param[in] mixed_derivs_min1_min2 - * The values of the the cross-derivatives at the lower boundary in the first dimension - * and the lower boundary in the second dimension. - * @param[in] mixed_derivs_max1_min2 - * The values of the the cross-derivatives at the upper boundary in the first dimension - * and the lower boundary in the second dimension. - * @param[in] mixed_derivs_min1_max2 - * The values of the the cross-derivatives at the lower boundary in the first dimension - * and the upper boundary in the second dimension. - * @param[in] mixed_derivs_max1_max2 - * The values of the the cross-derivatives at the upper boundary in the first dimension - * and the upper boundary in the second dimension. + * @param[in] derivs1 + * The values of the derivatives in the first dimension. + * @param[in] derivs2 + * The values of the derivatives in the second dimension. + * @param[in] mixed_derivs_1_2 + * The values of the the cross-derivatives in the first and the second dimension. */ - template + template void operator()( ddc::ChunkSpan< double, @@ -413,54 +559,21 @@ class SplineBuilder2D Layout, memory_space> spline, ddc::ChunkSpan vals, - std::optional, - Layout, - memory_space>> derivs_min1 - = std::nullopt, - std::optional, - Layout, - memory_space>> derivs_max1 - = std::nullopt, - std::optional, - Layout, - memory_space>> derivs_min2 - = std::nullopt, - std::optional, - Layout, - memory_space>> derivs_max2 - = std::nullopt, - std::optional, - Layout, - memory_space>> mixed_derivs_min1_min2 - = std::nullopt, - std::optional, - Layout, - memory_space>> mixed_derivs_max1_min2 - = std::nullopt, - std::optional, + DerivsLayout, + memory_space> derivs1, + ddc::ChunkSpan< double const, - batched_derivs_domain_type, - Layout, - memory_space>> mixed_derivs_min1_max2 - = std::nullopt, - std::optional, + DerivsLayout, + memory_space> derivs2, + ddc::ChunkSpan< double const, - batched_derivs_domain_type, - Layout, - memory_space>> mixed_derivs_max1_max2 - = std::nullopt) const; + whole_derivs_domain_type, + DerivsLayout, + memory_space> mixed_derivs_1_2) const; }; @@ -476,7 +589,7 @@ template < ddc::BoundCond BcLower2, ddc::BoundCond BcUpper2, ddc::SplineSolver Solver> -template +template void SplineBuilder2D< ExecSpace, MemorySpace, @@ -496,74 +609,53 @@ operator()( 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, + DerivsLayout, + memory_space> derivs1, + ddc::ChunkSpan< double const, - batched_derivs_domain_type, - Layout, - memory_space>> const mixed_derivs_min1_max2, - std::optional, + DerivsLayout, + memory_space> derivs2, + ddc::ChunkSpan< double const, - batched_derivs_domain_type, - Layout, - memory_space>> const mixed_derivs_max1_max2) const + whole_derivs_domain_type, + DerivsLayout, + memory_space> mixed_derivs_1_2) 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))); + // Spline1-approximate derivs_min2 (to spline1_deriv_min) + auto const spline_batched_deriv_domain = ddc::detail::get_whole_derivs_domain( + ddc::select(batched_interpolation_domain), + m_spline_builder_deriv1.batched_spline_domain(batched_interpolation_domain), + bsplines_type2::degree() / 2); - ddc::Chunk spline1_deriv_min_alloc( - m_spline_builder_deriv1.batched_spline_domain(batched_interpolation_deriv_domain), + ddc::Chunk spline1_deriv_alloc( + spline_batched_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()); + auto spline1_deriv = spline1_deriv_alloc.span_view(); + if constexpr (BcLower2 == ddc::BoundCond::HERMITE) { + constexpr auto min_deriv_2 = ddc::detail::dmin; + + auto const spline1_deriv_min_strided = detail::derivs(spline1_deriv, min_deriv_2); + auto const spline1_deriv_min + = detail::strided_to_discrete_domain_chunkspan(spline1_deriv_min_strided); + + auto const derivs_min2_strided = detail::derivs(derivs2, min_deriv_2); + auto const derivs_min2 = detail::strided_to_discrete_domain_chunkspan(derivs_min2_strided); + m_spline_builder_deriv1( spline1_deriv_min, - *derivs_min2, - mixed_derivs_min1_min2, - mixed_derivs_max1_min2); - } else { - spline1_deriv_min_opt = std::nullopt; + derivs_min2, + detail::derivs(mixed_derivs_1_2, min_deriv_2)); } // Spline1-approximate vals (to spline1) @@ -572,26 +664,27 @@ operator()( ddc::KokkosAllocator()); ddc::ChunkSpan const spline1 = spline1_alloc.span_view(); - m_spline_builder1(spline1, vals, derivs_min1, derivs_max1); + m_spline_builder1(spline1, vals, derivs1); // Spline1-approximate derivs_max2 (to spline1_deriv_max) - ddc::Chunk spline1_deriv_max_alloc( - 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()); if constexpr (BcUpper2 == ddc::BoundCond::HERMITE) { + constexpr auto max_deriv_2 = ddc::detail::dmax; + + auto const spline1_deriv_max_strided = detail::derivs(spline1_deriv, max_deriv_2); + auto const spline1_deriv_max + = detail::strided_to_discrete_domain_chunkspan(spline1_deriv_max_strided); + + auto const derivs_max2_strided = detail::derivs(derivs2, max_deriv_2); + auto const derivs_max2 = detail::strided_to_discrete_domain_chunkspan(derivs_max2_strided); + m_spline_builder_deriv1( spline1_deriv_max, - *derivs_max2, - mixed_derivs_min1_max2, - mixed_derivs_max1_max2); - } else { - spline1_deriv_max_opt = std::nullopt; + derivs_max2, + detail::derivs(mixed_derivs_1_2, max_deriv_2)); } // Spline2-approximate spline1 - m_spline_builder2(spline, spline1.span_cview(), spline1_deriv_min_opt, spline1_deriv_max_opt); + m_spline_builder2(spline, spline1.span_cview(), spline1_deriv.span_cview()); } } // namespace ddc diff --git a/include/ddc/kernels/splines/spline_builder_3d.hpp b/include/ddc/kernels/splines/spline_builder_3d.hpp index e20171457..cde6c2840 100644 --- a/include/ddc/kernels/splines/spline_builder_3d.hpp +++ b/include/ddc/kernels/splines/spline_builder_3d.hpp @@ -190,6 +190,24 @@ class SplineBuilder3D interpolation_discrete_dimension_type1, deriv_type1>; + /** + * @brief The type of the whole Derivs domain (1D dimension of interest and cartesian + * product of 1D Deriv domain and batch domain) in the first dimension, to be passed as + * argument to the builder, 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 dimensions of interest X, Y and Z, + * this is StridedDiscreteDomain, Y, Z, T> + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using whole_derivs_domain_type1 = detail::to_whole_derivs_domain_t< + ddc::detail::TypeSeq, + ddc::detail::TypeSeq, + ddc::to_type_seq_t>; + /** * @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. @@ -207,6 +225,24 @@ class SplineBuilder3D interpolation_discrete_dimension_type2, deriv_type2>; + /** + * @brief The type of the whole Derivs domain (1D dimension of interest and cartesian + * product of 1D Deriv domain and batch domain) in the second dimension, to be passed as + * argument to the builder, 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 dimensions of interest X, Y and Z, + * this is StridedDiscreteDomain, Z, T> + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using whole_derivs_domain_type2 = detail::to_whole_derivs_domain_t< + ddc::detail::TypeSeq, + ddc::detail::TypeSeq, + ddc::to_type_seq_t>; + /** * @brief The type of the whole Derivs domain (cartesian product of the 1D Deriv domain * and the associated batch domain) in the third dimension, preserving the order of dimensions. @@ -224,6 +260,24 @@ class SplineBuilder3D interpolation_discrete_dimension_type3, deriv_type3>; + /** + * @brief The type of the whole Derivs domain (1D dimension of interest and cartesian + * product of 1D Deriv domain and batch domain) in the third dimension, to be passed as + * argument to the builder, 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 dimensions of interest X, Y and Z, + * this is StridedDiscreteDomain, T> + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using whole_derivs_domain_type3 = detail::to_whole_derivs_domain_t< + ddc::detail::TypeSeq, + ddc::detail::TypeSeq, + ddc::to_type_seq_t>; + /** * @brief The type of the whole Derivs domain (cartesian product of the 2D Deriv domain * and the associated batch domain) in the first and second dimensions, preserving the order @@ -245,6 +299,26 @@ class SplineBuilder3D interpolation_domain_type2, deriv_type2>; + /** + * @brief The type of the whole Derivs domain (cartesian product of the 2D Deriv domain + * and the associated batch domain) in the first and second dimensions, to be passed + * as argument to the builder, preserving the order of dimensions. + * + * @tparam The batched discrete domain on which the interpolation points are defined. + * + * Example: For batched_interpolation_domain_type = DiscreteDomain and dimensions of interest X, Y and Z, + * this is StridedDiscreteDomain, Deriv, Z, T> + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using whole_derivs_domain_type1_2 = detail::to_whole_derivs_domain_t< + ddc::detail::TypeSeq< + interpolation_discrete_dimension_type1, + interpolation_discrete_dimension_type2>, + ddc::detail::TypeSeq, + ddc::to_type_seq_t>; + /** * @brief The type of the whole Derivs domain (cartesian product of the 2D Deriv domain * and the associated batch domain) in the second and third dimensions, preserving the order @@ -266,6 +340,26 @@ class SplineBuilder3D interpolation_domain_type3, deriv_type3>; + /** + * @brief The type of the whole Derivs domain (cartesian product of the 2D Deriv domain + * and the associated batch domain) in the second and third dimensions, to be passed + * as argument to the builder, preserving the order of dimensions. + * + * @tparam The batched discrete domain on which the interpolation points are defined. + * + * Example: For batched_interpolation_domain_type = DiscreteDomain and dimensions of interest X, Y and Z, + * this is StridedDiscreteDomain, Deriv, T> + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using whole_derivs_domain_type2_3 = detail::to_whole_derivs_domain_t< + ddc::detail::TypeSeq< + interpolation_discrete_dimension_type2, + interpolation_discrete_dimension_type3>, + ddc::detail::TypeSeq, + ddc::to_type_seq_t>; + /** * @brief The type of the whole Derivs domain (cartesian product of the 2D Deriv domain * and the associated batch domain) in the first and third dimensions, preserving the order @@ -287,6 +381,26 @@ class SplineBuilder3D interpolation_domain_type3, deriv_type3>; + /** + * @brief The type of the whole Derivs domain (cartesian product of the 2D Deriv domain + * and the associated batch domain) in the first and third dimensions, to be passed + * as argument to the builder, preserving the order of dimensions. + * + * @tparam The batched discrete domain on which the interpolation points are defined. + * + * Example: For batched_interpolation_domain_type = DiscreteDomain and dimensions of interest X, Y and Z, + * this is StridedDiscreteDomain, Y, Deriv, T> + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using whole_derivs_domain_type1_3 = detail::to_whole_derivs_domain_t< + ddc::detail::TypeSeq< + interpolation_discrete_dimension_type1, + interpolation_discrete_dimension_type3>, + ddc::detail::TypeSeq, + ddc::to_type_seq_t>; + /** * @brief The type of the whole Derivs domain (cartesian product of the 3D Deriv domain * and the batch domain), preserving the order of dimensions. @@ -308,6 +422,27 @@ class SplineBuilder3D interpolation_discrete_dimension_type3>, ddc::detail::TypeSeq>>; + /** + * @brief The type of the whole Derivs domain (cartesian product of the 3D Deriv domain + * and the associated batch domain) to be passed as argument to the builder, preserving + * the order of dimensions. + * + * @tparam The batched discrete domain on which the interpolation points are defined. + * + * Example: For batched_interpolation_domain_type = DiscreteDomain and dimensions of interest X, Y and Z, + * this is StridedDiscreteDomain, Deriv, Deriv, T> + */ + template < + class BatchedInterpolationDDom, + class = std::enable_if_t>> + using whole_derivs_domain_type = detail::to_whole_derivs_domain_t< + ddc::detail::TypeSeq< + interpolation_discrete_dimension_type1, + interpolation_discrete_dimension_type2, + interpolation_discrete_dimension_type3>, + ddc::detail::TypeSeq, + ddc::to_type_seq_t>; + private: builder_type1 m_spline_builder1; builder_type_2_3 m_spline_builder_2_3; @@ -503,60 +638,22 @@ class SplineBuilder3D * The coefficients of the spline computed by this SplineBuilder. * @param[in] vals * The values of the function at the interpolation mesh. - * @param[in] derivs_min1 - * The values of the derivatives at the lower boundary in the first dimension. - * @param[in] derivs_max1 - * The values of the derivatives at the upper boundary in the first dimension. - * @param[in] derivs_min2 - * The values of the derivatives at the lower boundary in the second dimension. - * @param[in] derivs_max2 - * The values of the derivatives at the upper boundary in the second dimension. - * @param[in] derivs_min3 - * The values of the derivatives at the lower boundary in the third dimension. - * @param[in] derivs_max3 - * The values of the derivatives at the upper boundary in the third dimension. - * @param[in] mixed_derivs_min1_min2 - * The values of the the cross-derivatives at the lower boundary in the first dimension and the lower boundary in the second dimension. - * @param[in] mixed_derivs_max1_min2 - * The values of the the cross-derivatives at the upper boundary in the first dimension and the lower boundary in the second dimension. - * @param[in] mixed_derivs_min1_max2 - * The values of the the cross-derivatives at the lower boundary in the first dimension and the upper boundary in the second dimension. - * @param[in] mixed_derivs_max1_max2 - * The values of the the cross-derivatives at the upper boundary in the first dimension and the upper boundary in the second dimension. - * @param[in] mixed_derivs_min2_min3 - * The values of the the cross-derivatives at the lower boundary in the second dimension and the lower boundary in the third dimension. - * @param[in] mixed_derivs_max2_min3 - * The values of the the cross-derivatives at the upper boundary in the second dimension and the lower boundary in the third dimension. - * @param[in] mixed_derivs_min2_max3 - * The values of the the cross-derivatives at the lower boundary in the second dimension and the upper boundary in the third dimension. - * @param[in] mixed_derivs_max2_max3 - * The values of the the cross-derivatives at the upper boundary in the second dimension and the upper boundary in the third dimension. - * @param[in] mixed_derivs_min1_min3 - * The values of the the cross-derivatives at the lower boundary in the first dimension and the lower boundary in the third dimension. - * @param[in] mixed_derivs_max1_min3 - * The values of the the cross-derivatives at the upper boundary in the first dimension and the lower boundary in the third dimension. - * @param[in] mixed_derivs_min1_max3 - * The values of the the cross-derivatives at the lower boundary in the first dimension and the upper boundary in the third dimension. - * @param[in] mixed_derivs_max1_max3 - * The values of the the cross-derivatives at the upper boundary in the first dimension and the upper boundary in the third dimension. - * @param[in] mixed_derivs_min1_min2_min3 - * The values of the the cross-derivatives at the lower boundary in the first dimension, the lower boundary in the second dimension and the lower boundary in the third dimension. - * @param[in] mixed_derivs_max1_min2_min3 - * The values of the the cross-derivatives at the upper boundary in the first dimension, the lower boundary in the second dimension and the lower boundary in the third dimension. - * @param[in] mixed_derivs_min1_max2_min3 - * The values of the the cross-derivatives at the lower boundary in the first dimension, the upper boundary in the second dimension and the lower boundary in the third dimension. - * @param[in] mixed_derivs_max1_max2_min3 - * The values of the the cross-derivatives at the upper boundary in the first dimension, the upper boundary in the second dimension and the lower boundary in the third dimension. - * @param[in] mixed_derivs_min1_min2_max3 - * The values of the the cross-derivatives at the lower boundary in the first dimension, the lower boundary in the second dimension and the upper boundary in the third dimension. - * @param[in] mixed_derivs_max1_min2_max3 - * The values of the the cross-derivatives at the upper boundary in the first dimension, the lower boundary in the second dimension and the upper boundary in the third dimension. - * @param[in] mixed_derivs_min1_max2_max3 - * The values of the the cross-derivatives at the lower boundary in the first dimension, the upper boundary in the second dimension and the upper boundary in the third dimension. - * @param[in] mixed_derivs_max1_max2_max3 - * The values of the the cross-derivatives at the upper boundary in the first dimension, the upper boundary in the second dimension and the upper boundary in the third dimension. + * @param[in] derivs1 + * The values of the derivatives in the first dimension. + * @param[in] derivs2 + * The values of the derivatives in the second dimension. + * @param[in] derivs3 + * The values of the derivatives in the third dimension. + * @param[in] mixed_derivs1_2 + * The values of the the cross-derivatives in the first dimension and in the second dimension. + * @param[in] mixed_derivs2_3 + * The values of the the cross-derivatives in the second dimension and in the third dimension. + * @param[in] mixed_derivs1_3 + * The values of the the cross-derivatives in the first dimension and in the third dimension. + * @param[in] mixed_derivs1_2_3 + * The values of the the cross-derivatives in the first dimension, in the second dimension and in the third dimension. */ - template + template void operator()( ddc::ChunkSpan< double, @@ -564,162 +661,41 @@ class SplineBuilder3D Layout, memory_space> spline, ddc::ChunkSpan vals, - std::optional, - Layout, - memory_space>> derivs_min1 - = std::nullopt, - std::optional, - Layout, - memory_space>> derivs_max1 - = std::nullopt, - std::optional, - Layout, - memory_space>> derivs_min2 - = std::nullopt, - std::optional, - Layout, - memory_space>> derivs_max2 - = std::nullopt, - std::optional, - Layout, - memory_space>> derivs_min3 - = std::nullopt, - std::optional, - Layout, - memory_space>> derivs_max3 - = std::nullopt, - std::optional, - Layout, - memory_space>> mixed_derivs_min1_min2 - = std::nullopt, - std::optional, - Layout, - memory_space>> mixed_derivs_max1_min2 - = std::nullopt, - std::optional, - Layout, - memory_space>> mixed_derivs_min1_max2 - = std::nullopt, - std::optional, - Layout, - memory_space>> mixed_derivs_max1_max2 - = std::nullopt, - std::optional, - Layout, - memory_space>> mixed_derivs_min2_min3 - = std::nullopt, - std::optional, - Layout, - memory_space>> mixed_derivs_max2_min3 - = std::nullopt, - std::optional, - Layout, - memory_space>> mixed_derivs_min2_max3 - = std::nullopt, - std::optional, - Layout, - memory_space>> mixed_derivs_max2_max3 - = std::nullopt, - std::optional, - Layout, - memory_space>> mixed_derivs_min1_min3 - = std::nullopt, - std::optional, - Layout, - memory_space>> mixed_derivs_max1_min3 - = std::nullopt, - std::optional, - Layout, - memory_space>> mixed_derivs_min1_max3 - = std::nullopt, - std::optional, - Layout, - memory_space>> mixed_derivs_max1_max3 - = std::nullopt, - std::optional, - Layout, - memory_space>> mixed_derivs_min1_min2_min3 - = std::nullopt, - std::optional, - Layout, - memory_space>> mixed_derivs_max1_min2_min3 - = std::nullopt, - std::optional, + DerivsLayout, + memory_space> derivs1, + ddc::ChunkSpan< double const, - batched_derivs_domain_type, - Layout, - memory_space>> mixed_derivs_min1_max2_min3 - = std::nullopt, - std::optional, + DerivsLayout, + memory_space> derivs2, + ddc::ChunkSpan< double const, - batched_derivs_domain_type, - Layout, - memory_space>> mixed_derivs_max1_max2_min3 - = std::nullopt, - std::optional, + DerivsLayout, + memory_space> derivs3, + ddc::ChunkSpan< double const, - batched_derivs_domain_type, - Layout, - memory_space>> mixed_derivs_min1_min2_max3 - = std::nullopt, - std::optional, + DerivsLayout, + memory_space> mixed_derivs1_2, + ddc::ChunkSpan< double const, - batched_derivs_domain_type, - Layout, - memory_space>> mixed_derivs_max1_min2_max3 - = std::nullopt, - std::optional, + DerivsLayout, + memory_space> mixed_derivs2_3, + ddc::ChunkSpan< double const, - batched_derivs_domain_type, - Layout, - memory_space>> mixed_derivs_min1_max2_max3 - = std::nullopt, - std::optional, + DerivsLayout, + memory_space> mixed_derivs1_3, + ddc::ChunkSpan< double const, - batched_derivs_domain_type, - Layout, - memory_space>> mixed_derivs_max1_max2_max3 - = std::nullopt) const; + whole_derivs_domain_type, + DerivsLayout, + memory_space> mixed_derivs1_2_3) const; }; @@ -739,7 +715,7 @@ template < ddc::BoundCond BcLower3, ddc::BoundCond BcUpper3, ddc::SplineSolver Solver> -template +template void SplineBuilder3D< ExecSpace, MemorySpace, @@ -763,283 +739,200 @@ operator()( Layout, memory_space> spline, ddc::ChunkSpan vals, - std::optional, - Layout, - memory_space>> derivs_min1, - std::optional, - Layout, - memory_space>> derivs_max1, - std::optional, - Layout, - memory_space>> derivs_min2, - std::optional, - Layout, - memory_space>> derivs_max2, - std::optional, - Layout, - memory_space>> derivs_min3, - std::optional, - Layout, - memory_space>> derivs_max3, - std::optional, - Layout, - memory_space>> mixed_derivs_min1_min2, - std::optional, - Layout, - memory_space>> mixed_derivs_max1_min2, - std::optional, - Layout, - memory_space>> mixed_derivs_min1_max2, - std::optional, - Layout, - memory_space>> mixed_derivs_max1_max2, - std::optional, - Layout, - memory_space>> mixed_derivs_min2_min3, - std::optional, - Layout, - memory_space>> mixed_derivs_max2_min3, - std::optional, - Layout, - memory_space>> mixed_derivs_min2_max3, - std::optional, - Layout, - memory_space>> mixed_derivs_max2_max3, - std::optional, - Layout, - memory_space>> mixed_derivs_min1_min3, - std::optional, - Layout, - memory_space>> mixed_derivs_max1_min3, - std::optional, - Layout, - memory_space>> mixed_derivs_min1_max3, - std::optional, - Layout, - memory_space>> mixed_derivs_max1_max3, - std::optional, - Layout, - memory_space>> mixed_derivs_min1_min2_min3, - std::optional, - Layout, - memory_space>> mixed_derivs_max1_min2_min3, - std::optional, + DerivsLayout, + memory_space> derivs1, + ddc::ChunkSpan< double const, - batched_derivs_domain_type, - Layout, - memory_space>> mixed_derivs_min1_max2_min3, - std::optional, + DerivsLayout, + memory_space> derivs2, + ddc::ChunkSpan< double const, - batched_derivs_domain_type, - Layout, - memory_space>> mixed_derivs_max1_max2_min3, - std::optional, + DerivsLayout, + memory_space> derivs3, + ddc::ChunkSpan< double const, - batched_derivs_domain_type, - Layout, - memory_space>> mixed_derivs_min1_min2_max3, - std::optional, + DerivsLayout, + memory_space> mixed_derivs1_2, + ddc::ChunkSpan< double const, - batched_derivs_domain_type, - Layout, - memory_space>> mixed_derivs_max1_min2_max3, - std::optional, + DerivsLayout, + memory_space> mixed_derivs2_3, + ddc::ChunkSpan< double const, - batched_derivs_domain_type, - Layout, - memory_space>> mixed_derivs_min1_max2_max3, - std::optional, + DerivsLayout, + memory_space> mixed_derivs1_3, + ddc::ChunkSpan< double const, - batched_derivs_domain_type, - Layout, - memory_space>> mixed_derivs_max1_max2_max3) const + whole_derivs_domain_type, + DerivsLayout, + memory_space> mixed_derivs1_2_3) const { auto const batched_interpolation_domain = vals.domain(); + using ddim2 = interpolation_discrete_dimension_type2; + using ddim3 = interpolation_discrete_dimension_type3; + using detail::dmax; + using detail::dmin; + assert(interpolation_domain() == interpolation_domain_type(batched_interpolation_domain)); // Build the derivatives along the second dimension - auto const batched_interpolation_deriv_domain2 - = ddc::replace_dim_of( - batched_interpolation_domain, - ddc::DiscreteDomain( - ddc::DiscreteElement(1), - ddc::DiscreteVector(bsplines_type2::degree() / 2))); - - ddc::Chunk spline_derivs_min2_alloc( - m_spline_builder1.batched_spline_domain(batched_interpolation_deriv_domain2), + auto const spline_batched_derivs2_domain = ddc::detail::get_whole_derivs_domain( + ddc::select(batched_interpolation_domain), + m_spline_builder1.batched_spline_domain(batched_interpolation_domain), + bsplines_type2::degree() / 2); + + ddc::Chunk spline_derivs2_alloc( + spline_batched_derivs2_domain, ddc::KokkosAllocator()); - auto spline_derivs_min2 = spline_derivs_min2_alloc.span_view(); - auto spline_derivs_min2_opt = std::optional(spline_derivs_min2.span_cview()); + auto spline_derivs2 = spline_derivs2_alloc.span_view(); + if constexpr (BcLower2 == ddc::BoundCond::HERMITE) { + auto const spline_derivs_min2_strided = detail::derivs(spline_derivs2, dmin); + auto const spline_derivs_min2 + = detail::strided_to_discrete_domain_chunkspan(spline_derivs_min2_strided); + + auto const derivs_min2_strided = detail::derivs(derivs2, dmin); + auto const derivs_min2 = detail::strided_to_discrete_domain_chunkspan(derivs_min2_strided); m_spline_builder1( spline_derivs_min2, - *derivs_min2, - mixed_derivs_min1_min2, - mixed_derivs_max1_min2); - } else { - spline_derivs_min2_opt = std::nullopt; + derivs_min2, + detail::derivs(mixed_derivs1_2, dmin)); } - ddc::Chunk spline_derivs_max2_alloc( - m_spline_builder1.batched_spline_domain(batched_interpolation_deriv_domain2), - ddc::KokkosAllocator()); - auto spline_derivs_max2 = spline_derivs_max2_alloc.span_view(); - auto spline_derivs_max2_opt = std::optional(spline_derivs_max2.span_cview()); if constexpr (BcUpper2 == ddc::BoundCond::HERMITE) { + auto const spline_derivs_max2_strided = detail::derivs(spline_derivs2, dmax); + auto const spline_derivs_max2 + = detail::strided_to_discrete_domain_chunkspan(spline_derivs_max2_strided); + + auto const derivs_max2_strided = detail::derivs(derivs2, dmax); + auto const derivs_max2 = detail::strided_to_discrete_domain_chunkspan(derivs_max2_strided); m_spline_builder1( spline_derivs_max2, - *derivs_max2, - mixed_derivs_min1_max2, - mixed_derivs_max1_max2); - } else { - spline_derivs_max2_opt = std::nullopt; + derivs_max2, + detail::derivs(mixed_derivs1_2, dmax)); } // Build the derivatives along the third dimension - auto const batched_interpolation_deriv_domain3 - = ddc::replace_dim_of( - batched_interpolation_domain, - ddc::DiscreteDomain( - ddc::DiscreteElement(1), - ddc::DiscreteVector(bsplines_type3::degree() / 2))); - - ddc::Chunk spline_derivs_min3_alloc( - m_spline_builder1.batched_spline_domain(batched_interpolation_deriv_domain3), + auto const spline_batched_derivs3_domain = ddc::detail::get_whole_derivs_domain( + ddc::select(batched_interpolation_domain), + m_spline_builder1.batched_spline_domain(batched_interpolation_domain), + bsplines_type3::degree() / 2); + + ddc::Chunk spline_derivs3_alloc( + spline_batched_derivs3_domain, ddc::KokkosAllocator()); - auto spline_derivs_min3 = spline_derivs_min3_alloc.span_view(); - auto spline_derivs_min3_opt = std::optional(spline_derivs_min3.span_cview()); + auto spline_derivs3 = spline_derivs3_alloc.span_view(); + if constexpr (BcLower3 == ddc::BoundCond::HERMITE) { + auto const spline_derivs_min3_strided = detail::derivs(spline_derivs3, dmin); + auto const spline_derivs_min3 + = detail::strided_to_discrete_domain_chunkspan(spline_derivs_min3_strided); + + auto const derivs_min3_strided = detail::derivs(derivs3, dmin); + auto const derivs_min3 = detail::strided_to_discrete_domain_chunkspan(derivs_min3_strided); m_spline_builder1( spline_derivs_min3, - *derivs_min3, - mixed_derivs_min1_min3, - mixed_derivs_max1_min3); - } else { - spline_derivs_min3_opt = std::nullopt; + derivs_min3, + detail::derivs(mixed_derivs1_3, dmin)); } - ddc::Chunk spline_derivs_max3_alloc( - m_spline_builder1.batched_spline_domain(batched_interpolation_deriv_domain3), - ddc::KokkosAllocator()); - auto spline_derivs_max3 = spline_derivs_max3_alloc.span_view(); - auto spline_derivs_max3_opt = std::optional(spline_derivs_max3.span_cview()); if constexpr (BcUpper3 == ddc::BoundCond::HERMITE) { + auto const spline_derivs_max3_strided = detail::derivs(spline_derivs3, dmax); + auto const spline_derivs_max3 + = detail::strided_to_discrete_domain_chunkspan(spline_derivs_max3_strided); + + auto const derivs_max3_strided = detail::derivs(derivs3, dmax); + auto const derivs_max3 = detail::strided_to_discrete_domain_chunkspan(derivs_max3_strided); m_spline_builder1( spline_derivs_max3, - *derivs_max3, - mixed_derivs_min1_max3, - mixed_derivs_max1_max3); - } else { - spline_derivs_max3_opt = std::nullopt; + derivs_max3, + detail::derivs(mixed_derivs1_3, dmax)); } // Build the cross derivatives along the second and third dimensions - auto batched_interpolation_deriv_domain2_3 - = ddc::replace_dim_of( - batched_interpolation_deriv_domain2, - ddc::DiscreteDomain( - ddc::DiscreteElement(1), - ddc::DiscreteVector(bsplines_type3::degree() / 2))); - - ddc::Chunk spline_derivs_min2_min3_alloc( - m_spline_builder1.batched_spline_domain(batched_interpolation_deriv_domain2_3), + auto const spline_batched_derivs2_3_domain + = ddc::detail::get_whole_derivs_domain( + ddc::select(batched_interpolation_domain), + m_spline_builder1.batched_spline_domain(batched_interpolation_domain), + bsplines_type2::degree() / 2, + bsplines_type3::degree() / 2); + + ddc::Chunk spline_derivs2_3_alloc( + spline_batched_derivs2_3_domain, ddc::KokkosAllocator()); - auto spline_derivs_min2_min3 = spline_derivs_min2_min3_alloc.span_view(); - auto spline_derivs_min2_min3_opt = std::optional(spline_derivs_min2_min3.span_cview()); + auto spline_derivs2_3 = spline_derivs2_3_alloc.span_view(); + if constexpr (BcLower2 == ddc::BoundCond::HERMITE || BcLower3 == ddc::BoundCond::HERMITE) { + auto const spline_derivs_min2_min3_strided + = detail::derivs(spline_derivs2_3, dmin, dmin); + auto const spline_derivs_min2_min3 + = detail::strided_to_discrete_domain_chunkspan(spline_derivs_min2_min3_strided); + + auto const derivs_min2_min3_strided + = detail::derivs(mixed_derivs2_3, dmin, dmin); + auto const derivs_min2_min3 + = detail::strided_to_discrete_domain_chunkspan(derivs_min2_min3_strided); m_spline_builder1( spline_derivs_min2_min3, - *mixed_derivs_min2_min3, - mixed_derivs_min1_min2_min3, - mixed_derivs_max1_min2_min3); - } else { - spline_derivs_min2_min3_opt = std::nullopt; + derivs_min2_min3, + detail::derivs(mixed_derivs1_2_3, dmin, dmin)); } - ddc::Chunk spline_derivs_min2_max3_alloc( - m_spline_builder1.batched_spline_domain(batched_interpolation_deriv_domain2_3), - ddc::KokkosAllocator()); - auto spline_derivs_min2_max3 = spline_derivs_min2_max3_alloc.span_view(); - auto spline_derivs_min2_max3_opt = std::optional(spline_derivs_min2_max3.span_cview()); if constexpr (BcLower2 == ddc::BoundCond::HERMITE || BcUpper3 == ddc::BoundCond::HERMITE) { + auto const spline_derivs_min2_max3_strided + = detail::derivs(spline_derivs2_3, dmin, dmax); + auto const spline_derivs_min2_max3 + = detail::strided_to_discrete_domain_chunkspan(spline_derivs_min2_max3_strided); + + auto const derivs_min2_max3_strided + = detail::derivs(mixed_derivs2_3, dmin, dmax); + auto const derivs_min2_max3 + = detail::strided_to_discrete_domain_chunkspan(derivs_min2_max3_strided); m_spline_builder1( spline_derivs_min2_max3, - *mixed_derivs_min2_max3, - mixed_derivs_min1_min2_max3, - mixed_derivs_max1_min2_max3); - } else { - spline_derivs_min2_max3_opt = std::nullopt; + derivs_min2_max3, + detail::derivs(mixed_derivs1_2_3, dmin, dmax)); } - ddc::Chunk spline_derivs_max2_min3_alloc( - m_spline_builder1.batched_spline_domain(batched_interpolation_deriv_domain2_3), - ddc::KokkosAllocator()); - auto spline_derivs_max2_min3 = spline_derivs_max2_min3_alloc.span_view(); - auto spline_derivs_max2_min3_opt = std::optional(spline_derivs_max2_min3.span_cview()); if constexpr (BcUpper2 == ddc::BoundCond::HERMITE || BcLower3 == ddc::BoundCond::HERMITE) { + auto const spline_derivs_max2_min3_strided + = detail::derivs(spline_derivs2_3, dmax, dmin); + auto const spline_derivs_max2_min3 + = detail::strided_to_discrete_domain_chunkspan(spline_derivs_max2_min3_strided); + + auto const derivs_max2_min3_strided + = detail::derivs(mixed_derivs2_3, dmax, dmin); + auto const derivs_max2_min3 + = detail::strided_to_discrete_domain_chunkspan(derivs_max2_min3_strided); m_spline_builder1( spline_derivs_max2_min3, - *mixed_derivs_max2_min3, - mixed_derivs_min1_max2_min3, - mixed_derivs_max1_max2_min3); - } else { - spline_derivs_max2_min3_opt = std::nullopt; + derivs_max2_min3, + detail::derivs(mixed_derivs1_2_3, dmax, dmin)); } - ddc::Chunk spline_derivs_max2_max3_alloc( - m_spline_builder1.batched_spline_domain(batched_interpolation_deriv_domain2_3), - ddc::KokkosAllocator()); - auto spline_derivs_max2_max3 = spline_derivs_max2_max3_alloc.span_view(); - auto spline_derivs_max2_max3_opt = std::optional(spline_derivs_max2_max3.span_cview()); if constexpr (BcUpper2 == ddc::BoundCond::HERMITE || BcUpper3 == ddc::BoundCond::HERMITE) { + auto const spline_derivs_max2_max3_strided + = detail::derivs(spline_derivs2_3, dmax, dmax); + auto const spline_derivs_max2_max3 + = detail::strided_to_discrete_domain_chunkspan(spline_derivs_max2_max3_strided); + + auto const derivs_max2_max3_strided + = detail::derivs(mixed_derivs2_3, dmax, dmax); + auto const derivs_max2_max3 + = detail::strided_to_discrete_domain_chunkspan(derivs_max2_max3_strided); m_spline_builder1( spline_derivs_max2_max3, - *mixed_derivs_max2_max3, - mixed_derivs_min1_max2_max3, - mixed_derivs_max1_max2_max3); - } else { - spline_derivs_max2_max3_opt = std::nullopt; + derivs_max2_max3, + detail::derivs(mixed_derivs1_2_3, dmax, dmax)); } // Spline1-approximate vals (to spline1) @@ -1048,19 +941,14 @@ operator()( ddc::KokkosAllocator()); ddc::ChunkSpan const spline1 = spline1_alloc.span_view(); - m_spline_builder1(spline1, vals, derivs_min1, derivs_max1); + m_spline_builder1(spline1, vals, derivs1); m_spline_builder_2_3( spline, spline1.span_cview(), - spline_derivs_min2_opt, - spline_derivs_max2_opt, - spline_derivs_min3_opt, - spline_derivs_max3_opt, - spline_derivs_min2_min3_opt, - spline_derivs_max2_min3_opt, - spline_derivs_min2_max3_opt, - spline_derivs_max2_max3_opt); + spline_derivs2.span_cview(), + spline_derivs3.span_cview(), + spline_derivs2_3.span_cview()); } } // namespace ddc diff --git a/tests/splines/batched_2d_spline_builder.cpp b/tests/splines/batched_2d_spline_builder.cpp index 1dab67ead..0e9efe1ad 100644 --- a/tests/splines/batched_2d_spline_builder.cpp +++ b/tests/splines/batched_2d_spline_builder.cpp @@ -4,9 +4,7 @@ #include #include -#if defined(BC_HERMITE) -# include -#endif + #if defined(BSPLINES_TYPE_UNIFORM) # include #endif @@ -184,21 +182,31 @@ void Batched2dSplineTest() ddc::DiscreteDomain const dom_vals(dom_vals_tmp, interpolation_domain1, interpolation_domain2); -#if defined(BC_HERMITE) // Create the derivs domain + int const nbc_i1 = s_degree / 2; + int const nbc_i2 = s_degree / 2; + ddc::DiscreteDomain> const - derivs_domain1(DElem>(1), DVect>(s_degree / 2)); + derivs_ddom1(DElem>(1), DVect>(nbc_i1)); ddc::DiscreteDomain> const - derivs_domain2(DElem>(1), DVect>(s_degree / 2)); + derivs_ddom2(DElem>(1), DVect>(nbc_i2)); ddc::DiscreteDomain, ddc::Deriv> const - derivs_domain(derivs_domain1, derivs_domain2); + derivs_ddom(derivs_ddom1, derivs_ddom2); - 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 + ddc::StridedDiscreteDomain, ddc::Deriv> const derivs_domain( + DElem, ddc::Deriv>(1, 1), + DVect, ddc::Deriv>(nbc_i1, nbc_i2), + DVect, ddc::Deriv>(1, 1)); + + auto const whole_derivs_domain1 = ddc::detail::get_whole_derivs_domain< + ddc::Deriv>(interpolation_domain1, dom_vals, nbc_i1); + auto const whole_derivs_domain2 = ddc::detail::get_whole_derivs_domain< + ddc::Deriv>(interpolation_domain2, dom_vals, nbc_i2); + auto const whole_derivs_domain = ddc::detail::get_whole_derivs_domain< + ddc::Deriv, + ddc::Deriv>(interpolation_domain, dom_vals, nbc_i1, nbc_i2); + + auto const dom_derivs = ddc::remove_dims_of(whole_derivs_domain); // Create a SplineBuilder over BSplines and batched along other dimensions using some boundary conditions ddc::SplineBuilder2D< @@ -236,141 +244,156 @@ void Batched2dSplineTest() 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(); + + ddc::Chunk derivs1_alloc(whole_derivs_domain1, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs1 = derivs1_alloc.span_view(); + if (s_bcl == ddc::BoundCond::HERMITE) { - ddc::Chunk derivs_1d_lhs1_host_alloc( - ddc::DiscreteDomain, DDimI2>(derivs_domain1, interpolation_domain2), + ddc::ChunkSpan const derivs1_lhs_view + = ddc::detail::derivs(derivs1, ddc::detail::dmin); + + ddc::Chunk derivs1_lhs_host_alloc( + ddc::DiscreteDomain, DDimI2>(derivs_ddom1, interpolation_domain2), ddc::HostAllocator()); - ddc::ChunkSpan const derivs_1d_lhs1_host = derivs_1d_lhs1_host_alloc.span_view(); + ddc::ChunkSpan const derivs1_lhs_host = derivs1_lhs_host_alloc.span_view(); ddc::for_each( - derivs_1d_lhs1_host.domain(), + derivs1_lhs_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); + derivs1_lhs_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(); + auto derivs1_lhs_alloc = ddc::create_mirror_view_and_copy(exec_space, derivs1_lhs_host); + ddc::ChunkSpan const derivs1_lhs = derivs1_lhs_alloc.span_view(); ddc::parallel_for_each( exec_space, - derivs_1d_lhs.domain(), + derivs1_lhs_view.domain(), KOKKOS_LAMBDA( - typename decltype(derivs_1d_lhs.domain())::discrete_element_type const e) { - derivs_1d_lhs(e) = derivs_1d_lhs1(DElem, DDimI2>(e)); + typename decltype(derivs1_lhs_view + .domain())::discrete_element_type const e) { + derivs1_lhs_view(e) = derivs1_lhs(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), + if (s_bcr == ddc::BoundCond::HERMITE) { + ddc::ChunkSpan const derivs1_rhs_view + = ddc::detail::derivs(derivs1, ddc::detail::dmax); + + ddc::Chunk derivs1_rhs_host_alloc( + ddc::DiscreteDomain, DDimI2>(derivs_ddom1, interpolation_domain2), ddc::HostAllocator()); - ddc::ChunkSpan const derivs_1d_rhs1_host = derivs_1d_rhs1_host_alloc.span_view(); + ddc::ChunkSpan const derivs1_rhs_host = derivs1_rhs_host_alloc.span_view(); ddc::for_each( - derivs_1d_rhs1_host.domain(), + derivs1_rhs_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); + derivs1_rhs_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(); + auto derivs1_rhs_alloc = ddc::create_mirror_view_and_copy(exec_space, derivs1_rhs_host); + ddc::ChunkSpan const derivs1_rhs = derivs1_rhs_alloc.span_view(); ddc::parallel_for_each( exec_space, - derivs_1d_rhs.domain(), + derivs1_rhs_view.domain(), KOKKOS_LAMBDA( - typename decltype(derivs_1d_rhs.domain())::discrete_element_type const e) { - derivs_1d_rhs(e) = derivs_1d_rhs1(DElem, DDimI2>(e)); + typename decltype(derivs1_rhs_view + .domain())::discrete_element_type const e) { + derivs1_rhs_view(e) = derivs1_rhs(DElem, DDimI2>(e)); }); } - ddc::Chunk derivs2_lhs_alloc(dom_derivs2, ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs2_lhs = derivs2_lhs_alloc.span_view(); + ddc::Chunk derivs2_alloc(whole_derivs_domain2, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs2 = derivs2_alloc.span_view(); + if (s_bcl == ddc::BoundCond::HERMITE) { - ddc::Chunk derivs2_lhs1_host_alloc( - ddc::DiscreteDomain>(interpolation_domain1, derivs_domain2), + ddc::ChunkSpan const derivs2_lhs_view + = ddc::detail::derivs(derivs2, ddc::detail::dmin); + + ddc::Chunk derivs2_lhs_host_alloc( + ddc::DiscreteDomain>(interpolation_domain1, derivs_ddom2), ddc::HostAllocator()); - ddc::ChunkSpan const derivs2_lhs1_host = derivs2_lhs1_host_alloc.span_view(); + ddc::ChunkSpan const derivs2_lhs_host = derivs2_lhs_host_alloc.span_view(); ddc::for_each( - derivs2_lhs1_host.domain(), + derivs2_lhs_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); + derivs2_lhs_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(); + auto derivs2_lhs_alloc = ddc::create_mirror_view_and_copy(exec_space, derivs2_lhs_host); + ddc::ChunkSpan const derivs2_lhs = derivs2_lhs_alloc.span_view(); ddc::parallel_for_each( exec_space, - derivs2_lhs.domain(), + derivs2_lhs_view.domain(), KOKKOS_LAMBDA( - typename decltype(derivs2_lhs.domain())::discrete_element_type const e) { - derivs2_lhs(e) = derivs2_lhs1(DElem>(e)); + typename decltype(derivs2_lhs_view + .domain())::discrete_element_type const e) { + derivs2_lhs_view(e) = derivs2_lhs(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), + if (s_bcr == ddc::BoundCond::HERMITE) { + ddc::ChunkSpan const derivs2_rhs_view + = ddc::detail::derivs(derivs2, ddc::detail::dmax); + + ddc::Chunk derivs2_rhs_host_alloc( + ddc::DiscreteDomain>(interpolation_domain1, derivs_ddom2), ddc::HostAllocator()); - ddc::ChunkSpan const derivs2_rhs1_host = derivs2_rhs1_host_alloc.span_view(); + ddc::ChunkSpan const derivs2_rhs_host = derivs2_rhs_host_alloc.span_view(); ddc::for_each( - derivs2_rhs1_host.domain(), + derivs2_rhs_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); + derivs2_rhs_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(); + auto derivs2_rhs_alloc = ddc::create_mirror_view_and_copy(exec_space, derivs2_rhs_host); + ddc::ChunkSpan const derivs2_rhs = derivs2_rhs_alloc.span_view(); ddc::parallel_for_each( exec_space, - derivs2_rhs.domain(), + derivs2_rhs_view.domain(), KOKKOS_LAMBDA( - typename decltype(derivs2_rhs.domain())::discrete_element_type const e) { - derivs2_rhs(e) = derivs2_rhs1(DElem>(e)); + typename decltype(derivs2_rhs_view + .domain())::discrete_element_type const e) { + derivs2_rhs_view(e) = derivs2_rhs(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(); + ddc::Chunk derivs_mixed_1_2_alloc( + whole_derivs_domain, + ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs_mixed_1_2 = derivs_mixed_1_2_alloc.span_view(); + + if (s_bcl == ddc::BoundCond::HERMITE || s_bcr == ddc::BoundCond::HERMITE) { + ddc::ChunkSpan const derivs_mixed_lhs1_lhs2_view = ddc::detail:: + derivs(derivs_mixed_1_2, ddc::detail::dmin, ddc::detail::dmin); + ddc::ChunkSpan const derivs_mixed_rhs1_lhs2_view = ddc::detail:: + derivs(derivs_mixed_1_2, ddc::detail::dmax, ddc::detail::dmin); + ddc::ChunkSpan const derivs_mixed_lhs1_rhs2_view = ddc::detail:: + derivs(derivs_mixed_1_2, ddc::detail::dmin, ddc::detail::dmax); + ddc::ChunkSpan const derivs_mixed_rhs1_rhs2_view = ddc::detail:: + derivs(derivs_mixed_1_2, ddc::detail::dmax, ddc::detail::dmax); + + ddc::Chunk derivs_mixed_lhs1_lhs2_host_alloc(derivs_ddom, ddc::HostAllocator()); + ddc::ChunkSpan const derivs_mixed_lhs1_lhs2_host + = derivs_mixed_lhs1_lhs2_host_alloc.span_view(); + ddc::Chunk derivs_mixed_rhs1_lhs2_host_alloc(derivs_ddom, ddc::HostAllocator()); + ddc::ChunkSpan const derivs_mixed_rhs1_lhs2_host + = derivs_mixed_rhs1_lhs2_host_alloc.span_view(); + ddc::Chunk derivs_mixed_lhs1_rhs2_host_alloc(derivs_ddom, ddc::HostAllocator()); + ddc::ChunkSpan const derivs_mixed_lhs1_rhs2_host + = derivs_mixed_lhs1_rhs2_host_alloc.span_view(); + ddc::Chunk derivs_mixed_rhs1_rhs2_host_alloc(derivs_ddom, ddc::HostAllocator()); + ddc::ChunkSpan const derivs_mixed_rhs1_rhs2_host + = derivs_mixed_rhs1_rhs2_host_alloc.span_view(); for (std::size_t ii = 1; ii < static_cast(derivs_domain.template extent>()) + 1; @@ -378,69 +401,59 @@ void Batched2dSplineTest() for (std::size_t jj = 1; jj < static_cast(derivs_domain.template extent>()) + 1; ++jj) { - derivs_mixed_lhs_lhs1_host( + derivs_mixed_lhs1_lhs2_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( + derivs_mixed_rhs1_lhs2_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( + derivs_mixed_lhs1_rhs2_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( + derivs_mixed_rhs1_rhs2_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(); + auto derivs_mixed_lhs1_lhs2_alloc + = ddc::create_mirror_view_and_copy(exec_space, derivs_mixed_lhs1_lhs2_host); + ddc::ChunkSpan const derivs_mixed_lhs1_lhs2 = derivs_mixed_lhs1_lhs2_alloc.span_view(); + auto derivs_mixed_rhs1_lhs2_alloc + = ddc::create_mirror_view_and_copy(exec_space, derivs_mixed_rhs1_lhs2_host); + ddc::ChunkSpan const derivs_mixed_rhs1_lhs2 = derivs_mixed_rhs1_lhs2_alloc.span_view(); + auto derivs_mixed_lhs1_rhs2_alloc + = ddc::create_mirror_view_and_copy(exec_space, derivs_mixed_lhs1_rhs2_host); + ddc::ChunkSpan const derivs_mixed_lhs1_rhs2 = derivs_mixed_lhs1_rhs2_alloc.span_view(); + auto derivs_mixed_rhs1_rhs2_alloc + = ddc::create_mirror_view_and_copy(exec_space, derivs_mixed_rhs1_rhs2_host); + ddc::ChunkSpan const derivs_mixed_rhs1_rhs2 = derivs_mixed_rhs1_rhs2_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)); + derivs_mixed_lhs1_lhs2_view(e) + = derivs_mixed_lhs1_lhs2(DElem, ddc::Deriv>(e)); + derivs_mixed_rhs1_lhs2_view(e) + = derivs_mixed_rhs1_lhs2(DElem, ddc::Deriv>(e)); + derivs_mixed_lhs1_rhs2_view(e) + = derivs_mixed_lhs1_rhs2(DElem, ddc::Deriv>(e)); + derivs_mixed_rhs1_rhs2_view(e) + = derivs_mixed_rhs1_rhs2(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 + derivs1.span_cview(), + derivs2.span_cview(), + derivs_mixed_1_2.span_cview()); // Instantiate a SplineEvaluator over interest dimension and batched along other dimensions #if defined(BC_PERIODIC) diff --git a/tests/splines/batched_3d_spline_builder.cpp b/tests/splines/batched_3d_spline_builder.cpp index e012fe2e4..2075272e1 100644 --- a/tests/splines/batched_3d_spline_builder.cpp +++ b/tests/splines/batched_3d_spline_builder.cpp @@ -208,14 +208,17 @@ void Batched3dSplineTest() interpolation_domain2, interpolation_domain3); -#if defined(BC_HERMITE) // Create the derivs domain + int const nbc_i1 = s_degree / 2; + int const nbc_i2 = s_degree / 2; + int const nbc_i3 = s_degree / 2; + ddc::DiscreteDomain> const - derivs_domain1(DElem>(1), DVect>(s_degree / 2)); + derivs_domain1(DElem>(1), DVect>(nbc_i1)); ddc::DiscreteDomain> const - derivs_domain2(DElem>(1), DVect>(s_degree / 2)); + derivs_domain2(DElem>(1), DVect>(nbc_i2)); ddc::DiscreteDomain> const - derivs_domain3(DElem>(1), DVect>(s_degree / 2)); + derivs_domain3(DElem>(1), DVect>(nbc_i3)); ddc::DiscreteDomain, ddc::Deriv, DDimI3> const derivs_domain12(derivs_domain1, derivs_domain2, interpolation_domain3); ddc::DiscreteDomain, ddc::Deriv> const @@ -233,10 +236,38 @@ void Batched3dSplineTest() auto const dom_derivs23 = ddc::replace_dim_of>(dom_derivs2, derivs_domain3); auto const dom_derivs13 - = ddc::replace_dim_of>(dom_derivs1, derivs_domain3); + = ddc::replace_dim_of>(dom_derivs3, derivs_domain1); auto const dom_derivs_all = ddc::replace_dim_of>(dom_derivs12, derivs_domain3); -#endif + + auto const whole_derivs_domain1 = ddc::detail::get_whole_derivs_domain< + ddc::Deriv>(interpolation_domain1, dom_vals, nbc_i1); + auto const whole_derivs_domain2 = ddc::detail::get_whole_derivs_domain< + ddc::Deriv>(interpolation_domain2, dom_vals, nbc_i2); + auto const whole_derivs_domain3 = ddc::detail::get_whole_derivs_domain< + ddc::Deriv>(interpolation_domain3, dom_vals, nbc_i3); + auto const whole_derivs_domain1_2 + = ddc::detail::get_whole_derivs_domain, ddc::Deriv>( + ddc::select(interpolation_domain), + dom_vals, + nbc_i1, + nbc_i2); + auto const whole_derivs_domain2_3 + = ddc::detail::get_whole_derivs_domain, ddc::Deriv>( + ddc::select(interpolation_domain), + dom_vals, + nbc_i2, + nbc_i3); + auto const whole_derivs_domain1_3 + = ddc::detail::get_whole_derivs_domain, ddc::Deriv>( + ddc::select(interpolation_domain), + dom_vals, + nbc_i1, + nbc_i3); + auto const whole_derivs_domain = ddc::detail::get_whole_derivs_domain< + ddc::Deriv, + ddc::Deriv, + ddc::Deriv>(interpolation_domain, dom_vals, nbc_i1, nbc_i2, nbc_i3); // Create a SplineBuilder over BSplines and batched along other dimensions using some boundary conditions ddc::SplineBuilder3D< @@ -278,13 +309,16 @@ void Batched3dSplineTest() 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 derivs1_lhs_view_alloc(dom_derivs1, ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs1_lhs_view = derivs1_lhs_view_alloc.span_view(); + ddc::Chunk derivs1_alloc(whole_derivs_domain1, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs1 = derivs1_alloc.span_view(); + if (s_bcl == ddc::BoundCond::HERMITE) { + ddc::ChunkSpan const derivs1_lhs_view + = ddc::detail::derivs(derivs1, ddc::detail::dmin); + ddc::Chunk derivs1_lhs_host_alloc( ddc::DiscreteDomain< ddc::Deriv, @@ -314,9 +348,10 @@ void Batched3dSplineTest() }); } - ddc::Chunk derivs1_rhs_view_alloc(dom_derivs1, ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs1_rhs_view = derivs1_rhs_view_alloc.span_view(); if (s_bcr == ddc::BoundCond::HERMITE) { + ddc::ChunkSpan const derivs1_rhs_view + = ddc::detail::derivs(derivs1, ddc::detail::dmax); + ddc::Chunk derivs1_rhs_host_alloc( ddc::DiscreteDomain< ddc::Deriv, @@ -346,9 +381,13 @@ void Batched3dSplineTest() }); } - ddc::Chunk derivs2_lhs_view_alloc(dom_derivs2, ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs2_lhs_view = derivs2_lhs_view_alloc.span_view(); + ddc::Chunk derivs2_alloc(whole_derivs_domain2, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs2 = derivs2_alloc.span_view(); + if (s_bcl == ddc::BoundCond::HERMITE) { + ddc::ChunkSpan const derivs2_lhs_view + = ddc::detail::derivs(derivs2, ddc::detail::dmin); + ddc::Chunk derivs2_lhs_host_alloc( ddc::DiscreteDomain< DDimI1, @@ -379,9 +418,10 @@ void Batched3dSplineTest() }); } - ddc::Chunk derivs2_rhs_view_alloc(dom_derivs2, ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs2_rhs_view = derivs2_rhs_view_alloc.span_view(); if (s_bcr == ddc::BoundCond::HERMITE) { + ddc::ChunkSpan const derivs2_rhs_view + = ddc::detail::derivs(derivs2, ddc::detail::dmax); + ddc::Chunk derivs2_rhs_host_alloc( ddc::DiscreteDomain< DDimI1, @@ -412,9 +452,13 @@ void Batched3dSplineTest() }); } - ddc::Chunk derivs3_lhs_view_alloc(dom_derivs3, ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs3_lhs_view = derivs3_lhs_view_alloc.span_view(); + ddc::Chunk derivs3_alloc(whole_derivs_domain3, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs3 = derivs3_alloc.span_view(); + if (s_bcl == ddc::BoundCond::HERMITE) { + ddc::ChunkSpan const derivs3_lhs_view + = ddc::detail::derivs(derivs3, ddc::detail::dmin); + ddc::Chunk derivs3_lhs_host_alloc( ddc::DiscreteDomain< DDimI1, @@ -446,9 +490,10 @@ void Batched3dSplineTest() }); } - ddc::Chunk derivs3_rhs_view_alloc(dom_derivs3, ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs3_rhs_view = derivs3_rhs_view_alloc.span_view(); if (s_bcr == ddc::BoundCond::HERMITE) { + ddc::ChunkSpan const derivs3_rhs_view + = ddc::detail::derivs(derivs3, ddc::detail::dmax); + ddc::Chunk derivs3_rhs_host_alloc( ddc::DiscreteDomain< DDimI1, @@ -480,28 +525,21 @@ void Batched3dSplineTest() }); } - ddc::Chunk derivs_mixed_lhs1_lhs2_view_alloc( - dom_derivs12, - ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_lhs1_lhs2_view - = derivs_mixed_lhs1_lhs2_view_alloc.span_view(); - ddc::Chunk derivs_mixed_rhs1_lhs2_view_alloc( - dom_derivs12, + ddc::Chunk derivs_mixed_1_2_alloc( + whole_derivs_domain1_2, ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_rhs1_lhs2_view - = derivs_mixed_rhs1_lhs2_view_alloc.span_view(); - ddc::Chunk derivs_mixed_lhs1_rhs2_view_alloc( - dom_derivs12, - ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_lhs1_rhs2_view - = derivs_mixed_lhs1_rhs2_view_alloc.span_view(); - ddc::Chunk derivs_mixed_rhs1_rhs2_view_alloc( - dom_derivs12, - ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_rhs1_rhs2_view - = derivs_mixed_rhs1_rhs2_view_alloc.span_view(); + ddc::ChunkSpan const derivs_mixed_1_2 = derivs_mixed_1_2_alloc.span_view(); if (s_bcl == ddc::BoundCond::HERMITE) { + ddc::ChunkSpan const derivs_mixed_lhs1_lhs2_view = ddc::detail:: + derivs(derivs_mixed_1_2, ddc::detail::dmin, ddc::detail::dmin); + ddc::ChunkSpan const derivs_mixed_rhs1_lhs2_view = ddc::detail:: + derivs(derivs_mixed_1_2, ddc::detail::dmax, ddc::detail::dmin); + ddc::ChunkSpan const derivs_mixed_lhs1_rhs2_view = ddc::detail:: + derivs(derivs_mixed_1_2, ddc::detail::dmin, ddc::detail::dmax); + ddc::ChunkSpan const derivs_mixed_rhs1_rhs2_view = ddc::detail:: + derivs(derivs_mixed_1_2, ddc::detail::dmax, ddc::detail::dmax); + ddc::Chunk derivs_mixed_lhs1_lhs2_host_alloc(derivs_domain12, ddc::HostAllocator()); ddc::ChunkSpan const derivs_mixed_lhs1_lhs2_host = derivs_mixed_lhs1_lhs2_host_alloc.span_view(); @@ -580,28 +618,21 @@ void Batched3dSplineTest() }); } - ddc::Chunk derivs_mixed_lhs2_lhs3_view_alloc( - dom_derivs23, - ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_lhs2_lhs3_view - = derivs_mixed_lhs2_lhs3_view_alloc.span_view(); - ddc::Chunk derivs_mixed_rhs2_lhs3_view_alloc( - dom_derivs23, - ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_rhs2_lhs3_view - = derivs_mixed_rhs2_lhs3_view_alloc.span_view(); - ddc::Chunk derivs_mixed_lhs2_rhs3_view_alloc( - dom_derivs23, + ddc::Chunk derivs_mixed_2_3_alloc( + whole_derivs_domain2_3, ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_lhs2_rhs3_view - = derivs_mixed_lhs2_rhs3_view_alloc.span_view(); - ddc::Chunk derivs_mixed_rhs2_rhs3_view_alloc( - dom_derivs23, - ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_rhs2_rhs3_view - = derivs_mixed_rhs2_rhs3_view_alloc.span_view(); + ddc::ChunkSpan const derivs_mixed_2_3 = derivs_mixed_2_3_alloc.span_view(); if (s_bcl == ddc::BoundCond::HERMITE) { + ddc::ChunkSpan const derivs_mixed_lhs2_lhs3_view = ddc::detail:: + derivs(derivs_mixed_2_3, ddc::detail::dmin, ddc::detail::dmin); + ddc::ChunkSpan const derivs_mixed_rhs2_lhs3_view = ddc::detail:: + derivs(derivs_mixed_2_3, ddc::detail::dmax, ddc::detail::dmin); + ddc::ChunkSpan const derivs_mixed_lhs2_rhs3_view = ddc::detail:: + derivs(derivs_mixed_2_3, ddc::detail::dmin, ddc::detail::dmax); + ddc::ChunkSpan const derivs_mixed_rhs2_rhs3_view = ddc::detail:: + derivs(derivs_mixed_2_3, ddc::detail::dmax, ddc::detail::dmax); + ddc::Chunk derivs_mixed_lhs2_lhs3_host_alloc(derivs_domain23, ddc::HostAllocator()); ddc::ChunkSpan const derivs_mixed_lhs2_lhs3_host = derivs_mixed_lhs2_lhs3_host_alloc.span_view(); @@ -680,28 +711,21 @@ void Batched3dSplineTest() }); } - ddc::Chunk derivs_mixed_lhs1_lhs3_view_alloc( - dom_derivs13, - ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_lhs1_lhs3_view - = derivs_mixed_lhs1_lhs3_view_alloc.span_view(); - ddc::Chunk derivs_mixed_rhs1_lhs3_view_alloc( - dom_derivs13, - ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_rhs1_lhs3_view - = derivs_mixed_rhs1_lhs3_view_alloc.span_view(); - ddc::Chunk derivs_mixed_lhs1_rhs3_view_alloc( - dom_derivs13, - ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_lhs1_rhs3_view - = derivs_mixed_lhs1_rhs3_view_alloc.span_view(); - ddc::Chunk derivs_mixed_rhs1_rhs3_view_alloc( - dom_derivs13, + ddc::Chunk derivs_mixed_1_3_alloc( + whole_derivs_domain1_3, ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_rhs1_rhs3_view - = derivs_mixed_rhs1_rhs3_view_alloc.span_view(); + ddc::ChunkSpan const derivs_mixed_1_3 = derivs_mixed_1_3_alloc.span_view(); if (s_bcl == ddc::BoundCond::HERMITE) { + ddc::ChunkSpan const derivs_mixed_lhs1_lhs3_view = ddc::detail:: + derivs(derivs_mixed_1_3, ddc::detail::dmin, ddc::detail::dmin); + ddc::ChunkSpan const derivs_mixed_rhs1_lhs3_view = ddc::detail:: + derivs(derivs_mixed_1_3, ddc::detail::dmax, ddc::detail::dmin); + ddc::ChunkSpan const derivs_mixed_lhs1_rhs3_view = ddc::detail:: + derivs(derivs_mixed_1_3, ddc::detail::dmin, ddc::detail::dmax); + ddc::ChunkSpan const derivs_mixed_rhs1_rhs3_view = ddc::detail:: + derivs(derivs_mixed_1_3, ddc::detail::dmax, ddc::detail::dmax); + ddc::Chunk derivs_mixed_lhs1_lhs3_host_alloc(derivs_domain13, ddc::HostAllocator()); ddc::ChunkSpan const derivs_mixed_lhs1_lhs3_host = derivs_mixed_lhs1_lhs3_host_alloc.span_view(); @@ -780,48 +804,53 @@ void Batched3dSplineTest() }); } - ddc::Chunk derivs_mixed_lhs1_lhs2_lhs3_view_alloc( - dom_derivs_all, - ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_lhs1_lhs2_lhs3_view - = derivs_mixed_lhs1_lhs2_lhs3_view_alloc.span_view(); - ddc::Chunk derivs_mixed_rhs1_lhs2_lhs3_view_alloc( - dom_derivs_all, - ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_rhs1_lhs2_lhs3_view - = derivs_mixed_rhs1_lhs2_lhs3_view_alloc.span_view(); - ddc::Chunk derivs_mixed_lhs1_rhs2_lhs3_view_alloc( - dom_derivs_all, - ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_lhs1_rhs2_lhs3_view - = derivs_mixed_lhs1_rhs2_lhs3_view_alloc.span_view(); - ddc::Chunk derivs_mixed_rhs1_rhs2_lhs3_view_alloc( - dom_derivs_all, - ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_rhs1_rhs2_lhs3_view - = derivs_mixed_rhs1_rhs2_lhs3_view_alloc.span_view(); - ddc::Chunk derivs_mixed_lhs1_lhs2_rhs3_view_alloc( - dom_derivs_all, + ddc::Chunk derivs_mixed_1_2_3_alloc( + whole_derivs_domain, ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_lhs1_lhs2_rhs3_view - = derivs_mixed_lhs1_lhs2_rhs3_view_alloc.span_view(); - ddc::Chunk derivs_mixed_rhs1_lhs2_rhs3_view_alloc( - dom_derivs_all, - ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_rhs1_lhs2_rhs3_view - = derivs_mixed_rhs1_lhs2_rhs3_view_alloc.span_view(); - ddc::Chunk derivs_mixed_lhs1_rhs2_rhs3_view_alloc( - dom_derivs_all, - ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_lhs1_rhs2_rhs3_view - = derivs_mixed_lhs1_rhs2_rhs3_view_alloc.span_view(); - ddc::Chunk derivs_mixed_rhs1_rhs2_rhs3_view_alloc( - dom_derivs_all, - ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_mixed_rhs1_rhs2_rhs3_view - = derivs_mixed_rhs1_rhs2_rhs3_view_alloc.span_view(); + ddc::ChunkSpan const derivs_mixed_1_2_3 = derivs_mixed_1_2_3_alloc.span_view(); + + if (s_bcl == ddc::BoundCond::HERMITE || s_bcr == ddc::BoundCond::HERMITE) { + ddc::ChunkSpan const derivs_mixed_lhs1_lhs2_lhs3_view = ddc::detail:: + derivs(derivs_mixed_1_2_3, + ddc::detail::dmin, + ddc::detail::dmin, + ddc::detail::dmin); + ddc::ChunkSpan const derivs_mixed_rhs1_lhs2_lhs3_view = ddc::detail:: + derivs(derivs_mixed_1_2_3, + ddc::detail::dmax, + ddc::detail::dmin, + ddc::detail::dmin); + ddc::ChunkSpan const derivs_mixed_lhs1_rhs2_lhs3_view = ddc::detail:: + derivs(derivs_mixed_1_2_3, + ddc::detail::dmin, + ddc::detail::dmax, + ddc::detail::dmin); + ddc::ChunkSpan const derivs_mixed_rhs1_rhs2_lhs3_view = ddc::detail:: + derivs(derivs_mixed_1_2_3, + ddc::detail::dmax, + ddc::detail::dmax, + ddc::detail::dmin); + ddc::ChunkSpan const derivs_mixed_lhs1_lhs2_rhs3_view = ddc::detail:: + derivs(derivs_mixed_1_2_3, + ddc::detail::dmin, + ddc::detail::dmin, + ddc::detail::dmax); + ddc::ChunkSpan const derivs_mixed_rhs1_lhs2_rhs3_view = ddc::detail:: + derivs(derivs_mixed_1_2_3, + ddc::detail::dmax, + ddc::detail::dmin, + ddc::detail::dmax); + ddc::ChunkSpan const derivs_mixed_lhs1_rhs2_rhs3_view = ddc::detail:: + derivs(derivs_mixed_1_2_3, + ddc::detail::dmin, + ddc::detail::dmax, + ddc::detail::dmax); + ddc::ChunkSpan const derivs_mixed_rhs1_rhs2_rhs3_view = ddc::detail:: + derivs(derivs_mixed_1_2_3, + ddc::detail::dmax, + ddc::detail::dmax, + ddc::detail::dmax); - if (s_bcl == ddc::BoundCond::HERMITE) { ddc::Chunk derivs_mixed_lhs1_lhs2_lhs3_host_alloc( derivs_domain_all, ddc::HostAllocator()); @@ -1005,48 +1034,24 @@ void Batched3dSplineTest() DElem, ddc::Deriv, 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(derivs1_lhs_view.span_cview()), - std::optional(derivs1_rhs_view.span_cview()), - std::optional(derivs2_lhs_view.span_cview()), - std::optional(derivs2_rhs_view.span_cview()), - std::optional(derivs3_lhs_view.span_cview()), - std::optional(derivs3_rhs_view.span_cview()), - std::optional(derivs_mixed_lhs1_lhs2_view.span_cview()), - std::optional(derivs_mixed_rhs1_lhs2_view.span_cview()), - std::optional(derivs_mixed_lhs1_rhs2_view.span_cview()), - std::optional(derivs_mixed_rhs1_rhs2_view.span_cview()), - std::optional(derivs_mixed_lhs2_lhs3_view.span_cview()), - std::optional(derivs_mixed_rhs2_lhs3_view.span_cview()), - std::optional(derivs_mixed_lhs2_rhs3_view.span_cview()), - std::optional(derivs_mixed_rhs2_rhs3_view.span_cview()), - std::optional(derivs_mixed_lhs1_lhs3_view.span_cview()), - std::optional(derivs_mixed_rhs1_lhs3_view.span_cview()), - std::optional(derivs_mixed_lhs1_rhs3_view.span_cview()), - std::optional(derivs_mixed_rhs1_rhs3_view.span_cview()), - std::optional(derivs_mixed_lhs1_lhs2_lhs3_view.span_cview()), - std::optional(derivs_mixed_rhs1_lhs2_lhs3_view.span_cview()), - std::optional(derivs_mixed_lhs1_rhs2_lhs3_view.span_cview()), - std::optional(derivs_mixed_rhs1_rhs2_lhs3_view.span_cview()), - std::optional(derivs_mixed_lhs1_lhs2_rhs3_view.span_cview()), - std::optional(derivs_mixed_rhs1_lhs2_rhs3_view.span_cview()), - std::optional(derivs_mixed_lhs1_rhs2_rhs3_view.span_cview()), - std::optional(derivs_mixed_rhs1_rhs2_rhs3_view.span_cview())); -#else - spline_builder(coef, vals.span_cview()); -#endif - - // Instantiate a SplineEvaluator over interest dimension and batched along other dimensions + derivs1.span_cview(), + derivs2.span_cview(), + derivs3.span_cview(), + derivs_mixed_1_2.span_cview(), + derivs_mixed_2_3.span_cview(), + derivs_mixed_1_3.span_cview(), + derivs_mixed_1_2_3.span_cview()); + +// Instantiate a SplineEvaluator over interest dimension and batched along other dimensions #if defined(BC_PERIODIC) using extrapolation_rule_1_type = ddc::PeriodicExtrapolationRule; using extrapolation_rule_2_type = ddc::PeriodicExtrapolationRule; diff --git a/tests/splines/batched_spline_builder.cpp b/tests/splines/batched_spline_builder.cpp index d98c21d69..6fa1098b7 100644 --- a/tests/splines/batched_spline_builder.cpp +++ b/tests/splines/batched_spline_builder.cpp @@ -1,4 +1,4 @@ -// Copyright (C) The DDC development team, see COPYRIGHT.md file +// Copyright (C) The DDC development team, see COPYRIGHT.md filebatched // // SPDX-License-Identifier: MIT @@ -160,12 +160,14 @@ void BatchedSplineTest() ddc::DiscreteDomain(DElem(0), DVect(ncells))...); ddc::DiscreteDomain const dom_vals(dom_vals_tmp, interpolation_domain); -#if defined(BC_HERMITE) // Create the derivs domain - 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 + ddc::StridedDiscreteDomain> const derivs_domain( + DElem>(1), + DVect>(s_degree_x / 2), + DVect>(1)); + + auto const whole_derivs_domain = ddc::detail::get_whole_derivs_domain< + ddc::Deriv>(interpolation_domain, dom_vals, s_degree_x / 2); // Create a SplineBuilder over BSplines and batched along other dimensions using some boundary conditions ddc::SplineBuilder< @@ -197,69 +199,56 @@ void BatchedSplineTest() 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(); + ddc::Chunk derivs_alloc(whole_derivs_domain, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs = derivs_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)) + ddc::ChunkSpan const derivs_lhs_view = derivs[interpolation_domain.front()]; + + ddc::Chunk derivs_lhs_host_alloc(derivs_domain, ddc::HostAllocator()); + ddc::ChunkSpan const derivs_lhs_host = derivs_lhs_host_alloc.span_view(); + for (int ii = 1; ii < derivs_lhs_host.domain().template extent>() + 1; ++ii) { + derivs_lhs_host(typename decltype(derivs_lhs_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(); + auto derivs_lhs_alloc = ddc::create_mirror_view_and_copy(exec_space, derivs_lhs_host); + ddc::ChunkSpan const derivs_lhs = derivs_lhs_alloc.span_view(); ddc::parallel_for_each( exec_space, - derivs_lhs.domain(), + derivs_lhs_view.domain(), KOKKOS_LAMBDA( - typename decltype(derivs_lhs.domain())::discrete_element_type const e) { - derivs_lhs(e) = derivs_lhs1(DElem>(e)); - }); + typename decltype(derivs_lhs_view.domain())::discrete_element_type const + e) { derivs_lhs_view(e) = derivs_lhs(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)) + ddc::ChunkSpan const derivs_rhs_view = derivs[interpolation_domain.back()]; + + ddc::Chunk derivs_rhs_host_alloc(derivs_domain, ddc::HostAllocator()); + ddc::ChunkSpan const derivs_rhs_host = derivs_rhs_host_alloc.span_view(); + for (int ii = 1; ii < derivs_rhs_host.domain().template extent>() + 1; ++ii) { + derivs_rhs_host(typename decltype(derivs_rhs_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(); + auto derivs_rhs_alloc = ddc::create_mirror_view_and_copy(exec_space, derivs_rhs_host); + ddc::ChunkSpan const derivs_rhs = derivs_rhs_alloc.span_view(); ddc::parallel_for_each( exec_space, - derivs_rhs.domain(), + derivs_rhs_view.domain(), KOKKOS_LAMBDA( - typename decltype(derivs_rhs.domain())::discrete_element_type const e) { - derivs_rhs(e) = derivs_rhs1(DElem>(e)); - }); + typename decltype(derivs_rhs_view.domain())::discrete_element_type const + e) { derivs_rhs_view(e) = derivs_rhs(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 + spline_builder(coef, vals.span_cview(), derivs.span_cview()); // Instantiate a SplineEvaluator over interest dimension and batched along other dimensions #if defined(BC_PERIODIC) diff --git a/tests/splines/extrapolation_rule.cpp b/tests/splines/extrapolation_rule.cpp index f1a218ea6..84dee0026 100644 --- a/tests/splines/extrapolation_rule.cpp +++ b/tests/splines/extrapolation_rule.cpp @@ -198,6 +198,23 @@ void ExtrapolationRuleSplineTest() = spline_builder.interpolation_domain(); auto const dom_spline = spline_builder.batched_spline_domain(dom_vals); + // Create empty ChunkSpans for the derivatives, as they are not used here + ddc::ChunkSpan< + double const, + ddc::StridedDiscreteDomain, DDimI2>, + Kokkos::layout_right, + MemorySpace> const derivs1 {}; + ddc::ChunkSpan< + double const, + ddc::StridedDiscreteDomain>, + Kokkos::layout_right, + MemorySpace> const derivs2 {}; + ddc::ChunkSpan< + double const, + ddc::StridedDiscreteDomain, ddc::Deriv>, + Kokkos::layout_right, + MemorySpace> const mixed_derivs1_2 {}; + // 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(); @@ -220,7 +237,12 @@ void ExtrapolationRuleSplineTest() ddc::ChunkSpan const coef = coef_alloc.span_view(); // Finally compute the spline by filling `coef` - spline_builder(coef, vals.span_cview()); + spline_builder( + coef, + vals.span_cview(), + derivs1.span_cview(), + derivs2.span_cview(), + mixed_derivs1_2.span_cview()); // Instantiate a SplineEvaluator over interest dimension and batched along other dimensions #if defined(ER_NULL) diff --git a/tests/splines/multiple_batch_domain_2d_spline_builder.cpp b/tests/splines/multiple_batch_domain_2d_spline_builder.cpp index a2316cabd..040c95a5b 100644 --- a/tests/splines/multiple_batch_domain_2d_spline_builder.cpp +++ b/tests/splines/multiple_batch_domain_2d_spline_builder.cpp @@ -5,9 +5,6 @@ #include #include #include -#if defined(BC_HERMITE) -# include -#endif #if defined(BSPLINES_TYPE_UNIFORM) # include #endif @@ -176,28 +173,40 @@ std::tuple ComputeEvaluationError( using I1 = typename DDimI1::continuous_dimension_type; using I2 = typename DDimI2::continuous_dimension_type; -#if defined(BC_HERMITE) + // Compute useful domains + ddc::DiscreteDomain const dom_interpolation + = spline_builder.interpolation_domain(); + auto const dom_spline = spline_builder.batched_spline_domain(dom_vals); + ddc::DiscreteDomain const interpolation_domain1(dom_vals); ddc::DiscreteDomain const interpolation_domain2(dom_vals); + ddc::DiscreteDomain const + interpolation_domain(interpolation_domain1, interpolation_domain2); + // Create the derivs domain + int const nbc_i1 = s_degree / 2; + int const nbc_i2 = s_degree / 2; + ddc::DiscreteDomain> const - derivs_domain1(DElem>(1), DVect>(s_degree / 2)); + derivs_ddom1(DElem>(1), DVect>(nbc_i1)); ddc::DiscreteDomain> const - derivs_domain2(DElem>(1), DVect>(s_degree / 2)); + derivs_ddom2(DElem>(1), DVect>(nbc_i2)); 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); + derivs_ddom(derivs_ddom1, derivs_ddom2); + + ddc::StridedDiscreteDomain, ddc::Deriv> const derivs_domain( + DElem, ddc::Deriv>(1, 1), + DVect, ddc::Deriv>(nbc_i1, nbc_i2), + DVect, ddc::Deriv>(1, 1)); + + auto const whole_derivs_domain1 = ddc::detail::get_whole_derivs_domain< + ddc::Deriv>(interpolation_domain1, dom_vals, nbc_i1); + auto const whole_derivs_domain2 = ddc::detail::get_whole_derivs_domain< + ddc::Deriv>(interpolation_domain2, dom_vals, nbc_i2); + auto const whole_derivs_domain = ddc::detail::get_whole_derivs_domain< + ddc::Deriv, + ddc::Deriv>(interpolation_domain, dom_vals, nbc_i1, nbc_i2); + auto const dom_derivs = ddc::remove_dims_of(whole_derivs_domain); // 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()); @@ -215,141 +224,156 @@ std::tuple ComputeEvaluationError( 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(); + + ddc::Chunk derivs1_alloc(whole_derivs_domain1, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs1 = derivs1_alloc.span_view(); + if (s_bcl == ddc::BoundCond::HERMITE) { - ddc::Chunk derivs_1d_lhs1_host_alloc( - ddc::DiscreteDomain, DDimI2>(derivs_domain1, interpolation_domain2), + ddc::ChunkSpan const derivs1_lhs_view + = ddc::detail::derivs(derivs1, ddc::detail::dmin); + + ddc::Chunk derivs1_lhs_host_alloc( + ddc::DiscreteDomain, DDimI2>(derivs_ddom1, interpolation_domain2), ddc::HostAllocator()); - ddc::ChunkSpan const derivs_1d_lhs1_host = derivs_1d_lhs1_host_alloc.span_view(); + ddc::ChunkSpan const derivs1_lhs_host = derivs1_lhs_host_alloc.span_view(); ddc::for_each( - derivs_1d_lhs1_host.domain(), + derivs1_lhs_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); + derivs1_lhs_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(); + auto derivs1_lhs_alloc = ddc::create_mirror_view_and_copy(exec_space, derivs1_lhs_host); + ddc::ChunkSpan const derivs1_lhs = derivs1_lhs_alloc.span_view(); ddc::parallel_for_each( exec_space, - derivs_1d_lhs.domain(), + derivs1_lhs_view.domain(), KOKKOS_LAMBDA( - typename decltype(derivs_1d_lhs.domain())::discrete_element_type const e) { - derivs_1d_lhs(e) = derivs_1d_lhs1(DElem, DDimI2>(e)); + typename decltype(derivs1_lhs_view + .domain())::discrete_element_type const e) { + derivs1_lhs_view(e) = derivs1_lhs(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::ChunkSpan const derivs1_rhs_view + = ddc::detail::derivs(derivs1, ddc::detail::dmax); + + ddc::Chunk derivs1_rhs_host_alloc( + ddc::DiscreteDomain, DDimI2>(derivs_ddom1, interpolation_domain2), ddc::HostAllocator()); - ddc::ChunkSpan const derivs_1d_rhs1_host = derivs_1d_rhs1_host_alloc.span_view(); + ddc::ChunkSpan const derivs1_rhs_host = derivs1_rhs_host_alloc.span_view(); ddc::for_each( - derivs_1d_rhs1_host.domain(), + derivs1_rhs_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); + derivs1_rhs_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(); + auto derivs1_rhs_alloc = ddc::create_mirror_view_and_copy(exec_space, derivs1_rhs_host); + ddc::ChunkSpan const derivs1_rhs = derivs1_rhs_alloc.span_view(); ddc::parallel_for_each( exec_space, - derivs_1d_rhs.domain(), + derivs1_rhs_view.domain(), KOKKOS_LAMBDA( - typename decltype(derivs_1d_rhs.domain())::discrete_element_type const e) { - derivs_1d_rhs(e) = derivs_1d_rhs1(DElem, DDimI2>(e)); + typename decltype(derivs1_rhs_view + .domain())::discrete_element_type const e) { + derivs1_rhs_view(e) = derivs1_rhs(DElem, DDimI2>(e)); }); } - ddc::Chunk derivs2_lhs_alloc(dom_derivs2, ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs2_lhs = derivs2_lhs_alloc.span_view(); + ddc::Chunk derivs2_alloc(whole_derivs_domain2, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs2 = derivs2_alloc.span_view(); + if (s_bcl == ddc::BoundCond::HERMITE) { - ddc::Chunk derivs2_lhs1_host_alloc( - ddc::DiscreteDomain>(interpolation_domain1, derivs_domain2), + ddc::ChunkSpan const derivs2_lhs_view + = ddc::detail::derivs(derivs2, ddc::detail::dmin); + + ddc::Chunk derivs2_lhs_host_alloc( + ddc::DiscreteDomain>(interpolation_domain1, derivs_ddom2), ddc::HostAllocator()); - ddc::ChunkSpan const derivs2_lhs1_host = derivs2_lhs1_host_alloc.span_view(); + ddc::ChunkSpan const derivs2_lhs_host = derivs2_lhs_host_alloc.span_view(); ddc::for_each( - derivs2_lhs1_host.domain(), + derivs2_lhs_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); + derivs2_lhs_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(); + auto derivs2_lhs_alloc = ddc::create_mirror_view_and_copy(exec_space, derivs2_lhs_host); + ddc::ChunkSpan const derivs2_lhs = derivs2_lhs_alloc.span_view(); ddc::parallel_for_each( exec_space, - derivs2_lhs.domain(), + derivs2_lhs_view.domain(), KOKKOS_LAMBDA( - typename decltype(derivs2_lhs.domain())::discrete_element_type const e) { - derivs2_lhs(e) = derivs2_lhs1(DElem>(e)); + typename decltype(derivs2_lhs_view + .domain())::discrete_element_type const e) { + derivs2_lhs_view(e) = derivs2_lhs(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::ChunkSpan const derivs2_rhs_view + = ddc::detail::derivs(derivs2, ddc::detail::dmax); + + ddc::Chunk derivs2_rhs_host_alloc( + ddc::DiscreteDomain>(interpolation_domain1, derivs_ddom2), ddc::HostAllocator()); - ddc::ChunkSpan const derivs2_rhs1_host = derivs2_rhs1_host_alloc.span_view(); + ddc::ChunkSpan const derivs2_rhs_host = derivs2_rhs_host_alloc.span_view(); ddc::for_each( - derivs2_rhs1_host.domain(), + derivs2_rhs_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); + derivs2_rhs_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(); + auto derivs2_rhs_alloc = ddc::create_mirror_view_and_copy(exec_space, derivs2_rhs_host); + ddc::ChunkSpan const derivs2_rhs = derivs2_rhs_alloc.span_view(); ddc::parallel_for_each( exec_space, - derivs2_rhs.domain(), + derivs2_rhs_view.domain(), KOKKOS_LAMBDA( - typename decltype(derivs2_rhs.domain())::discrete_element_type const e) { - derivs2_rhs(e) = derivs2_rhs1(DElem>(e)); + typename decltype(derivs2_rhs_view + .domain())::discrete_element_type const e) { + derivs2_rhs_view(e) = derivs2_rhs(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(); + ddc::Chunk derivs_mixed_1_2_alloc( + whole_derivs_domain, + ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs_mixed_1_2 = derivs_mixed_1_2_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(); + ddc::ChunkSpan const derivs_mixed_lhs1_lhs2_view = ddc::detail:: + derivs(derivs_mixed_1_2, ddc::detail::dmin, ddc::detail::dmin); + ddc::ChunkSpan const derivs_mixed_rhs1_lhs2_view = ddc::detail:: + derivs(derivs_mixed_1_2, ddc::detail::dmax, ddc::detail::dmin); + ddc::ChunkSpan const derivs_mixed_lhs1_rhs2_view = ddc::detail:: + derivs(derivs_mixed_1_2, ddc::detail::dmin, ddc::detail::dmax); + ddc::ChunkSpan const derivs_mixed_rhs1_rhs2_view = ddc::detail:: + derivs(derivs_mixed_1_2, ddc::detail::dmax, ddc::detail::dmax); + + ddc::Chunk derivs_mixed_lhs1_lhs2_host_alloc(derivs_ddom, ddc::HostAllocator()); + ddc::ChunkSpan const derivs_mixed_lhs1_lhs2_host + = derivs_mixed_lhs1_lhs2_host_alloc.span_view(); + ddc::Chunk derivs_mixed_rhs1_lhs2_host_alloc(derivs_ddom, ddc::HostAllocator()); + ddc::ChunkSpan const derivs_mixed_rhs1_lhs2_host + = derivs_mixed_rhs1_lhs2_host_alloc.span_view(); + ddc::Chunk derivs_mixed_lhs1_rhs2_host_alloc(derivs_ddom, ddc::HostAllocator()); + ddc::ChunkSpan const derivs_mixed_lhs1_rhs2_host + = derivs_mixed_lhs1_rhs2_host_alloc.span_view(); + ddc::Chunk derivs_mixed_rhs1_rhs2_host_alloc(derivs_ddom, ddc::HostAllocator()); + ddc::ChunkSpan const derivs_mixed_rhs1_rhs2_host + = derivs_mixed_rhs1_rhs2_host_alloc.span_view(); for (std::size_t ii = 1; ii < static_cast(derivs_domain.template extent>()) + 1; @@ -357,69 +381,59 @@ std::tuple ComputeEvaluationError( for (std::size_t jj = 1; jj < static_cast(derivs_domain.template extent>()) + 1; ++jj) { - derivs_mixed_lhs_lhs1_host( + derivs_mixed_lhs1_lhs2_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( + derivs_mixed_rhs1_lhs2_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( + derivs_mixed_lhs1_rhs2_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( + derivs_mixed_rhs1_rhs2_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(); + auto derivs_mixed_lhs1_lhs2_alloc + = ddc::create_mirror_view_and_copy(exec_space, derivs_mixed_lhs1_lhs2_host); + ddc::ChunkSpan const derivs_mixed_lhs1_lhs2 = derivs_mixed_lhs1_lhs2_alloc.span_view(); + auto derivs_mixed_rhs1_lhs2_alloc + = ddc::create_mirror_view_and_copy(exec_space, derivs_mixed_rhs1_lhs2_host); + ddc::ChunkSpan const derivs_mixed_rhs1_lhs2 = derivs_mixed_rhs1_lhs2_alloc.span_view(); + auto derivs_mixed_lhs1_rhs2_alloc + = ddc::create_mirror_view_and_copy(exec_space, derivs_mixed_lhs1_rhs2_host); + ddc::ChunkSpan const derivs_mixed_lhs1_rhs2 = derivs_mixed_lhs1_rhs2_alloc.span_view(); + auto derivs_mixed_rhs1_rhs2_alloc + = ddc::create_mirror_view_and_copy(exec_space, derivs_mixed_rhs1_rhs2_host); + ddc::ChunkSpan const derivs_mixed_rhs1_rhs2 = derivs_mixed_rhs1_rhs2_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)); + derivs_mixed_lhs1_lhs2_view(e) + = derivs_mixed_lhs1_lhs2(DElem, ddc::Deriv>(e)); + derivs_mixed_rhs1_lhs2_view(e) + = derivs_mixed_rhs1_lhs2(DElem, ddc::Deriv>(e)); + derivs_mixed_lhs1_rhs2_view(e) + = derivs_mixed_lhs1_rhs2(DElem, ddc::Deriv>(e)); + derivs_mixed_rhs1_rhs2_view(e) + = derivs_mixed_rhs1_rhs2(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 + derivs1.span_cview(), + derivs2.span_cview(), + derivs_mixed_1_2.span_cview()); // Instantiate chunk of coordinates of dom_interpolation ddc::Chunk coords_eval_alloc(dom_vals, ddc::KokkosAllocator, MemorySpace>()); diff --git a/tests/splines/multiple_batch_domain_spline_builder.cpp b/tests/splines/multiple_batch_domain_spline_builder.cpp index f42bbab20..74d3b5053 100644 --- a/tests/splines/multiple_batch_domain_spline_builder.cpp +++ b/tests/splines/multiple_batch_domain_spline_builder.cpp @@ -157,18 +157,20 @@ std::tuple ComputeEvaluationError( { 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); + // Create the derivs domains + ddc::StridedDiscreteDomain> const derivs_domain( + DElem>(1), + DVect>(s_degree_x / 2), + DVect>(1)); + + auto const whole_derivs_domain = ddc::detail::get_whole_derivs_domain< + ddc::Deriv>(dom_interpolation, dom_vals, s_degree_x / 2); + // 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(); @@ -183,69 +185,56 @@ std::tuple ComputeEvaluationError( 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(); + ddc::Chunk derivs_alloc(whole_derivs_domain, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs = derivs_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)) + ddc::ChunkSpan const derivs_lhs_view = derivs[dom_interpolation.front()]; + + ddc::Chunk derivs_lhs_host_alloc(derivs_domain, ddc::HostAllocator()); + ddc::ChunkSpan const derivs_lhs_host = derivs_lhs_host_alloc.span_view(); + for (int ii = 1; ii < derivs_lhs_host.domain().template extent>() + 1; ++ii) { + derivs_lhs_host(typename decltype(derivs_lhs_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(); + auto derivs_lhs_alloc = ddc::create_mirror_view_and_copy(exec_space, derivs_lhs_host); + ddc::ChunkSpan const derivs_lhs = derivs_lhs_alloc.span_view(); ddc::parallel_for_each( exec_space, - derivs_lhs.domain(), + derivs_lhs_view.domain(), KOKKOS_LAMBDA( - typename decltype(derivs_lhs.domain())::discrete_element_type const e) { - derivs_lhs(e) = derivs_lhs1(DElem>(e)); - }); + typename decltype(derivs_lhs_view.domain())::discrete_element_type const + e) { derivs_lhs_view(e) = derivs_lhs(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)) + ddc::ChunkSpan const derivs_rhs_view = derivs[dom_interpolation.back()]; + + ddc::Chunk derivs_rhs_host_alloc(derivs_domain, ddc::HostAllocator()); + ddc::ChunkSpan const derivs_rhs_host = derivs_rhs_host_alloc.span_view(); + for (int ii = 1; ii < derivs_rhs_host.domain().template extent>() + 1; ++ii) { + derivs_rhs_host(typename decltype(derivs_rhs_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(); + auto derivs_rhs_alloc = ddc::create_mirror_view_and_copy(exec_space, derivs_rhs_host); + ddc::ChunkSpan const derivs_rhs = derivs_rhs_alloc.span_view(); ddc::parallel_for_each( exec_space, - derivs_rhs.domain(), + derivs_rhs_view.domain(), KOKKOS_LAMBDA( - typename decltype(derivs_rhs.domain())::discrete_element_type const e) { - derivs_rhs(e) = derivs_rhs1(DElem>(e)); - }); + typename decltype(derivs_rhs_view.domain())::discrete_element_type const + e) { derivs_rhs_view(e) = derivs_rhs(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 + spline_builder(coef, vals.span_cview(), derivs.span_cview()); // Instantiate chunk of coordinates of dom_interpolation ddc::Chunk coords_eval_alloc(dom_vals, ddc::KokkosAllocator, MemorySpace>()); diff --git a/tests/splines/non_periodic_spline_builder.cpp b/tests/splines/non_periodic_spline_builder.cpp index 4b5755f12..c712845c1 100644 --- a/tests/splines/non_periodic_spline_builder.cpp +++ b/tests/splines/non_periodic_spline_builder.cpp @@ -107,10 +107,13 @@ void TestNonPeriodicSplineBuilderTestIdentity() // The chunk is filled with garbage data, we need to initialize it ddc::Chunk coef(dom_bsplines_x, ddc::KokkosAllocator()); - // 3. Create the interpolation domain + // 3. Create the interpolation and deriv domains ddc::init_discrete_space(GrevillePoints::get_sampling()); ddc::DiscreteDomain const interpolation_domain(GrevillePoints::get_domain()); + auto const whole_derivs_domain = ddc::detail::get_whole_derivs_domain< + ddc::Deriv>(interpolation_domain, s_degree_x / 2); + // 4. Create a SplineBuilder over BSplines using some boundary conditions ddc::SplineBuilder< execution_space, @@ -131,8 +134,11 @@ void TestNonPeriodicSplineBuilderTestIdentity() KOKKOS_LAMBDA(DElemX const ix) { yvals(ix) = evaluator(ddc::coordinate(ix)); }); int const shift = s_degree_x % 2; // shift = 0 for even order, 1 for odd order - ddc::Chunk derivs_lhs_alloc(derivs_domain, ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_lhs = derivs_lhs_alloc.span_view(); + + ddc::Chunk derivs_alloc(whole_derivs_domain, ddc::KokkosAllocator()); + ddc::ChunkSpan const derivs = derivs_alloc.span_view(); + + ddc::ChunkSpan const derivs_lhs = derivs[interpolation_domain.front()]; if (s_bcl == ddc::BoundCond::HERMITE) { ddc::parallel_for_each( execution_space(), @@ -142,8 +148,7 @@ void TestNonPeriodicSplineBuilderTestIdentity() }); } - ddc::Chunk derivs_rhs_alloc(derivs_domain, ddc::KokkosAllocator()); - ddc::ChunkSpan const derivs_rhs = derivs_rhs_alloc.span_view(); + ddc::ChunkSpan const derivs_rhs = derivs[interpolation_domain.back()]; if (s_bcr == ddc::BoundCond::HERMITE) { ddc::parallel_for_each( execution_space(), @@ -154,19 +159,7 @@ void TestNonPeriodicSplineBuilderTestIdentity() } // 6. Finally build the spline by filling `coef` -#if defined(BCL_HERMITE) - std::optional const deriv_l(derivs_lhs.span_cview()); -#else - decltype(std::optional(derivs_lhs.span_cview())) const deriv_l = std::nullopt; -#endif - -#if defined(BCR_HERMITE) - std::optional const deriv_r(derivs_rhs.span_cview()); -#else - decltype(std::optional(derivs_rhs.span_cview())) const deriv_r = std::nullopt; -#endif - - spline_builder(coef.span_view(), yvals.span_cview(), deriv_l, deriv_r); + spline_builder(coef.span_view(), yvals.span_cview(), derivs.span_cview()); // 7. Create a SplineEvaluator to evaluate the spline at any point in the domain of the BSplines ddc::NullExtrapolationRule const extrapolation_rule; diff --git a/tests/splines/periodic_spline_builder.cpp b/tests/splines/periodic_spline_builder.cpp index d508f7229..681e4160b 100644 --- a/tests/splines/periodic_spline_builder.cpp +++ b/tests/splines/periodic_spline_builder.cpp @@ -105,10 +105,17 @@ void TestPeriodicSplineBuilderTestIdentity() yvals.domain(), KOKKOS_LAMBDA(DElemX const ix) { yvals(ix) = evaluator(ddc::coordinate(ix)); }); - // 6. Finally build the spline by filling `coef` - spline_builder(coef.span_view(), yvals.span_cview()); + // 6. Create a empty chunk for the derivatives, as they are not used with periodic boundary conditions + ddc::ChunkSpan< + double const, + ddc::StridedDiscreteDomain>, + Kokkos::layout_right, + memory_space> const derivs {}; - // 7. Create a SplineEvaluator to evaluate the spline at any point in the domain of the BSplines + // 7. Finally build the spline by filling `coef` + spline_builder(coef.span_view(), yvals.span_cview(), derivs); + + // 8. Create a SplineEvaluator to evaluate the spline at any point in the domain of the BSplines ddc::PeriodicExtrapolationRule const periodic_extrapolation; ddc::SplineEvaluator< execution_space, diff --git a/tests/splines/periodicity_spline_builder.cpp b/tests/splines/periodicity_spline_builder.cpp index c4162ba0b..8084bceb2 100644 --- a/tests/splines/periodicity_spline_builder.cpp +++ b/tests/splines/periodicity_spline_builder.cpp @@ -148,8 +148,15 @@ void PeriodicitySplineBuilderTest() ddc::Chunk coef_alloc(dom_bsplines, ddc::KokkosAllocator()); ddc::ChunkSpan const coef = coef_alloc.span_view(); + // Instantiate empty chunk of derivatives + ddc::ChunkSpan< + double const, + ddc::StridedDiscreteDomain, ddc::Deriv>, + Kokkos::layout_right, + MemorySpace> const derivs {}; + // Finally compute the spline by filling `coef` - spline_builder(coef, vals.span_cview()); + spline_builder(coef, vals.span_cview(), derivs.span_cview()); // Instantiate a SplineEvaluator over interest dimension and batched along other dimensions ddc::PeriodicExtrapolationRule const extrapolation_rule;