Rounding error introduced by #519

[EmilyBourne]

PR #519 introduced tolerances on the evaluation of bsplines in order to avoid issues due to rounding errors near the domain boundaries. However the PR didn't change the evaluation algorithm. As a result if the rounding error occurs near the upper bound the value of the returned index is too large. This can lead to chunks being evaluated outside their domain.

Minimum reproducer

#include <iostream>

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

struct Vx
{
    static bool constexpr PERIODIC = false;
};

struct GridVx : ddc::UniformPointSampling<Vx>
{
};

struct BSplines : ddc::UniformBSplines<Vx, 3>
{
};

int main(int argc, char** argv)
{
    Kokkos::ScopeGuard kokkos_scope(argc, argv);
    ddc::ScopeGuard ddc_scope(argc, argv);

    ddc::Coordinate<Vx> min(-6.0);
    ddc::Coordinate<Vx> max(6.0);
    ddc::DiscreteVector<GridVx> ncells(34);
    ddc::init_discrete_space<BSplines>(min, max, ncells);

    using InterpPointInitMethod = ddc::GrevilleInterpolationPoints<BSplines, ddc::BoundCond::HERMITE, ddc::BoundCond::HERMITE>;
    ddc::init_discrete_space<GridVx>(InterpPointInitMethod::get_sampling<GridVx>());
    ddc::DiscreteDomain<GridVx> idx_range(InterpPointInitMethod::get_domain<GridVx>());

    // Shows that ddc::coordinate(idx_range.back()) > max
    std::cout << ddc::coordinate(idx_range.back()) << " " << max << " "
              << (ddc::coordinate(idx_range.back()) - max) << std::endl;

    using SplineVxBuilder = ddc::SplineBuilder<
            Kokkos::DefaultExecutionSpace,
            Kokkos::DefaultExecutionSpace::memory_space,
            BSplines,
            GridVx,
            ddc::BoundCond::HERMITE,
            ddc::BoundCond::HERMITE,
            ddc::SplineSolver::LAPACK,
            GridVx>;

    SplineVxBuilder builder(idx_range);

    ddc::DiscreteDomain<BSplines> spline_idx_range = builder.spline_domain();

    std::array<double, 4> basis_vals_ptr;
    ddc::DSpan1D basis_vals(basis_vals_ptr.data(), 4);

    ddc::DiscreteElement<BSplines> idx
            = ddc::discrete_space<BSplines>().eval_basis(basis_vals, ddc::coordinate(idx_range.back()));

    std::cout << "Iteration range : " << idx << " " << idx+3 << std::endl;
    std::cout << "Last valid index : " << spline_idx_range.back() << std::endl;

    return EXIT_SUCCESS;
}

Output:

(6) (6) (3.55271e-15)
Iteration range : (34) (37)
Last valid index : (36)

The "iteration range" refers to the range of bsplines which should be referenced by the 4 bspline values in the array.
Iteration over this range occurs here:

ddc/include/ddc/kernels/splines/spline_evaluator.hpp

Lines 450 to 452 in 93b822f

	for (std::size_t i = 0; i < bsplines_type::degree() + 1; ++i) {
	y += spline_coef(ddc::DiscreteElement<bsplines_type>(jmin + i)) * vals[i];
	}

(not critical as the too large coordinate leads to extrapolation conditions being evaluated instead)

and here:

ddc/include/ddc/kernels/splines/constant_extrapolation_rule.hpp

Lines 71 to 73 in 93b822f

	for (std::size_t i = 0; i < BSplines::degree() + 1; ++i) {
	y += spline_coef(idx + i) * vals[i];
	}

(can cause unexpected failures in simulations)

cc: @Geoflow

[tpadioleau]

I think in #519 we should have also updated the branching in get_icell_and_offset to take into account the new tolerance.