# HG changeset patch
# User paugier <pierre.augier@univ-grenoble-alpes.fr>
# Date 1545341748 -3600
#      Thu Dec 20 22:35:48 2018 +0100
# Node ID 1c0a407b8eb53b4bcf02fd630caf51f4992aa16f
# Parent  8968b988f0274a91270689c7d5dcecaeae2ee05a
Taylor-Green testcase + less py27 compatibility

diff --git a/.travis.yml b/.travis.yml
--- 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
--- 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
diff --git a/bench/spectraldns/bench_3d.py b/bench/spectraldns/bench_3d.py
--- 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
--- 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
--- 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
--- 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/test_cases/Taylor_Green_vortices/run_simul.py b/doc/test_cases/Taylor_Green_vortices/run_simul.py
new file mode 100644
--- /dev/null
+++ b/doc/test_cases/Taylor_Green_vortices/run_simul.py
@@ -0,0 +1,76 @@
+"""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_state_phys_file
+sim = load_state_phys_file()
+sim.output.phys_fields.set_equation_crosssection(f'x={{sim.oper.Lx/4}}')
+sim.output.phys_fields.animate('vx')
+
+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
--- a/fluidsim/__init__.py
+++ b/fluidsim/__init__.py
@@ -35,7 +35,7 @@
 
 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")
 ):
     from fluiddyn.util import mpi
 
diff --git a/fluidsim/base/basilisk/operators.py b/fluidsim/base/basilisk/operators.py
--- 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/basilisk/solver.py b/fluidsim/base/basilisk/solver.py
--- a/fluidsim/base/basilisk/solver.py
+++ b/fluidsim/base/basilisk/solver.py
@@ -12,7 +12,6 @@
    :private-members:
 
 """
-from __future__ import print_function
 
 from fluidsim.base.solvers.base import SimulBase
 from fluidsim.base.solvers.info_base import InfoSolverBase
diff --git a/fluidsim/base/basilisk/time_stepping.py b/fluidsim/base/basilisk/time_stepping.py
--- a/fluidsim/base/basilisk/time_stepping.py
+++ b/fluidsim/base/basilisk/time_stepping.py
@@ -9,7 +9,6 @@
    :private-members:
 
 """
-from __future__ import print_function
 
 from time import time
 from signal import signal
diff --git a/fluidsim/base/forcing/anisotropic.py b/fluidsim/base/forcing/anisotropic.py
--- 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/specific.py b/fluidsim/base/forcing/specific.py
--- 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
diff --git a/fluidsim/base/init_fields.py b/fluidsim/base/init_fields.py
--- a/fluidsim/base/init_fields.py
+++ b/fluidsim/base/init_fields.py
@@ -8,11 +8,7 @@
    :private-members:
 
 """
-from __future__ import division
-from __future__ import print_function
 
-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
--- 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
diff --git a/fluidsim/base/output/base.py b/fluidsim/base/output/base.py
--- 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
--- 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
diff --git a/fluidsim/base/output/phys_fields2d.py b/fluidsim/base/output/phys_fields2d.py
--- 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,10 +463,10 @@
         object corresponding to a 2D contour plot.
 
         """
-        if isinstance(vecx, basestring):
+        if isinstance(vecx, str):
             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/print_stdout.py b/fluidsim/base/output/print_stdout.py
--- 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/prob_dens_func.py b/fluidsim/base/output/prob_dens_func.py
--- a/fluidsim/base/output/prob_dens_func.py
+++ b/fluidsim/base/output/prob_dens_func.py
@@ -1,7 +1,3 @@
-from __future__ import division
-from __future__ import print_function
-
-from builtins import range
 from past.utils import old_div
 import h5py
 import numpy as np
diff --git a/fluidsim/base/output/spatial_means.py b/fluidsim/base/output/spatial_means.py
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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/base/time_stepping/base.py b/fluidsim/base/time_stepping/base.py
--- a/fluidsim/base/time_stepping/base.py
+++ b/fluidsim/base/time_stepping/base.py
@@ -8,10 +8,6 @@
    :private-members:
 
 """
