discrete_to_cartesian.hpp

// SPDX-License-Identifier: MIT
#pragma once
#include <iostream>
 
#include <ddc/ddc.hpp>
 
#include <sll/view.hpp>
 
#include "coordinate_converter.hpp"
#include "curvilinear2d_to_cartesian.hpp"
#include "jacobian.hpp"
#include "mapping_tools.hpp"
#include "pseudo_cartesian_compatible_mapping.hpp"
 
template <
        class X,
        class Y,
        class SplineEvaluator,
        class R = typename SplineEvaluator::continuous_dimension_type1,
        class Theta = typename SplineEvaluator::continuous_dimension_type2,
        class MemorySpace = typename SplineEvaluator::memory_space>
class DiscreteToCartesian
    : public CoordinateConverter<ddc::Coordinate<R, Theta>, ddc::Coordinate<X, Y>>
    , public NonAnalyticalJacobian<ddc::Coordinate<R, Theta>>
    , public PseudoCartesianCompatibleMapping
    , public Curvilinear2DToCartesian<X, Y, R, Theta>
{
    static_assert(std::is_same_v<MemorySpace, typename SplineEvaluator::memory_space>);
 
public:
    using BSplineR = typename SplineEvaluator::bsplines_type1;
    using BSplineTheta = typename SplineEvaluator::bsplines_type2;
 
    using cartesian_tag_x = typename Curvilinear2DToCartesian<X, Y, R, Theta>::cartesian_tag_x;
    using cartesian_tag_y = typename Curvilinear2DToCartesian<X, Y, R, Theta>::cartesian_tag_y;
    using curvilinear_tag_r = typename Curvilinear2DToCartesian<X, Y, R, Theta>::curvilinear_tag_r;
    using curvilinear_tag_theta =
            typename Curvilinear2DToCartesian<X, Y, R, Theta>::curvilinear_tag_theta;
 
    using CoordArg = ddc::Coordinate<R, Theta>;
    using CoordResult = ddc::Coordinate<X, Y>;
 
private:
    using spline_idx_range = ddc::DiscreteDomain<BSplineR, BSplineTheta>;
 
    using SplineType = ddc::
            ChunkView<double, spline_idx_range, std::experimental::layout_right, MemorySpace>;
 
    using IdxRangeTheta = typename SplineEvaluator::evaluation_domain_type2;
    using IdxTheta = typename IdxRangeTheta::discrete_element_type;
 
private:
    SplineType m_x_spline_representation;
    SplineType m_y_spline_representation;
    SplineEvaluator m_spline_evaluator;
    IdxRangeTheta m_idx_range_theta;
 
public:
    KOKKOS_FUNCTION DiscreteToCartesian(
            SplineType curvilinear_to_x,
            SplineType curvilinear_to_y,
            SplineEvaluator const& evaluator,
            IdxRangeTheta idx_range_theta)
        : m_x_spline_representation(curvilinear_to_x)
        , m_y_spline_representation(curvilinear_to_y)
        , m_spline_evaluator(evaluator)
        , m_idx_range_theta(idx_range_theta)
    {
    }
 
    KOKKOS_FUNCTION ddc::Coordinate<X, Y> operator()(
            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord) const final
    {
        const double x = m_spline_evaluator(coord, m_x_spline_representation.span_cview());
        const double y = m_spline_evaluator(coord, m_y_spline_representation.span_cview());
        return ddc::Coordinate<X, Y>(x, y);
    }
 
    KOKKOS_FUNCTION void jacobian_matrix(
            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord,
            Matrix_2x2& matrix) const final
    {
        matrix[0][0]
                = m_spline_evaluator.deriv_dim_1(coord, m_x_spline_representation.span_cview());
        matrix[0][1]
                = m_spline_evaluator.deriv_dim_2(coord, m_x_spline_representation.span_cview());
        matrix[1][0]
                = m_spline_evaluator.deriv_dim_1(coord, m_y_spline_representation.span_cview());
        matrix[1][1]
                = m_spline_evaluator.deriv_dim_2(coord, m_y_spline_representation.span_cview());
    }
 
    KOKKOS_FUNCTION double jacobian_11(
            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord) const final
    {
        return m_spline_evaluator.deriv_dim_1(coord, m_x_spline_representation.span_cview());
    }
 
    KOKKOS_FUNCTION double jacobian_12(
            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord) const final
    {
        return m_spline_evaluator.deriv_dim_2(coord, m_x_spline_representation.span_cview());
    }
 
    KOKKOS_FUNCTION double jacobian_21(
            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord) const final
    {
        return m_spline_evaluator.deriv_dim_1(coord, m_y_spline_representation.span_cview());
    }
 
    KOKKOS_FUNCTION double jacobian_22(
            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const& coord) const final
    {
        return m_spline_evaluator.deriv_dim_2(coord, m_y_spline_representation.span_cview());
    }
 
 
 
    KOKKOS_FUNCTION void to_pseudo_cartesian_jacobian_center_matrix(Matrix_2x2& matrix) const final
    {
        matrix[0][0] = 0;
        matrix[0][1] = 0;
        matrix[1][0] = 0;
        matrix[1][1] = 0;
 
        // Average the values at (r = 0, theta):
        for (IdxTheta ip : m_idx_range_theta) {
            const double th = ddc::coordinate(ip);
            ddc::Coordinate<curvilinear_tag_r, curvilinear_tag_theta> const coord(0, th);
            double const deriv_1_x
                    = m_spline_evaluator.deriv_dim_1(coord, m_x_spline_representation.span_cview());
            double const deriv_1_2_x
                    = m_spline_evaluator
                              .deriv_1_and_2(coord, m_x_spline_representation.span_cview());
            double const deriv_1_y
                    = m_spline_evaluator.deriv_dim_1(coord, m_y_spline_representation.span_cview());
            double const deriv_1_2_y
                    = m_spline_evaluator
                              .deriv_1_and_2(coord, m_y_spline_representation.span_cview());
 
            // Matrix from pseudo-Cart domain to physical domain by logical domain
            double const j11 = deriv_1_x * Kokkos::cos(th) - deriv_1_2_x * Kokkos::sin(th);
            double const j12 = deriv_1_x * Kokkos::sin(th) + deriv_1_2_x * Kokkos::cos(th);
            double const j21 = deriv_1_y * Kokkos::cos(th) - deriv_1_2_y * Kokkos::sin(th);
            double const j22 = deriv_1_y * Kokkos::sin(th) + deriv_1_2_y * Kokkos::cos(th);
 
            double const jacobian = j11 * j22 - j12 * j21;
            // Matrix from physical domain to pseudo_cart domain by logical domain
            assert(fabs(jacobian) >= 1e-16);
 
            matrix[0][0] += j22 / jacobian;
            matrix[0][1] += -j12 / jacobian;
            matrix[1][0] += -j21 / jacobian;
            matrix[1][1] += j11 / jacobian;
        }
 
        int const theta_size = m_idx_range_theta.size();
        matrix[0][0] /= theta_size;
        matrix[0][1] /= theta_size;
        matrix[1][0] /= theta_size;
        matrix[1][1] /= theta_size;
    }
 
 
    KOKKOS_FUNCTION double to_pseudo_cartesian_jacobian_11_center() const final
    {
        Matrix_2x2 jacobian;
        to_pseudo_cartesian_jacobian_center_matrix(jacobian);
        return jacobian[0][0];
    }
 
 
    KOKKOS_FUNCTION double to_pseudo_cartesian_jacobian_12_center() const final
    {
        Matrix_2x2 jacobian;
        to_pseudo_cartesian_jacobian_center_matrix(jacobian);
        return jacobian[0][1];
    }
 
 
    KOKKOS_FUNCTION double to_pseudo_cartesian_jacobian_21_center() const final
    {
        Matrix_2x2 jacobian;
        to_pseudo_cartesian_jacobian_center_matrix(jacobian);
        return jacobian[1][0];
    }
 
 
    KOKKOS_FUNCTION double to_pseudo_cartesian_jacobian_22_center() const final
    {
        Matrix_2x2 jacobian;
        to_pseudo_cartesian_jacobian_center_matrix(jacobian);
        return jacobian[1][1];
    }
 
    KOKKOS_INLINE_FUNCTION const ddc::Coordinate<X, Y> control_point(
            ddc::DiscreteElement<BSplineR, BSplineTheta> const& el) const
    {
        return ddc::Coordinate<X, Y>(m_x_spline_representation(el), m_y_spline_representation(el));
    }
};
 
 
namespace mapping_detail {
template <
        class X,
        class Y,
        class SplineEvaluator,
        class R,
        class Theta,
        class MemorySpace,
        class ExecSpace>
struct MappingAccessibility<
        ExecSpace,
        DiscreteToCartesian<X, Y, SplineEvaluator, R, Theta, MemorySpace>>
{
    static constexpr bool value = Kokkos::SpaceAccessibility<ExecSpace, MemorySpace>::accessible;
};
} // namespace mapping_detail