I suggest that we try an implementation of the for_each algorithm independent of the DiscreteDomain and DiscreteVector classes. I would prefer an implementation that should work for any class that implements the concept of a cartesian product, that remains to be defined.

This also applies to transform_reduce and deep_copy.

Should we have different public headers for the library, with a naming convention, or only a single header like #include <ddc>

See https://gitlab.maisondelasimulation.fr/gysela-developpers/voicexx/-/issues/40

We currently identify spline basis BSplines by a tag that is derived from the "real"/points discrete domain.

A spline basis does in fact depend on the domain rather than the dimension. We should see how we want to take that into account.

See https://gitlab.maisondelasimulation.fr/gysela-developpers/voicexx/-/issues/70

Could we take an AccessorPolicy template parameter in ChunkCommon that we would forward to internal_mdspan_type & allocation_mdspan_type or do we have specific constraints ton the accessor policy we support?

From https://gitlab.maisondelasimulation.fr/gysela-developpers/voicexx/-/issues/90

There is cmakelang that brings cmake-format that can help us to format the CMakeLists.txt. I think we should use it.

The link to cmake_format

• namespace ddc
• document the concepts
• architectures support: NVIDIA and AMD GPUs, CPUs
• respect naming conventions

In our parallel loops, it would be great to offer a way to the user to execute a per-thread preambule

From https://gitlab.maisondelasimulation.fr/gysela-developpers/voicexx/-/issues/71

ChunkCommon stores redundant information about the domain size in m_domain + m_internal_mdspan. We should see whether the advantages in terms of access performance overcome the drawback of increased copy cost each time we pass it as parameter.

I think it would be useful to allow expressions such as:

ddc::DiscreteElement<Dim1, Dim2> old_pos = ...;
ddc::DiscreteElement<Dim1, Dim2> new_pos = old_pos + ddc::DiscreteVector<Dim2>(3);

Currently it is only possible to do:

ddc::DiscreteElement<Dim1, Dim2> old_pos = ...;
ddc::DiscreteElement<Dim1, Dim2> new_pos = old_pos + ddc::DiscreteVector<Dim1, Dim2>(3, 0);

This addition would be useful in cases where we know that Dim2 is one of the dimensions but we don't want to specialise the function for all possible values of Dim1 (present or not present)

A step of a UniformPointSampling is a Coordinate. I think this is a mistake. A step is a length, but a Coordinate is a position. The return type should be double

From https://gitlab.maisondelasimulation.fr/gysela-developpers/voicexx/-/issues/61

We should run the bench and check that deepcopy perf is close to memcpy.

google benchmark can export a json file we can process in python for example.

See https://gitlab.maisondelasimulation.fr/gysela-developpers/voicexx/-/issues/20

One can overcome this issue by making a new continuous dimension. The drawback is that Coordinates do not belong to the same continuous space and need to be manually cast.

See https://gitlab.maisondelasimulation.fr/gysela-developpers/voicexx/-/issues/34

Currently and according to Voice++#26 , we only discretize points (vertices).
I would like to introduce edges discretization also.

For example, in 3D, this could look like so:

using IDimX = ×DiscretizationOf<RDimX>;
using IDimY = ×DiscretizationOf<RDimY>;
using IDimZ = ×DiscretizationOf<RDimZ>;

using IVerticesSpace = DiscreteSpace<IDimX,        IDimY,        IDimZ>;
using IXEdgesSpace   = DiscreteSpace<IDimX::Edges, IDimY,        IDimZ>;
using IYEdgesSpace   = DiscreteSpace<IDimX,        IDimY::Edges, IDimZ>;
using IXFacesSpace   = DiscreteSpace<IDimX,        IDimY::Edges, IDimZ::Edges>;
using IYFacesSpace   = DiscreteSpace<IDimX::Edges, IDimY,        IDimZ::Edges>;
using ICellsSpace    = DiscreteSpace<IDimX::Edges, IDimY::Edges, IDimZ::Edges>;

The real type would be a tensor product that would combine either Point<RDim> or Interval<RDim>

using RVertex = VectorElement<Point<RDimX>,    Point<RDimY>,    Point<RDimZ>>;
using RXEdge  = VectorElement<Interval<RDimX>, Point<RDimY>,    Point<RDimZ>>;
using RYEdge  = VectorElement<Point<RDimX>,    Interval<RDimY>, Point<RDimZ>>;
using RXFace  = VectorElement<Point<RDimX>,    Interval<RDimY>, Interval<RDimZ>>;
using RYFace  = VectorElement<Interval<RDimX>, Point<RDimY>,    Interval<RDimZ>>;
using RCell   = VectorElement<Interval<RDimX>, Interval<RDimY>, Interval<RDimZ>>;

