diff --git a/pyphare/pyphare/core/box.py b/pyphare/pyphare/core/box.py
index 15d5e571d..9c6201292 100644
--- a/pyphare/pyphare/core/box.py
+++ b/pyphare/pyphare/core/box.py
@@ -67,6 +67,9 @@ def __eq__(self, other):
         return isinstance(other, Box) and (self.lower == other.lower).all() and (self.upper == other.upper).all()
 
     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)
 
@@ -183,5 +186,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 dd8d9c379..1663af67e 100644
--- a/pyphare/pyphare/core/gridlayout.py
+++ b/pyphare/pyphare/core/gridlayout.py
@@ -276,6 +276,22 @@ 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):
         """
         from a qty and a direction, returns a 1d array containing
diff --git a/pyphare/pyphare/pharein/maxwellian_fluid_model.py b/pyphare/pyphare/pharein/maxwellian_fluid_model.py
index 9356d742e..4311b0679 100644
--- a/pyphare/pyphare/pharein/maxwellian_fluid_model.py
+++ b/pyphare/pyphare/pharein/maxwellian_fluid_model.py
@@ -142,104 +142,54 @@ def to_dict(self):
 #------------------------------------------------------------------------------
 
     def validate(self, sim):
-        import itertools
         import numpy as np
         from ..core.gridlayout import GridLayout, directions as all_directions
         from ..core.box import Box
         from pyphare.pharein import py_fn_wrapper
 
-        dl = np.array(sim.dl)
         directions = all_directions[:sim.ndim]
-
         domain_box = Box([0] * sim.ndim, np.asarray(sim.cells) - 1)
         assert len(sim.origin) == domain_box.ndim
+        ndim = domain_box.ndim
 
         def _primal_directions(qty):
             return [1 if yee_element_is_primal(qty, direction) else 0 for direction in directions]
 
-        def _nbrGhosts(qty, layout, primal_directions):
+        def _nbrGhosts(layout, qty, primal_directions):
             return [
-                layout.nbrGhosts(sim.interp_order, "primal" if primal_directions[dim] else "dual") for dim in range(sim.ndim)
+                layout.nbrGhosts(sim.interp_order, "primal" if primal_directions[dim] else "dual") for dim in range(ndim)
             ]
 
-        def qty_shifts(qty, nbrGhosts):
-          shifts = []
-          for dim in range(sim.ndim):
-              if yee_element_is_primal(qty, directions[dim]):
-                  shifts += [np.arange(-nbrGhosts[dim], nbrGhosts[dim] + 1) * dl[dim]]
-              else:
-                  shifts += [np.arange(-nbrGhosts[dim], nbrGhosts[dim]) *  dl[dim] + (.5 * dl[dim])]
-          return shifts
+        def split_to_columns(ar):
+            ar = ar.reshape(-1, ndim) # drop outer arrays if any
+            return [ar[:,dim] for dim in range(ndim)]
 
-        def _1d_check(layout, nbrGhosts, primal_directions, qty, fn):
-            fnX = py_fn_wrapper(fn)(layout.yeeCoordsFor(qty, "x"))
-            include_border = primal_directions[0]
-            np.testing.assert_allclose(fnX[:nbrGhosts * 2 + include_border] , fnX[-(nbrGhosts * 2) - include_border:], atol=1e-15)
-            return (
-                np.allclose(fnX[:nbrGhosts * 2 + include_border] , fnX[-(nbrGhosts * 2) - include_border:], atol=1e-15)
-            )
+        def _nd_check(layout, qty, nbrGhosts, fn):
+            from pyphare.pharesee.geometry import hierarchy_overlaps
+            from pyphare.pharesee.hierarchy import hierarchy_from_box
+            from ..core import box as boxm
+
+            hier = hierarchy_from_box(domain_box, nbrGhosts)
+            coords = layout.meshCoords(qty)
+
+            for ilvl, overlaps in hierarchy_overlaps(hier).items():
+                for overlap in overlaps:
+                    (pd1, pd2), box, offsets = overlap["pdatas"], overlap["box"], overlap["offset"]
+                    slice1 = boxm.select(coords, boxm.amr_to_local(box, boxm.shift(pd1.ghost_box, offsets[0])))
+                    slice2 = boxm.select(coords, boxm.amr_to_local(box, boxm.shift(pd2.ghost_box, offsets[1])))
+                    if not np.allclose(fn(*split_to_columns(slice1)) , fn(*split_to_columns(slice2)), atol=1e-15):
+                        return False
+            return True
 
 
-        def split_to_columns(ar):
-            ar = ar.reshape(-1, sim.ndim) # drop outer arrays if any
-            return [ar[:,dim] for dim in range(sim.ndim)]
-
-
-        def _2d_check(nbrGhosts, primal_directions, qty, fn, shift):
-            layout = GridLayout(domain_box, sim.origin + (shift * sim.dl), sim.dl, sim.interp_order)
-            x = layout.yeeCoordsFor(qty, "x")[nbrGhosts[0]:-(nbrGhosts[0]-1)-(primal_directions[0])]
-            y = layout.yeeCoordsFor(qty, "y")[nbrGhosts[1]:-(nbrGhosts[1]-1)-(primal_directions[1])]
-            X,Y = np.meshgrid(x,y,indexing="ij")
-            coords = np.array([X.flatten(),Y.flatten()]).T
-            top = coords[:len(x)]
-            bot = coords[-len(x):]
-            left  = coords[0 :: len(x)]
-            right = coords[len(x) - 1 :: len(x)]
-            return (
-                np.allclose(fn(*split_to_columns(left)) , fn(*split_to_columns(right)), atol=1e-15) and
-                np.allclose(fn(*split_to_columns(top))  , fn(*split_to_columns(bot)), atol=1e-15)
-            )
-
-        def _3d_check(nbrGhosts, primal_directions, qty, fn, shift):
-            def get_3d_faces(M):
-                def get_face(M, dim, front_side):
-                    return M[tuple((0 if front_side else -1) if i == dim else slice(None) for i in range(M.ndim))]
-                return [get_face(M, dim, front_side) for dim in range(sim.ndim) for front_side in (True, False)]
-
-            layout = GridLayout(domain_box, sim.origin + (shift * sim.dl), sim.dl, sim.interp_order)
-            x = layout.yeeCoordsFor(qty, "x")[nbrGhosts[0]:-(nbrGhosts[0]-1)-(primal_directions[0])]
-            y = layout.yeeCoordsFor(qty, "y")[nbrGhosts[1]:-(nbrGhosts[1]-1)-(primal_directions[1])]
-            z = layout.yeeCoordsFor(qty, "z")[nbrGhosts[2]:-(nbrGhosts[2]-1)-(primal_directions[2])]
-            coords = np.array([X.flatten(), Y.flatten(), Z.flatten()]).T.reshape((len(x), len(y), len(z), sim.ndim))
-            left, right, top, bot, front, back = get_3d_faces(coords)
-            return (
-                np.allclose(fn(split_to_columns(left)) , fn(split_to_columns(right)), atol=1e-15) and
-                np.allclose(fn(split_to_columns(top))  , fn(split_to_columns(bot)), atol=1e-15) and
-                np.allclose(fn(split_to_columns(front)), fn(split_to_columns(back)), atol=1e-15)
-            )
 
         def periodic_function_check(vec_field, dic):
             valid = True
             for xyz in ["x", "y", "z"]:
                 qty = vec_field + xyz
-                fn = dic[qty]
-
                 layout = GridLayout(domain_box, sim.origin, sim.dl, sim.interp_order)
-                primal_directions = _primal_directions(qty)
-                nbrGhosts = _nbrGhosts(qty, layout, primal_directions)
-
-                if sim.ndim == 1:
-                    valid &= _1d_check(layout, nbrGhosts[0], primal_directions, qty, fn)
-
-                if sim.ndim == 2:
-                    permutations = itertools.product(*qty_shifts(qty, nbrGhosts))
-                    valid &= all([_2d_check(nbrGhosts, primal_directions, qty, fn, np.asarray(perm)) for perm in permutations])
-
-                if sim.ndim == 3:
-                    # untested block - add in during 3d/splitting PR https://github.com/PHAREHUB/PHARE/pull/314
-                    permutations = itertools.product(*qty_shifts(qty, nbrGhosts))
-                    valid &= all([_3d_check(nbrGhosts, primal_directions, qty, fn, np.asarray(perm)) for perm in permutations])
-
+                nbrGhosts = _nbrGhosts(layout, qty, _primal_directions(qty))
+                valid &= _nd_check(layout, qty, nbrGhosts, fn=dic[qty])
             return valid
 
         if sim.boundary_types[0] == "periodic":
diff --git a/pyphare/pyphare/pharein/simulation.py b/pyphare/pyphare/pharein/simulation.py
index 1ea8349a7..c8b8ccf45 100644
--- a/pyphare/pyphare/pharein/simulation.py
+++ b/pyphare/pyphare/pharein/simulation.py
@@ -194,9 +194,13 @@ def check_boundaries(ndim, **kwargs):
     3: [4, 5, 8, 9, 25]
   },
   3: {
-    1: [6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 26, 27],
-    2: [6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 26, 27, 64],
-    3: [6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 26, 27, 125]
+    1: [6, 12],
+    2: [6, 12],
+    3: [6, 12]
+    # full below
+    # 1: [6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 26, 27],
+    # 2: [6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 26, 27, 64],
+    # 3: [6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 26, 27, 125]
   }
 } # Default refined_particle_nbr per dim/interp is considered index 0 of list
 def check_refined_particle_nbr(ndim, **kwargs):
diff --git a/pyphare/pyphare/pharesee/geometry.py b/pyphare/pyphare/pharesee/geometry.py
index d259fa8a6..fa75c7745 100644
--- a/pyphare/pyphare/pharesee/geometry.py
+++ b/pyphare/pyphare/pharesee/geometry.py
@@ -49,13 +49,25 @@ 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)),
+        }
+
+    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)),
         }
 
-    raise ValueError("Unhandeled dimension")
+
+
 
 def touch_domain_border(box, domain_box, border):
     if border == "upper":
@@ -70,13 +82,18 @@ def periodicity_shifts(domain_box):
 
     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)],
@@ -104,13 +121,33 @@ def periodicity_shifts(domain_box):
           "topleftright" : [
               *shifts["leftright"], *shifts["top"],
               shifts["topleft"][-1], shifts["topright"][-1]],
-          "bottomtopleftright" : [ # one patch covers domain
+          "leftrightbottomtop" : [ # one patch covers domain
               *shifts["bottomleft"], *shifts["topright"],
               shifts["bottomright"][-1], shifts["topleft"][-1]]
         })
 
     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"],
+            ]
+        })
+
+    shifts = {"".join(sorted(k)) : l for k,l in shifts.items()}
 
     return shifts
 
@@ -206,7 +243,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():
@@ -249,6 +286,7 @@ def hierarchy_overlaps(hierarchy, time=0):
 
 
 
+
 def get_periodic_list(patches, domain_box, n_ghosts):
     """
     given a list of patches and a domain box the function
