Skip to content

WIP: Use StridedDiscreteDomain with SplineBuilder #892

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 5 commits into
base: main
Choose a base branch
from
Draft

WIP: Use StridedDiscreteDomain with SplineBuilder #892

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
120e471
Allow derivative layout to be different to values layout
EmilyBourne Jul 18, 2025
799f757
Add missing strides() method. Allow conversion from DiscreteDomain to…
EmilyBourne Jul 18, 2025
02fb8e2
Fix slice of ChunkSpan on StridedDiscreteDomain. Fixes #891
EmilyBourne Jul 18, 2025
f3f1cce
Start adding test
EmilyBourne Jul 18, 2025
a29eb51
Merge remote-tracking branch 'origin/main' into devel-issue886
EmilyBourne Jul 23, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 6 additions & 6 deletions include/ddc/kernels/splines/spline_builder.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -468,7 +468,7 @@ class SplineBuilder
* @param[in] derivs_xmax The values of the derivatives at the upper boundary
* (used only with BoundCond::HERMITE upper boundary condition).
*/
template <class Layout, class BatchedInterpolationDDom>
template <class Layout, class BatchedInterpolationDDom, class LayoutDeriv = Kokkos::layout_right>
void operator()(
ddc::ChunkSpan<
double,
Expand All @@ -479,13 +479,13 @@ class SplineBuilder
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> derivs_xmin
= std::nullopt,
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> derivs_xmax
= std::nullopt) const;

