ABLForcing.H

Go to the documentation of this file.

#ifndef ABLFORCING_H
#define ABLFORCING_H
 
#include "src/equation_systems/icns/MomentumSource.H"
#include "src/core/SimTime.H"
#include "src/utilities/trig_ops.H"
#include "src/utilities/linear_interpolation.H"
#include "AMReX_REAL.H"
 
using namespace amrex::literals;
 
namespace kynema_sgf::pde::icns {

class ABLForcing : public MomentumSource::Register<ABLForcing>
{
public:
    static std::string identifier() { return "ABLForcing"; }
 
    explicit ABLForcing(const CFDSim& sim);
 
    ~ABLForcing() override;
 
    void operator()(
        int lev, FieldState fstate, amrex::MultiFab& src_term) const override;
 
    void set_target_velocities(amrex::Real ux, amrex::Real uy)
    {
        m_target_vel[0] = ux;
        m_target_vel[1] = uy;
    }
 
    void set_mean_velocities(amrex::Real ux, amrex::Real uy)
    {
        m_mean_vel[0] = ux;
        m_mean_vel[1] = uy;
 
        const auto& current_time = m_time.current_time();
        const auto& new_time = m_time.new_time();
        const auto& nph_time = 0.5_rt * (current_time + new_time);
        const auto& dt = m_time.delta_t();
        const auto& t_step = m_time.time_index();
 
        if (!m_vel_timetable.empty()) {
            // Forces should be applied at n+1/2. Because ABL forcing is a
            // delta, the difference between the target velocity (at n+1) and
            // the current velocity (at n) puts the force term at n+1/2
            const amrex::Real new_spd = ::kynema_sgf::interp::linear(
                m_time_table, m_speed_table, new_time);
            const amrex::Real new_dir = ::kynema_sgf::interp::linear_angle(
                m_time_table, m_direction_table, new_time,
                2.0_rt * std::numbers::pi_v<amrex::Real>);
 
            m_target_vel[0] = new_spd * std::cos(new_dir);
            m_target_vel[1] = new_spd * std::sin(new_dir);
        }
 
        m_abl_forcing[0] = (m_target_vel[0] - m_mean_vel[0]) / dt;
        m_abl_forcing[1] = (m_target_vel[1] - m_mean_vel[1]) / dt;
 
        if (m_write_force_timetable &&
            amrex::ParallelDescriptor::IOProcessor() &&
            (t_step % m_force_outfreq == 0) &&
            (current_time >= m_force_outstart)) {
            std::ofstream outfile;
            // Forces are recorded at n+1/2
            outfile.open(m_force_timetable, std::ios::out | std::ios::app);
            outfile << std::setprecision(17) << nph_time << "\t"
                    << m_abl_forcing[0] << "\t" << m_abl_forcing[1] << "\t"
                    << 0.0_rt << '\n';
        }
    }
 
    [[nodiscard]] amrex::RealArray abl_forcing() const { return m_abl_forcing; }
 
    [[nodiscard]] amrex::Real forcing_height() const
    {
        return m_forcing_height;
    }
 
private:
    const SimTime& m_time;
    const amrex::AmrCore& m_mesh;

    bool m_use_phase_ramp{false};

    int m_n_band{2};

    amrex::RealArray m_abl_forcing{0.0_rt};

    std::string m_vel_timetable;

    bool m_write_force_timetable{false};
    std::string m_force_timetable;
    int m_force_outfreq{1};
    amrex::Real m_force_outstart{0.0_rt};

    amrex::Vector<amrex::Real> m_time_table;

    amrex::Vector<amrex::Real> m_speed_table;

    amrex::Vector<amrex::Real> m_direction_table;

    amrex::Vector<amrex::Real> m_target_vel{0.0_rt, 0.0_rt, 0.0_rt};

    amrex::RealArray m_mean_vel{0.0_rt};

    amrex::Real m_forcing_height;

    amrex::Real m_forcing_mphase0;
    amrex::Real m_forcing_mphase1;

    amrex::Real m_water_level;

    const Field* m_vof;
};
 
} // namespace kynema_sgf::pde::icns
 
#endif /* ABLFORCING_H */