@@ -278,38 +316,50 @@ def get_periodic_list(patches, domain_box, n_ghosts):
             shift_patch(first_patch, domain_box.shape)
             sorted_patches.append(first_patch)
 
-    if dim == 2:
+        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:
+    def borders_per(box):
+        return "".join([key for key, side in sides.items() if box * side is not None])
 
-            in_sides = borders_per(boxm.grow(patch.box, n_ghosts))
+    for patch in patches:
 
-            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)
+        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
 
 
 
 
+
+
 def ghost_area_boxes(hierarchy, quantities, levelNbrs=[], time=0):
     """
     this function returns boxes representing ghost cell boxes for all levels
@@ -412,18 +462,7 @@ def level_ghost_boxes(hierarchy, quantities, levelNbrs=[], time=None):
 
                 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
index 7e109cee9..61be25d59 100644
--- a/pyphare/pyphare/pharesee/hierarchy.py
+++ b/pyphare/pyphare/pharesee/hierarchy.py
@@ -87,13 +87,7 @@ def select(self, box):
 
         overlap = box * gbox
         if overlap is not None:
-            lower = self.layout.AMRToLocal(overlap.lower)
-            upper = self.layout.AMRToLocal(overlap.upper)
-
-            if box.ndim == 1:
-                return self.dataset[lower[0] : upper[0] + 1]
-            if box.ndim == 2:
-                return self.dataset[lower[0]:upper[0] + 1 , lower[1] : upper[1] + 1]
+            return boxm.select(self.dataset, self.layout.AMRBoxToLocal(overlap))
         return np.array([])
 
     def __getitem__(self, box):
@@ -151,6 +145,19 @@ def __init__(self, layout, field_name, data, **kwargs):
 
 
 
+    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):
     """
@@ -1431,6 +1438,18 @@ def hierarchy_from_sim(simulator, qty, pop=""):
     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), [.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):
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/def.cmake b/res/cmake/def.cmake
index f7f2d3d47..ccbfeafa3 100644
--- a/res/cmake/def.cmake
+++ b/res/cmake/def.cmake
@@ -253,3 +253,19 @@ if (test AND ${PHARE_EXEC_LEVEL_MIN} GREATER 0) # 0 = no tests
 
 endif()
 
+# -DwithPGO_GEN
+if (withPGO_GEN)
+  set (PHARE_PGO_FLAGS -fprofile-generate)
+endif(withPGO_GEN)
+
+# -DwithPGO_USE
+if (withPGO_USE)
+  set (PHARE_PGO_FLAGS -fprofile-use)
+endif(withPGO_USE)
+
+if (withPGO_GEN OR withPGO_USE)
+  set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${PHARE_PGO_FLAGS}" )
+  set (CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} ${PHARE_PGO_FLAGS}"  )
+  set (CMAKE_SHARED_LINKER_FLAGS "${CMAKE_SHARED_LINKER_FLAGS} ${PHARE_PGO_FLAGS}"  )
+endif() # (withPGO_GEN OR withPGO_USE)
+
diff --git a/res/cmake/options.cmake b/res/cmake/options.cmake
index 66efd8fc5..c6a7c2b53 100644
--- a/res/cmake/options.cmake
+++ b/res/cmake/options.cmake
@@ -66,10 +66,21 @@ option(withCaliper "Use LLNL Caliper" 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" ON)
+
 # -DtestDuringBuild=ON enabled if devMode=ON, disabled if asan=ON (needs LD_PRELOAD)
 option(testDuringBuild "Runs C++ unit tests after they are built" OFF)
 
 
+
+
+# -DwithPGO_GEN
+option(withPGO_GEN "Activate generate step of PGO")
+# -DwithPGO_A
+option(withPGO_USE "Activate usage step of PGO")
+
+
 # Controlling the activation of tests
 if (NOT DEFINED PHARE_EXEC_LEVEL_MIN)
   set(PHARE_EXEC_LEVEL_MIN 1)
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 ce6b8dba0..da1f95cf6 100644
--- a/src/amr/data/particles/refine/splitter.hpp
+++ b/src/amr/data/particles/refine/splitter.hpp
@@ -55,6 +55,9 @@ class PatternDispatcher
     template<typename Particle, typename Particles>
     void dispatch(Particle const& particle, Particles& particles, size_t idx) const
     {
+        using Weight_t = std::decay_t<decltype(particle.weight)>;
+        using Delta_t  = std::decay_t<decltype(particle.delta[0])>;
+
         constexpr auto dimension = Particle::dimension;
         constexpr auto refRatio  = PHARE::amr::refinementRatio;
         constexpr std::array power{refRatio, refRatio * refRatio, refRatio * refRatio * refRatio};
@@ -64,19 +67,21 @@ 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++];
-                fineParticle.weight = particle.weight * pattern.weight_ * power[dimension - 1];
-                fineParticle.charge = particle.charge;
-                fineParticle.iCell  = particle.iCell;
-                fineParticle.delta  = particle.delta;
-                fineParticle.v      = particle.v;
+                fineParticle.weight       = particle.weight * weight * power[dimension - 1];
+                fineParticle.charge       = particle.charge;
+                fineParticle.iCell        = particle.iCell;
+                fineParticle.delta        = particle.delta;
+                fineParticle.v            = particle.v;
 
                 for (size_t iDim = 0; iDim < dimension; iDim++)
                 {
-                    fineParticle.delta[iDim] += pattern.deltas_[rpIndex][iDim];
-                    float integra = std::floor(fineParticle.delta[iDim]);
+                    fineParticle.delta[iDim]
+                        += static_cast<Delta_t>(pattern.deltas_[rpIndex][iDim]);
+                    Delta_t integra = std::floor(fineParticle.delta[iDim]);
                     fineParticle.delta[iDim] -= integra;
                     fineParticle.iCell[iDim] += static_cast<int32_t>(integra);
                 }
@@ -109,25 +114,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 PurplePattern
+{
+};
+template<typename dim>
+struct BrownPattern
 {
 };
 template<typename dim>
-struct BrownDispatcher
+struct LimePattern
 {
 };
 template<typename dim>
-struct PinkDispatcher
+struct WhitePattern
 {
 };
 
diff --git a/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp b/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp
index 0f562703d..812202716 100644
--- a/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp
+++ b/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp
@@ -138,6 +138,8 @@ namespace amr
         void registerLevel(std::shared_ptr<SAMRAI::hier::PatchHierarchy> const& hierarchy,
                            int const levelNumber) override
         {
+            PHARE_LOG_SCOPE("HybridHybridMessengerStrategy::registerLevel");
+
             auto const level = hierarchy->getPatchLevel(levelNumber);
 
             magneticSharedNodes_.registerLevel(hierarchy, level);
@@ -235,6 +237,8 @@ namespace amr
         void initLevel(IPhysicalModel& model, SAMRAI::hier::PatchLevel& level,
                        double const initDataTime) override
         {
+            PHARE_LOG_SCOPE("HybridHybridMessengerStrategy::initLevel");
+
             auto levelNumber = level.getLevelNumber();
 
             magneticInit_.fill(levelNumber, initDataTime);
diff --git a/src/core/data/grid/gridlayoutimplyee.hpp b/src/core/data/grid/gridlayoutimplyee.hpp
index bf9f52df9..f0394b648 100644
--- a/src/core/data/grid/gridlayoutimplyee.hpp
+++ b/src/core/data/grid/gridlayoutimplyee.hpp
@@ -310,7 +310,7 @@ namespace core
              * the following is only valid when dual and primal do not have the same number of
             ghosts
              * and that depends on the interp order
-             * It is commented out because ghosts are hard coded to 5 for now.
+             * It is commented out because ghosts are hard coded to 3 for now.
              *
             if constexpr (interp_order == 1 || interp_order == 2 || interp_order == 4)
                 return -1;
@@ -330,7 +330,7 @@ namespace core
              * the following is only valid when dual and primal do not have the same number of
             ghosts
              * and that depends on the interp order
-             * It is commented out because ghosts are hard coded to 5 for now.
+             * It is commented out because ghosts are hard coded to 3 for now.
              *
             if constexpr (interp_order == 1 || interp_order == 2 || interp_order == 4)
                 return 1;
diff --git a/src/core/data/ions/particle_initializers/maxwellian_particle_initializer.hpp b/src/core/data/ions/particle_initializers/maxwellian_particle_initializer.hpp
index 9c44ec623..82d9f05f4 100644
--- a/src/core/data/ions/particle_initializers/maxwellian_particle_initializer.hpp
+++ b/src/core/data/ions/particle_initializers/maxwellian_particle_initializer.hpp
@@ -9,10 +9,10 @@
 #include "core/data/grid/gridlayoutdefs.hpp"
 #include "core/hybrid/hybrid_quantities.hpp"
 #include "core/utilities/types.hpp"
+#include "core/utilities/point/point.hpp"
 #include "core/data/ions/particle_initializers/particle_initializer.hpp"
 #include "core/data/particles/particle.hpp"
 #include "initializer/data_provider.hpp"
-#include "core/utilities/point/point.hpp"
 
 
 namespace PHARE::core
diff --git a/src/core/data/ndarray/ndarray_vector.hpp b/src/core/data/ndarray/ndarray_vector.hpp
index c70d4f525..63428d9ac 100644
--- a/src/core/data/ndarray/ndarray_vector.hpp
+++ b/src/core/data/ndarray/ndarray_vector.hpp
@@ -2,6 +2,7 @@
 #define PHARE_CORE_DATA_NDARRAY_NDARRAY_VECTOR_HPP
 
 #include <stdexcept>
+#include <algorithm>
 #include <array>
 #include <cstdint>
 #include <vector>
@@ -109,15 +110,17 @@ class MaskedView
     auto operator=(data_type value) { mask_.fill(array_, value); }
 
     auto xstart() const { return mask_.min(); }
-
     auto xend() const { return shape_[0] - 1 - mask_.max(); }
 
 
     auto ystart() const { return mask_.min(); }
-
     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_;
@@ -308,10 +311,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;
     }
 
@@ -321,24 +324,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;
         }
     }
