Solver.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 "SolverConfig.hpp"
#include "shambackends/vec.hpp"
#include "shammodels/common/SolverLog.hpp"
#include "shammodels/gsph/modules/GSPHGhostHandler.hpp"
#include "shammodels/gsph/modules/SolverStorage.hpp"
#include "shamrock/patch/PatchDataLayerLayout.hpp"
#include "shamrock/scheduler/ComputeField.hpp"
#include "shamrock/scheduler/InterfacesUtility.hpp"
#include "shamrock/scheduler/SerialPatchTree.hpp"
#include "shamrock/scheduler/ShamrockCtx.hpp"
#include "shamsolvergraph/edge/IDataEdgeSerializable.hpp"
#include "shamsys/legacy/log.hpp"
#include "shamtree/TreeTraversalCache.hpp"
#include <algorithm>
#include <memory>
#include <stdexcept>
#include <variant>
 
namespace shammodels::gsph {
 
    struct TimestepLog {
        i32 rank;
        f64 rate;
        u64 npart;
        f64 tcompute;
 
        inline f64 rate_sum() { return shamalgs::collective::allreduce_sum(rate); }
        inline u64 npart_sum() { return shamalgs::collective::allreduce_sum(npart); }
        inline f64 tcompute_max() { return shamalgs::collective::allreduce_max(tcompute); }
    };

    template<class Tvec, template<class> class SPHKernel>
    class Solver {
        public:
        using Tscal              = shambase::VecComponent<Tvec>;
        static constexpr u32 dim = shambase::VectorProperties<Tvec>::dimension;
        using Kernel             = SPHKernel<Tscal>;
 
        using Config = SolverConfig<Tvec, SPHKernel>;
 
        using u_morton = u32;
 
        static constexpr Tscal Rkern = Kernel::Rkern;
 
        ShamrockCtx &context;
        inline PatchScheduler &scheduler() { return shambase::get_check_ref(context.sched); }
 
        SolverStorage<Tvec, u_morton> storage{};
 
        Config solver_config;
        SolverLog solve_logs;

        inline Tscal &time_edge_value() {
            return scheduler()
                .synchronized_data
                .template get_edge_ref<shamrock::solvergraph::IDataEdgeSerializable<Tscal>>("time")
                .data;
        }

        inline Tscal &dt_edge_value() {
            return scheduler()
                .synchronized_data
                .template get_edge_ref<shamrock::solvergraph::IDataEdgeSerializable<Tscal>>("dt")
                .data;
        }
 
        inline Tscal get_time() { return time_edge_value(); }
        inline void set_time(Tscal t) { time_edge_value() = t; }
        inline Tscal get_dt() { return dt_edge_value(); }
        inline void set_next_dt(Tscal dt) { dt_edge_value() = dt; }

        inline void ensure_time_state_edges() {
            auto &sync    = scheduler().synchronized_data;
            auto names    = sync.get_edge_names();
            auto has_edge = [&](const std::string &name) {
                return std::find(names.begin(), names.end(), name) != names.end();
            };
 
            if (!has_edge("time")) {
                auto edge = sync.register_edge(
                    "time", shamrock::solvergraph::IDataEdgeSerializable<Tscal>("time", "t"));
                edge->data = 0;
            }
            if (!has_edge("dt")) {
                auto edge = sync.register_edge(
                    "dt", shamrock::solvergraph::IDataEdgeSerializable<Tscal>("dt", "dt"));
                edge->data = 0;
            }
        }
 
        inline void init_required_fields() { solver_config.set_layout(context.get_pdl_write()); }
 
        // Serial patch tree control
        void gen_serial_patch_tree();
        inline void reset_serial_patch_tree() { storage.serial_patch_tree.reset(); }
 
        // Ghost handling - use GSPH ghost handler with Newtonian field names
        using GhostHandle      = GSPHGhostHandler<Tvec>;
        using GhostHandleCache = typename GhostHandle::CacheMap;
 
        void gen_ghost_handler(Tscal time_val);
        inline void reset_ghost_handler() {
            shambase::get_check_ref(storage.ghost_handler).free_alloc();
        }
 
        void build_ghost_cache();
        void clear_ghost_cache();
 
        void merge_position_ghost();
 
        // Tree operations
        using RTree = typename Config::RTree;
        void build_merged_pos_trees();
        void clear_merged_pos_trees();
 
        void compute_presteps_rint();
        void reset_presteps_rint();
 
        void start_neighbors_cache();
        void reset_neighbors_cache();
 
        void gsph_prestep(Tscal time_val, Tscal dt);
 
        void apply_position_boundary(Tscal time_val);
 
        void do_predictor_leapfrog(Tscal dt);
 
        void init_ghost_layout();
 
        void communicate_merge_ghosts_fields();
        void reset_merge_ghosts_fields();
 
        void compute_eos_fields();
        void reset_eos_fields();

        void copy_eos_to_patchdata();

        void compute_density();

        void compute_gradients();
 
        void prepare_corrector();

        void update_derivs();

        Tscal compute_dt_cfl();
 
        bool apply_corrector(Tscal dt, u64 Npart_all);
 
        void update_sync_load_values();
 
        Solver(ShamrockCtx &context) : context(context) {}
 
        void init_solver_graph();
 
        void vtk_do_dump(std::string filename, bool add_patch_world_id);
 
        inline void print_timestep_logs() {
            if (shamcomm::world_rank() == 0) {
                logger::info_ln(
                    "GSPH", "iteration since start :", solve_logs.get_iteration_count());
                logger::info_ln(
                    "GSPH", "time since start :", shambase::details::get_wtime(), "(s)");
            }
        }
 
        TimestepLog evolve_once();
 
        Tscal evolve_once_time_expl(Tscal t_current, Tscal dt_input) {
            set_time(t_current);
            set_next_dt(dt_input);
            evolve_once();
            return get_dt();
        }
 
        inline bool evolve_until(Tscal target_time, i32 niter_max = -1) {
            auto step = [&]() {
                Tscal dt = get_dt();
                Tscal t  = get_time();
 
                if (t > target_time) {
                    throw shambase::make_except_with_loc<std::invalid_argument>(
                        "the target time is higher than the current time");
                }
 
                if (t + dt > target_time) {
                    set_next_dt(target_time - t);
                }
                evolve_once();
            };
 
            i32 iter_count = 0;
 
            while (get_time() < target_time) {
                step();
                iter_count++;
 
                if ((iter_count >= niter_max) && (niter_max != -1)) {
                    logger::info_ln("GSPH", "stopping evolve until because of niter =", iter_count);
                    return false;
                }
            }
 
            print_timestep_logs();
 
            return true;
        }
    };
 
} // namespace shammodels::gsph