cartesian_to_circular.hpp

// SPDX-License-Identifier: MIT
#pragma once
 
#include <ddc/ddc.hpp>
 
#include <sll/view.hpp>
 
#include "mapping_tools.hpp"
 
// Pre-declaration of analytical inverse
template <class R, class Theta, class X, class Y>
class CircularToCartesian;
 
template <class X, class Y, class R, class Theta>
class CartesianToCircular
{
public:
    using cartesian_tag_x = X;
    using cartesian_tag_y = Y;
    using curvilinear_tag_r = R;
    using curvilinear_tag_theta = Theta;
 
    using CoordArg = ddc::Coordinate<X, Y>;
    using CoordResult = ddc::Coordinate<R, Theta>;
 
public:
    CartesianToCircular() = default;
 
    KOKKOS_FUNCTION CartesianToCircular(CartesianToCircular const& other) {}
 
    CartesianToCircular(CartesianToCircular&& x) = default;
 
    ~CartesianToCircular() = default;
 
    CartesianToCircular& operator=(CartesianToCircular const& x) = default;
 
    CartesianToCircular& operator=(CartesianToCircular&& x) = default;
 
    KOKKOS_FUNCTION ddc::Coordinate<R, Theta> operator()(ddc::Coordinate<X, Y> const& coord) const
    {
        const double x = ddc::get<X>(coord);
        const double y = ddc::get<Y>(coord);
        const double r = Kokkos::sqrt(x * x + y * y);
        const double theta = Kokkos::atan2(y, x);
        return ddc::Coordinate<R, Theta>(r, theta);
    }
 
    KOKKOS_FUNCTION double jacobian(ddc::Coordinate<X, Y> const& coord)
    {
        const double x = ddc::get<X>(coord);
        const double y = ddc::get<Y>(coord);
        return 2. * (x * x - y * y) / Kokkos::pow(x * x + y * y, 1.5);
    }
 
    KOKKOS_FUNCTION void jacobian_matrix(ddc::Coordinate<X, Y> const& coord, Matrix_2x2& matrix)
 
    {
        const double x = ddc::get<X>(coord);
        const double y = ddc::get<Y>(coord);
        matrix[0][0] = 2 * x / Kokkos::pow(x * x + y * y, 0.5);
        matrix[0][1] = 2 * y / Kokkos::pow(x * x + y * y, 0.5);
        matrix[1][0] = -y / Kokkos::pow(x * x + y * y, 2.);
        matrix[1][1] = x / Kokkos::pow(x * x + y * y, 2.);
    }
 
    KOKKOS_FUNCTION double jacobian_11(ddc::Coordinate<X, Y> const& coord)
    {
        const double x = ddc::get<X>(coord);
        const double y = ddc::get<Y>(coord);
        return 2 * x / Kokkos::pow(x * x + y * y, 0.5);
        ;
    }
 
    KOKKOS_FUNCTION double jacobian_12(ddc::Coordinate<X, Y> const& coord)
    {
        const double x = ddc::get<X>(coord);
        const double y = ddc::get<Y>(coord);
        return 2 * y / Kokkos::pow(x * x + y * y, 0.5);
    }
 
    KOKKOS_FUNCTION double jacobian_21(ddc::Coordinate<X, Y> const& coord)
    {
        const double x = ddc::get<X>(coord);
        const double y = ddc::get<Y>(coord);
        return -y / Kokkos::pow(x * x + y * y, 2.);
    }
 
    KOKKOS_FUNCTION double jacobian_22(ddc::Coordinate<X, Y> const& coord)
    {
        const double x = ddc::get<X>(coord);
        const double y = ddc::get<Y>(coord);
        return x / Kokkos::pow(x * x + y * y, 2.);
    }
 
    CircularToCartesian<R, Theta, X, Y> get_inverse_mapping() const
    {
        return CircularToCartesian<R, Theta, X, Y>();
    }
};
 
namespace mapping_detail {
template <class X, class Y, class R, class Theta, class ExecSpace>
struct MappingAccessibility<ExecSpace, CartesianToCircular<X, Y, R, Theta>> : std::true_type
{
};
} // namespace mapping_detail