diff --git a/.github/workflows/cmake_macos.yml b/.github/workflows/cmake_macos.yml
index 759ac2b8e..dd8141b34 100644
--- a/.github/workflows/cmake_macos.yml
+++ b/.github/workflows/cmake_macos.yml
@@ -87,7 +87,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: ${{runner.workspace}}/build
diff --git a/pyphare/pyphare/core/box.py b/pyphare/pyphare/core/box.py
index 9b59b3a1c..42da5084a 100644
--- a/pyphare/pyphare/core/box.py
+++ b/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
diff --git a/pyphare/pyphare/core/gridlayout.py b/pyphare/pyphare/core/gridlayout.py
index 7cc3b3ace..3a62b50ab 100644
--- a/pyphare/pyphare/core/gridlayout.py
+++ b/pyphare/pyphare/core/gridlayout.py
@@ -366,6 +366,20 @@ 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
diff --git a/pyphare/pyphare/pharesee/geometry.py b/pyphare/pyphare/pharesee/geometry.py
index a7cce9e1e..0d2a9e30b 100644
--- a/pyphare/pyphare/pharesee/geometry.py
+++ b/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
 
@@ -471,18 +517,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] = {}
diff --git a/pyphare/pyphare/pharesee/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy.py
new file mode 100644
index 000000000..794cf3161
--- /dev/null
+++ b/pyphare/pyphare/pharesee/hierarchy.py
@@ -0,0 +1,1834 @@
+import os
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from ..core import box as boxm
+from ..core.box import Box
+from ..core.gridlayout import GridLayout
+from ..core.phare_utilities import deep_copy, refinement_ratio
+from .particles import Particles
+from ..core.phare_utilities import listify
+
+
+class PatchData:
+    """
+    base class for FieldData and ParticleData
+    this class just factors common geometrical properties
+    """
+
+    def __init__(self, layout, quantity):
+        """
+        :param layout: a GridLayout representing the domain on which the data
+             is defined
+        :param quantity: ['field', 'particle']
+        """
+        self.quantity = quantity
+        self.box = layout.box
+        self.origin = layout.origin
+        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):
+    """
+    Concrete type of PatchData representing a physical quantity
+    defined on a grid.
+    """
+
+    @property
+    def x(self):
+        withGhosts = self.field_name != "tags"
+        if self._x is None:
+            self._x = self.layout.yeeCoordsFor(
+                self.field_name,
+                "x",
+                withGhosts=withGhosts,
+                centering=self.centerings[0],
+            )
+        return self._x
+
+    @property
+    def y(self):
+        withGhosts = self.field_name != "tags"
+        if self._y is None:
+            self._y = self.layout.yeeCoordsFor(
+                self.field_name,
+                "y",
+                withGhosts=withGhosts,
+                centering=self.centerings[1],
+            )
+        return self._y
+
+    @property
+    def z(self):
+        withGhosts = self.field_name != "tags"
+        if self._z is None:
+            self._z = self.layout.yeeCoordsFor(
+                self.field_name,
+                "z",
+                withGhosts=withGhosts,
+                centering=self.centerings[2],
+            )
+        return self._z
+
+    def primal_directions(self):
+        return self.size - self.ghost_box.shape
+
+    def __str__(self):
+        return "FieldData: (box=({}, {}), key={})".format(
+            self.layout.box, self.layout.box.shape, self.field_name
+        )
+
+    def __repr__(self):
+        return self.__str__()
+
+    def select(self, box):
+        """
+        return view of internal data based on overlap of input box
+           returns a view +1 in size in primal directions
+        """
+        assert isinstance(box, Box) and box.ndim == self.box.ndim
+
+        gbox = self.ghost_box.copy()
+        gbox.upper += self.primal_directions()
+
+        box = box.copy()
+        box.upper += self.primal_directions()
+
+        overlap = box * gbox
+        if overlap is not None:
+            return boxm.select(self.dataset, self.layout.AMRBoxToLocal(overlap))
+        return np.array([])
+
+    def __getitem__(self, box):
+        return self.select(box)
+
+    def __init__(self, layout, field_name, data, **kwargs):
+        """
+        :param layout: A GridLayout representing the domain on which data is defined
+        :param field_name: the name of the field (e.g. "Bx")
+        :param data: the dataset from which data can be accessed
+        """
+        super().__init__(layout, "field")
+        self._x = None
+        self._y = None
+        self._z = None
+
+        self.layout = layout
+        self.field_name = field_name
+        self.name = field_name
+        self.dl = np.asarray(layout.dl)
+        self.ndim = layout.box.ndim
+        self.ghosts_nbr = np.zeros(self.ndim, dtype=int)
+
+        if field_name in layout.centering["X"]:
+            directions = ["X", "Y", "Z"][: layout.box.ndim]  # drop unused directions
+            self.centerings = [
+                layout.qtyCentering(field_name, direction) for direction in directions
+            ]
+        elif "centering" in kwargs:
+            if isinstance(kwargs["centering"], list):
+                self.centerings = kwargs["centering"]
+                assert len(self.centerings) == self.ndim
+            else:
+                if self.ndim != 1:
+                    raise ValueError(
+                        "FieldData invalid dimenion for centering argument, expected list for dim > 1"
+                    )
+                self.centerings = [kwargs["centering"]]
+        else:
+            raise ValueError(
+                f"centering not specified and cannot be inferred from field name : {field_name}"
+            )
+
+        if self.field_name != "tags":
+            for i, centering in enumerate(self.centerings):
+                self.ghosts_nbr[i] = layout.nbrGhosts(layout.interp_order, centering)
+
+        self.ghost_box = boxm.grow(layout.box, self.ghosts_nbr)
+
+        self.size = np.copy(self.ghost_box.shape)
+        self.offset = np.zeros(self.ndim)
+
+        for i, centering in enumerate(self.centerings):
+            if centering == "primal":
+                self.size[i] = self.ghost_box.shape[i] + 1
+            else:
+                self.size[i] = self.ghost_box.shape[i]
+                self.offset[i] = 0.5 * self.dl[i]
+
+        self.dataset = data
+
+    def meshgrid(self, select=None):
+        def grid():
+            if self.ndim == 1:
+                return [self.x]
+            if self.ndim == 2:
+                return np.meshgrid(self.x, self.y, indexing="ij")
+            return np.meshgrid(self.x, self.y, self.z, indexing="ij")
+
+        mesh = grid()
+        if select is not None:
+            return tuple(g[select] for g in mesh)
+        return mesh
+
+
+class ParticleData(PatchData):
+    """
+    Concrete type of PatchData representing particles in a region
+    """
+
+    def __init__(self, layout, data, pop_name):
+        """
+        :param layout: A GridLayout object representing the domain in which particles are
+        :param data: dataset containing particles
+        """
+        super().__init__(layout, "particles")
+        self.dataset = data
+        self.pop_name = pop_name
+        self.name = pop_name
+        self.ndim = layout.box.ndim
+
+        self.pop_name = pop_name
+        if layout.interp_order == 1:
+            self.ghosts_nbr = np.array([1] * layout.box.ndim)
+        elif layout.interp_order == 2 or layout.interp_order == 3:
+            self.ghosts_nbr = np.array([2] * layout.box.ndim)
+        else:
+            raise RuntimeError(
+                "invalid interpolation order {}".format(layout.interp_order)
+            )
+
+        self.ghost_box = boxm.grow(layout.box, self.ghosts_nbr)
+        assert (self.box.lower == self.ghost_box.lower + self.ghosts_nbr).all()
+
+    def select(self, box):
+        return self.dataset[box]
+
+    def __getitem__(self, box):
+        return self.select(box)
+
+    def size(self):
+        return self.dataset.size()
+
+
+class Patch:
+    """
+    A patch represents a hyper-rectangular region of space
+    """
+
+    def __init__(self, patch_datas, patch_id=""):
+        """
+        :param patch_datas: a list of PatchData objects
+        these are assumed to "belong" to the Patch so to
+        share the same origin, mesh size and box.
+        """
+        pdata0 = list(patch_datas.values())[0]  # 0 represents all others
+        self.layout = pdata0.layout
+        self.box = pdata0.layout.box
+        self.origin = pdata0.layout.origin
+        self.dl = pdata0.layout.dl
+        self.patch_datas = patch_datas
+        self.id = patch_id
+
+    def __str__(self):
+        return f"Patch: box( {self.box}), id({self.id})"
+
+    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()
+
+    def __call__(self, qty, **kwargs):
+        # take slice/slab of 1/2d array from 2/3d array
+        if "x" in kwargs and len(kwargs) == 1:
+            cut = kwargs["x"]
+            idim = 0
+        elif "y" in kwargs and len(kwargs) == 1:
+            cut = kwargs["y"]
+            idim = 1
+        else:
+            raise ValueError("need to specify either x or y cut coordinate")
+        pd = self.patch_datas[qty]
+        origin = pd.origin[idim]
+        idx = int((cut - origin) / pd.layout.dl[idim])
+        nbrGhosts = pd.ghosts_nbr[idim]
+        if idim == 0:
+            return pd.dataset[idx + nbrGhosts, nbrGhosts:-nbrGhosts]
+        elif idim == 1:
+            return pd.dataset[nbrGhosts:-nbrGhosts, idx + nbrGhosts]
+
+
+class PatchLevel:
+    """is a collection of patches"""
+
+    def __init__(self, lvl_nbr, patches):
+        self.level_number = lvl_nbr
+        self.patches = patches
+
+    def __iter__(self):
+        return self.patches.__iter__()
+
+    def level_range(self):
+        name = list(self.patches[0].patch_datas.keys())[0]
+        return min([patch.patch_datas[name].x.min() for patch in self.patches]), max(
+            [patch.patch_datas[name].x.max() for patch in self.patches]
+        )
+
+
+def are_adjacent(lower, upper, atol=1e-6):
+    return np.abs(upper[0] - lower[-1]) < atol
+
+
+def overlap_mask_1d(x, dl, level, qty):
+    """
+    return the mask for x where x is overlaped by the qty patch datas
+    on the given level, assuming that this level is finer than the one of x
+
+    :param x: 1d array containing the [x] position
+    :param dl: list containing the grid steps where x is defined
+    :param level: a given level associated to a hierarchy
+    :param qty: ['Bx', 'By', 'Bz', 'Ex', 'Ey', 'Ez', 'Fx', 'Fy', 'Fz', 'Vx', 'Vy', 'Vz', 'rho']
+    """
+
+    is_overlaped = np.ones(x.shape[0], dtype=bool) * False
+
+    for patch in level.patches:
+        pdata = patch.patch_datas[qty]
+        ghosts_nbr = pdata.ghosts_nbr
+
+        fine_x = pdata.x[ghosts_nbr[0] - 1 : -ghosts_nbr[0] + 1]
+
+        fine_dl = pdata.dl
+        local_dl = dl
+
+        if fine_dl[0] < local_dl[0]:
+            xmin, xmax = fine_x.min(), fine_x.max()
+
+            overlaped_idx = np.where((x > xmin) & (x < xmax))[0]
+
+            is_overlaped[overlaped_idx] = True
+
+        else:
+            raise ValueError("level needs to have finer grid resolution than that of x")
+
+    return is_overlaped
+
+
+def overlap_mask_2d(x, y, dl, level, qty):
+    """
+    return the mask for x & y where ix & y are overlaped by the qty patch datas
+    on the given level, assuming that this level is finer than the one of x & y
+    important note : this mask is flatten
+
+    :param x: 1d array containing the [x] position
+    :param y: 1d array containing the [y] position
+    :param dl: list containing the grid steps where x and y are defined
+    :param level: a given level associated to a hierarchy
+    :param qty: ['Bx', 'By', 'Bz', 'Ex', 'Ey', 'Ez', 'Fx', 'Fy', 'Fz', 'Vx', 'Vy', 'Vz', 'rho']
+    """
+
+    is_overlaped = np.ones([x.shape[0] * y.shape[0]], dtype=bool) * False
+
+    for patch in level.patches:
+        pdata = patch.patch_datas[qty]
+        ghosts_nbr = pdata.ghosts_nbr
+
+        fine_x = pdata.x[ghosts_nbr[0] - 1 : -ghosts_nbr[0] + 1]
+        fine_y = pdata.y[ghosts_nbr[1] - 1 : -ghosts_nbr[1] + 1]
+
+        fine_dl = pdata.dl
+        local_dl = dl
+
+        if (fine_dl[0] < local_dl[0]) and (fine_dl[1] < local_dl[1]):
+            xmin, xmax = fine_x.min(), fine_x.max()
+            ymin, ymax = fine_y.min(), fine_y.max()
+
+            xv, yv = np.meshgrid(x, y, indexing="ij")
+            xf = xv.flatten()
+            yf = yv.flatten()
+
+            overlaped_idx = np.where(
+                (xf > xmin) & (xf < xmax) & (yf > ymin) & (yf < ymax)
+            )[0]
+
+            is_overlaped[overlaped_idx] = True
+
+        else:
+            raise ValueError(
+                "level needs to have finer grid resolution than that of x or y"
+            )
+
+    return is_overlaped
+
+
+def flat_finest_field(hierarchy, qty, time=None):
+    """
+    returns 2 flattened arrays containing the data (with shape [Npoints])
+    and the coordinates (with shape [Npoints, Ndim]) for the given
+    hierarchy of qty.
+
+    :param hierarchy: the hierarchy where qty is defined
+    :param qty: the field (eg "Bx") that we want
+    """
+
+    dim = hierarchy.ndim
+
+    if dim == 1:
+        return flat_finest_field_1d(hierarchy, qty, time)
+    elif dim == 2:
+        return flat_finest_field_2d(hierarchy, qty, time)
+    elif dim == 3:
+        raise RuntimeError("Not implemented")
+        # return flat_finest_field_3d(hierarchy, qty, time)
+    else:
+        raise ValueError("the dim of a hierarchy should be 1, 2 or 3")
+
+
+def flat_finest_field_1d(hierarchy, qty, time=None):
+    lvl = hierarchy.levels(time)
+
+    for ilvl in range(hierarchy.finest_level(time) + 1)[::-1]:
+        patches = lvl[ilvl].patches
+
+        for ip, patch in enumerate(patches):
+            pdata = patch.patch_datas[qty]
+
+            # all but 1 ghost nodes are removed in order to limit
+            # the overlapping, but to keep enough point to avoid
+            # any extrapolation for the interpolator
+            needed_points = pdata.ghosts_nbr - 1
+
+            # data = pdata.dataset[patch.box] # TODO : once PR 551 will be merged...
+            data = pdata.dataset[needed_points[0] : -needed_points[0]]
+            x = pdata.x[needed_points[0] : -needed_points[0]]
+
+            if ilvl == hierarchy.finest_level(time):
+                if ip == 0:
+                    final_data = data
+                    final_x = x
+                else:
+                    final_data = np.concatenate((final_data, data))
+                    final_x = np.concatenate((final_x, x))
+
+            else:
+                is_overlaped = overlap_mask_1d(
+                    x, pdata.dl, hierarchy.level(ilvl + 1, time), qty
+                )
+
+                finest_data = data[~is_overlaped]
+                finest_x = x[~is_overlaped]
+
+                final_data = np.concatenate((final_data, finest_data))
+                final_x = np.concatenate((final_x, finest_x))
+
+    return final_data, final_x
+
+
+def flat_finest_field_2d(hierarchy, qty, time=None):
+    lvl = hierarchy.levels(time)
+
+    for ilvl in range(hierarchy.finest_level(time) + 1)[::-1]:
+        patches = lvl[ilvl].patches
+
+        for ip, patch in enumerate(patches):
+            pdata = patch.patch_datas[qty]
+
+            # all but 1 ghost nodes are removed in order to limit
+            # the overlapping, but to keep enough point to avoid
+            # any extrapolation for the interpolator
+            needed_points = pdata.ghosts_nbr - 1
+
+            # data = pdata.dataset[patch.box] # TODO : once PR 551 will be merged...
+            data = pdata.dataset[
+                needed_points[0] : -needed_points[0],
+                needed_points[1] : -needed_points[1],
+            ]
+            x = pdata.x[needed_points[0] : -needed_points[0]]
+            y = pdata.y[needed_points[1] : -needed_points[1]]
+
+            xv, yv = np.meshgrid(x, y, indexing="ij")
+
+            data_f = data.flatten()
+            xv_f = xv.flatten()
+            yv_f = yv.flatten()
+
+            if ilvl == hierarchy.finest_level(time):
+                if ip == 0:
+                    final_data = data_f
+                    tmp_x = xv_f
+                    tmp_y = yv_f
+                else:
+                    final_data = np.concatenate((final_data, data_f))
+                    tmp_x = np.concatenate((tmp_x, xv_f))
+                    tmp_y = np.concatenate((tmp_y, yv_f))
+
+            else:
+                is_overlaped = overlap_mask_2d(
+                    x, y, pdata.dl, hierarchy.level(ilvl + 1, time), qty
+                )
+
+                finest_data = data_f[~is_overlaped]
+                finest_x = xv_f[~is_overlaped]
+                finest_y = yv_f[~is_overlaped]
+
+                final_data = np.concatenate((final_data, finest_data))
+                tmp_x = np.concatenate((tmp_x, finest_x))
+                tmp_y = np.concatenate((tmp_y, finest_y))
+
+    final_xy = np.stack((tmp_x, tmp_y), axis=1)
+
+    return final_data, final_xy
+
+
+def finest_part_data(hierarchy, time=None):
+    """
+    returns a dict {popname : Particles}
+    Particles contained in the dict are those from
+    the finest patches available at a given location
+    """
+    from copy import deepcopy
+
+    from .particles import remove
+
+    # we are going to return a dict {popname : Particles}
+    # we prepare it with population names
+    aPatch = hierarchy.level(0, time=time).patches[0]
+    particles = {popname: None for popname in aPatch.patch_datas.keys()}
+
+    # our strategy is to explore the hierarchy from the finest
+    # level to the coarsest. at Each level we keep only particles
+    # that are in cells that are not overlaped by finer boxes
+
+    # this dict keeps boxes for patches at each level
+    # each level will thus need this dict to see next finer boxes
+    lvlPatchBoxes = {ilvl: [] for ilvl in range(hierarchy.finest_level(time) + 1)}
+
+    for ilvl in range(hierarchy.finest_level(time) + 1)[::-1]:
+        plvl = hierarchy.level(ilvl, time=time)
+        for ip, patch in enumerate(plvl.patches):
+            lvlPatchBoxes[ilvl].append(patch.box)
+            for popname, pdata in patch.patch_datas.items():
+                # if we're at the finest level
+                # we need to keep all particles
+                if ilvl == hierarchy.finest_level(time):
+                    if particles[popname] is None:
+                        particles[popname] = deepcopy(pdata.dataset)
+                    else:
+                        particles[popname].add(deepcopy(pdata.dataset))
+
+                # if there is a finer level
+                # we need to keep only those of the current patch
+                # that are not in cells overlaped by finer patch boxes
+                else:
+                    icells = pdata.dataset.iCells
+                    parts = deepcopy(pdata.dataset)
+                    create = True
+                    for finerBox in lvlPatchBoxes[ilvl + 1]:
+                        coarseFinerBox = boxm.coarsen(finerBox, refinement_ratio)
+                        within = np.where(
+                            (icells >= coarseFinerBox.lower[0])
+                            & (icells <= coarseFinerBox.upper[0])
+                        )[0]
+                        if create:
+                            toRemove = within
+                            create = False
+                        else:
+                            toRemove = np.concatenate((toRemove, within))
+
+                    if toRemove.size != 0:
+                        parts = remove(parts, toRemove)
+                    if parts is not None:
+                        particles[popname].add(parts)
+    return particles
+
+
+class PatchHierarchy(object):
+    """is a collection of patch levels"""
+
+    def __init__(
+        self,
+        patch_levels,
+        domain_box,
+        refinement_ratio=2,
+        time=0.0,
+        data_files=None,
+        **kwargs,
+    ):
+        self.patch_levels = patch_levels
+        self.ndim = len(domain_box.lower)
+        self.time_hier = {}
+        self.time_hier.update({self.format_timestamp(time): patch_levels})
+
+        self.domain_box = domain_box
+        self.refinement_ratio = refinement_ratio
+
+        self.data_files = {}
+        self._sim = None
+
+        if data_files is not None:
+            self.data_files.update(data_files)
+
+        if len(self.quantities()) > 1:
+            for qty in self.quantities():
+                if qty in self.__dict__:
+                    continue
+                first = True
+                for time, levels in self.time_hier.items():
+                    new_lvls = {}
+                    for ilvl, level in levels.items():
+                        patches = []
+                        for patch in level.patches:
+                            patches += [Patch({qty: patch.patch_datas[qty]}, patch.id)]
+                        new_lvls[ilvl] = PatchLevel(ilvl, patches)
+                    if first:
+                        self.__dict__[qty] = PatchHierarchy(
+                            new_lvls, self.domain_box, time=time
+                        )
+                        first = False
+                    else:
+                        self.qty.time_hier[time] = new_lvls  # pylint: disable=E1101
+
+    def __getitem__(self, qty):
+        return self.__dict__[qty]
+
+    @property
+    def sim(self):
+        if self._sim:
+            return self._sim
+
+        if "py_attrs" not in self.data_files:
+            raise ValueError("Simulation is not available for deserialization")
+
+        from ..pharein.simulation import deserialize
+
+        try:
+            self._sim = deserialize(
+                self.data_files["py_attrs"].attrs["serialized_simulation"]
+            )
+        except Exception as e:
+            raise RuntimeError(f"Failed to deserialize simulation from data file : {e}")
+        return self._sim
+
+    def __call__(self, qty=None, **kwargs):
+        # take slice/slab of 1/2d array from 2/3d array
+        def cuts(c, coord):
+            return c > coord.min() and c < coord.max()
+
+        class Extractor:
+            def __init__(self):
+                self.exclusions = []
+
+            def extract(self, coord, data):
+                mask = coord == coord
+                for exclusion in self.exclusions:
+                    idx = np.where(
+                        (coord > exclusion[0] - 1e-6) & (coord < exclusion[1] + 1e-6)
+                    )[0]
+                    mask[idx] = False
+
+                self.exclusions += [(coord.min(), coord.max())]
+                return coord[mask], data[mask]
+
+        def domain_coords(patch, qty):
+            pd = patch.patch_datas[qty]
+            nbrGhosts = pd.ghosts_nbr[0]
+            return pd.x[nbrGhosts:-nbrGhosts], pd.y[nbrGhosts:-nbrGhosts]
+
+        if len(kwargs) < 1 or len(kwargs) > 3:
+            raise ValueError("Error - must provide coordinates")
+        if qty == None:
+            if len(self.quantities()) == 1:
+                qty = self.quantities()[0]
+            else:
+                raise ValueError(
+                    "The PatchHierarchy has several quantities but none is specified"
+                )
+
+        if len(kwargs) == 1:  # 1D cut
+            if "x" in kwargs:
+                c = kwargs["x"]
+                slice_dim = 1
+                cst_dim = 0
+            else:
+                c = kwargs["y"]
+                slice_dim = 0
+                cst_dim = 1
+
+        extractor = Extractor()
+        datas = []
+        coords = []
+        ilvls = list(self.levels().keys())[::-1]
+
+        for ilvl in ilvls:
+            lvl = self.patch_levels[ilvl]
+            for patch in lvl.patches:
+                slice_coord = domain_coords(patch, qty)[slice_dim]
+                cst_coord = domain_coords(patch, qty)[cst_dim]
+
+                if cuts(c, cst_coord):
+                    data = patch(qty, **kwargs)
+                    coord_keep, data_keep = extractor.extract(slice_coord, data)
+                    datas += [data_keep]
+                    coords += [coord_keep]
+
+        cut = np.concatenate(datas)
+        coords = np.concatenate(coords)
+        ic = np.argsort(coords)
+        coords = coords[ic]
+        cut = cut[ic]
+        return coords, cut
+
+    def _default_time(self):
+        return self.times()[0]
+
+    def finest_level(self, time=None):
+        if time is None:
+            time = self._default_time()
+        return max(list(self.levels(time=time).keys()))
+
+    def levels(self, time=None):
+        if time is None:
+            time = self._default_time()
+        return self.time_hier[self.format_timestamp(time)]
+
+    def level(self, level_number, time=None):
+        return self.levels(time)[level_number]
+
+    def levelNbr(self, time=None):
+        if time is None:
+            time = self._default_time()
+        return len(self.levels(time).items())
+
+    def levelNbrs(self, time=None):
+        if time is None:
+            time = self._default_time()
+        return list(self.levels(time).keys())
+
+    def is_homogeneous(self):
+        """
+        return True if all patches of all levels at all times
+        have the same patch data quantities
+        """
+        qties = self._quantities()
+        it_is = True
+        for time, levels in self.time_hier.items():
+            for ilvl, lvl in levels.items():
+                for patch in lvl.patches:
+                    it_is &= qties == list(patch.patch_datas.keys())
+        return it_is
+
+    def _quantities(self):
+        for ilvl, lvl in self.levels().items():
+            if len(lvl.patches) > 0:
+                return list(self.level(ilvl).patches[0].patch_datas.keys())
+        return []
+
+    def quantities(self):
+        if not self.is_homogeneous():
+            raise RuntimeError("Error - hierarchy is not homogeneous")
+        return self._quantities()
+
+    def global_min(self, qty, **kwargs):
+        time = kwargs.get("time", self._default_time())
+        first = True
+        for ilvl, lvl in self.levels(time).items():
+            for patch in lvl.patches:
+                pd = patch.patch_datas[qty]
+                if first:
+                    m = pd.dataset[:].min()
+                    first = False
+                else:
+                    m = min(m, pd.dataset[:].min())
+
+        return m
+
+    def global_max(self, qty, **kwargs):
+        time = kwargs.get("time", self._default_time())
+        first = True
+        for ilvl, lvl in self.levels(time).items():
+            for patch in lvl.patches:
+                pd = patch.patch_datas[qty]
+                if first:
+                    m = pd.dataset[:].max()
+                    first = False
+                else:
+                    m = max(m, pd.dataset[:].max())
+
+        return m
+
+    def refined_domain_box(self, level_number):
+        """
+        returns the domain box refined for a given level number
+        """
+        assert level_number >= 0
+        return boxm.refine(self.domain_box, self.refinement_ratio**level_number)
+
+    def format_timestamp(self, timestamp):
+        if isinstance(timestamp, str):
+            return timestamp
+        return "{:.10f}".format(timestamp)
+
+    def level_domain_box(self, level_number):
+        if level_number == 0:
+            return self.domain_box
+        return self.refined_domain_box(level_number)
+
+    def __str__(self):
+        s = "Hierarchy: \n"
+        for t, patch_levels in self.time_hier.items():
+            for ilvl, lvl in patch_levels.items():
+                s = s + "Level {}\n".format(ilvl)
+                for ip, patch in enumerate(lvl.patches):
+                    for qty_name, pd in patch.patch_datas.items():
+                        pdstr = "    P{ip} {type} {pdname} box is {box} and ghost box is {gbox}"
+                        s = s + pdstr.format(
+                            ip=ip,
+                            type=type(pd.dataset),
+                            pdname=qty_name,
+                            box=patch.box,
+                            gbox=pd.ghost_box,
+                        )
+                        s = s + "\n"
+        return s
+
+    def times(self):
+        return np.sort(np.asarray(list(self.time_hier.keys())))
+
+    def plot_patches(self, save=False):
+        fig, ax = plt.subplots(figsize=(10, 3))
+        for ilvl, lvl in self.levels(0.0).items():
+            lvl_offset = ilvl * 0.1
+            for patch in lvl.patches:
+                dx = patch.dl[0]
+                x0 = patch.box.lower * dx
+                x1 = patch.box.upper * dx
+                xcells = np.arange(x0, x1 + dx, dx)
+                y = lvl_offset + np.zeros_like(xcells)
+                ax.plot(xcells, y, marker=".")
+
+        if save:
+            fig.savefig("hierarchy.png")
+
+    def box_to_Rectangle(self, box):
+        from matplotlib.patches import Rectangle
+
+        return Rectangle(box.lower, *box.shape)
+
+    def plot_2d_patches(self, ilvl, collections, **kwargs):
+        if isinstance(collections, list) and all(
+            [isinstance(el, Box) for el in collections]
+        ):
+            collections = [{"boxes": collections}]
+
+        from matplotlib.collections import PatchCollection
+
+        level_domain_box = self.level_domain_box(ilvl)
+        mi, ma = level_domain_box.lower.min(), level_domain_box.upper.max()
+
+        fig, ax = kwargs.get("subplot", plt.subplots(figsize=(6, 6)))
+
+        for collection in collections:
+            facecolor = collection.get("facecolor", "none")
+            edgecolor = collection.get("edgecolor", "purple")
+            alpha = collection.get("alpha", 1)
+            rects = [self.box_to_Rectangle(box) for box in collection["boxes"]]
+
+            ax.add_collection(
+                PatchCollection(
+                    rects, facecolor=facecolor, alpha=alpha, edgecolor=edgecolor
+                )
+            )
+
+        if "title" in kwargs:
+            from textwrap import wrap
+
+            xfigsize = int(fig.get_size_inches()[0] * 10)  # 10 characters per inch
+            ax.set_title("\n".join(wrap(kwargs["title"], xfigsize)))
+
+        major_ticks = np.arange(mi - 5, ma + 5 + 5, 5)
+        ax.set_xticks(major_ticks)
+        ax.set_yticks(major_ticks)
+
+        minor_ticks = np.arange(mi - 5, ma + 5 + 5, 1)
+        ax.set_xticks(minor_ticks, minor=True)
+        ax.set_yticks(minor_ticks, minor=True)
+
+        ax.grid(which="both")
+
+        return fig
+
+    def plot1d(self, **kwargs):
+        """
+        plot
+        """
+        usr_lvls = kwargs.get("levels", (0,))
+        qty = kwargs.get("qty", None)
+        time = kwargs.get("time", self.times()[0])
+
+        if "ax" not in kwargs:
+            fig, ax = plt.subplots()
+        else:
+            ax = kwargs["ax"]
+            fig = ax.figure
+        for lvl_nbr, level in self.levels(time).items():
+            if lvl_nbr not in usr_lvls:
+                continue
+            for ip, patch in enumerate(level.patches):
+                pdata_nbr = len(patch.patch_datas)
+                pdata_names = list(patch.patch_datas.keys())
+                if qty is None and pdata_nbr != 1:
+                    multiple = "multiple quantities in patch, "
+                    err = (
+                        multiple
+                        + "please specify a quantity in "
+                        + " ".join(pdata_names)
+                    )
+                    raise ValueError(err)
+                if qty is None:
+                    qty = pdata_names[0]
+
+                layout = patch.patch_datas[qty].layout
+                nbrGhosts = layout.nbrGhostFor(qty)
+                val = patch.patch_datas[qty][patch.box]
+                x = patch.patch_datas[qty].x[nbrGhosts[0] : -nbrGhosts[0]]
+                label = "L{level}P{patch}".format(level=lvl_nbr, patch=ip)
+                marker = kwargs.get("marker", "")
+                ls = kwargs.get("ls", "--")
+                color = kwargs.get("color", "k")
+                ax.plot(x, val, label=label, marker=marker, ls=ls, color=color)
+
+        ax.set_title(kwargs.get("title", ""))
+        ax.set_xlabel(kwargs.get("xlabel", "x"))
+        ax.set_ylabel(kwargs.get("ylabel", qty))
+        if "xlim" in kwargs:
+            ax.set_xlim(kwargs["xlim"])
+        if "ylim" in kwargs:
+            ax.set_ylim(kwargs["ylim"])
+
+        if kwargs.get("legend", None) is not None:
+            ax.legend()
+
+        if "filename" in kwargs:
+            fig.savefig(kwargs["filename"])
+
+    def plot2d(self, **kwargs):
+        from matplotlib.patches import Rectangle
+        from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+        time = kwargs.get("time", self._default_time())
+        usr_lvls = kwargs.get("levels", self.levelNbrs(time))
+        default_qty = None
+        if len(self.quantities()) == 1:
+            default_qty = self.quantities()[0]
+        qty = kwargs.get("qty", default_qty)
+
+        if "ax" not in kwargs:
+            fig, ax = plt.subplots()
+        else:
+            ax = kwargs["ax"]
+            fig = ax.figure
+
+        glob_min = self.global_min(qty)
+        glob_max = self.global_max(qty)
+        # assumes max 5 levels...
+        patchcolors = {ilvl: "k" for ilvl in usr_lvls}
+        patchcolors = kwargs.get("patchcolors", patchcolors)
+        if not isinstance(patchcolors, dict):
+            patchcolors = dict(zip(usr_lvls, patchcolors))
+
+        linewidths = {ilvl: 1 for ilvl in usr_lvls}
+        linewidths = kwargs.get("lw", linewidths)
+        if not isinstance(linewidths, dict):
+            linewidths = dict(zip(usr_lvls, linewidths))
+        linestyles = {ilvl: "-" for ilvl in usr_lvls}
+        linestyles = kwargs.get("ls", linestyles)
+        if not isinstance(linestyles, dict):
+            linestyles = dict(zip(usr_lvls, linestyles))
+
+        for lvl_nbr, lvl in self.levels(time).items():
+            if lvl_nbr not in usr_lvls:
+                continue
+            for patch in self.level(lvl_nbr, time).patches:
+                pdat = patch.patch_datas[qty]
+                data = pdat.dataset[:]
+                nbrGhosts = pdat.ghosts_nbr
+                x = pdat.x
+                y = pdat.y
+
+                # if nbrGhosts is 0, we cannot do array[0,-0]
+                if np.all(nbrGhosts == np.zeros_like(nbrGhosts)):
+                    x = np.copy(x)
+                    y = np.copy(y)
+                else:
+                    data = pdat[patch.box]
+                    x = np.copy(x[nbrGhosts[0] : -nbrGhosts[0]])
+                    y = np.copy(y[nbrGhosts[1] : -nbrGhosts[1]])
+                dx, dy = pdat.layout.dl
+                x -= dx * 0.5
+                y -= dy * 0.5
+                x = np.append(x, x[-1] + dx)
+                y = np.append(y, y[-1] + dy)
+                im = ax.pcolormesh(
+                    x,
+                    y,
+                    data.T,
+                    cmap=kwargs.get("cmap", "Spectral_r"),
+                    vmin=kwargs.get("vmin", glob_min - 1e-6),
+                    vmax=kwargs.get("vmax", glob_max + 1e-6),
+                )
+
+                if kwargs.get("plot_patches", False) is True:
+                    r = Rectangle(
+                        (patch.box.lower[0] * dx, patch.box.lower[1] * dy),
+                        patch.box.shape[0] * dx,
+                        patch.box.shape[1] * dy,
+                        fc="none",
+                        ec=patchcolors[lvl_nbr],
+                        alpha=0.4,
+                        lw=linewidths[lvl_nbr],
+                        ls=linestyles[lvl_nbr],
+                    )
+                    ax.add_patch(r)
+
+        ax.set_aspect(kwargs.get("aspect", "equal"))
+        ax.set_title(kwargs.get("title", ""))
+        ax.set_xlabel(kwargs.get("xlabel", "x"))
+        ax.set_ylabel(kwargs.get("ylabel", "y"))
+        if "xlim" in kwargs:
+            ax.set_xlim(kwargs["xlim"])
+        if "ylim" in kwargs:
+            ax.set_ylim(kwargs["ylim"])
+
+        divider = make_axes_locatable(ax)
+        cax = divider.append_axes("right", size="5%", pad=0.08)
+        plt.colorbar(im, ax=ax, cax=cax)
+
+        if kwargs.get("legend", None) is not None:
+            ax.legend()
+
+        if "filename" in kwargs:
+            fig.savefig(kwargs["filename"], dpi=kwargs.get("dpi", 200))
+
+        return fig, ax
+
+    def plot3d(self, **kwargs):
+        """!HAX!"""
+        time = kwargs.get("time", self._default_time())
+        usr_lvls = kwargs.get("levels", self.levelNbrs(time))
+        default_qty = None
+        if len(self.quantities()) == 1:
+            default_qty = self.quantities()[0]
+        qty = kwargs.get("qty", default_qty)
+        for lvl_nbr, lvl in self.levels(time).items():
+            if lvl_nbr not in usr_lvls:
+                continue
+            for patch in self.level(lvl_nbr, time).patches:
+                pdat = patch.patch_datas[qty]
+                primals = pdat.primal_directions()
+                if primals[0]:
+                    pdat._x = pdat.x[:-1]
+                if primals[1]:
+                    pdat._y = pdat.y[:-1]
+                pdat.dataset = pdat.dataset[:, :, int(pdat.ghost_box.shape[2] / 2)]
+                patch.box.lower = patch.box.lower[:-1]
+                patch.box.upper = patch.box.upper[:-1]
+                patch.box.ndim = 2
+
+                pdat.ghost_box.lower = pdat.ghost_box.lower[:-1]
+                pdat.ghost_box.upper = pdat.ghost_box.upper[:-1]
+                pdat.ghost_box.ndim = 2
+                pdat.size = np.copy(pdat.ghost_box.shape)
+                pdat.layout.dl = pdat.layout.dl[:-1]
+
+        return self.plot2d(**kwargs)  # ¯\_(ツ)_/¯
+
+    def plot(self, **kwargs):
+        if self.ndim == 1:
+            return self.plot1d(**kwargs)
+        elif self.ndim == 2:
+            return self.plot2d(**kwargs)
+        elif self.ndim == 3:
+            return self.plot3d(**kwargs)
+
+    def dist_plot(self, **kwargs):
+        """
+        plot phase space of a particle hierarchy
+        """
+        import copy
+
+        from .plotting import dist_plot as dp
+
+        usr_lvls = kwargs.get("levels", (0,))
+        finest = kwargs.get("finest", False)
+        pops = kwargs.get("pop", [])
+        time = kwargs.get("time", self.times()[0])
+        axis = kwargs.get("axis", ("Vx", "Vy"))
+        all_pops = list(self.level(0, time).patches[0].patch_datas.keys())
+
+        vmin = kwargs.get("vmin", -2)
+        vmax = kwargs.get("vmax", 2)
+        dv = kwargs.get("dv", 0.05)
+        vbins = vmin + dv * np.arange(int((vmax - vmin) / dv))
+
+        if finest:
+            final = finest_part_data(self)
+            if axis[0] == "x":
+                xbins = amr_grid(self, time)
+                bins = (xbins, vbins)
+            else:
+                bins = (vbins, vbins)
+            kwargs["bins"] = bins
+
+        else:
+            final = {pop: None for pop in all_pops}
+            for lvl_nbr, level in self.levels(time).items():
+                if lvl_nbr not in usr_lvls:
+                    continue
+                for ip, patch in enumerate(level.patches):
+                    if len(pops) == 0:
+                        pops = list(patch.patch_datas.keys())
+
+                    for pop in pops:
+                        tmp = copy.copy(patch.patch_datas[pop].dataset)
+
+                        if final[pop] is None:
+                            final[pop] = tmp
+                        else:
+                            final[pop].add(tmp)
+
+        # select particles
+        if "select" in kwargs:
+            for pop, particles in final.items():
+                final[pop] = kwargs["select"](particles)
+
+        return final, dp(final, **kwargs)
+
+
+def amr_grid(hierarchy, time):
+    """returns a non-uniform contiguous primal grid
+    associated to the given hierarchy
+    """
+    lvlPatchBoxes = {ilvl: [] for ilvl in range(hierarchy.finest_level() + 1)}
+    finalCells = {ilvl: None for ilvl in range(hierarchy.finest_level() + 1)}
+    lvl = hierarchy.levels(time)
+
+    for ilvl in range(hierarchy.finest_level(time) + 1)[::-1]:
+        sorted_patches = sorted(lvl[ilvl].patches, key=lambda p: p.layout.box.lower[0])
+
+        for ip, patch in enumerate(sorted_patches):
+            box = patch.layout.box
+            lvlPatchBoxes[ilvl].append(box)
+
+            # we create a list of all cells in the current patch
+            # remember that if the box upper cell is, say = 40,
+            # it means that the upper node is the lower node of cell 41
+            # so to get all primal nodes of a patch we need to include
+            # one past the upper cell.
+            # this said we do not want to include that last primal nodes
+            # all the time because that would be a duplicate with the lower
+            # node of the next patch. We only want to add it for the LAST
+            # (because sorted) patch. We also do not want to do it on levels
+            # other than L0 because the last primal node of the last patch
+            # of L_i is the first primal node of a L_{i-1} node, so including it
+            # would also mean adding a duplicate.
+            last = 1 if ilvl == 0 and ip == len(sorted_patches) - 1 else 0
+            cells = np.arange(box.lower[0], box.upper[0] + 1 + last)
+
+            # finest level has no next finer so we take all cells
+            if ilvl == hierarchy.finest_level(time):
+                if finalCells[ilvl] is None:
+                    finalCells[ilvl] = cells
+                else:
+                    finalCells[ilvl] = np.concatenate((finalCells[ilvl], cells))
+
+            else:
+                # on other levels
+                # we take only grids not overlaped by next finer
+                coarsenedNextFinerBoxes = [
+                    boxm.coarsen(b, refinement_ratio) for b in lvlPatchBoxes[ilvl + 1]
+                ]
+                for coarseBox in coarsenedNextFinerBoxes:
+                    ccells = np.arange(coarseBox.lower[0], coarseBox.upper[0] + 1)
+                    inter, icells, iccells = np.intersect1d(
+                        cells, ccells, return_indices=True
+                    )
+                    cells = np.delete(cells, icells)
+                if len(cells):
+                    if finalCells[ilvl] is None:
+                        finalCells[ilvl] = cells
+                    else:
+                        finalCells[ilvl] = np.unique(
+                            np.concatenate((finalCells[ilvl], cells))
+                        )
+
+    # now we have all cells for each level we
+    # just need to compute the primal coordinates
+    # and concatenate in a single array
+    for ilvl in range(hierarchy.finest_level() + 1):
+        if ilvl == 0:
+            x = finalCells[ilvl] * hierarchy.level(ilvl).patches[0].layout.dl[0]
+        else:
+            xx = finalCells[ilvl] * hierarchy.level(ilvl).patches[0].layout.dl[0]
+            x = np.concatenate((x, xx))
+
+    return np.sort(x)
+
+
+def is_root_lvl(patch_level):
+    return patch_level.level_number == 0
+
+
+field_qties = {
+    "EM_B_x": "Bx",
+    "EM_B_y": "By",
+    "EM_B_z": "Bz",
+    "EM_E_x": "Ex",
+    "EM_E_y": "Ey",
+    "EM_E_z": "Ez",
+    "flux_x": "Fx",
+    "flux_y": "Fy",
+    "flux_z": "Fz",
+    "bulkVelocity_x": "Vx",
+    "bulkVelocity_y": "Vy",
+    "bulkVelocity_z": "Vz",
+    "momentum_tensor_xx": "Mxx",
+    "momentum_tensor_xy": "Mxy",
+    "momentum_tensor_xz": "Mxz",
+    "momentum_tensor_yy": "Myy",
+    "momentum_tensor_yz": "Myz",
+    "momentum_tensor_zz": "Mzz",
+    "density": "rho",
+    "mass_density": "rho",
+    "tags": "tags",
+}
+
+
+particle_files_patterns = ("domain", "patchGhost", "levelGhost")
+
+
+def is_particle_file(filename):
+    return any([pattern in filename for pattern in particle_files_patterns])
+
+
+def particle_dataset_name(basename):
+    """
+    return "alpha_domain" from ions_pop_alpha_domain.h5
+    """
+    popname = basename.strip(".h5").split("_")[-2]
+    part_type = basename.strip(".h5").split("_")[-1]
+    dataset_name = popname + "_" + part_type
+
+    return dataset_name
+
+
+def is_pop_fluid_file(basename):
+    return (is_particle_file(basename) is False) and "pop" in basename
+
+
+def make_layout(h5_patch_grp, cell_width, interp_order):
+    origin = h5_patch_grp.attrs["origin"]
+    upper = h5_patch_grp.attrs["upper"]
+    lower = h5_patch_grp.attrs["lower"]
+    return GridLayout(Box(lower, upper), origin, cell_width, interp_order=interp_order)
+
+
+def create_from_all_times(time, hier):
+    return time is None and hier is None
+
+
+def create_from_one_time(time, hier):
+    return time is not None and hier is None
+
+
+def load_all_times(time, hier):
+    return time is None and hier is not None
+
+
+def load_one_time(time, hier):
+    return time is not None and hier is not None
+
+
+def compute_hier_from_(h, compute):
+    """
+    returns a hierarchy resulting from calling 'compute'
+    on each patch of the given hierarchy 'h'
+
+    compute is a function taking a Patch and returning
+    a list of dicts with the following keys:
+
+        "data": ndarray containing the data
+        "name": str, name of the data living on that patch, must be unique
+        "centering": str, ["dual", "primal"]
+
+     caveat: routine only works in 1D so far.
+    """
+    assert len(h.time_hier) == 1  # only single time hierarchies now
+    patch_levels = {}
+    for ilvl, lvl in h.patch_levels.items():
+        patches = {}
+        for ip, patch in enumerate(lvl.patches):
+            new_patch_datas = {}
+            layout = patch.layout
+            datas = compute(patch)
+            for data in datas:
+                pd = FieldData(
+                    layout, data["name"], data["data"], centering=data["centering"]
+                )
+                new_patch_datas[data["name"]] = pd
+            if ilvl not in patches:
+                patches[ilvl] = []
+            patches[ilvl].append(Patch(new_patch_datas, patch.id))
+
+        patch_levels[ilvl] = PatchLevel(ilvl, patches[ilvl])
+
+    t = list(h.time_hier.keys())[0]
+    return PatchHierarchy(patch_levels, h.domain_box, refinement_ratio, time=t)
+
+
+def are_compatible_hierarchies(hierarchies):
+    return True
+
+
+def extract_patchdatas(hierarchies, ilvl, ipatch):
+    """
+    returns a dict {patchdata_name:patchdata} from a list of hierarchies for patch ipatch at level ilvl
+    """
+    patches = [h.patch_levels[ilvl].patches[ipatch] for h in hierarchies]
+    patch_datas = {
+        pdname: pd for p in patches for pdname, pd in list(p.patch_datas.items())
+    }
+    return patch_datas
+
+
+def new_patchdatas_from(compute, patchdatas, layout, id, **kwargs):
+    new_patch_datas = {}
+    datas = compute(patchdatas, id=id, **kwargs)
+    for data in datas:
+        pd = FieldData(layout, data["name"], data["data"], centering=data["centering"])
+        new_patch_datas[data["name"]] = pd
+    return new_patch_datas
+
+
+def new_patches_from(compute, hierarchies, ilvl, **kwargs):
+    reference_hier = hierarchies[0]
+    new_patches = []
+    patch_nbr = len(reference_hier.patch_levels[ilvl].patches)
+    for ip in range(patch_nbr):
+        current_patch = reference_hier.patch_levels[ilvl].patches[ip]
+        layout = current_patch.layout
+        patch_datas = extract_patchdatas(hierarchies, ilvl, ip)
+        new_patch_datas = new_patchdatas_from(
+            compute, patch_datas, layout, id=current_patch.id, **kwargs
+        )
+        new_patches.append(Patch(new_patch_datas, current_patch.id))
+    return new_patches
+
+
+def compute_hier_from(compute, hierarchies, **kwargs):
+    """return a new hierarchy using the callback 'compute' on all patchdatas
+    of the given hierarchies
+    """
+    if not are_compatible_hierarchies(hierarchies):
+        raise RuntimeError("hierarchies are not compatible")
+
+    hierarchies = listify(hierarchies)
+    reference_hier = hierarchies[0]
+    domain_box = reference_hier.domain_box
+    patch_levels = {}
+    for ilvl in range(reference_hier.levelNbr()):
+        patch_levels[ilvl] = PatchLevel(
+            ilvl, new_patches_from(compute, hierarchies, ilvl, **kwargs)
+        )
+
+    assert len(reference_hier.time_hier) == 1  # only single time hierarchies now
+    t = list(reference_hier.time_hier.keys())[0]
+    return PatchHierarchy(patch_levels, domain_box, refinement_ratio, time=t)
+
+
+def pop_name(basename):
+    return basename.strip(".h5").split("_")[2]
+
+
+def add_to_patchdata(patch_datas, h5_patch_grp, basename, layout):
+    """
+    adds data in the h5_patch_grp in the given PatchData dict
+    returns True if valid h5 patch found
+    """
+
+    if is_particle_file(basename):
+        v = np.asarray(h5_patch_grp["v"])
+        s = v.size
+        v = v[:].reshape(int(s / 3), 3)
+        nbrParts = v.shape[0]
+        dl = np.zeros((nbrParts, layout.ndim))
+        for i in range(layout.ndim):
+            dl[:, i] = layout.dl[i]
+
+        particles = Particles(
+            icells=h5_patch_grp["iCell"],
+            deltas=h5_patch_grp["delta"],
+            v=v,
+            weights=h5_patch_grp["weight"],
+            charges=h5_patch_grp["charge"],
+            dl=dl,
+        )
+
+        pdname = particle_dataset_name(basename)
+        if pdname in patch_datas:
+            raise ValueError("error - {} already in patchdata".format(pdname))
+
+        patch_datas[pdname] = ParticleData(layout, particles, pop_name(basename))
+
+    else:
+        for dataset_name in h5_patch_grp.keys():
+            dataset = h5_patch_grp[dataset_name]
+
+            if dataset_name not in field_qties:
+                raise RuntimeError(
+                    "invalid dataset name : {} is not in {}".format(
+                        dataset_name, field_qties
+                    )
+                )
+
+            pdata = FieldData(layout, field_qties[dataset_name], dataset)
+
+            pdata_name = field_qties[dataset_name]
+
+            if is_pop_fluid_file(basename):
+                pdata_name = pop_name(basename) + "_" + pdata_name
+
+            if dataset_name in patch_datas:
+                raise ValueError("error - {} already in patchdata".format(dataset_name))
+
+            patch_datas[pdata_name] = pdata
+
+    return True  # valid patch assumed
+
+
+def patch_has_datasets(h5_patch_grp):
+    return len(h5_patch_grp.keys()) > 0
+
+
+h5_time_grp_key = "t"
+
+
+def hierarchy_fromh5(h5_filename, time, hier, silent=True):
+    import h5py
+
+    data_file = h5py.File(h5_filename, "r")
+    basename = os.path.basename(h5_filename)
+
+    root_cell_width = np.asarray(data_file.attrs["cell_width"])
+    interp = data_file.attrs["interpOrder"]
+    domain_box = Box(
+        [0] * len(data_file.attrs["domain_box"]), data_file.attrs["domain_box"]
+    )
+
+    if create_from_all_times(time, hier):
+        # first create from first time
+        # then add all other times
+        if not silent:
+            print("creating hierarchy from all times in file")
+        times = list(data_file[h5_time_grp_key].keys())
+        hier = hierarchy_fromh5(h5_filename, time=times[0], hier=hier)
+        if len(times) > 1:
+            for t in times[1:]:
+                hierarchy_fromh5(h5_filename, t, hier)
+        return hier
+
+    if create_from_one_time(time, hier):
+        if not silent:
+            print("creating hierarchy from time {}".format(time))
+        t = time
+
+        h5_time_grp = data_file[h5_time_grp_key][time]
+        patch_levels = {}
+
+        for plvl_key in h5_time_grp.keys():
+            h5_patch_lvl_grp = h5_time_grp[plvl_key]
+            ilvl = int(plvl_key[2:])
+            lvl_cell_width = root_cell_width / refinement_ratio**ilvl
+            patches = {}
+
+            if ilvl not in patches:
+                patches[ilvl] = []
+
+            for pkey in h5_patch_lvl_grp.keys():
+                h5_patch_grp = h5_patch_lvl_grp[pkey]
+
+                if patch_has_datasets(h5_patch_grp):
+                    patch_datas = {}
+                    layout = make_layout(h5_patch_grp, lvl_cell_width, interp)
+                    add_to_patchdata(patch_datas, h5_patch_grp, basename, layout)
+
+                    patches[ilvl].append(
+                        Patch(patch_datas, h5_patch_grp.name.split("/")[-1])
+                    )
+
+                    patch_levels[ilvl] = PatchLevel(ilvl, patches[ilvl])
+
+        diag_hier = PatchHierarchy(
+            patch_levels, domain_box, refinement_ratio, t, data_file
+        )
+
+        return diag_hier
+
+    if load_one_time(time, hier):
+        if not silent:
+            print("loading data at time {} into existing hierarchy".format(time))
+        h5_time_grp = data_file[h5_time_grp_key][time]
+        t = time
+
+        if t in hier.time_hier:
+            if not silent:
+                print("time already exist, adding data...")
+
+            # time already exists in the hierarchy
+            # all we need to do is adding the data
+            # as patchDatas in the appropriate patches
+            # and levels, if data compatible with hierarchy
+
+            patch_levels = hier.time_hier[t]
+
+            for plvl_key in h5_time_grp.keys():
+                ilvl = int(plvl_key[2:])
+                lvl_cell_width = root_cell_width / refinement_ratio**ilvl
+
+                for ipatch, pkey in enumerate(h5_time_grp[plvl_key].keys()):
+                    h5_patch_grp = h5_time_grp[plvl_key][pkey]
+
+                    if patch_has_datasets(h5_patch_grp):
+                        hier_patch = patch_levels[ilvl].patches[ipatch]
+                        origin = h5_time_grp[plvl_key][pkey].attrs["origin"]
+                        upper = h5_time_grp[plvl_key][pkey].attrs["upper"]
+                        lower = h5_time_grp[plvl_key][pkey].attrs["lower"]
+                        file_patch_box = Box(lower, upper)
+
+                        assert file_patch_box == hier_patch.box
+                        assert (abs(origin - hier_patch.origin) < 1e-6).all()
+                        assert (abs(lvl_cell_width - hier_patch.dl) < 1e-6).all()
+
+                        layout = make_layout(h5_patch_grp, lvl_cell_width, interp)
+                        add_to_patchdata(
+                            hier_patch.patch_datas, h5_patch_grp, basename, layout
+                        )
+
+            return hier
+
+        if not silent:
+            print("adding data to new time")
+        # time does not exist in the hierarchy
+        # we have to create a brand new set of patchLevels
+        # containing patches, and load data in their patchdatas
+
+        patch_levels = {}
+
+        for plvl_key in h5_time_grp.keys():
+            ilvl = int(plvl_key[2:])
+
+            lvl_cell_width = root_cell_width / refinement_ratio**ilvl
+            lvl_patches = []
+
+            for ipatch, pkey in enumerate(h5_time_grp[plvl_key].keys()):
+                h5_patch_grp = h5_time_grp[plvl_key][pkey]
+
+                if patch_has_datasets(h5_patch_grp):
+                    layout = make_layout(h5_patch_grp, lvl_cell_width, interp)
+                    patch_datas = {}
+                    add_to_patchdata(patch_datas, h5_patch_grp, basename, layout)
+                    lvl_patches.append(Patch(patch_datas))
+
+            patch_levels[ilvl] = PatchLevel(ilvl, lvl_patches)
+
+        hier.time_hier[t] = patch_levels
+        return hier
+
+    if load_all_times(time, hier):
+        if not silent:
+            print("loading all times in existing hier")
+        for time in data_file[h5_time_grp_key].keys():
+            hier = hierarchy_fromh5(h5_filename, time, hier)
+
+        return hier
+
+
+def quantidic(ilvl, wrangler):
+    pl = wrangler.getPatchLevel(ilvl)
+
+    return {
+        "density": pl.getDensity,
+        "bulkVelocity_x": pl.getVix,
+        "bulkVelocity_y": pl.getViy,
+        "bulkVelocity_z": pl.getViz,
+        "EM_B_x": pl.getBx,
+        "EM_B_y": pl.getBy,
+        "EM_B_z": pl.getBz,
+        "EM_E_x": pl.getEx,
+        "EM_E_y": pl.getEy,
+        "EM_E_z": pl.getEz,
+        "flux_x": pl.getFx,
+        "flux_y": pl.getFy,
+        "flux_z": pl.getFz,
+        "particles": pl.getParticles,
+    }
+
+
+def isFieldQty(qty):
+    return qty in (
+        "density",
+        "bulkVelocity_x",
+        "bulkVelocity_y",
+        "bulkVelocity_z",
+        "EM_B_x",
+        "EM_B_y",
+        "EM_B_z",
+        "EM_E_x",
+        "EM_E_y",
+        "EM_E_z",
+        "flux_x",
+        "flux_y",
+        "flux_z",
+        "momentumTensor_xx",
+        "momentumTensor_xy",
+        "momentumTensor_xz",
+        "momentumTensor_yy",
+        "momentumTensor_yz",
+        "momentumTensor_zz",
+    )
+
+
+def hierarchy_from_sim(simulator, qty, pop=""):
+    dw = simulator.data_wrangler()
+    nbr_levels = dw.getNumberOfLevels()
+    patch_levels = {}
+
+    root_cell_width = simulator.cell_width()
+    domain_box = Box([0] * len(root_cell_width), simulator.domain_box())
+    assert len(domain_box.ndim) == len(simulator.domain_box().ndim)
+
+    for ilvl in range(nbr_levels):
+        lvl_cell_width = root_cell_width / refinement_ratio**ilvl
+
+        patches = {ilvl: [] for ilvl in range(nbr_levels)}
+        getters = quantidic(ilvl, dw)
+
+        if isFieldQty(qty):
+            wpatches = getters[qty]()
+            for patch in wpatches:
+                patch_datas = {}
+                lower = patch.lower
+                upper = patch.upper
+                origin = patch.origin
+                layout = GridLayout(
+                    Box(lower, upper),
+                    origin,
+                    lvl_cell_width,
+                    interp_order=simulator.interporder(),
+                )
+                pdata = FieldData(layout, field_qties[qty], patch.data)
+                patch_datas[qty] = pdata
+                patches[ilvl].append(Patch(patch_datas))
+
+        elif qty == "particles":
+            if pop == "":
+                raise ValueError("must specify pop argument for particles")
+            # here the getter returns a dict like this
+            # {'protons': {'patchGhost': [<pybindlibs.cpp.PatchDataContiguousParticles_1 at 0x119f78970>,
+            # <pybindlibs.cpp.PatchDataContiguousParticles_1 at 0x119f78f70>],
+            # 'domain': [<pybindlibs.cpp.PatchDataContiguousParticles_1 at 0x119f78d70>,
+            # <pybindlibs.cpp.PatchDataContiguousParticles_1 at 0x119f78770>]}}
+
+            # domain particles are assumed to always be here
+            # but patchGhost and levelGhost may not be, depending on the level
+
+            populationdict = getters[qty](pop)[pop]
+
+            dom_dw_patches = populationdict["domain"]
+            for patch in dom_dw_patches:
+                patch_datas = {}
+
+                lower = patch.lower
+                upper = patch.upper
+                origin = patch.origin
+                layout = GridLayout(
+                    Box(lower, upper),
+                    origin,
+                    lvl_cell_width,
+                    interp_order=simulator.interp_order(),
+                )
+                v = np.asarray(patch.data.v).reshape(int(len(patch.data.v) / 3), 3)
+
+                domain_particles = Particles(
+                    icells=np.asarray(patch.data.iCell),
+                    deltas=np.asarray(patch.data.delta),
+                    v=v,
+                    weights=np.asarray(patch.data.weight),
+                    charges=np.asarray(patch.data.charge),
+                )
+
+                patch_datas[pop + "_particles"] = ParticleData(
+                    layout, domain_particles, pop
+                )
+                patches[ilvl].append(Patch(patch_datas))
+
+            # ok now let's add the patchGhost if present
+            # note that patchGhost patches may not be the same list as the
+            # domain patches... since not all patches may not have patchGhost while they do have
+            # domain... while looping on the patchGhost items, we need to search in
+            # the already created patches which one to which add the patchGhost particles
+
+            for ghostParticles in ["patchGhost", "levelGhost"]:
+                if ghostParticles in populationdict:
+                    for dwpatch in populationdict[ghostParticles]:
+                        v = np.asarray(dwpatch.data.v)
+                        s = v.size
+                        v = v[:].reshape(int(s / 3), 3)
+
+                        patchGhost_part = Particles(
+                            icells=np.asarray(dwpatch.data.iCell),
+                            deltas=np.asarray(dwpatch.data.delta),
+                            v=v,
+                            weights=np.asarray(dwpatch.data.weight),
+                            charges=np.asarray(dwpatch.data.charge),
+                        )
+
+                        box = Box(dwpatch.lower, dwpatch.upper)
+
+                        # now search which of the already created patches has the same box
+                        # once found we add the new particles to the ones already present
+
+                        patch = [p for p in patches[ilvl] if p.box == box][0]
+                        patch.patch_datas[pop + "_particles"].dataset.add(
+                            patchGhost_part
+                        )
+
+        else:
+            raise ValueError("{} is not a valid quantity".format(qty))
+
+        patch_levels[ilvl] = PatchLevel(ilvl, patches[ilvl])
+
+    return PatchHierarchy(patch_levels, domain_box, time=simulator.currentTime())
+
+
+def hierarchy_from_box(domain_box, ghosts_nbr):
+    """
+    constructs a basic hierarchy with one patch for level 0 as the entire domain
+    """
+    layout = GridLayout(
+        domain_box, np.asarray([0] * domain_box.ndim), [0.1] * domain_box.ndim, 1
+    )
+    pdata = PatchData(layout, "qty")
+    object.__setattr__(pdata, "ghosts_nbr", np.asarray(ghosts_nbr))
+    object.__setattr__(pdata, "ghost_box", boxm.grow(layout.box, ghosts_nbr))
+    return PatchHierarchy({0: PatchLevel(0, [Patch({"qty": pdata})])}, domain_box)
+
+
+def hierarchy_from(
+    simulator=None, qty=None, pop="", h5_filename=None, time=None, hier=None
+):
+    """
+    this function reads an HDF5 PHARE file and returns a PatchHierarchy from
+    which data is accessible.
+    if 'time' is None, all times in the file will be read, if a time is given
+    then only that time will be read
+    if 'hier' is None, then a new hierarchy will be created, if not then the
+    given hierarchy 'hier' will be filled.
+
+    The function fails if the data is already in hierarchy
+    """
+
+    if simulator is not None and h5_filename is not None:
+        raise ValueError("cannot pass both a simulator and a h5 file")
+
+    if h5_filename is not None:
+        return hierarchy_fromh5(h5_filename, time, hier)
+
+    if simulator is not None and qty is not None:
+        return hierarchy_from_sim(simulator, qty, pop=pop)
+
+    raise ValueError("can't make hierarchy")
+
+
+def merge_particles(hierarchy):
+    for time, patch_levels in hierarchy.time_hier.items():
+        for ilvl, plvl in patch_levels.items():
+            for ip, patch in enumerate(plvl.patches):
+                pdatas = patch.patch_datas
+                domain_pdata = [
+                    (pdname, pd) for pdname, pd in pdatas.items() if "domain" in pdname
+                ][0]
+
+                pghost_pdatas = [
+                    (pdname, pd)
+                    for pdname, pd in pdatas.items()
+                    if "patchGhost" in pdname
+                ]
+                lghost_pdatas = [
+                    (pdname, pd)
+                    for pdname, pd in pdatas.items()
+                    if "levelGhost" in pdname
+                ]
+
+                pghost_pdata = pghost_pdatas[0] if pghost_pdatas else None
+                lghost_pdata = lghost_pdatas[0] if lghost_pdatas else None
+
+                if pghost_pdata is not None:
+                    domain_pdata[1].dataset.add(pghost_pdata[1].dataset)
+                    del pdatas[pghost_pdata[0]]
+
+                if lghost_pdata is not None:
+                    domain_pdata[1].dataset.add(lghost_pdata[1].dataset)
+                    del pdatas[lghost_pdata[0]]
+
+                popname = domain_pdata[0].split("_")[0]
+                pdatas[popname + "_particles"] = pdatas[domain_pdata[0]]
+                del pdatas[domain_pdata[0]]
+
+
+def h5_filename_from(diagInfo):
+    # diagInfo.quantity starts with a / , hence   [1:]
+    return (diagInfo.quantity + ".h5").replace("/", "_")[1:]
+
+
+def get_times_from_h5(filepath):
+    import h5py
+
+    f = h5py.File(filepath, "r")
+    times = np.array(sorted([float(s) for s in list(f["t"].keys())]))
+    f.close()
+    return times
+
+
+def getPatch(hier, point):
+    """
+    convenience function mainly for debug. returns a dict
+    {ilevel:patch}  for patches in which the given point is
+    """
+    patches = {}
+    counts = {ilvl: 0 for ilvl in hier.levels().keys()}
+    for ilvl, lvl in hier.levels().items():
+        for p in lvl.patches:
+            px, py = point
+            dx, dy = p.layout.dl
+            ix = int(px / dx)
+            iy = int(py / dy)
+            if (ix, iy) in p.box:
+                patches[ilvl] = p
+                counts[ilvl] += 1
+    for k, v in counts.items():
+        if v > 1:
+            print("error : ", k, v)
+            raise RuntimeError("more than one patch found for point")
+    return patches
diff --git a/pyphare/pyphare/pharesee/run/utils.py b/pyphare/pyphare/pharesee/run/utils.py
index e9cb5b57d..1b7a9c499 100644
--- a/pyphare/pyphare/pharesee/run/utils.py
+++ b/pyphare/pyphare/pharesee/run/utils.py
@@ -359,7 +359,7 @@ def _get_rank(patchdatas, patch_id, **kwargs):
     layout = reference_pd.layout
     centering = "dual"
     nbrGhosts = layout.nbrGhosts(layout.interp_order, centering)
