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 "config/AVConfig.hpp"
#include "config/BCConfig.hpp"
#include "shambackends/math.hpp"
#include "shambackends/typeAliasVec.hpp"
#include "shambackends/type_traits.hpp"
#include "shambackends/vec.hpp"
#include "shamcomm/worldInfo.hpp"
#include "shammath/sphkernels.hpp"
#include "shammodels/common/EOSConfig.hpp"
#include "shammodels/common/ExtForceConfig.hpp"
#include "shammodels/sph/config/MHDConfig.hpp"
#include "shamrock/experimental_features.hpp"
#include "shamrock/io/json_print_diff.hpp"
#include "shamrock/io/json_std_optional.hpp"
#include "shamrock/io/json_utils.hpp"
#include "shamrock/io/units_json.hpp"
#include "shamrock/patch/PatchDataLayerLayout.hpp"
#include "shamrock/scheduler/PatchScheduler.hpp"
#include "shamsys/NodeInstance.hpp"
#include "shamsys/legacy/log.hpp"
#include "shamtree/CompressedLeafBVH.hpp"
#include "shamtree/RadixTree.hpp"
#include <shamunits/Constants.hpp>
#include <shamunits/UnitSystem.hpp>
#include <stdexcept>
#include <variant>
#include <vector>
 
namespace shammodels::sph {

    template<class Tvec, template<class> class SPHKernel>
    struct SolverConfig;

    template<class Tvec>
    struct SolverStatusVar;

    template<class Tscal>
    struct CFLConfig {

        Tscal cfl_cour;

        Tscal cfl_force;

        Tscal cfl_multiplier_stiffness = 2;

        Tscal eta_sink = 0.05;
    };
 
    template<class Tvec>
    struct ParticleKillingConfig {
        using Tscal = shambase::VecComponent<Tvec>;
        struct Sphere {
            Tvec center;
            Tscal radius;
        };
 
        using kill_t = std::variant<Sphere>;
 
        std::vector<kill_t> kill_list;
 
        inline void add_kill_sphere(const Tvec &center, Tscal radius) {
            kill_list.push_back(Sphere{center, radius});
        }
    };
 
    template<class Tscal>
    struct DustConfig {
 
        struct None {};
 
        struct MonofluidTVI {
            u32 ndust;
            bool pure_diffusion_mode = false;
        };
 
        struct MonofluidComplete {
            u32 ndust;
        };

        using Variant = std::variant<None, MonofluidTVI, MonofluidComplete>;
 
        Variant current_mode = None{};
 
        inline void set_none() { current_mode = None{}; }
        inline void set_monofluid_tvi(u32 nvar, bool pure_diffusion_mode = false) {
            current_mode = MonofluidTVI{nvar, pure_diffusion_mode};
        }
        inline void set_monofluid_complete(u32 nvar) { current_mode = MonofluidComplete{nvar}; }
 
        inline bool is_none() { return std::holds_alternative<None>(current_mode); }
        inline bool is_monofluid_tvi() { return bool(std::get_if<MonofluidTVI>(&current_mode)); }
        inline bool is_monofluid_complete() {
            return bool(std::get_if<MonofluidComplete>(&current_mode));
        }
 
        inline void mode_to_json(nlohmann::json &j) const {
            if (const None *cfg = std::get_if<None>(&current_mode)) {
                j = {{"type", "none"}};
            } else if (const MonofluidTVI *cfg = std::get_if<MonofluidTVI>(&current_mode)) {
                j
                    = {{"type", "monofluid_tvi"},
                       {"ndust", cfg->ndust},
                       {"pure_diffusion_mode", cfg->pure_diffusion_mode}};
            } else if (
                const MonofluidComplete *cfg = std::get_if<MonofluidComplete>(&current_mode)) {
                j = {{"type", "monofluid_complete"}, {"ndust", cfg->ndust}};
            } else {
                shambase::throw_unimplemented();
            }
        }
 
        inline void mode_from_json(const nlohmann::json &j) {
            const std::string type = j.at("type").get<std::string>();
            if (type == "none") {
                set_none();
            } else if (type == "monofluid_tvi") {
                set_monofluid_tvi(
                    j.at("ndust").get<u32>(), j.at("pure_diffusion_mode").get<bool>());
            } else if (type == "monofluid_complete") {
                set_monofluid_complete(j.at("ndust").get<u32>());
            } else {
                shambase::throw_unimplemented();
            }
        }
 
        inline bool has_s_j_field() {
            return is_monofluid_tvi(); // S_j = sqrt(\rho \epsilon_j)
        }
 
        inline bool has_epsilon_field() {
            return bool(std::get_if<MonofluidComplete>(&current_mode));
        }
 
        inline bool has_deltav_field() {
            return bool(std::get_if<MonofluidComplete>(&current_mode));
        }
 
