advection_domain.hpp

// SPDX-License-Identifier: MIT
#pragma once
#include <cassert>
#include <typeinfo>
 
#include <sll/mapping/circular_to_cartesian.hpp>
#include <sll/mapping/curvilinear2d_to_cartesian.hpp>
 
#include "advection_domain.hpp"
#include "ddc_alias_inline_functions.hpp"
#include "ddc_aliases.hpp"
#include "ddc_helper.hpp"
#include "directional_tag.hpp"
#include "geometry.hpp"
#include "geometry_pseudo_cartesian.hpp"
#include "math_tools.hpp"
#include "vector_field.hpp"
#include "vector_field_mem.hpp"
 
template <class Mapping>
class AdvectionDomain
{
public:
    virtual ~AdvectionDomain() {};
};
 
 
template <class Mapping>
class AdvectionPhysicalDomain : public AdvectionDomain<Mapping>
{
public:
    using X_adv = X;
    using Y_adv = Y;
    using CoordXY_adv = Coord<X_adv, Y_adv>;
 
private:
    Mapping const& m_mapping;
 
public:
    AdvectionPhysicalDomain(Mapping const& mapping) : m_mapping(mapping) {};
    ~AdvectionPhysicalDomain() {};
 
    void advect_feet(
            host_t<FieldRTheta<CoordRTheta>> feet_coords_rp,
            host_t<DConstVectorFieldRTheta<X_adv, Y_adv>> advection_field,
            double dt) const
    {
        using namespace ddc;
 
        IdxRangeRTheta const idx_range_rp = get_idx_range<GridR, GridTheta>(feet_coords_rp);
        CoordXY coord_center(m_mapping(CoordRTheta(0, 0)));
 
        ddc::for_each(idx_range_rp, [&](IdxRTheta const irp) {
            CoordRTheta const coord_rp(feet_coords_rp(irp));
            CoordXY const coord_xy = m_mapping(coord_rp);
 
            CoordXY const feet_xy = coord_xy - dt * advection_field(irp);
 
            if (norm_inf(feet_xy - coord_center) < 1e-15) {
                feet_coords_rp(irp) = CoordRTheta(0, 0);
            } else {
                feet_coords_rp(irp) = m_mapping(feet_xy);
                ddc::select<Theta>(feet_coords_rp(irp)) = ddcHelper::restrict_to_idx_range(
                        ddc::select<Theta>(feet_coords_rp(irp)),
                        IdxRangeTheta(idx_range_rp));
            }
        });
    }
};
 
 
template <class Mapping>
class AdvectionPseudoCartesianDomain : public AdvectionDomain<Mapping>
{
public:
    using Matrix2x2 = std::array<std::array<double, 2>, 2>;
 
    using X_adv = DimX_pC;
    using Y_adv = DimY_pC;
    using CoordXY_adv = Coord<X_adv, Y_adv>;
 
 
private:
    Mapping const& m_mapping;
    double m_epsilon;
 
public:
    AdvectionPseudoCartesianDomain(Mapping const& mapping, double epsilon = 1e-12)
        : m_mapping(mapping)
        , m_epsilon(epsilon) {};
    ~AdvectionPseudoCartesianDomain() {};
 
    void advect_feet(
            host_t<FieldRTheta<CoordRTheta>> feet_coords_rp,
            host_t<DConstVectorFieldRTheta<X_adv, Y_adv>> const& advection_field,
            double const dt) const
    {
        static_assert(!std::is_same_v<Mapping, CircularToCartesian<X, Y, R, Theta>>);
        IdxRangeRTheta const idx_range_rp = get_idx_range(advection_field);
 
        CircularToCartesian<X_adv, Y_adv, R, Theta> const pseudo_Cartesian_mapping;
        CoordXY_adv const center_xy_pseudo_cart
                = CoordXY_adv(pseudo_Cartesian_mapping(CoordRTheta(0., 0.)));
 
        ddc::for_each(idx_range_rp, [&](IdxRTheta const irp) {
            CoordRTheta const coord_rp(feet_coords_rp(irp));
            CoordXY_adv const coord_xy_pseudo_cart = pseudo_Cartesian_mapping(coord_rp);
            CoordXY_adv const feet_xy_pseudo_cart
                    = coord_xy_pseudo_cart - dt * advection_field(irp);
 
            if (norm_inf(feet_xy_pseudo_cart - center_xy_pseudo_cart) < 1e-15) {
                feet_coords_rp(irp) = CoordRTheta(0, 0);
            } else {
                feet_coords_rp(irp) = pseudo_Cartesian_mapping(feet_xy_pseudo_cart);
                ddc::select<Theta>(feet_coords_rp(irp)) = ddcHelper::restrict_to_idx_range(
                        ddc::select<Theta>(feet_coords_rp(irp)),
                        IdxRangeTheta(idx_range_rp));
            }
        });
    }
};