[Stagging PR] Phantom emulation mode

doc/sphinx/examples/sph/run_phantom_run_emulation.py

@@ -38,4 +38,49 @@
 shamrock.utils.url.download_file(input_file_url, input_file_path)
 shamrock.utils.url.download_file(dump_file_url, dump_file_path)
 
-shamrock.utils.phantom.run_phantom_simulation(dump_folder, "disc")
+
+def plot_that_rippa(ctx, model, dump_number):
+    pixel_x = 1920
+    pixel_y = 1080
+    radius = 5
+    center = (0.0, 0.0, 0.0)
+
+    aspect = pixel_x / pixel_y
+    pic_range = [-radius * aspect, radius * aspect, -radius, radius]
+    delta_x = (radius * 2 * aspect, 0.0, 0.0)
+    delta_y = (0.0, radius * 2, 0.0)
+
+    arr_rho = model.render_cartesian_column_integ(
+        "rho",
+        "f64",
+        center=(0.0, 0.0, 0.0),
+        delta_x=delta_x,
+        delta_y=delta_y,
+        nx=pixel_x,
+        ny=pixel_y,
+    )
+
+    import copy
+
+    import matplotlib
+    import matplotlib.pyplot as plt
+
+    my_cmap = copy.copy(matplotlib.colormaps.get_cmap("gist_heat"))  # copy the default cmap
+    my_cmap.set_bad(color="black")
+
+    plt.figure(figsize=(16 / 2, 9 / 2))
+    res = plt.imshow(arr_rho, cmap=my_cmap, origin="lower", extent=pic_range, norm="log", vmin=1e-9)
+
+    cbar = plt.colorbar(res, extend="both")
+    cbar.set_label(r"$\int \rho \, \mathrm{d} z$ [code unit]")
+    # or r"$\rho$ [code unit]" for slices
+
+    plt.title("t = {:0.3f} [code unit]".format(model.get_time()))
+    plt.xlabel("x")
+    plt.ylabel("z")
+    plt.show()
+
+
+ctx, model, in_params = shamrock.utils.phantom.run_phantom_simulation(
+    dump_folder, "disc", callback=plot_that_rippa
+)

src/pylib/shamrock/utils/phantom/__init__.py

@@ -138,7 +138,7 @@ def load_simulation(simulation_path, dump_file_name=None, in_file_name=None, do_
     return ctx, model, in_params
 
 
-def run_phantom_simulation(simulation_folder, sim_name):
+def run_phantom_simulation(simulation_folder, sim_name, callback=None):
     """
     Run a Phantom simulation in Shamrock.
     """
@@ -191,6 +191,12 @@ def do_dump(dump_number):
         ):
             os.remove(os.path.join(simulation_folder, f"{sim_name}_00000.tmp"))
 
+        if callback is not None:
+            callback(ctx, model, dump_number)
+
+    model.set_next_dt(0)
+    model.timestep()
+
     do_dump(dump_number)
 
     dump_number += 1
@@ -204,7 +210,7 @@ def do_dump(dump_number):
             next_time = tmax
             last_step = True
 
-        last_step = model.evolve_until(next_time)
+        model.evolve_until(next_time)
 
         do_dump(dump_number)
         dump_number += 1

src/shambackends/include/shambackends/DeviceBufferRef.hpp

@@ -0,0 +1,26 @@
+// -------------------------------------------------------//
+//
+// SHAMROCK code for hydrodynamics
+// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
+// 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 DeviceBufferRef.hpp
+ * @author Timothée David--Cléris (tim.shamrock@proton.me)
+ * @brief
+ *
+ */
+
+#include "shambackends/DeviceBuffer.hpp"
+
+namespace sham {
+
+    template<class T, USMKindTarget target = device>
+    using DeviceBufferRef = std::reference_wrapper<DeviceBuffer<T, target>>;
+
+}

src/shammodels/gsph/include/shammodels/gsph/modules/SolverStorage.hpp

