3d

Author: PhilipDeegan

.github/workflows/cmake_macos.yml

@@ -78,7 +78,7 @@ jobs:
               -DCMAKE_BUILD_TYPE=RelWithDebInfo \
               -DENABLE_SAMRAI_TESTS=OFF -DCMAKE_C_COMPILER_LAUNCHER=ccache \
               -DCMAKE_CXX_COMPILER_LAUNCHER=ccache -DlowResourceTests=ON \
-              -DCMAKE_CXX_FLAGS="-DPHARE_DIAG_DOUBLES=1 "
+              -DCMAKE_CXX_FLAGS="-DPHARE_DIAG_DOUBLES=1 -DPHARE_SIMULATORS=2 "
 
     - name: Build
       working-directory: ${{github.workspace}}/build

pyphare/pyphare/core/box.py

@@ -62,6 +62,9 @@ def __eq__(self, other):
         )
 
     def __sub__(self, other):
+        if isinstance(other, (list, tuple)):
+            assert all([isinstance(item, Box) for item in other])
+            return remove_all(self, other)
         assert isinstance(other, Box)
         return remove(self, other)
 
@@ -179,5 +182,42 @@ def copy(arr, replace):
     return list(boxes.values())
 
 
+def remove_all(box, to_remove):
+    if len(to_remove) > 0:
+        remaining = box - to_remove[0]
+        for to_rm in to_remove[1:]:
+            tmp, remove = [], []
+            for i, rem in enumerate(remaining):
+                if rem * to_rm is not None:
+                    remove.append(i)
+                    tmp += rem - to_rm
+            for rm in reversed(remove):
+                del remaining[rm]
+            remaining += tmp
+        return remaining
+    return box
+
+
+
 def amr_to_local(box, ref_box):
     return Box(box.lower - ref_box.lower, box.upper - ref_box.lower)
+
+
+def select(data, box):
+    return data[tuple([slice(l, u + 1) for l,u in zip(box.lower, box.upper)])]
+
+class DataSelector:
+    """
+        can't assign return to function unless []
+        usage
+        DataSelector(data)[box] = val
+    """
+    def __init__(self, data):
+        self.data = data
+    def __getitem__(self, box_or_slice):
+        if isinstance(box_or_slice, Box):
+            return select(self.data, box_or_slice)
+        return self.data[box_or_slice]
+
+    def __setitem__(self, box_or_slice, val):
+        self.__getitem__(box_or_slice)[:] = val

pyphare/pyphare/core/gridlayout.py

@@ -359,6 +359,24 @@ def yeeCoords(self, knode, iStart, centering, direction, ds, origin, derivOrder)
 
         return x
 
