[Solvergraph][Ramses][PIC] Draft PR for PIC solvergraph node with empty kernel on grid

src/shammodels/ramses/CMakeLists.txt

@@ -44,6 +44,7 @@ set(Sources
     src/modules/TransformGhostLayer.cpp
     src/modules/FuseGhostLayer.cpp
     src/modules/render/GridRender.cpp
+    src/modules/ParticleInCell.cpp
 )
 
 if(SHAMROCK_USE_SHARED_LIB)

src/shammodels/ramses/include/shammodels/ramses/Solver.hpp

@@ -71,6 +71,11 @@ namespace shammodels::basegodunov {
                 context.pdata_layout_add_field<Tscal>(
                     "rho_gas_pscal", (npscal_gas * AMRBlock::block_size));
             }
+
+            if (solver_config.is_pic_enabled()) {
+                context.pdata_layout_add_field<Tscal>("mass_particles", AMRBlock::block_size);
+                context.pdata_layout_add_field<Tscal>("rho_pic", AMRBlock::block_size);
+            }
         }
 
         Solver(ShamrockCtx &context) : context(context) {}

src/shammodels/ramses/include/shammodels/ramses/SolverConfig.hpp

@@ -153,6 +153,10 @@ namespace shammodels::basegodunov {
     template<class Tvec, class TgridVec>
     struct SolverConfig;
 
+    struct PICConfig {
+        bool enabled = false;
+    };
+
 }; // namespace shammodels::basegodunov
 
 template<class Tvec>
@@ -258,6 +262,15 @@ struct shammodels::basegodunov::SolverConfig {
     // Units Config (END)
     //////////////////////////////////////////////////////////////////////////////////////////////
 
+    //////////////////////////////////////////////////////////////////////////////////////////////
+    // Particle-In-Cell config
+    //////////////////////////////////////////////////////////////////////////////////////////////
+    PICConfig pic_config{};
+    inline bool is_pic_enabled() { return pic_config.enabled; }
+    //////////////////////////////////////////////////////////////////////////////////////////////
+    // Particle-In-Cell config (END)
+    //////////////////////////////////////////////////////////////////////////////////////////////
+
     //////////////////////////////////////////////////////////////////////////////////////////////
     // Solver status variables
     //////////////////////////////////////////////////////////////////////////////////////////////

src/shammodels/ramses/include/shammodels/ramses/modules/ParticleInCell.hpp

@@ -0,0 +1,61 @@
+// -------------------------------------------------------//
+//
+// SHAMROCK code for hydrodynamics
+// Copyright (c) 2021-2025 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 ParticleInCell.hpp
+ * @author Anass Serhani (anass.serhani@cnrs.fr)
+ * @author Timothée David--Cléris (tim.shamrock@proton.me) --no git blame--
+ * @brief
+ */
+
+#include "shambackends/vec.hpp"
+#include "shamrock/solvergraph/IFieldSpan.hpp"
+#include "shamrock/solvergraph/INode.hpp"
+#include "shamrock/solvergraph/Indexes.hpp"
+
+namespace shammodels::basegodunov::modules {
+    template<class Tvec>
+    class NodePIC : public shamrock::solvergraph::INode {
+        using Tscal = shambase::VecComponent<Tvec>;
+        u32 block_size;
+
+        public:
+        NodePIC(u32 block_size) : block_size(block_size) {}
+
+        struct Edges {
+            const shamrock::solvergraph::Indexes<u32> &sizes;
+            const shamrock::solvergraph::IFieldSpan<Tscal> &spans_mass_particles;
+            shamrock::solvergraph::IFieldSpan<Tscal> &spans_rho_pic;
+        };
+
+        inline void set_edges(
+            std::shared_ptr<shamrock::solvergraph::Indexes<u32>> sizes,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> spans_mass_particles,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> spans_rho_pic) {
+            __internal_set_ro_edges({sizes, spans_mass_particles});
+            __internal_set_rw_edges({spans_rho_pic});
+        }
+
+        inline Edges get_edges() {
+            return Edges{
+                get_ro_edge<shamrock::solvergraph::Indexes<u32>>(0),
+                get_ro_edge<shamrock::solvergraph::IFieldSpan<Tscal>>(1),
+                get_rw_edge<shamrock::solvergraph::IFieldSpan<Tscal>>(0),
+            };
+        }
+
+        void _impl_evaluate_internal();
+
+        inline virtual std::string _impl_get_label() const { return "ParticleInCell"; };
+
+        virtual std::string _impl_get_tex() const;
+    };
+} // namespace shammodels::basegodunov::modules

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

@@ -82,6 +82,9 @@ namespace shammodels::basegodunov {
         std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> refs_rho_dust;
         std::shared_ptr<shamrock::solvergraph::FieldRefs<Tvec>> refs_rhov_dust;
 
+        std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> refs_mass_particles;
+        std::shared_ptr<shamrock::solvergraph::Field<Tscal>> rho_pic;
+
         std::shared_ptr<shamrock::solvergraph::Field<Tvec>> vel;
         std::shared_ptr<shamrock::solvergraph::Field<Tscal>> press;
         std::shared_ptr<shamrock::solvergraph::Field<Tvec>> vel_dust;

src/shammodels/ramses/src/Solver.cpp

@@ -38,6 +38,7 @@
 #include "shammodels/ramses/modules/FuseGhostLayer.hpp"
 #include "shammodels/ramses/modules/InterpolateToFace.hpp"
 #include "shammodels/ramses/modules/NodeComputeFlux.hpp"
+#include "shammodels/ramses/modules/ParticleInCell.hpp"
 #include "shammodels/ramses/modules/SlopeLimitedGradient.hpp"
 #include "shammodels/ramses/modules/TimeIntegrator.hpp"
 #include "shammodels/ramses/modules/TransformGhostLayer.hpp"
@@ -301,6 +302,11 @@ void shammodels::basegodunov::Solver<Tvec, TgridVec>::init_solver_graph() {
             u32 npscal_gas = solver_config.npscal_gas_config.npscal_gas;
             ghost_layout.add_field<Tscal>("rho_gas_pscal", npscal_gas * AMRBlock::block_size);
         }
+
+        if (solver_config.is_pic_enabled()) {
+            ghost_layout.add_field<Tscal>("mass_particles", AMRBlock::block_size);
+            ghost_layout.add_field<Tscal>("rho_pic", AMRBlock::block_size);
+        }
     }
 
     ////////////////////////////////////////////////////////////////////////////////
@@ -654,6 +660,14 @@ void shammodels::basegodunov::Solver<Tvec, TgridVec>::init_solver_graph() {
                 "flux_rhov_dust_face_zp", "flux_rhov_dust_face_zp", ndust);
     }
 