-from __future__ import division
-from __future__ import print_function
-
-from builtins import object
 
 from signal import signal
 from warnings import warn
diff --git a/fluidsim/magic.py b/fluidsim/magic.py
--- 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
--- 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
--- 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/test/test_operators2d.py b/fluidsim/operators/test/test_operators2d.py
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- a/fluidsim/solvers/ns2d/init_fields.py
+++ b/fluidsim/solvers/ns2d/init_fields.py
@@ -14,9 +14,7 @@
    :members:
 
 """
-from __future__ import division
 
-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
--- 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
--- 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
--- 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
--- a/fluidsim/solvers/ns2d/output/spect_energy_budget.py
+++ b/fluidsim/solvers/ns2d/output/spect_energy_budget.py
@@ -6,10 +6,7 @@
    :private-members:
 
 """
-from __future__ import division
-from __future__ import print_function
 
-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
--- a/fluidsim/solvers/ns2d/output/spectra.py
+++ b/fluidsim/solvers/ns2d/output/spectra.py
@@ -6,10 +6,7 @@
    :private-members:
 
 """
-from __future__ import division
-from __future__ import print_function
 
-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
--- 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
--- 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
--- a/fluidsim/solvers/ns2d/strat/diag/solver.py
+++ b/fluidsim/solvers/ns2d/strat/diag/solver.py
@@ -6,9 +6,7 @@
    :private-members:
 
 """
-from __future__ import division
 
-from past.utils import old_div
 from fluidsim.base.setofvariables import SetOfVariables
 
 from fluidsim.base.solvers.pseudo_spect import (
@@ -124,12 +122,9 @@
 
     # 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 = (
-        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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
 
diff --git a/fluidsim/solvers/ns2d/strat/output/spectra_multidim.py b/fluidsim/solvers/ns2d/strat/output/spectra_multidim.py
--- a/fluidsim/solvers/ns2d/strat/output/spectra_multidim.py
+++ b/fluidsim/solvers/ns2d/strat/output/spectra_multidim.py
@@ -29,20 +29,28 @@
         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
+        )
 
-        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
 
-        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,
-            "spectrumkykx_am": spectrumkykx_am
+            "spectrumkykx_am": spectrumkykx_am,
         }
 
         return dict_spectra
@@ -103,14 +111,16 @@
             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":
