czarny_to_cartesian.hpp

// SPDX-License-Identifier: MIT
#pragma once
 
#include <cmath>
 
#include <ddc/ddc.hpp>
 
#include <sll/view.hpp>
 
#include "coordinate_converter.hpp"
#include "curvilinear2d_to_cartesian.hpp"
#include "mapping_tools.hpp"
#include "pseudo_cartesian_compatible_mapping.hpp"
 
 
template <class X, class Y, class R, class Theta>
class CzarnyToCartesian
    : public CoordinateConverter<ddc::Coordinate<X, Y>, ddc::Coordinate<R, Theta>>
    , public PseudoCartesianCompatibleMapping
    , public Curvilinear2DToCartesian<X, Y, R, Theta>
{
public:
    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:
    double m_epsilon;
    double m_e;
 
public:
    CzarnyToCartesian(double epsilon, double e) : m_epsilon(epsilon), m_e(e) {}
 
    KOKKOS_FUNCTION CzarnyToCartesian(CzarnyToCartesian const& other)
        : m_epsilon(other.epsilon())
        , m_e(other.e())
    {
    }
 
    CzarnyToCartesian(CzarnyToCartesian&& x) = default;
 
    ~CzarnyToCartesian() = default;
 
    CzarnyToCartesian& operator=(CzarnyToCartesian const& x) = default;
 
    CzarnyToCartesian& operator=(CzarnyToCartesian&& x) = default;
 
    KOKKOS_FUNCTION double epsilon() const
    {
        return m_epsilon;
    }
 
    KOKKOS_FUNCTION double e() const
    {
        return m_e;
    }
 
    KOKKOS_FUNCTION ddc::Coordinate<X, Y> operator()(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
        const double tmp1
                = Kokkos::sqrt(m_epsilon * (m_epsilon + 2.0 * r * Kokkos::cos(theta)) + 1.0);
 
        const double x = (1.0 - tmp1) / m_epsilon;
        const double y = m_e * r * Kokkos::sin(theta)
                         / (Kokkos::sqrt(1.0 - 0.25 * m_epsilon * m_epsilon) * (2.0 - tmp1));
 
        return ddc::Coordinate<X, Y>(x, y);
    }
 
    KOKKOS_FUNCTION ddc::Coordinate<R, Theta> operator()(
            ddc::Coordinate<X, Y> const& coord) const final
    {
        const double x = ddc::get<X>(coord);
        const double y = ddc::get<Y>(coord);
        const double ex = 1. + m_epsilon * x;
        const double ex2 = (m_epsilon * x * x - 2. * x - m_epsilon);
        const double xi2 = 1. / (1. - m_epsilon * m_epsilon * 0.25);
        const double xi = Kokkos::sqrt(xi2);
        const double r = Kokkos::sqrt(y * y * ex * ex / (m_e * m_e * xi2) + ex2 * ex2 * 0.25);
        double theta
                = Kokkos::atan2(2. * y * ex, (m_e * xi * (m_epsilon * x * x - 2. * x - m_epsilon)));
        if (theta < 0) {
            theta = 2 * M_PI + theta;
        }
        return ddc::Coordinate<R, Theta>(r, theta);
    }
 
    KOKKOS_FUNCTION double jacobian(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
        const double xi = Kokkos::sqrt(1. / (1. - m_epsilon * m_epsilon * 0.25));
        return -r / Kokkos::sqrt(1 + m_epsilon * (m_epsilon + 2.0 * r * Kokkos::cos(theta))) * m_e
               * xi
               / (2 - Kokkos::sqrt(1 + m_epsilon * (m_epsilon + 2.0 * r * Kokkos::cos(theta))));
    }
 
 
    KOKKOS_FUNCTION void jacobian_matrix(ddc::Coordinate<R, Theta> const& coord, Matrix_2x2& matrix)
            const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
 
        const double sin_theta = Kokkos::sin(theta);
        const double cos_theta = Kokkos::cos(theta);
        const double xi2 = 1. / (1. - m_epsilon * m_epsilon * 0.25);
        const double xi = Kokkos::sqrt(xi2);
        const double sqrt_eps = Kokkos::sqrt(m_epsilon * (m_epsilon + 2.0 * r * cos_theta) + 1.0);
        const double sqrt_eps_2 = 2.0 - sqrt_eps;
 
        matrix[0][0] = -cos_theta / sqrt_eps;
        matrix[0][1] = r * sin_theta / sqrt_eps;
        matrix[1][0] = m_e * m_epsilon * r * sin_theta * cos_theta * xi
                               / (sqrt_eps_2 * sqrt_eps_2 * sqrt_eps)
                       + m_e * sin_theta * xi / sqrt_eps_2;
        matrix[1][1] = r
                       * (-m_e * m_epsilon * r * sin_theta * sin_theta * xi
                                  / (sqrt_eps_2 * sqrt_eps_2 * sqrt_eps)
                          + m_e * cos_theta * xi / sqrt_eps_2);
    }
 
    KOKKOS_FUNCTION double jacobian_11(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
        return -Kokkos::cos(theta)
               / Kokkos::sqrt(m_epsilon * (m_epsilon + 2.0 * r * Kokkos::cos(theta)) + 1.0);
    }
 
    KOKKOS_FUNCTION double jacobian_12(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
        return r * Kokkos::sin(theta)
               / Kokkos::sqrt(m_epsilon * (m_epsilon + 2.0 * r * Kokkos::cos(theta)) + 1.0);
    }
 
    KOKKOS_FUNCTION double jacobian_21(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
 
        const double sin_theta = Kokkos::sin(theta);
        const double cos_theta = Kokkos::cos(theta);
        const double xi2 = 1. / (1. - m_epsilon * m_epsilon * 0.25);
        const double xi = Kokkos::sqrt(xi2);
        const double tmp1 = Kokkos::sqrt(m_epsilon * (m_epsilon + 2.0 * r * cos_theta) + 1.0);
        const double tmp2 = 2.0 - tmp1;
        return m_e * m_epsilon * r * sin_theta * cos_theta * xi / (tmp2 * tmp2 * tmp1)
               + m_e * sin_theta * xi / tmp2;
    }
 
    KOKKOS_FUNCTION double jacobian_22(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
 
        const double sin_theta = Kokkos::sin(theta);
        const double cos_theta = Kokkos::cos(theta);
        const double xi2 = 1. / (1. - m_epsilon * m_epsilon * 0.25);
        const double xi = Kokkos::sqrt(xi2);
        const double tmp1 = Kokkos::sqrt(m_epsilon * (m_epsilon + 2.0 * r * cos_theta) + 1.0);
        const double tmp2 = 2.0 - tmp1;
        return r
               * (-m_e * m_epsilon * r * sin_theta * sin_theta * xi / (tmp2 * tmp2 * tmp1)
                  + m_e * cos_theta * xi / tmp2);
    }
 
    KOKKOS_FUNCTION void inv_jacobian_matrix(
            ddc::Coordinate<R, Theta> const& coord,
            Matrix_2x2& matrix) const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
 
        assert(r >= 1e-15);
 
        const double sin_theta = Kokkos::sin(theta);
        const double cos_theta = Kokkos::cos(theta);
        const double xi = Kokkos::sqrt(1. / (1. - m_epsilon * m_epsilon * 0.25));
        const double divisor = 2 - Kokkos::sqrt(1 + m_epsilon * (m_epsilon + 2.0 * r * cos_theta));
 
        const double fact_1 = 1 / Kokkos::sqrt(1 + m_epsilon * (m_epsilon + 2.0 * r * cos_theta));
        const double fact_2 = m_e * m_epsilon * xi * r * sin_theta * fact_1 / divisor / divisor;
        const double fact_3 = m_e * xi / divisor;
 
        matrix[0][0] = -1 / fact_1 * (-sin_theta * fact_2 + cos_theta * fact_3) / fact_3;
        matrix[0][1] = sin_theta / fact_3;
        matrix[1][0] = 1 / r / fact_1 * (cos_theta * fact_2 + sin_theta * fact_3) / fact_3;
        matrix[1][1] = 1 / r * cos_theta / fact_3;
    }
 
 
 
    KOKKOS_FUNCTION double inv_jacobian_11(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
 
        const double sin_theta = Kokkos::sin(theta);
        const double cos_theta = Kokkos::cos(theta);
        const double xi = Kokkos::sqrt(1. / (1. - m_epsilon * m_epsilon * 0.25));
        const double divisor = 2 - Kokkos::sqrt(1 + m_epsilon * (m_epsilon + 2.0 * r * cos_theta));
 
        const double fact_1 = 1 / Kokkos::sqrt(1 + m_epsilon * (m_epsilon + 2.0 * r * cos_theta));
        const double fact_2 = m_e * m_epsilon * xi * r * sin_theta * fact_1 / divisor / divisor;
        const double fact_3 = m_e * xi / divisor;
 
        return -1 / fact_1 * (-sin_theta * fact_2 + cos_theta * fact_3) / fact_3;
    }
 
    KOKKOS_FUNCTION double inv_jacobian_12(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
 
        const double sin_theta = Kokkos::sin(theta);
        const double cos_theta = Kokkos::cos(theta);
        const double xi = Kokkos::sqrt(1. / (1. - m_epsilon * m_epsilon * 0.25));
        const double divisor = 2 - Kokkos::sqrt(1 + m_epsilon * (m_epsilon + 2.0 * r * cos_theta));
 
        const double fact_3 = m_e * xi / divisor;
        return sin_theta / fact_3;
    }
 
    KOKKOS_FUNCTION double inv_jacobian_21(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
 
        assert(r >= 1e-15);
 
        const double sin_theta = Kokkos::sin(theta);
        const double cos_theta = Kokkos::cos(theta);
        const double xi = Kokkos::sqrt(1. / (1. - m_epsilon * m_epsilon * 0.25));
        const double divisor = 2 - Kokkos::sqrt(1 + m_epsilon * (m_epsilon + 2.0 * r * cos_theta));
 
        const double fact_1 = 1 / Kokkos::sqrt(1 + m_epsilon * (m_epsilon + 2.0 * r * cos_theta));
        const double fact_2 = m_e * m_epsilon * xi * r * sin_theta * fact_1 / divisor / divisor;
        const double fact_3 = m_e * xi / divisor;
 
        return 1 / r / fact_1 * (cos_theta * fact_2 + sin_theta * fact_3) / fact_3;
    }
 
    KOKKOS_FUNCTION double inv_jacobian_22(ddc::Coordinate<R, Theta> const& coord) const
    {
        const double r = ddc::get<R>(coord);
        const double theta = ddc::get<Theta>(coord);
 
        assert(r >= 1e-15);
 
        const double cos_theta = Kokkos::cos(theta);
        const double xi = Kokkos::sqrt(1. / (1. - m_epsilon * m_epsilon * 0.25));
        const double divisor = 2 - Kokkos::sqrt(1 + m_epsilon * (m_epsilon + 2.0 * r * cos_theta));
 
        const double fact_3 = m_e * xi / divisor;
        return 1 / r * cos_theta / fact_3;
    }
 
 
    KOKKOS_FUNCTION void to_pseudo_cartesian_jacobian_center_matrix(Matrix_2x2& matrix) const final
    {
        const double xi = Kokkos::sqrt(1. / (1. - m_epsilon * m_epsilon * 0.25));
        const double sqrt_eps_2 = Kokkos::sqrt(1. + m_epsilon * m_epsilon);
        matrix[0][0] = -sqrt_eps_2;
        matrix[0][1] = 0.;
        matrix[1][0] = 0.;
        matrix[1][1] = (2 - sqrt_eps_2) / m_e / xi;
    }
 
    KOKKOS_FUNCTION double to_pseudo_cartesian_jacobian_11_center() const final
    {
        return -Kokkos::sqrt(1 + m_epsilon * m_epsilon);
    }
 
    KOKKOS_FUNCTION double to_pseudo_cartesian_jacobian_12_center() const final
    {
        return 0;
    }
 
    KOKKOS_FUNCTION double to_pseudo_cartesian_jacobian_21_center() const final
    {
        return 0;
    }
 
    KOKKOS_FUNCTION double to_pseudo_cartesian_jacobian_22_center() const final
    {
        const double xi = Kokkos::sqrt(1. / (1. - m_epsilon * m_epsilon * 0.25));
        return (2 - Kokkos::sqrt(1 + m_epsilon * m_epsilon)) / m_e / xi;
    }
};
 
namespace mapping_detail {
template <class X, class Y, class R, class Theta, class ExecSpace>
struct MappingAccessibility<ExecSpace, CzarnyToCartesian<X, Y, R, Theta>> : std::true_type
{
};
} // namespace mapping_detail