The to_real() function would yield a type that depends on the discrete space:

IVerticesSpace vertices_space = /*[...]*/;
ICellsSpace cell_space = /*[...]*/;
RVertex = vertices_space.to_real(/*[...]*/);
RCell = cell_space.to_real(/*[...]*/);

We should offer a way to provide a scratch buffer in our parallel loops that would support efficient access when all iterations access the buffer with the same pattern (one buffer/thread on the right NUMA node for OMP, interleaved buffers on GPU, etc.)
I believe this can be built on top of #56 & #57 .

We currently use global variables to store some data related to discrete spaces. I do not see any implementation of a device global variable for now, only a discussion about it here, see also there.

As for now, discrete dimensions are represented by the discrete spaces themselves. For example this forces the definition of two classes to differentiate nodes from cells in a 1d mesh.
We want to decouple the two in order to be able to ease the definition of multiple discrete dimensions.

The file is full of untested templates that might not even compile when used

FFT branch.
Thomas suspect it is due to the compiler. I should check with another one to confirm that once the CPU version is available.

In the context of solving the Poisson equation one needs to identify the null frequency as a particular case. For now the user has to assume the null frequency is on the first position of the fourier domain. I think it can be improved.

With the recent possibility for mdspan to choose the index_type, I think we should a signed integer type common for DiscreteElementType, std::extents and iterators difference_type that is to say std::ptrdiff_t.

We should upgrade to jammy and lunar Ubuntu docker images.

          I don't know if this can be improved

Originally posted by @tpadioleau in #123 (comment)

Forked branches can't create pull requests which pass tests because they don't have the rights to push the docker.

This could be avoided if the docker is packed into an artifact instead of pushing it

Depending on the choice of normalization, the norm_coef can be 1. In these cases we could optimize out the loop normalizing the data.

It is an experiment. If it goes worse doing that we will revert to the current state.

See https://gitlab.maisondelasimulation.fr/gysela-developpers/voicexx/-/issues/67

For member accessors, we've gone with the no-get naming convention: data.value() rather than data.get_value().

For free functions though, it's not clear. We have for example auto get_domain(ChunkType const& chunk)

We should settle on a convention and apply it to the whole of DDC.

Our types in classes haven't been renamed in the big !73 renaming. For example, we have stuff like:

template <class... DDims>
class DiscreteDomain
{
public:
    using mcoord_type = DiscreteCoordinate<DDims...>;
};

This should read more like:

template <class... DDims>
class DiscreteDomain
{
public:
    using discrete_coordinate_type = DiscreteCoordinate<DDims...>;
};

I suggest to introduce a ddc namespace along with using namespace ddc; in the <ddc/ddc.hpp> header to not break application codes for now.

It would be nice to offer a multi-dimensional allocator that would return potentially non-contiguous data based on a requested domain

In the following code:

DiscreteDomain<BSplines> const dom_bsplines_singular(
            DiscreteElement<BSplines>(0),
            DiscreteVector<BSplines>(BSplines::n_singular_basis()));
Chunk<double, DiscreteDomain<BSplines>> singular_spline_coef(dom_bsplines_singular);
for_each(singular_spline_coef.domain(), [&](DiscreteElement<BSplines> const i) {
    singular_spline_coef(i) = 1.0;
});

where n_singular_basis() returns 0, the loop is entered resulting in a seg fault at the line singular_spline_coef(i) = 1.0;

See https://gitlab.maisondelasimulation.fr/gysela-developpers/voicexx/-/issues/69

We don't fully document our API currently. This should be improved

Cf. https://gitlab.maisondelasimulation.fr/gysela-developpers/voicexx/-/issues/23

I can do:

Coordinate<MyDim> get_coord(DiscreteElement<UniformPointSampling<MyDim>> el) {
    return coordinate(el);
}

but not:

Coordinate<MyDim1, MyDim2> get_coord(DiscreteElement<UniformPointSampling<MyDim1>,UniformPointSampling<MyDim2>> el) {
    return coordinate(el);
}

This should be possible

We currently apply the constexpr, inline & noexcept modifiers to functions somewhat randomly. This should be reviewed. We should come up with some guidelines and apply them all-around.

Issue following issue 41 in the gitlab.

I think we need three types similar to the manipulation of pointers in C++:

• T*
• std::size_t
• std::ptrdiff_t

In our case we would have:

• DiscreteCoordinate, identifier of a element of discretization
• DiscreteSize, multidimensional size
• DiscreteVector, multidimensional displacement