Add the possibility to compute 2nd derivatives or higher with the spline evaluators

[PaulineVidal]

Currently, the SplineEvaluator and SplineEvaluator2D operators offers methods to compute the derivative along the first dimension, (the second dimension and the cross-derivatives for 2D). However, to compute cross-derivatives involving a change of coordinates, we need the second derivatives in the two dimensions (2D case).
E.g. for a mapping $\mathcal{F}: (u,v) \mapsto (x,y)$ and a given function $f$,

$$\begin{aligned} \partial_{xy} f & = \partial_{xy}u \partial_u f + \partial_{xy}v \partial_v f \\\ & + \partial_x u \partial_y u \partial_{u^2} f \\\ & + \partial_x v \partial_y v \partial_{v^2} f \\\ & + (\partial_x u \partial_y v + \partial_x v \partial_y u) \partial_{uv} f \end{aligned}$$

• $\partial_{xy}u, \partial_{xy}v, \partial_x u, \partial_y u , \partial_x v,$ and $\partial_y v$ should be provided by the Jacobian and Hessian of the mapping,
• $\partial_u f , \partial_v f$ and $\partial_{uv} f$ are already provided by the spline evaluator 2D,
• only remain $\partial_{u^2} f$ and $\partial_{v^2} f$ .

I noticed that the computation of $n$-th derivatives is already possible in the definition of the B-splines: https://github.com/CExA-project/ddc/blob/c8aac4840109ceeede757d6f01090a3e30403338/include/ddc/kernels/splines/bsplines_non_uniform.hpp#L551C34-L551C57

I think that in https://github.com/CExA-project/ddc/blob/main/include/ddc/kernels/splines/spline_evaluator_2d.hpp the main changes are to create another tag, e.g.:

    struct eval_deriv2_type
    {
    };

and later in the code https://github.com/CExA-project/ddc/blob/main/include/ddc/kernels/splines/spline_evaluator_2d.hpp#L1271, add the new option, e.g.:

    template <class EvalType1, class EvalType2, class Layout, class... CoordsDims>
    KOKKOS_INLINE_FUNCTION double eval_no_bc(
            ddc::Coordinate<CoordsDims...> const& coord_eval,
            ddc::ChunkSpan<double const, spline_domain_type, Layout, memory_space> const
                    spline_coef) const
    {
        static_assert(
                std::is_same_v<EvalType1, eval_type> || std::is_same_v<EvalType1, eval_deriv_type>);
        static_assert(
                std::is_same_v<EvalType2, eval_type> || std::is_same_v<EvalType2, eval_deriv_type>);
        ddc::DiscreteElement<bsplines_type1> jmin1;
        ddc::DiscreteElement<bsplines_type2> jmin2;

        std::array<double, bsplines_type1::degree() + 1> vals1_ptr;
        Kokkos::mdspan<double, Kokkos::extents<std::size_t, bsplines_type1::degree() + 1>> const
                vals1(vals1_ptr.data());
        std::array<double, bsplines_type2::degree() + 1> vals2_ptr;
        Kokkos::mdspan<double, Kokkos::extents<std::size_t, bsplines_type2::degree() + 1>> const
                vals2(vals2_ptr.data());

        std::array<double, 3 * (bsplines_1::degree() + 1)> derivs1_ptr;
        ddc::DSpan2D const derivs1(derivs1_ptr.data(), bsplines_1::degree() + 1, 3);
        std::array<double, 3 * (bsplines_2::degree() + 1)> derivs2_ptr;
        ddc::DSpan2D const derivs2(derivs2_ptr.data(), bsplines_2::degree() + 1, 3);

        ddc::Coordinate<continuous_dimension_type1> const coord_eval_interest1(coord_eval);
        ddc::Coordinate<continuous_dimension_type2> const coord_eval_interest2(coord_eval);

        if constexpr (std::is_same_v<EvalType1, eval_type>) {
            jmin1 = ddc::discrete_space<bsplines_type1>().eval_basis(vals1, coord_eval_interest1);
        } else if constexpr (std::is_same_v<EvalType1, eval_deriv_type>) {
            jmin1 = ddc::discrete_space<bsplines_type1>().eval_deriv(vals1, coord_eval_interest1);
        } else if constexpr (std::is_same_v<EvalType1, eval_deriv2_type>) {
            jmin1 = ddc::discrete_space<bsplines_1>()
                            .eval_basis_and_n_derivs(derivs1, coord_eval_interest1, 2);
            std::size_t two = 2;
            for (std::size_t i = 0; i < bsplines_1::degree() + 1; ++i) {
                vals1[i] = derivs1(i, two);
            }
        }
        if constexpr (std::is_same_v<EvalType2, eval_type>) {
            jmin2 = ddc::discrete_space<bsplines_type2>().eval_basis(vals2, coord_eval_interest2);
        } else if constexpr (std::is_same_v<EvalType2, eval_deriv_type>) {
            jmin2 = ddc::discrete_space<bsplines_type2>().eval_deriv(vals2, coord_eval_interest2);
        } else if constexpr (std::is_same_v<EvalType2, eval_deriv2_type>) {
            jmin2 = ddc::discrete_space<bsplines_2>()
                            .eval_basis_and_n_derivs(derivs2, coord_eval_interest2, 2);
            std::size_t two = 2;
            for (std::size_t i = 0; i < bsplines_2::degree() + 1; ++i) {
                vals2[i] = derivs2(i, two);
            }
        }

        double y = 0.0;
        for (std::size_t i = 0; i < bsplines_type1::degree() + 1; ++i) {
            for (std::size_t j = 0; j < bsplines_type2::degree() + 1; ++j) {
                y += spline_coef(
                             ddc::DiscreteElement<
                                     bsplines_type1,
                                     bsplines_type2>(jmin1 + i, jmin2 + j))
                     * vals1[i] * vals2[j];
            }
        }
        return y;
    }

