SAMRAI: Algorithm/schedule per quantity component (#589)

pyphare/pyphare/core/gridlayout.py

@@ -3,7 +3,7 @@
 import math
 import numpy as np
 
-from .phare_utilities import listify
+from .phare_utilities import is_scalar, listify
 from .box import Box
 
 directions = ["x", "y", "z"]
@@ -204,7 +204,23 @@ def allocSizeDerived(self, interpOrder, centering, nbrCells):
         return size
 
 
+    def AMRPointToLocal(self, point):
+        return point - self.box.lower + self.physicalStartIndex(self.interp_order, "dual")
+
+    def AMRBoxToLocal(self, box):
+        local = box.copy()
+        local.lower = self.AMRPointToLocal(local.lower)
+        local.upper = self.AMRPointToLocal(local.upper)
+        return local
+
+    def AMRToLocal(self, toLocal):
+        assert not is_scalar(toLocal)
+        if isinstance(toLocal, Box):
+            return self.AMRBoxToLocal(toLocal)
+        return self.AMRPointToLocal(toLocal)
+
     def AMRIndexToLocal(self, dim, index):
+        assert is_scalar(index)
         dualStart = self.physicalStartIndex(self.interp_order, "dual")
         return index - self.box.lower[dim] + dualStart
 

pyphare/pyphare/core/phare_utilities.py

@@ -143,3 +143,19 @@ def print_trace():
     import sys, traceback
     _, _, tb = sys.exc_info()
     traceback.print_tb(tb)
+
+
+def deep_copy(item, memo, excludes=[]):
+    from copy import deepcopy
+
+    clazz = item.__class__
+    that = clazz.__new__(clazz)
+    memo[id(item)] = that
+    for key, value in item.__dict__.items():
+        # some objects are not picklable so cannot be deepcopied eg. h5py datasets
+        if key in excludes:
+            setattr(that, key, value)
+        else:
+            setattr(that, key, deepcopy(value, memo))
+    return that
+

pyphare/pyphare/pharesee/geometry.py

@@ -412,18 +412,18 @@ def level_ghost_boxes(hierarchy, quantities, levelNbrs=[], time=None):
 
                 for gabox in ghostAreaBoxes:
 
-                    remaining = boxm.remove(gabox, check_patches[0].box)
+                    remaining = gabox - check_patches[0].box
 
                     for patch in check_patches[1:]:
-                        tmp = remaining
+                        tmp = []
                         remove = []
                         for i, rem in enumerate(remaining):
                             if rem * patch.box is not None:
-                                tmp += boxm.remove(rem, patch.box)
                                 remove.append(i)
-                        for rm in remove:
-                            del tmp[rm]
-                        remaining = tmp
+                                tmp += rem - patch.box
+                        for rm in reversed(remove):
+                            del remaining[rm]
+                        remaining += tmp
 
                     if ilvl not in lvl_gaboxes:
                         lvl_gaboxes[ilvl] = {}

pyphare/pyphare/pharesee/hierarchy.py

@@ -9,7 +9,7 @@
 from ..core.box import Box
 from ..core.gridlayout import GridLayout
 import matplotlib.pyplot as plt
-from ..core.phare_utilities import refinement_ratio
+from ..core.phare_utilities import refinement_ratio, deep_copy
 
 
 
@@ -31,6 +31,10 @@ def __init__(self, layout, quantity):
         self.layout   = layout
 
 
+    def __deepcopy__(self, memo):
+        no_copy_keys = ["dataset"] # do not copy these things
+        return deep_copy(self, memo, no_copy_keys)
+
 
 
 class FieldData(PatchData):
@@ -80,8 +84,9 @@ def select(self, box):
 
         overlap = box * gbox
         if overlap is not None:
-            lower = self.layout.AMRIndexToLocal(dim=box.ndim - 1, index=overlap.lower)
-            upper  = self.layout.AMRIndexToLocal(dim=box.ndim - 1, index=overlap.upper)
+            lower = self.layout.AMRToLocal(overlap.lower)
+            upper = self.layout.AMRToLocal(overlap.upper)
+
             if box.ndim == 1:
                 return self.dataset[lower[0] : upper[0] + 1]
             if box.ndim == 2:
@@ -198,6 +203,16 @@ def __repr__(self):
         return self.__str__()
 
 
+    def copy(self):
+        """does not copy patchdatas.datasets (see class PatchData)"""
+        from copy import deepcopy
+        return deepcopy(self)
+
+    def __copy__(self):
+        return self.copy()
+
+
+
 class PatchLevel:
     """is a collection of patches """
 

pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py

@@ -366,14 +366,10 @@ class ParticleLevelGhostGeometryTest(AGeometryTest):
                     Box([10, -1], [10, 3]),
                     Box([10, 14], [10, 23]),
                     Box([10, 34], [10, 40]),
-                    Box([0, -1], [9, -1]),
-                    Box([0, 40], [9, 40]),
                     Box([19, -1], [19, 3]),
                     Box([19, 14], [19, 23]),
                     Box([19, 34], [19, 40]),
                     Box([30, -1], [30, 40]),
-                    Box([20, -1], [29, -1]),
-                    Box([20, 40], [29, 40]),
                     Box([10, 3], [19, 3]),
                     Box([10, 14], [19, 14]),
                     Box([10, 23], [19, 23]),
@@ -400,8 +396,6 @@ def test_level_ghost_boxes(self, refinement_boxes, expected):
                 [actual["boxes"] for actual in lvl_gaboxes[ilvl][key]], []
             )
 