        inline u32 get_dust_nvar() {
            if (None *cfg = std::get_if<None>(&current_mode)) {
                shambase::throw_with_loc<std::invalid_argument>(
                    "Querying a dust nvar with no dust as config is ... discutable ...");
                return 0;
            } else if (MonofluidTVI *cfg = std::get_if<MonofluidTVI>(&current_mode)) {
                return cfg->ndust;
            } else if (MonofluidComplete *cfg = std::get_if<MonofluidComplete>(&current_mode)) {
                return cfg->ndust;
            } else {
                shambase::throw_unimplemented("How did you get here ???");
            }
            return 0;
        }
 
        struct ConstantStoppingTimes {
            std::vector<Tscal> stopping_times;
        };
 
        struct EpsteinDrag {
            static constexpr bool supersonic_correction = false;
            Tscal gamma;
            std::vector<Tscal> grains_sizes;
            std::vector<Tscal> grains_densities;
        };
 
        std::variant<None, ConstantStoppingTimes, EpsteinDrag> dust_drag_mode = None{};
 
        inline void drag_mode_to_json(nlohmann::json &j) const {
            if (std::holds_alternative<None>(dust_drag_mode)) {
                j = {{"type", "none"}};
            } else if (
                const ConstantStoppingTimes *cfg
                = std::get_if<ConstantStoppingTimes>(&dust_drag_mode)) {
                j = {{"type", "constant_stopping_times"}, {"stopping_times", cfg->stopping_times}};
            } else if (const EpsteinDrag *cfg = std::get_if<EpsteinDrag>(&dust_drag_mode)) {
                j
                    = {{"type", "epstein_drag"},
                       {"gamma", cfg->gamma},
                       {"grains_sizes", cfg->grains_sizes},
                       {"grains_densities", cfg->grains_densities}};
            } else {
                shambase::throw_unimplemented();
            }
        }
 
        inline void drag_mode_from_json(const nlohmann::json &j) {
            if (j.at("type").get<std::string>() == "none") {
                dust_drag_mode = None{};
            } else if (j.at("type").get<std::string>() == "constant_stopping_times") {
                dust_drag_mode
                    = ConstantStoppingTimes{j.at("stopping_times").get<std::vector<Tscal>>()};
            } else if (j.at("type").get<std::string>() == "epstein_drag") {
                dust_drag_mode = EpsteinDrag{
                    j.at("gamma").get<Tscal>(),
                    j.at("grains_sizes").get<std::vector<Tscal>>(),
                    j.at("grains_densities").get<std::vector<Tscal>>()};
            } else {
                shambase::throw_unimplemented();
            }
        }
 
        inline void set_drag_constant(ConstantStoppingTimes in) { dust_drag_mode = std::move(in); }
 
        inline void set_drag_epstein(EpsteinDrag in) { dust_drag_mode = std::move(in); }
 
        inline void check_config() {
            bool is_not_none = !is_none();
            if (is_not_none) {
 
                if (!shamrock::are_experimental_features_allowed()) {
                    shambase::throw_with_loc<std::runtime_error>(
                        "Dust config != None is experimental");
                } else {
                    ON_RANK_0(
                        logger::warn_ln(
                            "SPH::config",
                            "Dust config != None is work in progress, use it at your own risk"));
                }
 
                if (std::holds_alternative<None>(dust_drag_mode)) {
                    throw shambase::make_except_with_loc<std::runtime_error>(
                        "you must select a drag mode for the dust if the dust is on !");
                } else if (
                    ConstantStoppingTimes *cfg
                    = std::get_if<ConstantStoppingTimes>(&dust_drag_mode)) {
                    if (get_dust_nvar() != cfg->stopping_times.size()) {
                        throw shambase::make_except_with_loc<std::invalid_argument>(
                            "stopping_times size does not match the number of dust bins");
                    }
                } else if (EpsteinDrag *cfg = std::get_if<EpsteinDrag>(&dust_drag_mode)) {
                    if (get_dust_nvar() != cfg->grains_densities.size()) {
                        throw shambase::make_except_with_loc<std::invalid_argument>(
                            "grains_densities size does not match the number of dust bins");
                    }
 
                    if (get_dust_nvar() != cfg->grains_sizes.size()) {
                        throw shambase::make_except_with_loc<std::invalid_argument>(
                            "grains_sizes size does not match the number of dust bins");
                    }
                }
            }
        }
    };
 
    struct SmoothingLengthConfig {
        struct DensityBased {};
        struct DensityBasedNeighLim {
            u32 max_neigh_count = 500;
        };
 
        using mode = std::variant<DensityBased, DensityBasedNeighLim>;
 
        mode config = DensityBased{};
 
        void set_density_based() { config = DensityBased{}; }
        void set_density_based_neigh_lim(u32 max_neigh_count) {
            config = DensityBasedNeighLim{max_neigh_count};
        }
 
        bool is_density_based_neigh_lim() const {
            return std::holds_alternative<DensityBasedNeighLim>(config);
        }
    };
 
    struct SelfGravConfig {
 
        struct SFMM {
            u32 order;
            f64 opening_angle;
            bool leaf_lowering;
            u32 reduction_level;
        };
 
        struct FMM {
            u32 order;
            f64 opening_angle;
            u32 reduction_level;
        };
 
