shared node value consistency

pyphare/pyphare/pharein/simulation.py

@@ -317,9 +317,11 @@ def get_max_ghosts():
         return max(grid.nbrGhosts(kwargs["interp_order"], x) for x in ['primal','dual'])
 
     max_ghosts = get_max_ghosts()
-    small_invalid_patch_size = phare_utilities.np_array_ify(max_ghosts - 1, ndim)
+    small_invalid_patch_size = phare_utilities.np_array_ify(max_ghosts, ndim)
     largest_patch_size =  kwargs.get("largest_patch_size", None)
-    smallest_patch_size = phare_utilities.np_array_ify(max_ghosts, ndim)
+
+    # to prevent primal ghost overlaps of non adjacent patches, we need smallest_patch_size+=1
+    smallest_patch_size = phare_utilities.np_array_ify(max_ghosts, ndim) + 1
     if "smallest_patch_size" in kwargs and kwargs["smallest_patch_size"] is not None:
         smallest_patch_size = phare_utilities.np_array_ify(kwargs["smallest_patch_size"], ndim)
 

src/amr/CMakeLists.txt

@@ -36,6 +36,7 @@ set( SOURCES_INC
      messengers/mhd_hybrid_messenger_strategy.h
      messengers/mhd_messenger.h
      messengers/messenger_factory.h
+     messengers/ghost_only_variable_fill_pattern.cpp
      messengers/messenger_info.h
      messengers/hybrid_messenger_info.h
      messengers/mhd_messenger_info.h

src/amr/data/field/field_geometry.h

@@ -9,6 +9,24 @@
 #include "core/data/grid/gridlayout_impl.h"
 #include "field_overlap.h"
 
+
+
+namespace PHARE::amr
+{
+class AFieldGeometry : public SAMRAI::hier::BoxGeometry
+{
+public:
+    virtual ~AFieldGeometry() {}
+    AFieldGeometry(SAMRAI::hier::Box const& box)
+        : patchBox{box}
+    {
+    }
+
+    SAMRAI::hier::Box const patchBox;
+};
+
+} // namespace PHARE::amr
+
 namespace PHARE
 {
 namespace amr
@@ -17,7 +35,7 @@ namespace amr
     /**
      * @brief The FieldGeometry class
      */
-    class FieldGeometry : public SAMRAI::hier::BoxGeometry
+    class FieldGeometry : public AFieldGeometry
     {
     public:
         static constexpr std::size_t dimension    = GridLayoutT::dimension;
@@ -27,7 +45,8 @@ namespace amr
          * with a temporary gridlayout
          */
         FieldGeometry(SAMRAI::hier::Box const& box, GridLayoutT layout, PhysicalQuantity qty)
-            : ghostBox_{toFieldBox(box, qty, layout)}
+            : AFieldGeometry{box}
+            , ghostBox_{toFieldBox(box, qty, layout)}
             , interiorBox_{toFieldBox(box, qty, layout, false)}
             , layout_{std::move(layout)}
             , quantity_{qty}
@@ -222,7 +241,6 @@ namespace amr
         PhysicalQuantity quantity_;
 
 
-
         /*** \brief Compute destination box representing the intersection of two geometry
          *
          *   \param destinationBoxes BoxContainer that will be filled of box
@@ -286,10 +304,10 @@ namespace amr
 
             SAMRAI::hier::Box const together(destinationBox * sourceBox * fillField);
 
-
             // if the interesection is not empty we either push it into the container
             // if we don't want to fill the interior we remove it from the intersection
             // which may add multiple boxes to the container.
+
             if (!together.empty())
             {
                 if (overwriteInterior)

src/amr/messengers/ghost_only_variable_fill_pattern.cpp

@@ -0,0 +1,10 @@
+
+
+#include "ghost_only_variable_fill_pattern.h"
+
+namespace PHARE::amr
+{
+// the value of this string cannot change, SAMRAI breaks :|
+const std::string GhostOnlyVariablyFillPattern::s_name_id = "BOX_GEOMETRY_FILL_PATTERN";
+
+} // namespace PHARE::amr

src/amr/messengers/ghost_only_variable_fill_pattern.h

@@ -0,0 +1,97 @@
+
+#ifndef GhostOnlyVariablyFillPattern_
+#define GhostOnlyVariablyFillPattern_
+
+#include <cassert>
+#include "SAMRAI/xfer/VariableFillPattern.h"
+#include "amr/data/field/field_geometry.h"
+
+namespace PHARE::amr
+{
+class GhostOnlyVariablyFillPattern : public SAMRAI::xfer::VariableFillPattern
+{
+public:
+    GhostOnlyVariablyFillPattern() {}
+    virtual ~GhostOnlyVariablyFillPattern() {}
+
+    std::shared_ptr<SAMRAI::hier::BoxOverlap>
+    calculateOverlap(const SAMRAI::hier::BoxGeometry& dst_geometry,
+                     const SAMRAI::hier::BoxGeometry& src_geometry,
+                     const SAMRAI::hier::Box& dst_patch_box, const SAMRAI::hier::Box& src_mask,
+                     const SAMRAI::hier::Box& fill_box, const bool overwrite_interior_,
+                     const SAMRAI::hier::Transformation& transformation) const
+    {
+#ifndef DEBUG_CHECK_DIM_ASSERTIONS
+        NULL_USE(dst_patch_box);
+#endif
+        TBOX_ASSERT_OBJDIM_EQUALITY2(dst_patch_box, src_mask);
+
+        bool overwrite_interior = true; // replace func param
+        assert(overwrite_interior_ == overwrite_interior);
+
+        auto& dst_cast = dynamic_cast<AFieldGeometry const&>(dst_geometry);
+        auto& src_cast = dynamic_cast<AFieldGeometry const&>(src_geometry);
+
+        // for shared border node value sync
+        if (src_cast.patchBox.getGlobalId().getOwnerRank()
+            != dst_cast.patchBox.getGlobalId().getOwnerRank())
+            overwrite_interior = src_cast.patchBox.getLocalId() > dst_cast.patchBox.getLocalId();
+
+        return dst_geometry.calculateOverlap(src_geometry, src_mask, fill_box, overwrite_interior,
+                                             transformation);
+    }
+
+    const std::string& getPatternName() const { return s_name_id; }
+
+private:
+    GhostOnlyVariablyFillPattern(const GhostOnlyVariablyFillPattern&);            // not implemented
+    GhostOnlyVariablyFillPattern& operator=(const GhostOnlyVariablyFillPattern&); // not implemented
+
+    static const std::string s_name_id; // = "GHOST_ONLY_FILL_PATTERN";
+
+    const SAMRAI::hier::IntVector& getStencilWidth()
+    {
+        TBOX_ERROR("getStencilWidth() should not be\n"
+                   << "called.  This pattern creates overlaps based on\n"
+                   << "the BoxGeometry objects and is not restricted to a\n"
+                   << "specific stencil.\n");
+
+        /*
+         * Dummy return value that will never get reached.
+         */
+        return SAMRAI::hier::IntVector::getZero(SAMRAI::tbox::Dimension(1));
+    }
+
+    /*
+     *************************************************************************
+     *
+     * Compute BoxOverlap that specifies data to be filled by refinement
+     * operator.
+     *
+     *************************************************************************
+     */
+    std::shared_ptr<SAMRAI::hier::BoxOverlap>
+    computeFillBoxesOverlap(const SAMRAI::hier::BoxContainer& fill_boxes,
+                            const SAMRAI::hier::BoxContainer& node_fill_boxes,
+                            const SAMRAI::hier::Box& patch_box, const SAMRAI::hier::Box& data_box,
+                            const SAMRAI::hier::PatchDataFactory& pdf) const
+    {
+        NULL_USE(node_fill_boxes);
+
+        /*
+         * For this (default) case, the overlap is simply the intersection of
+         * fill_boxes and data_box.
+         */
+        SAMRAI::hier::Transformation transformation(
+            SAMRAI::hier::IntVector::getZero(patch_box.getDim()));
+
+        SAMRAI::hier::BoxContainer overlap_boxes(fill_boxes);
+        overlap_boxes.intersectBoxes(data_box);
+
+        return pdf.getBoxGeometry(patch_box)->setUpOverlap(overlap_boxes, transformation);
+    }
+};
+
+} // namespace PHARE::amr
+
+#endif /*GhostOnlyVariablyFillPattern*/

