From f2121be0f783f17bf898b7970d4a58058b0cb124 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Thu, 30 Jul 2026 20:29:16 +0900 Subject: [PATCH 1/7] everything except sink update --- .../common/setup/GeneratorMCDisc.hpp | 173 ++++++++++++ .../src/modules/setup/GeneratorMCDisc.cpp | 185 +++++++++++++ src/shammodels/gsph/CMakeLists.txt | 3 + .../gsph/include/shammodels/gsph/Model.hpp | 17 ++ .../include/shammodels/gsph/SolverConfig.hpp | 25 ++ .../gsph/modules/ComputeLoadBalanceValue.hpp | 50 ++++ .../shammodels/gsph/modules/GSPHSetup.hpp | 92 +++++++ .../shammodels/gsph/modules/SolverStorage.hpp | 3 + .../gsph/modules/setup/GeneratorMCDisc.hpp | 174 ++++++++++++ .../gsph/modules/setup/IGSPHSetupNode.hpp | 143 ++++++++++ .../src/modules/ComputeLoadBalanceValue.cpp | 39 +++ src/shammodels/gsph/src/modules/GSPHSetup.cpp | 184 +++++++++++++ .../gsph/src/modules/GeneratorMCDisc.cpp | 181 ++++++++++++ src/shammodels/gsph/src/pyGSPHModel.cpp | 260 +++++++++++++++++- 14 files changed, 1526 insertions(+), 3 deletions(-) create mode 100644 src/shammodels/common/include/shammodels/common/setup/GeneratorMCDisc.hpp create mode 100644 src/shammodels/common/src/modules/setup/GeneratorMCDisc.cpp create mode 100644 src/shammodels/gsph/include/shammodels/gsph/modules/ComputeLoadBalanceValue.hpp create mode 100644 src/shammodels/gsph/include/shammodels/gsph/modules/GSPHSetup.hpp create mode 100644 src/shammodels/gsph/include/shammodels/gsph/modules/setup/GeneratorMCDisc.hpp create mode 100644 src/shammodels/gsph/include/shammodels/gsph/modules/setup/IGSPHSetupNode.hpp create mode 100644 src/shammodels/gsph/src/modules/ComputeLoadBalanceValue.cpp create mode 100644 src/shammodels/gsph/src/modules/GSPHSetup.cpp create mode 100644 src/shammodels/gsph/src/modules/GeneratorMCDisc.cpp diff --git a/src/shammodels/common/include/shammodels/common/setup/GeneratorMCDisc.hpp b/src/shammodels/common/include/shammodels/common/setup/GeneratorMCDisc.hpp new file mode 100644 index 0000000000..b7258c57a8 --- /dev/null +++ b/src/shammodels/common/include/shammodels/common/setup/GeneratorMCDisc.hpp @@ -0,0 +1,173 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2026 Timothée David--Cléris +// 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 + +/** + * @file GeneratorMCDisc.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief + * + */ + +#include "shambase/constants.hpp" +#include "shamalgs/collective/InvariantParallelGenerator.hpp" +#include "shamalgs/collective/indexing.hpp" +#include "shamrock/scheduler/ShamrockCtx.hpp" + +namespace shammodels::sph::modules { + + template class SPHKernel, class TSolverConfig, class TSetupNodeBase> + class GeneratorMCDisc : public TSetupNodeBase { + using Tscal = shambase::VecComponent; + static constexpr u32 dim = shambase::VectorProperties::dimension; + using Kernel = SPHKernel; + + using Config = TSolverConfig; + + ShamrockCtx &context; + Config &solver_config; + + struct DiscOutput { + sycl::vec pos; + Tscal rho; + }; + + Tscal pmass; + + class DiscIterator; + DiscIterator generator; + Tscal init_h_factor; + + std::function vel_profile; + std::function cs_profile; + + static DiscIterator make_generator( + Tscal part_mass, + Tscal disc_mass, + Tscal r_in, + Tscal r_out, + std::function sigma_profile, + std::function H_profile, + std::mt19937_64 eng) { + return DiscIterator(part_mass, disc_mass, r_in, r_out, sigma_profile, H_profile, eng); + } + + public: + GeneratorMCDisc( + ShamrockCtx &context, + Config &solver_config, + Tscal part_mass, + Tscal disc_mass, + Tscal r_in, + Tscal r_out, + std::function sigma_profile, + std::function H_profile, + std::function vel_profile, + std::function cs_profile, + std::mt19937_64 eng, + Tscal init_h_factor) + : context(context), solver_config(solver_config), + generator( + make_generator(part_mass, disc_mass, r_in, r_out, sigma_profile, H_profile, eng)), + init_h_factor(init_h_factor), pmass(part_mass), vel_profile(vel_profile), + cs_profile(cs_profile) {} + + bool is_done(); + + shamrock::patch::PatchDataLayer next_n(u32 nmax); + + std::string get_name() { return "GeneratorMCDisc"; } + ISPHSetupNode_Dot get_dot_subgraph() { return ISPHSetupNode_Dot{get_name(), 0, {}}; } + }; + +} // namespace shammodels::sph::modules + +template class SPHKernel, class TConfig, class TSetupNodeBase> +class shammodels::common::GeneratorMCDisc::DiscIterator { + + bool done = false; + u64 current_index = 0; + + Tscal part_mass; + Tscal disc_mass; + u64 Npart; + + Tscal r_in; + Tscal r_out; + std::function sigma_profile; + std::function H_profile; + + shamalgs::collective::InvariantParallelGenerator generator; + + static constexpr Tscal _2pi = 2 * shambase::constants::pi; + + Tscal f_func(Tscal r) { return r * sigma_profile(r); } + + DiscOutput next(u64 seed); + + public: + DiscIterator( + Tscal part_mass, + Tscal disc_mass, + Tscal r_in, + Tscal r_out, + std::function sigma_profile, + std::function H_profile, + std::mt19937_64 eng) + : DiscIterator( + part_mass, + disc_mass, + r_in, + r_out, + sigma_profile, + H_profile, + eng, + disc_mass / part_mass) {} + + DiscIterator( + Tscal part_mass, + Tscal disc_mass, + Tscal r_in, + Tscal r_out, + std::function sigma_profile, + std::function H_profile, + std::mt19937_64 eng, + u64 Npart) + : part_mass(part_mass), disc_mass(disc_mass), Npart(Npart), r_in(r_in), r_out(r_out), + sigma_profile(sigma_profile), H_profile(H_profile), generator(eng, Npart), + current_index(0) { + + shamlog_debug_ln( + "GeneratorMCDisc", + "part_mass", + part_mass, + "disc_mass", + disc_mass, + "r_in", + r_in, + "r_out", + r_out, + "Npart", + Npart); + } + + inline bool is_done() { + return generator.is_done(); + } // just to make sure the result is not tempered with + + inline std::vector next_n(u64 nmax) { + std::vector seeds = generator.next_n(nmax); + std::vector ret{}; + for (u64 seed : seeds) { + ret.push_back(next(seed)); + } + return ret; + } +}; diff --git a/src/shammodels/common/src/modules/setup/GeneratorMCDisc.cpp b/src/shammodels/common/src/modules/setup/GeneratorMCDisc.cpp new file mode 100644 index 0000000000..8144a2cd2c --- /dev/null +++ b/src/shammodels/common/src/modules/setup/GeneratorMCDisc.cpp @@ -0,0 +1,185 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2026 Timothée David--Cléris +// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1 +// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information +// +// -------------------------------------------------------// + +/** + * @file GeneratorMCDisc.cpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) + * @brief + * + */ + +#include "shambase/constants.hpp" +#include "shamalgs/collective/indexing.hpp" +#include "shamalgs/random.hpp" +#include "shammodels/common/setup/GeneratorMCDisc.hpp" +#include "shammodels/sph/math/density.hpp" + +template< + class Tvec, + template class SPHKernel, + class TConfig, + class TSetupNodeBase + >> auto shammodels::sph::modules::GeneratorMCDisc::DiscIterator::next( + u64 seed) -> DiscOutput { + + std::mt19937_64 eng_local(seed); // ensure that 1 part = 1 random draw + + Tscal fmax = f_func(r_out); + + auto find_r = [&]() { + while (true) { + Tscal u2 = shamalgs::primitives::mock_value(eng_local, 0, fmax); + Tscal r = shamalgs::primitives::mock_value(eng_local, r_in, r_out); + if (u2 < f_func(r)) { + return r; + } + } + }; + + auto theta = shamalgs::primitives::mock_value(eng_local, 0, _2pi); + auto Gauss = shamalgs::random::mock_gaussian(eng_local); + + // depends on sigma profile + Tscal r = find_r(); + Tscal sigma = sigma_profile(r); + + // depends on H profile & sigma profile (through r) + Tscal H = H_profile(r); + Tscal z = H * Gauss; + + auto pos = sycl::vec{r * sycl::cos(theta), r * sycl::sin(theta), z}; + + // extrapolate the density from sigma profile + Tscal fs = 1; + Tscal rho = (sigma * fs) * sycl::exp(-z * z / (2 * H * H)); + + DiscOutput out{.pos = pos, .rho = rho}; + + // increase counter + check if finished + current_index++; + if (current_index == Npart) { + done = true; + } + + return out; +} + +template class SPHKernel> +bool shammodels::sph::modules::GeneratorMCDisc::is_done() { + return generator.is_done(); +} + +template class SPHKernel> +shamrock::patch::PatchDataLayer shammodels::sph::modules::GeneratorMCDisc::next_n( + u32 nmax) { + + using namespace shamrock::patch; + PatchScheduler &sched = shambase::get_check_ref(context.sched); + + std::vector pos_data; + + // Fill pos_data if the scheduler has some patchdata in this rank + if (!generator.is_done()) { + u64 loc_gen_count = nmax; + pos_data = generator.next_n(loc_gen_count); + } + + // extract data from disc output + std::vector vec_pos; + std::vector vec_rho; + + vec_pos.reserve(pos_data.size()); + vec_rho.reserve(pos_data.size()); + + for (DiscOutput o : pos_data) { + vec_pos.push_back(o.pos); + vec_rho.push_back(o.rho); + } + + // compute the hpart from the rho + std::vector vec_h; + vec_h.reserve(pos_data.size()); + for (Tscal rho : vec_rho) { + vec_h.push_back(shamrock::sph::h_rho(pmass, rho, Kernel::hfactd) * init_h_factor); + } + + // compute velocities + std::vector vec_vel; + vec_vel.reserve(pos_data.size()); + for (size_t i = 0; i < vec_pos.size(); i++) { + Tvec vel = vel_profile(vec_pos[i]); + vec_vel.push_back(vel); + } + + // compute the cs + bool need_cs = solver_config.is_eos_locally_isothermal(); + + std::vector vec_cs; + if (need_cs) { + if (!cs_profile) { + throw shambase::make_except_with_loc( + "With this EOS you need to provide a cs_profile"); + } + vec_cs.reserve(pos_data.size()); + for (size_t i = 0; i < vec_pos.size(); i++) { + Tscal cs = cs_profile(vec_pos[i]); + vec_cs.push_back(cs); + } + } + + // Make a patchdata from pos_data + PatchDataLayer tmp(sched.get_layout_ptr_old()); + if (!pos_data.empty()) { + tmp.resize(pos_data.size()); + tmp.fields_raz(); + + { + u32 len = pos_data.size(); + PatchDataField &f + = tmp.get_field(sched.pdl_old().get_field_idx("xyz")); + sycl::buffer buf(vec_pos.data(), len); + f.override(buf, len); + } + + { + u32 len = pos_data.size(); + PatchDataField &f + = tmp.get_field(sched.pdl_old().get_field_idx("vxyz")); + sycl::buffer buf(vec_vel.data(), len); + f.override(buf, len); + } + { + u32 len = vec_pos.size(); + PatchDataField &f + = tmp.get_field(sched.pdl_old().get_field_idx("hpart")); + sycl::buffer buf(vec_h.data(), len); + f.override(buf, len); + } + + if (need_cs) { + u32 len = vec_pos.size(); + PatchDataField &f + = tmp.get_field(sched.pdl_old().get_field_idx("soundspeed")); + sycl::buffer buf(vec_cs.data(), len); + f.override(buf, len); + } + } + + return tmp; +} + +using namespace shammath; +template class shammodels::sph::modules::GeneratorMCDisc; +template class shammodels::sph::modules::GeneratorMCDisc; +template class shammodels::sph::modules::GeneratorMCDisc; + +template class shammodels::sph::modules::GeneratorMCDisc; +template class shammodels::sph::modules::GeneratorMCDisc; +template class shammodels::sph::modules::GeneratorMCDisc; diff --git a/src/shammodels/gsph/CMakeLists.txt b/src/shammodels/gsph/CMakeLists.txt index 10123c00f3..0a0960b89b 100644 --- a/src/shammodels/gsph/CMakeLists.txt +++ b/src/shammodels/gsph/CMakeLists.txt @@ -18,6 +18,9 @@ set(Sources src/modules/UpdateDerivs.cpp src/modules/GSPHGhostHandler.cpp src/modules/io/VTKDump.cpp + src/modules/GSPHSetup.cpp + src/modules/GeneratorMCDisc.cpp + src/modules/ComputeLoadBalanceValue.cpp ) if(SHAMROCK_USE_SHARED_LIB) diff --git a/src/shammodels/gsph/include/shammodels/gsph/Model.hpp b/src/shammodels/gsph/include/shammodels/gsph/Model.hpp index 17238ccf8a..aa3396b09e 100644 --- a/src/shammodels/gsph/include/shammodels/gsph/Model.hpp +++ b/src/shammodels/gsph/include/shammodels/gsph/Model.hpp @@ -35,6 +35,7 @@ #include "shamcomm/logs.hpp" #include "shammodels/common/setup/generators.hpp" #include "shammodels/gsph/Solver.hpp" +#include "shammodels/gsph/modules/GSPHSetup.hpp" #include "shammodels/sph/math/density.hpp" #include "shamrock/io/ShamrockDump.hpp" #include "shamrock/patch/PatchDataLayer.hpp" @@ -132,6 +133,22 @@ namespace shammodels::gsph { void add_cube_fcc_3d(Tscal dr, std::pair _box); void add_cube_hcp_3d(Tscal dr, std::pair _box); + inline std::unique_ptr> get_setup() { + return std::make_unique>( + ctx, solver.solver_config, solver.storage); + } + + inline void add_sink(Tscal mass, Tvec pos, Tvec velocity, Tscal accretion_radius) { + if (solver.storage.sinks.is_empty()) { + solver.storage.sinks.set({}); + } + + shamlog_debug_ln("SPH", "add sink :", mass, pos, velocity, accretion_radius); + + solver.storage.sinks.get().push_back( + {pos, velocity, {}, {}, mass, {}, accretion_radius}); + } + //////////////////////////////////////////////////////////////////////////////////////////// // Field manipulation //////////////////////////////////////////////////////////////////////////////////////////// diff --git a/src/shammodels/gsph/include/shammodels/gsph/SolverConfig.hpp b/src/shammodels/gsph/include/shammodels/gsph/SolverConfig.hpp index b9a0faa001..aa3238fe07 100644 --- a/src/shammodels/gsph/include/shammodels/gsph/SolverConfig.hpp +++ b/src/shammodels/gsph/include/shammodels/gsph/SolverConfig.hpp @@ -192,6 +192,31 @@ struct shammodels::gsph::SolverConfig { inline void set_eos_isothermal(Tscal cs) { eos_config.set_isothermal(cs); } + /** + * @brief Set the EOS configuration to a locally isothermal equation of state fromFarris 2014 + * + * @param cs0 Soundspeed at the reference radius + * @param q Power exponent of the soundspeed profile + * @param r0 Reference radius + */ + inline void set_eos_locally_isothermalFA2014(Tscal h_over_r) { + eos_config.set_locally_isothermalFA2014(h_over_r); + } + + /** + * @brief Set the EOS configuration to a locally isothermal equation of state from Farris 2014 + * extended to q != 1/2 + * + * @param cs0 Soundspeed at the reference radius + * @param q Power exponent of the soundspeed profile + * @param r0 Reference radius + * @param n_sinks Number of sinks to consider for the equation of state + */ + 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); + } + ////////////////////////////////////////////////////////////////////////////////////////////// // EOS Config (END) ////////////////////////////////////////////////////////////////////////////////////////////// diff --git a/src/shammodels/gsph/include/shammodels/gsph/modules/ComputeLoadBalanceValue.hpp b/src/shammodels/gsph/include/shammodels/gsph/modules/ComputeLoadBalanceValue.hpp new file mode 100644 index 0000000000..236c0f0abe --- /dev/null +++ b/src/shammodels/gsph/include/shammodels/gsph/modules/ComputeLoadBalanceValue.hpp @@ -0,0 +1,50 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2026 Timothée David--Cléris +// 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 + +/** + * @file ComputeLoadBalanceValue.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief + * + */ + +#include "shambackends/typeAliasVec.hpp" +#include "shambackends/vec.hpp" +#include "shammodels/gsph/SolverConfig.hpp" +#include "shammodels/gsph/modules/SolverStorage.hpp" +#include "shamrock/scheduler/ShamrockCtx.hpp" + +namespace shammodels::gsph::modules { + + template class SPHKernel> + class ComputeLoadBalanceValue { + public: + using Tscal = shambase::VecComponent; + static constexpr u32 dim = shambase::VectorProperties::dimension; + using Kernel = SPHKernel; + + using Config = SolverConfig; + using Storage = SolverStorage; + + ShamrockCtx &context; + Config &solver_config; + Storage &storage; + + ComputeLoadBalanceValue(ShamrockCtx &context, Config &solver_config, Storage &storage) + : context(context), solver_config(solver_config), storage(storage) {} + + void update_load_balancing(); + + private: + inline PatchScheduler &scheduler() { return shambase::get_check_ref(context.sched); } + }; + +} // namespace shammodels::gsph::modules diff --git a/src/shammodels/gsph/include/shammodels/gsph/modules/GSPHSetup.hpp b/src/shammodels/gsph/include/shammodels/gsph/modules/GSPHSetup.hpp new file mode 100644 index 0000000000..f59db0b6f6 --- /dev/null +++ b/src/shammodels/gsph/include/shammodels/gsph/modules/GSPHSetup.hpp @@ -0,0 +1,92 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2026 Timothée David--Cléris +// 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 + +/** + * @file GSPHSetup.hpp + * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) + * @brief + * + */ + +#include "shambackends/typeAliasVec.hpp" +#include "shambackends/vec.hpp" +#include "shammodels/gsph/SolverConfig.hpp" +#include "shammodels/gsph/modules/SolverStorage.hpp" +#include "shammodels/gsph/modules/setup/IGSPHSetupNode.hpp" +#include "shamrock/scheduler/ShamrockCtx.hpp" +#include + +namespace shammodels::gsph::modules { + + template class SPHKernel> + class GSPHSetup { + public: + using Tscal = shambase::VecComponent; + static constexpr u32 dim = shambase::VectorProperties::dimension; + using Kernel = SPHKernel; + + using Config = SolverConfig; + using Storage = SolverStorage; + + ShamrockCtx &context; + Config &solver_config; + Storage &storage; + + GSPHSetup(ShamrockCtx &context, Config &solver_config, Storage &storage) + : context(context), solver_config(solver_config), storage(storage) {} + + void apply_setup( + SetupNodePtr setup, + bool part_reordering, + std::optional insert_step = std::nullopt); + + std::shared_ptr make_generator_disc_mc( + Tscal part_mass, + Tscal disc_mass, + Tscal r_in, + Tscal r_out, + std::function sigma_profile, + std::function H_profile, + std::function vel_profile, + std::function cs_profile, + std::mt19937_64 eng, + Tscal init_h_factor); + + std::shared_ptr make_generator_from_context(ShamrockCtx &context_other); + + std::shared_ptr make_combiner_add( + SetupNodePtr parent1, SetupNodePtr parent2); + + std::shared_ptr make_modifier_warp_disc( + SetupNodePtr parent, Tscal Rwarp, Tscal Hwarp, Tscal inclination, Tscal posangle); + + std::shared_ptr make_modifier_custom_warp( + SetupNodePtr parent, + std::function inc_profile, + std::function psi_profile, + std::function k_profile); + + std::shared_ptr make_modifier_add_offset( + SetupNodePtr parent, Tvec offset_postion, Tvec offset_velocity); + + std::shared_ptr make_modifier_filter( + SetupNodePtr parent, std::function filter); + + std::shared_ptr make_modifier_split_part( + SetupNodePtr parent, u64 n_split, u64 seed, Tscal h_scaling); + + private: + inline PatchScheduler &scheduler() { return shambase::get_check_ref(context.sched); } + + u64 injected_parts = 0; + }; + +} // namespace shammodels::gsph::modules diff --git a/src/shammodels/gsph/include/shammodels/gsph/modules/SolverStorage.hpp b/src/shammodels/gsph/include/shammodels/gsph/modules/SolverStorage.hpp index 3217443eea..66c65276d1 100644 --- a/src/shammodels/gsph/include/shammodels/gsph/modules/SolverStorage.hpp +++ b/src/shammodels/gsph/include/shammodels/gsph/modules/SolverStorage.hpp @@ -29,6 +29,7 @@ #include "shambackends/vec.hpp" #include "shammodels/gsph/modules/GSPHGhostHandler.hpp" #include "shammodels/gsph/solvergraph/GhostHandlerEdge.hpp" +#include "shammodels/sph/SinkPartStruct.hpp" #include "shammodels/sph/solvergraph/NeighCache.hpp" #include "shamrock/scheduler/SerialPatchTree.hpp" #include "shamrock/scheduler/ShamrockCtx.hpp" @@ -136,6 +137,8 @@ namespace shammodels::gsph { Component> old_axyz; Component> old_duint; + Component>> sinks; + /// Timing statistics struct Timings { f64 interface = 0; diff --git a/src/shammodels/gsph/include/shammodels/gsph/modules/setup/GeneratorMCDisc.hpp b/src/shammodels/gsph/include/shammodels/gsph/modules/setup/GeneratorMCDisc.hpp new file mode 100644 index 0000000000..1e8587818e --- /dev/null +++ b/src/shammodels/gsph/include/shammodels/gsph/modules/setup/GeneratorMCDisc.hpp @@ -0,0 +1,174 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2026 Timothée David--Cléris +// 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 + +/** + * @file GeneratorMCDisc.hpp + * @brief + * + */ + +#include "shambase/constants.hpp" +#include "shamalgs/collective/InvariantParallelGenerator.hpp" +#include "shamalgs/collective/indexing.hpp" +#include "shammodels/gsph/SolverConfig.hpp" +#include "shammodels/gsph/modules/setup/IGSPHSetupNode.hpp" +#include "shamrock/scheduler/ShamrockCtx.hpp" + +namespace shammodels::gsph::modules { + + template class SPHKernel> + class GeneratorMCDisc : public IGSPHSetupNode { + using Tscal = shambase::VecComponent; + static constexpr u32 dim = shambase::VectorProperties::dimension; + using Kernel = SPHKernel; + + using Config = SolverConfig; + + ShamrockCtx &context; + Config &solver_config; + + struct DiscOutput { + sycl::vec pos; + Tscal rho; + }; + + Tscal pmass; + + class DiscIterator; + DiscIterator generator; + Tscal init_h_factor; + + std::function vel_profile; + std::function cs_profile; + + static DiscIterator make_generator( + Tscal part_mass, + Tscal disc_mass, + Tscal r_in, + Tscal r_out, + std::function sigma_profile, + std::function H_profile, + std::mt19937_64 eng) { + return DiscIterator(part_mass, disc_mass, r_in, r_out, sigma_profile, H_profile, eng); + } + + public: + GeneratorMCDisc( + ShamrockCtx &context, + Config &solver_config, + Tscal part_mass, + Tscal disc_mass, + Tscal r_in, + Tscal r_out, + std::function sigma_profile, + std::function H_profile, + std::function vel_profile, + std::function cs_profile, + std::mt19937_64 eng, + Tscal init_h_factor) + : context(context), solver_config(solver_config), + generator( + make_generator(part_mass, disc_mass, r_in, r_out, sigma_profile, H_profile, eng)), + init_h_factor(init_h_factor), pmass(part_mass), vel_profile(vel_profile), + cs_profile(cs_profile) {} + + bool is_done(); + + shamrock::patch::PatchDataLayer next_n(u32 nmax); + + std::string get_name() { return "GeneratorMCDisc"; } + IGSPHSetupNode_Dot get_dot_subgraph() { return IGSPHSetupNode_Dot{get_name(), 0, {}}; } + }; + +} // namespace shammodels::gsph::modules + +template class SPHKernel> +class shammodels::gsph::modules::GeneratorMCDisc::DiscIterator { + + bool done = false; + u64 current_index = 0; + + Tscal part_mass; + Tscal disc_mass; + u64 Npart; + + Tscal r_in; + Tscal r_out; + std::function sigma_profile; + std::function H_profile; + + shamalgs::collective::InvariantParallelGenerator generator; + + static constexpr Tscal _2pi = 2 * shambase::constants::pi; + + Tscal f_func(Tscal r) { return r * sigma_profile(r); } + + DiscOutput next(u64 seed); + + public: + DiscIterator( + Tscal part_mass, + Tscal disc_mass, + Tscal r_in, + Tscal r_out, + std::function sigma_profile, + std::function H_profile, + std::mt19937_64 eng) + : DiscIterator( + part_mass, + disc_mass, + r_in, + r_out, + sigma_profile, + H_profile, + eng, + disc_mass / part_mass) {} + + DiscIterator( + Tscal part_mass, + Tscal disc_mass, + Tscal r_in, + Tscal r_out, + std::function sigma_profile, + std::function H_profile, + std::mt19937_64 eng, + u64 Npart) + : part_mass(part_mass), disc_mass(disc_mass), Npart(Npart), r_in(r_in), r_out(r_out), + sigma_profile(sigma_profile), H_profile(H_profile), generator(eng, Npart), + current_index(0) { + + shamlog_debug_ln( + "GeneratorMCDisc", + "part_mass", + part_mass, + "disc_mass", + disc_mass, + "r_in", + r_in, + "r_out", + r_out, + "Npart", + Npart); + } + + inline bool is_done() { + return generator.is_done(); + } // just to make sure the result is not tempered with + + inline std::vector next_n(u64 nmax) { + std::vector seeds = generator.next_n(nmax); + std::vector ret{}; + for (u64 seed : seeds) { + ret.push_back(next(seed)); + } + return ret; + } +}; diff --git a/src/shammodels/gsph/include/shammodels/gsph/modules/setup/IGSPHSetupNode.hpp b/src/shammodels/gsph/include/shammodels/gsph/modules/setup/IGSPHSetupNode.hpp new file mode 100644 index 0000000000..62dcea693f --- /dev/null +++ b/src/shammodels/gsph/include/shammodels/gsph/modules/setup/IGSPHSetupNode.hpp @@ -0,0 +1,143 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2026 Timothée David--Cléris +// 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 + +/** + * @file IGSPHSetupNode.hpp + * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) + * @brief + * + */ + +#include "shambase/aliases_int.hpp" +#include "shamrock/patch/PatchDataLayer.hpp" +#include +#include + +namespace shammodels::gsph::modules { + + /** + * @brief This struct is used to generate a dot graph of the setup tree + * + * It is composed of a name, a type and a vector of inputs. + */ + struct IGSPHSetupNode_Dot { + std::string name; + u32 type; + std::vector inputs; + + /** + * @brief This function generate a dot graph for the setup tree + * + * This function is used to generate a dot graph that describes the + * setup tree. It takes a counter and a string as input, and update + * the counter and the string to generate the dot graph. + * + * @param counter a counter that is used to generate the node id in + * the dot graph + * @param out the string that will be updated to contain the dot graph + * @return the new value of the counter + */ + u32 add_node(u32 &counter, std::string &out) { + + std::vector inputs_id{}; + for (auto &in : inputs) { + inputs_id.push_back(in.add_node(counter, out)); + } + + u32 counter_val = counter; + counter++; + + out += "node_" + std::to_string(counter_val) + " [label=\"" + name + "\"];\n"; + + for (auto i : inputs_id) { + out += "node_" + std::to_string(i) + " -> node_" + std::to_string(counter_val) + + ";\n"; + } + + return counter_val; + } + }; + + /** + * @class IGSPHSetupNode + * @brief This class is an interface that all SPH setup nodes must implement. + * It describe an operation associated to a node in the setup tree. + */ + class IGSPHSetupNode { + public: + /** + * @brief This function return true if the setup is done + * + * @return true if done, false otherwise + */ + virtual bool is_done() = 0; + + /** + * @brief This function generate patchdata with at most nmax per MPI ranks + * This function is always assumed as called by every ranks simultaneously + * + * @param nmax + * @return shamrock::patch::PatchData + */ + virtual shamrock::patch::PatchDataLayer next_n(u32 nmax) = 0; + + /** + * @brief Get the name of the node + * @return The name of the node + */ + virtual std::string get_name() = 0; + + /** + * @brief Get a dot subgraph describing the node and its childrens (recursively) + * + * This function should return a IGSPHSetupNode_Dot object which contains + * all the information needed to generate a dot graph for the node and + * its children. + * + * @return A IGSPHSetupNode_Dot object + */ + virtual IGSPHSetupNode_Dot get_dot_subgraph() = 0; + + /** + * @brief Virtual destructor for the IGSPHSetupNode class + */ + virtual ~IGSPHSetupNode() = default; + + /** + * @brief Generate a dot graph for the setup tree + * + * This function returns a string containing a dot graph that describes the + * setup tree. + * + * @return A string containing a dot graph + */ + std::string get_dot() { + std::string out; + + out += "digraph G {\n"; + out += "rankdir=LR;\n"; + + u32 counter = 0; + u32 final_node = get_dot_subgraph().add_node(counter, out); + + out += "node_" + std::to_string(counter + 1) + " [label=\"Simulation\"];\n"; + out += "node_" + std::to_string(final_node) + " -> node_" + std::to_string(counter + 1) + + ";\n"; + + out += "}\n"; + return out; + } + }; + + /// Alias for a shared pointer to an IGSPHSetupNode + using SetupNodePtr = std::shared_ptr; + +} // namespace shammodels::gsph::modules diff --git a/src/shammodels/gsph/src/modules/ComputeLoadBalanceValue.cpp b/src/shammodels/gsph/src/modules/ComputeLoadBalanceValue.cpp new file mode 100644 index 0000000000..c98d32c8c4 --- /dev/null +++ b/src/shammodels/gsph/src/modules/ComputeLoadBalanceValue.cpp @@ -0,0 +1,39 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2026 Timothée David--Cléris +// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1 +// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information +// +// -------------------------------------------------------// + +/** + * @file ComputeLoadBalanceValue.cpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) + * @brief + * + */ + +#include "shammodels/gsph/modules/ComputeLoadBalanceValue.hpp" +#include "shammath/sphkernels.hpp" +#include "shamsys/legacy/log.hpp" + +template class SPHKernel> +void shammodels::gsph::modules::ComputeLoadBalanceValue::update_load_balancing() { + StackEntry stack_loc{}; + + shamlog_debug_ln("ComputeLoadBalanceValue", "update load balancing"); + scheduler().update_local_load_value([&](shamrock::patch::Patch p) { + return scheduler().patch_data.owned_data.get(p.id_patch).get_obj_cnt(); + }); +} + +using namespace shammath; +template class shammodels::gsph::modules::ComputeLoadBalanceValue; +template class shammodels::gsph::modules::ComputeLoadBalanceValue; +template class shammodels::gsph::modules::ComputeLoadBalanceValue; + +template class shammodels::gsph::modules::ComputeLoadBalanceValue; +template class shammodels::gsph::modules::ComputeLoadBalanceValue; +template class shammodels::gsph::modules::ComputeLoadBalanceValue; diff --git a/src/shammodels/gsph/src/modules/GSPHSetup.cpp b/src/shammodels/gsph/src/modules/GSPHSetup.cpp new file mode 100644 index 0000000000..23aa223df6 --- /dev/null +++ b/src/shammodels/gsph/src/modules/GSPHSetup.cpp @@ -0,0 +1,184 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2026 Timothée David--Cléris +// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1 +// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information +// +// -------------------------------------------------------// + +/** + * @file GSPHSetup.cpp + * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) + * @brief + * + */ + +#include "shambase/DistributedData.hpp" +#include "shambase/aliases_int.hpp" +#include "shambase/memory.hpp" +#include "shambase/string.hpp" +#include "shambase/tabulate.hpp" +#include "shamalgs/collective/are_all_rank_true.hpp" +#include "shamalgs/primitives/is_all_true.hpp" +#include "shambackends/DeviceBuffer.hpp" +#include "shambackends/SyclMpiTypes.hpp" +#include "shambackends/kernel_call.hpp" +#include "shamcomm/logs.hpp" +#include "shamcomm/worldInfo.hpp" +#include "shamcomm/wrapper.hpp" +#include "shammodels/gsph/modules/ComputeLoadBalanceValue.hpp" +#include "shammodels/gsph/modules/GSPHSetup.hpp" +#include "shammodels/gsph/modules/setup/GeneratorMCDisc.hpp" +#include "shammodels/sph/modules/ParticleReordering.hpp" +#include "shamrock/patch/PatchDataLayer.hpp" +#include "shamrock/scheduler/DataInserterUtility.hpp" +#include "shamsys/NodeInstance.hpp" +#include +#include + +template class SPHKernel> +inline std::shared_ptr shammodels::gsph::modules:: + GSPHSetup::make_generator_disc_mc( + Tscal part_mass, + Tscal disc_mass, + Tscal r_in, + Tscal r_out, + std::function sigma_profile, + std::function H_profile, + std::function vel_profile, + std::function cs_profile, + std::mt19937_64 eng, + Tscal init_h_factor) { + return std::shared_ptr(new gsph::modules::GeneratorMCDisc( + context, + solver_config, + part_mass, + disc_mass, + r_in, + r_out, + sigma_profile, + H_profile, + vel_profile, + cs_profile, + eng, + init_h_factor)); +} + +template class SPHKernel> +void shammodels::gsph::modules::GSPHSetup::apply_setup( + SetupNodePtr setup, bool part_reordering, std::optional insert_step) { + + if (!bool(setup)) { + shambase::throw_with_loc("The setup shared pointer is empty"); + } + + shambase::Timer time_setup; + time_setup.start(); + StackEntry stack_loc{}; + + PatchScheduler &sched = shambase::get_check_ref(context.sched); + + auto compute_load = [&]() { + modules::ComputeLoadBalanceValue(context, solver_config, storage) + .update_load_balancing(); + }; + + auto has_pdat = [&]() { + bool ret = false; + using namespace shamrock::patch; + sched.for_each_local_patchdata([&](const Patch &p, PatchDataLayer &pdat) { + ret = true; + }); + return ret; + }; + + shamrock::DataInserterUtility inserter(sched); + u32 _insert_step = sched.crit_patch_split * 8; + if (bool(insert_step)) { + _insert_step = insert_step.value(); + } + + while (!setup->is_done()) { + + shamrock::patch::PatchDataLayer pdat = setup->next_n((has_pdat()) ? _insert_step : 0); + + u64 injected + = inserter.push_patch_data(pdat, "xyz", sched.crit_patch_split * 8, compute_load); + + injected_parts += injected; + } + + u32 final_balancing_steps = 3; + for (u32 i = 0; i < final_balancing_steps; i++) { + ON_RANK_0( + logger::info_ln( + "SPH setup", "Final load balancing step", i, "of", final_balancing_steps)); + inserter.balance_load(compute_load); + } + + time_setup.stop(); + if (shamcomm::world_rank() == 0) { + logger::info_ln("SPH setup", "the setup took :", time_setup.elapsed_sec(), "s"); + } +} + +struct SetupLog { + struct State { + std::vector count_per_rank; + std::vector> msg_list; + } state; + + u64 step_counter = 0; + + nlohmann::json json_data = nlohmann::json::array(); + + void log_state() { + nlohmann::json step_data; + step_data["step_counter"] = step_counter; + step_data["count_per_rank"] = state.count_per_rank; + step_data["msg_list"] = state.msg_list; + json_data.push_back(step_data); + } + + void dump_state() { + std::string fname = "setup_log_step.json"; + if (shamcomm::world_rank() == 0) { + logger::normal_ln("SPH setup", "dumping setup log to ", fname); + } + + std::ofstream file(fname); + file << json_data.dump(4); + file.close(); + + step_counter++; + } + + void update_count_per_rank(u64 count) { + std::vector tmp{count}; + std::vector recv_count_per_rank; + shamalgs::collective::vector_allgatherv(tmp, recv_count_per_rank, MPI_COMM_WORLD); + state.count_per_rank = recv_count_per_rank; + log_state(); + if (step_counter % 20 == 0) + dump_state(); + } + + void update_msg_list(std::vector> &msg_list) { + state.msg_list = msg_list; + log_state(); + if (step_counter % 20 == 0) + dump_state(); + } +}; + +inline constexpr f64 golden_number = 1.61803398874989484820458683436563; + +using namespace shammath; +template class shammodels::gsph::modules::GSPHSetup; +template class shammodels::gsph::modules::GSPHSetup; +template class shammodels::gsph::modules::GSPHSetup; + +template class shammodels::gsph::modules::GSPHSetup; +template class shammodels::gsph::modules::GSPHSetup; +template class shammodels::gsph::modules::GSPHSetup; diff --git a/src/shammodels/gsph/src/modules/GeneratorMCDisc.cpp b/src/shammodels/gsph/src/modules/GeneratorMCDisc.cpp new file mode 100644 index 0000000000..6f65ec6eb9 --- /dev/null +++ b/src/shammodels/gsph/src/modules/GeneratorMCDisc.cpp @@ -0,0 +1,181 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2026 Timothée David--Cléris +// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1 +// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information +// +// -------------------------------------------------------// + +/** + * @file GeneratorMCDisc.cpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) + * @brief + * + */ + +#include "shambase/constants.hpp" +#include "shamalgs/collective/indexing.hpp" +#include "shamalgs/random.hpp" +#include "shammodels/gsph/modules/setup/GeneratorMCDisc.hpp" +#include "shammodels/sph/math/density.hpp" + +template class SPHKernel> +auto shammodels::gsph::modules::GeneratorMCDisc::DiscIterator::next(u64 seed) + -> DiscOutput { + + std::mt19937_64 eng_local(seed); // ensure that 1 part = 1 random draw + + Tscal fmax = f_func(r_out); + + auto find_r = [&]() { + while (true) { + Tscal u2 = shamalgs::primitives::mock_value(eng_local, 0, fmax); + Tscal r = shamalgs::primitives::mock_value(eng_local, r_in, r_out); + if (u2 < f_func(r)) { + return r; + } + } + }; + + auto theta = shamalgs::primitives::mock_value(eng_local, 0, _2pi); + auto Gauss = shamalgs::random::mock_gaussian(eng_local); + + // depends on sigma profile + Tscal r = find_r(); + Tscal sigma = sigma_profile(r); + + // depends on H profile & sigma profile (through r) + Tscal H = H_profile(r); + Tscal z = H * Gauss; + + auto pos = sycl::vec{r * sycl::cos(theta), r * sycl::sin(theta), z}; + + // extrapolate the density from sigma profile + Tscal fs = 1; + Tscal rho = (sigma * fs) * sycl::exp(-z * z / (2 * H * H)); + + DiscOutput out{.pos = pos, .rho = rho}; + + // increase counter + check if finished + current_index++; + if (current_index == Npart) { + done = true; + } + + return out; +} + +template class SPHKernel> +bool shammodels::gsph::modules::GeneratorMCDisc::is_done() { + return generator.is_done(); +} + +template class SPHKernel> +shamrock::patch::PatchDataLayer shammodels::gsph::modules::GeneratorMCDisc::next_n( + u32 nmax) { + + using namespace shamrock::patch; + PatchScheduler &sched = shambase::get_check_ref(context.sched); + + std::vector pos_data; + + // Fill pos_data if the scheduler has some patchdata in this rank + if (!generator.is_done()) { + u64 loc_gen_count = nmax; + pos_data = generator.next_n(loc_gen_count); + } + + // extract data from disc output + std::vector vec_pos; + std::vector vec_rho; + + vec_pos.reserve(pos_data.size()); + vec_rho.reserve(pos_data.size()); + + for (DiscOutput o : pos_data) { + vec_pos.push_back(o.pos); + vec_rho.push_back(o.rho); + } + + // compute the hpart from the rho + std::vector vec_h; + vec_h.reserve(pos_data.size()); + for (Tscal rho : vec_rho) { + vec_h.push_back(shamrock::sph::h_rho(pmass, rho, Kernel::hfactd) * init_h_factor); + } + + // compute velocities + std::vector vec_vel; + vec_vel.reserve(pos_data.size()); + for (size_t i = 0; i < vec_pos.size(); i++) { + Tvec vel = vel_profile(vec_pos[i]); + vec_vel.push_back(vel); + } + + // compute the cs + bool need_cs = false; // solver_config.is_eos_locally_isothermal(); + + std::vector vec_cs; + if (need_cs) { + if (!cs_profile) { + throw shambase::make_except_with_loc( + "With this EOS you need to provide a cs_profile"); + } + vec_cs.reserve(pos_data.size()); + for (size_t i = 0; i < vec_pos.size(); i++) { + Tscal cs = cs_profile(vec_pos[i]); + vec_cs.push_back(cs); + } + } + + // Make a patchdata from pos_data + PatchDataLayer tmp(sched.get_layout_ptr_old()); + if (!pos_data.empty()) { + tmp.resize(pos_data.size()); + tmp.fields_raz(); + + { + u32 len = pos_data.size(); + PatchDataField &f + = tmp.get_field(sched.pdl_old().get_field_idx("xyz")); + sycl::buffer buf(vec_pos.data(), len); + f.override(buf, len); + } + + { + u32 len = pos_data.size(); + PatchDataField &f + = tmp.get_field(sched.pdl_old().get_field_idx("vxyz")); + sycl::buffer buf(vec_vel.data(), len); + f.override(buf, len); + } + { + u32 len = vec_pos.size(); + PatchDataField &f + = tmp.get_field(sched.pdl_old().get_field_idx("hpart")); + sycl::buffer buf(vec_h.data(), len); + f.override(buf, len); + } + + if (need_cs) { + u32 len = vec_pos.size(); + PatchDataField &f + = tmp.get_field(sched.pdl_old().get_field_idx("soundspeed")); + sycl::buffer buf(vec_cs.data(), len); + f.override(buf, len); + } + } + + return tmp; +} + +using namespace shammath; +template class shammodels::gsph::modules::GeneratorMCDisc; +template class shammodels::gsph::modules::GeneratorMCDisc; +template class shammodels::gsph::modules::GeneratorMCDisc; + +template class shammodels::gsph::modules::GeneratorMCDisc; +template class shammodels::gsph::modules::GeneratorMCDisc; +template class shammodels::gsph::modules::GeneratorMCDisc; diff --git a/src/shammodels/gsph/src/pyGSPHModel.cpp b/src/shammodels/gsph/src/pyGSPHModel.cpp index 0a18e79836..d1858ffb97 100644 --- a/src/shammodels/gsph/src/pyGSPHModel.cpp +++ b/src/shammodels/gsph/src/pyGSPHModel.cpp @@ -42,8 +42,9 @@ void add_gsph_instance(py::module &m, std::string name_config, std::string name_ using Tscal = shambase::VecComponent; - using T = Model; - using TConfig = typename T::SolverConfig; + using T = Model; + using TSPHSetup = modules::GSPHSetup; + using TConfig = typename T::SolverConfig; shamlog_debug_ln("[Py]", "registering class :", name_config, typeid(T).name()); shamlog_debug_ln("[Py]", "registering class :", name_model, typeid(T).name()); @@ -129,6 +130,23 @@ void add_gsph_instance(py::module &m, std::string name_config, std::string name_ cs : float Sound speed )==") + .def( + "set_eos_locally_isothermalFA2014", + [](TConfig &self, Tscal h_over_r) { + self.set_eos_locally_isothermalFA2014(h_over_r); + }, + py::kw_only(), + py::arg("h_over_r")) + .def( + "set_eos_locally_isothermalFA2014_extended", + [](TConfig &self, Tscal cs0, Tscal q, Tscal r0, u32 n_sinks) { + self.set_eos_locally_isothermalFA2014_extended(cs0, q, r0, n_sinks); + }, + py::kw_only(), + py::arg("cs0"), + py::arg("q"), + py::arg("r0"), + py::arg("n_sinks")) // Boundary config .def("set_boundary_free", &TConfig::set_boundary_free) .def("set_boundary_periodic", &TConfig::set_boundary_periodic) @@ -169,6 +187,212 @@ void add_gsph_instance(py::module &m, std::string name_config, std::string name_ py::arg("split_load_value"), py::arg("merge_load_value")); + std::string setup_name = name_model + "_SPHSetup"; + py::class_(m, setup_name.c_str()) + .def( + "make_generator_disc_mc", + [](TSPHSetup &self, + Tscal part_mass, + Tscal disc_mass, + Tscal r_in, + Tscal r_out, + std::function sigma_profile, + std::function H_profile, + std::function rot_profile, + std::function cs_profile, + std::function velocity_field, + std::function cs_field, + u64 random_seed, + Tscal init_h_factor) { + auto build_vel_lambda = [&]() -> std::function { + if (!velocity_field && !rot_profile) { + throw shambase::make_except_with_loc( + "make_generator_disc_mc: either velocity_field or rot_profile must be " + "provided, you must provide one of them"); + } + + if (velocity_field && rot_profile) { + throw shambase::make_except_with_loc( + "make_generator_disc_mc: either velocity_field or rot_profile must be " + "provided, you cannot provide both"); + } + + if (velocity_field) { + return std::move(velocity_field); + } + return [vth_r = std::move(rot_profile)](Tvec pos) { + pos[2] = 0; // to get the cylindrical radius + Tscal r = sycl::length(pos); + + auto etheta = sycl::vec{-pos.y(), pos.x(), 0}; + etheta /= sycl::length(etheta); + + return vth_r(r) * etheta; + }; + }; + + auto build_cs_lambda = [&]() -> std::function { + bool need_cs = false; // self.solver_config.is_eos_locally_isothermal(); + + if (!need_cs) { + if (cs_field) { + if (shamcomm::world_rank() == 0) { + logger::warn_ln( + "SPHSetup", + "make_generator_disc_mc: with the current EOS, cs_field is " + "ignored"); + } + } + if (cs_profile) { + if (shamcomm::world_rank() == 0) { + logger::warn_ln( + "SPHSetup", + "make_generator_disc_mc: with the current EOS, cs_profile is " + "ignored"); + } + } + return std::function{}; + } + + if (!cs_field && !cs_profile) { + throw shambase::make_except_with_loc( + "make_generator_disc_mc: either cs_field or cs_profile must be " + "provided, you must provide one of them"); + } + + if (cs_field && cs_profile) { + throw shambase::make_except_with_loc( + "make_generator_disc_mc: either cs_field or cs_profile must be " + "provided, you cannot provide both"); + } + + if (cs_field) { + return std::move(cs_field); + } + + return [cs_r = std::move(cs_profile)](Tvec pos) { + pos[2] = 0; // to get the cylindrical radius + Tscal r = sycl::length(pos); + return cs_r(r); + }; + }; + + return self.make_generator_disc_mc( + part_mass, + disc_mass, + r_in, + r_out, + std::move(sigma_profile), + std::move(H_profile), + build_vel_lambda(), + build_cs_lambda(), + std::mt19937_64(random_seed), + init_h_factor); + }, + py::kw_only(), + py::arg("part_mass"), + py::arg("disc_mass"), + py::arg("r_in"), + py::arg("r_out"), + py::arg("sigma_profile"), + py::arg("H_profile"), + py::arg("rot_profile") = std::function{}, + py::arg("cs_profile") = std::function{}, + py::arg("velocity_field") = std::function{}, + py::arg("cs_field") = std::function{}, + py::arg("random_seed"), + py::arg("init_h_factor") = 0.8, + R"pbdoc( + Create a Monte Carlo disc particle generator. + + Particles are sampled in cylindrical coordinates: the radius is drawn + with rejection sampling from ``sigma_profile``, the azimuth is uniform, + and the vertical coordinate follows a Gaussian with scale ``H_profile(r)``. + The initial density is extrapolated from the surface density profile, and + smoothing lengths are set from that density. + + Args: + part_mass: Mass of each SPH particle. + disc_mass: Total disc mass. The particle count is ``disc_mass / part_mass``. + r_in: Inner disc radius. + r_out: Outer disc radius. + sigma_profile: Surface density profile ``sigma(r)``. + H_profile: Disc scale height profile ``H(r)``. + rot_profile: Azimuthal speed profile ``v_theta(r)``. The velocity is + projected along the cylindrical azimuthal direction at each + particle position. Mutually exclusive with ``velocity_field``. + cs_profile: Sound speed profile ``c_s(r)``. Evaluated at the cylindrical + radius of each particle. Required when the solver uses a locally + isothermal EOS. Mutually exclusive with ``cs_field``. + velocity_field: Velocity profile ``v(x, y, z)``. Mutually exclusive + with ``rot_profile``. + cs_field: Sound speed profile ``c_s(x, y, z)``. Required when the solver + uses a locally isothermal EOS. Mutually exclusive with ``cs_profile``. + random_seed: Seed for the Monte Carlo sampler. + init_h_factor: Multiplier applied to the smoothing length inferred from + the generated density. Defaults to ``0.8``. + + Notes: + Exactly one of ``velocity_field`` or ``rot_profile`` must be provided. + + If the solver uses a locally isothermal EOS, exactly one of ``cs_field`` + or ``cs_profile`` must be provided. Otherwise both sound-speed profiles + are ignored and a warning is emitted if either is supplied. + + Returns: + A setup node to pass to :py:meth:`apply_setup`. + )pbdoc") + .def( + "apply_setup", + [](TSPHSetup &self, + modules::SetupNodePtr setup, + bool part_reordering, + std::optional gen_step, + std::optional insert_step, + std::optional msg_count_limit, + std::optional msg_size_limit, + std::optional max_msg_size, + bool do_setup_log, + bool use_new_setup, + bool speculative_balancing) { + if (bool(gen_step)) { + ON_RANK_0( + logger::warn_ln("SPHSetup", "gen_step is ignored when using old setup")); + } + if (bool(msg_count_limit)) { + ON_RANK_0( + logger::warn_ln( + "SPHSetup", "msg_count_limit is ignored when using old setup")); + } + if (bool(msg_size_limit)) { + ON_RANK_0( + logger::warn_ln( + "SPHSetup", "msg_size_limit is ignored when using old setup")); + } + if (bool(max_msg_size)) { + ON_RANK_0( + logger::warn_ln( + "SPHSetup", "max_msg_size is ignored when using old setup")); + } + if (bool(do_setup_log)) { + ON_RANK_0( + logger::warn_ln( + "SPHSetup", "do_setup_log is ignored when using old setup")); + } + return self.apply_setup(setup, part_reordering, insert_step); + }, + py::arg("setup"), + py::kw_only(), + py::arg("part_reordering") = true, + py::arg("gen_step") = std::nullopt, + py::arg("insert_step") = std::nullopt, + py::arg("msg_count_limit") = std::nullopt, + py::arg("rank_comm_size_limit") = std::nullopt, + py::arg("max_msg_size") = std::nullopt, + py::arg("do_setup_log") = false, + py::arg("use_new_setup") = true, + py::arg("speculative_balancing") = false); + py::class_(m, name_model.c_str()) .def(py::init([](ShamrockCtx &ctx) { return std::make_unique(ctx); @@ -352,6 +576,28 @@ void add_gsph_instance(py::module &m, std::string name_config, std::string name_ return self.solver.solver_config; }) .def("set_solver_config", &T::set_solver_config) + .def("add_sink", &T::add_sink) + .def( + "get_sinks", + [](T &self) { + py::list list_out; + + if (!self.solver.storage.sinks.is_empty()) { + for (auto &sink : self.solver.storage.sinks.get()) { + py::dict sink_dic; + sink_dic["pos"] = sink.pos; + sink_dic["velocity"] = sink.velocity; + sink_dic["sph_acceleration"] = sink.sph_acceleration; + sink_dic["ext_acceleration"] = sink.ext_acceleration; + sink_dic["mass"] = sink.mass; + sink_dic["angular_momentum"] = sink.angular_momentum; + sink_dic["accretion_radius"] = sink.accretion_radius; + list_out.append(sink_dic); + } + } + + return list_out; + }) .def("do_vtk_dump", &T::do_vtk_dump) .def("solver_logs_last_rate", &T::solver_logs_last_rate) .def("solver_logs_last_obj_count", &T::solver_logs_last_obj_count) @@ -410,7 +656,8 @@ void add_gsph_instance(py::module &m, std::string name_config, std::string name_ Example ------- >>> model.dump("checkpoint.shamrock") -)=="); +)==") + .def("get_setup", &T::get_setup); } using namespace shammodels::gsph; @@ -422,6 +669,13 @@ ON_PYTHON_INIT { using namespace shammodels::gsph; + py::class_< + shammodels::gsph::modules::IGSPHSetupNode, + std::shared_ptr>(mgsph, "IGSPHSetupNode") + .def("get_dot", [](std::shared_ptr &self) { + return self->get_dot(); + }); + // Register GSPH models for different kernels add_gsph_instance( mgsph, "GSPHModel_f64_3_M4_SolverConfig", "GSPHModel_f64_3_M4"); From a079ee15060ba38d64ef572ba869ed66bba493e8 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Fri, 31 Jul 2026 16:28:55 +0900 Subject: [PATCH 2/7] sink update ok for single sink --- src/shammodels/gsph/CMakeLists.txt | 2 + .../include/shammodels/gsph/SolverConfig.hpp | 10 + .../gsph/modules/ExternalForces.hpp | 67 +++ .../gsph/modules/SinkParticlesUpdate.hpp | 57 +++ src/shammodels/gsph/src/Solver.cpp | 23 +- src/shammodels/gsph/src/SolverConfig.cpp | 2 + .../gsph/src/modules/ExternalForces.cpp | 450 ++++++++++++++++++ .../gsph/src/modules/SinkParticlesUpdate.cpp | 423 ++++++++++++++++ 8 files changed, 1032 insertions(+), 2 deletions(-) create mode 100644 src/shammodels/gsph/include/shammodels/gsph/modules/ExternalForces.hpp create mode 100644 src/shammodels/gsph/include/shammodels/gsph/modules/SinkParticlesUpdate.hpp create mode 100644 src/shammodels/gsph/src/modules/ExternalForces.cpp create mode 100644 src/shammodels/gsph/src/modules/SinkParticlesUpdate.cpp diff --git a/src/shammodels/gsph/CMakeLists.txt b/src/shammodels/gsph/CMakeLists.txt index 0a0960b89b..a8c2f0b79b 100644 --- a/src/shammodels/gsph/CMakeLists.txt +++ b/src/shammodels/gsph/CMakeLists.txt @@ -21,6 +21,8 @@ set(Sources src/modules/GSPHSetup.cpp src/modules/GeneratorMCDisc.cpp src/modules/ComputeLoadBalanceValue.cpp + src/modules/SinkParticlesUpdate.cpp + src/modules/ExternalForces.cpp ) if(SHAMROCK_USE_SHARED_LIB) diff --git a/src/shammodels/gsph/include/shammodels/gsph/SolverConfig.hpp b/src/shammodels/gsph/include/shammodels/gsph/SolverConfig.hpp index aa3238fe07..283c673c08 100644 --- a/src/shammodels/gsph/include/shammodels/gsph/SolverConfig.hpp +++ b/src/shammodels/gsph/include/shammodels/gsph/SolverConfig.hpp @@ -112,6 +112,16 @@ struct shammodels::gsph::SolverConfig { } } + inline Tscal get_constant_c() const { + if (!unit_sys) { + ON_RANK_0(logger::warn_ln("gsph::Config", "the unit system is not set")); + shamunits::Constants ctes{shamunits::UnitSystem{}}; + return ctes.c(); + } else { + return shamunits::Constants{*unit_sys}.c(); + } + } + ////////////////////////////////////////////////////////////////////////////////////////////// // Units Config (END) ////////////////////////////////////////////////////////////////////////////////////////////// diff --git a/src/shammodels/gsph/include/shammodels/gsph/modules/ExternalForces.hpp b/src/shammodels/gsph/include/shammodels/gsph/modules/ExternalForces.hpp new file mode 100644 index 0000000000..397d16258d --- /dev/null +++ b/src/shammodels/gsph/include/shammodels/gsph/modules/ExternalForces.hpp @@ -0,0 +1,67 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2026 Timothée David--Cléris +// 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 + +/** + * @file ExternalForces.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) + * @brief + * + */ + +#include "shambackends/typeAliasVec.hpp" +#include "shambackends/vec.hpp" +#include "shammodels/gsph/SolverConfig.hpp" +#include "shammodels/gsph/modules/SolverStorage.hpp" +#include "shamrock/scheduler/ShamrockCtx.hpp" + +namespace shammodels::gsph::modules { + + template class SPHKernel> + class ExternalForces { + public: + using Tscal = shambase::VecComponent; + static constexpr u32 dim = shambase::VectorProperties::dimension; + using Kernel = SPHKernel; + + using Config = SolverConfig; + using Storage = SolverStorage; + + ShamrockCtx &context; + Config &solver_config; + Storage &storage; + + ExternalForces(ShamrockCtx &context, Config &solver_config, Storage &storage) + : context(context), solver_config(solver_config), storage(storage) {} + + /** + * @brief is ran once per timestep, it computes the forces that are independant of velocity + * + */ + void compute_ext_forces_indep_v(); + + /** + * @brief add external forces to the particle acceleration, note that forces dependant on + * velocity shlould be added here + * + */ + void add_ext_forces(); + + void point_mass_accrete_particles(); + + private: + using SolverConfigExtForce = typename Config::ExtForceConfig; + using EF_PointMass = typename SolverConfigExtForce::PointMass; + + inline PatchScheduler &scheduler() { return shambase::get_check_ref(context.sched); } + }; + +} // namespace shammodels::gsph::modules diff --git a/src/shammodels/gsph/include/shammodels/gsph/modules/SinkParticlesUpdate.hpp b/src/shammodels/gsph/include/shammodels/gsph/modules/SinkParticlesUpdate.hpp new file mode 100644 index 0000000000..6a99c0b08a --- /dev/null +++ b/src/shammodels/gsph/include/shammodels/gsph/modules/SinkParticlesUpdate.hpp @@ -0,0 +1,57 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2026 Timothée David--Cléris +// 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 + +/** + * @file SinkParticlesUpdate.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief + * + */ + +#include "shambackends/typeAliasVec.hpp" +#include "shambackends/vec.hpp" +#include "shammodels/gsph/SolverConfig.hpp" +#include "shammodels/gsph/modules/SolverStorage.hpp" +#include "shammodels/sph/SinkPartStruct.hpp" +#include "shamrock/scheduler/ShamrockCtx.hpp" + +namespace shammodels::gsph::modules { + + template class SPHKernel> + class SinkParticlesUpdate { + public: + using Tscal = shambase::VecComponent; + static constexpr u32 dim = shambase::VectorProperties::dimension; + using Kernel = SPHKernel; + + using Config = SolverConfig; + using Storage = SolverStorage; + + ShamrockCtx &context; + Config &solver_config; + Storage &storage; + + using Sink = sph::SinkParticle; + + SinkParticlesUpdate(ShamrockCtx &context, Config &solver_config, Storage &storage) + : context(context), solver_config(solver_config), storage(storage) {} + + void accrete_particles(Tscal dt); + void predictor_step(Tscal dt); + void compute_sph_forces(); + void compute_ext_forces(); + void corrector_step(Tscal dt); + + private: + inline PatchScheduler &scheduler() { return shambase::get_check_ref(context.sched); } + }; + +} // namespace shammodels::gsph::modules diff --git a/src/shammodels/gsph/src/Solver.cpp b/src/shammodels/gsph/src/Solver.cpp index a713955538..a985d18986 100644 --- a/src/shammodels/gsph/src/Solver.cpp +++ b/src/shammodels/gsph/src/Solver.cpp @@ -11,7 +11,7 @@ * @file Solver.cpp * @author Guo Yansong (guo.yansong.ngy@gmail.com) * @author Timothée David--Cléris (tim.shamrock@proton.me) - * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) --no git blame-- + * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) --no gi´t blame-- * @brief GSPH Solver implementation * * The GSPH method originated from: @@ -38,7 +38,9 @@ #include "shammodels/gsph/Solver.hpp" #include "shammodels/gsph/SolverConfig.hpp" #include "shammodels/gsph/config/FieldNames.hpp" +#include "shammodels/gsph/modules/ExternalForces.hpp" #include "shammodels/gsph/modules/GSPHUtilities.hpp" +#include "shammodels/gsph/modules/SinkParticlesUpdate.hpp" #include "shammodels/gsph/modules/UpdateDerivs.hpp" #include "shammodels/gsph/modules/io/VTKDump.hpp" #include "shammodels/sph/modules/IterateSmoothingLengthDensity.hpp" @@ -1534,7 +1536,11 @@ template class Kern> void shammodels::gsph::Solver::update_derivs() { StackEntry stack_loc{}; // GSPH derivative update using Riemann solver - gsph::modules::UpdateDerivs(context, solver_config, storage).update_derivs(); + gsph::modules::UpdateDerivs derivs(context, solver_config, storage); + derivs.update_derivs(); + + modules::ExternalForces ext_forces(context, solver_config, storage); + ext_forces.add_ext_forces(); } template class Kern> @@ -1760,10 +1766,22 @@ shammodels::gsph::TimestepLog shammodels::gsph::Solver::evolve_once( // 7. CFL: compute next timestep // ========================================================================= + modules::SinkParticlesUpdate sink_update(context, solver_config, storage); + modules::ExternalForces ext_forces(context, solver_config, storage); + + // STEP 0: SINK PARTICLES + sink_update.accrete_particles(dt); + ext_forces.point_mass_accrete_particles(); + + sink_update.predictor_step(dt); + // STEP 1: PREDICTOR - move particles using OLD accelerations // (On first iteration, accelerations are zero, so this is just position drift) do_predictor_leapfrog(dt); + sink_update.compute_ext_forces(); + ext_forces.compute_ext_forces_indep_v(); + // STEP 2: BOUNDARY - apply boundary conditions to NEW positions // Build serial patch tree first (needed for boundary application) gen_serial_patch_tree(); @@ -1819,6 +1837,7 @@ shammodels::gsph::TimestepLog shammodels::gsph::Solver::evolve_once( // STEP 6: CORRECTOR - refine velocities apply_corrector(dt, Npart_all); + sink_update.corrector_step(dt); // STEP 7: CFL - compute next timestep Tscal dt_next = compute_dt_cfl(); diff --git a/src/shammodels/gsph/src/SolverConfig.cpp b/src/shammodels/gsph/src/SolverConfig.cpp index 2bd44afa95..bc1ccead3a 100644 --- a/src/shammodels/gsph/src/SolverConfig.cpp +++ b/src/shammodels/gsph/src/SolverConfig.cpp @@ -35,6 +35,8 @@ void shammodels::gsph::SolverConfig::set_layout( // Smoothing length pdl.add_field(names::common::hpart, 1); + pdl.add_field("axyz_ext", 1); + // Internal energy (for adiabatic EOS) if (has_field_uint()) { pdl.add_field(names::newtonian::uint, 1); diff --git a/src/shammodels/gsph/src/modules/ExternalForces.cpp b/src/shammodels/gsph/src/modules/ExternalForces.cpp new file mode 100644 index 0000000000..32c59ee54c --- /dev/null +++ b/src/shammodels/gsph/src/modules/ExternalForces.cpp @@ -0,0 +1,450 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2026 Timothée David--Cléris +// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1 +// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information +// +// -------------------------------------------------------// + +/** + * @file ExternalForces.cpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) + * @brief + * + */ + +#include "shambase/memory.hpp" +#include "shambackends/kernel_call.hpp" +#include "shambackends/kernel_call_distrib.hpp" +#include "shamcomm/logs.hpp" +#include "shammath/sphkernels.hpp" +#include "shammodels/common/modules/AddForceCentralGravPotential.hpp" +#include "shammodels/gsph/modules/ExternalForces.hpp" +#include "shammodels/gsph/modules/SinkParticlesUpdate.hpp" +#include "shamrock/solvergraph/IDataEdge.hpp" +#include "shamrock/solvergraph/INode.hpp" +#include "shamrock/solvergraph/NodeSetEdge.hpp" +#include "shamrock/solvergraph/OperationSequence.hpp" +#include "shamrock/solvergraph/SolverGraph.hpp" +#include "shamsys/legacy/log.hpp" +#include "shamunits/Constants.hpp" + +namespace shambase { + + template + std::shared_ptr to_shared(T &&t) { + return std::make_shared(std::forward(t)); + } +} // namespace shambase + +template class SPHKernel> +void shammodels::gsph::modules::ExternalForces::compute_ext_forces_indep_v() { + + StackEntry stack_loc{}; + + sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue(); + + Tscal gpart_mass = solver_config.gpart_mass; + + using namespace shamrock; + using namespace shamrock::patch; + + PatchDataLayerLayout &pdl = scheduler().pdl_old(); + + const u32 iaxyz_ext = pdl.get_field_idx("axyz_ext"); + modules::SinkParticlesUpdate sink_update(context, solver_config, storage); + + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + PatchDataField &field = pdat.get_field(iaxyz_ext); + field.field_raz(); + }); + + sink_update.compute_sph_forces(); + + if (solver_config.ext_force_config.ext_forces.empty()) { + return; + } + + auto field_xyz = shamrock::solvergraph::FieldRefs::make_shared("", ""); + + shamrock::solvergraph::NodeSetEdge> set_field_xyz( + [&](shamrock::solvergraph::FieldRefs &field_xyz_edge) { + shamrock::solvergraph::DDPatchDataFieldRef field_xyz_refs = {}; + scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) { + auto &field = pdat.get_field(0); + field_xyz_refs.add_obj(p.id_patch, std::ref(field)); + }); + field_xyz_edge.set_refs(field_xyz_refs); + }); + set_field_xyz.set_edges(field_xyz); + set_field_xyz.evaluate(); + + auto field_axyz_ext = shamrock::solvergraph::FieldRefs::make_shared("", ""); + + shamrock::solvergraph::NodeSetEdge> set_field_axyz_ext( + [&](shamrock::solvergraph::FieldRefs &field_axyz_ext_edge) { + shamrock::solvergraph::DDPatchDataFieldRef field_axyz_ext_refs = {}; + scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) { + auto &field = pdat.get_field(iaxyz_ext); + field_axyz_ext_refs.add_obj(p.id_patch, std::ref(field)); + }); + field_axyz_ext_edge.set_refs(field_axyz_ext_refs); + }); + set_field_axyz_ext.set_edges(field_axyz_ext); + set_field_axyz_ext.evaluate(); + + auto sizes = shamrock::solvergraph::Indexes::make_shared("", ""); + + shamrock::solvergraph::NodeSetEdge> set_sizes( + [&](shamrock::solvergraph::Indexes &sizes) { + sizes.indexes = {}; + scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) { + sizes.indexes.add_obj(p.id_patch, pdat.get_obj_cnt()); + }); + }); + set_sizes.set_edges(sizes); + set_sizes.evaluate(); + + auto constant_G = shamrock::solvergraph::IDataEdge::make_shared("", ""); + auto constant_c = shamrock::solvergraph::IDataEdge::make_shared("", ""); + + shamrock::solvergraph::NodeSetEdge> set_constant_G( + [&](shamrock::solvergraph::IDataEdge &constant_G) { + constant_G.data = solver_config.get_constant_G(); + }); + + shamrock::solvergraph::NodeSetEdge> set_constant_c( + [&](shamrock::solvergraph::IDataEdge &constant_c) { + constant_c.data = solver_config.get_constant_c(); + }); + + set_constant_G.set_edges(constant_G); + set_constant_c.set_edges(constant_c); + + std::vector> add_ext_forces_seq{}; + + for (auto var_force : solver_config.ext_force_config.ext_forces) { + if (EF_PointMass *ext_force = std::get_if(&var_force.val)) { + + auto central_mass = shamrock::solvergraph::IDataEdge::make_shared("", ""); + auto central_pos = shamrock::solvergraph::IDataEdge::make_shared("", ""); + + shamrock::solvergraph::NodeSetEdge> + set_central_mass([cmass = ext_force->central_mass]( + shamrock::solvergraph::IDataEdge ¢ral_mass) { + central_mass.data = cmass; + }); + set_central_mass.set_edges(central_mass); + + shamrock::solvergraph::NodeSetEdge> + set_central_pos([&](shamrock::solvergraph::IDataEdge ¢ral_pos) { + central_pos.data = {}; // no support for offset yet + }); + set_central_pos.set_edges(central_pos); + + common::modules::AddForceCentralGravPotential add_force_central_grav_potential; + add_force_central_grav_potential.set_edges( + constant_G, central_mass, central_pos, field_xyz, sizes, field_axyz_ext); + + add_ext_forces_seq.push_back( + std::make_shared( + "Point mass", + std::vector>{ + shambase::to_shared(std::move(set_central_pos)), + shambase::to_shared(std::move(set_central_mass)), + shambase::to_shared(std::move(add_force_central_grav_potential))})); + + } else { + shambase::throw_unimplemented("this force is not handled, yet ..."); + } + } + + set_constant_G.evaluate(); + set_constant_c.evaluate(); + + if (add_ext_forces_seq.size() > 0) { + shamrock::solvergraph::OperationSequence seq( + "Add external forces", std::move(add_ext_forces_seq)); + seq.evaluate(); + } +} + +template +std::shared_ptr register_constant_set( + shamrock::solvergraph::SolverGraph &solver_graph, std::string name, std::function getter) { + solver_graph.register_edge(name, shamrock::solvergraph::IDataEdge("", "")); + + solver_graph.register_node( + "set_" + name, + shamrock::solvergraph::NodeSetEdge>( + [getter](shamrock::solvergraph::IDataEdge &edge) { + edge.data = getter(); + })); + + solver_graph + .get_node_ref>>( + "set_" + name) + .set_edges(solver_graph.get_edge_ptr_base(name)); + + return solver_graph.get_node_ptr_base("set_" + name); +} + +template class SPHKernel> +void shammodels::gsph::modules::ExternalForces::add_ext_forces() { + + StackEntry stack_loc{}; + + sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue(); + + Tscal gpart_mass = solver_config.gpart_mass; + + using namespace shamrock; + using namespace shamrock::patch; + + PatchDataLayerLayout &pdl = scheduler().pdl_old(); + + const u32 iaxyz = pdl.get_field_idx("axyz"); + const u32 ivxyz = pdl.get_field_idx("vxyz"); + const u32 iaxyz_ext = pdl.get_field_idx("axyz_ext"); + + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + sham::DeviceBuffer &buf_axyz = pdat.get_field_buf_ref(iaxyz); + sham::DeviceBuffer &buf_axyz_ext = pdat.get_field_buf_ref(iaxyz_ext); + + sham::EventList depends_list; + auto axyz = buf_axyz.get_write_access(depends_list); + auto axyz_ext = buf_axyz_ext.get_read_access(depends_list); + + auto e = q.submit(depends_list, [&](sycl::handler &cgh) { + shambase::parallel_for( + cgh, pdat.get_obj_cnt(), "add ext force acc to acc", [=](u64 gid) { + axyz[gid] += axyz_ext[gid]; + }); + }); + + buf_axyz.complete_event_state(e); + buf_axyz_ext.complete_event_state(e); + }); + + if (solver_config.ext_force_config.ext_forces.empty()) { + return; // skip if no external forces + } + + using SolverConfigExtForce = typename Config::ExtForceConfig; + using EF_PointMass = typename SolverConfigExtForce::PointMass; + + using namespace shamrock::solvergraph; + SolverGraph solver_graph{}; + + auto set_constant_G = register_constant_set(solver_graph, "constant_G", [&]() { + return solver_config.get_constant_G(); + }); + auto set_constant_c = register_constant_set(solver_graph, "constant_c", [&]() { + return solver_config.get_constant_c(); + }); + + bool is_G_needed = false; + bool is_c_needed = false; + + for (auto var_force : solver_config.ext_force_config.ext_forces) { + if (EF_PointMass *ext_force = std::get_if(&var_force.val)) { + } else { + shambase::throw_unimplemented("this force is not handled, yet ..."); + } + } + + std::vector> add_ext_forces_seq{}; + + if (is_G_needed) { + add_ext_forces_seq.push_back(set_constant_G); + } + if (is_c_needed) { + add_ext_forces_seq.push_back(set_constant_c); + } + + auto field_xyz = solver_graph.register_edge("field_xyz", FieldRefs("", "")); + auto field_vxyz = solver_graph.register_edge("field_vxyz", FieldRefs("", "")); + auto field_axyz = solver_graph.register_edge("field_axyz", FieldRefs("", "")); + auto field_sizes = solver_graph.register_edge("field_sizes", Indexes("", "")); + + auto set_field_xyz = solver_graph.register_node( + "set_field_xyz", NodeSetEdge>([&](FieldRefs &field_xyz_edge) { + DDPatchDataFieldRef field_xyz_refs = {}; + scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) { + auto &field = pdat.get_field(0); + field_xyz_refs.add_obj(p.id_patch, std::ref(field)); + }); + field_xyz_edge.set_refs(field_xyz_refs); + })); + shambase::get_check_ref(set_field_xyz).set_edges(field_xyz); + + auto set_field_vxyz = solver_graph.register_node( + "set_field_vxyz", NodeSetEdge>([&](FieldRefs &field_vxyz_edge) { + DDPatchDataFieldRef field_vxyz_refs = {}; + scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) { + auto &field = pdat.get_field(ivxyz); + field_vxyz_refs.add_obj(p.id_patch, std::ref(field)); + }); + field_vxyz_edge.set_refs(field_vxyz_refs); + })); + shambase::get_check_ref(set_field_vxyz).set_edges(field_vxyz); + + auto set_field_axyz = solver_graph.register_node( + "set_field_axyz", NodeSetEdge>([&](FieldRefs &field_axyz_edge) { + DDPatchDataFieldRef field_axyz_refs = {}; + scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) { + auto &field = pdat.get_field(iaxyz); + field_axyz_refs.add_obj(p.id_patch, std::ref(field)); + }); + field_axyz_edge.set_refs(field_axyz_refs); + })); + shambase::get_check_ref(set_field_axyz).set_edges(field_axyz); + + auto set_field_sizes = solver_graph.register_node( + "set_field_sizes", NodeSetEdge>([&](Indexes &sizes) { + sizes.indexes = {}; + scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) { + sizes.indexes.add_obj(p.id_patch, pdat.get_obj_cnt()); + }); + })); + shambase::get_check_ref(set_field_sizes).set_edges(field_sizes); + + add_ext_forces_seq.push_back(set_field_xyz); + add_ext_forces_seq.push_back(set_field_vxyz); + add_ext_forces_seq.push_back(set_field_axyz); + add_ext_forces_seq.push_back(set_field_sizes); + + for (u32 i = 0; i < solver_config.ext_force_config.ext_forces.size(); i++) { + + auto &var_force = solver_config.ext_force_config.ext_forces[i]; + + std::string prefix = shambase::format("ext_force_{}_", i); + + if (EF_PointMass *ext_force = std::get_if(&var_force.val)) { + + } else { + shambase::throw_unimplemented("this force is not handled, yet ..."); + } + } + + if (add_ext_forces_seq.size() > 0) { + OperationSequence seq("Add external forces", std::move(add_ext_forces_seq)); + seq.evaluate(); + } +} + +template class SPHKernel> +void shammodels::gsph::modules::ExternalForces::point_mass_accrete_particles() { + + StackEntry stack_loc{}; + + Tscal gpart_mass = solver_config.gpart_mass; + + using namespace shamrock; + using namespace shamrock::patch; + + using SolverConfigExtForce = typename Config::ExtForceConfig; + using EF_PointMass = typename SolverConfigExtForce::PointMass; + + PatchDataLayerLayout &pdl = scheduler().pdl_old(); + const u32 ixyz = pdl.get_field_idx("xyz"); + const u32 ivxyz = pdl.get_field_idx("vxyz"); + + auto dev_sched = shamsys::instance::get_compute_scheduler_ptr(); + + sham::DeviceQueue &q = shambase::get_check_ref(dev_sched).get_queue(); + + for (auto var_force : solver_config.ext_force_config.ext_forces) { + + Tvec pos_accretion; + Tscal Racc; + + if (EF_PointMass *ext_force = std::get_if(&var_force.val)) { + pos_accretion = {0, 0, 0}; + Racc = ext_force->Racc; + } else { + continue; + } + + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + u32 Nobj = pdat.get_obj_cnt(); + + sham::DeviceBuffer &buf_xyz = pdat.get_field_buf_ref(ixyz); + sham::DeviceBuffer &buf_vxyz = pdat.get_field_buf_ref(ivxyz); + + sycl::buffer not_accreted(Nobj); + sycl::buffer accreted(Nobj); + + sham::EventList depends_list; + auto xyz = buf_xyz.get_read_access(depends_list); + + auto e = q.submit(depends_list, [&](sycl::handler &cgh) { + sycl::accessor not_acc{not_accreted, cgh, sycl::write_only, sycl::no_init}; + sycl::accessor acc{accreted, cgh, sycl::write_only, sycl::no_init}; + + Tvec r_sink = pos_accretion; + Tscal acc_rad2 = Racc * Racc; + + shambase::parallel_for(cgh, Nobj, "check accretion", [=](i32 id_a) { + Tvec r = xyz[id_a] - r_sink; + bool not_accreted = sycl::dot(r, r) > acc_rad2; + not_acc[id_a] = (not_accreted) ? 1 : 0; + acc[id_a] = (!not_accreted) ? 1 : 0; + }); + }); + + buf_xyz.complete_event_state(e); + + std::tuple>, u32> id_list_keep + = shamalgs::numeric::stream_compact(q.q, not_accreted, Nobj); + + std::tuple>, u32> id_list_accrete + = shamalgs::numeric::stream_compact(q.q, accreted, Nobj); + + // sum accreted values onto sink + + if (std::get<1>(id_list_accrete) > 0) { + + u32 Naccrete = std::get<1>(id_list_accrete); + + Tscal acc_mass = gpart_mass * Naccrete; + + sham::DeviceBuffer pxyz_acc(Naccrete, dev_sched); + + sham::EventList depends_list; + + auto vxyz = buf_vxyz.get_read_access(depends_list); + auto accretion_p = pxyz_acc.get_write_access(depends_list); + + auto e = q.submit(depends_list, [&, gpart_mass](sycl::handler &cgh) { + sycl::accessor id_acc{*std::get<0>(id_list_accrete), cgh, sycl::read_only}; + + shambase::parallel_for( + cgh, Naccrete, "compute sum momentum accretion", [=](i32 id_a) { + accretion_p[id_a] = gpart_mass * vxyz[id_acc[id_a]]; + }); + }); + + buf_vxyz.complete_event_state(e); + pxyz_acc.complete_event_state(e); + + Tvec acc_pxyz = shamalgs::primitives::sum(dev_sched, pxyz_acc, 0, Naccrete); + + logger::raw_ln("central potential accretion : += ", acc_mass); + + pdat.keep_ids(*std::get<0>(id_list_keep), std::get<1>(id_list_keep)); + } + }); + } +} + +using namespace shammath; +template class shammodels::gsph::modules::ExternalForces; +template class shammodels::gsph::modules::ExternalForces; +template class shammodels::gsph::modules::ExternalForces; + +template class shammodels::gsph::modules::ExternalForces; +template class shammodels::gsph::modules::ExternalForces; +template class shammodels::gsph::modules::ExternalForces; diff --git a/src/shammodels/gsph/src/modules/SinkParticlesUpdate.cpp b/src/shammodels/gsph/src/modules/SinkParticlesUpdate.cpp new file mode 100644 index 0000000000..70192f10e8 --- /dev/null +++ b/src/shammodels/gsph/src/modules/SinkParticlesUpdate.cpp @@ -0,0 +1,423 @@ +// -------------------------------------------------------// +// +// SHAMROCK code for hydrodynamics +// Copyright (c) 2021-2026 Timothée David--Cléris +// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1 +// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information +// +// -------------------------------------------------------// + +/** + * @file SinkParticlesUpdate.cpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) + * @brief + * + */ + +#include "shambase/DistributedData.hpp" +#include "shambase/narrowing.hpp" +#include "shamalgs/collective/reduction.hpp" +#include "shamalgs/details/numeric/numeric.hpp" +#include "shamalgs/primitives/reduction.hpp" +#include "shambackends/DeviceBuffer.hpp" +#include "shambackends/kernel_call.hpp" +#include "shamcomm/logs.hpp" +#include "shammath/sphkernels.hpp" +#include "shammodels/gsph/modules/SinkParticlesUpdate.hpp" +#include + +template class SPHKernel> +void shammodels::gsph::modules::SinkParticlesUpdate::accrete_particles(Tscal dt) { + StackEntry stack_loc{}; + + Tscal gpart_mass = solver_config.gpart_mass; + + if (storage.sinks.is_empty()) { + return; + } + + using namespace shamrock; + using namespace shamrock::patch; + + PatchDataLayerLayout &pdl = scheduler().pdl_old(); + const u32 ixyz = pdl.get_field_idx("xyz"); + const u32 ivxyz = pdl.get_field_idx("vxyz"); + const u32 iaxyz = pdl.get_field_idx("axyz"); + + auto dev_sched = shamsys::instance::get_compute_scheduler_ptr(); + sham::DeviceQueue &q = shambase::get_check_ref(dev_sched).get_queue(); + + std::vector &sink_parts = storage.sinks.get(); + + u32 sink_id = 0; + bool had_accretion = false; + std::string log = "sink accretion :"; + + struct AccretionFlagBufs { + sham::DeviceBuffer not_accreted; + sham::DeviceBuffer accreted; + }; + + for (size_t sink_id = 0; sink_id < sink_parts.size(); sink_id++) { + Sink &s = sink_parts[sink_id]; + + Tvec r_sink = s.pos; + Tvec v_sink = s.velocity; + Tscal acc_rad2 = s.accretion_radius * s.accretion_radius; + + // flags particles for accretion + shambase::DistributedData accretion_flag_bufs{}; + + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + u32 Nobj = pdat.get_obj_cnt(); + + sham::DeviceBuffer &buf_xyz = pdat.get_field_buf_ref(ixyz); + sham::DeviceBuffer &buf_vxyz = pdat.get_field_buf_ref(ivxyz); + + sham::DeviceBuffer not_accreted(Nobj, dev_sched); + sham::DeviceBuffer accreted(Nobj, dev_sched); + + sham::kernel_call( + q, + sham::MultiRef{buf_xyz}, + sham::MultiRef{not_accreted, accreted}, + Nobj, + [r_sink, acc_rad2]( + u32 id_a, + const Tvec *__restrict xyz, + u32 *__restrict not_acc, + u32 *__restrict acc) { + Tvec r = xyz[id_a] - r_sink; + bool not_accreted = sycl::dot(r, r) > acc_rad2; + not_acc[id_a] = (not_accreted) ? 1 : 0; + acc[id_a] = (!not_accreted) ? 1 : 0; + }); + + accretion_flag_bufs.add_obj( + cur_p.id_patch, AccretionFlagBufs{std::move(not_accreted), std::move(accreted)}); + }); + + // list the ids that will be accreted + shambase::DistributedData> bufs_id_list_accrete{}; + + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + u32 Nobj = pdat.get_obj_cnt(); + + sham::DeviceBuffer &accreted = accretion_flag_bufs.get(cur_p.id_patch).accreted; + + sham::DeviceBuffer id_list_accrete + = shamalgs::stream_compact(dev_sched, accreted, Nobj); + + bufs_id_list_accrete.add_obj(cur_p.id_patch, std::move(id_list_accrete)); + }); + + // compute the accreted mass, position moment and linear momentum + Tscal s_acc_mass = 0; + Tvec s_acc_mxyz = {0, 0, 0}; + Tvec s_acc_pxyz = {0, 0, 0}; + Tvec s_acc_maxyz = {0, 0, 0}; + Tvec s_acc_lxyz = {0, 0, 0}; + + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + u32 Nobj = pdat.get_obj_cnt(); + + sham::DeviceBuffer &buf_xyz = pdat.get_field_buf_ref(ixyz); + sham::DeviceBuffer &buf_vxyz = pdat.get_field_buf_ref(ivxyz); + sham::DeviceBuffer &buf_axyz = pdat.get_field_buf_ref(iaxyz); + + sham::DeviceBuffer &id_list_accrete = bufs_id_list_accrete.get(cur_p.id_patch); + + // sum accreted values onto sink + if (id_list_accrete.get_size() > 0) { + u32 Naccrete = shambase::narrow_or_throw(id_list_accrete.get_size()); + + Tscal acc_mass = gpart_mass * Naccrete; + + sham::DeviceBuffer pxyz_acc(Naccrete, dev_sched); + sham::DeviceBuffer maxyz_acc(Naccrete, dev_sched); + sham::DeviceBuffer mxyz_acc(Naccrete, dev_sched); + sham::DeviceBuffer lxyz_acc(Naccrete, dev_sched); + + sham::kernel_call( + q, + sham::MultiRef{buf_xyz, buf_vxyz, buf_axyz, id_list_accrete}, + sham::MultiRef{pxyz_acc, mxyz_acc, maxyz_acc, lxyz_acc}, + Naccrete, + [gpart_mass, r_sink, v_sink, dt]( + u32 id_a, + const Tvec *__restrict xyz, + const Tvec *__restrict vxyz, + const Tvec *__restrict axyz, + const u32 *__restrict id_acc, + Tvec *__restrict accretion_p, + Tvec *__restrict accretion_mr, + Tvec *__restrict accretion_ma, + Tvec *__restrict accretion_l) { + u32 i_a = id_acc[id_a]; + Tvec r = xyz[i_a]; + Tvec v = vxyz[i_a]; + Tvec a = axyz[i_a]; + accretion_p[id_a] = gpart_mass * v; + accretion_mr[id_a] = gpart_mass * r; + accretion_ma[id_a] = gpart_mass * a; + + // dirty trick to account for the residual acceleration in the spin. This + // allows us to maitain a much better angular momentum conservation. + v += a * dt / 2; + accretion_l[id_a] = gpart_mass * sycl::cross(r - r_sink, v - v_sink); + }); + + Tvec acc_pxyz = shamalgs::primitives::sum(dev_sched, pxyz_acc, 0, Naccrete); + Tvec acc_mxyz = shamalgs::primitives::sum(dev_sched, mxyz_acc, 0, Naccrete); + Tvec acc_maxyz = shamalgs::primitives::sum(dev_sched, maxyz_acc, 0, Naccrete); + Tvec acc_lxyz = shamalgs::primitives::sum(dev_sched, lxyz_acc, 0, Naccrete); + + s_acc_mass += acc_mass; + s_acc_pxyz += acc_pxyz; + s_acc_mxyz += acc_mxyz; + s_acc_maxyz += acc_maxyz; + s_acc_lxyz += acc_lxyz; + } + }); + + Tscal sum_acc_mass = shamalgs::collective::allreduce_sum(s_acc_mass); + + // if there is accretion continue otherwise skip that part + if (sum_acc_mass <= 0) { + continue; + } + + Tvec sum_acc_pxyz = shamalgs::collective::allreduce_sum(s_acc_pxyz); + Tvec sum_acc_mxyz = shamalgs::collective::allreduce_sum(s_acc_mxyz); + Tvec sum_acc_maxyz = shamalgs::collective::allreduce_sum(s_acc_maxyz); + Tvec sum_acc_lxyz = shamalgs::collective::allreduce_sum(s_acc_lxyz); + + // compute the new sink values + Tscal new_mass = s.mass + sum_acc_mass; + Tvec new_pos = (sum_acc_mxyz + s.pos * s.mass) / (s.mass + sum_acc_mass); + Tvec new_vel = (sum_acc_pxyz + s.velocity * s.mass) / (s.mass + sum_acc_mass); + Tvec new_acc = (sum_acc_maxyz + s.sph_acceleration * s.mass) / (s.mass + sum_acc_mass); + Tvec new_ang_mom = s.angular_momentum + sum_acc_lxyz + - new_mass * sycl::cross(new_pos - s.pos, new_vel - s.velocity); + + // write back the updated sink state + auto new_state = s; + new_state.mass = new_mass; + new_state.pos = new_pos; + new_state.velocity = new_vel; + new_state.angular_momentum = new_ang_mom; + new_state.sph_acceleration = new_acc; + + had_accretion = true; + log += shambase::format( + "\n id {} deltas : mass={} r={} v={} l={}", + sink_id, + new_state.mass - s.mass, + new_state.pos - s.pos, + new_state.velocity - s.velocity, + new_state.angular_momentum - s.angular_momentum); + + s = new_state; + + // evict accreted particles from patches + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + u32 Nobj = pdat.get_obj_cnt(); + + sham::DeviceBuffer ¬_accreted + = accretion_flag_bufs.get(cur_p.id_patch).not_accreted; + sham::DeviceBuffer &accreted = accretion_flag_bufs.get(cur_p.id_patch).accreted; + + sham::DeviceBuffer &id_list_accrete = bufs_id_list_accrete.get(cur_p.id_patch); + + if (id_list_accrete.get_size() > 0) { + + sham::DeviceBuffer id_list_keep + = shamalgs::stream_compact(dev_sched, not_accreted, Nobj); + + pdat.keep_ids( + id_list_keep, shambase::narrow_or_throw(id_list_keep.get_size())); + } + }); + } + + if (shamcomm::world_rank() == 0 && had_accretion) { + logger::info_ln("sph::Sink", log); + } +} + +template class SPHKernel> +void shammodels::gsph::modules::SinkParticlesUpdate::predictor_step(Tscal dt) { + + StackEntry stack_loc{}; + + if (storage.sinks.is_empty()) { + return; + } + + compute_ext_forces(); + + std::vector &sink_parts = storage.sinks.get(); + + for (Sink &s : sink_parts) { + s.velocity += (dt / 2) * (s.sph_acceleration + s.ext_acceleration); + } + + for (Sink &s : sink_parts) { + s.pos += (dt) *s.velocity; + } +} + +template class SPHKernel> +void shammodels::gsph::modules::SinkParticlesUpdate::corrector_step(Tscal dt) { + + StackEntry stack_loc{}; + + if (storage.sinks.is_empty()) { + return; + } + + std::vector &sink_parts = storage.sinks.get(); + + for (Sink &s : sink_parts) { + s.velocity += (dt / 2) * (s.sph_acceleration + s.ext_acceleration); + } +} + +template class SPHKernel> +void shammodels::gsph::modules::SinkParticlesUpdate::compute_sph_forces() { + + StackEntry stack_loc{}; + + Tscal gpart_mass = solver_config.gpart_mass; + + if (storage.sinks.is_empty()) { + return; + } + + std::vector &sink_parts = storage.sinks.get(); + + Tscal G = solver_config.get_constant_G(); + Tscal epsilon_grav = 1e-9; + + using namespace shamrock; + using namespace shamrock::patch; + + PatchDataLayerLayout &pdl = scheduler().pdl_old(); + const u32 ixyz = pdl.get_field_idx("xyz"); + const u32 iaxyz_ext = pdl.get_field_idx("axyz_ext"); + + auto dev_sched = shamsys::instance::get_compute_scheduler_ptr(); + sham::DeviceQueue &q = shambase::get_check_ref(dev_sched).get_queue(); + + std::vector result_acc_sinks{}; + + for (Sink &s : sink_parts) { + + Tvec sph_acc_sink = {}; + + scheduler().for_each_patchdata_nonempty( + [&, G, epsilon_grav, gpart_mass](Patch cur_p, PatchDataLayer &pdat) { + sham::DeviceBuffer &buf_xyz = pdat.get_field_buf_ref(ixyz); + sham::DeviceBuffer &buf_axyz_ext = pdat.get_field_buf_ref(iaxyz_ext); + + sham::DeviceBuffer buf_sync_axyz(pdat.get_obj_cnt(), dev_sched); + + Tscal sink_mass = s.mass; + Tscal sink_racc = s.accretion_radius; + Tvec sink_pos = s.pos; + + sham::EventList depends_list; + auto xyz = buf_xyz.get_read_access(depends_list); + auto axyz_ext = buf_axyz_ext.get_write_access(depends_list); + auto axyz_sync = buf_sync_axyz.get_write_access(depends_list); + + auto e = q.submit( + depends_list, + [&, G, epsilon_grav, sink_mass, sink_pos, sink_racc](sycl::handler &cgh) { + shambase::parallel_for( + cgh, pdat.get_obj_cnt(), "sink-sph forces", [=](i32 id_a) { + Tvec r_a = xyz[id_a]; + + Tvec delta = r_a - sink_pos; + Tscal d = sycl::length(delta); + + Tvec force = G * delta / (d * d * d); + + // This is a hack to avoid the sink kaboom effect + // when the particle is being advected close to the sink before + // being accreted + if (d < sink_racc) { + force = {0, 0, 0}; + } + + axyz_sync[id_a] = force * gpart_mass; + axyz_ext[id_a] += -force * sink_mass; + }); + }); + + buf_xyz.complete_event_state(e); + buf_axyz_ext.complete_event_state(e); + buf_sync_axyz.complete_event_state(e); + + sph_acc_sink + += shamalgs::primitives::sum(dev_sched, buf_sync_axyz, 0, pdat.get_obj_cnt()); + }); + + result_acc_sinks.push_back(sph_acc_sink); + } + + std::vector gathered_result_acc_sinks{}; + shamalgs::collective::vector_allgatherv( + result_acc_sinks, gathered_result_acc_sinks, MPI_COMM_WORLD); + + u32 id_s = 0; + for (Sink &s : sink_parts) { + + s.sph_acceleration = {}; + + for (u32 rid = 0; rid < shamcomm::world_size(); rid++) { + s.sph_acceleration += gathered_result_acc_sinks[rid * sink_parts.size() + id_s]; + } + + id_s++; + } +} + +template class SPHKernel> +void shammodels::gsph::modules::SinkParticlesUpdate::compute_ext_forces() { + + StackEntry stack_loc{}; + + if (storage.sinks.is_empty()) { + return; + } + + std::vector &sink_parts = storage.sinks.get(); + + for (Sink &s : sink_parts) { + s.ext_acceleration = Tvec{}; + } + + Tscal G = solver_config.get_constant_G(); + Tscal epsilon_grav_sink = 1e-9; + + for (Sink &s1 : sink_parts) { + Tvec sum{}; + for (Sink &s2 : sink_parts) { + Tvec rij = s1.pos - s2.pos; + Tscal rij_scal = sycl::length(rij); + sum -= G * s2.mass * rij / (rij_scal * rij_scal * rij_scal + epsilon_grav_sink); + } + s1.ext_acceleration = sum; + } +} + +using namespace shammath; +template class shammodels::gsph::modules::SinkParticlesUpdate; +template class shammodels::gsph::modules::SinkParticlesUpdate; +template class shammodels::gsph::modules::SinkParticlesUpdate; + +template class shammodels::gsph::modules::SinkParticlesUpdate; +template class shammodels::gsph::modules::SinkParticlesUpdate; +template class shammodels::gsph::modules::SinkParticlesUpdate; From b8165c0e65f7f2da78f9e47bbf7c4aed9419f728 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Fri, 31 Jul 2026 17:28:17 +0900 Subject: [PATCH 3/7] add sink-sink cfl condition --- .../gsph/include/shammodels/gsph/Solver.hpp | 7 +++ .../include/shammodels/gsph/SolverConfig.hpp | 1 + src/shammodels/gsph/src/Solver.cpp | 46 +++++++++++++++++++ 3 files changed, 54 insertions(+) diff --git a/src/shammodels/gsph/include/shammodels/gsph/Solver.hpp b/src/shammodels/gsph/include/shammodels/gsph/Solver.hpp index b03a434a29..d60591fe80 100644 --- a/src/shammodels/gsph/include/shammodels/gsph/Solver.hpp +++ b/src/shammodels/gsph/include/shammodels/gsph/Solver.hpp @@ -216,6 +216,13 @@ namespace shammodels::gsph { */ Tscal compute_dt_cfl(); + /** + * @brief Compute sink timestep constraint + * + * @return Minimum CFL timestep across all particles + */ + Tscal compute_sink_cfl(); + bool apply_corrector(Tscal dt, u64 Npart_all); void update_sync_load_values(); diff --git a/src/shammodels/gsph/include/shammodels/gsph/SolverConfig.hpp b/src/shammodels/gsph/include/shammodels/gsph/SolverConfig.hpp index 283c673c08..719e0488ce 100644 --- a/src/shammodels/gsph/include/shammodels/gsph/SolverConfig.hpp +++ b/src/shammodels/gsph/include/shammodels/gsph/SolverConfig.hpp @@ -72,6 +72,7 @@ namespace shammodels::gsph { struct CFLConfig { Tscal cfl_cour = 0.3; ///< CFL condition for the courant factor Tscal cfl_force = 0.25; ///< CFL condition for the force + Tscal eta_sink = 0.05; }; } // namespace shammodels::gsph diff --git a/src/shammodels/gsph/src/Solver.cpp b/src/shammodels/gsph/src/Solver.cpp index a985d18986..fdb6d875c1 100644 --- a/src/shammodels/gsph/src/Solver.cpp +++ b/src/shammodels/gsph/src/Solver.cpp @@ -1647,6 +1647,52 @@ typename shammodels::gsph::Solver::Tscal shammodels::gsph::Solver class Kern> +typename shammodels::gsph::Solver::Tscal shammodels::gsph::Solver:: + compute_sink_cfl() { + StackEntry stack_loc{}; + + Tscal C_cour = solver_config.cfl_config.cfl_cour; + Tscal C_force = solver_config.cfl_config.cfl_force; + Tscal eta_phi = solver_config.cfl_config.eta_sink; + + Tscal sink_sink_cfl = shambase::get_infty(); + + Tscal G = solver_config.get_constant_G(); + + std::vector> &sink_parts = storage.sinks.get(); + + for (u32 i = 0; i < sink_parts.size(); i++) { + sph::SinkParticle &s_i = sink_parts[i]; + Tscal sink_sink_cfl_i = shambase::get_infty(); + + Tvec f_i = s_i.ext_acceleration; + + Tscal grad_phi_i_sq = sham::dot(f_i, f_i); // m^2.s^-4 + + if (grad_phi_i_sq == 0) { + continue; + } + + for (u32 j = 0; j < sink_parts.size(); j++) { + sph::SinkParticle &s_j = sink_parts[j]; + if (i == j) { + continue; + } + Tvec rij = s_i.pos - s_j.pos; + Tscal rij_scal = sycl::length(rij); + Tscal phi_ij = G * s_j.mass / rij_scal; // J / kg = m^2.s^-2 + Tscal term_ij = sham::abs(phi_ij) / grad_phi_i_sq; // s^2 + Tscal dt_ij = C_force * eta_phi * sycl::sqrt(term_ij); // s + sink_sink_cfl_i = sham::min(sink_sink_cfl_i, dt_ij); + } + + sink_sink_cfl = sham::min(sink_sink_cfl, sink_sink_cfl_i); + } + + return sink_sink_cfl; +} + template class Kern> bool shammodels::gsph::Solver::apply_corrector(Tscal dt, u64 Npart_all) { StackEntry stack_loc{}; From f5168466ff2f8a9d842284eb24c9a651e006841a Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Mon, 3 Aug 2026 00:03:06 +0900 Subject: [PATCH 4/7] author update --- .../include/shammodels/gsph/modules/setup/GeneratorMCDisc.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/shammodels/gsph/include/shammodels/gsph/modules/setup/GeneratorMCDisc.hpp b/src/shammodels/gsph/include/shammodels/gsph/modules/setup/GeneratorMCDisc.hpp index 1e8587818e..48b3baa3ca 100644 --- a/src/shammodels/gsph/include/shammodels/gsph/modules/setup/GeneratorMCDisc.hpp +++ b/src/shammodels/gsph/include/shammodels/gsph/modules/setup/GeneratorMCDisc.hpp @@ -11,6 +11,7 @@ /** * @file GeneratorMCDisc.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) * @brief * */ From 928ffbf18713cccf14c52dcf6c43adf0b8727750 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Mon, 3 Aug 2026 17:38:53 +0900 Subject: [PATCH 5/7] minor fixes --- .../gsph/include/shammodels/gsph/Solver.hpp | 2 +- .../include/shammodels/gsph/SolverConfig.hpp | 6 ++-- .../shammodels/gsph/config/FieldNames.hpp | 3 ++ .../shammodels/gsph/modules/GSPHSetup.hpp | 5 +-- src/shammodels/gsph/src/Solver.cpp | 7 +++- src/shammodels/gsph/src/SolverConfig.cpp | 2 +- .../gsph/src/modules/ExternalForces.cpp | 24 ------------- src/shammodels/gsph/src/modules/GSPHSetup.cpp | 3 +- .../gsph/src/modules/GeneratorMCDisc.cpp | 2 +- src/shammodels/gsph/src/pyGSPHModel.cpp | 34 +++++++++---------- 10 files changed, 32 insertions(+), 56 deletions(-) diff --git a/src/shammodels/gsph/include/shammodels/gsph/Solver.hpp b/src/shammodels/gsph/include/shammodels/gsph/Solver.hpp index d60591fe80..54d38cfd73 100644 --- a/src/shammodels/gsph/include/shammodels/gsph/Solver.hpp +++ b/src/shammodels/gsph/include/shammodels/gsph/Solver.hpp @@ -212,7 +212,7 @@ namespace shammodels::gsph { * - Courant condition: dt_cour = C_cour * h / vsig * - Force condition: dt_force = C_force * sqrt(h / |a|) * - * @return Minimum CFL timestep across all particles + * @return Minimum CFL timestep across all sinks */ Tscal compute_dt_cfl(); diff --git a/src/shammodels/gsph/include/shammodels/gsph/SolverConfig.hpp b/src/shammodels/gsph/include/shammodels/gsph/SolverConfig.hpp index 968a94a498..b414574ff2 100644 --- a/src/shammodels/gsph/include/shammodels/gsph/SolverConfig.hpp +++ b/src/shammodels/gsph/include/shammodels/gsph/SolverConfig.hpp @@ -226,11 +226,9 @@ struct shammodels::gsph::SolverConfig { inline void set_eos_isothermal(Tscal cs) { eos_config.set_isothermal(cs); } /** - * @brief Set the EOS configuration to a locally isothermal equation of state fromFarris 2014 + * `@brief` Set the EOS configuration to a locally isothermal equation of state from Farris 2014 * - * @param cs0 Soundspeed at the reference radius - * @param q Power exponent of the soundspeed profile - * @param r0 Reference radius + * `@param` h_over_r Disc aspect ratio */ inline void set_eos_locally_isothermalFA2014(Tscal h_over_r) { eos_config.set_locally_isothermalFA2014(h_over_r); diff --git a/src/shammodels/gsph/include/shammodels/gsph/config/FieldNames.hpp b/src/shammodels/gsph/include/shammodels/gsph/config/FieldNames.hpp index f8af1dbe7c..d3501f8616 100644 --- a/src/shammodels/gsph/include/shammodels/gsph/config/FieldNames.hpp +++ b/src/shammodels/gsph/include/shammodels/gsph/config/FieldNames.hpp @@ -45,6 +45,9 @@ namespace shammodels::gsph::names { /// 3-acceleration field inline constexpr const char *axyz = "axyz"; + /// 3-acceleration field due to external forces + inline constexpr const char *axyz_ext = "axyz_ext"; + /// Specific internal energy u inline constexpr const char *uint = "uint"; diff --git a/src/shammodels/gsph/include/shammodels/gsph/modules/GSPHSetup.hpp b/src/shammodels/gsph/include/shammodels/gsph/modules/GSPHSetup.hpp index f59db0b6f6..af1f5b2192 100644 --- a/src/shammodels/gsph/include/shammodels/gsph/modules/GSPHSetup.hpp +++ b/src/shammodels/gsph/include/shammodels/gsph/modules/GSPHSetup.hpp @@ -43,10 +43,7 @@ namespace shammodels::gsph::modules { GSPHSetup(ShamrockCtx &context, Config &solver_config, Storage &storage) : context(context), solver_config(solver_config), storage(storage) {} - void apply_setup( - SetupNodePtr setup, - bool part_reordering, - std::optional insert_step = std::nullopt); + void apply_setup(SetupNodePtr setup, std::optional insert_step = std::nullopt); std::shared_ptr make_generator_disc_mc( Tscal part_mass, diff --git a/src/shammodels/gsph/src/Solver.cpp b/src/shammodels/gsph/src/Solver.cpp index b7be4f2658..4a96fa6da9 100644 --- a/src/shammodels/gsph/src/Solver.cpp +++ b/src/shammodels/gsph/src/Solver.cpp @@ -1682,7 +1682,6 @@ typename shammodels::gsph::Solver::Tscal shammodels::gsph::Solver::Tscal shammodels::gsph::Solver> &sink_parts = storage.sinks.get(); + if (storage.sinks.is_empty()) { + return sink_sink_cfl; + } + for (u32 i = 0; i < sink_parts.size(); i++) { sph::SinkParticle &s_i = sink_parts[i]; Tscal sink_sink_cfl_i = shambase::get_infty(); @@ -1917,9 +1920,11 @@ shammodels::gsph::TimestepLog shammodels::gsph::Solver::evolve_once( // STEP 7: CFL - compute next timestep Tscal dt_next = compute_dt_cfl(); + Tscal dt_cfl = compute_dt_cfl(); // Ensure dt doesn't grow too fast (max 2x per step), but allow any value if dt was 0 if (dt > Tscal(0)) { + dt_next = sham::min(dt_next, dt_cfl); dt_next = sham::min(dt_next, Tscal(2) * dt); } diff --git a/src/shammodels/gsph/src/SolverConfig.cpp b/src/shammodels/gsph/src/SolverConfig.cpp index bc1ccead3a..e4133d0d96 100644 --- a/src/shammodels/gsph/src/SolverConfig.cpp +++ b/src/shammodels/gsph/src/SolverConfig.cpp @@ -35,7 +35,7 @@ void shammodels::gsph::SolverConfig::set_layout( // Smoothing length pdl.add_field(names::common::hpart, 1); - pdl.add_field("axyz_ext", 1); + pdl.add_field(names::newtonian::axyz_ext, 1); // Internal energy (for adiabatic EOS) if (has_field_uint()) { diff --git a/src/shammodels/gsph/src/modules/ExternalForces.cpp b/src/shammodels/gsph/src/modules/ExternalForces.cpp index 32c59ee54c..c7c701013b 100644 --- a/src/shammodels/gsph/src/modules/ExternalForces.cpp +++ b/src/shammodels/gsph/src/modules/ExternalForces.cpp @@ -238,32 +238,8 @@ void shammodels::gsph::modules::ExternalForces::add_ext_forces( using namespace shamrock::solvergraph; SolverGraph solver_graph{}; - auto set_constant_G = register_constant_set(solver_graph, "constant_G", [&]() { - return solver_config.get_constant_G(); - }); - auto set_constant_c = register_constant_set(solver_graph, "constant_c", [&]() { - return solver_config.get_constant_c(); - }); - - bool is_G_needed = false; - bool is_c_needed = false; - - for (auto var_force : solver_config.ext_force_config.ext_forces) { - if (EF_PointMass *ext_force = std::get_if(&var_force.val)) { - } else { - shambase::throw_unimplemented("this force is not handled, yet ..."); - } - } - std::vector> add_ext_forces_seq{}; - if (is_G_needed) { - add_ext_forces_seq.push_back(set_constant_G); - } - if (is_c_needed) { - add_ext_forces_seq.push_back(set_constant_c); - } - auto field_xyz = solver_graph.register_edge("field_xyz", FieldRefs("", "")); auto field_vxyz = solver_graph.register_edge("field_vxyz", FieldRefs("", "")); auto field_axyz = solver_graph.register_edge("field_axyz", FieldRefs("", "")); diff --git a/src/shammodels/gsph/src/modules/GSPHSetup.cpp b/src/shammodels/gsph/src/modules/GSPHSetup.cpp index 23aa223df6..4cf5fa8367 100644 --- a/src/shammodels/gsph/src/modules/GSPHSetup.cpp +++ b/src/shammodels/gsph/src/modules/GSPHSetup.cpp @@ -30,7 +30,6 @@ #include "shammodels/gsph/modules/ComputeLoadBalanceValue.hpp" #include "shammodels/gsph/modules/GSPHSetup.hpp" #include "shammodels/gsph/modules/setup/GeneratorMCDisc.hpp" -#include "shammodels/sph/modules/ParticleReordering.hpp" #include "shamrock/patch/PatchDataLayer.hpp" #include "shamrock/scheduler/DataInserterUtility.hpp" #include "shamsys/NodeInstance.hpp" @@ -67,7 +66,7 @@ inline std::shared_ptr shammodels::gs template class SPHKernel> void shammodels::gsph::modules::GSPHSetup::apply_setup( - SetupNodePtr setup, bool part_reordering, std::optional insert_step) { + SetupNodePtr setup, std::optional insert_step) { if (!bool(setup)) { shambase::throw_with_loc("The setup shared pointer is empty"); diff --git a/src/shammodels/gsph/src/modules/GeneratorMCDisc.cpp b/src/shammodels/gsph/src/modules/GeneratorMCDisc.cpp index 6f65ec6eb9..b9b385b32e 100644 --- a/src/shammodels/gsph/src/modules/GeneratorMCDisc.cpp +++ b/src/shammodels/gsph/src/modules/GeneratorMCDisc.cpp @@ -115,7 +115,7 @@ shamrock::patch::PatchDataLayer shammodels::gsph::modules::GeneratorMCDisc vec_cs; if (need_cs) { diff --git a/src/shammodels/gsph/src/pyGSPHModel.cpp b/src/shammodels/gsph/src/pyGSPHModel.cpp index 6afc8c1afb..4b5dfcae87 100644 --- a/src/shammodels/gsph/src/pyGSPHModel.cpp +++ b/src/shammodels/gsph/src/pyGSPHModel.cpp @@ -42,9 +42,9 @@ void add_gsph_instance(py::module &m, std::string name_config, std::string name_ using Tscal = shambase::VecComponent; - using T = Model; - using TSPHSetup = modules::GSPHSetup; - using TConfig = typename T::SolverConfig; + using T = Model; + using TGSPHSetup = modules::GSPHSetup; + using TConfig = typename T::SolverConfig; shamlog_debug_ln("[Py]", "registering class :", name_config, typeid(T).name()); shamlog_debug_ln("[Py]", "registering class :", name_model, typeid(T).name()); @@ -235,11 +235,11 @@ void add_gsph_instance(py::module &m, std::string name_config, std::string name_ py::arg("split_load_value"), py::arg("merge_load_value")); - std::string setup_name = name_model + "_SPHSetup"; - py::class_(m, setup_name.c_str()) + std::string setup_name = name_model + "_GSPHSetup"; + py::class_(m, setup_name.c_str()) .def( "make_generator_disc_mc", - [](TSPHSetup &self, + [](TGSPHSetup &self, Tscal part_mass, Tscal disc_mass, Tscal r_in, @@ -286,7 +286,7 @@ void add_gsph_instance(py::module &m, std::string name_config, std::string name_ if (cs_field) { if (shamcomm::world_rank() == 0) { logger::warn_ln( - "SPHSetup", + "GSPHSetup", "make_generator_disc_mc: with the current EOS, cs_field is " "ignored"); } @@ -294,7 +294,7 @@ void add_gsph_instance(py::module &m, std::string name_config, std::string name_ if (cs_profile) { if (shamcomm::world_rank() == 0) { logger::warn_ln( - "SPHSetup", + "GSPHSetup", "make_generator_disc_mc: with the current EOS, cs_profile is " "ignored"); } @@ -360,7 +360,7 @@ void add_gsph_instance(py::module &m, std::string name_config, std::string name_ smoothing lengths are set from that density. Args: - part_mass: Mass of each SPH particle. + part_mass: Mass of each GSPH particle. disc_mass: Total disc mass. The particle count is ``disc_mass / part_mass``. r_in: Inner disc radius. r_out: Outer disc radius. @@ -392,9 +392,8 @@ void add_gsph_instance(py::module &m, std::string name_config, std::string name_ )pbdoc") .def( "apply_setup", - [](TSPHSetup &self, + [](TGSPHSetup &self, modules::SetupNodePtr setup, - bool part_reordering, std::optional gen_step, std::optional insert_step, std::optional msg_count_limit, @@ -405,33 +404,32 @@ void add_gsph_instance(py::module &m, std::string name_config, std::string name_ bool speculative_balancing) { if (bool(gen_step)) { ON_RANK_0( - logger::warn_ln("SPHSetup", "gen_step is ignored when using old setup")); + logger::warn_ln("GSPHSetup", "gen_step is ignored when using old setup")); } if (bool(msg_count_limit)) { ON_RANK_0( logger::warn_ln( - "SPHSetup", "msg_count_limit is ignored when using old setup")); + "GSPHSetup", "msg_count_limit is ignored when using old setup")); } if (bool(msg_size_limit)) { ON_RANK_0( logger::warn_ln( - "SPHSetup", "msg_size_limit is ignored when using old setup")); + "GSPHSetup", "msg_size_limit is ignored when using old setup")); } if (bool(max_msg_size)) { ON_RANK_0( logger::warn_ln( - "SPHSetup", "max_msg_size is ignored when using old setup")); + "GSPHSetup", "max_msg_size is ignored when using old setup")); } if (bool(do_setup_log)) { ON_RANK_0( logger::warn_ln( - "SPHSetup", "do_setup_log is ignored when using old setup")); + "GSPHSetup", "do_setup_log is ignored when using old setup")); } - return self.apply_setup(setup, part_reordering, insert_step); + return self.apply_setup(setup, insert_step); }, py::arg("setup"), py::kw_only(), - py::arg("part_reordering") = true, py::arg("gen_step") = std::nullopt, py::arg("insert_step") = std::nullopt, py::arg("msg_count_limit") = std::nullopt, From 32f873f5fc9c3047d4fc24dfc4e67bef17a10a84 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Tue, 4 Aug 2026 17:02:53 +0900 Subject: [PATCH 6/7] add sinks to json header --- src/shammodels/gsph/include/shammodels/gsph/Model.hpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/shammodels/gsph/include/shammodels/gsph/Model.hpp b/src/shammodels/gsph/include/shammodels/gsph/Model.hpp index aa3396b09e..4fd2a74f1f 100644 --- a/src/shammodels/gsph/include/shammodels/gsph/Model.hpp +++ b/src/shammodels/gsph/include/shammodels/gsph/Model.hpp @@ -429,6 +429,12 @@ namespace shammodels::gsph { nlohmann::json metadata; metadata["solver_config"] = solver.solver_config; + if (solver.storage.sinks.is_empty()) { + metadata["sinks"] = nlohmann::json{}; + } else { + metadata["sinks"] = solver.storage.sinks.get(); + } + shamrock::write_shamrock_dump( fname, metadata.dump(4), shambase::get_check_ref(ctx.sched)); } From 5ac9143a71e89b1a68d6f9ca1211cc4bb8476520 Mon Sep 17 00:00:00 2001 From: Yona Lapeyre Date: Tue, 4 Aug 2026 17:54:39 +0900 Subject: [PATCH 7/7] abort the merge --- .../common/setup/GeneratorMCDisc.hpp | 173 ---------------- .../src/modules/setup/GeneratorMCDisc.cpp | 185 ------------------ 2 files changed, 358 deletions(-) delete mode 100644 src/shammodels/common/include/shammodels/common/setup/GeneratorMCDisc.hpp delete mode 100644 src/shammodels/common/src/modules/setup/GeneratorMCDisc.cpp diff --git a/src/shammodels/common/include/shammodels/common/setup/GeneratorMCDisc.hpp b/src/shammodels/common/include/shammodels/common/setup/GeneratorMCDisc.hpp deleted file mode 100644 index b7258c57a8..0000000000 --- a/src/shammodels/common/include/shammodels/common/setup/GeneratorMCDisc.hpp +++ /dev/null @@ -1,173 +0,0 @@ -// -------------------------------------------------------// -// -// SHAMROCK code for hydrodynamics -// Copyright (c) 2021-2026 Timothée David--Cléris -// 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 - -/** - * @file GeneratorMCDisc.hpp - * @author Timothée David--Cléris (tim.shamrock@proton.me) - * @brief - * - */ - -#include "shambase/constants.hpp" -#include "shamalgs/collective/InvariantParallelGenerator.hpp" -#include "shamalgs/collective/indexing.hpp" -#include "shamrock/scheduler/ShamrockCtx.hpp" - -namespace shammodels::sph::modules { - - template class SPHKernel, class TSolverConfig, class TSetupNodeBase> - class GeneratorMCDisc : public TSetupNodeBase { - using Tscal = shambase::VecComponent; - static constexpr u32 dim = shambase::VectorProperties::dimension; - using Kernel = SPHKernel; - - using Config = TSolverConfig; - - ShamrockCtx &context; - Config &solver_config; - - struct DiscOutput { - sycl::vec pos; - Tscal rho; - }; - - Tscal pmass; - - class DiscIterator; - DiscIterator generator; - Tscal init_h_factor; - - std::function vel_profile; - std::function cs_profile; - - static DiscIterator make_generator( - Tscal part_mass, - Tscal disc_mass, - Tscal r_in, - Tscal r_out, - std::function sigma_profile, - std::function H_profile, - std::mt19937_64 eng) { - return DiscIterator(part_mass, disc_mass, r_in, r_out, sigma_profile, H_profile, eng); - } - - public: - GeneratorMCDisc( - ShamrockCtx &context, - Config &solver_config, - Tscal part_mass, - Tscal disc_mass, - Tscal r_in, - Tscal r_out, - std::function sigma_profile, - std::function H_profile, - std::function vel_profile, - std::function cs_profile, - std::mt19937_64 eng, - Tscal init_h_factor) - : context(context), solver_config(solver_config), - generator( - make_generator(part_mass, disc_mass, r_in, r_out, sigma_profile, H_profile, eng)), - init_h_factor(init_h_factor), pmass(part_mass), vel_profile(vel_profile), - cs_profile(cs_profile) {} - - bool is_done(); - - shamrock::patch::PatchDataLayer next_n(u32 nmax); - - std::string get_name() { return "GeneratorMCDisc"; } - ISPHSetupNode_Dot get_dot_subgraph() { return ISPHSetupNode_Dot{get_name(), 0, {}}; } - }; - -} // namespace shammodels::sph::modules - -template class SPHKernel, class TConfig, class TSetupNodeBase> -class shammodels::common::GeneratorMCDisc::DiscIterator { - - bool done = false; - u64 current_index = 0; - - Tscal part_mass; - Tscal disc_mass; - u64 Npart; - - Tscal r_in; - Tscal r_out; - std::function sigma_profile; - std::function H_profile; - - shamalgs::collective::InvariantParallelGenerator generator; - - static constexpr Tscal _2pi = 2 * shambase::constants::pi; - - Tscal f_func(Tscal r) { return r * sigma_profile(r); } - - DiscOutput next(u64 seed); - - public: - DiscIterator( - Tscal part_mass, - Tscal disc_mass, - Tscal r_in, - Tscal r_out, - std::function sigma_profile, - std::function H_profile, - std::mt19937_64 eng) - : DiscIterator( - part_mass, - disc_mass, - r_in, - r_out, - sigma_profile, - H_profile, - eng, - disc_mass / part_mass) {} - - DiscIterator( - Tscal part_mass, - Tscal disc_mass, - Tscal r_in, - Tscal r_out, - std::function sigma_profile, - std::function H_profile, - std::mt19937_64 eng, - u64 Npart) - : part_mass(part_mass), disc_mass(disc_mass), Npart(Npart), r_in(r_in), r_out(r_out), - sigma_profile(sigma_profile), H_profile(H_profile), generator(eng, Npart), - current_index(0) { - - shamlog_debug_ln( - "GeneratorMCDisc", - "part_mass", - part_mass, - "disc_mass", - disc_mass, - "r_in", - r_in, - "r_out", - r_out, - "Npart", - Npart); - } - - inline bool is_done() { - return generator.is_done(); - } // just to make sure the result is not tempered with - - inline std::vector next_n(u64 nmax) { - std::vector seeds = generator.next_n(nmax); - std::vector ret{}; - for (u64 seed : seeds) { - ret.push_back(next(seed)); - } - return ret; - } -}; diff --git a/src/shammodels/common/src/modules/setup/GeneratorMCDisc.cpp b/src/shammodels/common/src/modules/setup/GeneratorMCDisc.cpp deleted file mode 100644 index 8144a2cd2c..0000000000 --- a/src/shammodels/common/src/modules/setup/GeneratorMCDisc.cpp +++ /dev/null @@ -1,185 +0,0 @@ -// -------------------------------------------------------// -// -// SHAMROCK code for hydrodynamics -// Copyright (c) 2021-2026 Timothée David--Cléris -// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1 -// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information -// -// -------------------------------------------------------// - -/** - * @file GeneratorMCDisc.cpp - * @author Timothée David--Cléris (tim.shamrock@proton.me) - * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr) - * @brief - * - */ - -#include "shambase/constants.hpp" -#include "shamalgs/collective/indexing.hpp" -#include "shamalgs/random.hpp" -#include "shammodels/common/setup/GeneratorMCDisc.hpp" -#include "shammodels/sph/math/density.hpp" - -template< - class Tvec, - template class SPHKernel, - class TConfig, - class TSetupNodeBase - >> auto shammodels::sph::modules::GeneratorMCDisc::DiscIterator::next( - u64 seed) -> DiscOutput { - - std::mt19937_64 eng_local(seed); // ensure that 1 part = 1 random draw - - Tscal fmax = f_func(r_out); - - auto find_r = [&]() { - while (true) { - Tscal u2 = shamalgs::primitives::mock_value(eng_local, 0, fmax); - Tscal r = shamalgs::primitives::mock_value(eng_local, r_in, r_out); - if (u2 < f_func(r)) { - return r; - } - } - }; - - auto theta = shamalgs::primitives::mock_value(eng_local, 0, _2pi); - auto Gauss = shamalgs::random::mock_gaussian(eng_local); - - // depends on sigma profile - Tscal r = find_r(); - Tscal sigma = sigma_profile(r); - - // depends on H profile & sigma profile (through r) - Tscal H = H_profile(r); - Tscal z = H * Gauss; - - auto pos = sycl::vec{r * sycl::cos(theta), r * sycl::sin(theta), z}; - - // extrapolate the density from sigma profile - Tscal fs = 1; - Tscal rho = (sigma * fs) * sycl::exp(-z * z / (2 * H * H)); - - DiscOutput out{.pos = pos, .rho = rho}; - - // increase counter + check if finished - current_index++; - if (current_index == Npart) { - done = true; - } - - return out; -} - -template class SPHKernel> -bool shammodels::sph::modules::GeneratorMCDisc::is_done() { - return generator.is_done(); -} - -template class SPHKernel> -shamrock::patch::PatchDataLayer shammodels::sph::modules::GeneratorMCDisc::next_n( - u32 nmax) { - - using namespace shamrock::patch; - PatchScheduler &sched = shambase::get_check_ref(context.sched); - - std::vector pos_data; - - // Fill pos_data if the scheduler has some patchdata in this rank - if (!generator.is_done()) { - u64 loc_gen_count = nmax; - pos_data = generator.next_n(loc_gen_count); - } - - // extract data from disc output - std::vector vec_pos; - std::vector vec_rho; - - vec_pos.reserve(pos_data.size()); - vec_rho.reserve(pos_data.size()); - - for (DiscOutput o : pos_data) { - vec_pos.push_back(o.pos); - vec_rho.push_back(o.rho); - } - - // compute the hpart from the rho - std::vector vec_h; - vec_h.reserve(pos_data.size()); - for (Tscal rho : vec_rho) { - vec_h.push_back(shamrock::sph::h_rho(pmass, rho, Kernel::hfactd) * init_h_factor); - } - - // compute velocities - std::vector vec_vel; - vec_vel.reserve(pos_data.size()); - for (size_t i = 0; i < vec_pos.size(); i++) { - Tvec vel = vel_profile(vec_pos[i]); - vec_vel.push_back(vel); - } - - // compute the cs - bool need_cs = solver_config.is_eos_locally_isothermal(); - - std::vector vec_cs; - if (need_cs) { - if (!cs_profile) { - throw shambase::make_except_with_loc( - "With this EOS you need to provide a cs_profile"); - } - vec_cs.reserve(pos_data.size()); - for (size_t i = 0; i < vec_pos.size(); i++) { - Tscal cs = cs_profile(vec_pos[i]); - vec_cs.push_back(cs); - } - } - - // Make a patchdata from pos_data - PatchDataLayer tmp(sched.get_layout_ptr_old()); - if (!pos_data.empty()) { - tmp.resize(pos_data.size()); - tmp.fields_raz(); - - { - u32 len = pos_data.size(); - PatchDataField &f - = tmp.get_field(sched.pdl_old().get_field_idx("xyz")); - sycl::buffer buf(vec_pos.data(), len); - f.override(buf, len); - } - - { - u32 len = pos_data.size(); - PatchDataField &f - = tmp.get_field(sched.pdl_old().get_field_idx("vxyz")); - sycl::buffer buf(vec_vel.data(), len); - f.override(buf, len); - } - { - u32 len = vec_pos.size(); - PatchDataField &f - = tmp.get_field(sched.pdl_old().get_field_idx("hpart")); - sycl::buffer buf(vec_h.data(), len); - f.override(buf, len); - } - - if (need_cs) { - u32 len = vec_pos.size(); - PatchDataField &f - = tmp.get_field(sched.pdl_old().get_field_idx("soundspeed")); - sycl::buffer buf(vec_cs.data(), len); - f.override(buf, len); - } - } - - return tmp; -} - -using namespace shammath; -template class shammodels::sph::modules::GeneratorMCDisc; -template class shammodels::sph::modules::GeneratorMCDisc; -template class shammodels::sph::modules::GeneratorMCDisc; - -template class shammodels::sph::modules::GeneratorMCDisc; -template class shammodels::sph::modules::GeneratorMCDisc; -template class shammodels::sph::modules::GeneratorMCDisc;