spline_foot_finder.hpp

// SPDX-License-Identifier: MIT
#pragma once
#include <functional>
 
#include <sll/mapping/cartesian_to_pseudo_cartesian.hpp>
#include <sll/mapping/circular_to_cartesian.hpp>
 
#include "ddc_alias_inline_functions.hpp"
#include "ddc_aliases.hpp"
#include "geometry_pseudo_cartesian.hpp"
#include "ifoot_finder.hpp"
#include "vector_mapper.hpp"
 
template <class TimeStepper, class AdvectionDomain, class Mapping>
class SplineFootFinder : public IFootFinder
{
private:
    using X_adv = typename AdvectionDomain::X_adv;
    using Y_adv = typename AdvectionDomain::Y_adv;
 
    using CircularToPseudoCartesian = CircularToCartesian<X_adv, Y_adv, R, Theta>;
    using AdvectionMapping = CartesianToPseudoCartesian<Mapping, CircularToPseudoCartesian>;
 
    TimeStepper const& m_time_stepper;
 
    AdvectionDomain const& m_advection_domain;
    AdvectionMapping m_mapping;
 
    SplineRThetaBuilder const& m_builder_advection_field;
    SplineRThetaEvaluatorConstBound const& m_evaluator_advection_field;
 
 
public:
    SplineFootFinder(
            TimeStepper const& time_stepper,
            AdvectionDomain const& advection_domain,
            Mapping const& mapping,
            SplineRThetaBuilder const& builder_advection_field,
            SplineRThetaEvaluatorConstBound const& evaluator_advection_field,
            double epsilon = 1e-12)
        : m_time_stepper(time_stepper)
        , m_advection_domain(advection_domain)
        , m_mapping(CircularToPseudoCartesian(), mapping, epsilon)
        , m_builder_advection_field(builder_advection_field)
        , m_evaluator_advection_field(evaluator_advection_field)
    {
    }
 
    void operator()(
            host_t<FieldRTheta<CoordRTheta>> feet,
            host_t<DConstVectorFieldRTheta<X, Y>> advection_field,
            double dt) const final
    {
        host_t<VectorSplineCoeffsMem2D<X_adv, Y_adv>> advection_field_in_adv_domain_coefs(
                get_spline_idx_range(m_builder_advection_field));
 
        // Compute the advection field in the advection domain.
        auto advection_field_in_adv_domain = create_geometry_mirror_view(
                Kokkos::DefaultHostExecutionSpace(),
                advection_field,
                m_mapping);
 
        // Get the coefficients of the advection field in the advection domain.
        m_builder_advection_field(
                ddcHelper::get<X_adv>(advection_field_in_adv_domain_coefs),
                ddcHelper::get<X_adv>(get_const_field(advection_field_in_adv_domain)));
        m_builder_advection_field(
                ddcHelper::get<Y_adv>(advection_field_in_adv_domain_coefs),
                ddcHelper::get<Y_adv>(get_const_field(advection_field_in_adv_domain)));
 
 
        // The function describing how the derivative of the evolve function is calculated.
        std::function<
                void(host_t<DVectorFieldRTheta<X_adv, Y_adv>>,
                     host_t<ConstFieldRTheta<CoordRTheta>>)>
                dy = [&](host_t<DVectorFieldRTheta<X_adv, Y_adv>> updated_advection_field,
                         host_t<ConstFieldRTheta<CoordRTheta>> feet) {
                    m_evaluator_advection_field(
                            get_field(ddcHelper::get<X_adv>(updated_advection_field)),
                            get_const_field(feet),
                            get_const_field(
                                    ddcHelper::get<X_adv>(advection_field_in_adv_domain_coefs)));
                    m_evaluator_advection_field(
                            get_field(ddcHelper::get<Y_adv>(updated_advection_field)),
                            get_const_field(feet),
                            get_const_field(
                                    ddcHelper::get<Y_adv>(advection_field_in_adv_domain_coefs)));
                };
 
        // The function describing how the value(s) are updated using the derivative.
        std::function<
                void(host_t<FieldRTheta<CoordRTheta>>,
                     host_t<DConstVectorFieldRTheta<X_adv, Y_adv>>,
                     double)>
                update_function = [&](host_t<FieldRTheta<CoordRTheta>> feet,
                                      host_t<DConstVectorFieldRTheta<X_adv, Y_adv>> advection_field,
                                      double dt) {
                    // Compute the characteristic feet at t^n:
                    m_advection_domain.advect_feet(feet, advection_field, dt);
 
                    // Treatment to conserve the C0 property of the advected function:
                    unify_value_at_center_pt(feet);
                    // Test if the values are the same at the center point
                    is_unified(feet);
                };
 
 
        // Solve the characteristic equation
        m_time_stepper.update(Kokkos::DefaultHostExecutionSpace(), feet, dt, dy, update_function);
 
        is_unified(feet);
    }
 
 
 
private:
    template <class T>
    void is_unified(Field<T, IdxRangeRTheta, Kokkos::HostSpace> const& values) const
    {
        IdxRangeR const r_idx_range = get_idx_range<GridR>(values);
        IdxRangeTheta const theta_idx_range = get_idx_range<GridTheta>(values);
        if (std::fabs(ddc::coordinate(r_idx_range.front())) < 1e-15) {
            ddc::for_each(theta_idx_range, [&](const IdxTheta ip) {
                if (norm_inf(
                            values(r_idx_range.front(), ip)
                            - values(r_idx_range.front(), theta_idx_range.front()))
                    > 1e-15) {
                    std::cout << "WARNING ! -> Discontinous at the center point." << std::endl;
                }
                assert(values(r_idx_range.front(), ip)
                       == values(r_idx_range.front(), theta_idx_range.front()));
            });
        }
    }
 
 
    template <class T>
    void unify_value_at_center_pt(FieldRTheta<T> values) const
    {
        IdxRangeR const r_idx_range = get_idx_range<GridR>(values);
        IdxRangeTheta const theta_idx_range = get_idx_range<GridTheta>(values);
        if (std::fabs(ddc::coordinate(r_idx_range.front())) < 1e-15) {
            ddc::for_each(theta_idx_range, [&](const IdxTheta ip) {
                values(r_idx_range.front(), ip)
                        = values(r_idx_range.front(), theta_idx_range.front());
            });
        }
    }
};