init_ghosted gives rounding errors on domain boundary

[EmilyBourne]

When using init_ghosted on a uniform point sampling the addition of the ghosts leads to rounding errors on the domain boundary. This issue can be avoided using tolerances as was done in #519 so it may not be serious, but it should at least be documented so users are aware of this potential pitfall.

Minimum reprocible example

#include <iostream>

#include <ddc/ddc.hpp>

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

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

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

    ddc::Coordinate<Vx> vx_min(-8);
    ddc::Coordinate<Vx> vx_max(8);
    ddc::DiscreteVector<DDimVx> npoints(128);
    ddc::DiscreteVector<DDimVx> nghosts(3);

    ddc::DiscreteDomain<DDimVx> pre_ghost;
    ddc::DiscreteDomain<DDimVx> post_ghost;
    ddc::DiscreteDomain<DDimVx> main_domain;
    ddc::DiscreteDomain<DDimVx> full_domain;
    std::tie(main_domain, full_domain, pre_ghost, post_ghost) =
        ddc::init_discrete_space<DDimVx>(DDimVx::template init_ghosted<DDimVx>(
            vx_min, vx_max, npoints, nghosts, nghosts));

    std::cout << (vx_min - ddc::coordinate(main_domain.front())) << std::endl;

    return 0;
}

Output is:

(-8.88178e-16)

[tpadioleau]

It is true that you cannot expect exact equality between vx_min and ddc::coordinate(main_domain.front())/vx_xmax and ddc::coordinate(main_domain.back()). I don't think it is related to init_ghosted though, I guess init also has the same behavior. Are you relying on it in Gysela ?

[EmilyBourne]

I guess init also has the same behavior

I think for init you cannot expect exact equality between vx_xmax and ddc::coordinate(main_domain.back()) but you can expect exact equality between vx_xmin and ddc::coordinate(main_domain.front())

We are not relying on this for Gysela, but in DDC you were relying on it here:

ddc/include/ddc/kernels/splines/bsplines_uniform.hpp

Lines 557 to 558 in 0da9a03

	assert(x - rmin() >= -length() * 1e-14);
	assert(rmax() - x >= -length() * 1e-14);

if the extrapolation condition for a spline evaluator evaluates the value at the boundary.

As I say in the issue

it may not be serious, but it should at least be documented so users are aware of this potential pitfall.