@@ -346,7 +349,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>
@@ -356,17 +396,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;
     }
 
@@ -386,20 +432,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
@@ -416,7 +463,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());
 
 
@@ -445,7 +492,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 a43563764..82c533d01 100644
--- a/src/core/utilities/box/box.hpp
+++ b/src/core/utilities/box/box.hpp
@@ -6,6 +6,7 @@
 #include "core/utilities/point/point.hpp"
 #include "core/utilities/meta/meta_utilities.hpp"
 
+#include <vector>
 #include <cstddef>
 #include <algorithm>
 #include <optional>
@@ -135,6 +136,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>
@@ -189,7 +210,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..750f66890 100644
--- a/src/core/utilities/meta/meta_utilities.hpp
+++ b/src/core/utilities/meta/meta_utilities.hpp
@@ -78,7 +78,13 @@ namespace core
 
                           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>,
+
+                          // TODO add in the rest of 3d nbrParticles permutations
+                          // possibly consider compile time activation for uncommon cases
+                          SimulatorOption<DimConst<3>, InterpConst<1>, 6, 12 /*, 27*/>,
+                          SimulatorOption<DimConst<3>, InterpConst<2>, 6, 12>,
+                          SimulatorOption<DimConst<3>, InterpConst<3>, 6, 12>>{};
     }
 
 
diff --git a/src/core/utilities/point/point.hpp b/src/core/utilities/point/point.hpp
index b99f823c3..cf7359b3a 100644
--- a/src/core/utilities/point/point.hpp
+++ b/src/core/utilities/point/point.hpp
@@ -177,7 +177,7 @@ namespace core
 
     template<typename... Indexes>
     Point(Indexes... indexes)
-        -> Point<typename std::tuple_element<0, std::tuple<Indexes...>>::type, sizeof...(indexes)>;
+        ->Point<typename std::tuple_element<0, std::tuple<Indexes...>>::type, sizeof...(indexes)>;
 
     template<typename Type, std::size_t dim>
     auto& operator<<(std::ostream& os, Point<Type, dim> const& p)
@@ -203,6 +203,18 @@ PHARE::core::Point<Type, dim> abs(PHARE::core::Point<Type, dim> const& point)
     return postive;
 }
 
+template<typename Type, std::size_t dim>
+auto to_string(PHARE::core::Point<Type, dim> const& point)
+{
+    return point.str();
+}
+
+template<typename Type, std::size_t dim>
+auto to_string(PHARE::core::Point<Type, dim>&& point)
+{
+    return point.str();
+}
+
 } // namespace std
 
 
diff --git a/src/phare/phare.cpp b/src/phare/phare.cpp
index fec4d4c06..16dab12ba 100644
--- a/src/phare/phare.cpp
+++ b/src/phare/phare.cpp
@@ -32,7 +32,8 @@ std::unique_ptr<PHARE::initializer::DataProvider> fromCommandLine(int argc, char
     return nullptr;
 }
 
-int main(int argc, char** argv)
+
+int naim(int argc, char** argv)
 {
     std::string const welcome = R"~(
                   _____   _    _            _____   ______
@@ -76,4 +77,14 @@ int main(int argc, char** argv)
         std::cout << simulator->currentTime() << "\n";
         //    time += simulator.timeStep();
     }
+
+    return 0;
+}
+
+
+#if !defined(PHARE_EXE_NO_MAIN)
+int main(int argc, char** argv)
+{
+    return naim(argc, argv);
 }
+#endif
diff --git a/src/phare/phare_init.py b/src/phare/phare_init.py
index 41aa24932..cb863be3f 100644
--- a/src/phare/phare_init.py
+++ b/src/phare/phare_init.py
@@ -20,7 +20,7 @@
     max_nbr_levels=2,          # (default=1) max nbr of levels in the AMR hierarchy
     refinement = "tagging",
     #refinement_boxes={"L0": {"B0": [(10, ), (20, )]}},
-    diag_options={"format": "phareh5", "options": {"dir": "phare_outputs"}}
+    diag_options={"format": "phareh5", "options": {"dir": "phare_outputs", "mode":"overwrite"}}
 )
 
 
diff --git a/src/python3/CMakeLists.txt b/src/python3/CMakeLists.txt
index c858b68de..8f3f40495 100644
--- a/src/python3/CMakeLists.txt
+++ b/src/python3/CMakeLists.txt
@@ -2,23 +2,37 @@ cmake_minimum_required (VERSION 3.9)
 
 project(phare_python3)
 
+set(PHARE_PYBIND_CXX_FLAGS  -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE})
+set(PHARE_PYBIND_LD_FLAGS )
+
+# set(PHARE_PYBIND_CXX_FLAGS "-DPHARE_EXEC_PYBIND=1 -nostartfiles -Wl,--entry=naim")
+# set(PHARE_PYBIND_LD_FLAGS "${PHARE_PYBIND_CXX_FLAGS}")
+# set(PHARE_PYBIND_CXX_FLAGS "${PHARE_PYBIND_CXX_FLAGS} -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE}")
+
 pybind11_add_module(cpp cpp_simulator.cpp)
 target_link_libraries(cpp PUBLIC phare_simulator)
-target_compile_options(cpp PRIVATE ${PHARE_FLAGS} -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE}) # pybind fails with Werror
+target_compile_options(
+  cpp
+  PRIVATE ${PHARE_FLAGS} ${PHARE_PYBIND_CXX_FLAGS}) # pybind fails with Werror
 set_target_properties(cpp
     PROPERTIES
     LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/pybindlibs"
+    LINK_FLAGS "${PHARE_PYBIND_LD_FLAGS}"
 )
+
 # this is on by default "pybind11_add_module" but can interfere with coverage so we disable it if coverage is enabled
 set_property(TARGET cpp PROPERTY INTERPROCEDURAL_OPTIMIZATION ${PHARE_INTERPROCEDURAL_OPTIMIZATION})
 
 if (CMAKE_BUILD_TYPE STREQUAL "Debug")
   pybind11_add_module(cpp_dbg cpp_simulator.cpp)
   target_link_libraries(cpp_dbg PUBLIC phare_simulator)