-            self.assertEqual(expected[ilvl], ghost_area_box_list)
-
             fig = hierarchy.plot_2d_patches(
                 ilvl,
                 collections=[
@@ -419,6 +413,9 @@ def test_level_ghost_boxes(self, refinement_boxes, expected):
                 ),
             )
             fig.savefig(f"{type(self).__name__}_lvl_{ilvl}_{self._testMethodName}.png")
+            self.assertEqual(expected[ilvl], ghost_area_box_list)
+
+
 
     @data(
         (

readme.md

@@ -66,7 +66,7 @@ PhD, PostDoc, [contact us](mailto:[email protected]).
 ### Active core team
 
 - [Nicolas Aunai](https://github.com/nicolasaunai)
-- [Philip Deegan](https://github.com/Dekken)
+- [Philip Deegan](https://github.com/PhilipDeegan)
 - [Roch Smets](https://github.com/rochsmets)
 - [Andrea Ciardi](https://sites.google.com/site/andreaciardihomepage/home)
 

src/amr/data/field/field_variable_fill_pattern.h

@@ -29,22 +29,20 @@ namespace PHARE::amr
 // This class is mostly a copy of BoxGeometryVariableFillPattern
 class FieldFillPattern : public SAMRAI::xfer::VariableFillPattern
 {
-protected:
+public:
     FieldFillPattern(std::optional<bool> overwrite_interior)
         : opt_overwrite_interior_{overwrite_interior}
     {
     }
 
-public:
-    template<typename Subclass>
     static auto make_shared(std::shared_ptr<SAMRAI::hier::RefineOperator> const& samrai_op)
     {
         auto const& op = dynamic_cast<AFieldRefineOperator const&>(*samrai_op);
 
         if (op.node_only)
-            return std::make_shared<Subclass>(std::nullopt);
+            return std::make_shared<FieldFillPattern>(std::nullopt);
 
-        return std::make_shared<Subclass>(false);
+        return std::make_shared<FieldFillPattern>(false);
     }
 
 

src/amr/data/field/time_interpolate/field_linear_time_interpolate.h

@@ -21,6 +21,12 @@ using core::dirZ;
 template<typename GridLayoutT, typename FieldT>
 class FieldLinearTimeInterpolate : public SAMRAI::hier::TimeInterpolateOperator
 {
+    static std::size_t constexpr dim = GridLayoutT::dimension;
+    static_assert(dim > 0 && dim <= 3);
+
+    using PhysicalQuantity = decltype(std::declval<FieldT>().physicalQuantity());
+    using FieldDataT       = FieldData<GridLayoutT, FieldT>;
+
 public:
     using GridLayoutImpl = typename GridLayoutT::implT;
 
@@ -30,51 +36,34 @@ class FieldLinearTimeInterpolate : public SAMRAI::hier::TimeInterpolateOperator
     }
 
 
-
-
     virtual ~FieldLinearTimeInterpolate() = default;
 
 
-
-
     void timeInterpolate(SAMRAI::hier::PatchData& destData, SAMRAI::hier::Box const& where,
                          SAMRAI::hier::BoxOverlap const& /*overlap*/,
                          SAMRAI::hier::PatchData const& srcDataOld,
                          SAMRAI::hier::PatchData const& srcDataNew) const override
     {
-        //
-
         auto& fieldDataDest = dynamic_cast<FieldDataT&>(destData);
 
         auto const& fieldDataSrcOld = dynamic_cast<FieldDataT const&>(srcDataOld);
         auto const& fieldDataSrcNew = dynamic_cast<FieldDataT const&>(srcDataNew);
 
         double const interpTime = fieldDataDest.getTime();
-
-        double const oldTime = fieldDataSrcOld.getTime();
-        double const newTime = fieldDataSrcNew.getTime();
-
-
-        double const alpha = (interpTime - oldTime) / (newTime - oldTime);
-
-        auto const& layout = fieldDataDest.gridLayout;
-
-
-
-        auto& fieldDest = fieldDataDest.field;
+        double const oldTime    = fieldDataSrcOld.getTime();
+        double const newTime    = fieldDataSrcNew.getTime();
+        double const alpha      = (interpTime - oldTime) / (newTime - oldTime);
 
         auto const& fieldSrcOld = fieldDataSrcOld.field;
         auto const& fieldSrcNew = fieldDataSrcNew.field;
+        auto& fieldDest         = fieldDataDest.field;
 
-
-
-        auto qty = fieldDest.physicalQuantity();
-
-
+        auto const& layout = fieldDataDest.gridLayout;
         auto const whereLayout
             = FieldGeometry<GridLayoutT, PhysicalQuantity>::layoutFromBox(where, layout);
 
         bool const withGhost{true};
+        auto qty                  = fieldDest.physicalQuantity();
         auto const interpolateBox = FieldGeometry<GridLayoutT, PhysicalQuantity>::toFieldBox(
             where, qty, whereLayout, !withGhost);
 
@@ -86,12 +75,8 @@ class FieldLinearTimeInterpolate : public SAMRAI::hier::TimeInterpolateOperator
         auto srcGhostBox = FieldGeometry<GridLayoutT, PhysicalQuantity>::toFieldBox(
             fieldDataSrcNew.getBox(), qty, fieldDataSrcNew.gridLayout, withGhost);
 
-        auto const localDestBox
-            = AMRToLocal(static_cast<std::add_const_t<decltype(finalBox)>>(finalBox), ghostBox);
-
-        auto const localSrcBox
-            = AMRToLocal(static_cast<std::add_const_t<decltype(finalBox)>>(finalBox), srcGhostBox);
-
+        auto const localDestBox = AMRToLocal(finalBox, ghostBox);
+        auto const localSrcBox  = AMRToLocal(finalBox, srcGhostBox);
 
         if constexpr (dim == 1)
         {
@@ -149,20 +134,7 @@ class FieldLinearTimeInterpolate : public SAMRAI::hier::TimeInterpolateOperator
                 }
             }
         }
-        else
-        {
-            static_assert(dim > 0 && dim <= 3);
-        }
     }
-
-
-
-
-private:
-    static std::size_t constexpr dim = GridLayoutImpl::dimension;
-
-    using PhysicalQuantity = decltype(std::declval<FieldT>().physicalQuantity());
-    using FieldDataT       = FieldData<GridLayoutT, FieldT>;
 };
 
 } // namespace PHARE::amr