diff --git a/.github/workflows/cmake_ubuntu.yml b/.github/workflows/cmake_ubuntu.yml
index 0cc64efde..b5919ce7d 100644
--- a/.github/workflows/cmake_ubuntu.yml
+++ b/.github/workflows/cmake_ubuntu.yml
@@ -57,7 +57,7 @@ jobs:
               -DCMAKE_C_COMPILER_LAUNCHER=ccache \
               -DCMAKE_CXX_COMPILER_LAUNCHER=ccache \
               -DlowResourceTests=ON -DdevMode=ON -Dbench=ON \
-              -DCMAKE_CXX_FLAGS="-O3 -DPHARE_DIAG_DOUBLES=1 " -Dphare_configurator=ON
+              -DCMAKE_CXX_FLAGS="-O3 -DPHARE_DIAG_DOUBLES=1" -Dphare_configurator=ON
 
     - name: Build
       working-directory: ${{github.workspace}}/build
diff --git a/pyphare/pyphare/pharein/__init__.py b/pyphare/pyphare/pharein/__init__.py
index 5a2514275..65955f1aa 100644
--- a/pyphare/pyphare/pharein/__init__.py
+++ b/pyphare/pyphare/pharein/__init__.py
@@ -328,6 +328,9 @@ def as_paths(rb):
             )
 
     if len(simulation.diagnostics) > 0:
+        if simulation.diag_options is not None and "format" in simulation.diag_options:
+            add_string(diag_path + "format", simulation.diag_options["format"])
+
         if simulation.diag_options is not None and "options" in simulation.diag_options:
             add_string(
                 diag_path + "filePath", simulation.diag_options["options"]["dir"]
diff --git a/pyphare/pyphare/pharein/diagnostics.py b/pyphare/pyphare/pharein/diagnostics.py
index e11801513..1a63f8871 100644
--- a/pyphare/pyphare/pharein/diagnostics.py
+++ b/pyphare/pyphare/pharein/diagnostics.py
@@ -62,9 +62,9 @@ def validate_timestamps(clazz, key, **kwargs):
     timestamps = phare_utilities.np_array_ify(kwargs.get(key, []))
 
     if np.any(timestamps < init_time):
-        raise RuntimeError(
-            f"Error: timestamp({sim.time_step_nbr}) cannot be less than simulation.init_time({init_time}))"
-        )
+        timestamps = timestamps[timestamps >= init_time]
+        print(f"Warning: some timestamps below ({init_time}) are filtered")
+
     if np.any(timestamps > sim.final_time):
         raise RuntimeError(
             f"Error: timestamp({sim.time_step_nbr}) cannot be greater than simulation.final_time({sim.final_time}))"
diff --git a/pyphare/pyphare/pharein/simulation.py b/pyphare/pyphare/pharein/simulation.py
index 684c07084..b0db542b3 100644
--- a/pyphare/pyphare/pharein/simulation.py
+++ b/pyphare/pyphare/pharein/simulation.py
@@ -160,31 +160,39 @@ def check_time(**kwargs):
         and "time_step" not in kwargs
     )
 
+    if sum([final_and_dt, final_and_nsteps, nsteps_and_dt]) != 1:
+        raise ValueError(
+            "Error: Specify either 'final_time' and 'time_step' or 'time_step_nbr' and 'time_step'"
+            + " or 'final_time' and 'time_step_nbr'"
+        )
+
     start_time = kwargs.get("restart_options", {}).get("restart_time", 0)
 
+    def _final_time():
+        if "final_time" in kwargs:
+            return kwargs["final_time"]
+        return start_time + kwargs["time_step"] * kwargs["time_step_nbr"]
+
+    final_time = _final_time()
+    total_time = final_time - start_time
+    if total_time < 0:
+        raise RuntimeError("Simulation time cannot be negative - review inputs")
+    if phare_utilities.fp_equal(total_time, 0):
+        return 0, 0, final_time
+
     if final_and_dt:
-        time_step_nbr = int(kwargs["final_time"] / kwargs["time_step"])
-        time_step = kwargs["final_time"] / time_step_nbr
+        time_step_nbr = int(total_time / kwargs["time_step"])
+        time_step = total_time / time_step_nbr
 
     elif final_and_nsteps:
-        time_step = kwargs["final_time"] / kwargs["time_step_nbr"]
+        time_step = total_time / kwargs["time_step_nbr"]
         time_step_nbr = kwargs["time_step_nbr"]
 
-    elif nsteps_and_dt:
+    else:  # must be nsteps_and_dt
         time_step = kwargs["time_step"]
         time_step_nbr = kwargs["time_step_nbr"]
 
-    else:
-        raise ValueError(
-            "Error: Specify either 'final_time' and 'time_step' or 'time_step_nbr' and 'time_step'"
-            + " or 'final_time' and 'time_step_nbr'"
-        )
-
-    return (
-        time_step_nbr,
-        time_step,
-        kwargs.get("final_time", start_time + time_step * time_step_nbr),
-    )
+    return time_step_nbr, time_step, final_time
 
 
 # ------------------------------------------------------------------------------
@@ -485,7 +493,7 @@ def check_directory(directory, key):
 # diag_options = {"format":"phareh5", "options": {"dir": "phare_ouputs/"}}
 def check_diag_options(**kwargs):
     diag_options = kwargs.get("diag_options", None)
-    formats = ["phareh5"]
+    formats = ["phareh5", "pharevtkhdf"]
     if diag_options is not None and "format" in diag_options:
         if diag_options["format"] not in formats:
             raise ValueError("Error - diag_options format is invalid")
@@ -731,8 +739,10 @@ def wrapper(simulation_object, **kwargs):
                 kwargs["max_nbr_levels"],
             ) = check_refinement_boxes(ndim, **kwargs)
         else:
-            kwargs["max_nbr_levels"] = kwargs.get("max_nbr_levels", None)
-            assert kwargs["max_nbr_levels"] is not None  # this needs setting otherwise
+            if "max_nbr_levels" not in kwargs:
+                print("WARNING, 'max_nbr_levels' is not set, defaulting to 1")
+            kwargs["max_nbr_levels"] = kwargs.get("max_nbr_levels", 1)
+
             kwargs["refinement_boxes"] = None
             kwargs["tagging_threshold"] = kwargs.get("tagging_threshold", 0.1)
 
diff --git a/pyphare/pyphare/pharesee/phare_vtk/__init__.py b/pyphare/pyphare/pharesee/phare_vtk/__init__.py
new file mode 100644
index 000000000..9af411dae
--- /dev/null
+++ b/pyphare/pyphare/pharesee/phare_vtk/__init__.py
@@ -0,0 +1,7 @@
+#
+#
+#
+
+from .plot import plot
+
+__all__ = ["plot"]
diff --git a/pyphare/pyphare/pharesee/phare_vtk/base.py b/pyphare/pyphare/pharesee/phare_vtk/base.py
new file mode 100644
index 000000000..fcd62e771
--- /dev/null
+++ b/pyphare/pyphare/pharesee/phare_vtk/base.py
@@ -0,0 +1,110 @@
+#
+#
+#
+
+import vtk
+
+
+class PhaseOutput:
+    def __init__(self, **kwargs):
+        self.kwargs = {**kwargs}
+
+    def __iter__(self):
+        return self.kwargs.items().__iter__()
+
+
+def surface_filter(output, **kwargs):
+    surface = vtk.vtkDataSetSurfaceFilter()
+    surface.SetInputConnection(output.GetOutputPort())
+    return PhaseOutput(surface=surface)
+
+
+def composite_data_geometry_filter(output, **kwargs):
+    geom = vtk.vtkCompositeDataGeometryFilter()
+    geom.SetInputConnection(output.GetOutputPort())
+    return PhaseOutput(geom=geom)
+
+
+def poly_data_mapper(output, **kwargs):
+    array_name = kwargs.get("array_name", "data")
+    mapper = vtk.vtkPolyDataMapper()
+    mapper.SetInputConnection(output.GetOutputPort())
+    mapper.SelectColorArray(array_name)
+    mapper.SetScalarModeToUsePointData()
+    mapper.SetScalarRange(
+        output.GetOutput().GetPointData().GetArray(array_name).GetRange(0)
+    )
+    return PhaseOutput(mapper=mapper)
+
+
+def all_times_in(reader):
+    info = reader.GetOutputInformation(0)
+
+    time_key = vtk.vtkStreamingDemandDrivenPipeline.TIME_STEPS()
+
+    if not info.Has(time_key):
+        raise RuntimeError("VTK Error: cannot get times from file")
+
+    times = info.Get(time_key)
+    return times
+
+
+def bounds_in(reader):
+    amr = reader.GetOutput()
+    bbox = vtk.vtkBoundingBox()
+
+    for level in range(amr.GetNumberOfLevels()):
+        for idx in range(amr.GetNumberOfDataSets(level)):
+            ds = amr.GetDataSet(level, idx)
+            if ds is not None:
+                b = [0.0] * 6
+                ds.GetBounds(b)
+                bbox.AddBounds(b)
+
+    global_bounds = [0.0] * 6
+    bbox.GetBounds(global_bounds)
+    return global_bounds
+
+
+def _default_phases():
+    # Typical use case for PHARE fields
+    return [surface_filter, composite_data_geometry_filter, poly_data_mapper]
+
+
+class VtkFile:
+    def __init__(self, filename, time=None, array_name="data", phases=None):
+        phases = phases if phases is not None else _default_phases()
+        if len(phases) == 0:
+            raise RuntimeError("Error: Zero phases!")
+
+        self.array_name = array_name
+        self.reader = vtk.vtkHDFReader()
+        self.reader.SetFileName(filename)
+        self.reader.Update()
+        self.reader.UpdateInformation()
+
+        self.times = all_times_in(self.reader)
+        self.reader.UpdateTimeStep(time if time is not None else self.times[-1])
+
+        _in = self.reader
+        for i in range(len(phases)):
+            ret = phases[i](_in, array_name=array_name)
+            for key, val in ret:
+                if hasattr(val, "Update"):
+                    val.Update()
+                setattr(self, key, val)
+                _in = val
+
+        self.mapper = _in
+        self.bounds = bounds_in(self.reader)
+        self.spacing = self.reader.GetOutput().GetDataSet(0, 0).GetSpacing()
+
+
+class VtkFieldFile(VtkFile):
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+
+
+class VtkTensorFieldFile(VtkFile):
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
diff --git a/pyphare/pyphare/pharesee/phare_vtk/plot.py b/pyphare/pyphare/pharesee/phare_vtk/plot.py
new file mode 100644
index 000000000..b5a81aec8
--- /dev/null
+++ b/pyphare/pyphare/pharesee/phare_vtk/plot.py
@@ -0,0 +1,38 @@
+#
+#
+#
+
+
+import vtk
+from .base import VtkTensorFieldFile
+
+
+def plot(vtk_file, out_file="vtk.png"):
+    if isinstance(vtk_file, str):
+        vtk_file = VtkTensorFieldFile(vtk_file)
+    vtk_file.geom.GetOutput().GetPointData().SetActiveScalars("data")
+
+    colors = vtk.vtkNamedColors()
+
+    renderer = vtk.vtkRenderer()
+    renderer.AddActor(vtk.vtkActor(mapper=vtk_file.mapper))
+    renderer.SetBackground(colors.GetColor3d("Silver"))
+
+    ren_win = vtk.vtkRenderWindow()
+    ren_win.AddRenderer(renderer)
+    ren_win.SetSize(800, 600)
+    ren_win.OffScreenRenderingOn()  # do not flash image on screen
+
+    iren = vtk.vtkRenderWindowInteractor()
+    iren.SetRenderWindow(ren_win)
+
+    ren_win.Render()
+
+    w2if = vtk.vtkWindowToImageFilter()
+    w2if.SetInput(ren_win)
+    w2if.Update()
+
+    writer = vtk.vtkPNGWriter()
+    writer.SetFileName(out_file)
+    writer.SetInputConnection(w2if.GetOutputPort())
+    writer.Write()
diff --git a/pyphare/pyphare/pharesee/plotting.py b/pyphare/pyphare/pharesee/plotting.py
index 8f7614919..cef70607f 100644
--- a/pyphare/pyphare/pharesee/plotting.py
+++ b/pyphare/pyphare/pharesee/plotting.py
@@ -325,7 +325,7 @@ def finest_field_plot(run_path, qty, **kwargs):
         ax.set_aspect("equal")
         divider = make_axes_locatable(ax)
         cax = divider.append_axes("right", size="5%", pad=0.08)
-        cb = plt.colorbar(im, ax=ax, cax=cax)
+        plt.colorbar(im, ax=ax, cax=cax)
     else:
         raise ValueError("finest_field_plot not yet ready for 3d")
 
diff --git a/pyphare/pyphare_tests/test_pharesee/test_geometry.py b/pyphare/pyphare_tests/test_pharesee/test_geometry.py
index 6023f429f..c9ae8d1f5 100644
--- a/pyphare/pyphare_tests/test_pharesee/test_geometry.py
+++ b/pyphare/pyphare_tests/test_pharesee/test_geometry.py
@@ -116,7 +116,6 @@ def test_particle_ghost_area_boxes(self):
         }
 
         gaboxes = ghost_area_boxes(hierarchy, "particles")
-        particles = "particles"
 
         # same number of levels
         self.assertEqual(len(expected), len(gaboxes))
diff --git a/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py b/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py
index 908484607..edd8a7fbd 100644
--- a/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py
+++ b/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py
@@ -2,12 +2,12 @@
 from ddt import ddt
 import numpy as np
 
-from pyphare.simulator.simulator import Simulator
+from pyphare.core.box import Box
 from pyphare.pharesee.run import Run
+from pyphare.simulator.simulator import Simulator
 from pyphare.pharesee.hierarchy import PatchHierarchy
 from pyphare.pharesee.hierarchy import ScalarField, VectorField
 from pyphare.core.operators import dot, cross, sqrt, modulus, grad
-from pyphare.core.box import Box
 
 diag_outputs = "phare_outputs/"
 time_step_nbr = 20
@@ -152,12 +152,12 @@ def vthz(x, y):
 
         Simulator(config()).run()
 
-    def test_data_is_a_hierarchy(self):
+    def _test_data_is_a_hierarchy(self):
         r = Run(self.diag_dir())
         B = r.GetB(0.0)
         self.assertTrue(isinstance(B, PatchHierarchy))
 
-    def test_can_read_multiple_times(self):
+    def _test_can_read_multiple_times(self):
         r = Run(self.diag_dir())
         times = (0.0, 0.1)
         B = r.GetB(times)
@@ -168,7 +168,7 @@ def test_can_read_multiple_times(self):
             self.assertEqual(len(hier.times()), 2)
             self.assertTrue(np.allclose(hier.times().astype(np.float32), times))
 
-    def test_hierarchy_is_refined(self):
+    def _test_hierarchy_is_refined(self):
         r = Run(self.diag_dir())
         time = 0.0
         B = r.GetB(time)
@@ -177,19 +177,19 @@ def test_hierarchy_is_refined(self):
         self.assertEqual(len(B.levels(time)), 2)
         self.assertEqual(len(B.levels(time)), B.levelNbr(time))
 
-    def test_can_get_nbytes(self):
+    def _test_can_get_nbytes(self):
         r = Run(self.diag_dir())
         time = 0.0
         B = r.GetB(time)
         self.assertGreater(B.nbytes(), 0)
 
-    def test_hierarchy_has_patches(self):
+    def _test_hierarchy_has_patches(self):
         r = Run(self.diag_dir())
         time = 0.0
         B = r.GetB(time)
         self.assertGreater(B.nbrPatches(), 0)
 
-    def test_access_patchdatas_as_hierarchies(self):
+    def _test_access_patchdatas_as_hierarchies(self):
         r = Run(self.diag_dir())
         time = 0.0
         B = r.GetB(time)
@@ -197,7 +197,7 @@ def test_access_patchdatas_as_hierarchies(self):
         self.assertTrue(isinstance(B.y, PatchHierarchy))
         self.assertTrue(isinstance(B.z, PatchHierarchy))
 
-    def test_partial_domain_hierarchies(self):
+    def _test_partial_domain_hierarchies(self):
         import matplotlib.pyplot as plt
         from matplotlib.patches import Rectangle
 
@@ -226,7 +226,7 @@ def test_partial_domain_hierarchies(self):
         self.assertTrue(isinstance(Bpartial, PatchHierarchy))
         self.assertLess(Bpartial.nbrPatches(), B.nbrPatches())
 
-    def test_scalarfield_quantities(self):
+    def _test_scalarfield_quantities(self):
         r = Run(self.diag_dir())
         time = 0.0
         Ni = r.GetNi(time)
@@ -234,7 +234,7 @@ def test_scalarfield_quantities(self):
         self.assertTrue(isinstance(Ni, ScalarField))
         self.assertTrue(isinstance(Vi, VectorField))
 
-    def test_sum_two_scalarfields(self):
+    def _test_sum_two_scalarfields(self):
         r = Run(self.diag_dir())
         time = 0.0
         Ni = r.GetNi(time)
@@ -243,7 +243,7 @@ def test_sum_two_scalarfields(self):
         self.assertTrue(isinstance(s, ScalarField))
         self.assertEqual(s.quantities(), ["value"])
 
-    def test_sum_with_scalar(self):
+    def _test_sum_with_scalar(self):
         r = Run(self.diag_dir())
         time = 0.0
         Ni = r.GetNi(time)
@@ -253,7 +253,7 @@ def test_sum_with_scalar(self):
             self.assertTrue(isinstance(s, ScalarField))
             self.assertEqual(s.quantities(), ["value"])
 
-    def test_scalarfield_difference(self):
+    def _test_scalarfield_difference(self):
         r = Run(self.diag_dir())
         time = 0.0
         Ni = r.GetNi(time)
