advection_domain.hpp

// SPDX-License-Identifier: MIT
#pragma once
#include <cassert>
#include <typeinfo>
 
#include <sll/mapping/cartesian_to_circular.hpp>
#include <sll/mapping/circular_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 "l_norm_tools.hpp"
#include "vector_field.hpp"
#include "vector_field_mem.hpp"
 
template <class LogicalToPhysicalMapping>
class AdvectionDomain
{
    static_assert(is_analytical_mapping_v<LogicalToPhysicalMapping>);
 
private:
    using PhysicalToLogicalMapping = inverse_mapping_t<LogicalToPhysicalMapping>;
 
public:
    using X_adv = typename PhysicalToLogicalMapping::cartesian_tag_x;
    using Y_adv = typename PhysicalToLogicalMapping::cartesian_tag_y;
    using CoordXY_adv = Coord<X_adv, Y_adv>;
 
private:
    LogicalToPhysicalMapping m_to_cartesian_mapping;
    PhysicalToLogicalMapping m_to_curvilinear_mapping;
 
public:
    AdvectionDomain(LogicalToPhysicalMapping const& to_physical_mapping)
        : m_to_cartesian_mapping(to_physical_mapping)
        , m_to_curvilinear_mapping(to_physical_mapping.get_inverse_mapping())
    {
    }
 
    void advect_feet(
            host_t<FieldRTheta<CoordRTheta>> feet_coords_rp,
            host_t<DConstVectorFieldRTheta<X_adv, Y_adv>> advection_field,
            double dt) const
    {
        IdxRangeRTheta const idx_range_rp = get_idx_range<GridR, GridTheta>(feet_coords_rp);
 
        CoordXY_adv coord_center(m_to_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 = m_to_cartesian_mapping(coord_rp);
 
            CoordXY_adv 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_to_curvilinear_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));
            }
        });
    }
};