diff --git a/include/ddc/kernels/splines/bsplines_uniform.hpp b/include/ddc/kernels/splines/bsplines_uniform.hpp index 73c33ea49..ff2042a61 100644 --- a/include/ddc/kernels/splines/bsplines_uniform.hpp +++ b/include/ddc/kernels/splines/bsplines_uniform.hpp @@ -401,9 +401,9 @@ KOKKOS_INLINE_FUNCTION ddc::DiscreteElement UniformBSplines:: Impl::eval_basis( DSpan1D values, ddc::Coordinate const& x, - [[maybe_unused]] std::size_t const deg) const + [[maybe_unused]] std::size_t const degree) const { - assert(values.size() == deg + 1); + assert(values.size() == degree + 1); double offset; int jmin; diff --git a/include/ddc/kernels/splines/spline_builder.hpp b/include/ddc/kernels/splines/spline_builder.hpp index d10262eef..304f5dc20 100644 --- a/include/ddc/kernels/splines/spline_builder.hpp +++ b/include/ddc/kernels/splines/spline_builder.hpp @@ -165,7 +165,7 @@ class SplineBuilder private: batched_interpolation_domain_type m_batched_interpolation_domain; - int m_offset; + int m_offset = 0; double m_dx; // average cell size for normalization of derivatives @@ -173,7 +173,7 @@ class SplineBuilder std::unique_ptr> matrix; /// Calculate offset so that the matrix is diagonally dominant - int compute_offset(interpolation_domain_type const& interpolation_domain); + void compute_offset(interpolation_domain_type const& interpolation_domain, int& offset); public: /** @@ -197,13 +197,15 @@ class SplineBuilder std::optional cols_per_chunk = std::nullopt, std::optional preconditioner_max_block_size = std::nullopt) : m_batched_interpolation_domain(batched_interpolation_domain) - , m_offset(compute_offset(interpolation_domain())) , m_dx((ddc::discrete_space().rmax() - ddc::discrete_space().rmin()) / ddc::discrete_space().ncells()) { static_assert( ((BcLower == BoundCond::PERIODIC) == (BcUpper == BoundCond::PERIODIC)), "Incompatible boundary conditions"); + check_valid_grid(); + + compute_offset(interpolation_domain(), m_offset); // Calculate block sizes int lower_block_size; @@ -467,6 +469,11 @@ class SplineBuilder std::optional preconditioner_max_block_size = std::nullopt); void build_matrix_system(); + + void check_valid_grid(); + + template + static void check_n_points_in_cell(int n_points_in_cell, KnotElement current_cell_end_idx); }; template < @@ -478,7 +485,7 @@ template < ddc::BoundCond BcUpper, SplineSolver Solver, class... IDimX> -int SplineBuilder< +void SplineBuilder< ExecSpace, MemorySpace, BSplines, @@ -486,9 +493,9 @@ int SplineBuilder< BcLower, BcUpper, Solver, - IDimX...>::compute_offset(interpolation_domain_type const& interpolation_domain) + IDimX...>:: + compute_offset(interpolation_domain_type const& interpolation_domain, int& offset) { - int offset; if constexpr (bsplines_type::is_periodic()) { // Calculate offset so that the matrix is diagonally dominant std::array values_ptr; @@ -511,7 +518,6 @@ int SplineBuilder< } else { offset = 0; } - return offset; } template < @@ -1066,4 +1072,91 @@ SplineBuilder< std::move(coefficients_derivs_xmax_out)); } +template < + class ExecSpace, + class MemorySpace, + class BSplines, + class InterpolationDDim, + ddc::BoundCond BcLower, + ddc::BoundCond BcUpper, + SplineSolver Solver, + class... IDimX> +template +void SplineBuilder< + ExecSpace, + MemorySpace, + BSplines, + InterpolationDDim, + BcLower, + BcUpper, + Solver, + IDimX...>:: + check_n_points_in_cell(int const n_points_in_cell, KnotElement const current_cell_end_idx) +{ + if (n_points_in_cell > BSplines::degree() + 1) { + KnotElement const rmin_idx = ddc::discrete_space().break_point_domain().front(); + int const failed_cell = (current_cell_end_idx - rmin_idx).value(); + throw std::runtime_error( + "The spline problem is overconstrained. There are " + + std::to_string(n_points_in_cell) + " points in the " + std::to_string(failed_cell) + + "-th cell."); + } +} + +template < + class ExecSpace, + class MemorySpace, + class BSplines, + class InterpolationDDim, + ddc::BoundCond BcLower, + ddc::BoundCond BcUpper, + SplineSolver Solver, + class... IDimX> +void SplineBuilder< + ExecSpace, + MemorySpace, + BSplines, + InterpolationDDim, + BcLower, + BcUpper, + Solver, + IDimX...>::check_valid_grid() +{ + std::size_t const n_interp_points = interpolation_domain().size(); + std::size_t const expected_npoints + = ddc::discrete_space().nbasis() - s_nbc_xmin - s_nbc_xmax; + if (n_interp_points != expected_npoints) { + throw std::runtime_error( + "Incorrect number of points supplied to NonUniformInterpolationPoints. " + "(Received : " + + std::to_string(n_interp_points) + + ", expected : " + std::to_string(expected_npoints)); + } + int n_points_in_cell = 0; + auto current_cell_end_idx = ddc::discrete_space().break_point_domain().front() + 1; + ddc::for_each(interpolation_domain(), [&](auto idx) { + ddc::Coordinate const point = ddc::coordinate(idx); + if (point > ddc::coordinate(current_cell_end_idx)) { + // Check the points found in the previous cell + check_n_points_in_cell(n_points_in_cell, current_cell_end_idx); + // Initialise the number of points in the subsequent cell, including the new point + n_points_in_cell = 1; + // Move to the next cell + current_cell_end_idx += 1; + } else if (point == ddc::coordinate(current_cell_end_idx)) { + // Check the points found in the previous cell including the point on the boundary + check_n_points_in_cell(n_points_in_cell + 1, current_cell_end_idx); + // Initialise the number of points in the subsequent cell, including the point on the boundary + n_points_in_cell = 1; + // Move to the next cell + current_cell_end_idx += 1; + } else { + // Indicate that the point is in the cell + n_points_in_cell += 1; + } + }); + // Check the number of points in the final cell + check_n_points_in_cell(n_points_in_cell, current_cell_end_idx); +} + } // namespace ddc diff --git a/tests/splines/CMakeLists.txt b/tests/splines/CMakeLists.txt index 12d99dc48..c11def96a 100644 --- a/tests/splines/CMakeLists.txt +++ b/tests/splines/CMakeLists.txt @@ -9,7 +9,13 @@ include(GoogleTest) set(SPLINES_TEST_DEGREE_MIN 3 CACHE STRING "Minimum degree to test splines.") set(SPLINES_TEST_DEGREE_MAX 3 CACHE STRING "Maximum degree to test splines.") -add_executable(splines_tests ../main.cpp knots_as_interpolation_points.cpp view.cpp) +add_executable( + splines_tests + ../main.cpp + knots_as_interpolation_points.cpp + spline_builder.cpp + view.cpp +) target_compile_features(splines_tests PUBLIC cxx_std_17) target_link_libraries(splines_tests PUBLIC DDC::core DDC::splines GTest::gtest) diff --git a/tests/splines/spline_builder.cpp b/tests/splines/spline_builder.cpp new file mode 100644 index 000000000..8a1d60ffd --- /dev/null +++ b/tests/splines/spline_builder.cpp @@ -0,0 +1,139 @@ +// Copyright (C) The DDC development team, see COPYRIGHT.md file +// +// SPDX-License-Identifier: MIT + +#include +#include +#include + +#include +#include + +#include + +#include + +struct DimX +{ + static constexpr bool PERIODIC = true; +}; + +using CoordX = ddc::Coordinate; + +static constexpr std::size_t s_degree_x = 2; + +struct BSplinesX : ddc::UniformBSplines +{ +}; + +struct IDimX : ddc::NonUniformPointSampling +{ +}; + +using execution_space = Kokkos::DefaultHostExecutionSpace; +using memory_space = Kokkos::HostSpace; + +TEST(SplineBuilder, ShortInterpolationGrid) +{ + CoordX const x0(0.); + CoordX const xN(1.); + std::size_t const ncells = 5; + + ddc::init_discrete_space(x0, xN, ncells); + + // One point missing + std::vector const range {0.1, 0.3, 0.5, 0.7}; + + ddc::DiscreteDomain const interpolation_domain + = ddc::init_discrete_space(IDimX::init(range)); + + EXPECT_THROW( + (ddc::SplineBuilder< + execution_space, + memory_space, + BSplinesX, + IDimX, + ddc::BoundCond::PERIODIC, + ddc::BoundCond::PERIODIC, + ddc::SplineSolver::GINKGO, + IDimX>(interpolation_domain)), + std::runtime_error); +} + +TEST(SplineBuilder, LongInterpolationGrid) +{ + CoordX const x0(0.); + CoordX const xN(1.); + std::size_t const ncells = 5; + + ddc::init_discrete_space(x0, xN, ncells); + + // One point too much + std::vector const range {0.1, 0.3, 0.5, 0.7, 0.9, 0.95}; + + ddc::DiscreteDomain const interpolation_domain + = ddc::init_discrete_space(IDimX::init(range)); + + EXPECT_THROW( + (ddc::SplineBuilder< + execution_space, + memory_space, + BSplinesX, + IDimX, + ddc::BoundCond::PERIODIC, + ddc::BoundCond::PERIODIC, + ddc::SplineSolver::GINKGO, + IDimX>(interpolation_domain)), + std::runtime_error); +} + +TEST(SplineBuilder, BadShapeInterpolationGrid) +{ + CoordX const x0(0.); + CoordX const xN(1.); + std::size_t const ncells = 5; + + ddc::init_discrete_space(x0, xN, ncells); + + // All points end up in the first cell ]0, 0.2[ + std::vector const range {0.1, 0.11, 0.12, 0.13, 0.14}; + + ddc::DiscreteDomain const interpolation_domain + = ddc::init_discrete_space(IDimX::init(range)); + + EXPECT_THROW( + (ddc::SplineBuilder< + execution_space, + memory_space, + BSplinesX, + IDimX, + ddc::BoundCond::PERIODIC, + ddc::BoundCond::PERIODIC, + ddc::SplineSolver::GINKGO, + IDimX>(interpolation_domain)), + std::runtime_error); +} + +TEST(SplineBuilder, CorrectInterpolationGrid) +{ + CoordX const x0(0.); + CoordX const xN(1.); + std::size_t const ncells = 5; + + ddc::init_discrete_space(x0, xN, ncells); + + std::vector const range {0.05, 0.15, 0.5, 0.85, 0.95}; + + ddc::DiscreteDomain const interpolation_domain + = ddc::init_discrete_space(IDimX::init(range)); + + EXPECT_NO_THROW((ddc::SplineBuilder< + execution_space, + memory_space, + BSplinesX, + IDimX, + ddc::BoundCond::PERIODIC, + ddc::BoundCond::PERIODIC, + ddc::SplineSolver::GINKGO, + IDimX>(interpolation_domain))); +}