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/__init__.py

@@ -64,12 +64,6 @@ def NO_GUI():
     mpl.use("Agg")
 
 
-def getSimulation():
-    from .global_vars import sim
-
-    return sim
-
-
 def _patch_data_ids(restart_file_dir):
     """
     for restarts we save samrai patch data ids to the restart files, which we access from here
@@ -125,10 +119,15 @@ def clearDict():
     pp.stop()
 
 
-def populateDict():
-    from .global_vars import sim as simulation
+def populateDict(simulation=None):
+    # from .global_vars import sim as simulation
     import pybindlibs.dictator as pp
 
+    if simulation is None:
+        from .global_vars import sim
+
+        simulation = sim
+
     # pybind complains if receiving wrong type
     def add_int(path, val):
         pp.add_int(path, int(val))

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")