@@ -265,7 +265,7 @@ def test_scalarfield_difference(self):
             self.assertTrue(isinstance(s, ScalarField))
             self.assertEqual(s.quantities(), ["value"])
 
-    def test_scalarfield_product(self):
+    def _test_scalarfield_product(self):
         r = Run(self.diag_dir())
         time = 0.0
         Ni = r.GetNi(time)
@@ -276,7 +276,7 @@ def test_scalarfield_product(self):
             self.assertTrue(isinstance(s, ScalarField))
             self.assertEqual(s.quantities(), ["value"])
 
-    def test_scalarfield_division(self):
+    def _test_scalarfield_division(self):
         r = Run(self.diag_dir())
         time = 0.0
         Ni = r.GetNi(time)
@@ -289,14 +289,14 @@ def test_scalarfield_division(self):
             self.assertTrue(isinstance(s, ScalarField))
             self.assertEqual(s.quantities(), ["value"])
 
-    def test_scalarfield_sqrt(self):
+    def _test_scalarfield_sqrt(self):
         r = Run(self.diag_dir())
         time = 0.0
         Ni = r.GetNi(time)
         nisqrt = sqrt(Ni)
         self.assertTrue(isinstance(nisqrt, ScalarField))
 
-    def test_vectorfield_dot_product(self):
+    def _test_vectorfield_dot_product(self):
         r = Run(self.diag_dir())
         time = 0.0
         B = r.GetB(time)
@@ -306,7 +306,7 @@ def test_vectorfield_dot_product(self):
             self.assertTrue(isinstance(s, ScalarField))
             self.assertEqual(s.quantities(), ["value"])
 
-    def test_vectorfield_binary_ops(self):
+    def _test_vectorfield_binary_ops(self):
         r = Run(self.diag_dir())
         time = 0.0
         B = r.GetB(time)
@@ -321,7 +321,7 @@ def test_vectorfield_binary_ops(self):
             self.assertTrue(isinstance(s, VectorField))
             self.assertEqual(s.quantities(), ["x", "y", "z"])
 
-    def test_scalarfield_to_vectorfield_ops(self):
+    def _test_scalarfield_to_vectorfield_ops(self):
         r = Run(self.diag_dir())
         time = 0.0
         Ni = r.GetNi(time)
@@ -330,7 +330,7 @@ def test_scalarfield_to_vectorfield_ops(self):
             self.assertTrue(isinstance(s, VectorField))
             self.assertEqual(s.quantities(), ["x", "y", "z"])
 
-    def test_vectorfield_to_vectorfield_ops(self):
+    def _test_vectorfield_to_vectorfield_ops(self):
         r = Run(self.diag_dir())
         time = 0.0
         B = r.GetB(time)
@@ -340,6 +340,17 @@ def test_vectorfield_to_vectorfield_ops(self):
             self.assertTrue(isinstance(s, VectorField))
             self.assertEqual(s.quantities(), ["x", "y", "z"])
 
+    def test_all(self):
+        """
+        DO NOT RUN MULTIPLE SIMULATIONS!
+        """
+
+        checks = 0
+        for test in [method for method in dir(self) if method.startswith("_test_")]:
+            getattr(self, test)()
+            checks += 1
+        self.assertEqual(checks, 18)  # update if you add new tests
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/src/amr/amr_constants.hpp b/src/amr/amr_constants.hpp
index a7c43570d..00751b414 100644
--- a/src/amr/amr_constants.hpp
+++ b/src/amr/amr_constants.hpp
@@ -1,10 +1,16 @@
-#ifndef AMR_CONSTANTS_HPP
-#define AMR_CONSTANTS_HPP
+#ifndef PHARE_AMR_CONSTANTS_HPP
+#define PHARE_AMR_CONSTANTS_HPP
 
 #include <cstddef>
-namespace PHARE::amr {
+
+namespace PHARE::amr
+{
+
 static std::size_t constexpr refinementRatio = 2;
-}
+
+static std::size_t constexpr MAX_LEVEL_IDX = 9; // 10 levels
+
+} // namespace PHARE::amr
 
 
-#endif // AMR_CONSTANTS_HPP
+#endif // PHARE_AMR_CONSTANTS_HPP
diff --git a/src/amr/physical_models/hybrid_model.hpp b/src/amr/physical_models/hybrid_model.hpp
index 08dfcc878..a485d0a1f 100644
--- a/src/amr/physical_models/hybrid_model.hpp
+++ b/src/amr/physical_models/hybrid_model.hpp
@@ -28,15 +28,15 @@ class HybridModel : public IPhysicalModel<AMR_Types>
     using Interface              = IPhysicalModel<AMR_Types>;
     using amr_types              = AMR_Types;
     using electrons_t            = Electrons;
-    using patch_t                = typename AMR_Types::patch_t;
-    using level_t                = typename AMR_Types::level_t;
+    using patch_t                = AMR_Types::patch_t;
+    using level_t                = AMR_Types::level_t;
     using gridlayout_type        = GridLayoutT;
     using electromag_type        = Electromag;
-    using vecfield_type          = typename Electromag::vecfield_type;
-    using field_type             = typename vecfield_type::field_type;
+    using vecfield_type          = Electromag::vecfield_type;
+    using field_type             = vecfield_type::field_type;
     using grid_type              = Grid_t;
     using ions_type              = Ions;
-    using particle_array_type    = typename Ions::particle_array_type;
+    using particle_array_type    = Ions::particle_array_type;
     using resources_manager_type = amr::ResourcesManager<gridlayout_type, grid_type>;
     using ParticleInitializerFactory
         = core::ParticleInitializerFactory<particle_array_type, gridlayout_type>;
diff --git a/src/amr/resources_manager/amr_utils.hpp b/src/amr/resources_manager/amr_utils.hpp
index f3e3da60c..feaebf92d 100644
--- a/src/amr/resources_manager/amr_utils.hpp
+++ b/src/amr/resources_manager/amr_utils.hpp
@@ -4,8 +4,10 @@
 
 #include "core/def.hpp"
 #include "core/def/phare_mpi.hpp" // IWYU pragma: keep
+#include "core/utilities/mpi_utils.hpp"
 #include "core/utilities/constants.hpp"
 
+#include "amr/amr_constants.hpp"
 #include "amr/types/amr_types.hpp"
 #include "amr/utilities/box/amr_box.hpp"
 
@@ -16,6 +18,8 @@
 #include <SAMRAI/hier/BoxOverlap.h>
 #include <SAMRAI/hier/HierarchyNeighbors.h>
 #include <SAMRAI/geom/CartesianPatchGeometry.h>
+
+#include <algorithm>
 #include <stdexcept>
 
 namespace PHARE
@@ -275,10 +279,69 @@ namespace amr
              ++iLevel)
         {
             visitLevel<GridLayout>(*hierarchy.getPatchLevel(iLevel), resman,
-                                   std::forward<Action>(action), std::forward<Args...>(args...));
+                                   std::forward<Action>(action), std::forward<Args>(args)...);
         }
     }
 
+
+    template<std::size_t dim>
+    auto boxesPerRankOn(auto const& level)
+    {
+        auto const& mapping      = level.getProcessorMapping();
+        auto const& global_boxes = level.getBoxes();
+
+        std::vector<std::vector<core::Box<int, dim>>> boxes_per_rank(core::mpi::size());
+        assert(global_boxes.size() == level.getGlobalNumberOfPatches());
+
+        auto gbox_iter = global_boxes.begin();
+        for (int i = 0; i < level.getGlobalNumberOfPatches(); ++i, ++gbox_iter)
+            boxes_per_rank[mapping.getProcessorAssignment(i)].emplace_back(
+                phare_box_from<dim>(*gbox_iter));
+
+        return boxes_per_rank;
+    }
+
+
+
+    template<typename Action>
+    void onLevels(auto& hierarchy, Action&& action, std::size_t const minlvl = 0,
+                  std::size_t const maxlvl = MAX_LEVEL_IDX)
+    {
+        if (hierarchy.getNumberOfLevels() < 1)
+            throw std::runtime_error("Hierarchy must have a level");
+
+        std::size_t const hier_levels = hierarchy.getNumberOfLevels() - 1; // size vs index
+        std::size_t const max         = hier_levels < maxlvl ? hier_levels : maxlvl;
+
+        for (std::size_t ilvl = minlvl; ilvl <= max; ++ilvl)
+            if (auto lvl = hierarchy.getPatchLevel(ilvl))
+                action(*lvl);
+    }
+
+
+    /* This function differs from the one above in that the "orMissing(int ilvl)" function will be
+       executed when the level does not exist in the hierarchy at the current time, but we want to
+       do something anyway for that "index".
+    */
+    void onLevels(auto& hierarchy, auto&& onLevel, auto&& orMissing, std::size_t const minlvl,
+                  std::size_t const maxlvl)
+    {
+        if (hierarchy.getNumberOfLevels() < 1)
+            throw std::runtime_error("Hierarchy must have a level");
+
+        std::size_t const hier_levels   = hierarchy.getNumberOfLevels();
+        std::size_t const last_existing = hier_levels - 1;
+        std::size_t const max_existing  = std::min(maxlvl, last_existing);
+
+        for (std::size_t ilvl = minlvl; ilvl <= max_existing; ++ilvl)
+            onLevel(*hierarchy.getPatchLevel(ilvl));
+
+        for (std::size_t ilvl = std::max(minlvl, hier_levels); ilvl <= maxlvl; ++ilvl)
+            orMissing(ilvl);
+    }
+
+
+
 } // namespace amr
 } // namespace PHARE
 
diff --git a/src/amr/wrappers/hierarchy.hpp b/src/amr/wrappers/hierarchy.hpp
index 201cdcddf..662473ab0 100644
--- a/src/amr/wrappers/hierarchy.hpp
+++ b/src/amr/wrappers/hierarchy.hpp
@@ -2,7 +2,6 @@
 #define PHARE_AMR_HIERARCHY_HPP
 
 
-
 #include "core/def.hpp"
 #include "core/logger.hpp"
 #include "core/def/phare_mpi.hpp" // IWYU pragma: keep
@@ -12,6 +11,7 @@
 #include "initializer/data_provider.hpp"
 
 #include "amr/samrai.hpp"
+#include "amr/amr_constants.hpp"
 
 
 #include <SAMRAI/algs/TimeRefinementIntegrator.h>
@@ -89,9 +89,7 @@ class HierarchyRestarter
     void closeRestartFile() { SamraiLifeCycle::getRestartManager()->closeRestartFile(); }
 
     NO_DISCARD bool isFromRestart() const
-    {
-        return SamraiLifeCycle::getRestartManager()->isFromRestart();
-    }
+    { return SamraiLifeCycle::getRestartManager()->isFromRestart(); }
 
 private:
     std::optional<std::string> static restartFilePath(auto const& dict)
@@ -119,6 +117,7 @@ class Hierarchy : public HierarchyRestarter, public SAMRAI::hier::PatchHierarchy
     NO_DISCARD auto const& boundaryConditions() const { return boundaryConditions_; }
     NO_DISCARD auto const& cellWidth() const { return cellWidth_; }
     NO_DISCARD auto const& domainBox() const { return domainBox_; }
+    NO_DISCARD auto const& maxLevel() const { return maxLevel_; }
 
 
 
@@ -145,6 +144,7 @@ class Hierarchy : public HierarchyRestarter, public SAMRAI::hier::PatchHierarchy
     std::vector<double> const cellWidth_;
     std::vector<int> const domainBox_;
     std::vector<std::string> boundaryConditions_;
+    std::size_t maxLevel_ = 0;
 };
 
 
@@ -241,7 +241,17 @@ Hierarchy::Hierarchy(initializer::PHAREDict const& dict,
     , cellWidth_(cellWidth.data(), cellWidth.data() + dimension)
     , domainBox_(domainBox.data(), domainBox.data() + dimension)
     , boundaryConditions_(boundaryConditions.data(), boundaryConditions.data() + dimension)
+
 {
+    auto const max_nbr_levels = dict["simulation"]["AMR"]["max_nbr_levels"].template to<int>();
+    if (max_nbr_levels < 1)
+        throw std::runtime_error("Invalid max_nbr_levels, must be >= 1");
+
+    maxLevel_ = max_nbr_levels - 1;
+
+    if (maxLevel_ > MAX_LEVEL_IDX)
+        throw std::runtime_error("Invalid max_nbr_levels, must be <= "
+                                 + std::to_string(MAX_LEVEL_IDX + 1));
 }
 
 
diff --git a/src/core/data/grid/gridlayout.hpp b/src/core/data/grid/gridlayout.hpp
index 56df86052..655d1f38b 100644
--- a/src/core/data/grid/gridlayout.hpp
+++ b/src/core/data/grid/gridlayout.hpp
@@ -2,16 +2,15 @@
 #define PHARE_CORE_GRID_GridLayout_HPP
 
 
-#include "core/hybrid/hybrid_quantities.hpp"
+#include "core/def.hpp"
 #include "core/utilities/types.hpp"
-#include "core/data/field/field.hpp"
-#include "gridlayoutdefs.hpp"
-#include "core/utilities/algorithm.hpp"
 #include "core/utilities/box/box.hpp"
 #include "core/utilities/constants.hpp"
 #include "core/utilities/index/index.hpp"
 #include "core/utilities/point/point.hpp"
-#include "core/def.hpp"
+#include "core/hybrid/hybrid_quantities.hpp"
+
+#include "gridlayoutdefs.hpp"
 
 #include <array>
 #include <tuple>
@@ -878,7 +877,6 @@ namespace core
             return GridLayoutImpl::centering(hybridQuantity);
         }
 