-  target_compile_options(cpp_dbg PRIVATE ${PHARE_FLAGS} -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE} -DPHARE_DIAG_DOUBLES=1 -DPHARE_CPP_MOD_NAME=cpp_dbg)
+  target_compile_options(
+    cpp_dbg
+    PRIVATE ${PHARE_FLAGS} ${PHARE_PYBIND_CXX_FLAGS} -DPHARE_DIAG_DOUBLES=1 -DPHARE_CPP_MOD_NAME=cpp_dbg)
   set_target_properties(cpp_dbg
       PROPERTIES
       LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/pybindlibs"
+      LINK_FLAGS "${PHARE_PYBIND_LD_FLAGS}"
   )
   set_property(TARGET cpp_dbg PROPERTY INTERPROCEDURAL_OPTIMIZATION ${PHARE_INTERPROCEDURAL_OPTIMIZATION})
 endif (CMAKE_BUILD_TYPE STREQUAL "Debug")
diff --git a/src/python3/cpp_simulator.cpp b/src/python3/cpp_simulator.cpp
index a5c0a5d1c..03e55b582 100644
--- a/src/python3/cpp_simulator.cpp
+++ b/src/python3/cpp_simulator.cpp
@@ -24,3 +24,9 @@ PYBIND11_MODULE(PHARE_CPP_MOD_NAME, m)
     declarePatchData<std::vector<double>, 3>(m, "PatchDataVectorDouble_3D");
 }
 } // namespace PHARE::pydata
+
+#if defined(PHARE_EXEC_PYBIND)
+#define PHARE_EXE_NO_MAIN 1
+#include "phare/phare.cpp"
+#undef PHARE_EXE_NO_MAIN
+#endif
diff --git a/src/simulator/simulator.hpp b/src/simulator/simulator.hpp
index 95d9242af..25183ecf6 100644
--- a/src/simulator/simulator.hpp
+++ b/src/simulator/simulator.hpp
@@ -78,6 +78,7 @@ class Simulator : public ISimulator
 
     Simulator(PHARE::initializer::PHAREDict const& dict,
               std::shared_ptr<PHARE::amr::Hierarchy> const& hierarchy);
+
     ~Simulator()
     {
         if (coutbuf != nullptr)
diff --git a/tests/amr/data/field/refine/test_refine_field.py b/tests/amr/data/field/refine/test_refine_field.py
index 308825c81..5716d6e3f 100644
--- a/tests/amr/data/field/refine/test_refine_field.py
+++ b/tests/amr/data/field/refine/test_refine_field.py
@@ -151,6 +151,10 @@ def refine(field, **kwargs):
             )
 
 
+    if fine_box.ndim == 3:
+        assert False
+
+
     # the fine data is ghost_nbr larger than expected to include ghost values from the coarser, so we drop the outer values
     if field.box.ndim == 1:
         fine_data = fine_data[field.ghosts_nbr[0] : -field.ghosts_nbr[0]]
