From 79a333b858db7daaf762114167871dc85302403e Mon Sep 17 00:00:00 2001 From: Jacksen Narvaez Date: Mon, 22 Jun 2026 22:40:20 -0230 Subject: [PATCH 1/4] fix: move sink loading outside particle insertion loop to prevent duplicate sinks being added once per range chunk --- src/shammodels/sph/src/Model.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/shammodels/sph/src/Model.cpp b/src/shammodels/sph/src/Model.cpp index 1e7d080283..cf57f9bfb9 100644 --- a/src/shammodels/sph/src/Model.cpp +++ b/src/shammodels/sph/src/Model.cpp @@ -1388,8 +1388,10 @@ void shammodels::sph::Model::init_from_phantom_dump( post_insert_data(sched); - // add sinks + } + // add sinks + if (phdump.blocks.size() > 1) { PhantomDumpBlock &sink_block = phdump.blocks[1]; { std::vector xsink, ysink, zsink; From 26ea3abaf118357a1ddd2564ed2155b6401f4f00 Mon Sep 17 00:00:00 2001 From: Jacksen Narvaez Date: Tue, 23 Jun 2026 10:40:16 -0230 Subject: [PATCH 2/4] Read MHD block (B/rho and psi) from dump files --- .../shammodels/sph/config/MHDConfig.hpp | 5 ++ .../shammodels/sph/io/Phantom2Shamrock.hpp | 11 +++++ src/shammodels/sph/src/Model.cpp | 46 ++++++++++++++++++- .../sph/src/io/Phantom2Shamrock.cpp | 15 ++++++ 4 files changed, 75 insertions(+), 2 deletions(-) diff --git a/src/shammodels/sph/include/shammodels/sph/config/MHDConfig.hpp b/src/shammodels/sph/include/shammodels/sph/config/MHDConfig.hpp index 33daf37ee0..4d3bb38584 100644 --- a/src/shammodels/sph/include/shammodels/sph/config/MHDConfig.hpp +++ b/src/shammodels/sph/include/shammodels/sph/config/MHDConfig.hpp @@ -56,6 +56,11 @@ struct shammodels::sph::MHDConfig { void set(Variant v) { config = v; } + void set_ideal_mhd_constrained_hyper_para( + Tscal sigma_mhd, Tscal alpha_u) { + set(IdealMHD_constrained_hyper_para{sigma_mhd, alpha_u}); + } + inline bool has_B_field() { bool is_B = bool(std::get_if(&config)) || bool(std::get_if(&config)); diff --git a/src/shammodels/sph/include/shammodels/sph/io/Phantom2Shamrock.hpp b/src/shammodels/sph/include/shammodels/sph/io/Phantom2Shamrock.hpp index 3581e1c8e0..4d0a1deb7e 100644 --- a/src/shammodels/sph/include/shammodels/sph/io/Phantom2Shamrock.hpp +++ b/src/shammodels/sph/include/shammodels/sph/io/Phantom2Shamrock.hpp @@ -18,6 +18,7 @@ #include "shammodels/common/EOSConfig.hpp" #include "shammodels/sph/config/AVConfig.hpp" +#include "shammodels/sph/config/MHDConfig.hpp" #include "shammodels/sph/config/BCConfig.hpp" #include "shammodels/sph/io/PhantomDump.hpp" #include "shamunits/UnitSystem.hpp" @@ -54,6 +55,16 @@ namespace shammodels::sph { template AVConfig get_shamrock_avconfig(PhantomDump &phdump); + /** + * @brief Generate an Shamrock MHD configuration from a PhantomDump object. + * + * @param phdump Reference to the PhantomDump object. + * + * @return The MHD configuration corresponding to the given PhantomDump object. + */ + template + MHDConfig get_shamrock_mhdconfig(PhantomDump &phdump); + /** * @brief Get the shamrock units object * \todo load also magfd diff --git a/src/shammodels/sph/src/Model.cpp b/src/shammodels/sph/src/Model.cpp index cf57f9bfb9..31224797bf 100644 --- a/src/shammodels/sph/src/Model.cpp +++ b/src/shammodels/sph/src/Model.cpp @@ -1202,6 +1202,7 @@ auto shammodels::sph::Model::gen_config_from_phantom_dump( conf.eos_config = get_shamrock_eosconfig(phdump, bypass_error); conf.artif_viscosity = get_shamrock_avconfig(phdump); + conf.mhd_config = get_shamrock_mhdconfig(phdump); conf.set_units(get_shamrock_units(phdump)); @@ -1224,9 +1225,12 @@ void shammodels::sph::Model::init_from_phantom_dump( std::vector xyz, vxyz; std::vector h, u, alpha; + std::vector Brhoxyz; + std::vector psich; { std::vector x, y, z, vx, vy, vz; + std::vector Brhox, Brhoy, Brhoz; phdump.blocks[0].fill_vec("x", x); phdump.blocks[0].fill_vec("y", y); @@ -1261,6 +1265,13 @@ void shammodels::sph::Model::init_from_phantom_dump( // expand the box d *= box_tolerance; + // Log the original calculated box dimensions for reference + if (shamcomm::world_rank() == 0) { + logger::info_ln("Model", + "Box domain (from particles): center = ", center, + ", half-size = ", d); + } + resize_simulation_box({center - d, center + d}); } @@ -1273,12 +1284,25 @@ void shammodels::sph::Model::init_from_phantom_dump( phdump.blocks[0].fill_vec("u", u); phdump.blocks[0].fill_vec("alpha", alpha); + // MHD fields + for (auto &block : phdump.blocks) { + block.fill_vec("B/rhox", Brhox); + block.fill_vec("B/rhoy", Brhoy); + block.fill_vec("B/rhoz", Brhoz); + block.fill_vec("psi/ch", psich); + } + for (u32 i = 0; i < x.size(); i++) { xyz.push_back({x[i], y[i], z[i]}); } + for (u32 i = 0; i < vx.size(); i++) { vxyz.push_back({vx[i], vy[i], vz[i]}); } + + for (u32 i = 0; i < Brhox.size(); i++) { + Brhoxyz.push_back({Brhox[i], Brhoy[i], Brhoz[i]}); + } } // Load time infos @@ -1332,8 +1356,8 @@ void shammodels::sph::Model::init_from_phantom_dump( patch_coord.lower, patch_coord.upper); - std::vector ins_xyz, ins_vxyz; - std::vector ins_h, ins_u, ins_alpha; + std::vector ins_xyz, ins_vxyz, ins_Brhoxyz; + std::vector ins_h, ins_u, ins_alpha, ins_psich; for (u64 i : sel_index) { ins_xyz.push_back(xyz[i]); } @@ -1353,6 +1377,16 @@ void shammodels::sph::Model::init_from_phantom_dump( ins_alpha.push_back(alpha[i]); } } + if (Brhoxyz.size() > 0) { + for (u64 i : sel_index) { + ins_Brhoxyz.push_back(Brhoxyz[i]); + } + } + if (psich.size() > 0) { + for (u64 i : sel_index) { + ins_psich.push_back(psich[i]); + } + } PatchDataLayer ptmp(sched.get_layout_ptr_old()); ptmp.resize(sel_index.size()); @@ -1370,6 +1404,14 @@ void shammodels::sph::Model::init_from_phantom_dump( ptmp.override_patch_field("uint", ins_u); } + if (ins_Brhoxyz.size() > 0) { + ptmp.override_patch_field("B/rho", ins_Brhoxyz); + } + + if (ins_psich.size() > 0) { + ptmp.override_patch_field("psi/ch", ins_psich); + } + pdat.insert_elements(ptmp); }); diff --git a/src/shammodels/sph/src/io/Phantom2Shamrock.cpp b/src/shammodels/sph/src/io/Phantom2Shamrock.cpp index 29e73ab6e5..e8c637cf8e 100644 --- a/src/shammodels/sph/src/io/Phantom2Shamrock.cpp +++ b/src/shammodels/sph/src/io/Phantom2Shamrock.cpp @@ -143,6 +143,21 @@ namespace shammodels::sph { /// explicit instanciation for f64_3 template AVConfig get_shamrock_avconfig(PhantomDump &phdump); + template + MHDConfig get_shamrock_mhdconfig(PhantomDump &phdump) { + MHDConfig cfg{}; + + // Defaut values + cfg.set_ideal_mhd_constrained_hyper_para(0.1, 1); + + return cfg; + } + + /// explicit instanciation for f32_3 + template MHDConfig get_shamrock_mhdconfig(PhantomDump &phdump); + /// explicit instanciation for f64_3 + template MHDConfig get_shamrock_mhdconfig(PhantomDump &phdump); + template shamunits::UnitSystem get_shamrock_units(PhantomDump &phdump) { From 34377c1f979d1bd3e08478df80b6e748211bbb09 Mon Sep 17 00:00:00 2001 From: Jacksen Narvaez Date: Tue, 23 Jun 2026 12:20:57 -0230 Subject: [PATCH 3/4] Write MHD block (B/rho and psi) in dump files --- .../sph/include/shammodels/sph/Model.hpp | 2 + src/shammodels/sph/src/Model.cpp | 60 +++++++++++++++++++ 2 files changed, 62 insertions(+) diff --git a/src/shammodels/sph/include/shammodels/sph/Model.hpp b/src/shammodels/sph/include/shammodels/sph/Model.hpp index 2fcfcae5e4..640e6b75a0 100644 --- a/src/shammodels/sph/include/shammodels/sph/Model.hpp +++ b/src/shammodels/sph/include/shammodels/sph/Model.hpp @@ -964,6 +964,8 @@ namespace shammodels::sph { private: void add_pdat_to_phantom_block( PhantomDumpBlock &block, shamrock::patch::PatchDataLayer &pdat); + void add_pdat_to_phantom_block_mhd( + PhantomDumpBlock &block, shamrock::patch::PatchDataLayer &pdat); template inline void warp_disc( diff --git a/src/shammodels/sph/src/Model.cpp b/src/shammodels/sph/src/Model.cpp index 31224797bf..febfe77893 100644 --- a/src/shammodels/sph/src/Model.cpp +++ b/src/shammodels/sph/src/Model.cpp @@ -1520,6 +1520,57 @@ void shammodels::sph::Model::add_pdat_to_phantom_block( block.tot_count = block.blocks_fort_real[xid].vals.size(); } +template class SPHKernel> +void shammodels::sph::Model::add_pdat_to_phantom_block_mhd( + PhantomDumpBlock &block, shamrock::patch::PatchDataLayer &pdat) { + + if (solver.solver_config.has_field_B_on_rho()) { + std::vector Brhoxyz = pdat.fetch_data("B/rho"); + + u64 Brhoxid = block.get_ref_fort_real("B/rhox"); + u64 Brhoyid = block.get_ref_fort_real("B/rhoy"); + u64 Brhozid = block.get_ref_fort_real("B/rhoz"); + + for (auto vec : Brhoxyz) { + block.blocks_fort_real[Brhoxid].vals.push_back(vec.x()); + block.blocks_fort_real[Brhoyid].vals.push_back(vec.y()); + block.blocks_fort_real[Brhozid].vals.push_back(vec.z()); + } + + block.tot_count = block.blocks_fort_real[Brhoxid].vals.size(); + } + + if (solver.solver_config.has_field_psi_on_ch()) { + std::vector psich = pdat.fetch_data("psi/ch"); + u64 psid = block.get_ref_fort_real("psi/ch"); + for (auto ps_ : psich) { + block.blocks_fort_real[psid].vals.push_back(ps_); + } + } + + if (solver.solver_config.has_field_divB()) { + std::vector divB = pdat.fetch_data("divB"); + u64 divBid = block.get_ref_f32("divB"); + for (auto d_ : divB) { + block.blocks_f32[divBid].vals.push_back(d_); + } + } + + if (solver.solver_config.has_field_curlB()) { + std::vector curlB = pdat.fetch_data("curlB"); + + u64 curlBxid = block.get_ref_f32("curlBx"); + u64 curlByid = block.get_ref_f32("curlBy"); + u64 curlBzid = block.get_ref_f32("curlBz"); + + for (auto vec : curlB) { + block.blocks_f32[curlBxid].vals.push_back(vec.x()); + block.blocks_f32[curlByid].vals.push_back(vec.y()); + block.blocks_f32[curlBzid].vals.push_back(vec.z()); + } + } +} + template class SPHKernel> shammodels::sph::PhantomDump shammodels::sph::Model::make_phantom_dump() { StackEntry stack_loc{}; @@ -1621,6 +1672,7 @@ shammodels::sph::PhantomDump shammodels::sph::Model::make_phant write_shamrock_units_in_phantom_dump(solver.solver_config.unit_sys, dump, bypass_error_check); PhantomDumpBlock block_part; + PhantomDumpBlock block_mhd; { NamedStackEntry stack_loc{"gather data"}; @@ -1630,6 +1682,10 @@ shammodels::sph::PhantomDump shammodels::sph::Model::make_phant for (auto &dat : gathered) { add_pdat_to_phantom_block(block_part, shambase::get_check_ref(dat)); } + + for (auto &dat : gathered) { + add_pdat_to_phantom_block_mhd(block_mhd, shambase::get_check_ref(dat)); + } } dump.blocks.push_back(std::move(block_part)); @@ -1665,6 +1721,10 @@ shammodels::sph::PhantomDump shammodels::sph::Model::make_phant dump.blocks.push_back(std::move(sink_block)); } + PhantomDumpBlock block_3rd; + dump.blocks.push_back(std::move(block_3rd)); + dump.blocks.push_back(std::move(block_mhd)); + return dump; } From e0feea1144f0c8509d7928241023e3983bbea520 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 24 Jun 2026 09:21:30 +0000 Subject: [PATCH 4/4] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../sph/include/shammodels/sph/config/MHDConfig.hpp | 3 +-- .../sph/include/shammodels/sph/io/Phantom2Shamrock.hpp | 2 +- src/shammodels/sph/src/Model.cpp | 10 ++++------ 3 files changed, 6 insertions(+), 9 deletions(-) diff --git a/src/shammodels/sph/include/shammodels/sph/config/MHDConfig.hpp b/src/shammodels/sph/include/shammodels/sph/config/MHDConfig.hpp index 4d3bb38584..76e90cfa3c 100644 --- a/src/shammodels/sph/include/shammodels/sph/config/MHDConfig.hpp +++ b/src/shammodels/sph/include/shammodels/sph/config/MHDConfig.hpp @@ -56,8 +56,7 @@ struct shammodels::sph::MHDConfig { void set(Variant v) { config = v; } - void set_ideal_mhd_constrained_hyper_para( - Tscal sigma_mhd, Tscal alpha_u) { + void set_ideal_mhd_constrained_hyper_para(Tscal sigma_mhd, Tscal alpha_u) { set(IdealMHD_constrained_hyper_para{sigma_mhd, alpha_u}); } diff --git a/src/shammodels/sph/include/shammodels/sph/io/Phantom2Shamrock.hpp b/src/shammodels/sph/include/shammodels/sph/io/Phantom2Shamrock.hpp index 4d0a1deb7e..c533374d22 100644 --- a/src/shammodels/sph/include/shammodels/sph/io/Phantom2Shamrock.hpp +++ b/src/shammodels/sph/include/shammodels/sph/io/Phantom2Shamrock.hpp @@ -18,8 +18,8 @@ #include "shammodels/common/EOSConfig.hpp" #include "shammodels/sph/config/AVConfig.hpp" -#include "shammodels/sph/config/MHDConfig.hpp" #include "shammodels/sph/config/BCConfig.hpp" +#include "shammodels/sph/config/MHDConfig.hpp" #include "shammodels/sph/io/PhantomDump.hpp" #include "shamunits/UnitSystem.hpp" #include diff --git a/src/shammodels/sph/src/Model.cpp b/src/shammodels/sph/src/Model.cpp index febfe77893..d4cb0e0e0b 100644 --- a/src/shammodels/sph/src/Model.cpp +++ b/src/shammodels/sph/src/Model.cpp @@ -1266,11 +1266,10 @@ void shammodels::sph::Model::init_from_phantom_dump( d *= box_tolerance; // Log the original calculated box dimensions for reference - if (shamcomm::world_rank() == 0) { - logger::info_ln("Model", - "Box domain (from particles): center = ", center, - ", half-size = ", d); - } + if (shamcomm::world_rank() == 0) { + logger::info_ln( + "Model", "Box domain (from particles): center = ", center, ", half-size = ", d); + } resize_simulation_box({center - d, center + d}); } @@ -1429,7 +1428,6 @@ void shammodels::sph::Model::init_from_phantom_dump( .update_load_balancing(); post_insert_data(sched); - } // add sinks