Provide a method to combine types such as DiscreteDomains

[EmilyBourne]

It is quite common to want to combine 2 existing DiscreteDomain types into an extended DiscreteDomain. E.g. from DiscreteDomain<GridBatch1, GridBatch2> and DiscreteDomain<OperatorDim1, OperatorDim2> I would like to obtain a DiscreteDomain<GridBatch1, GridBatch2, OperatorDim1, OperatorDim2>. Currently to achieve this we need to transit via TypeSeq objects:

using BatchedDomain = ddc::detail::convert_type_seq_to_discrete_domain_t<type_seq_cat_t<ddc::to_type_seq_t<BatchDomain>, ddc::to_type_seq_t<OperatorDomain>>>;

It would be nice to have an operator which can do this. This could be implemented as above or more generally

Example implementation

template <class... Containers>
struct Combine;

template <class Container>
struct Combine<Container>
{
    using type = Container;
};

template <
        template <typename...>
        class Container,
        class... Tags,
        class... OTags,
        class... TailContainers>
struct Combine<Container<Tags...>, Container<OTags...>, TailContainers...>
{
    using type = Combine<Container<Tags..., OTags...>, TailContainers...>::type;
};

/**
 * @brief A helper structure to determine the type when combining tags across containers.
 *
 * E.g. combine_t<Idx<Tag1>, Idx<Tag2>> == Idx<Tag1, Tag2>
 *
 * @tparam Containers The containers that should be combined. They should differ only in the tags.
 */
template <class... Containers>
using combine_t = typename detail::Combine<Containers...>::type;

The advantage of the general method is that it can also be used for DiscreteElement and DiscreteVector. It is rare to need to combine these without also needing to combine a DiscreteDomain (from which the element and vector can be extracted) but this does happen.

[tpadioleau]

We do have ddc::cartesian_prod_t that works only on DiscreteDomain. From the DiscreteDomain you can retrieve discrete_element_type and discrete_vector_type. Is it enough ?

[EmilyBourne]

I am using this here:
https://github.com/gyselax/gyselalibxx/blob/0abf78de106fa6fca755a00f9b6fded1f16d235e/src/data_types/derivative_field.hpp#L366-L373

So the code:

    template <class... DElem>
    KOKKOS_FUNCTION constexpr reference operator()(DElem... elems) const noexcept
    {
        static_assert((ddc::is_discrete_element_v<DElem> && ...));
        using full_index_type = detail::combine_t<DElem...>;
        full_index_type elem(elems...);
        return base_type::get_internal_field(elem)();
    }

would become:

    template <class... Grid1D>
    KOKKOS_FUNCTION constexpr reference operator()(ddc::DiscreteElement<Grid1D>... elems) const noexcept
    {
        using full_index_range_type = ddc::cartesian_prod_t<ddc::DiscreteDomain<Grid1D>...>;
        using full_index_type = typename full_index_range_type::discrete_element_type;
        full_index_type elem(elems...);
        return base_type::get_internal_field(elem)();
    }

?
This is probably enough. But is there a reason to not let ddc::cartesian_prod_t work for elements and vectors?

[tpadioleau]

Do you need to keep the ordering of labels ? This piece of code could become

template <class... DElem>
KOKKOS_FUNCTION constexpr reference operator()(DElem... elems) const noexcept
{
    static_assert((ddc::is_discrete_element_v<DElem> && ...));
    return base_type::get_internal_field(discrete_element_type(elems...))();
}

But is there a reason to not let ddc::cartesian_prod_t work for elements and vectors?

I think we only needed it for DiscreteDomain when it was first introduced. I don't feel this is the correct name if we also want to apply it on DiscreteElement and DiscreteVector, I think the name refers to sets rather than elements of sets.

[EmilyBourne]

Do you need to keep the ordering of labels ? This piece of code could become

template <class... DElem>
KOKKOS_FUNCTION constexpr reference operator()(DElem... elems) const noexcept
{
static_assert((ddc::is_discrete_element_v && ...));
return base_type::get_internal_field(discrete_element_type(elems...))();
}

But is there a reason to not let ddc::cartesian_prod_t work for elements and vectors?

How would discrete_element_type be defined? You seem to have removed my calculation of this type without putting anything in to replace it

[EmilyBourne]

It can probably be built from my get_internal_field though

[tpadioleau]