-    shape = grow(reference_pd.box, [nbrGhosts] * 2).shape
+    shape = grow(reference_pd.box, [nbrGhosts] * ndim).shape
 
     if ndim == 1:
         pass
diff --git a/res/amr/splitting.yml b/res/amr/splitting.yml
index e71aa05e8..14ecdb96d 100644
--- a/res/amr/splitting.yml
+++ b/res/amr/splitting.yml
@@ -99,3 +99,43 @@ dimension_2:
       delta:  1 2
       weight: .140625 .09375 .0625 .0234375 .015625 .00390625
 
+
+dimension_3:
+  interp_1:
+    N_particles_6:
+      delta:  .966431
+      weight: .166666
+
+    N_particles_12:
+      delta:  .74823
+      weight: .083333
+
+    N_particles_27:
+      delta:  1
+      weight: .125 .0625
+
+  interp_2:
+    N_particles_6:
+      delta:  1.149658
+      weight: .166666
+
+    N_particles_12:
+      delta:  .888184
+      weight: .083333
+
+    N_particles_27:
+      delta:  1.111333
+      weight: .099995 .055301
+
+  interp_3:
+    N_particles_6:
+      delta:  1.312622
+      weight: .166666
+
+    N_particles_12:
+      delta:  1.012756
+      weight: .083333
+
+    N_particles_27:
+      delta:  1.276815
+      weight: .104047 .05564
diff --git a/res/cmake/funcs.cmake b/res/cmake/funcs.cmake
new file mode 100644
index 000000000..922aca10a
--- /dev/null
+++ b/res/cmake/funcs.cmake
@@ -0,0 +1,192 @@
+# public test functions
+#
+#   add_phare_test($binary $directory)
+#    execute binary in target directory, with mpirun when -DtestMPI=ON
+#
+#   add_python3_test($name $file $directory)
+#    launch python3 file described by name in target directory, with mpirun when -DtestMPI=ON
+#
+#   add_no_mpi_phare_test($binary $directory)
+#    execute binary in target directory, does not run when -DtestMPI=ON
+#
+#   add_no_mpi_python3_test($name $file $directory)
+#    launch python3 file described by name in target directory, does not run when -DtestMPI=ON
+#
+#  phare_exec(level target exe directory)
+#    execute exe identified by target in directory
+#      if level >= PHARE_EXEC_LEVEL_MIN AND level <= PHARE_EXEC_LEVEL_MAX
+#
+#  phare_python3_exec(level target file directory $ARGS)
+#    execute file identified by target in directory
+#      if level >= PHARE_EXEC_LEVEL_MIN AND level <= PHARE_EXEC_LEVEL_MAX
+#
+#  Note to developers - do not use cmake variable function arguments for functions
+#      phare_python3_exec
+#      phare_mpi_python3_exec
+#   if these function calls are to files executing python unit tests as they will interfere
+
+if (test AND ${PHARE_EXEC_LEVEL_MIN} GREATER 0) # 0 = no tests
+
+  if (NOT DEFINED PHARE_MPI_PROCS)
+    set(PHARE_MPI_PROCS 1)
+    if(testMPI)
+      set(PHARE_MPI_PROCS 2)
+    endif()
+  endif()
+
+  function(set_exe_paths_ binary)
+    set_property(TEST ${binary}        PROPERTY ENVIRONMENT "PYTHONPATH=${PHARE_PYTHONPATH}")
+    # ASAN detects leaks by default, even in system/third party libraries
+    set_property(TEST ${binary} APPEND PROPERTY ENVIRONMENT "ASAN_OPTIONS=detect_leaks=0")
+    set_property(TEST ${binary} APPEND PROPERTY ENVIRONMENT PHARE_SKIP_CLI=1 )
+  endfunction(set_exe_paths_)
+
+  function(add_phare_test_ binary directory)
+    target_compile_options(${binary} PRIVATE ${PHARE_WERROR_FLAGS} -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE})
+    set_exe_paths_(${binary})
+    set_property(TEST ${binary} APPEND PROPERTY ENVIRONMENT GMON_OUT_PREFIX=gprof.${binary})
+    set_property(TEST ${binary} APPEND PROPERTY ENVIRONMENT PHARE_MPI_PROCS=${PHARE_MPI_PROCS})
+    if(testDuringBuild)
+      add_custom_command(
+           TARGET ${binary}
+           POST_BUILD
+           WORKING_DIRECTORY ${directory}
+           COMMAND ${CMAKE_CTEST_COMMAND} -C $<CONFIGURATION> -R \"^${binary}$$\" --output-on-failures
+      )
+    endif(testDuringBuild)
+  endfunction(add_phare_test_)
+
+  function(add_no_mpi_phare_test binary directory)
+    if(NOT testMPI OR (testMPI AND forceSerialTests))
+      add_test(NAME ${binary} COMMAND ./${binary} WORKING_DIRECTORY ${directory})
+      add_phare_test_(${binary} ${directory})
+    else()
+      # this prevents building targets even when added via "add_executable"
+      set_target_properties(${binary} PROPERTIES EXCLUDE_FROM_ALL 1 EXCLUDE_FROM_DEFAULT_BUILD 1)
+    endif()
+  endfunction(add_no_mpi_phare_test)
+
+  function(add_no_mpi_python3_test name file directory)
+    if(NOT testMPI OR (testMPI AND forceSerialTests))
+      add_test(NAME py3_${name} COMMAND python3 -u ${file} WORKING_DIRECTORY ${directory})
+      set_exe_paths_(py3_${name})
+    endif()
+  endfunction(add_no_mpi_python3_test)
+
+  if(testMPI)
+    function(add_phare_test binary directory)
+      add_test(NAME ${binary} COMMAND mpirun -n ${PHARE_MPI_PROCS} ${PHARE_MPIRUN_POSTFIX} ./${binary} WORKING_DIRECTORY ${directory})
+      add_phare_test_(${binary} ${directory})
+    endfunction(add_phare_test)
+
+    function(add_python3_test name file directory)
+      add_test(NAME py3_${name} COMMAND mpirun -n ${PHARE_MPI_PROCS} ${PHARE_MPIRUN_POSTFIX} python3 -u ${file} WORKING_DIRECTORY ${directory})
+      set_exe_paths_(py3_${name})
+    endfunction(add_python3_test)
+
+    function(add_mpi_python3_test N name file directory)
+      add_test(NAME py3_${name}_mpi_n_${N} COMMAND mpirun -n ${N} ${PHARE_MPIRUN_POSTFIX} python3 ${file} WORKING_DIRECTORY ${directory})
+      set_exe_paths_(py3_${name}_mpi_n_${N})
+    endfunction(add_mpi_python3_test)
+
+  else()
+    function(add_phare_test binary directory)
+      add_no_mpi_phare_test(${binary} ${directory})
+    endfunction(add_phare_test)
+
+    function(add_python3_test name file directory)
+      add_no_mpi_python3_test(${name} ${file} ${directory})
+    endfunction(add_python3_test)
+
+    function(add_mpi_python3_test N name file directory)
+      # do nothing
+    endfunction(add_mpi_python3_test)
+  endif(testMPI)
+
+
+  function(add_phare_build_test binary directory)
+  endfunction(add_phare_build_test)
+
+
+  if(DEFINED GTEST_ROOT)
+    set(GTEST_ROOT ${GTEST_ROOT} CACHE PATH "Path to googletest")
+    find_package(GTest REQUIRED)
+    set(GTEST_LIBS GTest::GTest GTest::Main)
+  else()
+    set(GTEST_ROOT ${CMAKE_CURRENT_SOURCE_DIR}/subprojects/googletest)
+
+    if (NOT EXISTS ${GTEST_ROOT})
+      execute_process(COMMAND ${Git} clone https://github.com/google/googletest ${GTEST_ROOT})
+    endif()
+
+    add_subdirectory(subprojects/googletest)
+    set(GTEST_INCLUDE_DIRS
+      $<BUILD_INTERFACE:${gtest_SOURCE_DIR}/include>
+      $<BUILD_INTERFACE:${gmock_SOURCE_DIR}/include>)
+    set(GTEST_LIBS gtest gmock)
+
+  endif()
+
+  function(phare_exec level target exe directory)
+    if(${level} GREATER_EQUAL ${PHARE_EXEC_LEVEL_MIN} AND ${level} LESS_EQUAL ${PHARE_EXEC_LEVEL_MAX})
+      add_test(NAME ${target} COMMAND ${exe} WORKING_DIRECTORY ${directory})
+    endif()
+  endfunction(phare_exec)
+  # use
+  #  phare_exec(1 test_id ./binary ${CMAKE_CURRENT_BINARY_DIR})
+
+
+  function(phare_mpi_python3_exec level N target file directory)
+    if(${level} GREATER_EQUAL ${PHARE_EXEC_LEVEL_MIN} AND ${level} LESS_EQUAL ${PHARE_EXEC_LEVEL_MAX})
+      string (REPLACE ";" " " CLI_ARGS "${ARGN}")
+      if(${N} EQUAL 1)
+        add_test(
+            NAME py3_${target}
+            COMMAND python3 -u ${file} ${CLI_ARGS}
+            WORKING_DIRECTORY ${directory})
+        set_exe_paths_(py3_${target})
+      else()
+        add_test(
+            NAME py3_${target}_mpi_n_${N}
+            COMMAND mpirun -n ${N} ${PHARE_MPIRUN_POSTFIX} python3 -u ${file} ${CLI_ARGS}
+            WORKING_DIRECTORY ${directory})
+        set_exe_paths_(py3_${target}_mpi_n_${N})
+      endif()
+    endif()
+  endfunction(phare_mpi_python3_exec)
+  # use
+  #  phare_mpi_python3_exec(1 2 test_id script.py ${CMAKE_CURRENT_BINARY_DIR} $ARGS)
+
+
+  if(testMPI)
+    function(phare_python3_exec level target file directory)
+      phare_mpi_python3_exec(${level} 2 ${target} ${file} ${directory})
+    endfunction(phare_python3_exec)
+  else()
+    function(phare_python3_exec level target file directory)
+      if(${level} GREATER_EQUAL ${PHARE_EXEC_LEVEL_MIN} AND ${level} LESS_EQUAL ${PHARE_EXEC_LEVEL_MAX})
+        string (REPLACE ";" " " CLI_ARGS "${ARGN}")
+        add_test(NAME py3_${target} COMMAND python3 -u ${file} ${CLI_ARGS} WORKING_DIRECTORY ${directory})
+        set_exe_paths_(py3_${target})
+      endif()
+    endfunction(phare_python3_exec)
+  endif(testMPI)
+  # use
+  #  phare_python3_exec(1 test_id script.py ${CMAKE_CURRENT_BINARY_DIR} $ARGS)
+
+
+  set(GTEST_INCLUDE_DIRS ${GTEST_INCLUDE_DIRS} ${PHARE_PROJECT_DIR})
+
+  enable_testing()
+
+endif()
+
+# useful to see what's available after importing a package
+function(phare_print_all_vars)
+  get_cmake_property(_variableNames VARIABLES)
+  list (SORT _variableNames)
+  foreach (_variableName ${_variableNames})
+      message(STATUS "${_variableName}=${${_variableName}}")
+  endforeach()
+endfunction(phare_print_all_vars)
+
diff --git a/res/cmake/options.cmake b/res/cmake/options.cmake
index 20266c84a..29552db0c 100644
--- a/res/cmake/options.cmake
+++ b/res/cmake/options.cmake
@@ -70,10 +70,14 @@ option(withPhlop "Use phlop" OFF)
 # -DlowResourceTests=ON
 option(lowResourceTests "Disable heavy tests for CI (2d/3d/etc" OFF)
 
+# -DhighResourceTests=ON
+option(highResourceTests "Enable heavy tests for CI (3d/etc" OFF)
+
 # -DtestDuringBuild=ON enabled if devMode=ON, disabled if asan=ON (needs LD_PRELOAD)
 option(testDuringBuild "Runs C++ unit tests after they are built" OFF)
 
 
+
 # -DPGO_GEN=OFF profile guided optimization generate
 option(PGO_GEN "profile guided optimization generate" OFF)
 # -DPGO_USE=OFF
@@ -84,6 +88,7 @@ option(PGO_USE "profile guided optimization use" OFF)
 option(PHARE_COMPILER_WORKAROUNDS "Attempt silence compiler false positives" ON) # on by default
 
 
+
 # Controlling the activation of tests
 if (NOT DEFINED PHARE_EXEC_LEVEL_MIN)
   set(PHARE_EXEC_LEVEL_MIN 1)
diff --git a/src/amr/data/field/coarsening/magnetic_field_coarsener.hpp b/src/amr/data/field/coarsening/magnetic_field_coarsener.hpp
index 39d816413..e46abb1e9 100644
--- a/src/amr/data/field/coarsening/magnetic_field_coarsener.hpp
+++ b/src/amr/data/field/coarsening/magnetic_field_coarsener.hpp
@@ -31,110 +31,216 @@ using core::dirZ;
  * of the enclosed fine faces
  *
  */
-template<std::size_t dimension>
+template<std::size_t dim>
 class MagneticFieldCoarsener
 {
+    using Point_t = core::Point<int, dim>;
+
 public:
-    MagneticFieldCoarsener(std::array<core::QtyCentering, dimension> const centering,
+    MagneticFieldCoarsener(std::array<core::QtyCentering, dim> const centering,
                            SAMRAI::hier::Box const& sourceBox,
                            SAMRAI::hier::Box const& destinationBox,
                            SAMRAI::hier::IntVector const& ratio)
         : centering_{centering}
         , sourceBox_{sourceBox}
         , destinationBox_{destinationBox}
-
     {
     }
 
+
     template<typename FieldT>
-    void operator()(FieldT const& fineField, FieldT& coarseField,
-                    core::Point<int, dimension> coarseIndex)
+    void operator()(FieldT const& fineField, FieldT& coarseField, core::Point<int, dim> coarseIndex)
     {
         // For the moment we only take the case of field with the same centering
         TBOX_ASSERT(fineField.physicalQuantity() == coarseField.physicalQuantity());
 
-        core::Point<int, dimension> fineStartIndex;
+        coarsen<dim>(fine_start_index(coarseIndex), fineField, coarseField,
+                     AMRToLocal(coarseIndex, destinationBox_));
+    }
 
+private:
+    auto fine_start_index(core::Point<int, dim> const coarseIndex) const
+    {
+        core::Point<int, dim> fineStartIndex;
         fineStartIndex[dirX] = coarseIndex[dirX] * this->ratio_;
-
-        if constexpr (dimension > 1)
+        if constexpr (dim > 1)
         {
             fineStartIndex[dirY] = coarseIndex[dirY] * this->ratio_;
-            if constexpr (dimension > 2)
+            if constexpr (dim > 2)
             {
                 fineStartIndex[dirZ] = coarseIndex[dirZ] * this->ratio_;
             }
         }
+        return AMRToLocal(fineStartIndex, sourceBox_);
+    }
 
-        fineStartIndex = AMRToLocal(fineStartIndex, sourceBox_);
-        coarseIndex    = AMRToLocal(coarseIndex, destinationBox_);
 
-        // the following kinda assumes where B is, i.e. Yee layout centering
-        // as it only does faces pirmal-dual, dual-primal and dual-dual
+    template<std::size_t D, typename FieldT>
+    typename std::enable_if<D == 1, void>::type
+    coarsen(Point_t const fineStartIndex, FieldT const& fineField, FieldT& coarseField,
+            Point_t const coarseIndex);
 
-        if constexpr (dimension == 1)
-        {
-            // in 1D div(B) is automatically satisfied so using this coarsening
-            // opertor is probably not better than the default one, but we do that
-            // for a kind of consistency...
-            // coarse flux is equal to fine flux and we're 1D so there is flux partitioned
-            // only for By and Bz, Bx is equal to the fine value
+    template<std::size_t D, typename FieldT>
+    typename std::enable_if<D == 2, void>::type
+    coarsen(Point_t const fineStartIndex, FieldT const& fineField, FieldT& coarseField,
+            Point_t const coarseIndex);
+    template<std::size_t D, typename FieldT>
+    typename std::enable_if<D == 3, void>::type
+    coarsen(Point_t const fineStartIndex, FieldT const& fineField, FieldT& coarseField,
+            Point_t const coarseIndex);
 
-            if (centering_[dirX] == core::QtyCentering::primal) // bx
-            {
-                coarseField(coarseIndex[dirX]) = fineField(fineStartIndex[dirX]);
-            }
-            else if (centering_[dirX] == core::QtyCentering::dual) // by and bz
-            {
-                coarseField(coarseIndex[dirX])
-                    = 0.5 * (fineField(fineStartIndex[dirX] + 1) + fineField(fineStartIndex[dirX]));
-            }
-        }
-
-        if constexpr (dimension == 2)
-        {
-            if (centering_[dirX] == core::QtyCentering::primal
-                and centering_[dirY] == core::QtyCentering::dual)
-            {
-                coarseField(coarseIndex[dirX], coarseIndex[dirY])
-                    = 0.5
-                      * (fineField(fineStartIndex[dirX], fineStartIndex[dirY])
-                         + fineField(fineStartIndex[dirX], fineStartIndex[dirY] + 1));
-            }
-            else if (centering_[dirX] == core::QtyCentering::dual
-                     and centering_[dirY] == core::QtyCentering::primal)
-            {
-                coarseField(coarseIndex[dirX], coarseIndex[dirY])
-                    = 0.5
-                      * (fineField(fineStartIndex[dirX], fineStartIndex[dirY])
-                         + fineField(fineStartIndex[dirX] + 1, fineStartIndex[dirY]));
-            }
-            else if (centering_[dirX] == core::QtyCentering::dual
-                     and centering_[dirY] == core::QtyCentering::dual)
-            {
-                coarseField(coarseIndex[dirX], coarseIndex[dirY])
-                    = 0.25
-                      * (fineField(fineStartIndex[dirX], fineStartIndex[dirY])
-                         + fineField(fineStartIndex[dirX] + 1, fineStartIndex[dirY])
-                         + fineField(fineStartIndex[dirX], fineStartIndex[dirY] + 1)
-                         + fineField(fineStartIndex[dirX] + 1, fineStartIndex[dirY] + 1));
-            }
-            else
-            {
-                throw std::runtime_error("no magnetic field should end up here");
-            }
-        }
-        else if constexpr (dimension == 3)
-        {
-            throw std::runtime_error("Not Implemented yet");
-        }
-    }
 
-private:
-    std::array<core::QtyCentering, dimension> const centering_;
+    std::array<core::QtyCentering, dim> const centering_;
     SAMRAI::hier::Box const sourceBox_;
     SAMRAI::hier::Box const destinationBox_;
     static int constexpr ratio_ = 2;
 };
+
+template<std::size_t dim>
+template<std::size_t D, typename FieldT>
+typename std::enable_if<D == 1, void>::type
+MagneticFieldCoarsener<dim>::coarsen(Point_t const fineStartIndex, FieldT const& fineField,
+                                     FieldT& coarseField, Point_t const coarseIndex)
+{
+    if (centering_[dirX] == core::QtyCentering::primal) // bx
+    {
+        coarseField(coarseIndex[dirX]) = fineField(fineStartIndex[dirX]);
+    }
+    else if (centering_[dirX] == core::QtyCentering::dual) // by and bz
+    {
+        coarseField(coarseIndex[dirX])
+            = 0.5 * (fineField(fineStartIndex[dirX] + 1) + fineField(fineStartIndex[dirX]));
+    }
+}
+
+template<std::size_t dim>
+template<std::size_t D, typename FieldT>
+typename std::enable_if<D == 2, void>::type
+MagneticFieldCoarsener<dim>::coarsen(Point_t const fineStartIndex, FieldT const& fineField,
+                                     FieldT& coarseField, Point_t const coarseIndex)
+{
+    if (centering_[dirX] == core::QtyCentering::primal
+        and centering_[dirY] == core::QtyCentering::dual)
+    {
+        coarseField(coarseIndex[dirX], coarseIndex[dirY])
+            = 0.5
+              * (fineField(fineStartIndex[dirX], fineStartIndex[dirY])
+                 + fineField(fineStartIndex[dirX], fineStartIndex[dirY] + 1));
+    }
+    else if (centering_[dirX] == core::QtyCentering::dual
+             and centering_[dirY] == core::QtyCentering::primal)
+    {
+        coarseField(coarseIndex[dirX], coarseIndex[dirY])
+            = 0.5
+              * (fineField(fineStartIndex[dirX], fineStartIndex[dirY])
+                 + fineField(fineStartIndex[dirX] + 1, fineStartIndex[dirY]));
+    }
+    else if (centering_[dirX] == core::QtyCentering::dual
+             and centering_[dirY] == core::QtyCentering::dual)
+    {
+        coarseField(coarseIndex[dirX], coarseIndex[dirY])
+            = 0.25
+              * (fineField(fineStartIndex[dirX], fineStartIndex[dirY])
+                 + fineField(fineStartIndex[dirX] + 1, fineStartIndex[dirY])
+                 + fineField(fineStartIndex[dirX], fineStartIndex[dirY] + 1)
+                 + fineField(fineStartIndex[dirX] + 1, fineStartIndex[dirY] + 1));
+    }
+    else
+    {
+        throw std::runtime_error("no magnetic field should end up here");
+    }
+}
+
+// clang-format off
+
+// clang-format on
+
+template<std::size_t dim>
+template<std::size_t D, typename FieldT>
+typename std::enable_if<D == 3, void>::type
+MagneticFieldCoarsener<dim>::coarsen(Point_t const fineStartIndex, FieldT const& fineField,
+                                     FieldT& coarseField, Point_t const coarseIndex)
+{
+    auto& ci  = coarseIndex;
+    auto& cf  = coarseField;
+    auto& fsi = fineStartIndex;
+    auto& ff  = fineField;
+
+    using Fn     = std::function<void()>;
+    using Fn_arr = std::array<Fn, 2>;
+    using A0_arr = std::array<Fn_arr, 2>;
+
+    // X = 0 = primal
+    // Y = 1 = dual
+    std::array<A0_arr, 2> const fns{
+        A0_arr{Fn_arr{[]() { // XXX
+                          throw std::runtime_error("no magnetic field should end up here");
+                      },
+                      [&]() { // XXY
+                          double constexpr static V = .5;
+
+                          cf(ci) = (ff(fsi[dirX], fsi[dirY], fsi[dirZ]) + //
+                                    ff(fsi[dirX], fsi[dirY], fsi[dirZ] + 1))
+                                   * V;
+                      }},
+               Fn_arr{[&]() { // XYX
+                          double constexpr static V = .5;
+
+                          cf(ci) = (ff(fsi[dirX], fsi[dirY], fsi[dirZ]) + //
+                                    ff(fsi[dirX], fsi[dirY] + 1, fsi[dirZ]))
+                                   * V;
+                      },
+                      [&]() { // XYY
+                          double constexpr static V = .5;
+
+                          cf(ci) = (ff(fsi[dirX], fsi[dirY], fsi[dirZ]) + //
+                                    ff(fsi[dirX], fsi[dirY] + 1, fsi[dirZ] + 1))
+                                   * V;
+                      }}},
+        A0_arr{Fn_arr{[&]() {
+                          // YXX
+                          double constexpr static V = .5;
+
+                          cf(ci) = (ff(fsi[dirX], fsi[dirY], fsi[dirZ]) + //
+                                    ff(fsi[dirX] + 1, fsi[dirY], fsi[dirZ]))
+                                   * V;
+                      },
+                      [&]() {
+                          // YXY
+                          double constexpr static V = .5;
+
+                          cf(ci) = (ff(fsi[dirX], fsi[dirY], fsi[dirZ]) + //
+                                    ff(fsi[dirX] + 1, fsi[dirY], fsi[dirZ] + 1))
+                                   * V;
+                      }},
+               Fn_arr{[&]() {
+                          // YYX
+                          double constexpr static V = .5;
+
+                          cf(ci) = (ff(fsi[dirX], fsi[dirY], fsi[dirZ]) + //
+                                    ff(fsi[dirX] + 1, fsi[dirY] + 1, fsi[dirZ]))
+                                   * V;
+                      },
+                      [&]() {
+                          // YYY
+                          double constexpr static V = .125;
+
+                          cf(ci) = (ff(fsi[dirX], fsi[dirY], fsi[dirZ]) +         //
+                                    ff(fsi[dirX] + 1, fsi[dirY], fsi[dirZ]) +     //
+                                    ff(fsi[dirX], fsi[dirY] + 1, fsi[dirZ]) +     //
+                                    ff(fsi[dirX], fsi[dirY], fsi[dirZ] + 1) +     //
+                                    ff(fsi[dirX] + 1, fsi[dirY] + 1, fsi[dirZ]) + //
+                                    ff(fsi[dirX] + 1, fsi[dirY], fsi[dirZ] + 1) + //
+                                    ff(fsi[dirX], fsi[dirY] + 1, fsi[dirZ] + 1) + //
+                                    ff(fsi[dirX] + 1, fsi[dirY] + 1, fsi[dirZ] + 1))
+                                   * V;
+                      }}}};
+
+    auto cast = [](auto const& qty) {
+        return static_cast<std::underlying_type_t<core::QtyCentering>>(qty);
+    };
+    fns[cast(centering_[0])][cast(centering_[1])][cast(centering_[2])]();
+}
+
 } // namespace PHARE::amr
 #endif
