SolverConfig.hpp

Go to the documentation of this file.

// -------------------------------------------------------//
//
// SHAMROCK code for hydrodynamics
// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
//
// -------------------------------------------------------//
 
#pragma once
 
#include "shambase/exception.hpp"
#include "shambase/string.hpp"
#include "shambackends/vec.hpp"
#include "shamcomm/logs.hpp"
#include "shammodels/common/amr/AMRBlock.hpp"
#include "shammodels/ramses/config/enum_DragSolverMode.hpp"
#include "shammodels/ramses/config/enum_DustRiemannSolverMode.hpp"
#include "shammodels/ramses/config/enum_GravityMode.hpp"
#include "shammodels/ramses/config/enum_RiemannSolverMode.hpp"
#include "shammodels/ramses/config/enum_SlopeMode.hpp"
#include "shamrock/experimental_features.hpp"
#include "shamrock/patch/PatchDataLayerLayout.hpp"
#include "shamrock/scheduler/PatchScheduler.hpp"
#include <shamrock/io/json_std_optional.hpp>
#include <shamunits/Constants.hpp>
#include <shamunits/UnitSystem.hpp>
#include <stdexcept>
#include <variant>
 
namespace shammodels::basegodunov {
 
    struct DragConfig {
        DragSolverMode drag_solver_config = NoDrag;
        std::vector<f64> alphas;
        bool enable_frictional_heating
            = false; // 0 to turn off and 1 when all dissipation is deposited to the gas
    };
 
    struct DustConfig {
        DustRiemannSolverMode dust_riemann_config = NoDust;
        u32 ndust                                 = 0;
 
        inline bool is_dust_on() {
            if (dust_riemann_config != NoDust) {
 
                if (ndust == 0) {
                    throw shambase::make_except_with_loc<std::runtime_error>(
                        "Dust is on with ndust == 0");
                }
                return true;
            }
            return false;
        }
    };
    struct PassiveScalarGasConfig {
        u32 npscal_gas = 0;
 
        inline bool is_gas_passive_scalar_on() { return npscal_gas > 0; }
    };
 
    template<class Tvec>
    struct GravityConfig {
        using Tscal              = shambase::VecComponent<Tvec>;
        GravityMode gravity_mode = NoGravity;
        bool analytical_gravity  = false; // whether to use an external analytical gravity
        Tscal tol                = 1e-6;
        inline Tscal get_tolerance() { return tol; }
        inline bool is_gravity_on() { return gravity_mode != NoGravity; }
    };
 
    template<class Tvec>
    struct SolverStatusVar;
 
    template<class Tvec, class TgridVec>
    struct AMRMode {
 
        using Tscal = shambase::VecComponent<Tvec>;
 
        struct None {};
        struct DensityBased {
            Tscal crit_mass;
        };
 
        using mode = std::variant<None, DensityBased>;
 
        mode config = None{};
        void set_refine_none() { config = None{}; }
        void set_refine_density_based(Tscal crit_mass) { config = DensityBased{crit_mass}; }
 
        bool need_level_zero_compute() { return false; }
        bool need_amr_level_compute() { return false; }
    };
 
    struct BCConfig {
        enum class GhostType { Periodic = 0, Reflective = 1, Outflow = 2 };
 
        GhostType ghost_type_x = GhostType::Periodic;
        GhostType ghost_type_y = GhostType::Periodic;
        GhostType ghost_type_z = GhostType::Periodic;
 
        GhostType get_x() const { return ghost_type_x; }
        GhostType get_y() const { return ghost_type_y; }
        GhostType get_z() const { return ghost_type_z; }
 
        void set_x(GhostType ghost_type) { ghost_type_x = ghost_type; }
        void set_y(GhostType ghost_type) { ghost_type_y = ghost_type; }
        void set_z(GhostType ghost_type) { ghost_type_z = ghost_type; }
    };
 
    template<class Tvec, class TgridVec>
    struct SolverConfig;
 
}; // namespace shammodels::basegodunov
 
template<class Tvec>
struct shammodels::basegodunov::SolverStatusVar {
 
    using Tscal = shambase::VecComponent<Tvec>;
 
    Tscal time = 0; 
    Tscal dt   = 0; 
};
 
template<class Tvec, class TgridVec>
struct shammodels::basegodunov::SolverConfig {
 
    using Tscal = shambase::VecComponent<Tvec>;
 
    Tscal eos_gamma = 5. / 3.;
 
    Tscal grid_coord_to_pos_fact = 1;
 
    static constexpr u32 NsideBlockPow = 1;
    using AMRBlock                     = amr::AMRBlock<Tvec, TgridVec, NsideBlockPow>;
 
    inline void set_eos_gamma(Tscal gamma) { eos_gamma = gamma; }
 
    RiemannSolverMode riemann_config  = HLL;
    SlopeMode slope_config            = VanLeer_sym;
    bool face_half_time_interpolation = true;
 
    inline bool should_compute_rho_mean() { return is_gravity_on() && is_boundary_periodic(); }
 
