[SPH] Hard wall boundary condition

src/shammodels/sph/CMakeLists.txt

@@ -45,6 +45,7 @@ set(Sources
     src/modules/SPHSetup.cpp
     src/modules/KillParticles.cpp
     src/modules/GetParticlesOutsideSphere.cpp
+    src/modules/GetParticlesInWall.cpp
     src/modules/IterateSmoothingLengthDensity.cpp
     src/modules/IterateSmoothingLengthDensityNeighLim.cpp
     src/modules/LoopSmoothingLengthIter.cpp

src/shammodels/sph/include/shammodels/sph/SolverConfig.hpp

@@ -105,6 +105,26 @@ namespace shammodels::sph {
         }
     };
 
+    template<class Tvec>
+    struct ParticleDisableConfig {
+        using Tscal = shambase::VecComponent<Tvec>;
+
+        struct Wall {
+            Tvec pos;
+            Tscal length;
+            Tscal width;
+            Tscal thickness;
+        };
+
+        using disable_t = std::variant<Wall>;
+
+        std::vector<disable_t> disable_list;
+
+        inline void add_disable_wall(const Tvec &pos, Tscal length, Tscal width, Tscal thickness) {
+            disable_list.push_back(Wall{pos, length, width, thickness});
+        }
+    };
+
     template<class Tscal>
     struct DustConfig {
 
@@ -340,6 +360,7 @@ struct shammodels::sph::SolverConfig {
     //////////////////////////////////////////////////////////////////////////////////////////////
 
     ParticleKillingConfig<Tvec> particle_killing;
+    ParticleDisableConfig<Tvec> particle_disable;
 
     //////////////////////////////////////////////////////////////////////////////////////////////
     // Particle killing config (END)
@@ -950,6 +971,41 @@ namespace shammodels::sph {
         }
     }
 
+    // JSON serialisation for ParticleDisableConfig
+    template<class Tvec>
+    inline void to_json(nlohmann::json &j, const ParticleDisableConfig<Tvec> &p) {
+        j = nlohmann::json::array();
+        for (const auto &disable : p.disable_list) {
+            if (std::holds_alternative<typename ParticleDisableConfig<Tvec>::Wall>(disable)) {
+                const auto &wall = std::get<typename ParticleDisableConfig<Tvec>::Wall>(disable);
+                j.push_back(
+                    {{"type", "wall"},
+                     {"pos", wall.pos},
+                     {"length", wall.length},
+                     {"width", wall.width},
+                     {"thickness", wall.thickness}});
+            }
+            // If more types are added to disable_t, handle them here
+        }
+    }
+
+    template<class Tvec>
+    inline void from_json(const nlohmann::json &j, ParticleDisableConfig<Tvec> &p) {
+        p.disable_list.clear();
+        for (const auto &item : j) {
+            std::string type = item.at("type").get<std::string>();
+            if (type == "wall") {
+                typename ParticleDisableConfig<Tvec>::Wall wall;
+                item.at("pos").get_to(wall.pos);
+                item.at("length").get_to(wall.length);
+                item.at("width").get_to(wall.width);
+                item.at("thickness").get_to(wall.thickness);
+                p.disable_list.push_back(wall);
+            }
+            // If more types are added to disable_t, handle them here
+        }
+    }
+
     // JSON serialization for SmoothingLengthConfig
     inline void to_json(nlohmann::json &j, const SmoothingLengthConfig &p) {
         if (const SmoothingLengthConfig::DensityBased *conf

src/shammodels/sph/include/shammodels/sph/modules/GetParticlesInWall.hpp

@@ -0,0 +1,70 @@
+// -------------------------------------------------------//
+//
+// 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 GetParticlesInWall.hpp
+ * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr)
+ * @brief Declares the GetParticlesInWall module for disabling particles update.
+ *
+ */
+
+#include "shamrock/solvergraph/DistributedBuffers.hpp"
+#include "shamrock/solvergraph/IFieldRefs.hpp"
+#include "shamrock/solvergraph/INode.hpp"
+#include "shamrock/solvergraph/Indexes.hpp"
+#include "shamrock/solvergraph/PatchDataLayerRefs.hpp"
+
+namespace shammodels::sph::modules {
+
+    template<typename Tvec>
+    class GetParticlesInWall : public shamrock::solvergraph::INode {
+
+        using Tscal = shambase::VecComponent<Tvec>;
+
+        Tvec wall_pos;
+        Tscal wall_length;
+        Tscal wall_width;
+        Tscal wall_thickness;
+
+        public:
+        GetParticlesInWall(
+            const Tvec &wall_pos, Tscal wall_length, Tscal wall_width, Tscal wall_thickness)
+            : wall_pos(wall_pos), wall_length(wall_length), wall_width(wall_width),
+              wall_thickness(wall_thickness) {}
+
+        struct Edges {
+            const shamrock::solvergraph::IFieldRefs<Tvec> &pos;
+            const shamrock::solvergraph::Indexes<u32> &sizes;
+            shamrock::solvergraph::IFieldSpan<u32> &part_ids_in_wall;
+        };
+
+        inline void set_edges(
+            std::shared_ptr<shamrock::solvergraph::IFieldRefs<Tvec>> pos,
+            std::shared_ptr<shamrock::solvergraph::Indexes<u32>> sizes,
+            std::shared_ptr<shamrock::solvergraph::IFieldSpan<u32>> part_ids_in_wall) {
+            __internal_set_ro_edges({pos, sizes});
+            __internal_set_rw_edges({part_ids_in_wall});
+        }
+
+        inline Edges get_edges() {
+            return Edges{
+                get_ro_edge<shamrock::solvergraph::IFieldRefs<Tvec>>(0),
+                get_ro_edge<shamrock::solvergraph::Indexes<u32>>(1),
+                get_rw_edge<shamrock::solvergraph::IFieldSpan<u32>>(0)};
+        }
+
+        void _impl_evaluate_internal();
+
+        inline virtual std::string _impl_get_label() const { return "GetParticlesInWall"; };
+
+        virtual std::string _impl_get_tex() const;
+    };
+} // namespace shammodels::sph::modules

src/shammodels/sph/src/Solver.cpp

@@ -52,6 +52,7 @@
 #include "shammodels/sph/modules/DiffOperator.hpp"
 #include "shammodels/sph/modules/DiffOperatorDtDivv.hpp"
 #include "shammodels/sph/modules/ExternalForces.hpp"
+#include "shammodels/sph/modules/GetParticlesInWall.hpp"
 #include "shammodels/sph/modules/GetParticlesOutsideSphere.hpp"
 #include "shammodels/sph/modules/IterateSmoothingLengthDensity.hpp"
 #include "shammodels/sph/modules/IterateSmoothingLengthDensityNeighLim.hpp"
@@ -97,6 +98,7 @@
 #include "shamsys/legacy/log.hpp"
 #include "shamtree/KarrasRadixTreeField.hpp"
 #include "shamtree/TreeTraversalCache.hpp"
+#include <shambackends/sycl.hpp>
 #include <memory>
 #include <stdexcept>
 #include <vector>
@@ -2017,6 +2019,81 @@ shammodels::sph::TimestepLog shammodels::sph::Solver<Tvec, Kern>::evolve_once()
 
         update_derivs();
 
+        bool has_walls = solver_config.particle_disable.disable_list.empty();
+        logger::raw_ln("value of has_wallas: ", has_walls);
+        if (has_walls) {
+            logger::raw_ln("@@@@@@@@@@@@ in wall disabling routine @@@@@@@@@@@@@@@");
+
+            // select phantom particles (create a mask of the ids of the particles to disable)
+            using namespace shamrock::solvergraph;
+            SolverGraph &solver_graph = storage.solver_graph;
+
+            auto xyz_edge = solver_graph.get_edge_ptr<FieldRefs<Tvec>>("xyz");
+
+            // get number of particle per patch: array of sizes of each patch
+            auto sizes = std::make_shared<shamrock::solvergraph::Indexes<u32>>("sizes", "Np");
+            scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) {
+                sizes->indexes.add_obj(p.id_patch, pdat.get_obj_cnt());
+            });
+
+            // create the actual mask field
+            auto part_to_disable_field = std::make_shared<shamrock::solvergraph::Field<u32>>(
+                1, "part_to_disable", "part_to_disable");
+
+            // properly initialize it for all patches
+            scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) {
+                auto &field = part_to_disable_field->get_field(p.id_patch);
+                field.resize(pdat.get_obj_cnt());
+            });
+
+            // now cast to the interface type through inheritance
+            auto part_to_disable = part_to_disable_field;
+
+            std::vector<std::shared_ptr<shamrock::solvergraph::INode>> part_disable_sequence{};
+
+            using disable_t    = typename ParticleDisableConfig<Tvec>::disable_t; // the types
+            using disable_wall = typename ParticleDisableConfig<Tvec>::Wall;
+
+            for (disable_t &disable_obj : solver_config.particle_disable.disable_list) {
+                if (disable_wall *disable_info = std::get_if<disable_wall>(&disable_obj)) {
+
+                    modules::GetParticlesInWall<Tvec> node_selector(
+                        disable_info->pos,
+                        disable_info->length,
+                        disable_info->width,
+                        disable_info->thickness);
+                    node_selector.set_edges(xyz_edge, sizes, part_to_disable);
+
+                    part_disable_sequence.push_back(
+                        std::make_shared<decltype(node_selector)>(std::move(node_selector)));
+                }
+            }
+
+            if (!part_disable_sequence.empty()) {
+                auto disable_sequence_node = solver_graph.register_node(
+                    "part_disable_selectors",
+                    OperationSequence("part disable selectors", std::move(part_disable_sequence)));
+
+                disable_sequence_node->evaluate();
+            }
+
+            // now set accelerations to 0 for these particles
+            // for this, let's reuse _impl_evaluate_internal() twice form Node Compute Omega
+            auto axyz_edge
+                = solver_graph.get_edge_ptr<shamrock::solvergraph::IFieldSpan<Tvec>>("axyz");
+            auto du_edge = solver_graph.get_edge_ptr<FieldRefs<Tscal>>("duint");
+
+            modules::SetWhenMask<Tscal> set_omega_mask{0};
+            set_omega_mask.set_edges(storage.part_counts, part_to_disable, du_edge);
+
+            set_omega_mask.evaluate();
+
+            modules::SetWhenMask<Tvec> set_omega_mask_vec({0.0, 0.0, 0.0});
+            set_omega_mask_vec.set_edges(storage.part_counts, part_to_disable, axyz_edge);
+            set_omega_mask_vec.evaluate();
+        }
+        // done disabling, now do the integration
+
         modules::ConservativeCheck<Tvec, Kern> cv_check(context, solver_config, storage);
         cv_check.check_conservation();
 

