Skip to content
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.

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

Add functions to check if the provided interpolation points are valid #686

Merged
merged 10 commits into from
Dec 4, 2024
Merged
4 changes: 2 additions & 2 deletions include/ddc/kernels/splines/bsplines_uniform.hpp
Original file line number Diff line number Diff line change
@@ -401,9 +401,9 @@ KOKKOS_INLINE_FUNCTION ddc::DiscreteElement<DDim> UniformBSplines<CDim, D>::
Impl<DDim, MemorySpace>::eval_basis(
DSpan1D values,
ddc::Coordinate<CDim> 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;
107 changes: 100 additions & 7 deletions include/ddc/kernels/splines/spline_builder.hpp
Copy link
Member

Choose a reason for hiding this comment

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

Sorry but I don't get the logic of the check. Are we checking the number of points in each cell ? If so shouldn't current_cell_end_idx increment in the loop as soon as point > ddc::coordinate(current_cell_end_idx) is true ?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Hmm, I agree. Either I missed something while copy/pasting or the error output isn't correct in Gysela either. I don't see where current_cell_idx is updated

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Fixed 26840d5

Original file line number Diff line number Diff line change
@@ -165,15 +165,15 @@ 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

// interpolator specific
std::unique_ptr<ddc::detail::SplinesLinearProblem<exec_space>> 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<std::size_t> cols_per_chunk = std::nullopt,
std::optional<unsigned int> preconditioner_max_block_size = std::nullopt)
: m_batched_interpolation_domain(batched_interpolation_domain)
, m_offset(compute_offset(interpolation_domain()))
, m_dx((ddc::discrete_space<BSplines>().rmax() - ddc::discrete_space<BSplines>().rmin())
/ ddc::discrete_space<BSplines>().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<unsigned int> preconditioner_max_block_size = std::nullopt);

void build_matrix_system();

void check_valid_grid();

template <class KnotElement>
static void check_n_points_in_cell(int n_points_in_cell, KnotElement current_cell_end_idx);
};

template <
@@ -478,17 +485,17 @@ template <
ddc::BoundCond BcUpper,
SplineSolver Solver,
class... IDimX>
int SplineBuilder<
void SplineBuilder<
ExecSpace,
MemorySpace,
BSplines,
InterpolationDDim,
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<double, bsplines_type::degree() + 1> 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 <class KnotElement>
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<BSplines>().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<BSplines>().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<BSplines>().break_point_domain().front() + 1;
ddc::for_each(interpolation_domain(), [&](auto idx) {
ddc::Coordinate<continuous_dimension_type> 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
8 changes: 7 additions & 1 deletion tests/splines/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -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)
139 changes: 139 additions & 0 deletions tests/splines/spline_builder.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,139 @@
// Copyright (C) The DDC development team, see COPYRIGHT.md file
//
// SPDX-License-Identifier: MIT

#include <cstddef>
#include <stdexcept>
#include <vector>

#include <ddc/ddc.hpp>
#include <ddc/kernels/splines.hpp>

#include <gtest/gtest.h>

#include <Kokkos_Core.hpp>

struct DimX
{
static constexpr bool PERIODIC = true;
};

using CoordX = ddc::Coordinate<DimX>;

static constexpr std::size_t s_degree_x = 2;

struct BSplinesX : ddc::UniformBSplines<DimX, s_degree_x>
{
};

struct IDimX : ddc::NonUniformPointSampling<DimX>
{
};

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<BSplinesX>(x0, xN, ncells);

// One point missing
std::vector<double> const range {0.1, 0.3, 0.5, 0.7};

ddc::DiscreteDomain<IDimX> const interpolation_domain
= ddc::init_discrete_space<IDimX>(IDimX::init<IDimX>(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<BSplinesX>(x0, xN, ncells);

// One point too much
std::vector<double> const range {0.1, 0.3, 0.5, 0.7, 0.9, 0.95};

ddc::DiscreteDomain<IDimX> const interpolation_domain
= ddc::init_discrete_space<IDimX>(IDimX::init<IDimX>(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<BSplinesX>(x0, xN, ncells);

// All points end up in the first cell ]0, 0.2[
std::vector<double> const range {0.1, 0.11, 0.12, 0.13, 0.14};

ddc::DiscreteDomain<IDimX> const interpolation_domain
= ddc::init_discrete_space<IDimX>(IDimX::init<IDimX>(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<BSplinesX>(x0, xN, ncells);

std::vector<double> const range {0.05, 0.15, 0.5, 0.85, 0.95};

ddc::DiscreteDomain<IDimX> const interpolation_domain
= ddc::init_discrete_space<IDimX>(IDimX::init<IDimX>(range));

EXPECT_NO_THROW((ddc::SplineBuilder<
execution_space,
memory_space,
BSplinesX,
IDimX,
ddc::BoundCond::PERIODIC,
ddc::BoundCond::PERIODIC,
ddc::SplineSolver::GINKGO,
IDimX>(interpolation_domain)));
}