Expand Down Expand Up @@ -793,7 +793,7 @@ template <
ddc::BoundCond BcLower,
ddc::BoundCond BcUpper,
SplineSolver Solver>
template <class Layout, class BatchedInterpolationDDom>
template <class Layout, class BatchedInterpolationDDom, class LayoutDeriv>
void SplineBuilder<ExecSpace, MemorySpace, BSplines, InterpolationDDim, BcLower, BcUpper, Solver>::
operator()(
ddc::ChunkSpan<
Expand All @@ -805,12 +805,12 @@ operator()(
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> const derivs_xmin,
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> const derivs_xmax) const
{
auto const batched_interpolation_domain = vals.domain();
Expand Down
36 changes: 18 additions & 18 deletions include/ddc/kernels/splines/spline_builder_2d.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -405,7 +405,7 @@ class SplineBuilder2D
* The values of the the cross-derivatives at the upper boundary in the first dimension
* and the upper boundary in the second dimension.
*/
template <class Layout, class BatchedInterpolationDDom>
template <class Layout, class BatchedInterpolationDDom, class LayoutDeriv = Kokkos::layout_right>
void operator()(
ddc::ChunkSpan<
double,
Expand All @@ -416,49 +416,49 @@ class SplineBuilder2D
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type1<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> derivs_min1
= std::nullopt,
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type1<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> derivs_max1
= std::nullopt,
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type2<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> derivs_min2
= std::nullopt,
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type2<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> derivs_max2
= std::nullopt,
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> mixed_derivs_min1_min2
= std::nullopt,
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> mixed_derivs_max1_min2
= std::nullopt,
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> mixed_derivs_min1_max2
= std::nullopt,
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> mixed_derivs_max1_max2
= std::nullopt) const;
};
Expand All @@ -476,7 +476,7 @@ template <
ddc::BoundCond BcLower2,
ddc::BoundCond BcUpper2,
ddc::SplineSolver Solver>
template <class Layout, class BatchedInterpolationDDom>
template <class Layout, class BatchedInterpolationDDom, class LayoutDeriv>
void SplineBuilder2D<
ExecSpace,
MemorySpace,
Expand All @@ -499,42 +499,42 @@ operator()(
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type1<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> const derivs_min1,
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type1<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> const derivs_max1,
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type2<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> const derivs_min2,
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type2<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> const derivs_max2,
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> const mixed_derivs_min1_min2,
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> const mixed_derivs_max1_min2,
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> const mixed_derivs_min1_max2,
std::optional<ddc::ChunkSpan<
double const,
batched_derivs_domain_type<BatchedInterpolationDDom>,
Layout,
LayoutDeriv,
memory_space>> const mixed_derivs_max1_max2) const
{
auto const batched_interpolation_domain = vals.domain();
Expand Down
23 changes: 22 additions & 1 deletion include/ddc/strided_discrete_domain.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@

#include "detail/type_seq.hpp"

#include "discrete_domain.hpp"
#include "discrete_element.hpp"
#include "discrete_vector.hpp"

Expand Down Expand Up @@ -95,14 +96,22 @@ class StridedDiscreteDomain
/// Construct a StridedDiscreteDomain by copies and merge of domains
template <
class... DDoms,
class = std::enable_if_t<(is_strided_discrete_domain_v<DDoms> && ...)>>
std::enable_if_t<(is_strided_discrete_domain_v<DDoms> && ...), bool> = true>
KOKKOS_FUNCTION constexpr explicit StridedDiscreteDomain(DDoms const&... domains)
: m_element_begin(domains.front()...)
, m_extents(domains.extents()...)
, m_strides(domains.strides()...)
{
}

/// Construct a StridedDiscreteDomain from a DiscreteDomain
KOKKOS_FUNCTION constexpr explicit StridedDiscreteDomain(DiscreteDomain<DDims...> const& domain)
: m_element_begin(domain.front())
, m_extents(domain.extents())
, m_strides((DiscreteVector<DDims>{1})...)
{
}

/** Construct a StridedDiscreteDomain starting from element_begin with size points.
* @param element_begin the lower bound in each direction
* @param extents the number of points in each direction
Expand Down Expand Up @@ -340,6 +349,13 @@ class StridedDiscreteDomain<>
{
}

// Construct a StridedDiscreteDomain from a reordered copy of `domain`
template <class... ODDims>
KOKKOS_FUNCTION constexpr explicit StridedDiscreteDomain(
[[maybe_unused]] DiscreteDomain<ODDims...> const& domain)
{
}

/** Construct a StridedDiscreteDomain starting from element_begin with size points.
* @param element_begin the lower bound in each direction
* @param size the number of points in each direction
Expand Down Expand Up @@ -387,6 +403,11 @@ class StridedDiscreteDomain<>
return {};
}

static KOKKOS_FUNCTION constexpr discrete_vector_type strides() noexcept
{
return {};
}

static KOKKOS_FUNCTION constexpr discrete_element_type front() noexcept
{
return {};
Expand Down
31 changes: 31 additions & 0 deletions tests/splines/batched_spline_builder.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -165,6 +165,13 @@ void BatchedSplineTest()
ddc::DiscreteDomain<ddc::Deriv<I>> const
derivs_domain(DElem<ddc::Deriv<I>>(1), DVect<ddc::Deriv<I>>(s_degree_x / 2));
auto const dom_derivs = ddc::replace_dim_of<DDimI, ddc::Deriv<I>>(dom_vals, derivs_domain);
// Create the derivs domain
ddc::StridedDiscreteDomain<DDimI, ddc::Deriv<I>> const
derivs_domain_strided(DElem<DDimI, ddc::Deriv<I>>(interpolation_domain.front(), derivs_domain.front()),
DVect<DDimI, ddc::Deriv<I>>(DVect<DDimI>(1), derivs_domain.extents()),
DVect<DDimI, ddc::Deriv<I>>(interpolation_domain.extents() - 1, DVect<ddc::Deriv<I>>(1)));
ddc::remove_dims_of_t<ddc::StridedDiscreteDomain<DDims...>, DDimI> const dom_vals_tmp_strided(dom_vals_tmp);
ddc::StridedDiscreteDomain<DDims..., ddc::Deriv<I>> const dom_derivs_strided(dom_vals_tmp_strided, derivs_domain_strided);
#endif

// Create a SplineBuilder over BSplines<I> and batched along other dimensions using some boundary conditions
Expand Down Expand Up @@ -200,6 +207,8 @@ void BatchedSplineTest()
#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_strided_alloc(dom_derivs_strided, ddc::KokkosAllocator<double, MemorySpace>());
ddc::ChunkSpan const derivs_strided = derivs_strided_alloc.span_view();
ddc::Chunk derivs_lhs_alloc(dom_derivs, ddc::KokkosAllocator<double, MemorySpace>());
ddc::ChunkSpan const derivs_lhs = derivs_lhs_alloc.span_view();
if (s_bcl == ddc::BoundCond::HERMITE) {
Expand All @@ -220,6 +229,8 @@ void BatchedSplineTest()
typename decltype(derivs_lhs.domain())::discrete_element_type const e) {
derivs_lhs(e) = derivs_lhs1(DElem<ddc::Deriv<I>>(e));
});

Kokkos::deep_copy(derivs_strided[interpolation_domain.front()].allocation_kokkos_view(), derivs_lhs.allocation_kokkos_view());
}

ddc::Chunk derivs_rhs_alloc(dom_derivs, ddc::KokkosAllocator<double, MemorySpace>());
Expand All @@ -243,6 +254,7 @@ void BatchedSplineTest()
typename decltype(derivs_rhs.domain())::discrete_element_type const e) {
derivs_rhs(e) = derivs_rhs1(DElem<ddc::Deriv<I>>(e));
});
Kokkos::deep_copy(derivs_strided[interpolation_domain.back()].allocation_kokkos_view(), derivs_rhs.allocation_kokkos_view());
}
#endif

Expand All @@ -257,6 +269,25 @@ void BatchedSplineTest()
vals.span_cview(),
std::optional(derivs_lhs.span_cview()),
std::optional(derivs_rhs.span_cview()));

{
ddc::Chunk coef_2_alloc(dom_spline, ddc::KokkosAllocator<double, MemorySpace>());
ddc::ChunkSpan const coef_2 = coef_2_alloc.span_view();
spline_builder(
coef_2,
vals.span_cview(),
std::optional(derivs_strided[interpolation_domain.front()].span_cview()),
std::optional(derivs_strided[interpolation_domain.back()].span_cview()));
Comment on lines +279 to +280
Copy link
Member

@tpadioleau tpadioleau Jul 23, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think a way to properly do this would be to do a cartesian product between a 1D strided domain and a 1D contiguous domain. Once you slice by a DiscreteElement in the strided dimension you would recover the contiguous domain over derivatives.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We could allow strided domains in the spline operator but it feels artificial, the operator still requires all derivatives, so it would work only when the stride is unitary.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think a way to properly do this would be to do a cartesian product between a 1D strided domain and a 1D contiguous domain

I agree. I couldn't find a way to do this though. I guess it's not possible right now?

Copy link
Collaborator Author

@EmilyBourne EmilyBourne Jul 23, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We could allow strided domains in the spline operator but it feels artificial, the operator still requires all derivatives, so it would work only when the stride is unitary.

Yes. If a cross product of strided and unstrided domains was possible then there would be no need to allow strides domains here.

However I wonder if it is worth considering keeping this in mind for the handling of ND splines. If I had to pass derivs_strided here instead of derivs_strided[front] and derivs_strided[back] then it would be much easier to calculate the number of arguments required for the ND case (N! N(N+1)/2 derivatives+ 1 argument for the values)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess it's not possible right now?

No indeed, but this is something to work on to be able to release v1 in my opinion.

it would be much easier to calculate the number of arguments required for the ND case (N! derivatives+ 1 argument for the values)

Do you claim that it is possible to pass all derivatives in one chunkspan ? Still we would need to check at runtime that the dimension of derivatives has unitary stride right ?

If we had this product feature, in 2D we could consider passing a single

ChunkSpan<double, Product<StridedDiscreteDomain<DDimI1>, DiscreteDomain<Deriv<I1>>>> derivs_dim1;

instead of two arrays derivs_min1 and derivs_max1. Do you think we can do better than this with only Chunk/ChunkSpan data structures ? This would also allow to check that the user passes derivatives at the correct location ?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you claim that it is possible to pass all derivatives in one chunkspan ?

Not all derivatives, but all derivatives of the same type. So 3 chunks for 2D:
derivs1 (domain = $\textbraceleft x_{1,0}, x_{1,N} \textbraceright \times [x_{2,0}, x_{2,N}]$), derivs2 (domain = $[ x_{1,0}, x_{1,N} ] \times \textbraceleft x_{2,0}, x_{2,N}\textbraceright$), cross_derivs (domain = $\textbraceleft x_{1,0}, x_{1,N} \textbraceright \times \textbraceleft x_{2,0}, x_{2,N}\textbraceright$)

Still we would need to check at runtime that the dimension of derivatives has unitary stride right ?

Yes probably. This could be in an assert I guess?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The formula for the number of derivatives required is:
$$\sum_{i=1}^{N} \binom{N}{i}$$
if we pass all derivatives of the same type in one chunkspan.
If we pass explicitly as we have done up to now the formula is:
$$\sum_{i=1}^{N} 2^i \binom{N}{i}$$

N N derivs with strided N derivs without strided
1 1 2
2 3 8
3 7 26
4 15 80
5 31 242

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Still we would need to check at runtime that the dimension of derivatives has unitary stride right ?

Yes probably. This could be in an assert I guess?

Yes we can but my point is to say that we need to develop this product feature. We cannot propagate too many unnecessary StridedDiscreteDomain, there could be a performance penalty.

double const max_norm_error = ddc::parallel_transform_reduce(
exec_space,
coef.domain(),
0.,
ddc::reducer::max<double>(),
KOKKOS_LAMBDA(DElem<DDims...> const e) {
return Kokkos::abs(coef(e) - coef_2(e));
});
EXPECT_LE(max_norm_error, 1e-14);
}
#else
spline_builder(coef, vals.span_cview());
#endif
Expand Down
Loading