The addition of public method such as deriv_dim_1, deriv_dim_2 etc should be pretty similar for the second derivatives by changing the tag eval_deriv_type' by eval_deriv2_type'. Suggestion of name second_deriv_dim_1, second_deriv_dim_2. But maybe as it seems to be general for the B-splines, it may be possible to generalise it to the evaluator too, so put the order as an input?

[EmilyBourne]

The 2Devaluator has a second derivative:

ddc/include/ddc/kernels/splines/spline_evaluator_2d.hpp

Lines 573 to 593 in 23760f3

	/**
	* @brief Double-differentiate 2D spline function (described by its spline coefficients) at a given coordinate along specified dimensions of interest.
	*
	* The spline coefficients represent a 2D spline function defined on a B-splines (basis splines). They can be
	* obtained via various methods, such as using a SplineBuilder2D.
	*
	* Note: double-differentiation other than cross-differentiation is not supported atm. See #440
	*
	* @tparam InterestDim1 First dimension along which differentiation is performed.
	* @tparam InterestDim2 Second dimension along which differentiation is performed.
	*
	* @param coord_eval The coordinate where the spline is double-differentiated. Note that only the components along the dimensions of interest are used.
	* @param spline_coef A ChunkSpan storing the 2D spline coefficients.
	*
	* @return The derivative of the spline function at the desired coordinate.
	*/
	template <class InterestDim1, class InterestDim2, class Layout, class... CoordsDims>
	KOKKOS_FUNCTION double deriv2(
	ddc::Coordinate<CoordsDims...> const& coord_eval,
	ddc::ChunkSpan<double const, spline_domain_type, Layout, memory_space> const
	spline_coef) const

In the future I think it would be better to not have multiple differentiating functions with names like deriv, deriv2, deriv3 etc, but rather to accept an argument ddc::DiscreteElement<ddc::Deriv<Dim>>/ddc::DiscreteDomain<ddc::Deriv<Dim>> specifying which derivatives should be provided in the output