+    def meshCoords(self, qty):
+        ndim = self.ndim
+        assert ndim > 0 and ndim < 4
+        x = self.yeeCoordsFor(qty, "x")
+        if ndim == 1:
+            return x
+        y = self.yeeCoordsFor(qty, "y")
+        if ndim == 2:
+            X, Y = np.meshgrid(x, y, indexing="ij")
+            return np.array([X.flatten(), Y.flatten()]).T.reshape(
+                (len(x), len(y), ndim)
+            )
+        z = self.yeeCoordsFor(qty, "z")
+        X, Y, Z = np.meshgrid(x, y, z, indexing="ij")
+        return np.array([X.flatten(), Y.flatten(), Z.flatten()]).T.reshape(
+            (len(x), len(y), len(z), ndim)
+        )
+
     def yeeCoordsFor(self, qty, direction, withGhosts=True, **kwargs):
         """
         from a qty and a direction, returns a 1d array containing

pyphare/pyphare/cpp/__init__.py

@@ -3,7 +3,7 @@
 #
 
 
-def cpp_lib(override=None):
+def cpp_lib():
     import importlib
 
     return importlib.import_module("pybindlibs.cpp")

pyphare/pyphare/pharein/maxwellian_fluid_model.py

@@ -193,6 +193,10 @@ def validate(self, sim, atol=1e-15):
             self.validate1d(sim, atol)
         elif sim.ndim == 2:
             self.validate2d(sim, atol)
+        elif sim.ndim == 3:
+            self.validate3d(sim, atol)
+        else:
+            raise ValueError("Unknown dimension")
 
     def validate1d(self, sim, atol):
         domain_box = Box([0] * sim.ndim, sim.cells)
@@ -301,3 +305,79 @@ def getCoord(L, R, idir):
             )
             if sim.strict:
                 raise RuntimeError("Simulation is not periodic")
+
+    def validate3d(self, sim, atol):
+        domain_box = Box([0] * sim.ndim, sim.cells)
+        layout = GridLayout(domain_box, domain_box.lower, sim.dl, sim.interp_order)
+        nbrDualGhosts = layout.nbrGhostsPrimal(sim.interp_order)
+        nbrPrimalGhosts = layout.nbrGhostsPrimal(sim.interp_order)
+        directions = ["X", "Y", "Z"]
+        domain = sim.simulation_domain()
+        bx = self.model_dict["bx"]
+        by = self.model_dict["by"]
+        bz = self.model_dict["bz"]
+        is_periodic = True
+        not_periodic = []
+
+        def getCoord(L, R, idir):
+            if idir == 0:
+                return (L, np.zeros_like(L), np.zeros_like(L)), (
+                    R,
+                    np.zeros_like(R),
+                    np.zeros_like(R),
+                )
+            else:
+                return (np.zeros_like(L), np.zeros_like(L), L), (
+                    np.zeros_like(R),
+                    np.zeros_like(R),
+                    R,
+                )
+
+        phare_utilities.debug_print("3d periodic validation")
+        for idir in np.arange(sim.ndim):
+            phare_utilities.debug_print("validating direction ...", idir)
+            if sim.boundary_types[idir] == "periodic":
+                phare_utilities.debug_print(f"direction {idir} is periodic?")
+                dual_left = (np.arange(-nbrDualGhosts, nbrDualGhosts) + 0.5) * sim.dl[0]
+                dual_right = dual_left + domain[0]
+                primal_left = np.arange(-nbrPrimalGhosts, nbrPrimalGhosts) * sim.dl[0]
+                primal_right = primal_left + domain[0]
+
+                direction = directions[idir]
+
+                for b_i, b_name in zip((bx, by, bz), ("Bx", "By", "Bz")):
+                    if layout.qtyIsDual(b_name, direction):
+                        L, R = dual_left, dual_right
+                    else:
+                        L, R = primal_left, primal_right
+
+                coordsL, coordsR = getCoord(L, R, idir)
+                check = np.allclose(b_i(*coordsL), b_i(*coordsR), atol=atol, rtol=0)
+                if not check:
+                    not_periodic += [(b_name, idir)]
+                is_periodic &= check
+
+                for pop in self.populations:
+                    functions = ("vx", "vy", "vz", "vthx", "vthy", "vthz")
+                    L, R = primal_left, primal_right
+                    coordsL, coordsR = getCoord(L, R, idir)
+
+                    for fn in functions:
+                        f = self.model_dict[pop][fn]
+                        fL = f(*coordsL)
+                        fR = f(*coordsR)
+                        check = np.allclose(fL, fR, atol=atol, rtol=0)
+                        phare_utilities.debug_print(
+                            f"checked {fn} : fL = {fL} and fR = {fR} and check = {check}"
+                        )
+                        if not check:
+                            not_periodic += [(fn, idir)]
+                        is_periodic &= check
+
+        if not is_periodic:
+            print(
+                "Warning: Simulation is periodic but some functions are not : ",
+                not_periodic,
+            )
+            if sim.strict:
+                raise RuntimeError("Simulation is not periodic")

pyphare/pyphare/pharesee/geometry.py

@@ -54,19 +54,24 @@ def domain_border_ghost_boxes(domain_box, patches):
     elif domain_box.ndim == 2:
         upper_x, upper_y = domain_box.upper
         return {
-            "bottom": Box(
-                (
-                    0,
-                    0,
-                ),
-                (upper_x, ghost_box_width),
-            ),
-            "top": Box((0, upper_y - ghost_box_width), (upper_x, upper_y)),
             "left": Box((0, 0), (ghost_box_width, upper_y)),
             "right": Box((upper_x - ghost_box_width, 0), (upper_x, upper_y)),
+            "bottom": Box((0, 0), (upper_x, ghost_box_width)),
+            "top": Box((0, upper_y - ghost_box_width), (upper_x, upper_y)),
         }
 
-    raise ValueError("Unhandeled dimension")
+    else:
+        upper_x, upper_y, upper_z = domain_box.upper
+        return {
+            "left": Box((0, 0, 0), (ghost_box_width, upper_y, upper_z)),
+            "right": Box(
+                (upper_x - ghost_box_width, 0, 0), (upper_x, upper_y, upper_z)
+            ),
+            "bottom": Box((0, 0, 0), (upper_x, ghost_box_width, upper_z)),
+            "top": Box((0, upper_y - ghost_box_width, 0), (upper_x, upper_y, upper_z)),
+            "front": Box((0, 0, 0), (upper_x, upper_y, ghost_box_width)),
+            "back": Box((0, 0, upper_z - ghost_box_width), (upper_x, upper_y, upper_z)),
+        }
 
 
 def touch_domain_border(box, domain_box, border):
@@ -79,21 +84,29 @@ def touch_domain_border(box, domain_box, border):
 
 
 def periodicity_shifts(domain_box):
+    shifts = {}
+
     if domain_box.ndim == 1:
         shape_x = domain_box.shape
-        return {
+        shifts = {
             "left": shape_x,
             "right": -shape_x,
         }
+        shifts.update({"leftright": [shifts["left"], shifts["right"]]})
 
     if domain_box.ndim == 2:
         shape_x, shape_y = domain_box.shape
+    if domain_box.ndim == 3:
+        shape_x, shape_y, shape_z = domain_box.shape
+
+    if domain_box.ndim > 1:
         shifts = {
             "left": [(shape_x, 0)],
             "right": [(-shape_x, 0)],
             "bottom": [(0, shape_y)],
             "top": [(0, -shape_y)],
         }
+
         shifts.update(
             {
                 "bottomleft": [*shifts["left"], *shifts["bottom"], (shape_x, shape_y)],
@@ -134,7 +147,7 @@ def periodicity_shifts(domain_box):
                     shifts["topleft"][-1],
                     shifts["topright"][-1],
                 ],
-                "bottomtopleftright": [  # one patch covers domain
+                "leftrightbottomtop": [  # one patch covers domain
                     *shifts["bottomleft"],
                     *shifts["topright"],
                     shifts["bottomright"][-1],
@@ -144,7 +157,35 @@ def periodicity_shifts(domain_box):
         )
 
     if domain_box.ndim == 3:
-        raise ValueError("Unhandeled dimension")
+        front = {
+            f"{k}front": [(v[0], v[1], shape_z) for v in l] for k, l in shifts.items()
+        }
+        back = {
+            f"{k}back": [([v[0], v[1], -shape_z]) for v in l] for k, l in shifts.items()
+        }
+
+        shifts = {k: [([v[0], v[1], 0]) for v in l] for k, l in shifts.items()}
+
+        shifts.update(front)
+        shifts.update(back)
+        shifts.update(
+            {
+                "back": [(0, 0, -shape_z)],
+                "front": [(0, 0, shape_z)],
+                "leftrightbottomtopfrontback": [
+                    *shifts["bottomleftfront"],
+                    *shifts["bottomrightback"],
+                    *shifts["topleftfront"],
+                    *shifts["toprightback"],
+                ],
+            }
+        )
+
+    assert len(list(shifts.keys())) == len(
+        set(["".join(sorted(k)) for k in list(shifts.keys())])
+    )
+
+    shifts = {"".join(sorted(k)): l for k, l in shifts.items()}
 
     return shifts
 
@@ -246,7 +287,7 @@ def borders_per(patch):
     ]
 
     for patch_i, ref_patch in enumerate(border_patches):
-        in_sides = borders_per_patch[ref_patch]
+        in_sides = "".join(sorted(borders_per_patch[ref_patch]))
         assert in_sides in shifts
 
         for ref_pdname, ref_pd in ref_patch.patch_datas.items():
@@ -336,36 +377,41 @@ def get_periodic_list(patches, domain_box, n_ghosts):
             shift_patch(first_patch, domain_box.shape)
             sorted_patches.append(first_patch)
 
+        return sorted_patches
+
+    dbu = domain_box.upper
+
     if dim == 2:
         sides = {
-            "bottom": Box([0, 0], [domain_box.upper[0], 0]),
-            "top": Box(
-                [0, domain_box.upper[1]], [domain_box.upper[0], domain_box.upper[1]]
-            ),
-            "left": Box([0, 0], [0, domain_box.upper[1]]),
-            "right": Box(
-                [domain_box.upper[0], 0], [domain_box.upper[0], domain_box.upper[1]]
-            ),
+            "left": Box([0, 0], [0, dbu[1]]),
+            "right": Box([dbu[0], 0], [dbu[0], dbu[1]]),
+            "bottom": Box([0, 0], [dbu[0], 0]),
+            "top": Box([0, dbu[1]], [dbu[0], dbu[1]]),
         }
 
-        shifts = periodicity_shifts(domain_box)
+    else:
+        sides = {
+            "left": Box([0, 0, 0], [0, dbu[1], dbu[2]]),
+            "right": Box([dbu[0], 0, 0], [dbu[0], dbu[1], dbu[2]]),
+            "bottom": Box([0, 0, 0], [dbu[0], 0, dbu[2]]),
+            "top": Box([0, dbu[1], 0], [dbu[0], dbu[1], dbu[2]]),
+            "front": Box([0, 0, 0], [dbu[0], dbu[1], 0]),
+            "back": Box([0, 0, dbu[2]], [dbu[0], dbu[1], dbu[2]]),
+        }
 
-        def borders_per(box):
-            return "".join(
-                [key for key, side in sides.items() if box * side is not None]
-            )
+    shifts = periodicity_shifts(domain_box)
 
-        for patch in patches:
-            in_sides = borders_per(boxm.grow(patch.box, n_ghosts))
+    def borders_per(box):
+        return "".join([key for key, side in sides.items() if box * side is not None])
 
-            if in_sides in shifts:  # in_sides might be empty, so no borders
-                for shift in shifts[in_sides]:
-                    patch_copy = copy(patch)
-                    shift_patch(patch_copy, shift)
-                    sorted_patches.append(patch_copy)
+    for patch in patches:
+        in_sides = "".join(sorted(borders_per(boxm.grow(patch.box, n_ghosts))))
 
-    if dim == 3:
-        raise ValueError("not yet implemented")
+        if in_sides in shifts:  # in_sides might be empty, so no borders
+            for shift in shifts[in_sides]:
+                patch_copy = copy(patch)
+                shift_patch(patch_copy, shift)
+                sorted_patches.append(patch_copy)
 
     return sorted_patches
 
@@ -470,18 +516,7 @@ def level_ghost_boxes(hierarchy, quantities, levelNbrs=[], time=None):
                     check_patches = patches
 
                 for gabox in ghostAreaBoxes:
-                    remaining = gabox - check_patches[0].box
-
-                    for patch in check_patches[1:]:
-                        tmp = []
-                        remove = []
-                        for i, rem in enumerate(remaining):
-                            if rem * patch.box is not None:
-                                remove.append(i)
-                                tmp += rem - patch.box
-                        for rm in reversed(remove):
-                            del remaining[rm]
-                        remaining += tmp
+                    remaining = gabox - [p.box for p in check_patches]
 
                     if ilvl not in lvl_gaboxes:
                         lvl_gaboxes[ilvl] = {}