-
         NO_DISCARD constexpr static std::array<std::array<QtyCentering, dimension>, 6>
         centering(HybridQuantity::Tensor hybridQuantity)
         {
@@ -1026,6 +1024,29 @@ namespace core
 
 
 
+        /**
+         * @brief BxToMoments return the indexes and associated coef to compute the linear
+         * interpolation necessary to project Bx onto moments.
+         */
+        NO_DISCARD auto static constexpr BxToMoments() { return GridLayoutImpl::BxToMoments(); }
+
+
+        /**
+         * @brief ByToMoments return the indexes and associated coef to compute the linear
+         * interpolation necessary to project By onto moments.
+         */
+        NO_DISCARD auto static constexpr ByToMoments() { return GridLayoutImpl::ByToMoments(); }
+
+
+        /**
+         * @brief BzToMoments return the indexes and associated coef to compute the linear
+         * interpolation necessary to project Bz onto moments.
+         */
+        NO_DISCARD auto static constexpr BzToMoments() { return GridLayoutImpl::BzToMoments(); }
+
+
+
+
         /**
          * @brief ExToMoments return the indexes and associated coef to compute the linear
          * interpolation necessary to project Ex onto moments.
@@ -1156,30 +1177,23 @@ namespace core
 
 
 
-
-        auto AMRGhostBoxFor(auto const& field) const
+        auto AMRBoxFor(auto const& field) const
         {
             auto const centerings = centering(field);
-            auto const growBy     = [&]() {
-                std::array<int, dimension> arr;
-                for (std::uint8_t i = 0; i < dimension; ++i)
-                    arr[i] = nbrGhosts(centerings[i]);
-                return arr;
-            }();
-            auto ghostBox = grow(AMRBox_, growBy);
+            auto box              = AMRBox_;
             for (std::uint8_t i = 0; i < dimension; ++i)
-                ghostBox.upper[i] += (centerings[i] == QtyCentering::primal) ? 1 : 0;
-            return ghostBox;
+                box.upper[i] += (centerings[i] == QtyCentering::primal) ? 1 : 0;
+            return box;
         }
 
-        auto AMRBoxFor(auto const& field) const
+        auto AMRGhostBoxFor(auto const& field) const
         {
             auto const centerings = centering(field);
-            return grow(AMRGhostBoxFor(field), for_N_make_array<dimension>([&](auto i) {
-                            return -1 * nbrGhosts(centerings[i]);
-                        }));
+            return grow(AMRBoxFor(field), for_N_make_array<dimension>(
+                                              [&](auto i) { return nbrGhosts(centerings[i]); }));
         }
 
+
         template<typename Field, typename Fn>
         void evalOnBox(Field& field, Fn&& fn) const
         {
@@ -1202,6 +1216,15 @@ namespace core
 
         auto levelNumber() const { return levelNumber_; }
 
+
+        auto amr_lcl_idx(auto const& box) const { return boxes_iterator{box, AMRToLocal(box)}; }
+        auto amr_lcl_idx() const { return amr_lcl_idx(AMRBox()); }
+        auto domain_amr_lcl_idx(auto const& field) const { return amr_lcl_idx(AMRBoxFor(field)); }
+        auto ghost_amr_lcl_idx(auto const& field) const
+        {
+            return amr_lcl_idx(AMRGhostBoxFor(field));
+        }
+
     private:
         template<typename Field, typename IndicesFn, typename Fn>
         static void evalOnBox_(Field& field, Fn& fn, IndicesFn& startToEnd)
@@ -1514,51 +1537,6 @@ namespace core
         }
 
 
-
-        struct AMRLocalIndexer //
-        {
-            GridLayout const* layout;
-            Box<int, dimension> amr_box;
-            Box<std::uint32_t, dimension> lcl_box = layout->AMRToLocal(amr_box);
-
-            struct iterator
-            {
-                void operator++()
-                {
-                    ++amr_it;
-                    ++lcl_it;
-                }
-                auto operator*() { return std::forward_as_tuple(*amr_it, *lcl_it); }
-                auto operator==(auto const& that) const
-                {
-                    return amr_it == that.amr_it and lcl_it == that.lcl_it;
-                }
-                auto operator!=(auto const& that) const
-                {
-                    return amr_it != that.amr_it or lcl_it != that.lcl_it;
-                }
-
-                AMRLocalIndexer const* amr_lcl_indexer;
-                box_iterator<int, dimension> amr_it;
-                box_iterator<std::uint32_t, dimension> lcl_it;
-            };
-
-            auto begin() const { return iterator{this, amr_box.begin(), lcl_box.begin()}; }
-            auto end() const { return iterator{this, amr_box.end(), lcl_box.end()}; }
-        };
-
-    public:
-        auto amr_lcl_idx(auto const& box) const { return AMRLocalIndexer{this, box}; }
-        auto amr_lcl_idx() const { return amr_lcl_idx(AMRBox()); }
-        auto domain_amr_lcl_idx(auto const& field) const
-        {
-            return AMRLocalIndexer{this, AMRBoxFor(field)};
-        }
-        auto ghost_amr_lcl_idx(auto const& field) const
-        {
-            return AMRLocalIndexer{this, AMRGhostBoxFor(field)};
-        }
-
     private:
         std::array<double, dimension> meshSize_;
         Point<double, dimension> origin_;
diff --git a/src/core/data/grid/gridlayoutdefs.hpp b/src/core/data/grid/gridlayoutdefs.hpp
index 23920396e..1b3cb02a0 100644
--- a/src/core/data/grid/gridlayoutdefs.hpp
+++ b/src/core/data/grid/gridlayoutdefs.hpp
@@ -1,11 +1,11 @@
 #ifndef PHARE_CORE_GRID_GRIDLAYOUTDEFS_HPP
 #define PHARE_CORE_GRID_GRIDLAYOUTDEFS_HPP
 
-#include <cstddef>
 
-#include "core/hybrid/hybrid_quantities.hpp"
-#include "core/utilities/types.hpp"
 #include "core/utilities/point/point.hpp"
+#include "core/hybrid/hybrid_quantities.hpp"
+
+#include <cstddef>
 
 namespace PHARE
 {
diff --git a/src/core/data/grid/gridlayoutimplyee.hpp b/src/core/data/grid/gridlayoutimplyee.hpp
index d60537315..9a3746fef 100644
--- a/src/core/data/grid/gridlayoutimplyee.hpp
+++ b/src/core/data/grid/gridlayoutimplyee.hpp
@@ -424,19 +424,19 @@ namespace core
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
             else if constexpr (dimension == 2)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift, 0}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
             else if constexpr (dimension == 3)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift, 0, 0}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
         }
 
@@ -457,19 +457,19 @@ namespace core
                 // in 1D the moment is already on Ey so return 1 point with no shift
                 // with coef 1.
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 1.};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             else if constexpr (dimension == 2)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, iShift}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
             else if constexpr (dimension == 3)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, iShift, 0}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
         }
 
@@ -490,19 +490,116 @@ namespace core
                 // in 1D or 2D the moment is already on Ez so return 1 point with no shift
                 // with coef 1.
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 1.};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             else if constexpr (dimension == 2)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 1.};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             else if constexpr (dimension == 3)
             {
                 // in 3D we need two points, the second with a primalToDual shift along Z
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, 0, iShift}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
+            }
+        }
+
+
+
+        NO_DISCARD auto static constexpr BxToMoments()
+        {
+            // Bx is primal dual dual
+            // moments are primal primal primal
+            // operation is thus Pdd to Ppp
+            [[maybe_unused]] auto constexpr iShift = dualToPrimal();
+
+            if constexpr (dimension == 1)
+            {
+                constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 1.};
+                return std::array{P1};
+            }
+            if constexpr (dimension == 2)
+            {
+                constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 0.5};
+                constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, iShift}, 0.5};
+                return std::array{P1, P2};
+            }
+            else if constexpr (dimension == 3)
+            {
+                constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.25};
+                constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, iShift, 0}, 0.25};
+                constexpr WeightPoint<dimension> P3{Point<int, dimension>{0, 0, iShift}, 0.25};
+                constexpr WeightPoint<dimension> P4{Point<int, dimension>{0, iShift, iShift}, 0.25};
+                return std::array{P1, P2, P3, P4};
+            }
+        }
+
+
+
+        NO_DISCARD auto static constexpr ByToMoments()
+        {
+            // By is dual primal dual
+            // moments are primal primal primal
+            // operation is thus Dpd to Ppp
+
+            auto constexpr iShift = dualToPrimal();
+
+            if constexpr (dimension == 1)
+            {
+                constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 0.5};
+                constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift}, 0.5};
+                return std::array{P1, P2};
+            }
+            if constexpr (dimension == 2)
+            {
+                constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 0.5};
+                constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift, 0}, 0.5};
+                return std::array{P1, P2};
+            }
+            else if constexpr (dimension == 3)
+            {
+                constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.25};
+                constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift, 0, 0}, 0.25};
+                constexpr WeightPoint<dimension> P3{Point<int, dimension>{0, 0, iShift}, 0.25};
+                constexpr WeightPoint<dimension> P4{Point<int, dimension>{iShift, 0, iShift}, 0.25};
+                return std::array{P1, P2, P3, P4};
+            }
+        }
+
+
+
+
+        NO_DISCARD auto static constexpr BzToMoments()
+        {
+            // Bz is dual dual primal
+            // moments are primal primal primal
+            // operation is thus Ddp to Ppp
+
+            auto constexpr iShift = dualToPrimal();
+
+            if constexpr (dimension == 1)
+            {
+                constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 0.5};
+                constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift}, 0.5};
+                return std::array{P1, P2};
+            }
+            if constexpr (dimension == 2)
+            {
+                constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 0.25};
+                constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift, 0}, 0.25};
+                constexpr WeightPoint<dimension> P3{Point<int, dimension>{0, iShift}, 0.25};
+                constexpr WeightPoint<dimension> P4{Point<int, dimension>{iShift, iShift}, 0.25};
+                return std::array{P1, P2, P3, P4};
+            }
+            else if constexpr (dimension == 3)
+            {
+                constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.25};
+                constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift, 0, 0}, 0.25};
+                constexpr WeightPoint<dimension> P3{Point<int, dimension>{0, iShift, 0}, 0.25};
+                constexpr WeightPoint<dimension> P4{Point<int, dimension>{iShift, iShift, 0}, 0.25};
+                return std::array{P1, P2, P3, P4};
             }
         }
 
@@ -521,19 +618,19 @@ namespace core
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
             if constexpr (dimension == 2)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift, 0}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
             else if constexpr (dimension == 3)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift, 0, 0}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
         }
 
@@ -551,19 +648,19 @@ namespace core
             if constexpr (dimension == 1)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 1.0};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             if constexpr (dimension == 2)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, iShift}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
             else if constexpr (dimension == 3)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, iShift, 0}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
         }
 
@@ -581,18 +678,18 @@ namespace core
             if constexpr (dimension == 1)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 1.0};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             if constexpr (dimension == 2)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 1.0};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             else if constexpr (dimension == 3)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, 0, iShift}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
         }
 
@@ -610,19 +707,19 @@ namespace core
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
             if constexpr (dimension == 2)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift, 0}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
             else if constexpr (dimension == 3)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift, 0, 0}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
         }
 
@@ -639,19 +736,19 @@ namespace core
             if constexpr (dimension == 1)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 1.0};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             if constexpr (dimension == 2)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, iShift}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
             else if constexpr (dimension == 3)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, iShift, 0}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
         }
 
@@ -669,18 +766,18 @@ namespace core
             if constexpr (dimension == 1)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 1.0};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             if constexpr (dimension == 2)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 1.0};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             else if constexpr (dimension == 3)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, 0, iShift}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
         }
 
@@ -696,7 +793,7 @@ namespace core
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{p2dShift}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
             else if constexpr (dimension == 2)
             {
@@ -705,7 +802,7 @@ namespace core
                 constexpr WeightPoint<dimension> P3{Point<int, dimension>{p2dShift, 0}, 0.25};
                 constexpr WeightPoint<dimension> P4{Point<int, dimension>{p2dShift, d2pShift},
                                                     0.25};
-                return std::array<WeightPoint<dimension>, 4>{P1, P2, P3, P4};
+                return std::array{P1, P2, P3, P4};
             }
             else if constexpr (dimension == 3)
             {
@@ -722,13 +819,14 @@ namespace core
                                                     0.125};
                 constexpr WeightPoint<dimension> P8{
                     Point<int, dimension>{p2dShift, d2pShift, d2pShift}, 0.125};
-                return std::array<WeightPoint<dimension>, 8>{P1, P2, P3, P4, P5, P6, P7, P8};
+                return std::array{P1, P2, P3, P4, P5, P6, P7, P8};
             }
         }
 
 
         NO_DISCARD auto static constexpr ByToEx()
-        { // By is dual primal dual
+        {
+            // By is dual primal dual
             // Ex is dual primal primal
             // operation is thus dpD to dpP
             // shift only in the Z direction
@@ -737,18 +835,18 @@ namespace core
             if constexpr (dimension == 1)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 1.0};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             if constexpr (dimension == 2)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 1.0};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             else if constexpr (dimension == 3)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, 0, iShift}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
         }
 
@@ -766,20 +864,20 @@ namespace core
             if constexpr (dimension == 1)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 1.};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             else if constexpr (dimension == 2)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, iShift}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
 
             else if constexpr (dimension == 3)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, iShift, 0}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
         }
 
@@ -796,7 +894,7 @@ namespace core
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{d2pShift}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
 
             else if constexpr (dimension == 2)
@@ -806,7 +904,7 @@ namespace core
                 constexpr WeightPoint<dimension> P3{Point<int, dimension>{0, d2pShift}, 0.25};
                 constexpr WeightPoint<dimension> P4{Point<int, dimension>{d2pShift, d2pShift},
                                                     0.25};
-                return std::array<WeightPoint<dimension>, 4>{P1, P2, P3, P4};
+                return std::array{P1, P2, P3, P4};
             }
             else if constexpr (dimension == 3)
             {
@@ -822,7 +920,7 @@ namespace core
                                                     0.25};
                 constexpr WeightPoint<dimension> P8{
                     Point<int, dimension>{d2pShift, d2pShift, p2dShift}, 0.25};
-                return std::array<WeightPoint<dimension>, 8>{P1, P2, P3, P4, P5, P6, P7, P8};
+                return std::array{P1, P2, P3, P4, P5, P6, P7, P8};
             }
         }
 
@@ -841,19 +939,19 @@ namespace core
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
             if constexpr (dimension == 2)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift, 0}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
             else if constexpr (dimension == 3)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift, 0, 0}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
         }
 
@@ -871,20 +969,20 @@ namespace core
             if constexpr (dimension == 1)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 1.};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             else if constexpr (dimension == 2)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, iShift}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
 
             else if constexpr (dimension == 3)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, iShift, 0}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
         }
 
@@ -902,19 +1000,19 @@ namespace core
             if constexpr (dimension == 1)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 1};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             if constexpr (dimension == 2)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 1.};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             else if constexpr (dimension == 3)
             {
                 // in 3D we need two points, the second with a dualToPrimal shift along Z
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, 0, iShift}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
         }
 
@@ -931,7 +1029,7 @@ namespace core
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{d2pShift}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
             else if constexpr (dimension == 2)
             {
@@ -940,7 +1038,7 @@ namespace core
                 constexpr WeightPoint<dimension> P3{Point<int, dimension>{0, p2dShift}, 0.25};
                 constexpr WeightPoint<dimension> P4{Point<int, dimension>{d2pShift, p2dShift},
                                                     0.25};
-                return std::array<WeightPoint<dimension>, 4>{P1, P2, P3, P4};
+                return std::array{P1, P2, P3, P4};
             }
             else if constexpr (dimension == 3)
             {
@@ -956,7 +1054,7 @@ namespace core
                                                     0.25};
                 constexpr WeightPoint<dimension> P8{
                     Point<int, dimension>{d2pShift, p2dShift, d2pShift}, 0.25};
-                return std::array<WeightPoint<dimension>, 8>{P1, P2, P3, P4, P5, P6, P7, P8};
+                return std::array{P1, P2, P3, P4, P5, P6, P7, P8};
             }
         }
 
@@ -973,19 +1071,19 @@ namespace core
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
             if constexpr (dimension == 2)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift, 0}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
             else if constexpr (dimension == 3)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.5};
                 constexpr WeightPoint<dimension> P2{Point<int, dimension>{iShift, 0, 0}, 0.5};
-                return std::array<WeightPoint<dimension>, 2>{P1, P2};
+                return std::array{P1, P2};
             }
         }
 
@@ -1001,17 +1099,17 @@ namespace core
             if constexpr (dimension == 1)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 1};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             if constexpr (dimension == 2)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 1.0};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             else if constexpr (dimension == 3)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 1.0};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
         }
 
@@ -1027,17 +1125,17 @@ namespace core
             if constexpr (dimension == 1)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 1};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             if constexpr (dimension == 2)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 1.0};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             else if constexpr (dimension == 3)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 1.0};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
         }
 
@@ -1053,17 +1151,17 @@ namespace core
             if constexpr (dimension == 1)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 1};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             if constexpr (dimension == 2)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 1.0};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
             else if constexpr (dimension == 3)
             {
                 constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 1.0};
-                return std::array<WeightPoint<dimension>, 1>{P1};
+                return std::array{P1};
             }
         }
     }; // namespace core
diff --git a/src/core/data/ndarray/ndarray_vector.hpp b/src/core/data/ndarray/ndarray_vector.hpp
index a008c76ef..5c8b370c3 100644
--- a/src/core/data/ndarray/ndarray_vector.hpp
+++ b/src/core/data/ndarray/ndarray_vector.hpp
@@ -2,81 +2,73 @@
 #define PHARE_CORE_DATA_NDARRAY_NDARRAY_VECTOR_HPP
 
 #include "core/def.hpp"
-#include <stdexcept>
+#include "core/utilities/types.hpp"
+
 #include <array>
-#include <cstdint>
 #include <vector>
-#include <tuple>
+#include <cstdint>
 #include <numeric>
-#include <iostream>
-
-
-#include "core/utilities/types.hpp"
+#include <stdexcept>
 
 
 namespace PHARE::core
 {
-template<std::size_t dim, bool c_ordering = true, typename DataType = double>
+template<std::size_t dim, bool c_ordering = true>
 struct NdArrayViewer
 {
-    template<typename NCells, typename... Indexes>
-    NO_DISCARD static DataType const& at(DataType const* data, NCells const& nCells,
-                                         Indexes const&... indexes)
+    template<typename NCells, template<typename, std::size_t> typename Indexes, typename Index>
+    static inline std::uint32_t idx(NCells const& nCells, Indexes<Index, dim> const& indexes)
     {
-        auto params = std::forward_as_tuple(indexes...);
-        static_assert(sizeof...(Indexes) == dim);
-        // static_assert((... && std::is_unsigned_v<decltype(indexes)>)); TODO : manage later if
-        // this test should be included
-
         if constexpr (dim == 1)
-        {
-            auto i = std::get<0>(params);
+            return idx(nCells, indexes[0]);
 
-            return data[i];
-        }
+        else if constexpr (dim == 2)
+            return idx(nCells, indexes[0], indexes[1]);
 
-        if constexpr (dim == 2)
-        {
-            auto i = std::get<0>(params);
-            auto j = std::get<1>(params);
+        else if constexpr (dim == 3)
+            return idx(nCells, indexes[0], indexes[1], indexes[2]);
+    }
 
-            if constexpr (c_ordering)
-                return data[j + i * nCells[1]];
-            else
-                return data[i + j * nCells[0]];
-        }
+    static inline std::uint32_t idx(auto const& /*nCells*/, std::uint32_t const i) { return i; }
+    static inline std::uint32_t idx(auto const& nCells, std::uint32_t const i,
+                                    std::uint32_t const j)
+    {
+        if constexpr (c_ordering)
+            return j + i * nCells[1];
+        else
+            return i + j * nCells[0];
+    }
 
-        if constexpr (dim == 3)
-        {
-            auto i = std::get<0>(params);
-            auto j = std::get<1>(params);
-            auto k = std::get<2>(params);
-
-            if constexpr (c_ordering)
-                return data[k + j * nCells[2] + i * nCells[1] * nCells[2]];
-            else
-                return data[i + j * nCells[0] + k * nCells[1] * nCells[0]];
-        }
+    static inline std::uint32_t idx(auto const& nCells, std::uint32_t const i,
+                                    std::uint32_t const j, std::uint32_t const k)
+    {
+        if constexpr (c_ordering)
+            return k + j * nCells[2] + i * nCells[1] * nCells[2];
+        else
+            return i + j * nCells[0] + k * nCells[1] * nCells[0];
     }
 
-    template<typename NCells, template<typename, std::size_t> typename Indexes, typename Index>
-    NO_DISCARD static DataType const& at(DataType const* data, NCells const& nCells,
-                                         Indexes<Index, dim> const& indexes)
+    template<template<typename, std::size_t> typename Indexes, typename Index>
+    NO_DISCARD static inline auto& at(auto* data, auto const& nCells,
+                                      Indexes<Index, dim> const& indexes)
 
     {
-        if constexpr (dim == 1)
-            return at(data, nCells, indexes[0]);
-
-        else if constexpr (dim == 2)
-            return at(data, nCells, indexes[0], indexes[1]);
+        auto const& i = idx(nCells, indexes);
+        assert(i < product(nCells, std::uint32_t{1}));
+        return data[i];
+    }
 
-        else if constexpr (dim == 3)
-            return at(data, nCells, indexes[0], indexes[1], indexes[2]);
+    static inline auto& at(auto* data, auto const nCells, auto const... indexes)
+    {
+        auto const& i = idx(nCells, indexes...);
+        assert(i < product(nCells, std::uint32_t{1}));
+        return data[i];
     }
 };
 
 
 
