diff --git a/.travis.yml b/.travis.yml
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_LnRyYXZpcy55bWw=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_LnRyYXZpcy55bWw= 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -50,7 +50,7 @@
 install:
     - pip install -U pip setuptools wheel six colorlog
     - pip install -U pipenv tox-pipenv coverage fluiddevops
-    - pip install -U numpy cython mako mpi4py "fluidpythran>=0.1.4"  # To be removed when pip==10 is out and this step can be bootstrapped
+    - pip install -U numpy cython mako mpi4py fluidpythran  # To be removed when pip==10 is out and this step can be bootstrapped
 
 before_script:
     - |
diff --git a/Makefile b/Makefile
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_TWFrZWZpbGU=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_TWFrZWZpbGU= 100644
--- a/Makefile
+++ b/Makefile
@@ -32,7 +32,8 @@
 _tests_coverage:
 	mkdir -p .coverage
 	coverage run -p -m fluidsim.util.testing -v
-	mpirun -np 2 coverage run -p -m fluidsim.util.testing -v
+	FLUIDPYTHRAN_NO_REPLACE=1 coverage run -p -m fluidsim.util.testing -v
+	FLUIDPYTHRAN_NO_REPLACE=1 mpirun -np 2 coverage run -p -m fluidsim.util.testing -v
 
 _report_coverage:
 	coverage combine
@@ -44,6 +45,11 @@
 
 coverage: _tests_coverage _report_coverage
 
+coverage_short:
+	mkdir -p .coverage
+	FLUIDPYTHRAN_NO_REPLACE=1 coverage run -p -m fluidsim.util.testing -v
+	make _report_coverage
+
 lint:
 	pylint -rn --rcfile=pylintrc --jobs=$(shell nproc) fluidsim --exit-zero
 
diff --git a/bench/spectraldns/bench_3d.py b/bench/spectraldns/bench_3d.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_YmVuY2gvc3BlY3RyYWxkbnMvYmVuY2hfM2QucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_YmVuY2gvc3BlY3RyYWxkbnMvYmVuY2hfM2QucHk= 100644
--- a/bench/spectraldns/bench_3d.py
+++ b/bench/spectraldns/bench_3d.py
@@ -16,7 +16,6 @@
   mpirun -np 2 fluidsim-bench 128 -d 3 -s ns3d -it 10
 
 """
-from __future__ import print_function
 import time
 import numpy as np
 from numpy import zeros, sum, float64, sin, cos, prod, asscalar
diff --git a/bench/spectraldns/profile_3d.py b/bench/spectraldns/profile_3d.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_YmVuY2gvc3BlY3RyYWxkbnMvcHJvZmlsZV8zZC5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_YmVuY2gvc3BlY3RyYWxkbnMvcHJvZmlsZV8zZC5weQ== 100644
--- a/bench/spectraldns/profile_3d.py
+++ b/bench/spectraldns/profile_3d.py
@@ -13,7 +13,6 @@
   fluidsim-bench 128 -d 3 -s ns3d -it 10
 
 """
-from __future__ import print_function
 import time
 import numpy as np
 from numpy import zeros, sum, float64, sin, cos, prod, asscalar
diff --git a/bitbucket-pipelines.yml b/bitbucket-pipelines.yml
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Yml0YnVja2V0LXBpcGVsaW5lcy55bWw=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Yml0YnVja2V0LXBpcGVsaW5lcy55bWw= 100644
--- a/bitbucket-pipelines.yml
+++ b/bitbucket-pipelines.yml
@@ -14,5 +14,5 @@
             caches:
               - pip
             script:
-              - pip install -U tox "fluidpythran>=0.1.4" --user
+              - pip install -U tox fluidpythran --user
               - tox -e py36,codecov
diff --git a/doc/example_from_numpy_doc.py b/doc/example_from_numpy_doc.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_ZG9jL2V4YW1wbGVfZnJvbV9udW1weV9kb2MucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_ZG9jL2V4YW1wbGVfZnJvbV9udW1weV9kb2MucHk= 100644
--- a/doc/example_from_numpy_doc.py
+++ b/doc/example_from_numpy_doc.py
@@ -8,7 +8,6 @@
 a line by itself, preferably preceeded by a blank line.
 
 """
-from __future__ import division, absolute_import, print_function
 
 import os # standard library imports first
 
@@ -121,5 +120,3 @@
     """
 
     pass
-
-
diff --git a/doc/examples/simul_ns3dbouss_initfields_in_script.py b/doc/examples/simul_ns3dbouss_initfields_in_script.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_ZG9jL2V4YW1wbGVzL3NpbXVsX25zM2Rib3Vzc19pbml0ZmllbGRzX2luX3NjcmlwdC5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_ZG9jL2V4YW1wbGVzL3NpbXVsX25zM2Rib3Vzc19pbml0ZmllbGRzX2luX3NjcmlwdC5weQ== 100644
--- a/doc/examples/simul_ns3dbouss_initfields_in_script.py
+++ b/doc/examples/simul_ns3dbouss_initfields_in_script.py
@@ -16,5 +16,5 @@
 
 params = Simul.create_default_params()
 
-params.output.sub_directory = 'examples'
+params.output.sub_directory = "examples"
 
@@ -20,9 +20,9 @@
 
-nx = 48//2
-ny = 64//2
-nz = 96//2
+nx = 48 // 2
+ny = 64 // 2
+nz = 96 // 2
 Lx = 3
 params.oper.nx = nx
 params.oper.ny = ny
 params.oper.nz = nz
 params.oper.Lx = Lx
@@ -24,13 +24,13 @@
 Lx = 3
 params.oper.nx = nx
 params.oper.ny = ny
 params.oper.nz = nz
 params.oper.Lx = Lx
-params.oper.Ly = Ly = Lx/nx*ny
-params.oper.Lz = Lz = Lx/nx*nz
-fft = 'fftwmpi3d'
+params.oper.Ly = Ly = Lx / nx * ny
+params.oper.Lz = Lz = Lx / nx * nz
+fft = "fftwmpi3d"
 # fft = 'fftw1d'
 # fft = 'pfft'
 # fft = 'p3dfft'
 
 # for sequential runs, just comment these 2 lines
@@ -32,9 +32,9 @@
 # fft = 'fftw1d'
 # fft = 'pfft'
 # fft = 'p3dfft'
 
 # for sequential runs, just comment these 2 lines
-params.oper.type_fft = 'fluidfft.fft3d.mpi_with_' + fft
+params.oper.type_fft = "fluidfft.fft3d.mpi_with_" + fft
 params.short_name_type_run = fft
 
 
@@ -62,7 +62,7 @@
 
 """
 n = 8
-C = 1.
-dx = Lx/nx
+C = 1.0
+dx = Lx / nx
 B = 1
 D = 1
@@ -67,5 +67,5 @@
 B = 1
 D = 1
-eps = 1e-2*B**(3/2)*D**(1/2)
-params.nu_8 = (dx/C)**((3*n-2)/3) * eps**(1/3)
+eps = 1e-2 * B ** (3 / 2) * D ** (1 / 2)
+params.nu_8 = (dx / C) ** ((3 * n - 2) / 3) * eps ** (1 / 3)
 
@@ -71,4 +71,4 @@
 
-printby0(f'nu_8 = {params.nu_8:.3e}')
+printby0(f"nu_8 = {params.nu_8:.3e}")
 
 params.time_stepping.USE_T_END = True
@@ -73,4 +73,4 @@
 
 params.time_stepping.USE_T_END = True
-params.time_stepping.t_end = 10.
+params.time_stepping.t_end = 10.0
 
@@ -76,5 +76,5 @@
 
-params.init_fields.type = 'in_script'
+params.init_fields.type = "in_script"
 
 params.output.periods_print.print_stdout = 1e-1
 
@@ -86,8 +86,9 @@
 
 # here we have to initialize the flow fields
 
-variables = {k: 1e-2 * sim.oper.create_arrayX_random()
-             for k in ('vx', 'vy', 'vz')}
+variables = {
+    k: 1e-2 * sim.oper.create_arrayX_random() for k in ("vx", "vy", "vz")
+}
 
 X, Y, Z = sim.oper.get_XYZ_loc()
 
@@ -91,8 +92,8 @@
 
 X, Y, Z = sim.oper.get_XYZ_loc()
 
-x0 = Lx/2.
-y0 = Ly/2.
-z0 = Lz/2.
-R2 = (X-x0)**2 + (Y-y0)**2 + (Z-z0)**2
+x0 = Lx / 2.0
+y0 = Ly / 2.0
+z0 = Lz / 2.0
+R2 = (X - x0) ** 2 + (Y - y0) ** 2 + (Z - z0) ** 2
 r0 = 0.5
@@ -98,6 +99,6 @@
 r0 = 0.5
-b = -0.5*(1-np.tanh((R2 - r0**2)/0.2**2))
-variables['b'] = b
+b = -0.5 * (1 - np.tanh((R2 - r0 ** 2) / 0.2 ** 2))
+variables["b"] = b
 
 sim.state.init_statephys_from(**variables)
 
@@ -106,8 +107,9 @@
 
 sim.time_stepping.start()
 
-printby0(f"""
+printby0(
+    f"""
 To visualize the output with Paraview, create a file states_phys.xmf with:
 
 fluidsim-create-xml-description {sim.output.path_run}
 
@@ -110,5 +112,6 @@
 To visualize the output with Paraview, create a file states_phys.xmf with:
 
 fluidsim-create-xml-description {sim.output.path_run}
 
-""")
+"""
+)
diff --git a/doc/test_cases/Taylor_Green_vortices/run_simul.py b/doc/test_cases/Taylor_Green_vortices/run_simul.py
new file mode 100644
index 0000000000000000000000000000000000000000..e7e3275ab7675c5fc429ecdfa056b44c05b78136_ZG9jL3Rlc3RfY2FzZXMvVGF5bG9yX0dyZWVuX3ZvcnRpY2VzL3J1bl9zaW11bC5weQ==
--- /dev/null
+++ b/doc/test_cases/Taylor_Green_vortices/run_simul.py
@@ -0,0 +1,79 @@
+"""Taylor-Green Vortex at Re = 1600
+===================================
+
+https://www.grc.nasa.gov/hiocfd/wp-content/uploads/sites/22/case_c3.3.pdf
+
+
+"""
+
+import numpy as np
+
+from fluidsim.solvers.ns3d.solver import Simul
+
+params = Simul.create_default_params()
+
+Re = 1600
+V0 = 1.0
+L = 1
+params.nu_2 = V0 * L / Re
+
+params.init_fields.type = "in_script"
+
+params.time_stepping.t_end = 20.0 * L / V0
+
+nx = 128 // 4
+params.oper.nx = params.oper.ny = params.oper.nz = nx
+lx = params.oper.Lx = params.oper.Ly = params.oper.Lz = 2 * np.pi * L
+
+params.output.periods_save.phys_fields = 4
+params.output.periods_save.spatial_means = 0.2
+params.output.periods_save.spectra = 0.5
+
+sim = Simul(params)
+
+X, Y, Z = sim.oper.get_XYZ_loc()
+
+vx = V0 * np.sin(X / L) * np.cos(Y / L) * np.cos(Z / L)
+vy = -V0 * np.cos(X / L) * np.sin(Y / L) * np.cos(Z / L)
+vz = sim.oper.create_arrayX(value=0)
+
+sim.state.init_statephys_from(vx=vx, vy=vy, vz=vz)
+
+sim.state.statespect_from_statephys()
+sim.state.statephys_from_statespect()
+
+# only useful to plot fields before time_stepping
+# sim.output.init_with_initialized_state()
+
+if __name__ == "__main__":
+
+    sim.time_stepping.start()
+
+    from fluiddyn.util.mpi import printby0
+
+    printby0(
+        f"""
+To visualize the output with Paraview, create a file states_phys.xmf with:
+
+fluidsim-create-xml-description {sim.output.path_run}
+
+# To visualize with fluidsim:
+
+cd {sim.output.path_run}
+ipython
+
+# in ipython:
+
+from fluidsim import load_sim_for_plot
+sim = load_sim_for_plot()
+
+sim.output.phys_fields.set_equation_crosssection(f'x={{sim.oper.Lx/4}}')
+sim.output.phys_fields.animate('vx')
+
+sim.output.phys_fields.plot(field="vx", time=10)
+
+sim.output.spatial_means.plot()
+sim.output.spectra.plot1d(tmin=12, tmax=16, coef_compensate=5/3)
+
+"""
+    )
diff --git a/fluidsim/__init__.py b/fluidsim/__init__.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vX19pbml0X18ucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vX19pbml0X18ucHk= 100644
--- a/fluidsim/__init__.py
+++ b/fluidsim/__init__.py
@@ -33,5 +33,7 @@
         Path(os.environ["FLUIDSIM_PATH"]) / ".fluidpythran"
     )
 
+_is_testing = False
+
 if any(
     any(test_tool in arg for arg in sys.argv)
@@ -36,4 +38,4 @@
 if any(
     any(test_tool in arg for arg in sys.argv)
-    for test_tool in ("pytest", "unittest", "fluidsim-test")
+    for test_tool in ("pytest", "unittest", "fluidsim-test", "coverage")
 ):
@@ -39,4 +41,5 @@
 ):
+    _is_testing = True
     from fluiddyn.util import mpi
 
     mpi.printby0(
diff --git a/fluidsim/base/basilisk/operators.py b/fluidsim/base/basilisk/operators.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay9vcGVyYXRvcnMucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay9vcGVyYXRvcnMucHk= 100644
--- a/fluidsim/base/basilisk/operators.py
+++ b/fluidsim/base/basilisk/operators.py
@@ -8,9 +8,6 @@
    :private-members:
 
 """
-from __future__ import division
-
-from builtins import object
 
 import numpy as np
 
diff --git a/fluidsim/base/forcing/anisotropic.py b/fluidsim/base/forcing/anisotropic.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9mb3JjaW5nL2FuaXNvdHJvcGljLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9mb3JjaW5nL2FuaXNvdHJvcGljLnB5 100644
--- a/fluidsim/base/forcing/anisotropic.py
+++ b/fluidsim/base/forcing/anisotropic.py
@@ -7,8 +7,6 @@
    :private-members:
 
 """
-from __future__ import division
-from __future__ import print_function
 
 from math import degrees
 
diff --git a/fluidsim/base/forcing/base.py b/fluidsim/base/forcing/base.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9mb3JjaW5nL2Jhc2UucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9mb3JjaW5nL2Jhc2UucHk= 100644
--- a/fluidsim/base/forcing/base.py
+++ b/fluidsim/base/forcing/base.py
@@ -118,7 +118,10 @@
         ClassForcing = dict_classes[self.type_forcing]
 
         self.sim = sim
-        self.forcing_maker = ClassForcing(sim)
+        if not sim.params.ONLY_COARSE_OPER:
+            self.forcing_maker = ClassForcing(sim)
+        else:
+            self.forcing_maker = None
 
         self._t_last_computed = -np.inf
 
diff --git a/fluidsim/base/forcing/specific.py b/fluidsim/base/forcing/specific.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9mb3JjaW5nL3NwZWNpZmljLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9mb3JjaW5nL3NwZWNpZmljLnB5 100644
--- a/fluidsim/base/forcing/specific.py
+++ b/fluidsim/base/forcing/specific.py
@@ -41,10 +41,6 @@
    :private-members:
 
 """
