[SPH] fix potential interaction asymmetry in ghosts

exemples/benchmarks/sph_weak_scale_test.py

@@ -51,7 +51,8 @@
 )
 cfg.set_boundary_periodic()
 cfg.set_eos_adiabatic(gamma)
-cfg.set_max_neigh_cache_size(int(100e9))
+cfg.set_filter_empty_patch_gz(False)
+
 cfg.print_status()
 model.set_solver_config(cfg)
 model.init_scheduler(scheduler_split_val, scheduler_merge_val)

src/shammodels/sph/include/shammodels/sph/BasicSPHGhosts.hpp

@@ -105,7 +105,8 @@ namespace shammodels::sph {
         GeneratorMap find_interfaces(
             SerialPatchTree<vec> &sptree,
             shamrock::patch::PatchtreeField<flt> &int_range_max_tree,
-            shamrock::patch::PatchField<flt> &int_range_max);
+            shamrock::patch::PatchField<flt> &int_range_max,
+            bool filter_empty_patch_gz);
 
         /**
          * @brief precompute interfaces members and cache result in the return
@@ -131,11 +132,12 @@ namespace shammodels::sph {
         CacheMap make_interface_cache(
             SerialPatchTree<vec> &sptree,
             shamrock::patch::PatchtreeField<flt> &int_range_max_tree,
-            shamrock::patch::PatchField<flt> &int_range_max) {
+            shamrock::patch::PatchField<flt> &int_range_max,
+            bool filter_empty_patch_gz) {
             StackEntry stack_loc{};
 
             return gen_id_table_interfaces(
-                find_interfaces(sptree, int_range_max_tree, int_range_max));
+                find_interfaces(sptree, int_range_max_tree, int_range_max, filter_empty_patch_gz));
         }
 
         /**

src/shammodels/sph/include/shammodels/sph/SPHUtilities.hpp

@@ -77,7 +77,10 @@ namespace shammodels::sph {
         SPHUtilities(PatchScheduler &sched) : sched(sched) {}
 
         inline InterfBuildCache build_interf_cache(
-            GhostHndl &interf_handle, SerialPatchTree<vec> &sptree, flt h_evol_max) {
+            GhostHndl &interf_handle,
+            SerialPatchTree<vec> &sptree,
+            flt h_evol_max,
+            bool filter_empty_patch_gz) {
 
             using namespace shamrock::patch;
 
@@ -100,7 +103,8 @@ namespace shammodels::sph {
                     return sham::max_8points(h0, h1, h2, h3, h4, h5, h6, h7);
                 });
 
-            return interf_handle.make_interface_cache(sptree, interactR_mpi_tree, interactR_patch);
+            return interf_handle.make_interface_cache(
+                sptree, interactR_mpi_tree, interactR_patch, filter_empty_patch_gz);
         }
 
         static void iterate_smoothing_length_cache(

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

@@ -277,6 +277,11 @@ struct shammodels::sph::SolverConfig {
 
     inline void set_particle_tracking(bool state) { track_particles_id = state; }
 
+    bool filter_empty_patch_gz = false;
+    inline void set_filter_empty_patch_gz(bool state) { filter_empty_patch_gz = state; }
+
+    inline bool has_filter_empty_patch_gz() { return filter_empty_patch_gz; }
+
     //////////////////////////////////////////////////////////////////////////////////////////////
     // Units Config
     //////////////////////////////////////////////////////////////////////////////////////////////

src/shammodels/sph/src/BasicSPHGhosts.cpp

@@ -254,11 +254,14 @@ inline void for_each_patch_shift(
 
 using namespace shammodels::sph;
 
+#define NEW_GZ_INTERACT_CRIT
+
 template<class vec>
 auto BasicSPHGhostHandler<vec>::find_interfaces(
     SerialPatchTree<vec> &sptree,
     shamrock::patch::PatchtreeField<flt> &int_range_max_tree,
-    shamrock::patch::PatchField<flt> &int_range_max) -> GeneratorMap {
+    shamrock::patch::PatchField<flt> &int_range_max,
+    bool filter_empty_patch_gz) -> GeneratorMap {
 
     StackEntry stack_loc{};
 
@@ -306,21 +309,42 @@ auto BasicSPHGhostHandler<vec>::find_interfaces(
 
                         sptree.host_for_each_leafs(
                             [&](u64 tree_id, PtNode n) {
+#ifdef NEW_GZ_INTERACT_CRIT
+                                // The ghost interaction was not symmetric before
+                                // now it is since we check for the max of the two h_max
+                                flt receiv_h_max = acc_tf[tree_id];
+                                shammath::AABB<vec> tree_cell{n.box_min, n.box_max};
+                                tree_cell
+                                    = tree_cell.expand_all(sham::max(sender_h_max, receiv_h_max));
+
+                                return tree_cell
+                                    .get_intersect(
+                                        shammath::AABB<vec>{
+                                            sender_bsize_off.lower, sender_bsize_off.upper})
+                                    .is_not_empty();
+#else
                                 flt receiv_h_max = acc_tf[tree_id];
                                 CoordRange<vec> receiv_exp{
                                     n.box_min - receiv_h_max, n.box_max + receiv_h_max};
 
                                 return receiv_exp.get_intersect(sender_bsize_off).is_not_empty();
+#endif
                             },
                             [&](u64 id_found, PtNode n) {
                                 if ((id_found == psender.id_patch) && (xoff == 0) && (yoff == 0)
                                     && (zoff == 0)) {
                                     return;
                                 }
 
+#ifdef NEW_GZ_INTERACT_CRIT
+                                CoordRange<vec> receiv_exp
+                                    = CoordRange<vec>{n.box_min, n.box_max}.expand_all(
+                                        sham::max(sender_h_max, int_range_max.get(id_found)));
+#else
                                 CoordRange<vec> receiv_exp
                                     = CoordRange<vec>{n.box_min, n.box_max}.expand_all(
                                         int_range_max.get(id_found));
+#endif
 
                                 CoordRange<vec> interf_volume = sender_bsize.get_intersect(
                                     receiv_exp.add_offset(-periodic_offset));
@@ -361,21 +385,40 @@ auto BasicSPHGhostHandler<vec>::find_interfaces(
 
                 sptree.host_for_each_leafs(
                     [&](u64 tree_id, PtNode n) {
+#ifdef NEW_GZ_INTERACT_CRIT
+                        // The ghost interaction was not symmetric before
+                        // now it is since we check for the max of the two h_max
+                        flt receiv_h_max = acc_tf[tree_id];
+                        shammath::AABB<vec> tree_cell{n.box_min, n.box_max};
+                        tree_cell = tree_cell.expand_all(sham::max(sender_h_max, receiv_h_max));
+
+                        return tree_cell
+                            .get_intersect(
+                                shammath::AABB<vec>{sender_bsize_off.lower, sender_bsize_off.upper})
+                            .is_not_empty();
+#else
                         flt receiv_h_max = acc_tf[tree_id];
                         CoordRange<vec> receiv_exp{
                             n.box_min - receiv_h_max, n.box_max + receiv_h_max};
 
                         return receiv_exp.get_intersect(sender_bsize_off).is_not_empty();
+#endif
                     },
                     [&](u64 id_found, PtNode n) {
                         if ((id_found == psender.id_patch) && (xoff == 0) && (yoff == 0)
                             && (zoff == 0)) {
                             return;
                         }
 
+#ifdef NEW_GZ_INTERACT_CRIT
+                        CoordRange<vec> receiv_exp
+                            = CoordRange<vec>{n.box_min, n.box_max}.expand_all(
+                                sham::max(sender_h_max, int_range_max.get(id_found)));
+#else
                         CoordRange<vec> receiv_exp
                             = CoordRange<vec>{n.box_min, n.box_max}.expand_all(
                                 int_range_max.get(id_found));
+#endif
 
                         CoordRange<vec> interf_volume
                             = sender_bsize.get_intersect(receiv_exp.add_offset(-offset));
@@ -413,18 +456,36 @@ auto BasicSPHGhostHandler<vec>::find_interfaces(
 
             sptree.host_for_each_leafs(
                 [&](u64 tree_id, PtNode n) {
+#ifdef NEW_GZ_INTERACT_CRIT
+                    // The ghost interaction was not symmetric before
+                    // now it is since we check for the max of the two h_max
+                    flt receiv_h_max = acc_tf[tree_id];
+                    shammath::AABB<vec> tree_cell{n.box_min, n.box_max};
+                    tree_cell = tree_cell.expand_all(sham::max(sender_h_max, receiv_h_max));
+
+                    return tree_cell
+                        .get_intersect(
+                            shammath::AABB<vec>{sender_bsize_off.lower, sender_bsize_off.upper})
+                        .is_not_empty();
+#else
                     flt receiv_h_max = acc_tf[tree_id];
                     CoordRange<vec> receiv_exp{n.box_min - receiv_h_max, n.box_max + receiv_h_max};
 
                     return receiv_exp.get_intersect(sender_bsize_off).is_not_empty();
+#endif
                 },
                 [&](u64 id_found, PtNode n) {
                     if (id_found == psender.id_patch) {
                         return;
                     }
 
+#ifdef NEW_GZ_INTERACT_CRIT
+                    CoordRange<vec> receiv_exp = CoordRange<vec>{n.box_min, n.box_max}.expand_all(
+                        sham::max(sender_h_max, int_range_max.get(id_found)));
+#else
                     CoordRange<vec> receiv_exp = CoordRange<vec>{n.box_min, n.box_max}.expand_all(
                         int_range_max.get(id_found));
+#endif
 
                     CoordRange<vec> interf_volume
                         = sender_bsize.get_intersect(receiv_exp.add_offset(-periodic_offset));
@@ -441,11 +502,57 @@ auto BasicSPHGhostHandler<vec>::find_interfaces(
         });
     }
 
-    // interf_map.for_each([](u64 sender, u64 receiver, InterfaceBuildInfos build){
-    //     logger::raw_ln("found interface
-    //     :",sender,"->",receiver,"ratio:",build.volume_ratio,
-    //     "volume:",build.cut_volume.lower,build.cut_volume.upper);
-    // });
+    if (filter_empty_patch_gz) {
+        PatchField<u32> patchdata_obj_cnt = sched.map_owned_to_patch_field_simple<u32>(
+            [&](const Patch p, PatchDataLayer &pdat) -> u32 {
+                return pdat.get_obj_cnt();
+            });
+
+        GeneratorMap interf_map2 = {};
+
+        interf_map.for_each([&](u64 sender, u64 receiver, InterfaceBuildInfos build) {
+            u32 sender_obj_cnt   = patchdata_obj_cnt.get(sender);
+            u32 receiver_obj_cnt = patchdata_obj_cnt.get(receiver);
+
+            if (sender_obj_cnt == 0 || receiver_obj_cnt == 0) {
+                logger::debug_ln("Ghost", "patch ", sender, "->", receiver, "is empty, skipping");
+                return;
+            }
+
+            interf_map2.add_obj(sender, receiver, InterfaceBuildInfos{build});
+        });
+        interf_map = interf_map2;
+    }
+
+    u32 total          = 0;
+    u32 empty_receiver = 0;
+    u32 empty_sender   = 0;
+    u32 empty_either   = 0;
+    interf_map.for_each([&](u64 sender, u64 receiver, InterfaceBuildInfos build) {
+        // logger::raw_ln("found interface :",sender,"->",receiver,"ratio:",build.volume_ratio,
+        //"volume:",build.cut_volume.lower,build.cut_volume.upper,
+        // sched.get_patch(receiver).load_value);
+        total++;
+        if (sched.get_patch(receiver).load_value == 0) {
+            empty_receiver++;
+        }
+        if (sched.get_patch(sender).load_value == 0) {
+            empty_sender++;
+        }
+        if (sched.get_patch(receiver).load_value == 0 || sched.get_patch(sender).load_value == 0) {
+            empty_either++;
+        }
+    });
+    logger::debug_ln(
+        "Ghost",
+        "total:",
+        total,
+        "empty_receiver:",
+        empty_receiver,
+        "empty_sender:",
+        empty_sender,
+        "empty_either:",
+        empty_either);
 
     return interf_map;
 }

src/shammodels/sph/src/Solver.cpp

@@ -362,7 +362,8 @@ void shammodels::sph::Solver<Tvec, Kern>::build_ghost_cache() {
     storage.ghost_patch_cache.set(sph_utils.build_interf_cache(
         storage.ghost_handler.get(),
         storage.serial_patch_tree.get(),
-        solver_config.htol_up_coarse_cycle));
+        solver_config.htol_up_coarse_cycle,
+        solver_config.has_filter_empty_patch_gz()));
 
     // storage.ghost_handler.get().gen_debug_patch_ghost(storage.ghost_patch_cache.get());
 }

src/shammodels/sph/src/pySPHModel.cpp

@@ -69,6 +69,7 @@ void add_instance(py::module &m, std::string name_config, std::string name_model
                               "    -> calling this is a no-op,\n"
                               "    -> you can remove the call to that function"););
             })
+        .def("set_filter_empty_patch_gz", &TConfig::set_filter_empty_patch_gz)
         .def("set_smoothing_length_density_based", &TConfig::set_smoothing_length_density_based)
         .def(
             "set_smoothing_length_density_based_neigh_lim",

src/shamrock/include/shamrock/scheduler/PatchScheduler.hpp

@@ -287,6 +287,10 @@ class PatchScheduler {
         });
     }
 
+    inline const shamrock::patch::Patch &get_patch(u64 patch_id) const {
+        return patch_list.global[patch_list.id_patch_to_global_idx.at(patch_id)];
+    }
+
     inline u32 get_patch_rank_owner(u64 patch_id) {
         shamrock::patch::Patch &cur_p
             = patch_list.global[patch_list.id_patch_to_global_idx.at(patch_id)];