How to manage tensors

[EmilyBourne]

I am not sure if this is out of scope for DDC or not.
I need to add tensors to my code (e.g. metric tensor, gradient vector). Doing this cleanly requires several of the concepts that DDC introduces. However this work is not about discontinuous dimensions hence why it may be out of scope.

I think this could be managed with a Tensor type. E.g:

template <class ElementType, class... ValidIndexSet>
class Tensor
{
    static_assert(std::conjunction_v<ddc::is_type_seq<ValidIndexSet>...>);

private:
    using AllIndexSets = ddc::TypeSeq<ValidIndexSet...>;

private:
    static constexpr std::size_t n_elements = (ddc::type_seq_size_v<ValidIndexSet> * ...);
    std::array<ElementType, n_elements> data;

public:
    KOKKOS_FUNCTION Tensor();

    template<class... QueryIndexTag>
    KOKKOS_FUNCTION ElementType& get() {
        static_assert(sizeof...(QueryIndexTag) == sizeof...(ValidIndexSet));
        ...
    }

    template<class... QueryIndexTag>
    KOKKOS_FUNCTION ElementType get() const;

    KOKKOS_FUNCTION std::size_t size() const {return n_elements;}

    static constexpr std::size_t rank() const {return sizeof...(ValidIndexSet);}

    template<std::size_t i>
    using index_set_t = type_seq_element_t<i, AllIndexSets>;
};

I am happy to provide a PoC for this if you think it is in scope.

It may be helpful to bring TypeSeq out of detail/ in order to implement this.

In the places where I would use this type I would also be interested in typing my dimensions depending on whether they are covariant/contravariant. This could probably be avoided in DDC if it is out of scope but I mention it already as it will play a much more important role in any implementation of a tensor_mul function.

E.g. a dot product may be implemented like this:

struct R_cov;
struct Theta_cov;

struct R
{
    static bool constexpr PERIODIC = false;
    static bool constexpr IS_COVARIANT = false;
    static bool constexpr IS_CONTRAVARIANT = true;
    using Dual = R_cov;
};
struct Theta
{
    static bool constexpr PERIODIC = true;
    static bool constexpr IS_COVARIANT = false;
    static bool constexpr IS_CONTRAVARIANT = true;
    using Dual = Theta_cov;
};

struct R_cov
{
    static bool constexpr PERIODIC = false;
    static bool constexpr IS_COVARIANT = true;
    static bool constexpr IS_CONTRAVARIANT = false;
    using Dual = R;
};
struct Theta_cov
{
    static bool constexpr PERIODIC = true;
    static bool constexpr IS_COVARIANT = true;
    static bool constexpr IS_CONTRAVARIANT = false;
    using Dual = Theta;
};

double dot_product(Tensor<double, ddc::TypeSeq<R_cov, Theta_cov>> a, Tensor<double, ddc::TypeSeq<R, Theta>> b) {
    return a.get<R_cov>()*b.get<R>() + a.get<Theta_cov>()*b.get<Theta>();
}

Therefore if co/contravariant concepts are not added along with tensor concepts any tensor_mul functions may be unusable.
This addition wouldn't necessarily force all users to care about co/contravariance.
For codes that don't care about co/contravariance (e.g. because they use cartesian coordinates) it is possible to write:

struct X
{
    static bool constexpr PERIODIC = false;
    static bool constexpr IS_COVARIANT = true;
    static bool constexpr IS_CONTRAVARIANT = true;
    using Dual = X;
};

[EmilyBourne]

cc: @tpadioleau @jbigot

[tpadioleau]

What kind of operation you want to support ?

[EmilyBourne]

The Tensor class would simply be a storage type. I think I put in the description all the methods that I think would be useful.

Beyond that some kind of tensor_mul function would be very useful. E.g. to go from a contravariant vector to a covariant vector.

However the exact form of this method will require some thought. It would be fantastic to be able to do something like:

tensor_mul<"i", "ij", "j">(cov_vec, metric_tensor, contra_vec);

but this is not trivial to implement.

As I mentioned in the description though, if you don't want to include covariant/contravariant concepts then I am not sure it is worth implementing this funciton.

The Tensor class is still useful with the tensor_mul function as dedicated routines can be written in user code.

[EmilyBourne]

Maybe some operators would also be useful. E.g. Tensor::operator*(double) Tensor::operator+(Tensor)

[blegouix]

I worked a lot on this topic in my project https://blegouix.github.io/similie/tensor_module.html

[tpadioleau]

Thanks for the reference, I let @EmilyBourne see if it can help.

[EmilyBourne]

Thanks for the link @blegouix . I had a look, it looks like a useful way to describe Chunks of tensors. I didn't understand how you handle co/contravariance though. Are you assuming that everything is stored in contravariant form and that the metric tensor is available?

I noticed that you chose the GPL licence. This is more restrictive than the MIT licence. I don't know if this means that we can't include things from your library.

It looks like we have gone in similar directions on some points. I think this could be useful if you want to incorporate it into DDC. For Gyselalib I think it is missing some tools and is more complex than needed on other points so I think I will keep going with my implementation for now.

[blegouix]

Maybe take a look at https://github.com/blegouix/similie/blob/ee2b3034e5ac270519cd9fd98a190468733c1c76/examples/2d_vector_laplacian.cpp#L121 or in the tests for practical usage. Yes the metric is available and used to lower or upper indices. I know this is almost undocumented but you can pm me on the DDC slack to ask me anything you want.

As you may know I have to indefinitely stall my project and I cannot do any contribution in it or DDC until at least several months. And I cannot reasonably consider working on such a big topic in DDC directly without a position in the team... I'll relax the license to MIT though so feel free to extract part of the code if you want.