From 9dc5b6d0fb0298655f66aa9575314e8ad480cc73 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Wed, 1 Apr 2026 12:20:45 +0200 Subject: [PATCH 01/23] [SPH] modernize sink accretion routine --- .../sph/src/modules/SinkParticlesUpdate.cpp | 75 +++++++++---------- 1 file changed, 35 insertions(+), 40 deletions(-) diff --git a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp index e72b6c095c..8e0e9445fc 100644 --- a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp +++ b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp @@ -15,10 +15,14 @@ * */ -#include "shammodels/sph/modules/SinkParticlesUpdate.hpp" +#include "shambase/narrowing.hpp" +#include "shamalgs/details/numeric/numeric.hpp" #include "shamalgs/numeric.hpp" #include "shamalgs/reduction.hpp" +#include "shambackends/DeviceBuffer.hpp" +#include "shambackends/kernel_call.hpp" #include "shammath/sphkernels.hpp" +#include "shammodels/sph/modules/SinkParticlesUpdate.hpp" #include "shamsys/legacy/log.hpp" template class SPHKernel> @@ -57,67 +61,58 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par sham::DeviceBuffer &buf_xyz = pdat.get_field_buf_ref(ixyz); sham::DeviceBuffer &buf_vxyz = pdat.get_field_buf_ref(ivxyz); - sycl::buffer not_accreted(Nobj); - sycl::buffer accreted(Nobj); - - sham::EventList depends_list1; - auto xyz = buf_xyz.get_read_access(depends_list1); - - auto e1 = q.submit(depends_list1, [&](sycl::handler &cgh) { - sycl::accessor not_acc{not_accreted, cgh, sycl::write_only, sycl::no_init}; - sycl::accessor acc{accreted, cgh, sycl::write_only, sycl::no_init}; + sham::DeviceBuffer not_accreted(Nobj, dev_sched); + sham::DeviceBuffer accreted(Nobj, dev_sched); - Tvec r_sink = s.pos; - Tscal acc_rad2 = s.accretion_radius * s.accretion_radius; + Tvec r_sink = s.pos; + Tscal acc_rad2 = s.accretion_radius * s.accretion_radius; - shambase::parallel_for(cgh, Nobj, "check accretion", [=](i32 id_a) { + sham::kernel_call( + q, + sham::MultiRef{buf_xyz}, + sham::MultiRef{not_accreted, accreted}, + Nobj, + [r_sink, acc_rad2](u32 id_a, const Tvec *__restrict xyz, u32 *not_acc, u32 *acc) { Tvec r = xyz[id_a] - r_sink; bool not_accreted = sycl::dot(r, r) > acc_rad2; not_acc[id_a] = (not_accreted) ? 1 : 0; acc[id_a] = (!not_accreted) ? 1 : 0; }); - }); - - buf_xyz.complete_event_state(e1); - - std::tuple>, u32> id_list_keep - = shamalgs::numeric::stream_compact(q.q, not_accreted, Nobj); - std::tuple>, u32> id_list_accrete - = shamalgs::numeric::stream_compact(q.q, accreted, Nobj); + sham::DeviceBuffer id_list_keep + = shamalgs::stream_compact(dev_sched, not_accreted, Nobj); + sham::DeviceBuffer id_list_accrete + = shamalgs::stream_compact(dev_sched, accreted, Nobj); // sum accreted values onto sink + if (id_list_accrete.get_size() > 0) { - if (std::get<1>(id_list_accrete) > 0) { - - u32 Naccrete = std::get<1>(id_list_accrete); + u32 Naccrete = shambase::narrow_or_throw(id_list_accrete.get_size()); Tscal acc_mass = gpart_mass * Naccrete; sham::DeviceBuffer pxyz_acc(Naccrete, dev_sched); - sham::EventList depends_list2; - auto vxyz = buf_vxyz.get_read_access(depends_list2); - auto accretion_p = pxyz_acc.get_write_access(depends_list2); - - auto e = q.submit(depends_list2, [&, gpart_mass](sycl::handler &cgh) { - sycl::accessor id_acc{*std::get<0>(id_list_accrete), cgh, sycl::read_only}; - - shambase::parallel_for( - cgh, Naccrete, "compute sum momentum accretion", [=](i32 id_a) { - accretion_p[id_a] = gpart_mass * vxyz[id_acc[id_a]]; - }); - }); - - buf_vxyz.complete_event_state(e); - pxyz_acc.complete_event_state(e); + sham::kernel_call( + q, + sham::MultiRef{buf_vxyz, id_list_accrete}, + sham::MultiRef{pxyz_acc}, + Naccrete, + [gpart_mass]( + u32 id_a, + const Tvec *__restrict vxyz, + const u32 *__restrict id_acc, + Tvec *accretion_p) { + accretion_p[id_a] = gpart_mass * vxyz[id_acc[id_a]]; + }); Tvec acc_pxyz = shamalgs::primitives::sum(dev_sched, pxyz_acc, 0, Naccrete); s_acc_mass += acc_mass; s_acc_pxyz += acc_pxyz; - pdat.keep_ids(*std::get<0>(id_list_keep), std::get<1>(id_list_keep)); + pdat.keep_ids( + id_list_keep, shambase::narrow_or_throw(id_list_keep.get_size())); } }); From b0a59e733016e0019d2c2a6fbdfc9f7916c01309 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Thu, 2 Apr 2026 00:46:23 +0200 Subject: [PATCH 02/23] start implementing everything --- src/shammath/include/shammath/matrix_op.hpp | 29 ++ .../include/shammodels/sph/SinkPartStruct.hpp | 4 + .../sph/src/modules/SinkParticlesUpdate.cpp | 271 ++++++++++++++++-- 3 files changed, 284 insertions(+), 20 deletions(-) diff --git a/src/shammath/include/shammath/matrix_op.hpp b/src/shammath/include/shammath/matrix_op.hpp index dde113f365..1d22ffeeeb 100644 --- a/src/shammath/include/shammath/matrix_op.hpp +++ b/src/shammath/include/shammath/matrix_op.hpp @@ -334,6 +334,35 @@ namespace shammath { } } + /** + * @brief Compute the determinant of a 3x3 matrix. + * + * @param input The input matrix to invert. + * @return the determinant + */ + template + inline T mat_det_33( + const std::mdspan, Layout, Accessor> &input) { + + T &a00 = input(0, 0); + T &a10 = input(1, 0); + T &a20 = input(2, 0); + + T &a01 = input(0, 1); + T &a11 = input(1, 1); + T &a21 = input(2, 1); + + T &a02 = input(0, 2); + T &a12 = input(1, 2); + T &a22 = input(2, 2); + + T det + = (-a02 * a11 * a20 + a01 * a12 * a20 + a02 * a10 * a21 - a00 * a12 * a21 + - a01 * a10 * a22 + a00 * a11 * a22); + + return det; + } + /** * @brief Compute the inverse of a 3x3 matrix. * diff --git a/src/shammodels/sph/include/shammodels/sph/SinkPartStruct.hpp b/src/shammodels/sph/include/shammodels/sph/SinkPartStruct.hpp index 0a58f2232d..0c64835d79 100644 --- a/src/shammodels/sph/include/shammodels/sph/SinkPartStruct.hpp +++ b/src/shammodels/sph/include/shammodels/sph/SinkPartStruct.hpp @@ -34,6 +34,8 @@ namespace shammodels::sph { Tscal mass; Tvec angular_momentum; Tscal accretion_radius; + + bool is_torque_free = false; }; template @@ -50,6 +52,7 @@ namespace shammodels::sph { {"mass", p.mass}, {"angular_momentum", p.angular_momentum}, {"accretion_radius", p.accretion_radius}, + {"is_torque_free", p.is_torque_free}, }; } @@ -65,6 +68,7 @@ namespace shammodels::sph { j.at("mass").get_to(p.mass); j.at("angular_momentum").get_to(p.angular_momentum); j.at("accretion_radius").get_to(p.accretion_radius); + j.at("is_torque_free").get_to(p.is_torque_free); } } // namespace shammodels::sph diff --git a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp index 8e0e9445fc..2e1bf44399 100644 --- a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp +++ b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp @@ -15,15 +15,24 @@ * */ +#include "shambase/DistributedData.hpp" #include "shambase/narrowing.hpp" +#include "shamalgs/collective/reduction.hpp" #include "shamalgs/details/numeric/numeric.hpp" #include "shamalgs/numeric.hpp" +#include "shamalgs/primitives/reduction.hpp" #include "shamalgs/reduction.hpp" #include "shambackends/DeviceBuffer.hpp" #include "shambackends/kernel_call.hpp" +#include "shambackends/math.hpp" +#include "shamcomm/logs.hpp" +#include "shammath/matrix.hpp" +#include "shammath/matrix_op.hpp" #include "shammath/sphkernels.hpp" #include "shammodels/sph/modules/SinkParticlesUpdate.hpp" #include "shamsys/legacy/log.hpp" +#include +#include template class SPHKernel> void shammodels::sph::modules::SinkParticlesUpdate::accrete_particles() { @@ -50,10 +59,20 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par u32 sink_id = 0; bool had_accretion = false; std::string log = "sink accretion :"; + + struct AccretionFlagBufs { + sham::DeviceBuffer not_accreted; + sham::DeviceBuffer accreted; + }; + for (Sink &s : sink_parts) { - Tscal s_acc_mass = 0; - Tvec s_acc_pxyz = {0, 0, 0}; + Tvec r_sink = s.pos; + Tvec v_sink = s.velocity; + Tscal acc_rad2 = s.accretion_radius * s.accretion_radius; + + // flags particles for accretion + shambase::DistributedData accretion_flag_bufs{}; scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { u32 Nobj = pdat.get_obj_cnt(); @@ -64,67 +83,279 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par sham::DeviceBuffer not_accreted(Nobj, dev_sched); sham::DeviceBuffer accreted(Nobj, dev_sched); - Tvec r_sink = s.pos; - Tscal acc_rad2 = s.accretion_radius * s.accretion_radius; - sham::kernel_call( q, sham::MultiRef{buf_xyz}, sham::MultiRef{not_accreted, accreted}, Nobj, - [r_sink, acc_rad2](u32 id_a, const Tvec *__restrict xyz, u32 *not_acc, u32 *acc) { + [r_sink, acc_rad2]( + u32 id_a, + const Tvec *__restrict xyz, + u32 *__restrict not_acc, + u32 *__restrict acc) { Tvec r = xyz[id_a] - r_sink; bool not_accreted = sycl::dot(r, r) > acc_rad2; not_acc[id_a] = (not_accreted) ? 1 : 0; acc[id_a] = (!not_accreted) ? 1 : 0; }); - sham::DeviceBuffer id_list_keep - = shamalgs::stream_compact(dev_sched, not_accreted, Nobj); + accretion_flag_bufs.add_obj( + cur_p.id_patch, AccretionFlagBufs{std::move(not_accreted), std::move(accreted)}); + }); + + // list the ids that will be accreted + shambase::DistributedData> bufs_id_list_accrete{}; + + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + u32 Nobj = pdat.get_obj_cnt(); + + sham::DeviceBuffer &accreted = accretion_flag_bufs.get(cur_p.id_patch).accreted; + sham::DeviceBuffer id_list_accrete = shamalgs::stream_compact(dev_sched, accreted, Nobj); + bufs_id_list_accrete.add_obj(cur_p.id_patch, std::move(id_list_accrete)); + }); + + // compute the accreted mass, position moment and linear momentum + Tscal s_acc_mass = 0; + Tvec s_acc_mxyz = {0, 0, 0}; + Tvec s_acc_pxyz = {0, 0, 0}; + + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + u32 Nobj = pdat.get_obj_cnt(); + + sham::DeviceBuffer &buf_xyz = pdat.get_field_buf_ref(ixyz); + sham::DeviceBuffer &buf_vxyz = pdat.get_field_buf_ref(ivxyz); + + sham::DeviceBuffer ¬_accreted + = accretion_flag_bufs.get(cur_p.id_patch).not_accreted; + sham::DeviceBuffer &accreted = accretion_flag_bufs.get(cur_p.id_patch).accreted; + + sham::DeviceBuffer &id_list_accrete = bufs_id_list_accrete.get(cur_p.id_patch); + // sum accreted values onto sink if (id_list_accrete.get_size() > 0) { - u32 Naccrete = shambase::narrow_or_throw(id_list_accrete.get_size()); Tscal acc_mass = gpart_mass * Naccrete; sham::DeviceBuffer pxyz_acc(Naccrete, dev_sched); + sham::DeviceBuffer mxyz_acc(Naccrete, dev_sched); sham::kernel_call( q, - sham::MultiRef{buf_vxyz, id_list_accrete}, - sham::MultiRef{pxyz_acc}, + sham::MultiRef{buf_xyz, buf_vxyz, id_list_accrete}, + sham::MultiRef{pxyz_acc, mxyz_acc}, Naccrete, - [gpart_mass]( + [gpart_mass, r_sink, v_sink]( u32 id_a, + const Tvec *__restrict xyz, const Tvec *__restrict vxyz, const u32 *__restrict id_acc, - Tvec *accretion_p) { - accretion_p[id_a] = gpart_mass * vxyz[id_acc[id_a]]; + Tvec *accretion_p, + Tvec *accretion_mr) { + u32 i_a = id_acc[id_a]; + Tvec r = xyz[i_a]; + Tvec v = vxyz[i_a]; + accretion_p[id_a] = gpart_mass * v; + accretion_mr[id_a] = gpart_mass * r; }); Tvec acc_pxyz = shamalgs::primitives::sum(dev_sched, pxyz_acc, 0, Naccrete); + Tvec acc_mxyz = shamalgs::primitives::sum(dev_sched, mxyz_acc, 0, Naccrete); s_acc_mass += acc_mass; s_acc_pxyz += acc_pxyz; - - pdat.keep_ids( - id_list_keep, shambase::narrow_or_throw(id_list_keep.get_size())); + s_acc_mxyz += acc_mxyz; } }); Tscal sum_acc_mass = shamalgs::collective::allreduce_sum(s_acc_mass); Tvec sum_acc_pxyz = shamalgs::collective::allreduce_sum(s_acc_pxyz); + Tvec sum_acc_mxyz = shamalgs::collective::allreduce_sum(s_acc_mxyz); + + // compute the new sink values + Tscal new_mass = s.mass + sum_acc_mass; + Tvec new_pos = (sum_acc_mxyz + s.pos * s.mass) / (s.mass + sum_acc_mass); + Tvec new_vel = (sum_acc_pxyz + s.velocity * s.mass) / (s.mass + sum_acc_mass); + + // use them to compute the angular_momentum updated + Tvec s_acc_lxyz = {0, 0, 0}; + Tvec sum_acc_lxyz = shamalgs::collective::allreduce_sum(s_acc_lxyz); + + Tvec new_ang_mom = s.angular_momentum + sum_acc_lxyz; + + // write back the updated sink state + auto new_state = s; + new_state.mass = new_mass; + new_state.pos = new_pos; + new_state.velocity = new_vel; + new_state.angular_momentum = new_ang_mom; + + s = new_state; - s.mass += sum_acc_mass; - s.velocity += sum_acc_pxyz / s.mass; if (sum_acc_mass > 0) { had_accretion = true; log += shambase::format( - "\n id {} mass {} vel {}", sink_id, sum_acc_mass, sum_acc_pxyz / s.mass); + "\n id {} mass {} vel {} l {}", + sink_id, + sum_acc_mass, + sum_acc_pxyz / s.mass, + sum_acc_lxyz); + } + + // evict accreted particles from patches + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + u32 Nobj = pdat.get_obj_cnt(); + + sham::DeviceBuffer ¬_accreted + = accretion_flag_bufs.get(cur_p.id_patch).not_accreted; + sham::DeviceBuffer &accreted = accretion_flag_bufs.get(cur_p.id_patch).accreted; + + sham::DeviceBuffer &id_list_accrete = bufs_id_list_accrete.get(cur_p.id_patch); + + if (id_list_accrete.get_size() > 0) { + + sham::DeviceBuffer id_list_keep + = shamalgs::stream_compact(dev_sched, not_accreted, Nobj); + + pdat.keep_ids( + id_list_keep, shambase::narrow_or_throw(id_list_keep.get_size())); + } + }); + + if (s.is_torque_free && sham::dot(s.angular_momentum, s.angular_momentum) > 0) { + // boost particles in a radius of 4 racc to empty the sink angular_momentum + + struct BoostWeight { + Tscal radius; + Tscal weight(Tscal r) const { + if (r < radius) { + return 1; + } else { + return 0; + } + } + }; + + BoostWeight weight_func{s.accretion_radius * 10}; + + Tvec r_sink = s.pos; + + using mat = shammath::mat; + + mat I_rank{}; + + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + u32 Nobj = pdat.get_obj_cnt(); + + sham::DeviceBuffer &buf_xyz = pdat.get_field_buf_ref(ixyz); + + sham::DeviceBuffer weighted_intertia(9 * Nobj, dev_sched); + sham::kernel_call( + q, + sham::MultiRef{buf_xyz}, + sham::MultiRef{weighted_intertia}, + Nobj, + [r_sink, gpart_mass, Nobj, weight_func]( + u32 id_a, const Tvec *__restrict xyz, Tscal *weighted_intertia) { + Tvec r_a = xyz[id_a] - r_sink; + Tscal r_a2 = sham::dot(r_a, r_a); + + Tscal w = gpart_mass * weight_func.weight(sycl::sqrt(r_a2)); + + mat I{}; + + { // I = w [Id r_a^2 - r_a r_a^T] + auto Imd = I.get_mdspan(); + + Imd(0, 0) += w * (r_a.y() * r_a.y() + r_a.z() * r_a.z()); + Imd(1, 1) += w * (r_a.x() * r_a.x() + r_a.z() * r_a.z()); + Imd(2, 2) += w * (r_a.x() * r_a.x() + r_a.y() * r_a.y()); + + Imd(0, 1) -= w * r_a.x() * r_a.y(); + Imd(0, 2) -= w * r_a.x() * r_a.z(); + Imd(1, 2) -= w * r_a.y() * r_a.z(); + + // symmetry time ! + Imd(1, 0) = Imd(0, 1); + Imd(2, 0) = Imd(0, 2); + Imd(2, 1) = Imd(1, 2); + } + +#pragma unroll + for (u32 i = 0; i < 9; i++) { + weighted_intertia[Nobj * i + id_a] = I.data[i]; + } + }); + + for (u32 i = 0; i < 9; i++) { + I_rank.data[i] += shamalgs::primitives::sum( + dev_sched, weighted_intertia, Nobj * i, Nobj * (i + 1)); + } + }); + + mat weighted_intertia_sum{}; + + for (u32 i = 0; i < 9; i++) { + weighted_intertia_sum.data[i] = shamalgs::collective::allreduce_sum(I_rank.data[i]); + } + + mat I_inv{}; + if (shammath::mat_det_33(weighted_intertia_sum.get_mdspan()) != 0) { + shammath::mat_inv_33(weighted_intertia_sum.get_mdspan(), I_inv.get_mdspan()); + } + + shammath::vec delta_l{}; + delta_l.data[0] = s.angular_momentum.x(); + delta_l.data[1] = s.angular_momentum.y(); + delta_l.data[2] = s.angular_momentum.z(); + + shammath::vec delta_w{}; + shammath::mat_prod( + I_inv.get_mdspan(), delta_l.get_mdspan_mat_col(), delta_w.get_mdspan_mat_col()); + + Tvec vec_delta_w = {delta_w.data[0], delta_w.data[1], delta_w[2]}; + + logger::raw_ln("delta l =", delta_l.data[0], delta_l.data[1], delta_l.data[2]); + logger::raw_ln("delta w =", delta_w.data[0], delta_w.data[1], delta_w.data[2]); + + Tvec lboost_sum_rank{}; + scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { + u32 Nobj = pdat.get_obj_cnt(); + + sham::DeviceBuffer &buf_xyz = pdat.get_field_buf_ref(ixyz); + sham::DeviceBuffer &buf_vxyz = pdat.get_field_buf_ref(ivxyz); + + sham::DeviceBuffer l_boost(Nobj, dev_sched); + + sham::kernel_call( + q, + sham::MultiRef{buf_xyz}, + sham::MultiRef{buf_vxyz, l_boost}, + Nobj, + [r_sink, delta_w = vec_delta_w, weight_func, gpart_mass]( + u32 id_a, const Tvec *__restrict xyz, Tvec *vxyz, Tvec *lboost) { + Tvec r_a = xyz[id_a] - r_sink; + Tscal r_a2 = sham::dot(r_a, r_a); + + Tscal W = weight_func.weight(sycl::sqrt(r_a2)); + + Tvec v_boost = W * sycl::cross(delta_w, r_a); + + vxyz[id_a] += v_boost; + lboost[id_a] += gpart_mass * sycl::cross(r_a, v_boost); + }); + + lboost_sum_rank = shamalgs::primitives::sum(dev_sched, l_boost, 0, Nobj); + }); + + Tvec lboost = shamalgs::collective::allreduce_sum(lboost_sum_rank); + + logger::raw_ln("lboost =", lboost); + + s.angular_momentum -= lboost; } sink_id++; From 2be5d2a462c0988566ba6c69fbd369f25a5568f0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Thu, 2 Apr 2026 16:41:10 +0200 Subject: [PATCH 03/23] clean sping accretion --- .../sph/src/modules/SinkParticlesUpdate.cpp | 122 ++++++++++++++---- 1 file changed, 96 insertions(+), 26 deletions(-) diff --git a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp index 2e1bf44399..e125be5c8e 100644 --- a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp +++ b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp @@ -31,8 +31,7 @@ #include "shammath/sphkernels.hpp" #include "shammodels/sph/modules/SinkParticlesUpdate.hpp" #include "shamsys/legacy/log.hpp" -#include -#include +#include template class SPHKernel> void shammodels::sph::modules::SinkParticlesUpdate::accrete_particles() { @@ -50,12 +49,50 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par PatchDataLayerLayout &pdl = scheduler().pdl_old(); const u32 ixyz = pdl.get_field_idx("xyz"); const u32 ivxyz = pdl.get_field_idx("vxyz"); + const u32 iaxyz = pdl.get_field_idx("axyz"); auto dev_sched = shamsys::instance::get_compute_scheduler_ptr(); sham::DeviceQueue &q = shambase::get_check_ref(dev_sched).get_queue(); std::vector &sink_parts = storage.sinks.get(); + /* // Uncomment this if you want to compute total momentum to check some stuff + auto get_tot_mom = [&]() -> Tvec { + Tvec total_momentum = {}; + + sham::DeviceBuffer total_momentum_part(0, dev_sched); + + scheduler().for_each_patchdata_nonempty( + [&](const shamrock::patch::Patch p, shamrock::patch::PatchDataLayer &pdat) { + u32 len = pdat.get_obj_cnt(); + + total_momentum_part.resize(len); + + sham::DeviceBuffer &vxyz_buf = pdat.get_field_buf_ref(ivxyz); + + sham::kernel_call( + q, + sham::MultiRef{vxyz_buf}, + sham::MultiRef{total_momentum_part}, + len, + [gpart_mass]( + u32 i, const Tvec *__restrict vxyz, Tvec *__restrict total_momentum_part) { + total_momentum_part[i] = gpart_mass * vxyz[i]; + }); + + total_momentum += shamalgs::primitives::sum(dev_sched, total_momentum_part, 0, len); + }); + + Tvec tot_total_momentum = shamalgs::collective::allreduce_sum(total_momentum); + + for (auto &sink : sink_parts) { + tot_total_momentum += sink.mass * sink.velocity; + } + + return tot_total_momentum; + }; +*/ + u32 sink_id = 0; bool had_accretion = false; std::string log = "sink accretion :"; @@ -121,16 +158,15 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par Tscal s_acc_mass = 0; Tvec s_acc_mxyz = {0, 0, 0}; Tvec s_acc_pxyz = {0, 0, 0}; + Tvec s_acc_maxyz = {0, 0, 0}; + Tvec s_acc_lxyz = {0, 0, 0}; scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { u32 Nobj = pdat.get_obj_cnt(); sham::DeviceBuffer &buf_xyz = pdat.get_field_buf_ref(ixyz); sham::DeviceBuffer &buf_vxyz = pdat.get_field_buf_ref(ivxyz); - - sham::DeviceBuffer ¬_accreted - = accretion_flag_bufs.get(cur_p.id_patch).not_accreted; - sham::DeviceBuffer &accreted = accretion_flag_bufs.get(cur_p.id_patch).accreted; + sham::DeviceBuffer &buf_axyz = pdat.get_field_buf_ref(iaxyz); sham::DeviceBuffer &id_list_accrete = bufs_id_list_accrete.get(cur_p.id_patch); @@ -141,50 +177,61 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par Tscal acc_mass = gpart_mass * Naccrete; sham::DeviceBuffer pxyz_acc(Naccrete, dev_sched); + sham::DeviceBuffer maxyz_acc(Naccrete, dev_sched); sham::DeviceBuffer mxyz_acc(Naccrete, dev_sched); + sham::DeviceBuffer lxyz_acc(Naccrete, dev_sched); sham::kernel_call( q, - sham::MultiRef{buf_xyz, buf_vxyz, id_list_accrete}, - sham::MultiRef{pxyz_acc, mxyz_acc}, + sham::MultiRef{buf_xyz, buf_vxyz, buf_axyz, id_list_accrete}, + sham::MultiRef{pxyz_acc, mxyz_acc, maxyz_acc, lxyz_acc}, Naccrete, [gpart_mass, r_sink, v_sink]( u32 id_a, const Tvec *__restrict xyz, const Tvec *__restrict vxyz, + const Tvec *__restrict axyz, const u32 *__restrict id_acc, Tvec *accretion_p, - Tvec *accretion_mr) { + Tvec *accretion_mr, + Tvec *accretion_ma, + Tvec *accretion_l) { u32 i_a = id_acc[id_a]; Tvec r = xyz[i_a]; Tvec v = vxyz[i_a]; + Tvec a = axyz[i_a]; accretion_p[id_a] = gpart_mass * v; accretion_mr[id_a] = gpart_mass * r; + accretion_ma[id_a] = gpart_mass * a; + accretion_l[id_a] = gpart_mass * sycl::cross(r - r_sink, v - v_sink); }); - Tvec acc_pxyz = shamalgs::primitives::sum(dev_sched, pxyz_acc, 0, Naccrete); - Tvec acc_mxyz = shamalgs::primitives::sum(dev_sched, mxyz_acc, 0, Naccrete); + Tvec acc_pxyz = shamalgs::primitives::sum(dev_sched, pxyz_acc, 0, Naccrete); + Tvec acc_mxyz = shamalgs::primitives::sum(dev_sched, mxyz_acc, 0, Naccrete); + Tvec acc_maxyz = shamalgs::primitives::sum(dev_sched, maxyz_acc, 0, Naccrete); + Tvec acc_lxyz = shamalgs::primitives::sum(dev_sched, lxyz_acc, 0, Naccrete); s_acc_mass += acc_mass; s_acc_pxyz += acc_pxyz; s_acc_mxyz += acc_mxyz; + s_acc_maxyz += acc_maxyz; + s_acc_lxyz += acc_lxyz; } }); Tscal sum_acc_mass = shamalgs::collective::allreduce_sum(s_acc_mass); Tvec sum_acc_pxyz = shamalgs::collective::allreduce_sum(s_acc_pxyz); Tvec sum_acc_mxyz = shamalgs::collective::allreduce_sum(s_acc_mxyz); + Tvec sum_acc_maxyz = shamalgs::collective::allreduce_sum(s_acc_maxyz); + Tvec sum_acc_lxyz = shamalgs::collective::allreduce_sum(s_acc_lxyz); // compute the new sink values - Tscal new_mass = s.mass + sum_acc_mass; - Tvec new_pos = (sum_acc_mxyz + s.pos * s.mass) / (s.mass + sum_acc_mass); - Tvec new_vel = (sum_acc_pxyz + s.velocity * s.mass) / (s.mass + sum_acc_mass); - - // use them to compute the angular_momentum updated - Tvec s_acc_lxyz = {0, 0, 0}; - Tvec sum_acc_lxyz = shamalgs::collective::allreduce_sum(s_acc_lxyz); - - Tvec new_ang_mom = s.angular_momentum + sum_acc_lxyz; + Tscal new_mass = s.mass + sum_acc_mass; + Tvec new_pos = (sum_acc_mxyz + s.pos * s.mass) / (s.mass + sum_acc_mass); + Tvec new_vel = (sum_acc_pxyz + s.velocity * s.mass) / (s.mass + sum_acc_mass); + Tvec new_acc = (sum_acc_maxyz + s.sph_acceleration * s.mass) / (s.mass + sum_acc_mass); + Tvec new_ang_mom = s.angular_momentum + sum_acc_lxyz + - new_mass * sycl::cross(new_pos - s.pos, new_vel - s.velocity); // write back the updated sink state auto new_state = s; @@ -192,19 +239,23 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par new_state.pos = new_pos; new_state.velocity = new_vel; new_state.angular_momentum = new_ang_mom; - - s = new_state; + new_state.sph_acceleration = new_acc; if (sum_acc_mass > 0) { had_accretion = true; log += shambase::format( - "\n id {} mass {} vel {} l {}", + "\n id {} deltas : mass={} r={} v={} l={}", sink_id, - sum_acc_mass, - sum_acc_pxyz / s.mass, - sum_acc_lxyz); + new_state.mass - s.mass, + new_state.pos - s.pos, + new_state.velocity - s.velocity, + new_state.angular_momentum - s.angular_momentum); } + // Tvec tot_mom_before = get_tot_mom(); + + s = new_state; + // evict accreted particles from patches scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { u32 Nobj = pdat.get_obj_cnt(); @@ -225,6 +276,25 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par } }); + /* + Tvec tot_mom_after = get_tot_mom(); + Tvec delta = tot_mom_after - tot_mom_before; + + logger::raw_ln(tot_mom_before, tot_mom_after, tot_mom_after - tot_mom_before); + logger::raw_ln(sycl::length(delta), sycl::length(tot_mom_after)); + + if (sycl::length(delta) > 1e-12) { + throw ""; + } + + if (sycl::length(tot_mom_after) > 1e-12) { + throw ""; + } + */ + } + + // now time for torque free torque restitution + for (Sink &s : sink_parts) { if (s.is_torque_free && sham::dot(s.angular_momentum, s.angular_momentum) > 0) { // boost particles in a radius of 4 racc to empty the sink angular_momentum From e861974df5743ebab2f9d5ebd0d448a235d4e3b4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Thu, 2 Apr 2026 17:58:52 +0200 Subject: [PATCH 04/23] niiiiice --- .../sph/modules/AnalysisAngularMomentum.hpp | 96 +++++++++++++++++++ .../sph/modules/SinkParticlesUpdate.hpp | 2 +- src/shammodels/sph/src/Solver.cpp | 2 +- .../sph/src/modules/SinkParticlesUpdate.cpp | 59 +++++++++++- src/shammodels/sph/src/pySPHModel.cpp | 32 +++++++ 5 files changed, 186 insertions(+), 5 deletions(-) create mode 100644 src/shammodels/sph/include/shammodels/sph/modules/AnalysisAngularMomentum.hpp diff --git a/src/shammodels/sph/include/shammodels/sph/modules/AnalysisAngularMomentum.hpp b/src/shammodels/sph/include/shammodels/sph/modules/AnalysisAngularMomentum.hpp new file mode 100644 index 0000000000..2f2d56d03e --- /dev/null +++ b/src/shammodels/sph/include/shammodels/sph/modules/AnalysisAngularMomentum.hpp @@ -0,0 +1,96 @@ +// -------------------------------------------------------// +// +// 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 AnalysisAngularMomentum.hpp + * @author Timothée David--Cléris (tim.shamrock@proton.me) + * @brief AnalysisTotalMomentum class with one method AnalysisTotalMomentum.get_total_momentum() + * + */ + +#include "shambase/memory.hpp" +#include "shamalgs/primitives/reduction.hpp" +#include "shambackends/DeviceQueue.hpp" +#include "shambackends/DeviceScheduler.hpp" +#include "shammodels/sph/Model.hpp" +#include "shamrock/scheduler/PatchScheduler.hpp" +#include "shamrock/scheduler/ShamrockCtx.hpp" +#include + +namespace shammodels::sph::modules { + + template class SPHKernel> + class AnalysisAngularMomentum { + public: + using Tscal = shambase::VecComponent; + static constexpr u32 dim = shambase::VectorProperties::dimension; + + using Solver = Solver; + + Model &model; + Solver &solver; + ShamrockCtx &ctx; + + AnalysisAngularMomentum(Model &model) + : model(model), ctx(model.ctx), solver(model.solver) {}; + + auto get_angular_momentum() -> Tvec { + PatchScheduler &sched = shambase::get_check_ref(ctx.sched); + auto dev_sched_ptr = shamsys::instance::get_compute_scheduler_ptr(); + sham::DeviceQueue &q = shambase::get_check_ref(dev_sched_ptr).get_queue(); + + const u32 ivxyz = sched.pdl_old().template get_field_idx("vxyz"); + const u32 ixyz = sched.pdl_old().template get_field_idx("xyz"); + const Tscal pmass = solver.solver_config.gpart_mass; + + Tvec angular_momentum = {}; + + sham::DeviceBuffer total_momentum_part(0, dev_sched_ptr); + + sched.for_each_patchdata_nonempty( + [&](const shamrock::patch::Patch p, shamrock::patch::PatchDataLayer &pdat) { + u32 len = pdat.get_obj_cnt(); + + total_momentum_part.resize(len); + + sham::DeviceBuffer &xyz_buf = pdat.get_field_buf_ref(ixyz); + sham::DeviceBuffer &vxyz_buf = pdat.get_field_buf_ref(ivxyz); + + sham::kernel_call( + q, + sham::MultiRef{xyz_buf, vxyz_buf}, + sham::MultiRef{total_momentum_part}, + len, + [pmass]( + u32 i, + const Tvec *__restrict xyz, + const Tvec *__restrict vxyz, + Tvec *__restrict total_momentum_part) { + total_momentum_part[i] = pmass * sycl::cross(xyz[i], vxyz[i]); + }); + + angular_momentum + += shamalgs::primitives::sum(dev_sched_ptr, total_momentum_part, 0, len); + }); + + Tvec tot_angular_momentum = shamalgs::collective::allreduce_sum(angular_momentum); + + if (!solver.storage.sinks.is_empty()) { + for (auto &sink : solver.storage.sinks.get()) { + tot_angular_momentum + += sink.mass * sycl::cross(sink.pos, sink.velocity) + sink.angular_momentum; + } + } + + return tot_angular_momentum; + } + }; +} // namespace shammodels::sph::modules diff --git a/src/shammodels/sph/include/shammodels/sph/modules/SinkParticlesUpdate.hpp b/src/shammodels/sph/include/shammodels/sph/modules/SinkParticlesUpdate.hpp index b046983fdd..b7c1c240b8 100644 --- a/src/shammodels/sph/include/shammodels/sph/modules/SinkParticlesUpdate.hpp +++ b/src/shammodels/sph/include/shammodels/sph/modules/SinkParticlesUpdate.hpp @@ -44,7 +44,7 @@ namespace shammodels::sph::modules { SinkParticlesUpdate(ShamrockCtx &context, Config &solver_config, Storage &storage) : context(context), solver_config(solver_config), storage(storage) {} - void accrete_particles(); + void accrete_particles(Tscal dt); void predictor_step(Tscal dt); void compute_sph_forces(); void compute_ext_forces(); diff --git a/src/shammodels/sph/src/Solver.cpp b/src/shammodels/sph/src/Solver.cpp index 0ef3a5ad05..96ab7e0c57 100644 --- a/src/shammodels/sph/src/Solver.cpp +++ b/src/shammodels/sph/src/Solver.cpp @@ -1638,7 +1638,7 @@ shammodels::sph::TimestepLog shammodels::sph::Solver::evolve_once() modules::SinkParticlesUpdate sink_update(context, solver_config, storage); modules::ExternalForces ext_forces(context, solver_config, storage); - sink_update.accrete_particles(); + sink_update.accrete_particles(dt); ext_forces.point_mass_accrete_particles(); sink_update.predictor_step(dt); diff --git a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp index e125be5c8e..3706a3ca87 100644 --- a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp +++ b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp @@ -34,7 +34,7 @@ #include template class SPHKernel> -void shammodels::sph::modules::SinkParticlesUpdate::accrete_particles() { +void shammodels::sph::modules::SinkParticlesUpdate::accrete_particles(Tscal dt) { StackEntry stack_loc{}; Tscal gpart_mass = solver_config.gpart_mass; @@ -93,6 +93,44 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par }; */ + auto get_ang_mom = [&]() -> Tvec { + Tvec angular_momentum = {}; + + sham::DeviceBuffer total_momentum_part(0, dev_sched); + + scheduler().for_each_patchdata_nonempty([&](const shamrock::patch::Patch p, + shamrock::patch::PatchDataLayer &pdat) { + u32 len = pdat.get_obj_cnt(); + + total_momentum_part.resize(len); + + sham::DeviceBuffer &xyz_buf = pdat.get_field_buf_ref(ixyz); + sham::DeviceBuffer &vxyz_buf = pdat.get_field_buf_ref(ivxyz); + sham::kernel_call( + q, + sham::MultiRef{xyz_buf, vxyz_buf}, + sham::MultiRef{total_momentum_part}, + len, + [gpart_mass]( + u32 i, + const Tvec *__restrict xyz, + const Tvec *__restrict vxyz, + Tvec *__restrict total_momentum_part) { + total_momentum_part[i] = gpart_mass * sycl::cross(xyz[i], vxyz[i]); + }); + + angular_momentum += shamalgs::primitives::sum(dev_sched, total_momentum_part, 0, len); + }); + + Tvec tot_angular_momentum = shamalgs::collective::allreduce_sum(angular_momentum); + + for (auto &sink : sink_parts) { + tot_angular_momentum + += sink.mass * sycl::cross(sink.pos, sink.velocity) + sink.angular_momentum; + } + return tot_angular_momentum; + }; + u32 sink_id = 0; bool had_accretion = false; std::string log = "sink accretion :"; @@ -186,7 +224,7 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par sham::MultiRef{buf_xyz, buf_vxyz, buf_axyz, id_list_accrete}, sham::MultiRef{pxyz_acc, mxyz_acc, maxyz_acc, lxyz_acc}, Naccrete, - [gpart_mass, r_sink, v_sink]( + [gpart_mass, r_sink, v_sink, dt]( u32 id_a, const Tvec *__restrict xyz, const Tvec *__restrict vxyz, @@ -203,7 +241,8 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par accretion_p[id_a] = gpart_mass * v; accretion_mr[id_a] = gpart_mass * r; accretion_ma[id_a] = gpart_mass * a; - accretion_l[id_a] = gpart_mass * sycl::cross(r - r_sink, v - v_sink); + accretion_l[id_a] + = gpart_mass * sycl::cross(r - r_sink, v + ((dt / 2) * a) - v_sink); }); Tvec acc_pxyz = shamalgs::primitives::sum(dev_sched, pxyz_acc, 0, Naccrete); @@ -253,6 +292,7 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par } // Tvec tot_mom_before = get_tot_mom(); + Tvec ang_mom_before = get_ang_mom(); s = new_state; @@ -291,6 +331,19 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par throw ""; } */ + + /* + Tvec ang_mom_after = get_ang_mom(); + Tvec delta = ang_mom_after - ang_mom_before; + + logger::raw_ln(ang_mom_before, ang_mom_after, ang_mom_after - ang_mom_before); + logger::raw_ln(sycl::length(delta), sycl::length(ang_mom_after)); + + if (sycl::length(delta) > 1e-12) { + // throw ""; + } + + // */ } // now time for torque free torque restitution diff --git a/src/shammodels/sph/src/pySPHModel.cpp b/src/shammodels/sph/src/pySPHModel.cpp index 74da727ebb..2c30117f76 100644 --- a/src/shammodels/sph/src/pySPHModel.cpp +++ b/src/shammodels/sph/src/pySPHModel.cpp @@ -25,6 +25,7 @@ #include "shammodels/common/shamrock_json_to_py_json.hpp" #include "shammodels/sph/Model.hpp" #include "shammodels/sph/io/PhantomDump.hpp" +#include "shammodels/sph/modules/AnalysisAngularMomentum.hpp" #include "shammodels/sph/modules/AnalysisBarycenter.hpp" #include "shammodels/sph/modules/AnalysisDisc.hpp" #include "shammodels/sph/modules/AnalysisEnergyKinetic.hpp" @@ -1315,6 +1316,21 @@ void add_analysisTotalMomentum_instance(py::module &m, const std::string &name_m }); } +template class SPHKernel> +void add_analysisAngularMomentum_instance(py::module &m, const std::string &name_model) { + using namespace shammodels::sph; + + using Tscal = shambase::VecComponent; + using T = Model; + + py::class_>(m, name_model.c_str()) + .def(py::init([](T &model) { + return std::make_unique>(model); + })) + .def("get_angular_momentum", [](modules::AnalysisAngularMomentum &self) { + return self.get_angular_momentum(); + }); +} using namespace shammodels::sph; template class SPHKernel> @@ -1491,6 +1507,20 @@ Register_pymod(pysphmodel) { add_analysisTotalMomentum_instance(msph, "AnalysisTotalMomentum_f64_3_C4"); add_analysisTotalMomentum_instance(msph, "AnalysisTotalMomentum_f64_3_C6"); + add_analysisAngularMomentum_instance( + msph, "AnalysisAngularMomentum_f64_3_M4"); + add_analysisAngularMomentum_instance( + msph, "AnalysisAngularMomentum_f64_3_M6"); + add_analysisAngularMomentum_instance( + msph, "AnalysisAngularMomentum_f64_3_M8"); + + add_analysisAngularMomentum_instance( + msph, "AnalysisAngularMomentum_f64_3_C2"); + add_analysisAngularMomentum_instance( + msph, "AnalysisAngularMomentum_f64_3_C4"); + add_analysisAngularMomentum_instance( + msph, "AnalysisAngularMomentum_f64_3_C6"); + register_analysis_impl_for_each_kernel(msph, "analysisBarycenter"); register_analysis_impl_for_each_kernel( msph, "analysisEnergyKinetic"); @@ -1498,4 +1528,6 @@ Register_pymod(pysphmodel) { msph, "analysisEnergyPotential"); register_analysis_impl_for_each_kernel( msph, "analysisTotalMomentum"); + register_analysis_impl_for_each_kernel( + msph, "analysisAngularMomentum"); } From 6e6581636bb8f5ccb685ea542a273b00f770b19a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Thu, 2 Apr 2026 18:04:06 +0200 Subject: [PATCH 05/23] yeaaaaah --- .../sph/src/modules/SinkParticlesUpdate.cpp | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp index 3706a3ca87..07f0cd04bf 100644 --- a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp +++ b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp @@ -230,10 +230,10 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par const Tvec *__restrict vxyz, const Tvec *__restrict axyz, const u32 *__restrict id_acc, - Tvec *accretion_p, - Tvec *accretion_mr, - Tvec *accretion_ma, - Tvec *accretion_l) { + Tvec *__restrict accretion_p, + Tvec *__restrict accretion_mr, + Tvec *__restrict accretion_ma, + Tvec *__restrict accretion_l) { u32 i_a = id_acc[id_a]; Tvec r = xyz[i_a]; Tvec v = vxyz[i_a]; @@ -241,8 +241,11 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par accretion_p[id_a] = gpart_mass * v; accretion_mr[id_a] = gpart_mass * r; accretion_ma[id_a] = gpart_mass * a; - accretion_l[id_a] - = gpart_mass * sycl::cross(r - r_sink, v + ((dt / 2) * a) - v_sink); + + // dirty trick to account for the residual acceleration in the spin. This + // allows us to maitain a much better angular momentum conservation. + v += a * dt / 2; + accretion_l[id_a] = gpart_mass * sycl::cross(r - r_sink, v - v_sink); }); Tvec acc_pxyz = shamalgs::primitives::sum(dev_sched, pxyz_acc, 0, Naccrete); From 888c68a9c9e9c28ff736f4a11c002951a9d29315 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Thu, 2 Apr 2026 18:04:18 +0200 Subject: [PATCH 06/23] tmp do not commit this one --- examples/sph/run_circular_disc_sink.py | 67 ++++++++++++++++++++------ 1 file changed, 52 insertions(+), 15 deletions(-) diff --git a/examples/sph/run_circular_disc_sink.py b/examples/sph/run_circular_disc_sink.py index e0ba5c9318..9d411f4040 100644 --- a/examples/sph/run_circular_disc_sink.py +++ b/examples/sph/run_circular_disc_sink.py @@ -83,7 +83,7 @@ plot_freq_stop = 1 dt_stop = 0.01 -nstop = 30 +nstop = 120 # The list of times at which the simulation will pause for analysis / dumping t_stop = [i * dt_stop for i in range(nstop + 1)] @@ -462,18 +462,21 @@ def save_analysis_data(filename, key, value, ianalysis): def analysis(ianalysis): - column_density_plot.analysis_save(ianalysis) - column_density_plot_hollywood.analysis_save(ianalysis) - vertical_density_plot.analysis_save(ianalysis) - v_z_slice_plot.analysis_save(ianalysis) - relative_azy_velocity_slice_plot.analysis_save(ianalysis) - vertical_shear_gradient_slice_plot.analysis_save(ianalysis) - dt_part_slice_plot.analysis_save(ianalysis) - column_particle_count_plot.analysis_save(ianalysis) + # column_density_plot.analysis_save(ianalysis) + # column_density_plot_hollywood.analysis_save(ianalysis) + # vertical_density_plot.analysis_save(ianalysis) + # v_z_slice_plot.analysis_save(ianalysis) + # relative_azy_velocity_slice_plot.analysis_save(ianalysis) + # vertical_shear_gradient_slice_plot.analysis_save(ianalysis) + # dt_part_slice_plot.analysis_save(ianalysis) + # column_particle_count_plot.analysis_save(ianalysis) barycenter, disc_mass = shamrock.model_sph.analysisBarycenter(model=model).get_barycenter() total_momentum = shamrock.model_sph.analysisTotalMomentum(model=model).get_total_momentum() + angular_momentum = shamrock.model_sph.analysisAngularMomentum( + model=model + ).get_angular_momentum() potential_energy = shamrock.model_sph.analysisEnergyPotential( model=model @@ -484,6 +487,7 @@ def analysis(ianalysis): save_analysis_data("barycenter.json", "barycenter", barycenter, ianalysis) save_analysis_data("disc_mass.json", "disc_mass", disc_mass, ianalysis) save_analysis_data("total_momentum.json", "total_momentum", total_momentum, ianalysis) + save_analysis_data("angular_momentum.json", "angular_momentum", angular_momentum, ianalysis) save_analysis_data("potential_energy.json", "potential_energy", potential_energy, ianalysis) save_analysis_data("kinetic_energy.json", "kinetic_energy", kinetic_energy, ianalysis) @@ -534,6 +538,7 @@ def analysis(ianalysis): import matplotlib import matplotlib.pyplot as plt +""" face_on_render_kwargs = { "x_unit": "au", "y_unit": "au", @@ -634,6 +639,7 @@ def analysis(ianalysis): contour_list=[1, 10, 100, 1000], **sink_params, ) +""" # %% # Make gif for the doc (plot_to_gif.py) # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -642,7 +648,7 @@ def analysis(ianalysis): # sphinx_gallery_multi_image = "single" -render_gif = True +render_gif = False # %% # Do it for rho integ @@ -739,7 +745,6 @@ def load_data_from_json(filename, key): plt.ylabel("barycenter") plt.legend(["x", "y", "z"]) plt.savefig(analysis_folder + "barycenter.png") -plt.show() # %% # load the json file for disc_mass @@ -751,7 +756,6 @@ def load_data_from_json(filename, key): plt.xlabel("t") plt.ylabel("disc_mass") plt.savefig(analysis_folder + "disc_mass.png") -plt.show() # %% # load the json file for total_momentum @@ -769,7 +773,24 @@ def load_data_from_json(filename, key): plt.ylabel("total_momentum") plt.legend(["x", "y", "z"]) plt.savefig(analysis_folder + "total_momentum.png") -plt.show() + +# %% +# load the json file for total_momentum +t, angular_momentum = load_data_from_json("angular_momentum.json", "angular_momentum") +angular_momentum_x = [d[0] - angular_momentum[0][0] for d in angular_momentum] +angular_momentum_y = [d[1] - angular_momentum[0][1] for d in angular_momentum] +angular_momentum_z = [d[2] - angular_momentum[0][2] for d in angular_momentum] + + +plt.figure(figsize=(8, 5), dpi=200) + +plt.plot(t, angular_momentum_x) +plt.plot(t, angular_momentum_y) +plt.plot(t, angular_momentum_z) +plt.xlabel("t") +plt.ylabel("total angular_momentum") +plt.legend(["x", "y", "z"]) +plt.savefig(analysis_folder + "angular_momentum.png") # %% # load the json file for energies @@ -786,7 +807,6 @@ def load_data_from_json(filename, key): plt.ylabel("energy") plt.legend(["potential_energy", "kinetic_energy", "total_energy"]) plt.savefig(analysis_folder + "energies.png") -plt.show() # %% # load the json file for sinks @@ -804,7 +824,6 @@ def load_data_from_json(filename, key): plt.ylabel("sink position") plt.legend() plt.savefig(analysis_folder + "sinks.png") -plt.show() # %% # Sink to barycenter distance @@ -828,6 +847,24 @@ def load_data_from_json(filename, key): plt.ylabel("sink pos - barycenter pos") plt.legend() plt.savefig(analysis_folder + "sink_to_barycenter_distance.png") + +# %% +# Sink angular momentum +t, sinks = load_data_from_json("sinks.json", "sinks") + +sinks_lx = np.array([d[0]["angular_momentum"][0] for d in sinks]) +sinks_ly = np.array([d[0]["angular_momentum"][1] for d in sinks]) +sinks_lz = np.array([d[0]["angular_momentum"][2] for d in sinks]) + + +plt.figure(figsize=(8, 5), dpi=200) +plt.plot(t, sinks_lx, label="sink 0 (l_x)") +plt.plot(t, sinks_ly, label="sink 0 (l_y)") +plt.plot(t, sinks_lz, label="sink 0 (l_z)") +plt.xlabel("t") +plt.ylabel("sink spin") +plt.legend() +plt.savefig(analysis_folder + "sink_angmom.png") plt.show() From 394be15151a1ca2de359cb83a5814ff5ed71dfc4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Fri, 3 Apr 2026 10:14:27 +0200 Subject: [PATCH 07/23] sync with main --- examples/sph/run_circular_disc_sink.py | 66 ++++++-------------------- 1 file changed, 14 insertions(+), 52 deletions(-) diff --git a/examples/sph/run_circular_disc_sink.py b/examples/sph/run_circular_disc_sink.py index 2000aae1fe..0ccfef2c99 100644 --- a/examples/sph/run_circular_disc_sink.py +++ b/examples/sph/run_circular_disc_sink.py @@ -82,13 +82,8 @@ dump_freq_stop = 2 plot_freq_stop = 1 -<<<<<<< torque-free-sink -dt_stop = 0.01 -nstop = 120 -======= dt_stop = 0.02 nstop = 30 ->>>>>>> main # The list of times at which the simulation will pause for analysis / dumping t_stop = [i * dt_stop for i in range(nstop + 1)] @@ -467,14 +462,14 @@ def save_analysis_data(filename, key, value, ianalysis): def analysis(ianalysis): - # column_density_plot.analysis_save(ianalysis) - # column_density_plot_hollywood.analysis_save(ianalysis) - # vertical_density_plot.analysis_save(ianalysis) - # v_z_slice_plot.analysis_save(ianalysis) - # relative_azy_velocity_slice_plot.analysis_save(ianalysis) - # vertical_shear_gradient_slice_plot.analysis_save(ianalysis) - # dt_part_slice_plot.analysis_save(ianalysis) - # column_particle_count_plot.analysis_save(ianalysis) + column_density_plot.analysis_save(ianalysis) + column_density_plot_hollywood.analysis_save(ianalysis) + vertical_density_plot.analysis_save(ianalysis) + v_z_slice_plot.analysis_save(ianalysis) + relative_azy_velocity_slice_plot.analysis_save(ianalysis) + vertical_shear_gradient_slice_plot.analysis_save(ianalysis) + dt_part_slice_plot.analysis_save(ianalysis) + column_particle_count_plot.analysis_save(ianalysis) barycenter, disc_mass = shamrock.model_sph.analysisBarycenter(model=model).get_barycenter() @@ -543,7 +538,6 @@ def analysis(ianalysis): import matplotlib import matplotlib.pyplot as plt -""" face_on_render_kwargs = { "x_unit": "au", "y_unit": "au", @@ -644,7 +638,6 @@ def analysis(ianalysis): contour_list=[1, 10, 100, 1000], **sink_params, ) -""" # %% # Make gif for the doc (plot_to_gif.py) # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -653,7 +646,7 @@ def analysis(ianalysis): # sphinx_gallery_multi_image = "single" -render_gif = False +render_gif = True # %% # Do it for rho integ @@ -750,6 +743,7 @@ def load_data_from_json(filename, key): plt.ylabel("barycenter") plt.legend(["x", "y", "z"]) plt.savefig(analysis_folder + "barycenter.png") +plt.show() # %% # load the json file for disc_mass @@ -761,6 +755,7 @@ def load_data_from_json(filename, key): plt.xlabel("t") plt.ylabel("disc_mass") plt.savefig(analysis_folder + "disc_mass.png") +plt.show() # %% # load the json file for total_momentum @@ -778,24 +773,7 @@ def load_data_from_json(filename, key): plt.ylabel("total_momentum") plt.legend(["x", "y", "z"]) plt.savefig(analysis_folder + "total_momentum.png") - -# %% -# load the json file for total_momentum -t, angular_momentum = load_data_from_json("angular_momentum.json", "angular_momentum") -angular_momentum_x = [d[0] - angular_momentum[0][0] for d in angular_momentum] -angular_momentum_y = [d[1] - angular_momentum[0][1] for d in angular_momentum] -angular_momentum_z = [d[2] - angular_momentum[0][2] for d in angular_momentum] - - -plt.figure(figsize=(8, 5), dpi=200) - -plt.plot(t, angular_momentum_x) -plt.plot(t, angular_momentum_y) -plt.plot(t, angular_momentum_z) -plt.xlabel("t") -plt.ylabel("total angular_momentum") -plt.legend(["x", "y", "z"]) -plt.savefig(analysis_folder + "angular_momentum.png") +plt.show() # %% @@ -831,6 +809,7 @@ def load_data_from_json(filename, key): plt.ylabel("energy") plt.legend(["potential_energy", "kinetic_energy", "total_energy"]) plt.savefig(analysis_folder + "energies.png") +plt.show() # %% # load the json file for sinks @@ -848,6 +827,7 @@ def load_data_from_json(filename, key): plt.ylabel("sink position") plt.legend() plt.savefig(analysis_folder + "sinks.png") +plt.show() # %% # Sink angular momentum @@ -890,24 +870,6 @@ def load_data_from_json(filename, key): plt.ylabel("sink pos - barycenter pos") plt.legend() plt.savefig(analysis_folder + "sink_to_barycenter_distance.png") - -# %% -# Sink angular momentum -t, sinks = load_data_from_json("sinks.json", "sinks") - -sinks_lx = np.array([d[0]["angular_momentum"][0] for d in sinks]) -sinks_ly = np.array([d[0]["angular_momentum"][1] for d in sinks]) -sinks_lz = np.array([d[0]["angular_momentum"][2] for d in sinks]) - - -plt.figure(figsize=(8, 5), dpi=200) -plt.plot(t, sinks_lx, label="sink 0 (l_x)") -plt.plot(t, sinks_ly, label="sink 0 (l_y)") -plt.plot(t, sinks_lz, label="sink 0 (l_z)") -plt.xlabel("t") -plt.ylabel("sink spin") -plt.legend() -plt.savefig(analysis_folder + "sink_angmom.png") plt.show() From a5b075f2d15076b114b61c59e81354a4f11efd0a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Sat, 4 Apr 2026 18:48:38 +0200 Subject: [PATCH 08/23] progress --- .../utils/analysis/AngularMomentumPlots.py | 165 ++++++++++++++++++ src/pylib/shamrock/utils/analysis/__init__.py | 5 + .../sph/include/shammodels/sph/Model.hpp | 8 +- .../include/shammodels/sph/SinkPartStruct.hpp | 2 +- src/shammodels/sph/src/Model.cpp | 3 +- .../sph/src/modules/SinkParticlesUpdate.cpp | 60 ++++--- src/shammodels/sph/src/pySPHModel.cpp | 10 +- 7 files changed, 225 insertions(+), 28 deletions(-) create mode 100644 src/pylib/shamrock/utils/analysis/AngularMomentumPlots.py diff --git a/src/pylib/shamrock/utils/analysis/AngularMomentumPlots.py b/src/pylib/shamrock/utils/analysis/AngularMomentumPlots.py new file mode 100644 index 0000000000..955c4338e0 --- /dev/null +++ b/src/pylib/shamrock/utils/analysis/AngularMomentumPlots.py @@ -0,0 +1,165 @@ +import numpy as np + +import shamrock.sys + +from .StandardPlotHelper import StandardPlotHelper + +try: + from numba import njit + + _HAS_NUMBA = True +except ImportError: + _HAS_NUMBA = False + + + + +def SliceAngularMomentum( + model, + ext_r, + nx, + ny, + ex, + ey, + center, + analysis_folder, + analysis_prefix, + do_normalization=True, + min_normalization=1e-9, + Lprojection = [0.,0.,1.] +): + def compute_angular_mom(helper): + if _HAS_NUMBA: + if shamrock.sys.world_rank() == 0: + print("Using numba for angular momentum in SliceAngularMomentum") + + pmass = model.get_particle_mass() + hfact = model.get_hfact() + + def internal( + size: int, + hpart: np.array, x: np.array, y: np.array, z: np.array, vx: np.array, vy: np.array, vz: np.array + ) -> np.array: + + rho = pmass * (hfact / hpart) ** 3 + + r = np.stack([x, y, z], axis=-1) + v = np.stack([vx, vy, vz], axis=-1) + L = rho[:, None] * np.cross(r, v) + + # project L onto the Lprojection vector + L = np.sum(L * Lprojection, axis=-1) + + return L + + if _HAS_NUMBA: + internal = njit(internal) + + def custom_getter(size: int, dic_out: dict) -> np.array: + return internal( + size, + dic_out["hpart"], + dic_out["xyz"][:, 0], + dic_out["xyz"][:, 1], + dic_out["xyz"][:, 2], + dic_out["vxyz"][:, 0], + dic_out["vxyz"][:, 1], + dic_out["vxyz"][:, 2], + ) + + arr_L = helper.slice_render( + "custom", + "f64", + do_normalization, + min_normalization, + custom_getter=custom_getter, + ) + + return arr_L + + return StandardPlotHelper( + model, + ext_r, + nx, + ny, + ex, + ey, + center, + analysis_folder, + analysis_prefix, + compute_function=compute_angular_mom, + ) + +def ColumnAverageAngularMomentum( + model, + ext_r, + nx, + ny, + ex, + ey, + center, + analysis_folder, + analysis_prefix, + min_normalization=1e-9, + Lprojection = [0.,0.,1.] +): + def compute_angular_mom(helper): + if _HAS_NUMBA: + if shamrock.sys.world_rank() == 0: + print("Using numba for angular momentum in SliceAngularMomentum") + + pmass = model.get_particle_mass() + hfact = model.get_hfact() + + def internal( + size: int, + hpart: np.array, x: np.array, y: np.array, z: np.array, vx: np.array, vy: np.array, vz: np.array + ) -> np.array: + + rho = pmass * (hfact / hpart) ** 3 + + r = np.stack([x, y, z], axis=-1) + v = np.stack([vx, vy, vz], axis=-1) + L = rho[:, None] * np.cross(r, v) + + # project L onto the Lprojection vector + L = np.sum(L * Lprojection, axis=-1) + + return L + + if _HAS_NUMBA: + internal = njit(internal) + + def custom_getter(size: int, dic_out: dict) -> np.array: + return internal( + size, + dic_out["hpart"], + dic_out["xyz"][:, 0], + dic_out["xyz"][:, 1], + dic_out["xyz"][:, 2], + dic_out["vxyz"][:, 0], + dic_out["vxyz"][:, 1], + dic_out["vxyz"][:, 2], + ) + + arr_L = helper.column_average_render( + "custom", + "f64", + min_normalization, + custom_getter=custom_getter, + ) + + return arr_L + + return StandardPlotHelper( + model, + ext_r, + nx, + ny, + ex, + ey, + center, + analysis_folder, + analysis_prefix, + compute_function=compute_angular_mom, + ) diff --git a/src/pylib/shamrock/utils/analysis/__init__.py b/src/pylib/shamrock/utils/analysis/__init__.py index 75cd2ed055..ed8679bb50 100644 --- a/src/pylib/shamrock/utils/analysis/__init__.py +++ b/src/pylib/shamrock/utils/analysis/__init__.py @@ -21,5 +21,10 @@ SliceBVerticalShearGradient, ) +from .AngularMomentumPlots import ( + SliceAngularMomentum, + ColumnAverageAngularMomentum, +) + # Performance analysis from .PerfHistory import PerfHistory diff --git a/src/shammodels/sph/include/shammodels/sph/Model.hpp b/src/shammodels/sph/include/shammodels/sph/Model.hpp index eb78e1feee..4fa9296dc5 100644 --- a/src/shammodels/sph/include/shammodels/sph/Model.hpp +++ b/src/shammodels/sph/include/shammodels/sph/Model.hpp @@ -173,15 +173,17 @@ namespace shammodels::sph { Tscal q, std::mt19937 eng); - inline void add_sink(Tscal mass, Tvec pos, Tvec velocity, Tscal accretion_radius) { + inline void add_sink( + Tscal mass, Tvec pos, Tvec velocity, Tscal accretion_radius, bool is_torque_free) { if (solver.storage.sinks.is_empty()) { solver.storage.sinks.set({}); } - shamlog_debug_ln("SPH", "add sink :", mass, pos, velocity, accretion_radius); + shamlog_debug_ln( + "SPH", "add sink :", mass, pos, velocity, accretion_radius, is_torque_free); solver.storage.sinks.get().push_back( - {pos, velocity, {}, {}, mass, {}, accretion_radius}); + {pos, velocity, {}, {}, mass, {}, accretion_radius, is_torque_free}); } template diff --git a/src/shammodels/sph/include/shammodels/sph/SinkPartStruct.hpp b/src/shammodels/sph/include/shammodels/sph/SinkPartStruct.hpp index 0c64835d79..23492deb0d 100644 --- a/src/shammodels/sph/include/shammodels/sph/SinkPartStruct.hpp +++ b/src/shammodels/sph/include/shammodels/sph/SinkPartStruct.hpp @@ -35,7 +35,7 @@ namespace shammodels::sph { Tvec angular_momentum; Tscal accretion_radius; - bool is_torque_free = false; + bool is_torque_free; }; template diff --git a/src/shammodels/sph/src/Model.cpp b/src/shammodels/sph/src/Model.cpp index 80b2473925..bae04a108d 100644 --- a/src/shammodels/sph/src/Model.cpp +++ b/src/shammodels/sph/src/Model.cpp @@ -1437,7 +1437,8 @@ void shammodels::sph::Model::init_from_phantom_dump( mass[i], {xsink[i], ysink[i], zsink[i]}, {vxsink[i], vysink[i], vzsink[i]}, - Racc[i]); + Racc[i], + false); } } } diff --git a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp index d4af1fa93c..8473d79fc4 100644 --- a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp +++ b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp @@ -63,12 +63,13 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par sham::DeviceBuffer total_momentum_part(0, dev_sched); scheduler().for_each_patchdata_nonempty( - [&](const shamrock::patch::Patch p, shamrock::patch::PatchDataLayer &pdat) { - u32 len = pdat.get_obj_cnt(); + [&](const shamrock::patch::Patch p, shamrock::patch::PatchDataLayer + &pdat) { u32 len = pdat.get_obj_cnt(); total_momentum_part.resize(len); - sham::DeviceBuffer &vxyz_buf = pdat.get_field_buf_ref(ivxyz); + sham::DeviceBuffer &vxyz_buf = + pdat.get_field_buf_ref(ivxyz); sham::kernel_call( q, @@ -76,14 +77,16 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par sham::MultiRef{total_momentum_part}, len, [gpart_mass]( - u32 i, const Tvec *__restrict vxyz, Tvec *__restrict total_momentum_part) { - total_momentum_part[i] = gpart_mass * vxyz[i]; + u32 i, const Tvec *__restrict vxyz, Tvec *__restrict + total_momentum_part) { total_momentum_part[i] = gpart_mass * vxyz[i]; }); - total_momentum += shamalgs::primitives::sum(dev_sched, total_momentum_part, 0, len); + total_momentum += shamalgs::primitives::sum(dev_sched, + total_momentum_part, 0, len); }); - Tvec tot_total_momentum = shamalgs::collective::allreduce_sum(total_momentum); + Tvec tot_total_momentum = + shamalgs::collective::allreduce_sum(total_momentum); for (auto &sink : sink_parts) { tot_total_momentum += sink.mass * sink.velocity; @@ -91,7 +94,7 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par return tot_total_momentum; }; -*/ + */ auto get_ang_mom = [&]() -> Tvec { Tvec angular_momentum = {}; @@ -243,8 +246,9 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par accretion_mr[id_a] = gpart_mass * r; accretion_ma[id_a] = gpart_mass * a; - // dirty trick to account for the residual acceleration in the spin. This - // allows us to maitain a much better angular momentum conservation. + // dirty trick to account for the residual acceleration in the + // spin. This allows us to maitain a much better angular momentum + // conservation. v += a * dt / 2; accretion_l[id_a] = gpart_mass * sycl::cross(r - r_sink, v - v_sink); }); @@ -321,7 +325,7 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par } }); } - + // now time for torque free torque restitution for (Sink &s : sink_parts) { if (s.is_torque_free && sham::dot(s.angular_momentum, s.angular_momentum) > 0) { @@ -330,15 +334,15 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par struct BoostWeight { Tscal radius; Tscal weight(Tscal r) const { - if (r < radius) { - return 1; + if (r < radius * Kernel::Rkern) { + return Kernel::W_3d(r, radius); } else { return 0; } } }; - BoostWeight weight_func{s.accretion_radius * 10}; + BoostWeight weight_func{s.accretion_radius * 2}; Tvec r_sink = s.pos; @@ -421,6 +425,7 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par logger::raw_ln("delta w =", delta_w.data[0], delta_w.data[1], delta_w.data[2]); Tvec lboost_sum_rank{}; + Tvec pboost_sum_rank{}; scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) { u32 Nobj = pdat.get_obj_cnt(); @@ -428,14 +433,19 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par sham::DeviceBuffer &buf_vxyz = pdat.get_field_buf_ref(ivxyz); sham::DeviceBuffer l_boost(Nobj, dev_sched); + sham::DeviceBuffer p_boost(Nobj, dev_sched); sham::kernel_call( q, sham::MultiRef{buf_xyz}, - sham::MultiRef{buf_vxyz, l_boost}, + sham::MultiRef{buf_vxyz, l_boost, p_boost}, Nobj, [r_sink, delta_w = vec_delta_w, weight_func, gpart_mass]( - u32 id_a, const Tvec *__restrict xyz, Tvec *vxyz, Tvec *lboost) { + u32 id_a, + const Tvec *__restrict xyz, + Tvec *__restrict vxyz, + Tvec *__restrict lboost, + Tvec *__restrict pboost) { Tvec r_a = xyz[id_a] - r_sink; Tscal r_a2 = sham::dot(r_a, r_a); @@ -445,17 +455,22 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par vxyz[id_a] += v_boost; lboost[id_a] += gpart_mass * sycl::cross(r_a, v_boost); + pboost[id_a] += gpart_mass * v_boost; }); - lboost_sum_rank = shamalgs::primitives::sum(dev_sched, l_boost, 0, Nobj); + lboost_sum_rank += shamalgs::primitives::sum(dev_sched, l_boost, 0, Nobj); + pboost_sum_rank += shamalgs::primitives::sum(dev_sched, p_boost, 0, Nobj); }); Tvec lboost = shamalgs::collective::allreduce_sum(lboost_sum_rank); + Tvec pboost = shamalgs::collective::allreduce_sum(pboost_sum_rank); - logger::raw_ln("lboost =", lboost); + logger::raw_ln("lboost =", lboost, "pboost =", pboost); s.angular_momentum -= lboost; -} + s.velocity -= pboost / s.mass; + } + } if (shamcomm::world_rank() == 0 && had_accretion) { logger::info_ln("sph::Sink", log); @@ -568,9 +583,10 @@ void shammodels::sph::modules::SinkParticlesUpdate::compute_sph Tvec force = G * delta / (d * d * d); - // This is a hack to avoid the sink kaboom effect - // when the particle is being advected close to the sink before - // being accreted + // This is a hack to avoid the sink + // kaboom effect when the particle is + // being advected close to the sink + // before being accreted if (d < sink_racc) { force = {0, 0, 0}; } diff --git a/src/shammodels/sph/src/pySPHModel.cpp b/src/shammodels/sph/src/pySPHModel.cpp index 2c30117f76..a9c5e1f9b0 100644 --- a/src/shammodels/sph/src/pySPHModel.cpp +++ b/src/shammodels/sph/src/pySPHModel.cpp @@ -824,7 +824,14 @@ void add_instance(py::module &m, std::string name_config, std::string name_model return self.solver.solver_config; }) .def("set_solver_config", &T::set_solver_config) - .def("add_sink", &T::add_sink) + .def( + "add_sink", + &T::add_sink, + py::arg("mass"), + py::arg("pos"), + py::arg("velocity"), + py::arg("accretion_radius"), + py::arg("is_torque_free") = false) .def( "get_sinks", [](T &self) { @@ -840,6 +847,7 @@ void add_instance(py::module &m, std::string name_config, std::string name_model sink_dic["mass"] = sink.mass; sink_dic["angular_momentum"] = sink.angular_momentum; sink_dic["accretion_radius"] = sink.accretion_radius; + sink_dic["is_torque_free"] = sink.is_torque_free; list_out.append(sink_dic); } } From 8f92e4c0049c09a3b2cf3074382a6f72dd33d433 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Sat, 4 Apr 2026 18:49:35 +0200 Subject: [PATCH 09/23] whoops --- .../utils/analysis/AngularMomentumPlots.py | 27 +++++++++++++------ 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/src/pylib/shamrock/utils/analysis/AngularMomentumPlots.py b/src/pylib/shamrock/utils/analysis/AngularMomentumPlots.py index 955c4338e0..220a47dce4 100644 --- a/src/pylib/shamrock/utils/analysis/AngularMomentumPlots.py +++ b/src/pylib/shamrock/utils/analysis/AngularMomentumPlots.py @@ -12,8 +12,6 @@ _HAS_NUMBA = False - - def SliceAngularMomentum( model, ext_r, @@ -26,7 +24,7 @@ def SliceAngularMomentum( analysis_prefix, do_normalization=True, min_normalization=1e-9, - Lprojection = [0.,0.,1.] + Lprojection=[0.0, 0.0, 1.0], ): def compute_angular_mom(helper): if _HAS_NUMBA: @@ -38,9 +36,15 @@ def compute_angular_mom(helper): def internal( size: int, - hpart: np.array, x: np.array, y: np.array, z: np.array, vx: np.array, vy: np.array, vz: np.array + hpart: np.array, + x: np.array, + y: np.array, + z: np.array, + vx: np.array, + vy: np.array, + vz: np.array, ) -> np.array: - + rho = pmass * (hfact / hpart) ** 3 r = np.stack([x, y, z], axis=-1) @@ -90,6 +94,7 @@ def custom_getter(size: int, dic_out: dict) -> np.array: compute_function=compute_angular_mom, ) + def ColumnAverageAngularMomentum( model, ext_r, @@ -101,7 +106,7 @@ def ColumnAverageAngularMomentum( analysis_folder, analysis_prefix, min_normalization=1e-9, - Lprojection = [0.,0.,1.] + Lprojection=[0.0, 0.0, 1.0], ): def compute_angular_mom(helper): if _HAS_NUMBA: @@ -113,9 +118,15 @@ def compute_angular_mom(helper): def internal( size: int, - hpart: np.array, x: np.array, y: np.array, z: np.array, vx: np.array, vy: np.array, vz: np.array + hpart: np.array, + x: np.array, + y: np.array, + z: np.array, + vx: np.array, + vy: np.array, + vz: np.array, ) -> np.array: - + rho = pmass * (hfact / hpart) ** 3 r = np.stack([x, y, z], axis=-1) From ccc0d07d08f424fd4d16b84b101abf2476be6cc0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Sat, 4 Apr 2026 18:51:43 +0200 Subject: [PATCH 10/23] early exit --- src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp index 8473d79fc4..27da533a12 100644 --- a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp +++ b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp @@ -368,6 +368,10 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par Tscal w = gpart_mass * weight_func.weight(sycl::sqrt(r_a2)); + if (w == 0) { + return; + } + mat I{}; { // I = w [Id r_a^2 - r_a r_a^T] From 6d4c0f3f8b5b769256b296b1ea93453d42e8fbf5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Tue, 7 Apr 2026 01:09:53 +0200 Subject: [PATCH 11/23] noicer --- examples/sph/run_compare_torque_free_sink.py | 579 ++++++++++++++++++ .../shamalgs/primitives/compute_histogram.hpp | 2 +- src/shammodels/common/src/pyCommonUtils.cpp | 4 + 3 files changed, 584 insertions(+), 1 deletion(-) create mode 100644 examples/sph/run_compare_torque_free_sink.py diff --git a/examples/sph/run_compare_torque_free_sink.py b/examples/sph/run_compare_torque_free_sink.py new file mode 100644 index 0000000000..a0c952708c --- /dev/null +++ b/examples/sph/run_compare_torque_free_sink.py @@ -0,0 +1,579 @@ +""" +Production run: Circular disc & central sink particle +===================================================== + +This example demonstrates how to run a smoothed particle hydrodynamics (SPH) +simulation of a circular disc orbiting around a central sink. + +The simulation models: + +- A central star with a given mass and accretion radius +- A gaseous disc with specified mass, inner/outer radii, and vertical structure +- Artificial viscosity for angular momentum transport +- Locally isothermal equation of state + +Also this simulation feature rolling dumps (see `purge_old_dumps` function) to save disk space. + +This example is the accumulation of 3 files in a single one to showcase the complete workflow. + +- The actual run script (runscript.py) +- Plot generation (make_plots.py) +- Animation from the plots (plot_to_gif.py) + +On a cluster or laptop, one can run the code as follows: + +.. code-block:: bash + + mpirun ./shamrock --sycl-cfg 0:0 --loglevel 1 --rscript runscript.py + + +then after the run is done (or while it is running), one can run the following to generate the plots: + +.. code-block:: bash + + python make_plots.py + + +""" + +# %% +# Runscript (runscript.py) +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# The runscript is the actual simulation with on the fly analysis & rolling dumps + +import glob +import json +import os # for makedirs + +import numpy as np + +try: + from numba import njit + + _HAS_NUMBA = True +except ImportError: + _HAS_NUMBA = False + + +from shamrock.utils.analysis import ( + AnalysisHelper, + ColumnDensityPlot, + ColumnParticleCount, + PerfHistory, + SliceDensityPlot, + SliceDiffVthetaProfile, + SliceDtPart, + SliceVzPlot, + VerticalShearGradient, +) + +import shamrock + +# If we use the shamrock executable to run this script instead of the python interpreter, +# we should not initialize the system as the shamrock executable needs to handle specific MPI logic +if not shamrock.sys.is_initialized(): + shamrock.change_loglevel(1) + shamrock.sys.init("0:0") + +# %% +# Setup units + +si = shamrock.UnitSystem() +sicte = shamrock.Constants(si) +codeu = shamrock.UnitSystem( + unit_time=sicte.year(), + unit_length=sicte.au(), + unit_mass=sicte.sol_mass(), +) +ucte = shamrock.Constants(codeu) +G = ucte.G() + +# %% +# List parameters + + +# Resolution +Npart = 200000 + +# Domain decomposition parameters +scheduler_split_val = int(1.0e7) # split patches with more than 1e7 particles +scheduler_merge_val = scheduler_split_val // 16 + +# Dump and plot frequency and duration of the simulation +dump_freq_stop = 2 + +dt_stop = 1 +nstop = 300 + +# The list of times at which the simulation will pause for analysis / dumping +t_stop = [i * dt_stop for i in range(nstop + 1)] + + +# Disc parameter +disc_mass = 0.01 # sol mass +rout = 20.0 # au +rin = 1.0 # au +H_r_0 = 0.05 +q = 0.5 +p = 3.0 / 2.0 +r0 = 1.0 + +# Viscosity parameter +alpha_AV = 1.0e-3 / 0.08 +alpha_u = 1.0 +beta_AV = 2.0 + +# Integrator parameters +C_cour = 0.3 +C_force = 0.25 + + +# Sink parameters +center_mass = 1.0 + + +# Disc profiles +def sigma_profile(r): + sigma_0 = 1.0 # We do not care as it will be renormalized + return sigma_0 * (r / r0) ** (-p) + + +def kep_profile(r): + return (G * center_mass / r) ** 0.5 + + +def omega_k(r): + return kep_profile(r) / r + + +def cs_profile(r): + cs_in = (H_r_0 * r0) * omega_k(r0) + return ((r / r0) ** (-q)) * cs_in + + +# Deduced quantities +pmass = disc_mass / Npart + +bsize = rout * 2 +bmin = (-bsize, -bsize, -bsize) +bmax = (bsize, bsize, bsize) + +cs0 = cs_profile(r0) + + +def rot_profile(r): + return ((kep_profile(r) ** 2) - (2 * p + q) * cs_profile(r) ** 2) ** 0.5 + + +def H_profile(r): + H = cs_profile(r) / omega_k(r) + # fact = (2.**0.5) * 3. # factor taken from phantom, to fasten thermalizing + fact = 1.0 + return fact * H + + +def simulate_disc(center_racc, is_torque_free, sim_folder): + + dump_folder = sim_folder + "dump/" + analysis_folder = sim_folder + "analysis/" + plot_folder = analysis_folder + "plots/" + + dump_prefix = dump_folder + "dump_" + + # %% + # Create the dump directory if it does not exist + if shamrock.sys.world_rank() == 0: + os.makedirs(sim_folder, exist_ok=True) + os.makedirs(dump_folder, exist_ok=True) + os.makedirs(analysis_folder, exist_ok=True) + os.makedirs(plot_folder, exist_ok=True) + + # %% + # Utility functions and quantities deduced from the base one + + # %% + # Start the context + # The context holds the data of the code + # We then init the layout of the field (e.g. the list of fields used by the solver) + + ctx = shamrock.Context() + ctx.pdata_layout_new() + + # %% + # Attach a SPH model to the context + + model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M4") + + # %% + # Dump handling + dump_helper = shamrock.utils.dump.ShamrockDumpHandleHelper(model, dump_prefix) + + # %% + # Load the last dump if it exists, setup otherwise + def setup_model(): + global disc_mass + + # Generate the default config + cfg = model.gen_default_config() + cfg.set_artif_viscosity_ConstantDisc(alpha_u=alpha_u, alpha_AV=alpha_AV, beta_AV=beta_AV) + cfg.set_eos_locally_isothermalLP07(cs0=cs0, q=q, r0=r0) + + cfg.add_kill_sphere( + center=(0, 0, 0), radius=bsize + ) # kill particles outside the simulation box + + cfg.set_units(codeu) + cfg.set_particle_mass(pmass) + # Set the CFL + cfg.set_cfl_cour(C_cour) + cfg.set_cfl_force(C_force) + + # Enable this to debug the neighbor counts + # cfg.set_show_neigh_stats(True) + + # Standard way to set the smoothing length (e.g. Price et al. 2018) + cfg.set_smoothing_length_density_based() + + cfg.set_save_dt_to_fields(True) + + # Standard density based smoothing lenght but with a neighbor count limit + # Use it if you have large slowdowns due to giant particles + # I recommend to use it if you have a circumbinary discs as the issue is very likely to happen + # cfg.set_smoothing_length_density_based_neigh_lim(500) + + cfg.set_save_dt_to_fields(True) + + # Set the solver config to be the one stored in cfg + model.set_solver_config(cfg) + + # Print the solver config + model.get_current_config().print_status() + + # Init the scheduler & fields + model.init_scheduler(scheduler_split_val, scheduler_merge_val) + + # Set the simulation box size + model.resize_simulation_box(bmin, bmax) + + # Create the setup + + setup = model.get_setup() + gen_disc = setup.make_generator_disc_mc( + part_mass=pmass, + disc_mass=disc_mass, + r_in=rin, + r_out=rout, + sigma_profile=sigma_profile, + H_profile=H_profile, + rot_profile=rot_profile, + cs_profile=cs_profile, + random_seed=666, + ) + + # Print the dot graph of the setup + print(gen_disc.get_dot()) + + # Apply the setup + setup.apply_setup(gen_disc) + + # correct the momentum and barycenter of the disc to 0 + analysis_momentum = shamrock.model_sph.analysisTotalMomentum(model=model) + total_momentum = analysis_momentum.get_total_momentum() + + if shamrock.sys.world_rank() == 0: + print(f"disc momentum = {total_momentum}") + + model.apply_momentum_offset((-total_momentum[0], -total_momentum[1], -total_momentum[2])) + + # Correct the barycenter before adding the sink + 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])) + + total_momentum = shamrock.model_sph.analysisTotalMomentum(model=model).get_total_momentum() + + if shamrock.sys.world_rank() == 0: + print(f"disc momentum after correction = {total_momentum}") + + barycenter, disc_mass = shamrock.model_sph.analysisBarycenter(model=model).get_barycenter() + + if shamrock.sys.world_rank() == 0: + print(f"disc barycenter after correction = {barycenter}") + + if not np.allclose(total_momentum, 0.0): + raise RuntimeError("disc momentum is not 0") + if not np.allclose(barycenter, 0.0): + raise RuntimeError("disc barycenter is not 0") + + # now that the barycenter & momentum are 0, we can add the sink + model.add_sink( + center_mass, (0, 0, 0), (0, 0, 0), center_racc, is_torque_free=is_torque_free + ) + + # Run a single step to init the integrator and smoothing length of the particles + # Here the htolerance is the maximum factor of evolution of the smoothing length in each + # Smoothing length iterations, increasing it affect the performance negatively but increse the + # convergence rate of the smoothing length + # this is why we increase it temporely to 1.3 before lowering it back to 1.1 (default value) + # Note that both ``change_htolerances`` can be removed and it will work the same but would converge + # more slowly at the first timestep + + model.change_htolerances(coarse=1.3, fine=1.1) + model.timestep() + model.change_htolerances(coarse=1.1, fine=1.1) + + dump_helper.load_last_dump_or(setup_model) + + def save_analysis_data(filename, key, value, ianalysis): + """Helper to save analysis data to a JSON file.""" + if shamrock.sys.world_rank() == 0: + filepath = os.path.join(analysis_folder, filename) + try: + with open(filepath, "r") as fp: + data = json.load(fp) + except (FileNotFoundError, json.JSONDecodeError): + data = {key: []} + data[key] = data[key][:ianalysis] + data[key].append({"t": model.get_time(), key: value}) + with open(filepath, "w") as fp: + json.dump(data, fp, indent=4) + + perf_analysis = PerfHistory(model, analysis_folder, "perf_history") + + column_density_plot = ColumnDensityPlot( + model, + ext_r=rout * 1.5, + nx=1024, + ny=1024, + ex=(1, 0, 0), + ey=(0, 1, 0), + center=(0, 0, 0), + analysis_folder=analysis_folder, + analysis_prefix="rho_integ_normal", + ) + + profile_plot = AnalysisHelper( + analysis_folder=os.path.join(analysis_folder, "plots"), + analysis_prefix="density_profile", + ) + + def analysis(ianalysis): + column_density_plot.analysis_save(ianalysis) + + barycenter, disc_mass = shamrock.model_sph.analysisBarycenter(model=model).get_barycenter() + + total_momentum = shamrock.model_sph.analysisTotalMomentum(model=model).get_total_momentum() + angular_momentum = shamrock.model_sph.analysisAngularMomentum( + model=model + ).get_angular_momentum() + + potential_energy = shamrock.model_sph.analysisEnergyPotential( + model=model + ).get_potential_energy() + + kinetic_energy = shamrock.model_sph.analysisEnergyKinetic(model=model).get_kinetic_energy() + + save_analysis_data("barycenter.json", "barycenter", barycenter, ianalysis) + save_analysis_data("disc_mass.json", "disc_mass", disc_mass, ianalysis) + save_analysis_data("total_momentum.json", "total_momentum", total_momentum, ianalysis) + save_analysis_data("angular_momentum.json", "angular_momentum", angular_momentum, ianalysis) + save_analysis_data("potential_energy.json", "potential_energy", potential_energy, ianalysis) + save_analysis_data("kinetic_energy.json", "kinetic_energy", kinetic_energy, ianalysis) + + sinks = model.get_sinks() + save_analysis_data("sinks.json", "sinks", sinks, ianalysis) + + perf_analysis.analysis_save(ianalysis) + + #''' + rho_field = model.compute_field("rho", "f64") + hpart_field = model.compute_field("hpart", "f64") + + def internal(size: int, x: np.array, y: np.array, z: np.array) -> np.array: + r = np.sqrt(x**2 + y**2 + z**2) + return r + + if _HAS_NUMBA: + internal = njit(internal) + + def custom_getter(size: int, dic_out: dict) -> np.array: + return internal( + size, + dic_out["xyz"][:, 0], + dic_out["xyz"][:, 1], + dic_out["xyz"][:, 2], + ) + + r_field = model.compute_field("custom", "f64", custom_getter) + + print(rho_field, r_field) + + x_min = 0.25 + x_max = rout * 1.5 + x_min_log = np.log10(x_min) + x_max_log = np.log10(x_max) + + bin_edges_x1d = np.logspace(x_min_log, x_max_log, 2049) + + histo = shamrock.compute_histogram( + bin_edges=bin_edges_x1d, + x_field=r_field, + y_field=rho_field, + do_average=True, + ) + + histo_convolve = shamrock.compute_histogram_convolve_x( + bin_edges=bin_edges_x1d, + x_field=r_field, + y_field=rho_field, + size_field=hpart_field, + do_average=True, + ) + + bin_edges_x = np.logspace(x_min_log, x_max_log, 1025) + bin_edges_y = np.logspace(-6, -3, 1025) + histo_top = shamrock.compute_histogram_2d( + bin_edges_x=bin_edges_x, + bin_edges_y=bin_edges_y, + x_field=r_field, + y_field=rho_field, + ) + histo_2d = np.array(histo_top).reshape(len(bin_edges_x) - 1, len(bin_edges_y) - 1) + + data = { + "bin_edges_x1d": bin_edges_x1d, + "bin_edges_x": bin_edges_x, + "bin_edges_y": bin_edges_y, + "histo": histo, + "histo_convolve": histo_convolve, + "histo_2d": histo_2d, + "time": model.get_time(), + } + + profile_plot.analysis_save(ianalysis, data) + + # %% + # Evolve the simulation + model.solver_logs_reset_cumulated_step_time() + model.solver_logs_reset_step_count() + + t_start = model.get_time() + + idump = 0 + iplot = 0 + istop = 0 + for ttarg in t_stop: + if ttarg >= t_start: + model.evolve_until(ttarg) + + if istop % dump_freq_stop == 0: + dump_helper.write_dump(idump, purge_old_dumps=True, keep_first=1, keep_last=3) + + # dump = model.make_phantom_dump() + # dump.save_dump(get_ph_dump_name(idump)) + + analysis(iplot) + + if istop % dump_freq_stop == 0: + idump += 1 + + iplot += 1 + istop += 1 + + +center_racc = 0.8 +sim_folder1 = f"_to_trash/circular_disc_sink_no_torque_free_{Npart}/" +simulate_disc(center_racc, False, sim_folder1) + +center_racc = 0.8 +sim_folder2 = f"_to_trash/circular_disc_sink_torque_free_{Npart}/" +simulate_disc(center_racc, True, sim_folder2) + +center_racc = 0.25 +sim_folder3 = f"_to_trash/circular_disc_sink_no_torque_free_0.25_{Npart}/" +simulate_disc(center_racc, False, sim_folder3) + +center_racc = 0.25 +sim_folder4 = f"_to_trash/circular_disc_sink_torque_free_0.25_{Npart}/" +simulate_disc(center_racc, True, sim_folder4) + +import matplotlib.pyplot as plt + + +def render_profiles(sim_folder1, sim_folder2, sim_folder3, sim_folder4): + profile_plot1 = AnalysisHelper( + analysis_folder=os.path.join(sim_folder1, "analysis", "plots"), + analysis_prefix="density_profile", + ) + profile_plot2 = AnalysisHelper( + analysis_folder=os.path.join(sim_folder2, "analysis", "plots"), + analysis_prefix="density_profile", + ) + profile_plot3 = AnalysisHelper( + analysis_folder=os.path.join(sim_folder3, "analysis", "plots"), + analysis_prefix="density_profile", + ) + profile_plot4 = AnalysisHelper( + analysis_folder=os.path.join(sim_folder4, "analysis", "plots"), + analysis_prefix="density_profile", + ) + profile_list_analysis_id1 = profile_plot1.get_list_analysis_id() + profile_list_analysis_id2 = profile_plot2.get_list_analysis_id() + profile_list_analysis_id3 = profile_plot3.get_list_analysis_id() + profile_list_analysis_id4 = profile_plot4.get_list_analysis_id() + + for iplot in profile_list_analysis_id1: + has_iplot_data_1 = iplot in profile_list_analysis_id1 + has_iplot_data_2 = iplot in profile_list_analysis_id2 + has_iplot_data_3 = iplot in profile_list_analysis_id3 + has_iplot_data_4 = iplot in profile_list_analysis_id4 + + if ( + not has_iplot_data_1 + or not has_iplot_data_2 + or not has_iplot_data_3 + or not has_iplot_data_4 + ): + continue + + data1 = profile_plot1.load_analysis(iplot).item() + data2 = profile_plot2.load_analysis(iplot).item() + data3 = profile_plot3.load_analysis(iplot).item() + data4 = profile_plot4.load_analysis(iplot).item() + + time = data1["time"] + + histo_convolve_torque_free = data1["histo_convolve"] + histo_convolve_not_torque_free = data2["histo_convolve"] + histo_convolve_torque_free_025 = data3["histo_convolve"] + histo_convolve_not_torque_free_025 = data4["histo_convolve"] + + bin_edges_x1d = data1["bin_edges_x1d"] + bin_center = (bin_edges_x1d[:-1] + bin_edges_x1d[1:]) / 2 + + plt.figure(dpi=150) + plt.plot(bin_center, histo_convolve_torque_free, label="racc=0.8") + plt.plot(bin_center, histo_convolve_not_torque_free, label="racc=0.8 torque free") + plt.plot(bin_center, histo_convolve_torque_free_025, label="racc=0.25") + plt.plot(bin_center, histo_convolve_not_torque_free_025, label="racc=0.25 torque free") + + text = f"t = {time:0.3f}" + from matplotlib.offsetbox import AnchoredText + + anchored_text = AnchoredText(text, loc=1) + plt.gca().add_artist(anchored_text) + + plt.xscale("log") + plt.yscale("log") + plt.ylim(1e-6, 1e-3) + plt.xlim(0.1, 10) + plt.legend(loc="upper left") + plt.savefig(f"_to_trash/compare_torque_free_sink_{iplot:07}.png") + plt.close() + + +render_profiles(sim_folder1, sim_folder2, sim_folder3, sim_folder4) diff --git a/src/shamalgs/include/shamalgs/primitives/compute_histogram.hpp b/src/shamalgs/include/shamalgs/primitives/compute_histogram.hpp index 77f76e495a..580ba229b8 100644 --- a/src/shamalgs/include/shamalgs/primitives/compute_histogram.hpp +++ b/src/shamalgs/include/shamalgs/primitives/compute_histogram.hpp @@ -312,7 +312,7 @@ namespace shamalgs::primitives { Tbins edge_sup = is_valid_point ? bin_edge_sup[ibin] : Tbins{}; // for each thread this will the sum of all the - // "func(in_data[group_size*i + local_data]) for all i" + // "func(in_data[group_size*i + local_id]) for all i" T local_sum = 0; for (size_t i = 0; i < element_count; i += group_size) { diff --git a/src/shammodels/common/src/pyCommonUtils.cpp b/src/shammodels/common/src/pyCommonUtils.cpp index 4549f81278..5eded60715 100644 --- a/src/shammodels/common/src/pyCommonUtils.cpp +++ b/src/shammodels/common/src/pyCommonUtils.cpp @@ -301,6 +301,10 @@ Register_pymod(shammodelcommonlibinit) { u32 iy = shamalgs::primitives::binary_search_upper_bound( y_bins, 0, ny + 1, y_val); + if (ix >= nx || iy >= ny) { + return; + } + using atomic_ref_T = sycl::atomic_ref< u64, sycl::memory_order_relaxed, From 4445c0b303e54ef2dcc3050d1c6de824cf9b3417 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Tue, 7 Apr 2026 14:03:21 +0200 Subject: [PATCH 12/23] noicer --- examples/production/compare_torque_free.py | 834 +++++++++++++++++++++ 1 file changed, 834 insertions(+) create mode 100644 examples/production/compare_torque_free.py diff --git a/examples/production/compare_torque_free.py b/examples/production/compare_torque_free.py new file mode 100644 index 0000000000..b58ffcb5d8 --- /dev/null +++ b/examples/production/compare_torque_free.py @@ -0,0 +1,834 @@ +""" +Production run: Circular disc & central sink particle +===================================================== + +This example demonstrates how to run a smoothed particle hydrodynamics (SPH) +simulation of a circular disc orbiting around a central sink. + +The simulation models: + +- A central star with a given mass and accretion radius +- A gaseous disc with specified mass, inner/outer radii, and vertical structure +- Artificial viscosity for angular momentum transport +- Locally isothermal equation of state + +Also this simulation feature rolling dumps (see `purge_old_dumps` function) to save disk space. + +This example is the accumulation of 3 files in a single one to showcase the complete workflow. + +- The actual run script (runscript.py) +- Plot generation (make_plots.py) +- Animation from the plots (plot_to_gif.py) + +On a cluster or laptop, one can run the code as follows: + +.. code-block:: bash + + mpirun ./shamrock --sycl-cfg 0:0 --loglevel 1 --rscript runscript.py + + +then after the run is done (or while it is running), one can run the following to generate the plots: + +.. code-block:: bash + + python make_plots.py + + +""" + +# %% +# Runscript (runscript.py) +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# The runscript is the actual simulation with on the fly analysis & rolling dumps + +import glob +import json +import os # for makedirs + +import numpy as np + +try: + from numba import njit + + _HAS_NUMBA = True +except ImportError: + _HAS_NUMBA = False + + +import shamrock + +# If we use the shamrock executable to run this script instead of the python interpreter, +# we should not initialize the system as the shamrock executable needs to handle specific MPI logic +if not shamrock.sys.is_initialized(): + shamrock.change_loglevel(1) + shamrock.sys.init("0:0") + +# %% +# Setup units + +si = shamrock.UnitSystem() +sicte = shamrock.Constants(si) +codeu = shamrock.UnitSystem( + unit_time=sicte.year(), + unit_length=sicte.au(), + unit_mass=sicte.sol_mass(), +) +ucte = shamrock.Constants(codeu) +G = ucte.G() + +# %% +# List parameters + +# Resolution +Npart = int(os.environ.get("NPART")) + +# Domain decomposition parameters +scheduler_split_val = int(1.0e7) # split patches with more than 1e7 particles +scheduler_merge_val = scheduler_split_val // 16 + +# Dump and plot frequency and duration of the simulation +dump_freq_stop = 2 +plot_freq_stop = 1 + +dt_stop = 0.1 +nstop = 10 + +# The list of times at which the simulation will pause for analysis / dumping +t_stop = [i * dt_stop for i in range(nstop + 1)] + + +# Sink parameters +center_mass = 1.0 +center_racc = float(os.environ.get("CENTER_RACC")) +center_is_torque_free = bool(os.environ.get("CENTER_IS_TORQUE_FREE")) + +# Disc parameter +disc_mass = 0.01 # sol mass +rout = 20.0 # au +rin = 1.0 # au +H_r_0 = 0.05 +q = 0.5 +p = 3.0 / 2.0 +r0 = 1.0 + +# Viscosity parameter +alpha_AV = 1.0e-3 / 0.08 +alpha_u = 1.0 +beta_AV = 2.0 + +# Integrator parameters +C_cour = 0.3 +C_force = 0.25 + + +sim_folder = f"_to_trash/circular_disc_sink_{center_is_torque_free}_{center_racc}_{Npart}/" + +dump_folder = sim_folder + "dump/" +analysis_folder = sim_folder + "analysis/" +plot_folder = analysis_folder + "plots/" + +dump_prefix = dump_folder + "dump_" + + +# Disc profiles +def sigma_profile(r): + sigma_0 = 1.0 # We do not care as it will be renormalized + return sigma_0 * (r / r0) ** (-p) + + +def kep_profile(r): + return (G * center_mass / r) ** 0.5 + + +def omega_k(r): + return kep_profile(r) / r + + +def cs_profile(r): + cs_in = (H_r_0 * r0) * omega_k(r0) + return ((r / r0) ** (-q)) * cs_in + + +# %% +# Create the dump directory if it does not exist +if shamrock.sys.world_rank() == 0: + os.makedirs(sim_folder, exist_ok=True) + os.makedirs(dump_folder, exist_ok=True) + os.makedirs(analysis_folder, exist_ok=True) + os.makedirs(plot_folder, exist_ok=True) + +# %% +# Utility functions and quantities deduced from the base one + +# Deduced quantities +pmass = disc_mass / Npart + +bsize = rout * 2 +bmin = (-bsize, -bsize, -bsize) +bmax = (bsize, bsize, bsize) + +cs0 = cs_profile(r0) + + +def rot_profile(r): + return ((kep_profile(r) ** 2) - (2 * p + q) * cs_profile(r) ** 2) ** 0.5 + + +def H_profile(r): + H = cs_profile(r) / omega_k(r) + # fact = (2.**0.5) * 3. # factor taken from phantom, to fasten thermalizing + fact = 1.0 + return fact * H + + +# %% +# Start the context +# The context holds the data of the code +# We then init the layout of the field (e.g. the list of fields used by the solver) + +ctx = shamrock.Context() +ctx.pdata_layout_new() + +# %% +# Attach a SPH model to the context + +model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M4") + + +# %% +# Dump handling + + +def get_vtk_dump_name(idump): + return dump_prefix + f"{idump:07}" + ".vtk" + + +def get_ph_dump_name(idump): + return dump_prefix + f"{idump:07}" + ".phdump" + + +dump_helper = shamrock.utils.dump.ShamrockDumpHandleHelper(model, dump_prefix) + +# %% +# Load the last dump if it exists, setup otherwise + + +def setup_model(): + global disc_mass + + # Generate the default config + cfg = model.gen_default_config() + cfg.set_artif_viscosity_ConstantDisc(alpha_u=alpha_u, alpha_AV=alpha_AV, beta_AV=beta_AV) + cfg.set_eos_locally_isothermalLP07(cs0=cs0, q=q, r0=r0) + + cfg.add_kill_sphere(center=(0, 0, 0), radius=bsize) # kill particles outside the simulation box + + cfg.set_units(codeu) + cfg.set_particle_mass(pmass) + # Set the CFL + cfg.set_cfl_cour(C_cour) + cfg.set_cfl_force(C_force) + + # Enable this to debug the neighbor counts + # cfg.set_show_neigh_stats(True) + + # Standard way to set the smoothing length (e.g. Price et al. 2018) + cfg.set_smoothing_length_density_based() + + cfg.set_save_dt_to_fields(True) + + # Standard density based smoothing lenght but with a neighbor count limit + # Use it if you have large slowdowns due to giant particles + # I recommend to use it if you have a circumbinary discs as the issue is very likely to happen + # cfg.set_smoothing_length_density_based_neigh_lim(500) + + cfg.set_save_dt_to_fields(True) + + # Set the solver config to be the one stored in cfg + model.set_solver_config(cfg) + + # Print the solver config + model.get_current_config().print_status() + + # Init the scheduler & fields + model.init_scheduler(scheduler_split_val, scheduler_merge_val) + + # Set the simulation box size + model.resize_simulation_box(bmin, bmax) + + # Create the setup + + setup = model.get_setup() + gen_disc = setup.make_generator_disc_mc( + part_mass=pmass, + disc_mass=disc_mass, + r_in=rin, + r_out=rout, + sigma_profile=sigma_profile, + H_profile=H_profile, + rot_profile=rot_profile, + cs_profile=cs_profile, + random_seed=666, + ) + + # Print the dot graph of the setup + print(gen_disc.get_dot()) + + # Apply the setup + setup.apply_setup(gen_disc) + + # correct the momentum and barycenter of the disc to 0 + analysis_momentum = shamrock.model_sph.analysisTotalMomentum(model=model) + total_momentum = analysis_momentum.get_total_momentum() + + if shamrock.sys.world_rank() == 0: + print(f"disc momentum = {total_momentum}") + + model.apply_momentum_offset((-total_momentum[0], -total_momentum[1], -total_momentum[2])) + + # Correct the barycenter before adding the sink + 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])) + + total_momentum = shamrock.model_sph.analysisTotalMomentum(model=model).get_total_momentum() + + if shamrock.sys.world_rank() == 0: + print(f"disc momentum after correction = {total_momentum}") + + barycenter, disc_mass = shamrock.model_sph.analysisBarycenter(model=model).get_barycenter() + + if shamrock.sys.world_rank() == 0: + print(f"disc barycenter after correction = {barycenter}") + + if not np.allclose(total_momentum, 0.0): + raise RuntimeError("disc momentum is not 0") + if not np.allclose(barycenter, 0.0): + raise RuntimeError("disc barycenter is not 0") + + # now that the barycenter & momentum are 0, we can add the sink + model.add_sink( + center_mass, (0, 0, 0), (0, 0, 0), center_racc, is_torque_free=center_is_torque_free + ) + + # Run a single step to init the integrator and smoothing length of the particles + # Here the htolerance is the maximum factor of evolution of the smoothing length in each + # Smoothing length iterations, increasing it affect the performance negatively but increse the + # convergence rate of the smoothing length + # this is why we increase it temporely to 1.3 before lowering it back to 1.1 (default value) + # Note that both ``change_htolerances`` can be removed and it will work the same but would converge + # more slowly at the first timestep + + model.change_htolerances(coarse=1.3, fine=1.1) + model.timestep() + model.change_htolerances(coarse=1.1, fine=1.1) + + +dump_helper.load_last_dump_or(setup_model) + + +# %% +# On the fly analysis +def save_rho_integ(ext, arr_rho, iplot): + if shamrock.sys.world_rank() == 0: + metadata = {"extent": [-ext, ext, -ext, ext], "time": model.get_time()} + np.save(plot_folder + f"rho_integ_{iplot:07}.npy", arr_rho) + + with open(plot_folder + f"rho_integ_{iplot:07}.json", "w") as fp: + json.dump(metadata, fp) + + +def save_analysis_data(filename, key, value, ianalysis): + """Helper to save analysis data to a JSON file.""" + if shamrock.sys.world_rank() == 0: + filepath = os.path.join(analysis_folder, filename) + try: + with open(filepath, "r") as fp: + data = json.load(fp) + except (FileNotFoundError, json.JSONDecodeError): + data = {key: []} + data[key] = data[key][:ianalysis] + data[key].append({"t": model.get_time(), key: value}) + with open(filepath, "w") as fp: + json.dump(data, fp, indent=4) + + +from shamrock.utils.analysis import ( + AnalysisHelper, + ColumnDensityPlot, + PerfHistory, + SliceDtPart, +) + +perf_analysis = PerfHistory(model, analysis_folder, "perf_history") + +column_density_plot = ColumnDensityPlot( + model, + ext_r=rout * 1.5, + nx=1024, + ny=1024, + ex=(1, 0, 0), + ey=(0, 1, 0), + center=(0, 0, 0), + analysis_folder=analysis_folder, + analysis_prefix="rho_integ_normal", +) + +dt_part_slice_plot = SliceDtPart( + model, + ext_r=rout * 0.5 / (16.0 / 9.0), # aspect ratio of 16:9 + nx=1920, + ny=1080, + ex=(1, 0, 0), + ey=(0, 0, 1), + center=((rin + rout) / 2, 0, 0), + analysis_folder=analysis_folder, + analysis_prefix="dt_part_slice", +) + +profile_plots = AnalysisHelper( + analysis_folder=os.path.join(analysis_folder, "plots"), + analysis_prefix="density_profile", +) + + +def analysis(ianalysis): + column_density_plot.analysis_save(ianalysis) + dt_part_slice_plot.analysis_save(ianalysis) + + barycenter, disc_mass = shamrock.model_sph.analysisBarycenter(model=model).get_barycenter() + + total_momentum = shamrock.model_sph.analysisTotalMomentum(model=model).get_total_momentum() + angular_momentum = shamrock.model_sph.analysisAngularMomentum( + model=model + ).get_angular_momentum() + + potential_energy = shamrock.model_sph.analysisEnergyPotential( + model=model + ).get_potential_energy() + + kinetic_energy = shamrock.model_sph.analysisEnergyKinetic(model=model).get_kinetic_energy() + + save_analysis_data("barycenter.json", "barycenter", barycenter, ianalysis) + save_analysis_data("disc_mass.json", "disc_mass", disc_mass, ianalysis) + save_analysis_data("total_momentum.json", "total_momentum", total_momentum, ianalysis) + save_analysis_data("angular_momentum.json", "angular_momentum", angular_momentum, ianalysis) + save_analysis_data("potential_energy.json", "potential_energy", potential_energy, ianalysis) + save_analysis_data("kinetic_energy.json", "kinetic_energy", kinetic_energy, ianalysis) + + sinks = model.get_sinks() + save_analysis_data("sinks.json", "sinks", sinks, ianalysis) + + perf_analysis.analysis_save(ianalysis) + + rho_field = model.compute_field("rho", "f64") + hpart_field = model.compute_field("hpart", "f64") + + def compute_r_field(): + + def internal(size: int, x: np.array, y: np.array, z: np.array) -> np.array: + r = np.sqrt(x**2 + y**2) + return r + + if _HAS_NUMBA: + internal = njit(internal) + + def custom_getter(size: int, dic_out: dict) -> np.array: + return internal( + size, + dic_out["xyz"][:, 0], + dic_out["xyz"][:, 1], + dic_out["xyz"][:, 2], + ) + + return model.compute_field("custom", "f64", custom_getter) + + def compute_Lz_field(): + pmass = model.get_particle_mass() + hfact = model.get_hfact() + + def internal( + size: int, + x: np.array, + y: np.array, + z: np.array, + hpart: np.array, + vx: np.array, + vy: np.array, + vz: np.array, + ) -> np.array: + rho = pmass * (hfact / hpart) ** 3 + Lz = rho * (x * vy - y * vx) + return Lz + + if _HAS_NUMBA: + internal = njit(internal) + + def custom_getter(size: int, dic_out: dict) -> np.array: + return internal( + size, + dic_out["xyz"][:, 0], + dic_out["xyz"][:, 1], + dic_out["xyz"][:, 2], + dic_out["hpart"], + dic_out["vxyz"][:, 0], + dic_out["vxyz"][:, 1], + dic_out["vxyz"][:, 2], + ) + + return model.compute_field("custom", "f64", custom_getter) + + r_field = compute_r_field() + Lz_field = compute_Lz_field() + + print(rho_field, r_field, Lz_field) + + x_min = 0.25 + x_max = rout * 1.5 + x_min_log = np.log10(x_min) + x_max_log = np.log10(x_max) + + bin_edges_x1d = np.logspace(x_min_log, x_max_log, 2049) + + rho_profile = shamrock.compute_histogram_convolve_x( + bin_edges=bin_edges_x1d, + x_field=r_field, + y_field=rho_field, + size_field=hpart_field, + do_average=True, + ) + + Lz_profile = shamrock.compute_histogram_convolve_x( + bin_edges=bin_edges_x1d, + x_field=r_field, + y_field=Lz_field, + size_field=hpart_field, + do_average=True, + ) + + data = { + "bin_edges_x1d": bin_edges_x1d, + "rho_profile": rho_profile, + "Lz_profile": Lz_profile, + "time": model.get_time(), + } + + profile_plots.analysis_save(ianalysis, data) + + +# %% +# Evolve the simulation +model.solver_logs_reset_cumulated_step_time() +model.solver_logs_reset_step_count() + +t_start = model.get_time() + +idump = 0 +iplot = 0 +istop = 0 +for ttarg in t_stop: + if ttarg >= t_start: + model.evolve_until(ttarg) + + if istop % dump_freq_stop == 0: + dump_helper.write_dump(idump, purge_old_dumps=True, keep_first=1, keep_last=3) + + # dump = model.make_phantom_dump() + # dump.save_dump(get_ph_dump_name(idump)) + + if istop % plot_freq_stop == 0: + analysis(iplot) + + if istop % dump_freq_stop == 0: + idump += 1 + + if istop % plot_freq_stop == 0: + iplot += 1 + + istop += 1 + +# %% +# Plot generation +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# Load the on-the-fly analysis after the run to make the plots +# (everything in this section can be in another file) + +import matplotlib +import matplotlib.pyplot as plt + +face_on_render_kwargs = { + "x_unit": "au", + "y_unit": "au", + "time_unit": "year", + "x_label": "x", + "y_label": "y", +} + +sink_params = { + "sink_scale_factor": 1, + "sink_color": "green", + "sink_linewidth": 1, + "sink_fill": False, +} + +column_density_plot.render_all( + **face_on_render_kwargs, + field_unit="kg.m^-2", + field_label="$\\int \\rho \\, \\mathrm{{d}} z$", + vmin=1, + vmax=1e4, + norm="log", + **sink_params, +) + +dt_part_slice_plot.render_all( + **face_on_render_kwargs, + field_unit="year", + field_label="$\\Delta t$", + vmin=1e-4, + vmax=1, + norm="log", + contour_list=[1e-4, 1e-3, 1e-2, 1e-1, 1], + **sink_params, +) + + +def profile_plots_func(iplot, data): + + data = data.item() + + bin_edges_x1d = data["bin_edges_x1d"] + rho_profile = data["rho_profile"] + Lz_profile = data["Lz_profile"] + time = data["time"] + + bin_center = (bin_edges_x1d[:-1] + bin_edges_x1d[1:]) / 2 + + plt.figure(dpi=150) + plt.plot(bin_center, rho_profile) + plt.xlabel("r") + plt.ylabel("density") + + plt.xscale("log") + plt.yscale("log") + + text = f"t = {time:0.3f}" + from matplotlib.offsetbox import AnchoredText + + anchored_text = AnchoredText(text, loc=2) + plt.gca().add_artist(anchored_text) + + plt.savefig(profile_plots.analysis_prefix + f"rho_{iplot:07}.png") + plt.close() + + plt.figure(dpi=150) + plt.plot(bin_center, Lz_profile) + plt.xlabel("r") + plt.ylabel("Lz") + + plt.xscale("log") + plt.yscale("log") + + text = f"t = {time:0.3f}" + from matplotlib.offsetbox import AnchoredText + + anchored_text = AnchoredText(text, loc=2) + plt.gca().add_artist(anchored_text) + + plt.savefig(profile_plots.analysis_prefix + f"Lz_{iplot:07}.png") + plt.close() + + +profile_plots.render_all(profile_plots_func) + +# %% +# Make gif for the doc (plot_to_gif.py) +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# Convert PNG sequence to Image sequence in mpl + +# sphinx_gallery_multi_image = "single" + + +render_gif = True + +# %% +# Do it for rho integ +if render_gif: + ani = column_density_plot.render_gif(gif_filename="rho_integ.gif", save_animation=True) + if ani is not None: + plt.show() + + +# %% +# Make a gif from the plots +if render_gif and shamrock.sys.world_rank() == 0: + ani = dt_part_slice_plot.render_gif(gif_filename="dt_part_slice.gif", save_animation=True) + if ani is not None: + plt.show() + + +# %% +# helper function to load data from JSON files +def load_data_from_json(filename, key): + filepath = os.path.join(analysis_folder, filename) + with open(filepath, "r") as fp: + data = json.load(fp)[key] + t = [d["t"] for d in data] + values = [d[key] for d in data] + return t, values + + +# %% +# load the json file for barycenter +t, barycenter = load_data_from_json("barycenter.json", "barycenter") +barycenter_x = [d[0] for d in barycenter] +barycenter_y = [d[1] for d in barycenter] +barycenter_z = [d[2] for d in barycenter] + +plt.figure(figsize=(8, 5), dpi=200) + +plt.plot(t, barycenter_x) +plt.plot(t, barycenter_y) +plt.plot(t, barycenter_z) +plt.xlabel("t") +plt.ylabel("barycenter") +plt.legend(["x", "y", "z"]) +plt.savefig(analysis_folder + "barycenter.png") +plt.show() + +# %% +# load the json file for disc_mass +t, disc_mass = load_data_from_json("disc_mass.json", "disc_mass") + +plt.figure(figsize=(8, 5), dpi=200) + +plt.plot(t, disc_mass) +plt.xlabel("t") +plt.ylabel("disc_mass") +plt.savefig(analysis_folder + "disc_mass.png") +plt.show() + +# %% +# load the json file for total_momentum +t, total_momentum = load_data_from_json("total_momentum.json", "total_momentum") +total_momentum_x = [d[0] for d in total_momentum] +total_momentum_y = [d[1] for d in total_momentum] +total_momentum_z = [d[2] for d in total_momentum] + +plt.figure(figsize=(8, 5), dpi=200) + +plt.plot(t, total_momentum_x) +plt.plot(t, total_momentum_y) +plt.plot(t, total_momentum_z) +plt.xlabel("t") +plt.ylabel("total_momentum") +plt.legend(["x", "y", "z"]) +plt.savefig(analysis_folder + "total_momentum.png") +plt.show() + + +# %% +# load the json file for total_momentum +t, angular_momentum = load_data_from_json("angular_momentum.json", "angular_momentum") +angular_momentum_x = [d[0] - angular_momentum[0][0] for d in angular_momentum] +angular_momentum_y = [d[1] - angular_momentum[0][1] for d in angular_momentum] +angular_momentum_z = [d[2] - angular_momentum[0][2] for d in angular_momentum] + + +plt.figure(figsize=(8, 5), dpi=200) + +plt.plot(t, angular_momentum_x) +plt.plot(t, angular_momentum_y) +plt.plot(t, angular_momentum_z) +plt.xlabel("t") +plt.ylabel(r"$\mathrm{L} - \mathrm{L}(t=0)$") +plt.legend(["x", "y", "z"]) +plt.savefig(analysis_folder + "angular_momentum.png") + +# %% +# load the json file for energies +t, potential_energy = load_data_from_json("potential_energy.json", "potential_energy") +_, kinetic_energy = load_data_from_json("kinetic_energy.json", "kinetic_energy") + +total_energy = [p + k for p, k in zip(potential_energy, kinetic_energy)] + +plt.figure(figsize=(8, 5), dpi=200) +plt.plot(t, potential_energy) +plt.plot(t, kinetic_energy) +plt.plot(t, total_energy) +plt.xlabel("t") +plt.ylabel("energy") +plt.legend(["potential_energy", "kinetic_energy", "total_energy"]) +plt.savefig(analysis_folder + "energies.png") +plt.show() + +# %% +# load the json file for sinks +t, sinks = load_data_from_json("sinks.json", "sinks") + +sinks_x = [d[0]["pos"][0] for d in sinks] +sinks_y = [d[0]["pos"][1] for d in sinks] +sinks_z = [d[0]["pos"][2] for d in sinks] + +plt.figure(figsize=(8, 5), dpi=200) +plt.plot(t, sinks_x, label="sink 0 (x)") +plt.plot(t, sinks_y, label="sink 0 (y)") +plt.plot(t, sinks_z, label="sink 0 (z)") +plt.xlabel("t") +plt.ylabel("sink position") +plt.legend() +plt.savefig(analysis_folder + "sinks.png") +plt.show() + +# %% +# Sink angular momentum +t, sinks = load_data_from_json("sinks.json", "sinks") + +sinks_lx = np.array([d[0]["angular_momentum"][0] for d in sinks]) +sinks_ly = np.array([d[0]["angular_momentum"][1] for d in sinks]) +sinks_lz = np.array([d[0]["angular_momentum"][2] for d in sinks]) + + +plt.figure(figsize=(8, 5), dpi=200) +plt.plot(t, sinks_lx, label="sink 0 (l_x)") +plt.plot(t, sinks_ly, label="sink 0 (l_y)") +plt.plot(t, sinks_lz, label="sink 0 (l_z)") +plt.xlabel("t") +plt.ylabel("sink spin") +plt.legend() +plt.savefig(analysis_folder + "sink_angular_momentum.png") +plt.show() + +# %% +# Sink to barycenter distance +t, sinks = load_data_from_json("sinks.json", "sinks") +_, barycenter = load_data_from_json("barycenter.json", "barycenter") + +barycenter_x = np.array([d[0] for d in barycenter]) +barycenter_y = np.array([d[1] for d in barycenter]) +barycenter_z = np.array([d[2] for d in barycenter]) + +sinks_x = np.array([d[0]["pos"][0] for d in sinks]) +sinks_y = np.array([d[0]["pos"][1] for d in sinks]) +sinks_z = np.array([d[0]["pos"][2] for d in sinks]) + + +plt.figure(figsize=(8, 5), dpi=200) +plt.plot(t, sinks_x - barycenter_x, label="sink 0 (x)") +plt.plot(t, sinks_y - barycenter_y, label="sink 0 (y)") +plt.plot(t, sinks_z - barycenter_z, label="sink 0 (z)") +plt.xlabel("t") +plt.ylabel("sink pos - barycenter pos") +plt.legend() +plt.savefig(analysis_folder + "sink_to_barycenter_distance.png") +plt.show() + + +# %% +# Plot the performance history (Switch close_plots to True if doing a long run) +perf_analysis.plot_perf_history(close_plots=False) +plt.show() From 86550172b7eba84e58bab074f927c68676207951 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Tue, 7 Apr 2026 14:11:08 +0200 Subject: [PATCH 13/23] more orbits --- examples/production/compare_torque_free.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/production/compare_torque_free.py b/examples/production/compare_torque_free.py index b58ffcb5d8..349f337e3d 100644 --- a/examples/production/compare_torque_free.py +++ b/examples/production/compare_torque_free.py @@ -90,8 +90,8 @@ dump_freq_stop = 2 plot_freq_stop = 1 -dt_stop = 0.1 -nstop = 10 +dt_stop = 1 +nstop = 2000 # The list of times at which the simulation will pause for analysis / dumping t_stop = [i * dt_stop for i in range(nstop + 1)] From 1c495c140c4bb6a8077eb496602713af3edef546 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Tue, 7 Apr 2026 18:22:16 +0200 Subject: [PATCH 14/23] more orbits --- examples/production/compare_torque_free.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/production/compare_torque_free.py b/examples/production/compare_torque_free.py index 349f337e3d..95f7f47402 100644 --- a/examples/production/compare_torque_free.py +++ b/examples/production/compare_torque_free.py @@ -80,7 +80,7 @@ # List parameters # Resolution -Npart = int(os.environ.get("NPART")) +Npart = int(float(os.environ.get("NPART"))) # Domain decomposition parameters scheduler_split_val = int(1.0e7) # split patches with more than 1e7 particles From 4be91bcc320b2fde05da89cf4150c1cc72b84666 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Tue, 7 Apr 2026 18:50:19 +0200 Subject: [PATCH 15/23] fix --- examples/production/compare_torque_free.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/production/compare_torque_free.py b/examples/production/compare_torque_free.py index 95f7f47402..944cb38024 100644 --- a/examples/production/compare_torque_free.py +++ b/examples/production/compare_torque_free.py @@ -100,7 +100,7 @@ # Sink parameters center_mass = 1.0 center_racc = float(os.environ.get("CENTER_RACC")) -center_is_torque_free = bool(os.environ.get("CENTER_IS_TORQUE_FREE")) +center_is_torque_free = bool(os.environ.get("CENTER_IS_TORQUE_FREE").lower() == "true") # Disc parameter disc_mass = 0.01 # sol mass @@ -123,6 +123,8 @@ sim_folder = f"_to_trash/circular_disc_sink_{center_is_torque_free}_{center_racc}_{Npart}/" +print(sim_folder) + dump_folder = sim_folder + "dump/" analysis_folder = sim_folder + "analysis/" plot_folder = analysis_folder + "plots/" From 573a5b7095cadd7cc908a900d60992d00f4437db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Thu, 9 Apr 2026 14:42:18 +0200 Subject: [PATCH 16/23] add analysis --- examples/production/make_plots.py | 80 +++++++++++++++++++++++++++++++ 1 file changed, 80 insertions(+) create mode 100644 examples/production/make_plots.py diff --git a/examples/production/make_plots.py b/examples/production/make_plots.py new file mode 100644 index 0000000000..66ae041503 --- /dev/null +++ b/examples/production/make_plots.py @@ -0,0 +1,80 @@ +cases = { + "racc=0.8": { + "path": "_to_trash/circular_disc_sink_False_0.8_10000000", + "is_torque_free": False, + "racc": 0.8, + }, + "racc=0.25": { + "path": "_to_trash/circular_disc_sink_False_0.25_10000000", + "is_torque_free": False, + "racc": 0.25, + }, + "racc=0.8 (torque-free)": { + "path": "_to_trash/circular_disc_sink_True_0.8_10000000", + "is_torque_free": True, + "racc": 0.8, + }, + "racc=0.25 (torque-free)": { + "path": "_to_trash/circular_disc_sink_True_0.25_10000000", + "is_torque_free": True, + "racc": 0.25, + }, +} + +import os + +import matplotlib.pyplot as plt +from shamrock.utils.analysis import ( + AnalysisHelper, +) + + +def render_profiles(cases): + helpers = {} + for key, case in cases.items(): + helpers[case] = {} + helpers[case]["helper"] = AnalysisHelper( + analysis_folder=os.path.join(cases[case]["path"], "analysis", "plots"), + analysis_prefix="density_profile", + ) + + for k in helpers: + helpers[k]["list_analysis_id"] = helpers[k]["helper"].get_list_analysis_id() + + profile_list_analysis_id1 = helpers["racc=0.8"]["list_analysis_id"] + + for iplot in profile_list_analysis_id1: + has_iplot_data = True + for k in helpers: + has_iplot_data = has_iplot_data and iplot in helpers[k]["list_analysis_id"] + if not has_iplot_data: + continue + + data = {} + for k in helpers: + data[k] = helpers[k]["helper"].load_analysis(iplot).item() + + time = data["racc=0.8"]["time"] + bin_edges_x1d = data["racc=0.8"]["bin_edges_x1d"] + bin_center = (bin_edges_x1d[:-1] + bin_edges_x1d[1:]) / 2 + + plt.figure(dpi=150) + for k, v in data.items(): + plt.plot(bin_center, v["histo_convolve"], label=k) + + text = f"t = {time:0.3f}" + from matplotlib.offsetbox import AnchoredText + + anchored_text = AnchoredText(text, loc=1) + plt.gca().add_artist(anchored_text) + + plt.xscale("log") + plt.yscale("log") + plt.ylim(1e-6, 1e-3) + plt.xlim(0.1, 10) + plt.legend(loc="upper left") + plt.savefig(f"_to_trash/compare_torque_free_sink_{iplot:07}.png") + plt.close() + + +render_profiles(cases) From 598c9b27352bf3bdcca680dcc644dc898f8131be Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Mon, 13 Apr 2026 15:44:46 +0200 Subject: [PATCH 17/23] update --- examples/production/make_plots.py | 171 ++++++++++++++++++++++++++---- 1 file changed, 153 insertions(+), 18 deletions(-) diff --git a/examples/production/make_plots.py b/examples/production/make_plots.py index 66ae041503..c858d620fd 100644 --- a/examples/production/make_plots.py +++ b/examples/production/make_plots.py @@ -4,6 +4,11 @@ "is_torque_free": False, "racc": 0.8, }, + "racc=0.5": { + "path": "_to_trash/circular_disc_sink_False_0.5_10000000", + "is_torque_free": False, + "racc": 0.5, + }, "racc=0.25": { "path": "_to_trash/circular_disc_sink_False_0.25_10000000", "is_torque_free": False, @@ -14,6 +19,11 @@ "is_torque_free": True, "racc": 0.8, }, + "racc=0.5 (torque-free)": { + "path": "_to_trash/circular_disc_sink_True_0.5_10000000", + "is_torque_free": True, + "racc": 0.5, + }, "racc=0.25 (torque-free)": { "path": "_to_trash/circular_disc_sink_True_0.25_10000000", "is_torque_free": True, @@ -21,20 +31,21 @@ }, } +import json import os - +import numpy as np import matplotlib.pyplot as plt from shamrock.utils.analysis import ( AnalysisHelper, ) -def render_profiles(cases): +def render_rho_profiles(cases): helpers = {} for key, case in cases.items(): - helpers[case] = {} - helpers[case]["helper"] = AnalysisHelper( - analysis_folder=os.path.join(cases[case]["path"], "analysis", "plots"), + helpers[key] = {} + helpers[key]["helper"] = AnalysisHelper( + analysis_folder=os.path.join(cases[key]["path"], "analysis", "plots"), analysis_prefix="density_profile", ) @@ -44,23 +55,69 @@ def render_profiles(cases): profile_list_analysis_id1 = helpers["racc=0.8"]["list_analysis_id"] for iplot in profile_list_analysis_id1: - has_iplot_data = True + + ref_key = None + data = {} for k in helpers: - has_iplot_data = has_iplot_data and iplot in helpers[k]["list_analysis_id"] - if not has_iplot_data: - continue + if iplot in helpers[k]["list_analysis_id"]: + data[k] = helpers[k]["helper"].load_analysis(iplot).item() + if ref_key is None: + ref_key = k + + time = data[ref_key]["time"] + bin_edges_x1d = data[ref_key]["bin_edges_x1d"] + bin_center = (bin_edges_x1d[:-1] + bin_edges_x1d[1:]) / 2 + + plt.figure(dpi=150) + for k, v in data.items(): + plt.plot(bin_center, v["rho_profile"], label=k) + + text = f"t = {time:0.3f}" + from matplotlib.offsetbox import AnchoredText + + anchored_text = AnchoredText(text, loc=1) + plt.gca().add_artist(anchored_text) + + plt.xscale("log") + plt.yscale("log") + #plt.ylim(1e-6, 1e-3) + #plt.xlim(0.1, 21) + plt.legend(loc="lower left") + plt.savefig(f"_to_trash/compare_torque_free_sink_rho_profile_{iplot:07}.png") + plt.close() + + +def render_Lz_profiles(cases): + helpers = {} + for key, case in cases.items(): + helpers[key] = {} + helpers[key]["helper"] = AnalysisHelper( + analysis_folder=os.path.join(cases[key]["path"], "analysis", "plots"), + analysis_prefix="density_profile", + ) + + for k in helpers: + helpers[k]["list_analysis_id"] = helpers[k]["helper"].get_list_analysis_id() + + profile_list_analysis_id1 = helpers["racc=0.8"]["list_analysis_id"] + for iplot in profile_list_analysis_id1: + + ref_key = None data = {} for k in helpers: - data[k] = helpers[k]["helper"].load_analysis(iplot).item() + if iplot in helpers[k]["list_analysis_id"]: + data[k] = helpers[k]["helper"].load_analysis(iplot).item() + if ref_key is None: + ref_key = k - time = data["racc=0.8"]["time"] - bin_edges_x1d = data["racc=0.8"]["bin_edges_x1d"] + time = data[ref_key]["time"] + bin_edges_x1d = data[ref_key]["bin_edges_x1d"] bin_center = (bin_edges_x1d[:-1] + bin_edges_x1d[1:]) / 2 plt.figure(dpi=150) for k, v in data.items(): - plt.plot(bin_center, v["histo_convolve"], label=k) + plt.plot(bin_center, v["Lz_profile"], label=k) text = f"t = {time:0.3f}" from matplotlib.offsetbox import AnchoredText @@ -70,11 +127,89 @@ def render_profiles(cases): plt.xscale("log") plt.yscale("log") - plt.ylim(1e-6, 1e-3) - plt.xlim(0.1, 10) - plt.legend(loc="upper left") - plt.savefig(f"_to_trash/compare_torque_free_sink_{iplot:07}.png") + #plt.ylim(1e-6, 1e-3) + #plt.xlim(0.1, 21) + plt.legend(loc="lower left") + plt.savefig(f"_to_trash/compare_torque_free_sink_Lz_profile_{iplot:07}.png") plt.close() -render_profiles(cases) +def compare_acc_rate(cases): + curves = {} + for key, case in cases.items(): + disc_mass_filename = os.path.join(cases[key]["path"], "analysis", "disc_mass.json") + + ret = {} + with open(disc_mass_filename, "r") as fp: + disc_mass_data = json.load(fp)["disc_mass"] + ret["t"] = [d["t"] for d in disc_mass_data] + ret["disc_mass"] = [d["disc_mass"] for d in disc_mass_data] + + # time derivative of the disc mass + ret["disc_mass_rate"] = np.diff(ret["disc_mass"]) / np.diff(ret["t"]) + curves[key] = ret + + + + plt.figure(dpi=150) + for k, v in curves.items(): + plt.plot(v["t"], v["disc_mass"], label=k) + + #plt.xscale("log") + #plt.yscale("log") + #plt.ylim(1e-6, 1e-3) + #plt.xlim(0.1, 21) + plt.legend(loc="lower left") + plt.savefig(f"_to_trash/compare_torque_free_sink_acc_rate_disc_mass.png") + plt.close() + + plt.figure(dpi=150) + for k, v in curves.items(): + plt.plot(v["t"][1:], v["disc_mass_rate"], label=k) + + #plt.xscale("log") + #plt.yscale("log") + #plt.ylim(1e-6, 1e-3) + #plt.xlim(0.1, 21) + plt.legend(loc="lower right") + plt.savefig(f"_to_trash/compare_torque_free_sink_acc_rate_diff.png") + plt.close() + +def compare_dsink_rate(cases): + + def get_r(pos): + return np.linalg.norm(pos) + + curves = {} + for key, case in cases.items(): + disc_mass_filename = os.path.join(cases[key]["path"], "analysis", "sinks.json") + + ret = {} + with open(disc_mass_filename, "r") as fp: + sinks_dat = json.load(fp)["sinks"] + ret["t"] = [d["t"] for d in sinks_dat] + ret["r"] = [get_r(d["sinks"][0]["pos"]) for d in sinks_dat] + + curves[key] = ret + + + + plt.figure(dpi=150) + for k, v in curves.items(): + plt.plot(v["t"], v["r"], label=k) + + #plt.xscale("log") + #plt.yscale("log") + #plt.ylim(1e-6, 1e-3) + #plt.xlim(0.1, 21) + plt.legend(loc="lower left") + plt.savefig(f"_to_trash/compare_torque_free_sink_r.png") + plt.close() + + + + +compare_dsink_rate(cases) +compare_acc_rate(cases) +render_rho_profiles(cases) +render_Lz_profiles(cases) \ No newline at end of file From 4c44fa6783d6d7c69ecd0373940b62c7281ff26c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Sat, 18 Apr 2026 10:45:23 +0200 Subject: [PATCH 18/23] boost radisu is now a param --- examples/production/compare_torque_free.py | 31 ++++++++-- examples/production/make_plots.py | 58 +++++++++---------- .../sph/include/shammodels/sph/Model.hpp | 26 ++++++++- .../include/shammodels/sph/SinkPartStruct.hpp | 3 + src/shammodels/sph/src/Model.cpp | 3 +- .../sph/src/modules/SinkParticlesUpdate.cpp | 3 +- src/shammodels/sph/src/pySPHModel.cpp | 20 ++++--- 7 files changed, 93 insertions(+), 51 deletions(-) diff --git a/examples/production/compare_torque_free.py b/examples/production/compare_torque_free.py index 944cb38024..61aad5e2bd 100644 --- a/examples/production/compare_torque_free.py +++ b/examples/production/compare_torque_free.py @@ -91,16 +91,32 @@ plot_freq_stop = 1 dt_stop = 1 -nstop = 2000 +nstop = 100 # The list of times at which the simulation will pause for analysis / dumping t_stop = [i * dt_stop for i in range(nstop + 1)] +run_params = [] + +for center_is_torque_free in [True, False]: + for center_racc in [0.25, 0.5, 0.8]: + if center_is_torque_free: + for center_torque_boost_radius_fact in [1.5, 2.0, 3.0]: + run_params.append( + (center_racc, center_is_torque_free, center_torque_boost_radius_fact) + ) + else: + run_params.append((center_racc, center_is_torque_free, 0.0)) + +print(run_params, len(run_params)) + +run_id = os.environ.get("RUN_ID") # Sink parameters center_mass = 1.0 -center_racc = float(os.environ.get("CENTER_RACC")) -center_is_torque_free = bool(os.environ.get("CENTER_IS_TORQUE_FREE").lower() == "true") +center_racc = run_params[run_id][0] +center_is_torque_free = run_params[run_id][1] +center_torque_boost_radius_fact = run_params[run_id][2] # Disc parameter disc_mass = 0.01 # sol mass @@ -121,7 +137,7 @@ C_force = 0.25 -sim_folder = f"_to_trash/circular_disc_sink_{center_is_torque_free}_{center_racc}_{Npart}/" +sim_folder = f"_to_trash/circular_disc_sink_{center_is_torque_free}_{center_racc}_{center_torque_boost_radius_fact}_{Npart}/" print(sim_folder) @@ -314,7 +330,12 @@ def setup_model(): # now that the barycenter & momentum are 0, we can add the sink model.add_sink( - center_mass, (0, 0, 0), (0, 0, 0), center_racc, is_torque_free=center_is_torque_free + center_mass, + (0, 0, 0), + (0, 0, 0), + center_racc, + is_torque_free=center_is_torque_free, + torque_boost_radius_fact=center_torque_boost_radius_fact, ) # Run a single step to init the integrator and smoothing length of the particles diff --git a/examples/production/make_plots.py b/examples/production/make_plots.py index c858d620fd..f4736b6896 100644 --- a/examples/production/make_plots.py +++ b/examples/production/make_plots.py @@ -33,8 +33,9 @@ import json import os -import numpy as np + import matplotlib.pyplot as plt +import numpy as np from shamrock.utils.analysis import ( AnalysisHelper, ) @@ -55,7 +56,6 @@ def render_rho_profiles(cases): profile_list_analysis_id1 = helpers["racc=0.8"]["list_analysis_id"] for iplot in profile_list_analysis_id1: - ref_key = None data = {} for k in helpers: @@ -80,8 +80,8 @@ def render_rho_profiles(cases): plt.xscale("log") plt.yscale("log") - #plt.ylim(1e-6, 1e-3) - #plt.xlim(0.1, 21) + # plt.ylim(1e-6, 1e-3) + # plt.xlim(0.1, 21) plt.legend(loc="lower left") plt.savefig(f"_to_trash/compare_torque_free_sink_rho_profile_{iplot:07}.png") plt.close() @@ -102,7 +102,6 @@ def render_Lz_profiles(cases): profile_list_analysis_id1 = helpers["racc=0.8"]["list_analysis_id"] for iplot in profile_list_analysis_id1: - ref_key = None data = {} for k in helpers: @@ -127,8 +126,8 @@ def render_Lz_profiles(cases): plt.xscale("log") plt.yscale("log") - #plt.ylim(1e-6, 1e-3) - #plt.xlim(0.1, 21) + # plt.ylim(1e-6, 1e-3) + # plt.xlim(0.1, 21) plt.legend(loc="lower left") plt.savefig(f"_to_trash/compare_torque_free_sink_Lz_profile_{iplot:07}.png") plt.close() @@ -138,7 +137,7 @@ def compare_acc_rate(cases): curves = {} for key, case in cases.items(): disc_mass_filename = os.path.join(cases[key]["path"], "analysis", "disc_mass.json") - + ret = {} with open(disc_mass_filename, "r") as fp: disc_mass_data = json.load(fp)["disc_mass"] @@ -149,32 +148,31 @@ def compare_acc_rate(cases): ret["disc_mass_rate"] = np.diff(ret["disc_mass"]) / np.diff(ret["t"]) curves[key] = ret - - plt.figure(dpi=150) for k, v in curves.items(): plt.plot(v["t"], v["disc_mass"], label=k) - #plt.xscale("log") - #plt.yscale("log") - #plt.ylim(1e-6, 1e-3) - #plt.xlim(0.1, 21) + # plt.xscale("log") + # plt.yscale("log") + # plt.ylim(1e-6, 1e-3) + # plt.xlim(0.1, 21) plt.legend(loc="lower left") - plt.savefig(f"_to_trash/compare_torque_free_sink_acc_rate_disc_mass.png") + plt.savefig("_to_trash/compare_torque_free_sink_acc_rate_disc_mass.png") plt.close() - + plt.figure(dpi=150) for k, v in curves.items(): plt.plot(v["t"][1:], v["disc_mass_rate"], label=k) - #plt.xscale("log") - #plt.yscale("log") - #plt.ylim(1e-6, 1e-3) - #plt.xlim(0.1, 21) + # plt.xscale("log") + # plt.yscale("log") + # plt.ylim(1e-6, 1e-3) + # plt.xlim(0.1, 21) plt.legend(loc="lower right") - plt.savefig(f"_to_trash/compare_torque_free_sink_acc_rate_diff.png") + plt.savefig("_to_trash/compare_torque_free_sink_acc_rate_diff.png") plt.close() + def compare_dsink_rate(cases): def get_r(pos): @@ -183,7 +181,7 @@ def get_r(pos): curves = {} for key, case in cases.items(): disc_mass_filename = os.path.join(cases[key]["path"], "analysis", "sinks.json") - + ret = {} with open(disc_mass_filename, "r") as fp: sinks_dat = json.load(fp)["sinks"] @@ -192,24 +190,20 @@ def get_r(pos): curves[key] = ret - - plt.figure(dpi=150) for k, v in curves.items(): plt.plot(v["t"], v["r"], label=k) - #plt.xscale("log") - #plt.yscale("log") - #plt.ylim(1e-6, 1e-3) - #plt.xlim(0.1, 21) + # plt.xscale("log") + # plt.yscale("log") + # plt.ylim(1e-6, 1e-3) + # plt.xlim(0.1, 21) plt.legend(loc="lower left") - plt.savefig(f"_to_trash/compare_torque_free_sink_r.png") + plt.savefig("_to_trash/compare_torque_free_sink_r.png") plt.close() - - compare_dsink_rate(cases) compare_acc_rate(cases) render_rho_profiles(cases) -render_Lz_profiles(cases) \ No newline at end of file +render_Lz_profiles(cases) diff --git a/src/shammodels/sph/include/shammodels/sph/Model.hpp b/src/shammodels/sph/include/shammodels/sph/Model.hpp index 4fa9296dc5..061c6d0002 100644 --- a/src/shammodels/sph/include/shammodels/sph/Model.hpp +++ b/src/shammodels/sph/include/shammodels/sph/Model.hpp @@ -174,16 +174,36 @@ namespace shammodels::sph { std::mt19937 eng); inline void add_sink( - Tscal mass, Tvec pos, Tvec velocity, Tscal accretion_radius, bool is_torque_free) { + Tscal mass, + Tvec pos, + Tvec velocity, + Tscal accretion_radius, + bool is_torque_free, + Tscal torque_boost_radius_fact) { if (solver.storage.sinks.is_empty()) { solver.storage.sinks.set({}); } shamlog_debug_ln( - "SPH", "add sink :", mass, pos, velocity, accretion_radius, is_torque_free); + "SPH", + "add sink :", + mass, + pos, + velocity, + accretion_radius, + is_torque_free, + torque_boost_radius_fact); solver.storage.sinks.get().push_back( - {pos, velocity, {}, {}, mass, {}, accretion_radius, is_torque_free}); + {pos, + velocity, + {}, + {}, + mass, + {}, + accretion_radius, + is_torque_free, + torque_boost_radius_fact}); } template diff --git a/src/shammodels/sph/include/shammodels/sph/SinkPartStruct.hpp b/src/shammodels/sph/include/shammodels/sph/SinkPartStruct.hpp index 23492deb0d..9eb0e95be7 100644 --- a/src/shammodels/sph/include/shammodels/sph/SinkPartStruct.hpp +++ b/src/shammodels/sph/include/shammodels/sph/SinkPartStruct.hpp @@ -36,6 +36,7 @@ namespace shammodels::sph { Tscal accretion_radius; bool is_torque_free; + Tscal torque_boost_radius_fact; }; template @@ -53,6 +54,7 @@ namespace shammodels::sph { {"angular_momentum", p.angular_momentum}, {"accretion_radius", p.accretion_radius}, {"is_torque_free", p.is_torque_free}, + {"torque_boost_radius_fact", p.torque_boost_radius_fact}, }; } @@ -69,6 +71,7 @@ namespace shammodels::sph { j.at("angular_momentum").get_to(p.angular_momentum); j.at("accretion_radius").get_to(p.accretion_radius); j.at("is_torque_free").get_to(p.is_torque_free); + j.at("torque_boost_radius_fact").get_to(p.torque_boost_radius_fact); } } // namespace shammodels::sph diff --git a/src/shammodels/sph/src/Model.cpp b/src/shammodels/sph/src/Model.cpp index bae04a108d..589515eb83 100644 --- a/src/shammodels/sph/src/Model.cpp +++ b/src/shammodels/sph/src/Model.cpp @@ -1438,7 +1438,8 @@ void shammodels::sph::Model::init_from_phantom_dump( {xsink[i], ysink[i], zsink[i]}, {vxsink[i], vysink[i], vzsink[i]}, Racc[i], - false); + false, + 2.0); } } } diff --git a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp index 27da533a12..29fe1457e4 100644 --- a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp +++ b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp @@ -342,7 +342,8 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par } }; - BoostWeight weight_func{s.accretion_radius * 2}; + BoostWeight weight_func{ + s.accretion_radius * s.torque_boost_radius_fact / Kernel::Rkern}; Tvec r_sink = s.pos; diff --git a/src/shammodels/sph/src/pySPHModel.cpp b/src/shammodels/sph/src/pySPHModel.cpp index 855579d0e7..6ac27434eb 100644 --- a/src/shammodels/sph/src/pySPHModel.cpp +++ b/src/shammodels/sph/src/pySPHModel.cpp @@ -832,7 +832,8 @@ void add_instance(py::module &m, std::string name_config, std::string name_model py::arg("pos"), py::arg("velocity"), py::arg("accretion_radius"), - py::arg("is_torque_free") = false) + py::arg("is_torque_free") = false, + py::arg("torque_boost_radius_fact") = 2.0) .def( "get_sinks", [](T &self) { @@ -841,14 +842,15 @@ void add_instance(py::module &m, std::string name_config, std::string name_model if (!self.solver.storage.sinks.is_empty()) { for (auto &sink : self.solver.storage.sinks.get()) { py::dict sink_dic; - sink_dic["pos"] = sink.pos; - sink_dic["velocity"] = sink.velocity; - sink_dic["sph_acceleration"] = sink.sph_acceleration; - sink_dic["ext_acceleration"] = sink.ext_acceleration; - sink_dic["mass"] = sink.mass; - sink_dic["angular_momentum"] = sink.angular_momentum; - sink_dic["accretion_radius"] = sink.accretion_radius; - sink_dic["is_torque_free"] = sink.is_torque_free; + sink_dic["pos"] = sink.pos; + sink_dic["velocity"] = sink.velocity; + sink_dic["sph_acceleration"] = sink.sph_acceleration; + sink_dic["ext_acceleration"] = sink.ext_acceleration; + sink_dic["mass"] = sink.mass; + sink_dic["angular_momentum"] = sink.angular_momentum; + sink_dic["accretion_radius"] = sink.accretion_radius; + sink_dic["is_torque_free"] = sink.is_torque_free; + sink_dic["torque_boost_radius_fact"] = sink.torque_boost_radius_fact; list_out.append(sink_dic); } } From ffb463385a7cc8b71e6dbb994dc2037bdc4a67e9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Sat, 18 Apr 2026 11:11:14 +0200 Subject: [PATCH 19/23] whoops --- examples/production/compare_torque_free.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/production/compare_torque_free.py b/examples/production/compare_torque_free.py index 61aad5e2bd..6fa64d4662 100644 --- a/examples/production/compare_torque_free.py +++ b/examples/production/compare_torque_free.py @@ -110,7 +110,7 @@ print(run_params, len(run_params)) -run_id = os.environ.get("RUN_ID") +run_id = int(os.environ.get("RUN_ID")) # Sink parameters center_mass = 1.0 From 81eafeedd7dfee3ce748d9dbc5c9c033d42f5360 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Sun, 19 Apr 2026 13:00:36 +0200 Subject: [PATCH 20/23] ensure validity of the boost vector --- .../sph/src/modules/SinkParticlesUpdate.cpp | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp index 29fe1457e4..72fad4c72b 100644 --- a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp +++ b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp @@ -17,6 +17,7 @@ #include "shambase/DistributedData.hpp" #include "shambase/narrowing.hpp" +#include "shambase/string.hpp" #include "shamalgs/collective/reduction.hpp" #include "shamalgs/details/numeric/numeric.hpp" #include "shamalgs/numeric.hpp" @@ -426,6 +427,24 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par Tvec vec_delta_w = {delta_w.data[0], delta_w.data[1], delta_w[2]}; + if (sham::has_nan_or_inf(vec_delta_w)) { + logger::warn_ln( + "Sink", + shambase::format( + "Nan or inf detected after inertia matrix inversion\n" + " delta_l = {}\n" + " delta_w = {}\n" + " weighted_intertia_sum = {}\n" + " I_inv = {}\n" + "I will therefore skip the boost step (delta_w = 0)", + delta_l, + delta_w, + weighted_intertia_sum.data, + I_inv.data)); + + delta_w = {}; + } + logger::raw_ln("delta l =", delta_l.data[0], delta_l.data[1], delta_l.data[2]); logger::raw_ln("delta w =", delta_w.data[0], delta_w.data[1], delta_w.data[2]); From 966e1f36144c3c527529f312c81814dc13619c14 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Sun, 19 Apr 2026 13:01:43 +0200 Subject: [PATCH 21/23] upload --- examples/production/compare_torque_free.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/production/compare_torque_free.py b/examples/production/compare_torque_free.py index 6fa64d4662..d46bbdc6fa 100644 --- a/examples/production/compare_torque_free.py +++ b/examples/production/compare_torque_free.py @@ -99,7 +99,7 @@ run_params = [] for center_is_torque_free in [True, False]: - for center_racc in [0.25, 0.5, 0.8]: + for center_racc in [0.1, 0.25, 0.5, 0.8]: if center_is_torque_free: for center_torque_boost_radius_fact in [1.5, 2.0, 3.0]: run_params.append( From 6bd43464c0f984522ccf58b8e313c1b0982e5852 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Sun, 19 Apr 2026 13:14:36 +0200 Subject: [PATCH 22/23] fix --- src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp index 72fad4c72b..a1c6a43119 100644 --- a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp +++ b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp @@ -24,6 +24,7 @@ #include "shamalgs/primitives/reduction.hpp" #include "shamalgs/reduction.hpp" #include "shambackends/DeviceBuffer.hpp" +#include "shambackends/fmt_bindings/fmt_defs.hpp" #include "shambackends/kernel_call.hpp" #include "shambackends/math.hpp" #include "shamcomm/logs.hpp" @@ -437,8 +438,8 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par " weighted_intertia_sum = {}\n" " I_inv = {}\n" "I will therefore skip the boost step (delta_w = 0)", - delta_l, - delta_w, + delta_l.data, + delta_w.data, weighted_intertia_sum.data, I_inv.data)); From 464fca14e68cff5d1b72f69e81274cb07582bbf4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Sun, 19 Apr 2026 13:47:49 +0200 Subject: [PATCH 23/23] fix --- src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp index a1c6a43119..f5b08d9fe1 100644 --- a/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp +++ b/src/shammodels/sph/src/modules/SinkParticlesUpdate.cpp @@ -443,7 +443,8 @@ void shammodels::sph::modules::SinkParticlesUpdate::accrete_par weighted_intertia_sum.data, I_inv.data)); - delta_w = {}; + delta_w = {}; + vec_delta_w = {delta_w.data[0], delta_w.data[1], delta_w[2]}; } logger::raw_ln("delta l =", delta_l.data[0], delta_l.data[1], delta_l.data[2]);