@@ -70,6 +70,9 @@ namespace shammodels::gsph {
 
         using RTree = shamtree::CompressedLeafBVH<Tmorton, Tvec, 3>;
 
+        std::shared_ptr<shamrock::solvergraph::ExchangeGhostLayer> exchange_gz_node;
+        std::shared_ptr<shamrock::solvergraph::ExchangeGhostLayer> exchange_gz_positions;
+
         /// Particle counts per patch
         std::shared_ptr<shamrock::solvergraph::Indexes<u32>> part_counts;
         std::shared_ptr<shamrock::solvergraph::Indexes<u32>> part_counts_with_ghost;

src/shammodels/gsph/src/Solver.cpp

@@ -93,6 +93,10 @@ void shammodels::gsph::Solver<Tvec, Kern>::init_solver_graph() {
     storage.soundspeed
         = std::make_shared<shamrock::solvergraph::Field<Tscal>>(1, "soundspeed", "c_s");
 
+    storage.exchange_gz_node
+        = std::make_shared<shamrock::solvergraph::ExchangeGhostLayer>(storage.ghost_layout);
+    storage.exchange_gz_positions
+        = std::make_shared<shamrock::solvergraph::ExchangeGhostLayer>(storage.xyzh_ghost_layout);
     // Initialize gradient fields for MUSCL reconstruction
     // These are only used when reconstruct_config.is_muscl() == true
     storage.grad_density
@@ -190,6 +194,9 @@ template<class Tvec, template<class> class Kern>
 void shammodels::gsph::Solver<Tvec, Kern>::merge_position_ghost() {
     StackEntry stack_loc{};
 
+    std::shared_ptr<shamrock::solvergraph::ExchangeGhostLayer> exchange_gz_node
+        = std::make_shared<shamrock::solvergraph::ExchangeGhostLayer>(storage.xyzh_ghost_layout);
+
     storage.merged_xyzh.set(
         storage.ghost_handler.get().build_comm_merge_positions(storage.ghost_patch_cache.get()));
 
@@ -618,7 +625,7 @@ void shammodels::gsph::Solver<Tvec, Kern>::init_ghost_layout() {
     storage.ghost_layout = std::make_shared<shamrock::patch::PatchDataLayerLayout>();
 
     shamrock::patch::PatchDataLayerLayout &ghost_layout
-        = shambase::get_check_ref(storage.ghost_layout.get());
+        = shambase::get_check_ref(storage.ghost_layout);
 
     solver_config.set_ghost_layout(ghost_layout);
 }
@@ -746,6 +753,9 @@ void shammodels::gsph::Solver<Tvec, Kern>::communicate_merge_ghosts_fields() {
             }
         });
 
+    std::shared_ptr<shamrock::solvergraph::ExchangeGhostLayer> exchange_gz_node
+        = std::make_shared<shamrock::solvergraph::ExchangeGhostLayer>(storage.ghost_layout);
+
     // Communicate ghost data across MPI ranks
     shambase::DistributedDataShared<PatchDataLayer> interf_pdat
         = ghost_handle.communicate_pdat(ghost_layout_ptr, std::move(pdat_interf));

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

@@ -151,6 +151,20 @@ namespace shammodels::sph {
         f64 utime  = phdump.read_header_float<f64>("utime");
         f64 umagfd = phdump.read_header_float<f64>("umagfd");
 
+        if (udist == 1.0 || umass == 1.0 || utime == 1.0 || umagfd == 3.54491) {
+            logger::warn_ln("SPH", "phantom dump units are not set, defaulting to SI");
+            logger::warn_ln("SPH", "udist =", udist);
+            logger::warn_ln("SPH", "umass =", umass);
+            logger::warn_ln("SPH", "utime =", utime);
+            logger::warn_ln("SPH", "umagfd =", umagfd);
+
+            return shamunits::UnitSystem<Tscal>();
+        }
+
+        // convert from CGS to SI
+        udist /= 100.0;
+        umass /= 1000.0;
+
         return shamunits::UnitSystem<Tscal>(
             utime, udist, umass
             // unit_current = 1 ,
@@ -169,9 +183,10 @@ namespace shammodels::sph {
             dump.table_header_f64.add("umass", units->kg_inv);
             dump.table_header_f64.add("utime", units->s_inv);
 
-            f64 umass = units->template to<shamunits::units::kg>();
+            // Back to freakin CGS (worst units system ever, well no ... there is imperial)
+            f64 umass = units->template to<shamunits::units::kg>() / 1000.0;
             f64 utime = units->template to<shamunits::units::s>();
-            f64 udist = units->template to<shamunits::units::m>();
+            f64 udist = units->template to<shamunits::units::m>() / 100.0;
 
             shamunits::Constants<Tscal> ctes{*units};
             f64 ccst    = ctes.c();

src/shammodels/sph/src/modules/SPHSetup.cpp

@@ -444,6 +444,9 @@ void shammodels::sph::modules::SPHSetup<Tvec, SPHKernel>::apply_setup_new(
                 has_been_limited = false;
                 using namespace shamrock::patch;
 
+                auto dev_sched = shamsys::instance::get_compute_scheduler_ptr();
+                sham::DeviceBuffer<u32> mask_get_ids_where(0, dev_sched);
+
                 // inject in local domains first
                 PatchCoordTransform<Tvec> ptransf = sched.get_sim_box().get_patch_transform<Tvec>();
                 sched.for_each_local_patchdata([&](const Patch p, PatchDataLayer &pdat) {