+
 template<typename Array, typename Mask>
 class MaskedView
 {
@@ -102,7 +94,7 @@ class MaskedView
     template<typename... Indexes>
     NO_DISCARD DataType const& operator()(Indexes... indexes) const
     {
-        return NdArrayViewer<dimension, true, DataType>::at(array_.data(), shape_, indexes...);
+        return NdArrayViewer<dimension, true>::at(array_.data(), shape_, indexes...);
     }
 
     template<typename... Indexes>
@@ -139,7 +131,7 @@ class NdArrayView
     static std::size_t const dimension  = dim;
     using type                          = DataType;
     using pointer_type                  = DataType*;
-    using viewer                        = NdArrayViewer<dim, c_ordering, DataType>;
+    using viewer                        = NdArrayViewer<dim, c_ordering>;
 
     explicit NdArrayView(pointer_type ptr, std::array<std::uint32_t, dim> const& nCells)
         : ptr_{ptr}
@@ -292,7 +284,7 @@ class NdArrayVector
     template<typename... Indexes>
     NO_DISCARD DataType const& operator()(Indexes... indexes) const
     {
-        return NdArrayViewer<dim, c_ordering, DataType>::at(data_.data(), nCells_, indexes...);
+        return NdArrayViewer<dim, c_ordering>::at(data_.data(), nCells_, indexes...);
     }
 
     template<typename... Indexes>
@@ -304,7 +296,7 @@ class NdArrayVector
     template<typename Index>
     NO_DISCARD DataType const& operator()(std::array<Index, dim> const& indexes) const
     {
-        return NdArrayViewer<dim, c_ordering, DataType>::at(data_.data(), nCells_, indexes);
+        return NdArrayViewer<dim, c_ordering>::at(data_.data(), nCells_, indexes);
     }
 
     template<typename Index>
diff --git a/src/core/data/tensorfield/tensorfield.hpp b/src/core/data/tensorfield/tensorfield.hpp
index 48a6ca2d8..0a8afd270 100644
--- a/src/core/data/tensorfield/tensorfield.hpp
+++ b/src/core/data/tensorfield/tensorfield.hpp
@@ -205,7 +205,7 @@ class TensorField
 
     NO_DISCARD auto& componentNames() const { return componentNames_; }
     NO_DISCARD auto& physicalQuantity() const { return qty_; }
-    NO_DISCARD auto constexpr static size() { return N; }
+    NO_DISCARD auto static constexpr size() { return N; }
 
 private:
     auto static _get_index_for(Component component)
diff --git a/src/core/utilities/algorithm.hpp b/src/core/utilities/algorithm.hpp
index ae1aba764..0ba51317a 100644
--- a/src/core/utilities/algorithm.hpp
+++ b/src/core/utilities/algorithm.hpp
@@ -1,13 +1,21 @@
-#ifndef PHARE_ALGORITHM_HPP
-#define PHARE_ALGORITHM_HPP
+#ifndef PHARE_CORE_UTILITIES_ALGORITHM_HPP
+#define PHARE_CORE_UTILITIES_ALGORITHM_HPP
+
 
 #include "core/def.hpp"
-#include "core/utilities/span.hpp"
 
+#include "core/utilities/span.hpp"
+#include "core/utilities/types.hpp"
+#include "core/data/field/field.hpp"
+#include "core/utilities/box/box.hpp"
+#include "core/data/grid/gridlayoutdefs.hpp"
+#include "core/data/ndarray/ndarray_vector.hpp"
+#include "core/data/tensorfield/tensorfield.hpp"
 
 
-#include <algorithm>
 #include <string>
+#include <algorithm>
+#include <stdexcept>
 
 namespace PHARE
 {
@@ -64,6 +72,84 @@ void average(Span const& f1, Span const& f2, Span& avg)
         av[i] = (d1[i] + d2[i]) * .5;
 }
 
+
+template<typename PhysicalQuantity>
+auto convert_to_primal(        //
+    auto const& src,           //
+    auto const& layout,        //
+    auto const lix,            //
+    PhysicalQuantity const qty //
+)
+{
+    using PQ = PhysicalQuantity;
+
+    if (qty == PQ::Bx)
+        return layout.project(src, lix, layout.BxToMoments());
+    else if (qty == PQ::By)
+        return layout.project(src, lix, layout.ByToMoments());
+    else if (qty == PQ::Bz)
+        return layout.project(src, lix, layout.BzToMoments());
+
+    else if (qty == PQ::Ex)
+        return layout.project(src, lix, layout.ExToMoments());
+    else if (qty == PQ::Ey)
+        return layout.project(src, lix, layout.EyToMoments());
+    else if (qty == PQ::Ez)
+        return layout.project(src, lix, layout.EzToMoments());
+
+    throw std::runtime_error("Quantity not supported for conversion to primal.");
+}
+
+template<std::size_t dim, typename... Ts>
+auto& convert_to_fortran_primal(  // DOES NOT WORK ON GHOST BOX!
+    Field<dim, Ts...>& dst,       //
+    Field<dim, Ts...> const& src, //
+    auto const& layout            //
+)
+{
+    bool static constexpr c_ordering = false;
+
+    if (not all(layout.centering(dst), [](auto const c) { return c == QtyCentering::primal; }))
+        throw std::runtime_error("Invalid operation, all directions must be primal");
+
+    auto const all_primal
+        = all(layout.centering(src), [](auto const c) { return c == QtyCentering::primal; });
+
+    auto const lcl_box = layout.AMRToLocal(layout.AMRBoxFor(dst));
+    auto dst_view      = make_array_view<c_ordering>(dst.data(), *lcl_box.shape());
+
+    auto const lcl_zero_box = [&]() {
+        auto box = lcl_box;
+        box.lower -= lcl_box.lower; // 0
+        box.upper -= lcl_box.lower; // shift
+        return box;
+    }();
+
+    if (all_primal)
+        for (auto const [lix, lix0] : boxes_iterator(lcl_box, lcl_zero_box))
+            dst_view(lix0) = src(lix);
+
+    else
+        for (auto const [lix, lix0] : boxes_iterator(lcl_box, lcl_zero_box))
+            dst_view(lix0) = convert_to_primal(src, layout, lix, src.physicalQuantity());
+
+    return dst;
+}
+
+template<typename Field_t, typename PQ, std::size_t rank>
+auto& convert_to_fortran_primal(               //
+    TensorField<Field_t, PQ, rank>& dst,       //
+    TensorField<Field_t, PQ, rank> const& src, //
+    auto const& layout                         //
+)
+{
+    for (std::size_t ci = 0; ci < src.size(); ++ci)
+        convert_to_fortran_primal(dst[ci], src[ci], layout);
+    return dst;
+}
+
+
+
 } // namespace PHARE::core
 
-#endif
+#endif // PHARE_CORE_UTILITIES_ALGORITHM_HPP
diff --git a/src/core/utilities/box/box.hpp b/src/core/utilities/box/box.hpp
index 6e55d5196..57bdd9ecf 100644
--- a/src/core/utilities/box/box.hpp
+++ b/src/core/utilities/box/box.hpp
@@ -2,15 +2,16 @@
 #define PHARE_CORE_UTILITIES_BOX_BOX_HPP
 
 
+#include "core/def.hpp"
 #include "core/utilities/types.hpp"
 #include "core/utilities/point/point.hpp"
 #include "core/utilities/meta/meta_utilities.hpp"
-#include "core/def.hpp"
 
+#include <tuple>
 #include <cstddef>
-#include <algorithm>
 #include <optional>
 #include <iostream>
+#include <algorithm>
 
 namespace PHARE::core
 {
@@ -26,6 +27,7 @@ template<typename Type, std::size_t dim>
 struct Box
 {
     using value_type                = Type;
+    using iterator                  = box_iterator<Type, dim>;
     static constexpr auto dimension = dim;
 
 
@@ -93,11 +95,10 @@ struct Box
         return *this;
     }
 
-    NO_DISCARD auto shape() const { return upper - lower + 1; }
-    NO_DISCARD auto size() const { return core::product(shape()); }
 
+    NO_DISCARD auto shape() const { return upper - lower + 1; }
+    NO_DISCARD auto size() const { return core::product(shape(), std::size_t{1}); }
 
-    using iterator = box_iterator<Type, dim>;
     NO_DISCARD auto begin() { return iterator{this, lower}; }
 
     //   // since the 1D scan of the multidimensional box is done assuming C ordering
@@ -188,6 +189,7 @@ class box_iterator
         return *this;
     }
 
