diff --git a/fluidsim/base/output/phys_fields.py b/fluidsim/base/output/phys_fields.py
index c6bb9b9b3aeafd19ad71c511a2a701e42847673a_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMucHk=..8afa4454ac9fbe9d917927230a1b34d94ba9d943_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMucHk= 100644
--- a/fluidsim/base/output/phys_fields.py
+++ b/fluidsim/base/output/phys_fields.py
@@ -205,5 +205,6 @@
                     _create_variable(group_state_phys, k, field_seq)
         else:
             h5file.atomic = False
-            if len(self.oper.shapeX_loc) == 2:
+            ndim = len(self.oper.shapeX_loc)
+            if ndim == 2:
                 xstart, ystart = self.oper.seq_indices_first_X
@@ -209,3 +210,3 @@
                 xstart, ystart = self.oper.seq_indices_first_X
-            if len(self.oper.shapeX_loc) == 3:
+            elif ndim == 3:
                 xstart, ystart, zstart = self.oper.seq_indices_first_X
@@ -211,4 +212,6 @@
                 xstart, ystart, zstart = self.oper.seq_indices_first_X
+            else:
+                raise NotImplementedError
             xend = xstart + self.oper.shapeX_loc[0]
             yend = ystart + self.oper.shapeX_loc[1]
             for k in state_phys.keys:
diff --git a/fluidsim/base/test/test_base_solver.py b/fluidsim/base/test/test_base_solver.py
index c6bb9b9b3aeafd19ad71c511a2a701e42847673a_Zmx1aWRzaW0vYmFzZS90ZXN0L3Rlc3RfYmFzZV9zb2x2ZXIucHk=..8afa4454ac9fbe9d917927230a1b34d94ba9d943_Zmx1aWRzaW0vYmFzZS90ZXN0L3Rlc3RfYmFzZV9zb2x2ZXIucHk= 100644
--- a/fluidsim/base/test/test_base_solver.py
+++ b/fluidsim/base/test/test_base_solver.py
@@ -2,8 +2,7 @@
 import shutil
 
 import fluiddyn as fld
-from fluiddyn.io import stdout_redirected
-import fluiddyn.util.mpi as mpi
+from fluiddyn.util import mpi
 
 # to get fld.show
 import fluiddyn.output
@@ -21,8 +20,7 @@
         params.time_stepping.it_end = 4
         params.time_stepping.deltat0 = 0.1
 
-        with stdout_redirected():
-            self.sim = Simul(params)
+        self.sim = Simul(params)
 
     def tearDown(self):
         # clean by removing the directory
@@ -30,5 +28,6 @@
             if hasattr(self, "sim"):
                 shutil.rmtree(self.sim.output.path_run)
 
+    @unittest.skipIf(mpi.nb_proc > 1, "Od solvers do not work with mpi")
     def test_simul(self):
         """Should be able to run a base experiment."""
@@ -33,7 +32,6 @@
     def test_simul(self):
         """Should be able to run a base experiment."""
-        with stdout_redirected():
-            self.sim.time_stepping.start()
+        self.sim.time_stepping.start()
 
         fld.show()
 
diff --git a/fluidsim/operators/operators0d.py b/fluidsim/operators/operators0d.py
index c6bb9b9b3aeafd19ad71c511a2a701e42847673a_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczBkLnB5..8afa4454ac9fbe9d917927230a1b34d94ba9d943_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczBkLnB5 100644
--- a/fluidsim/operators/operators0d.py
+++ b/fluidsim/operators/operators0d.py
@@ -9,9 +9,7 @@
 
 """
 
-from __future__ import division
-
-from builtins import object
+from fluiddyn.util import mpi
 
 
 class Operators0D(object):
@@ -26,6 +24,8 @@
         params._set_child("oper")
 
     def __init__(self, params=None, SEQUENTIAL=None):
+        if mpi.nb_proc > 1:
+            raise ValueError
 
         self.params = params
         self.axes = tuple()