diff --git a/tests/amr/data/particles/refine/input/input_1d_ratio_2.txt b/tests/amr/data/particles/refine/input/input_1d_ratio_2.txt
index 86e8b0c0c..496be4bf0 100644
--- a/tests/amr/data/particles/refine/input/input_1d_ratio_2.txt
+++ b/tests/amr/data/particles/refine/input/input_1d_ratio_2.txt
@@ -1,6 +1,6 @@
 CartesianGridGeometry
 {
-    domain_boxes = [ (0) , (64) ]
+    domain_boxes = [ (0) , (20) ]
     x_lo = 0.0
     x_up = 1.0
     periodic_dimension = 1
diff --git a/tests/amr/data/particles/refine/input/input_2d_ratio_2.txt b/tests/amr/data/particles/refine/input/input_2d_ratio_2.txt
index 1b0d9ef37..7d7d5d6bb 100644
--- a/tests/amr/data/particles/refine/input/input_2d_ratio_2.txt
+++ b/tests/amr/data/particles/refine/input/input_2d_ratio_2.txt
@@ -1,6 +1,6 @@
 CartesianGridGeometry
 {
-    domain_boxes = [ (0, 0) , (64, 64) ]
+    domain_boxes = [ (0, 0) , (20, 20) ]
     x_lo = 0.0, 0.0
     x_up = 1.0, 1.0
     periodic_dimension = 1, 1
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 687cd5670..f1cee33b5 100644
--- a/tests/amr/data/particles/refine/test_split.cpp
+++ b/tests/amr/data/particles/refine/test_split.cpp
@@ -20,7 +20,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 59d9c8a50..8c3021d57 100644
--- a/tests/core/numerics/interpolator/test_main.cpp
+++ b/tests/core/numerics/interpolator/test_main.cpp
@@ -836,6 +836,80 @@ 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;
+
+    using PHARE_TYPES               = PHARE::core::PHARE_Types<dim, interp_order>;
+    using NdArray_t                 = typename PHARE_TYPES::Array_t;
+    using ParticleArray_t           = typename PHARE_TYPES::ParticleArray_t;
+    using GridLayout_t              = typename PHARE_TYPES::GridLayout_t;
+    constexpr static auto safeLayer = static_cast<int>(1 + ghostWidthForParticles<interp_order>());
+
+    ACollectionOfParticles_3d()
+        : particles{grow(layout.AMRBox(), safeLayer)}
+        , rho{"field", HybridQuantity::Scalar::rho, nx, ny, nz}
+        , vx{"v_x", HybridQuantity::Scalar::Vx, nx, ny, nz}
+        , vy{"v_y", HybridQuantity::Scalar::Vy, nx, ny, nz}
+        , vz{"v_z", HybridQuantity::Scalar::Vz, nx, ny, nz}
+        , v{"v", HybridQuantity::Vector::V}
+    {
+        v.setBuffer("v_x", &vx);
+        v.setBuffer("v_y", &vy);
+        v.setBuffer("v_z", &vz);
+
+        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);
+    }
+
+    GridLayout<GridLayoutImplYee<dim, Interpolator::interp_order>> layout{
+        ConstArray<double, dim>(.1), {nx, ny}, ConstArray<double, dim>(0)};
+
+    ParticleArray_t particles;
+    Field<NdArray_t, typename HybridQuantity::Scalar> rho, vx, vy, vz;
+    VecField<NdArray_t, HybridQuantity> v;
+    Interpolator interpolator;
+};
+TYPED_TEST_SUITE_P(ACollectionOfParticles_3d);
+
+
+TYPED_TEST_P(ACollectionOfParticles_3d, DepositCorrectlyTheirWeight_3d)
+{
+    // EXPECT_DOUBLE_EQ(this->rho(7, 7, 7), 1.0);
+    // EXPECT_DOUBLE_EQ(this->vx(7, 7, 7), 2.0);
+    // EXPECT_DOUBLE_EQ(this->vy(7, 7, 7), -1.0);
+    // EXPECT_DOUBLE_EQ(this->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/core/numerics/pusher/test_pusher.cpp b/tests/core/numerics/pusher/test_pusher.cpp
index 154b598c3..a17073588 100644
--- a/tests/core/numerics/pusher/test_pusher.cpp
+++ b/tests/core/numerics/pusher/test_pusher.cpp
@@ -176,12 +176,9 @@ TEST_F(APusher3D, trajectoryIsOk)
 
     for (decltype(nt) i = 0; i < nt; ++i)
     {
-        actual[0][i]
-            = (particlesOut[0].iCell[0] + particlesOut[0].delta[0]) * static_cast<float>(dxyz[0]);
-        actual[1][i]
-            = (particlesOut[0].iCell[1] + particlesOut[0].delta[1]) * static_cast<float>(dxyz[1]);
-        actual[2][i]
-            = (particlesOut[0].iCell[2] + particlesOut[0].delta[2]) * static_cast<float>(dxyz[2]);
+        actual[0][i] = (particlesOut[0].iCell[0] + particlesOut[0].delta[0]) * dxyz[0];
+        actual[1][i] = (particlesOut[0].iCell[1] + particlesOut[0].delta[1]) * dxyz[1];
+        actual[2][i] = (particlesOut[0].iCell[2] + particlesOut[0].delta[2]) * dxyz[2];
 
         pusher->move(
             rangeIn, rangeOut, em, mass, interpolator, layout,
@@ -206,10 +203,8 @@ TEST_F(APusher2D, trajectoryIsOk)
 
     for (decltype(nt) i = 0; i < nt; ++i)
     {
-        actual[0][i]
-            = (particlesOut[0].iCell[0] + particlesOut[0].delta[0]) * static_cast<float>(dxyz[0]);
-        actual[1][i]
-            = (particlesOut[0].iCell[1] + particlesOut[0].delta[1]) * static_cast<float>(dxyz[1]);
+        actual[0][i] = (particlesOut[0].iCell[0] + particlesOut[0].delta[0]) * dxyz[0];
+        actual[1][i] = (particlesOut[0].iCell[1] + particlesOut[0].delta[1]) * dxyz[1];
 
         pusher->move(
             rangeIn, rangeOut, em, mass, interpolator, layout, [](auto& rge) { return rge; },
@@ -232,8 +227,7 @@ TEST_F(APusher1D, trajectoryIsOk)
 
     for (decltype(nt) i = 0; i < nt; ++i)
     {
-        actual[0][i]
-            = (particlesOut[0].iCell[0] + particlesOut[0].delta[0]) * static_cast<float>(dxyz[0]);
+        actual[0][i] = (particlesOut[0].iCell[0] + particlesOut[0].delta[0]) * dxyz[0];
 
         pusher->move(rangeIn, rangeOut, em, mass, interpolator, layout, selector, selector);
 
diff --git a/tests/diagnostic/CMakeLists.txt b/tests/diagnostic/CMakeLists.txt
index 3a5964a5f..b7b5378c1 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..5df839152
--- /dev/null
+++ b/tests/diagnostic/test-diagnostics_3d.cpp
@@ -0,0 +1,35 @@
+
+#include <mpi.h>
+
+#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/diagnostic/test_diagnostics.hpp b/tests/diagnostic/test_diagnostics.hpp
index b2b2cf50a..19ee34a01 100644
--- a/tests/diagnostic/test_diagnostics.hpp
+++ b/tests/diagnostic/test_diagnostics.hpp
@@ -45,7 +45,7 @@ void validateFluidGhosts(Data const& data, GridLayout const& layout, Field const
         std::size_t nans = 0;
         for (auto const& v : data)
             if (std::isnan(v))
-                nans++;
+                ++nans;
 
         auto phyStartIdx = layout.physicalStartIndex(field.physicalQuantity(), core::Direction::X);
         core::NdArrayMask mask{0, phyStartIdx - 2};
@@ -85,7 +85,6 @@ void validateFluidGhosts(Data const& data, GridLayout const& layout, Field const
     }
     else if constexpr (dim == 3)
     {
-        throw std::runtime_error("not implemented");
     }
 }
 
diff --git a/tests/simulator/__init__.py b/tests/simulator/__init__.py
index 3f0f5c928..6034831ca 100644
--- a/tests/simulator/__init__.py
+++ b/tests/simulator/__init__.py
@@ -12,6 +12,8 @@ def parse_cli_args(pop_from_sys = True):
         sys.argv = [sys.argv[0]]
     return r
 
+
+
 # Block accidental dictionary key rewrites
 class NoOverwriteDict(dict):
     def __init__(self, dict):
@@ -82,12 +84,13 @@ def density_2d_periodic(sim, x, y):
     xx, yy = meshify(x, y)
     return np.exp(-(xx-0.5*xmax)**2)*np.exp(-(yy-ymax/2.)**2) + background_particles
 
-# 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.)**2)*np.exp(-(zz-zmax/2.)**2) + background_particles
+    return r
 
 
 def defaultPopulationSettings(sim, density_fn, vbulk_fn):
@@ -250,3 +253,23 @@ def clean_up_diags_dirs(self):
             for diag_dir in self.diag_dirs:
                 if os.path.exists(diag_dir):
                     shutil.rmtree(diag_dir)
+
+
+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 eb51cca8c..d6f1c4367 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)
+        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_1d.py b/tests/simulator/advance/test_fields_advance_1d.py
index fc9579ea4..474fd8691 100644
--- a/tests/simulator/advance/test_fields_advance_1d.py
+++ b/tests/simulator/advance/test_fields_advance_1d.py
@@ -5,7 +5,7 @@
 
 import unittest
 from ddt import ddt, data, unpack
-from pyphare.core.box import Box, Box1D, nDBox
+from pyphare.core.box import Box1D
 from tests.simulator.test_advance import AdvanceTestBase
 
 import matplotlib
diff --git a/tests/simulator/advance/test_fields_advance_2d.py b/tests/simulator/advance/test_fields_advance_2d.py
index 7ccfbc387..4eae863ca 100644
--- a/tests/simulator/advance/test_fields_advance_2d.py
+++ b/tests/simulator/advance/test_fields_advance_2d.py
@@ -5,7 +5,7 @@
 
 import unittest
 from ddt import ddt, data, unpack
-from pyphare.core.box import Box, Box2D, nDBox
+from pyphare.core.box import Box2D
 from tests.simulator.test_advance import AdvanceTestBase
 
 import matplotlib
@@ -50,6 +50,7 @@ def test_overlaped_fields_are_equal_with_min_max_patch_size_of_max_ghosts(self,
         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}/{interp_order}/{self.ddt_test_id()}"
         largest_patch_size, smallest_patch_size = check_patch_size(ndim, interp_order=interp_order, cells=[60] * ndim)
+        print("largest_patch_size, smallest_patch_size", largest_patch_size, smallest_patch_size)
         datahier = self.getHierarchy(interp_order, refinement_boxes, "eb", diag_outputs=diag_outputs,
                                   smallest_patch_size=smallest_patch_size, largest_patch_size=smallest_patch_size,
                                   time_step=time_step, time_step_nbr=time_step_nbr, ndim=ndim, nbr_part_per_cell=ppc)
@@ -62,7 +63,7 @@ def test_overlaped_fields_are_equal_with_min_max_patch_size_of_max_ghosts(self,
        *per_interp(({"L0": {"B0": Box2D(6, 23)}})),
        *per_interp(({"L0": {"B0": Box2D( 2, 12), "B1": Box2D(13, 25)}})),
        *per_interp(({"L0": {"B0": Box2D( 5, 20)}, "L1": {"B0": Box2D(15, 19)}})),
-       *per_interp(({"L0": {"B0": Box2D( 5, 20)}, "L1": {"B0": Box2D(12, 38)}, "L2": {"B0": Box2D(30, 52)} })),
+       *per_interp(({"L0": {"B0": Box2D( 5, 20)}, "L1": {"B0": Box2D(12, 37)}, "L2": {"B0": Box2D(30, 51)} })),
     )
     @unpack
     def test_field_coarsening_via_subcycles(self, interp_order, refinement_boxes):
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..81f0c238e
--- /dev/null
+++ b/tests/simulator/advance/test_fields_advance_3d.py
@@ -0,0 +1,88 @@
+"""
+  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
+from ddt import ddt, data, unpack
+from pyphare.core.box import Box3D
+from tests.simulator.test_advance import AdvanceTestBase
+
+import matplotlib
+
+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
+        diag_outputs=f"phare_overlaped_fields_are_equal_{ndim}_{self.ddt_test_id()}"
+        datahier = self.getHierarchy(interp_order, refinement_boxes, "eb", diag_outputs=diag_outputs, cells=cells,
+                                  time_step=time_step, time_step_nbr=time_step_nbr, ndim=ndim, 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
+        diag_outputs=f"phare_overlaped_fields_are_equal_with_min_max_patch_size_of_max_ghosts/{ndim}/{interp_order}/{self.ddt_test_id()}"
+        largest_patch_size, smallest_patch_size = check_patch_size(ndim, interp_order=interp_order, cells=[cells] * ndim)
+        datahier = self.getHierarchy(interp_order, refinement_boxes, "eb", diag_outputs=diag_outputs, cells=cells,
+                                  smallest_patch_size=smallest_patch_size, largest_patch_size=smallest_patch_size,
+                                  time_step=time_step, time_step_nbr=time_step_nbr, ndim=ndim, 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)
+
+
+    # @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_1d.py b/tests/simulator/advance/test_particles_advance_1d.py
index d770d3aa6..c9ba702a0 100644
--- a/tests/simulator/advance/test_particles_advance_1d.py
+++ b/tests/simulator/advance/test_particles_advance_1d.py
@@ -5,7 +5,7 @@
 
 import unittest
 from ddt import ddt, data, unpack
-from pyphare.core.box import Box, Box1D, nDBox
+from pyphare.core.box import Box1D
 from tests.simulator.test_advance import AdvanceTestBase
 
 import matplotlib
diff --git a/tests/simulator/advance/test_particles_advance_2d.py b/tests/simulator/advance/test_particles_advance_2d.py
index 969ad568e..2a9bd66c9 100644
--- a/tests/simulator/advance/test_particles_advance_2d.py
+++ b/tests/simulator/advance/test_particles_advance_2d.py
@@ -5,7 +5,7 @@
 
 import unittest
 from ddt import ddt, data, unpack
-from pyphare.core.box import Box, Box2D, nDBox
+from pyphare.core.box import Box2D
 from tests.simulator.test_advance import AdvanceTestBase
 
 import matplotlib
@@ -33,6 +33,7 @@ def test_overlapped_particledatas_have_identical_particles(self, interp_order, r
         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)
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..53f7dc062
--- /dev/null
+++ b/tests/simulator/advance/test_particles_advance_3d.py
@@ -0,0 +1,44 @@
+"""
+  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
+from ddt import ddt, data, unpack
+from pyphare.core.box import Box3D
+from tests.simulator.test_advance import AdvanceTestBase
+
+import matplotlib
+
+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 51808cbdc..38ee5de50 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_1d.py b/tests/simulator/initialize/test_fields_init_1d.py
index e2f4f5311..2380e1702 100644
--- a/tests/simulator/initialize/test_fields_init_1d.py
+++ b/tests/simulator/initialize/test_fields_init_1d.py
@@ -5,7 +5,7 @@
 
 import unittest
 from ddt import ddt, data, unpack
-from pyphare.core.box import Box, Box1D, nDBox
+
 from tests.simulator.test_initialization import InitializationTest
 
 import matplotlib
diff --git a/tests/simulator/initialize/test_fields_init_2d.py b/tests/simulator/initialize/test_fields_init_2d.py
index a37c47599..a5b72fc9d 100644
--- a/tests/simulator/initialize/test_fields_init_2d.py
+++ b/tests/simulator/initialize/test_fields_init_2d.py
@@ -5,7 +5,7 @@
 
 import unittest
 from ddt import ddt, data, unpack
-from pyphare.core.box import Box, Box2D, nDBox
+
 from tests.simulator.test_initialization import InitializationTest
 
 import matplotlib
@@ -22,7 +22,7 @@ class InitializationTest(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..b4059c8c6
--- /dev/null
+++ b/tests/simulator/initialize/test_fields_init_3d.py
@@ -0,0 +1,44 @@
+"""
+  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 numpy as np
+import unittest
+from ddt import ddt, data, unpack
+
+from tests.simulator.test_initialization import InitializationTest
+
+import matplotlib
+
+matplotlib.use("Agg")  # for systems without GUI
+
+ndim = 3
+interp_orders = [1, 2, 3]
+ppc, cells = 10, 20
+
+@ddt
+class InitializationTest(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 622d9b789..ac8308090 100644
--- a/tests/simulator/initialize/test_particles_init_1d.py
+++ b/tests/simulator/initialize/test_particles_init_1d.py
@@ -5,7 +5,7 @@
 
 import unittest
 from ddt import ddt, data, unpack
-from pyphare.core.box import Box, Box1D, nDBox
+from pyphare.core.box import Box1D
 from tests.simulator.test_initialization import InitializationTest
 
 import matplotlib
diff --git a/tests/simulator/initialize/test_particles_init_2d.py b/tests/simulator/initialize/test_particles_init_2d.py
index 023de7eea..2fbd0f133 100644
--- a/tests/simulator/initialize/test_particles_init_2d.py
+++ b/tests/simulator/initialize/test_particles_init_2d.py
@@ -5,7 +5,7 @@
 
 import unittest
 from ddt import ddt, data, unpack
-from pyphare.core.box import Box, Box2D, nDBox
+from pyphare.core.box import Box2D
 from tests.simulator.test_initialization import InitializationTest
 
 import matplotlib
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..7f45f60a0
--- /dev/null
+++ b/tests/simulator/initialize/test_particles_init_3d.py
@@ -0,0 +1,94 @@
+"""
+  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
+from ddt import ddt, data, unpack
+from pyphare.core.box import Box, Box3D, nDBox
+from tests.simulator.test_initialization import InitializationTest
+
+import matplotlib
+
+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 InitializationTest(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(
+                interp_order,
+                refinement_boxes,
+                "particles",
+                ndim=ndim,
+                cells=cells, nbr_part_per_cell=ppc,
+                diag_outputs=f"phare_outputs/test_levelghost/{ndim}/{interp_order}/{self.ddt_test_id()}",
+            )
+        )
+        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 91d428d09..d06d4f43e 100644
--- a/tests/simulator/per_test.hpp
+++ b/tests/simulator/per_test.hpp
@@ -122,6 +122,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 ff1078a09..4f1bc047f 100644
--- a/tests/simulator/refined_particle_nbr.py
+++ b/tests/simulator/refined_particle_nbr.py
@@ -9,7 +9,6 @@
 import numpy as np
 import pyphare.pharein as ph
 from pyphare.simulator.simulator import Simulator
-from tests.simulator import NoOverwriteDict
 from tests.simulator import populate_simulation
 from pyphare.cpp import splitter_type
 
@@ -28,41 +27,32 @@ 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
 
         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()
@@ -76,54 +66,76 @@ 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)
-                )
+                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
                     self.assertTrue(max_per_pop < prev_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)
+
+    # 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 tearDown(self):
-        # needed in case exception is raised in test and Simulator
-        # not reset properly
-        if self.simulator is not None:
-            self.simulator.reset()
 
+    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 0657470d5..7dff6838d 100644
--- a/tests/simulator/test_advance.py
+++ b/tests/simulator/test_advance.py
@@ -192,13 +192,14 @@ def vthz(*xyz):
             return mom_hier
 
 
+
+
+
     def base_test_overlaped_fields_are_equal(self, datahier, coarsest_time):
         checks = 0
         for ilvl, overlaps in hierarchy_overlaps(datahier, coarsest_time).items():
             for overlap in overlaps:
-                pd1, pd2 = overlap["pdatas"]
-                box      = overlap["box"]
-                offsets  = overlap["offset"]
+                (pd1, pd2), box, offsets = overlap["pdatas"], overlap["box"], overlap["offset"]
 
                 self.assertEqual(pd1.quantity, pd2.quantity)
 
@@ -214,22 +215,11 @@ 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]))
-                    loc_b2 = 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]
+                    slice1 = boxm.select(pd1.dataset, boxm.amr_to_local(box, boxm.shift(pd1.ghost_box, offsets[0])))
+                    slice2 = boxm.select(pd2.dataset, boxm.amr_to_local(box, boxm.shift(pd2.ghost_box, offsets[1])))
+                    assert slice1.dtype == np.float64
 
                     try:
