From 5a10f3cef9bd05aa38593325335fa71cbdd29fb3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Thu, 5 Feb 2026 15:23:52 +0100 Subject: [PATCH 01/24] [SPH] ensure that dust can only be enabled in experimental mode --- .../sph/include/shammodels/sph/SolverConfig.hpp | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp index 7012dc42f7..15d3c3f802 100644 --- a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp +++ b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp @@ -155,10 +155,15 @@ namespace shammodels::sph { bool is_not_none = bool(std::get_if(¤t_mode)) || bool(std::get_if(¤t_mode)); if (is_not_none) { - ON_RANK_0( - logger::warn_ln( - "SPH::config", - "Dust config != None is work in progress, use it at your own risk")); + if (!shamrock::are_experimental_features_allowed()) { + shambase::throw_with_loc( + "Dust config != None is experimental"); + } else { + ON_RANK_0( + logger::warn_ln( + "SPH::config", + "Dust config != None is work in progress, use it at your own risk")); + } } } }; From bb0d28fc67c9b0a2af12776dfda53a3b6a003a05 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Thu, 5 Feb 2026 15:33:14 +0100 Subject: [PATCH 02/24] add S_j to config --- .../sph/include/shammodels/sph/SolverConfig.hpp | 9 +++++++++ src/shammodels/sph/src/SolverConfig.cpp | 5 +++++ 2 files changed, 14 insertions(+) diff --git a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp index 15d3c3f802..daf2f57bfb 100644 --- a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp +++ b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp @@ -127,6 +127,15 @@ namespace shammodels::sph { inline void set_monofluid_tvi(u32 nvar) { current_mode = MonofluidTVI{nvar}; } inline void set_monofluid_complete(u32 nvar) { current_mode = MonofluidComplete{nvar}; } + inline bool is_monofluid_tvi() { return bool(std::get_if(¤t_mode)); } + inline bool is_monofluid_complete() { + return bool(std::get_if(¤t_mode)); + } + + inline bool has_Sj_field() { + return is_monofluid_tvi(); // S_j = sqrt(\rho \epsilon_j) + } + inline bool has_epsilon_field() { return bool(std::get_if(¤t_mode)) || bool(std::get_if(¤t_mode)); diff --git a/src/shammodels/sph/src/SolverConfig.cpp b/src/shammodels/sph/src/SolverConfig.cpp index 4259af596e..e7f13e2d46 100644 --- a/src/shammodels/sph/src/SolverConfig.cpp +++ b/src/shammodels/sph/src/SolverConfig.cpp @@ -111,6 +111,11 @@ namespace shammodels::sph { if (should_save_dt_to_fields()) { pdl.add_field("dt_part", 1); } + + if (dust_config.has_Sj_field()) { + u32 ndust = dust_config.get_dust_nvar(); + pdl.add_field("Sj", ndust); + } } template class SPHKernel> From 8676e83974271ffcb51738cd40b0e1b118b701f0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Thu, 5 Feb 2026 16:01:10 +0100 Subject: [PATCH 03/24] monofluid formula --- .../shammodels/sph/modules/UpdateDerivs.hpp | 4 + src/shammodels/sph/src/SolverConfig.cpp | 8 +- .../sph/src/modules/UpdateDerivs.cpp | 185 +++++++++++++++++- 3 files changed, 194 insertions(+), 3 deletions(-) diff --git a/src/shammodels/sph/include/shammodels/sph/modules/UpdateDerivs.hpp b/src/shammodels/sph/include/shammodels/sph/modules/UpdateDerivs.hpp index 3075248eb5..78fc5a54b7 100644 --- a/src/shammodels/sph/include/shammodels/sph/modules/UpdateDerivs.hpp +++ b/src/shammodels/sph/include/shammodels/sph/modules/UpdateDerivs.hpp @@ -60,6 +60,10 @@ namespace shammodels::sph::modules { void update_derivs_cd10(VaryingCD10 cfg); void update_derivs_disc_visco(ConstantDisc cfg); + using DustConfig = typename Config::DustConfig; + + void update_derivs_dust_monofluid_tvi_Sj(DustConfig cfg); + using Cfg_MHD = typename Config::MHDConfig; using NoneMHD = typename Cfg_MHD::None; diff --git a/src/shammodels/sph/src/SolverConfig.cpp b/src/shammodels/sph/src/SolverConfig.cpp index e7f13e2d46..85ba030a2e 100644 --- a/src/shammodels/sph/src/SolverConfig.cpp +++ b/src/shammodels/sph/src/SolverConfig.cpp @@ -114,7 +114,8 @@ namespace shammodels::sph { if (dust_config.has_Sj_field()) { u32 ndust = dust_config.get_dust_nvar(); - pdl.add_field("Sj", ndust); + pdl.add_field("s_j", ndust); + pdl.add_field("ds_j_dt", ndust); } } @@ -156,6 +157,11 @@ namespace shammodels::sph { u32 ndust = dust_config.get_dust_nvar(); ghost_layout.add_field("deltav", ndust); } + + if (dust_config.has_Sj_field()) { + u32 ndust = dust_config.get_dust_nvar(); + ghost_layout.add_field("s_j", ndust); + } } }; // namespace shammodels::sph diff --git a/src/shammodels/sph/src/modules/UpdateDerivs.cpp b/src/shammodels/sph/src/modules/UpdateDerivs.cpp index b9cef7ce21..e59f1eb729 100644 --- a/src/shammodels/sph/src/modules/UpdateDerivs.cpp +++ b/src/shammodels/sph/src/modules/UpdateDerivs.cpp @@ -27,8 +27,9 @@ template class SPHKernel> void shammodels::sph::modules::UpdateDerivs::update_derivs() { - Cfg_AV cfg_av = solver_config.artif_viscosity; - Cfg_MHD cfg_mhd = solver_config.mhd_config; + Cfg_AV cfg_av = solver_config.artif_viscosity; + Cfg_MHD cfg_mhd = solver_config.mhd_config; + DustConfig cfg_dust = solver_config.dust_config; if (Constant *v = std::get_if(&cfg_av.config)) { update_derivs_constantAV(*v); @@ -49,6 +50,11 @@ void shammodels::sph::modules::UpdateDerivs::update_derivs() { } else { shambase::throw_unimplemented(); } + + if (cfg_dust.has_Sj_field()) { + // we can do it separately because the backreaction is done only through the pressure + update_derivs_dust_monofluid_tvi_Sj(cfg_dust); + } } template class SPHKernel> @@ -1207,6 +1213,181 @@ void shammodels::sph::modules::UpdateDerivs::update_derivs_MHD( }); } +template class SPHKernel> +void shammodels::sph::modules::UpdateDerivs::update_derivs_dust_monofluid_tvi_Sj( + DustConfig cfg) { + StackEntry stack_loc{}; + + using namespace shamrock; + using namespace shamrock::patch; + + PatchDataLayerLayout &pdl = scheduler().pdl(); + + const u32 ixyz = pdl.get_field_idx("xyz"); + const u32 ivxyz = pdl.get_field_idx("vxyz"); + const u32 iaxyz = pdl.get_field_idx("axyz"); + const u32 ihpart = pdl.get_field_idx("hpart"); + const u32 is_j = pdl.get_field_idx("s_j"); + const u32 ids_j_dt = pdl.get_field_idx("ds_j_dt"); + + shamrock::patch::PatchDataLayerLayout &ghost_layout + = shambase::get_check_ref(storage.ghost_layout.get()); + u32 ihpart_interf = ghost_layout.get_field_idx("hpart"); + u32 ivxyz_interf = ghost_layout.get_field_idx("vxyz"); + u32 iomega_interf = ghost_layout.get_field_idx("omega"); + u32 is_j_interf = ghost_layout.get_field_idx("s_j"); + + u32 ndust = cfg.get_dust_nvar(); + + auto &merged_xyzh = storage.merged_xyzh.get(); + shamrock::solvergraph::Field &omega = shambase::get_check_ref(storage.omega); + shambase::DistributedData &mpdats = storage.merged_patchdata_ghost.get(); + + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + PatchDataLayer &mpdat = mpdats.get(cur_p.id_patch); + sham::DeviceBuffer &buf_xyz + = merged_xyzh.get(cur_p.id_patch).template get_field_buf_ref(0); + sham::DeviceBuffer &buf_axyz = pdat.get_field_buf_ref(iaxyz); + sham::DeviceBuffer &buf_vxyz = mpdat.get_field_buf_ref(ivxyz_interf); + sham::DeviceBuffer &buf_hpart = mpdat.get_field_buf_ref(ihpart_interf); + sham::DeviceBuffer &buf_omega = mpdat.get_field_buf_ref(iomega_interf); + sham::DeviceBuffer &buf_pressure + = shambase::get_check_ref(storage.pressure).get_field(cur_p.id_patch).get_buf(); + sham::DeviceBuffer &buf_alpha_AV + = storage.alpha_av_ghost.get().get(cur_p.id_patch).get_buf(); + sham::DeviceBuffer &buf_cs + = shambase::get_check_ref(storage.soundspeed).get_field(cur_p.id_patch).get_buf(); + + sham::DeviceBuffer &buf_s_j = mpdat.get_field_buf_ref(is_j_interf); + sham::DeviceBuffer &buf_ds_j_dt = pdat.get_field_buf_ref(ids_j_dt); + + sycl::range range_npart{pdat.get_obj_cnt()}; + + tree::ObjectCache &pcache + = shambase::get_check_ref(storage.neigh_cache).get_cache(cur_p.id_patch); + + ///////////////////////////////////////////// + + sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue(); + sham::EventList depends_list; + + auto xyz = buf_xyz.get_read_access(depends_list); + auto axyz = buf_axyz.get_write_access(depends_list); + auto vxyz = buf_vxyz.get_read_access(depends_list); + auto hpart = buf_hpart.get_read_access(depends_list); + auto omega = buf_omega.get_read_access(depends_list); + auto pressure = buf_pressure.get_read_access(depends_list); + auto alpha_AV = buf_alpha_AV.get_read_access(depends_list); + auto cs = buf_cs.get_read_access(depends_list); + auto ploop_ptrs = pcache.get_read_access(depends_list); + + auto s_j = buf_s_j.get_read_access(depends_list); + auto ds_j_dt = buf_ds_j_dt.get_write_access(depends_list); + + auto e = q.submit(depends_list, [&](sycl::handler &cgh) { + const Tscal pmass = solver_config.gpart_mass; + + tree::ObjectCacheIterator particle_looper(ploop_ptrs); + + // sycl::accessor hmax_tree{tree_field_hmax, cgh, sycl::read_only}; + + // sycl::stream out {4096,1024,cgh}; + + constexpr Tscal Rker2 = Kernel::Rkern * Kernel::Rkern; + + shambase::parallel_for( + cgh, pdat.get_obj_cnt() * ndust, "compute force CD10 AV", [=](u64 gid) { + u32 id_a = (u32) gid / ndust; + u32 jdust = (u32) gid % ndust; + + using namespace shamrock::sph; + + Tscal term1 = 0; + Tscal term2 = 0; + + Tscal h_a = hpart[id_a]; + Tvec xyz_a = xyz[id_a]; + Tvec vxyz_a = vxyz[id_a]; + Tscal P_a = pressure[id_a]; + Tscal cs_a = cs[id_a]; + Tscal omega_a = omega[id_a]; + const Tscal alpha_a = alpha_AV[id_a]; + Tscal rho_a = rho_h(pmass, h_a, Kernel::hfactd); + Tscal rho_a_sq = rho_a * rho_a; + Tscal rho_a_inv = 1. / rho_a; + + Tscal s_j_a = s_j[gid]; + Tscal Ttilde_sj_a = 0; // TODO + + // f32 P_a = cs * cs * rho_a; + + Tscal omega_a_rho_a_inv = 1 / (omega_a * rho_a); + + Tvec force_pressure{0, 0, 0}; + Tscal tmpdU_pressure = 0; + + particle_looper.for_each_object(id_a, [&](u32 id_b) { + // compute only omega_a + Tvec dr = xyz_a - xyz[id_b]; + Tscal rab2 = sycl::dot(dr, dr); + Tscal h_b = hpart[id_b]; + + if (rab2 > h_a * h_a * Rker2 && rab2 > h_b * h_b * Rker2) { + return; + } + + Tvec vxyz_b = vxyz[id_b]; + Tscal P_b = pressure[id_b]; + Tscal omega_b = omega[id_b]; + const Tscal alpha_b = alpha_AV[id_b]; + Tscal cs_b = cs[id_b]; + + Tscal s_j_b = s_j[id_b * ndust + jdust]; + Tscal Ttilde_sj_b = 0; // TODO + + Tscal rab = sycl::sqrt(rab2); + + Tscal rho_b = rho_h(pmass, h_b, Kernel::hfactd); + + Tscal Fab_a = Kernel::dW_3d(rab, h_a); + Tscal Fab_b = Kernel::dW_3d(rab, h_b); + + Tvec v_ab = vxyz_a - vxyz_b; + + Tvec r_ab_unit = dr * sham::inv_sat_positive(rab); + + Tscal F_ab_bar = (Fab_a + Fab_b) / 2; + Tscal delta_P = P_a - P_b; + Tscal Ts_weighted = (Ttilde_sj_a / rho_a + Ttilde_sj_b / rho_b); + + term1 += (pmass * s_j_b / rho_b) * Ts_weighted * delta_P * F_ab_bar; + term2 += pmass * sham::dot(v_ab, r_ab_unit * Fab_a); + }); + + // eq 51, Hutchison 2018 + Tscal ds_j_dt_a = Tscal{-0.5} * term1 + (s_j_a / (2 * rho_a * omega_a)) * term2; + + ds_j_dt[gid] = ds_j_dt_a; + }); + }); + + buf_xyz.complete_event_state(e); + buf_axyz.complete_event_state(e); + buf_vxyz.complete_event_state(e); + buf_hpart.complete_event_state(e); + buf_omega.complete_event_state(e); + buf_pressure.complete_event_state(e); + buf_alpha_AV.complete_event_state(e); + buf_cs.complete_event_state(e); + buf_s_j.complete_event_state(e); + buf_ds_j_dt.complete_event_state(e); + + sham::EventList resulting_events; + resulting_events.add_event(e); + pcache.complete_event_state(resulting_events); + }); +} + using namespace shammath; template class shammodels::sph::modules::UpdateDerivs; template class shammodels::sph::modules::UpdateDerivs; From 8acd1504a3a7cc66bc2935f7847ae82a00605b7f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Mon, 9 Mar 2026 11:16:24 +0100 Subject: [PATCH 04/24] fix compile --- src/shammodels/sph/src/modules/UpdateDerivs.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/shammodels/sph/src/modules/UpdateDerivs.cpp b/src/shammodels/sph/src/modules/UpdateDerivs.cpp index 83b7c639e7..52942d3955 100644 --- a/src/shammodels/sph/src/modules/UpdateDerivs.cpp +++ b/src/shammodels/sph/src/modules/UpdateDerivs.cpp @@ -1065,7 +1065,7 @@ void shammodels::sph::modules::UpdateDerivs::update_derivs_dust using namespace shamrock; using namespace shamrock::patch; - PatchDataLayerLayout &pdl = scheduler().pdl(); + PatchDataLayerLayout &pdl = scheduler().pdl_old(); const u32 ixyz = pdl.get_field_idx("xyz"); const u32 ivxyz = pdl.get_field_idx("vxyz"); From ccdc245c60d813d6d59758a9122d438d9b763c37 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Mon, 9 Mar 2026 11:52:06 +0100 Subject: [PATCH 05/24] move to standalone node --- src/shammodels/sph/CMakeLists.txt | 1 + .../modules/NodeUpdateDerivsMonofluidTVI.hpp | 76 ++++++++ .../modules/NodeUpdateDerivsMonofluidTVI.cpp | 162 ++++++++++++++++++ 3 files changed, 239 insertions(+) create mode 100644 src/shammodels/sph/include/shammodels/sph/modules/NodeUpdateDerivsMonofluidTVI.hpp create mode 100644 src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp diff --git a/src/shammodels/sph/CMakeLists.txt b/src/shammodels/sph/CMakeLists.txt index 295ff990ca..795e093011 100644 --- a/src/shammodels/sph/CMakeLists.txt +++ b/src/shammodels/sph/CMakeLists.txt @@ -51,6 +51,7 @@ set(Sources src/modules/LoopSmoothingLengthIter.cpp src/modules/ComputeNeighStats.cpp src/modules/NodeUpdateDerivsVaryingAlphaAV.cpp + src/modules/NodeUpdateDerivsMonofluidTVI.cpp src/modules/setup/GeneratorMCDisc.cpp src/modules/setup/ModifierApplyDiscWarp.cpp src/modules/setup/ModifierApplyCustomWarp.cpp diff --git a/src/shammodels/sph/include/shammodels/sph/modules/NodeUpdateDerivsMonofluidTVI.hpp b/src/shammodels/sph/include/shammodels/sph/modules/NodeUpdateDerivsMonofluidTVI.hpp new file mode 100644 index 0000000000..21fd02865f --- /dev/null +++ b/src/shammodels/sph/include/shammodels/sph/modules/NodeUpdateDerivsMonofluidTVI.hpp @@ -0,0 +1,76 @@ +// -------------------------------------------------------// +// +// 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 NodeUpdateDerivsMonofluidTVI.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief + * + */ + +#include "shambackends/vec.hpp" +#include "shammodels/sph/solvergraph/NeighCache.hpp" +#include "shamrock/solvergraph/IFieldSpan.hpp" +#include "shamrock/solvergraph/INode.hpp" +#include "shamrock/solvergraph/Indexes.hpp" +#include "shamrock/solvergraph/ScalarEdge.hpp" + +#define NODE_UPDATE_DERIVS_MONOFLUID_TVI_EDGES(X_RO, X_RW) \ + /* scalars */ \ + X_RO(shamrock::solvergraph::ScalarEdge, gpart_mass) \ + X_RO(shamrock::solvergraph::ScalarEdge, alpha_u) \ + X_RO(shamrock::solvergraph::ScalarEdge, beta_AV) \ + \ + /* counts */ \ + X_RO(shamrock::solvergraph::Indexes, part_counts) \ + X_RO(shamrock::solvergraph::Indexes, part_counts_with_ghost) \ + \ + /* fields */ \ + X_RO(shamrock::solvergraph::IFieldSpan, xyz) \ + X_RO(shamrock::solvergraph::IFieldSpan, hpart) \ + X_RO(shamrock::solvergraph::IFieldSpan, vxyz) \ + X_RO(shamrock::solvergraph::IFieldSpan, omega) \ + X_RO(shamrock::solvergraph::IFieldSpan, pressure) \ + X_RO(shamrock::solvergraph::IFieldSpan, cs) \ + X_RO(shamrock::solvergraph::IFieldSpan, s_j) \ + X_RO(shamrock::solvergraph::IFieldSpan, Ttilde_sj) \ + \ + /* neigh */ \ + X_RO(shammodels::sph::solvergraph::NeighCache, neigh_cache) \ + \ + /* outputs */ \ + X_RW(shamrock::solvergraph::IFieldSpan, ds_j_dt) + +namespace shammodels::sph::modules { + + template class SPHKernel> + class NodeUpdateDerivsMonofluidTVI : public shamrock::solvergraph::INode { + + using Tscal = shambase::VecComponent; + + static constexpr Tscal kernel_radius = SPHKernel::Rkern; + + u32 ndust; + + public: + NodeUpdateDerivsMonofluidTVI(u32 ndust) : ndust(ndust) {} + + EXPAND_NODE_EDGES(NODE_UPDATE_DERIVS_MONOFLUID_TVI_EDGES) + + void _impl_evaluate_internal(); + + inline virtual std::string _impl_get_label() const { + return "NodeUpdateDerivsMonofluidTVI"; + }; + + inline virtual std::string _impl_get_tex() const { return "TODO"; }; + }; +} // namespace shammodels::sph::modules diff --git a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp new file mode 100644 index 0000000000..92bdafdf0a --- /dev/null +++ b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp @@ -0,0 +1,162 @@ +// -------------------------------------------------------// +// +// 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 NodeUpdateDerivsMonofluidTVI.cpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief + * + */ + +#include "shammodels/sph/modules/NodeUpdateDerivsMonofluidTVI.hpp" +#include "shambackends/kernel_call_distrib.hpp" +#include "shammath/sphkernels.hpp" +#include "shammodels/sph/math/density.hpp" +#include "shamrock/patch/PatchDataField.hpp" + +template class SPHKernel> +struct KernelUpdateDerivsMonofluidTVI { + using Tscal = shambase::VecComponent; + using Kernel = SPHKernel; + static constexpr Tscal hfactd = Kernel::hfactd; + static constexpr Tscal Rkern = Kernel::Rkern; + static constexpr Tscal Rker2 = Rkern * Rkern; + + Tscal pmass; + u32 ndust; + + inline void operator()( + u32 thread_id, + const Tvec *__restrict xyz, + const Tscal *__restrict hpart, + const Tvec *__restrict vxyz, + const Tscal *__restrict omega, + const Tscal *__restrict pressure, + const Tscal *__restrict s_j, + const Tscal *__restrict Ttilde_sj, + shamrock::tree::ObjectCache::ptrs_read ploop_ptrs, + Tscal *__restrict ds_j_dt) const { + + u32 id_a = (u32) thread_id / ndust; + u32 jdust = (u32) thread_id % ndust; + + Tscal h_a = hpart[id_a]; + Tvec xyz_a = xyz[id_a]; + Tvec vxyz_a = vxyz[id_a]; + Tscal P_a = pressure[id_a]; + Tscal omega_a = omega[id_a]; + Tscal s_j_a = s_j[thread_id]; + Tscal Ttilde_sj_a = Ttilde_sj[thread_id]; + + using namespace shamrock::sph; + Tscal rho_a = rho_h(pmass, h_a, Kernel::hfactd); + Tscal rho_a_sq = rho_a * rho_a; + Tscal rho_a_inv = 1. / rho_a; + Tscal omega_a_rho_a_inv = 1 / (omega_a * rho_a); + + Tscal term1 = 0; + Tscal term2 = 0; + + shamrock::tree::ObjectCacheIterator particle_looper(ploop_ptrs); + particle_looper.for_each_object(id_a, [&](u32 id_b) { + Tvec dr = xyz_a - xyz[id_b]; + Tscal rab2 = sycl::dot(dr, dr); + Tscal h_b = hpart[id_b]; + + if (rab2 > h_a * h_a * Rker2 && rab2 > h_b * h_b * Rker2) { + return; + } + + Tvec vxyz_b = vxyz[id_b]; + Tscal P_b = pressure[id_b]; + Tscal omega_b = omega[id_b]; + + Tscal s_j_b = s_j[id_b * ndust + jdust]; + Tscal Ttilde_sj_b = Ttilde_sj[id_b * ndust + jdust]; + + Tscal rab = sycl::sqrt(rab2); + + Tscal rho_b = rho_h(pmass, h_b, Kernel::hfactd); + + Tscal Fab_a = Kernel::dW_3d(rab, h_a); + Tscal Fab_b = Kernel::dW_3d(rab, h_b); + + Tvec v_ab = vxyz_a - vxyz_b; + + Tvec r_ab_unit = dr * sham::inv_sat_positive(rab); + + Tscal F_ab_bar = (Fab_a + Fab_b) / 2; + Tscal delta_P = P_a - P_b; + Tscal Ts_weighted = (Ttilde_sj_a / rho_a + Ttilde_sj_b / rho_b); + + term1 += (pmass * s_j_b / rho_b) * Ts_weighted * delta_P * F_ab_bar; + term2 += pmass * sham::dot(v_ab, r_ab_unit * Fab_a); + }); + + // eq 51, Hutchison 2018 + Tscal ds_j_dt_a = Tscal{-0.5} * term1 + (s_j_a / (2 * rho_a * omega_a)) * term2; + + ds_j_dt[thread_id] = ds_j_dt_a; + } +}; + +template class SPHKernel> +void shammodels::sph::modules::NodeUpdateDerivsMonofluidTVI:: + _impl_evaluate_internal() { + + __shamrock_stack_entry(); + + auto edges = get_edges(); + + auto &part_counts_with_ghost = edges.part_counts_with_ghost.indexes; + auto &part_counts = edges.part_counts.indexes; + + // check that all input edges have the particles with ghosts zones + edges.xyz.check_sizes(part_counts_with_ghost); + edges.hpart.check_sizes(part_counts_with_ghost); + edges.vxyz.check_sizes(part_counts_with_ghost); + edges.omega.check_sizes(part_counts_with_ghost); + edges.pressure.check_sizes(part_counts_with_ghost); + edges.s_j.check_sizes(part_counts_with_ghost); + edges.Ttilde_sj.check_sizes(part_counts_with_ghost); + + // ensure that the output edges are of size part_counts (output without ghosts zones) + edges.ds_j_dt.ensure_sizes(part_counts); + + const Tscal pmass = edges.gpart_mass.value; + const Tscal alpha_u = edges.alpha_u.value; + const Tscal beta_AV = edges.beta_AV.value; + + using ComputeKernel = KernelUpdateDerivsMonofluidTVI; + + // call the kernel for each patches with part_counts.get(id_patch) threads of patch id_patch + sham::distributed_data_kernel_call( + shamsys::instance::get_compute_scheduler_ptr(), + sham::DDMultiRef{ + edges.xyz.get_spans(), + edges.hpart.get_spans(), + edges.vxyz.get_spans(), + edges.omega.get_spans(), + edges.pressure.get_spans(), + edges.s_j.get_spans(), + edges.Ttilde_sj.get_spans(), + edges.neigh_cache}, + sham::DDMultiRef{edges.ds_j_dt.get_spans()}, + part_counts, + ComputeKernel{pmass, ndust}); +} + +using namespace shammath; +template class shammodels::sph::modules::NodeUpdateDerivsMonofluidTVI; +template class shammodels::sph::modules::NodeUpdateDerivsMonofluidTVI; +template class shammodels::sph::modules::NodeUpdateDerivsMonofluidTVI; + +template class shammodels::sph::modules::NodeUpdateDerivsMonofluidTVI; +template class shammodels::sph::modules::NodeUpdateDerivsMonofluidTVI; +template class shammodels::sph::modules::NodeUpdateDerivsMonofluidTVI; From 215653548a12f09c514e0a3bfec7100eaac71d24 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Mon, 9 Mar 2026 13:36:21 +0100 Subject: [PATCH 06/24] finish move --- .../include/shammodels/sph/SolverConfig.hpp | 2 +- .../modules/NodeUpdateDerivsMonofluidTVI.hpp | 3 - src/shammodels/sph/src/Solver.cpp | 36 +++ src/shammodels/sph/src/SolverConfig.cpp | 4 +- .../modules/NodeUpdateDerivsMonofluidTVI.cpp | 4 +- .../sph/src/modules/UpdateDerivs.cpp | 205 ++++++------------ 6 files changed, 109 insertions(+), 145 deletions(-) diff --git a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp index 66105bd0f6..9d30cffb88 100644 --- a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp +++ b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp @@ -135,7 +135,7 @@ namespace shammodels::sph { return bool(std::get_if(¤t_mode)); } - inline bool has_Sj_field() { + inline bool has_s_j_field() { return is_monofluid_tvi(); // S_j = sqrt(\rho \epsilon_j) } diff --git a/src/shammodels/sph/include/shammodels/sph/modules/NodeUpdateDerivsMonofluidTVI.hpp b/src/shammodels/sph/include/shammodels/sph/modules/NodeUpdateDerivsMonofluidTVI.hpp index 21fd02865f..0fbc75410c 100644 --- a/src/shammodels/sph/include/shammodels/sph/modules/NodeUpdateDerivsMonofluidTVI.hpp +++ b/src/shammodels/sph/include/shammodels/sph/modules/NodeUpdateDerivsMonofluidTVI.hpp @@ -26,8 +26,6 @@ #define NODE_UPDATE_DERIVS_MONOFLUID_TVI_EDGES(X_RO, X_RW) \ /* scalars */ \ X_RO(shamrock::solvergraph::ScalarEdge, gpart_mass) \ - X_RO(shamrock::solvergraph::ScalarEdge, alpha_u) \ - X_RO(shamrock::solvergraph::ScalarEdge, beta_AV) \ \ /* counts */ \ X_RO(shamrock::solvergraph::Indexes, part_counts) \ @@ -39,7 +37,6 @@ X_RO(shamrock::solvergraph::IFieldSpan, vxyz) \ X_RO(shamrock::solvergraph::IFieldSpan, omega) \ X_RO(shamrock::solvergraph::IFieldSpan, pressure) \ - X_RO(shamrock::solvergraph::IFieldSpan, cs) \ X_RO(shamrock::solvergraph::IFieldSpan, s_j) \ X_RO(shamrock::solvergraph::IFieldSpan, Ttilde_sj) \ \ diff --git a/src/shammodels/sph/src/Solver.cpp b/src/shammodels/sph/src/Solver.cpp index b3e8130c04..84f637bbb3 100644 --- a/src/shammodels/sph/src/Solver.cpp +++ b/src/shammodels/sph/src/Solver.cpp @@ -110,6 +110,7 @@ void shammodels::sph::Solver::init_solver_graph() { bool has_psi_field = solver_config.has_field_psi_on_ch(); bool has_epsilon_field = solver_config.dust_config.has_epsilon_field(); bool has_deltav_field = solver_config.dust_config.has_deltav_field(); + bool has_s_j_field = solver_config.dust_config.has_s_j_field(); using namespace shamrock::solvergraph; @@ -146,6 +147,10 @@ void shammodels::sph::Solver::init_solver_graph() { solver_graph.register_edge("deltav", FieldRefs("deltav", "\\Delta v")); solver_graph.register_edge("dtdeltav", FieldRefs("dtdeltav", "d\\Delta v")); } + if (has_s_j_field) { + solver_graph.register_edge("s_j", FieldRefs("s_j", "S_j")); + solver_graph.register_edge("ds_j_dt", FieldRefs("ds_j_dt", "dS_j/dt")); + } { auto set_gpart_mass = solver_graph.register_node( @@ -328,6 +333,25 @@ void shammodels::sph::Solver::init_solver_graph() { attach_field_sequence.push_back(attach_dtdeltav); } + if (has_s_j_field) { + auto attach_s_j + = solver_graph.register_node("attach_s_j", GetFieldRefFromLayer(pdl, "s_j")); + shambase::get_check_ref(attach_s_j) + .set_edges( + solver_graph.get_edge_ptr("scheduler_patchdata"), + solver_graph.get_edge_ptr>("s_j")); + attach_field_sequence.push_back(attach_s_j); + } + + if (has_s_j_field) { + auto attach_ds_j_dt = solver_graph.register_node( + "attach_ds_j_dt", GetFieldRefFromLayer(pdl, "ds_j_dt")); + shambase::get_check_ref(attach_ds_j_dt) + .set_edges( + solver_graph.get_edge_ptr("scheduler_patchdata"), + solver_graph.get_edge_ptr>("ds_j_dt")); + attach_field_sequence.push_back(attach_ds_j_dt); + } solver_graph.register_node( "attach fields to scheduler", OperationSequence("attach fields", std::move(attach_field_sequence))); @@ -414,6 +438,18 @@ void shammodels::sph::Solver::init_solver_graph() { half_step_sequence.push_back(half_step_deltav); } + if (has_s_j_field) { + auto half_step_s_j = solver_graph.register_node( + prefix + "_s_j", shammodels::common::modules::ForwardEuler{}); + shambase::get_check_ref(half_step_s_j) + .set_edges( + solver_graph.get_edge_ptr>("dt_half"), + solver_graph.get_edge_ptr>("ds_j_dt"), + solver_graph.get_edge_ptr>("part_counts"), + solver_graph.get_edge_ptr>("s_j")); + half_step_sequence.push_back(half_step_s_j); + } + return OperationSequence("half step", std::move(half_step_sequence)); }; diff --git a/src/shammodels/sph/src/SolverConfig.cpp b/src/shammodels/sph/src/SolverConfig.cpp index 85ba030a2e..b5d3c8d957 100644 --- a/src/shammodels/sph/src/SolverConfig.cpp +++ b/src/shammodels/sph/src/SolverConfig.cpp @@ -112,7 +112,7 @@ namespace shammodels::sph { pdl.add_field("dt_part", 1); } - if (dust_config.has_Sj_field()) { + if (dust_config.has_s_j_field()) { u32 ndust = dust_config.get_dust_nvar(); pdl.add_field("s_j", ndust); pdl.add_field("ds_j_dt", ndust); @@ -158,7 +158,7 @@ namespace shammodels::sph { ghost_layout.add_field("deltav", ndust); } - if (dust_config.has_Sj_field()) { + if (dust_config.has_s_j_field()) { u32 ndust = dust_config.get_dust_nvar(); ghost_layout.add_field("s_j", ndust); } diff --git a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp index 92bdafdf0a..3c3bd72ed7 100644 --- a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp +++ b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp @@ -129,9 +129,7 @@ void shammodels::sph::modules::NodeUpdateDerivsMonofluidTVI:: // ensure that the output edges are of size part_counts (output without ghosts zones) edges.ds_j_dt.ensure_sizes(part_counts); - const Tscal pmass = edges.gpart_mass.value; - const Tscal alpha_u = edges.alpha_u.value; - const Tscal beta_AV = edges.beta_AV.value; + const Tscal pmass = edges.gpart_mass.value; using ComputeKernel = KernelUpdateDerivsMonofluidTVI; diff --git a/src/shammodels/sph/src/modules/UpdateDerivs.cpp b/src/shammodels/sph/src/modules/UpdateDerivs.cpp index 52942d3955..b94495f2ee 100644 --- a/src/shammodels/sph/src/modules/UpdateDerivs.cpp +++ b/src/shammodels/sph/src/modules/UpdateDerivs.cpp @@ -26,6 +26,7 @@ #include "shammodels/sph/math/forces.hpp" #include "shammodels/sph/math/mhd.hpp" #include "shammodels/sph/math/q_ab.hpp" +#include "shammodels/sph/modules/NodeUpdateDerivsMonofluidTVI.hpp" #include "shammodels/sph/modules/NodeUpdateDerivsVaryingAlphaAV.hpp" #include "shammodels/sph/modules/UpdateDerivs.hpp" #include "shamphys/mhd.hpp" @@ -62,7 +63,7 @@ void shammodels::sph::modules::UpdateDerivs::update_derivs() { shambase::throw_unimplemented(); } - if (cfg_dust.has_Sj_field()) { + if (cfg_dust.has_s_j_field()) { // we can do it separately because the backreaction is done only through the pressure update_derivs_dust_monofluid_tvi_Sj(cfg_dust); } @@ -1087,149 +1088,81 @@ void shammodels::sph::modules::UpdateDerivs::update_derivs_dust shamrock::solvergraph::Field &omega = shambase::get_check_ref(storage.omega); shambase::DistributedData &mpdats = storage.merged_patchdata_ghost.get(); - scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { - PatchDataLayer &mpdat = mpdats.get(cur_p.id_patch); - sham::DeviceBuffer &buf_xyz - = merged_xyzh.get(cur_p.id_patch).template get_field_buf_ref(0); - sham::DeviceBuffer &buf_axyz = pdat.get_field_buf_ref(iaxyz); - sham::DeviceBuffer &buf_vxyz = mpdat.get_field_buf_ref(ivxyz_interf); - sham::DeviceBuffer &buf_hpart = mpdat.get_field_buf_ref(ihpart_interf); - sham::DeviceBuffer &buf_omega = mpdat.get_field_buf_ref(iomega_interf); - sham::DeviceBuffer &buf_pressure - = shambase::get_check_ref(storage.pressure).get_field(cur_p.id_patch).get_buf(); - sham::DeviceBuffer &buf_alpha_AV - = storage.alpha_av_ghost.get().get(cur_p.id_patch).get_buf(); - sham::DeviceBuffer &buf_cs - = shambase::get_check_ref(storage.soundspeed).get_field(cur_p.id_patch).get_buf(); - - sham::DeviceBuffer &buf_s_j = mpdat.get_field_buf_ref(is_j_interf); - sham::DeviceBuffer &buf_ds_j_dt = pdat.get_field_buf_ref(ids_j_dt); - - sycl::range range_npart{pdat.get_obj_cnt()}; - - tree::ObjectCache &pcache - = shambase::get_check_ref(storage.neigh_cache).get_cache(cur_p.id_patch); - - ///////////////////////////////////////////// - - sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue(); - sham::EventList depends_list; - - auto xyz = buf_xyz.get_read_access(depends_list); - auto axyz = buf_axyz.get_write_access(depends_list); - auto vxyz = buf_vxyz.get_read_access(depends_list); - auto hpart = buf_hpart.get_read_access(depends_list); - auto omega = buf_omega.get_read_access(depends_list); - auto pressure = buf_pressure.get_read_access(depends_list); - auto alpha_AV = buf_alpha_AV.get_read_access(depends_list); - auto cs = buf_cs.get_read_access(depends_list); - auto ploop_ptrs = pcache.get_read_access(depends_list); - - auto s_j = buf_s_j.get_read_access(depends_list); - auto ds_j_dt = buf_ds_j_dt.get_write_access(depends_list); - - auto e = q.submit(depends_list, [&](sycl::handler &cgh) { - const Tscal pmass = solver_config.gpart_mass; - - tree::ObjectCacheIterator particle_looper(ploop_ptrs); - - // sycl::accessor hmax_tree{tree_field_hmax, cgh, sycl::read_only}; - - // sycl::stream out {4096,1024,cgh}; - - constexpr Tscal Rker2 = Kernel::Rkern * Kernel::Rkern; - - shambase::parallel_for( - cgh, pdat.get_obj_cnt() * ndust, "compute force CD10 AV", [=](u64 gid) { - u32 id_a = (u32) gid / ndust; - u32 jdust = (u32) gid % ndust; - - using namespace shamrock::sph; - - Tscal term1 = 0; - Tscal term2 = 0; - - Tscal h_a = hpart[id_a]; - Tvec xyz_a = xyz[id_a]; - Tvec vxyz_a = vxyz[id_a]; - Tscal P_a = pressure[id_a]; - Tscal cs_a = cs[id_a]; - Tscal omega_a = omega[id_a]; - const Tscal alpha_a = alpha_AV[id_a]; - Tscal rho_a = rho_h(pmass, h_a, Kernel::hfactd); - Tscal rho_a_sq = rho_a * rho_a; - Tscal rho_a_inv = 1. / rho_a; - - Tscal s_j_a = s_j[gid]; - Tscal Ttilde_sj_a = 0; // TODO - - // f32 P_a = cs * cs * rho_a; - - Tscal omega_a_rho_a_inv = 1 / (omega_a * rho_a); - - Tvec force_pressure{0, 0, 0}; - Tscal tmpdU_pressure = 0; - - particle_looper.for_each_object(id_a, [&](u32 id_b) { - // compute only omega_a - Tvec dr = xyz_a - xyz[id_b]; - Tscal rab2 = sycl::dot(dr, dr); - Tscal h_b = hpart[id_b]; - - if (rab2 > h_a * h_a * Rker2 && rab2 > h_b * h_b * Rker2) { - return; - } - - Tvec vxyz_b = vxyz[id_b]; - Tscal P_b = pressure[id_b]; - Tscal omega_b = omega[id_b]; - const Tscal alpha_b = alpha_AV[id_b]; - Tscal cs_b = cs[id_b]; - - Tscal s_j_b = s_j[id_b * ndust + jdust]; - Tscal Ttilde_sj_b = 0; // TODO - - Tscal rab = sycl::sqrt(rab2); - - Tscal rho_b = rho_h(pmass, h_b, Kernel::hfactd); - - Tscal Fab_a = Kernel::dW_3d(rab, h_a); - Tscal Fab_b = Kernel::dW_3d(rab, h_b); + auto &part_counts = storage.part_counts; + auto &part_counts_with_ghost = storage.part_counts_with_ghost; + auto &xyz_refs = storage.positions_with_ghosts; + auto &pressure_field = storage.pressure; - Tvec v_ab = vxyz_a - vxyz_b; + std::shared_ptr> vxyz_refs + = std::make_shared>("vxyz", "v"); + { + shambase::get_check_ref(vxyz_refs).set_refs( + mpdats.map>>( + [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) { + return std::ref(mpdat.get_field(ivxyz_interf)); + })); + } - Tvec r_ab_unit = dr * sham::inv_sat_positive(rab); + std::shared_ptr> hpart_refs + = std::make_shared>("hpart", "h"); + { // if was just reset before this call + shambase::get_check_ref(hpart_refs) + .set_refs(mpdats.map>>( + [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) { + return std::ref(mpdat.get_field(ihpart_interf)); + })); + } - Tscal F_ab_bar = (Fab_a + Fab_b) / 2; - Tscal delta_P = P_a - P_b; - Tscal Ts_weighted = (Ttilde_sj_a / rho_a + Ttilde_sj_b / rho_b); + std::shared_ptr> omega_refs + = std::make_shared>("omega", "omega"); + { + shambase::get_check_ref(omega_refs) + .set_refs(mpdats.map>>( + [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) { + return std::ref(mpdat.get_field(iomega_interf)); + })); + } - term1 += (pmass * s_j_b / rho_b) * Ts_weighted * delta_P * F_ab_bar; - term2 += pmass * sham::dot(v_ab, r_ab_unit * Fab_a); - }); + // s_j_interf + std::shared_ptr> s_j_refs + = std::make_shared>("s_j", "s_j"); + { + shambase::get_check_ref(s_j_refs).set_refs( + mpdats.map>>( + [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) { + return std::ref(mpdat.get_field(is_j_interf)); + })); + } - // eq 51, Hutchison 2018 - Tscal ds_j_dt_a = Tscal{-0.5} * term1 + (s_j_a / (2 * rho_a * omega_a)) * term2; + std::shared_ptr> Ttilde_sj_field + = std::make_shared>(ndust, "Ttilde_sj", "Ttilde_sj"); + // TODO: fill the field with the correct values - ds_j_dt[gid] = ds_j_dt_a; - }); - }); + shamrock::solvergraph::SolverGraph &solver_graph = storage.solver_graph; + auto ds_j_dt_refs + = solver_graph.get_edge_ptr>("ds_j_dt"); - buf_xyz.complete_event_state(e); - buf_axyz.complete_event_state(e); - buf_vxyz.complete_event_state(e); - buf_hpart.complete_event_state(e); - buf_omega.complete_event_state(e); - buf_pressure.complete_event_state(e); - buf_alpha_AV.complete_event_state(e); - buf_cs.complete_event_state(e); - buf_s_j.complete_event_state(e); - buf_ds_j_dt.complete_event_state(e); + auto gpart_mass + = solver_graph.get_edge_ptr>("gpart_mass"); - sham::EventList resulting_events; - resulting_events.add_event(e); - pcache.complete_event_state(resulting_events); - }); + std::shared_ptr> node + = std::make_shared>(ndust); + { + node->set_edges( + gpart_mass, + part_counts, + part_counts_with_ghost, + xyz_refs, + hpart_refs, + vxyz_refs, + omega_refs, + pressure_field, + s_j_refs, + Ttilde_sj_field, + storage.neigh_cache, + ds_j_dt_refs); + } + node->evaluate(); } using namespace shammath; From 190f083d8335614c9f9ede5587db90983c93ad16 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Mon, 9 Mar 2026 13:50:50 +0100 Subject: [PATCH 07/24] cleaner --- .../modules/NodeUpdateDerivsMonofluidTVI.cpp | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp index 3c3bd72ed7..2d64601e34 100644 --- a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp +++ b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp @@ -33,6 +33,7 @@ struct KernelUpdateDerivsMonofluidTVI { inline void operator()( u32 thread_id, + // input const Tvec *__restrict xyz, const Tscal *__restrict hpart, const Tvec *__restrict vxyz, @@ -41,10 +42,11 @@ struct KernelUpdateDerivsMonofluidTVI { const Tscal *__restrict s_j, const Tscal *__restrict Ttilde_sj, shamrock::tree::ObjectCache::ptrs_read ploop_ptrs, + // output Tscal *__restrict ds_j_dt) const { - u32 id_a = (u32) thread_id / ndust; - u32 jdust = (u32) thread_id % ndust; + u32 id_a = thread_id / ndust; + u32 jdust = thread_id % ndust; Tscal h_a = hpart[id_a]; Tvec xyz_a = xyz[id_a]; @@ -73,10 +75,9 @@ struct KernelUpdateDerivsMonofluidTVI { return; } - Tvec vxyz_b = vxyz[id_b]; - Tscal P_b = pressure[id_b]; - Tscal omega_b = omega[id_b]; - + Tvec vxyz_b = vxyz[id_b]; + Tscal P_b = pressure[id_b]; + Tscal omega_b = omega[id_b]; Tscal s_j_b = s_j[id_b * ndust + jdust]; Tscal Ttilde_sj_b = Ttilde_sj[id_b * ndust + jdust]; @@ -100,9 +101,7 @@ struct KernelUpdateDerivsMonofluidTVI { }); // eq 51, Hutchison 2018 - Tscal ds_j_dt_a = Tscal{-0.5} * term1 + (s_j_a / (2 * rho_a * omega_a)) * term2; - - ds_j_dt[thread_id] = ds_j_dt_a; + ds_j_dt[thread_id] = Tscal{-0.5} * term1 + (s_j_a / (2 * rho_a * omega_a)) * term2; } }; From 7f07ef926117513f1e706dbf31c772c5d4fbd1c8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Mon, 9 Mar 2026 15:45:48 +0100 Subject: [PATCH 08/24] ts compute --- .../sph/modules/ComputeDustTtilde.hpp | 129 ++++++++++++++++++ .../sph/src/modules/UpdateDerivs.cpp | 14 +- 2 files changed, 142 insertions(+), 1 deletion(-) create mode 100644 src/shammodels/sph/include/shammodels/sph/modules/ComputeDustTtilde.hpp diff --git a/src/shammodels/sph/include/shammodels/sph/modules/ComputeDustTtilde.hpp b/src/shammodels/sph/include/shammodels/sph/modules/ComputeDustTtilde.hpp new file mode 100644 index 0000000000..086ad66507 --- /dev/null +++ b/src/shammodels/sph/include/shammodels/sph/modules/ComputeDustTtilde.hpp @@ -0,0 +1,129 @@ +// -------------------------------------------------------// +// +// 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 ComputeDustTtilde.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief + * + */ + +#include "shambackends/kernel_call_distrib.hpp" +#include "shambackends/vec.hpp" +#include "shammodels/sph/math/density.hpp" +#include "shammodels/sph/solvergraph/NeighCache.hpp" +#include "shamrock/solvergraph/IFieldSpan.hpp" +#include "shamrock/solvergraph/INode.hpp" +#include "shamrock/solvergraph/Indexes.hpp" +#include "shamrock/solvergraph/ScalarEdge.hpp" + +#define NODE_COMPUTE_DUST_TTILDE_EDGES(X_RO, X_RW) \ + /* scalars */ \ + X_RO(shamrock::solvergraph::ScalarEdge, gpart_mass) \ + \ + /* counts */ \ + X_RO(shamrock::solvergraph::Indexes, part_counts) \ + \ + /* fields */ \ + X_RO(shamrock::solvergraph::IFieldSpan, hpart) \ + X_RO(shamrock::solvergraph::IFieldSpan, s_j) \ + X_RO(shamrock::solvergraph::IFieldSpan, t_j) \ + \ + /* outputs */ \ + X_RW(shamrock::solvergraph::IFieldSpan, Ttilde_sj) + +namespace shammodels::sph::modules { + + template class SPHKernel> + class ComputeDustTtilde : public shamrock::solvergraph::INode { + + using Tscal = shambase::VecComponent; + using Kernel = SPHKernel; + + static constexpr Tscal kernel_radius = SPHKernel::Rkern; + + u32 ndust; + + public: + ComputeDustTtilde(u32 ndust) : ndust(ndust) {} + + EXPAND_NODE_EDGES(NODE_COMPUTE_DUST_TTILDE_EDGES) + + inline void _impl_evaluate_internal() { + + __shamrock_stack_entry(); + + auto edges = get_edges(); + + auto &part_counts = edges.part_counts.indexes; + + // check that all input edges have the particles with ghosts zones + edges.hpart.check_sizes(part_counts); + edges.s_j.check_sizes(part_counts); + edges.t_j.check_sizes(part_counts); + + // ensure that the output edges are of size part_counts (output without ghosts zones) + edges.Ttilde_sj.ensure_sizes(part_counts); + + const Tscal pmass = edges.gpart_mass.value; + + // call the kernel for each patches with part_counts.get(id_patch) threads of patch + // id_patch + sham::distributed_data_kernel_call( + shamsys::instance::get_compute_scheduler_ptr(), + sham::DDMultiRef{ + edges.hpart.get_spans(), edges.s_j.get_spans(), edges.t_j.get_spans()}, + sham::DDMultiRef{edges.Ttilde_sj.get_spans()}, + part_counts, + [pmass, ndust = ndust]( + u32 thread_id, + const Tscal *__restrict hpart, + const Tscal *__restrict s_j, + const Tscal *__restrict t_j, + Tscal *__restrict Ttilde_sj) { + u32 id_a = thread_id / ndust; + u32 jdust = thread_id % ndust; + + Tscal h_a = hpart[id_a]; + Tscal sj_a = s_j[thread_id]; + Tscal tj_a = t_j[thread_id]; + + using namespace shamrock::sph; + Tscal rho_a = rho_h(pmass, hpart, Kernel::hfactd); + + auto epsilon = [&](Tscal sj) { + return sj * sj / rho_a; + }; + /* + * Hutchison 2018 eq 15 + * \tilde{T}_{sj} = \epsilon_j t_j - \sum_{k} \epsilon_k^2 t_k + */ + + Tscal eps_j_a = epsilon(sj_a); + Tscal Ttilde_sj_a = eps_j_a * tj_a; + + for (u32 k = 0; k < ndust; k++) { + Tscal sk_a = s_j[id_a * ndust + k]; + Tscal tk_a = t_j[id_a * ndust + k]; + + Tscal eps_k_a = epsilon(sk_a); + Ttilde_sj_a -= eps_k_a * eps_k_a * tk_a; + } + + Ttilde_sj[thread_id] = Ttilde_sj_a; + }); + } + + inline virtual std::string _impl_get_label() const { return "ComputeDustTtilde"; }; + + inline virtual std::string _impl_get_tex() const { return "TODO"; }; + }; +} // namespace shammodels::sph::modules diff --git a/src/shammodels/sph/src/modules/UpdateDerivs.cpp b/src/shammodels/sph/src/modules/UpdateDerivs.cpp index b94495f2ee..5f3dd82897 100644 --- a/src/shammodels/sph/src/modules/UpdateDerivs.cpp +++ b/src/shammodels/sph/src/modules/UpdateDerivs.cpp @@ -26,6 +26,7 @@ #include "shammodels/sph/math/forces.hpp" #include "shammodels/sph/math/mhd.hpp" #include "shammodels/sph/math/q_ab.hpp" +#include "shammodels/sph/modules/ComputeDustTtilde.hpp" #include "shammodels/sph/modules/NodeUpdateDerivsMonofluidTVI.hpp" #include "shammodels/sph/modules/NodeUpdateDerivsVaryingAlphaAV.hpp" #include "shammodels/sph/modules/UpdateDerivs.hpp" @@ -1134,9 +1135,12 @@ void shammodels::sph::modules::UpdateDerivs::update_derivs_dust })); } + std::shared_ptr> t_j_field + = std::make_shared>(ndust, "t_j", "t_j"); + // TODO: fill the field with the correct values + std::shared_ptr> Ttilde_sj_field = std::make_shared>(ndust, "Ttilde_sj", "Ttilde_sj"); - // TODO: fill the field with the correct values shamrock::solvergraph::SolverGraph &solver_graph = storage.solver_graph; auto ds_j_dt_refs @@ -1145,6 +1149,14 @@ void shammodels::sph::modules::UpdateDerivs::update_derivs_dust auto gpart_mass = solver_graph.get_edge_ptr>("gpart_mass"); + std::shared_ptr> node_tj + = std::make_shared>(ndust); + { + node_tj->set_edges( + gpart_mass, part_counts_with_ghost, hpart_refs, s_j_refs, t_j_field, Ttilde_sj_field); + } + node_tj->evaluate(); + std::shared_ptr> node = std::make_shared>(ndust); { From 7c75bae4a877cca3de24ce924efcbe9a4c2ca325 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Mon, 9 Mar 2026 16:07:36 +0100 Subject: [PATCH 09/24] yeah boy --- .../sph/modules/ComputeDustTtilde.hpp | 4 +- .../modules/SetDustStoppingTimeConstant.hpp | 94 +++++++++++++++++++ .../sph/src/modules/UpdateDerivs.cpp | 15 ++- 3 files changed, 110 insertions(+), 3 deletions(-) create mode 100644 src/shammodels/sph/include/shammodels/sph/modules/SetDustStoppingTimeConstant.hpp diff --git a/src/shammodels/sph/include/shammodels/sph/modules/ComputeDustTtilde.hpp b/src/shammodels/sph/include/shammodels/sph/modules/ComputeDustTtilde.hpp index 086ad66507..122796ccec 100644 --- a/src/shammodels/sph/include/shammodels/sph/modules/ComputeDustTtilde.hpp +++ b/src/shammodels/sph/include/shammodels/sph/modules/ComputeDustTtilde.hpp @@ -19,11 +19,11 @@ #include "shambackends/kernel_call_distrib.hpp" #include "shambackends/vec.hpp" #include "shammodels/sph/math/density.hpp" -#include "shammodels/sph/solvergraph/NeighCache.hpp" #include "shamrock/solvergraph/IFieldSpan.hpp" #include "shamrock/solvergraph/INode.hpp" #include "shamrock/solvergraph/Indexes.hpp" #include "shamrock/solvergraph/ScalarEdge.hpp" +#include "shamsys/NodeInstance.hpp" #define NODE_COMPUTE_DUST_TTILDE_EDGES(X_RO, X_RW) \ /* scalars */ \ @@ -97,7 +97,7 @@ namespace shammodels::sph::modules { Tscal tj_a = t_j[thread_id]; using namespace shamrock::sph; - Tscal rho_a = rho_h(pmass, hpart, Kernel::hfactd); + Tscal rho_a = rho_h(pmass, h_a, Kernel::hfactd); auto epsilon = [&](Tscal sj) { return sj * sj / rho_a; diff --git a/src/shammodels/sph/include/shammodels/sph/modules/SetDustStoppingTimeConstant.hpp b/src/shammodels/sph/include/shammodels/sph/modules/SetDustStoppingTimeConstant.hpp new file mode 100644 index 0000000000..fb19373968 --- /dev/null +++ b/src/shammodels/sph/include/shammodels/sph/modules/SetDustStoppingTimeConstant.hpp @@ -0,0 +1,94 @@ +// -------------------------------------------------------// +// +// 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 SetDustStoppingTimeConstant.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief + * + */ + +#include "shambackends/DeviceBuffer.hpp" +#include "shambackends/kernel_call_distrib.hpp" +#include "shambackends/vec.hpp" +#include "shammodels/sph/math/density.hpp" +#include "shamrock/solvergraph/IFieldSpan.hpp" +#include "shamrock/solvergraph/INode.hpp" +#include "shamrock/solvergraph/Indexes.hpp" +#include "shamrock/solvergraph/ScalarEdge.hpp" +#include "shamsys/NodeInstance.hpp" +#include + +#define NODE_SET_DUST_STOPPING_TIME_CTE_EDGES(X_RO, X_RW) \ + /* scalars */ \ + X_RO(shamrock::solvergraph::ScalarEdge>, t_j_0) \ + \ + /* counts */ \ + X_RO(shamrock::solvergraph::Indexes, part_counts) \ + \ + /* fields */ \ + X_RW(shamrock::solvergraph::IFieldSpan, t_j) + +namespace shammodels::sph::modules { + + template class SPHKernel> + class SetDustStoppingTimeConstant : public shamrock::solvergraph::INode { + + using Tscal = shambase::VecComponent; + using Kernel = SPHKernel; + + static constexpr Tscal kernel_radius = SPHKernel::Rkern; + + u32 ndust; + + public: + SetDustStoppingTimeConstant(u32 ndust) : ndust(ndust) {} + + EXPAND_NODE_EDGES(NODE_SET_DUST_STOPPING_TIME_CTE_EDGES) + + inline void _impl_evaluate_internal() { + + __shamrock_stack_entry(); + + auto edges = get_edges(); + + auto &part_counts = edges.part_counts.indexes; + const std::vector &inputs_tj = edges.t_j_0.value; + + // ensure that the output edges are of size part_counts + edges.t_j.ensure_sizes(part_counts); + + sham::DeviceBuffer t_j_0(ndust, shamsys::instance::get_compute_scheduler_ptr()); + t_j_0.copy_from_stdvec(inputs_tj); + + auto &q = shamsys::instance::get_compute_scheduler().get_queue(); + + part_counts.for_each([&](u64 id, u32 count) { + // call the kernel for each patches with part_counts.get(id_patch) threads of patch + // id_patch + sham::kernel_call( + q, + sham::MultiRef{t_j_0}, + sham::MultiRef{edges.t_j.get_spans().get(id)}, + part_counts.get(id), + [ndust + = ndust](u32 thread_id, const Tscal *__restrict t_j_0, Tscal *__restrict t_j) { + u32 jdust = thread_id % ndust; + t_j[thread_id] = t_j_0[jdust]; + }); + }); + } + + inline virtual std::string _impl_get_label() const { return "ComputeDustTtilde"; }; + + inline virtual std::string _impl_get_tex() const { return "TODO"; }; + }; +} // namespace shammodels::sph::modules diff --git a/src/shammodels/sph/src/modules/UpdateDerivs.cpp b/src/shammodels/sph/src/modules/UpdateDerivs.cpp index 5f3dd82897..8660ebcadb 100644 --- a/src/shammodels/sph/src/modules/UpdateDerivs.cpp +++ b/src/shammodels/sph/src/modules/UpdateDerivs.cpp @@ -29,6 +29,7 @@ #include "shammodels/sph/modules/ComputeDustTtilde.hpp" #include "shammodels/sph/modules/NodeUpdateDerivsMonofluidTVI.hpp" #include "shammodels/sph/modules/NodeUpdateDerivsVaryingAlphaAV.hpp" +#include "shammodels/sph/modules/SetDustStoppingTimeConstant.hpp" #include "shammodels/sph/modules/UpdateDerivs.hpp" #include "shamphys/mhd.hpp" #include "shamrock/patch/PatchDataFieldSpan.hpp" @@ -36,6 +37,8 @@ #include "shamrock/solvergraph/IFieldSpan.hpp" #include "shamrock/solvergraph/Indexes.hpp" #include "shamrock/solvergraph/ScalarEdge.hpp" +#include +#include template class SPHKernel> void shammodels::sph::modules::UpdateDerivs::update_derivs() { @@ -1137,7 +1140,17 @@ void shammodels::sph::modules::UpdateDerivs::update_derivs_dust std::shared_ptr> t_j_field = std::make_shared>(ndust, "t_j", "t_j"); - // TODO: fill the field with the correct values + + std::shared_ptr>> input_t_j + = std::make_shared>>("", ""); + input_t_j->value = std::vector(ndust, Tscal{0}); + + std::shared_ptr> node_set_tj + = std::make_shared>(ndust); + { + node_set_tj->set_edges(input_t_j, part_counts_with_ghost, t_j_field); + } + node_set_tj->evaluate(); std::shared_ptr> Ttilde_sj_field = std::make_shared>(ndust, "Ttilde_sj", "Ttilde_sj"); From 7d7383503157dc1281ef4c1665d910c5e899a3ea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Mon, 9 Mar 2026 16:19:51 +0100 Subject: [PATCH 10/24] connect to config --- src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp | 6 ++++++ src/shammodels/sph/src/modules/UpdateDerivs.cpp | 2 +- src/shammodels/sph/src/pySPHModel.cpp | 5 +++++ 3 files changed, 12 insertions(+), 1 deletion(-) diff --git a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp index 9d30cffb88..b56b001798 100644 --- a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp +++ b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp @@ -121,6 +121,8 @@ namespace shammodels::sph { u32 ndust; }; + std::vector stopping_times; + /// Variant type to store the EOS configuration using Variant = std::variant; @@ -176,6 +178,10 @@ namespace shammodels::sph { "SPH::config", "Dust config != None is work in progress, use it at your own risk")); } + + if (get_dust_nvar() != stopping_times.size()) { + throw ""; + } } } }; diff --git a/src/shammodels/sph/src/modules/UpdateDerivs.cpp b/src/shammodels/sph/src/modules/UpdateDerivs.cpp index 8660ebcadb..6f69865fb0 100644 --- a/src/shammodels/sph/src/modules/UpdateDerivs.cpp +++ b/src/shammodels/sph/src/modules/UpdateDerivs.cpp @@ -1143,7 +1143,7 @@ void shammodels::sph::modules::UpdateDerivs::update_derivs_dust std::shared_ptr>> input_t_j = std::make_shared>>("", ""); - input_t_j->value = std::vector(ndust, Tscal{0}); + input_t_j->value = cfg.stopping_times; std::shared_ptr> node_set_tj = std::make_shared>(ndust); diff --git a/src/shammodels/sph/src/pySPHModel.cpp b/src/shammodels/sph/src/pySPHModel.cpp index fb9cbd9831..d5877253dc 100644 --- a/src/shammodels/sph/src/pySPHModel.cpp +++ b/src/shammodels/sph/src/pySPHModel.cpp @@ -260,6 +260,11 @@ void add_instance(py::module &m, std::string name_config, std::string name_model }, py::kw_only(), py::arg("ndust")) + .def( + "set_dust_stopping_times", + [](TConfig &self, std::vector ts) { + self.dust_config.stopping_times = ts; + }) .def("add_ext_force_point_mass", &TConfig::add_ext_force_point_mass) .def( "add_ext_force_lense_thirring", From 149695ff929aec24a82f5283b3ed9c7371f0a8c7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Mon, 9 Mar 2026 16:37:57 +0100 Subject: [PATCH 11/24] prepare tests --- examples/sph/run_dustybox.py | 82 ++++++++++++++++++++++++++++++++++++ 1 file changed, 82 insertions(+) create mode 100644 examples/sph/run_dustybox.py diff --git a/examples/sph/run_dustybox.py b/examples/sph/run_dustybox.py new file mode 100644 index 0000000000..d2437874c6 --- /dev/null +++ b/examples/sph/run_dustybox.py @@ -0,0 +1,82 @@ +""" +Testing Sod tube with SPH +========================= + +CI test for Sod tube with SPH +""" + +import matplotlib.pyplot as plt + +import shamrock +import numpy as np + +gamma = 1.4 +rho_g = 1 +P_g = 1 + +u_g = P_g / ((gamma - 1) * rho_g) + + +bmin = (-0.6, -0.6, -0.6) +bmax = (0.6, 0.6, 0.6) + +N_target = 1e4 +scheduler_split_val = int(2e7) +scheduler_merge_val = int(1) + +xm, ym, zm = bmin +xM, yM, zM = bmax +vol_b = (xM - xm) * (yM - ym) * (zM - zm) + +part_vol = vol_b / N_target + +# lattice volume +HCP_PACKING_DENSITY = 0.74 +part_vol_lattice = HCP_PACKING_DENSITY * part_vol + +dr = (part_vol_lattice / ((4.0 / 3.0) * np.pi)) ** (1.0 / 3.0) + +pmass = -1 + +ctx = shamrock.Context() +ctx.pdata_layout_new() + +model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M6") + +cfg = model.gen_default_config() +# cfg.set_artif_viscosity_Constant(alpha_u = 1, alpha_AV = 1, beta_AV = 2) +# cfg.set_artif_viscosity_VaryingMM97(alpha_min = 0.1,alpha_max = 1,sigma_decay = 0.1, alpha_u = 1, beta_AV = 2) +cfg.set_artif_viscosity_VaryingCD10( + alpha_min=0.0, alpha_max=1, sigma_decay=0.1, alpha_u=1, beta_AV=2 +) +cfg.set_boundary_periodic() +cfg.set_eos_adiabatic(gamma) +cfg.print_status() +model.set_solver_config(cfg) + +model.init_scheduler(int(1e8), 1) + + +bmin, bmax = model.get_ideal_hcp_box(dr, bmin, bmax) +xm, ym, zm = bmin +xM, yM, zM = bmax + +model.resize_simulation_box(bmin, bmax) + +setup = model.get_setup() +gen = setup.make_generator_lattice_hcp(dr, bmin, bmax) +setup.apply_setup(gen, insert_step=scheduler_split_val) + + + +vol_b = (xM - xm) * (yM - ym) * (zM - zm) +totmass = rho_g * vol_b +print("Total mass :", totmass) + +pmass = model.total_mass_to_part_mass(totmass) +model.set_particle_mass(pmass) +print("Current part mass :", pmass) + + +model.set_cfl_cour(0.1) +model.set_cfl_force(0.1) From f890fb1705b351c6eca4ee62ba0020f659ea9350 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Mon, 9 Mar 2026 23:45:22 +0100 Subject: [PATCH 12/24] prepare what will follow --- examples/sph/run_dustybox.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/examples/sph/run_dustybox.py b/examples/sph/run_dustybox.py index d2437874c6..3e801929c6 100644 --- a/examples/sph/run_dustybox.py +++ b/examples/sph/run_dustybox.py @@ -80,3 +80,18 @@ model.set_cfl_cour(0.1) model.set_cfl_force(0.1) + +dic = ctx.collect_data() + +print(dic["xyz"].shape) + +x = dic["xyz"][:,0] +y = dic["xyz"][:,1] +z = dic["xyz"][:,2] + +r = np.sqrt(x*x + y*y + z*z) + +hpart = dic["hpart"] + +plt.scatter(r,hpart) +plt.show() From b548d50475290b7dfc9c16e89eeba33ec7d0364a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Tue, 10 Mar 2026 14:55:42 +0100 Subject: [PATCH 13/24] run but unstable --- examples/sph/run_dusty.py | 166 ++++++++++++++++++ examples/sph/run_dustybox.py | 97 ---------- .../sph/include/shammodels/sph/Model.hpp | 38 ++-- .../include/shammodels/sph/SolverConfig.hpp | 3 +- src/shammodels/sph/src/Solver.cpp | 12 ++ src/shammodels/sph/src/modules/ComputeEos.cpp | 63 +++++++ src/shammodels/sph/src/pySPHModel.cpp | 27 ++- 7 files changed, 283 insertions(+), 123 deletions(-) create mode 100644 examples/sph/run_dusty.py delete mode 100644 examples/sph/run_dustybox.py diff --git a/examples/sph/run_dusty.py b/examples/sph/run_dusty.py new file mode 100644 index 0000000000..4110cc8892 --- /dev/null +++ b/examples/sph/run_dusty.py @@ -0,0 +1,166 @@ +""" +Testing Sod tube with SPH +========================= + +CI test for Sod tube with SPH +""" + +import matplotlib.pyplot as plt + +import shamrock + +shamrock.enable_experimental_features() +import numpy as np + +gamma = 1.4 +rho = 1 + + +def func_rho_t(r): + return rho + + +def func_rho_d_j(r, idust): + r_ = np.sqrt(r[0] ** 2 + r[1] ** 2 + r[2] ** 2) + return (0.1 / ndust) * max(1 - (r_ / rc) ** 2, 0) + + +def func_rho_g(r): + return rho - sum([func_rho_d_j(r, i) for i in range(ndust)]) + + +def func_s_j(r, idust): + rho_t = func_rho_t(r) + rho_d_j = [func_rho_d_j(r, i) for i in range(ndust)] + eps_j = rho_d_j[idust] / rho_t + return np.sqrt(rho_t * eps_j) + + +cs_g = 1 + + +def uint_g(r): + rho_g = func_rho_g(r) + P = rho_g * cs_g * cs_g / gamma + return P / ((gamma - 1) * rho_g) + + +ndust = 3 +rc = 0.3 +stopping_times = [1e-4, 1e-3, 1e-2] + + +bmin = (-0.6, -0.6, -0.6) +bmax = (0.6, 0.6, 0.6) + +N_target = 1e4 +scheduler_split_val = int(2e7) +scheduler_merge_val = int(1) + +xm, ym, zm = bmin +xM, yM, zM = bmax +vol_b = (xM - xm) * (yM - ym) * (zM - zm) + +part_vol = vol_b / N_target + +# lattice volume +HCP_PACKING_DENSITY = 0.74 +part_vol_lattice = HCP_PACKING_DENSITY * part_vol + +dr = (part_vol_lattice / ((4.0 / 3.0) * np.pi)) ** (1.0 / 3.0) + +pmass = -1 + +ctx = shamrock.Context() +ctx.pdata_layout_new() + +model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M6") + +cfg = model.gen_default_config() +# cfg.set_artif_viscosity_Constant(alpha_u = 1, alpha_AV = 1, beta_AV = 2) +# cfg.set_artif_viscosity_VaryingMM97(alpha_min = 0.1,alpha_max = 1,sigma_decay = 0.1, alpha_u = 1, beta_AV = 2) +cfg.set_artif_viscosity_VaryingCD10( + alpha_min=0.0, alpha_max=1, sigma_decay=0.1, alpha_u=1, beta_AV=2 +) +cfg.set_dust_mode_monofluid_tvi(ndust) +cfg.set_dust_stopping_times(stopping_times) +cfg.set_boundary_periodic() +cfg.set_eos_adiabatic(gamma) +cfg.print_status() +model.set_solver_config(cfg) + +model.init_scheduler(int(1e8), 1) + + +bmin, bmax = model.get_ideal_hcp_box(dr, bmin, bmax) +xm, ym, zm = bmin +xM, yM, zM = bmax + +model.resize_simulation_box(bmin, bmax) + +setup = model.get_setup() +gen = setup.make_generator_lattice_hcp(dr, bmin, bmax) +setup.apply_setup(gen, insert_step=scheduler_split_val) + +for i in range(ndust): + + def func_s(r): + return func_s_j(r, i) + + model.set_field_value_lambda_f64("s_j", func_s, i) + +model.set_field_value_lambda_f64("uint", uint_g) + +vol_b = (xM - xm) * (yM - ym) * (zM - zm) +totmass = rho * vol_b +print("Total mass :", totmass) + +pmass = model.total_mass_to_part_mass(totmass) +model.set_particle_mass(pmass) +print("Current part mass :", pmass) + + +model.set_cfl_cour(0.1) +model.set_cfl_force(0.1) + +model.timestep() + +for j in range(10): + if j > 0: + for i in range(10): + model.timestep() + + dic = ctx.collect_data() + print(dic["s_j"]) + + print(dic["xyz"].shape) + + x = dic["xyz"][:, 0] + y = dic["xyz"][:, 1] + z = dic["xyz"][:, 2] + s_j = dic["s_j"].reshape(-1, ndust) + ds_j_dt = dic["ds_j_dt"].reshape(-1, ndust) + cs = dic["soundspeed"] + + print(s_j) + + r = np.sqrt(x * x + y * y + z * z) + + hpart = dic["hpart"] + rho = pmass * (model.get_hfact() / np.array(hpart)) ** 3 + + fig, axs = plt.subplots(nrows=1, ncols=4, figsize=(15, 5)) + axs[0].scatter(r, rho, label="rho") + axs[0].legend() + for i in range(ndust): + axs[1].scatter(r, s_j[:, i] ** 2 / rho, label=f"eps_j_{i}") + axs[1].legend() + + axs[2].scatter(r, cs, label="soundspeed") + axs[2].legend() + + axs[3].scatter(r, ds_j_dt[:, 0], label=f"ds_j_dt_{0}") + axs[3].scatter(r, ds_j_dt[:, 1], label=f"ds_j_dt_{1}") + axs[3].scatter(r, ds_j_dt[:, 2], label=f"ds_j_dt_{2}") + axs[3].legend() + plt.show() diff --git a/examples/sph/run_dustybox.py b/examples/sph/run_dustybox.py deleted file mode 100644 index 3e801929c6..0000000000 --- a/examples/sph/run_dustybox.py +++ /dev/null @@ -1,97 +0,0 @@ -""" -Testing Sod tube with SPH -========================= - -CI test for Sod tube with SPH -""" - -import matplotlib.pyplot as plt - -import shamrock -import numpy as np - -gamma = 1.4 -rho_g = 1 -P_g = 1 - -u_g = P_g / ((gamma - 1) * rho_g) - - -bmin = (-0.6, -0.6, -0.6) -bmax = (0.6, 0.6, 0.6) - -N_target = 1e4 -scheduler_split_val = int(2e7) -scheduler_merge_val = int(1) - -xm, ym, zm = bmin -xM, yM, zM = bmax -vol_b = (xM - xm) * (yM - ym) * (zM - zm) - -part_vol = vol_b / N_target - -# lattice volume -HCP_PACKING_DENSITY = 0.74 -part_vol_lattice = HCP_PACKING_DENSITY * part_vol - -dr = (part_vol_lattice / ((4.0 / 3.0) * np.pi)) ** (1.0 / 3.0) - -pmass = -1 - -ctx = shamrock.Context() -ctx.pdata_layout_new() - -model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M6") - -cfg = model.gen_default_config() -# cfg.set_artif_viscosity_Constant(alpha_u = 1, alpha_AV = 1, beta_AV = 2) -# cfg.set_artif_viscosity_VaryingMM97(alpha_min = 0.1,alpha_max = 1,sigma_decay = 0.1, alpha_u = 1, beta_AV = 2) -cfg.set_artif_viscosity_VaryingCD10( - alpha_min=0.0, alpha_max=1, sigma_decay=0.1, alpha_u=1, beta_AV=2 -) -cfg.set_boundary_periodic() -cfg.set_eos_adiabatic(gamma) -cfg.print_status() -model.set_solver_config(cfg) - -model.init_scheduler(int(1e8), 1) - - -bmin, bmax = model.get_ideal_hcp_box(dr, bmin, bmax) -xm, ym, zm = bmin -xM, yM, zM = bmax - -model.resize_simulation_box(bmin, bmax) - -setup = model.get_setup() -gen = setup.make_generator_lattice_hcp(dr, bmin, bmax) -setup.apply_setup(gen, insert_step=scheduler_split_val) - - - -vol_b = (xM - xm) * (yM - ym) * (zM - zm) -totmass = rho_g * vol_b -print("Total mass :", totmass) - -pmass = model.total_mass_to_part_mass(totmass) -model.set_particle_mass(pmass) -print("Current part mass :", pmass) - - -model.set_cfl_cour(0.1) -model.set_cfl_force(0.1) - -dic = ctx.collect_data() - -print(dic["xyz"].shape) - -x = dic["xyz"][:,0] -y = dic["xyz"][:,1] -z = dic["xyz"][:,2] - -r = np.sqrt(x*x + y*y + z*z) - -hpart = dic["hpart"] - -plt.scatter(r,hpart) -plt.show() diff --git a/src/shammodels/sph/include/shammodels/sph/Model.hpp b/src/shammodels/sph/include/shammodels/sph/Model.hpp index 97f514113a..a105473a6a 100644 --- a/src/shammodels/sph/include/shammodels/sph/Model.hpp +++ b/src/shammodels/sph/include/shammodels/sph/Model.hpp @@ -177,38 +177,32 @@ namespace shammodels::sph { template inline void set_field_value_lambda( - std::string field_name, const std::function pos_to_val) { + std::string field_name, const std::function pos_to_val, const i32 offset) { StackEntry stack_loc{}; + PatchScheduler &sched = shambase::get_check_ref(ctx.sched); - sched.patch_data.for_each_patchdata( - [&](u64 patch_id, shamrock::patch::PatchDataLayer &pdat) { - PatchDataField &xyz - = pdat.template get_field(sched.pdl_old().get_field_idx("xyz")); - PatchDataField &f - = pdat.template get_field(sched.pdl_old().get_field_idx(field_name)); + u32 ixyz = sched.pdl_old().get_field_idx("xyz"); + u32 ifield = sched.pdl_old().get_field_idx(field_name); - if (f.get_nvar() != 1) { - shambase::throw_unimplemented(); - } + sched.patch_data.for_each_patchdata( + [&](u64 patch_id, shamrock::patch::PatchDataLayer &pdat) { + PatchDataField &xyz = pdat.template get_field(ixyz); + PatchDataField &f = pdat.template get_field(ifield); - { - auto &buf = f.get_buf(); - auto acc = buf.copy_to_stdvec(); + auto f_nvar = f.get_nvar(); - auto &buf_xyz = xyz.get_buf(); - auto acc_xyz = buf_xyz.copy_to_stdvec(); + auto acc = f.get_buf().copy_to_stdvec(); + auto acc_xyz = xyz.get_buf().copy_to_stdvec(); - for (u32 i = 0; i < f.get_obj_cnt(); i++) { - Tvec r = acc_xyz[i]; + for (u32 i = 0; i < pdat.get_obj_cnt(); i++) { + acc[i * f_nvar + offset] = pos_to_val(acc_xyz[i]); + } - acc[i] = pos_to_val(r); - } + logger::raw_ln(acc); - buf.copy_from_stdvec(acc); - buf_xyz.copy_from_stdvec(acc_xyz); - } + f.get_buf().copy_from_stdvec(acc); }); } diff --git a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp index b56b001798..8c5e73338f 100644 --- a/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp +++ b/src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp @@ -142,8 +142,7 @@ namespace shammodels::sph { } inline bool has_epsilon_field() { - return bool(std::get_if(¤t_mode)) - || bool(std::get_if(¤t_mode)); + return bool(std::get_if(¤t_mode)); } inline bool has_deltav_field() { diff --git a/src/shammodels/sph/src/Solver.cpp b/src/shammodels/sph/src/Solver.cpp index 84f637bbb3..b1eed9277a 100644 --- a/src/shammodels/sph/src/Solver.cpp +++ b/src/shammodels/sph/src/Solver.cpp @@ -1299,6 +1299,7 @@ void shammodels::sph::Solver::communicate_merge_ghosts_fields() { bool has_curlB_field = solver_config.has_field_curlB(); bool has_epsilon_field = solver_config.dust_config.has_epsilon_field(); bool has_deltav_field = solver_config.dust_config.has_deltav_field(); + bool has_s_j_field = solver_config.dust_config.has_s_j_field(); PatchDataLayerLayout &pdl = scheduler().pdl_old(); const u32 ixyz = pdl.get_field_idx("xyz"); @@ -1329,6 +1330,7 @@ void shammodels::sph::Solver::communicate_merge_ghosts_fields() { const u32 iepsilon = (has_epsilon_field) ? pdl.get_field_idx("epsilon") : 0; const u32 ideltav = (has_deltav_field) ? pdl.get_field_idx("deltav") : 0; + const u32 is_j = (has_s_j_field) ? pdl.get_field_idx("s_j") : 0; auto &ghost_layout_ptr = storage.ghost_layout; shamrock::patch::PatchDataLayerLayout &ghost_layout = shambase::get_check_ref(ghost_layout_ptr); @@ -1350,6 +1352,7 @@ void shammodels::sph::Solver::communicate_merge_ghosts_fields() { const u32 iepsilon_interf = (has_epsilon_field) ? ghost_layout.get_field_idx("epsilon") : 0; const u32 ideltav_interf = (has_deltav_field) ? ghost_layout.get_field_idx("deltav") : 0; + const u32 is_j_interf = (has_s_j_field) ? ghost_layout.get_field_idx("s_j") : 0; using InterfaceBuildInfos = typename sph::BasicSPHGhostHandler::InterfaceBuildInfos; @@ -1422,6 +1425,11 @@ void shammodels::sph::Solver::communicate_merge_ghosts_fields() { sender_patch.get_field(ideltav).append_subset_to( buf_idx, cnt, pdat.get_field(ideltav_interf)); } + + if (has_s_j_field) { + sender_patch.get_field(is_j).append_subset_to( + buf_idx, cnt, pdat.get_field(is_j_interf)); + } }); ghost_handle.template modify_interface_native( @@ -1498,6 +1506,10 @@ void shammodels::sph::Solver::communicate_merge_ghosts_fields() { pdat_new.get_field(ideltav_interf).insert(pdat.get_field(ideltav)); } + if (has_s_j_field) { + pdat_new.get_field(is_j_interf).insert(pdat.get_field(is_j)); + } + pdat_new.check_field_obj_cnt_match(); return pdat_new; diff --git a/src/shammodels/sph/src/modules/ComputeEos.cpp b/src/shammodels/sph/src/modules/ComputeEos.cpp index ded3f9b76a..26dcb34fbc 100644 --- a/src/shammodels/sph/src/modules/ComputeEos.cpp +++ b/src/shammodels/sph/src/modules/ComputeEos.cpp @@ -19,6 +19,7 @@ #include "shambase/DistributedData.hpp" #include "shambase/exception.hpp" #include "shambackends/kernel_call.hpp" +#include "shamcomm/logs.hpp" #include "shammath/sphkernels.hpp" #include "shammodels/sph/math/density.hpp" #include "shammodels/sph/modules/ComputeEos.hpp" @@ -89,6 +90,50 @@ struct RhoGetterMonofluid { void complete_event_state(sycl::event e) { buf_h.complete_event_state(e); } }; +template +struct RhoGetterSJ { + sham::DeviceBuffer &buf_h; + sham::DeviceBuffer &buf_s_j; + u32 nvar_dust; + Tscal pmass; + Tscal hfact; + + struct accessed { + const Tscal *h; + const Tscal *buf_s_j; + u32 nvar_dust; + Tscal pmass; + Tscal hfact; + + Tscal operator()(u32 i) const { + using namespace shamrock::sph; + Tscal rho = rho_h(pmass, h[i], hfact); + + Tscal epsilon_sum = 0; + for (u32 j = 0; j < nvar_dust; j++) { + Tscal s_j = buf_s_j[i * nvar_dust + j]; + epsilon_sum += s_j * s_j / rho; + } + + logger::raw_ln(i, rho, epsilon_sum, rho * (1 - epsilon_sum)); + + return rho * (1 - epsilon_sum); + } + }; + + accessed get_read_access(sham::EventList &depends_list) { + auto h = buf_h.get_read_access(depends_list); + auto s_j = buf_s_j.get_read_access(depends_list); + + return accessed{h, s_j, nvar_dust, pmass, hfact}; + } + + void complete_event_state(sycl::event e) { + buf_h.complete_event_state(e); + buf_s_j.complete_event_state(e); + } +}; + template class SPHKernel> template void shammodels::sph::modules::ComputeEos::compute_eos_internal( @@ -561,6 +606,9 @@ void shammodels::sph::modules::ComputeEos::compute_eos() { shambase::get_check_ref(storage.pressure).ensure_sizes(counts_with_ghosts); shambase::get_check_ref(storage.soundspeed).ensure_sizes(counts_with_ghosts); + logger::raw_ln( + solver_config.dust_config.has_epsilon_field(), solver_config.dust_config.has_s_j_field()); + if (solver_config.dust_config.has_epsilon_field()) { u32 iepsilon_interf = ghost_layout.get_field_idx("epsilon"); @@ -574,6 +622,21 @@ void shammodels::sph::modules::ComputeEos::compute_eos() { gpart_mass, Kernel::hfactd}; }); + } else if (solver_config.dust_config.has_s_j_field()) { + + u32 is_j_interf = ghost_layout.get_field_idx("s_j"); + u32 nvar_dust = solver_config.dust_config.get_dust_nvar(); + + compute_eos_internal([&](PatchDataLayer &mpdat) { + return RhoGetterSJ{ + mpdat.get_field_buf_ref(ihpart_interf), + mpdat.get_field_buf_ref(is_j_interf), + nvar_dust, + gpart_mass, + Kernel::hfactd}; + }); + + logger::info_ln("ComputeEos", "SJ field has been computed"); } else { compute_eos_internal([&](PatchDataLayer &mpdat) { return RhoGetterBase{ diff --git a/src/shammodels/sph/src/pySPHModel.cpp b/src/shammodels/sph/src/pySPHModel.cpp index d5877253dc..2154aa4860 100644 --- a/src/shammodels/sph/src/pySPHModel.cpp +++ b/src/shammodels/sph/src/pySPHModel.cpp @@ -38,6 +38,7 @@ #include #include #include +#include template class SPHKernel> void add_instance(py::module &m, std::string name_config, std::string name_model) { @@ -746,8 +747,30 @@ void add_instance(py::module &m, std::string name_config, std::string name_model "unknown field type"); } }) - .def("set_field_value_lambda_f64_3", &T::template set_field_value_lambda) - .def("set_field_value_lambda_f64", &T::template set_field_value_lambda) + .def( + "set_field_value_lambda_f64", + [](T &self, + std::string field_name, + const std::function pos_to_val, + const i32 offset) { + return self.template set_field_value_lambda( + std::move(field_name), pos_to_val, offset); + }, + py::arg("field_name"), + py::arg("pos_to_val"), + py::arg("offset") = 0) + .def( + "set_field_value_lambda_f64_3", + [](T &self, + std::string field_name, + const std::function pos_to_val, + const i32 offset) { + return self.template set_field_value_lambda( + std::move(field_name), pos_to_val, offset); + }, + py::arg("field_name"), + py::arg("pos_to_val"), + py::arg("offset") = 0) .def("remap_positions", &T::remap_positions) //.def("set_field_value_lambda_f64_3",[](T&self,std::string field_name, const // std::function pos_to_val){ From 65c375d5de72d98274cb8f7e25a2ffb0536c36e2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Tue, 10 Mar 2026 15:59:02 +0100 Subject: [PATCH 14/24] clean --- src/shammodels/sph/src/modules/ComputeEos.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/shammodels/sph/src/modules/ComputeEos.cpp b/src/shammodels/sph/src/modules/ComputeEos.cpp index 26dcb34fbc..81de5b5275 100644 --- a/src/shammodels/sph/src/modules/ComputeEos.cpp +++ b/src/shammodels/sph/src/modules/ComputeEos.cpp @@ -115,8 +115,6 @@ struct RhoGetterSJ { epsilon_sum += s_j * s_j / rho; } - logger::raw_ln(i, rho, epsilon_sum, rho * (1 - epsilon_sum)); - return rho * (1 - epsilon_sum); } }; From daca6fa50d4eb38945f1466579120eb8f2a756f0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Tue, 10 Mar 2026 18:23:20 +0100 Subject: [PATCH 15/24] clean --- examples/sph/run_dusty.py | 2 +- .../sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp | 7 ++++++- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/examples/sph/run_dusty.py b/examples/sph/run_dusty.py index 4110cc8892..8989e986c4 100644 --- a/examples/sph/run_dusty.py +++ b/examples/sph/run_dusty.py @@ -46,7 +46,7 @@ def uint_g(r): ndust = 3 -rc = 0.3 +rc = 0.6 stopping_times = [1e-4, 1e-3, 1e-2] diff --git a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp index 2d64601e34..795d75f593 100644 --- a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp +++ b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp @@ -102,6 +102,7 @@ struct KernelUpdateDerivsMonofluidTVI { // eq 51, Hutchison 2018 ds_j_dt[thread_id] = Tscal{-0.5} * term1 + (s_j_a / (2 * rho_a * omega_a)) * term2; + } }; @@ -132,6 +133,10 @@ void shammodels::sph::modules::NodeUpdateDerivsMonofluidTVI:: using ComputeKernel = KernelUpdateDerivsMonofluidTVI; + auto total_specie_count = part_counts.template map([&](u64 id, u32 count) { + return count * ndust; + }); + // call the kernel for each patches with part_counts.get(id_patch) threads of patch id_patch sham::distributed_data_kernel_call( shamsys::instance::get_compute_scheduler_ptr(), @@ -145,7 +150,7 @@ void shammodels::sph::modules::NodeUpdateDerivsMonofluidTVI:: edges.Ttilde_sj.get_spans(), edges.neigh_cache}, sham::DDMultiRef{edges.ds_j_dt.get_spans()}, - part_counts, + total_specie_count, ComputeKernel{pmass, ndust}); } From d1b722cd2d5f78ba456181d30392573483f56414 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Wed, 11 Mar 2026 00:40:21 +0100 Subject: [PATCH 16/24] progress --- examples/sph/run_dusty.py | 24 +- examples/sph/run_dusty_settle.py | 235 ++++++++++++++++++ .../shammodels/common/ExtForceConfig.hpp | 44 +++- .../modules/AddForceVerticalDiscPotential.hpp | 97 ++++++++ .../common/modules/ForwardEuler.hpp | 36 ++- .../sph/modules/ComputeDustTtilde.hpp | 9 +- .../shammodels/sph/modules/ExternalForces.hpp | 9 +- .../modules/SetDustStoppingTimeConstant.hpp | 2 +- src/shammodels/sph/src/Solver.cpp | 7 +- .../sph/src/modules/ExternalForces.cpp | 40 +++ .../modules/NodeUpdateDerivsMonofluidTVI.cpp | 4 +- src/shammodels/sph/src/pySPHModel.cpp | 8 + 12 files changed, 477 insertions(+), 38 deletions(-) create mode 100644 examples/sph/run_dusty_settle.py create mode 100644 src/shammodels/common/include/shammodels/common/modules/AddForceVerticalDiscPotential.hpp diff --git a/examples/sph/run_dusty.py b/examples/sph/run_dusty.py index 8989e986c4..c1caf9e782 100644 --- a/examples/sph/run_dusty.py +++ b/examples/sph/run_dusty.py @@ -45,9 +45,9 @@ def uint_g(r): return P / ((gamma - 1) * rho_g) -ndust = 3 -rc = 0.6 -stopping_times = [1e-4, 1e-3, 1e-2] +ndust = 1 +rc = 0.25 +stopping_times = [ 1e-1] bmin = (-0.6, -0.6, -0.6) @@ -85,7 +85,7 @@ def uint_g(r): cfg.set_dust_mode_monofluid_tvi(ndust) cfg.set_dust_stopping_times(stopping_times) cfg.set_boundary_periodic() -cfg.set_eos_adiabatic(gamma) +cfg.set_eos_isothermal(1) cfg.print_status() model.set_solver_config(cfg) @@ -125,10 +125,11 @@ def func_s(r): model.timestep() +tnext = 0 for j in range(10): if j > 0: - for i in range(10): - model.timestep() + tnext += 0.1 + model.evolve_until(tnext) dic = ctx.collect_data() print(dic["s_j"]) @@ -159,8 +160,11 @@ def func_s(r): axs[2].scatter(r, cs, label="soundspeed") axs[2].legend() - axs[3].scatter(r, ds_j_dt[:, 0], label=f"ds_j_dt_{0}") - axs[3].scatter(r, ds_j_dt[:, 1], label=f"ds_j_dt_{1}") - axs[3].scatter(r, ds_j_dt[:, 2], label=f"ds_j_dt_{2}") + for i in range(ndust): + axs[3].scatter(r, ds_j_dt[:, i], label=f"ds_j_dt_{i}") axs[3].legend() - plt.show() + + for k in range(4): + axs[k].set_xlim(-0.1,1.1) + plt.savefig(f"mono_{j}.png") + plt.close() diff --git a/examples/sph/run_dusty_settle.py b/examples/sph/run_dusty_settle.py new file mode 100644 index 0000000000..d019c71e1a --- /dev/null +++ b/examples/sph/run_dusty_settle.py @@ -0,0 +1,235 @@ +""" +Testing Sod tube with SPH +========================= + +CI test for Sod tube with SPH +""" + +import matplotlib.pyplot as plt + +import shamrock + +shamrock.enable_experimental_features() +import numpy as np + +gamma = 1.4 +rho_i = 1 +central_mass = 1 +R0 = 1 +H_r_0 = 0.05 + +si = shamrock.UnitSystem() +sicte = shamrock.Constants(si) +codeu = shamrock.UnitSystem( + unit_time=3600 * 24 * 365, + unit_length=sicte.au(), + unit_mass=sicte.sol_mass(), +) +ucte = shamrock.Constants(codeu) +G = ucte.G() + +def kep_profile(r): + return (G * central_mass / r) ** 0.5 + + +def omega_k(r): + return kep_profile(r) / r + +H = H_r_0 * R0 +cs = H * omega_k(R0) +box = 8*H + +print(f"cs = {cs}") +print(f"H = {H}") + + +def scaling_rho(r): + x,y,z = r + + loc_h = H/1.44 + gaussian = np.exp(-(y**2) / (2*loc_h*loc_h)) / (loc_h * np.sqrt(2 * np.pi)) + return gaussian + +def func_rho_t(r): + return rho_i * scaling_rho(r) + +def func_rho_d_j(r, idust): + return (0.1 / ndust) * rho_i * scaling_rho(r) + + +def func_rho_g(r): + return rho_i*scaling_rho(r) - sum([func_rho_d_j(r, i) for i in range(ndust)]) + + +def func_s_j(r, idust): + rho_t = func_rho_t(r) + rho_d_j = [func_rho_d_j(r, i) for i in range(ndust)] + eps_j = rho_d_j[idust] / rho_t + return np.sqrt(rho_t * eps_j) + + +cs_g = 1 + + +def uint_g(r): + rho_g = func_rho_g(r) + P = rho_g * cs_g * cs_g / gamma + return P / ((gamma - 1) * rho_g) + + +ndust = 4 +rc = 0.25 +stopping_times = np.logspace(-4,-1,ndust)*omega_k(R0) +from scipy.special import erfinv + +bmin = (-box/4, -box, -box/4) +bmax = (box/4, box, box/4) + +N_target = 1e5 +scheduler_split_val = int(2e7) +scheduler_merge_val = int(1) + +xm, ym, zm = bmin +xM, yM, zM = bmax +vol_b = (xM - xm) * (yM - ym) * (zM - zm) + +part_vol = vol_b / N_target + +# lattice volume +HCP_PACKING_DENSITY = 0.74 +part_vol_lattice = HCP_PACKING_DENSITY * part_vol + +dr = (part_vol_lattice / ((4.0 / 3.0) * np.pi)) ** (1.0 / 3.0) + +pmass = -1 + +ctx = shamrock.Context() +ctx.pdata_layout_new() + +model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M6") + +cfg = model.gen_default_config() +# cfg.set_artif_viscosity_Constant(alpha_u = 1, alpha_AV = 1, beta_AV = 2) +# cfg.set_artif_viscosity_VaryingMM97(alpha_min = 0.1,alpha_max = 1,sigma_decay = 0.1, alpha_u = 1, beta_AV = 2) +cfg.set_artif_viscosity_VaryingCD10( + alpha_min=0.0, alpha_max=1, sigma_decay=0.1, alpha_u=1, beta_AV=2 +) +cfg.set_dust_mode_monofluid_tvi(ndust) +cfg.set_dust_stopping_times(stopping_times) +cfg.add_ext_force_vertical_disc_potential(central_mass = 1, R0 = 1) +cfg.set_boundary_periodic() +cfg.set_units(codeu) +cfg.set_eos_isothermal(cs) +cfg.print_status() +model.set_solver_config(cfg) + +model.init_scheduler(int(1e8), 1) + + +bmin, bmax = model.get_ideal_hcp_box(dr, bmin, bmax) +xm, ym, zm = bmin +xM, yM, zM = bmax + +model.resize_simulation_box(bmin, bmax) + +setup = model.get_setup() +gen = setup.make_generator_lattice_hcp(dr, bmin, bmax) +setup.apply_setup(gen, insert_step=scheduler_split_val) + +vol_b = (xM - xm) * (yM - ym) * (zM - zm) +totmass = rho_i * vol_b +print("Total mass :", totmass) + +pmass = model.total_mass_to_part_mass(totmass) +model.set_particle_mass(pmass) +print("Current part mass :", pmass) + +# Correct the barycenter +analysis_barycenter = shamrock.model_sph.analysisBarycenter(model=model) +barycenter, disc_mass = analysis_barycenter.get_barycenter() + +if shamrock.sys.world_rank() == 0: + print(f"disc barycenter = {barycenter}") + +model.apply_position_offset((-barycenter[0], -barycenter[1], -barycenter[2])) + +def f_remap(r): + x,y,z = r + + rn = max(abs(yM),abs(ym)) + print(y,H, H*erfinv(y/rn) ) + y = H*erfinv(y/rn) + + y = min(y,yM) + y = max(y,ym) + return (x,y,z) + +model.remap_positions(f_remap) + + +for i in range(ndust): + + def func_s(r): + return func_s_j(r, i) + + model.set_field_value_lambda_f64("s_j", func_s, i) + +model.set_field_value_lambda_f64("uint", uint_g) + + + +model.set_cfl_cour(0.1) +model.set_cfl_force(0.1) + +model.timestep() + +tnext = 0 +for j in range(1000): + if j > 0: + tnext += 0.02 + model.evolve_until(tnext) + + dic = ctx.collect_data() + print(dic["s_j"]) + + print(dic["xyz"].shape) + + x = dic["xyz"][:, 0] + y = dic["xyz"][:, 1] + z = dic["xyz"][:, 2] + s_j = dic["s_j"].reshape(-1, ndust) + ds_j_dt = dic["ds_j_dt"].reshape(-1, ndust) + cs = dic["soundspeed"] + + print(s_j) + + r = y + + hpart = dic["hpart"] + rho = pmass * (model.get_hfact() / np.array(hpart)) ** 3 + + print("compute original rho") + estimated_rho = [func_rho_t(dic["xyz"][kk]) for kk in range(len(dic["xyz"]))] + + sz = 2 + + fig, axs = plt.subplots(nrows=1, ncols=4, figsize=(15, 5)) + axs[0].scatter(y, rho, label="rho",s=sz) + for i in range(ndust): + axs[0].scatter(y, s_j[:, i] ** 2, label=f"rho_j_{i}",s=sz) + #axs[0].scatter(y,estimated_rho) + + axs[0].legend() + for i in range(ndust): + axs[1].scatter(y, s_j[:, i] ** 2 / rho, label=f"eps_j_{i}",s=sz) + axs[1].legend() + + axs[2].scatter(y, cs, label="soundspeed",s=sz) + axs[2].legend() + + for i in range(ndust): + axs[3].scatter(y, ds_j_dt[:, i], label=f"ds_j_dt_{i}",s=sz) + axs[3].legend() + + plt.savefig(f"mono_{j}.png") + plt.close() diff --git a/src/shammodels/common/include/shammodels/common/ExtForceConfig.hpp b/src/shammodels/common/include/shammodels/common/ExtForceConfig.hpp index 21308286fc..fea3439460 100644 --- a/src/shammodels/common/include/shammodels/common/ExtForceConfig.hpp +++ b/src/shammodels/common/include/shammodels/common/ExtForceConfig.hpp @@ -68,7 +68,14 @@ namespace shammodels { : shear_base(shear_base), shear_dir(shear_dir), Omega_0(Omega_0), eta(eta), q(q) {}; }; - using VariantForce = std::variant; + /// -GMy / sqrt(R0^2 + y^2) + struct VerticalDiscPotential { + Tscal central_mass; + Tscal R0; + }; + + using VariantForce + = std::variant; VariantForce val; }; @@ -78,9 +85,10 @@ namespace shammodels { using Tscal = shambase::VecComponent; static constexpr u32 dim = shambase::VectorProperties::dimension; - using PointMass = typename ExtForceVariant::PointMass; - using LenseThirring = typename ExtForceVariant::LenseThirring; - using ShearingBoxForce = typename ExtForceVariant::ShearingBoxForce; + using PointMass = typename ExtForceVariant::PointMass; + using LenseThirring = typename ExtForceVariant::LenseThirring; + using ShearingBoxForce = typename ExtForceVariant::ShearingBoxForce; + using VerticalDiscPotential = typename ExtForceVariant::VerticalDiscPotential; using VariantForce = std::variant; @@ -108,6 +116,10 @@ namespace shammodels { ext_forces.push_back(ExtForceVariant{ShearingBoxForce{Omega_0, eta, q}}); } + + inline void add_vertical_disc_potential(Tscal central_mass, Tscal R0) { + ext_forces.push_back(ExtForceVariant{VerticalDiscPotential{central_mass, R0}}); + } }; } // namespace shammodels @@ -117,9 +129,10 @@ namespace shammodels { inline void to_json(nlohmann::json &j, const ExtForceVariant &p) { using T = ExtForceVariant; - using PointMass = typename T::PointMass; - using LenseThirring = typename T::LenseThirring; - using ShearingBoxForce = typename T::ShearingBoxForce; + using PointMass = typename T::PointMass; + using LenseThirring = typename T::LenseThirring; + using ShearingBoxForce = typename T::ShearingBoxForce; + using VerticalDiscPotential = typename T::VerticalDiscPotential; if (const PointMass *v = std::get_if(&p.val)) { j = { @@ -141,6 +154,11 @@ namespace shammodels { {"eta", v->eta}, {"q", v->q}, }; + } else if (const VerticalDiscPotential *v = std::get_if(&p.val)) { + j + = {{"force_type", "vertical_disc_potential"}, + {"central_mass", v->central_mass}, + {"R0", v->R0}}; } else { shambase::throw_unimplemented(); } @@ -158,9 +176,10 @@ namespace shammodels { std::string force_type; j.at("force_type").get_to(force_type); - using PointMass = typename T::PointMass; - using LenseThirring = typename T::LenseThirring; - using ShearingBoxForce = typename T::ShearingBoxForce; + using PointMass = typename T::PointMass; + using LenseThirring = typename T::LenseThirring; + using ShearingBoxForce = typename T::ShearingBoxForce; + using VerticalDiscPotential = typename T::VerticalDiscPotential; if (force_type == "point_mass") { p.val = PointMass{ @@ -182,6 +201,11 @@ namespace shammodels { j.at("eta").get(), j.at("q").get(), }; + } else if (force_type == "vertical_disc_potential") { + p.val = VerticalDiscPotential{ + j.at("central_mass").get(), + j.at("R0").get(), + }; } else { shambase::throw_unimplemented("wtf !"); } diff --git a/src/shammodels/common/include/shammodels/common/modules/AddForceVerticalDiscPotential.hpp b/src/shammodels/common/include/shammodels/common/modules/AddForceVerticalDiscPotential.hpp new file mode 100644 index 0000000000..9fe1c612c0 --- /dev/null +++ b/src/shammodels/common/include/shammodels/common/modules/AddForceVerticalDiscPotential.hpp @@ -0,0 +1,97 @@ +// -------------------------------------------------------// +// +// 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 AddForceCentralGravPotential.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief Adds the acceleration from a central gravitational potential (point mass). + * + */ + +#include "shambackends/kernel_call_distrib.hpp" +#include "shambackends/vec.hpp" +#include "shamrock/solvergraph/IDataEdge.hpp" +#include "shamrock/solvergraph/IFieldSpan.hpp" +#include "shamrock/solvergraph/INode.hpp" +#include "shamrock/solvergraph/Indexes.hpp" +#include "shamsys/NodeInstance.hpp" + +namespace shammodels::common::modules { + template + class AddForceVerticalDiscPotential : public shamrock::solvergraph::INode { + + using Tscal = shambase::VecComponent; + + public: + AddForceVerticalDiscPotential() = default; + + struct Edges { + const shamrock::solvergraph::IDataEdge &constant_G; + const shamrock::solvergraph::IDataEdge ¢ral_mass; + const shamrock::solvergraph::IDataEdge &R0; + const shamrock::solvergraph::IFieldSpan &spans_positions; + const shamrock::solvergraph::Indexes &sizes; + shamrock::solvergraph::IFieldSpan &spans_accel_ext; + }; + + inline void set_edges( + std::shared_ptr> constant_G, + std::shared_ptr> central_mass, + std::shared_ptr> R0, + std::shared_ptr> spans_positions, + std::shared_ptr> sizes, + std::shared_ptr> spans_accel_ext) { + __internal_set_ro_edges({constant_G, central_mass, R0, spans_positions, sizes}); + __internal_set_rw_edges({spans_accel_ext}); + } + + inline Edges get_edges() { + return Edges{ + get_ro_edge>(0), + get_ro_edge>(1), + get_ro_edge>(2), + get_ro_edge>(3), + get_ro_edge>(4), + get_rw_edge>(0)}; + } + + inline void _impl_evaluate_internal() { + + __shamrock_stack_entry(); + + auto edges = get_edges(); + + edges.spans_positions.check_sizes(edges.sizes.indexes); + edges.spans_accel_ext.ensure_sizes(edges.sizes.indexes); + + Tscal cmass = edges.central_mass.data; + Tscal G = edges.constant_G.data; + Tscal R0 = edges.R0.data; + + sham::distributed_data_kernel_call( + shamsys::instance::get_compute_scheduler_ptr(), + sham::DDMultiRef{edges.spans_positions.get_spans()}, + sham::DDMultiRef{edges.spans_accel_ext.get_spans()}, + edges.sizes.indexes, + [mGM = -cmass * G, R02 = R0 * R0](u32 gid, const Tvec *xyz, Tvec *axyz_ext) { + Tscal y_a = xyz[gid].y(); + axyz_ext[gid] += mGM * y_a / sycl::sqrt(R02 + y_a * y_a); + }); + } + + inline virtual std::string _impl_get_label() const { + return "AddForceVerticalDiscPotential"; + }; + + inline virtual std::string _impl_get_tex() const { return "TODO"; }; + }; + +} // namespace shammodels::common::modules diff --git a/src/shammodels/common/include/shammodels/common/modules/ForwardEuler.hpp b/src/shammodels/common/include/shammodels/common/modules/ForwardEuler.hpp index b0f90caa1f..6450cdeef2 100644 --- a/src/shammodels/common/include/shammodels/common/modules/ForwardEuler.hpp +++ b/src/shammodels/common/include/shammodels/common/modules/ForwardEuler.hpp @@ -31,8 +31,10 @@ namespace shammodels::common::modules { using Tscal = shambase::VecComponent; + u32 nvar; + public: - ForwardEuler() = default; + ForwardEuler(u32 nvar = 1) : nvar(nvar) {} struct Edges { const shamrock::solvergraph::IDataEdge &dt; @@ -68,14 +70,32 @@ namespace shammodels::common::modules { Tscal dt = edges.dt.data; - sham::distributed_data_kernel_call( - shamsys::instance::get_compute_scheduler_ptr(), - sham::DDMultiRef{edges.time_derivative.get_spans()}, - sham::DDMultiRef{edges.field.get_spans()}, - edges.sizes.indexes, - [dt](u32 gid, const T *time_derivative, T *field) { - field[gid] = field[gid] + dt * time_derivative[gid]; + if (nvar == 1) { + + sham::distributed_data_kernel_call( + shamsys::instance::get_compute_scheduler_ptr(), + sham::DDMultiRef{edges.time_derivative.get_spans()}, + sham::DDMultiRef{edges.field.get_spans()}, + edges.sizes.indexes, + [dt](u32 gid, const T *time_derivative, T *field) { + field[gid] = field[gid] + dt * time_derivative[gid]; + }); + + } else { + + auto var_count = edges.sizes.indexes.template map([&](u64 id, u32 count) { + return count * nvar; }); + + sham::distributed_data_kernel_call( + shamsys::instance::get_compute_scheduler_ptr(), + sham::DDMultiRef{edges.time_derivative.get_spans()}, + sham::DDMultiRef{edges.field.get_spans()}, + var_count, + [dt](u32 gid, const T *time_derivative, T *field) { + field[gid] = field[gid] + dt * time_derivative[gid]; + }); + } } inline virtual std::string _impl_get_label() const { return "ForwardEuler"; }; diff --git a/src/shammodels/sph/include/shammodels/sph/modules/ComputeDustTtilde.hpp b/src/shammodels/sph/include/shammodels/sph/modules/ComputeDustTtilde.hpp index 122796ccec..085dcd6652 100644 --- a/src/shammodels/sph/include/shammodels/sph/modules/ComputeDustTtilde.hpp +++ b/src/shammodels/sph/include/shammodels/sph/modules/ComputeDustTtilde.hpp @@ -18,6 +18,7 @@ #include "shambackends/kernel_call_distrib.hpp" #include "shambackends/vec.hpp" +#include "shamcomm/logs.hpp" #include "shammodels/sph/math/density.hpp" #include "shamrock/solvergraph/IFieldSpan.hpp" #include "shamrock/solvergraph/INode.hpp" @@ -75,6 +76,10 @@ namespace shammodels::sph::modules { const Tscal pmass = edges.gpart_mass.value; + auto total_specie_count = part_counts.template map([&](u64 id, u32 count) { + return count * ndust; + }); + // call the kernel for each patches with part_counts.get(id_patch) threads of patch // id_patch sham::distributed_data_kernel_call( @@ -82,7 +87,7 @@ namespace shammodels::sph::modules { sham::DDMultiRef{ edges.hpart.get_spans(), edges.s_j.get_spans(), edges.t_j.get_spans()}, sham::DDMultiRef{edges.Ttilde_sj.get_spans()}, - part_counts, + total_specie_count, [pmass, ndust = ndust]( u32 thread_id, const Tscal *__restrict hpart, @@ -118,6 +123,8 @@ namespace shammodels::sph::modules { Ttilde_sj_a -= eps_k_a * eps_k_a * tk_a; } + // logger::raw_ln("Ttilde_sj_a", jdust, Ttilde_sj_a, eps_j_a, tj_a); + Ttilde_sj[thread_id] = Ttilde_sj_a; }); } diff --git a/src/shammodels/sph/include/shammodels/sph/modules/ExternalForces.hpp b/src/shammodels/sph/include/shammodels/sph/modules/ExternalForces.hpp index 94588b99ff..4918f7c105 100644 --- a/src/shammodels/sph/include/shammodels/sph/modules/ExternalForces.hpp +++ b/src/shammodels/sph/include/shammodels/sph/modules/ExternalForces.hpp @@ -57,10 +57,11 @@ namespace shammodels::sph::modules { void point_mass_accrete_particles(); private: - using SolverConfigExtForce = typename Config::ExtForceConfig; - using EF_PointMass = typename SolverConfigExtForce::PointMass; - using EF_LenseThirring = typename SolverConfigExtForce::LenseThirring; - using EF_ShearingBoxForce = typename SolverConfigExtForce::ShearingBoxForce; + using SolverConfigExtForce = typename Config::ExtForceConfig; + using EF_PointMass = typename SolverConfigExtForce::PointMass; + using EF_LenseThirring = typename SolverConfigExtForce::LenseThirring; + using EF_ShearingBoxForce = typename SolverConfigExtForce::ShearingBoxForce; + using EF_VerticalDiscPotential = typename SolverConfigExtForce::VerticalDiscPotential; inline PatchScheduler &scheduler() { return shambase::get_check_ref(context.sched); } }; diff --git a/src/shammodels/sph/include/shammodels/sph/modules/SetDustStoppingTimeConstant.hpp b/src/shammodels/sph/include/shammodels/sph/modules/SetDustStoppingTimeConstant.hpp index fb19373968..cc5e70c863 100644 --- a/src/shammodels/sph/include/shammodels/sph/modules/SetDustStoppingTimeConstant.hpp +++ b/src/shammodels/sph/include/shammodels/sph/modules/SetDustStoppingTimeConstant.hpp @@ -78,7 +78,7 @@ namespace shammodels::sph::modules { q, sham::MultiRef{t_j_0}, sham::MultiRef{edges.t_j.get_spans().get(id)}, - part_counts.get(id), + part_counts.get(id) * ndust, [ndust = ndust](u32 thread_id, const Tscal *__restrict t_j_0, Tscal *__restrict t_j) { u32 jdust = thread_id % ndust; diff --git a/src/shammodels/sph/src/Solver.cpp b/src/shammodels/sph/src/Solver.cpp index b1eed9277a..40fcda89b9 100644 --- a/src/shammodels/sph/src/Solver.cpp +++ b/src/shammodels/sph/src/Solver.cpp @@ -439,8 +439,9 @@ void shammodels::sph::Solver::init_solver_graph() { } if (has_s_j_field) { + u32 ndust = solver_config.dust_config.get_dust_nvar(); auto half_step_s_j = solver_graph.register_node( - prefix + "_s_j", shammodels::common::modules::ForwardEuler{}); + prefix + "_s_j", shammodels::common::modules::ForwardEuler(ndust)); shambase::get_check_ref(half_step_s_j) .set_edges( solver_graph.get_edge_ptr>("dt_half"), @@ -2514,8 +2515,8 @@ shammodels::sph::TimestepLog shammodels::sph::Solver::evolve_once() auto cfl_dt = cfl_dt_buf.get_write_access(depends_list); auto e = q.submit(depends_list, [&](sycl::handler &cgh) { - Tscal C_cour = solver_config.cfl_config.cfl_cour - * solver_config.time_state.cfl_multiplier; + Tscal C_cour = solver_config.cfl_config.cfl_cour + * solver_config.time_state.cfl_multiplier; Tscal C_force = solver_config.cfl_config.cfl_force * solver_config.time_state.cfl_multiplier; diff --git a/src/shammodels/sph/src/modules/ExternalForces.cpp b/src/shammodels/sph/src/modules/ExternalForces.cpp index e906430e27..70392490e7 100644 --- a/src/shammodels/sph/src/modules/ExternalForces.cpp +++ b/src/shammodels/sph/src/modules/ExternalForces.cpp @@ -18,11 +18,13 @@ #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/common/modules/AddForceLenseThirring.hpp" #include "shammodels/common/modules/AddForceShearingBoxInertialPart.hpp" #include "shammodels/common/modules/AddForceShearingBoxNonInertial.hpp" +#include "shammodels/common/modules/AddForceVerticalDiscPotential.hpp" #include "shammodels/sph/modules/ExternalForces.hpp" #include "shammodels/sph/modules/SinkParticlesUpdate.hpp" #include "shamrock/solvergraph/IDataEdge.hpp" @@ -200,6 +202,38 @@ void shammodels::sph::modules::ExternalForces::compute_ext_forc shambase::to_shared(std::move(set_eta)), shambase::to_shared(std::move(add_force_shearing_box_inertial_part))})); + } else if ( + EF_VerticalDiscPotential *ext_force + = std::get_if(&var_force.val)) { + + auto central_mass = shamrock::solvergraph::IDataEdge::make_shared("", ""); + auto R0 = 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_R0( + [r = ext_force->R0](shamrock::solvergraph::IDataEdge &R0) { + R0.data = r; // no support for offset yet + }); + set_R0.set_edges(R0); + + common::modules::AddForceVerticalDiscPotential add_force_central_grav_potential; + add_force_central_grav_potential.set_edges( + constant_G, central_mass, R0, 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_R0)), + 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 ..."); } @@ -300,6 +334,9 @@ void shammodels::sph::modules::ExternalForces::add_ext_forces() is_c_needed = true; } else if ( EF_ShearingBoxForce *ext_force = std::get_if(&var_force.val)) { + } else if ( + EF_VerticalDiscPotential *ext_force + = std::get_if(&var_force.val)) { } else { shambase::throw_unimplemented("this force is not handled, yet ..."); } @@ -451,6 +488,9 @@ void shammodels::sph::modules::ExternalForces::add_ext_forces() add_ext_forces_seq.push_back(set_q); add_ext_forces_seq.push_back(solver_graph.get_node_ptr_base(prefix_shearing_box)); + } else if ( + EF_VerticalDiscPotential *ext_force + = std::get_if(&var_force.val)) { } else { shambase::throw_unimplemented("this force is not handled, yet ..."); } diff --git a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp index 795d75f593..1ba1832063 100644 --- a/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp +++ b/src/shammodels/sph/src/modules/NodeUpdateDerivsMonofluidTVI.cpp @@ -16,6 +16,7 @@ #include "shammodels/sph/modules/NodeUpdateDerivsMonofluidTVI.hpp" #include "shambackends/kernel_call_distrib.hpp" +#include "shamcomm/logs.hpp" #include "shammath/sphkernels.hpp" #include "shammodels/sph/math/density.hpp" #include "shamrock/patch/PatchDataField.hpp" @@ -96,13 +97,14 @@ struct KernelUpdateDerivsMonofluidTVI { Tscal delta_P = P_a - P_b; Tscal Ts_weighted = (Ttilde_sj_a / rho_a + Ttilde_sj_b / rho_b); + // logger::raw_ln("Ts_weighted", Ts_weighted); + term1 += (pmass * s_j_b / rho_b) * Ts_weighted * delta_P * F_ab_bar; term2 += pmass * sham::dot(v_ab, r_ab_unit * Fab_a); }); // eq 51, Hutchison 2018 ds_j_dt[thread_id] = Tscal{-0.5} * term1 + (s_j_a / (2 * rho_a * omega_a)) * term2; - } }; diff --git a/src/shammodels/sph/src/pySPHModel.cpp b/src/shammodels/sph/src/pySPHModel.cpp index 2154aa4860..3035d4d127 100644 --- a/src/shammodels/sph/src/pySPHModel.cpp +++ b/src/shammodels/sph/src/pySPHModel.cpp @@ -286,6 +286,14 @@ void add_instance(py::module &m, std::string name_config, std::string name_model py::arg("Omega_0"), py::arg("eta"), py::arg("q")) + .def( + "add_ext_force_vertical_disc_potential", + [](TConfig &self, Tscal central_mass, Tscal R0) { + self.ext_force_config.add_vertical_disc_potential(central_mass, R0); + }, + py::kw_only(), + py::arg("central_mass"), + py::arg("R0")) .def("set_units", &TConfig::set_units) .def( "get_units", From ec0cff259a2b477ed820e53aea8c4c37ab9823ec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Wed, 11 Mar 2026 00:45:59 +0100 Subject: [PATCH 17/24] [SPH] add disc vertical potential for stratified boxes --- .../shammodels/common/ExtForceConfig.hpp | 44 +++++++-- .../modules/AddForceVerticalDiscPotential.hpp | 97 +++++++++++++++++++ .../shammodels/sph/modules/ExternalForces.hpp | 9 +- .../sph/src/modules/ExternalForces.cpp | 40 ++++++++ src/shammodels/sph/src/pySPHModel.cpp | 8 ++ 5 files changed, 184 insertions(+), 14 deletions(-) create mode 100644 src/shammodels/common/include/shammodels/common/modules/AddForceVerticalDiscPotential.hpp diff --git a/src/shammodels/common/include/shammodels/common/ExtForceConfig.hpp b/src/shammodels/common/include/shammodels/common/ExtForceConfig.hpp index 21308286fc..fea3439460 100644 --- a/src/shammodels/common/include/shammodels/common/ExtForceConfig.hpp +++ b/src/shammodels/common/include/shammodels/common/ExtForceConfig.hpp @@ -68,7 +68,14 @@ namespace shammodels { : shear_base(shear_base), shear_dir(shear_dir), Omega_0(Omega_0), eta(eta), q(q) {}; }; - using VariantForce = std::variant; + /// -GMy / sqrt(R0^2 + y^2) + struct VerticalDiscPotential { + Tscal central_mass; + Tscal R0; + }; + + using VariantForce + = std::variant; VariantForce val; }; @@ -78,9 +85,10 @@ namespace shammodels { using Tscal = shambase::VecComponent; static constexpr u32 dim = shambase::VectorProperties::dimension; - using PointMass = typename ExtForceVariant::PointMass; - using LenseThirring = typename ExtForceVariant::LenseThirring; - using ShearingBoxForce = typename ExtForceVariant::ShearingBoxForce; + using PointMass = typename ExtForceVariant::PointMass; + using LenseThirring = typename ExtForceVariant::LenseThirring; + using ShearingBoxForce = typename ExtForceVariant::ShearingBoxForce; + using VerticalDiscPotential = typename ExtForceVariant::VerticalDiscPotential; using VariantForce = std::variant; @@ -108,6 +116,10 @@ namespace shammodels { ext_forces.push_back(ExtForceVariant{ShearingBoxForce{Omega_0, eta, q}}); } + + inline void add_vertical_disc_potential(Tscal central_mass, Tscal R0) { + ext_forces.push_back(ExtForceVariant{VerticalDiscPotential{central_mass, R0}}); + } }; } // namespace shammodels @@ -117,9 +129,10 @@ namespace shammodels { inline void to_json(nlohmann::json &j, const ExtForceVariant &p) { using T = ExtForceVariant; - using PointMass = typename T::PointMass; - using LenseThirring = typename T::LenseThirring; - using ShearingBoxForce = typename T::ShearingBoxForce; + using PointMass = typename T::PointMass; + using LenseThirring = typename T::LenseThirring; + using ShearingBoxForce = typename T::ShearingBoxForce; + using VerticalDiscPotential = typename T::VerticalDiscPotential; if (const PointMass *v = std::get_if(&p.val)) { j = { @@ -141,6 +154,11 @@ namespace shammodels { {"eta", v->eta}, {"q", v->q}, }; + } else if (const VerticalDiscPotential *v = std::get_if(&p.val)) { + j + = {{"force_type", "vertical_disc_potential"}, + {"central_mass", v->central_mass}, + {"R0", v->R0}}; } else { shambase::throw_unimplemented(); } @@ -158,9 +176,10 @@ namespace shammodels { std::string force_type; j.at("force_type").get_to(force_type); - using PointMass = typename T::PointMass; - using LenseThirring = typename T::LenseThirring; - using ShearingBoxForce = typename T::ShearingBoxForce; + using PointMass = typename T::PointMass; + using LenseThirring = typename T::LenseThirring; + using ShearingBoxForce = typename T::ShearingBoxForce; + using VerticalDiscPotential = typename T::VerticalDiscPotential; if (force_type == "point_mass") { p.val = PointMass{ @@ -182,6 +201,11 @@ namespace shammodels { j.at("eta").get(), j.at("q").get(), }; + } else if (force_type == "vertical_disc_potential") { + p.val = VerticalDiscPotential{ + j.at("central_mass").get(), + j.at("R0").get(), + }; } else { shambase::throw_unimplemented("wtf !"); } diff --git a/src/shammodels/common/include/shammodels/common/modules/AddForceVerticalDiscPotential.hpp b/src/shammodels/common/include/shammodels/common/modules/AddForceVerticalDiscPotential.hpp new file mode 100644 index 0000000000..c1c8f709fe --- /dev/null +++ b/src/shammodels/common/include/shammodels/common/modules/AddForceVerticalDiscPotential.hpp @@ -0,0 +1,97 @@ +// -------------------------------------------------------// +// +// 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 AddForceVerticalDiscPotential.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief Adds the acceleration from a central gravitational potential (point mass). + * + */ + +#include "shambackends/kernel_call_distrib.hpp" +#include "shambackends/vec.hpp" +#include "shamrock/solvergraph/IDataEdge.hpp" +#include "shamrock/solvergraph/IFieldSpan.hpp" +#include "shamrock/solvergraph/INode.hpp" +#include "shamrock/solvergraph/Indexes.hpp" +#include "shamsys/NodeInstance.hpp" + +namespace shammodels::common::modules { + template + class AddForceVerticalDiscPotential : public shamrock::solvergraph::INode { + + using Tscal = shambase::VecComponent; + + public: + AddForceVerticalDiscPotential() = default; + + struct Edges { + const shamrock::solvergraph::IDataEdge &constant_G; + const shamrock::solvergraph::IDataEdge ¢ral_mass; + const shamrock::solvergraph::IDataEdge &R0; + const shamrock::solvergraph::IFieldSpan &spans_positions; + const shamrock::solvergraph::Indexes &sizes; + shamrock::solvergraph::IFieldSpan &spans_accel_ext; + }; + + inline void set_edges( + std::shared_ptr> constant_G, + std::shared_ptr> central_mass, + std::shared_ptr> R0, + std::shared_ptr> spans_positions, + std::shared_ptr> sizes, + std::shared_ptr> spans_accel_ext) { + __internal_set_ro_edges({constant_G, central_mass, R0, spans_positions, sizes}); + __internal_set_rw_edges({spans_accel_ext}); + } + + inline Edges get_edges() { + return Edges{ + get_ro_edge>(0), + get_ro_edge>(1), + get_ro_edge>(2), + get_ro_edge>(3), + get_ro_edge>(4), + get_rw_edge>(0)}; + } + + inline void _impl_evaluate_internal() { + + __shamrock_stack_entry(); + + auto edges = get_edges(); + + edges.spans_positions.check_sizes(edges.sizes.indexes); + edges.spans_accel_ext.ensure_sizes(edges.sizes.indexes); + + Tscal cmass = edges.central_mass.data; + Tscal G = edges.constant_G.data; + Tscal R0 = edges.R0.data; + + sham::distributed_data_kernel_call( + shamsys::instance::get_compute_scheduler_ptr(), + sham::DDMultiRef{edges.spans_positions.get_spans()}, + sham::DDMultiRef{edges.spans_accel_ext.get_spans()}, + edges.sizes.indexes, + [mGM = -cmass * G, R02 = R0 * R0](u32 gid, const Tvec *xyz, Tvec *axyz_ext) { + Tscal y_a = xyz[gid].y(); + axyz_ext[gid] += mGM * y_a / sycl::sqrt(R02 + y_a * y_a); + }); + } + + inline virtual std::string _impl_get_label() const { + return "AddForceVerticalDiscPotential"; + }; + + inline virtual std::string _impl_get_tex() const { return "TODO"; }; + }; + +} // namespace shammodels::common::modules diff --git a/src/shammodels/sph/include/shammodels/sph/modules/ExternalForces.hpp b/src/shammodels/sph/include/shammodels/sph/modules/ExternalForces.hpp index 94588b99ff..4918f7c105 100644 --- a/src/shammodels/sph/include/shammodels/sph/modules/ExternalForces.hpp +++ b/src/shammodels/sph/include/shammodels/sph/modules/ExternalForces.hpp @@ -57,10 +57,11 @@ namespace shammodels::sph::modules { void point_mass_accrete_particles(); private: - using SolverConfigExtForce = typename Config::ExtForceConfig; - using EF_PointMass = typename SolverConfigExtForce::PointMass; - using EF_LenseThirring = typename SolverConfigExtForce::LenseThirring; - using EF_ShearingBoxForce = typename SolverConfigExtForce::ShearingBoxForce; + using SolverConfigExtForce = typename Config::ExtForceConfig; + using EF_PointMass = typename SolverConfigExtForce::PointMass; + using EF_LenseThirring = typename SolverConfigExtForce::LenseThirring; + using EF_ShearingBoxForce = typename SolverConfigExtForce::ShearingBoxForce; + using EF_VerticalDiscPotential = typename SolverConfigExtForce::VerticalDiscPotential; inline PatchScheduler &scheduler() { return shambase::get_check_ref(context.sched); } }; diff --git a/src/shammodels/sph/src/modules/ExternalForces.cpp b/src/shammodels/sph/src/modules/ExternalForces.cpp index e906430e27..70392490e7 100644 --- a/src/shammodels/sph/src/modules/ExternalForces.cpp +++ b/src/shammodels/sph/src/modules/ExternalForces.cpp @@ -18,11 +18,13 @@ #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/common/modules/AddForceLenseThirring.hpp" #include "shammodels/common/modules/AddForceShearingBoxInertialPart.hpp" #include "shammodels/common/modules/AddForceShearingBoxNonInertial.hpp" +#include "shammodels/common/modules/AddForceVerticalDiscPotential.hpp" #include "shammodels/sph/modules/ExternalForces.hpp" #include "shammodels/sph/modules/SinkParticlesUpdate.hpp" #include "shamrock/solvergraph/IDataEdge.hpp" @@ -200,6 +202,38 @@ void shammodels::sph::modules::ExternalForces::compute_ext_forc shambase::to_shared(std::move(set_eta)), shambase::to_shared(std::move(add_force_shearing_box_inertial_part))})); + } else if ( + EF_VerticalDiscPotential *ext_force + = std::get_if(&var_force.val)) { + + auto central_mass = shamrock::solvergraph::IDataEdge::make_shared("", ""); + auto R0 = 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_R0( + [r = ext_force->R0](shamrock::solvergraph::IDataEdge &R0) { + R0.data = r; // no support for offset yet + }); + set_R0.set_edges(R0); + + common::modules::AddForceVerticalDiscPotential add_force_central_grav_potential; + add_force_central_grav_potential.set_edges( + constant_G, central_mass, R0, 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_R0)), + 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 ..."); } @@ -300,6 +334,9 @@ void shammodels::sph::modules::ExternalForces::add_ext_forces() is_c_needed = true; } else if ( EF_ShearingBoxForce *ext_force = std::get_if(&var_force.val)) { + } else if ( + EF_VerticalDiscPotential *ext_force + = std::get_if(&var_force.val)) { } else { shambase::throw_unimplemented("this force is not handled, yet ..."); } @@ -451,6 +488,9 @@ void shammodels::sph::modules::ExternalForces::add_ext_forces() add_ext_forces_seq.push_back(set_q); add_ext_forces_seq.push_back(solver_graph.get_node_ptr_base(prefix_shearing_box)); + } else if ( + EF_VerticalDiscPotential *ext_force + = std::get_if(&var_force.val)) { } else { shambase::throw_unimplemented("this force is not handled, yet ..."); } diff --git a/src/shammodels/sph/src/pySPHModel.cpp b/src/shammodels/sph/src/pySPHModel.cpp index 507cade058..5c78a4c195 100644 --- a/src/shammodels/sph/src/pySPHModel.cpp +++ b/src/shammodels/sph/src/pySPHModel.cpp @@ -290,6 +290,14 @@ void add_instance(py::module &m, std::string name_config, std::string name_model py::arg("Omega_0"), py::arg("eta"), py::arg("q")) + .def( + "add_ext_force_vertical_disc_potential", + [](TConfig &self, Tscal central_mass, Tscal R0) { + self.ext_force_config.add_vertical_disc_potential(central_mass, R0); + }, + py::kw_only(), + py::arg("central_mass"), + py::arg("R0")) .def("set_units", &TConfig::set_units) .def( "get_units", From 883ef983fab3ea911396771ff6a8e9fab8f07f52 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David--Cl=C3=A9ris=20Timoth=C3=A9e?= Date: Wed, 11 Mar 2026 09:49:12 +0100 Subject: [PATCH 18/24] Apply suggestions from code review Co-authored-by: gemini-code-assist[bot] <176961590+gemini-code-assist[bot]@users.noreply.github.com> --- .../common/modules/AddForceVerticalDiscPotential.hpp | 4 ++-- src/shammodels/sph/src/modules/ExternalForces.cpp | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/shammodels/common/include/shammodels/common/modules/AddForceVerticalDiscPotential.hpp b/src/shammodels/common/include/shammodels/common/modules/AddForceVerticalDiscPotential.hpp index c1c8f709fe..59eb75ec8c 100644 --- a/src/shammodels/common/include/shammodels/common/modules/AddForceVerticalDiscPotential.hpp +++ b/src/shammodels/common/include/shammodels/common/modules/AddForceVerticalDiscPotential.hpp @@ -12,7 +12,7 @@ /** * @file AddForceVerticalDiscPotential.hpp * @author Timothée David--Cléris (tim.shamrock@proton.me) - * @brief Adds the acceleration from a central gravitational potential (point mass). + * @brief Adds the acceleration from a vertical disc potential. * */ @@ -83,7 +83,7 @@ namespace shammodels::common::modules { edges.sizes.indexes, [mGM = -cmass * G, R02 = R0 * R0](u32 gid, const Tvec *xyz, Tvec *axyz_ext) { Tscal y_a = xyz[gid].y(); - axyz_ext[gid] += mGM * y_a / sycl::sqrt(R02 + y_a * y_a); + axyz_ext[gid].y() += mGM * y_a / sycl::sqrt(R02 + y_a * y_a); }); } diff --git a/src/shammodels/sph/src/modules/ExternalForces.cpp b/src/shammodels/sph/src/modules/ExternalForces.cpp index 70392490e7..86c2d886cc 100644 --- a/src/shammodels/sph/src/modules/ExternalForces.cpp +++ b/src/shammodels/sph/src/modules/ExternalForces.cpp @@ -222,17 +222,17 @@ void shammodels::sph::modules::ExternalForces::compute_ext_forc }); set_R0.set_edges(R0); - common::modules::AddForceVerticalDiscPotential add_force_central_grav_potential; - add_force_central_grav_potential.set_edges( + common::modules::AddForceVerticalDiscPotential add_force_vertical_disc_potential; + add_force_vertical_disc_potential.set_edges( constant_G, central_mass, R0, field_xyz, sizes, field_axyz_ext); add_ext_forces_seq.push_back( std::make_shared( - "Point mass", + "Vertical disc potential", std::vector>{ shambase::to_shared(std::move(set_R0)), shambase::to_shared(std::move(set_central_mass)), - shambase::to_shared(std::move(add_force_central_grav_potential))})); + shambase::to_shared(std::move(add_force_vertical_disc_potential))})); } else { shambase::throw_unimplemented("this force is not handled, yet ..."); From 0c094298ca565e13674e62c116687bf539bde644 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Wed, 11 Mar 2026 09:49:41 +0100 Subject: [PATCH 19/24] tmp --- examples/sph/run_dusty_settle.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/sph/run_dusty_settle.py b/examples/sph/run_dusty_settle.py index d019c71e1a..068cb20e5b 100644 --- a/examples/sph/run_dusty_settle.py +++ b/examples/sph/run_dusty_settle.py @@ -183,6 +183,9 @@ def func_s(r): model.timestep() + +#TODO: add function to modify fields e.g. get rho and do stuff according to it + tnext = 0 for j in range(1000): if j > 0: From 238cc6bcf3e9dbccc55eabdec35ff369d2ef3f00 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Wed, 11 Mar 2026 10:06:37 +0100 Subject: [PATCH 20/24] cleaner --- examples/sph/run_dusty.py | 4 +- examples/sph/run_dusty_settle.py | 53 ++++++++++--------- .../modules/AddForceVerticalDiscPotential.hpp | 46 ++++++---------- src/shammodels/sph/src/Solver.cpp | 4 +- 4 files changed, 50 insertions(+), 57 deletions(-) diff --git a/examples/sph/run_dusty.py b/examples/sph/run_dusty.py index c1caf9e782..96e0a57a57 100644 --- a/examples/sph/run_dusty.py +++ b/examples/sph/run_dusty.py @@ -47,7 +47,7 @@ def uint_g(r): ndust = 1 rc = 0.25 -stopping_times = [ 1e-1] +stopping_times = [1e-1] bmin = (-0.6, -0.6, -0.6) @@ -165,6 +165,6 @@ def func_s(r): axs[3].legend() for k in range(4): - axs[k].set_xlim(-0.1,1.1) + axs[k].set_xlim(-0.1, 1.1) plt.savefig(f"mono_{j}.png") plt.close() diff --git a/examples/sph/run_dusty_settle.py b/examples/sph/run_dusty_settle.py index 068cb20e5b..0521e38f97 100644 --- a/examples/sph/run_dusty_settle.py +++ b/examples/sph/run_dusty_settle.py @@ -28,6 +28,7 @@ ucte = shamrock.Constants(codeu) G = ucte.G() + def kep_profile(r): return (G * central_mass / r) ** 0.5 @@ -35,30 +36,33 @@ def kep_profile(r): def omega_k(r): return kep_profile(r) / r + H = H_r_0 * R0 cs = H * omega_k(R0) -box = 8*H +box = 8 * H print(f"cs = {cs}") print(f"H = {H}") def scaling_rho(r): - x,y,z = r + x, y, z = r - loc_h = H/1.44 - gaussian = np.exp(-(y**2) / (2*loc_h*loc_h)) / (loc_h * np.sqrt(2 * np.pi)) + loc_h = H / (2**0.5) + gaussian = np.exp(-(y**2) / (2 * loc_h * loc_h)) / (loc_h * np.sqrt(2 * np.pi)) return gaussian + def func_rho_t(r): return rho_i * scaling_rho(r) + def func_rho_d_j(r, idust): return (0.1 / ndust) * rho_i * scaling_rho(r) def func_rho_g(r): - return rho_i*scaling_rho(r) - sum([func_rho_d_j(r, i) for i in range(ndust)]) + return rho_i * scaling_rho(r) - sum([func_rho_d_j(r, i) for i in range(ndust)]) def func_s_j(r, idust): @@ -79,11 +83,11 @@ def uint_g(r): ndust = 4 rc = 0.25 -stopping_times = np.logspace(-4,-1,ndust)*omega_k(R0) +stopping_times = np.logspace(-4, -1, ndust) * omega_k(R0) from scipy.special import erfinv -bmin = (-box/4, -box, -box/4) -bmax = (box/4, box, box/4) +bmin = (-box / 4, -box, -box / 4) +bmax = (box / 4, box, box / 4) N_target = 1e5 scheduler_split_val = int(2e7) @@ -116,7 +120,7 @@ def uint_g(r): ) cfg.set_dust_mode_monofluid_tvi(ndust) cfg.set_dust_stopping_times(stopping_times) -cfg.add_ext_force_vertical_disc_potential(central_mass = 1, R0 = 1) +cfg.add_ext_force_vertical_disc_potential(central_mass=1, R0=1) cfg.set_boundary_periodic() cfg.set_units(codeu) cfg.set_eos_isothermal(cs) @@ -153,16 +157,18 @@ def uint_g(r): model.apply_position_offset((-barycenter[0], -barycenter[1], -barycenter[2])) + def f_remap(r): - x,y,z = r + x, y, z = r + + rn = max(abs(yM), abs(ym)) + print(y, H, H * erfinv(y / rn)) + y = H * erfinv(y / rn) - rn = max(abs(yM),abs(ym)) - print(y,H, H*erfinv(y/rn) ) - y = H*erfinv(y/rn) + y = min(y, yM) + y = max(y, ym) + return (x, y, z) - y = min(y,yM) - y = max(y,ym) - return (x,y,z) model.remap_positions(f_remap) @@ -177,14 +183,13 @@ def func_s(r): model.set_field_value_lambda_f64("uint", uint_g) - model.set_cfl_cour(0.1) model.set_cfl_force(0.1) model.timestep() -#TODO: add function to modify fields e.g. get rho and do stuff according to it +# TODO: add function to modify fields e.g. get rho and do stuff according to it tnext = 0 for j in range(1000): @@ -217,21 +222,21 @@ def func_s(r): sz = 2 fig, axs = plt.subplots(nrows=1, ncols=4, figsize=(15, 5)) - axs[0].scatter(y, rho, label="rho",s=sz) + axs[0].scatter(y, rho, label="rho", s=sz) for i in range(ndust): - axs[0].scatter(y, s_j[:, i] ** 2, label=f"rho_j_{i}",s=sz) - #axs[0].scatter(y,estimated_rho) + axs[0].scatter(y, s_j[:, i] ** 2, label=f"rho_j_{i}", s=sz) + # axs[0].scatter(y,estimated_rho) axs[0].legend() for i in range(ndust): - axs[1].scatter(y, s_j[:, i] ** 2 / rho, label=f"eps_j_{i}",s=sz) + axs[1].scatter(y, s_j[:, i] ** 2 / rho, label=f"eps_j_{i}", s=sz) axs[1].legend() - axs[2].scatter(y, cs, label="soundspeed",s=sz) + axs[2].scatter(y, cs, label="soundspeed", s=sz) axs[2].legend() for i in range(ndust): - axs[3].scatter(y, ds_j_dt[:, i], label=f"ds_j_dt_{i}",s=sz) + axs[3].scatter(y, ds_j_dt[:, i], label=f"ds_j_dt_{i}", s=sz) axs[3].legend() plt.savefig(f"mono_{j}.png") diff --git a/src/shammodels/common/include/shammodels/common/modules/AddForceVerticalDiscPotential.hpp b/src/shammodels/common/include/shammodels/common/modules/AddForceVerticalDiscPotential.hpp index 59eb75ec8c..5415a3def7 100644 --- a/src/shammodels/common/include/shammodels/common/modules/AddForceVerticalDiscPotential.hpp +++ b/src/shammodels/common/include/shammodels/common/modules/AddForceVerticalDiscPotential.hpp @@ -24,6 +24,18 @@ #include "shamrock/solvergraph/Indexes.hpp" #include "shamsys/NodeInstance.hpp" +/// declare the list of edges for this node +#define NODE_EDGES(X_RO, X_RW) \ + /* inputs */ \ + X_RO(shamrock::solvergraph::IDataEdge, constant_G) \ + X_RO(shamrock::solvergraph::IDataEdge, central_mass) \ + X_RO(shamrock::solvergraph::IDataEdge, R0) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_positions) \ + X_RO(shamrock::solvergraph::Indexes, sizes) \ + \ + /* outputs */ \ + X_RW(shamrock::solvergraph::IFieldSpan, spans_accel_ext) + namespace shammodels::common::modules { template class AddForceVerticalDiscPotential : public shamrock::solvergraph::INode { @@ -33,35 +45,7 @@ namespace shammodels::common::modules { public: AddForceVerticalDiscPotential() = default; - struct Edges { - const shamrock::solvergraph::IDataEdge &constant_G; - const shamrock::solvergraph::IDataEdge ¢ral_mass; - const shamrock::solvergraph::IDataEdge &R0; - const shamrock::solvergraph::IFieldSpan &spans_positions; - const shamrock::solvergraph::Indexes &sizes; - shamrock::solvergraph::IFieldSpan &spans_accel_ext; - }; - - inline void set_edges( - std::shared_ptr> constant_G, - std::shared_ptr> central_mass, - std::shared_ptr> R0, - std::shared_ptr> spans_positions, - std::shared_ptr> sizes, - std::shared_ptr> spans_accel_ext) { - __internal_set_ro_edges({constant_G, central_mass, R0, spans_positions, sizes}); - __internal_set_rw_edges({spans_accel_ext}); - } - - inline Edges get_edges() { - return Edges{ - get_ro_edge>(0), - get_ro_edge>(1), - get_ro_edge>(2), - get_ro_edge>(3), - get_ro_edge>(4), - get_rw_edge>(0)}; - } + EXPAND_NODE_EDGES(NODE_EDGES) inline void _impl_evaluate_internal() { @@ -76,6 +60,7 @@ namespace shammodels::common::modules { Tscal G = edges.constant_G.data; Tscal R0 = edges.R0.data; + // call the GPU kernel for each patches sham::distributed_data_kernel_call( shamsys::instance::get_compute_scheduler_ptr(), sham::DDMultiRef{edges.spans_positions.get_spans()}, @@ -95,3 +80,6 @@ namespace shammodels::common::modules { }; } // namespace shammodels::common::modules + +// remove the macro before exiting +#undef NODE_EDGES diff --git a/src/shammodels/sph/src/Solver.cpp b/src/shammodels/sph/src/Solver.cpp index 40fcda89b9..b71e23c1b0 100644 --- a/src/shammodels/sph/src/Solver.cpp +++ b/src/shammodels/sph/src/Solver.cpp @@ -2515,8 +2515,8 @@ shammodels::sph::TimestepLog shammodels::sph::Solver::evolve_once() auto cfl_dt = cfl_dt_buf.get_write_access(depends_list); auto e = q.submit(depends_list, [&](sycl::handler &cgh) { - Tscal C_cour = solver_config.cfl_config.cfl_cour - * solver_config.time_state.cfl_multiplier; + Tscal C_cour = solver_config.cfl_config.cfl_cour + * solver_config.time_state.cfl_multiplier; Tscal C_force = solver_config.cfl_config.cfl_force * solver_config.time_state.cfl_multiplier; From 6d1c1dc18fd29258d524c7aa3472c06472f23ec0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Wed, 11 Mar 2026 10:55:33 +0100 Subject: [PATCH 21/24] add velocity dissipation --- examples/sph/run_dusty_settle.py | 2 + .../shammodels/common/ExtForceConfig.hpp | 28 +++++-- .../modules/AddForceVelocityDissipation.hpp | 80 +++++++++++++++++++ .../shammodels/sph/modules/ExternalForces.hpp | 1 + .../sph/src/modules/ExternalForces.cpp | 33 ++++++++ src/shammodels/sph/src/pySPHModel.cpp | 7 ++ 6 files changed, 146 insertions(+), 5 deletions(-) create mode 100644 src/shammodels/common/include/shammodels/common/modules/AddForceVelocityDissipation.hpp diff --git a/examples/sph/run_dusty_settle.py b/examples/sph/run_dusty_settle.py index 0521e38f97..fcfca5681b 100644 --- a/examples/sph/run_dusty_settle.py +++ b/examples/sph/run_dusty_settle.py @@ -121,6 +121,7 @@ def uint_g(r): cfg.set_dust_mode_monofluid_tvi(ndust) cfg.set_dust_stopping_times(stopping_times) cfg.add_ext_force_vertical_disc_potential(central_mass=1, R0=1) +cfg.add_ext_force_velocity_dissipation(eta=10) cfg.set_boundary_periodic() cfg.set_units(codeu) cfg.set_eos_isothermal(cs) @@ -240,4 +241,5 @@ def func_s(r): axs[3].legend() plt.savefig(f"mono_{j}.png") + model.do_vtk_dump(f"dump_stratif_{j}.vtk", True) plt.close() diff --git a/src/shammodels/common/include/shammodels/common/ExtForceConfig.hpp b/src/shammodels/common/include/shammodels/common/ExtForceConfig.hpp index fea3439460..30bae02dec 100644 --- a/src/shammodels/common/include/shammodels/common/ExtForceConfig.hpp +++ b/src/shammodels/common/include/shammodels/common/ExtForceConfig.hpp @@ -68,14 +68,23 @@ namespace shammodels { : shear_base(shear_base), shear_dir(shear_dir), Omega_0(Omega_0), eta(eta), q(q) {}; }; - /// -GMy / sqrt(R0^2 + y^2) + /// f = -GMy / sqrt(R0^2 + y^2) struct VerticalDiscPotential { Tscal central_mass; Tscal R0; }; - using VariantForce - = std::variant; + /// f = -eta v + struct VelocityDissipation { + Tscal eta; + }; + + using VariantForce = std::variant< + PointMass, + LenseThirring, + ShearingBoxForce, + VerticalDiscPotential, + VelocityDissipation>; VariantForce val; }; @@ -89,8 +98,7 @@ namespace shammodels { using LenseThirring = typename ExtForceVariant::LenseThirring; using ShearingBoxForce = typename ExtForceVariant::ShearingBoxForce; using VerticalDiscPotential = typename ExtForceVariant::VerticalDiscPotential; - - using VariantForce = std::variant; + using VelocityDissipation = typename ExtForceVariant::VelocityDissipation; std::vector> ext_forces; @@ -120,6 +128,10 @@ namespace shammodels { inline void add_vertical_disc_potential(Tscal central_mass, Tscal R0) { ext_forces.push_back(ExtForceVariant{VerticalDiscPotential{central_mass, R0}}); } + + inline void add_velocity_dissipation(Tscal eta) { + ext_forces.push_back(ExtForceVariant{VelocityDissipation{eta}}); + } }; } // namespace shammodels @@ -133,6 +145,7 @@ namespace shammodels { using LenseThirring = typename T::LenseThirring; using ShearingBoxForce = typename T::ShearingBoxForce; using VerticalDiscPotential = typename T::VerticalDiscPotential; + using VelocityDissipation = typename T::VelocityDissipation; if (const PointMass *v = std::get_if(&p.val)) { j = { @@ -159,6 +172,8 @@ namespace shammodels { = {{"force_type", "vertical_disc_potential"}, {"central_mass", v->central_mass}, {"R0", v->R0}}; + } else if (const VelocityDissipation *v = std::get_if(&p.val)) { + j = {{"force_type", "velocity_dissipation"}, {"eta", v->eta}}; } else { shambase::throw_unimplemented(); } @@ -180,6 +195,7 @@ namespace shammodels { using LenseThirring = typename T::LenseThirring; using ShearingBoxForce = typename T::ShearingBoxForce; using VerticalDiscPotential = typename T::VerticalDiscPotential; + using VelocityDissipation = typename T::VelocityDissipation; if (force_type == "point_mass") { p.val = PointMass{ @@ -206,6 +222,8 @@ namespace shammodels { j.at("central_mass").get(), j.at("R0").get(), }; + } else if (force_type == "velocity_dissipation") { + p.val = VelocityDissipation{j.at("eta").get()}; } else { shambase::throw_unimplemented("wtf !"); } diff --git a/src/shammodels/common/include/shammodels/common/modules/AddForceVelocityDissipation.hpp b/src/shammodels/common/include/shammodels/common/modules/AddForceVelocityDissipation.hpp new file mode 100644 index 0000000000..05106eda72 --- /dev/null +++ b/src/shammodels/common/include/shammodels/common/modules/AddForceVelocityDissipation.hpp @@ -0,0 +1,80 @@ +// -------------------------------------------------------// +// +// 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 AddForceVelocityDissipation.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief Adds the acceleration from a vertical disc potential. + * + */ + +#include "shambackends/kernel_call_distrib.hpp" +#include "shambackends/vec.hpp" +#include "shamrock/solvergraph/IDataEdge.hpp" +#include "shamrock/solvergraph/IFieldSpan.hpp" +#include "shamrock/solvergraph/INode.hpp" +#include "shamrock/solvergraph/Indexes.hpp" +#include "shamsys/NodeInstance.hpp" + +/// declare the list of edges for this node +#define NODE_EDGES(X_RO, X_RW) \ + /* inputs */ \ + X_RO(shamrock::solvergraph::IDataEdge, eta) \ + X_RO(shamrock::solvergraph::IFieldSpan, spans_velocity) \ + X_RO(shamrock::solvergraph::Indexes, sizes) \ + \ + /* outputs */ \ + X_RW(shamrock::solvergraph::IFieldSpan, spans_accel_ext) + +namespace shammodels::common::modules { + template + class AddForceVelocityDissipation : public shamrock::solvergraph::INode { + + using Tscal = shambase::VecComponent; + + public: + AddForceVelocityDissipation() = default; + + EXPAND_NODE_EDGES(NODE_EDGES) + + inline void _impl_evaluate_internal() { + + __shamrock_stack_entry(); + + auto edges = get_edges(); + + edges.spans_velocity.check_sizes(edges.sizes.indexes); + edges.spans_accel_ext.ensure_sizes(edges.sizes.indexes); + + Tscal eta = edges.eta.data; + + // call the GPU kernel for each patches + sham::distributed_data_kernel_call( + shamsys::instance::get_compute_scheduler_ptr(), + sham::DDMultiRef{edges.spans_velocity.get_spans()}, + sham::DDMultiRef{edges.spans_accel_ext.get_spans()}, + edges.sizes.indexes, + [eta](u32 gid, const Tvec *vxyz, Tvec *axyz_ext) { + axyz_ext[gid] -= eta * vxyz[gid]; + }); + } + + inline virtual std::string _impl_get_label() const { + return "AddForceVelocityDissipation"; + }; + + inline virtual std::string _impl_get_tex() const { return "TODO"; }; + }; + +} // namespace shammodels::common::modules + +// remove the macro before exiting +#undef NODE_EDGES diff --git a/src/shammodels/sph/include/shammodels/sph/modules/ExternalForces.hpp b/src/shammodels/sph/include/shammodels/sph/modules/ExternalForces.hpp index 4918f7c105..70b7ad5824 100644 --- a/src/shammodels/sph/include/shammodels/sph/modules/ExternalForces.hpp +++ b/src/shammodels/sph/include/shammodels/sph/modules/ExternalForces.hpp @@ -62,6 +62,7 @@ namespace shammodels::sph::modules { using EF_LenseThirring = typename SolverConfigExtForce::LenseThirring; using EF_ShearingBoxForce = typename SolverConfigExtForce::ShearingBoxForce; using EF_VerticalDiscPotential = typename SolverConfigExtForce::VerticalDiscPotential; + using EF_VelocityDissipation = typename SolverConfigExtForce::VelocityDissipation; inline PatchScheduler &scheduler() { return shambase::get_check_ref(context.sched); } }; diff --git a/src/shammodels/sph/src/modules/ExternalForces.cpp b/src/shammodels/sph/src/modules/ExternalForces.cpp index 86c2d886cc..30271ffa6b 100644 --- a/src/shammodels/sph/src/modules/ExternalForces.cpp +++ b/src/shammodels/sph/src/modules/ExternalForces.cpp @@ -24,6 +24,7 @@ #include "shammodels/common/modules/AddForceLenseThirring.hpp" #include "shammodels/common/modules/AddForceShearingBoxInertialPart.hpp" #include "shammodels/common/modules/AddForceShearingBoxNonInertial.hpp" +#include "shammodels/common/modules/AddForceVelocityDissipation.hpp" #include "shammodels/common/modules/AddForceVerticalDiscPotential.hpp" #include "shammodels/sph/modules/ExternalForces.hpp" #include "shammodels/sph/modules/SinkParticlesUpdate.hpp" @@ -234,6 +235,10 @@ void shammodels::sph::modules::ExternalForces::compute_ext_forc shambase::to_shared(std::move(set_central_mass)), shambase::to_shared(std::move(add_force_vertical_disc_potential))})); + } else if ( + EF_VelocityDissipation *ext_force + = std::get_if(&var_force.val)) { + } else { shambase::throw_unimplemented("this force is not handled, yet ..."); } @@ -337,6 +342,9 @@ void shammodels::sph::modules::ExternalForces::add_ext_forces() } else if ( EF_VerticalDiscPotential *ext_force = std::get_if(&var_force.val)) { + } else if ( + EF_VelocityDissipation *ext_force + = std::get_if(&var_force.val)) { } else { shambase::throw_unimplemented("this force is not handled, yet ..."); } @@ -491,6 +499,31 @@ void shammodels::sph::modules::ExternalForces::add_ext_forces() } else if ( EF_VerticalDiscPotential *ext_force = std::get_if(&var_force.val)) { + } else if ( + EF_VelocityDissipation *ext_force + = std::get_if(&var_force.val)) { + std::string prefix_eta = prefix + "eta_"; + std::string prefix_velocity_dissipation = prefix + "velocity_dissipation_"; + + auto set_eta + = register_constant_set(solver_graph, prefix_eta, [eta = ext_force->eta]() { + return eta; + }); + + auto add_force_shearing_box_non_inertial = solver_graph.register_node( + prefix_velocity_dissipation, + shammodels::common::modules::AddForceVelocityDissipation()); + shambase::get_check_ref(add_force_shearing_box_non_inertial) + .set_edges( + solver_graph.get_edge_ptr>(prefix_eta), + solver_graph.get_edge_ptr>("field_vxyz"), + solver_graph.get_edge_ptr>("field_sizes"), + solver_graph.get_edge_ptr>("field_axyz")); + + add_ext_forces_seq.push_back(set_eta); + add_ext_forces_seq.push_back( + solver_graph.get_node_ptr_base(prefix_velocity_dissipation)); + } else { shambase::throw_unimplemented("this force is not handled, yet ..."); } diff --git a/src/shammodels/sph/src/pySPHModel.cpp b/src/shammodels/sph/src/pySPHModel.cpp index ecb02a2329..051de89df2 100644 --- a/src/shammodels/sph/src/pySPHModel.cpp +++ b/src/shammodels/sph/src/pySPHModel.cpp @@ -296,6 +296,13 @@ void add_instance(py::module &m, std::string name_config, std::string name_model py::arg("Omega_0"), py::arg("eta"), py::arg("q")) + .def( + "add_ext_force_velocity_dissipation", + [](TConfig &self, Tscal eta) { + self.ext_force_config.add_velocity_dissipation(eta); + }, + py::kw_only(), + py::arg("eta")) .def( "add_ext_force_vertical_disc_potential", [](TConfig &self, Tscal central_mass, Tscal R0) { From 5f3f0dc91da00fc22bf91f0070c0dd371df5b6e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Thu, 12 Mar 2026 11:40:48 +0100 Subject: [PATCH 22/24] fix --- examples/sph/run_dusty_settle.py | 33 +++++++------- .../shammodels/common/ExtForceConfig.hpp | 29 ------------ .../sph/include/shammodels/sph/Model.hpp | 45 +++++++++++++++++-- .../sph/src/modules/ExternalForces.cpp | 6 --- src/shammodels/sph/src/pySPHModel.cpp | 5 +-- 5 files changed, 61 insertions(+), 57 deletions(-) diff --git a/examples/sph/run_dusty_settle.py b/examples/sph/run_dusty_settle.py index fcfca5681b..5e595f7163 100644 --- a/examples/sph/run_dusty_settle.py +++ b/examples/sph/run_dusty_settle.py @@ -65,13 +65,6 @@ def func_rho_g(r): return rho_i * scaling_rho(r) - sum([func_rho_d_j(r, i) for i in range(ndust)]) -def func_s_j(r, idust): - rho_t = func_rho_t(r) - rho_d_j = [func_rho_d_j(r, i) for i in range(ndust)] - eps_j = rho_d_j[idust] / rho_t - return np.sqrt(rho_t * eps_j) - - cs_g = 1 @@ -83,7 +76,7 @@ def uint_g(r): ndust = 4 rc = 0.25 -stopping_times = np.logspace(-4, -1, ndust) * omega_k(R0) +stopping_times = np.logspace(-2, -1, ndust) * omega_k(R0) from scipy.special import erfinv bmin = (-box / 4, -box, -box / 4) @@ -172,22 +165,26 @@ def f_remap(r): model.remap_positions(f_remap) +model.set_field_value_lambda_f64("uint", uint_g) + + +model.set_cfl_cour(0.1) +model.set_cfl_force(0.1) +model.timestep() -for i in range(ndust): +def compute_sj_new(patchdata): + hpart = patchdata["hpart"] + rho = pmass * (model.get_hfact() / np.array(hpart)) ** 3 - def func_s(r): - return func_s_j(r, i) + epsilon_target = 0.1 / ndust + s = np.sqrt(rho * epsilon_target) - model.set_field_value_lambda_f64("s_j", func_s, i) + return s -model.set_field_value_lambda_f64("uint", uint_g) -model.set_cfl_cour(0.1) -model.set_cfl_force(0.1) -model.timestep() # TODO: add function to modify fields e.g. get rho and do stuff according to it @@ -198,6 +195,10 @@ def func_s(r): tnext += 0.02 model.evolve_until(tnext) + if(j == 70): + for k in range(ndust): + model.overwrite_field_value_f64("s_j",compute_sj_new,k) + dic = ctx.collect_data() print(dic["s_j"]) diff --git a/src/shammodels/common/include/shammodels/common/ExtForceConfig.hpp b/src/shammodels/common/include/shammodels/common/ExtForceConfig.hpp index e19b2a40a2..30bae02dec 100644 --- a/src/shammodels/common/include/shammodels/common/ExtForceConfig.hpp +++ b/src/shammodels/common/include/shammodels/common/ExtForceConfig.hpp @@ -68,17 +68,12 @@ namespace shammodels { : shear_base(shear_base), shear_dir(shear_dir), Omega_0(Omega_0), eta(eta), q(q) {}; }; -<<<<<<< patch-2026-02-05-15-21 /// f = -GMy / sqrt(R0^2 + y^2) -======= - /// -GMy / sqrt(R0^2 + y^2) ->>>>>>> main struct VerticalDiscPotential { Tscal central_mass; Tscal R0; }; -<<<<<<< patch-2026-02-05-15-21 /// f = -eta v struct VelocityDissipation { Tscal eta; @@ -90,10 +85,6 @@ namespace shammodels { ShearingBoxForce, VerticalDiscPotential, VelocityDissipation>; -======= - using VariantForce - = std::variant; ->>>>>>> main VariantForce val; }; @@ -107,12 +98,7 @@ namespace shammodels { using LenseThirring = typename ExtForceVariant::LenseThirring; using ShearingBoxForce = typename ExtForceVariant::ShearingBoxForce; using VerticalDiscPotential = typename ExtForceVariant::VerticalDiscPotential; -<<<<<<< patch-2026-02-05-15-21 using VelocityDissipation = typename ExtForceVariant::VelocityDissipation; -======= - - using VariantForce = std::variant; ->>>>>>> main std::vector> ext_forces; @@ -142,13 +128,10 @@ namespace shammodels { inline void add_vertical_disc_potential(Tscal central_mass, Tscal R0) { ext_forces.push_back(ExtForceVariant{VerticalDiscPotential{central_mass, R0}}); } -<<<<<<< patch-2026-02-05-15-21 inline void add_velocity_dissipation(Tscal eta) { ext_forces.push_back(ExtForceVariant{VelocityDissipation{eta}}); } -======= ->>>>>>> main }; } // namespace shammodels @@ -162,10 +145,7 @@ namespace shammodels { using LenseThirring = typename T::LenseThirring; using ShearingBoxForce = typename T::ShearingBoxForce; using VerticalDiscPotential = typename T::VerticalDiscPotential; -<<<<<<< patch-2026-02-05-15-21 using VelocityDissipation = typename T::VelocityDissipation; -======= ->>>>>>> main if (const PointMass *v = std::get_if(&p.val)) { j = { @@ -192,11 +172,8 @@ namespace shammodels { = {{"force_type", "vertical_disc_potential"}, {"central_mass", v->central_mass}, {"R0", v->R0}}; -<<<<<<< patch-2026-02-05-15-21 } else if (const VelocityDissipation *v = std::get_if(&p.val)) { j = {{"force_type", "velocity_dissipation"}, {"eta", v->eta}}; -======= ->>>>>>> main } else { shambase::throw_unimplemented(); } @@ -218,10 +195,7 @@ namespace shammodels { using LenseThirring = typename T::LenseThirring; using ShearingBoxForce = typename T::ShearingBoxForce; using VerticalDiscPotential = typename T::VerticalDiscPotential; -<<<<<<< patch-2026-02-05-15-21 using VelocityDissipation = typename T::VelocityDissipation; -======= ->>>>>>> main if (force_type == "point_mass") { p.val = PointMass{ @@ -248,11 +222,8 @@ namespace shammodels { j.at("central_mass").get(), j.at("R0").get(), }; -<<<<<<< patch-2026-02-05-15-21 } else if (force_type == "velocity_dissipation") { p.val = VelocityDissipation{j.at("eta").get()}; -======= ->>>>>>> main } else { shambase::throw_unimplemented("wtf !"); } diff --git a/src/shammodels/sph/include/shammodels/sph/Model.hpp b/src/shammodels/sph/include/shammodels/sph/Model.hpp index dbbe364ce6..6e43a7bb57 100644 --- a/src/shammodels/sph/include/shammodels/sph/Model.hpp +++ b/src/shammodels/sph/include/shammodels/sph/Model.hpp @@ -31,6 +31,7 @@ #include "shammodels/sph/math/density.hpp" #include "shammodels/sph/modules/ComputeLoadBalanceValue.hpp" #include "shammodels/sph/modules/SPHSetup.hpp" +#include "shampylib/PatchDataToPy.hpp" #include "shamrock/io/ShamrockDump.hpp" #include "shamrock/patch/PatchDataLayer.hpp" #include "shamrock/scheduler/ReattributeDataUtility.hpp" @@ -215,6 +216,44 @@ namespace shammodels::sph { }); } + template + inline void overwrite_field_value( + std::string field_name, + const std::function(py::dict)> field_compute, + const i32 offset) { + + StackEntry stack_loc{}; + + PatchScheduler &sched = shambase::get_check_ref(ctx.sched); + + u32 ifield = sched.pdl_old().get_field_idx(field_name); + + sched.patch_data.for_each_patchdata( + [&](u64 patch_id, shamrock::patch::PatchDataLayer &pdat) { + PatchDataField &f = pdat.template get_field(ifield); + + auto result = field_compute(shamrock::pdat_to_dic(pdat)); + + if (result.size() != f.get_obj_cnt()) { + throw shambase::make_except_with_loc(shambase::format( + "result.size() != f.get_obj_cnt() ({} != {})", + result.size(), + f.get_obj_cnt())); + } + + auto f_nvar = f.get_nvar(); + auto acc = f.get_buf().copy_to_stdvec(); + + for (u32 i = 0; i < pdat.get_obj_cnt(); i++) { + acc[i * f_nvar + offset] = result[i]; + } + + logger::raw_ln(acc); + + f.get_buf().copy_from_stdvec(acc); + }); + } + /** * @brief Add a disc distribution * @@ -942,9 +981,9 @@ namespace shammodels::sph { * (1. + sycl::sin(shambase::constants::pi / (2. * Hwarp) * (R - Rwarp))) * sycl::sin(incl_rad)); - psi = shambase::constants::pi - * Rwarp / (4. * Hwarp) * sycl::sin(incl_rad) - / sycl::sqrt(1. - (0.5 * sycl::pow(sycl::sin(incl_rad), 2))); + psi = shambase::constants::pi + * Rwarp / (4. * Hwarp) * sycl::sin(incl_rad) + / sycl::sqrt(1. - (0.5 * sycl::pow(sycl::sin(incl_rad), 2))); Tscal psimax = sycl::max(psimax, psi); Tscal x = pos[i].x(); Tscal y = pos[i].y(); diff --git a/src/shammodels/sph/src/modules/ExternalForces.cpp b/src/shammodels/sph/src/modules/ExternalForces.cpp index deacbd5b93..30271ffa6b 100644 --- a/src/shammodels/sph/src/modules/ExternalForces.cpp +++ b/src/shammodels/sph/src/modules/ExternalForces.cpp @@ -24,10 +24,7 @@ #include "shammodels/common/modules/AddForceLenseThirring.hpp" #include "shammodels/common/modules/AddForceShearingBoxInertialPart.hpp" #include "shammodels/common/modules/AddForceShearingBoxNonInertial.hpp" -<<<<<<< patch-2026-02-05-15-21 #include "shammodels/common/modules/AddForceVelocityDissipation.hpp" -======= ->>>>>>> main #include "shammodels/common/modules/AddForceVerticalDiscPotential.hpp" #include "shammodels/sph/modules/ExternalForces.hpp" #include "shammodels/sph/modules/SinkParticlesUpdate.hpp" @@ -238,13 +235,10 @@ void shammodels::sph::modules::ExternalForces::compute_ext_forc shambase::to_shared(std::move(set_central_mass)), shambase::to_shared(std::move(add_force_vertical_disc_potential))})); -<<<<<<< patch-2026-02-05-15-21 } else if ( EF_VelocityDissipation *ext_force = std::get_if(&var_force.val)) { -======= ->>>>>>> main } else { shambase::throw_unimplemented("this force is not handled, yet ..."); } diff --git a/src/shammodels/sph/src/pySPHModel.cpp b/src/shammodels/sph/src/pySPHModel.cpp index 09366c012c..0a178c190f 100644 --- a/src/shammodels/sph/src/pySPHModel.cpp +++ b/src/shammodels/sph/src/pySPHModel.cpp @@ -297,7 +297,6 @@ void add_instance(py::module &m, std::string name_config, std::string name_model py::arg("eta"), py::arg("q")) .def( -<<<<<<< patch-2026-02-05-15-21 "add_ext_force_velocity_dissipation", [](TConfig &self, Tscal eta) { self.ext_force_config.add_velocity_dissipation(eta); @@ -305,8 +304,6 @@ void add_instance(py::module &m, std::string name_config, std::string name_model py::kw_only(), py::arg("eta")) .def( -======= ->>>>>>> main "add_ext_force_vertical_disc_potential", [](TConfig &self, Tscal central_mass, Tscal R0) { self.ext_force_config.add_vertical_disc_potential(central_mass, R0); @@ -800,6 +797,8 @@ void add_instance(py::module &m, std::string name_config, std::string name_model py::arg("field_name"), py::arg("pos_to_val"), py::arg("offset") = 0) + .def("overwrite_field_value_f64", &T::template overwrite_field_value) + .def("overwrite_field_value_f64_3", &T::template overwrite_field_value) .def("remap_positions", &T::remap_positions) //.def("set_field_value_lambda_f64_3",[](T&self,std::string field_name, const // std::function pos_to_val){ From 2018540d4fce46a8602cfd706e69eec63de3d64d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Thu, 12 Mar 2026 20:40:08 +0100 Subject: [PATCH 23/24] better --- examples/sph/run_dusty_settle.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/examples/sph/run_dusty_settle.py b/examples/sph/run_dusty_settle.py index 5e595f7163..f775ece302 100644 --- a/examples/sph/run_dusty_settle.py +++ b/examples/sph/run_dusty_settle.py @@ -76,13 +76,14 @@ def uint_g(r): ndust = 4 rc = 0.25 -stopping_times = np.logspace(-2, -1, ndust) * omega_k(R0) +stopping_times = np.logspace(-3, 0, ndust) * omega_k(R0) +print(stopping_times, omega_k(R0)) from scipy.special import erfinv bmin = (-box / 4, -box, -box / 4) bmax = (box / 4, box, box / 4) -N_target = 1e5 +N_target = 1e4 scheduler_split_val = int(2e7) scheduler_merge_val = int(1) @@ -195,7 +196,7 @@ def compute_sj_new(patchdata): tnext += 0.02 model.evolve_until(tnext) - if(j == 70): + if(j == 100): for k in range(ndust): model.overwrite_field_value_f64("s_j",compute_sj_new,k) @@ -229,6 +230,7 @@ def compute_sj_new(patchdata): axs[0].scatter(y, s_j[:, i] ** 2, label=f"rho_j_{i}", s=sz) # axs[0].scatter(y,estimated_rho) + axs[0].set_yscale('log') axs[0].legend() for i in range(ndust): axs[1].scatter(y, s_j[:, i] ** 2 / rho, label=f"eps_j_{i}", s=sz) From d6745fb1098aefa5c0b6ba0807840b9e9e610a27 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Fri, 13 Mar 2026 14:59:48 +0100 Subject: [PATCH 24/24] update setup --- examples/sph/run_dusty_settle.py | 35 ++++++++++++++++++++------------ 1 file changed, 22 insertions(+), 13 deletions(-) diff --git a/examples/sph/run_dusty_settle.py b/examples/sph/run_dusty_settle.py index f775ece302..e92741c997 100644 --- a/examples/sph/run_dusty_settle.py +++ b/examples/sph/run_dusty_settle.py @@ -76,14 +76,15 @@ def uint_g(r): ndust = 4 rc = 0.25 -stopping_times = np.logspace(-3, 0, ndust) * omega_k(R0) +stokes = np.logspace(-3, 0, ndust) +stopping_times = stokes / omega_k(R0) print(stopping_times, omega_k(R0)) from scipy.special import erfinv bmin = (-box / 4, -box, -box / 4) bmax = (box / 4, box, box / 4) -N_target = 1e4 +N_target = 1e5 scheduler_split_val = int(2e7) scheduler_merge_val = int(1) @@ -193,10 +194,10 @@ def compute_sj_new(patchdata): tnext = 0 for j in range(1000): if j > 0: - tnext += 0.02 + tnext += 0.1 model.evolve_until(tnext) - if(j == 100): + if(j == 20): for k in range(ndust): model.overwrite_field_value_f64("s_j",compute_sj_new,k) @@ -224,24 +225,32 @@ def compute_sj_new(patchdata): sz = 2 - fig, axs = plt.subplots(nrows=1, ncols=4, figsize=(15, 5)) - axs[0].scatter(y, rho, label="rho", s=sz) + fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(12, 5)) + time = model.get_time() + fig.suptitle(f"t = {time:.2f}") + + fig.subplots_adjust(left=0.07, right=0.98, wspace=0.4) + + axs[0].scatter(y, rho, label="gas", s=sz) for i in range(ndust): - axs[0].scatter(y, s_j[:, i] ** 2, label=f"rho_j_{i}", s=sz) + axs[0].scatter(y, s_j[:, i] ** 2, label=f"dust {i} ts = {stopping_times[i]:.2f}", s=sz) # axs[0].scatter(y,estimated_rho) + axs[0].set_ylabel(r"$\rho$") + axs[0].set_xlabel(r"$y$") axs[0].set_yscale('log') axs[0].legend() for i in range(ndust): - axs[1].scatter(y, s_j[:, i] ** 2 / rho, label=f"eps_j_{i}", s=sz) + axs[1].scatter(y, s_j[:, i] ** 2 / rho, label=f"dust {i} ts = {stopping_times[i]:.2f}", s=sz) + axs[1].set_ylabel(r"$\epsilon_j$") + axs[1].set_xlabel(r"$y$") axs[1].legend() - axs[2].scatter(y, cs, label="soundspeed", s=sz) - axs[2].legend() - for i in range(ndust): - axs[3].scatter(y, ds_j_dt[:, i], label=f"ds_j_dt_{i}", s=sz) - axs[3].legend() + axs[2].scatter(y, ds_j_dt[:, i], label=f"dust {i} ts = {stopping_times[i]:.2f}", s=sz) + axs[2].set_ylabel(r"$\frac{d s_j}{dt}$") + axs[2].set_xlabel(r"$y$") + axs[2].legend() plt.savefig(f"mono_{j}.png") model.do_vtk_dump(f"dump_stratif_{j}.vtk", True)