-from __future__ import division
-from __future__ import print_function
-
-from builtins import object, str
 
 from copy import deepcopy
 from math import radians
@@ -89,6 +85,8 @@
         if isinstance(obj, dict):
             kwargs = obj
         else:
+            if self.sim.params.forcing.key_forced is None:
+                raise ValueError("params.forcing.key_forced must be initialized.")
             kwargs = {self.sim.params.forcing.key_forced: obj}
         self.fstate.init_statespect_from(**kwargs)
 
@@ -115,6 +113,8 @@
         if isinstance(obj, dict):
             kwargs = {key: self.sim.oper.fft(value) for key, value in obj.items()}
         else:
+            if self.sim.params.forcing.key_forced is None:
+                raise ValueError("params.forcing.key_forced must be initialized.")
             kwargs = {self.sim.params.forcing.key_forced: self.sim.oper.fft(obj)}
         return kwargs
 
diff --git a/fluidsim/base/init_fields.py b/fluidsim/base/init_fields.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9pbml0X2ZpZWxkcy5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9pbml0X2ZpZWxkcy5weQ== 100644
--- a/fluidsim/base/init_fields.py
+++ b/fluidsim/base/init_fields.py
@@ -8,6 +8,4 @@
    :private-members:
 
 """
-from __future__ import division
-from __future__ import print_function
 
@@ -13,6 +11,4 @@
 
-from builtins import range
-from builtins import object
 import os
 from copy import deepcopy
 
diff --git a/fluidsim/base/output/__init__.py b/fluidsim/base/output/__init__.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvX19pbml0X18ucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvX19pbml0X18ucHk= 100644
--- a/fluidsim/base/output/__init__.py
+++ b/fluidsim/base/output/__init__.py
@@ -30,7 +30,6 @@
    :private-members:
 
 """
-from __future__ import division
 
 import os
 from glob import glob
@@ -49,101 +48,10 @@
 def create_description_xmf_files(path=None):
     """Create description xmf files for Paraview"""
 
-    if path is None:
-        path = os.getcwd()
-
-    if os.path.isdir(path):
-        paths = glob(path + "/state_phys*.nc")
-    else:
-        paths = glob(path)
-
-    if len(paths) == 0:
-        raise ValueError("No file corresponds to this path.")
-
-    for path in paths:
-        path_out = path.split(".nc")[0] + ".xmf"
-
-        base_name = os.path.basename(path)
-
-        with h5netcdf.File(path) as f:
-            ndim = 3
-
-            nx = f["/info_simul/params/oper"].attrs["nx"]
-            Lx = f["/info_simul/params/oper"].attrs["Lx"]
-            deltax = Lx / nx
-            try:
-                ny = f["/info_simul/params/oper"].attrs["ny"]
-                Ly = f["/info_simul/params/oper"].attrs["Ly"]
-                deltay = Ly / ny
-            except KeyError:
-                ndim = 1
-            try:
-                nz = f["/info_simul/params/oper"].attrs["nz"]
-                Lz = f["/info_simul/params/oper"].attrs["Lz"]
-                deltaz = Lz / nz
-            except KeyError:
-                ndim = 2
-
-            keys = list(f["/state_phys"].keys())
-
-        if ndim == 1:
-            geometry_type = "Origin_Dx"
-            dims_data = "{}".format(nx)
-            origins = "0"
-            deltaxs = "{}".format(deltax)
-        elif ndim == 2:
-            geometry_type = "Origin_DxDy"
-            dims_data = "{} {}".format(ny, nx)
-            origins = "0 0"
-            deltaxs = "{} {}".format(deltay, deltax)
-        elif ndim == 3:
-            geometry_type = "Origin_DxDyDz"
-            dims_data = "{} {} {}".format(nz, ny, nx)
-            origins = "0 0 0"
-            deltaxs = "{} {} {}".format(deltaz, deltay, deltax)
-
-        txt = """<?xml version="1.0" ?>
-<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>
-<Xdmf>
-  <Domain>
-    <Grid Name="my_uniform_grid" GridType="Uniform">
-      <Time Value="0." />
-      <Topology TopologyType="{ndim}DCoRectMesh" Dimensions="{dims_data}">
-      </Topology>
-      <Geometry GeometryType="{geometry_type}">
-        <DataItem Dimensions="{ndim}" NumberType="Float" Format="XML">
-        {origins}
-        </DataItem>
-        <DataItem Dimensions="{ndim}" NumberType="Float" Format="XML">
-        {deltaxs}
-        </DataItem>
-      </Geometry>
-""".format(
-            ndim=ndim,
-            geometry_type=geometry_type,
-            dims_data=dims_data,
-            origins=origins,
-            deltaxs=deltaxs,
-        )
-
-        for key in keys:
-            txt += """
-      <Attribute Name="{key}" AttributeType="Scalar" Center="Node">
-        <DataItem Dimensions="{dims_data}" NumberType="Float" Format="HDF">
-          {file_name}:/state_phys/{key}
-        </DataItem>
-      </Attribute>
-""".format(
-                key=key, dims_data=dims_data, file_name=base_name
-            )
-
-        txt += """    </Grid>
-  </Domain>
-</Xdmf>
-"""
-
-        with open(path_out, "w") as f:
-            f.write(txt)
+    raise DeprecationWarning(
+        "create_description_xmf_files is deprecated. "
+        "Use create_description_xmf_file instead."
+    )
 
 
 def create_description_xmf_file(path=None):
diff --git a/fluidsim/base/output/base.py b/fluidsim/base/output/base.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvYmFzZS5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvYmFzZS5weQ== 100644
--- a/fluidsim/base/output/base.py
+++ b/fluidsim/base/output/base.py
@@ -17,10 +17,6 @@
 
 """
 
-from __future__ import print_function
-
-from builtins import str
-from builtins import object
 import datetime
 import os
 import shutil
diff --git a/fluidsim/base/output/movies.py b/fluidsim/base/output/movies.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvbW92aWVzLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvbW92aWVzLnB5 100644
--- a/fluidsim/base/output/movies.py
+++ b/fluidsim/base/output/movies.py
@@ -21,9 +21,6 @@
 
 """
 
-from __future__ import division
-from __future__ import print_function
-from builtins import object
 import os
 import numpy as np
 import matplotlib.pyplot as plt
@@ -454,6 +451,6 @@
             or equation.startswith("iz=")
             or equation.startswith("z=")
         ):
-            x = self.oper.x_seq
-            y = self.oper.y_seq
+            x = self.oper.get_grid1d_seq("x")
+            y = self.oper.get_grid1d_seq("y")
         elif equation.startswith("iy=") or equation.startswith("y="):
@@ -459,4 +456,4 @@
         elif equation.startswith("iy=") or equation.startswith("y="):
-            x = self.oper.x_seq
-            y = self.oper.z_seq
+            x = self.oper.get_grid1d_seq("x")
+            y = self.oper.get_grid1d_seq("z")
         elif equation.startswith("ix=") or equation.startswith("x="):
@@ -462,6 +459,6 @@
         elif equation.startswith("ix=") or equation.startswith("x="):
-            x = self.oper.y_seq
-            y = self.oper.z_seq
+            x = self.oper.get_grid1d_seq("y")
+            y = self.oper.get_grid1d_seq("z")
         else:
             raise NotImplementedError
 
diff --git a/fluidsim/base/output/phys_fields.py b/fluidsim/base/output/phys_fields.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMucHk= 100644
--- a/fluidsim/base/output/phys_fields.py
+++ b/fluidsim/base/output/phys_fields.py
@@ -13,8 +13,6 @@
 
 """
 
-from builtins import str
-
 import re
 import os
 import datetime
@@ -258,7 +256,11 @@
         if equation is None:
             equation = self._equation
 
-        if time is None and idx_time is None:
+        if (
+            not self.sim.params.ONLY_COARSE_OPER
+            and time is None
+            and idx_time is None
+        ):
             # we get the field from the state
             field, key = self.get_field_to_plot_from_state(
                 field=key, equation=equation
@@ -318,10 +320,10 @@
             field = field[iz, ...]
         elif equation.startswith("z="):
             z = eval(equation[len("z=") :])
-            iz = abs(self.output.sim.oper.z_seq - z).argmin()
+            iz = abs(self.output.sim.oper.get_grid1d_seq("z") - z).argmin()
             field = field[iz, ...]
         elif equation.startswith("iy="):
             iy = eval(equation[len("iy=") :])
             field = field[:, iy, :]
         elif equation.startswith("y="):
             y = eval(equation[len("y=") :])
@@ -322,13 +324,13 @@
             field = field[iz, ...]
         elif equation.startswith("iy="):
             iy = eval(equation[len("iy=") :])
             field = field[:, iy, :]
         elif equation.startswith("y="):
             y = eval(equation[len("y=") :])
-            iy = abs(self.output.sim.oper.y_seq - y).argmin()
+            iy = abs(self.output.sim.oper.get_grid1d_seq("y") - y).argmin()
             field = field[:, iy, :]
         elif equation.startswith("ix="):
             ix = eval(equation[len("ix=") :])
             field = field[..., ix]
         elif equation.startswith("x="):
             x = eval(equation[len("x=") :])
@@ -329,10 +331,10 @@
             field = field[:, iy, :]
         elif equation.startswith("ix="):
             ix = eval(equation[len("ix=") :])
             field = field[..., ix]
         elif equation.startswith("x="):
             x = eval(equation[len("x=") :])
-            ix = abs(self.output.sim.oper.x_seq - x).argmin()
+            ix = abs(self.output.sim.oper.get_grid1d_seq("x") - x).argmin()
             field = field[..., ix]
         else:
             raise NotImplementedError
@@ -354,7 +356,7 @@
     return time
 
 
-class SetOfPhysFieldFiles(object):
+class SetOfPhysFieldFiles:
     """A set of physical field files.
 
     """
@@ -386,7 +388,8 @@
     ):
 
         if time is None and idx_time is None:
-            raise ValueError()
+            self.update_times()
+            time = self.times[-1]
 
         if not interpolate_time and time is not None:
             idx, time_closest = self.get_closest_time_file(time)
@@ -444,7 +447,7 @@
 
             elif equation.startswith("z="):
                 z = eval(equation[len("z=") :])
-                iz = abs(self.output.sim.oper.z_seq - z).argmin()
+                iz = abs(self.output.sim.oper.get_grid1d_seq("z") - z).argmin()
                 return dset[iz, ...]
 
             elif equation.startswith("iy="):
@@ -453,7 +456,7 @@
 
             elif equation.startswith("y="):
                 y = eval(equation[len("y=") :])
-                iy = abs(self.output.sim.oper.y_seq - y).argmin()
+                iy = abs(self.output.sim.oper.get_grid1d_seq("y") - y).argmin()
                 return dset[:, iy, :]
 
             elif equation.startswith("ix="):
@@ -462,7 +465,7 @@
 
             elif equation.startswith("x="):
                 x = eval(equation[len("x=") :])
-                ix = abs(self.output.sim.oper.x_seq - x).argmin()
+                ix = abs(self.output.sim.oper.get_grid1d_seq("x") - x).argmin()
                 return dset[..., ix]
 
             else:
diff --git a/fluidsim/base/output/phys_fields1d.py b/fluidsim/base/output/phys_fields1d.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMxZC5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMxZC5weQ== 100644
--- a/fluidsim/base/output/phys_fields1d.py
+++ b/fluidsim/base/output/phys_fields1d.py
@@ -14,8 +14,6 @@
 
 """
 
-from past.builtins import basestring
-
 import numpy as np
 
 from fluiddyn.util import mpi
@@ -45,7 +43,7 @@
         elif isinstance(field, np.ndarray):
             key_field = "given array"
             is_field_ready = True
-        elif isinstance(field, basestring):
+        elif isinstance(field, str):
             key_field = field
 
         assert key_field is not None
diff --git a/fluidsim/base/output/phys_fields2d.py b/fluidsim/base/output/phys_fields2d.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMyZC5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMyZC5weQ== 100644
--- a/fluidsim/base/output/phys_fields2d.py
+++ b/fluidsim/base/output/phys_fields2d.py
@@ -13,11 +13,6 @@
 
 """
 
-from __future__ import division
-from __future__ import print_function
-
-from past.builtins import basestring
-
 import os
 
 import numpy as np
@@ -384,7 +379,7 @@
         elif isinstance(field, np.ndarray):
             key_field = "given array"
             is_field_ready = True
-        elif isinstance(field, basestring):
+        elif isinstance(field, str):
             key_field = field
 
         assert key_field is not None
@@ -468,6 +463,6 @@
         object corresponding to a 2D contour plot.
 
         """
-        if isinstance(vecx, basestring):
+        if isinstance(vecx, str):
             vecx = self.get_field_to_plot(vecx)
 
@@ -472,6 +467,6 @@
             vecx = self.get_field_to_plot(vecx)
 
-        if isinstance(vecy, basestring):
+        if isinstance(vecy, str):
             vecy = self.get_field_to_plot(vecy)
 
         if XX is None and YY is None:
diff --git a/fluidsim/base/output/phys_fields3d.py b/fluidsim/base/output/phys_fields3d.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMzZC5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMzZC5weQ== 100644
--- a/fluidsim/base/output/phys_fields3d.py
+++ b/fluidsim/base/output/phys_fields3d.py
@@ -13,8 +13,6 @@
 
 """
 
-from past.builtins import basestring
-
 import numpy as np
 import matplotlib.pyplot as plt
 
@@ -49,6 +47,22 @@
         self._equation = equation
         self.movies._equation = equation
 
+    def _get_grid1d(self, equation):
+
+        if equation.startswith("iz=") or equation.startswith("z="):
+            x_seq = self.oper.get_grid1d_seq("x")
+            y_seq = self.oper.get_grid1d_seq("y")
+        elif equation.startswith("iy=") or equation.startswith("y="):
+            x_seq = self.oper.get_grid1d_seq("x")
+            y_seq = self.oper.get_grid1d_seq("z")
+        elif equation.startswith("ix=") or equation.startswith("x="):
+            x_seq = self.oper.get_grid1d_seq("y")
+            y_seq = self.oper.get_grid1d_seq("z")
+        else:
+            raise NotImplementedError
+
+        return x_seq, y_seq
+
     def plot(
         self,
         field=None,
@@ -106,7 +120,7 @@
         elif isinstance(field, np.ndarray):
             key_field = "given array"
             is_field_ready = True
-        elif isinstance(field, basestring):
+        elif isinstance(field, str):
             key_field = field
 
         assert key_field is not None
@@ -144,7 +158,10 @@
             # we have to get the field from a file
             self.set_of_phys_files.update_times()
             if key_field not in self.sim.state.keys_state_phys:
-                raise ValueError("key not in state.keys_state_phys")
+                raise ValueError(
+                    f'key "{key_field}" not in state.keys_state_phys '
+                    f"({self.sim.state.keys_state_phys})"
+                )
 
             field = self.get_field_to_plot(
                 key=key_field, time=time, equation=equation
@@ -163,17 +180,7 @@
             else:
                 fig, ax = self.output.figure_axe(numfig=numfig)
 
-            if equation.startswith("iz=") or equation.startswith("z="):
-                x_seq = self.oper.x_seq
-                y_seq = self.oper.y_seq
-            elif equation.startswith("iy=") or equation.startswith("y="):
-                x_seq = self.oper.x_seq
-                y_seq = self.oper.z_seq
-            elif equation.startswith("ix=") or equation.startswith("x="):
-                x_seq = self.oper.y_seq
-                y_seq = self.oper.z_seq
-            else:
-                raise NotImplementedError
+            x_seq, y_seq = self._get_grid1d(equation)
 
             [XX_seq, YY_seq] = np.meshgrid(x_seq, y_seq)
             try:
@@ -225,6 +232,6 @@
         object corresponding to a 2D contour plot.
 
         """