    // Dust config
 
    DustConfig dust_config{};
    DragConfig drag_config{};
 
    inline bool is_dust_on() { return dust_config.is_dust_on(); }
    // get alpha values from user
    // alphas is the dust collision rate (the inverse of the stopping time)
    inline void set_alphas_static(f32 alpha_values) {
        StackEntry stack_lock{};
        drag_config.alphas.push_back(alpha_values);
    }
 
    // Dust config (END)
 
    BCConfig bc_config{};
 
    // Gas passive scalars config
 
    PassiveScalarGasConfig npscal_gas_config{};
 
    inline bool is_gas_passive_scalar_on() { return npscal_gas_config.is_gas_passive_scalar_on(); }
    // Gas passive scalars config (END)
 
    // Gravity config
    inline Tscal get_constant_G() {
        if (!unit_sys) {
            ON_RANK_0(logger::warn_ln("amr::Config", "the unit system is not set"));
            shamunits::Constants<Tscal> ctes{shamunits::UnitSystem<Tscal>{}};
            return ctes.G();
        } else {
            return shamunits::Constants<Tscal>{*unit_sys}.G();
        }
    }
    inline bool is_boundary_periodic() { return true; }
    GravityConfig<Tvec> gravity_config{};
    inline Tscal get_constant_4piG() {
        auto scal_G = get_constant_G();
        return 4 * M_PI * scal_G;
    }
    inline Tscal get_grav_tol() { return gravity_config.get_tolerance(); }
    inline bool is_gravity_on() { return gravity_config.is_gravity_on(); }
    inline bool is_coordinate_field_required() { return gravity_config.analytical_gravity; }
 
    // Gravity config (END)
 
    AMRMode<Tvec, TgridVec> amr_mode = {};
 
    // Units Config
 
    std::optional<shamunits::UnitSystem<Tscal>> unit_sys = {};
 
    inline void set_units(shamunits::UnitSystem<Tscal> new_sys) { unit_sys = new_sys; }
    // Units Config (END)
 
    PatchSchedulerConfig scheduler_conf = {};
 
    // Solver status variables
 
    using SolverStatusVar = SolverStatusVar<Tvec>;
    SolverStatusVar time_state;
    inline void set_time(Tscal t) { time_state.time = t; }
    inline void set_next_dt(Tscal dt) { time_state.dt = dt; }
    inline Tscal get_time() { return time_state.time; }
    inline Tscal get_dt() { return time_state.dt; }
 
    Tscal Csafe = 0.9;
    // Solver status variables (END)
 
    inline void check_config() {
        if (grid_coord_to_pos_fact <= 0) {
            shambase::throw_with_loc<std::runtime_error>(shambase::format(
                "grid_coord_to_pos_fact must be > 0, got {}", grid_coord_to_pos_fact));
        }
 
        if (is_dust_on()) {
            ON_RANK_0(logger::warn_ln("Ramses::SolverConfig", "Dust is experimental"));
        }
 
        if (is_gravity_on()) {
            ON_RANK_0(logger::warn_ln("Ramses::SolverConfig", "Self gravity is experimental"));
            u32 mode = gravity_config.gravity_mode;
 
            shamrock::experimental_feature_check(
                shambase::format(
                    "self gravity mode is not enabled but gravity mode is set to {} (> 0 whith 0 "
                    "== "
                    "NoGravity mode)",
                    mode));
        }
 
        if (!(eos_gamma > 1.0)) {
            shambase::throw_with_loc<std::invalid_argument>(
                shambase::format("Gamma must be > 1, currently Gamma = {}", eos_gamma));
        }
 
        if (is_gas_passive_scalar_on()) {
            ON_RANK_0(logger::warn_ln("Ramses::SolverConfig", "Passive scalars are experimental"));
            shamrock::experimental_feature_check(
                shambase::format(
                    "gas passive scalars mode is not enabled but gas passive scalars mode is set "
                    "to {}"
                    "> 0",
                    npscal_gas_config.npscal_gas));
        }
    }
 
    void set_layout(shamrock::patch::PatchDataLayerLayout &pdl);
};
 
namespace shammodels::basegodunov {
 
    template<class Tvec>
    inline void to_json(nlohmann::json &j, const SolverStatusVar<Tvec> &p) {
        j = nlohmann::json{{"time", p.time}, {"dt", p.dt}};
    }
 
    template<class Tvec>
    inline void from_json(const nlohmann::json &j, SolverStatusVar<Tvec> &p) {
        using Tscal = typename SolverStatusVar<Tvec>::Tscal;
        j.at("time").get_to<Tscal>(p.time);
        j.at("dt").get_to<Tscal>(p.dt);
    }
 
    template<class Tvec, class TgridVec>
    void to_json(nlohmann::json &j, const SolverConfig<Tvec, TgridVec> &p);
    template<class Tvec, class TgridVec>
    void from_json(const nlohmann::json &j, SolverConfig<Tvec, TgridVec> &p);
 
} // namespace shammodels::basegodunov