-                        assert slice1.dtype == np.float64
                         np.testing.assert_equal(slice1, slice2)
                         checks += 1
                     except AssertionError as e:
@@ -300,8 +290,8 @@ def _test_overlapped_particledatas_have_identical_particles(self, ndim, interp_o
                         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
 
@@ -324,7 +314,7 @@ def _test_L0_particle_number_conservation(self, ndim, interp_order, ppc=100):
 
 
 
-    def _test_field_coarsening_via_subcycles(self, dim, interp_order, refinement_boxes, **kwargs):
+    def _test_field_coarsening_via_subcycles(self, dim, interp_order, refinement_boxes, cells=60, **kwargs):
         print("test_field_coarsening_via_subcycles for dim/interp : {}/{}".format(dim, interp_order))
 
         from tests.amr.data.field.coarsening.test_coarsen_field import coarsen
@@ -333,7 +323,7 @@ def _test_field_coarsening_via_subcycles(self, dim, interp_order, refinement_box
         time_step_nbr=3
 
         diag_outputs=f"subcycle_coarsening/{dim}/{interp_order}/{self.ddt_test_id()}"
-        datahier = self.getHierarchy(interp_order, refinement_boxes, "eb", cells=60,
+        datahier = self.getHierarchy(interp_order, refinement_boxes, "eb", cells=cells,
                                       diag_outputs=diag_outputs, time_step=0.001,
                                       extra_diag_options={"fine_dump_lvl_max": 10},
                                       time_step_nbr=time_step_nbr,
@@ -396,12 +386,7 @@ def _test_field_coarsening_via_subcycles(self, dim, interp_order, refinement_box
                                     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, coarse_pd, fine_pd, coarseBox, fine_pdDataset, afterCoarse)
                                     np.testing.assert_allclose(coarse_pdDataset, afterCoarse, atol=1e-16, rtol=0)
diff --git a/tests/simulator/test_initialization.py b/tests/simulator/test_initialization.py
index 47322f44e..f113755c2 100644
--- a/tests/simulator/test_initialization.py
+++ b/tests/simulator/test_initialization.py
@@ -7,7 +7,7 @@
 from pyphare.pharein.diagnostics import ParticleDiagnostics, FluidDiagnostics, ElectromagDiagnostics
 from pyphare.pharein import ElectronModel
 from pyphare.pharein.simulation import Simulation
-from pyphare.pharesee.geometry import level_ghost_boxes, hierarchy_overlaps, touch_domain_border
+from pyphare.pharesee.geometry import level_ghost_boxes
 from pyphare.pharesee.particles import aggregate as aggregate_particles
 import pyphare.core.box as boxm
 from pyphare.core.box import nDBox
@@ -211,11 +211,14 @@ def vthz(*xyz):
 
 
 
-    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))
-        hier = self.getHierarchy(interp_order, refinement_boxes=None, qty="b", ndim=dim,
+        now = self.datetime_now()
+        hier = self.getHierarchy(interp_order, refinement_boxes=None, qty="b", ndim=dim, 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
         model = global_vars.sim.model
@@ -261,16 +264,28 @@ def _test_B_is_as_provided_by_user(self, dim, interp_order, **kwargs):
                     np.testing.assert_allclose(bz, bz_fn(xbz, ybz).reshape(bz.shape), atol=1e-16, rtol=0)
 
                 if dim == 3:
-                    raise ValueError("Unsupported dimension")
+                    zbx   = bx_pd.z[:]
+                    zby   = by_pd.z[:]
+                    zbz   = bz_pd.z[:]
 
+                    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):
+
+
+
+    def _test_bulkvel_is_as_provided_by_user(self, dim, interp_order, ppc=100, **kwargs):
         hier = self.getHierarchy(interp_order, {"L0": {"B0": nDBox(dim, 10, 19)}},
-                                 "moments", nbr_part_per_cell=100, beam=True, ndim=dim,  # ppc needs to be 10000?
-                                  diag_outputs=f"test_bulkV/{dim}/{interp_order}/{self.ddt_test_id()}")
+                                 "moments", nbr_part_per_cell=ppc, beam=True, ndim=dim,  # ppc needs to be 10000?
+                                  diag_outputs=f"test_bulkV/{dim}/{interp_order}/{self.ddt_test_id()}", **kwargs)
 
         from pyphare.pharein import global_vars
         model = global_vars.sim.model
@@ -286,92 +301,46 @@ 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):
-                    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)
+                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)])
 