-        if isinstance(vecx, basestring):
+        if isinstance(vecx, str):
             vecx = self.get_field_to_plot(vecx)
 
@@ -229,7 +236,7 @@
             vecx = self.get_field_to_plot(vecx)
 
-        if isinstance(vecy, basestring):
+        if isinstance(vecy, str):
             vecy = self.get_field_to_plot(vecy)
 
         if XX is None and YY is None:
             equation = self._equation
@@ -232,18 +239,9 @@
             vecy = self.get_field_to_plot(vecy)
 
         if XX is None and YY is None:
             equation = self._equation
-            if equation.startswith("iz=") or equation.startswith("z="):
-                x_seq = self.oper.x_seq
-                y_seq = self.oper.y_seq
-            elif equation.startswith("iy=") or equation.startswith("y="):
-                x_seq = self.oper.x_seq
-                y_seq = self.oper.z_seq
-            elif equation.startswith("ix=") or equation.startswith("x="):
-                x_seq = self.oper.y_seq
-                y_seq = self.oper.z_seq
-            else:
-                raise NotImplementedError
+
+            x_seq, y_seq = self._get_grid1d(equation)
 
             [XX, YY] = np.meshgrid(x_seq, y_seq)
 
diff --git a/fluidsim/base/output/print_stdout.py b/fluidsim/base/output/print_stdout.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcHJpbnRfc3Rkb3V0LnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcHJpbnRfc3Rkb3V0LnB5 100644
--- a/fluidsim/base/output/print_stdout.py
+++ b/fluidsim/base/output/print_stdout.py
@@ -1,6 +1,3 @@
-from __future__ import print_function
-
-from builtins import object
 from time import time
 import os
 import sys
diff --git a/fluidsim/base/output/spatial_means.py b/fluidsim/base/output/spatial_means.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvc3BhdGlhbF9tZWFucy5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvc3BhdGlhbF9tZWFucy5weQ== 100644
--- a/fluidsim/base/output/spatial_means.py
+++ b/fluidsim/base/output/spatial_means.py
@@ -1,5 +1,3 @@
-from __future__ import division, print_function
-
 import os
 import numpy as np
 import json
diff --git a/fluidsim/base/output/spectra.py b/fluidsim/base/output/spectra.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvc3BlY3RyYS5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvc3BlY3RyYS5weQ== 100644
--- a/fluidsim/base/output/spectra.py
+++ b/fluidsim/base/output/spectra.py
@@ -1,5 +1,3 @@
-from __future__ import print_function
-
 import os
 
 import numpy as np
diff --git a/fluidsim/base/output/spectra3d.py b/fluidsim/base/output/spectra3d.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvc3BlY3RyYTNkLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvc3BlY3RyYTNkLnB5 100644
--- a/fluidsim/base/output/spectra3d.py
+++ b/fluidsim/base/output/spectra3d.py
@@ -1,5 +1,3 @@
-from __future__ import print_function
-
 import os
 
 import numpy as np
diff --git a/fluidsim/base/output/spectra_multidim.py b/fluidsim/base/output/spectra_multidim.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvc3BlY3RyYV9tdWx0aWRpbS5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvc3BlY3RyYV9tdWx0aWRpbS5weQ== 100644
--- a/fluidsim/base/output/spectra_multidim.py
+++ b/fluidsim/base/output/spectra_multidim.py
@@ -1,5 +1,3 @@
-from __future__ import print_function
-
 import os
 
 import numpy as np
diff --git a/fluidsim/base/output/time_signals_fft.py b/fluidsim/base/output/time_signals_fft.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvdGltZV9zaWduYWxzX2ZmdC5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvdGltZV9zaWduYWxzX2ZmdC5weQ== 100644
--- a/fluidsim/base/output/time_signals_fft.py
+++ b/fluidsim/base/output/time_signals_fft.py
@@ -1,5 +1,3 @@
-from __future__ import division
-from builtins import range
 import h5py
 
 import os
diff --git a/fluidsim/base/preprocess/pseudo_spect.py b/fluidsim/base/preprocess/pseudo_spect.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9wcmVwcm9jZXNzL3BzZXVkb19zcGVjdC5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9wcmVwcm9jZXNzL3BzZXVkb19zcGVjdC5weQ== 100644
--- a/fluidsim/base/preprocess/pseudo_spect.py
+++ b/fluidsim/base/preprocess/pseudo_spect.py
@@ -9,7 +9,6 @@
 
 
 """
-from __future__ import division, print_function
 
 import numpy as np
 from .base import PreprocessBase
diff --git a/fluidsim/base/solvers/pseudo_spect.py b/fluidsim/base/solvers/pseudo_spect.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vYmFzZS9zb2x2ZXJzL3BzZXVkb19zcGVjdC5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vYmFzZS9zb2x2ZXJzL3BzZXVkb19zcGVjdC5weQ== 100644
--- a/fluidsim/base/solvers/pseudo_spect.py
+++ b/fluidsim/base/solvers/pseudo_spect.py
@@ -17,7 +17,6 @@
    :private-members:
 
 """
-from __future__ import division
 
 import os
 
diff --git a/fluidsim/magic.py b/fluidsim/magic.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vbWFnaWMucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vbWFnaWMucHk= 100644
--- a/fluidsim/magic.py
+++ b/fluidsim/magic.py
@@ -6,7 +6,6 @@
    :private-members:
 
 """
-from __future__ import print_function
 
 from IPython.core.magic import Magics, magics_class, line_magic
 from IPython.core import magic_arguments
diff --git a/fluidsim/operators/op_finitediff.py b/fluidsim/operators/op_finitediff.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vb3BlcmF0b3JzL29wX2Zpbml0ZWRpZmYucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vb3BlcmF0b3JzL29wX2Zpbml0ZWRpZmYucHk= 100644
--- a/fluidsim/operators/op_finitediff.py
+++ b/fluidsim/operators/op_finitediff.py
@@ -13,9 +13,7 @@
    :private-members:
 
 """
-from __future__ import division
-from __future__ import print_function
-from builtins import range
+
 from past.utils import old_div
 
 import numpy as np