        struct MM {
            u32 order;
            f64 opening_angle;
            u32 reduction_level;
        };
 
        struct Direct {
            bool reference_mode = false;
        };
 
        struct None {};
 
        using mode = std::variant<SFMM, FMM, MM, Direct, None>;
 
        mode config = None{};
 
        void set_none() { config = None{}; }
        void set_direct(bool reference_mode = false) { config = Direct{reference_mode}; }
        void set_mm(u32 mm_order, f64 opening_angle, u32 reduction_level) {
            config = MM{
                .order           = mm_order,
                .opening_angle   = opening_angle,
                .reduction_level = reduction_level};
        }
        void set_fmm(u32 order, f64 opening_angle, u32 reduction_level) {
            config = FMM{
                .order = order, .opening_angle = opening_angle, .reduction_level = reduction_level};
        }
        void set_sfmm(u32 order, f64 opening_angle, bool leaf_lowering, u32 reduction_level) {
            config = SFMM{
                .order           = order,
                .opening_angle   = opening_angle,
                .leaf_lowering   = leaf_lowering,
                .reduction_level = reduction_level};
        }
 
        bool is_none() const { return std::holds_alternative<None>(config); }
        bool is_direct() const { return std::holds_alternative<Direct>(config); }
        bool is_mm() const { return std::holds_alternative<MM>(config); }
        bool is_fmm() const { return std::holds_alternative<FMM>(config); }
        bool is_sfmm() const { return std::holds_alternative<SFMM>(config); }
 
        bool is_sg_on() const { return !is_none(); }
        bool is_sg_off() const { return is_none(); }
 
        struct SofteningPlummer {
            f64 epsilon;
        };
 
        using mode_soft          = std::variant<SofteningPlummer>;
        mode_soft softening_mode = SofteningPlummer{1e-9};
 
        void set_softening_plummer(f64 epsilon) { softening_mode = SofteningPlummer{epsilon}; }
        void set_softening_none() { set_softening_plummer(0.); }
 
        bool is_softening_plummer() const {
            return std::holds_alternative<SofteningPlummer>(softening_mode);
        }
    };
 
} // namespace shammodels::sph
 
template<class Tvec>
struct shammodels::sph::SolverStatusVar {

    using Tscal = shambase::VecComponent<Tvec>;
 
    Tscal time   = 0; 
    Tscal dt_sph = 0; 
 
    Tscal cfl_multiplier = 1e-2; 
};
 
template<class Tvec, template<class> class SPHKernel>
struct shammodels::sph::SolverConfig {

    using Tscal = shambase::VecComponent<Tvec>;
    static constexpr u32 dim = shambase::VectorProperties<Tvec>::dimension;
    using Kernel = SPHKernel<Tscal>;
    using u_morton = u32;
 
    using RTree = shamtree::CompressedLeafBVH<u_morton, Tvec, 3>;

    static constexpr Tscal Rkern = Kernel::Rkern;
 
    Tscal gpart_mass;            
    CFLConfig<Tscal> cfl_config; 
 
    bool track_particles_id = false;
 
    inline void set_particle_tracking(bool state) { track_particles_id = state; }
 
    PatchSchedulerConfig scheduler_conf = {};

    // Units Config

    std::optional<shamunits::UnitSystem<Tscal>> unit_sys = {};

    inline void set_units(shamunits::UnitSystem<Tscal> new_sys) { unit_sys = new_sys; }