+    // will be filled by NodePIC
+    if (solver_config.is_pic_enabled()) {
+        storage.refs_mass_particles = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>(
+            "mass_particles", "m_{particles}");
+        storage.rho_pic = std::make_shared<shamrock::solvergraph::Field<Tscal>>(
+            AMRBlock::block_size, "rho_pic", "\\rho_{pic}");
+    }
+
     ////////////////////////////////////////////////////////////////////////////////
     /// Nodes
     ////////////////////////////////////////////////////////////////////////////////
@@ -847,6 +861,17 @@ void shammodels::basegodunov::Solver<Tvec, TgridVec>::init_solver_graph() {
             solver_sequence.push_back(
                 std::make_shared<decltype(attach_rhov_dust)>(std::move(attach_rhov_dust)));
         }
+
+        if (solver_config.is_pic_enabled()) { // attach spans to PIC field with ghosts (temporary)
+            shamrock::solvergraph::GetFieldRefFromLayer<Tscal> attach_mass_particles
+                = shamrock::solvergraph::GetFieldRefFromLayer<Tscal>(
+                    storage.ghost_layout, "mass_particles");
+            attach_mass_particles.set_edges(
+                storage.merged_patchdata_ghost, storage.refs_mass_particles);
+            solver_sequence.push_back(
+                std::make_shared<decltype(attach_mass_particles)>(
+                    std::move(attach_mass_particles)));
+        }
     }
 
     { // build trees
@@ -918,6 +943,22 @@ void shammodels::basegodunov::Solver<Tvec, TgridVec>::init_solver_graph() {
         solver_sequence.push_back(std::make_shared<decltype(node2)>(std::move(node2)));
     }
 
+    // Build PIC (Particle in Cell) node
+    if (solver_config.is_pic_enabled()) {
+        std::vector<std::shared_ptr<shamrock::solvergraph::INode>> pic_sequence;
+
+        {
+            modules::NodePIC<Tvec> node{AMRBlock::block_size};
+            node.set_edges(
+                storage.block_counts_with_ghost, storage.refs_mass_particles, storage.rho_pic);
+
+            pic_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
+        }
+
+        shamrock::solvergraph::OperationSequence seq("Particle in Cell", std::move(pic_sequence));
+        solver_sequence.push_back(std::make_shared<decltype(seq)>(std::move(seq)));
+    }
+
     { // Build ConsToPrim node
         std::vector<std::shared_ptr<shamrock::solvergraph::INode>> const_to_prim_sequence;
 
@@ -1353,10 +1394,19 @@ void shammodels::basegodunov::Solver<Tvec, TgridVec>::init_solver_graph() {
     shamrock::solvergraph::OperationSequence seq("Solver", std::move(solver_sequence));
     storage.solver_sequence = std::make_shared<decltype(seq)>(std::move(seq));
 
-    if (false) {
-        logger::raw_ln(" -- tex:\n" + shambase::get_check_ref(storage.solver_sequence).get_tex());
-        logger::raw_ln(
-            " -- dot:\n" + shambase::get_check_ref(storage.solver_sequence).get_dot_graph());
+    if (true) {
+        // logger::raw_ln(" -- tex:\n" +
+        // shambase::get_check_ref(storage.solver_sequence).get_tex()); logger::raw_ln(
+        //     " -- dot:\n" + shambase::get_check_ref(storage.solver_sequence).get_dot_graph());
+
+        if (shamcomm::world_rank() == 0) {
+            std::ofstream graph_dot_file;
+            graph_dot_file.open("./solvergraph.dot");
+            graph_dot_file << "digraph G {\n";
+            graph_dot_file << shambase::get_check_ref(storage.solver_sequence).get_dot_graph();
+            graph_dot_file << "}\n";
+            graph_dot_file.close();
+        }
     }
 }
 

src/shammodels/ramses/src/modules/ParticleInCell.cpp

@@ -0,0 +1,95 @@
+// -------------------------------------------------------//
+//
+// SHAMROCK code for hydrodynamics
+// Copyright (c) 2021-2025 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
+//
+// -------------------------------------------------------//
+
+/**
+ * @file ParticleInCell.cpp
+ * @author Anass Serhani (anass.serhani@cnrs.fr)
+ * @author Timothée David--Cléris (tim.shamrock@proton.me) --no git blame--
+ * @brief
+ */
+
+#include "shammodels/ramses/modules/ParticleInCell.hpp"
+#include "shambackends/kernel_call_distrib.hpp"
+#include "shammath/riemann.hpp"
+#include "shamrock/patch/PatchDataField.hpp"
+#include "shamsys/NodeInstance.hpp"
+
+namespace {
+
+    template<class Tvec>
+    struct KernelPIC {
+        using Tscal = shambase::VecComponent<Tvec>;
+
+        inline static void kernel(
+            const shambase::DistributedData<shamrock::PatchDataFieldSpanPointer<Tscal>>
+                &spans_mass_particles,
+            shambase::DistributedData<shamrock::PatchDataFieldSpanPointer<Tscal>> &spans_rho_pic,
+            const shambase::DistributedData<u32> &sizes,
+            u32 block_size) {
+
+            shambase::DistributedData<u32> cell_counts
+                = sizes.map<u32>([&](u64 id, u32 block_count) {
+                      u32 cell_count = block_count * block_size;
+                      return cell_count;
+                  });
+
+            sham::distributed_data_kernel_call(
+                shamsys::instance::get_compute_scheduler_ptr(),
+                sham::DDMultiRef{spans_mass_particles},
+                sham::DDMultiRef{spans_rho_pic},
+                cell_counts,
+                [](u32 i, const Tscal *__restrict mass_particles, Tscal *__restrict rho_pic) {
+                    /*
+
+                    PIC KERNELS IN PROGRESS
+
+                    */
+                });
+        }
+    };
+
+} // namespace
+
+namespace shammodels::basegodunov::modules {
+
+    template<class Tvec>
+    void NodePIC<Tvec>::_impl_evaluate_internal() {
+        auto edges = get_edges();
+
+        edges.spans_mass_particles.check_sizes(edges.sizes.indexes);
+        edges.spans_rho_pic.ensure_sizes(edges.sizes.indexes);
+
+        KernelPIC<Tvec>::kernel(
+            edges.spans_mass_particles.get_spans(),
+            edges.spans_rho_pic.get_spans(),
+            edges.sizes.indexes,
+            block_size);
+    }
+
+    template<class Tvec>
+    std::string NodePIC<Tvec>::_impl_get_tex() const {
+        auto block_count    = get_ro_edge_base(0).get_tex_symbol();
+        auto mass_particles = get_ro_edge_base(1).get_tex_symbol();
+        auto rho_pic        = get_rw_edge_base(0).get_tex_symbol();
+
+        std::string tex = R"tex(
+            // TODO: Add TeX description here
+        )tex";
+
+        shambase::replace_all(tex, "{rho_pic}", rho_pic);
+        shambase::replace_all(tex, "{mass_particles}", mass_particles);
+        shambase::replace_all(tex, "{block_count}", block_count);
+        shambase::replace_all(tex, "{block_size}", shambase::format("{}", block_size));
+
+        return tex;
+    }
+
+} // namespace shammodels::basegodunov::modules
+
+template class shammodels::basegodunov::modules::NodePIC<f64_3>;

src/shammodels/ramses/src/pyRamsesModel.cpp

@@ -205,8 +205,13 @@ namespace shammodels::basegodunov {
                 [](TConfig &self) {
                     self.gravity_config.gravity_mode = BICGSTAB;
                 })
-            .def("set_npscal_gas", [](TConfig &self, u32 npscal_gas) {
-                self.npscal_gas_config.npscal_gas = npscal_gas;
+            .def(
+                "set_npscal_gas",
+                [](TConfig &self, u32 npscal_gas) {
+                    self.npscal_gas_config.npscal_gas = npscal_gas;
+                })
+            .def("set_pic", [](TConfig &self) {
+                self.pic_config.enabled = true;
             });
 
         std::string sod_tube_analysis_name = name_model + "_AnalysisSodTube";