diff --git a/fluidsim/operators/operators1d.py b/fluidsim/operators/operators1d.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczFkLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczFkLnB5 100644
--- a/fluidsim/operators/operators1d.py
+++ b/fluidsim/operators/operators1d.py
@@ -9,8 +9,6 @@
 
 """
 
-from __future__ import division
-
 import numpy as np
 from fluiddyn.calcul.easypyfft import FFTW1DReal2Complex
 from .base import OperatorsBase1D
diff --git a/fluidsim/operators/operators2d.py b/fluidsim/operators/operators2d.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczJkLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczJkLnB5 100644
--- a/fluidsim/operators/operators2d.py
+++ b/fluidsim/operators/operators2d.py
@@ -18,6 +18,7 @@
 from fluidfft.fft2d.operators import OperatorsPseudoSpectral2D as _Operators
 
 from ..base.setofvariables import SetOfVariables
+from .. import _is_testing
 
 fp = FluidPythran()
 
@@ -21,7 +22,7 @@
 
 fp = FluidPythran()
 
-if not fp.is_transpiling and not fp.is_compiled:
+if not fp.is_transpiling and not fp.is_compiled and not _is_testing:
     warn(
         "operators2d.py has to be pythranized to be efficient! "
         "Install pythran and recompile."
@@ -797,6 +798,18 @@
                 for ikxc in range(nKxc):
                     arr[iky, ikxc] = arr_coarse[ikyc, ikxc]
 
+    def get_grid1d_seq(self, axe="x"):
+
+        if axe not in ("x", "y"):
+            raise ValueError
+
+        if self.params.ONLY_COARSE_OPER:
+            number_points = getattr(self.params.oper, "n" + axe)
+            length = getattr(self, "L" + axe)
+            return np.linspace(0, length, number_points)
+        else:
+            return getattr(self, axe + "_seq")
+
 
 # energy_arr = self.sum_wavenumbers(abs(arr)**2)
 # energy_array_coarse_after = oper_coarse.sum_wavenumbers(
diff --git a/fluidsim/operators/operators3d.py b/fluidsim/operators/operators3d.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczNkLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczNkLnB5 100644
--- a/fluidsim/operators/operators3d.py
+++ b/fluidsim/operators/operators3d.py
@@ -21,6 +21,7 @@
 from fluidsim.base.setofvariables import SetOfVariables
 
 from .operators2d import OperatorsPseudoSpectral2D as OpPseudoSpectral2D
+from .. import _is_testing
 
 fp = FluidPythran()
 
@@ -75,7 +76,7 @@
     ff_fft[np.nonzero(where_dealiased)] = 0.0
 
 
-if not fp.is_transpiling and not fp.is_compiled:
+if not fp.is_transpiling and not fp.is_compiled and not _is_testing:
     # for example if Pythran is not available
     dealiasing_variable = dealiasing_variable_numpy
     dealiasing_setofvar = dealiasing_setofvar_numpy
@@ -106,7 +107,7 @@
 
     Kx: Af
     Ky: Af
-    inv_Ksquare_nozero: Af
+    inv_K_square_nozero: Af
 
     @staticmethod
     def _complete_params_with_default(params):
@@ -114,7 +115,7 @@
         """
         attribs = {
             "type_fft": "default",
-            "type_fft2d": "default",
+            "type_fft2d": "sequential",
             "coef_dealiasing": 2.0 / 3,
             "nx": 48,
             "ny": 48,
@@ -166,5 +167,10 @@
 
         self.params = params
         self.axes = ("z", "y", "x")
+
+        params.oper.nx = int(params.oper.nx)
+        params.oper.ny = int(params.oper.ny)
+        params.oper.nz = int(params.oper.nz)
+
         if params.ONLY_COARSE_OPER:
             nx = ny = nz = 4
@@ -169,7 +175,3 @@
         if params.ONLY_COARSE_OPER:
             nx = ny = nz = 4
-            # Unchanged, but type cast into integers
-            params.oper.nx = int(params.oper.nx)
-            params.oper.ny = int(params.oper.ny)
-            params.oper.nz = int(params.oper.nz)
         else:
@@ -175,7 +177,7 @@
         else:
-            nx = params.oper.nx = int(params.oper.nx)
-            ny = params.oper.ny = int(params.oper.ny)
-            nz = params.oper.nz = int(params.oper.nz)
+            nx = params.oper.nx
+            ny = params.oper.ny
+            nz = params.oper.nz
 
         super(OperatorsPseudoSpectral3D, self).__init__(
             nx,
@@ -195,7 +197,7 @@
 
         if (
             any([fft.startswith(s) for s in ["fluidfft.fft2d.", "fft2d."]])
-            or fft == "default"
+            or fft in ("default", "sequential")
             or fft is None
         ):
             self.oper2d = OpPseudoSpectral2D(params2d)
@@ -286,11 +288,11 @@
         """
 
         divh_fft = 1j * (self.Kx * vx_fft + self.Ky * vy_fft)
-        urx_fft = vx_fft - divh_fft * self.Kx * self.inv_Ksquare_nozero
-        ury_fft = vy_fft - divh_fft * self.Ky * self.inv_Ksquare_nozero
+        urx_fft = vx_fft - divh_fft * self.Kx * self.inv_K_square_nozero
+        ury_fft = vy_fft - divh_fft * self.Ky * self.inv_K_square_nozero
 
         udx_fft = vx_fft - urx_fft
         udy_fft = vy_fft - ury_fft
 
         return urx_fft, ury_fft, udx_fft, udy_fft
 
@@ -291,9 +293,21 @@
 
         udx_fft = vx_fft - urx_fft
         udy_fft = vy_fft - ury_fft
 
         return urx_fft, ury_fft, udx_fft, udy_fft
 
+    def get_grid1d_seq(self, axe="x"):
+
+        if axe not in ("x", "y", "z"):
+            raise ValueError
+
+        if self.params.ONLY_COARSE_OPER:
+            number_points = getattr(self.params.oper, "n" + axe)
+            length = getattr(self, "L" + axe)
+            return np.linspace(0, length, number_points)
+        else:
+            return getattr(self, axe + "_seq")
+
 
 def _ik_from_ikc(ikc, nkc, nk):
     if ikc <= nkc / 2.0:
diff --git a/fluidsim/operators/test/test_operators2d.py b/fluidsim/operators/test/test_operators2d.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vb3BlcmF0b3JzL3Rlc3QvdGVzdF9vcGVyYXRvcnMyZC5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vb3BlcmF0b3JzL3Rlc3QvdGVzdF9vcGVyYXRvcnMyZC5weQ== 100644
--- a/fluidsim/operators/test/test_operators2d.py
+++ b/fluidsim/operators/test/test_operators2d.py
@@ -1,5 +1,3 @@
-from __future__ import division
-
 import unittest
 import numpy as np
 import sys
diff --git a/fluidsim/solvers/ad1d/init_fields.py b/fluidsim/solvers/ad1d/init_fields.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9hZDFkL2luaXRfZmllbGRzLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9hZDFkL2luaXRfZmllbGRzLnB5 100644
--- a/fluidsim/solvers/ad1d/init_fields.py
+++ b/fluidsim/solvers/ad1d/init_fields.py
@@ -11,7 +11,6 @@
    :members:
 
 """
-from __future__ import division
 from past.utils import old_div
 import numpy as np
 
diff --git a/fluidsim/solvers/ad1d/pseudo_spect/state.py b/fluidsim/solvers/ad1d/pseudo_spect/state.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9hZDFkL3BzZXVkb19zcGVjdC9zdGF0ZS5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9hZDFkL3BzZXVkb19zcGVjdC9zdGF0ZS5weQ== 100644
--- a/fluidsim/solvers/ad1d/pseudo_spect/state.py
+++ b/fluidsim/solvers/ad1d/pseudo_spect/state.py
@@ -1,8 +1,6 @@
 """State for the AD1D solver (:mod:`fluidsim.solvers.ad1d.pseudo_spect.state`)
 ==============================================================================
 """
-from __future__ import print_function
-
 
 from fluidsim.base.state import StatePseudoSpectral
 
diff --git a/fluidsim/solvers/ad1d/solver.py b/fluidsim/solvers/ad1d/solver.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9hZDFkL3NvbHZlci5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9hZDFkL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/ad1d/solver.py
+++ b/fluidsim/solvers/ad1d/solver.py
@@ -8,7 +8,6 @@
    :private-members:
 
 """
-from __future__ import print_function
 
 from fluidsim.base.solvers.base import SimulBase
 from fluidsim.base.setofvariables import SetOfVariables
diff --git a/fluidsim/solvers/ad1d/state.py b/fluidsim/solvers/ad1d/state.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9hZDFkL3N0YXRlLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9hZDFkL3N0YXRlLnB5 100644
--- a/fluidsim/solvers/ad1d/state.py
+++ b/fluidsim/solvers/ad1d/state.py
@@ -1,8 +1,6 @@
 """State for the AD1D solver (:mod:`fluidsim.solvers.ad1d.state`)
 =======================================================================
 """
-from __future__ import print_function
-
 
 from fluidsim.base.state import StateBase
 
diff --git a/fluidsim/solvers/models0d/lorenz/output/__init__.py b/fluidsim/solvers/models0d/lorenz/output/__init__.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9tb2RlbHMwZC9sb3Jlbnovb3V0cHV0L19faW5pdF9fLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9tb2RlbHMwZC9sb3Jlbnovb3V0cHV0L19faW5pdF9fLnB5 100644
--- a/fluidsim/solvers/models0d/lorenz/output/__init__.py
+++ b/fluidsim/solvers/models0d/lorenz/output/__init__.py
@@ -15,7 +15,6 @@
    :private-members:
 
 """
-from __future__ import division
 
 from fluidsim.base.output import OutputBase
 
diff --git a/fluidsim/solvers/models0d/lorenz/output/print_stdout.py b/fluidsim/solvers/models0d/lorenz/output/print_stdout.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9tb2RlbHMwZC9sb3Jlbnovb3V0cHV0L3ByaW50X3N0ZG91dC5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9tb2RlbHMwZC9sb3Jlbnovb3V0cHV0L3ByaW50X3N0ZG91dC5weQ== 100644
--- a/fluidsim/solvers/models0d/lorenz/output/print_stdout.py
+++ b/fluidsim/solvers/models0d/lorenz/output/print_stdout.py
@@ -7,9 +7,6 @@
 
 """
 
-from __future__ import print_function, division
-
-from builtins import range
 import numpy as np
 import matplotlib.pyplot as plt
 from mpl_toolkits.mplot3d import Axes3D
diff --git a/fluidsim/solvers/models0d/lorenz/solver.py b/fluidsim/solvers/models0d/lorenz/solver.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9tb2RlbHMwZC9sb3Jlbnovc29sdmVyLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9tb2RlbHMwZC9sb3Jlbnovc29sdmVyLnB5 100644
--- a/fluidsim/solvers/models0d/lorenz/solver.py
+++ b/fluidsim/solvers/models0d/lorenz/solver.py
@@ -17,8 +17,6 @@
 
 """
 
-from __future__ import division
-
 from math import sqrt
 
 from fluidsim.base.setofvariables import SetOfVariables
diff --git a/fluidsim/solvers/models0d/predaprey/output/__init__.py b/fluidsim/solvers/models0d/predaprey/output/__init__.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9tb2RlbHMwZC9wcmVkYXByZXkvb3V0cHV0L19faW5pdF9fLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9tb2RlbHMwZC9wcmVkYXByZXkvb3V0cHV0L19faW5pdF9fLnB5 100644
--- a/fluidsim/solvers/models0d/predaprey/output/__init__.py
+++ b/fluidsim/solvers/models0d/predaprey/output/__init__.py
@@ -15,7 +15,6 @@
    :private-members:
 
 """
-from __future__ import division
 
 from math import log
 
diff --git a/fluidsim/solvers/models0d/predaprey/output/print_stdout.py b/fluidsim/solvers/models0d/predaprey/output/print_stdout.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9tb2RlbHMwZC9wcmVkYXByZXkvb3V0cHV0L3ByaW50X3N0ZG91dC5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9tb2RlbHMwZC9wcmVkYXByZXkvb3V0cHV0L3ByaW50X3N0ZG91dC5weQ== 100644
--- a/fluidsim/solvers/models0d/predaprey/output/print_stdout.py
+++ b/fluidsim/solvers/models0d/predaprey/output/print_stdout.py
@@ -7,9 +7,6 @@
 
 """
 
-from __future__ import print_function, division
-
-from builtins import range
 import numpy as np
 
 from fluidsim.base.output.print_stdout import PrintStdOutBase
diff --git a/fluidsim/solvers/models0d/predaprey/solver.py b/fluidsim/solvers/models0d/predaprey/solver.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9tb2RlbHMwZC9wcmVkYXByZXkvc29sdmVyLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9tb2RlbHMwZC9wcmVkYXByZXkvc29sdmVyLnB5 100644
--- a/fluidsim/solvers/models0d/predaprey/solver.py
+++ b/fluidsim/solvers/models0d/predaprey/solver.py
@@ -17,8 +17,6 @@
 
 """
 
-from __future__ import division
-
 from fluidsim.base.setofvariables import SetOfVariables
 
 from fluidsim.base.solvers.base import InfoSolverBase, SimulBase
diff --git a/fluidsim/solvers/ns2d/forcing.py b/fluidsim/solvers/ns2d/forcing.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL2ZvcmNpbmcucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL2ZvcmNpbmcucHk= 100644
--- a/fluidsim/solvers/ns2d/forcing.py
+++ b/fluidsim/solvers/ns2d/forcing.py
@@ -5,8 +5,6 @@
    :members:
 
 """
-from __future__ import division
-from __future__ import print_function
 
 from fluidsim.base.forcing import ForcingBasePseudoSpectral
 from fluidsim.base.forcing.anisotropic import (
diff --git a/fluidsim/solvers/ns2d/init_fields.py b/fluidsim/solvers/ns2d/init_fields.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL2luaXRfZmllbGRzLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL2luaXRfZmllbGRzLnB5 100644
--- a/fluidsim/solvers/ns2d/init_fields.py
+++ b/fluidsim/solvers/ns2d/init_fields.py
@@ -14,5 +14,4 @@
    :members:
 
 """
-from __future__ import division
 
@@ -18,5 +17,4 @@
 
-from builtins import range
 import numpy as np
 
 from fluiddyn.util import mpi
diff --git a/fluidsim/solvers/ns2d/output/__init__.py b/fluidsim/solvers/ns2d/output/__init__.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9fX2luaXRfXy5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9fX2luaXRfXy5weQ== 100644
--- a/fluidsim/solvers/ns2d/output/__init__.py
+++ b/fluidsim/solvers/ns2d/output/__init__.py
@@ -18,7 +18,6 @@
    :private-members:
 
 """
-from __future__ import division
 
 import numpy as np
 
diff --git a/fluidsim/solvers/ns2d/output/print_stdout.py b/fluidsim/solvers/ns2d/output/print_stdout.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9wcmludF9zdGRvdXQucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9wcmludF9zdGRvdXQucHk= 100644
--- a/fluidsim/solvers/ns2d/output/print_stdout.py
+++ b/fluidsim/solvers/ns2d/output/print_stdout.py
@@ -7,9 +7,6 @@
 
 """
 
-from __future__ import print_function, division
-
-from builtins import range
 import numpy as np
 
 from fluidsim.base.output.print_stdout import PrintStdOutBase
diff --git a/fluidsim/solvers/ns2d/output/spatial_means.py b/fluidsim/solvers/ns2d/output/spatial_means.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9zcGF0aWFsX21lYW5zLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9zcGF0aWFsX21lYW5zLnB5 100644
--- a/fluidsim/solvers/ns2d/output/spatial_means.py
+++ b/fluidsim/solvers/ns2d/output/spatial_means.py
@@ -7,9 +7,6 @@
 
 """
 
-from __future__ import division, print_function
-
-from builtins import range
 import os
 import numpy as np
 import matplotlib.pyplot as plt
diff --git a/fluidsim/solvers/ns2d/output/spect_energy_budget.py b/fluidsim/solvers/ns2d/output/spect_energy_budget.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9zcGVjdF9lbmVyZ3lfYnVkZ2V0LnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9zcGVjdF9lbmVyZ3lfYnVkZ2V0LnB5 100644
--- a/fluidsim/solvers/ns2d/output/spect_energy_budget.py
+++ b/fluidsim/solvers/ns2d/output/spect_energy_budget.py
@@ -6,6 +6,4 @@
    :private-members:
 
 """
-from __future__ import division
-from __future__ import print_function
 
@@ -11,5 +9,4 @@
 
-from builtins import range
 import numpy as np
 import h5py
 
diff --git a/fluidsim/solvers/ns2d/output/spectra.py b/fluidsim/solvers/ns2d/output/spectra.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9zcGVjdHJhLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9zcGVjdHJhLnB5 100644
--- a/fluidsim/solvers/ns2d/output/spectra.py
+++ b/fluidsim/solvers/ns2d/output/spectra.py
@@ -6,6 +6,4 @@
    :private-members:
 
 """
-from __future__ import division
-from __future__ import print_function
 
@@ -11,5 +9,4 @@
 
-from builtins import range
 import h5py
 
 import numpy as np
diff --git a/fluidsim/solvers/ns2d/output/spectra_multidim.py b/fluidsim/solvers/ns2d/output/spectra_multidim.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9zcGVjdHJhX211bHRpZGltLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9zcGVjdHJhX211bHRpZGltLnB5 100644
--- a/fluidsim/solvers/ns2d/output/spectra_multidim.py
+++ b/fluidsim/solvers/ns2d/output/spectra_multidim.py
@@ -6,8 +6,6 @@
    :private-members:
 
 """
-from __future__ import division
-from __future__ import print_function
 
 from fluidsim.base.output.spectra_multidim import SpectraMultiDim
 
diff --git a/fluidsim/solvers/ns2d/state.py b/fluidsim/solvers/ns2d/state.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0YXRlLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0YXRlLnB5 100644
--- a/fluidsim/solvers/ns2d/state.py
+++ b/fluidsim/solvers/ns2d/state.py
@@ -6,7 +6,6 @@
    :private-members:
 
 """
-from __future__ import print_function
 
 import numpy as np
 
diff --git a/fluidsim/solvers/ns2d/strat/diag/solver.py b/fluidsim/solvers/ns2d/strat/diag/solver.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L2RpYWcvc29sdmVyLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L2RpYWcvc29sdmVyLnB5 100644
--- a/fluidsim/solvers/ns2d/strat/diag/solver.py
+++ b/fluidsim/solvers/ns2d/strat/diag/solver.py
@@ -6,5 +6,4 @@
    :private-members:
 
 """
-from __future__ import division
 
@@ -10,5 +9,4 @@
 
-from past.utils import old_div
 from fluidsim.base.setofvariables import SetOfVariables
 
 from fluidsim.base.solvers.pseudo_spect import (
@@ -124,5 +122,5 @@
 
     # params.oper.type_fft = 'FFTWPY'
 
-    delta_x = old_div(params.oper.Lx, params.oper.nx)
+    delta_x = params.oper.Lx / params.oper.nx
     params.nu_8 = (
@@ -128,8 +126,5 @@
     params.nu_8 = (
-        2.0
-        * 10e-1
-        * params.forcing.forcing_rate ** (old_div(1.0, 3))
-        * delta_x ** 8
+        2.0 * 10e-1 * params.forcing.forcing_rate ** (1.0 / 3) * delta_x ** 8
     )
 
     params.time_stepping.t_end = 1.0
diff --git a/fluidsim/solvers/ns2d/strat/init_fields.py b/fluidsim/solvers/ns2d/strat/init_fields.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L2luaXRfZmllbGRzLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L2luaXRfZmllbGRzLnB5 100644
--- a/fluidsim/solvers/ns2d/strat/init_fields.py
+++ b/fluidsim/solvers/ns2d/strat/init_fields.py
@@ -15,8 +15,6 @@
 
 """
 
-from __future__ import print_function
-
 import random
 import numpy as np
 
diff --git a/fluidsim/solvers/ns2d/strat/output/phys_fields.py b/fluidsim/solvers/ns2d/strat/output/phys_fields.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9waHlzX2ZpZWxkcy5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9waHlzX2ZpZWxkcy5weQ== 100644
--- a/fluidsim/solvers/ns2d/strat/output/phys_fields.py
+++ b/fluidsim/solvers/ns2d/strat/output/phys_fields.py
@@ -8,7 +8,6 @@
    :private-members:
 
 """
-from __future__ import print_function
 
 from fluidsim.base.output.phys_fields2d import PhysFieldsBase2D
 
diff --git a/fluidsim/solvers/ns2d/strat/output/print_stdout.py b/fluidsim/solvers/ns2d/strat/output/print_stdout.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9wcmludF9zdGRvdXQucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9wcmludF9zdGRvdXQucHk= 100644
--- a/fluidsim/solvers/ns2d/strat/output/print_stdout.py
+++ b/fluidsim/solvers/ns2d/strat/output/print_stdout.py
@@ -7,8 +7,6 @@
 
 """
 
-from __future__ import print_function, division
-
 import os
 
 import numpy as np
diff --git a/fluidsim/solvers/ns2d/strat/output/spatial_means.py b/fluidsim/solvers/ns2d/strat/output/spatial_means.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9zcGF0aWFsX21lYW5zLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9zcGF0aWFsX21lYW5zLnB5 100644
--- a/fluidsim/solvers/ns2d/strat/output/spatial_means.py
+++ b/fluidsim/solvers/ns2d/strat/output/spatial_means.py
@@ -7,8 +7,6 @@
 
 """
 
-from __future__ import division, print_function
-
 import os
 import numpy as np
 import matplotlib.pyplot as plt
diff --git a/fluidsim/solvers/ns2d/strat/output/spect_energy_budget.py b/fluidsim/solvers/ns2d/strat/output/spect_energy_budget.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9zcGVjdF9lbmVyZ3lfYnVkZ2V0LnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9zcGVjdF9lbmVyZ3lfYnVkZ2V0LnB5 100644
--- a/fluidsim/solvers/ns2d/strat/output/spect_energy_budget.py
+++ b/fluidsim/solvers/ns2d/strat/output/spect_energy_budget.py
@@ -6,8 +6,6 @@
    :private-members:
 
 """
-from __future__ import division
-from __future__ import print_function
 
 import numpy as np
 import h5py
diff --git a/fluidsim/solvers/ns2d/strat/output/spectra.py b/fluidsim/solvers/ns2d/strat/output/spectra.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9zcGVjdHJhLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9zcGVjdHJhLnB5 100644
--- a/fluidsim/solvers/ns2d/strat/output/spectra.py
+++ b/fluidsim/solvers/ns2d/strat/output/spectra.py
@@ -6,7 +6,6 @@
    :private-members:
 
 """
-from __future__ import print_function
 
 import h5py
 
@@ -229,4 +228,5 @@
 
             # Set limits axis y
             ymin = E_kx_plot[0:50] * kx_plot[0:50] ** coef_compensate_kx
+            ymin = ymin[ymin > 0]
             ymax = E_kx_plot * kx_plot ** coef_compensate_kx
@@ -232,5 +232,5 @@
             ymax = E_kx_plot * kx_plot ** coef_compensate_kx
-            ax1.set_ylim(ymin=1e-1 * ymin.min(), ymax=ymax.max())
+            ax1.set_ylim((1e-1 * ymin.min(), ymax.max()))
 
             # Parameters figure E(k_y)
             fig, ax2 = self.output.figure_axe()
diff --git a/fluidsim/solvers/ns2d/strat/output/spectra_multidim.py b/fluidsim/solvers/ns2d/strat/output/spectra_multidim.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9zcGVjdHJhX211bHRpZGltLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9zcGVjdHJhX211bHRpZGltLnB5 100644
--- a/fluidsim/solvers/ns2d/strat/output/spectra_multidim.py
+++ b/fluidsim/solvers/ns2d/strat/output/spectra_multidim.py
@@ -29,6 +29,10 @@
         am_fft = self.sim.state.compute("am_fft")
 
         # Computes multidimensional spectra
-        spectrumkykx_EK = self.oper.compute_spectrum_kykx(energyK_fft, folded=False)
-        spectrumkykx_EA = self.oper.compute_spectrum_kykx(energyA_fft, folded=False)
+        spectrumkykx_EK = self.oper.compute_spectrum_kykx(
+            energyK_fft, folded=False
+        )
+        spectrumkykx_EA = self.oper.compute_spectrum_kykx(
+            energyA_fft, folded=False
+        )
 
@@ -34,4 +38,4 @@
 
-        energy_ap_fft = abs(ap_fft)**2
-        energy_am_fft = abs(am_fft)**2
+        energy_ap_fft = abs(ap_fft) ** 2
+        energy_am_fft = abs(am_fft) ** 2
 
@@ -37,8 +41,12 @@
 
-        spectrumkykx_ap = self.oper.compute_spectrum_kykx(energy_ap_fft, folded=False)
-        spectrumkykx_am = self.oper.compute_spectrum_kykx(energy_am_fft, folded=False)
+        spectrumkykx_ap = self.oper.compute_spectrum_kykx(
+            energy_ap_fft, folded=False
+        )
+        spectrumkykx_am = self.oper.compute_spectrum_kykx(
+            energy_am_fft, folded=False
+        )
 
         dict_spectra = {
             "spectrumkykx_EK": spectrumkykx_EK,
             "spectrumkykx_EA": spectrumkykx_EA,
             "spectrumkykx_ap": spectrumkykx_ap,
@@ -40,9 +48,9 @@
 
         dict_spectra = {
             "spectrumkykx_EK": spectrumkykx_EK,
             "spectrumkykx_EA": spectrumkykx_EA,
             "spectrumkykx_ap": spectrumkykx_ap,
-            "spectrumkykx_am": spectrumkykx_am
+            "spectrumkykx_am": spectrumkykx_am,
         }
 
         return dict_spectra
@@ -87,6 +95,6 @@
                 data = f["spectrumkykx_EA"]
             elif key == "ap_fft" or not key:
                 data = f["spectrumkykx_ap"]
-                text_plot = "$\hat{a}_+$"
+                text_plot = r"$\hat{a}_+$"
             elif key == "am_fft":
                 data = f["spectrumkykx_am"]
@@ -91,6 +99,6 @@
             elif key == "am_fft":
                 data = f["spectrumkykx_am"]
-                text_plot = "$\hat{a}_-$"
+                text_plot = r"$\hat{a}_-$"
             else:
                 raise ValueError("Key unknown.")
 
@@ -103,6 +111,6 @@
             itmin = np.argmin(abs(times - tmin))
             itmax = np.argmin(abs(times - tmax))
 
-            data_plot = (data[itmin : itmax, :, :]).mean(0)
+            data_plot = (data[itmin:itmax, :, :]).mean(0)
 
             if key != "ap_fft":
@@ -107,6 +115,8 @@
 
             if key != "ap_fft":
-                vmax_modified =  ((ap_fft_spectrum[itmin : itmax, :, :]).mean(0)).max()
+                vmax_modified = (
+                    (ap_fft_spectrum[itmin:itmax, :, :]).mean(0)
+                ).max()
 
         # Create array kz with negative values
         kz = 2 * np.pi * np.fft.fftfreq(oper.ny, oper.Ly / oper.ny)
@@ -110,7 +120,7 @@
 
         # Create array kz with negative values
         kz = 2 * np.pi * np.fft.fftfreq(oper.ny, oper.Ly / oper.ny)
-        kz[kz.shape[0]//2] *= -1
+        kz[kz.shape[0] // 2] *= -1
 
         # Create mesh of wave-numbers
         KX, KZ = np.meshgrid(kx, kz)
@@ -136,9 +146,13 @@
 
         # Modify grid
         kz_modified = np.empty_like(kz)
-        kz_modified[0:kz_modified.shape[0]//2 - 1] = kz[kz_modified.shape[0]//2 + 1:]
-        kz_modified[kz_modified.shape[0]//2 - 1:] = kz[0:kz_modified.shape[0]//2 + 1]
+        kz_modified[0 : kz_modified.shape[0] // 2 - 1] = kz[
+            kz_modified.shape[0] // 2 + 1 :
+        ]
+        kz_modified[kz_modified.shape[0] // 2 - 1 :] = kz[
+            0 : kz_modified.shape[0] // 2 + 1
+        ]
 
         KX, KZ = np.meshgrid(kx, kz_modified)
 
         data_plot_modified = np.empty_like(data_plot)
@@ -141,9 +155,13 @@
 
         KX, KZ = np.meshgrid(kx, kz_modified)
 
         data_plot_modified = np.empty_like(data_plot)
-        data_plot_modified[0:kz_modified.shape[0]//2 - 1, :] = data_plot[kz_modified.shape[0]//2 + 1:, :]
-        data_plot_modified[kz_modified.shape[0]//2 - 1:, :] = data_plot[0:kz_modified.shape[0]//2 + 1, :]
+        data_plot_modified[0 : kz_modified.shape[0] // 2 - 1, :] = data_plot[
+            kz_modified.shape[0] // 2 + 1 :, :
+        ]
+        data_plot_modified[kz_modified.shape[0] // 2 - 1 :, :] = data_plot[
+            0 : kz_modified.shape[0] // 2 + 1, :
+        ]
 
         # Vmax is 1% of the maximum value.
         vmin = 0
@@ -160,4 +178,5 @@
 
         if isinstance(pforcing.tcrandom_anisotropic.angle, str):
             from math import radians
+
             angle = radians(
@@ -163,5 +182,6 @@
             angle = radians(
-                float((pforcing.tcrandom_anisotropic.angle).split("°")[0]))
+                float((pforcing.tcrandom_anisotropic.angle).split("°")[0])
+            )
         else:
             angle = pforcing.tcrandom_anisotropic.angle
 
@@ -165,7 +185,6 @@
         else:
             angle = pforcing.tcrandom_anisotropic.angle
 
-
         x_rect = np.sin(angle) * deltak * pforcing.nkmin_forcing
 
         z_rect = np.cos(angle) * deltak * pforcing.nkmin_forcing
@@ -174,9 +193,16 @@
 
         height = abs(z_rect - np.cos(angle) * deltak * pforcing.nkmax_forcing)
 
-        rect1 = patches.Rectangle((x_rect,z_rect),width,height,linewidth=1,edgecolor='r',facecolor='none')
+        rect1 = patches.Rectangle(
+            (x_rect, z_rect),
+            width,
+            height,
+            linewidth=1,
+            edgecolor="r",
+            facecolor="none",
+        )
 
         ax.add_patch(rect1)
 
         if pforcing.tcrandom_anisotropic.kz_negative_enable:
             rect2 = patches.Rectangle(
@@ -178,10 +204,15 @@
 
         ax.add_patch(rect1)
 
         if pforcing.tcrandom_anisotropic.kz_negative_enable:
             rect2 = patches.Rectangle(
-                (x_rect,-(z_rect + height)), width, height, linewidth=1,
-                edgecolor='r',facecolor='none')
+                (x_rect, -(z_rect + height)),
+                width,
+                height,
+                linewidth=1,
+                edgecolor="r",
+                facecolor="none",
+            )
 
             ax.add_patch(rect2)
 
@@ -192,6 +223,6 @@
                 width=2 * pforcing.nkmin_forcing * deltak,
                 height=2 * pforcing.nkmin_forcing * deltak,
                 angle=0,
-                theta1=-90.,
-                theta2=90.,
+                theta1=-90.0,
+                theta2=90.0,
                 linestyle="-.",
@@ -197,5 +228,5 @@
                 linestyle="-.",
-                color="red"
+                color="red",
             )
         )
         ax.add_patch(
@@ -205,5 +236,5 @@
                 height=2 * pforcing.nkmax_forcing * deltak,
                 angle=0,
                 theta1=-90,
-                theta2=90.,
+                theta2=90.0,
                 linestyle="-.",
@@ -209,5 +240,5 @@
                 linestyle="-.",
-                color="red"
+                color="red",
             )
         )
 
@@ -215,11 +246,13 @@
         ikx_text = np.argmin(abs(kx - kx[ikx] * 0.7))
         ikz_text = np.argmin(abs(kz - kz[ikz] * 0.7))
 
-        ax.text(kx[ikx_text],
-                kz[ikz_text],
-                "{}".format(text_plot),
-                color="white",
-                fontsize=15)
+        ax.text(
+            kx[ikx_text],
+            kz[ikz_text],
+            "{}".format(text_plot),
+            color="white",
+            fontsize=15,
+        )
 
         # Colorbar
         fig.colorbar(spectrum)
diff --git a/fluidsim/solvers/ns2d/strat/output/test/__init__.py b/fluidsim/solvers/ns2d/strat/output/test/__init__.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC90ZXN0L19faW5pdF9fLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC90ZXN0L19faW5pdF9fLnB5 100644
--- a/fluidsim/solvers/ns2d/strat/output/test/__init__.py
+++ b/fluidsim/solvers/ns2d/strat/output/test/__init__.py
@@ -1,5 +1,3 @@
-from __future__ import print_function
-
 import unittest
 import numpy as np
 
diff --git a/fluidsim/solvers/ns2d/strat/output/test/test_print_stdout.py b/fluidsim/solvers/ns2d/strat/output/test/test_print_stdout.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC90ZXN0L3Rlc3RfcHJpbnRfc3Rkb3V0LnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC90ZXN0L3Rlc3RfcHJpbnRfc3Rkb3V0LnB5 100644
--- a/fluidsim/solvers/ns2d/strat/output/test/test_print_stdout.py
+++ b/fluidsim/solvers/ns2d/strat/output/test/test_print_stdout.py
@@ -1,5 +1,3 @@
-from __future__ import print_function
-
 import unittest
 import numpy as np
 
diff --git a/fluidsim/solvers/ns2d/strat/output/test/test_spatial_means.py b/fluidsim/solvers/ns2d/strat/output/test/test_spatial_means.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC90ZXN0L3Rlc3Rfc3BhdGlhbF9tZWFucy5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC90ZXN0L3Rlc3Rfc3BhdGlhbF9tZWFucy5weQ== 100644
--- a/fluidsim/solvers/ns2d/strat/output/test/test_spatial_means.py
+++ b/fluidsim/solvers/ns2d/strat/output/test/test_spatial_means.py
@@ -1,5 +1,3 @@
-from __future__ import print_function
-
 import unittest
 import numpy as np
 
@@ -3,7 +1,6 @@
 import unittest
 import numpy as np
 
-from fluidsim.solvers.test.test_ns import run_mini_simul
 from fluidsim.solvers.ns2d.strat.output.test import TestNS2DStrat, mpi
 
 
diff --git a/fluidsim/solvers/ns2d/strat/output/test/test_spect_energy_budg.py b/fluidsim/solvers/ns2d/strat/output/test/test_spect_energy_budg.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC90ZXN0L3Rlc3Rfc3BlY3RfZW5lcmd5X2J1ZGcucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC90ZXN0L3Rlc3Rfc3BlY3RfZW5lcmd5X2J1ZGcucHk= 100644
--- a/fluidsim/solvers/ns2d/strat/output/test/test_spect_energy_budg.py
+++ b/fluidsim/solvers/ns2d/strat/output/test/test_spect_energy_budg.py
@@ -1,5 +1,3 @@
-from __future__ import print_function
-
 import unittest
 import numpy as np
 
diff --git a/fluidsim/solvers/ns2d/strat/output/test/test_spectra.py b/fluidsim/solvers/ns2d/strat/output/test/test_spectra.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC90ZXN0L3Rlc3Rfc3BlY3RyYS5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC90ZXN0L3Rlc3Rfc3BlY3RyYS5weQ== 100644
--- a/fluidsim/solvers/ns2d/strat/output/test/test_spectra.py
+++ b/fluidsim/solvers/ns2d/strat/output/test/test_spectra.py
@@ -1,5 +1,3 @@
-from __future__ import print_function
-
 import unittest
 import numpy as np
 
diff --git a/fluidsim/solvers/ns2d/strat/output/test/test_spectra_multidim.py b/fluidsim/solvers/ns2d/strat/output/test/test_spectra_multidim.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC90ZXN0L3Rlc3Rfc3BlY3RyYV9tdWx0aWRpbS5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC90ZXN0L3Rlc3Rfc3BlY3RyYV9tdWx0aWRpbS5weQ== 100644
--- a/fluidsim/solvers/ns2d/strat/output/test/test_spectra_multidim.py
+++ b/fluidsim/solvers/ns2d/strat/output/test/test_spectra_multidim.py
@@ -1,5 +1,3 @@
-from __future__ import print_function
-
 import unittest
 import numpy as np
 
diff --git a/fluidsim/solvers/ns2d/strat/test_solver.py b/fluidsim/solvers/ns2d/strat/test_solver.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L3Rlc3Rfc29sdmVyLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L3Rlc3Rfc29sdmVyLnB5 100644
--- a/fluidsim/solvers/ns2d/strat/test_solver.py
+++ b/fluidsim/solvers/ns2d/strat/test_solver.py
@@ -1,5 +1,4 @@
 # -*- coding: utf-8 -*-
-from __future__ import division
 
 import unittest
 import shutil
diff --git a/fluidsim/solvers/ns2d/strat/time_stepping.py b/fluidsim/solvers/ns2d/strat/time_stepping.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L3RpbWVfc3RlcHBpbmcucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L3RpbWVfc3RlcHBpbmcucHk= 100644
--- a/fluidsim/solvers/ns2d/strat/time_stepping.py
+++ b/fluidsim/solvers/ns2d/strat/time_stepping.py
@@ -10,8 +10,6 @@
 
 """
 
-from __future__ import print_function
-
 import numpy as np
 
 from math import pi
diff --git a/fluidsim/solvers/ns2d/test/test_solver.py b/fluidsim/solvers/ns2d/test/test_solver.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczJkL3Rlc3QvdGVzdF9zb2x2ZXIucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczJkL3Rlc3QvdGVzdF9zb2x2ZXIucHk= 100644
--- a/fluidsim/solvers/ns2d/test/test_solver.py
+++ b/fluidsim/solvers/ns2d/test/test_solver.py
@@ -1,6 +1,4 @@
-from __future__ import division
-
 import unittest
 import shutil
 
 import numpy as np
@@ -3,7 +1,8 @@
 import unittest
 import shutil
 
 import numpy as np
+import matplotlib.pyplot as plt
 
 import fluidsim as fls
 
@@ -109,6 +108,17 @@
                 sim2.output.increments.plot()
                 sim2.output.increments.load_pdf_from_file()
 
+                sim2.output.phys_fields.animate(
+                    "ux",
+                    dt_frame_in_sec=1e-6,
+                    dt_equations=0.3,
+                    repeat=False,
+                    clim=(-1, 1),
+                    save_file=False,
+                    numfig=1,
+                )
+                sim2.output.phys_fields.plot()
+
             # `compute('q')` two times for better coverage...
             sim.state.get_var("q")
             sim.state.get_var("q")
@@ -132,10 +142,7 @@
                     save_file=False,
                     numfig=1,
                 )
-
-
-# class TestSolverNS2DFluidfft(TestSolverNS2D):
-#     Simul = Simul2
+        plt.close("all")
 
 
 if __name__ == "__main__":
diff --git a/fluidsim/solvers/ns3d/bouss/solver.py b/fluidsim/solvers/ns3d/bouss/solver.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczNkL2JvdXNzL3NvbHZlci5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczNkL2JvdXNzL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/ns3d/bouss/solver.py
+++ b/fluidsim/solvers/ns3d/bouss/solver.py
@@ -12,7 +12,6 @@
    :private-members:
 
 """
-from __future__ import division
 
 from fluidfft.fft3d.operators import vector_product
 
diff --git a/fluidsim/solvers/ns3d/bouss/test_ns3d_bouss.py b/fluidsim/solvers/ns3d/bouss/test_ns3d_bouss.py
new file mode 100644
index 0000000000000000000000000000000000000000..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczNkL2JvdXNzL3Rlc3RfbnMzZF9ib3Vzcy5weQ==
--- /dev/null
+++ b/fluidsim/solvers/ns3d/bouss/test_ns3d_bouss.py
@@ -0,0 +1,100 @@
+import unittest
+
+import matplotlib.pyplot as plt
+
+from fluiddyn.io import stdout_redirected
+import fluiddyn.util.mpi as mpi
+
+import fluidsim as fls
+
+from ..test_solver import TestSolverBase
+
+from .solver import Simul
+
+
+class TestSolverNS3DBouss(TestSolverBase):
+    Simul = Simul
+
+    def test_output(self):
+
+        params = self.get_params()
+
+        params.init_fields.type = "dipole"
+
+        # save all outputs!
+        periods = params.output.periods_save
+        for key in periods._key_attribs:
+            periods[key] = 0.2
+
+        with stdout_redirected():
+            self.sim = sim = self.Simul(params)
+            sim.time_stepping.start()
+
+            if mpi.nb_proc == 1:
+
+                sim2 = fls.load_sim_for_plot(sim.output.path_run)
+                sim2.output.print_stdout.load()
+                sim2.output.print_stdout.plot()
+                sim2.output.spatial_means.load()
+                sim2.output.spatial_means.load_dataset()
+                sim2.output.spatial_means.plot()
+                sim2.output.spectra.load1d_mean()
+                sim2.output.spectra.load3d_mean()
+                sim2.output.spectra.plot1d(
+                    tmin=0.1,
+                    tmax=10,
+                    delta_t=0.01,
+                    coef_compensate=5 / 3,
+                    coef_plot_k3=1.0,
+                    coef_plot_k53=1.0,
+                )
+                sim2.output.spectra.plot3d(
+                    tmin=0.1,
+                    tmax=10,
+                    delta_t=0.01,
+                    coef_compensate=5 / 3,
+                    coef_plot_k3=1.0,
+                    coef_plot_k53=1.0,
+                )
+
+                sim2.output.phys_fields.set_equation_crosssection(
+                    f"x={sim.oper.Lx/4}"
+                )
+                sim2.output.phys_fields.animate("vx")
+
+                sim2.output.phys_fields.plot(field="vx", time=10)
+
+        plt.close("all")
+
+    def test_init_in_script(self):
+
+        params = self.get_params()
+
+        params.init_fields.type = "in_script"
+
+        with stdout_redirected():
+            self.sim = sim = Simul(params)
+
+        # here we have to initialize the flow fields
+
+        variables = {
+            k: 1e-2 * sim.oper.create_arrayX_random() for k in ("vx", "vy", "vz")
+        }
+
+        X, Y, Z = sim.oper.get_XYZ_loc()
+
+        sim.state.init_statephys_from(**variables)
+
+        sim.state.statespect_from_statephys()
+        sim.state.statephys_from_statespect()
+
+        vx_fft = sim.state.get_var("vx_fft")
+        vy_fft = sim.state.get_var("vy_fft")
+        vz_fft = sim.state.get_var("vz_fft")
+
+        sim.state.init_from_vxvyfft(vx_fft, vy_fft)
+        sim.state.init_from_vxvyvzfft(vx_fft, vy_fft, vz_fft)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/fluidsim/solvers/ns3d/init_fields.py b/fluidsim/solvers/ns3d/init_fields.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczNkL2luaXRfZmllbGRzLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczNkL2luaXRfZmllbGRzLnB5 100644
--- a/fluidsim/solvers/ns3d/init_fields.py
+++ b/fluidsim/solvers/ns3d/init_fields.py
@@ -11,5 +11,4 @@
    :members:
 
 """
-from __future__ import division
 
@@ -15,5 +14,4 @@
 
-from builtins import range
 import numpy as np
 
 from fluiddyn.util import mpi
diff --git a/fluidsim/solvers/ns3d/output/__init__.py b/fluidsim/solvers/ns3d/output/__init__.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczNkL291dHB1dC9fX2luaXRfXy5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczNkL291dHB1dC9fX2luaXRfXy5weQ== 100644
--- a/fluidsim/solvers/ns3d/output/__init__.py
+++ b/fluidsim/solvers/ns3d/output/__init__.py
@@ -1,6 +1,3 @@
-from __future__ import division
-
-
 import numpy as np
 
 from fluidsim.base.output import OutputBasePseudoSpectral
@@ -68,6 +65,5 @@
         return nrj_x_fft + nrj_y_fft + nrj_z_fft
 
     def compute_energy(self):
-        print("compute_energy, ns3d")
         energy_fft = self.compute_energy_fft()
         return self.sum_wavenumbers(energy_fft)
diff --git a/fluidsim/solvers/ns3d/output/print_stdout.py b/fluidsim/solvers/ns3d/output/print_stdout.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczNkL291dHB1dC9wcmludF9zdGRvdXQucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczNkL291dHB1dC9wcmludF9zdGRvdXQucHk= 100644
--- a/fluidsim/solvers/ns3d/output/print_stdout.py
+++ b/fluidsim/solvers/ns3d/output/print_stdout.py
@@ -1,6 +1,3 @@
-from __future__ import print_function, division
-
-from builtins import range
 import numpy as np
 
 from fluidsim.base.output.print_stdout import PrintStdOutBase
@@ -99,5 +96,5 @@
         ax1.set_title(
             "info stdout, solver "
             + self.output.name_solver
-            + ", nh = {0:5d}".format(self.nx)
+            + ", nx = {0:5d}".format(self.params.oper.nx)
         )
@@ -103,8 +100,7 @@
         )
-        ax1.hold(True)
         ax1.plot(t, deltat, "k", linewidth=2)
 
         size_axe[1] = 0.08
         ax2 = fig.add_axes(size_axe)
         ax2.set_xlabel("t")
         ax2.set_ylabel("E(t), deltaE(t)")
@@ -105,9 +101,8 @@
         ax1.plot(t, deltat, "k", linewidth=2)
 
         size_axe[1] = 0.08
         ax2 = fig.add_axes(size_axe)
         ax2.set_xlabel("t")
         ax2.set_ylabel("E(t), deltaE(t)")
-        ax2.hold(True)
         ax2.plot(t, E, "k", linewidth=2)
         ax2.plot(t, deltaE, "b", linewidth=2)
diff --git a/fluidsim/solvers/ns3d/output/spatial_means.py b/fluidsim/solvers/ns3d/output/spatial_means.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczNkL291dHB1dC9zcGF0aWFsX21lYW5zLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczNkL291dHB1dC9zcGF0aWFsX21lYW5zLnB5 100644
--- a/fluidsim/solvers/ns3d/output/spatial_means.py
+++ b/fluidsim/solvers/ns3d/output/spatial_means.py
@@ -7,9 +7,6 @@
 
 """
 
-from __future__ import division, print_function
-
-from builtins import range
 import os
 import numpy as np
 
@@ -13,7 +10,6 @@
 import os
 import numpy as np
 
-
 from fluiddyn.util import mpi
 
 from fluidsim.base.output.spatial_means import SpatialMeansBase
@@ -48,10 +44,12 @@
             vy_fft = self.sim.state.state_spect.get_var("vy_fft")
             vz_fft = self.sim.state.state_spect.get_var("vz_fft")
 
-            PK1_fft = np.real(
-                vx_fft.conj() * fx_fft
-                + vy_fft.conj() * fy_fft
-                + vz_fft.conj() * fz_fft
+            PK1_fft = np.ascontiguousarray(
+                np.real(
+                    vx_fft.conj() * fx_fft
+                    + vy_fft.conj() * fy_fft
+                    + vz_fft.conj() * fz_fft
+                )
             )
             PK2_fft = (
                 (abs(fx_fft) ** 2 + abs(fy_fft) ** 2 + abs(fz_fft) ** 2)
@@ -168,6 +166,9 @@
 
         dict_results["t"] = t
         dict_results["E"] = E
+        dict_results["Ex"] = Ex
+        dict_results["Ey"] = Ey
+        dict_results["Ez"] = Ez
 
         dict_results["PK1"] = PK1
         dict_results["PK2"] = PK2
@@ -192,11 +193,5 @@
         epsK_hypo = dict_results["epsK_hypo"]
         epsK_tot = dict_results["epsK_tot"]
 
-        width_axe = 0.85
-        height_axe = 0.39
-        x_left_axe = 0.12
-        z_bottom_axe = 0.55
-
-        size_axe = [x_left_axe, z_bottom_axe, width_axe, height_axe]
-        fig, ax1 = self.output.figure_axe(size_axe=size_axe)
+        fig, ax = self.output.figure_axe()
         fig.suptitle("Energy and enstrophy")
@@ -202,7 +197,7 @@
         fig.suptitle("Energy and enstrophy")
-        ax1.set_ylabel("$E(t)$")
-        ax1.plot(t, E, "k", linewidth=2)
-        ax1.plot(t, Ex, "b")
-        ax1.plot(t, Ey, "r")
-        ax1.plot(t, Ez, "c")
+        ax.set_ylabel("$E(t)$")
+        ax.plot(t, E, "k", linewidth=2)
+        ax.plot(t, Ex, "b")
+        ax.plot(t, Ey, "r")
+        ax.plot(t, Ez, "c")
 
@@ -208,5 +203,3 @@
 
-        z_bottom_axe = 0.54
-        size_axe[1] = z_bottom_axe
-        fig, ax1 = self.output.figure_axe(size_axe=size_axe)
+        fig, ax = self.output.figure_axe()
         fig.suptitle("Dissipation of energy and enstrophy")
@@ -212,3 +205,3 @@
         fig.suptitle("Dissipation of energy and enstrophy")
-        ax1.set_ylabel(r"$\epsilon_K(t)$")
+        ax.set_ylabel(r"$\epsilon_K(t)$")
 
@@ -214,4 +207,5 @@
 
-        ax1.plot(t, epsK, "r", linewidth=2)
-        ax1.plot(t, epsK_hypo, "g", linewidth=2)
-        ax1.plot(t, epsK_tot, "k", linewidth=2)
+        ax.plot(t, epsK, "r", linewidth=2)
+        if self.sim.params.nu_m4 != 0:
+            ax.plot(t, epsK_hypo, "g", linewidth=2)
+        ax.plot(t, epsK_tot, "k", linewidth=2)
diff --git a/fluidsim/solvers/ns3d/solver.py b/fluidsim/solvers/ns3d/solver.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczNkL3NvbHZlci5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczNkL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/ns3d/solver.py
+++ b/fluidsim/solvers/ns3d/solver.py
@@ -9,7 +9,6 @@
 
 
 """
-from __future__ import division
 
 import sys
 
diff --git a/fluidsim/solvers/ns3d/strat/output/spatial_means.py b/fluidsim/solvers/ns3d/strat/output/spatial_means.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczNkL3N0cmF0L291dHB1dC9zcGF0aWFsX21lYW5zLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczNkL3N0cmF0L291dHB1dC9zcGF0aWFsX21lYW5zLnB5 100644
--- a/fluidsim/solvers/ns3d/strat/output/spatial_means.py
+++ b/fluidsim/solvers/ns3d/strat/output/spatial_means.py
@@ -7,9 +7,6 @@
 
 """
 
-from __future__ import division, print_function
-
-from builtins import range
 import os
 import numpy as np
 
diff --git a/fluidsim/solvers/ns3d/strat/solver.py b/fluidsim/solvers/ns3d/strat/solver.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczNkL3N0cmF0L3NvbHZlci5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczNkL3N0cmF0L3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/ns3d/strat/solver.py
+++ b/fluidsim/solvers/ns3d/strat/solver.py
@@ -12,7 +12,6 @@
    :private-members:
 
 """
-from __future__ import division
 
 from fluidfft.fft3d.operators import vector_product
 
diff --git a/fluidsim/solvers/ns3d/strat/test_ns3d_strat.py b/fluidsim/solvers/ns3d/strat/test_ns3d_strat.py
new file mode 100644
index 0000000000000000000000000000000000000000..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczNkL3N0cmF0L3Rlc3RfbnMzZF9zdHJhdC5weQ==
--- /dev/null
+++ b/fluidsim/solvers/ns3d/strat/test_ns3d_strat.py
@@ -0,0 +1,146 @@
+import unittest
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from fluiddyn.io import stdout_redirected
+import fluiddyn.util.mpi as mpi
+
+import fluidsim as fls
+
+from ..test_solver import TestSolverBase
+
+from .solver import Simul
+
+
+class TestSolverNS3DStrat(TestSolverBase):
+    Simul = Simul
+
+    def test_tendency(self):
+
+        params = self.get_params()
+
+        params.init_fields.type = "noise"
+        params.output.HAS_TO_SAVE = False
+
+        with stdout_redirected():
+            self.sim = sim = self.Simul(params)
+
+        tend = sim.tendencies_nonlin(state_spect=sim.state.state_spect)
+
+        T_tot = np.ascontiguousarray(sim.oper.create_arrayK(value=0.0).real)
+
+        for axis in ("x", "y", "z"):
+            key = f"v{axis}_fft"
+            vi_fft = sim.state.get_var(key)
+            Fi_fft = tend.get_var(key)
+            T_tot += np.real(Fi_fft.conj() * vi_fft)
+
+        key = "b_fft"
+        b_fft = sim.state.get_var(key)
+        F_fft = tend.get_var(key)
+
+        T_tot += 1.0 / sim.params.N ** 2 * np.real(F_fft.conj() * b_fft)
+
+        ratio = sim.oper.sum_wavenumbers(T_tot) / sim.oper.sum_wavenumbers(
+            abs(T_tot)
+        )
+
+        self.assertGreater(1e-15, abs(ratio))
+
+    def test_output(self):
+
+        params = self.get_params()
+
+        params.init_fields.type = "dipole"
+
+        params.forcing.enable = True
+        params.forcing.type = "in_script"
+
+        # save all outputs!
+        periods = params.output.periods_save
+        for key in periods._key_attribs:
+            periods[key] = 0.2
+
+        with stdout_redirected():
+            self.sim = sim = self.Simul(params)
+
+            oper = sim.oper
+            X, Y, Z = oper.get_XYZ_loc()
+
+            def compute_forcing_fft_each_time(self):
+                return {'vx_fft': oper.create_arrayK(value=0)}
+
+            sim.forcing.forcing_maker.monkeypatch_compute_forcing_fft_each_time(
+                compute_forcing_fft_each_time)
+
+            sim.time_stepping.start()
+
+            if mpi.nb_proc == 1:
+
+                sim2 = fls.load_sim_for_plot(sim.output.path_run)
+                sim2.output.print_stdout.load()
+                sim2.output.print_stdout.plot()
+                sim2.output.spatial_means.load()
+                sim2.output.spatial_means.load_dataset()
+                sim2.output.spatial_means.plot()
+                sim2.output.spectra.load1d_mean()
+                sim2.output.spectra.load3d_mean()
+                sim2.output.spectra.plot1d(
+                    tmin=0.1,
+                    tmax=10,
+                    delta_t=0.01,
+                    coef_compensate=5 / 3,
+                    coef_plot_k3=1.0,
+                    coef_plot_k53=1.0,
+                )
+                sim2.output.spectra.plot3d(
+                    tmin=0.1,
+                    tmax=10,
+                    delta_t=0.01,
+                    coef_compensate=5 / 3,
+                    coef_plot_k3=1.0,
+                    coef_plot_k53=1.0,
+                )
+
+                sim2.output.phys_fields.set_equation_crosssection(
+                    f"x={sim.oper.Lx/4}"
+                )
+                sim2.output.phys_fields.animate("vx")
+
+                sim2.output.phys_fields.plot(field="vx", time=10)
+
+        plt.close("all")
+
+    def test_init_in_script(self):
+
+        params = self.get_params()
+
+        params.init_fields.type = "in_script"
+
+        with stdout_redirected():
+            self.sim = sim = Simul(params)
+
+        # here we have to initialize the flow fields
+
+        variables = {
+            k: 1e-2 * sim.oper.create_arrayX_random() for k in ("vx", "vy", "vz")
+        }
+
+        X, Y, Z = sim.oper.get_XYZ_loc()
+
+        sim.state.init_statephys_from(**variables)
+
+        sim.state.statespect_from_statephys()
+        sim.state.statephys_from_statespect()
+
+        vx_fft = sim.state.get_var("vx_fft")
+        vy_fft = sim.state.get_var("vy_fft")
+        vz_fft = sim.state.get_var("vz_fft")
+
+        sim.state.init_from_vxvyfft(vx_fft, vy_fft)
+        sim.state.init_from_vxvyvzfft(vx_fft, vy_fft, vz_fft)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/fluidsim/solvers/ns3d/test_solver.py b/fluidsim/solvers/ns3d/test_solver.py
new file mode 100644
index 0000000000000000000000000000000000000000..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczNkL3Rlc3Rfc29sdmVyLnB5
--- /dev/null
+++ b/fluidsim/solvers/ns3d/test_solver.py
@@ -0,0 +1,189 @@
+import unittest
+import shutil
+import sys
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import fluidsim as fls
+
+import fluiddyn.util.mpi as mpi
+from fluiddyn.io import stdout_redirected
+
+from fluidsim.solvers.ns3d.solver import Simul
+
+from fluidsim.base.output import run
+
+
+class TestSolverBase(unittest.TestCase):
+    Simul = Simul
+
+    def tearDown(self):
+        # clean by removing the directory
+        if mpi.rank == 0:
+            if hasattr(self, "sim"):
+                shutil.rmtree(self.sim.output.path_run)
+
+    def get_params(self):
+
+        params = self.Simul.create_default_params()
+
+        params.short_name_type_run = "test"
+
+        nx = 32
+        params.oper.nx = nx
+        params.oper.ny = nx
+        params.oper.nz = nx
+
+        Lx = 6.0
+        params.oper.Lx = Lx
+        params.oper.Ly = Lx
+        params.oper.Lz = Lx
+
+        params.oper.coef_dealiasing = 2.0 / 3
+        params.nu_8 = 2.0
+
+        params.time_stepping.t_end = 0.5
+
+        return params
+
+
+class TestSolverNS3D(TestSolverBase):
+    def test_tendency(self):
+
+        params = self.get_params()
+
+        params.init_fields.type = "noise"
+        params.output.HAS_TO_SAVE = False
+
+        with stdout_redirected():
+            self.sim = sim = self.Simul(params)
+
+        tend = sim.tendencies_nonlin(state_spect=sim.state.state_spect)
+
+        T_tot = np.ascontiguousarray(sim.oper.create_arrayK(value=0.0).real)
+
+        for axis in ("x", "y", "z"):
+            key = f"v{axis}_fft"
+            vi_fft = sim.state.get_var(key)
+            Fi_fft = tend.get_var(key)
+            T_tot += np.real(Fi_fft.conj() * vi_fft)
+
+        ratio = sim.oper.sum_wavenumbers(T_tot) / sim.oper.sum_wavenumbers(
+            abs(T_tot)
+        )
+
+        self.assertGreater(1e-15, abs(ratio))
+
+    def test_output(self):
+
+        params = self.get_params()
+
+        params.init_fields.type = "dipole"
+
+        params.forcing.enable = True
+        params.forcing.type = "in_script"
+        params.forcing.key_forced = "vx_fft"
+
+        # save all outputs!
+        periods = params.output.periods_save
+        for key in periods._key_attribs:
+            periods[key] = 0.2
+
+        with stdout_redirected():
+            self.sim = sim = self.Simul(params)
+            sim.time_stepping.start()
+
+            # compute twice for better coverage
+            sim.state.compute("rotz")
+            sim.state.compute("rotz")
+
+            if mpi.nb_proc == 1:
+
+                sim2 = fls.load_sim_for_plot(sim.output.path_run)
+                sim2.output.print_stdout.load()
+                sim2.output.print_stdout.plot()
+                sim2.output.spatial_means.load()
+                sim2.output.spatial_means.load_dataset()
+                sim2.output.spatial_means.plot()
+                sim2.output.spectra.load1d_mean()
+                sim2.output.spectra.load3d_mean()
+                sim2.output.spectra.plot1d(
+                    tmin=0.1,
+                    tmax=10,
+                    delta_t=0.01,
+                    coef_compensate=5 / 3,
+                    coef_plot_k3=1.0,
+                    coef_plot_k53=1.0,
+                )
+                sim2.output.spectra.plot3d(
+                    tmin=0.1,
+                    tmax=10,
+                    delta_t=0.01,
+                    coef_compensate=5 / 3,
+                    coef_plot_k3=1.0,
+                    coef_plot_k53=1.0,
+                )
+
+                sim2.output.phys_fields.set_equation_crosssection(
+                    f"x={sim.oper.Lx/4}"
+                )
+                sim2.output.phys_fields.animate("vx")
+
+                sim2.output.phys_fields.plot(field="vx", time=10)
+
+            path_run = sim.output.path_run
+            if mpi.nb_proc > 1:
+                path_run = mpi.comm.bcast(path_run)
+
+            sim3 = fls.load_state_phys_file(path_run, modif_save_params=False)
+            sim3.params.time_stepping.t_end += 0.2
+            sim3.time_stepping.start()
+
+            if mpi.nb_proc == 1:
+                sim3.output.phys_fields.animate(
+                    "vx",
+                    dt_frame_in_sec=1e-6,
+                    dt_equations=0.3,
+                    repeat=False,
+                    clim=(-1, 1),
+                    save_file=False,
+                    numfig=1,
+                )
+
+                sys.argv = ["fluidsim-create-xml-description", path_run]
+                run()
+        plt.close("all")
+
+    def test_init_in_script(self):
+
+        params = self.get_params()
+
+        params.init_fields.type = "in_script"
+
+        with stdout_redirected():
+            self.sim = sim = Simul(params)
+
+        # here we have to initialize the flow fields
+
+        variables = {
+            k: 1e-2 * sim.oper.create_arrayX_random() for k in ("vx", "vy", "vz")
+        }
+
+        X, Y, Z = sim.oper.get_XYZ_loc()
+
+        sim.state.init_statephys_from(**variables)
+
+        sim.state.statespect_from_statephys()
+        sim.state.statephys_from_statespect()
+
+        vx_fft = sim.state.get_var("vx_fft")
+        vy_fft = sim.state.get_var("vy_fft")
+        vz_fft = sim.state.get_var("vz_fft")
+
+        sim.state.init_from_vxvyfft(vx_fft, vy_fft)
+        sim.state.init_from_vxvyvzfft(vx_fft, vy_fft, vz_fft)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/fluidsim/solvers/ns3d/try_load.py b/fluidsim/solvers/ns3d/try_load.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczNkL3RyeV9sb2FkLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczNkL3RyeV9sb2FkLnB5 100644
--- a/fluidsim/solvers/ns3d/try_load.py
+++ b/fluidsim/solvers/ns3d/try_load.py
@@ -1,4 +1,3 @@
-from __future__ import division
 from past.utils import old_div
 from fluidsim.solvers.ns3d.solver import Simul
 import numpy as np
diff --git a/fluidsim/solvers/ns3d/try_save.py b/fluidsim/solvers/ns3d/try_save.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9uczNkL3RyeV9zYXZlLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9uczNkL3RyeV9zYXZlLnB5 100644
--- a/fluidsim/solvers/ns3d/try_save.py
+++ b/fluidsim/solvers/ns3d/try_save.py
@@ -1,4 +1,3 @@
-from __future__ import division
 from past.utils import old_div
 from fluidsim.solvers.ns3d.solver import Simul
 import numpy as np
diff --git a/fluidsim/solvers/plate2d/diag/solver.py b/fluidsim/solvers/plate2d/diag/solver.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL2RpYWcvc29sdmVyLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL2RpYWcvc29sdmVyLnB5 100644
--- a/fluidsim/solvers/plate2d/diag/solver.py
+++ b/fluidsim/solvers/plate2d/diag/solver.py
@@ -12,9 +12,6 @@
 
 """
 
-from __future__ import print_function
-from __future__ import division
-
 from past.utils import old_div
 import numpy as np
 
diff --git a/fluidsim/solvers/plate2d/dimensional.py b/fluidsim/solvers/plate2d/dimensional.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL2RpbWVuc2lvbmFsLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL2RpbWVuc2lvbmFsLnB5 100644
--- a/fluidsim/solvers/plate2d/dimensional.py
+++ b/fluidsim/solvers/plate2d/dimensional.py
@@ -12,10 +12,6 @@
 
 """
 
-from __future__ import print_function
-from __future__ import division
-
-from builtins import object
 from past.utils import old_div
 import numpy as np
 
diff --git a/fluidsim/solvers/plate2d/init_fields.py b/fluidsim/solvers/plate2d/init_fields.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL2luaXRfZmllbGRzLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL2luaXRfZmllbGRzLnB5 100644
--- a/fluidsim/solvers/plate2d/init_fields.py
+++ b/fluidsim/solvers/plate2d/init_fields.py
@@ -4,7 +4,6 @@
 
 
 """
-from __future__ import division
 
 import numpy as np
 
diff --git a/fluidsim/solvers/plate2d/output/spatial_means.py b/fluidsim/solvers/plate2d/output/spatial_means.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL291dHB1dC9zcGF0aWFsX21lYW5zLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL291dHB1dC9zcGF0aWFsX21lYW5zLnB5 100644
--- a/fluidsim/solvers/plate2d/output/spatial_means.py
+++ b/fluidsim/solvers/plate2d/output/spatial_means.py
@@ -10,5 +10,4 @@
    :private-members:
 
 """
-from __future__ import division, print_function
 
@@ -14,5 +13,4 @@
 
-from builtins import range
 import os
 import numpy as np
 
diff --git a/fluidsim/solvers/plate2d/output/spectra.py b/fluidsim/solvers/plate2d/output/spectra.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL291dHB1dC9zcGVjdHJhLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL291dHB1dC9zcGVjdHJhLnB5 100644
--- a/fluidsim/solvers/plate2d/output/spectra.py
+++ b/fluidsim/solvers/plate2d/output/spectra.py
@@ -9,6 +9,4 @@
    :private-members:
 
 """
-from __future__ import division
-from __future__ import print_function
 
@@ -14,5 +12,4 @@
 
-from builtins import range
 from past.utils import old_div
 import h5py
 
diff --git a/fluidsim/solvers/plate2d/solver.py b/fluidsim/solvers/plate2d/solver.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL3NvbHZlci5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/plate2d/solver.py
+++ b/fluidsim/solvers/plate2d/solver.py
@@ -16,9 +16,6 @@
 
 """
 
-from __future__ import print_function
-from __future__ import division
-
 import numpy as np
 
 from fluidsim.base.setofvariables import SetOfVariables
diff --git a/fluidsim/solvers/plate2d/state.py b/fluidsim/solvers/plate2d/state.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL3N0YXRlLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL3N0YXRlLnB5 100644
--- a/fluidsim/solvers/plate2d/state.py
+++ b/fluidsim/solvers/plate2d/state.py
@@ -1,7 +1,6 @@
 """Plate2d state (:mod:`fluidsim.solvers.plate2d.state`)
 ==============================================================
 """
-from __future__ import print_function
 
 import numpy as np
 
diff --git a/fluidsim/solvers/plate2d/test/test_solver.py b/fluidsim/solvers/plate2d/test/test_solver.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL3Rlc3QvdGVzdF9zb2x2ZXIucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL3Rlc3QvdGVzdF9zb2x2ZXIucHk= 100644
--- a/fluidsim/solvers/plate2d/test/test_solver.py
+++ b/fluidsim/solvers/plate2d/test/test_solver.py
@@ -1,5 +1,3 @@
-from __future__ import division
-
 import unittest
 import shutil
 
diff --git a/fluidsim/solvers/sw1l/exactlin/modified/solver.py b/fluidsim/solvers/sw1l/exactlin/modified/solver.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2V4YWN0bGluL21vZGlmaWVkL3NvbHZlci5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2V4YWN0bGluL21vZGlmaWVkL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/sw1l/exactlin/modified/solver.py
+++ b/fluidsim/solvers/sw1l/exactlin/modified/solver.py
@@ -7,8 +7,6 @@
 
 """
 
-from __future__ import division, print_function
-
 import numpy as np
 
 from fluidsim.base.setofvariables import SetOfVariables
diff --git a/fluidsim/solvers/sw1l/exactlin/solver.py b/fluidsim/solvers/sw1l/exactlin/solver.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2V4YWN0bGluL3NvbHZlci5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2V4YWN0bGluL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/sw1l/exactlin/solver.py
+++ b/fluidsim/solvers/sw1l/exactlin/solver.py
@@ -7,8 +7,6 @@
 
 """
 
-from __future__ import division, print_function
-
 import numpy as np
 
 
diff --git a/fluidsim/solvers/sw1l/forcing.py b/fluidsim/solvers/sw1l/forcing.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2ZvcmNpbmcucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2ZvcmNpbmcucHk= 100644
--- a/fluidsim/solvers/sw1l/forcing.py
+++ b/fluidsim/solvers/sw1l/forcing.py
@@ -3,6 +3,4 @@
 ===================================================
 
 """
-from __future__ import division
-from __future__ import print_function
 
@@ -8,5 +6,4 @@
 
-from builtins import object
 import numpy as np
 
 from fluiddyn.util import mpi
diff --git a/fluidsim/solvers/sw1l/init_fields.py b/fluidsim/solvers/sw1l/init_fields.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2luaXRfZmllbGRzLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2luaXRfZmllbGRzLnB5 100644
--- a/fluidsim/solvers/sw1l/init_fields.py
+++ b/fluidsim/solvers/sw1l/init_fields.py
@@ -20,11 +20,6 @@
    :private-members:
 
 """
-from __future__ import division
-from __future__ import print_function
-
-
-from builtins import range
 
 import numpy as np
 
diff --git a/fluidsim/solvers/sw1l/modified/output.py b/fluidsim/solvers/sw1l/modified/output.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL21vZGlmaWVkL291dHB1dC5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9zdzFsL21vZGlmaWVkL291dHB1dC5weQ== 100644
--- a/fluidsim/solvers/sw1l/modified/output.py
+++ b/fluidsim/solvers/sw1l/modified/output.py
@@ -1,5 +1,4 @@
 """ """
-from __future__ import division
 
 import numpy as np
 
diff --git a/fluidsim/solvers/sw1l/modified/solver.py b/fluidsim/solvers/sw1l/modified/solver.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL21vZGlmaWVkL3NvbHZlci5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9zdzFsL21vZGlmaWVkL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/sw1l/modified/solver.py
+++ b/fluidsim/solvers/sw1l/modified/solver.py
@@ -8,8 +8,6 @@
 due to the rotational velocity.
 """
 
-from __future__ import division, print_function
-
 import numpy as np
 
 from fluidsim.base.setofvariables import SetOfVariables
diff --git a/fluidsim/solvers/sw1l/modified/state.py b/fluidsim/solvers/sw1l/modified/state.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL21vZGlmaWVkL3N0YXRlLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9zdzFsL21vZGlmaWVkL3N0YXRlLnB5 100644
--- a/fluidsim/solvers/sw1l/modified/state.py
+++ b/fluidsim/solvers/sw1l/modified/state.py
@@ -10,7 +10,6 @@
    :private-members:
 
 """
-from __future__ import print_function
 
 from fluidsim.base.setofvariables import SetOfVariables
 
diff --git a/fluidsim/solvers/sw1l/onlywaves/solver.py b/fluidsim/solvers/sw1l/onlywaves/solver.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL29ubHl3YXZlcy9zb2x2ZXIucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9zdzFsL29ubHl3YXZlcy9zb2x2ZXIucHk= 100644
--- a/fluidsim/solvers/sw1l/onlywaves/solver.py
+++ b/fluidsim/solvers/sw1l/onlywaves/solver.py
@@ -7,8 +7,6 @@
 equations with zeros QG PV.
 """
 
-from __future__ import division, print_function
-
 import numpy as np
 
 from fluidsim.base.setofvariables import SetOfVariables
diff --git a/fluidsim/solvers/sw1l/output/print_stdout.py b/fluidsim/solvers/sw1l/output/print_stdout.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9wcmludF9zdGRvdXQucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9wcmludF9zdGRvdXQucHk= 100644
--- a/fluidsim/solvers/sw1l/output/print_stdout.py
+++ b/fluidsim/solvers/sw1l/output/print_stdout.py
@@ -1,6 +1,3 @@
-from __future__ import print_function, division
-
-from builtins import range
 import numpy as np
 
 from fluidsim.base.output.print_stdout import PrintStdOutBase
diff --git a/fluidsim/solvers/sw1l/output/spatial_means.py b/fluidsim/solvers/sw1l/output/spatial_means.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9zcGF0aWFsX21lYW5zLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9zcGF0aWFsX21lYW5zLnB5 100644
--- a/fluidsim/solvers/sw1l/output/spatial_means.py
+++ b/fluidsim/solvers/sw1l/output/spatial_means.py
@@ -1,5 +1,3 @@
-from __future__ import division
-
 import os
 
 import numpy as np
diff --git a/fluidsim/solvers/sw1l/output/spect_energy_budget.py b/fluidsim/solvers/sw1l/output/spect_energy_budget.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9zcGVjdF9lbmVyZ3lfYnVkZ2V0LnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9zcGVjdF9lbmVyZ3lfYnVkZ2V0LnB5 100644
--- a/fluidsim/solvers/sw1l/output/spect_energy_budget.py
+++ b/fluidsim/solvers/sw1l/output/spect_energy_budget.py
@@ -1,7 +1,3 @@
-from __future__ import division
-from __future__ import print_function
-
-from builtins import range
 import numpy as np
 import h5py
 
diff --git a/fluidsim/solvers/sw1l/output/test/__init__.py b/fluidsim/solvers/sw1l/output/test/__init__.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC90ZXN0L19faW5pdF9fLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC90ZXN0L19faW5pdF9fLnB5 100644
--- a/fluidsim/solvers/sw1l/output/test/__init__.py
+++ b/fluidsim/solvers/sw1l/output/test/__init__.py
@@ -1,5 +1,3 @@
-from __future__ import print_function
-
 import unittest
 from shutil import rmtree
 
diff --git a/fluidsim/solvers/sw1l/output/test/test_phys_fields.py b/fluidsim/solvers/sw1l/output/test/test_phys_fields.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC90ZXN0L3Rlc3RfcGh5c19maWVsZHMucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC90ZXN0L3Rlc3RfcGh5c19maWVsZHMucHk= 100644
--- a/fluidsim/solvers/sw1l/output/test/test_phys_fields.py
+++ b/fluidsim/solvers/sw1l/output/test/test_phys_fields.py
@@ -1,5 +1,3 @@
-from __future__ import print_function
-
 import unittest
 
 from fluidsim.solvers.test.test_ns import run_mini_simul
diff --git a/fluidsim/solvers/sw1l/output/test/test_print_stdout.py b/fluidsim/solvers/sw1l/output/test/test_print_stdout.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC90ZXN0L3Rlc3RfcHJpbnRfc3Rkb3V0LnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC90ZXN0L3Rlc3RfcHJpbnRfc3Rkb3V0LnB5 100644
--- a/fluidsim/solvers/sw1l/output/test/test_print_stdout.py
+++ b/fluidsim/solvers/sw1l/output/test/test_print_stdout.py
@@ -1,5 +1,3 @@
-from __future__ import print_function
-
 import unittest
 
 import numpy as np
diff --git a/fluidsim/solvers/sw1l/output/test/test_spatial_means.py b/fluidsim/solvers/sw1l/output/test/test_spatial_means.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC90ZXN0L3Rlc3Rfc3BhdGlhbF9tZWFucy5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC90ZXN0L3Rlc3Rfc3BhdGlhbF9tZWFucy5weQ== 100644
--- a/fluidsim/solvers/sw1l/output/test/test_spatial_means.py
+++ b/fluidsim/solvers/sw1l/output/test/test_spatial_means.py
@@ -1,5 +1,3 @@
-from __future__ import print_function
-
 import unittest
 
 from fluidsim.solvers.test.test_ns import run_mini_simul
diff --git a/fluidsim/solvers/sw1l/output/test/test_spect_energy_budg.py b/fluidsim/solvers/sw1l/output/test/test_spect_energy_budg.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC90ZXN0L3Rlc3Rfc3BlY3RfZW5lcmd5X2J1ZGcucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC90ZXN0L3Rlc3Rfc3BlY3RfZW5lcmd5X2J1ZGcucHk= 100644
--- a/fluidsim/solvers/sw1l/output/test/test_spect_energy_budg.py
+++ b/fluidsim/solvers/sw1l/output/test/test_spect_energy_budg.py
@@ -1,5 +1,3 @@
-from __future__ import print_function
-
 import unittest
 import numpy as np
 
diff --git a/fluidsim/solvers/sw1l/output/test/test_spectra.py b/fluidsim/solvers/sw1l/output/test/test_spectra.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC90ZXN0L3Rlc3Rfc3BlY3RyYS5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC90ZXN0L3Rlc3Rfc3BlY3RyYS5weQ== 100644
--- a/fluidsim/solvers/sw1l/output/test/test_spectra.py
+++ b/fluidsim/solvers/sw1l/output/test/test_spectra.py
@@ -1,5 +1,3 @@
-from __future__ import print_function
-
 import unittest
 
 from fluidsim.solvers.sw1l.output.test import BaseTestCase, mpi
diff --git a/fluidsim/solvers/sw1l/test_operators.py b/fluidsim/solvers/sw1l/test_operators.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy9zdzFsL3Rlc3Rfb3BlcmF0b3JzLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy9zdzFsL3Rlc3Rfb3BlcmF0b3JzLnB5 100644
--- a/fluidsim/solvers/sw1l/test_operators.py
+++ b/fluidsim/solvers/sw1l/test_operators.py
@@ -1,5 +1,3 @@
-from __future__ import division
-
 import unittest
 import numpy as np
 import sys
diff --git a/fluidsim/solvers/test/test_ns.py b/fluidsim/solvers/test/test_ns.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy90ZXN0L3Rlc3RfbnMucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy90ZXN0L3Rlc3RfbnMucHk= 100644
--- a/fluidsim/solvers/test/test_ns.py
+++ b/fluidsim/solvers/test/test_ns.py
@@ -1,5 +1,3 @@
-from __future__ import division
-
 import unittest
 from shutil import rmtree
 
diff --git a/fluidsim/solvers/waves2d/state.py b/fluidsim/solvers/waves2d/state.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vc29sdmVycy93YXZlczJkL3N0YXRlLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vc29sdmVycy93YXZlczJkL3N0YXRlLnB5 100644
--- a/fluidsim/solvers/waves2d/state.py
+++ b/fluidsim/solvers/waves2d/state.py
@@ -1,8 +1,6 @@
 """State for the Waves solver (:mod:`fluidsim.solvers.waves2d.state`)
 =======================================================================
 """
-from __future__ import print_function
-
 
 from fluidsim.base.state import StatePseudoSpectral
 
diff --git a/fluidsim/util/console/bench.py b/fluidsim/util/console/bench.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vdXRpbC9jb25zb2xlL2JlbmNoLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vdXRpbC9jb25zb2xlL2JlbmNoLnB5 100644
--- a/fluidsim/util/console/bench.py
+++ b/fluidsim/util/console/bench.py
@@ -2,7 +2,6 @@
 ======================================================
 
 """
-from __future__ import print_function, division
 
 import gc
 from collections import OrderedDict
diff --git a/fluidsim/util/console/profile.py b/fluidsim/util/console/profile.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vdXRpbC9jb25zb2xlL3Byb2ZpbGUucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vdXRpbC9jb25zb2xlL3Byb2ZpbGUucHk= 100644
--- a/fluidsim/util/console/profile.py
+++ b/fluidsim/util/console/profile.py
@@ -2,7 +2,6 @@
 ======================================================
 
 """
-from __future__ import print_function, division
 
 import gc
 import os
diff --git a/fluidsim/util/console/util.py b/fluidsim/util/console/util.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vdXRpbC9jb25zb2xlL3V0aWwucHk=..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vdXRpbC9jb25zb2xlL3V0aWwucHk= 100644
--- a/fluidsim/util/console/util.py
+++ b/fluidsim/util/console/util.py
@@ -2,8 +2,6 @@
 =============================================================================
 
 """
-from __future__ import print_function
-from __future__ import division
 import os
 from time import time
 
diff --git a/fluidsim/util/util.py b/fluidsim/util/util.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_Zmx1aWRzaW0vdXRpbC91dGlsLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_Zmx1aWRzaW0vdXRpbC91dGlsLnB5 100644
--- a/fluidsim/util/util.py
+++ b/fluidsim/util/util.py
@@ -246,7 +246,6 @@
         pass
 
     fix_old_params(params)
-    params.forcing.enable = False
 
     sim = solver.Simul(params)
     return sim
diff --git a/scripts/bugs/no_energy_conservation_ns2dstrat.py b/scripts/bugs/no_energy_conservation_ns2dstrat.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_c2NyaXB0cy9idWdzL25vX2VuZXJneV9jb25zZXJ2YXRpb25fbnMyZHN0cmF0LnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_c2NyaXB0cy9idWdzL25vX2VuZXJneV9jb25zZXJ2YXRpb25fbnMyZHN0cmF0LnB5 100644
--- a/scripts/bugs/no_energy_conservation_ns2dstrat.py
+++ b/scripts/bugs/no_energy_conservation_ns2dstrat.py
@@ -21,7 +21,6 @@
 python no_energy_conservation_ns2dstrat.py
 
 """
-from __future__ import print_function
 
 # import numpy as np
 import matplotlib.pyplot as plt
diff --git a/scripts/bugs/wrong_output_it0_initfield.py b/scripts/bugs/wrong_output_it0_initfield.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_c2NyaXB0cy9idWdzL3dyb25nX291dHB1dF9pdDBfaW5pdGZpZWxkLnB5..e7e3275ab7675c5fc429ecdfa056b44c05b78136_c2NyaXB0cy9idWdzL3dyb25nX291dHB1dF9pdDBfaW5pdGZpZWxkLnB5 100644
--- a/scripts/bugs/wrong_output_it0_initfield.py
+++ b/scripts/bugs/wrong_output_it0_initfield.py
@@ -19,7 +19,6 @@
 python wrong_output_it0_initfield.py
 
 """
-from __future__ import print_function
 
 import matplotlib.pyplot as plt
 
diff --git a/scripts/ns2d.strat/ns2dstrat_lmode.py b/scripts/ns2d.strat/ns2dstrat_lmode.py
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_c2NyaXB0cy9uczJkLnN0cmF0L25zMmRzdHJhdF9sbW9kZS5weQ==..e7e3275ab7675c5fc429ecdfa056b44c05b78136_c2NyaXB0cy9uczJkLnN0cmF0L25zMmRzdHJhdF9sbW9kZS5weQ== 100644
--- a/scripts/ns2d.strat/ns2dstrat_lmode.py
+++ b/scripts/ns2d.strat/ns2dstrat_lmode.py
@@ -4,8 +4,6 @@
 
 """
 
-from __future__ import print_function
-
 import numpy as np
 
 from math import pi
diff --git a/setup.cfg b/setup.cfg
index d9f8ed29f9f9e210d703baa5fdd56c7505f86c44_c2V0dXAuY2Zn..e7e3275ab7675c5fc429ecdfa056b44c05b78136_c2V0dXAuY2Zn 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -4,7 +4,10 @@
 [coverage:run]
 source = ./fluidsim
 data_file = .coverage/coverage
-omit = */test*.py
+omit =
+    */test*.py
+    */try_*.py
+    */_old_*.py
 plugins = Cython.Coverage
 
 [coverage:report]
@@ -12,6 +15,11 @@
 exclude_lines =
     raise ValueError
     if __name__ == .__main__.:
+    if "sphinx" in sys.modules:
+    except KeyError:
+    except ImportError:
+    except AttributeError:
+    except NotImplementedError:
 
 [coverage:html]
 directory = .coverage