-                    x = xx[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
-                    y = yy[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
+                def domain(patch_data):
+                    if dim == 1:
+                        return patch_data.dataset[select]
+                    return patch_data.dataset[:].reshape(patch.box.shape + (nbrGhosts * 2) + 1)[select]
 
-                    fpx = protons_Fx[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
-                    fpy = protons_Fy[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
-                    fpz = protons_Fz[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
+                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
 
-                    fbx = beam_Fx[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
-                    fby = beam_Fy[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
-                    fbz = beam_Fz[nbrGhosts:-nbrGhosts, nbrGhosts:-nbrGhosts]
+                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)
 
-                    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)
 
 
-
-
-
-    def _test_density_is_as_provided_by_user(self, dim, interp_order):
+    def _test_density_is_as_provided_by_user(self, dim, interp_order, **kwargs):
+        print(f"test_density_is_as_provided_by_user : dim {dim} interp_order {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))
-        hier = self.getHierarchy(interp_order, {"L0": {"B0": nDBox(dim, 10, 20)}},
+        nbParts = {1 : 10000, 2: 1000, 3: 20}
+        hier = self.getHierarchy(interp_order, {"L0": {"B0": nDBox(dim, 5, 14)}},
                                  qty="moments", nbr_part_per_cell=nbParts[dim], beam=True, ndim=dim,
-                                 diag_outputs=f"test_density/{dim}/{interp_order}/{self.ddt_test_id()}")
+                                 diag_outputs=f"test_density/{dim}/{interp_order}/{self.ddt_test_id()}", **kwargs)
 
         from pyphare.pharein import global_vars
         model = global_vars.sim.model
@@ -391,31 +360,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(dim)])
 
-                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
-
-                    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]
+                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
 
+                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)
@@ -430,11 +384,13 @@ def _test_density_is_as_provided_by_user(self, dim, interp_order):
 
 
 
-    def _test_density_decreases_as_1overSqrtN(self, dim, interp_order):
+    def _test_density_decreases_as_1overSqrtN(self, dim, interp_order, nbr_particles=None, cells=960):
         import matplotlib.pyplot as plt
         print(f"test_density_decreases_as_1overSqrtN, interp_order = {interp_order}")
 
-        nbr_particles = np.asarray([100, 1000, 5000, 10000])
+        if nbr_particles is None:
+            nbr_particles = np.asarray([100, 1000, 5000, 10000])
+
         noise = np.zeros(len(nbr_particles))
 
         for inbr,nbrpart in enumerate(nbr_particles):
@@ -443,9 +399,9 @@ def _test_density_decreases_as_1overSqrtN(self, dim, interp_order):
                                      nbr_part_per_cell=nbrpart,
                                      diag_outputs=f"1overSqrtN/{dim}/{interp_order}/{nbrpart}",
                                      density=lambda x:np.zeros_like(x)+1.,
-                                     smallest_patch_size=480,
-                                     largest_patch_size=480,
-                                     cells=960,
+                                     smallest_patch_size=int(cells/2),
+                                     largest_patch_size=int(cells/2),
+                                     cells=cells,
                                      dl=0.0125)
 
             from pyphare.pharein import global_vars
@@ -475,7 +431,6 @@ def _test_density_decreases_as_1overSqrtN(self, dim, interp_order):
             plt.close("all")
 
 
-
         plt.figure()
         plt.plot(nbr_particles, noise/noise[0], label=r"$\sigma/\sigma_0$")
         plt.plot(nbr_particles, 1/np.sqrt(nbr_particles/nbr_particles[0]), label=r"$1/sqrt(nppc/nppc0)$")
@@ -498,18 +453,18 @@ def _test_density_decreases_as_1overSqrtN(self, dim, interp_order):
 
 
 
-    def _test_nbr_particles_per_cell_is_as_provided(self, dim, interp_order, default_ppc=100):
+    def _test_nbr_particles_per_cell_is_as_provided(self, dim, interp_order, ppc=100, **kwargs):
         ddt_test_id = self.ddt_test_id()
-        datahier = self.getHierarchy(interp_order, {"L0": {"B0": nDBox(dim, 10, 20)}}, "particles", ndim=dim,
-                      diag_outputs=f"ppc/{dim}/{interp_order}/{ddt_test_id}")
+        datahier = self.getHierarchy(interp_order, {}, "particles", ndim=dim,
+                      diag_outputs=f"ppc/{dim}/{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, edges = np.histogramdd(icells)
-            self.assertTrue((H == default_ppc).all())
-
-
+            H, edges = np.histogramdd(icells, bins=patch.box.shape)
+            self.assertTrue((H == ppc).all())
 
 
 
@@ -536,6 +491,8 @@ def _test_domainparticles_have_correct_split_from_coarser_particle(self, ndim, i
         datahier = self.getHierarchy(interp_order, refinement_boxes, "particles", ndim=ndim,
            diag_outputs=f"coarser_split/{ndim}/{interp_order}/{ddt_test_id}", cells=30, **kwargs)
 
+        if cpp.mpi_rank() > 0: return
+
         from pyphare.pharein.global_vars import sim
         assert sim is not None and len(sim.cells) == ndim
 
diff --git a/tests/simulator/test_python_concurrent.py b/tests/simulator/test_python_concurrent.py
index 89e0d9fa4..c12589955 100644
--- a/tests/simulator/test_python_concurrent.py
+++ b/tests/simulator/test_python_concurrent.py
@@ -4,54 +4,137 @@
 """
 
 import os
+import time
 import unittest
-import multiprocessing
+from multiprocessing import Process, Queue, cpu_count
 
 from tests.simulator.test_validation import SimulatorValidation
-
+from tests.simulator.test_diagnostics import DiagnosticsTest  # mpirun -n 1/2/3/4
 from tests.simulator.initialize.test_fields_init_1d import InitializationTest as InitField1d
 from tests.simulator.initialize.test_particles_init_1d import InitializationTest as InitParticles1d
-
 from tests.simulator.advance.test_fields_advance_1d import AdvanceTest as AdvanceField1d
 from tests.simulator.advance.test_particles_advance_1d import AdvanceTest as AdvanceParticles1d
-
 from tests.simulator.initialize.test_fields_init_2d import InitializationTest as InitField2d
 from tests.simulator.initialize.test_particles_init_2d import InitializationTest as InitParticles2d
-
 from tests.simulator.advance.test_fields_advance_2d import AdvanceTest as AdvanceField2d
 from tests.simulator.advance.test_particles_advance_2d import AdvanceTest as AdvanceParticles2d
+from tests.simulator.initialize.test_fields_init_3d import InitializationTest as InitField3d
+from tests.simulator.initialize.test_particles_init_3d import InitializationTest as InitParticles3d
+from tests.simulator.advance.test_fields_advance_3d import AdvanceTest as AdvanceField3d
+from tests.simulator.advance.test_particles_advance_3d import AdvanceTest as AdvanceParticles3d
 
+from tools.python3 import run, find_on_path
 
-N_CORES = int(os.environ["N_CORES"]) if "N_CORES" in os.environ else multiprocessing.cpu_count()
-MPI_RUN = int(os.environ["MPI_RUN"]) if "MPI_RUN" in os.environ else 1
+N_CORES = int(os.environ["N_CORES"]) if "N_CORES" in os.environ else cpu_count()
+MPI_RUN = os.environ["MPI_RUN"] if "MPI_RUN" in os.environ else 1
 PRINT   = int(os.environ["PRINT"]) if "PRINT" in os.environ else 0
 
-def test_cmd(clazz, test_id):
-    return f"mpirun -n {MPI_RUN} python3 -m {clazz.__module__} {clazz.__name__}.{test_id}"
+MPI_RUN_EXTRA = os.environ["MPI_RUN_EXTRA"] if "MPI_RUN_EXTRA" in os.environ else ""
 
-if __name__ == "__main__":
+def test_cmd(clazz, test_id, mpi_run):
+    return f"mpirun {MPI_RUN_EXTRA} -n {mpi_run} python3 -m {clazz.__module__} {clazz.__name__}.{test_id}"
+
+test_classes_to_run = [
+  SimulatorValidation, DiagnosticsTest,
+  InitField1d,         InitParticles1d,
+  AdvanceField1d,      AdvanceParticles1d,
+  InitField2d,         InitParticles2d,
+  AdvanceField2d,      AdvanceParticles2d,
+  InitField3d,         InitParticles3d,
+  AdvanceField3d,      AdvanceParticles3d
+]
+
+class TestBatch:
+    def __init__(self, tests, mpi_run = 1):
+        self.tests = tests
+        self.mpi_run = mpi_run
 
-    test_classes_to_run = [
-      SimulatorValidation,
-      InitField1d,
-      InitParticles1d,
-      AdvanceField1d,
-      AdvanceParticles1d,
-      InitField2d,
-      InitParticles2d,
-      AdvanceField2d,
-      AdvanceParticles2d
-    ]
-
-    tests = []
-    loader = unittest.TestLoader()
-    for test_class in test_classes_to_run:
+def load_test_cases_in(classes, mpi_run = 1):
+    tests, loader = [], unittest.TestLoader()
+    for test_class in classes:
         for suite in loader.loadTestsFromTestCase(test_class):
-            tests += [test_cmd(type(suite), suite._testMethodName)]
+            tests += [test_cmd(type(suite), suite._testMethodName, mpi_run)]
+    return TestBatch(tests, mpi_run)
 
-    from tools.python3 import run_mp
-    if PRINT:
-        for test in tests:
-            print(test)
+def build_batches():
+    batches = []
+    if MPI_RUN=="cmake":
+        batches += [load_test_cases_in(test_classes_to_run, 1)]
+        batches += [load_test_cases_in(test_classes_to_run, 2)]
+        batches += [load_test_cases_in([DiagnosticsTest], 3)]
+        batches += [load_test_cases_in([DiagnosticsTest], 4)]
     else:
-        run_mp(tests, N_CORES, check=True)
+        batches += [load_test_cases_in(test_classes_to_run, int(MPI_RUN))]
+    return batches
+
+class CallableTest:
+    def __init__(self, batch_index, cmd):
+        self.batch_index = batch_index
+        self.cmd = cmd
+        self.run = None
+
+    def __call__(self, **kwargs):
+        self.run = run(self.cmd.split(), shell=False, capture_output=True, check=True, print_cmd=False)
+        return self
+
+class CoreCount:
+  def __init__(self, cores_avail):
+      self.cores_avail = cores_avail
+      self.proces = []
+      self.fin = []
+
+def runner(runnable, queue):
+    queue.put(runnable())
+
+def print_tests(batches):
+    for batch in batches:
+        for test in batch.tests:
+            print(test)
+
+def main():
+    batches = build_batches()
+
+    if PRINT:
+        print_tests(batches)
+        return
+
+    cc = CoreCount(N_CORES)
+    assert cc.cores_avail >= max([batch.mpi_run for batch in batches])
+    cc.procs = [[] for batch in batches]
+    cc.fin = [0 for batch in batches]
+    pqueue = Queue()
+
+    def launch_tests():
+        for batch_index, batch in enumerate(batches):
+            offset = len(cc.procs[batch_index])
+            for test_index, test in enumerate(batch.tests[offset:]):
+                test_index += offset
+                if batch.mpi_run <= cc.cores_avail:
+                    test = CallableTest(batch_index, batches[batch_index].tests[test_index])
+                    cc.cores_avail -= batch.mpi_run
+                    cc.procs[batch_index] += [Process(target=runner, args=(test, (pqueue),))]
+                    cc.procs[batch_index][-1].daemon = True
+                    cc.procs[batch_index][-1].start()
+
+    def finished():
+        b = True
+        for batch_index, batch in enumerate(batches):
+            b &= cc.fin[batch_index] == len(batch.tests)
+        return b
+
+    def waiter(queue):
+        while True:
+            proc = queue.get()
+            time.sleep(.01) # don't throttle!
+            if (isinstance(proc, CallableTest)):
+                print(proc.cmd, f"finished in {proc.run.t:.2f} seconds")
+                cc.cores_avail += batches[proc.batch_index].mpi_run
+                cc.fin[proc.batch_index] += 1
+                launch_tests()
+                if finished(): break
+
+    launch_tests()
+    waiter(pqueue)
+
+if __name__ == "__main__":
+    main()
diff --git a/tools/ctest_exec_mp.py b/tools/ctest_exec_mp.py
deleted file mode 100644
index c14ae09ec..000000000
--- a/tools/ctest_exec_mp.py
+++ /dev/null
@@ -1,47 +0,0 @@
-
-
-import os
-import sys
-from pathlib import Path
-
-not_root_error = "script must be run from project root directory"
-assert all([os.path.exists(d) for d in ["tests", "tools", "CMakeLists.txt"]]), not_root_error
-root = Path(os.getcwd())
-sys.path.insert(0, ".")
-
-from tools.python3 import run, pushd
-import tools.python3.cmake as cmake
-import tools.python3.git as git
-
-
-def run_tests_log_to_file(data_dir, n_cores, tests):
-    """ ctest seems to merge stdout and stderr so we only need one output file """
-    import concurrent.futures
-
-    with concurrent.futures.ThreadPoolExecutor(max_workers=n_cores) as executor:
-        jobs = [
-          executor.submit(
-            run, cmake.test_cmd(test, verbose=True),
-            shell=True, capture_output=False, check=True,
-            stdout=open(os.path.join(str(data_dir), test), "w"))
-            for i, test in enumerate(tests)
-        ]
-        results = []
-        for future in concurrent.futures.as_completed(jobs):
-            try:
-                results += [future.result()]
-            except Exception as exc:
-                print("run_mp generated an exception: %s" % exc)
-        return results
-
-def main():
-    with pushd(os.path.join(str(root), "build")):
-        current_git_hash = git.hashes(1)[0]
-        data_dir = Path(os.path.join(str(root), "data_out", current_git_hash, "ctest_logs"))
-        os.makedirs(str(data_dir), exist_ok=True)
-        N_CORES= int(os.environ["N_CORES"]) if "N_CORES" in os.environ else 1
-        print(f"Launching ctests with N_CORES {N_CORES}")
-        run_tests_log_to_file(data_dir, N_CORES, cmake.list_tests())
-
-if __name__ == "__main__":
-    main()
diff --git a/tools/python3/__init__.py b/tools/python3/__init__.py
index 16db65902..6e8ecf55d 100644
--- a/tools/python3/__init__.py
+++ b/tools/python3/__init__.py
@@ -1,6 +1,18 @@
 def decode_bytes(input):
     return input.decode("ascii", errors="ignore")
 
+class RunTimer:
+    def __init__(self, cmd, shell=True, capture_output=True, check=False, print_cmd=True, **kwargs):
+        import time
+        import subprocess
+        self.cmd = cmd
+        start = time.time()
+        self.run = subprocess.run(self.cmd, shell=shell, capture_output=capture_output, check=check, **kwargs)
+        self.t = time.time() - start
+        self.stdout = self.run.stdout
+        self.stderr = self.run.stderr
+
+
 
 def run(cmd, shell=True, capture_output=True, check=False, print_cmd=True, **kwargs):
     """
@@ -11,10 +23,13 @@ def run(cmd, shell=True, capture_output=True, check=False, print_cmd=True, **kwa
     if print_cmd:
         print(f"running: {cmd}")
     try:
-        return subprocess.run(cmd, shell=shell, capture_output=capture_output, check=check, **kwargs)
+        return RunTimer(cmd, shell=shell, capture_output=capture_output, check=check, **kwargs)
     except subprocess.CalledProcessError as e: # only triggers on failure if check=True
-        print(f"run failed with error: {e}\n\t{e.stdout}\n\t{e.stderr} ")
-        raise RuntimeError(decode_bytes(e.stderr))
+        what = f"run failed with error: {e}"
+        print(what)
+        if capture_output:
+            raise RuntimeError(decode_bytes(e.stderr))
+        raise RuntimeError(what)
 
 def run_mp(cmds, N_CORES=None, **kwargs):
     """
@@ -31,9 +46,11 @@ def run_mp(cmds, N_CORES=None, **kwargs):
         results = []
         for future in concurrent.futures.as_completed(jobs):
             try:
-                results += [future.result()]
+                proc = future.result()
+                results += [proc]
                 if future.exception() is not None:
                     raise future.exception()
+                print(proc.cmd, f"finished in {proc.t:.2f} seconds")
             except Exception as exc:
                 if kwargs.get("check", False):
                     executor.shutdown(wait=False, cancel_futures=True)
@@ -43,11 +60,20 @@ def run_mp(cmds, N_CORES=None, **kwargs):
         return results
 
 
+def find_on_path(file):
+    import os
+    for dir in os.environ["PATH"].split(os.pathsep):
+        full = os.path.join(dir, file)
+        if os.path.exists(full):
+            return full
+    return ""
+
+
 def binary_exists_on_path(bin):
     """
     https://linux.die.net/man/1/which
     """
-    return run(f"which {bin}").returncode == 0
+    return len(find_on_path(bin))
 
 
 def scan_dir(path, files_only=False, dirs_only=False, drop=[]):
@@ -75,3 +101,4 @@ def pushd(new_cwd):
         yield
     finally:
         os.chdir(cwd)
+