src/amr/messengers/quantity_communicator.h

@@ -23,6 +23,7 @@
 #include <optional>
 
 
+#include "ghost_only_variable_fill_pattern.h"
 
 namespace PHARE
 {
@@ -48,6 +49,18 @@ namespace amr
     {
     };
 
+    class XGhostOnlyVariablyFillPattern : public GhostOnlyVariablyFillPattern
+    {
+    };
+
+    class YGhostOnlyVariablyFillPattern : public GhostOnlyVariablyFillPattern
+    {
+    };
+
+    class ZGhostOnlyVariablyFillPattern : public GhostOnlyVariablyFillPattern
+    {
+    };
+
     struct Refiner
     {
         using schedule_type  = SAMRAI::xfer::RefineSchedule;
@@ -125,6 +138,9 @@ namespace amr
          */
         void add(std::shared_ptr<Schedule> schedule, int levelNumber)
         {
+            // for shared border node value sync
+            schedule->setDeterministicUnpackOrderingFlag(true);
+
             schedules_[levelNumber] = std::move(schedule);
         }
     };
@@ -160,11 +176,12 @@ namespace amr
                 std::shared_ptr<SAMRAI::hier::TimeInterpolateOperator> timeOp)
     {
         std::shared_ptr<SAMRAI::xfer::VariableFillPattern> xVariableFillPattern
-            = std::make_shared<XVariableFillPattern>();
+            = std::make_shared<XGhostOnlyVariablyFillPattern>();
         std::shared_ptr<SAMRAI::xfer::VariableFillPattern> yVariableFillPattern
-            = std::make_shared<YVariableFillPattern>();
+            = std::make_shared<YGhostOnlyVariablyFillPattern>();
         std::shared_ptr<SAMRAI::xfer::VariableFillPattern> zVariableFillPattern
-            = std::make_shared<ZVariableFillPattern>();
+            = std::make_shared<ZGhostOnlyVariablyFillPattern>();
+
         Communicator<Refiner> com;
 
         auto registerRefine

tests/simulator/__init__.py

@@ -26,7 +26,7 @@ def basicSimulatorArgs(dim: int, interp: int, **kwargs):
     b0 = [[3] * dim, [8] * dim]
     args = {
         "interp_order": interp,
-        "smallest_patch_size": [5] * dim,
+        "smallest_patch_size": [6] * dim,
         "largest_patch_size": [20] * dim,
         "time_step_nbr": 1000,
         "final_time": 1.0,

tests/simulator/refined_particle_nbr.py

@@ -20,6 +20,7 @@ class SimulatorRefinedParticleNbr(unittest.TestCase):
 
     def __init__(self, *args, **kwargs):
         super(SimulatorRefinedParticleNbr, self).__init__(*args, **kwargs)
+        self.simulator = None
         with open(os.path.join(project_root, "res/amr/splitting.yml"), 'r') as stream:
             try:
                 self.yaml_root = yaml.safe_load(stream)

tests/simulator/test_advance.py

@@ -2,6 +2,7 @@
 cpp = cpp_lib()
 
 from pyphare.simulator.simulator import Simulator, startMPI
+from pyphare.core.phare_utilities import np_array_ify
 from pyphare.pharesee.hierarchy import hierarchy_from, merge_particles
 from pyphare.pharein import MaxwellianFluidModel
 from pyphare.pharein.diagnostics import ParticleDiagnostics, FluidDiagnostics, ElectromagDiagnostics
@@ -26,9 +27,9 @@ def ddt_test_id(self):
 
     def getHierarchy(self, interp_order, refinement_boxes, qty, nbr_part_per_cell=100,
                      diag_outputs="phare_outputs",
-                     smallest_patch_size=5, largest_patch_size=20,
+                     smallest_patch_size=6, largest_patch_size=20,
                      cells=120, time_step=0.001, model_init={},
-                     dl=0.1, extra_diag_options={}, time_step_nbr=1, timestamps=None):
+                     dl=0.1, extra_diag_options={}, time_step_nbr=1, timestamps=None, ndim=1):
 
         from pyphare.pharein import global_vars
         global_vars.sim = None
@@ -41,8 +42,8 @@ def getHierarchy(self, interp_order, refinement_boxes, qty, nbr_part_per_cell=10
             time_step_nbr=time_step_nbr,
             time_step=time_step,
             boundary_types="periodic",
-            cells=cells,
-            dl=dl,
+            cells=np_array_ify(cells, ndim),
+            dl=np_array_ify(dl, ndim),
             interp_order=interp_order,
             refinement_boxes=refinement_boxes,
             diag_options={"format": "phareh5",
@@ -210,9 +211,9 @@ def _test_overlaped_fields_are_equal(self, time_step, time_step_nbr, datahier):
                         slice2 = data2[loc_b2.lower[0]:loc_b2.upper[0] + 1]
 
                         try:
-                            np.testing.assert_allclose(slice1, slice2, atol=1e-6)
+                            np.testing.assert_equal(slice1, slice2)
                         except AssertionError as e:
-                            print("error", coarsest_time, overlap)
+                            print("error", pd1.name, coarsest_time, overlap)
                             raise e
 
         self.assertGreater(check, time_step_nbr)
@@ -235,30 +236,28 @@ def test_overlaped_fields_are_equal(self, refinement_boxes):
 
 
     @data(
+        {},
         {"L0": [Box1D(10, 19)]},
-        # {"L0": [Box2D(10, 19)]},
-        # {"L0": [Box3D(10, 19)]},
     )
     def test_overlaped_fields_are_equal_with_min_max_patch_size_of_max_ghosts(self, refinement_boxes):
         from pyphare.pharein.simulation import check_patch_size
 
-        dim = refinement_boxes["L0"][0].ndim
+        dim = 1
 
-        cells = [30] * dim
+        cells = [60] * dim
         time_step_nbr=3
         time_step=0.001
         diag_outputs=f"phare_overlaped_fields_are_equal_with_min_max_patch_size_of_max_ghosts{self.ddt_test_id()}"
         for interp_order in [1, 2, 3]:
             largest_patch_size, smallest_patch_size = check_patch_size(dim, interp_order=interp_order, cells=cells)
             datahier = self.getHierarchy(interp_order, refinement_boxes, "eb", diag_outputs=diag_outputs,
                                       smallest_patch_size=smallest_patch_size, largest_patch_size=smallest_patch_size,
-                                      time_step=time_step, time_step_nbr=time_step_nbr)
+                                      time_step=time_step, time_step_nbr=time_step_nbr, cells=cells)
             self._test_overlaped_fields_are_equal(time_step, time_step_nbr, datahier)
 
 
 
 
-
     def _test_patch_ghost_particle_are_clone_of_overlaped_patch_domain_particles(self, dim, interp_order, refinement_boxes):
         print("test_patch_ghost_particle_are_clone_of_overlaped_patch_domain_particles")
         print("interporder : {}".format(interp_order))
@@ -452,7 +451,7 @@ def _test_field_coarsening_via_subcycles(self, dim, interp_order, refinement_box
         datahier = self.getHierarchy(interp_order, refinement_boxes, "eb", cells=30,
                                       diag_outputs=diag_outputs, time_step=0.001,
                                       extra_diag_options={"fine_dump_lvl_max": 10},
-                                      time_step_nbr=time_step_nbr, smallest_patch_size=5,
+                                      time_step_nbr=time_step_nbr,
                                       largest_patch_size=30)
 
         lvl_steps = global_vars.sim.level_time_steps
@@ -556,7 +555,7 @@ def _test_field_level_ghosts_via_subcycles_and_coarser_interpolation(self, ndim,
 
         def _getHier(diag_dir, boxes=[]):
             return self.getHierarchy(interp_order, boxes, "eb", cells=30,
-                time_step_nbr=1, smallest_patch_size=5, largest_patch_size=30,
+                time_step_nbr=1,  largest_patch_size=30,
                 diag_outputs=diag_dir, extra_diag_options={"fine_dump_lvl_max": 10}, time_step=0.001,
                 model_init={"seed": rando}
             )
@@ -657,7 +656,7 @@ def test_hierarchy_timestamp_cadence(self, refinement_boxes):
                 diag_outputs=f"phare_outputs_hierarchy_timestamp_cadence_{self.ddt_test_id()}_{i}"
                 hier = self.getHierarchy(interp_order=1, refinement_boxes=refinement_boxes, qty="eb", cells=30,
                                               diag_outputs=diag_outputs, time_step=time_step,
-                                              time_step_nbr=time_step_nbr, smallest_patch_size=5,
+                                              time_step_nbr=time_step_nbr,
                                               largest_patch_size=30, timestamps=timestamps)
 
                 time_hier_keys = list(hier.time_hier.keys())

tests/simulator/test_initialization.py

@@ -551,8 +551,8 @@ def test_no_patch_ghost_on_refined_level_case(self, simInput):
     _has_patch_ghost_on_refined_level_case = (
       {
         "cells": 40,
-        "smallest_patch_size": 5,
-        "largest_patch_size": 5},
+        "smallest_patch_size": 6,
+        "largest_patch_size": 6},
     )
     @data(*_has_patch_ghost_on_refined_level_case)
     def test_has_patch_ghost_on_refined_level_case(self, simInput):