-                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)
-        kz[kz.shape[0]//2] *= -1
+        kz[kz.shape[0] // 2] *= -1
 
         # Create mesh of wave-numbers
         KX, KZ = np.meshgrid(kx, kz)
@@ -136,14 +146,22 @@
 
         # 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)
-        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,12 +178,13 @@
 
         if isinstance(pforcing.tcrandom_anisotropic.angle, str):
             from math import radians
+
             angle = radians(
-                float((pforcing.tcrandom_anisotropic.angle).split("°")[0]))
+                float((pforcing.tcrandom_anisotropic.angle).split("°")[0])
+            )
         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,14 +193,26 @@
 
         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(
-                (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,10 +223,10 @@
                 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="-.",
-                color="red"
+                color="red",
             )
         )
         ax.add_patch(
@@ -205,9 +236,9 @@
                 height=2 * pforcing.nkmax_forcing * deltak,
                 angle=0,
                 theta1=-90,
-                theta2=90.,
+                theta2=90.0,
                 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
--- 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
--- 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
--- a/fluidsim/solvers/ns2d/strat/output/test/test_spatial_means.py
+++ b/fluidsim/solvers/ns2d/strat/output/test/test_spatial_means.py
@@ -1,9 +1,6 @@
-from __future__ import print_function
-
 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- a/fluidsim/solvers/ns2d/test/test_solver.py
+++ b/fluidsim/solvers/ns2d/test/test_solver.py
@@ -1,5 +1,3 @@
-from __future__ import division
-
 import unittest
 import shutil
 
diff --git a/fluidsim/solvers/ns3d/bouss/solver.py b/fluidsim/solvers/ns3d/bouss/solver.py
--- 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/init_fields.py b/fluidsim/solvers/ns3d/init_fields.py
--- a/fluidsim/solvers/ns3d/init_fields.py
+++ b/fluidsim/solvers/ns3d/init_fields.py
@@ -11,9 +11,7 @@
    :members:
 
 """
-from __future__ import division
 
-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
--- 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
--- 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
diff --git a/fluidsim/solvers/ns3d/output/spatial_means.py b/fluidsim/solvers/ns3d/output/spatial_means.py
--- a/fluidsim/solvers/ns3d/output/spatial_means.py
+++ b/fluidsim/solvers/ns3d/output/spatial_means.py
@@ -7,13 +7,9 @@
 
 """
 
-from __future__ import division, print_function
-
-from builtins import range
 import os
 import numpy as np
 
-
 from fluiddyn.util import mpi
 
 from fluidsim.base.output.spatial_means import SpatialMeansBase
@@ -168,6 +164,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,26 +191,19 @@
         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")
-        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")
 
-        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")
-        ax1.set_ylabel(r"$\epsilon_K(t)$")
+        ax.set_ylabel(r"$\epsilon_K(t)$")
 
-        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
--- 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
--- 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
--- 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/try_load.py b/fluidsim/solvers/ns3d/try_load.py
--- 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
--- 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
--- 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
--- 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
--- 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
--- a/fluidsim/solvers/plate2d/output/spatial_means.py
+++ b/fluidsim/solvers/plate2d/output/spatial_means.py
@@ -10,9 +10,7 @@
    :private-members:
 
 """
-from __future__ import division, print_function
 
-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
--- a/fluidsim/solvers/plate2d/output/spectra.py
+++ b/fluidsim/solvers/plate2d/output/spectra.py
@@ -9,10 +9,7 @@
    :private-members:
 
 """
-from __future__ import division
-from __future__ import print_function
 
-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
--- 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
--- 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
--- 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
--- 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
--- 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
--- a/fluidsim/solvers/sw1l/forcing.py
+++ b/fluidsim/solvers/sw1l/forcing.py
@@ -3,10 +3,7 @@
 ===================================================
 
 """
-from __future__ import division
-from __future__ import print_function
 
-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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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/state.py b/fluidsim/solvers/sw1l/state.py
--- a/fluidsim/solvers/sw1l/state.py
+++ b/fluidsim/solvers/sw1l/state.py
@@ -6,8 +6,6 @@
    :private-members:
 
 """
-from __future__ import division
-from __future__ import print_function
 
 from past.utils import old_div
 import numpy as np
diff --git a/fluidsim/solvers/sw1l/test_operators.py b/fluidsim/solvers/sw1l/test_operators.py
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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/scripts/bugs/no_energy_conservation_ns2dstrat.py b/scripts/bugs/no_energy_conservation_ns2dstrat.py
--- 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
--- 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
--- 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/tox.ini b/tox.ini
--- a/tox.ini
+++ b/tox.ini
@@ -22,7 +22,6 @@
 deps =
     pip==18.0
     coverage
-    fluidpythran>=0.1.4
     mako
     numpy
     matplotlib
@@ -40,6 +39,7 @@
     pandas
     xarray
     setuptools_scm
+    hg+https://bitbucket.org/fluiddyn/fluidpythran#egg=fluidpythran
     hg+https://bitbucket.org/fluiddyn/fluiddyn#egg=fluiddyn
     hg+https://bitbucket.org/fluiddyn/fluidfft#egg=fluidfft
     pythran: git+https://github.com/serge-sans-paille/pythran#egg=pythran[deps]