[SPH] Add MHD block to Phantom I/O

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<class Tscal>
         inline void warp_disc(

src/shammodels/sph/include/shammodels/sph/config/MHDConfig.hpp

@@ -56,6 +56,10 @@ 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<IdealMHD_constrained_hyper_para>(&config))
                     || bool(std::get_if<NonIdealMHD>(&config));

src/shammodels/sph/include/shammodels/sph/io/Phantom2Shamrock.hpp

@@ -19,6 +19,7 @@
 #include "shammodels/common/EOSConfig.hpp"
 #include "shammodels/sph/config/AVConfig.hpp"
 #include "shammodels/sph/config/BCConfig.hpp"
+#include "shammodels/sph/config/MHDConfig.hpp"
 #include "shammodels/sph/io/PhantomDump.hpp"
 #include "shamunits/UnitSystem.hpp"
 #include <cstdlib>
@@ -54,6 +55,16 @@ namespace shammodels::sph {
     template<class Tvec>
     AVConfig<Tvec> 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<class Tvec>
+    MHDConfig<Tvec> get_shamrock_mhdconfig(PhantomDump &phdump);
+
     /**
      * @brief Get the shamrock units object
      * \todo load also magfd

src/shammodels/sph/src/Model.cpp

@@ -1202,6 +1202,7 @@ auto shammodels::sph::Model<Tvec, SPHKernel>::gen_config_from_phantom_dump(
 
     conf.eos_config      = get_shamrock_eosconfig<Tvec>(phdump, bypass_error);
     conf.artif_viscosity = get_shamrock_avconfig<Tvec>(phdump);
+    conf.mhd_config      = get_shamrock_mhdconfig<Tvec>(phdump);
 
     conf.set_units(get_shamrock_units<Tscal>(phdump));
 
@@ -1224,9 +1225,12 @@ void shammodels::sph::Model<Tvec, SPHKernel>::init_from_phantom_dump(
 
     std::vector<Tvec> xyz, vxyz;
     std::vector<Tscal> h, u, alpha;
+    std::vector<Tvec> Brhoxyz;
+    std::vector<Tscal> psich;
 
     {
         std::vector<Tscal> x, y, z, vx, vy, vz;
+        std::vector<Tscal> Brhox, Brhoy, Brhoz;
 
         phdump.blocks[0].fill_vec("x", x);
         phdump.blocks[0].fill_vec("y", y);
@@ -1261,6 +1265,12 @@ void shammodels::sph::Model<Tvec, SPHKernel>::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 +1283,25 @@ void shammodels::sph::Model<Tvec, SPHKernel>::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 +1355,8 @@ void shammodels::sph::Model<Tvec, SPHKernel>::init_from_phantom_dump(
                 patch_coord.lower,
                 patch_coord.upper);
 
-            std::vector<Tvec> ins_xyz, ins_vxyz;
-            std::vector<Tscal> ins_h, ins_u, ins_alpha;
+            std::vector<Tvec> ins_xyz, ins_vxyz, ins_Brhoxyz;
+            std::vector<Tscal> ins_h, ins_u, ins_alpha, ins_psich;
             for (u64 i : sel_index) {
                 ins_xyz.push_back(xyz[i]);
             }
@@ -1353,6 +1376,16 @@ void shammodels::sph::Model<Tvec, SPHKernel>::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 +1403,14 @@ void shammodels::sph::Model<Tvec, SPHKernel>::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);
         });
 
@@ -1387,9 +1428,10 @@ void shammodels::sph::Model<Tvec, SPHKernel>::init_from_phantom_dump(
             .update_load_balancing();
 
         post_insert_data<Tvec>(sched);
+    }
+    // add sinks
 
-        // add sinks
-
+    if (phdump.blocks.size() > 1) {
         PhantomDumpBlock &sink_block = phdump.blocks[1];
         {
             std::vector<Tscal> xsink, ysink, zsink;
@@ -1476,6 +1518,57 @@ void shammodels::sph::Model<Tvec, SPHKernel>::add_pdat_to_phantom_block(
     block.tot_count = block.blocks_fort_real[xid].vals.size();
 }
 
+template<class Tvec, template<class> class SPHKernel>
+void shammodels::sph::Model<Tvec, SPHKernel>::add_pdat_to_phantom_block_mhd(
+    PhantomDumpBlock &block, shamrock::patch::PatchDataLayer &pdat) {
+
+    if (solver.solver_config.has_field_B_on_rho()) {
+        std::vector<Tvec> Brhoxyz = pdat.fetch_data<Tvec>("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<Tscal> psich = pdat.fetch_data<Tscal>("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<Tscal> divB = pdat.fetch_data<Tscal>("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<Tvec> curlB = pdat.fetch_data<Tvec>("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 Tvec, template<class> class SPHKernel>
 shammodels::sph::PhantomDump shammodels::sph::Model<Tvec, SPHKernel>::make_phantom_dump() {
     StackEntry stack_loc{};
@@ -1577,6 +1670,7 @@ shammodels::sph::PhantomDump shammodels::sph::Model<Tvec, SPHKernel>::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"};
@@ -1586,6 +1680,10 @@ shammodels::sph::PhantomDump shammodels::sph::Model<Tvec, SPHKernel>::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));
@@ -1621,6 +1719,10 @@ shammodels::sph::PhantomDump shammodels::sph::Model<Tvec, SPHKernel>::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;
 }
 

src/shammodels/sph/src/io/Phantom2Shamrock.cpp

@@ -143,6 +143,21 @@ namespace shammodels::sph {
     /// explicit instanciation for f64_3
     template AVConfig<f64_3> get_shamrock_avconfig<f64_3>(PhantomDump &phdump);
 
+    template<class Tvec>
+    MHDConfig<Tvec> get_shamrock_mhdconfig(PhantomDump &phdump) {
+        MHDConfig<Tvec> cfg{};
+
+        // Defaut values
+        cfg.set_ideal_mhd_constrained_hyper_para(0.1, 1);
+
+        return cfg;
+    }
+
+    /// explicit instanciation for f32_3
+    template MHDConfig<f32_3> get_shamrock_mhdconfig<f32_3>(PhantomDump &phdump);
+    /// explicit instanciation for f64_3
+    template MHDConfig<f64_3> get_shamrock_mhdconfig<f64_3>(PhantomDump &phdump);
+
     template<class Tscal>
     shamunits::UnitSystem<Tscal> get_shamrock_units(PhantomDump &phdump) {