src/shammodels/sph/src/modules/ComputeOmega.cpp

@@ -107,6 +107,7 @@ std::string shammodels::sph::modules::SetWhenMask<T>::_impl_get_tex() const {
 }
 
 template class shammodels::sph::modules::SetWhenMask<f64>;
+template class shammodels::sph::modules::SetWhenMask<f64_3>;
 
 using namespace shammath;
 template class shammodels::sph::modules::NodeComputeOmega<f64_3, M4>;

src/shammodels/sph/src/modules/GetParticlesInWall.cpp

@@ -0,0 +1,103 @@
+// -------------------------------------------------------//
+//
+// 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
+//
+// -------------------------------------------------------//
+
+/**
+ * @file GetParticlesInWall.cpp
+ * @author Yona Lapeyre (yona.lapeyre@ens-lyon.fr)
+ * @brief Implements the GetParticlesInWall module, which identifies particles in a wall.
+ *
+ */
+
+#include "shambase/stacktrace.hpp"
+#include "shambase/string.hpp"
+#include "shambackends/kernel_call_distrib.hpp"
+#include "shammodels/sph/modules/GetParticlesInWall.hpp"
+
+namespace shammodels::sph::modules {
+
+    template<typename Tvec>
+    void GetParticlesInWall<Tvec>::_impl_evaluate_internal() {
+        StackEntry stack_loc{};
+
+        auto edges = get_edges();
+
+        auto &thread_counts = edges.sizes.indexes;
+
+        edges.pos.check_sizes(thread_counts);
+        edges.part_ids_in_wall.ensure_sizes(thread_counts);
+
+        auto &positions        = edges.pos.get_spans();
+        auto &part_ids_in_wall = edges.part_ids_in_wall.get_spans();
+
+        auto dev_sched = shamsys::instance::get_compute_scheduler_ptr();
+
+        sham::distributed_data_kernel_call(
+            dev_sched,
+            sham::DDMultiRef{positions},
+            sham::DDMultiRef{part_ids_in_wall},
+            thread_counts,
+            [wall_pos       = this->wall_pos,
+             wall_length    = this->wall_length,
+             wall_width     = this->wall_width,
+             wall_thickness = this->wall_thickness](
+                u32 i, const Tvec *__restrict pos, u32 *__restrict part_ids_in_wall) {
+                Tscal x = pos[i][0];
+                Tscal y = pos[i][1];
+                Tscal z = pos[i][2];
+
+                Tscal x0 = wall_pos[0];
+                Tscal y0 = wall_pos[1];
+                Tscal z0 = wall_pos[2];
+
+                bool in_wall = (x - x0 < wall_length) && (x - x0 > 0) && (y - y0 < wall_width)
+                               && (y - y0 > 0) && (z - z0 < wall_thickness) && (z - z0 > 0);
+
+                if (in_wall) {
+                    part_ids_in_wall[i] = 1;
+                } else {
+                    part_ids_in_wall[i] = 0;
+                }
+            });
+    }
+
+    /**
+     * @brief Returns the tex string for the GetParticlesInWall module.
+     *
+     * This module identifies particles inside a rectangular wall.
+     *
+     * @param[in] pos The position field.
+     * @param[in] part_ids_in_wall The particle id field.
+     * @return The tex string.
+     */
+    template<typename Tvec>
+    std::string GetParticlesInWall<Tvec>::_impl_get_tex() const {
+        auto pos              = get_ro_edge_base(0).get_tex_symbol();
+        auto part_ids_in_wall = get_rw_edge_base(0).get_tex_symbol();
+
+        std::string tex = R"tex(
+    Get particles inside the rectangular wall
+
+    \begin{align}
+    {part_ids_in_wall} &= \{i \text{ where } \vert\vert{x}_i - {x0}\vert\vert < {wall_length} \text{ and } \vert\vert{y}_i - {y0}\vert\vert < {wall_width} \text{ and } \vert\vert{z}_i - {z0}\vert\vert < {wall_thickness}\}\\
+    \end{align}
+    )tex";
+
+        shambase::replace_all(tex, "{x0}", std::to_string(wall_pos[0]));
+        shambase::replace_all(tex, "{y0}", std::to_string(wall_pos[1]));
+        shambase::replace_all(tex, "{z0}", std::to_string(wall_pos[2]));
+        shambase::replace_all(tex, "{wall_length}", std::to_string(wall_length));
+        shambase::replace_all(tex, "{wall_width}", std::to_string(wall_width));
+        shambase::replace_all(tex, "{wall_thickness}", std::to_string(wall_thickness));
+
+        return tex;
+    }
+
+} // namespace shammodels::sph::modules
+
+template class shammodels::sph::modules::GetParticlesInWall<f64_3>;

src/shammodels/sph/src/pySPHModel.cpp

@@ -745,6 +745,13 @@ void add_instance(py::module &m, std::string name_config, std::string name_model
             [](T &self) {
                 return self.solver.solver_config;
             })
+        .def(
+            "add_wall",
+            [](T &self, Tvec pos, Tscal length, Tscal width, Tscal thickness) {
+                ParticleDisableConfig<Tvec> config;
+                self.solver.solver_config.particle_disable.add_disable_wall(
+                    pos, length, width, thickness);
+            })
         .def("set_solver_config", &T::set_solver_config)
         .def("add_sink", &T::add_sink)
         .def(