How would discrete_element_type be defined? You seem to have removed my calculation of this type without putting anything in to replace it

We can take the one defined in DerivField

[EmilyBourne]

Actually the docstring of this method says:

 * @brief Get a modifiable reference to an element from a constant field.
* A Idx describes the element of interest. If information about the derivatives is
* missing then it is assumed that the 0-th order derivative is requested.     

So there are 2 possible definitions of discrete_element_type depending on what is provided to the method

[tpadioleau]

What about creating a discrete_element_type from elems... without the derivative and one DiscreteElement with the derivative ?

[EmilyBourne]

I'm not sure how that helps?

Suppose a 2D DerivativeField defined on (x,y), we have 4 possible types that we want for discrete_element_type:

1. DiscreteElement<GridX, GridY>
2. DiscreteElement<GridX, GridY, ddc::Deriv<X>>
3. DiscreteElement<GridX, GridY, ddc::Deriv<Y>>
4. DiscreteElement<GridX, GridY, ddc::Deriv<X>, ddc::Deriv<Y>>

The method get_internal_field contains if constexpr blocks to know which internal field is of interest depending on which of these 4 types is given:

1. Implies values internal_array[0]
2. Implies x derivatives internal_array[1]
3. Implies y derivatives internal_array[2]
4. Implies cross derivatives internal_array[3]

[EmilyBourne]

If I create 1 without and 1 with then one of the 2 will fail to compile. If you mean create multiple methods then this is not possible as the number of methods required depends on the number of dimensions (1D = 2 methods, 2D = 4 methods, 3D = 8 methods, etc)

[tpadioleau]

Ok what I had in mind does not work, I forgot you have multiple derivatives. Then maybe I would try a get_deriv_or ? Looping over the expected derivatives, if it is present take it else return a default value ? Then you can construct two separate DiscreteElement, one for the derivatives and one for the rest ?

[EmilyBourne]

if it is present take it else return a default value ?

Again I don't think this would work. The final type must not be DiscreteElement<GridX, GridY, ddc::Deriv<X>, ddc::Deriv<Y>> as the type is used to determine which internal_array. Each set of derivatives has a different domain so the internal arrays are:

• Chunk<DiscreteDomain<GridX, GridY>> (full domain)
• Chunk<DiscreteDomain<GridX, GridY, ddc::Deriv<X>>> (y domain times first and last elements from the x domain)
• Chunk<DiscreteDomain<GridX, GridY, ddc::Deriv<Y>>> (x domain times first and last elements from the y domain)
• Chunk<DiscreteDomain<GridX, GridY, ddc::Deriv<X>, ddc::Deriv<Y>>> (corners of the x,y domain)

Then you can construct two separate DiscreteElement, one for the derivatives and one for the rest ?

Again it is not two separate DiscreteElements but 2^N separate DiscreteElements with N a template parameter. Even if I managed to create these 2^N DiscreteElements, there would still be the question of identifying which of these is relevant to retrieve the internal_array

[tpadioleau]

Sorry I am failing to understand why we are not using a single DiscreteDomain<GridX, GridY, ddc::Deriv<X>, ddc::Deriv<Y>> to represent them all ? We would need a convention that the first ddc::DiscreteElement<ddc::Deriv<X>, ddc::Deriv<Y>> is the function, what you called Chunk<DiscreteDomain<GridX, GridY>>.

[EmilyBourne]

There are a few constraints to consider:

1. The 4 chunks have different domains on gridx,gridy so they cannot be saved in the same Chunk
2. Externally it would be strange to get back a Chunk<DiscreteDomain<GridX, GridY, ddc::Deriv<X>, ddc::Deriv<Y>>> to represent the values of a function at every point on a grid. This object would always need to be sliced with DiscreteElement<ddc::Deriv<X>, ddc::Deriv<Y>>(0,0) in order to retrieve the Chunk<DiscreteDomain<GridX, GridY>> which is expected by all operators which may work on this object.
3. Given 2 we need to be able to determine at compile time which of the 4 chunks should be returned
4. It is not possible to determine at compile time if the domain of a Chunk<DiscreteDomain<GridX, GridY, ddc::Deriv<X>, ddc::Deriv<Y>>> is an extension of a DiscreteDomain<ddc::Deriv<X>, ddc::Deriv<Y>> covering only the point DiscreteElement<ddc::Deriv<X>, ddc::Deriv<Y>>(0,0)