inv_jacobian_o_point.hpp

// SPDX-License-Identifier: MIT
#pragma once
 
#include <sll/view.hpp>
 
#include "cartesian_to_circular.hpp"
#include "circular_to_cartesian.hpp"
#include "czarny_to_cartesian.hpp"
#include "discrete_to_cartesian.hpp"
#include "mapping_tools.hpp"
 
template <class Mapping, class CoordRTheta>
class InvJacobianOPoint;
 
// Pre-declaration of CombinedMapping.
template <class Mapping1, class Mapping2>
class CombinedMapping;
 
template <class X, class Y, class R, class Theta, class Xpc, class Ypc>
class InvJacobianOPoint<
        CombinedMapping<
                CircularToCartesian<R, Theta, X, Y>,
                CartesianToCircular<Xpc, Ypc, R, Theta>>,
        ddc::Coordinate<R, Theta>>
{
    using CoordRTheta = ddc::Coordinate<R, Theta>;
 
private:
    using Mapping = CombinedMapping<
            CircularToCartesian<R, Theta, X, Y>,
            CartesianToCircular<Xpc, Ypc, R, Theta>>;
 
private:
    Mapping m_mapping;
 
public:
    KOKKOS_FUNCTION explicit InvJacobianOPoint(Mapping const& mapping) : m_mapping(mapping) {}
 
    KOKKOS_INLINE_FUNCTION Matrix_2x2 operator()() const
    {
        Matrix_2x2 J;
        J[0][0] = 1.;
        J[0][1] = 0.;
        J[1][0] = 0.;
        J[1][1] = 1.;
        return J;
    }
};
 
template <class X, class Y, class R, class Theta, class Xpc, class Ypc>
class InvJacobianOPoint<
        CombinedMapping<CzarnyToCartesian<R, Theta, X, Y>, CartesianToCircular<Xpc, Ypc, R, Theta>>,
        ddc::Coordinate<R, Theta>>
{
    using CoordRTheta = ddc::Coordinate<R, Theta>;
 
private:
    using Mapping = CombinedMapping<
            CzarnyToCartesian<R, Theta, X, Y>,
            CartesianToCircular<Xpc, Ypc, R, Theta>>;
 
private:
    Mapping m_mapping;
 
public:
    KOKKOS_FUNCTION explicit InvJacobianOPoint(Mapping const& mapping) : m_mapping(mapping) {}
 
    KOKKOS_INLINE_FUNCTION Matrix_2x2 operator()() const
    {
        Matrix_2x2 J;
        const double epsilon
                = m_mapping.template get<CzarnyToCartesian<R, Theta, X, Y>>().epsilon();
        const double e = m_mapping.template get<CzarnyToCartesian<R, Theta, X, Y>>().e();
        const double xi = Kokkos::sqrt(1. / (1. - epsilon * epsilon * 0.25));
        const double sqrt_eps_2 = Kokkos::sqrt(1. + epsilon * epsilon);
        J[0][0] = -sqrt_eps_2;
        J[0][1] = 0.;
        J[1][0] = 0.;
        J[1][1] = (2 - sqrt_eps_2) / e / xi;
        return J;
    }
};
 
template <
        class X,
        class Y,
        class SplineEvaluator,
        class R,
        class Theta,
        class MemorySpace,
        class Xpc,
        class Ypc>
class InvJacobianOPoint<
        CombinedMapping<
                DiscreteToCartesian<X, Y, SplineEvaluator, R, Theta, MemorySpace>,
                CartesianToCircular<Xpc, Ypc, R, Theta>>,
        ddc::Coordinate<R, Theta>>
{
    using CoordRTheta = ddc::Coordinate<R, Theta>;
 
private:
    using Mapping = CombinedMapping<
            DiscreteToCartesian<X, Y, SplineEvaluator, R, Theta, MemorySpace>,
            CartesianToCircular<Xpc, Ypc, R, Theta>>;
 
    using IdxRangeR = typename SplineEvaluator::evaluation_domain_type1;
    using IdxRangeTheta = typename SplineEvaluator::evaluation_domain_type2;
    using IdxR = typename IdxRangeR::discrete_element_type;
    using IdxTheta = typename IdxRangeTheta::discrete_element_type;
 
private:
    Mapping m_mapping;
 
public:
    KOKKOS_FUNCTION explicit InvJacobianOPoint(Mapping const& mapping) : m_mapping(mapping) {}
 
    KOKKOS_FUNCTION Matrix_2x2 operator()() const
    {
        DiscreteToCartesian<X, Y, SplineEvaluator, R, Theta, MemorySpace> const& discrete_mapping
                = m_mapping.template get<
                        DiscreteToCartesian<X, Y, SplineEvaluator, R, Theta, MemorySpace>>();
        Matrix_2x2 J;
        J[0][0] = 0;
        J[0][1] = 0;
        J[1][0] = 0;
        J[1][1] = 0;
        ddc::DiscreteDomain idx_range_singular_point = discrete_mapping.idx_range_singular_point();
        // Average the values at (r = 0, theta):
        IdxR ir(idx_range_singular_point.front());
        for (IdxTheta ip : IdxRangeTheta(idx_range_singular_point)) {
            ddc::Coordinate<R, Theta> coord(ddc::coordinate(ir), ddc::coordinate(ip));
            Matrix_2x2 J_first_order = discrete_mapping.first_order_jacobian_matrix_r_rtheta(coord);
 
            double th = ddc::get<Theta>(coord);
            Matrix_2x2 J_circ_r_rtheta;
            J_circ_r_rtheta[0][0] = Kokkos::cos(th);
            J_circ_r_rtheta[0][1] = Kokkos::sin(th);
            J_circ_r_rtheta[1][0] = -Kokkos::sin(th);
            J_circ_r_rtheta[1][1] = Kokkos::cos(th);
 
            Matrix_2x2 J_theta;
            J_theta = inverse(mat_mul(J_first_order, J_circ_r_rtheta));
 
            J[0][0] += J_theta[0][0];
            J[0][1] += J_theta[0][1];
            J[1][0] += J_theta[1][0];
            J[1][1] += J_theta[1][1];
        }
 
        int const theta_size = idx_range_singular_point.size();
        J[0][0] /= theta_size;
        J[0][1] /= theta_size;
        J[1][0] /= theta_size;
        J[1][1] /= theta_size;
        return J;
    }
};