diff --git a/src/amr/data/particles/refine/split.hpp b/src/amr/data/particles/refine/split.hpp
index 1d70a7689..73358bf77 100644
--- a/src/amr/data/particles/refine/split.hpp
+++ b/src/amr/data/particles/refine/split.hpp
@@ -3,6 +3,7 @@
 
 #include "amr/data/particles/refine/split_1d.hpp"
 #include "amr/data/particles/refine/split_2d.hpp"
+#include "amr/data/particles/refine/split_3d.hpp"
 
 
 #endif // endif PHARE_SPLIT_HPP
diff --git a/src/amr/data/particles/refine/split_1d.hpp b/src/amr/data/particles/refine/split_1d.hpp
index beef8ff26..c0aade255 100644
--- a/src/amr/data/particles/refine/split_1d.hpp
+++ b/src/amr/data/particles/refine/split_1d.hpp
@@ -17,11 +17,11 @@ using namespace PHARE::core;
 
 /**************************************************************************/
 template<>
-struct PinkDispatcher<DimConst<1>> : SplitPattern<DimConst<1>, RefinedParticlesConst<2>>
+struct PinkPattern<DimConst<1>> : SplitPattern<DimConst<1>, RefinedParticlesConst<2>>
 {
     using Super = SplitPattern<DimConst<1>, RefinedParticlesConst<2>>;
 
-    constexpr PinkDispatcher(float const weight, float const delta)
+    constexpr PinkPattern(float const weight, float const delta)
         : Super{weight}
     {
         Super::deltas_[0] = {delta};
@@ -31,7 +31,7 @@ struct PinkDispatcher<DimConst<1>> : SplitPattern<DimConst<1>, RefinedParticlesC
 
 
 /**************************************************************************/
-using SplitPattern_1_1_2_Dispatcher = PatternDispatcher<PinkDispatcher<DimConst<1>>>;
+using SplitPattern_1_1_2_Dispatcher = PatternDispatcher<PinkPattern<DimConst<1>>>;
 
 template<>
 struct Splitter<DimConst<1>, InterpConst<1>, RefinedParticlesConst<2>>
@@ -50,7 +50,7 @@ struct Splitter<DimConst<1>, InterpConst<1>, RefinedParticlesConst<2>>
 
 /**************************************************************************/
 using SplitPattern_1_1_3_Dispatcher
-    = PatternDispatcher<BlackDispatcher<DimConst<1>>, PinkDispatcher<DimConst<1>>>;
+    = PatternDispatcher<BlackPattern<DimConst<1>>, PinkPattern<DimConst<1>>>;
 
 template<>
 struct Splitter<DimConst<1>, InterpConst<1>, RefinedParticlesConst<3>>
@@ -68,7 +68,7 @@ struct Splitter<DimConst<1>, InterpConst<1>, RefinedParticlesConst<3>>
 
 
 /**************************************************************************/
-using SplitPattern_1_2_2_Dispatcher = PatternDispatcher<PinkDispatcher<DimConst<1>>>;
+using SplitPattern_1_2_2_Dispatcher = PatternDispatcher<PinkPattern<DimConst<1>>>;
 
 template<>
 struct Splitter<DimConst<1>, InterpConst<2>, RefinedParticlesConst<2>>
@@ -87,7 +87,7 @@ struct Splitter<DimConst<1>, InterpConst<2>, RefinedParticlesConst<2>>
 
 /**************************************************************************/
 using SplitPattern_1_2_3_Dispatcher
-    = PatternDispatcher<BlackDispatcher<DimConst<1>>, PinkDispatcher<DimConst<1>>>;
+    = PatternDispatcher<BlackPattern<DimConst<1>>, PinkPattern<DimConst<1>>>;
 
 template<>
 struct Splitter<DimConst<1>, InterpConst<2>, RefinedParticlesConst<3>>
@@ -106,7 +106,7 @@ struct Splitter<DimConst<1>, InterpConst<2>, RefinedParticlesConst<3>>
 
 /**************************************************************************/
 using SplitPattern_1_2_4_Dispatcher
-    = PatternDispatcher<PinkDispatcher<DimConst<1>>, PinkDispatcher<DimConst<1>>>;
+    = PatternDispatcher<PinkPattern<DimConst<1>>, PinkPattern<DimConst<1>>>;
 
 template<>
 struct Splitter<DimConst<1>, InterpConst<2>, RefinedParticlesConst<4>>
@@ -124,7 +124,7 @@ struct Splitter<DimConst<1>, InterpConst<2>, RefinedParticlesConst<4>>
 
 
 /**************************************************************************/
-using SplitPattern_1_3_2_Dispatcher = PatternDispatcher<PinkDispatcher<DimConst<1>>>;
+using SplitPattern_1_3_2_Dispatcher = PatternDispatcher<PinkPattern<DimConst<1>>>;
 
 template<>
 struct Splitter<DimConst<1>, InterpConst<3>, RefinedParticlesConst<2>>
@@ -143,7 +143,7 @@ struct Splitter<DimConst<1>, InterpConst<3>, RefinedParticlesConst<2>>
 
 /**************************************************************************/
 using SplitPattern_1_3_3_Dispatcher
-    = PatternDispatcher<BlackDispatcher<DimConst<1>>, PinkDispatcher<DimConst<1>>>;
+    = PatternDispatcher<BlackPattern<DimConst<1>>, PinkPattern<DimConst<1>>>;
 
 template<>
 struct Splitter<DimConst<1>, InterpConst<3>, RefinedParticlesConst<3>>
@@ -162,7 +162,7 @@ struct Splitter<DimConst<1>, InterpConst<3>, RefinedParticlesConst<3>>
 
 /**************************************************************************/
 using SplitPattern_1_3_4_Dispatcher
-    = PatternDispatcher<PinkDispatcher<DimConst<1>>, PinkDispatcher<DimConst<1>>>;
+    = PatternDispatcher<PinkPattern<DimConst<1>>, PinkPattern<DimConst<1>>>;
 
 template<>
 struct Splitter<DimConst<1>, InterpConst<3>, RefinedParticlesConst<4>>
@@ -181,8 +181,8 @@ struct Splitter<DimConst<1>, InterpConst<3>, RefinedParticlesConst<4>>
 
 /**************************************************************************/
 using SplitPattern_1_3_5_Dispatcher
-    = PatternDispatcher<BlackDispatcher<DimConst<1>>, PinkDispatcher<DimConst<1>>,
-                        PinkDispatcher<DimConst<1>>>;
+    = PatternDispatcher<BlackPattern<DimConst<1>>, PinkPattern<DimConst<1>>,
+                        PinkPattern<DimConst<1>>>;
 
 template<>
 struct Splitter<DimConst<1>, InterpConst<3>, RefinedParticlesConst<5>>
diff --git a/src/amr/data/particles/refine/split_2d.hpp b/src/amr/data/particles/refine/split_2d.hpp
index 53a0bfa3f..f082b98ab 100644
--- a/src/amr/data/particles/refine/split_2d.hpp
+++ b/src/amr/data/particles/refine/split_2d.hpp
@@ -18,11 +18,11 @@ using namespace PHARE::core;
 
 /**************************************************************************/
 template<>
-struct PurpleDispatcher<DimConst<2>> : SplitPattern<DimConst<2>, RefinedParticlesConst<4>>
+struct PurplePattern<DimConst<2>> : SplitPattern<DimConst<2>, RefinedParticlesConst<4>>
 {
     using Super = SplitPattern<DimConst<2>, RefinedParticlesConst<4>>;
 
-    constexpr PurpleDispatcher(float const weight, float const delta)
+    constexpr PurplePattern(float const weight, float const delta)
         : Super{weight}
     {
         Super::deltas_[0] = {-delta, -delta};
@@ -34,12 +34,12 @@ struct PurpleDispatcher<DimConst<2>> : SplitPattern<DimConst<2>, RefinedParticle
 
 
 template<>
-struct BrownDispatcher<DimConst<2>> : SplitPattern<DimConst<2>, RefinedParticlesConst<8>>
+struct BrownPattern<DimConst<2>> : SplitPattern<DimConst<2>, RefinedParticlesConst<8>>
 {
     using Super = SplitPattern<DimConst<2>, RefinedParticlesConst<8>>;
 
     template<typename Delta>
-    constexpr BrownDispatcher(float const weight, Delta const& delta)
+    constexpr BrownPattern(float const weight, Delta const& delta)
         : Super{weight}
     {
         for (std::size_t i = 0; i < 2; i++)
@@ -55,11 +55,11 @@ struct BrownDispatcher<DimConst<2>> : SplitPattern<DimConst<2>, RefinedParticles
 };
 
 template<>
-struct PinkDispatcher<DimConst<2>> : SplitPattern<DimConst<2>, RefinedParticlesConst<4>>
+struct PinkPattern<DimConst<2>> : SplitPattern<DimConst<2>, RefinedParticlesConst<4>>
 {
     using Super = SplitPattern<DimConst<2>, RefinedParticlesConst<4>>;
 
-    constexpr PinkDispatcher(float const weight, float const delta)
+    constexpr PinkPattern(float const weight, float const delta)
         : Super{weight}
     {
         Super::deltas_[0] = {0.0f, -delta};
@@ -71,7 +71,7 @@ struct PinkDispatcher<DimConst<2>> : SplitPattern<DimConst<2>, RefinedParticlesC
 
 
 /**************************************************************************/
-using SplitPattern_2_1_4_Dispatcher = PatternDispatcher<PurpleDispatcher<DimConst<2>>>;
+using SplitPattern_2_1_4_Dispatcher = PatternDispatcher<PurplePattern<DimConst<2>>>;
 
 template<>
 struct Splitter<DimConst<2>, InterpConst<1>, RefinedParticlesConst<4>>
@@ -90,7 +90,7 @@ struct Splitter<DimConst<2>, InterpConst<1>, RefinedParticlesConst<4>>
 
 /**************************************************************************/
 using SplitPattern_2_1_5_Dispatcher
-    = PatternDispatcher<BlackDispatcher<DimConst<2>>, PurpleDispatcher<DimConst<2>>>;
+    = PatternDispatcher<BlackPattern<DimConst<2>>, PurplePattern<DimConst<2>>>;
 
 template<>
 struct Splitter<DimConst<2>, InterpConst<1>, RefinedParticlesConst<5>>
@@ -109,7 +109,7 @@ struct Splitter<DimConst<2>, InterpConst<1>, RefinedParticlesConst<5>>
 
 /**************************************************************************/
 using SplitPattern_2_1_8_Dispatcher
-    = PatternDispatcher<PurpleDispatcher<DimConst<2>>, PinkDispatcher<DimConst<2>>>;
+    = PatternDispatcher<PurplePattern<DimConst<2>>, PinkPattern<DimConst<2>>>;
 
 template<>
 struct Splitter<DimConst<2>, InterpConst<1>, RefinedParticlesConst<8>>
@@ -128,8 +128,8 @@ struct Splitter<DimConst<2>, InterpConst<1>, RefinedParticlesConst<8>>
 
 /**************************************************************************/
 using SplitPattern_2_1_9_Dispatcher
-    = PatternDispatcher<BlackDispatcher<DimConst<2>>, PinkDispatcher<DimConst<2>>,
-                        PurpleDispatcher<DimConst<2>>>;
+    = PatternDispatcher<BlackPattern<DimConst<2>>, PinkPattern<DimConst<2>>,
+                        PurplePattern<DimConst<2>>>;
 
 template<>
 struct Splitter<DimConst<2>, InterpConst<1>, RefinedParticlesConst<9>>
@@ -147,7 +147,7 @@ struct Splitter<DimConst<2>, InterpConst<1>, RefinedParticlesConst<9>>
 
 
 /**************************************************************************/
-using SplitPattern_2_2_4_Dispatcher = PatternDispatcher<PinkDispatcher<DimConst<2>>>;
+using SplitPattern_2_2_4_Dispatcher = PatternDispatcher<PinkPattern<DimConst<2>>>;
 
 template<>
 struct Splitter<DimConst<2>, InterpConst<2>, RefinedParticlesConst<4>>
@@ -166,7 +166,7 @@ struct Splitter<DimConst<2>, InterpConst<2>, RefinedParticlesConst<4>>
 
 /**************************************************************************/
 using SplitPattern_2_2_5_Dispatcher
-    = PatternDispatcher<BlackDispatcher<DimConst<2>>, PinkDispatcher<DimConst<2>>>;
+    = PatternDispatcher<BlackPattern<DimConst<2>>, PinkPattern<DimConst<2>>>;
 
 template<>
 struct Splitter<DimConst<2>, InterpConst<2>, RefinedParticlesConst<5>>
@@ -185,7 +185,7 @@ struct Splitter<DimConst<2>, InterpConst<2>, RefinedParticlesConst<5>>
 
 /**************************************************************************/
 using SplitPattern_2_2_8_Dispatcher
-    = PatternDispatcher<PinkDispatcher<DimConst<2>>, PurpleDispatcher<DimConst<2>>>;
+    = PatternDispatcher<PinkPattern<DimConst<2>>, PurplePattern<DimConst<2>>>;
 
 template<>
 struct Splitter<DimConst<2>, InterpConst<2>, RefinedParticlesConst<8>>
@@ -204,8 +204,8 @@ struct Splitter<DimConst<2>, InterpConst<2>, RefinedParticlesConst<8>>
 
 /**************************************************************************/
 using SplitPattern_2_2_9_Dispatcher
-    = PatternDispatcher<BlackDispatcher<DimConst<2>>, PinkDispatcher<DimConst<2>>,
-                        PurpleDispatcher<DimConst<2>>>;
+    = PatternDispatcher<BlackPattern<DimConst<2>>, PinkPattern<DimConst<2>>,
+                        PurplePattern<DimConst<2>>>;
 
 template<>
 struct Splitter<DimConst<2>, InterpConst<2>, RefinedParticlesConst<9>>
@@ -224,8 +224,8 @@ struct Splitter<DimConst<2>, InterpConst<2>, RefinedParticlesConst<9>>
 
 /**************************************************************************/
 using SplitPattern_2_2_16_Dispatcher
-    = PatternDispatcher<PurpleDispatcher<DimConst<2>>, BrownDispatcher<DimConst<2>>,
-                        PurpleDispatcher<DimConst<2>>>;
+    = PatternDispatcher<PurplePattern<DimConst<2>>, BrownPattern<DimConst<2>>,
+                        PurplePattern<DimConst<2>>>;
 
 template<>
 struct Splitter<DimConst<2>, InterpConst<2>, RefinedParticlesConst<16>>
@@ -244,7 +244,7 @@ struct Splitter<DimConst<2>, InterpConst<2>, RefinedParticlesConst<16>>
 
 
 /**************************************************************************/
-using SplitPattern_2_3_4_Dispatcher = PatternDispatcher<PurpleDispatcher<DimConst<2>>>;
+using SplitPattern_2_3_4_Dispatcher = PatternDispatcher<PurplePattern<DimConst<2>>>;
 
 template<>
 struct Splitter<DimConst<2>, InterpConst<3>, RefinedParticlesConst<4>>
@@ -263,7 +263,7 @@ struct Splitter<DimConst<2>, InterpConst<3>, RefinedParticlesConst<4>>
 
 /**************************************************************************/
 using SplitPattern_2_3_5_Dispatcher
-    = PatternDispatcher<BlackDispatcher<DimConst<2>>, PinkDispatcher<DimConst<2>>>;
+    = PatternDispatcher<BlackPattern<DimConst<2>>, PinkPattern<DimConst<2>>>;
 
 template<>
 struct Splitter<DimConst<2>, InterpConst<3>, RefinedParticlesConst<5>>
@@ -282,7 +282,7 @@ struct Splitter<DimConst<2>, InterpConst<3>, RefinedParticlesConst<5>>
 
 /**************************************************************************/
 using SplitPattern_2_3_8_Dispatcher
-    = PatternDispatcher<PinkDispatcher<DimConst<2>>, PurpleDispatcher<DimConst<2>>>;
+    = PatternDispatcher<PinkPattern<DimConst<2>>, PurplePattern<DimConst<2>>>;
 
 template<>
 struct Splitter<DimConst<2>, InterpConst<3>, RefinedParticlesConst<8>>
@@ -301,8 +301,8 @@ struct Splitter<DimConst<2>, InterpConst<3>, RefinedParticlesConst<8>>
 
 /**************************************************************************/
 using SplitPattern_2_3_9_Dispatcher
-    = PatternDispatcher<BlackDispatcher<DimConst<2>>, PinkDispatcher<DimConst<2>>,
-                        PurpleDispatcher<DimConst<2>>>;
+    = PatternDispatcher<BlackPattern<DimConst<2>>, PinkPattern<DimConst<2>>,
+                        PurplePattern<DimConst<2>>>;
 
 template<>
 struct Splitter<DimConst<2>, InterpConst<3>, RefinedParticlesConst<9>>
@@ -321,9 +321,9 @@ struct Splitter<DimConst<2>, InterpConst<3>, RefinedParticlesConst<9>>
 
 /**************************************************************************/
 using SplitPattern_2_3_25_Dispatcher
-    = PatternDispatcher<BlackDispatcher<DimConst<2>>, PinkDispatcher<DimConst<2>>,
-                        PurpleDispatcher<DimConst<2>>, PinkDispatcher<DimConst<2>>,
-                        BrownDispatcher<DimConst<2>>, PurpleDispatcher<DimConst<2>>>;
+    = PatternDispatcher<BlackPattern<DimConst<2>>, PinkPattern<DimConst<2>>,
+                        PurplePattern<DimConst<2>>, PinkPattern<DimConst<2>>,
+                        BrownPattern<DimConst<2>>, PurplePattern<DimConst<2>>>;
 
 template<>
 struct Splitter<DimConst<2>, InterpConst<3>, RefinedParticlesConst<25>>
diff --git a/src/amr/data/particles/refine/split_3d.hpp b/src/amr/data/particles/refine/split_3d.hpp
new file mode 100644
index 000000000..508907abc
--- /dev/null
+++ b/src/amr/data/particles/refine/split_3d.hpp
@@ -0,0 +1,302 @@
+/*
+Splitting reference material can be found @
+  https://github.com/PHAREHUB/PHARE/wiki/SplitPattern
+*/
+
+#ifndef PHARE_SPLIT_3D_HPP
+#define PHARE_SPLIT_3D_HPP
+
+#include <array>
+#include <cstddef>
+#include "core/utilities/point/point.hpp"
+#include "core/utilities/types.hpp"
+#include "splitter.hpp"
+
+namespace PHARE::amr
+{
+using namespace PHARE::core;
+
+/**************************************************************************/
+template<> // 1 per face - centered
+struct PinkPattern<DimConst<3>> : SplitPattern<DimConst<3>, RefinedParticlesConst<6>>
+{
+    using Super = SplitPattern<DimConst<3>, RefinedParticlesConst<6>>;
+
+    constexpr PinkPattern(float const weight, float const delta)
+        : Super{weight}
+    {
+        constexpr float zero = 0;
+
+        for (std::size_t i = 0; i < 2; i++)
+        {
+            std::size_t offset = i * 3;
+            float sign         = i % 2 ? -1 : 1;
+            auto mode          = delta * sign;
+
+            Super::deltas_[0 + offset] = {mode, zero, zero};
+            Super::deltas_[1 + offset] = {zero, zero, mode};
+            Super::deltas_[2 + offset] = {zero, mode, zero};
+        }
+    }
+};
+
+template<> // 1 per corner
+struct PurplePattern<DimConst<3>> : SplitPattern<DimConst<3>, RefinedParticlesConst<8>>
+{
+    using Super = SplitPattern<DimConst<3>, RefinedParticlesConst<8>>;
+
+    constexpr PurplePattern(float const weight, float const delta)
+        : Super{weight}
+    {
+        for (std::size_t i = 0; i < 2; i++)
+        {
+            std::size_t offset = i * 4;
+            float sign         = i % 2 ? -1 : 1;
+            auto mode          = delta * sign;
+
+            Super::deltas_[0 + offset] = {mode, mode, mode};
+            Super::deltas_[1 + offset] = {mode, mode, -mode};
+            Super::deltas_[2 + offset] = {mode, -mode, mode};
+            Super::deltas_[3 + offset] = {mode, -mode, -mode};
+        }
+    }
+};
+
+
+template<> // 1 per edge - centered
+struct LimePattern<DimConst<3>> : SplitPattern<DimConst<3>, RefinedParticlesConst<12>>
+{
+    using Super = SplitPattern<DimConst<3>, RefinedParticlesConst<12>>;
+
+    constexpr LimePattern(float const weight, float const delta)
+        : Super{weight}
+    {
+        constexpr float zero = 0;
+
+        auto addSquare = [&delta, this](size_t offset, float y) {
+            Super::deltas_[0 + offset] = {zero, y, delta};
+            Super::deltas_[1 + offset] = {zero, y, -delta};
+            Super::deltas_[2 + offset] = {delta, y, zero};
+            Super::deltas_[3 + offset] = {-delta, y, zero};
+        };
+
+        addSquare(0, delta);  // top
+        addSquare(4, -delta); // bottom
+        addSquare(8, 0);      // middle
+    }
+};
+
+
+template<> // 2 per edge - equidistant from center
+class WhitePattern<DimConst<3>> : SplitPattern<DimConst<3>, RefinedParticlesConst<24>>
+{
+    using Super = SplitPattern<DimConst<3>, RefinedParticlesConst<24>>;
+
+    template<typename Delta>
+    constexpr WhitePattern(float const weight, Delta const& delta)
+        : Super{weight}
+    {
+        auto addSquare = [&delta, this](size_t offset, float x, float y, float z) {
+            Super::deltas_[0 + offset] = {x, y, z};
+            Super::deltas_[1 + offset] = {x, -y, z};
+            Super::deltas_[2 + offset] = {-x, y, z};
+            Super::deltas_[3 + offset] = {-x, -y, z};
+        };
+
+        auto addSquares = [addSquare, &delta, this](size_t offset, float x, float y, float z) {
+            addSquare(offset, x, y, z);
+            addSquare(offset + 4, x, y, -z);
+        };
+
+        addSquares(0, delta[0], delta[0], delta[1]);
+        addSquares(8, delta[0], delta[1], delta[0]);
+        addSquares(16, delta[1], delta[0], delta[0]);
+    }
+};
+
+
+/**************************************************************************/
+/****************************** INTERP == 1 *******************************/
+/**************************************************************************/
+using SplitPattern_3_1_6_Dispatcher = PatternDispatcher<PinkPattern<DimConst<3>>>;
+
+template<>
+struct Splitter<DimConst<3>, InterpConst<1>, RefinedParticlesConst<6>>
+    : public ASplitter<DimConst<3>, InterpConst<1>, RefinedParticlesConst<6>>,
+      SplitPattern_3_1_6_Dispatcher
+{
+    constexpr Splitter()
+        : SplitPattern_3_1_6_Dispatcher{{weight[0], delta[0]}}
+    {
+    }
+
+    static constexpr std::array<float, 1> delta  = {.966431};
+    static constexpr std::array<float, 1> weight = {.166666};
+};
+
+
+/**************************************************************************/
+using SplitPattern_3_1_12_Dispatcher = PatternDispatcher<LimePattern<DimConst<3>>>;
+
+template<>
+struct Splitter<DimConst<3>, InterpConst<1>, RefinedParticlesConst<12>>
+    : public ASplitter<DimConst<3>, InterpConst<1>, RefinedParticlesConst<12>>,
+      SplitPattern_3_1_12_Dispatcher
+{
+    constexpr Splitter()
+        : SplitPattern_3_1_12_Dispatcher{{weight[0], delta[0]}}
+    {
+    }
+
+    static constexpr std::array<float, 1> delta  = {.74823};
+    static constexpr std::array<float, 1> weight = {.083333};
+};
+
+
+/**************************************************************************/
+using SplitPattern_3_1_27_Dispatcher
+    = PatternDispatcher<BlackPattern<DimConst<3>>, PinkPattern<DimConst<3>>,
+                        PurplePattern<DimConst<3>>, LimePattern<DimConst<3>>>;
+
+template<>
+struct Splitter<DimConst<3>, InterpConst<1>, RefinedParticlesConst<27>>
+    : public ASplitter<DimConst<3>, InterpConst<1>, RefinedParticlesConst<27>>,
+      SplitPattern_3_1_27_Dispatcher
+{
+    constexpr Splitter()
+        : SplitPattern_3_1_27_Dispatcher{
+            {weight[0]}, {weight[1], delta[0]}, {weight[1], delta[0]}, {weight[1], delta[0]}}
+    {
+    }
+
+    static constexpr std::array<float, 1> delta  = {1};
+    static constexpr std::array<float, 2> weight = {0.125, 0.0625};
+};
+
+
+
+/**************************************************************************/
+/****************************** INTERP == 2 *******************************/
+/**************************************************************************/
+using SplitPattern_3_2_6_Dispatcher = PatternDispatcher<PinkPattern<DimConst<3>>>;
+
+template<>
+struct Splitter<DimConst<3>, InterpConst<2>, RefinedParticlesConst<6>>
+    : public ASplitter<DimConst<3>, InterpConst<2>, RefinedParticlesConst<6>>,
+      SplitPattern_3_2_6_Dispatcher
+{
+    constexpr Splitter()
+        : SplitPattern_3_2_6_Dispatcher{{weight[0], delta[0]}}
+    {
+    }
+
+    static constexpr std::array<float, 1> delta  = {1.149658};
+    static constexpr std::array<float, 1> weight = {.166666};
+};
+
+
+/**************************************************************************/
+using SplitPattern_3_2_12_Dispatcher = PatternDispatcher<LimePattern<DimConst<3>>>;
+
+template<>
+struct Splitter<DimConst<3>, InterpConst<2>, RefinedParticlesConst<12>>
+    : public ASplitter<DimConst<3>, InterpConst<2>, RefinedParticlesConst<12>>,
+      SplitPattern_3_2_12_Dispatcher
+{
+    constexpr Splitter()
+        : SplitPattern_3_2_12_Dispatcher{{weight[0], delta[0]}}
+    {
+    }
+
+    static constexpr std::array<float, 1> delta  = {.888184};
+    static constexpr std::array<float, 1> weight = {.083333};
+};
+
+
+/**************************************************************************/
+using SplitPattern_3_2_27_Dispatcher
+    = PatternDispatcher<BlackPattern<DimConst<3>>, PinkPattern<DimConst<3>>,
+                        PurplePattern<DimConst<3>>, LimePattern<DimConst<3>>>;
+
+template<>
+struct Splitter<DimConst<3>, InterpConst<2>, RefinedParticlesConst<27>>
+    : public ASplitter<DimConst<3>, InterpConst<2>, RefinedParticlesConst<27>>,
+      SplitPattern_3_2_27_Dispatcher
+{
+    constexpr Splitter()
+        : SplitPattern_3_2_27_Dispatcher{
+            {weight[0]}, {weight[1], delta[0]}, {weight[1], delta[0]}, {weight[1], delta[0]}}
+    {
+    }
+
+    static constexpr std::array<float, 1> delta  = {1.111333};
+    static constexpr std::array<float, 2> weight = {.099995, .055301};
+};
+
+
+/**************************************************************************/
+/****************************** INTERP == 3 *******************************/
+/**************************************************************************/
+using SplitPattern_3_3_6_Dispatcher = PatternDispatcher<PinkPattern<DimConst<3>>>;
+
+template<>
+struct Splitter<DimConst<3>, InterpConst<3>, RefinedParticlesConst<6>>
+    : public ASplitter<DimConst<3>, InterpConst<3>, RefinedParticlesConst<6>>,
+      SplitPattern_3_3_6_Dispatcher
+{
+    constexpr Splitter()
+        : SplitPattern_3_3_6_Dispatcher{{weight[0], delta[0]}}
+    {
+    }
+
+    static constexpr std::array<float, 1> delta  = {1.312622};
+    static constexpr std::array<float, 1> weight = {.166666};
+};
+
+
+/**************************************************************************/
+using SplitPattern_3_3_12_Dispatcher = PatternDispatcher<LimePattern<DimConst<3>>>;
+
+template<>
+struct Splitter<DimConst<3>, InterpConst<3>, RefinedParticlesConst<12>>
+    : public ASplitter<DimConst<3>, InterpConst<3>, RefinedParticlesConst<12>>,
+      SplitPattern_3_3_12_Dispatcher
+{
+    constexpr Splitter()
+        : SplitPattern_3_3_12_Dispatcher{{weight[0], delta[0]}}
+    {
+    }
+
+    static constexpr std::array<float, 1> delta  = {1.012756};
+    static constexpr std::array<float, 1> weight = {.083333};
+};
+
+
+/**************************************************************************/
+using SplitPattern_3_3_27_Dispatcher
+    = PatternDispatcher<BlackPattern<DimConst<3>>, PinkPattern<DimConst<3>>,
+                        PurplePattern<DimConst<3>>, LimePattern<DimConst<3>>>;
+
+template<>
+struct Splitter<DimConst<3>, InterpConst<3>, RefinedParticlesConst<27>>
+    : public ASplitter<DimConst<3>, InterpConst<3>, RefinedParticlesConst<27>>,
+      SplitPattern_3_3_27_Dispatcher
+{
+    constexpr Splitter()
+        : SplitPattern_3_3_27_Dispatcher{
+            {weight[0]}, {weight[1], delta[0]}, {weight[1], delta[0]}, {weight[1], delta[0]}}
+    {
+    }
+
+    static constexpr std::array<float, 1> delta  = {1.276815};
+    static constexpr std::array<float, 2> weight = {.104047, .05564};
+};
+
+
+/**************************************************************************/
+
+
+} // namespace PHARE::amr
+
+
+#endif /* PHARE_SPLIT_3D_HPP */
diff --git a/src/amr/data/particles/refine/splitter.hpp b/src/amr/data/particles/refine/splitter.hpp
index accd4d387..d2bd416a8 100644
--- a/src/amr/data/particles/refine/splitter.hpp
+++ b/src/amr/data/particles/refine/splitter.hpp
@@ -67,6 +67,7 @@ class PatternDispatcher
         using FineParticle = decltype(particles[0]); // may be a reference
 
         core::apply(patterns, [&](auto const& pattern) {
+            auto weight = static_cast<Weight_t>(pattern.weight_);
             for (size_t rpIndex = 0; rpIndex < pattern.deltas_.size(); rpIndex++)
             {
                 FineParticle fineParticle = particles[idx++];
@@ -114,25 +115,33 @@ class Splitter : public ASplitter<dimension, interp_order, nbRefinedPart>
 };
 
 template<typename dim>
-struct BlackDispatcher : SplitPattern<dim, core::RefinedParticlesConst<1>>
+struct BlackPattern : SplitPattern<dim, core::RefinedParticlesConst<1>>
 {
     using Super = SplitPattern<dim, core::RefinedParticlesConst<1>>;
-    constexpr BlackDispatcher(float const weight)
+    constexpr BlackPattern(float const weight)
         : Super{weight}
     {
     }
 };
 
 template<typename dim>
-struct PurpleDispatcher
+struct PinkPattern
 {
 };
 template<typename dim>
-struct BrownDispatcher
+struct PurplePattern
 {
 };
 template<typename dim>
-struct PinkDispatcher
+struct BrownPattern
+{
+};
+template<typename dim>
+struct LimePattern
+{
+};
+template<typename dim>
+struct WhitePattern
 {
 };
 
diff --git a/src/core/data/grid/gridlayoutdefs.hpp b/src/core/data/grid/gridlayoutdefs.hpp
index 23920396e..f5cef3bba 100644
--- a/src/core/data/grid/gridlayoutdefs.hpp
+++ b/src/core/data/grid/gridlayoutdefs.hpp
@@ -2,6 +2,7 @@
 #define PHARE_CORE_GRID_GRIDLAYOUTDEFS_HPP
 
 #include <cstddef>
+#include <cstdint>
 
 #include "core/hybrid/hybrid_quantities.hpp"
 #include "core/utilities/types.hpp"
@@ -14,7 +15,7 @@ namespace core
     enum class Direction { X, Y, Z };
 
 
-    enum class QtyCentering { primal = 0, dual = 1 };
+    enum class QtyCentering : std::uint16_t { primal = 0, dual = 1 };
 
 
     template<std::size_t dim>
diff --git a/src/core/data/ndarray/ndarray_vector.hpp b/src/core/data/ndarray/ndarray_vector.hpp
index 7d52d1d01..98898b796 100644
--- a/src/core/data/ndarray/ndarray_vector.hpp
+++ b/src/core/data/ndarray/ndarray_vector.hpp
@@ -3,6 +3,7 @@
 
 #include "core/def.hpp"
 #include <stdexcept>
+#include <algorithm>
 #include <array>
 #include <cstdint>
 #include <vector>
@@ -123,6 +124,10 @@ class MaskedView
     NO_DISCARD auto yend() const { return shape_[1] - 1 - mask_.max(); }
 
 
+    auto zstart() const { return mask_.min(); }
+    auto zend() const { return shape_[2] - 1 - mask_.max(); }
+
+
 private:
     Array& array_;
     std::array<std::uint32_t, dimension> shape_;
@@ -344,10 +349,10 @@ class NdArrayMask
     {
         auto shape = array.shape();
 
-        for (std::size_t i = min_; i <= max_; ++i)
+        for (auto i = min_; i <= max_; ++i)
             array(i) = val;
 
-        for (std::size_t i = shape[0] - 1 - max_; i <= shape[0] - 1 - min_; ++i)
+        for (auto i = shape[0] - 1 - max_; i <= shape[0] - 1 - min_; ++i)
             array(i) = val;
     }
 
@@ -357,24 +362,24 @@ class NdArrayMask
         auto shape = array.shape();
 
         // left border
-        for (std::size_t i = min_; i <= max_; ++i)
-            for (std::size_t j = min_; j <= shape[1] - 1 - max_; ++j)
+        for (auto i = min_; i <= max_; ++i)
+            for (auto j = min_; j <= shape[1] - 1 - max_; ++j)
                 array(i, j) = val;
 
         // right border
-        for (std::size_t i = shape[0] - 1 - max_; i <= shape[0] - 1 - min_; ++i)
-            for (std::size_t j = min_; j <= shape[1] - 1 - max_; ++j)
+        for (auto i = shape[0] - 1 - max_; i <= shape[0] - 1 - min_; ++i)
+            for (auto j = min_; j <= shape[1] - 1 - max_; ++j)
                 array(i, j) = val;
 
 
-        for (std::size_t i = min_; i <= shape[0] - 1 - min_; ++i)
+        for (auto i = min_; i <= shape[0] - 1 - min_; ++i)
         {
             // bottom border
-            for (std::size_t j = min_; j <= max_; ++j)
+            for (auto j = min_; j <= max_; ++j)
                 array(i, j) = val;
 
             // top border
-            for (std::size_t j = shape[1] - 1 - max_; j <= shape[1] - 1 - min_; ++j)
+            for (auto j = shape[1] - 1 - max_; j <= shape[1] - 1 - min_; ++j)
                 array(i, j) = val;
         }
     }
@@ -382,7 +387,44 @@ class NdArrayMask
     template<typename Array>
     void fill3D(Array& array, typename Array::type val) const
     {
-        throw std::runtime_error("3d not implemented");
+        auto shape = array.shape();
+
+        // left border
+        for (auto i = min_; i <= shape[0] - 1 - max_; ++i)
+            for (auto j = min_; j <= shape[1] - 1 - max_; ++j)
+                for (auto k = min_; k <= max_; ++k)
+                    array(i, j, k) = val;
+
+        // // right border
+        for (auto i = min_; i <= shape[0] - 1 - max_; ++i)
+            for (auto j = min_; j <= shape[1] - 1 - max_; ++j)
+                for (auto k = shape[2] - 1 - max_; k <= shape[2] - 1 - min_; ++k)
+                    array(i, j, k) = val;
+
+        for (auto i = min_; i <= shape[0] - 1 - min_; ++i)
+        {
+            // bottom border
+            for (auto j = min_; j <= max_; ++j)
+                for (auto k = min_; k <= shape[2] - 1 - min_; ++k)
+                    array(i, j, k) = val;
+
+            // top border
+            for (auto j = shape[1] - 1 - max_; j <= shape[1] - 1 - min_; ++j)
+                for (auto k = min_; k <= shape[2] - 1 - min_; ++k)
+                    array(i, j, k) = val;
+        }
+
+        // front
+        for (auto i = min_; i <= max_; ++i)
+            for (auto j = min_; j <= shape[1] - 1 - max_; ++j)
+                for (auto k = min_; k <= shape[2] - 1 - min_; ++k)
+                    array(i, j, k) = val;
+
+        // back
+        for (auto i = shape[0] - 1 - max_; i <= shape[0] - 1 - min_; ++i)
+            for (auto j = min_; j <= shape[1] - 1 - max_; ++j)
+                for (auto k = min_; k <= shape[2] - 1 - min_; ++k)
+                    array(i, j, k) = val;
     }
 
     template<typename Array>
@@ -392,17 +434,23 @@ class NdArrayMask
 
         std::size_t cells = 0;
 
-        if constexpr (Array::dimension == 1)
-            for (std::size_t i = min_; i <= max_; ++i)
+        for (auto i = min_; i <= max_; ++i)
+        {
+            if constexpr (Array::dimension == 1)
                 cells += 2;
 
-        if constexpr (Array::dimension == 2)
-            for (std::size_t i = min_; i <= max_; ++i)
+            if constexpr (Array::dimension == 2)
                 cells += (shape[0] - (i * 2) - 2) * 2 + (shape[1] - (i * 2) - 2) * 2 + 4;
 
-        if constexpr (Array::dimension == 3)
-            throw std::runtime_error("Not implemented dimension");
-
+            if constexpr (Array::dimension == 3)
+            {
+                auto [x, y, z] = shape;
+                x -= i * 2;
+                y -= i * 2;
+                z -= i * 2;
+                cells += (x * y * 2) + (y * (z - 2) * 2) + ((z - 2) * (x - 2) * 2);
+            }
+        }
         return cells;
     }
 
@@ -422,20 +470,21 @@ void operator>>(MaskedView<Array, Mask>&& inner, MaskedView<Array, Mask>&& outer
 {
     using MaskedView_t = MaskedView<Array, Mask>;
 
+    assert(inner.xstart() > outer.xstart() and inner.xend() < outer.xend());
+
     if constexpr (MaskedView_t::dimension == 1)
     {
-        assert(inner.xstart() > outer.xstart());
-        assert(inner.xend() < outer.xend());
         outer(outer.xstart()) = inner(inner.xstart());
         outer(outer.xend())   = inner(inner.xend());
     }
 
+    if constexpr (MaskedView_t::dimension > 1)
+    {
+        assert(inner.ystart() > outer.ystart() and inner.yend() < outer.yend());
+    }
 
     if constexpr (MaskedView_t::dimension == 2)
     {
-        assert(inner.xstart() > outer.xstart() and inner.xend() < outer.xend()
-               and inner.ystart() > outer.ystart() and inner.yend() < outer.yend());
-
         for (auto ix = inner.xstart(); ix <= inner.xend(); ++ix)
         {
             outer(ix, outer.ystart()) = inner(ix, inner.ystart()); // bottom
@@ -452,7 +501,7 @@ void operator>>(MaskedView<Array, Mask>&& inner, MaskedView<Array, Mask>&& outer
         for (auto ix = outer.xstart(); ix < inner.xstart(); ++ix)
             outer(ix, outer.ystart()) = inner(inner.xstart(), inner.ystart());
 
-        for (std::size_t iy = outer.ystart(); iy < inner.ystart(); ++iy)
+        for (auto iy = outer.ystart(); iy < inner.ystart(); ++iy)
             outer(outer.xstart(), iy) = inner(inner.xstart(), inner.ystart());
 
 
@@ -481,7 +530,72 @@ void operator>>(MaskedView<Array, Mask>&& inner, MaskedView<Array, Mask>&& outer
 
     if constexpr (MaskedView_t::dimension == 3)
     {
-        throw std::runtime_error("3d not implemented");
+        assert(inner.zstart() > outer.zstart() and inner.zend() < outer.zend());
+
+        for (auto i = inner.xstart(); i <= inner.xend(); ++i)
+        {
+            for (auto k = inner.zstart(); k <= inner.zend(); ++k)
+            {
+                outer(i, outer.ystart(), k) = inner(i, inner.ystart(), k);
+                outer(i, outer.yend(), k)   = inner(i, inner.yend(), k);
+            }
+            for (auto j = inner.ystart(); j <= inner.yend(); ++j)
+            {
+                outer(i, j, outer.zstart()) = inner(i, j, inner.zstart());
+                outer(i, j, outer.zend())   = inner(i, j, inner.zend());
+            }
+        }
+
+        for (auto j = inner.ystart(); j <= inner.yend(); ++j)
+        {
+            for (auto k = inner.zstart(); k <= inner.zend(); ++k)
+            {
+                outer(outer.xstart(), j, k) = inner(inner.xstart(), j, k);
+                outer(outer.xend(), j, k)   = inner(inner.xend(), j, k);
+            }
+        }
+
+        for (auto k = inner.zstart(); k <= inner.zend(); ++k)
+        {
+            outer(outer.xstart(), outer.ystart(), k) = inner(outer.xstart(), inner.ystart(), k);
+            outer(outer.xstart(), outer.yend(), k)   = inner(outer.xstart(), inner.yend(), k);
+
+            outer(outer.xend(), outer.ystart(), k) = inner(outer.xend(), inner.ystart(), k);
+            outer(outer.xend(), outer.yend(), k)   = inner(outer.xend(), inner.yend(), k);
+        }
+
+        for (auto j = inner.ystart(); j <= inner.yend(); ++j)
+        {
+            outer(outer.xstart(), j, outer.zstart()) = inner(outer.xstart(), j, inner.zstart());
+            outer(outer.xstart(), j, outer.zend())   = inner(outer.xstart(), j, inner.zend());
+
+            outer(outer.xend(), j, outer.zstart()) = inner(outer.xend(), j, inner.zstart());
+            outer(outer.xend(), j, outer.zend())   = inner(outer.xend(), j, inner.zend());
+        }
+
+        for (auto i = inner.xstart(); i <= inner.xend(); ++i)
+        {
+            outer(i, outer.ystart(), outer.zstart()) = inner(i, outer.ystart(), inner.zstart());
+            outer(i, outer.ystart(), outer.zend())   = inner(i, outer.ystart(), inner.zend());
+
+            outer(i, outer.yend(), outer.zstart()) = inner(i, outer.yend(), inner.zstart());
+            outer(i, outer.yend(), outer.zend())   = inner(i, outer.yend(), inner.zend());
+        }
+
+        auto corner = [&](auto xouter, auto xinner) {
+            outer(xouter, outer.ystart(), outer.zstart())
+                = inner(xinner, outer.ystart(), inner.zstart());
+
+            outer(xouter, outer.ystart(), outer.zend())
+                = inner(xinner, outer.ystart(), inner.zend());
+
+            outer(xouter, outer.yend(), outer.zstart())
+                = inner(xinner, outer.yend(), inner.zstart());
+            outer(xouter, outer.yend(), outer.zend()) = inner(xinner, outer.yend(), inner.zend());
+        };
+
+        corner(outer.xstart(), inner.xstart());
+        corner(outer.xend(), inner.xend());
     }
 }
 
diff --git a/src/core/utilities/box/box.hpp b/src/core/utilities/box/box.hpp
index 0ff5d1a8a..493cf307a 100644
--- a/src/core/utilities/box/box.hpp
+++ b/src/core/utilities/box/box.hpp
@@ -7,6 +7,7 @@
 #include "core/utilities/meta/meta_utilities.hpp"
 #include "core/def.hpp"
 
+#include <vector>
 #include <cstddef>
 #include <algorithm>
 #include <optional>
@@ -146,6 +147,26 @@ struct Box
         else
             return 6;
     }
+
+    auto as_points() const
+    {
+        std::vector<Point<Type, dim>> points;
+        if constexpr (dim == 1)
+            for (Type i = lower[0]; i <= upper[0]; ++i)
+                points.emplace_back(Point{i});
+
+        else if constexpr (dim == 2)
+            for (Type i = lower[0]; i <= upper[0]; ++i)
+                for (Type j = lower[1]; j <= upper[1]; ++j)
+                    points.emplace_back(Point{i, j});
+        else
+            for (Type i = lower[0]; i <= upper[0]; ++i)
+                for (Type j = lower[1]; j <= upper[1]; ++j)
+                    for (Type k = lower[2]; k <= upper[2]; ++k)
+                        points.emplace_back(Point{i, j, k});
+
+        return points;
+    }
 };
 
 template<typename Type, std::size_t dim>
@@ -194,7 +215,7 @@ class box_iterator
 
 
 template<typename T, std::size_t s>
-Box(Point<T, s> lower, Point<T, s> upper) -> Box<T, s>;
+Box(Point<T, s> lower, Point<T, s> upper)->Box<T, s>;
 
 
 
diff --git a/src/core/utilities/meta/meta_utilities.hpp b/src/core/utilities/meta/meta_utilities.hpp
index 3bea76ba3..00c13dbc3 100644
--- a/src/core/utilities/meta/meta_utilities.hpp
+++ b/src/core/utilities/meta/meta_utilities.hpp
@@ -6,6 +6,10 @@
 
 #include "core/utilities/types.hpp"
 
+#if !defined(PHARE_SIMULATORS)
+#define PHARE_SIMULATORS 3
+#endif
+
 namespace PHARE
 {
 namespace core
@@ -74,11 +78,24 @@ namespace core
         // inner tuple = dim, interp, list[possible nbrParticles for dim/interp]
         return std::tuple<SimulatorOption<DimConst<1>, InterpConst<1>, 2, 3>,
                           SimulatorOption<DimConst<1>, InterpConst<2>, 2, 3, 4>,
-                          SimulatorOption<DimConst<1>, InterpConst<3>, 2, 3, 4, 5>,
+                          SimulatorOption<DimConst<1>, InterpConst<3>, 2, 3, 4, 5>
 
+#if PHARE_SIMULATORS > 1
+                          ,
                           SimulatorOption<DimConst<2>, InterpConst<1>, 4, 5, 8, 9>,
                           SimulatorOption<DimConst<2>, InterpConst<2>, 4, 5, 8, 9, 16>,
-                          SimulatorOption<DimConst<2>, InterpConst<3>, 4, 5, 8, 9, 25>>{};
+                          SimulatorOption<DimConst<2>, InterpConst<3>, 4, 5, 8, 9, 25>
+#endif
+
+                          // TODO add in the rest of 3d nbrParticles permutations
+                          // possibly consider compile time activation for uncommon cases
+#if PHARE_SIMULATORS > 2
+                          ,
+                          SimulatorOption<DimConst<3>, InterpConst<1>, 6, 12 /*, 27*/>,
+                          SimulatorOption<DimConst<3>, InterpConst<2>, 6, 12>,
+                          SimulatorOption<DimConst<3>, InterpConst<3>, 6, 12>
+#endif
+>{};
     }
 
 
diff --git a/tests/amr/data/field/refine/test_refine_field.py b/tests/amr/data/field/refine/test_refine_field.py
index 188f80191..a3e66f756 100644
--- a/tests/amr/data/field/refine/test_refine_field.py
+++ b/tests/amr/data/field/refine/test_refine_field.py
@@ -102,6 +102,9 @@ def refine_electric(field, **kwargs):
             fine_data[::2, 1::2] = 0.5 * (coarse_data[:, 1:] + coarse_data[:, :-1])
             fine_data[1::2, 1::2] = 0.5 * (coarse_data[:, 1:] + coarse_data[:, :-1])
 
+    elif field.box.ndim == 3:
+        assert False # fix
+
     return cropToFieldData(fine_data, field)
 
 
diff --git a/tests/amr/data/particles/refine/input/input_3d_ratio_2.txt b/tests/amr/data/particles/refine/input/input_3d_ratio_2.txt
new file mode 100644
index 000000000..064cf8fb5
--- /dev/null
+++ b/tests/amr/data/particles/refine/input/input_3d_ratio_2.txt
@@ -0,0 +1,51 @@
+CartesianGridGeometry
+{
+    domain_boxes = [ (0, 0, 0) , (8, 8, 8) ]
+    x_lo = 0.0, 0.0, 0.0
+    x_up = 1.0, 1.0, 1.0
+    periodic_dimension = 1, 1, 1
+}
+Main
+{
+    dim = 1
+}
+PatchHierarchy
+{
+    max_levels  =   2
+    proper_nesting_buffer = 2
+    // vector of coarse ratio with dim dimension
+    ratio_to_coarser
+    {
+        level_1 =   2, 2, 2
+    }
+    smallest_patch_size
+    {
+        level_0 =   5, 5, 5
+        // All finer level will use same values in as level_0
+     }
+}
+ChopAndPackLoadBalancer
+{
+    bin_pack_method = "SPATIAL"
+}
+StandardTagAndInitialize
+{
+    at_0
+    {
+        cycle = 0
+        tag_0
+        {
+            tagging_method = "REFINE_BOXES"
+            level_0
+            {
+                boxes = [ (1, 1, 1) , (3, 3, 3)]
+            }
+        }
+    }
+}
+TileClustering
+{
+}
+GriddingAlgorithm
+{
+}
diff --git a/tests/amr/data/particles/refine/test_split.cpp b/tests/amr/data/particles/refine/test_split.cpp
index 8027c35a2..77651dd75 100644
--- a/tests/amr/data/particles/refine/test_split.cpp
+++ b/tests/amr/data/particles/refine/test_split.cpp
@@ -22,7 +22,7 @@ struct SplitterTest : public ::testing::Test
     SplitterTest() { Splitter splitter; }
 };
 
-using Splitters = testing::Types<Splitter<1, 1, 2>, Splitter<2, 1, 8> /*, Splitter<3, 1, 27>*/>;
+using Splitters = testing::Types<Splitter<1, 1, 2>, Splitter<2, 1, 8>, Splitter<3, 1, 27>>;
 
 TYPED_TEST_SUITE(SplitterTest, Splitters);
 
diff --git a/tests/core/data/ndarray/test_main.cpp b/tests/core/data/ndarray/test_main.cpp
index be0bebd9f..b034d1b08 100644
--- a/tests/core/data/ndarray/test_main.cpp
+++ b/tests/core/data/ndarray/test_main.cpp
@@ -392,6 +392,60 @@ TEST(MaskedView2d, maskOps2)
                                                                   + Mask{0u}.nCells(array));
 }
 
+TEST(MaskedView3d, maskOps3)
+{
+    constexpr std::size_t dim      = 3;
+    constexpr std::uint32_t size0  = 10;
+    constexpr std::uint32_t sizeCu = size0 * size0 * size0;
+    using Mask                     = PHARE::core::NdArrayMask;
+
+    auto sum = [](auto const& array) { return std::accumulate(array.begin(), array.end(), 0); };
+
+    {
+        NdArrayVector<dim> array{size0, size0, size0};
+        EXPECT_EQ(sum(array), 0);
+        std::fill(array.begin(), array.end(), 1);
+        EXPECT_EQ(sum(array), sizeCu);
+    }
+
+    {
+        NdArrayVector<dim> array{size0, size0, size0};
+        EXPECT_EQ(std::accumulate(array.begin(), array.end(), 0), 0);
+        array[Mask{0}] = 1;
+        EXPECT_EQ(std::accumulate(array.begin(), array.end(), 0), 488);
+
+        // outter cells of a 10**3 cube =
+        // (10 * 10 * 2) + (10 * 8 * 2) + (8 * 8 * 2);
+        // or
+        // (8 * 8 * 6) + (10 * 4) + (8 * 8);
+        // = 488
+    }
+
+    std::uint32_t ten = 10;
+    PHARE::core::NdArrayVector<3> array(ten, ten, ten);
+
+    array[Mask{0}] = 1;
+    EXPECT_EQ(sum(array), 488);
+    array[Mask{1}] >> array[Mask{0}];
+    EXPECT_EQ(sum(array), 0);
+
+    array[Mask{2}] = 1;
+    EXPECT_EQ(sum(array), 152);
+    array[Mask{1}] = 1;
+    EXPECT_EQ(sum(array), 448);
+    array[Mask{1}] = 0;
+    EXPECT_EQ(sum(array), 152);
+
+    array[Mask{2}] >> array[Mask{1}];
+    EXPECT_EQ(sum(array), 448);
+    array[Mask{2}] = 0;
+    EXPECT_EQ(sum(array), 296);
+
+    EXPECT_EQ(Mask{1}.nCells(array), 296);
+    EXPECT_EQ(Mask{2}.nCells(array), 152);
+}
+
+
 int main(int argc, char** argv)
 {
     ::testing::InitGoogleTest(&argc, argv);
diff --git a/tests/core/numerics/interpolator/test_main.cpp b/tests/core/numerics/interpolator/test_main.cpp
index 28e2fff14..b69025e7a 100644
--- a/tests/core/numerics/interpolator/test_main.cpp
+++ b/tests/core/numerics/interpolator/test_main.cpp
@@ -721,6 +721,73 @@ using My2dTypes = ::testing::Types<Interpolator<2, 1>, Interpolator<2, 2>, Inter
 INSTANTIATE_TYPED_TEST_SUITE_P(testInterpolator, ACollectionOfParticles_2d, My2dTypes);
 
 
+
+/*********************************************************************************************/
+template<typename Interpolator>
+struct ACollectionOfParticles_3d : public ::testing::Test
+{
+    static constexpr auto interp_order = Interpolator::interp_order;
+    static constexpr std::size_t dim   = 3;
+    static constexpr std::uint32_t nx = 15, ny = 15, nz = 15;
+    static constexpr int start = 0, end = 5;
+    constexpr static auto safeLayer = static_cast<int>(1 + ghostWidthForParticles<interp_order>());
+
+    using PHARE_TYPES      = PHARE::core::PHARE_Types<dim, interp_order>;
+    using ParticleArray_t  = typename PHARE_TYPES::ParticleArray_t;
+    using GridLayout_t     = typename PHARE_TYPES::GridLayout_t;
+    using Grid_t           = typename PHARE_TYPES::Grid_t;
+    using UsableVecFieldND = UsableVecField<dim>;
+
+    GridLayout_t layout{ConstArray<double, dim>(.1), {nx, ny}, ConstArray<double, dim>(0)};
+    ParticleArray_t particles;
+    Grid_t rho;
+    UsableVecFieldND v;
+    Interpolator interpolator;
+
+    ACollectionOfParticles_3d()
+        : particles{grow(layout.AMRBox(), safeLayer)}
+        , rho{"field", HybridQuantity::Scalar::rho, nx, ny, nz}
+        , v{"v", layout, HybridQuantity::Vector::V}
+    {
+        double weight = [](auto const& meshSize) {
+            return std::accumulate(meshSize.begin(), meshSize.end(), 1.0,
+                                   std::multiplies<double>());
+        }(layout.meshSize());
+
+        for (int i = start; i < end; i++)
+            for (int j = start; j < end; j++)
+                for (int k = start; k < end; k++)
+                {
+                    auto& part  = particles.emplace_back();
+                    part.iCell  = {i, j, k};
+                    part.delta  = ConstArray<double, dim>(.5);
+                    part.weight = weight;
+                    part.v[0]   = +2.;
+                    part.v[1]   = -1.;
+                    part.v[2]   = +1.;
+                }
+
+        interpolator(particles, rho, v, layout);
+    }
+};
+TYPED_TEST_SUITE_P(ACollectionOfParticles_3d);
+
+
+TYPED_TEST_P(ACollectionOfParticles_3d, DepositCorrectlyTheirWeight_3d)
+{
+    // auto const& [vx, vy, vz] = this->v();
+    // EXPECT_DOUBLE_EQ(this->rho(7, 7, 7), 1.0);
+    // EXPECT_DOUBLE_EQ(vx(7, 7, 7), 2.0);
+    // EXPECT_DOUBLE_EQ(vy(7, 7, 7), -1.0);
+    // EXPECT_DOUBLE_EQ(vz(7, 7, 7), 1.0);
+}
+REGISTER_TYPED_TEST_SUITE_P(ACollectionOfParticles_3d, DepositCorrectlyTheirWeight_3d);
+
+
+using My3dTypes = ::testing::Types<Interpolator<3, 1>, Interpolator<3, 2>, Interpolator<3, 3>>;
+INSTANTIATE_TYPED_TEST_SUITE_P(testInterpolator, ACollectionOfParticles_3d, My3dTypes);
+/*********************************************************************************************/
+
 int main(int argc, char** argv)
 {
     ::testing::InitGoogleTest(&argc, argv);
diff --git a/tests/diagnostic/CMakeLists.txt b/tests/diagnostic/CMakeLists.txt
index 1554c1020..8ce79e4be 100644
--- a/tests/diagnostic/CMakeLists.txt
+++ b/tests/diagnostic/CMakeLists.txt
@@ -34,9 +34,11 @@ if(HighFive)
 
   _add_diagnostics_test(test-diagnostics_1d)
   _add_diagnostics_test(test-diagnostics_2d)
+  _add_diagnostics_test(test-diagnostics_3d)
 
   configure_file(${CMAKE_CURRENT_SOURCE_DIR}/job_1d.py.in ${CMAKE_CURRENT_BINARY_DIR}/job_1d.py @ONLY)
   configure_file(${CMAKE_CURRENT_SOURCE_DIR}/job_2d.py.in ${CMAKE_CURRENT_BINARY_DIR}/job_2d.py @ONLY)
+  configure_file(${CMAKE_CURRENT_SOURCE_DIR}/job_3d.py.in ${CMAKE_CURRENT_BINARY_DIR}/job_3d.py @ONLY)
   configure_file(${CMAKE_CURRENT_SOURCE_DIR}/__init__.py ${CMAKE_CURRENT_BINARY_DIR}/__init__.py @ONLY)
 
   message(STATUS "diagnostic working directory " ${PHARE_PROJECT_DIR})
diff --git a/tests/diagnostic/job_3d.py.in b/tests/diagnostic/job_3d.py.in
new file mode 100644
index 000000000..27e45c041
--- /dev/null
+++ b/tests/diagnostic/job_3d.py.in
@@ -0,0 +1,14 @@
+#!/usr/bin/env python3
+
+import pyphare.pharein as ph
+from pyphare.pharein import ElectronModel
+from tests.simulator import basicSimulatorArgs, makeBasicModel
+from tests.diagnostic import dump_all_diags
+
+out = "phare_outputs/diags_3d/"
+simInput = {"diag_options": {"format": "phareh5", "options": {"dir": out, "mode" : "overwrite"}}}
+
+ph.Simulation(**basicSimulatorArgs(dim = 3, interp = 1, **simInput))
+model = makeBasicModel()
+ElectronModel(closure="isothermal",Te = 0.12)
+dump_all_diags(model.populations)
diff --git a/tests/diagnostic/test-diagnostics_3d.cpp b/tests/diagnostic/test-diagnostics_3d.cpp
new file mode 100644
index 000000000..54e5741eb
--- /dev/null
+++ b/tests/diagnostic/test-diagnostics_3d.cpp
@@ -0,0 +1,35 @@
+
+#include "core/def/phare_mpi.hpp"
+
+#include "test_diagnostics.ipp"
+
+static std::string const job_file = "job_3d";
+static std::string const out_dir  = "phare_outputs/diags_3d/";
+
+TYPED_TEST(Simulator3dTest, fluid)
+{
+    fluid_test(TypeParam{job_file}, out_dir);
+}
+
+TYPED_TEST(Simulator3dTest, particles)
+{
+    particles_test(TypeParam{job_file}, out_dir);
+}
+
+TYPED_TEST(Simulator3dTest, electromag)
+{
+    electromag_test(TypeParam{job_file}, out_dir);
+}
+
+TYPED_TEST(Simulator3dTest, allFromPython)
+{
+    allFromPython_test(TypeParam{job_file}, out_dir);
+}
+
+
+int main(int argc, char** argv)
+{
+    ::testing::InitGoogleTest(&argc, argv);
+    PHARE::SamraiLifeCycle samsam(argc, argv);
+    return RUN_ALL_TESTS();
+}
diff --git a/tests/functional/harris/harris_3d.py b/tests/functional/harris/harris_3d.py
new file mode 100644
index 000000000..ed6cecf9b
--- /dev/null
+++ b/tests/functional/harris/harris_3d.py
@@ -0,0 +1,188 @@
+#!/usr/bin/env python3
+import os
+import numpy as np
+from pathlib import Path
+
+import pyphare.pharein as ph
+from pyphare.cpp import cpp_lib
+from pyphare.pharesee.run import Run
+from pyphare.simulator.simulator import Simulator
+from pyphare.simulator.simulator import startMPI
+
+os.environ["PHARE_SCOPE_TIMING"] = "0"  # turn on scope timing
+
+
+ph.NO_GUI()
+cpp = cpp_lib()
+startMPI()
+
+cells = (40, 40, 40)
+dl = (0.2, 0.2, 0.2)
+
+diag_outputs = "phare_outputs/test/harris/3d"
+time_step_nbr = 1000
+time_step = 0.001
+final_time = time_step * time_step_nbr
+
+timestamps = [0, final_time]
+if time_step_nbr > 100:
+    dt = 100 * time_step
+    nt = final_time / dt + 1
+    timestamps = dt * np.arange(nt)
+
+plot_dir = Path(f"{diag_outputs}_plots")
+plot_dir.mkdir(parents=True, exist_ok=True)
+
+
+def config():
+    sim = ph.Simulation(
+        time_step=time_step,
+        time_step_nbr=time_step_nbr,
+        dl=dl,
+        cells=cells,
+        refinement="tagging",
+        max_nbr_levels=1,
+        hyper_resistivity=0.001,
+        resistivity=0.001,
+        diag_options={
+            "format": "phareh5",
+            "options": {"dir": diag_outputs, "mode": "overwrite"},
+        },
+        # strict=True,
+    )
+
+    def density(x, y, z):
+        L = sim.simulation_domain()[1]
+        return (
+            0.2
+            + 1.0 / np.cosh((y - L * 0.3) / 0.5) ** 2
+            + 1.0 / np.cosh((y - L * 0.7) / 0.5) ** 2
+        )
+
+    def S(y, y0, l):
+        return 0.5 * (1.0 + np.tanh((y - y0) / l))
+
+    def by(x, y, z):
+        Lx = sim.simulation_domain()[0]
+        Ly = sim.simulation_domain()[1]
+        w1 = 0.2
+        w2 = 1.0
+        x0 = x - 0.5 * Lx
+        y1 = y - 0.3 * Ly
+        y2 = y - 0.7 * Ly
+        w3 = np.exp(-(x0 * x0 + y1 * y1) / (w2 * w2))
+        w4 = np.exp(-(x0 * x0 + y2 * y2) / (w2 * w2))
+        w5 = 2.0 * w1 / w2
+        return (w5 * x0 * w3) + (-w5 * x0 * w4)
+
+    def bx(x, y, z):
+        Lx = sim.simulation_domain()[0]
+        Ly = sim.simulation_domain()[1]
+        w1 = 0.2
+        w2 = 1.0
+        x0 = x - 0.5 * Lx
+        y1 = y - 0.3 * Ly
+        y2 = y - 0.7 * Ly
+        w3 = np.exp(-(x0 * x0 + y1 * y1) / (w2 * w2))
+        w4 = np.exp(-(x0 * x0 + y2 * y2) / (w2 * w2))
+        w5 = 2.0 * w1 / w2
+        v1 = -1
+        v2 = 1.0
+        return (
+            v1
+            + (v2 - v1) * (S(y, Ly * 0.3, 0.5) - S(y, Ly * 0.7, 0.5))
+            + (-w5 * y1 * w3)
+            + (+w5 * y2 * w4)
+        )
+
+    def bz(x, y, z):
+        return bx(x, y, z) * -1
+
+    def b2(x, y, z):
+        return bx(x, y, z) ** 2 + by(x, y, z) ** 2 + bz(x, y, z) ** 2
+
+    def T(x, y, z):
+        K = 1
+        temp = 1.0 / density(x, y, z) * (K - b2(x, y, z) * 0.5)
+        assert np.all(temp > 0)
+        return temp
+
+    def vxyz(x, y, z):
+        return 0.0
+
+    def vthxyz(x, y, z):
+        return np.sqrt(T(x, y, z))
+
+    C = "xyz"
+    vvv = {
+        **{f"vbulk{c}": vxyz for c in C},
+        **{f"vth{c}": vthxyz for c in C},
+        "nbr_part_per_cell": 100,
+    }
+    protons = {"charge": 1, "density": density, **vvv, "init": {"seed": 12334}}
+    ph.MaxwellianFluidModel(bx=bx, by=by, bz=bz, protons=protons)
+    ph.ElectronModel(closure="isothermal", Te=0.0)
+
+    for quantity in ["density", "bulkVelocity"]:
+        ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
+    ph.FluidDiagnostics(
+        quantity="density", write_timestamps=timestamps, population_name="protons"
+    )
+    for quantity in ["E", "B"]:
+        ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
+    ph.InfoDiagnostics(quantity="particle_count")  # defaults all coarse time steps
+
+    return sim
+
+
+def plot_file_for_qty(qty, time):
+    return f"{plot_dir}/harris_{qty}_t{time}.png"
+
+
+def plot(diag_dir):
+    run = Run(diag_dir)
+    for time in timestamps:
+        # run.GetDivB(time).plot(
+        #     filename=plot_file_for_qty("divb", time),
+        #     plot_patches=True,
+        #     vmin=1e-11,
+        #     vmax=2e-10,
+        # )
+        # run.GetRanks(time).plot(
+        #     filename=plot_file_for_qty("Ranks", time),
+        #     plot_patches=True,
+        # )
+        run.GetN(time, pop_name="protons").plot(
+            filename=plot_file_for_qty("N", time),
+            plot_patches=True,
+        )
+        for c in ["x", "y", "z"]:
+            run.GetB(time).plot(
+                filename=plot_file_for_qty(f"b{c}", time),
+                qty=f"B{c}",
+                plot_patches=True,
+            )
+        # run.GetJ(time).plot(
+        #     filename=plot_file_for_qty("jz", time),
+        #     qty="Jz",
+        #     plot_patches=True,
+        #     vmin=-2,
+        #     vmax=2,
+        # )
+
+
+def main():
+    Simulator(config()).run()
+    if cpp.mpi_rank() == 0:
+        plot(diag_outputs)
+    # try:
+    #     from tools.python3 import plotting as m_plotting
+
+    #     m_plotting.plot_run_timer_data(diag_outputs, cpp.mpi_rank())
+    # except ImportError:
+    #     print("Phlop not found - install with: `pip install phlop`")
+    cpp.mpi_barrier()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/tests/simulator/__init__.py b/tests/simulator/__init__.py
index 5c100634c..5125ba623 100644
--- a/tests/simulator/__init__.py
+++ b/tests/simulator/__init__.py
@@ -94,12 +94,17 @@ def density_2d_periodic(sim, x, y):
     )
 
 
-# def density_3d_periodic(sim, x, y, z):
-#     xmax, ymax, zmax = sim.simulation_domain()
-#     background_particles = 0.3  # avoids 0 density
-#     xx, yy, zz = meshify(x, y, z)
-#     r = np.exp(-(xx-0.5*xmax)**2)*np.exp(-(yy-ymax/2.)**2)*np.exp(-(zz-zmax/2.)**2) + background_particles
-#     return r
+def density_3d_periodic(sim, x, y, z):
+    xmax, ymax, zmax = sim.simulation_domain()
+    background_particles = 0.3  # avoids 0 density
+    xx, yy, zz = meshify(x, y, z)
+    r = (
+        np.exp(-((xx - 0.5 * xmax) ** 2))
+        * np.exp(-((yy - ymax / 2.0) ** 2))
+        * np.exp(-((zz - zmax / 2.0) ** 2))
+        + background_particles
+    )
+    return r
 
 
 def defaultPopulationSettings(sim, density_fn, vbulk_fn):
@@ -308,3 +313,33 @@ def clean_up_diags_dirs(self):
                 if os.path.exists(diag_dir):
                     shutil.rmtree(diag_dir)
         cpp_lib().mpi_barrier()
+
+
+def debug_tracer():
+    """
+    print live stack trace during execution
+    """
+
+    import sys, os
+
+    def tracefunc(frame, event, arg, indent=[0]):
+        filename = os.path.basename(frame.f_code.co_filename)
+        line_number = frame.f_lineno
+        if event == "call":
+            indent[0] += 2
+            print(
+                "-" * indent[0] + "> enter function",
+                frame.f_code.co_name,
+                f"{filename} {line_number}",
+            )
+        elif event == "return":
+            print(
+                "<" + "-" * indent[0],
+                "exit function",
+                frame.f_code.co_name,
+                f"{filename} {line_number}",
+            )
+            indent[0] -= 2
+        return tracefunc
+
+    sys.setprofile(tracefunc)
diff --git a/tests/simulator/advance/CMakeLists.txt b/tests/simulator/advance/CMakeLists.txt
index 0b4ad77a5..4304d8aa0 100644
--- a/tests/simulator/advance/CMakeLists.txt
+++ b/tests/simulator/advance/CMakeLists.txt
@@ -18,6 +18,11 @@ if(HighFive)
     phare_mpi_python3_exec(9 ${PHARE_MPI_PROCS} advance-2d-fields     test_fields_advance_2d.py     ${CMAKE_CURRENT_BINARY_DIR})
     phare_mpi_python3_exec(9 ${PHARE_MPI_PROCS} advance-2d-particles  test_particles_advance_2d.py  ${CMAKE_CURRENT_BINARY_DIR})
 
+    if(highResourceTests) # off by default as it's quite intensive
+        phare_mpi_python3_exec(9 ${PHARE_MPI_PROCS} advance-3d-fields     test_fields_advance_3d.py     ${CMAKE_CURRENT_BINARY_DIR})
+        phare_mpi_python3_exec(9 ${PHARE_MPI_PROCS} advance-3d-particles  test_particles_advance_3d.py  ${CMAKE_CURRENT_BINARY_DIR})
+    endif(highResourceTests)
+
   endif()
 
 endif()
diff --git a/tests/simulator/advance/test_fields_advance_3d.py b/tests/simulator/advance/test_fields_advance_3d.py
new file mode 100644
index 000000000..272ff40a8
--- /dev/null
+++ b/tests/simulator/advance/test_fields_advance_3d.py
@@ -0,0 +1,107 @@
+"""
+  This file exists independently from test_advance.py to isolate dimension
+    test cases and allow each to be overridden in some way if required.
+"""
+
+import unittest
+
+import matplotlib
+from ddt import data, ddt, unpack
+from pyphare.core.box import Box3D
+
+from tests.simulator.test_advance import AdvanceTestBase
+
+matplotlib.use("Agg")  # for systems without GUI
+
+ndim = 3
+interp_orders = [1, 2, 3]
+ppc, cells = 10, 30
+
+
+def per_interp(dic):
+    return [(interp, dic) for interp in interp_orders]
+
+
+@ddt
+class AdvanceTest(AdvanceTestBase):
+    @data(
+        *per_interp({}),
+        *per_interp({"L0": [Box3D(10, 19)]}),
+        *per_interp({"L0": [Box3D(8, 20)]}),
+    )
+    @unpack
+    def test_overlaped_fields_are_equal(self, interp_order, refinement_boxes):
+        print(f"{self._testMethodName}_{ndim}d")
+        time_step_nbr = 3
+        time_step = 0.001
+        datahier = self.getHierarchy(
+            ndim,
+            interp_order,
+            refinement_boxes,
+            "eb",
+            cells=cells,
+            time_step=time_step,
+            time_step_nbr=time_step_nbr,
+            nbr_part_per_cell=ppc,
+        )
+        self._test_overlaped_fields_are_equal(datahier, time_step_nbr, time_step)
+
+    @data(
+        *per_interp({}),
+        *per_interp({"L0": [Box3D(5, 14)]}),
+    )
+    @unpack
+    def test_overlaped_fields_are_equal_with_min_max_patch_size_of_max_ghosts(
+        self, interp_order, refinement_boxes
+    ):
+        print(f"{self._testMethodName}_{ndim}d")
+        time_step_nbr = 3
+        time_step = 0.001
+        from pyphare.pharein.simulation import check_patch_size
+
+        largest_patch_size, smallest_patch_size = check_patch_size(
+            ndim, interp_order=interp_order, cells=[cells] * ndim
+        )
+        datahier = self.getHierarchy(
+            ndim,
+            interp_order,
+            refinement_boxes,
+            "eb",
+            cells=cells,
+            smallest_patch_size=smallest_patch_size,
+            largest_patch_size=smallest_patch_size,
+            time_step=time_step,
+            time_step_nbr=time_step_nbr,
+            nbr_part_per_cell=ppc,
+        )
+        self._test_overlaped_fields_are_equal(datahier, time_step_nbr, time_step)
+
+    # @data(
+    #    *per_interp(({"L0": {"B0": Box3D(10, 14)}})),
+    #    *per_interp(({"L0": {"B0": Box3D(10, 14), "B1": Box3D(15, 19)}})),
+    #    *per_interp(({"L0": {"B0": Box3D(6, 23)}})),
+    #    *per_interp(({"L0": {"B0": Box3D( 2, 12), "B1": Box3D(13, 25)}})),
+    #    *per_interp(({"L0": {"B0": Box3D( 5, 20)}, "L1": {"B0": Box3D(15, 19)}})),
+    #    *per_interp(({"L0": {"B0": Box3D( 5, 20)}, "L1": {"B0": Box3D(12, 38)}, "L2": {"B0": Box3D(30, 52)} })),
+    # )
+    # @unpack
+    # def test_field_coarsening_via_subcycles(self, interp_order, refinement_boxes):
+    #     print(f"{self._testMethodName}_{ndim}d")
+    #     self._test_field_coarsening_via_subcycles(ndim, interp_order, refinement_boxes, dl=.3, cells=cells)
+
+    # @unittest.skip("should change to work on moments")
+    # @data( # only supports a hierarchy with 2 levels
+    #    *per_interp(({"L0": [Box3D(0, 4)]})),
+    #    *per_interp(({"L0": [Box3D(10, 14)]})),
+    #    *per_interp(({"L0": [Box3D(0, 4), Box3D(10, 14)]})),
+    #    *per_interp(({"L0": [Box3D(0, 4), Box3D(5, 9), Box3D(10, 14)]})),
+    #    *per_interp(({"L0": [Box3D(20, 24)]})),
+    # )
+    # @unpack
+    # def test_field_level_ghosts_via_subcycles_and_coarser_interpolation(self, interp_order, refinement_boxes):
+    #     print(f"{self._testMethodName}_{ndim}d")
+    #     self._test_field_level_ghosts_via_subcycles_and_coarser_interpolation(ndim, interp_order, refinement_boxes)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/simulator/advance/test_particles_advance_3d.py b/tests/simulator/advance/test_particles_advance_3d.py
new file mode 100644
index 000000000..a095986d9
--- /dev/null
+++ b/tests/simulator/advance/test_particles_advance_3d.py
@@ -0,0 +1,51 @@
+"""
+  This file exists independently from test_advance.py to isolate dimension
+    test cases and allow each to be overridden in some way if required.
+"""
+
+import unittest
+
+import matplotlib
+from ddt import data, ddt, unpack
+from pyphare.core.box import Box3D
+
+from tests.simulator.test_advance import AdvanceTestBase
+
+matplotlib.use("Agg")  # for systems without GUI
+
+ndim = 3
+interp_orders = [1, 2, 3]
+ppc = 5
+
+
+def per_interp(dic):
+    return [(interp, dic) for interp in interp_orders]
+
+
+@ddt
+class AdvanceTest(AdvanceTestBase):
+    @data(
+        *per_interp({}),
+        *per_interp({"L0": [Box3D(10, 19)]}),
+        *per_interp({"L0": [Box3D(5, 9), Box3D(10, 14)]}),
+    )
+    @unpack
+    def test_overlapped_particledatas_have_identical_particles(
+        self, interp_order, refinement_boxes
+    ):
+        self._test_overlapped_particledatas_have_identical_particles(
+            ndim,
+            interp_order,
+            refinement_boxes,
+            ppc=ppc,
+            cells=40,
+            largest_patch_size=20,
+        )
+
+    @data(*interp_orders)
+    def test_L0_particle_number_conservation(self, interp):
+        self._test_L0_particle_number_conservation(ndim, interp, ppc=ppc, cells=30)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/simulator/initialize/CMakeLists.txt b/tests/simulator/initialize/CMakeLists.txt
index cbaf44d02..1f1c93827 100644
--- a/tests/simulator/initialize/CMakeLists.txt
+++ b/tests/simulator/initialize/CMakeLists.txt
@@ -18,6 +18,11 @@ if(HighFive)
     phare_mpi_python3_exec(9 ${PHARE_MPI_PROCS} init-2d-fields     test_fields_init_2d.py     ${CMAKE_CURRENT_BINARY_DIR})
     phare_mpi_python3_exec(9 ${PHARE_MPI_PROCS} init-2d-particles  test_particles_init_2d.py  ${CMAKE_CURRENT_BINARY_DIR})
 
+    if(highResourceTests)
+        phare_mpi_python3_exec(9 ${PHARE_MPI_PROCS} init-3d-fields     test_fields_init_3d.py     ${CMAKE_CURRENT_BINARY_DIR})
+        phare_mpi_python3_exec(9 ${PHARE_MPI_PROCS} init-3d-particles  test_particles_init_3d.py  ${CMAKE_CURRENT_BINARY_DIR})
+    endif()
+
   endif()
 
 endif()
diff --git a/tests/simulator/initialize/test_fields_init_2d.py b/tests/simulator/initialize/test_fields_init_2d.py
index af47a5d99..22e4587d1 100644
--- a/tests/simulator/initialize/test_fields_init_2d.py
+++ b/tests/simulator/initialize/test_fields_init_2d.py
@@ -23,7 +23,7 @@ class Initialization2DTest(InitializationTest):
     @data(*interp_orders)
     def test_B_is_as_provided_by_user(self, interp_order):
         print(f"\n{self._testMethodName}_{ndim}d")
-        self._test_B_is_as_provided_by_user(ndim, interp_order, nbr_part_per_cell=ppc)
+        self._test_B_is_as_provided_by_user(ndim, interp_order, ppc=ppc)
 
     @data(*interp_orders)
     def test_bulkvel_is_as_provided_by_user(self, interp_order):
diff --git a/tests/simulator/initialize/test_fields_init_3d.py b/tests/simulator/initialize/test_fields_init_3d.py
new file mode 100644
index 000000000..05d8f5a53
--- /dev/null
+++ b/tests/simulator/initialize/test_fields_init_3d.py
@@ -0,0 +1,49 @@
+"""
+  This file exists independently from test_initialization.py to isolate dimension
+    test cases and allow each to be overridden in some way if required.
+"""
+
+import unittest
+
+import numpy as np
+import matplotlib
+from ddt import data, ddt
+
+from tests.simulator.test_initialization import InitializationTest
+
+matplotlib.use("Agg")  # for systems without GUI
+
+ndim = 3
+interp_orders = [1, 2, 3]
+ppc, cells = 10, 20
+
+
+@ddt
+class Initialization3DTest(InitializationTest):
+    @data(*interp_orders)
+    def test_B_is_as_provided_by_user(self, interp_order):
+        print(f"\n{self._testMethodName}_{ndim}d")
+        self._test_B_is_as_provided_by_user(ndim, interp_order, ppc=ppc, cells=cells)
+
+    @data(*interp_orders)
+    def test_bulkvel_is_as_provided_by_user(self, interp_order):
+        print(f"\n{self._testMethodName}_{ndim}d")
+        self._test_bulkvel_is_as_provided_by_user(
+            ndim, interp_order, ppc=ppc, cells=cells
+        )
+
+    @data(*interp_orders)
+    def test_density_is_as_provided_by_user(self, interp_order):
+        print(f"\n{self._testMethodName}_{ndim}d")
+        self._test_density_is_as_provided_by_user(ndim, interp_order, cells=cells)
+
+    @data(*interp_orders)  # uses too much RAM - to isolate somehow
+    def test_density_decreases_as_1overSqrtN(self, interp_order):
+        print(f"\n{self._testMethodName}_{ndim}d")
+        self._test_density_decreases_as_1overSqrtN(
+            ndim, interp_order, np.asarray([10, 25, 50, 75]), cells=cells
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/simulator/initialize/test_particles_init_1d.py b/tests/simulator/initialize/test_particles_init_1d.py
index 60c44d692..708995b33 100644
--- a/tests/simulator/initialize/test_particles_init_1d.py
+++ b/tests/simulator/initialize/test_particles_init_1d.py
@@ -8,12 +8,9 @@
 import matplotlib
 from ddt import data, ddt, unpack
 from pyphare.core.box import Box1D
-from pyphare.cpp import cpp_lib
 
 from tests.simulator.test_initialization import InitializationTest
 
-cpp = cpp_lib()
-
 matplotlib.use("Agg")  # for systems without GUI
 
 ndim = 1
diff --git a/tests/simulator/initialize/test_particles_init_2d.py b/tests/simulator/initialize/test_particles_init_2d.py
index cc56392f8..758d3b0a2 100644
--- a/tests/simulator/initialize/test_particles_init_2d.py
+++ b/tests/simulator/initialize/test_particles_init_2d.py
@@ -23,7 +23,7 @@ def per_interp(dic):
 
 
 @ddt
-class Initialization1DTest(InitializationTest):
+class Initialization2DTest(InitializationTest):
     @data(*interp_orders)
     def test_nbr_particles_per_cell_is_as_provided(self, interp_order):
         print(f"{self._testMethodName}_{ndim}d")
diff --git a/tests/simulator/initialize/test_particles_init_3d.py b/tests/simulator/initialize/test_particles_init_3d.py
new file mode 100644
index 000000000..a87c5328b
--- /dev/null
+++ b/tests/simulator/initialize/test_particles_init_3d.py
@@ -0,0 +1,98 @@
+"""
+  This file exists independently from test_initialization.py to isolate dimension
+    test cases and allow each to be overridden in some way if required.
+"""
+
+import unittest
+
+import matplotlib
+from ddt import data, ddt, unpack
+from pyphare.core.box import Box3D
+
+from tests.simulator.test_initialization import InitializationTest
+
+matplotlib.use("Agg")  # for systems without GUI
+
+ndim = 3
+interp_orders = [1, 2, 3]
+ppc, cells = 10, 30
+
+
+def per_interp(dic):
+    return [(interp, dic) for interp in interp_orders]
+
+
+@ddt
+class Initialization3DTest(InitializationTest):
+    @data(*interp_orders)
+    def test_nbr_particles_per_cell_is_as_provided(self, interp_order):
+        print(f"{self._testMethodName}_{ndim}d")
+        self._test_nbr_particles_per_cell_is_as_provided(
+            ndim, interp_order, ppc, cells=cells
+        )
+
+    @data(
+        *per_interp(({"L0": {"B0": Box3D(10, 14)}})),
+        *per_interp(({"L0": {"B0": Box3D(10, 14)}, "L1": {"B0": Box3D(22, 26)}})),
+        *per_interp(({"L0": {"B0": Box3D(2, 6), "B1": Box3D(7, 11)}})),
+    )
+    @unpack
+    def test_levelghostparticles_have_correct_split_from_coarser_particle(
+        self, interp_order, refinement_boxes
+    ):
+        print(f"\n{self._testMethodName}_{ndim}d")
+        now = self.datetime_now()
+        self._test_levelghostparticles_have_correct_split_from_coarser_particle(
+            self.getHierarchy(
+                ndim,
+                interp_order,
+                refinement_boxes,
+                "particles",
+                cells=cells,
+                nbr_part_per_cell=ppc,
+            )
+        )
+        print(
+            f"\n{self._testMethodName}_{ndim}d took {self.datetime_diff(now)} seconds"
+        )
+
+    @data(
+        *per_interp(({"L0": {"B0": Box3D(10, 14)}})),
+        *per_interp(({"L0": {"B0": Box3D(5, 14)}, "L1": {"B0": Box3D(15, 19)}})),
+        *per_interp(({"L0": {"B0": Box3D(2, 12), "B1": Box3D(13, 25)}})),
+    )
+    @unpack
+    def test_domainparticles_have_correct_split_from_coarser_particle(
+        self, interp_order, refinement_boxes
+    ):
+        print(f"\n{self._testMethodName}_{ndim}d")
+        now = self.datetime_now()
+        self._test_domainparticles_have_correct_split_from_coarser_particle(
+            ndim, interp_order, refinement_boxes, nbr_part_per_cell=ppc
+        )
+        print(
+            f"\n{self._testMethodName}_{ndim}d took {self.datetime_diff(now)} seconds"
+        )
+
+    # @data({"cells": 40, "smallest_patch_size": 20, "largest_patch_size": 20, "nbr_part_per_cell" : ppc})
+    # def test_no_patch_ghost_on_refined_level_case(self, simInput):
+    #     print(f"\n{self._testMethodName}_{ndim}d")
+    #     now = self.datetime_now()
+    #     self._test_patch_ghost_on_refined_level_case(ndim, False, **simInput)
+    #     print(f"\n{self._testMethodName}_{ndim}d took {self.datetime_diff(now)} seconds")
+
+    # @data({"cells": 40, "interp_order": 1, "nbr_part_per_cell" : ppc})
+    # def test_has_patch_ghost_on_refined_level_case(self, simInput):
+    #     print(f"\n{self._testMethodName}_{ndim}d")
+    #     from pyphare.pharein.simulation import check_patch_size
+    #     diag_outputs=f"phare_overlaped_fields_are_equal_with_min_max_patch_size_of_max_ghosts_{ndim}_{self.ddt_test_id()}"
+    #     _, smallest_patch_size = check_patch_size(ndim, **simInput)
+    #     simInput["smallest_patch_size"] = smallest_patch_size
+    #     simInput["largest_patch_size"] = smallest_patch_size
+    #     now = self.datetime_now()
+    #     self._test_patch_ghost_on_refined_level_case(ndim, True, **simInput)
+    #     print(f"\n{self._testMethodName}_{ndim}d took {self.datetime_diff(now)} seconds")
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/simulator/per_test.hpp b/tests/simulator/per_test.hpp
index d4a8be4e0..901c7c428 100644
--- a/tests/simulator/per_test.hpp
+++ b/tests/simulator/per_test.hpp
@@ -108,6 +108,14 @@ using Simulators2d = testing::Types<
 TYPED_TEST_SUITE(Simulator2dTest, Simulators2d);
 
 
+template<typename Simulator>
+struct Simulator3dTest : public ::testing::Test
+{
+};
+using Simulator3d = testing::Types<SimulatorTestParam<3, 1, 6>, SimulatorTestParam<3, 2, 6>,
+                                   SimulatorTestParam<3, 3, 6>>;
+TYPED_TEST_SUITE(Simulator3dTest, Simulator3d);
+
 
 
 #endif /* PHARE_TEST_SIMULATOR_PER_TEST_H */
diff --git a/tests/simulator/refined_particle_nbr.py b/tests/simulator/refined_particle_nbr.py
index 95e21dda9..c111738f0 100644
--- a/tests/simulator/refined_particle_nbr.py
+++ b/tests/simulator/refined_particle_nbr.py
@@ -30,27 +30,21 @@ def __init__(self, *args, **kwargs):
                 print(exc)
                 sys.exit(1)
 
+    # needed in case exception is raised in test and Simulator not reset properly
+    def tearDown(self):
+        if self.simulator is not None:
+            self.simulator.reset()
+
     # while splitting particles may leave the domain area
     #  so we remove the particles from the border cells of each patch
     def _less_per_dim(self, dim, refined_particle_nbr, patch):
         if dim == 1:
             return refined_particle_nbr * 2
         cellNbr = patch.upper - patch.lower + 1
+        dim2 = refined_particle_nbr * ((cellNbr[0] * 2 + (cellNbr[1] * 2)))
         if dim == 2:
-            return refined_particle_nbr * ((cellNbr[0] * 2 + (cellNbr[1] * 2)))
-        raise ValueError("Unhandled dimension for function")
-
-    def _check_deltas_and_weights(self, dim, interp, refined_particle_nbr):
-        yaml_dim = self.yaml_root["dimension_" + str(dim)]
-        yaml_interp = yaml_dim["interp_" + str(interp)]
-        yaml_n_particles = yaml_interp["N_particles_" + str(refined_particle_nbr)]
-
-        yaml_delta = [float(s) for s in str(yaml_n_particles["delta"]).split(" ")]
-        yaml_weight = [float(s) for s in str(yaml_n_particles["weight"]).split(" ")]
-
-        splitter_t = splitter_type(dim, interp, refined_particle_nbr)
-        np.testing.assert_allclose(yaml_delta, splitter_t.delta)
-        np.testing.assert_allclose(yaml_weight, splitter_t.weight)
+            return dim2
+        return dim2 * (cellNbr[2] * 2)
 
     def _do_dim(self, dim, min_diff, max_diff):
         from pyphare.pharein.simulation import valid_refined_particle_nbr
@@ -58,11 +52,7 @@ def _do_dim(self, dim, min_diff, max_diff):
         for interp in range(1, 4):
             prev_split_particle_max = 0
             for refined_particle_nbr in valid_refined_particle_nbr[dim][interp]:
-                self._check_deltas_and_weights(dim, interp, refined_particle_nbr)
-
-                simInput = NoOverwriteDict(
-                    {"refined_particle_nbr": refined_particle_nbr}
-                )
+                simInput = {"refined_particle_nbr": refined_particle_nbr}
                 self.simulator = Simulator(populate_simulation(dim, interp, **simInput))
                 self.simulator.initialize()
                 dw = self.simulator.data_wrangler()
@@ -78,9 +68,6 @@ def _do_dim(self, dim, min_diff, max_diff):
                             per_pop += patch.data.size()
                     max_per_pop = max(max_per_pop, per_pop)
                 prev_min_diff = prev_split_particle_max * min_diff
-
-                # while splitting particles may leave the domain area
-                #  so we remove the particles from the border cells of each patch
                 self.assertTrue(max_per_pop > prev_min_diff - (leaving_particles))
                 if prev_split_particle_max > 0:
                     prev_max_diff = prev_min_diff * dim * max_diff
@@ -88,45 +75,70 @@ def _do_dim(self, dim, min_diff, max_diff):
                 prev_split_particle_max = max_per_pop
                 self.simulator = None
 
-    """ 1d
-      refine 10 cells in 1d, ppc 100
-        10 * 2 * ppc = 200
-        10 * 3 * ppc = 300 300 / 200 = 1.5
-        10 * 4 * ppc = 400 500 / 400 = 1.33
-        10 * 5 * ppc = 500 500 / 400 = 1.25
-      taking the minimul diff across permutations
-      current to previous should be at least this times bigger
-    """
-    PREVIOUS_ITERATION_MIN_DIFF_1d = 1.25
-    PREVIOUS_ITERATION_MAX_DIFF_1d = 1.50
+    #  1d
+    #  refine 10 cells in 1d, ppc 100
+    #    10 * 2 * ppc = 200
+    #    10 * 3 * ppc = 300 300 / 200 = 1.5
+    #    10 * 4 * ppc = 400 500 / 400 = 1.33
+    #    10 * 5 * ppc = 500 500 / 400 = 1.25
+    #  taking the minimul diff across permutations
+    #  current to previous should be at least this times bigger
+    PRIOR_MIN_DIFF_1d = 1.25
+    PRIOR_MAX_DIFF_1d = 1.50
 
     def test_1d(self):
         This = type(self)
-        self._do_dim(
-            1, This.PREVIOUS_ITERATION_MIN_DIFF_1d, This.PREVIOUS_ITERATION_MAX_DIFF_1d
-        )
-
-    """ 2d
-      refine 10x10 cells in 2d, ppc 100
-      10 * 10 * 4 * ppc = 400
-      10 * 10 * 8 * ppc = 800 800 / 400 = 1.5
-      10 * 10 * 9 * ppc = 900 900 / 800 = 1.125
-    """
-    PREVIOUS_ITERATION_MIN_DIFF_2d = 1.125
-    PREVIOUS_ITERATION_MAX_DIFF_2d = 1.50
+        self._do_dim(1, This.PRIOR_MIN_DIFF_1d, This.PRIOR_MAX_DIFF_1d)
+
+    #  2d
+    #  refine 10x10 cells in 2d, ppc 100
+    #    10 * 10 * 4 * ppc = 400
+    #    10 * 10 * 8 * ppc = 800 800 / 400 = 1.5
+    #    10 * 10 * 9 * ppc = 900 900 / 800 = 1.125
+    PRIOR_MIN_DIFF_2d = 1.125
+    PRIOR_MAX_DIFF_2d = 1.50
 
     def test_2d(self):
         This = type(self)
-        self._do_dim(
-            2, This.PREVIOUS_ITERATION_MIN_DIFF_2d, This.PREVIOUS_ITERATION_MAX_DIFF_2d
-        )
+        self._do_dim(2, This.PRIOR_MIN_DIFF_2d, This.PRIOR_MAX_DIFF_2d)
 
-    def tearDown(self):
-        # needed in case exception is raised in test and Simulator
-        # not reset properly
-        if self.simulator is not None:
-            self.simulator.reset()
+    # 3d
+    #  refine 10x10x10 cells in 3d, ppc 100
+    #    10 * 10 * 10 * 6 *  ppc = 6000
+    #    10 * 10 * 10 * 12 * ppc = 12000 - 12000 / 6000  = 2
+    #    10 * 10 * 10 * 27 * ppc = 27000 - 27000 / 12000 = 2.25
+    PRIOR_MIN_DIFF_3d = 2
+    PRIOR_MAX_DIFF_3d = 2.25
+
+    def test_3d(self):
+        This = type(self)
+        self._do_dim(3, This.PRIOR_MIN_DIFF_3d, This.PRIOR_MAX_DIFF_3d)
+
+    def _check_deltas_and_weights(self, dim, interp, refined_particle_nbr):
+        yaml_dim = self.yaml_root["dimension_" + str(dim)]
+        yaml_interp = yaml_dim["interp_" + str(interp)]
+        yaml_n_particles = yaml_interp["N_particles_" + str(refined_particle_nbr)]
+
+        yaml_delta = [float(s) for s in str(yaml_n_particles["delta"]).split(" ")]
+        yaml_weight = [float(s) for s in str(yaml_n_particles["weight"]).split(" ")]
+
+        splitter_t = splitter_type(dim, interp, refined_particle_nbr)
+        np.testing.assert_allclose(yaml_delta, splitter_t.delta)
+        np.testing.assert_allclose(yaml_weight, splitter_t.weight)
+
+    def test_values(self):
+        from pyphare.pharein.simulation import valid_refined_particle_nbr
+
+        for dim in range(1, 4):
+            for interp in range(1, 4):
+                for refined_particle_nbr in valid_refined_particle_nbr[dim][interp]:
+                    self._check_deltas_and_weights(dim, interp, refined_particle_nbr)
 
 
 if __name__ == "__main__":
-    unittest.main()
+    # the following ensures the order of execution so delta/weights are verified first
+    suite = unittest.TestSuite()
+    suite.addTest(SimulatorRefinedParticleNbr("test_values"))
+    for dim in range(1, 4):
+        suite.addTest(SimulatorRefinedParticleNbr("test_" + str(dim) + "d"))
+    unittest.TextTestRunner(failfast=True).run(suite)
diff --git a/tests/simulator/test_advance.py b/tests/simulator/test_advance.py
index 569f2c873..ecb8db623 100644
--- a/tests/simulator/test_advance.py
+++ b/tests/simulator/test_advance.py
@@ -86,6 +86,9 @@ def getHierarchy(
             strict=True,
         )
 
+        def S(x, x0, l):
+            return 0.5 * (1 + np.tanh((x - x0) / l))
+
         def bx(*xyz):
             return 1.0
 
@@ -259,29 +262,15 @@ def base_test_overlaped_fields_are_equal(self, datahier, coarsest_time):
                     # this is because the overlap box has been calculated from
                     # the intersection of possibly shifted patch data ghost boxes
 
-                    loc_b1 = boxm.amr_to_local(
-                        box, boxm.shift(pd1.ghost_box, offsets[0])
+                    slice1 = boxm.select(
+                        pd1.dataset,
+                        boxm.amr_to_local(box, boxm.shift(pd1.ghost_box, offsets[0])),
                     )
-                    loc_b2 = boxm.amr_to_local(
-                        box, boxm.shift(pd2.ghost_box, offsets[1])
+                    slice2 = boxm.select(
+                        pd2.dataset,
+                        boxm.amr_to_local(box, boxm.shift(pd2.ghost_box, offsets[1])),
                     )
-
-                    data1 = pd1.dataset
-                    data2 = pd2.dataset
-
-                    if box.ndim == 1:
-                        slice1 = data1[loc_b1.lower[0] : loc_b1.upper[0] + 1]
-                        slice2 = data2[loc_b2.lower[0] : loc_b2.upper[0] + 1]
-
-                    if box.ndim == 2:
-                        slice1 = data1[
-                            loc_b1.lower[0] : loc_b1.upper[0] + 1,
-                            loc_b1.lower[1] : loc_b1.upper[1] + 1,
-                        ]
-                        slice2 = data2[
-                            loc_b2.lower[0] : loc_b2.upper[0] + 1,
-                            loc_b2.lower[1] : loc_b2.upper[1] + 1,
-                        ]
+                    assert slice1.dtype == np.float64
 
                     try:
                         # empirical max absolute observed 5.2e-15
@@ -298,12 +287,9 @@ def base_test_overlaped_fields_are_equal(self, datahier, coarsest_time):
                         print(pd1.y.mean())
                         print(pd2.x.mean())
                         print(pd2.y.mean())
-                        print(loc_b1)
-                        print(loc_b2)
                         print(coarsest_time)
                         print(slice1)
                         print(slice2)
-                        print(data1[:])
                         if self.rethrow_:
                             raise e
                         return diff_boxes(slice1, slice2, box)
@@ -381,8 +367,9 @@ def _test_overlapped_particledatas_have_identical_particles(
                         part2.iCells = part2.iCells + offsets[1]
                         self.assertEqual(part1, part2)
 
-    def _test_L0_particle_number_conservation(self, ndim, interp_order, ppc=100):
-        cells = 120
+    def _test_L0_particle_number_conservation(
+        self, ndim, interp_order, ppc=100, cells=120
+    ):
         time_step_nbr = 10
         time_step = 0.001
 
@@ -409,7 +396,7 @@ def _test_L0_particle_number_conservation(self, ndim, interp_order, ppc=100):
             self.assertEqual(n_particles, n_particles_at_t)
 
     def _test_field_coarsening_via_subcycles(
-        self, dim, interp_order, refinement_boxes, **kwargs
+        self, dim, interp_order, refinement_boxes, cells=60, **kwargs
     ):
         print(
             "test_field_coarsening_via_subcycles for dim/interp : {}/{}".format(
@@ -423,12 +410,14 @@ def _test_field_coarsening_via_subcycles(
 
         time_step_nbr = 3
 
+        diag_outputs = f"subcycle_coarsening/{dim}/{interp_order}/{self.ddt_test_id()}"
         datahier = self.getHierarchy(
             dim,
             interp_order,
             refinement_boxes,
             "fields",
-            cells=60,
+            cells=cells,
+            diag_outputs=diag_outputs,
             time_step=0.001,
             extra_diag_options={"fine_dump_lvl_max": 10},
             time_step_nbr=time_step_nbr,
@@ -501,16 +490,7 @@ def _test_field_coarsening_via_subcycles(
                                 afterCoarse = np.copy(coarse_pdDataset)
 
                                 # change values that should be updated to make failure obvious
-                                assert dim < 3  # update
-                                if dim == 1:
-                                    afterCoarse[
-                                        dataBox.lower[0] : dataBox.upper[0] + 1
-                                    ] = -144123
-                                if dim == 2:
-                                    afterCoarse[
-                                        dataBox.lower[0] : dataBox.upper[0] + 1,
-                                        dataBox.lower[1] : dataBox.upper[1] + 1,
-                                    ] = -144123
+                                boxm.DataSelector(afterCoarse)[dataBox] = -144123
 
                                 coarsen(
                                     qty,
diff --git a/tests/simulator/test_initialization.py b/tests/simulator/test_initialization.py
index 2ea135909..70937fbf6 100644
--- a/tests/simulator/test_initialization.py
+++ b/tests/simulator/test_initialization.py
@@ -54,7 +54,6 @@ def getHierarchy(
         diag_outputs = self.unique_diag_dir_for_test_case(
             "phare_outputs/init", ndim, interp_order, diag_outputs
         )
-
         from pyphare.pharein import global_vars
 
         global_vars.sim = None
@@ -252,19 +251,26 @@ def vthz(*xyz):
                 )
             return mom_hier
 
-    def _test_B_is_as_provided_by_user(self, dim, interp_order, **kwargs):
+    def _test_B_is_as_provided_by_user(self, dim, interp_order, ppc=100, **kwargs):
         print(
             "test_B_is_as_provided_by_user : dim  {} interp_order : {}".format(
                 dim, interp_order
             )
         )
+        now = self.datetime_now()
         hier = self.getHierarchy(
             dim,
             interp_order,
             refinement_boxes=None,
             qty="b",
+            nbr_part_per_cell=ppc,
+            diag_outputs=f"test_b/{dim}/{interp_order}/{self.ddt_test_id()}",
             **kwargs,
         )
+        print(
+            f"\n{self._testMethodName}_{dim}d init took {self.datetime_diff(now)} seconds"
+        )
+        now = self.datetime_now()
 
         from pyphare.pharein import global_vars
 
@@ -320,16 +326,44 @@ def _test_B_is_as_provided_by_user(self, dim, interp_order, **kwargs):
                     )
 
                 if dim == 3:
-                    raise ValueError("Unsupported dimension")
+                    zbx = bx_pd.z[:]
+                    zby = by_pd.z[:]
+                    zbz = bz_pd.z[:]
 
-    def _test_bulkvel_is_as_provided_by_user(self, dim, interp_order):
+                    xbx, ybx, zbx = [
+                        a.flatten() for a in np.meshgrid(xbx, ybx, zbx, indexing="ij")
+                    ]
+                    xby, yby, zby = [
+                        a.flatten() for a in np.meshgrid(xby, yby, zby, indexing="ij")
+                    ]
+                    xbz, ybz, zbz = [
+                        a.flatten() for a in np.meshgrid(xbz, ybz, zbz, indexing="ij")
+                    ]
+
+                    np.testing.assert_allclose(
+                        bx, bx_fn(xbx, ybx, zbx), atol=1e-16, rtol=0
+                    )
+                    np.testing.assert_allclose(
+                        by, by_fn(xby, yby, zby).reshape(by.shape), atol=1e-16, rtol=0
+                    )
+                    np.testing.assert_allclose(
+                        bz, bz_fn(xbz, ybz, zbz).reshape(bz.shape), atol=1e-16, rtol=0
+                    )
+
+        print(f"\n{self._testMethodName}_{dim}d took {self.datetime_diff(now)} seconds")
+
+    def _test_bulkvel_is_as_provided_by_user(
+        self, dim, interp_order, ppc=100, **kwargs
+    ):
         hier = self.getHierarchy(
             dim,
             interp_order,
             {"L0": {"B0": nDBox(dim, 10, 19)}},
             "moments",
-            nbr_part_per_cell=100,
+            nbr_part_per_cell=ppc,
             beam=True,
+            diag_outputs=f"test_bulkV/{dim}/{interp_order}/{self.ddt_test_id()}",
+            **kwargs,
         )
 
         from pyphare.pharein import global_vars
@@ -347,117 +381,58 @@ def _test_bulkvel_is_as_provided_by_user(self, dim, interp_order):
             for ip, patch in enumerate(level.patches):
                 print("patch {}".format(ip))
 
-                layout = patch.patch_datas["protons_Fx"].layout
-                centering = layout.centering["X"][
-                    patch.patch_datas["protons_Fx"].field_name
-                ]
-                nbrGhosts = layout.nbrGhosts(interp_order, centering)
-
-                if dim == 1:
-                    x = patch.patch_datas["protons_Fx"].x[nbrGhosts:-nbrGhosts]
-
-                    fpx = patch.patch_datas["protons_Fx"].dataset[nbrGhosts:-nbrGhosts]
-                    fpy = patch.patch_datas["protons_Fy"].dataset[nbrGhosts:-nbrGhosts]
-                    fpz = patch.patch_datas["protons_Fz"].dataset[nbrGhosts:-nbrGhosts]
-
-                    fbx = patch.patch_datas["beam_Fx"].dataset[nbrGhosts:-nbrGhosts]
-                    fby = patch.patch_datas["beam_Fy"].dataset[nbrGhosts:-nbrGhosts]
-                    fbz = patch.patch_datas["beam_Fz"].dataset[nbrGhosts:-nbrGhosts]
-
-                    ni = patch.patch_datas["rho"].dataset[nbrGhosts:-nbrGhosts]
-
-                    vxact = (fpx + fbx) / ni
-                    vyact = (fpy + fby) / ni
-                    vzact = (fpz + fbz) / ni
-
-                    vxexp = (nprot(x) * vx_fn(x) + nbeam(x) * vx_fn(x)) / (
-                        nprot(x) + nbeam(x)
-                    )
-                    vyexp = (nprot(x) * vy_fn(x) + nbeam(x) * vy_fn(x)) / (
-                        nprot(x) + nbeam(x)
-                    )
-                    vzexp = (nprot(x) * vz_fn(x) + nbeam(x) * vz_fn(x)) / (
-                        nprot(x) + nbeam(x)
-                    )
-
-                    for vexp, vact in zip((vxexp, vyexp, vzexp), (vxact, vyact, vzact)):
-                        std = np.std(vexp - vact)
-                        print("sigma(user v - actual v) = {}".format(std))
-                        self.assertTrue(
-                            std < 1e-2
-                        )  # empirical value obtained from print just above
-
-                def reshape(patch_data, nGhosts):
+                pdatas = patch.patch_datas
+                layout = pdatas["protons_Fx"].layout
+                centering = layout.centering["X"][pdatas["protons_Fx"].field_name]
+                nbrGhosts = layout.nbrGhosts(
+                    interp_order, centering
+                )  # primal in all directions
+                select = tuple([slice(nbrGhosts, -nbrGhosts) for i in range(dim)])
+
+                def domain(patch_data):
+                    if dim == 1:
+                        return patch_data.dataset[select]
                     return patch_data.dataset[:].reshape(
-                        patch.box.shape + (nGhosts * 2) + 1
-                    )
-
-                if dim == 2:
-                    xx, yy = np.meshgrid(
-                        patch.patch_datas["protons_Fx"].x,
-                        patch.patch_datas["protons_Fx"].y,
-                        indexing="ij",
-                    )
-
-                    density = reshape(patch.patch_datas["rho"], nbrGhosts)
-
-                    protons_Fx = reshape(patch.patch_datas["protons_Fx"], nbrGhosts)
-                    protons_Fy = reshape(patch.patch_datas["protons_Fy"], nbrGhosts)
-                    protons_Fz = reshape(patch.patch_datas["protons_Fz"], nbrGhosts)
-
-                    beam_Fx = reshape(patch.patch_datas["beam_Fx"], nbrGhosts)
-                    beam_Fy = reshape(patch.patch_datas["beam_Fy"], nbrGhosts)
-                    beam_Fz = reshape(patch.patch_datas["beam_Fz"], nbrGhosts)
-
-                    x = xx[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
-                    y = yy[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
-
-                    fpx = protons_Fx[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
-                    fpy = protons_Fy[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
-                    fpz = protons_Fz[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
-
-                    fbx = beam_Fx[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
-                    fby = beam_Fy[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
-                    fbz = beam_Fz[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
-
-                    ni = density[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
-
-                    vxact = (fpx + fbx) / ni
-                    vyact = (fpy + fby) / ni
-                    vzact = (fpz + fbz) / ni
-
-                    vxexp = (nprot(x, y) * vx_fn(x, y) + nbeam(x, y) * vx_fn(x, y)) / (
-                        nprot(x, y) + nbeam(x, y)
-                    )
-                    vyexp = (nprot(x, y) * vy_fn(x, y) + nbeam(x, y) * vy_fn(x, y)) / (
-                        nprot(x, y) + nbeam(x, y)
-                    )
-                    vzexp = (nprot(x, y) * vz_fn(x, y) + nbeam(x, y) * vz_fn(x, y)) / (
-                        nprot(x, y) + nbeam(x, y)
-                    )
-
-                    for vexp, vact in zip((vxexp, vyexp, vzexp), (vxact, vyact, vzact)):
-                        self.assertTrue(np.std(vexp - vact) < 1e-2)
+                        patch.box.shape + (nbrGhosts * 2) + 1
+                    )[select]
+
+                ni = domain(pdatas["rho"])
+                vxact = (domain(pdatas["protons_Fx"]) + domain(pdatas["beam_Fx"])) / ni
+                vyact = (domain(pdatas["protons_Fy"]) + domain(pdatas["beam_Fy"])) / ni
+                vzact = (domain(pdatas["protons_Fz"]) + domain(pdatas["beam_Fz"])) / ni
+
+                select = pdatas["protons_Fx"].meshgrid(select)
+                vxexp = (
+                    nprot(*select) * vx_fn(*select) + nbeam(*select) * vx_fn(*select)
+                ) / (nprot(*select) + nbeam(*select))
+                vyexp = (
+                    nprot(*select) * vy_fn(*select) + nbeam(*select) * vy_fn(*select)
+                ) / (nprot(*select) + nbeam(*select))
+                vzexp = (
+                    nprot(*select) * vz_fn(*select) + nbeam(*select) * vz_fn(*select)
+                ) / (nprot(*select) + nbeam(*select))
+                for vexp, vact in zip((vxexp, vyexp, vzexp), (vxact, vyact, vzact)):
+                    self.assertTrue(np.std(vexp - vact) < 1e-2)
+
+    def _test_density_is_as_provided_by_user(self, ndim, interp_order, **kwargs):
+        print(
+            f"test_density_is_as_provided_by_user : dim {ndim} interp_order {interp_order}"
+        )
 
-    def _test_density_is_as_provided_by_user(self, dim, interp_order):
         empirical_dim_devs = {
             1: 6e-3,
             2: 3e-2,
+            3: 2e-1,
         }
-
-        nbParts = {1: 10000, 2: 1000}
-        print(
-            "test_density_is_as_provided_by_user : interp_order : {}".format(
-                interp_order
-            )
-        )
+        nbParts = {1: 10000, 2: 1000, 3: 20}
         hier = self.getHierarchy(
-            dim,
+            ndim,
             interp_order,
-            {"L0": {"B0": nDBox(dim, 10, 20)}},
+            {"L0": {"B0": nDBox(ndim, 5, 14)}},
             qty="moments",
-            nbr_part_per_cell=nbParts[dim],
+            nbr_part_per_cell=nbParts[ndim],
             beam=True,
+            **kwargs,
         )
 
         from pyphare.pharein import global_vars
@@ -479,34 +454,16 @@ def _test_density_is_as_provided_by_user(self, dim, interp_order):
                 layout = patch.patch_datas["rho"].layout
                 centering = layout.centering["X"][patch.patch_datas["rho"].field_name]
                 nbrGhosts = layout.nbrGhosts(interp_order, centering)
+                select = tuple([slice(nbrGhosts, -nbrGhosts) for i in range(ndim)])
 
-                if dim == 1:
-                    protons_expected = proton_density_fn(x[nbrGhosts:-nbrGhosts])
-                    beam_expected = beam_density_fn(x[nbrGhosts:-nbrGhosts])
-                    ion_expected = protons_expected + beam_expected
+                mesh = patch.patch_datas["rho"].meshgrid(select)
+                protons_expected = proton_density_fn(*mesh)
+                beam_expected = beam_density_fn(*mesh)
+                ion_expected = protons_expected + beam_expected
 
-                    ion_actual = ion_density[nbrGhosts:-nbrGhosts]
-                    beam_actual = beam_density[nbrGhosts:-nbrGhosts]
-                    protons_actual = proton_density[nbrGhosts:-nbrGhosts]
-
-                if dim == 2:
-                    y = patch.patch_datas["rho"].y
-                    xx, yy = np.meshgrid(x, y, indexing="ij")
-
-                    x0 = xx[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
-                    y0 = yy[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
-
-                    protons_expected = proton_density_fn(x0, y0)
-                    beam_expected = beam_density_fn(x0, y0)
-                    ion_expected = protons_expected + beam_expected
-
-                    ion_actual = ion_density[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
-                    beam_actual = beam_density[
-                        nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts
-                    ]
-                    protons_actual = proton_density[
-                        nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts
-                    ]
+                protons_actual = proton_density[select]
+                beam_actual = beam_density[select]
+                ion_actual = ion_density[select]
 
                 names = ("ions", "protons", "beam")
                 expected = (ion_expected, protons_expected, beam_expected)
@@ -519,11 +476,11 @@ def _test_density_is_as_provided_by_user(self, dim, interp_order):
                 for name, dev in devs.items():
                     print(f"sigma(user density - {name} density) = {dev}")
                     self.assertLess(
-                        dev, empirical_dim_devs[dim], f"{name} has dev = {dev}"
+                        dev, empirical_dim_devs[ndim], f"{name} has dev = {dev}"
                     )
 
     def _test_density_decreases_as_1overSqrtN(
-        self, dim, interp_order, nbr_particles=None, cells=960
+        self, ndim, interp_order, nbr_particles=None, cells=960
     ):
         import matplotlib.pyplot as plt
 
@@ -536,12 +493,12 @@ def _test_density_decreases_as_1overSqrtN(
 
         for inbr, nbrpart in enumerate(nbr_particles):
             hier = self.getHierarchy(
-                dim,
+                ndim,
                 interp_order,
                 None,
                 "moments",
                 nbr_part_per_cell=nbrpart,
-                diag_outputs=f"{nbrpart}",
+                diag_outputs=f"1overSqrtN/{ndim}/{interp_order}/{nbrpart}",
                 density=lambda *xyz: np.zeros(tuple(_.shape[0] for _ in xyz)) + 1.0,
                 smallest_patch_size=int(cells / 2),
                 largest_patch_size=int(cells / 2),
@@ -559,7 +516,7 @@ def _test_density_decreases_as_1overSqrtN(
 
             centering = layout.centering["X"][patch.patch_datas["rho"].field_name]
             nbrGhosts = layout.nbrGhosts(interp_order, centering)
-            select = tuple([slice(nbrGhosts, -nbrGhosts) for i in range(dim)])
+            select = tuple([slice(nbrGhosts, -nbrGhosts) for i in range(ndim)])
             ion_density = patch.patch_datas["rho"].dataset[:]
             mesh = patch.patch_datas["rho"].meshgrid(select)
 
@@ -568,14 +525,14 @@ def _test_density_decreases_as_1overSqrtN(
             noise[inbr] = np.std(expected - actual)
             print(f"noise is {noise[inbr]} for {nbrpart} particles per cell")
 
-            if dim == 1:
+            if ndim == 1:
                 x = patch.patch_datas["rho"].x
                 plt.figure()
-                plt.plot(x[nbrGhosts:-nbrGhosts], actual, label="actual")
-                plt.plot(x[nbrGhosts:-nbrGhosts], expected, label="expected")
+                plt.plot(x[select], actual, label="actual")
+                plt.plot(x[select], expected, label="expected")
                 plt.legend()
                 plt.title(r"$\sigma =$ {}".format(noise[inbr]))
-                plt.savefig(f"noise_{nbrpart}_interp_{dim}_{interp_order}.png")
+                plt.savefig(f"noise_{nbrpart}_interp_{ndim}_{interp_order}.png")
                 plt.close("all")
 
         plt.figure()
@@ -587,7 +544,7 @@ def _test_density_decreases_as_1overSqrtN(
         )
         plt.xlabel("nbr_particles")
         plt.legend()
-        plt.savefig(f"noise_nppc_interp_{dim}_{interp_order}.png")
+        plt.savefig(f"noise_nppc_interp_{ndim}_{interp_order}.png")
         plt.close("all")
 
         noiseMinusTheory = noise / noise[0] - 1 / np.sqrt(
@@ -601,25 +558,32 @@ def _test_density_decreases_as_1overSqrtN(
         )
         plt.xlabel("nbr_particles")
         plt.legend()
-        plt.savefig(f"noise_nppc_minus_theory_interp_{dim}_{interp_order}.png")
+        plt.savefig(f"noise_nppc_minus_theory_interp_{ndim}_{interp_order}.png")
         plt.close("all")
         self.assertGreater(3e-2, noiseMinusTheory[1:].mean())
 
     def _test_nbr_particles_per_cell_is_as_provided(
-        self, dim, interp_order, default_ppc=100
+        self, ndim, interp_order, ppc=100, **kwargs
     ):
+        ddt_test_id = self.ddt_test_id()
         datahier = self.getHierarchy(
-            dim,
+            ndim,
             interp_order,
-            {"L0": {"B0": nDBox(dim, 10, 20)}},
+            {},
             "particles",
+            diag_outputs=f"ppc/{ndim}/{interp_order}/{ddt_test_id}",
+            nbr_part_per_cell=ppc,
+            **kwargs,
         )
 
-        for patch in datahier.level(0).patches:
+        if cpp.mpi_rank() > 0:
+            return
+
+        for pi, patch in enumerate(datahier.level(0).patches):
             pd = patch.patch_datas["protons_particles"]
             icells = pd.dataset[patch.box].iCells
-            H, _ = np.histogramdd(icells)
-            self.assertTrue((H == default_ppc).all())
+            H, edges = np.histogramdd(icells, bins=patch.box.shape)
+            self.assertTrue((H == ppc).all())
 
     def _domainParticles_for(self, datahier, ilvl):
         patch0 = datahier.levels()[ilvl].patches[0]
@@ -654,6 +618,9 @@ def _test_domainparticles_have_correct_split_from_coarser_particle(
             **kwargs,
         )
 
+        if cpp.mpi_rank() > 0:
+            return
+
         from pyphare.pharein.global_vars import sim
 
         assert sim is not None and len(sim.cells) == ndim
@@ -691,15 +658,15 @@ def _test_domainparticles_have_correct_split_from_coarser_particle(
                     part2 = coarse_split_particles[pop_name].select(patch.box)
                     self.assertEqual(part1, part2)
 
-    def _test_patch_ghost_on_refined_level_case(self, dim, has_patch_ghost, **kwargs):
+    def _test_patch_ghost_on_refined_level_case(self, ndim, has_patch_ghost, **kwargs):
         import pyphare.pharein as ph
 
         out = "phare_outputs"
-        refinement_boxes = {"L0": [nDBox(dim, 10, 19)]}
+        refinement_boxes = {"L0": [nDBox(ndim, 10, 19)]}
         kwargs["interp_order"] = kwargs.get("interp_order", 1)
         kwargs["diag_outputs"] = f"{has_patch_ghost}"
         datahier = self.getHierarchy(
-            ndim=dim,
+            ndim,
             refinement_boxes=refinement_boxes,
             qty="particles_patch_ghost",
             **kwargs,
@@ -722,12 +689,12 @@ def _test_patch_ghost_on_refined_level_case(self, dim, has_patch_ghost, **kwargs
     def _test_levelghostparticles_have_correct_split_from_coarser_particle(
         self, datahier
     ):
-        dim = datahier.level(0).patches[0].box.ndim
+        ndim = datahier.level(0).patches[0].box.ndim
 
         from pyphare.pharein.global_vars import sim
 
         assert sim is not None
-        assert len(sim.cells) == dim
+        assert len(sim.cells) == ndim
 
         particle_level_ghost_boxes_per_level = level_ghost_boxes(datahier, "particles")
 
diff --git a/tools/python3/cmake.py b/tools/python3/cmake.py
new file mode 100644
index 000000000..3cd20a65a
--- /dev/null
+++ b/tools/python3/cmake.py
@@ -0,0 +1,48 @@
+from tools.python3 import decode_bytes, run
+
+
+def version():
+    pass
+
+
+def make_config_str(
+    path, samrai_dir=None, cxx_flags=None, use_ninja=False, use_ccache=False, extra=""
+):
+    """
+    FULL UNPORTABLE OPTIMIZATIONS = cxx_flags="-O3 -march=native -mtune=native"
+    """
+
+    samrai_dir = "" if samrai_dir is None else f"-DSAMRAI_ROOT={samrai_dir}"
+    cxx_flags = "" if cxx_flags is None else f'-DCMAKE_CXX_FLAGS="{cxx_flags}"'
+    ccache = "" if use_ccache is False else "-DCMAKE_CXX_COMPILER_LAUNCHER=ccache"
+    ninja = "" if not use_ninja else "-G Ninja"
+
+    return f"cmake {path} {samrai_dir} {cxx_flags} {ninja} {ccache} {extra}"
+
+
+def config(
+    path, samrai_dir=None, cxx_flags=None, use_ninja=False, use_ccache=False, extra=""
+):
+    cmd = make_config_str(path, samrai_dir, cxx_flags, use_ninja, extra)
+    run(cmd, capture_output=False)
+
+
+def build(use_ninja=False, threads=1):
+    run("ninja" if use_ninja else f"make -j{threads}", capture_output=False)
+
+
+def list_tests():
+    proc = run("ctest -N", capture_output=True)
+    out = decode_bytes(proc.stdout)
+    return [line.split(" ")[-1] for line in out.splitlines()[1:-2]]
+
+
+def test_cmd(test, verbose=False):
+    cmd = f"ctest -R {test}"
+    if verbose:
+        cmd = f"{cmd} -V"
+    return cmd
+
+
+def run_test(test, verbose=False, capture_output=False):
+    run(test_cmd(cmd, verbose=verbose), capture_output=capture_output)
diff --git a/tools/python3/git.py b/tools/python3/git.py
new file mode 100644
index 000000000..9e081d57c
--- /dev/null
+++ b/tools/python3/git.py
@@ -0,0 +1,35 @@
+from tools.python3 import decode_bytes, run
+import subprocess
+
+
+def current_branch():
+    return decode_bytes(run("git branch --show-current").stdout)
+
+
+def hashes(N=1):
+    return decode_bytes(run(f"git log -{N} --pretty=format:%h").stdout).splitlines()
+
+
+def log(N=10, use_short=False):
+    # https://git-scm.com/docs/pretty-formats
+    short = '--pretty=format:"%h%x09%an%x09%ad%x09%s"' if use_short else ""
+    return decode_bytes(run(f"git log -{N} {short}").stdout)
+
+
+def branch_exists(branch):
+    try:
+        run(f"git show-branch {branch}", check=True)
+    except subprocess.CalledProcessError as e:
+        return False  # exit failure means branch does not exist
+    return True
+
+
+def checkout(branch, create=False, recreate=False):
+    if recreate:
+        delete_branch(branch)
+        create = True
+
+    if create and not branch_exists(branch):
+        run(f"git checkout -b {branch}")
+    else:
+        run(f"git checkout {branch}")