+
     bool operator!=(box_iterator const& other) const
     {
         return box_ != other.box_ or index_ != other.index_;
@@ -203,6 +205,59 @@ class box_iterator
 template<typename T, std::size_t s>
 Box(Point<T, s> lower, Point<T, s> upper) -> Box<T, s>;
 
+template<typename... Boxes>
+struct boxes_iterator
+{
+    auto constexpr static N = sizeof...(Boxes);
+
+    boxes_iterator(Boxes const&... boxes)
+        : boxes{std::make_tuple(boxes...)}
+    {
+    }
+
+    struct iterator
+    {
+        using Tuple_t = std::tuple<typename Boxes::iterator...>;
+        static_assert(N == std::tuple_size_v<Tuple_t>);
+
+        iterator(std::tuple<typename Boxes::iterator...> iterators)
+            : its{iterators}
+        {
+        }
+
+        iterator& operator++()
+        {
+            for_N<N>([&](auto i) { ++std::get<i>(its); });
+            return *this;
+        }
+
+        auto operator*()
+        {
+            return for_N<N>([&](auto i) { return *std::get<i>(its); });
+        }
+
+        auto operator!=(iterator const& that) const
+        {
+            return for_N_any<N>([&](auto i) { return std::get<i>(its) != std::get<i>(that.its); });
+        }
+
+        Tuple_t its;
+    };
+
+
+
+    auto begin()
+    {
+        return iterator{for_N<N>([&](auto i) { return std::get<i>(boxes).begin(); })};
+    }
+    auto end()
+    {
+        return iterator{for_N<N>([&](auto i) { return std::get<i>(boxes).end(); })};
+    }
+
+    std::tuple<Boxes...> boxes;
+};
+
 
 
 /** this overload of isIn takes a Point and a Container of boxes
@@ -268,6 +323,7 @@ NO_DISCARD bool isIn(Point<Type, SIZE> const& point, Box<Type, SIZE> const& box)
 
 
 
+
 template<typename Type, std::size_t dim, typename OType>
 Box<Type, dim> grow(Box<Type, dim> const& box, OType const& size)
 {
diff --git a/src/diagnostic/detail/h5typewriter.hpp b/src/diagnostic/detail/h5typewriter.hpp
index 57d450eb0..a05aea449 100644
--- a/src/diagnostic/detail/h5typewriter.hpp
+++ b/src/diagnostic/detail/h5typewriter.hpp
@@ -20,7 +20,7 @@ template<typename Writer>
 class H5TypeWriter : public PHARE::diagnostic::TypeWriter
 {
 public:
-    using Attributes                = typename Writer::Attributes;
+    using Attributes                = Writer::Attributes;
     static constexpr auto dimension = Writer::dimension;
 
     H5TypeWriter(Writer& h5Writer)
diff --git a/src/diagnostic/detail/h5writer.hpp b/src/diagnostic/detail/h5writer.hpp
index 9e636e0b4..e6675d0d3 100644
--- a/src/diagnostic/detail/h5writer.hpp
+++ b/src/diagnostic/detail/h5writer.hpp
@@ -1,17 +1,22 @@
 #ifndef PHARE_DETAIL_DIAGNOSTIC_HIGHFIVE_HPP
 #define PHARE_DETAIL_DIAGNOSTIC_HIGHFIVE_HPP
 
-
-#include "core/data/vecfield/vecfield_component.hpp"
-#include "core/utilities/mpi_utils.hpp"
 #include "core/utilities/types.hpp"
-#include "core/utilities/meta/meta_utilities.hpp"
+#include "core/utilities/mpi_utils.hpp"
+#include "core/data/vecfield/vecfield_component.hpp"
+
+#include "initializer/data_provider.hpp"
 
 #include "hdf5/detail/h5/h5_file.hpp"
 
-#include "diagnostic/detail/h5typewriter.hpp"
-#include "diagnostic/diagnostic_manager.hpp"
+
 #include "diagnostic/diagnostic_props.hpp"
+#include "diagnostic/detail/h5typewriter.hpp"
+#include "diagnostic/detail/types/info.hpp"
+#include "diagnostic/detail/types/meta.hpp"
+#include "diagnostic/detail/types/fluid.hpp"
+#include "diagnostic/detail/types/particle.hpp"
+#include "diagnostic/detail/types/electromag.hpp"
 
 
 #if !defined(PHARE_DIAG_DOUBLES)
@@ -23,17 +28,6 @@ namespace PHARE::diagnostic::h5
 {
 using namespace hdf5::h5;
 
-template<typename Writer>
-class ElectromagDiagnosticWriter;
-template<typename Writer>
-class FluidDiagnosticWriter;
-template<typename Writer>
-class ParticlesDiagnosticWriter;
-template<typename Writer>
-class MetaDiagnosticWriter;
-template<typename Writer>
-class InfoDiagnosticWriter;
-
 
 
 template<typename ModelView>
@@ -57,9 +51,9 @@ class H5Writer
 
     template<typename Hierarchy, typename Model>
     H5Writer(Hierarchy& hier, Model& model, std::string const hifivePath, HiFile::AccessMode _flags)
-        : flags{_flags}
+        : modelView_{hier, model}
+        , flags{_flags}
         , filePath_{hifivePath}
-        , modelView_{hier, model}
     {
     }
 
@@ -168,7 +162,11 @@ class H5Writer
     auto& modelView() { return modelView_; }
     auto timestamp() const { return timestamp_; }
 
-    std::size_t minLevel = 0, maxLevel = 10; // TODO hard-coded to be parametrized somehow
+private:
+    ModelView modelView_;
+
+public:
+    std::size_t minLevel = 0, maxLevel = modelView_.maxLevel();
     HiFile::AccessMode flags;
 
 
@@ -176,7 +174,6 @@ class H5Writer
     double timestamp_ = 0;
     std::string filePath_;
     std::string patchPath_; // is passed around as "virtual write()" has no parameters
-    ModelView modelView_;
     Attributes fileAttributes_;
 
     std::unordered_map<std::string, HiFile::AccessMode> file_flags;
@@ -201,8 +198,8 @@ class H5Writer
 
     H5Writer(H5Writer const&)            = delete;
     H5Writer(H5Writer&&)                 = delete;
-    H5Writer& operator&(H5Writer const&) = delete;
-    H5Writer& operator&(H5Writer&&)      = delete;
+    H5Writer& operator=(H5Writer const&) = delete;
+    H5Writer& operator=(H5Writer&&)      = delete;
 
 
     //  State of this class is controlled via "dump()"
diff --git a/src/diagnostic/detail/vtk_types/electromag.hpp b/src/diagnostic/detail/vtk_types/electromag.hpp
new file mode 100644
index 000000000..00154913d
--- /dev/null
+++ b/src/diagnostic/detail/vtk_types/electromag.hpp
@@ -0,0 +1,103 @@
+#ifndef PHARE_DIAGNOSTIC_DETAIL_VTK_TYPES_ELECTROMAG_HPP
+#define PHARE_DIAGNOSTIC_DETAIL_VTK_TYPES_ELECTROMAG_HPP
+
+#include "diagnostic/detail/vtkh5_type_writer.hpp"
+
+#include <string>
+#include <vector>
+#include <optional>
+#include <unordered_map>
+
+namespace PHARE::diagnostic::vtkh5
+{
+
+template<typename H5Writer>
+class ElectromagDiagnosticWriter : public H5TypeWriter<H5Writer>
+{
+    using Super              = H5TypeWriter<H5Writer>;
+    using VTKFileWriter      = Super::VTKFileWriter;
+    using VTKFileInitializer = Super::VTKFileInitializer;
+
+public:
+    ElectromagDiagnosticWriter(H5Writer& h5Writer)
+        : Super{h5Writer}
+    {
+    }
+
+    void setup(DiagnosticProperties&) override;
+    void write(DiagnosticProperties&) override;
+    void compute(DiagnosticProperties&) override {}
+
+private:
+    struct Info
+    {
+        std::vector<std::size_t> offset_per_level
+            = std::vector<std::size_t>(amr::MAX_LEVEL_IDX + 1);
+    };
+
+    std::unordered_map<std::string, Info> mem;
+};
+
+
+template<typename H5Writer>
+void ElectromagDiagnosticWriter<H5Writer>::setup(DiagnosticProperties& diagnostic)
+{
+    auto& modelView = this->h5Writer_.modelView();
+    VTKFileInitializer initializer{diagnostic, this};
+
+    if (mem.count(diagnostic.quantity) == 0)
+        mem.try_emplace(diagnostic.quantity);
+    auto& info = mem[diagnostic.quantity];
+
+    // assumes exists for all models
+    auto const init = [&](auto const ilvl) -> std::optional<std::size_t> {
+        for (auto* vecField : this->h5Writer_.modelView().getElectromagFields())
+            if (diagnostic.quantity == "/" + vecField->name())
+                return initializer.template initTensorFieldFileLevel<1>(ilvl);
+
+        return std::nullopt;
+    };
+
+    modelView.onLevels(
+        [&](auto const& level) {
+            auto const ilvl = level.getLevelNumber();
+            if (auto const offset = init(ilvl))
+                info.offset_per_level[ilvl] = *offset;
+        },
+        [&](int const ilvl) { // missing level
+            init(ilvl);
+        },
+        this->h5Writer_.minLevel, this->h5Writer_.maxLevel);
+}
+
+
+
+template<typename H5Writer>
+void ElectromagDiagnosticWriter<H5Writer>::write(DiagnosticProperties& diagnostic)
+{
+    auto& modelView = this->h5Writer_.modelView();
+    auto& info      = mem[diagnostic.quantity];
+
+    modelView.onLevels(
+        [&](auto const& level) {
+            auto const ilvl = level.getLevelNumber();
+
+            VTKFileWriter writer{diagnostic, this, info.offset_per_level[ilvl]};
+
+            auto const write_quantity = [&](auto& layout, auto const&, auto const) {
+                PHARE_LOG_SCOPE(3, "ElectromagDiagnosticWriter<H5Writer>::write_quantity");
+
+                for (auto* vecField : this->h5Writer_.modelView().getElectromagFields())
+                    if (diagnostic.quantity == "/" + vecField->name())
+                        writer.template writeTensorField<1>(*vecField, layout);
+            };
+
+            modelView.visitHierarchy(write_quantity, ilvl, ilvl);
+        },
+        this->h5Writer_.minLevel, this->h5Writer_.maxLevel);
+}
+
+
+} // namespace PHARE::diagnostic::vtkh5
+
+#endif /* PHARE_DIAGNOSTIC_DETAIL_VTK_TYPES_ELECTROMAG_HPP */
diff --git a/src/diagnostic/detail/vtk_types/fluid.hpp b/src/diagnostic/detail/vtk_types/fluid.hpp
new file mode 100644
index 000000000..2b0ea8f1c
--- /dev/null
+++ b/src/diagnostic/detail/vtk_types/fluid.hpp
@@ -0,0 +1,246 @@
+#ifndef PHARE_DIAGNOSTIC_DETAIL_VTK_TYPES_FLUID_HPP
+#define PHARE_DIAGNOSTIC_DETAIL_VTK_TYPES_FLUID_HPP
+
+#include "core/logger.hpp"
+
+#include "amr/physical_models/mhd_model.hpp"
+#include "amr/physical_models/hybrid_model.hpp"
+#include "diagnostic/detail/vtkh5_type_writer.hpp"
+
+#include <string>
+#include <vector>
+#include <optional>
+#include <stdexcept>
+#include <unordered_map>
+
+namespace PHARE::diagnostic::vtkh5
+{
+
+template<typename H5Writer>
+class FluidDiagnosticWriter : public H5TypeWriter<H5Writer>
+{
+    using Super              = H5TypeWriter<H5Writer>;
+    using VTKFileWriter      = Super::VTKFileWriter;
+    using VTKFileInitializer = Super::VTKFileInitializer;
+
+public:
+    FluidDiagnosticWriter(H5Writer& h5Writer)
+        : Super{h5Writer}
+    {
+    }
+
+    void setup(DiagnosticProperties&) override;
+    void write(DiagnosticProperties&) override;
+    void compute(DiagnosticProperties&) override {};
+
+private:
+    struct Info
+    {
+        std::vector<std::size_t> offset_per_level
+            = std::vector<std::size_t>(amr::MAX_LEVEL_IDX + 1);
+    };
+
+    struct HybridFluidInitializer
+    {
+        std::optional<std::size_t> operator()(auto const ilvl);
+
+        FluidDiagnosticWriter* writer;
+        DiagnosticProperties& diagnostic;
+        VTKFileInitializer& file_initializer;
+    };
+
+    struct MhdFluidInitializer
+    {
+        std::optional<std::size_t> operator()(auto const ilvl);
+
+        FluidDiagnosticWriter* writer;
+        DiagnosticProperties& diagnostic;
+        VTKFileInitializer& file_initializer;
+    };
+
+    struct HybridFluidWriter
+    {
+        void operator()(auto const& layout);
+
+        FluidDiagnosticWriter* writer;
+        DiagnosticProperties& diagnostic;
+        VTKFileWriter& file_writer;
+    };
+
+    struct MhdFluidWriter
+    {
+        void operator()(auto const& layout);
+
+        FluidDiagnosticWriter* writer;
+        DiagnosticProperties& diagnostic;
+        VTKFileWriter& file_writer;
+    };
+
+    auto static isActiveDiag(DiagnosticProperties const& diagnostic, std::string const& tree,
+                             std::string const& var)
+    {
+        return diagnostic.quantity == tree + var;
+    };
+
+    std::unordered_map<std::string, Info> mem;
+};
+
+
+
+template<typename H5Writer>
+std::optional<std::size_t>
+FluidDiagnosticWriter<H5Writer>::HybridFluidInitializer::operator()(auto const ilvl)
+{
+    auto& modelView = writer->h5Writer_.modelView();
+    auto& ions      = modelView.getIons();
+    std::string const tree{"/ions/"};
+
+    if (isActiveDiag(diagnostic, tree, "charge_density"))
+        return file_initializer.initFieldFileLevel(ilvl);
+    if (isActiveDiag(diagnostic, tree, "mass_density"))
+        return file_initializer.initFieldFileLevel(ilvl);
+    if (isActiveDiag(diagnostic, tree, "bulkVelocity"))
+        return file_initializer.template initTensorFieldFileLevel<1>(ilvl);
+
+    for (auto& pop : ions)
+    {
+        auto const pop_tree = tree + "pop/" + pop.name() + "/";
+        if (isActiveDiag(diagnostic, pop_tree, "density"))
+            return file_initializer.initFieldFileLevel(ilvl);
+        else if (isActiveDiag(diagnostic, pop_tree, "charge_density"))
+            return file_initializer.initFieldFileLevel(ilvl);
+        else if (isActiveDiag(diagnostic, pop_tree, "flux"))
+            return file_initializer.template initTensorFieldFileLevel<1>(ilvl);
+    }
+
+    return std::nullopt;
+}
+
+
+
+template<typename H5Writer>
+std::optional<std::size_t>
+FluidDiagnosticWriter<H5Writer>::MhdFluidInitializer::operator()(auto const ilvl)
+{
+    return std::nullopt;
+}
+
+
+template<typename H5Writer>
+void FluidDiagnosticWriter<H5Writer>::setup(DiagnosticProperties& diagnostic)
+{
+    PHARE_LOG_SCOPE(3, "FluidDiagnosticWriter<H5Writer>::setup");
+
+    using Model_t   = H5Writer::ModelView::Model_t;
+    auto& modelView = this->h5Writer_.modelView();
+
+    VTKFileInitializer initializer{diagnostic, this};
+
+    if (mem.count(diagnostic.quantity) == 0)
+        mem.try_emplace(diagnostic.quantity);
+    auto& info = mem[diagnostic.quantity];
+
+    auto const init = [&](auto const ilvl) -> std::optional<std::size_t> {
+        //
+
+        if constexpr (solver::is_hybrid_model_v<Model_t>)
+            if (auto ret = HybridFluidInitializer{this, diagnostic, initializer}(ilvl))
+                return ret;
+
+        if constexpr (solver::is_mhd_model_v<Model_t>)
+            if (auto ret = MhdFluidInitializer{this, diagnostic, initializer}(ilvl))
+                return ret;
+
+        return std::nullopt;
+    };
+
+    modelView.onLevels(
+        [&](auto const& level) {
+            PHARE_LOG_SCOPE(3, "FluidDiagnosticWriter<H5Writer>::setup_level");
+
+            auto const ilvl = level.getLevelNumber();
+            if (auto const offset = init(ilvl))
+                info.offset_per_level[ilvl] = *offset;
+        },
+        [&](int const ilvl) {
+            PHARE_LOG_SCOPE(3, "FluidDiagnosticWriter<H5Writer>::setup_missing_level");
+
+            init(ilvl);
+        },
+        this->h5Writer_.minLevel, this->h5Writer_.maxLevel);
+}
+
+
+
+template<typename H5Writer>
+void FluidDiagnosticWriter<H5Writer>::HybridFluidWriter::operator()(auto const& layout)
+{
+    auto& modelView = writer->h5Writer_.modelView();
+    auto& ions      = modelView.getIons();
+    std::string const tree{"/ions/"};
+
+    if (isActiveDiag(diagnostic, tree, "charge_density"))
+        file_writer.writeField(ions.chargeDensity(), layout);
+    else if (isActiveDiag(diagnostic, tree, "mass_density"))
+        file_writer.writeField(ions.massDensity(), layout);
+    else if (isActiveDiag(diagnostic, tree, "bulkVelocity"))
+        file_writer.template writeTensorField<1>(ions.velocity(), layout);
+    else
+    {
+        for (auto& pop : ions)
+        {
+            auto const pop_tree = tree + "pop/" + pop.name() + "/";
+            if (isActiveDiag(diagnostic, pop_tree, "density"))
+                file_writer.writeField(pop.particleDensity(), layout);
+            else if (isActiveDiag(diagnostic, pop_tree, "charge_density"))
+                file_writer.writeField(pop.chargeDensity(), layout);
+            else if (isActiveDiag(diagnostic, pop_tree, "flux"))
+                file_writer.template writeTensorField<1>(pop.flux(), layout);
+        }
+    }
+}
+
+
+
+template<typename H5Writer>
+void FluidDiagnosticWriter<H5Writer>::MhdFluidWriter::operator()(auto const& layout)
+{
+    throw std::runtime_error("not implemented");
+}
+
+
+
+template<typename H5Writer>
+void FluidDiagnosticWriter<H5Writer>::write(DiagnosticProperties& diagnostic)
+{
+    PHARE_LOG_SCOPE(3, "FluidDiagnosticWriter<H5Writer>::write");
+
+    using Model_t    = H5Writer::ModelView::Model_t;
+    auto& modelView  = this->h5Writer_.modelView();
+    auto const& info = mem[diagnostic.quantity];
+
+    modelView.onLevels(
+        [&](auto const& level) {
+            auto const ilvl = level.getLevelNumber();
+
+            VTKFileWriter writer{diagnostic, this, info.offset_per_level[ilvl]};
+
+            auto const write_quantity = [&](auto& layout, auto const&, auto const) {
+                if constexpr (solver::is_hybrid_model_v<Model_t>)
+                    HybridFluidWriter{this, diagnostic, writer}(layout);
+
+                if constexpr (solver::is_mhd_model_v<Model_t>)
+                    MhdFluidWriter{this, diagnostic, writer}(layout);
+            };
+
+            modelView.visitHierarchy(write_quantity, ilvl, ilvl);
+        },
+        this->h5Writer_.minLevel, this->h5Writer_.maxLevel);
+}
+
+
+} // namespace PHARE::diagnostic::vtkh5
+
+
+
+#endif /* PHARE_DIAGNOSTIC_DETAIL_VTK_TYPES_FLUID_HPP */
diff --git a/src/diagnostic/detail/vtkh5_type_writer.hpp b/src/diagnostic/detail/vtkh5_type_writer.hpp
new file mode 100644
index 000000000..32a08f1c9
--- /dev/null
+++ b/src/diagnostic/detail/vtkh5_type_writer.hpp
@@ -0,0 +1,483 @@
+#ifndef PHARE_DIAGNOSTIC_DETAIL_VTK_H5_TYPE_WRITER_HPP
+#define PHARE_DIAGNOSTIC_DETAIL_VTK_H5_TYPE_WRITER_HPP
+
+
+#include "core/logger.hpp"
+#include "core/utilities/types.hpp"
+#include "core/utilities/algorithm.hpp"
+#include "core/utilities/mpi_utils.hpp"
+#include "core/data/tensorfield/tensorfield.hpp"
+
+#include "amr/resources_manager/amr_utils.hpp"
+
+#include "diagnostic/diagnostic_writer.hpp"
+
+#include "hdf5/detail/h5/h5_file.hpp"
+
+#include <string>
+#include <unordered_map>
+
+#if !defined(PHARE_DIAG_DOUBLES)
+#error // PHARE_DIAG_DOUBLES not defined
+#endif
+
+
+namespace PHARE::diagnostic::vtkh5
+{
+
+using namespace hdf5::h5;
+
+
+template<std::size_t dim>
+struct HierarchyData
+{
+    // data is duplicated in lower dimensions to fit a 3d view
+    //  so 2d data is * 2
+    //  and 1d data is * 4
+    constexpr static auto X_TIMES = std::array{4, 2, 1}[dim - 1];
+
+    static auto& INSTANCE()
+    {
+        static HierarchyData data;
+        return data;
+    }
+
+    static void reset(auto& h5Writer) // called once per dump
+    {
+        PHARE_LOG_SCOPE(3, "HierarchyData<H5Writer::reset");
+
+        auto const maxLevels = h5Writer.modelView().maxLevel() + 1;
+        auto& data           = INSTANCE();
+
+        data.level_data_size.resize(maxLevels);
+        data.n_boxes_per_level.resize(maxLevels);
+        data.level_boxes_per_rank.resize(maxLevels);
+        data.level_rank_data_size.resize(maxLevels);
+        data.flattened_lcl_level_boxes.resize(maxLevels);
+
+        h5Writer.modelView().onLevels(
+            [&](auto const& level) {
+                auto const ilvl = level.getLevelNumber();
+
+                data.n_boxes_per_level[ilvl]    = 0;
+                data.level_boxes_per_rank[ilvl] = amr::boxesPerRankOn<dim>(level);
+                data.flattened_lcl_level_boxes[ilvl]
+                    = flatten_boxes(data.level_boxes_per_rank[ilvl][core::mpi::rank()]);
+
+                for (std::size_t i = 0; i < data.level_boxes_per_rank[ilvl].size(); ++i)
+                {
+                    data.level_rank_data_size[ilvl].resize(core::mpi::size());
+
+                    data.level_rank_data_size[ilvl][i] = 0;
+                    for (auto box : data.level_boxes_per_rank[ilvl][i])
+                    {
+                        box.upper += 1;                                             // all primal
+                        data.level_rank_data_size[ilvl][i] += box.size() * X_TIMES; // no ghosts
+                    }
+                    data.n_boxes_per_level[ilvl] += data.level_boxes_per_rank[ilvl][i].size();
+                }
+
+                data.level_data_size[ilvl] = core::sum(data.level_rank_data_size[ilvl]);
+            },
+            [&](int const ilvl) { // ilvl does not exist currently
+                data.level_data_size[ilvl]   = 0;
+                data.n_boxes_per_level[ilvl] = 0;
+                data.level_rank_data_size[ilvl].clear();
+                data.level_rank_data_size[ilvl].resize(core::mpi::size());
+                data.level_boxes_per_rank[ilvl].clear();
+                data.level_boxes_per_rank[ilvl].resize(core::mpi::size());
+                data.flattened_lcl_level_boxes[ilvl].clear();
+            },
+            h5Writer.minLevel, h5Writer.maxLevel);
+    }
+
+
+    auto static flatten_boxes(auto const& boxes)
+    {
+        std::vector<std::array<int, dim * 2>> data(boxes.size());
+        std::size_t pos = 0;
+        for (auto const& box : boxes)
+        {
+            for (std::size_t i = 0; i < dim; ++i)
+            {
+                data[pos][0 + i * 2] = box.lower[i];
+                data[pos][1 + i * 2] = box.upper[i];
+            }
+            ++pos;
+        }
+        return data;
+    }
+
+    std::vector<std::vector<std::array<int, dim * 2>>> flattened_lcl_level_boxes;
+    std::vector<std::vector<std::vector<core::Box<int, dim>>>> level_boxes_per_rank;
+    std::vector<std::vector<int>> level_rank_data_size;
+    std::vector<int> n_boxes_per_level, level_data_size;
+};
+
+
+template<typename Writer>
+class H5TypeWriter : public PHARE::diagnostic::TypeWriter
+{
+    using FloatType                      = std::conditional_t<PHARE_DIAG_DOUBLES, double, float>;
+    using HierData                       = HierarchyData<Writer::dimension>;
+    using ModelView                      = Writer::ModelView;
+    using physical_quantity_type         = ModelView::Field::physical_quantity_type;
+    using Attributes                     = Writer::Attributes;
+    std::string static inline const base = "/VTKHDF";
+    std::string static inline const level_base = base + "/Level";
+    std::string static inline const step_level = base + "/Steps/Level";
+    auto static inline const level_data_path
+        = [](auto const ilvl) { return level_base + std::to_string(ilvl) + "/PointData/data"; };
+
+
+public:
+    static constexpr auto dimension = Writer::dimension;
+
+    H5TypeWriter(Writer& h5Writer)
+        : h5Writer_{h5Writer}
+    {
+    }
+
+    virtual void setup(DiagnosticProperties&) = 0;
+
+    void writeFileAttributes(DiagnosticProperties const&, Attributes const&);
+
+    //------ valid for all h5type writers -------------------------------------
+    void finalize(DiagnosticProperties& diagnostic)
+    {
+        // we close the file by removing the associated file
+        // from the map. This is done only at flush time otherwise
+        ++diagnostic.dumpIdx;
+
+        assert(diagnostic.params.contains("flush_every"));
+
+        std::size_t flushEvery = diagnostic.param<std::size_t>("flush_every");
+
+        if (flushEvery != Writer::flush_never and diagnostic.dumpIdx % flushEvery == 0)
+        {
+            fileData_.erase(diagnostic.quantity);
+            assert(fileData_.count(diagnostic.quantity) == 0);
+        }
+    }
+
+protected:
+    class VTKFileInitializer;
+
+    class VTKFileWriter;
+
+
+    auto& getOrCreateH5File(DiagnosticProperties const& diagnostic)
+    {
+        if (!fileExistsFor(diagnostic))
+            fileData_.emplace(diagnostic.quantity, this->h5Writer_.makeFile(diagnostic));
+        return *fileData_.at(diagnostic.quantity);
+    }
+
+    bool fileExistsFor(DiagnosticProperties const& diagnostic) const
+    {
+        return fileData_.count(diagnostic.quantity);
+    }
+
+
+    Writer& h5Writer_;
+    std::unordered_map<std::string, std::unique_ptr<HighFiveFile>> fileData_;
+};
+
+
+template<typename Writer>
+void H5TypeWriter<Writer>::writeFileAttributes(DiagnosticProperties const& prop,
+                                               Attributes const& attrs)
+{
+    if (fileExistsFor(prop)) // otherwise qty not supported (yet)
+        h5Writer_.writeGlobalAttributeDict(*fileData_.at(prop.quantity), attrs, "/");
+}
+
+
+template<typename Writer>
+class H5TypeWriter<Writer>::VTKFileInitializer
+{
+    auto static constexpr primal_qty         = physical_quantity_type::rho;
+    std::size_t static constexpr boxValsIn3D = 6; // lo0, up0, lo1, up1, lo2, up2
+public:
+    VTKFileInitializer(DiagnosticProperties const& prop, H5TypeWriter<Writer>* const typewriter);
+
+    void initFileLevel(int const ilvl);
+
+    std::size_t initFieldFileLevel(auto const ilvl) { return initAnyFieldLevel(ilvl); }
+
+    template<std::size_t rank = 2>
+    std::size_t initTensorFieldFileLevel(auto const ilvl)
+    {
+        return initAnyFieldLevel<core::detail::tensor_field_dim_from_rank<rank>()>(ilvl);
+    }
+
+private:
+    auto level_spacing(int const lvl) const
+    {
+        auto const mesh_size = typewriter->h5Writer_.modelView().cellWidth();
+        return core::for_N_make_array<dimension>(
+            [&](auto i) { return static_cast<float>(mesh_size[i] / std::pow(2, lvl)); });
+    }
+
+    template<std::size_t N = 1>
+    std::size_t initAnyFieldLevel(auto const ilvl)
+    {
+        h5file.create_resizable_2d_data_set<FloatType, N>(level_data_path(ilvl));
+        initFileLevel(ilvl);
+        resize_boxes(ilvl);
+        resize_data<N>(ilvl);
+        return data_offset;
+    }
+
+
+    template<std::size_t N = 1>
+    void resize_data(int const ilvl);
+
+    void resize_boxes(int const ilvl);
+
+    bool const newFile;
+    DiagnosticProperties const& diagnostic;
+    H5TypeWriter<Writer>* const typewriter;
+    HighFiveFile& h5file;
+    std::size_t data_offset;
+    ModelView& modelView = typewriter->h5Writer_.modelView();
+};
+
+
+template<typename Writer>
+class H5TypeWriter<Writer>::VTKFileWriter
+{
+    auto static constexpr primal_qty         = physical_quantity_type::rho;
+    std::size_t static constexpr boxValsIn3D = 6; // lo0, up0, lo1, up1, lo2, up2
+
+
+public:
+    VTKFileWriter(DiagnosticProperties const& prop, H5TypeWriter<Writer>* const typewriter,
+                  std::size_t const offset);
+
+    void writeField(auto const& field, auto const& layout);
+
+    template<std::size_t rank = 2>
+    void writeTensorField(auto const& tf, auto const& layout);
+
+
+
+private:
+    auto local_box(auto const& layout) const
+    {
+        return layout.AMRToLocal(layout.AMRBoxFor(primal_qty));
+    }
+
+    DiagnosticProperties const& diagnostic;
+    H5TypeWriter<Writer>* const typewriter;
+    HighFiveFile& h5file;
+    std::size_t data_offset = 0;
+    ModelView& modelView    = typewriter->h5Writer_.modelView();
+};
+
+
+template<typename Writer>
+H5TypeWriter<Writer>::VTKFileInitializer::VTKFileInitializer(DiagnosticProperties const& prop,
+                                                             H5TypeWriter<Writer>* const tw)
+    : newFile{!tw->fileExistsFor(prop)}
+    , diagnostic{prop}
+    , typewriter{tw}
+    , h5file{tw->getOrCreateH5File(prop)}
+{
+    {
+        PHARE_LOG_SCOPE(3, "VTKFileInitializer::VTKFileInitializer::0");
+        h5file.create_resizable_1d_data_set<FloatType>(base + "/Steps/Values");
+        auto steps_group = h5file.file().getGroup(base + "/Steps");
+        if (!steps_group.hasAttribute("NSteps"))
+            steps_group.template createAttribute<int>("NSteps", HighFive::DataSpace::From(0))
+                .write(0);
+        auto steps_attr = steps_group.getAttribute("NSteps");
+        steps_attr.write(steps_attr.template read<int>() + 1);
+    }
+
+    {
+        PHARE_LOG_SCOPE(3, "VTKFileInitializer::VTKFileInitializer::1");
+        auto const& timestamp = typewriter->h5Writer_.timestamp();
+        auto ds               = h5file.getDataSet(base + "/Steps/Values");
+        auto const old_size   = ds.getDimensions()[0];
+        ds.resize({old_size + 1});
+        ds.select({old_size}, {1}).write(timestamp);
+    }
+
+    if (newFile)
+    {
+        PHARE_LOG_SCOPE(3, "VTKFileInitializer::VTKFileInitializer::2");
+        auto root = h5file.file().getGroup(base);
+
+        if (!root.hasAttribute("Version"))
+            root.template createAttribute<std::array<int, 2>>("Version", std::array<int, 2>{2, 2});
+
+        if (!root.hasAttribute("Origin"))
+            root.template createAttribute<std::array<float, 3>>("Origin",
+                                                                std::array<float, 3>{0, 0, 0});
+
+        if (!root.hasAttribute("Type"))
+            root.template createAttribute<std::string>("Type", "OverlappingAMR");
+
+        if (!root.hasAttribute("GridDescription"))
+            root.template createAttribute<std::string>("GridDescription", "XYZ");
+    }
+}
+
+
+template<typename Writer>
+H5TypeWriter<Writer>::VTKFileWriter::VTKFileWriter(DiagnosticProperties const& prop,
+                                                   H5TypeWriter<Writer>* const tw,
+                                                   std::size_t const offset)
+    : diagnostic{prop}
+    , typewriter{tw}
+    , h5file{tw->getOrCreateH5File(prop)}
+    , data_offset{offset}
+{
+}
+
+
+template<typename Writer>
+void H5TypeWriter<Writer>::VTKFileWriter::writeField(auto const& field, auto const& layout)
+{
+    PHARE_LOG_SCOPE(3, "VTKFileWriter::writeField");
+
+    auto const& frimal = core::convert_to_fortran_primal(modelView.tmpField(), field, layout);
+    auto const size    = local_box(layout).size();
+    auto ds            = h5file.getDataSet(level_data_path(layout.levelNumber()));
+
+    PHARE_LOG_SCOPE(3, "VTKFileWriter::writeField::0");
+    for (std::uint16_t i = 0; i < HierData::X_TIMES; ++i)
+    {
+        ds.select({data_offset, 0}, {size, 1}).write_raw(frimal.data());
+        data_offset += size;
+    }
+}
+
+
+template<typename Writer>
+template<std::size_t rank>
+void H5TypeWriter<Writer>::VTKFileWriter::writeTensorField(auto const& tf, auto const& layout)
+{
+    PHARE_LOG_SCOPE(3, "VTKFileWriter::writeTensorField");
+
+    auto static constexpr N = core::detail::tensor_field_dim_from_rank<rank>();
+    auto const& frimal
+        = core::convert_to_fortran_primal(modelView.template tmpTensorField<rank>(), tf, layout);
+    auto const size = local_box(layout).size();
+    auto ds         = h5file.getDataSet(level_data_path(layout.levelNumber()));
+
+    PHARE_LOG_SCOPE(3, "VTKFileWriter::writeTensorField::0");
+    for (std::uint16_t i = 0; i < HierData::X_TIMES; ++i)
+    {
+        for (std::uint32_t c = 0; c < N; ++c)
+            ds.select({data_offset, c}, {size, 1}).write_raw(frimal[c].data());
+        data_offset += size;
+    }
+}
+
+
+template<typename Writer>
+void H5TypeWriter<Writer>::VTKFileInitializer::initFileLevel(int const ilvl)
+{
+    PHARE_LOG_SCOPE(3, "VTKFileInitializer::initFileLevel");
+
+    auto const lvl = std::to_string(ilvl);
+
+    h5file.create_resizable_2d_data_set<int, boxValsIn3D>(level_base + lvl + "/AMRBox");
+    h5file.create_resizable_1d_data_set<int>(step_level + lvl + "/AMRBoxOffset");
+    h5file.create_resizable_1d_data_set<int>(step_level + lvl + "/NumberOfAMRBox");
+    h5file.create_resizable_1d_data_set<int>(step_level + lvl + "/PointDataOffset/data");
+
+    auto level_group = h5file.file().getGroup(level_base + lvl);
+    if (!level_group.hasAttribute("Spacing"))
+    {
+        PHARE_LOG_SCOPE(3, "VTKFileInitializer::initFileLevel::0");
+
+        level_group.template createAttribute<std::array<float, dimension>>("Spacing",
+                                                                           level_spacing(ilvl));
+
+        level_group.createGroup("CellData");
+        level_group.createGroup("FieldData");
+
+        auto steps_group = h5file.file().getGroup(step_level + lvl);
+        steps_group.createGroup("CellDataOffset");
+        steps_group.createGroup("FieldDataOffset");
+    }
+}
+
+
+template<typename Writer>
+template<std::size_t N>
+void H5TypeWriter<Writer>::VTKFileInitializer::resize_data(int const ilvl)
+{
+    PHARE_LOG_SCOPE(3, "VTKFileInitializer::resize_data");
+
+    auto const& hier_data = HierData::INSTANCE();
+    auto const lvl        = std::to_string(ilvl);
+    auto const data_path  = level_data_path(ilvl);
+    auto point_data_ds    = h5file.getDataSet(data_path);
+    data_offset           = point_data_ds.getDimensions()[0];
+
+    {
+        PHARE_LOG_SCOPE(3, "VTKFileInitializer::resize_data::0");
+        auto ds             = h5file.getDataSet(step_level + lvl + "/PointDataOffset/data");
+        auto const old_size = ds.getDimensions()[0];
+        ds.resize({old_size + 1});
+        ds.select({old_size}, {1}).write(data_offset);
+    }
+
+    PHARE_LOG_SCOPE(3, "VTKFileInitializer::resize_data::1");
+    auto const& rank_data_sizes = hier_data.level_rank_data_size[ilvl];
+    auto const& level_data_size = hier_data.level_data_size[ilvl];
+    auto const new_size         = data_offset + level_data_size;
+
+    point_data_ds.resize({new_size, N});
+    for (int i = 0; i < core::mpi::rank(); ++i)
+        data_offset += rank_data_sizes[i];
+}
+
+
+template<typename Writer>
+void H5TypeWriter<Writer>::VTKFileInitializer::resize_boxes(int const ilvl)
+{
+    PHARE_LOG_SCOPE(3, "VTKFileInitializer::resize_boxes");
+
+    auto const lvl          = std::to_string(ilvl);
+    auto const& hier_data   = HierData::INSTANCE();
+    auto const& rank_boxes  = hier_data.level_boxes_per_rank[ilvl];
+    auto const& total_boxes = hier_data.n_boxes_per_level[ilvl];
+
+    {
+        PHARE_LOG_SCOPE(3, "VTKFileInitializer::resize_boxes::0");
+        auto ds             = h5file.getDataSet(step_level + lvl + "/NumberOfAMRBox");
+        auto const old_size = ds.getDimensions()[0];
+        ds.resize({old_size + 1});
+        ds.select({old_size}, {1}).write(total_boxes);
+    }
+
+    auto amrbox_ds  = h5file.getDataSet(level_base + lvl + "/AMRBox");
+    auto box_offset = amrbox_ds.getDimensions()[0];
+
+    {
+        PHARE_LOG_SCOPE(3, "VTKFileInitializer::resize_boxes::1");
+        auto ds             = h5file.getDataSet(step_level + lvl + "/AMRBoxOffset");
+        auto const old_size = ds.getDimensions()[0];
+        ds.resize({old_size + 1});
+        ds.select({old_size}, {1}).write(box_offset);
+    }
+
+    PHARE_LOG_SCOPE(3, "VTKFileInitializer::resize_boxes::2");
+    amrbox_ds.resize({box_offset + total_boxes, boxValsIn3D});
+    for (int i = 0; i < core::mpi::rank(); ++i)
+        box_offset += rank_boxes[i].size();
+
+    PHARE_LOG_SCOPE(3, "VTKFileInitializer::resize_boxes::3");
+    amrbox_ds.select({box_offset, 0}, {rank_boxes[core::mpi::rank()].size(), dimension * 2})
+        .write(hier_data.flattened_lcl_level_boxes[ilvl]);
+}
+
+
+} // namespace PHARE::diagnostic::vtkh5
+
+#endif // PHARE_DIAGNOSTIC_DETAIL_VTK_H5_TYPE_WRITER_HPP
diff --git a/src/diagnostic/detail/vtkh5_writer.hpp b/src/diagnostic/detail/vtkh5_writer.hpp
new file mode 100644
index 000000000..352096bbc
--- /dev/null
+++ b/src/diagnostic/detail/vtkh5_writer.hpp
@@ -0,0 +1,221 @@
+#ifndef PHARE_DIAGNOSTIC_DETAIL_VTK_H5_WRITER_HPP
+#define PHARE_DIAGNOSTIC_DETAIL_VTK_H5_WRITER_HPP
+
+
+#include "core/logger.hpp"
+#include "core/utilities/mpi_utils.hpp"
+
+#include "initializer/data_provider.hpp"
+
+#include "diagnostic/diagnostic_props.hpp"
+#include "diagnostic/detail/vtkh5_type_writer.hpp"
+
+#include "diagnostic/detail/vtk_types/fluid.hpp"
+#include "diagnostic/detail/vtk_types/electromag.hpp"
+
+#include "hdf5/detail/h5/h5_file.hpp"
+
+#include <memory>
+#include <string>
+#include <algorithm>
+#include <unordered_map>
+
+namespace PHARE::diagnostic::vtkh5
+{
+using namespace hdf5::h5;
+
+
+
+template<typename _ModelView>
+class H5Writer
+{
+    struct NullTypeWriter : public H5TypeWriter<H5Writer<_ModelView>>
+    {
+        NullTypeWriter(auto& h5Writer)
+            : H5TypeWriter<H5Writer<_ModelView>>{h5Writer}
+        {
+        }
+
+        void setup(DiagnosticProperties& prop) {}
+        void write(DiagnosticProperties& prop)
+        {
+            if (core::mpi::rank() == 0)
+            {
+                PHARE_LOG_LINE_SS( //
+                    "No diagnostic writer found for " + prop.type + ":" + prop.quantity);
+            }
+        }
+        void compute(DiagnosticProperties&) {}
+    };
+
+public:
+    using ModelView  = _ModelView;
+    using This       = H5Writer<ModelView>;
+    using GridLayout = ModelView::GridLayout;
+    using Attributes = ModelView::PatchProperties;
+
+    static constexpr auto dimension   = GridLayout::dimension;
+    static constexpr auto interpOrder = GridLayout::interp_order;
+    static constexpr auto READ_WRITE  = HiFile::AccessMode::OpenOrCreate;
+
+    // flush_never: disables manual file closing, but still occurrs via RAII
+    static constexpr std::size_t flush_never = 0;
+
+    template<typename Hierarchy, typename Model>
+    H5Writer(Hierarchy& hier, Model& model, std::string const hifivePath, HiFile::AccessMode _flags)
+        : modelView_{hier, model}
+        , flags{_flags}
+        , filePath_{hifivePath}
+    {
+    }
+
+    ~H5Writer() {}
+    H5Writer(H5Writer const&)            = delete;
+    H5Writer(H5Writer&&)                 = delete;
+    H5Writer& operator=(H5Writer const&) = delete;
+    H5Writer& operator=(H5Writer&&)      = delete;
+
+
+    template<typename Hierarchy, typename Model>
+    static auto make_unique(Hierarchy& hier, Model& model, initializer::PHAREDict const& dict)
+    {
+        std::string filePath     = dict["filePath"].template to<std::string>();
+        HiFile::AccessMode flags = READ_WRITE;
+        if (dict.contains("mode") and dict["mode"].template to<std::string>() == "overwrite")
+            flags |= HiFile::Truncate;
+        return std::make_unique<This>(hier, model, filePath, flags);
+    }
+
+
+    void dump(std::vector<DiagnosticProperties*> const&, double current_timestamp);
+    void dump_level(std::size_t level, std::vector<DiagnosticProperties*> const& diagnostics,
+                    double timestamp);
+
+    template<typename String>
+    auto getDiagnosticWriterForType(String& type)
+    {
+        return typeWriters_.at(type);
+    }
+
+
+    auto makeFile(DiagnosticProperties const& diagnostic)
+    {
+        return std::make_unique<HighFiveFile>(filePath_ + "/" + fileString(diagnostic.quantity),
+                                              file_flags[diagnostic.type + diagnostic.quantity]);
+    }
+
+    template<typename Dict>
+    static void writeGlobalAttributeDict(HighFiveFile& h5, Dict const& dict, std::string path)
+    {
+        dict.visit(
+            [&](std::string const& key, auto const& val) { h5.write_attribute(path, key, val); });
+    }
+
+
+    auto& modelView() { return modelView_; }
+    auto timestamp() const { return timestamp_; }
+
+private:
+    ModelView modelView_;
+
+public:
+    std::size_t minLevel = 0, maxLevel = modelView_.maxLevel();
+    HiFile::AccessMode flags;
+
+
+private:
+    double timestamp_ = 0;
+    std::string filePath_;
+    Attributes fileAttributes_;
+
+    std::unordered_map<std::string, HiFile::AccessMode> file_flags;
+
+    std::unordered_map<std::string, std::shared_ptr<H5TypeWriter<This>>> typeWriters_{
+        {"info", make_writer<NullTypeWriter>()},
+        {"meta", make_writer<NullTypeWriter>()},
+        {"fluid", make_writer<FluidDiagnosticWriter<This>>()},
+        {"electromag", make_writer<ElectromagDiagnosticWriter<This>>()},
+        {"particle", make_writer<NullTypeWriter>()} //
+    };
+
+    template<typename Writer>
+    std::shared_ptr<H5TypeWriter<This>> make_writer()
+    {
+        return std::make_shared<Writer>(*this);
+    }
+
+
+    static std::string fileString(std::string fileStr)
+    {
+        if (fileStr[0] == '/')
+            fileStr = fileStr.substr(1);
+        std::replace(fileStr.begin(), fileStr.end(), '/', '_');
+        return fileStr + ".vtkhdf";
+    }
+
+
+    //  State of this class is controlled via "dump()"
+    //  block public access to internal state
+    friend class H5TypeWriter<This>;
+    friend class FluidDiagnosticWriter<This>;
+    friend class ElectromagDiagnosticWriter<This>;
+};
+
+
+
+template<typename ModelView>
+void H5Writer<ModelView>::dump(std::vector<DiagnosticProperties*> const& diagnostics,
+                               double timestamp)
+{
+    timestamp_                             = timestamp;
+    fileAttributes_["dimension"]           = dimension;
+    fileAttributes_["interpOrder"]         = interpOrder;
+    fileAttributes_["domain_box"]          = modelView_.domainBox();
+    fileAttributes_["boundary_conditions"] = modelView_.boundaryConditions();
+
+    HierarchyData<dimension>::reset(*this);
+
+    for (auto* diagnostic : diagnostics)
+        if (!file_flags.count(diagnostic->type + diagnostic->quantity))
+            file_flags[diagnostic->type + diagnostic->quantity] = this->flags;
+
+    for (auto* diagnostic : diagnostics) // all collective calls first!
+    {
+        auto& type_writer = *typeWriters_.at(diagnostic->type);
+        type_writer.setup(*diagnostic);
+        type_writer.writeFileAttributes(*diagnostic, fileAttributes_);
+    }
+
+    for (auto* diagnostic : diagnostics)
+        typeWriters_.at(diagnostic->type)->write(*diagnostic);
+
+    for (auto* diagnostic : diagnostics)
+    {
+        typeWriters_.at(diagnostic->type)->finalize(*diagnostic);
+        // don't truncate past first dump
+        file_flags[diagnostic->type + diagnostic->quantity] = READ_WRITE;
+    }
+}
+
+template<typename ModelView>
+void H5Writer<ModelView>::dump_level(std::size_t level,
+                                     std::vector<DiagnosticProperties*> const& diagnostics,
+                                     double timestamp)
+{
+    std::size_t _minLevel = this->minLevel;
+    std::size_t _maxLevel = this->maxLevel;
+
+    this->minLevel = level;
+    this->maxLevel = level;
+
+    this->dump(diagnostics, timestamp);
+
+    this->minLevel = _minLevel;
+    this->maxLevel = _maxLevel;
+}
+
+
+
+} // namespace PHARE::diagnostic::vtkh5
+
+#endif /* PHARE_DIAGNOSTIC_DETAIL_VTK_H5_WRITER_HPP */
diff --git a/src/diagnostic/diagnostic_manager.hpp b/src/diagnostic/diagnostic_manager.hpp
index a7d0f0ea6..14c6560b2 100644
--- a/src/diagnostic/diagnostic_manager.hpp
+++ b/src/diagnostic/diagnostic_manager.hpp
@@ -31,7 +31,7 @@ void registerDiagnostics(DiagManager& dMan, initializer::PHAREDict const& diagsP
         while (diagsParams.contains(diagType)
                && diagsParams[diagType].contains(diagType + std::to_string(diagBlockID)))
         {
-            const std::string diagName = diagType + std::to_string(diagBlockID);
+            std::string const diagName = diagType + std::to_string(diagBlockID);
             dMan.addDiagDict(diagsParams[diagType][diagName]);
             ++diagBlockID;
         }
diff --git a/src/diagnostic/diagnostic_model_view.hpp b/src/diagnostic/diagnostic_model_view.hpp
index 8db43b3ff..f7a555e1a 100644
--- a/src/diagnostic/diagnostic_model_view.hpp
+++ b/src/diagnostic/diagnostic_model_view.hpp
@@ -1,6 +1,7 @@
 #ifndef DIAGNOSTIC_MODEL_VIEW_HPP
 #define DIAGNOSTIC_MODEL_VIEW_HPP
 
+#include "amr/amr_constants.hpp"
 #include "core/def.hpp"
 #include "core/utilities/mpi_utils.hpp"
 
@@ -33,8 +34,8 @@ class BaseModelView : public IModelView
     using GridLayout        = Model::gridlayout_type;
     using VecField          = Model::vecfield_type;
     using TensorFieldT      = Model::ions_type::tensorfield_type;
-    using GridLayoutT       = Model::gridlayout_type;
     using ResMan            = Model::resources_manager_type;
+    using Field             = Model::field_type;
     using TensorFieldData_t = ResMan::template UserTensorField_t</*rank=*/2>::patch_data_type;
     static constexpr auto dimension = Model::dimension;
 
@@ -67,41 +68,50 @@ class BaseModelView : public IModelView
         auto& rm   = *model_.resourcesManager;
         auto& ions = model_.state.ions;
 
-        for (auto patch : rm.enumerate(lvl, ions, sumTensor_))
+        for (auto patch : rm.enumerate(lvl, ions, tmpTensor_))
             for (std::uint8_t c = 0; c < N; ++c)
-                std::memcpy(sumTensor_[c].data(), ions[popidx].momentumTensor()[c].data(),
+                std::memcpy(tmpTensor_[c].data(), ions[popidx].momentumTensor()[c].data(),
                             ions[popidx].momentumTensor()[c].size() * sizeof(value_type));
 
         MTAlgos[popidx].getOrCreateSchedule(hierarchy_, lvl.getLevelNumber()).fillData(time);
 
-        for (auto patch : rm.enumerate(lvl, ions, sumTensor_))
+        for (auto patch : rm.enumerate(lvl, ions, tmpTensor_))
             for (std::uint8_t c = 0; c < N; ++c)
-                std::memcpy(ions[popidx].momentumTensor()[c].data(), sumTensor_[c].data(),
+                std::memcpy(ions[popidx].momentumTensor()[c].data(), tmpTensor_[c].data(),
                             ions[popidx].momentumTensor()[c].size() * sizeof(value_type));
     }
 
 
     template<typename Action>
-    void onLevels(Action&& action, int minlvl = 0, int maxlvl = 0)
+    void onLevels(Action&& action, std::size_t const minlvl = 0,
+                  std::size_t const maxlvl = amr::MAX_LEVEL_IDX)
     {
-        for (int ilvl = minlvl; ilvl < hierarchy_.getNumberOfLevels() && ilvl <= maxlvl; ++ilvl)
-            if (auto lvl = hierarchy_.getPatchLevel(ilvl))
-                action(*lvl);
+        amr::onLevels(hierarchy_, std::forward<Action>(action), minlvl, maxlvl);
+    }
+
+
+    template<typename OnLevel, typename OrMissing>
+    void onLevels(OnLevel&& onLevel, OrMissing&& orMissing, std::size_t const minlvl,
+                  std::size_t const maxlvl)
+    {
+        amr::onLevels(hierarchy_, std::forward<OnLevel>(onLevel),
+                      std::forward<OrMissing>(orMissing), minlvl, maxlvl);
     }
 
 
     template<typename Action>
     void visitHierarchy(Action&& action, int minLevel = 0, int maxLevel = 0)
     {
-        PHARE::amr::visitHierarchy<GridLayout>(hierarchy_, *model_.resourcesManager,
-                                               std::forward<Action>(action), minLevel, maxLevel,
-                                               model_);
+        amr::visitHierarchy<GridLayout>(hierarchy_, *model_.resourcesManager,
+                                        std::forward<Action>(action), minLevel, maxLevel, *this,
+                                        model_);
     }
 
     NO_DISCARD auto boundaryConditions() const { return hierarchy_.boundaryConditions(); }
     NO_DISCARD auto domainBox() const { return hierarchy_.domainBox(); }
     NO_DISCARD auto origin() const { return std::vector<double>(dimension, 0); }
     NO_DISCARD auto cellWidth() const { return hierarchy_.cellWidth(); }
+    NO_DISCARD auto maxLevel() const { return hierarchy_.maxLevel(); }
 
     NO_DISCARD std::string getLayoutTypeString() const
     {
@@ -141,6 +151,26 @@ class BaseModelView : public IModelView
         return model_.tags.at(key);
     }
 
+
+    auto& tmpField() { return tmpField_; }
+    auto& tmpVecField() { return tmpVec_; }
+
+    template<std::size_t rank = 2>
+    auto& tmpTensorField()
+    {
+        static_assert(rank > 0 and rank < 3);
+        if constexpr (rank == 1)
+            return tmpVec_;
+        else
+            return tmpTensor_;
+    }
+
+    NO_DISCARD auto getCompileTimeResourcesViewList()
+    {
+        return std::forward_as_tuple(tmpField_, tmpVec_, tmpTensor_);
+    }
+
+
 protected:
     Model& model_;
     Hierarchy& hierarchy_;
@@ -149,7 +179,7 @@ class BaseModelView : public IModelView
     {
         auto& rm = *model_.resourcesManager;
 
-        auto const dst_name = sumTensor_.name();
+        auto const dst_name = tmpTensor_.name();
 
         for (auto& pop : model_.state.ions)
         {
@@ -160,7 +190,7 @@ class BaseModelView : public IModelView
             MTAlgo.MTalgo->registerRefine(
                 idDst, idSrc, idDst, nullptr,
                 std::make_shared<
-                    amr::TensorFieldGhostInterpOverlapFillPattern<GridLayoutT, /*rank_=*/2>>());
+                    amr::TensorFieldGhostInterpOverlapFillPattern<GridLayout, /*rank_=*/2>>());
         }
 
         // can't create schedules here as the hierarchy has no levels yet
@@ -187,7 +217,9 @@ class BaseModelView : public IModelView
     };
 
     std::vector<MTAlgo> MTAlgos;
-    TensorFieldT sumTensor_{"PHARE_sumTensor", core::HybridQuantity::Tensor::M};
+    Field tmpField_{"PHARE_sumField", core::HybridQuantity::Scalar::rho};
+    VecField tmpVec_{"PHARE_sumVec", core::HybridQuantity::Vector::V};
+    TensorFieldT tmpTensor_{"PHARE_sumTensor", core::HybridQuantity::Tensor::M};
 };
 
 
diff --git a/src/diagnostic/diagnostic_writer.hpp b/src/diagnostic/diagnostic_writer.hpp
index a33c0a403..1d0b22344 100644
--- a/src/diagnostic/diagnostic_writer.hpp
+++ b/src/diagnostic/diagnostic_writer.hpp
@@ -3,7 +3,7 @@
 
 #include "diagnostic_props.hpp"
 
-#include "dict.hpp"
+
 
 namespace PHARE::diagnostic
 {
diff --git a/src/diagnostic/diagnostics.hpp b/src/diagnostic/diagnostics.hpp
index 74111ff37..041efc374 100644
--- a/src/diagnostic/diagnostics.hpp
+++ b/src/diagnostic/diagnostics.hpp
@@ -25,11 +25,7 @@
 #include "diagnostic_model_view.hpp"
 
 #include "diagnostic/detail/h5writer.hpp"
-#include "diagnostic/detail/types/electromag.hpp"
-#include "diagnostic/detail/types/particle.hpp"
-#include "diagnostic/detail/types/fluid.hpp"
-#include "diagnostic/detail/types/meta.hpp"
-#include "diagnostic/detail/types/info.hpp"
+#include "diagnostic/detail/vtkh5_writer.hpp"
 
 #endif
 
@@ -56,8 +52,13 @@ struct DiagnosticsManagerResolver
     {
 #if PHARE_HAS_HIGHFIVE
         using ModelView_t = ModelView<Hierarchy, Model>;
-        using Writer_t    = h5::H5Writer<ModelView_t>;
-        return DiagnosticsManager<Writer_t>::make_unique(hier, model, dict);
+        auto const format = cppdict::get_value(dict, "format", std::string{"phareh5"});
+
+        if (format == "phareh5")
+            return DiagnosticsManager<h5::H5Writer<ModelView_t>>::make_unique(hier, model, dict);
+        if (format == "pharevtkhdf")
+            return DiagnosticsManager<vtkh5::H5Writer<ModelView_t>>::make_unique(hier, model, dict);
+        throw std::runtime_error("DiagnosticsManagerResolver - unknown format " + format);
 #else
         return std::make_unique<NullOpDiagnosticsManager>();
 #endif
diff --git a/src/hdf5/detail/h5/h5_file.hpp b/src/hdf5/detail/h5/h5_file.hpp
index 510c0786b..f8208d75c 100644
--- a/src/hdf5/detail/h5/h5_file.hpp
+++ b/src/hdf5/detail/h5/h5_file.hpp
@@ -3,30 +3,29 @@
 
 #include "core/def.hpp"
 #include "core/def/phare_mpi.hpp" // IWYU pragma: keep
-#include "highfive/H5File.hpp"
-#include "highfive/H5Easy.hpp"
-
 #include "core/utilities/types.hpp"
 #include "core/utilities/mpi_utils.hpp"
 #include "core/utilities/meta/meta_utilities.hpp"
 
-namespace PHARE::hdf5::h5
+
+#include "highfive/H5File.hpp"
+#include "highfive/H5Easy.hpp"
+
+
+
+namespace PHARE::hdf5::h5::detail
 {
-using HiFile = HighFive::File;
-using FileOp = HighFive::File::AccessMode;
 
+// https://support.hdfgroup.org/documentation/hdf5/latest/hdf5_chunking.html
+static inline auto const CHUNK_SIZE = core::get_env_as("PHARE_H5_CHUNK_SIZE", std::size_t{1024});
 
+} // namespace PHARE::hdf5::h5::detail
 
-template<std::size_t dim, typename Data>
-NO_DISCARD auto decay_to_pointer(Data& data)
+namespace PHARE::hdf5::h5
 {
-    if constexpr (dim == 1)
-        return data.data();
-    if constexpr (dim == 2)
-        return data[0].data();
-    if constexpr (dim == 3)
-        return data[0][0].data();
-}
+using HiFile = HighFive::File;
+using FileOp = HighFive::File::AccessMode;
+
 
 template<typename Data, std::size_t dim>
 NO_DISCARD auto vector_for_dim()
@@ -107,11 +106,41 @@ class HighFiveFile
     template<typename Type, typename Size>
     auto create_data_set(std::string const& path, Size const& dataSetSize)
     {
+        if (exist(path))
+            return h5file_.getDataSet(path);
         createGroupsToDataSet(path);
         return h5file_.createDataSet<Type>(path, HighFive::DataSpace(dataSetSize));
     }
 
 
+    template<typename Type>
+    auto create_chunked_data_set(auto const& path, auto const chunk, auto const& dataspace)
+    {
+        if (exist(path))
+            return h5file_.getDataSet(path);
+        createGroupsToDataSet(path);
+        HighFive::DataSetCreateProps props;
+        props.add(HighFive::Chunking{chunk});
+        return h5file_.createDataSet(path, dataspace, HighFive::create_datatype<Type>(), props);
+    }
+
+    template<typename Type, std::size_t cols>
+    auto create_resizable_2d_data_set(auto const& path, hsize_t const x_chunk = detail::CHUNK_SIZE,
+                                      hsize_t const y_chunk = 1)
+    {
+        return create_chunked_data_set<Type>(
+            path, std::vector<hsize_t>{x_chunk, y_chunk},
+            HighFive::DataSpace({0, cols}, {HighFive::DataSpace::UNLIMITED, cols}));
+    }
+
+    template<typename Type>
+    auto create_resizable_1d_data_set(auto const& path, hsize_t const chunk = detail::CHUNK_SIZE)
+    {
+        return create_chunked_data_set<Type>(
+            path, std::vector<hsize_t>{chunk},
+            HighFive::DataSpace({0}, {HighFive::DataSpace::UNLIMITED}));
+    }
+
 
     /*
      * Communicate all dataset paths and sizes to all MPI process to allow each to create all
@@ -214,7 +243,7 @@ class HighFiveFile
         }();
         auto const paths = core::mpi::collect(path, mpi_size);
 
-        for (int i = 0; i < mpi_size; i++)
+        for (int i = 0; i < mpi_size; ++i)
         {
             std::string const keyPath = paths[i] == "null" ? "" : paths[i];
             if (keyPath.empty())
@@ -244,11 +273,18 @@ class HighFiveFile
         return attr;
     }
 
+    bool exist(std::string const& s) const { return h5file_.exist(s); }
+    auto getDataSet(std::string const& s)
+    {
+        if (!exist(s))
+            throw std::runtime_error("Dataset does not exist: " + s);
+        return h5file_.getDataSet(s);
+    }
 
-    HighFiveFile(HighFiveFile const&)             = delete;
-    HighFiveFile(HighFiveFile const&&)            = delete;
-    HighFiveFile& operator=(HighFiveFile const&)  = delete;
-    HighFiveFile& operator=(HighFiveFile const&&) = delete;
+    HighFiveFile(HighFiveFile const&)            = delete;
+    HighFiveFile(HighFiveFile&&)                 = delete;
+    HighFiveFile& operator=(HighFiveFile const&) = delete;
+    HighFiveFile& operator=(HighFiveFile&&)      = delete;
 
 private:
     HighFive::FileAccessProps fapl_;
diff --git a/src/phare/CMakeLists.txt b/src/phare/CMakeLists.txt
index f034996a0..2e56c3a08 100644
--- a/src/phare/CMakeLists.txt
+++ b/src/phare/CMakeLists.txt
@@ -6,7 +6,6 @@ if(PHARE_EXE)
 
   add_executable(phare-exe ${SOURCES_INC} phare.cpp)
   target_compile_options(phare-exe PRIVATE ${PHARE_WERROR_FLAGS} -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE})
-  configure_file(${CMAKE_CURRENT_SOURCE_DIR}/phare_init.py ${CMAKE_CURRENT_BINARY_DIR}/phare_init.py @ONLY)
   target_link_libraries(phare-exe PUBLIC
       ${PHARE_BASE_LIBS}
       phare_simulator
diff --git a/src/phare/phare_init.py b/src/phare/phare_init.py
deleted file mode 100644
index d79de42b7..000000000
--- a/src/phare/phare_init.py
+++ /dev/null
@@ -1,109 +0,0 @@
-#!/usr/bin/env python3
-
-import numpy as np
-
-import pyphare.pharein as ph
-
-# configure the simulation
-
-ph.Simulation(
-    smallest_patch_size=20,
-    largest_patch_size=20,
-    time_step_nbr=30000,  # number of time steps (not specified if time_step and final_time provided)
-    final_time=30.0,  # simulation final time (not specified if time_step and time_step_nbr provided)
-    boundary_types="periodic",  # boundary condition, string or tuple, length == len(cell) == len(dl)
-    cells=40,  # integer or tuple length == dimension
-    dl=0.3,  # mesh size of the root level, float or tuple
-    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"}},
-)
-
-
-def density(x):
-    return 1.0
-
-
-def by(x):
-    from pyphare.pharein.global_vars import sim
-
-    L = sim.simulation_domain()
-    return 0.1 * np.cos(2 * np.pi * x / L[0])
-
-
-def bz(x):
-    from pyphare.pharein.global_vars import sim
-
-    L = sim.simulation_domain()
-    return 0.1 * np.sin(2 * np.pi * x / L[0])
-
-
-def bx(x):
-    return 1.0
-
-
-def vx(x):
-    return 0.0
-
-
-def vy(x):
-    from pyphare.pharein.global_vars import sim
-
-    L = sim.simulation_domain()
-    return 0.1 * np.cos(2 * np.pi * x / L[0])
-
-
-def vz(x):
-    from pyphare.pharein.global_vars import sim
-
-    L = sim.simulation_domain()
-    return 0.1 * np.sin(2 * np.pi * x / L[0])
-
-
-def vthx(x):
-    return 0.01
-
-
-def vthy(x):
-    return 0.01
-
-
-def vthz(x):
-    return 0.01
-
-
-vvv = {
-    "vbulkx": vx,
-    "vbulky": vy,
-    "vbulkz": vz,
-    "vthx": vthx,
-    "vthy": vthy,
-    "vthz": vthz,
-}
-
-ph.MaxwellianFluidModel(
-    bx=bx,
-    by=by,
-    bz=bz,
-    protons={"charge": 1, "density": density, **vvv, "init": {"seed": 1337}},
-)
-
-ph.ElectronModel(closure="isothermal", Te=0.12)
-
-
-sim = ph.global_vars.sim
-
-timestamps = np.arange(0, sim.final_time + sim.time_step, 100 * sim.time_step)
-
-
-for quantity in ["E", "B"]:
-    ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
-
-
-for quantity in ["charge_density", "bulkVelocity"]:
-    ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
-
-pops = [
-    "protons",
-]
diff --git a/src/phare/phare_init_small.py b/src/phare/phare_init_small.py
deleted file mode 100644
index f3f8e76e1..000000000
--- a/src/phare/phare_init_small.py
+++ /dev/null
@@ -1,112 +0,0 @@
-#!/usr/bin/env python3
-
-import numpy as np
-
-import pyphare.pharein as ph
-
-# configure the simulation
-
-ph.Simulation(
-    smallest_patch_size=20,
-    largest_patch_size=20,
-    time_step_nbr=300,  # number of time steps (not specified if time_step and final_time provided)
-    time_step=0.001,  # simulation final time (not specified if time_step and time_step_nbr provided)
-    boundary_types="periodic",  # boundary condition, string or tuple, length == len(cell) == len(dl)
-    cells=40,  # integer or tuple length == dimension
-    dl=0.3,  # mesh size of the root level, float or tuple
-    # 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", "mode": "overwrite"},
-    },
-)
-
-
-def density(x):
-    return 1.0
-
-
-def by(x):
-    from pyphare.pharein.global_vars import sim
-
-    L = sim.simulation_domain()
-    return 0.1 * np.cos(2 * np.pi * x / L[0])
-
-
-def bz(x):
-    from pyphare.pharein.global_vars import sim
-
-    L = sim.simulation_domain()
-    return 0.1 * np.sin(2 * np.pi * x / L[0])
-
-
-def bx(x):
-    return 1.0
-
-
-def vx(x):
-    return 0.0
-
-
-def vy(x):
-    from pyphare.pharein.global_vars import sim
-
-    L = sim.simulation_domain()
-    return 0.1 * np.cos(2 * np.pi * x / L[0])
-
-
-def vz(x):
-    from pyphare.pharein.global_vars import sim
-
-    L = sim.simulation_domain()
-    return 0.1 * np.sin(2 * np.pi * x / L[0])
-
-
-def vthx(x):
-    return 0.01
-
-
-def vthy(x):
-    return 0.01
-
-
-def vthz(x):
-    return 0.01
-
-
-vvv = {
-    "vbulkx": vx,
-    "vbulky": vy,
-    "vbulkz": vz,
-    "vthx": vthx,
-    "vthy": vthy,
-    "vthz": vthz,
-}
-
-ph.MaxwellianFluidModel(
-    bx=bx,
-    by=by,
-    bz=bz,
-    protons={"charge": 1, "density": density, **vvv, "init": {"seed": 1337}},
-)
-
-ph.ElectronModel(closure="isothermal", Te=0.12)
-
-
-sim = ph.global_vars.sim
-
-timestamps = np.arange(0, sim.final_time + sim.time_step, 10 * sim.time_step)
-
-
-for quantity in ["E", "B"]:
-    ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
-
-
-for quantity in ["charge_density", "bulkVelocity"]:
-    ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
-
-pops = [
-    "protons",
-]
diff --git a/src/simulator/simulator.hpp b/src/simulator/simulator.hpp
index ae92bc2c5..5883d4912 100644
--- a/src/simulator/simulator.hpp
+++ b/src/simulator/simulator.hpp
@@ -359,7 +359,6 @@ Simulator<opts>::Simulator(PHARE::initializer::PHAREDict const& dict,
 {
     resman_ptr   = std::make_shared<ResourceManager_t>();
     currentTime_ = restart_time(dict);
-    finalTime_ += currentTime_;
 
     if (dict["simulation"].contains("restarts"))
         rMan = restarts::RestartsManagerResolver::make_unique(*hierarchy_, *resman_ptr,
diff --git a/tests/simulator/CMakeLists.txt b/tests/simulator/CMakeLists.txt
index dde92f091..23da715ac 100644
--- a/tests/simulator/CMakeLists.txt
+++ b/tests/simulator/CMakeLists.txt
@@ -18,6 +18,7 @@ if(HighFive)
   phare_python3_exec(9       diagnostics test_diagnostics.py  ${CMAKE_CURRENT_BINARY_DIR}) # serial or n = 2
   phare_python3_exec(9       restarts    test_restarts.py     ${CMAKE_CURRENT_BINARY_DIR}) # serial or n = 2
   phare_python3_exec(9       run         test_run.py          ${CMAKE_CURRENT_BINARY_DIR}) # serial or n = 2
+  phare_python3_exec(9       vtk_diagnostics test_vtk_diagnostics.py  ${CMAKE_CURRENT_BINARY_DIR}) # serial or n = 2
   # phare_python3_exec(9       tagging     test_tagging.py      ${CMAKE_CURRENT_BINARY_DIR}) # serial or n = 2
 
   if(testMPI)
@@ -27,10 +28,11 @@ if(HighFive)
     # doesn't make sense in serial
     phare_mpi_python3_exec(9 3 load_balancing test_load_balancing.py ${CMAKE_CURRENT_BINARY_DIR})
 
-  endif(testMPI)
+    phare_mpi_python3_exec(9 4 vtk_diagnostics test_vtk_diagnostics.py ${CMAKE_CURRENT_BINARY_DIR})
 
+  endif(testMPI)
 
-  # phare_python3_exec(11 test_diagnostic_timestamps test_diagnostic_timestamps.py ${CMAKE_CURRENT_BINARY_DIR})
+  phare_python3_exec(11 test_diagnostic_timestamps test_diagnostic_timestamps.py ${CMAKE_CURRENT_BINARY_DIR})
 
 endif()
 
diff --git a/tests/simulator/__init__.py b/tests/simulator/__init__.py
index 47aa1f8b7..dc4986426 100644
--- a/tests/simulator/__init__.py
+++ b/tests/simulator/__init__.py
@@ -250,7 +250,10 @@ def run(self, result=None):
     def unique_diag_dir_for_test_case(self, base_path, ndim, interp, post_path=""):
         from pyphare import cpp
 
-        return f"{base_path}/{self._testMethodName}/{cpp.mpi_size()}/{ndim}/{interp}/{post_path}"
+        base = f"{base_path}/{self._testMethodName}/{cpp.mpi_size()}/{ndim}/{interp}"
+        if post_path:
+            return base + "/" + post_path
+        return base
 
     def clean_up_diags_dirs(self):
         from pyphare import cpp
diff --git a/tests/simulator/test_restarts.py b/tests/simulator/test_restarts.py
index 8d36aabfa..fcba1fe8c 100644
--- a/tests/simulator/test_restarts.py
+++ b/tests/simulator/test_restarts.py
@@ -2,10 +2,10 @@
 #
 
 import copy
-from time import sleep
 import datetime
 import unittest
 import numpy as np
+from time import sleep
 from pathlib import Path
 
 from datetime import timedelta
@@ -428,7 +428,7 @@ def test_advanced_restarts_options(self):
         Dim / interp / etc are not relevant here
         """
         ndim, interp = 1, 1
-        print(f"test_advanced_restarts_options")
+        print("test_advanced_restarts_options")
 
         simput = copy.deepcopy(
             dup(
@@ -453,10 +453,9 @@ def test_advanced_restarts_options(self):
 
         ph.global_vars.sim = None
         ph.global_vars.sim = ph.Simulation(**simput)
-        model = setup_model()
+        setup_model()
         Simulator(ph.global_vars.sim).run().reset()
         self.register_diag_dir_for_cleanup(local_out)
-        diag_dir0 = local_out
 
         simput["restart_options"]["restart_time"] = "auto"
         self.assertEqual(0.01, ph.restarts.restart_time(simput["restart_options"]))
diff --git a/tests/simulator/test_vtk_diagnostics.py b/tests/simulator/test_vtk_diagnostics.py
new file mode 100644
index 000000000..e2348759b
--- /dev/null
+++ b/tests/simulator/test_vtk_diagnostics.py
@@ -0,0 +1,209 @@
+#
+
+import unittest
+import itertools
+import numpy as np
+from copy import deepcopy
+from ddt import data, ddt, unpack
+
+from pyphare import cpp
+import pyphare.pharein as ph
+from pyphare.pharesee.run import Run
+
+from pyphare.core import phare_utilities as phut
+from pyphare.simulator.simulator import startMPI
+from pyphare.simulator.simulator import Simulator
+from pyphare.pharesee.hierarchy import hierarchy_utils as hootils
+
+from tests.simulator import SimulatorTest
+from tests.diagnostic import dump_all_diags
+
+# only in 2d for now
+ppc_per_dim = [100, 25, 10]
+
+
+def config(ndim, interp, **simInput):
+    ppc = ppc_per_dim[ndim - 1]
+    sim = ph.Simulation(**simInput)
+
+    L = 0.5
+
+    def density(x, y):
+        Ly = sim.simulation_domain()[1]
+        return (
+            0.4
+            + 1.0 / np.cosh((y - Ly * 0.3) / L) ** 2
+            + 1.0 / np.cosh((y - Ly * 0.7) / L) ** 2
+        )
+
+    def S(y, y0, l):
+        return 0.5 * (1.0 + np.tanh((y - y0) / l))
+
+    def by(x, y):
+        Lx = sim.simulation_domain()[0]
+        Ly = sim.simulation_domain()[1]
+        sigma = 1.0
+        dB = 0.1
+
+        x0 = x - 0.5 * Lx
+        y1 = y - 0.3 * Ly
+        y2 = y - 0.7 * Ly
+
+        dBy1 = 2 * dB * x0 * np.exp(-(x0**2 + y1**2) / (sigma) ** 2)
+        dBy2 = -2 * dB * x0 * np.exp(-(x0**2 + y2**2) / (sigma) ** 2)
+
+        return dBy1 + dBy2
+
+    def bx(x, y):
+        Lx = sim.simulation_domain()[0]
+        Ly = sim.simulation_domain()[1]
+        sigma = 1.0
+        dB = 0.1
+
+        x0 = x - 0.5 * Lx
+        y1 = y - 0.3 * Ly
+        y2 = y - 0.7 * Ly
+
+        dBx1 = -2 * dB * y1 * np.exp(-(x0**2 + y1**2) / (sigma) ** 2)
+        dBx2 = 2 * dB * y2 * np.exp(-(x0**2 + y2**2) / (sigma) ** 2)
+
+        v1 = -1
+        v2 = 1.0
+        return v1 + (v2 - v1) * (S(y, Ly * 0.3, L) - S(y, Ly * 0.7, L)) + dBx1 + dBx2
+
+    def bz(x, y):
+        return 0.0
+
+    def b2(x, y):
+        return bx(x, y) ** 2 + by(x, y) ** 2 + bz(x, y) ** 2
+
+    def T(x, y):
+        K = 0.7
+        temp = 1.0 / density(x, y) * (K - b2(x, y) * 0.5)
+        assert np.all(temp > 0)
+        return temp
+
+    def vx(x, y):
+        return 0.0
+
+    def vy(x, y):
+        return 0.0
+
+    def vz(x, y):
+        return 0.0
+
+    def vthx(x, y):
+        return np.sqrt(T(x, y))
+
+    def vthy(x, y):
+        return np.sqrt(T(x, y))
+
+    def vthz(x, y):
+        return np.sqrt(T(x, y))
+
+    vvv = {
+        "vbulkx": vx,
+        "vbulky": vy,
+        "vbulkz": vz,
+        "vthx": vthx,
+        "vthy": vthy,
+        "vthz": vthz,
+    }
+
+    ph.MaxwellianFluidModel(
+        bx=bx,
+        by=by,
+        bz=bz,
+        protons={
+            "mass": 1,
+            "charge": 1,
+            "density": density,
+            **vvv,
+            "nbr_part_per_cell": ppc,
+            "init": {"seed": 1337},
+        },
+        alpha={
+            "mass": 4,
+            "charge": 1,
+            "density": density,
+            **vvv,
+            "nbr_part_per_cell": ppc,
+            "init": {"seed": 2334},
+        },
+    )
+    ph.ElectronModel(closure="isothermal", Te=0.12)
+    return sim
+
+
+out = "phare_outputs/vtk_diagnostic_test"
+simArgs = {
+    "time_step_nbr": 1,
+    "final_time": 0.001,
+    "boundary_types": "periodic",
+    "cells": 40,
+    "dl": 0.3,
+    "diag_options": {
+        "format": "pharevtkhdf",
+        "options": {"dir": out, "mode": "overwrite"},
+    },
+}
+
+
+def permute(dic):
+    ndims = [2]
+    interp_orders = [1]
+    dic.update(simArgs.copy())
+    return [
+        dict(
+            ndim=ndim,
+            interp=interp_order,
+            simInput=deepcopy(dic),
+        )
+        for ndim, interp_order in itertools.product(ndims, interp_orders)
+    ]
+
+
+@ddt
+class VTKDiagnosticsTest(SimulatorTest):
+    def __init__(self, *args, **kwargs):
+        super(VTKDiagnosticsTest, self).__init__(*args, **kwargs)
+        self.simulator = None
+        ph.global_vars.sim = None
+
+    def _run(self, ndim, interp, simInput, diag_dir="", **kwargs):
+        for key in ["cells", "dl", "boundary_types"]:
+            simInput[key] = list(phut.np_array_ify(simInput[key], ndim))
+        local_out = self.unique_diag_dir_for_test_case(
+            f"{out}{'/'+diag_dir if diag_dir else ''}", ndim, interp
+        )
+        self.register_diag_dir_for_cleanup(local_out)
+        simInput["diag_options"]["options"]["dir"] = local_out
+        simulation = config(ndim, interp, **simInput)
+        self.assertTrue(len(simulation.cells) == ndim)
+        dump_all_diags(simulation.model.populations)
+        Simulator(simulation).run().reset()
+        ph.global_vars.sim = None
+        return local_out
+
+    @data(*permute({}))
+    @unpack
+    def test_dump_diags(self, ndim, interp, simInput):
+        print("test_dump_diags dim/interp:{}/{}".format(ndim, interp))
+
+        b0 = [[10 for i in range(ndim)], [19 for i in range(ndim)]]
+        simInput["refinement_boxes"] = {"L0": {"B0": b0}}
+        local_out = self._run(ndim, interp, simInput)
+
+        if cpp.mpi_rank() == 0:
+            try:
+                from pyphare.pharesee.phare_vtk import plot as plot_vtk
+
+                plot_vtk(local_out + "/EM_B.vtkhdf", f"B{ndim}d_interp{interp}.vtk.png")
+                plot_vtk(local_out + "/EM_E.vtkhdf", f"E{ndim}d_interp{interp}.vtk.png")
+            except ModuleNotFoundError:
+                print("WARNING: vtk python module not found - cannot make plots")
+
+
+if __name__ == "__main__":
+    startMPI()
+    unittest.main()