    inline Tscal get_constant_G() {
        if (!unit_sys) {
            ON_RANK_0(logger::warn_ln("sph::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 Tscal get_constant_c() {
        if (!unit_sys) {
            ON_RANK_0(logger::warn_ln("sph::Config", "the unit system is not set"));
            shamunits::Constants<Tscal> ctes{shamunits::UnitSystem<Tscal>{}};
            return ctes.c();
        } else {
            return shamunits::Constants<Tscal>{*unit_sys}.c();
        }
    }

    inline Tscal get_constant_mu_0() {
        if (!unit_sys) {
            ON_RANK_0(logger::warn_ln("sph::Config", "the unit system is not set"));
            shamunits::Constants<Tscal> ctes{shamunits::UnitSystem<Tscal>{}};
            return ctes.mu_0();
        } else {
            return shamunits::Constants<Tscal>{*unit_sys}.mu_0();
        }
    }

    // Units Config (END)

    // Particle killing config
 
    ParticleKillingConfig<Tvec> particle_killing;

    // Particle killing config (END)

    // 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_sph = dt; }

    inline Tscal get_time() { return time_state.time; }

    inline Tscal get_dt_sph() { return time_state.dt_sph; }

    inline void set_cfl_multipler(Tscal lambda) { time_state.cfl_multiplier = lambda; }

    inline Tscal get_cfl_multipler() { return time_state.cfl_multiplier; }

    inline void set_cfl_mult_stiffness(Tscal cstiff) {
        cfl_config.cfl_multiplier_stiffness = cstiff;
    }

    inline Tscal get_cfl_mult_stiffness() { return cfl_config.cfl_multiplier_stiffness; }
 
    bool show_cfl_detail = false;

    // Solver status variables (END)

    // MHD Config
 
    using MHDConfig      = MHDConfig<Tvec>;
    MHDConfig mhd_config = {};

    inline void set_noMHD() {
        using Tmp = typename MHDConfig::None;
        mhd_config.set(Tmp{});
    }

    inline void set_IdealMHD(typename MHDConfig::IdealMHD_constrained_hyper_para v) {
        mhd_config.set(v);
    }
 
    inline void set_NonIdealMHD(typename MHDConfig::NonIdealMHD v) { mhd_config.set(v); }

    // MHD Config (END)

    // Dust config
 
    using DustConfig       = DustConfig<Tscal>;
    DustConfig dust_config = {};

    // Dust config (END)

    // Self gravity config
 
    SelfGravConfig self_grav_config = SelfGravConfig{};

    // Self gravity config (END)

    // Tree config
 
    u32 tree_reduction_level  = 3;    
    bool use_two_stage_search = true; 

    inline void set_tree_reduction_level(u32 level) { tree_reduction_level = level; }
    inline void set_two_stage_search(bool enable) { use_two_stage_search = enable; }
 
    bool show_neigh_stats = false;
    inline void set_show_neigh_stats(bool enable) { show_neigh_stats = enable; }
    // Tree config (END)

    // Solver behavior config
 
    bool combined_dtdiv_divcurlv_compute = false; 
    Tscal htol_up_coarse_cycle = 1.1;
    Tscal htol_up_fine_cycle  = 1.1;
    Tscal epsilon_h           = 1e-6; 
    u32 h_iter_per_subcycles  = 50;   
    u32 h_max_subcycles_count = 100;  
 
    SmoothingLengthConfig smoothing_length_config;
 
    inline void set_smoothing_length_density_based() {
        smoothing_length_config.set_density_based();
    }
    inline void set_smoothing_length_density_based_neigh_lim(u32 max_neigh_count) {
        smoothing_length_config.set_density_based_neigh_lim(max_neigh_count);
    }
 
    bool enable_particle_reordering = false;
    inline void set_enable_particle_reordering(bool enable) { enable_particle_reordering = enable; }
    u64 particle_reordering_step_freq = 1000;
    inline void set_particle_reordering_step_freq(u64 freq) {
        if (freq == 0) {
            shambase::throw_with_loc<std::invalid_argument>(
                "particle_reordering_step_freq cannot be zero");
        }
        particle_reordering_step_freq = freq;
    }
 
    bool save_dt_to_fields = false;
    inline void set_save_dt_to_fields(bool enable) { save_dt_to_fields = enable; }
    inline bool should_save_dt_to_fields() const { return save_dt_to_fields; }
 
    bool show_ghost_zone_graph = false;
    inline void set_show_ghost_zone_graph(bool enable) { show_ghost_zone_graph = enable; }

    // Solver behavior config (END)

    // EOS Config

    using EOSConfig = shammodels::EOSConfig<Tvec>;

    EOSConfig eos_config;

    inline bool is_eos_locally_isothermal() {
        using T = typename EOSConfig::LocallyIsothermal;
        return bool(std::get_if<T>(&eos_config.config));
    }

    inline bool is_eos_adiabatic() {
        using T = typename EOSConfig::Adiabatic;
        return bool(std::get_if<T>(&eos_config.config));
    }

    inline bool is_eos_polytropic() {
        using T = typename EOSConfig::Polytropic;
        return bool(std::get_if<T>(&eos_config.config));
    }

    inline bool is_eos_isothermal() {
        using T = typename EOSConfig::Isothermal;
        return bool(std::get_if<T>(&eos_config.config));
    }

    inline bool is_eos_fermi() {
        using T = typename EOSConfig::Fermi;
        return bool(std::get_if<T>(&eos_config.config));
    }

    inline void set_eos_isothermal(Tscal cs) { eos_config.set_isothermal(cs); }

    inline void set_eos_adiabatic(Tscal gamma) { eos_config.set_adiabatic(gamma); }

    inline void set_eos_polytropic(Tscal K, Tscal gamma) { eos_config.set_polytropic(K, gamma); }

    inline void set_eos_locally_isothermal() { eos_config.set_locally_isothermal(); }

    inline void set_eos_locally_isothermalLP07(Tscal cs0, Tscal q, Tscal r0) {
        eos_config.set_locally_isothermalLP07(cs0, q, r0);
    }

    inline void set_eos_locally_isothermalFA2014(Tscal h_over_r) {
        eos_config.set_locally_isothermalFA2014(h_over_r);
    }

    inline void set_eos_locally_isothermalFA2014_extended(
        Tscal cs0, Tscal q, Tscal r0, u32 n_sinks) {
        eos_config.set_locally_isothermalFA2014_extended(cs0, q, r0, n_sinks);
    }

    inline void set_eos_fermi(Tscal mu_e) { eos_config.set_fermi(mu_e); }

    // EOS Config (END)

    // Artificial viscosity Config

    using AVConfig = AVConfig<Tvec>;

    AVConfig artif_viscosity;

    inline void set_artif_viscosity_None() {
        using Tmp = typename AVConfig::None;
        artif_viscosity.set(Tmp{});
    }

    inline void set_artif_viscosity_Constant(typename AVConfig::Constant v) {
        artif_viscosity.set(v);
    }

    inline void set_artif_viscosity_VaryingMM97(typename AVConfig::VaryingMM97 v) {
        artif_viscosity.set(v);
    }

    inline void set_artif_viscosity_VaryingCD10(typename AVConfig::VaryingCD10 v) {
        artif_viscosity.set(v);
    }

    inline void set_artif_viscosity_ConstantDisc(typename AVConfig::ConstantDisc v) {
        artif_viscosity.set(v);
    }

    // Artificial viscosity Config (END)

    // Boundary Config

    using BCConfig = BCConfig<Tvec>;

    BCConfig boundary_config;

    inline void set_boundary_free() { boundary_config.set_free(); }

    inline void set_boundary_periodic() { boundary_config.set_periodic(); }

    inline void set_boundary_shearing_periodic(i32_3 shear_base, i32_3 shear_dir, Tscal speed) {
        boundary_config.set_shearing_periodic(shear_base, shear_dir, speed);
    }

    // Boundary Config (END)

    // Ext force Config

    using ExtForceConfig = shammodels::ExtForceConfig<Tvec>;

    ExtForceConfig ext_force_config{};

    inline void add_ext_force_point_mass(Tscal central_mass, Tscal Racc) {
        ext_force_config.add_point_mass(central_mass, Racc);
    }

    inline void add_ext_force_paczynski_wiita(Tscal central_mass, Tvec central_pos, Tscal Racc) {
        ext_force_config.add_paczynski_wiita(central_mass, central_pos, Racc);
    }

    inline void add_ext_force_lense_thirring(
        Tscal central_mass, Tscal Racc, Tscal a_spin, Tvec dir_spin) {
        ext_force_config.add_lense_thirring(central_mass, Racc, a_spin, dir_spin);
    }

    inline void add_ext_force_shearing_box(Tscal Omega_0, Tscal eta, Tscal q) {
        ext_force_config.add_shearing_box(Omega_0, eta, q);
    }

    // Ext force Config (END)

    // Debug dump config

    bool do_debug_dump = false;

    std::string debug_dump_filename = "";

    inline void set_debug_dump(bool _do_debug_dump, std::string _debug_dump_filename) {
        this->do_debug_dump       = _do_debug_dump;
        this->debug_dump_filename = _debug_dump_filename;
    }

    inline constexpr bool do_MHD_debug() { return false; }

    // Debug dump config (END)

    inline bool ghost_has_soundspeed() { return is_eos_locally_isothermal(); }

    inline bool has_field_uint() {
        // no barotropic for now
        return true;
    }

    inline bool has_field_alphaAV() { return artif_viscosity.has_alphaAV_field(); }

    inline bool has_field_divv() { return artif_viscosity.has_alphaAV_field(); }

    inline bool has_field_dtdivv() { return artif_viscosity.has_dtdivv_field(); }

    inline bool has_field_curlv() { return artif_viscosity.has_curlv_field() && (dim == 3); }

    inline bool has_axyz_in_ghost() { return has_field_dtdivv(); }

    inline bool has_field_soundspeed() {
        return artif_viscosity.has_field_soundspeed() || is_eos_locally_isothermal();
    }

    inline bool has_field_B_on_rho() { return mhd_config.has_B_field() && (dim == 3); }

    inline bool has_field_psi_on_ch() { return mhd_config.has_psi_field(); }

    inline bool has_field_divB() { return mhd_config.has_divB_field(); }

    inline bool has_field_curlB() { return mhd_config.has_curlB_field() && (dim == 3); }

    inline bool has_field_dtdivB() { return mhd_config.has_dtdivB_field(); }

    bool compute_luminosity = false;
    inline void use_luminosity(bool enable) { compute_luminosity = enable; }

    inline void print_status() {
        if (shamcomm::world_rank() != 0) {
            return;
        }
        logger::raw_ln("----- SPH Solver configuration -----");
        logger::raw_ln(nlohmann::json{*this}.dump(4));
        logger::raw_ln("------------------------------------");
    }
 
    inline void check_config() {
        dust_config.check_config();
 
        if (track_particles_id && false /*particle injection when added*/) {
            shamrock::experimental_feature_check(
                "particle injection is not yet compatible with particle id tracking");
        }
 
        if (track_particles_id) {
            shamrock::experimental_feature_check("Particle tracking is experimental");
        }
 
        if (!self_grav_config.is_none()) {
            shamrock::experimental_feature_check(
                "Self gravity is experimental, please enable experimental features to use it");
        }
    }
 
    void set_layout(shamrock::patch::PatchDataLayerLayout &pdl);
    void set_ghost_layout(shamrock::patch::PatchDataLayerLayout &ghost_layout);
};
 
namespace shammodels::sph {

    template<class Tscal>
    inline void to_json(nlohmann::json &j, const CFLConfig<Tscal> &p) {
        j = nlohmann::json{
            {"cfl_cour", p.cfl_cour},
            {"cfl_force", p.cfl_force},
            {"cfl_multiplier_stiffness", p.cfl_multiplier_stiffness},
            {"eta_sink", p.eta_sink}};
    }

    template<class Tscal>
    inline void from_json(const nlohmann::json &j, CFLConfig<Tscal> &p) {
        j.at("cfl_cour").get_to<Tscal>(p.cfl_cour);
        j.at("cfl_force").get_to<Tscal>(p.cfl_force);
        j.at("cfl_multiplier_stiffness").get_to<Tscal>(p.cfl_multiplier_stiffness);
 
        if (j.contains("eta_sink")) {
            j.at("eta_sink").get_to<Tscal>(p.eta_sink);
        } else {
            // Already set to default value
            ON_RANK_0(shamlog_warn_ln(
                "SPHConfig", "eta_sink not found when deserializing, defaulting to", p.eta_sink));
        }
    }

    template<class Tvec>
    inline void to_json(nlohmann::json &j, const SolverStatusVar<Tvec> &p) {
        j = nlohmann::json{
            {"time", p.time}, {"dt_sph", p.dt_sph}, {"cfl_multiplier", p.cfl_multiplier}};
    }

    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_sph").get_to<Tscal>(p.dt_sph);
        j.at("cfl_multiplier").get_to<Tscal>(p.cfl_multiplier);
    }
 
    // JSON serialization for ParticleKillingConfig
    template<class Tvec>
    inline void to_json(nlohmann::json &j, const ParticleKillingConfig<Tvec> &p) {
        j = nlohmann::json::array();
        for (const auto &kill : p.kill_list) {
            if (std::holds_alternative<typename ParticleKillingConfig<Tvec>::Sphere>(kill)) {
                const auto &sphere = std::get<typename ParticleKillingConfig<Tvec>::Sphere>(kill);
                j.push_back(
                    {{"type", "sphere"}, {"center", sphere.center}, {"radius", sphere.radius}});
            }
            // If more types are added to kill_t, handle them here
        }
    }
 
    template<class Tvec>
    inline void from_json(const nlohmann::json &j, ParticleKillingConfig<Tvec> &p) {
        p.kill_list.clear();
        for (const auto &item : j) {
            std::string type = item.at("type").get<std::string>();
            if (type == "sphere") {
                typename ParticleKillingConfig<Tvec>::Sphere sphere;
                item.at("center").get_to(sphere.center);
                item.at("radius").get_to(sphere.radius);
                p.kill_list.push_back(sphere);
            }
            // If more types are added to kill_t, handle them here
        }
    }
 
    // JSON serialization for SmoothingLengthConfig
    inline void to_json(nlohmann::json &j, const SmoothingLengthConfig &p) {
        if (const SmoothingLengthConfig::DensityBased *conf
            = std::get_if<SmoothingLengthConfig::DensityBased>(&p.config)) {
            j = {
                {"type", "density_based"},
            };
 
        } else if (
            const SmoothingLengthConfig::DensityBasedNeighLim *conf
            = std::get_if<SmoothingLengthConfig::DensityBasedNeighLim>(&p.config)) {
 
            j = {
                {"type", "density_based_neigh_lim"},
                {"max_neigh_count", conf->max_neigh_count},
            };
        } else {
            shambase::throw_unimplemented();
        }
    }
 
    inline void from_json(const nlohmann::json &j, SmoothingLengthConfig &p) {
        if (j.at("type").get<std::string>() == "density_based") {
            p.config = SmoothingLengthConfig::DensityBased{};
        } else if (j.at("type").get<std::string>() == "density_based_neigh_lim") {
            p.config
                = SmoothingLengthConfig::DensityBasedNeighLim{j.at("max_neigh_count").get<u32>()};
        } else {
            shambase::throw_unimplemented();
        }
    }

    inline void to_json(nlohmann::json &j, const SelfGravConfig &p) {
        if (const SelfGravConfig::SFMM *conf = std::get_if<SelfGravConfig::SFMM>(&p.config)) {
            j = {
                {"type", "sfmm"},
                {"order", conf->order},
                {"opening_angle", conf->opening_angle},
                {"reduction_level", conf->reduction_level},
                {"leaf_lowering", conf->leaf_lowering},
            };
        } else if (const SelfGravConfig::FMM *conf = std::get_if<SelfGravConfig::FMM>(&p.config)) {
            j = {
                {"type", "fmm"},
                {"order", conf->order},
                {"opening_angle", conf->opening_angle},
                {"reduction_level", conf->reduction_level},
            };
        } else if (const SelfGravConfig::MM *conf = std::get_if<SelfGravConfig::MM>(&p.config)) {
            j = {
                {"type", "mm"},
                {"order", conf->order},
                {"opening_angle", conf->opening_angle},
                {"reduction_level", conf->reduction_level},
            };
        } else if (
            const SelfGravConfig::Direct *conf = std::get_if<SelfGravConfig::Direct>(&p.config)) {
            j = {
                {"type", "direct"},
                {"reference_mode", conf->reference_mode},
            };
        } else if (
            const SelfGravConfig::None *conf = std::get_if<SelfGravConfig::None>(&p.config)) {
            j = {
                {"type", "none"},
            };
        }
 
        if (const SelfGravConfig::SofteningPlummer *conf
            = std::get_if<SelfGravConfig::SofteningPlummer>(&p.softening_mode)) {
            j["softening_mode"]   = "plummer";
            j["softening_length"] = conf->epsilon;
        } else {
            shambase::throw_unimplemented();
        }
    }

    inline void from_json(const nlohmann::json &j, SelfGravConfig &p) {
        if (j.at("type").get<std::string>() == "sfmm") {
            p.config = SelfGravConfig::SFMM{
                .order           = j.at("order").get<u32>(),
                .opening_angle   = j.at("opening_angle").get<f64>(),
                .leaf_lowering   = j.at("leaf_lowering").get<bool>(),
                .reduction_level = j.at("reduction_level").get<u32>()};
        } else if (j.at("type").get<std::string>() == "fmm") {
            p.config = SelfGravConfig::FMM{
                .order           = j.at("order").get<u32>(),
                .opening_angle   = j.at("opening_angle").get<f64>(),
                .reduction_level = j.at("reduction_level").get<u32>()};
        } else if (j.at("type").get<std::string>() == "mm") {
            p.config = SelfGravConfig::MM{
                .order           = j.at("order").get<u32>(),
                .opening_angle   = j.at("opening_angle").get<f64>(),
                .reduction_level = j.at("reduction_level").get<u32>()};
        } else if (j.at("type").get<std::string>() == "direct") {
            p.config = SelfGravConfig::Direct{j.at("reference_mode").get<bool>()};
        } else if (j.at("type").get<std::string>() == "none") {
            p.config = SelfGravConfig::None{};
        } else {
            throw shambase::make_except_with_loc<std::runtime_error>(
                "Invalid self gravity type: " + j.at("type").get<std::string>());
        }
 
        if (j.contains("softening_mode")) {
            std::string softening_mode = j.at("softening_mode").get<std::string>();
            if (softening_mode == "plummer") {
                p.softening_mode
                    = SelfGravConfig::SofteningPlummer{j.at("softening_length").get<f64>()};
            } else {
                throw shambase::make_except_with_loc<std::runtime_error>(
                    "Invalid softening mode: " + softening_mode);
            }
        }
    }
 
    // JSON serialization for DustConfig
    template<class Tvec>
    inline void to_json(nlohmann::json &j, const DustConfig<Tvec> &p) {
        j = {};
 
        p.mode_to_json(j["mode"]);
        p.drag_mode_to_json(j["drag_mode"]);
    }
 
    template<class Tvec>
    inline void from_json(const nlohmann::json &j, DustConfig<Tvec> &p) {
        p.mode_from_json(j.at("mode"));
        p.drag_mode_from_json(j.at("drag_mode"));
    }

    template<class Tvec, template<class> class SPHKernel>
    inline void to_json(nlohmann::json &j, const SolverConfig<Tvec, SPHKernel> &p) {
        using T       = SolverConfig<Tvec, SPHKernel>;
        using Tkernel = typename T::Kernel;
 
        std::string kernel_id = shambase::get_type_name<Tkernel>();
        std::string type_id   = shambase::get_type_name<Tvec>();
 
        j = nlohmann::json{
            // used for type checking
            {"kernel_id", kernel_id},
            {"type_id", type_id},
            // scheduler config
            {"scheduler_config", p.scheduler_conf},
            // actual data stored in the json
            {"gpart_mass", p.gpart_mass},
            {"cfl_config", p.cfl_config},
            {"unit_sys", p.unit_sys},
            {"time_state", p.time_state},
            {"show_cfl_detail", p.show_cfl_detail},
            // mhd config
            {"mhd_config", p.mhd_config},
            // dust config
            {"dust_config", p.dust_config},
            // self gravity config
            {"self_grav_config", p.self_grav_config},
            // tree config
            {"tree_reduction_level", p.tree_reduction_level},
            {"use_two_stage_search", p.use_two_stage_search},
            {"show_neigh_stats", p.show_neigh_stats},
            // solver behavior config
            {"combined_dtdiv_divcurlv_compute", p.combined_dtdiv_divcurlv_compute},
            {"htol_up_coarse_cycle", p.htol_up_coarse_cycle},
            {"htol_up_fine_cycle", p.htol_up_fine_cycle},
            {"epsilon_h", p.epsilon_h},
            {"smoothing_length_config", p.smoothing_length_config},
            {"h_iter_per_subcycles", p.h_iter_per_subcycles},
            {"h_max_subcycles_count", p.h_max_subcycles_count},
 
            {"enable_particle_reordering", p.enable_particle_reordering},
            {"particle_reordering_step_freq", p.particle_reordering_step_freq},
 
            {"save_dt_to_fields", p.save_dt_to_fields},
            {"show_ghost_zone_graph", p.show_ghost_zone_graph},
 
            {"eos_config", p.eos_config},
 
            {"artif_viscosity", p.artif_viscosity},
            {"boundary_config", p.boundary_config},
            {"ext_force_config", p.ext_force_config},
 
            {"do_debug_dump", p.do_debug_dump},
            {"debug_dump_filename", p.debug_dump_filename},
            // particle killing config
            {"particle_killing", p.particle_killing},
        };
    }

    template<class Tvec, template<class> class SPHKernel>
    inline void from_json(const nlohmann::json &j, SolverConfig<Tvec, SPHKernel> &p) {
        using T       = SolverConfig<Tvec, SPHKernel>;
        using Tkernel = typename T::Kernel;
 
        // type checking
        if (j.contains("kernel_id")) {
 
            std::string kernel_id = j.at("kernel_id").get<std::string>();
 
            if (kernel_id != shambase::get_type_name<Tkernel>()) {
                shambase::throw_with_loc<std::runtime_error>(
                    "Invalid type to deserialize, wanted " + shambase::get_type_name<Tvec>()
                    + " but got " + kernel_id);
            }
        }
 
        if (j.contains("type_id")) {
 
            std::string type_id = j.at("type_id").get<std::string>();
 
            if (type_id != shambase::get_type_name<Tvec>()) {
                shambase::throw_with_loc<std::runtime_error>(
                    "Invalid type to deserialize, wanted " + shambase::get_type_name<Tvec>()
                    + " but got " + type_id);
            }
        }
 
        bool has_used_defaults  = false;
        bool has_updated_config = false;
 
        auto _get_to_if_contains = [&](const std::string &key, auto &value) {
            shamrock::get_to_if_contains(j, key, value, has_used_defaults);
        };
 
        auto _get_to_if_contains_fallbacks = [&](const std::string &key,
                                                 auto &value,
                                                 std::initializer_list<const char *> fallbacks) {
            shamrock::get_to_if_contains_fallbacks(
                j, key, value, fallbacks, has_used_defaults, has_updated_config);
        };
 
        _get_to_if_contains("scheduler_config", p.scheduler_conf);
 
        // actual data stored in the json
        _get_to_if_contains("gpart_mass", p.gpart_mass);
        _get_to_if_contains("cfl_config", p.cfl_config);
        _get_to_if_contains("unit_sys", p.unit_sys);
        _get_to_if_contains("time_state", p.time_state);
        _get_to_if_contains("show_cfl_detail", p.show_cfl_detail);
        _get_to_if_contains("mhd_config", p.mhd_config);
        _get_to_if_contains("dust_config", p.dust_config);
        _get_to_if_contains("self_grav_config", p.self_grav_config);
        _get_to_if_contains("tree_reduction_level", p.tree_reduction_level);
        _get_to_if_contains("use_two_stage_search", p.use_two_stage_search);
        _get_to_if_contains("show_neigh_stats", p.show_neigh_stats);
        _get_to_if_contains("combined_dtdiv_divcurlv_compute", p.combined_dtdiv_divcurlv_compute);
 
        // Try new names first, fall back to old names for backward compatibility
        _get_to_if_contains_fallbacks(
            "htol_up_coarse_cycle", p.htol_up_coarse_cycle, {"htol_up_tol"});
        _get_to_if_contains_fallbacks("htol_up_fine_cycle", p.htol_up_fine_cycle, {"htol_up_iter"});
 
        _get_to_if_contains("epsilon_h", p.epsilon_h);
        _get_to_if_contains("smoothing_length_config", p.smoothing_length_config);
        _get_to_if_contains("h_iter_per_subcycles", p.h_iter_per_subcycles);
        _get_to_if_contains("h_max_subcycles_count", p.h_max_subcycles_count);
        _get_to_if_contains("enable_particle_reordering", p.enable_particle_reordering);
        _get_to_if_contains("particle_reordering_step_freq", p.particle_reordering_step_freq);
        _get_to_if_contains("save_dt_to_fields", p.save_dt_to_fields);
        _get_to_if_contains("show_ghost_zone_graph", p.show_ghost_zone_graph);
        _get_to_if_contains("eos_config", p.eos_config);
        _get_to_if_contains("artif_viscosity", p.artif_viscosity);
        _get_to_if_contains("boundary_config", p.boundary_config);
        _get_to_if_contains("ext_force_config", p.ext_force_config);
        _get_to_if_contains("do_debug_dump", p.do_debug_dump);
        _get_to_if_contains("debug_dump_filename", p.debug_dump_filename);
        _get_to_if_contains("particle_killing", p.particle_killing);
 
        if (has_used_defaults || has_updated_config) {
            if (shamcomm::world_rank() == 0) {
                logger::info_ln(
                    "SPH::SolverConfig",
                    shamrock::log_json_changes(p, j, has_used_defaults, has_updated_config));
            }
        }
    }
 
} // namespace shammodels::sph