diff --git a/bench/profiling/util_bench.py b/bench/profiling/util_bench.py
index dcfee750993158507a19b585c48795ecd9fe3f98_YmVuY2gvcHJvZmlsaW5nL3V0aWxfYmVuY2gucHk=..6360550698ec6db213709bb2f849713d9c584550_YmVuY2gvcHJvZmlsaW5nL3V0aWxfYmVuY2gucHk= 100644
--- a/bench/profiling/util_bench.py
+++ b/bench/profiling/util_bench.py
@@ -128,7 +128,7 @@
     nfile = (
         'result_bench_' + key_solver + '_' +
         sim.oper.produce_str_describing_grid() +
-        '_' + t_as_str + '_{}'.format(pid) + '.json')
+        '_' + t_as_str + f'_{pid}' + '.json')
 
     path = os.path.join(path_results, nfile)
 
diff --git a/bench/submit_benchmarks_legi.py b/bench/submit_benchmarks_legi.py
index dcfee750993158507a19b585c48795ecd9fe3f98_YmVuY2gvc3VibWl0X2JlbmNobWFya3NfbGVnaS5weQ==..6360550698ec6db213709bb2f849713d9c584550_YmVuY2gvc3VibWl0X2JlbmNobWFya3NfbGVnaS5weQ== 100644
--- a/bench/submit_benchmarks_legi.py
+++ b/bench/submit_benchmarks_legi.py
@@ -17,7 +17,7 @@
         ('fluidsim bench 1024 -d 2 -s '
          '{} -o /.fsnet/data/legi/calcul9/home/augier3pi/fluidsim_bench'
         ).format(solver),
-        name_run='fluidsim-bench_{}'.format(nb_mpi),
+        name_run=f'fluidsim-bench_{nb_mpi}',
         nb_nodes=nb_nodes,
         # nb_cores_per_node=nb_cores_per_node,
         nb_cores_per_node=cluster.nb_cores_per_node,
diff --git a/bench/util.py b/bench/util.py
index dcfee750993158507a19b585c48795ecd9fe3f98_YmVuY2gvdXRpbC5weQ==..6360550698ec6db213709bb2f849713d9c584550_YmVuY2gvdXRpbC5weQ== 100644
--- a/bench/util.py
+++ b/bench/util.py
@@ -28,7 +28,7 @@
         n1 = n0
 
     params._set_child('two_d', dict(
-        shape='{} {}'.format(n0, n1), time='00:20:00',
+        shape=f'{n0} {n1}', time='00:20:00',
         solver='ns2d',
         fft_seq=['fft2d.with_fftw1d',
                 'fft2d.with_fftw2d'],
@@ -41,7 +41,7 @@
         n2 = n0
 
     params._set_child('three_d', dict(
-        shape='{} {} {}'.format(n0, n1, n2), time='00:30:00',
+        shape=f'{n0} {n1} {n2}', time='00:30:00',
         solver='ns3d',
         fft_seq=['fft3d.with_fftw3d'],
         fft=['fft3d.mpi_with_fftw1d',
@@ -155,7 +155,7 @@
     else:
         cluster.submit_command(
             cmd,
-            name_run='{}_{}'.format(params_dim.solver, nb_mpi),
+            name_run=f'{params_dim.solver}_{nb_mpi}',
             nb_nodes=nb_nodes,
             nb_cores_per_node=nb_cores_per_node,
             walltime=params_dim.time,
diff --git a/doc/conf.py b/doc/conf.py
index dcfee750993158507a19b585c48795ecd9fe3f98_ZG9jL2NvbmYucHk=..6360550698ec6db213709bb2f849713d9c584550_ZG9jL2NvbmYucHk= 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -91,8 +91,8 @@
 master_doc = "index"
 
 # General information about the project.
-project = u"FluidSim"
-copyright = u"2015, Pierre Augier"
+project = "FluidSim"
+copyright = "2015, Pierre Augier"
 
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
@@ -239,8 +239,8 @@
     (
         "index",
         "fluidsim.tex",
-        u"fluidsim Documentation",
-        u"Pierre Augier",
+        "fluidsim Documentation",
+        "Pierre Augier",
         "manual",
     )
 ]
@@ -271,7 +271,7 @@
 # One entry per manual page. List of tuples
 # (source start file, name, description, authors, manual section).
 man_pages = [
-    ("index", "FluidSim", u"FluidSim Documentation", [u"Pierre Augier"], 1)
+    ("index", "FluidSim", "FluidSim Documentation", ["Pierre Augier"], 1)
 ]
 
 # If true, show URL addresses after external links.
@@ -287,8 +287,8 @@
     (
         "index",
         "FluidDyn",
-        u"FluidDyn Documentation",
-        u"Pierre Augier",
+        "FluidDyn Documentation",
+        "Pierre Augier",
         "FluidDyn",
         "One line description of project.",
         "Miscellaneous",
diff --git a/fluidsim/base/basilisk/operators.py b/fluidsim/base/basilisk/operators.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay9vcGVyYXRvcnMucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay9vcGVyYXRvcnMucHk= 100644
--- a/fluidsim/base/basilisk/operators.py
+++ b/fluidsim/base/basilisk/operators.py
@@ -14,7 +14,7 @@
 import basilisk.stream as basilisk
 
 
-class OperatorsBasilisk2D(object):
+class OperatorsBasilisk2D:
     """2D operators.
 
     """
diff --git a/fluidsim/base/basilisk/solver.py b/fluidsim/base/basilisk/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay9zb2x2ZXIucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay9zb2x2ZXIucHk= 100644
--- a/fluidsim/base/basilisk/solver.py
+++ b/fluidsim/base/basilisk/solver.py
@@ -24,7 +24,7 @@
 
     def _init_root(self):
 
-        super(InfoSolverBasilisk, self)._init_root()
+        super()._init_root()
 
         mod = "fluidsim.base.basilisk"
 
@@ -55,7 +55,7 @@
     def __init__(self, params):
         """Initialize parameters, state fields, and event loop."""
         bas = self.basilisk = basilisk
-        super(SimulBasilisk, self).__init__(params)
+        super().__init__(params)
 
         def init(i, t):
             bas.omega.f = bas.noise
diff --git a/fluidsim/base/basilisk/time_stepping.py b/fluidsim/base/basilisk/time_stepping.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay90aW1lX3N0ZXBwaW5nLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay90aW1lX3N0ZXBwaW5nLnB5 100644
--- a/fluidsim/base/basilisk/time_stepping.py
+++ b/fluidsim/base/basilisk/time_stepping.py
@@ -16,7 +16,7 @@
 import numpy as np
 
 
-class TimeSteppingBasilisk(object):
+class TimeSteppingBasilisk:
     """Time stepping class to handle Basilisk's event loop and FluidSim output.
 
     """
@@ -38,7 +38,7 @@
         self._has_to_stop = False
 
         def handler_signals(signal_number, stack):
-            print("signal {} received.".format(signal_number))
+            print(f"signal {signal_number} received.")
             self._has_to_stop = True
 
         signal(12, handler_signals)
diff --git a/fluidsim/base/forcing/anisotropic.py b/fluidsim/base/forcing/anisotropic.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9mb3JjaW5nL2FuaXNvdHJvcGljLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9mb3JjaW5nL2FuaXNvdHJvcGljLnB5 100644
--- a/fluidsim/base/forcing/anisotropic.py
+++ b/fluidsim/base/forcing/anisotropic.py
@@ -118,6 +118,6 @@
         title = (
             pforcing.type
             + "; "
-            + r"$nk_{{min}} = {} \delta k_x$; ".format(pforcing.nkmin_forcing)
-            + r"$nk_{{max}} = {} \delta k_z$; ".format(pforcing.nkmax_forcing)
+            + fr"$nk_{{min}} = {pforcing.nkmin_forcing} \delta k_x$; "
+            + fr"$nk_{{max}} = {pforcing.nkmax_forcing} \delta k_z$; "
             + r"$\theta = {:.0f}^\circ$; ".format(degrees(self.angle))
@@ -123,5 +123,5 @@
             + r"$\theta = {:.0f}^\circ$; ".format(degrees(self.angle))
-            + r"Forced modes = {}".format(self.nb_forced_modes)
+            + fr"Forced modes = {self.nb_forced_modes}"
         )
 
         ax.set_title(title)
diff --git a/fluidsim/base/forcing/base.py b/fluidsim/base/forcing/base.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9mb3JjaW5nL2Jhc2UucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9mb3JjaW5nL2Jhc2UucHk= 100644
--- a/fluidsim/base/forcing/base.py
+++ b/fluidsim/base/forcing/base.py
@@ -28,7 +28,7 @@
 )
 
 
-class ForcingBase(object):
+class ForcingBase:
     """Organize the forcing schemes (base class)"""
 
     @staticmethod
diff --git a/fluidsim/base/forcing/specific.py b/fluidsim/base/forcing/specific.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9mb3JjaW5nL3NwZWNpZmljLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9mb3JjaW5nL3NwZWNpZmljLnB5 100644
--- a/fluidsim/base/forcing/specific.py
+++ b/fluidsim/base/forcing/specific.py
@@ -53,7 +53,7 @@
 from fluidsim.base.setofvariables import SetOfVariables
 
 
-class SpecificForcing(object):
+class SpecificForcing:
     """Base class for specific forcing"""
 
     tag = "specific"
@@ -202,7 +202,7 @@
             pass
         else:
             if isinstance(angle, str):
-                if angle.endswith(u"°"):
+                if angle.endswith("°"):
                     angle = radians(float(angle[:-1]))
                 else:
                     raise ValueError(
@@ -468,7 +468,7 @@
     def _complete_params_with_default(cls, params):
         """This static method is used to complete the *params* container.
         """
-        super(NormalizedForcing, cls)._complete_params_with_default(params)
+        super()._complete_params_with_default(params)
         try:
             params.forcing.normalized
         except AttributeError:
@@ -696,7 +696,7 @@
     def _complete_params_with_default(cls, params):
         """This static method is used to complete the *params* container.
         """
-        super(RandomSimplePseudoSpectral, cls)._complete_params_with_default(
+        super()._complete_params_with_default(
             params
         )
 
@@ -743,9 +743,7 @@
     def _complete_params_with_default(cls, params):
         """This static method is used to complete the *params* container.
         """
-        super(
-            TimeCorrelatedRandomPseudoSpectral, cls
-        )._complete_params_with_default(params)
+        super()._complete_params_with_default(params)
 
         try:
             params.forcing.tcrandom
diff --git a/fluidsim/base/init_fields.py b/fluidsim/base/init_fields.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9pbml0X2ZpZWxkcy5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9pbml0X2ZpZWxkcy5weQ== 100644
--- a/fluidsim/base/init_fields.py
+++ b/fluidsim/base/init_fields.py
@@ -21,7 +21,7 @@
 from fluidsim.base.setofvariables import SetOfVariables
 
 
-class InitFieldsBase(object):
+class InitFieldsBase:
     """Initialization of the fields (base class)."""
 
     @staticmethod
@@ -138,7 +138,7 @@
         self._specific_init_fields()
 
 
-class SpecificInitFields(object):
+class SpecificInitFields:
     tag = "specific"
 
     @classmethod
@@ -155,7 +155,7 @@
 
     @classmethod
     def _complete_params_with_default(cls, params):
-        super(InitFieldsFromFile, cls)._complete_params_with_default(params)
+        super()._complete_params_with_default(params)
         params.init_fields._set_child(cls.tag, attribs={"path": ""})
         params.init_fields.from_file._set_doc(
             """
@@ -206,10 +206,10 @@
                 for r in axes:
                     # for example r can be: 'z', 'y', 'x'
                     r = r.decode("utf-8")
-                    nr = "n{}".format(r)
+                    nr = f"n{r}"
                     nr_file = group_oper.attrs[nr]
                     if params.oper[nr] != nr_file:
                         raise ValueError(
                             "this is not a correct state for this simulation\n"
                             "self.{0} != params_file.{0}".format(nr)
                         )
@@ -210,10 +210,10 @@
                     nr_file = group_oper.attrs[nr]
                     if params.oper[nr] != nr_file:
                         raise ValueError(
                             "this is not a correct state for this simulation\n"
                             "self.{0} != params_file.{0}".format(nr)
                         )
-                    Lr = "L{}".format(r)
+                    Lr = f"L{r}"
                     try:
                         Lr_file = group_oper.attrs[Lr]
                     except KeyError:
@@ -400,7 +400,7 @@
 
     @classmethod
     def _complete_params_with_default(cls, params):
-        super(InitFieldsConstant, cls)._complete_params_with_default(params)
+        super()._complete_params_with_default(params)
         params.init_fields._set_child(cls.tag, attribs={"value": 1.0})
         params.init_fields.constant._set_doc(
             """
@@ -424,7 +424,7 @@
     @classmethod
     def _complete_params_with_default(cls, params):
         """Complete the `params` container (class method)."""
-        super(InitFieldsNoise, cls)._complete_params_with_default(params)
+        super()._complete_params_with_default(params)
         params.init_fields._set_child(cls.tag, attribs={"max": 1.0})
         params.init_fields.noise._set_doc(
             """
diff --git a/fluidsim/base/output/__init__.py b/fluidsim/base/output/__init__.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvX19pbml0X18ucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvX19pbml0X18ucHk= 100644
--- a/fluidsim/base/output/__init__.py
+++ b/fluidsim/base/output/__init__.py
@@ -98,5 +98,5 @@
 
     if ndim == 1:
         geometry_type = "Origin_Dx"
-        dims_data = "{}".format(nx)
+        dims_data = f"{nx}"
         origins = "0"
@@ -102,4 +102,4 @@
         origins = "0"
-        deltaxs = "{}".format(deltax)
+        deltaxs = f"{deltax}"
     elif ndim == 2:
         geometry_type = "Origin_DxDy"
@@ -104,4 +104,4 @@
     elif ndim == 2:
         geometry_type = "Origin_DxDy"
-        dims_data = "{} {}".format(ny, nx)
+        dims_data = f"{ny} {nx}"
         origins = "0 0"
@@ -107,4 +107,4 @@
         origins = "0 0"
-        deltaxs = "{} {}".format(deltay, deltax)
+        deltaxs = f"{deltay} {deltax}"
     elif ndim == 3:
         geometry_type = "Origin_DxDyDz"
@@ -109,4 +109,4 @@
     elif ndim == 3:
         geometry_type = "Origin_DxDyDz"
-        dims_data = "{} {} {}".format(nz, ny, nx)
+        dims_data = f"{nz} {ny} {nx}"
         origins = "0 0 0"
@@ -112,5 +112,5 @@
         origins = "0 0 0"
-        deltaxs = "{} {} {}".format(deltaz, deltay, deltax)
+        deltaxs = f"{deltaz} {deltay} {deltax}"
 
     if ndim in (2, 3):
         vectors = []
@@ -124,7 +124,7 @@
         for key in keys:
             if key.endswith("x"):
                 vector = key[:-1]
-                vector_components = set([vector + compo for compo in components])
+                vector_components = {vector + compo for compo in components}
                 if vector_components.issubset(set(keys)):
                     vectors.append(vector)
 
diff --git a/fluidsim/base/output/base.py b/fluidsim/base/output/base.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvYmFzZS5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvYmFzZS5weQ== 100644
--- a/fluidsim/base/output/base.py
+++ b/fluidsim/base/output/base.py
@@ -36,7 +36,7 @@
 from fluidsim.util.util import load_params_simul
 
 
-class OutputBase(object):
+class OutputBase:
     """Handle the output."""
 
     @staticmethod
@@ -306,7 +306,7 @@
             if mpi.nb_proc == 1:
                 specifications += "sequential,\n"
             else:
-                specifications += "parallel ({} proc.)\n".format(mpi.nb_proc)
+                specifications += f"parallel ({mpi.nb_proc} proc.)\n"
             self.print_stdout(
                 "\nsolver "
                 + self.name_solver
@@ -514,7 +514,7 @@
         mem = arr.nbytes * 1.0e-6
         if mpi.nb_proc > 1:
             mem = mpi.comm.allreduce(mem, op=mpi.MPI.SUM)
-        self.print_stdout(string.ljust(30) + ": {0} Mo".format(mem))
+        self.print_stdout(string.ljust(30) + f": {mem} Mo")
 
 
 class OutputBasePseudoSpectral(OutputBase):
@@ -522,7 +522,7 @@
 
         oper = self.oper
         self.sum_wavenumbers = oper.sum_wavenumbers
-        super(OutputBasePseudoSpectral, self).init_with_oper_and_state()
+        super().init_with_oper_and_state()
 
     def compute_energy_fft(self):
         """Compute energy(k)"""
@@ -540,7 +540,7 @@
         return self.sum_wavenumbers(energy_fft)
 
 
-class SpecificOutput(object):
+class SpecificOutput:
     """Small class for features useful for specific outputs"""
 
     def __init__(
diff --git a/fluidsim/base/output/increments.py b/fluidsim/base/output/increments.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvaW5jcmVtZW50cy5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvaW5jcmVtZW50cy5weQ== 100644
--- a/fluidsim/base/output/increments.py
+++ b/fluidsim/base/output/increments.py
@@ -81,7 +81,7 @@
         self._bins = np.arange(0.5, self.nbins, dtype=float) / self.nbins
         self.keys_vars_to_compute = list(output.sim.state.state_phys.keys)
 
-        super(Increments, self).__init__(
+        super().__init__(
             output,
             period_save=params.output.periods_save.increments,
             has_to_plot_saved=params.output.increments.HAS_TO_PLOT_SAVED,
@@ -97,7 +97,7 @@
             axe.set_title(
                 r"pdf $\delta u_x (x)$, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
+                + f", nh = {self.nx:5d}"
             )
 
     def _online_plot_saving(self, dict_results, key="rot"):
@@ -192,7 +192,7 @@
         tmin_plot = times[imin_plot]
         tmax_plot = times[imax_plot]
 
-        to_print = "plot(tmin={0}, tmax={1}, delta_t={2:.2f})".format(
+        to_print = "plot(tmin={}, tmax={}, delta_t={:.2f})".format(
             tmin, tmax, delta_t
         )
         print(to_print)
@@ -225,9 +225,9 @@
         fig, ax1 = self.output.figure_axe(size_axe=size_axe)
         ax1.set_xlabel("$r_x$")
         ax1.set_ylabel(
-            r"$\langle \delta u^{" + "{0}".format(order) + "} \\rangle$"
+            r"$\langle \delta u^{" + f"{order}" + "} \\rangle$"
         )
 
         ax1.set_title(
             "struct. functions, solver "
             + self.output.name_solver
@@ -229,9 +229,9 @@
         )
 
         ax1.set_title(
             "struct. functions, solver "
             + self.output.name_solver
-            + ", nh = {0:5d}".format(self.nx)
+            + f", nh = {self.nx:5d}"
         )
         ax1.set_xscale("log")
         ax1.set_yscale(yscale)
@@ -463,10 +463,10 @@
             tmin=tmin, tmax=tmax, key_var=key_var, irx_to_plot=irx_to_plot
         )
 
-        to_print = "plot_pdf(tmin={0}, tmax={1})".format(tmin, tmax)
+        to_print = f"plot_pdf(tmin={tmin}, tmax={tmax})"
         print(to_print)
 
         fig, ax1 = self.output.figure_axe()
         ax1.set_title(
             "pdf increments, solver "
             + self.output.name_solver
@@ -467,10 +467,10 @@
         print(to_print)
 
         fig, ax1 = self.output.figure_axe()
         ax1.set_title(
             "pdf increments, solver "
             + self.output.name_solver
-            + ", nh = {0:5d}".format(self.nx)
+            + f", nh = {self.nx:5d}"
         )
         # +', c2 = {0:.4g}, f = {1:.4g}'.format(self.c2, self.f))
 
@@ -484,7 +484,7 @@
 
         for irxp, irx in enumerate(irx_to_plot):
 
-            print("color = {0}, rx = {1}".format(colors[irxp], self.rxs[irx]))
+            print("color = {}, rx = {}".format(colors[irxp], self.rxs[irx]))
 
             val_inc = values_inc_timemean[irxp]
 
diff --git a/fluidsim/base/output/movies.py b/fluidsim/base/output/movies.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvbW92aWVzLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvbW92aWVzLnB5 100644
--- a/fluidsim/base/output/movies.py
+++ b/fluidsim/base/output/movies.py
@@ -28,7 +28,7 @@
 from fluiddyn.util import mpi, is_run_from_jupyter
 
 
-class MoviesBase(object):
+class MoviesBase:
     """Base class defining most generic functions for movies."""
 
     def __init__(self, output):
@@ -374,7 +374,7 @@
     ):
         """Initializes animated figure."""
 
-        super(MoviesBase1D, self).init_animation(
+        super().init_animation(
             key_field, numfig, dt_equations, tmin, tmax, fig_kw, **kwargs
         )
 
@@ -400,5 +400,5 @@
         self._ani_line.set_data(x, y)
         self.ax.set_title(
             self.key_field
-            + ", $t = {0:.3f}$, ".format(time)
+            + f", $t = {time:.3f}$, "
             + self.output.name_solver
@@ -404,5 +404,5 @@
             + self.output.name_solver
-            + ", $n_x = {0:d}$".format(self.output.sim.oper.nx_seq)
+            + f", $n_x = {self.output.sim.oper.nx_seq:d}$"
         )
 
         return self._ani_line
diff --git a/fluidsim/base/output/phys_fields.py b/fluidsim/base/output/phys_fields.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMucHk= 100644
--- a/fluidsim/base/output/phys_fields.py
+++ b/fluidsim/base/output/phys_fields.py
@@ -91,7 +91,7 @@
             self.animate = self.movies.animate
             self.interact = self.movies.interact
 
-        super(PhysFieldsBase, self).__init__(
+        super().__init__(
             output,
             period_save=params.output.periods_save.phys_fields,
             period_plot=params.output.periods_plot.phys_fields,
@@ -155,7 +155,7 @@
                 time, self.sim.time_stepping.it, ext
             )
         else:
-            name_save = "state_phys_t{:07.3f}.{}".format(time, ext)
+            name_save = f"state_phys_t{time:07.3f}.{ext}"
 
         path_file = os.path.join(path_run, name_save)
         if os.path.exists(path_file):
diff --git a/fluidsim/base/output/phys_fields1d.py b/fluidsim/base/output/phys_fields1d.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMxZC5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMxZC5weQ== 100644
--- a/fluidsim/base/output/phys_fields1d.py
+++ b/fluidsim/base/output/phys_fields1d.py
@@ -26,7 +26,7 @@
 class MoviesBasePhysFields1D(MoviesBase1D):
     def __init__(self, output, phys_fields):
         self.phys_fields = phys_fields
-        super(MoviesBasePhysFields1D, self).__init__(output)
+        super().__init__(output)
         self._equation = None
 
 
@@ -75,5 +75,5 @@
     def _set_title(self, ax, key, time):
         title = (
             key
-            + ", $t = {0:.3f}$, ".format(time)
+            + f", $t = {time:.3f}$, "
             + self.output.name_solver
@@ -79,4 +79,4 @@
             + self.output.name_solver
-            + ", $n_x = {0:d}$".format(self.params.oper.nx)
+            + f", $n_x = {self.params.oper.nx:d}$"
         )
         ax.set_title(title)
diff --git a/fluidsim/base/output/phys_fields2d.py b/fluidsim/base/output/phys_fields2d.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMyZC5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMyZC5weQ== 100644
--- a/fluidsim/base/output/phys_fields2d.py
+++ b/fluidsim/base/output/phys_fields2d.py
@@ -32,7 +32,7 @@
 
     def __init__(self, output, phys_fields):
         self.phys_fields = phys_fields
-        super(MoviesBasePhysFields2D, self).__init__(output)
+        super().__init__(output)
         self._equation = None
 
     def init_animation(
@@ -45,8 +45,8 @@
 
         if dt_equations is None:
             dt_equations = np.median(np.diff(self.time_files))
-            print("dt_equations = {:.4f}".format(dt_equations))
+            print(f"dt_equations = {dt_equations:.4f}")
 
         if tmax is None:
             tmax = self.time_files.max()
 
@@ -49,8 +49,8 @@
 
         if tmax is None:
             tmax = self.time_files.max()
 
-        super(MoviesBasePhysFields2D, self).init_animation(
+        super().init_animation(
             key_field, numfig, dt_equations, tmin, tmax, fig_kw, **kwargs
         )
 
@@ -154,7 +154,7 @@
     def _get_axis_data(self, shape=None):
         """Get 1D arrays for setting the axes."""
 
-        x, y = super(MoviesBasePhysFields2D, self)._get_axis_data()
+        x, y = super()._get_axis_data()
         if shape is not None and (y.shape[0], x.shape[0]) != shape:
             path_file = os.path.join(self.output.path_run, "params_simul.xml")
             params = Parameters(path_file=path_file)
@@ -267,5 +267,5 @@
     def _set_title(self, ax, key, time, vmax=None):
         title = (
             key
-            + ", $t = {0:.3f}$, ".format(time)
+            + f", $t = {time:.3f}$, "
             + self.output.name_solver
@@ -271,4 +271,4 @@
             + self.output.name_solver
-            + ", $n_x = {0:d}$".format(self.params.oper.nx)
+            + f", $n_x = {self.params.oper.nx:d}$"
         )
         if vmax is not None:
@@ -273,6 +273,6 @@
         )
         if vmax is not None:
-            title += r", $|\vec{v}|_{max} = $" + "{0:.3f}".format(vmax)
+            title += r", $|\vec{v}|_{max} = $" + f"{vmax:.3f}"
         ax.set_title(title)
 
     def _init_online_plot(self):
diff --git a/fluidsim/base/output/phys_fields3d.py b/fluidsim/base/output/phys_fields3d.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMzZC5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcGh5c19maWVsZHMzZC5weQ== 100644
--- a/fluidsim/base/output/phys_fields3d.py
+++ b/fluidsim/base/output/phys_fields3d.py
@@ -27,7 +27,7 @@
 
 class PhysFieldsBase3D(PhysFieldsBase2D):
     def __init__(self, output):
-        super(PhysFieldsBase3D, self).__init__(output)
+        super().__init__(output)
         self.set_equation_crosssection("iz=0")
 
     def _init_movies(self):
diff --git a/fluidsim/base/output/print_stdout.py b/fluidsim/base/output/print_stdout.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcHJpbnRfc3Rkb3V0LnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcHJpbnRfc3Rkb3V0LnB5 100644
--- a/fluidsim/base/output/print_stdout.py
+++ b/fluidsim/base/output/print_stdout.py
@@ -6,7 +6,7 @@
 from fluiddyn.util import mpi, print_memory_usage
 
 
-class PrintStdOutBase(object):
+class PrintStdOutBase:
     """A :class:`PrintStdOutBase` object is used to print in both the
     stdout and the stdout.txt file, and also to print simple info on
     the current state of the simulation."""
@@ -77,7 +77,7 @@
         print_memory_usage("MEMORY_USAGE:")
 
     def _make_str_info(self):
-        return "it = {0:6d} ; t = {1:10.6g} ; deltat  = {2:10.5g}\n".format(
+        return "it = {:6d} ; t = {:10.6g} ; deltat  = {:10.5g}\n".format(
             self.sim.time_stepping.it,
             self.sim.time_stepping.t,
             self.sim.time_stepping.deltat,
diff --git a/fluidsim/base/output/prob_dens_func.py b/fluidsim/base/output/prob_dens_func.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcHJvYl9kZW5zX2Z1bmMucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvcHJvYl9kZW5zX2Z1bmMucHk= 100644
--- a/fluidsim/base/output/prob_dens_func.py
+++ b/fluidsim/base/output/prob_dens_func.py
@@ -30,7 +30,7 @@
         self.f = params.f
         self.nx = params.oper.nx
 
-        super(ProbaDensityFunc, self).__init__(
+        super().__init__(
             output,
             period_save=params.output.periods_save.pdf,
             has_to_plot_saved=params.output.pdf.HAS_TO_PLOT_SAVED,
@@ -45,8 +45,8 @@
             title = (
                 r"pdf $\eta$, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
-                + ", c = {0:.4g}, f = {1:.4g}".format(np.sqrt(self.c2), self.f)
+                + f", nh = {self.nx:5d}"
+                + ", c = {:.4g}, f = {:.4g}".format(np.sqrt(self.c2), self.f)
             )
             axe.set_title(title)
 
@@ -115,8 +115,8 @@
         ax1.set_title(
             "PDF, solver "
             + self.output.name_solver
-            + ", nh = {0:5d}".format(self.nx)
-            + ", c = {0:.4g}, f = {1:.4g}".format(np.sqrt(self.c2), self.f)
+            + f", nh = {self.nx:5d}"
+            + ", c = {:.4g}, f = {:.4g}".format(np.sqrt(self.c2), self.f)
         )
         ax1.set_xscale("linear")
         ax1.set_yscale("linear")
@@ -156,7 +156,7 @@
             tmin_plot = times[imin_plot]
             tmax_plot = times[imax_plot]
 
-            to_print = "plot(tmin={0}, tmax={1}, delta_t={2:.2f})".format(
+            to_print = "plot(tmin={}, tmax={}, delta_t={:.2f})".format(
                 tmin, tmax, delta_t
             )
             print(to_print)
diff --git a/fluidsim/base/output/spatial_means.py b/fluidsim/base/output/spatial_means.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvc3BhdGlhbF9tZWFucy5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvc3BhdGlhbF9tZWFucy5weQ== 100644
--- a/fluidsim/base/output/spatial_means.py
+++ b/fluidsim/base/output/spatial_means.py
@@ -89,7 +89,7 @@
             title = (
                 "mean quantities, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
+                + f", nh = {self.nx:5d}"
             )
             axe.set_title(title)
 
diff --git a/fluidsim/base/output/spect_energy_budget.py b/fluidsim/base/output/spect_energy_budget.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvc3BlY3RfZW5lcmd5X2J1ZGdldC5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvc3BlY3RfZW5lcmd5X2J1ZGdldC5weQ== 100644
--- a/fluidsim/base/output/spect_energy_budget.py
+++ b/fluidsim/base/output/spect_energy_budget.py
@@ -42,7 +42,7 @@
         self.sum_wavenumbers = output.sim.oper.sum_wavenumbers
 
         HAS_TO_PLOT_SAVED = params.output.spect_energy_budg.HAS_TO_PLOT_SAVED
-        super(SpectralEnergyBudgetBase, self).__init__(
+        super().__init__(
             output,
             period_save=params.output.periods_save.spect_energy_budg,
             has_to_plot_saved=HAS_TO_PLOT_SAVED,
@@ -76,7 +76,7 @@
             axe_a.set_title(
                 "energy flux, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
+                + f", nh = {self.nx:5d}"
             )
             axe_a.set_xscale("log")
 
diff --git a/fluidsim/base/output/spectra.py b/fluidsim/base/output/spectra.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvc3BlY3RyYS5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvc3BlY3RyYS5weQ== 100644
--- a/fluidsim/base/output/spectra.py
+++ b/fluidsim/base/output/spectra.py
@@ -12,7 +12,7 @@
 class MoviesSpectra(MoviesBase1D):
     def __init__(self, output, spectra):
         self.spectra = spectra
-        super(MoviesSpectra, self).__init__(output)
+        super().__init__(output)
 
     def init_animation(self, *args, **kwargs):
         if "xmax" not in kwargs:
@@ -23,7 +23,7 @@
         with h5py.File(self.spectra.path_file2D) as f:
             self.times = f["times"][...]
 
-        super(MoviesSpectra, self).init_animation(*args, **kwargs)
+        super().init_animation(*args, **kwargs)
 
     def get_field_to_plot(self, time, key=None):
         if key is None:
@@ -96,7 +96,7 @@
         self.spectrum2D_from_fft = output.sim.oper.spectrum2D_from_fft
         self.spectra1D_from_fft = output.sim.oper.spectra1D_from_fft
 
-        super(Spectra, self).__init__(
+        super().__init__(
             output,
             period_save=params.output.periods_save.spectra,
             has_to_plot_saved=params.output.spectra.HAS_TO_PLOT_SAVED,
@@ -168,7 +168,7 @@
             axe.set_title(
                 "spectra, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
+                + f", nh = {self.nx:5d}"
             )
 
     def _online_plot_saving(self):
diff --git a/fluidsim/base/output/spectra3d.py b/fluidsim/base/output/spectra3d.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvc3BlY3RyYTNkLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvc3BlY3RyYTNkLnB5 100644
--- a/fluidsim/base/output/spectra3d.py
+++ b/fluidsim/base/output/spectra3d.py
@@ -12,7 +12,7 @@
 class MoviesSpectra(MoviesBase1D):
     def __init__(self, output, spectra):
         self.spectra = spectra
-        super(MoviesSpectra, self).__init__(output)
+        super().__init__(output)
 
     def init_animation(self, *args, **kwargs):
         if "xmax" not in kwargs:
@@ -23,7 +23,7 @@
         with h5py.File(self.spectra.path_file3d) as f:
             self.times = f["times"][...]
 
-        super(MoviesSpectra, self).init_animation(*args, **kwargs)
+        super().init_animation(*args, **kwargs)
 
     def get_field_to_plot(self, time, key=None):
         if key is None:
@@ -93,7 +93,7 @@
         params = output.sim.params
         self.nx = int(params.oper.nx)
 
-        super(Spectra, self).__init__(
+        super().__init__(
             output,
             period_save=params.output.periods_save.spectra,
             has_to_plot_saved=params.output.spectra.HAS_TO_PLOT_SAVED,
@@ -166,7 +166,7 @@
             axe.set_title(
                 "spectra, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
+                + f", nh = {self.nx:5d}"
             )
 
     def _online_plot_saving(self, dict_spectra1d, dict_spectra3d):
diff --git a/fluidsim/base/output/spectra_multidim.py b/fluidsim/base/output/spectra_multidim.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvc3BlY3RyYV9tdWx0aWRpbS5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvc3BlY3RyYV9tdWx0aWRpbS5weQ== 100644
--- a/fluidsim/base/output/spectra_multidim.py
+++ b/fluidsim/base/output/spectra_multidim.py
@@ -29,7 +29,7 @@
         if not params.output.HAS_TO_SAVE:
             params.output.periods_save.spectra_multidim = False
 
-        super(SpectraMultiDim, self).__init__(
+        super().__init__(
             output,
             period_save=params.output.periods_save.spectra_multidim,
             has_to_plot_saved=params.output.spectra_multidim.HAS_TO_PLOT_SAVED,
@@ -92,7 +92,7 @@
             axe.set_title(
                 "Multidimensional spectra, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
+                + f", nh = {self.nx:5d}"
             )
 
     def _online_plot_saving(self, arg):
diff --git a/fluidsim/base/output/time_signals_fft.py b/fluidsim/base/output/time_signals_fft.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvdGltZV9zaWduYWxzX2ZmdC5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvdGltZV9zaWduYWxzX2ZmdC5weQ== 100644
--- a/fluidsim/base/output/time_signals_fft.py
+++ b/fluidsim/base/output/time_signals_fft.py
@@ -48,7 +48,7 @@
         if params.output.periods_save.time_signals_fft:
             self._init_save(sim)
 
-        super(TimeSignalsK, self).__init__(
+        super().__init__(
             output,
             period_save=params.output.periods_save.time_signals_fft,
             period_plot=params.output.periods_plot.time_signals_fft,
@@ -326,7 +326,7 @@
             ax1.set_xlabel("$t/T$")
             ax1.set_ylabel("signals (s$^{-1}$)")
             title = (
-                "signals eigenmodes, ikh = {0:.2f}, solver ".format(
+                "signals eigenmodes, ikh = {:.2f}, solver ".format(
                     (kh_shell[ish] / self.sim.oper.deltak)
                 )
                 + self.output.name_solver
@@ -330,8 +330,8 @@
                     (kh_shell[ish] / self.sim.oper.deltak)
                 )
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
-                + ", c2 = {0:.4g}, f = {1:.4g}".format(self.c2, self.f)
+                + f", nh = {self.nx:5d}"
+                + f", c2 = {self.c2:.4g}, f = {self.f:.4g}"
             )
             ax1.set_title(title)
 
@@ -457,8 +457,8 @@
         title = (
             "time spectra, solver "
             + self.output.name_solver
-            + ", nh = {0:5d}".format(self.nx)
-            + ", c = {0:.4g}, f = {1:.4g}".format(np.sqrt(self.c2), self.f)
+            + f", nh = {self.nx:5d}"
+            + ", c = {:.4g}, f = {:.4g}".format(np.sqrt(self.c2), self.f)
         )
         ax1.set_title(title)
 
diff --git a/fluidsim/base/preprocess/base.py b/fluidsim/base/preprocess/base.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9wcmVwcm9jZXNzL2Jhc2UucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9wcmVwcm9jZXNzL2Jhc2UucHk= 100644
--- a/fluidsim/base/preprocess/base.py
+++ b/fluidsim/base/preprocess/base.py
@@ -11,7 +11,7 @@
 from builtins import object
 
 
-class PreprocessBase(object):
+class PreprocessBase:
     _tag = "preprocess"
 
     @staticmethod
diff --git a/fluidsim/base/preprocess/pseudo_spect.py b/fluidsim/base/preprocess/pseudo_spect.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9wcmVwcm9jZXNzL3BzZXVkb19zcGVjdC5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9wcmVwcm9jZXNzL3BzZXVkb19zcGVjdC5weQ== 100644
--- a/fluidsim/base/preprocess/pseudo_spect.py
+++ b/fluidsim/base/preprocess/pseudo_spect.py
@@ -20,7 +20,7 @@
     def __call__(self):
         """Preprocesses if enabled."""
 
-        super(PreprocessPseudoSpectral, self).__call__()
+        super().__call__()
         if self.params.enable:
             if self.sim.params.forcing.enable:
                 if "forcing" in self.params.init_field_scale:
@@ -98,7 +98,7 @@
             raise ValueError("Unknown initial fields scaling: ", scale)
 
     def set_viscosity(self):
-        """Based on
+        r"""Based on
 
         - the initial total enstrophy, \Omega_0, or
 
@@ -356,7 +356,7 @@
     if verbose:
         length_scale = np.mean(kolmo_len)
         k_diss = 1.0 / length_scale
-        print("\nCALCULATED (eps={})".format(epsilon))
+        print(f"\nCALCULATED (eps={epsilon})")
         print(
             f"Dissipation wavenumber, k_d = {k_diss}; k_d / k_f = {k_diss / k_f}"
         )
diff --git a/fluidsim/base/solvers/finite_diff.py b/fluidsim/base/solvers/finite_diff.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9zb2x2ZXJzL2Zpbml0ZV9kaWZmLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9zb2x2ZXJzL2Zpbml0ZV9kaWZmLnB5 100644
--- a/fluidsim/base/solvers/finite_diff.py
+++ b/fluidsim/base/solvers/finite_diff.py
@@ -6,7 +6,7 @@
 
     def _init_root(self):
 
-        super(InfoSolverFiniteDiff, self)._init_root()
+        super()._init_root()
 
         # self.classes.State.module_name = 'fluidsim.base.state'
         # self.classes.State.class_name = 'StateBase'
diff --git a/fluidsim/base/solvers/info_base.py b/fluidsim/base/solvers/info_base.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9zb2x2ZXJzL2luZm9fYmFzZS5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9zb2x2ZXJzL2luZm9fYmFzZS5weQ== 100644
--- a/fluidsim/base/solvers/info_base.py
+++ b/fluidsim/base/solvers/info_base.py
@@ -26,7 +26,7 @@
         if len(kargs) == 0 or ("path_file" in kargs and "tag" not in kargs):
             kargs["tag"] = "solver"
 
-        super(InfoSolverBase, self).__init__(**kargs)
+        super().__init__(**kargs)
 
         if (
             "tag" in kargs
diff --git a/fluidsim/base/solvers/pseudo_spect.py b/fluidsim/base/solvers/pseudo_spect.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9zb2x2ZXJzL3BzZXVkb19zcGVjdC5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9zb2x2ZXJzL3BzZXVkb19zcGVjdC5weQ== 100644
--- a/fluidsim/base/solvers/pseudo_spect.py
+++ b/fluidsim/base/solvers/pseudo_spect.py
@@ -45,7 +45,7 @@
         - :class:`fluidsim.operators.operators.OperatorsPseudoSpectral2D`
 
         """
-        super(InfoSolverPseudoSpectral, self)._init_root()
+        super()._init_root()
 
         self.module_name = "fluidsim.base.solvers.pseudo_spect"
         self.class_name = "SimulBasePseudoSpectral"
@@ -92,7 +92,7 @@
 
         """
 
-        super(InfoSolverPseudoSpectral3D, self)._init_root()
+        super()._init_root()
 
         self.classes.Operators.module_name = "fluidsim.operators.operators3d"
         self.classes.Operators.class_name = "OperatorsPseudoSpectral3D"
diff --git a/fluidsim/base/sphericalharmo/phys_fields.py b/fluidsim/base/sphericalharmo/phys_fields.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9zcGhlcmljYWxoYXJtby9waHlzX2ZpZWxkcy5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9zcGhlcmljYWxoYXJtby9waHlzX2ZpZWxkcy5weQ== 100644
--- a/fluidsim/base/sphericalharmo/phys_fields.py
+++ b/fluidsim/base/sphericalharmo/phys_fields.py
@@ -34,5 +34,5 @@
     def _set_title(self, ax, key, time, vmax=None):
         title = (
             key
-            + ", $t = {0:.3f}$, ".format(time)
+            + f", $t = {time:.3f}$, "
             + self.output.name_solver
@@ -38,4 +38,4 @@
             + self.output.name_solver
-            + r", $l_\max = {0:d}$".format(self.params.oper.lmax)
+            + fr", $l_\max = {self.params.oper.lmax:d}$"
         )
         if vmax is not None:
@@ -40,6 +40,6 @@
         )
         if vmax is not None:
-            title += r", $|\vec{v}|_{max} = $" + "{0:.3f}".format(vmax)
+            title += r", $|\vec{v}|_{max} = $" + f"{vmax:.3f}"
         ax.set_title(title)
 
     def _compute_skip_quiver(self):
diff --git a/fluidsim/base/sphericalharmo/solver.py b/fluidsim/base/sphericalharmo/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9zcGhlcmljYWxoYXJtby9zb2x2ZXIucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9zcGhlcmljYWxoYXJtby9zb2x2ZXIucHk= 100644
--- a/fluidsim/base/sphericalharmo/solver.py
+++ b/fluidsim/base/sphericalharmo/solver.py
@@ -31,7 +31,7 @@
 
         """
 
-        super(InfoSolverSphericalHarmo, self)._init_root()
+        super()._init_root()
 
         here = "fluidsim.base.sphericalharmo"
 
diff --git a/fluidsim/base/sphericalharmo/state.py b/fluidsim/base/sphericalharmo/state.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9zcGhlcmljYWxoYXJtby9zdGF0ZS5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9zcGhlcmljYWxoYXJtby9zdGF0ZS5weQ== 100644
--- a/fluidsim/base/sphericalharmo/state.py
+++ b/fluidsim/base/sphericalharmo/state.py
@@ -111,4 +111,4 @@
             if "rot_sh" in kwargs:
                 self.init_from_rotsh(kwargs["rot_sh"])
         else:
-            super(StateSphericalHarmo, self).init_statespect_from(**kwargs)
+            super().init_statespect_from(**kwargs)
diff --git a/fluidsim/base/state.py b/fluidsim/base/state.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS9zdGF0ZS5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS9zdGF0ZS5weQ== 100644
--- a/fluidsim/base/state.py
+++ b/fluidsim/base/state.py
@@ -20,7 +20,7 @@
 from fluidsim.base.setofvariables import SetOfVariables
 
 
-class StateBase(object):
+class StateBase:
     """Contains the state variables and handles the access to fields.
 
     Parameters
@@ -202,7 +202,7 @@
         for key, value in list(kwargs.items()):
             if key not in self.keys_state_phys:
                 raise ValueError(
-                    'Do not know how to initialize with key "{}".'.format(key)
+                    f'Do not know how to initialize with key "{key}".'
                 )
 
             self.state_phys.set_var(key, value)
@@ -240,7 +240,7 @@
 
     def __init__(self, sim, oper=None):
 
-        super(StatePseudoSpectral, self).__init__(sim, oper)
+        super().__init__(sim, oper)
 
         self.keys_state_spect = sim.info.solver.classes.State.keys_state_spect
         self.state_spect = SetOfVariables(
@@ -380,7 +380,7 @@
         for key, value in list(kwargs.items()):
             if key not in self.keys_state_spect:
                 raise ValueError(
-                    'Do not know how to initialize with key "{}".'.format(key)
+                    f'Do not know how to initialize with key "{key}".'
                 )
 
             self.state_spect.set_var(key, value)
diff --git a/fluidsim/base/time_stepping/base.py b/fluidsim/base/time_stepping/base.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS90aW1lX3N0ZXBwaW5nL2Jhc2UucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS90aW1lX3N0ZXBwaW5nL2Jhc2UucHk= 100644
--- a/fluidsim/base/time_stepping/base.py
+++ b/fluidsim/base/time_stepping/base.py
@@ -16,7 +16,7 @@
 from fluiddyn.util import mpi
 
 
-class TimeSteppingBase(object):
+class TimeSteppingBase:
     """Universal time stepping class used for all solvers.
 
 
@@ -90,7 +90,7 @@
         self._has_to_stop = False
 
         def handler_signals(signal_number, stack):
-            print("signal {} received.".format(signal_number))
+            print(f"signal {signal_number} received.")
             self._has_to_stop = True
 
         try:
@@ -182,7 +182,7 @@
             self.sim.output.phys_fields.save()
         if self.params.time_stepping.USE_T_END:
             print_stdout(
-                "    compute until t = {0:10.6g}".format(
+                "    compute until t = {:10.6g}".format(
                     self.params.time_stepping.t_end
                 )
             )
@@ -192,7 +192,7 @@
                 self.one_time_step()
         else:
             print_stdout(
-                "    compute until it = {0:8d}".format(
+                "    compute until it = {:8d}".format(
                     self.params.time_stepping.it_end
                 )
             )
diff --git a/fluidsim/base/time_stepping/finite_diff.py b/fluidsim/base/time_stepping/finite_diff.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS90aW1lX3N0ZXBwaW5nL2Zpbml0ZV9kaWZmLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS90aW1lX3N0ZXBwaW5nL2Zpbml0ZV9kaWZmLnB5 100644
--- a/fluidsim/base/time_stepping/finite_diff.py
+++ b/fluidsim/base/time_stepping/finite_diff.py
@@ -32,7 +32,7 @@
         params.time_stepping.type_time_scheme = "RK2"
 
     def __init__(self, sim):
-        super(TimeSteppingFiniteDiffCrankNicolson, self).__init__(sim)
+        super().__init__(sim)
 
         self._init_compute_time_step()
         self._init_time_scheme()
@@ -44,7 +44,7 @@
         self._time_step_RK()
         if np.isnan(np.min(self.sim.state.state_phys)):
             raise ValueError(
-                "nan at it = {0}, t = {1:.4f}".format(self.it, self.t)
+                f"nan at it = {self.it}, t = {self.t:.4f}"
             )
 
     def _time_step_RK2(self):
diff --git a/fluidsim/base/time_stepping/pseudo_spect.py b/fluidsim/base/time_stepping/pseudo_spect.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS90aW1lX3N0ZXBwaW5nL3BzZXVkb19zcGVjdC5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS90aW1lX3N0ZXBwaW5nL3BzZXVkb19zcGVjdC5weQ== 100644
--- a/fluidsim/base/time_stepping/pseudo_spect.py
+++ b/fluidsim/base/time_stepping/pseudo_spect.py
@@ -47,7 +47,7 @@
     use_cython = True
 
 
-class ExactLinearCoefs(object):
+class ExactLinearCoefs:
     """Handle the computation of the exact coefficient for the RK4."""
 
     def __init__(self, time_stepping):
@@ -98,7 +98,7 @@
         params.time_stepping.USE_CFL = True
 
     def __init__(self, sim):
-        super(TimeSteppingPseudoSpectral, self).__init__(sim)
+        super().__init__(sim)
 
         self._init_freq_lin()
         self._init_compute_time_step()
@@ -188,7 +188,7 @@
         # np.isnan(np.sum seems to be really fast
         if np.isnan(np.sum(self.sim.state.state_spect[0])):
             raise ValueError(
-                "nan at it = {0}, t = {1:.4f}".format(self.it, self.t)
+                f"nan at it = {self.it}, t = {self.t:.4f}"
             )
 
     def _time_step_RK2(self):
diff --git a/fluidsim/base/time_stepping/simple.py b/fluidsim/base/time_stepping/simple.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vYmFzZS90aW1lX3N0ZXBwaW5nL3NpbXBsZS5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vYmFzZS90aW1lX3N0ZXBwaW5nL3NpbXBsZS5weQ== 100644
--- a/fluidsim/base/time_stepping/simple.py
+++ b/fluidsim/base/time_stepping/simple.py
@@ -16,7 +16,7 @@
 from .base import TimeSteppingBase
 
 
-class ExactLinearCoefs(object):
+class ExactLinearCoefs:
     """Handle the computation of the exact coefficient for the RK4."""
 
     def __init__(self, time_stepping):
@@ -59,7 +59,7 @@
     """
 
     def __init__(self, sim):
-        super(TimeSteppingSimple, self).__init__(sim)
+        super().__init__(sim)
 
         self._init_freq_lin()
         self._init_compute_time_step()
@@ -80,7 +80,7 @@
         self._time_step_RK()
         if np.isnan(np.sum(self.sim.state.state_phys[0])):
             raise ValueError(
-                "nan at it = {0}, t = {1:.4f}".format(self.it, self.t)
+                f"nan at it = {self.it}, t = {self.t:.4f}"
             )
 
     def _time_step_RK2(self):
diff --git a/fluidsim/magic.py b/fluidsim/magic.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vbWFnaWMucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vbWFnaWMucHk= 100644
--- a/fluidsim/magic.py
+++ b/fluidsim/magic.py
@@ -77,7 +77,7 @@
     """
 
     def __init__(self, shell, package_solvers):
-        super(FluidsimMagics, self).__init__(shell)
+        super().__init__(shell)
         self.package = package_solvers
 
     def is_defined(self, varname):
diff --git a/fluidsim/operators/base.py b/fluidsim/operators/base.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vb3BlcmF0b3JzL2Jhc2UucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vb3BlcmF0b3JzL2Jhc2UucHk= 100644
--- a/fluidsim/operators/base.py
+++ b/fluidsim/operators/base.py
@@ -14,7 +14,7 @@
 import numpy as np
 
 
-class OperatorsBase1D(object):
+class OperatorsBase1D:
     @staticmethod
     def _complete_params_with_default(params):
         """This static method is used to complete the *params* container.
@@ -49,7 +49,7 @@
         if (self.Lx / np.pi).is_integer():
             str_Lx = repr(int(self.Lx / np.pi)) + "pi"
         else:
-            str_Lx = "{:.3f}".format(self.Lx).rstrip("0")
+            str_Lx = f"{self.Lx:.3f}".rstrip("0")
 
         return str_Lx
 
@@ -62,8 +62,8 @@
         """Produce a string describing the operator."""
         str_Lx = self._str_describing_oper()
         return (
-            "{} operator 1D,\n".format(oper_method)
-            + "nx = {0:6d}\n".format(self.nx)
+            f"{oper_method} operator 1D,\n"
+            + f"nx = {self.nx:6d}\n"
             + "Lx = "
             + str_Lx
             + "\n"
diff --git a/fluidsim/operators/op_finitediff1d.py b/fluidsim/operators/op_finitediff1d.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vb3BlcmF0b3JzL29wX2Zpbml0ZWRpZmYxZC5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vb3BlcmF0b3JzL29wX2Zpbml0ZWRpZmYxZC5weQ== 100644
--- a/fluidsim/operators/op_finitediff1d.py
+++ b/fluidsim/operators/op_finitediff1d.py
@@ -17,7 +17,7 @@
 
 class OperatorFiniteDiff1DPeriodic(OperatorsBase1D):
     def __init__(self, params=None):
-        super(OperatorFiniteDiff1DPeriodic, self).__init__(params)
+        super().__init__(params)
         nx = self.nx
         self.nx_seq = nx
         dx = self.deltax
@@ -45,6 +45,4 @@
         return sparse.identity(self.size)
 
     def produce_long_str_describing_oper(self):
-        return super(
-            OperatorFiniteDiff1DPeriodic, self
-        ).produce_long_str_describing_oper("Finite difference")
+        return super().produce_long_str_describing_oper("Finite difference")
diff --git a/fluidsim/operators/op_finitediff2d.py b/fluidsim/operators/op_finitediff2d.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vb3BlcmF0b3JzL29wX2Zpbml0ZWRpZmYyZC5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vb3BlcmF0b3JzL29wX2Zpbml0ZWRpZmYyZC5weQ== 100644
--- a/fluidsim/operators/op_finitediff2d.py
+++ b/fluidsim/operators/op_finitediff2d.py
@@ -137,7 +137,7 @@
         if (self.Lx / np.pi).is_integer():
             str_Lx = repr(int(self.Lx // np.pi)) + "pi"
         else:
-            str_Lx = "{:.3f}".format(self.Lx).rstrip("0")
+            str_Lx = f"{self.Lx:.3f}".rstrip("0")
         if (self.Ly / np.pi).is_integer():
             str_Ly = repr(int(self.Ly // np.pi)) + "pi"
         else:
@@ -141,7 +141,7 @@
         if (self.Ly / np.pi).is_integer():
             str_Ly = repr(int(self.Ly // np.pi)) + "pi"
         else:
-            str_Ly = "{:.3f}".format(self.Ly).rstrip("0")
+            str_Ly = f"{self.Ly:.3f}".rstrip("0")
         return ("L=" + str_Lx + "x" + str_Ly + "_{}x{}").format(self.nx, self.ny)
 
     def produce_long_str_describing_oper(self):
@@ -149,7 +149,7 @@
         if (self.Lx / np.pi).is_integer():
             str_Lx = repr(int(self.Lx // np.pi)) + "pi"
         else:
-            str_Lx = "{:.3f}".format(self.Lx).rstrip("0")
+            str_Lx = f"{self.Lx:.3f}".rstrip("0")
         if (self.Ly / np.pi).is_integer():
             str_Ly = repr(int(self.Ly // np.pi)) + "pi"
         else:
@@ -153,6 +153,6 @@
         if (self.Ly / np.pi).is_integer():
             str_Ly = repr(int(self.Ly // np.pi)) + "pi"
         else:
-            str_Ly = "{:.3f}".format(self.Ly).rstrip("0")
+            str_Ly = f"{self.Ly:.3f}".rstrip("0")
         return (
             "Finite difference operator 2D,\n"
@@ -157,6 +157,6 @@
         return (
             "Finite difference operator 2D,\n"
-            "nx = {0:6d} ; ny = {1:6d}\n".format(self.nx, self.ny)
+            "nx = {:6d} ; ny = {:6d}\n".format(self.nx, self.ny)
             + "Lx = "
             + str_Lx
             + " ; Ly = "
diff --git a/fluidsim/operators/operators0d.py b/fluidsim/operators/operators0d.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczBkLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczBkLnB5 100644
--- a/fluidsim/operators/operators0d.py
+++ b/fluidsim/operators/operators0d.py
@@ -12,7 +12,7 @@
 from fluiddyn.util import mpi
 
 
-class Operators0D(object):
+class Operators0D:
     """0D operators.
 
     """
diff --git a/fluidsim/operators/operators1d.py b/fluidsim/operators/operators1d.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczFkLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczFkLnB5 100644
--- a/fluidsim/operators/operators1d.py
+++ b/fluidsim/operators/operators1d.py
@@ -73,9 +73,7 @@
         self.K = self.KX
 
     def produce_long_str_describing_oper(self):
-        return super(
-            OperatorsPseudoSpectral1D, self
-        ).produce_long_str_describing_oper("Pseudospectral")
+        return super().produce_long_str_describing_oper("Pseudospectral")
 
     def dealiasing(self, *args):
         for thing in args:
diff --git a/fluidsim/operators/operators2d.py b/fluidsim/operators/operators2d.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczJkLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczJkLnB5 100644
--- a/fluidsim/operators/operators2d.py
+++ b/fluidsim/operators/operators2d.py
@@ -117,7 +117,7 @@
         if params.ONLY_COARSE_OPER:
             nx = ny = 4
 
-        super(OperatorsPseudoSpectral2D, self).__init__(
+        super().__init__(
             nx,
             ny,
             params.oper.Lx,
@@ -455,7 +455,7 @@
         """
         sign = 1j ** order
         if sign.imag != 0:
-            raise ValueError("Order={} should be even!".format(order))
+            raise ValueError(f"Order={order} should be even!")
 
         if negative:
             sign *= -1
@@ -483,7 +483,7 @@
         """
         sign = 1.0 / 1j ** order
         if sign.imag != 0:
-            raise ValueError("Order={} should be even!".format(order))
+            raise ValueError(f"Order={order} should be even!")
 
         if negative:
             sign *= -1
diff --git a/fluidsim/operators/operators3d.py b/fluidsim/operators/operators3d.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczNkLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vb3BlcmF0b3JzL29wZXJhdG9yczNkLnB5 100644
--- a/fluidsim/operators/operators3d.py
+++ b/fluidsim/operators/operators3d.py
@@ -179,7 +179,7 @@
             ny = params.oper.ny
             nz = params.oper.nz
 
-        super(OperatorsPseudoSpectral3D, self).__init__(
+        super().__init__(
             nx,
             ny,
             nz,
diff --git a/fluidsim/solvers/ad1d/init_fields.py b/fluidsim/solvers/ad1d/init_fields.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9hZDFkL2luaXRfZmllbGRzLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9hZDFkL2luaXRfZmllbGRzLnB5 100644
--- a/fluidsim/solvers/ad1d/init_fields.py
+++ b/fluidsim/solvers/ad1d/init_fields.py
@@ -21,7 +21,7 @@
 
     @classmethod
     def _complete_params_with_default(cls, params):
-        super(InitFieldsCos, cls)._complete_params_with_default(params)
+        super()._complete_params_with_default(params)
 
     # params.init_fields._set_child(cls.tag, attribs={})
 
@@ -36,7 +36,7 @@
 
     @classmethod
     def _complete_params_with_default(cls, params):
-        super(InitFieldsGaussian, cls)._complete_params_with_default(params)
+        super()._complete_params_with_default(params)
 
     # params.init_fields._set_child(cls.tag, attribs={})
 
diff --git a/fluidsim/solvers/ad1d/pseudo_spect/solver.py b/fluidsim/solvers/ad1d/pseudo_spect/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9hZDFkL3BzZXVkb19zcGVjdC9zb2x2ZXIucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9hZDFkL3BzZXVkb19zcGVjdC9zb2x2ZXIucHk= 100644
--- a/fluidsim/solvers/ad1d/pseudo_spect/solver.py
+++ b/fluidsim/solvers/ad1d/pseudo_spect/solver.py
@@ -19,7 +19,7 @@
 
 class InfoSolverAD1DPseudoSpect(InfoSolverPseudoSpectral):
     def _init_root(self):
-        super(InfoSolverAD1DPseudoSpect, self)._init_root()
+        super()._init_root()
 
         package = "fluidsim.solvers.ad1d.pseudo_spect"
         self.module_name = package + ".solver"
diff --git a/fluidsim/solvers/ad1d/solver.py b/fluidsim/solvers/ad1d/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9hZDFkL3NvbHZlci5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9hZDFkL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/ad1d/solver.py
+++ b/fluidsim/solvers/ad1d/solver.py
@@ -17,7 +17,7 @@
 class InfoSolverAD1D(InfoSolverFiniteDiff):
     def _init_root(self):
 
-        super(InfoSolverAD1D, self)._init_root()
+        super()._init_root()
 
         package = "fluidsim.solvers.ad1d"
         self.module_name = package + ".solver"
diff --git a/fluidsim/solvers/models0d/lorenz/output/print_stdout.py b/fluidsim/solvers/models0d/lorenz/output/print_stdout.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9tb2RlbHMwZC9sb3Jlbnovb3V0cHV0L3ByaW50X3N0ZG91dC5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9tb2RlbHMwZC9sb3Jlbnovb3V0cHV0L3ByaW50X3N0ZG91dC5weQ== 100644
--- a/fluidsim/solvers/models0d/lorenz/output/print_stdout.py
+++ b/fluidsim/solvers/models0d/lorenz/output/print_stdout.py
@@ -34,7 +34,7 @@
         self.print_stdout = self.__call__
 
     def _make_str_info(self):
-        to_print = super(PrintStdOutLorenz, self)._make_str_info()
+        to_print = super()._make_str_info()
 
         if mpi.rank == 0:
             to_print += (
diff --git a/fluidsim/solvers/models0d/predaprey/output/print_stdout.py b/fluidsim/solvers/models0d/predaprey/output/print_stdout.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9tb2RlbHMwZC9wcmVkYXByZXkvb3V0cHV0L3ByaW50X3N0ZG91dC5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9tb2RlbHMwZC9wcmVkYXByZXkvb3V0cHV0L3ByaW50X3N0ZG91dC5weQ== 100644
--- a/fluidsim/solvers/models0d/predaprey/output/print_stdout.py
+++ b/fluidsim/solvers/models0d/predaprey/output/print_stdout.py
@@ -35,7 +35,7 @@
         self.print_stdout = self.__call__
 
     def _make_str_info(self):
-        to_print = super(PrintStdOutPredaPrey, self)._make_str_info()
+        to_print = super()._make_str_info()
 
         potential = self.output.compute_potential()
         if mpi.rank == 0:
diff --git a/fluidsim/solvers/models0d/predaprey/solver.py b/fluidsim/solvers/models0d/predaprey/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9tb2RlbHMwZC9wcmVkYXByZXkvc29sdmVyLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9tb2RlbHMwZC9wcmVkYXByZXkvc29sdmVyLnB5 100644
--- a/fluidsim/solvers/models0d/predaprey/solver.py
+++ b/fluidsim/solvers/models0d/predaprey/solver.py
@@ -41,7 +41,7 @@
     """Contain the information on the solver predaprey."""
 
     def _init_root(self):
-        super(InfoSolverPredaPrey, self)._init_root()
+        super()._init_root()
 
         package = "fluidsim.solvers.models0d.predaprey"
         self.module_name = package + ".solver"
@@ -72,7 +72,7 @@
         params._set_attribs(attribs)
 
     def __init__(self, *args, **kargs):
-        super(Simul, self).__init__(*args, **kargs)
+        super().__init__(*args, **kargs)
         p = self.params
         self.Xs = p.C / p.D
         self.Ys = p.A / p.B
diff --git a/fluidsim/solvers/ns2d/bouss/solver.py b/fluidsim/solvers/ns2d/bouss/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL2JvdXNzL3NvbHZlci5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL2JvdXNzL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/ns2d/bouss/solver.py
+++ b/fluidsim/solvers/ns2d/bouss/solver.py
@@ -29,7 +29,7 @@
 class InfoSolverNS2DBouss(InfoSolverNS2D):
     def _init_root(self):
 
-        super(InfoSolverNS2DBouss, self)._init_root()
+        super()._init_root()
 
         package = "fluidsim.solvers.ns2d.bouss"
         self.module_name = package + ".solver"
diff --git a/fluidsim/solvers/ns2d/bouss/state.py b/fluidsim/solvers/ns2d/bouss/state.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL2JvdXNzL3N0YXRlLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL2JvdXNzL3N0YXRlLnB5 100644
--- a/fluidsim/solvers/ns2d/bouss/state.py
+++ b/fluidsim/solvers/ns2d/bouss/state.py
@@ -40,7 +40,7 @@
 
     def __init__(self, sim, oper=None):
 
-        super(StateNS2DBouss, self).__init__(sim, oper)
+        super().__init__(sim, oper)
 
         self.field_tmp4 = np.empty_like(self.state_phys[0])
         self.field_tmp5 = np.empty_like(self.state_phys[0])
diff --git a/fluidsim/solvers/ns2d/init_fields.py b/fluidsim/solvers/ns2d/init_fields.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL2luaXRfZmllbGRzLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL2luaXRfZmllbGRzLnB5 100644
--- a/fluidsim/solvers/ns2d/init_fields.py
+++ b/fluidsim/solvers/ns2d/init_fields.py
@@ -30,7 +30,7 @@
     @classmethod
     def _complete_params_with_default(cls, params):
         """Complete the `params` container (class method)."""
-        super(InitFieldsNoise, cls)._complete_params_with_default(params)
+        super()._complete_params_with_default(params)
 
         params.init_fields._set_child(
             cls.tag, attribs={"velo_max": 1.0, "length": 0.0}
@@ -117,7 +117,7 @@
     @classmethod
     def _complete_params_with_default(cls, params):
         """Complete the `params` container (class method)."""
-        super(InitFieldsJet, cls)._complete_params_with_default(params)
+        super()._complete_params_with_default(params)
 
     # params.init_fields._set_child(cls.tag, attribs={})
 
@@ -153,7 +153,7 @@
     @classmethod
     def _complete_params_with_default(cls, params):
         """Complete the `params` container (class method)."""
-        super(InitFieldsDipole, cls)._complete_params_with_default(params)
+        super()._complete_params_with_default(params)
 
     # params.init_fields._set_child(cls.tag, attribs={})
 
diff --git a/fluidsim/solvers/ns2d/output/print_stdout.py b/fluidsim/solvers/ns2d/output/print_stdout.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9wcmludF9zdGRvdXQucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9wcmludF9zdGRvdXQucHk= 100644
--- a/fluidsim/solvers/ns2d/output/print_stdout.py
+++ b/fluidsim/solvers/ns2d/output/print_stdout.py
@@ -23,7 +23,7 @@
     """
 
     def _make_str_info(self):
-        to_print = super(PrintStdOutNS2D, self)._make_str_info()
+        to_print = super()._make_str_info()
 
         energy = self.output.compute_energy()
 
@@ -102,7 +102,7 @@
         ax.set_title(
             "info stdout, solver "
             + self.output.name_solver
-            + ", nh = {0:5d}".format(self.sim.oper.nx_seq)
+            + f", nh = {self.sim.oper.nx_seq:5d}"
         )
         ax.plot(t, deltat, "k", linewidth=2)
         fig.tight_layout()
@@ -119,7 +119,7 @@
         ax.set_title(
             "info stdout, solver "
             + self.output.name_solver
-            + ", nh = {0:5d}".format(self.sim.oper.nx_seq)
+            + f", nh = {self.sim.oper.nx_seq:5d}"
         )
 
         ax.set_xlabel("t")
diff --git a/fluidsim/solvers/ns2d/output/spatial_means.py b/fluidsim/solvers/ns2d/output/spatial_means.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9zcGF0aWFsX21lYW5zLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9zcGF0aWFsX21lYW5zLnB5 100644
--- a/fluidsim/solvers/ns2d/output/spatial_means.py
+++ b/fluidsim/solvers/ns2d/output/spatial_means.py
@@ -66,7 +66,7 @@
         if mpi.rank == 0:
             epsK_tot = epsK + epsK_hypo
 
-            self.file.write("####\ntime = {:11.5e}\n".format(tsim))
+            self.file.write(f"####\ntime = {tsim:11.5e}\n")
             to_print = (
                 "E    = {:11.5e} ; Z         = {:11.5e} \n"
                 "epsK = {:11.5e} ; epsK_hypo = {:11.5e} ; epsK_tot = {:11.5e} \n"
diff --git a/fluidsim/solvers/ns2d/output/spect_energy_budget.py b/fluidsim/solvers/ns2d/output/spect_energy_budget.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9zcGVjdF9lbmVyZ3lfYnVkZ2V0LnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9zcGVjdF9lbmVyZ3lfYnVkZ2V0LnB5 100644
--- a/fluidsim/solvers/ns2d/output/spect_energy_budget.py
+++ b/fluidsim/solvers/ns2d/output/spect_energy_budget.py
@@ -117,7 +117,7 @@
             imin_plot = np.argmin(abs(times - tmin))
             imax_plot = np.argmin(abs(times - tmax))
 
-            to_print = "plot(tmin={0}, tmax={1}, delta_t={2:.2f})".format(
+            to_print = "plot(tmin={}, tmax={}, delta_t={:.2f})".format(
                 tmin, tmax, delta_t
             )
             print(to_print)
@@ -126,8 +126,8 @@
             tmax_plot = times[imax_plot]
             print(
                 """plot spectral energy budget
-    tmin = {0:8.6g} ; tmax = {1:8.6g} ; delta_t = {2:8.6g}
-    imin = {3:8d} ; imax = {4:8d} ; delta_i = {5:8d}""".format(
+    tmin = {:8.6g} ; tmax = {:8.6g} ; delta_t = {:8.6g}
+    imin = {:8d} ; imax = {:8d} ; delta_i = {:8d}""".format(
                     tmin_plot,
                     tmax_plot,
                     delta_t,
diff --git a/fluidsim/solvers/ns2d/output/spectra.py b/fluidsim/solvers/ns2d/output/spectra.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9zcGVjdHJhLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL291dHB1dC9zcGVjdHJhLnB5 100644
--- a/fluidsim/solvers/ns2d/output/spectra.py
+++ b/fluidsim/solvers/ns2d/output/spectra.py
@@ -77,6 +77,6 @@
             tmax_plot = times[imax_plot]
 
             print(
-                "plot1d(tmin={0}, tmax={1}, delta_t={2:.2f},".format(
+                "plot1d(tmin={}, tmax={}, delta_t={:.2f},".format(
                     tmin, tmax, delta_t
                 )
@@ -81,7 +81,7 @@
                     tmin, tmax, delta_t
                 )
-                + " coef_compensate={0:.3f})".format(coef_compensate)
+                + f" coef_compensate={coef_compensate:.3f})"
             )
 
             print(
                 """plot 1D spectra
@@ -84,9 +84,9 @@
             )
 
             print(
                 """plot 1D spectra
-    tmin = {0:8.6g} ; tmax = {1:8.6g} ; delta_t = {2:8.6g}
-    imin = {3:8d} ; imax = {4:8d} ; delta_i = {5:8d}""".format(
+    tmin = {:8.6g} ; tmax = {:8.6g} ; delta_t = {:8.6g}
+    imin = {:8d} ; imax = {:8d} ; delta_i = {:8d}""".format(
                     tmin_plot,
                     tmax_plot,
                     delta_t,
@@ -102,7 +102,7 @@
             ax1.set_title(
                 "1D spectra, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
+                + f", nh = {self.nx:5d}"
             )
             ax1.set_xscale("log")
             ax1.set_yscale("log")
@@ -157,6 +157,6 @@
             tmax_plot = times[imax_plot]
 
             print(
-                "plot2s(tmin={0}, tmax={1}, delta_t={2:.2f},".format(
+                "plot2s(tmin={}, tmax={}, delta_t={:.2f},".format(
                     tmin, tmax, delta_t
                 )
@@ -161,7 +161,7 @@
                     tmin, tmax, delta_t
                 )
-                + " coef_compensate={0:.3f})".format(coef_compensate)
+                + f" coef_compensate={coef_compensate:.3f})"
             )
 
             print(
                 """plot 2D spectra
@@ -164,9 +164,9 @@
             )
 
             print(
                 """plot 2D spectra
-    tmin = {0:8.6g} ; tmax = {1:8.6g} ; delta_t = {2:8.6g}
-    imin = {3:8d} ; imax = {4:8d} ; delta_i = {5:8d}""".format(
+    tmin = {:8.6g} ; tmax = {:8.6g} ; delta_t = {:8.6g}
+    imin = {:8d} ; imax = {:8d} ; delta_i = {:8d}""".format(
                     tmin_plot,
                     tmax_plot,
                     delta_t,
@@ -182,7 +182,7 @@
             ax1.set_title(
                 "2D spectra, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
+                + f", nh = {self.nx:5d}"
             )
             ax1.set_xscale("log")
             ax1.set_yscale("log")
diff --git a/fluidsim/solvers/ns2d/solver.py b/fluidsim/solvers/ns2d/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL3NvbHZlci5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/ns2d/solver.py
+++ b/fluidsim/solvers/ns2d/solver.py
@@ -70,7 +70,7 @@
         - :class:`fluidsim.solvers.ns2d.forcing.ForcingNS2D`
 
         """
-        super(InfoSolverNS2D, self)._init_root()
+        super()._init_root()
 
         package = "fluidsim.solvers.ns2d"
         self.module_name = package + ".solver"
diff --git a/fluidsim/solvers/ns2d/state.py b/fluidsim/solvers/ns2d/state.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0YXRlLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0YXRlLnB5 100644
--- a/fluidsim/solvers/ns2d/state.py
+++ b/fluidsim/solvers/ns2d/state.py
@@ -39,7 +39,7 @@
 
     def __init__(self, sim, oper=None):
 
-        super(StateNS2D, self).__init__(sim, oper)
+        super().__init__(sim, oper)
 
         self.field_tmp0 = np.empty_like(self.state_phys[0])
         self.field_tmp1 = np.empty_like(self.state_phys[0])
@@ -141,5 +141,5 @@
             elif key == "uy_fft":
                 self.init_from_uyfft(arr)
             else:
-                super(StateNS2D, self).init_statespect_from(**kwargs)
+                super().init_statespect_from(**kwargs)
         else:
@@ -145,2 +145,2 @@
         else:
-            super(StateNS2D, self).init_statespect_from(**kwargs)
+            super().init_statespect_from(**kwargs)
diff --git a/fluidsim/solvers/ns2d/strat/init_fields.py b/fluidsim/solvers/ns2d/strat/init_fields.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L2luaXRfZmllbGRzLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L2luaXRfZmllbGRzLnB5 100644
--- a/fluidsim/solvers/ns2d/strat/init_fields.py
+++ b/fluidsim/solvers/ns2d/strat/init_fields.py
@@ -51,7 +51,7 @@
     @classmethod
     def _complete_params_with_default(cls, params):
         """Complete the `params` container (class method)."""
-        super(InitFieldsLinearMode, cls)._complete_params_with_default(params)
+        super()._complete_params_with_default(params)
         params.init_fields._set_child(
             cls.tag,
             attribs={
diff --git a/fluidsim/solvers/ns2d/strat/output/__init__.py b/fluidsim/solvers/ns2d/strat/output/__init__.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9fX2luaXRfXy5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9fX2luaXRfXy5weQ== 100644
--- a/fluidsim/solvers/ns2d/strat/output/__init__.py
+++ b/fluidsim/solvers/ns2d/strat/output/__init__.py
@@ -36,7 +36,7 @@
             self._init_froude_number()
             self.ratio_omegas = self._compute_ratio_omegas()
 
-        super(OutputStrat, self).__init__(sim)
+        super().__init__(sim)
 
     @staticmethod
     def _complete_info_solver(info_solver):
@@ -146,7 +146,7 @@
         """Compute froude number ONLY for anisotropic forcing."""
         angle = self.sim.params.forcing.tcrandom_anisotropic.angle
         if isinstance(angle, str):
-            if angle.endswith(u"°"):
+            if angle.endswith("°"):
                 angle = radians(float(angle[:-1]))
             else:
                 raise ValueError(
@@ -158,7 +158,7 @@
         self.froude_number = np.sin(angle)
 
     def _compute_ratio_omegas(self):
-        """Compute ratio omegas; gamma = \omega_l / \omega_{af}"""
+        r"""Compute ratio omegas; gamma = \omega_l / \omega_{af}"""
         params = self.sim.params
 
         omega_l = params.N * self.froude_number
@@ -208,7 +208,7 @@
 
     def _create_list_for_name_run(self):
         """Creates new name_run for the simulation."""
-        list_for_name_run = super(OutputStrat, self)._create_list_for_name_run()
+        list_for_name_run = super()._create_list_for_name_run()
         if self.sim.params.forcing.type.endswith("anisotropic"):
             str_describing_attribs_strat = (
                 self._produce_str_describing_attribs_strat()
diff --git a/fluidsim/solvers/ns2d/strat/output/phys_fields.py b/fluidsim/solvers/ns2d/strat/output/phys_fields.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9waHlzX2ZpZWxkcy5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9waHlzX2ZpZWxkcy5weQ== 100644
--- a/fluidsim/solvers/ns2d/strat/output/phys_fields.py
+++ b/fluidsim/solvers/ns2d/strat/output/phys_fields.py
@@ -17,8 +17,8 @@
 
     def update_animation(self, frame, **fargs):
 
-        super(PhysFields2DStrat, self).update_animation(frame, **fargs)
+        super().update_animation(frame, **fargs)
 
         if "ratio_omegas" in self.__dict__:
             title = (
                 self.key_field
@@ -21,9 +21,9 @@
 
         if "ratio_omegas" in self.__dict__:
             title = (
                 self.key_field
-                + ", R = {0:.2f}".format(self.output.ratio_omegas)
-                + ", F = {0:.1f}".format(self.output.froude_number)
-                + ", $t = {0:.3f}$".format(self.ani_times[frame])
+                + f", R = {self.output.ratio_omegas:.2f}"
+                + f", F = {self.output.froude_number:.1f}"
+                + ", $t = {:.3f}$".format(self.ani_times[frame])
             )
         else:
@@ -28,6 +28,6 @@
             )
         else:
-            title = self.key_field + ", $t = {0:.3f}$".format(
+            title = self.key_field + ", $t = {:.3f}$".format(
                 self.ani_times[frame]
             )
 
diff --git a/fluidsim/solvers/ns2d/strat/output/print_stdout.py b/fluidsim/solvers/ns2d/strat/output/print_stdout.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9wcmludF9zdGRvdXQucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9wcmludF9zdGRvdXQucHk= 100644
--- a/fluidsim/solvers/ns2d/strat/output/print_stdout.py
+++ b/fluidsim/solvers/ns2d/strat/output/print_stdout.py
@@ -26,7 +26,7 @@
     """
 
     def __init__(self, output):
-        super(PrintStdOutNS2DStrat, self).__init__(output)
+        super().__init__(output)
         self.path_memory = self.output.path_run + "/memory_out.txt"
         if mpi.rank == 0 and self.output._has_to_save:
             if not os.path.exists(self.path_memory):
@@ -36,7 +36,7 @@
                 self.file_memory.seek(0, 2)  # go to the end of the file
 
     def _make_str_info(self):
-        to_print = super(PrintStdOutNS2DStrat, self)._make_str_info()
+        to_print = super()._make_str_info()
 
         energyK, energyA, energyK_ux = self.output.compute_energies()
         energy = energyK + energyA
@@ -52,7 +52,7 @@
                 memory = get_memory_usage()
                 self._write_memory_txt()
 
-                to_print += "              memory  = {:9.3f} Mo.\n".format(memory)
+                to_print += f"              memory  = {memory:9.3f} Mo.\n"
 
             duration_left = self._evaluate_duration_left()
             if duration_left is not None:
@@ -65,7 +65,7 @@
         """Write memory .txt"""
         it = self.sim.time_stepping.it
         mem = get_memory_usage()
-        self.file_memory.write("{:.3f},{:.3f}\n".format(it, mem))
+        self.file_memory.write(f"{it:.3f},{mem:.3f}\n")
         self.file_memory.flush()
         os.fsync(self.file_memory.fileno())
 
@@ -172,7 +172,7 @@
         ax1.set_title(
             "info stdout, solver "
             + self.output.name_solver
-            + ", nh = {0:5d}".format(self.params.oper.nx)
+            + f", nh = {self.params.oper.nx:5d}"
         )
         ax1.plot(t, deltat, "k", linewidth=2)
 
@@ -188,7 +188,7 @@
         ax2.grid(True)
 
     def close(self):
-        super(PrintStdOutNS2DStrat, self).close()
+        super().close()
         try:
             self.file_memory.close()
         except AttributeError:
diff --git a/fluidsim/solvers/ns2d/strat/output/spatial_means.py b/fluidsim/solvers/ns2d/strat/output/spatial_means.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9zcGF0aWFsX21lYW5zLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9zcGF0aWFsX21lYW5zLnB5 100644
--- a/fluidsim/solvers/ns2d/strat/output/spatial_means.py
+++ b/fluidsim/solvers/ns2d/strat/output/spatial_means.py
@@ -87,7 +87,7 @@
         if mpi.rank == 0:
             epsK_tot = epsK + epsK_hypo
 
-            self.file.write("####\ntime = {:11.5e}\n".format(tsim))
+            self.file.write(f"####\ntime = {tsim:11.5e}\n")
             to_print = (
                 "Z    = {:11.5e} \n"
                 "E    = {:11.5e} ; EK         = {:11.5e} ; EA        = {:11.5e} \n"
diff --git a/fluidsim/solvers/ns2d/strat/output/spect_energy_budget.py b/fluidsim/solvers/ns2d/strat/output/spect_energy_budget.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9zcGVjdF9lbmVyZ3lfYnVkZ2V0LnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9zcGVjdF9lbmVyZ3lfYnVkZ2V0LnB5 100644
--- a/fluidsim/solvers/ns2d/strat/output/spect_energy_budget.py
+++ b/fluidsim/solvers/ns2d/strat/output/spect_energy_budget.py
@@ -235,7 +235,7 @@
         imin_plot = np.argmin(abs(times - tmin))
         imax_plot = np.argmin(abs(times - tmax))
 
-        to_print = "plot(tmin={0}, tmax={1}, delta_t={2:.2f})".format(
+        to_print = "plot(tmin={}, tmax={}, delta_t={:.2f})".format(
             tmin, tmax, delta_t
         )
         print(to_print)
@@ -244,8 +244,8 @@
         tmax_plot = times[imax_plot]
         print(
             """plot spectral energy budget
-            tmin = {0:8.6g} ; tmax = {1:8.6g} ; delta_t = {2:8.6g}
-            imin = {3:8d} ; imax = {4:8d} ; delta_i = {5:8d}""".format(
+            tmin = {:8.6g} ; tmax = {:8.6g} ; delta_t = {:8.6g}
+            imin = {:8d} ; imax = {:8d} ; delta_i = {:8d}""".format(
                 tmin_plot, tmax_plot, delta_t, imin_plot, imax_plot, delta_i_plot
             )
         )
@@ -259,7 +259,7 @@
         ax1.set_title(
             "2D spectra, solver "
             + self.output.name_solver
-            + ", nh = {0:5d}".format(self.nx)
+            + f", nh = {self.nx:5d}"
         )
 
         transferEK_kx = dset_transferEK_kx[imin_plot : imax_plot + 1].mean(0)
@@ -298,7 +298,7 @@
         ax2.set_title(
             "2D spectra, solver "
             + self.output.name_solver
-            + ", nh = {0:5d}".format(self.nx)
+            + f", nh = {self.nx:5d}"
         )
 
         transferEK_ky = dset_transferEK_ky[imin_plot : imax_plot + 1].mean(0)
diff --git a/fluidsim/solvers/ns2d/strat/output/spectra.py b/fluidsim/solvers/ns2d/strat/output/spectra.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9zcGVjdHJhLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9zcGVjdHJhLnB5 100644
--- a/fluidsim/solvers/ns2d/strat/output/spectra.py
+++ b/fluidsim/solvers/ns2d/strat/output/spectra.py
@@ -137,6 +137,6 @@
             tmax_plot = times[imax_plot]
 
             print(
-                "plot1d(tmin={0}, tmax={1}, delta_t={2:.2f},".format(
+                "plot1d(tmin={}, tmax={}, delta_t={:.2f},".format(
                     tmin, tmax, delta_t
                 )
@@ -141,7 +141,7 @@
                     tmin, tmax, delta_t
                 )
-                + " coef_compensate_kx={0:.3f})".format(coef_compensate_kx)
+                + f" coef_compensate_kx={coef_compensate_kx:.3f})"
             )
 
             print(
                 """plot 1D spectra
@@ -144,9 +144,9 @@
             )
 
             print(
                 """plot 1D spectra
-            tmin = {0:8.6g} ; tmax = {1:8.6g} ; delta_t = {2:8.6g}
-            imin = {3:8d} ; imax = {4:8d} ; delta_i = {5:8d}""".format(
+            tmin = {:8.6g} ; tmax = {:8.6g} ; delta_t = {:8.6g}
+            imin = {:8d} ; imax = {:8d} ; delta_i = {:8d}""".format(
                     tmin_plot,
                     tmax_plot,
                     delta_t,
@@ -165,7 +165,7 @@
             ax1.set_title(
                 "1D spectra, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
+                + f", nh = {self.nx:5d}"
             )
             ax1.set_xscale("log")
             ax1.set_yscale("log")
@@ -235,7 +235,7 @@
             # Parameters figure E(k_y)
             fig, ax2 = self.output.figure_axe()
             ax2.set_xlabel(r"$k_z$")
-            ax2.set_ylabel(r"$E(k_z)k_z^{{{}}}$".format(coef_compensate_kz))
+            ax2.set_ylabel(fr"$E(k_z)k_z^{{{coef_compensate_kz}}}$")
             ax2.set_title(
                 "1D spectra, solver "
                 + self.output.name_solver
@@ -239,7 +239,7 @@
             ax2.set_title(
                 "1D spectra, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
+                + f", nh = {self.nx:5d}"
             )
 
             ax2.set_xscale("log")
@@ -328,6 +328,6 @@
             tmax_plot = times[imax_plot]
 
             print(
-                "plot2s(tmin={0}, tmax={1}, delta_t={2:.2f},".format(
+                "plot2s(tmin={}, tmax={}, delta_t={:.2f},".format(
                     tmin, tmax, delta_t
                 )
@@ -332,7 +332,7 @@
                     tmin, tmax, delta_t
                 )
-                + " coef_compensate={0:.3f})".format(coef_compensate)
+                + f" coef_compensate={coef_compensate:.3f})"
             )
 
             print(
                 """plot 2D spectra
@@ -335,9 +335,9 @@
             )
 
             print(
                 """plot 2D spectra
-            tmin = {0:8.6g} ; tmax = {1:8.6g} ; delta_t = {2:8.6g}
-            imin = {3:8d} ; imax = {4:8d} ; delta_i = {5:8d}""".format(
+            tmin = {:8.6g} ; tmax = {:8.6g} ; delta_t = {:8.6g}
+            imin = {:8d} ; imax = {:8d} ; delta_i = {:8d}""".format(
                     tmin_plot,
                     tmax_plot,
                     delta_t,
@@ -354,7 +354,7 @@
             ax1.set_title(
                 "2D spectra, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
+                + f", nh = {self.nx:5d}"
             )
 
             ax1.set_xscale("log")
diff --git a/fluidsim/solvers/ns2d/strat/output/spectra_multidim.py b/fluidsim/solvers/ns2d/strat/output/spectra_multidim.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9zcGVjdHJhX211bHRpZGltLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L291dHB1dC9zcGVjdHJhX211bHRpZGltLnB5 100644
--- a/fluidsim/solvers/ns2d/strat/output/spectra_multidim.py
+++ b/fluidsim/solvers/ns2d/strat/output/spectra_multidim.py
@@ -249,7 +249,7 @@
         ax.text(
             kx[ikx_text],
             kz[ikz_text],
-            "{}".format(text_plot),
+            f"{text_plot}",
             color="white",
             fontsize=15,
         )
diff --git a/fluidsim/solvers/ns2d/strat/solver.py b/fluidsim/solvers/ns2d/strat/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L3NvbHZlci5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/ns2d/strat/solver.py
+++ b/fluidsim/solvers/ns2d/strat/solver.py
@@ -30,7 +30,7 @@
 class InfoSolverNS2DStrat(InfoSolverNS2D):
     def _init_root(self):
 
-        super(InfoSolverNS2DStrat, self)._init_root()
+        super()._init_root()
 
         package = "fluidsim.solvers.ns2d.strat"
         self.module_name = package + ".solver"
@@ -212,7 +212,7 @@
         energy_conserved = ratio < epsilon
         if not energy_conserved:
             print("Energy is not conserved!")
-            print("ratio = {}".format(ratio))
+            print(f"ratio = {ratio}")
             return False
         return True
 
diff --git a/fluidsim/solvers/ns2d/strat/state.py b/fluidsim/solvers/ns2d/strat/state.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L3N0YXRlLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L3N0YXRlLnB5 100644
--- a/fluidsim/solvers/ns2d/strat/state.py
+++ b/fluidsim/solvers/ns2d/strat/state.py
@@ -49,7 +49,7 @@
 
     def __init__(self, sim, oper=None):
 
-        super(StateNS2DStrat, self).__init__(sim, oper)
+        super().__init__(sim, oper)
 
         self.field_tmp4 = np.empty_like(self.state_phys[0])
         self.field_tmp5 = np.empty_like(self.state_phys[0])
diff --git a/fluidsim/solvers/ns2d/strat/time_stepping.py b/fluidsim/solvers/ns2d/strat/time_stepping.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L3RpbWVfc3RlcHBpbmcucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L3RpbWVfc3RlcHBpbmcucHk= 100644
--- a/fluidsim/solvers/ns2d/strat/time_stepping.py
+++ b/fluidsim/solvers/ns2d/strat/time_stepping.py
@@ -25,7 +25,7 @@
 
     @classmethod
     def _complete_params_with_default(cls, params):
-        super(TimeSteppingPseudoSpectralStrat, cls)._complete_params_with_default(
+        super()._complete_params_with_default(
             params
         )
 
@@ -36,7 +36,7 @@
         """
         Initialization compute time step solver ns2d.strat.
         """
-        super(TimeSteppingPseudoSpectralStrat, self)._init_compute_time_step()
+        super()._init_compute_time_step()
 
         # Coefficients dt
         self.coef_deltat_dispersion_relation = 1.0
@@ -204,5 +204,5 @@
         # np.isnan(np.sum seems to be really fast
         if np.isnan(np.sum(self.sim.state.state_spect[0])):
             raise ValueError(
-                "nan at it = {0}, t = {1:.4f}".format(self.it, self.t)
+                f"nan at it = {self.it}, t = {self.t:.4f}"
             )
diff --git a/fluidsim/solvers/ns3d/bouss/solver.py b/fluidsim/solvers/ns3d/bouss/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczNkL2JvdXNzL3NvbHZlci5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczNkL2JvdXNzL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/ns3d/bouss/solver.py
+++ b/fluidsim/solvers/ns3d/bouss/solver.py
@@ -23,7 +23,7 @@
 class InfoSolverNS3DBouss(InfoSolverNS3DStrat):
     def _init_root(self):
 
-        super(InfoSolverNS3DBouss, self)._init_root()
+        super()._init_root()
 
         package = "fluidsim.solvers.ns3d.bouss"
         self.module_name = package + ".solver"
diff --git a/fluidsim/solvers/ns3d/init_fields.py b/fluidsim/solvers/ns3d/init_fields.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczNkL2luaXRfZmllbGRzLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczNkL2luaXRfZmllbGRzLnB5 100644
--- a/fluidsim/solvers/ns3d/init_fields.py
+++ b/fluidsim/solvers/ns3d/init_fields.py
@@ -24,7 +24,7 @@
 
     @classmethod
     def _complete_params_with_default(cls, params):
-        super(InitFieldsDipole, cls)._complete_params_with_default(params)
+        super()._complete_params_with_default(params)
 
     # params.init_fields._set_child(cls.tag, attribs={'U': 1.})
 
@@ -78,7 +78,7 @@
     @classmethod
     def _complete_params_with_default(cls, params):
         """Complete the `params` container (class method)."""
-        super(InitFieldsNoise, cls)._complete_params_with_default(params)
+        super()._complete_params_with_default(params)
 
         params.init_fields._set_child(
             cls.tag, attribs={"velo_max": 1.0, "length": None}
diff --git a/fluidsim/solvers/ns3d/output/print_stdout.py b/fluidsim/solvers/ns3d/output/print_stdout.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczNkL291dHB1dC9wcmludF9zdGRvdXQucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczNkL291dHB1dC9wcmludF9zdGRvdXQucHk= 100644
--- a/fluidsim/solvers/ns3d/output/print_stdout.py
+++ b/fluidsim/solvers/ns3d/output/print_stdout.py
@@ -12,7 +12,7 @@
     """
 
     def _make_str_info(self):
-        to_print = super(PrintStdOutNS3D, self)._make_str_info()
+        to_print = super()._make_str_info()
 
         energy = self.output.compute_energy()
         if mpi.rank == 0:
@@ -93,7 +93,7 @@
         ax1.set_title(
             "info stdout, solver "
             + self.output.name_solver
-            + ", nx = {0:5d}".format(self.params.oper.nx)
+            + f", nx = {self.params.oper.nx:5d}"
         )
         ax1.plot(t, deltat, "k", linewidth=2)
 
diff --git a/fluidsim/solvers/ns3d/output/spatial_means.py b/fluidsim/solvers/ns3d/output/spatial_means.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczNkL291dHB1dC9zcGF0aWFsX21lYW5zLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczNkL291dHB1dC9zcGF0aWFsX21lYW5zLnB5 100644
--- a/fluidsim/solvers/ns3d/output/spatial_means.py
+++ b/fluidsim/solvers/ns3d/output/spatial_means.py
@@ -62,7 +62,7 @@
 
         if mpi.rank == 0:
 
-            self.file.write("####\ntime = {:11.5e}\n".format(tsim))
+            self.file.write(f"####\ntime = {tsim:11.5e}\n")
             to_print = (
                 "E    = {:11.5e}\n"
                 "Ex   = {:11.5e} ; Ey   = {:11.5e} ; Ez   = {:11.5e}\n"
diff --git a/fluidsim/solvers/ns3d/output/spectra.py b/fluidsim/solvers/ns3d/output/spectra.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczNkL291dHB1dC9zcGVjdHJhLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczNkL291dHB1dC9zcGVjdHJhLnB5 100644
--- a/fluidsim/solvers/ns3d/output/spectra.py
+++ b/fluidsim/solvers/ns3d/output/spectra.py
@@ -93,8 +93,8 @@
             tmax_plot = times[imax_plot]
 
             print(
-                "plot1d(tmin={}, tmax={}, delta_t={},".format(tmin, tmax, delta_t)
-                + " coef_compensate={:.3f})".format(coef_compensate)
+                f"plot1d(tmin={tmin}, tmax={tmax}, delta_t={delta_t},"
+                + f" coef_compensate={coef_compensate:.3f})"
             )
 
             print(
@@ -111,8 +111,8 @@
             )
 
             fig, ax = self.output.figure_axe()
-            ax.set_xlabel("${}$".format(key_k))
+            ax.set_xlabel(f"${key_k}$")
             ax.set_ylabel("spectra " + key)
             ax.set_title(
                 "1D spectra, solver "
                 + self.output.name_solver
@@ -115,8 +115,8 @@
             ax.set_ylabel("spectra " + key)
             ax.set_title(
                 "1D spectra, solver "
                 + self.output.name_solver
-                + ", nx = {0:5d}".format(self.nx)
+                + f", nx = {self.nx:5d}"
             )
             ax.set_xscale("log")
             ax.set_yscale("log")
@@ -178,6 +178,6 @@
             tmax_plot = times[imax_plot]
 
             print(
-                "plot3d(tmin={0}, tmax={1}, delta_t={2:.2f},".format(
+                "plot3d(tmin={}, tmax={}, delta_t={:.2f},".format(
                     tmin, tmax, delta_t
                 )
@@ -182,7 +182,7 @@
                     tmin, tmax, delta_t
                 )
-                + " coef_compensate={0:.3f})".format(coef_compensate)
+                + f" coef_compensate={coef_compensate:.3f})"
             )
 
             print(
                 """plot 3d spectra
@@ -185,9 +185,9 @@
             )
 
             print(
                 """plot 3d spectra
-    tmin = {0:8.6g} ; tmax = {1:8.6g} ; delta_t = {2:8.6g}
-    imin = {3:8d} ; imax = {4:8d} ; delta_i = {5:8d}""".format(
+    tmin = {:8.6g} ; tmax = {:8.6g} ; delta_t = {:8.6g}
+    imin = {:8d} ; imax = {:8d} ; delta_i = {:8d}""".format(
                     tmin_plot,
                     tmax_plot,
                     delta_t,
@@ -203,7 +203,7 @@
             ax1.set_title(
                 "3D spectra, solver "
                 + self.output.name_solver
-                + ", nx = {0:5d}".format(self.nx)
+                + f", nx = {self.nx:5d}"
             )
             ax1.set_xscale("log")
             ax1.set_yscale("log")
diff --git a/fluidsim/solvers/ns3d/solver.py b/fluidsim/solvers/ns3d/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczNkL3NvbHZlci5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczNkL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/ns3d/solver.py
+++ b/fluidsim/solvers/ns3d/solver.py
@@ -27,7 +27,7 @@
 class InfoSolverNS3D(InfoSolverPseudoSpectral3D):
     def _init_root(self):
 
-        super(InfoSolverNS3D, self)._init_root()
+        super()._init_root()
 
         package = "fluidsim.solvers.ns3d"
         self.module_name = package + ".solver"
diff --git a/fluidsim/solvers/ns3d/state.py b/fluidsim/solvers/ns3d/state.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczNkL3N0YXRlLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczNkL3N0YXRlLnB5 100644
--- a/fluidsim/solvers/ns3d/state.py
+++ b/fluidsim/solvers/ns3d/state.py
@@ -36,7 +36,7 @@
 
     def __init__(self, sim, oper=None):
 
-        super(StateNS3D, self).__init__(sim, oper)
+        super().__init__(sim, oper)
 
         self.fields_tmp = tuple(
             np.empty_like(self.state_phys[0]) for n in range(6)
diff --git a/fluidsim/solvers/ns3d/strat/output/spatial_means.py b/fluidsim/solvers/ns3d/strat/output/spatial_means.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczNkL3N0cmF0L291dHB1dC9zcGF0aWFsX21lYW5zLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczNkL3N0cmF0L291dHB1dC9zcGF0aWFsX21lYW5zLnB5 100644
--- a/fluidsim/solvers/ns3d/strat/output/spatial_means.py
+++ b/fluidsim/solvers/ns3d/strat/output/spatial_means.py
@@ -21,7 +21,7 @@
 
     def __init__(self, output):
         self.one_over_N2 = 1.0 / output.sim.params.N ** 2
-        super(SpatialMeansNS3DStrat, self).__init__(output)
+        super().__init__(output)
 
     def _save_one_time(self):
         tsim = self.sim.time_stepping.t
@@ -77,7 +77,7 @@
 
         if mpi.rank == 0:
 
-            self.file.write("####\ntime = {:11.5e}\n".format(tsim))
+            self.file.write(f"####\ntime = {tsim:11.5e}\n")
             to_print = (
                 "E    = {:11.5e}\n"
                 "EA   = {:11.5e} ; EKz   = {:11.5e} ; "
diff --git a/fluidsim/solvers/ns3d/strat/solver.py b/fluidsim/solvers/ns3d/strat/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczNkL3N0cmF0L3NvbHZlci5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczNkL3N0cmF0L3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/ns3d/strat/solver.py
+++ b/fluidsim/solvers/ns3d/strat/solver.py
@@ -23,7 +23,7 @@
 class InfoSolverNS3DStrat(InfoSolverNS3D):
     def _init_root(self):
 
-        super(InfoSolverNS3DStrat, self)._init_root()
+        super()._init_root()
 
         package = "fluidsim.solvers.ns3d.strat"
         self.module_name = package + ".solver"
diff --git a/fluidsim/solvers/ns3d/time_stepping.py b/fluidsim/solvers/ns3d/time_stepping.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9uczNkL3RpbWVfc3RlcHBpbmcucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9uczNkL3RpbWVfc3RlcHBpbmcucHk= 100644
--- a/fluidsim/solvers/ns3d/time_stepping.py
+++ b/fluidsim/solvers/ns3d/time_stepping.py
@@ -20,5 +20,5 @@
         # np.isnan(np.sum seems to be really fast
         if np.isnan(np.sum(state_spect[0])):
             raise ValueError(
-                "nan at it = {0}, t = {1:.4f}".format(self.it, self.t)
+                f"nan at it = {self.it}, t = {self.t:.4f}"
             )
diff --git a/fluidsim/solvers/plate2d/dimensional.py b/fluidsim/solvers/plate2d/dimensional.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL2RpbWVuc2lvbmFsLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL2RpbWVuc2lvbmFsLnB5 100644
--- a/fluidsim/solvers/plate2d/dimensional.py
+++ b/fluidsim/solvers/plate2d/dimensional.py
@@ -14,7 +14,7 @@
 import numpy as np
 
 
-class Converter(object):
+class Converter:
     r"""Converter between dimensional and adimensional variables.
 
     Parameters
diff --git a/fluidsim/solvers/plate2d/operators.py b/fluidsim/solvers/plate2d/operators.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL29wZXJhdG9ycy5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL29wZXJhdG9ycy5weQ== 100644
--- a/fluidsim/solvers/plate2d/operators.py
+++ b/fluidsim/solvers/plate2d/operators.py
@@ -43,7 +43,7 @@
     """
 
     def __init__(self, params):
-        super(OperatorsPseudoSpectralPlate2D, self).__init__(params)
+        super().__init__(params)
 
         self.KX2 = self.KX ** 2
         self.KY2 = self.KY ** 2
diff --git a/fluidsim/solvers/plate2d/output/__init__.py b/fluidsim/solvers/plate2d/output/__init__.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL291dHB1dC9fX2luaXRfXy5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL291dHB1dC9fX2luaXRfXy5weQ== 100644
--- a/fluidsim/solvers/plate2d/output/__init__.py
+++ b/fluidsim/solvers/plate2d/output/__init__.py
@@ -84,6 +84,6 @@
         params.output.phys_fields.field_to_plot = "z"
 
     def _create_list_for_name_run(self):
-        list_for_name_run = super(Output, self)._create_list_for_name_run()
+        list_for_name_run = super()._create_list_for_name_run()
 
         if self.sim.params.forcing.enable:
@@ -88,6 +88,6 @@
 
         if self.sim.params.forcing.enable:
-            str_P = "P{:5.0e}".format(self.sim.params.forcing.forcing_rate)
+            str_P = f"P{self.sim.params.forcing.forcing_rate:5.0e}"
             str_P = str_P.replace("+", "")
             list_for_name_run.insert(2, str_P)
 
diff --git a/fluidsim/solvers/plate2d/output/correlations_freq.py b/fluidsim/solvers/plate2d/output/correlations_freq.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL291dHB1dC9jb3JyZWxhdGlvbnNfZnJlcS5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL291dHB1dC9jb3JyZWxhdGlvbnNfZnJlcS5weQ== 100644
--- a/fluidsim/solvers/plate2d/output/correlations_freq.py
+++ b/fluidsim/solvers/plate2d/output/correlations_freq.py
@@ -177,7 +177,7 @@
     def __init__(self, output):
         params = output.sim.params
         pcorrel_freq = params.output.correl_freq
-        super(CorrelationsFreq, self).__init__(
+        super().__init__(
             output,
             period_save=params.output.periods_save.correl_freq,
             has_to_plot_saved=pcorrel_freq.HAS_TO_PLOT_SAVED,
@@ -347,7 +347,7 @@
             ax.set_title(
                 "Correlation, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.params.oper.nx)
+                + f", nh = {self.params.oper.nx:5d}"
             )
 
     def _online_plot_saving(self, dict_results):
diff --git a/fluidsim/solvers/plate2d/output/spatial_means.py b/fluidsim/solvers/plate2d/output/spatial_means.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL291dHB1dC9zcGF0aWFsX21lYW5zLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL291dHB1dC9zcGF0aWFsX21lYW5zLnB5 100644
--- a/fluidsim/solvers/plate2d/output/spatial_means.py
+++ b/fluidsim/solvers/plate2d/output/spatial_means.py
@@ -103,7 +103,7 @@
             energy = energy_k + energy_l + energy_e
             epsK_tot = epsK + epsK_hypo
 
-            self.file.write("####\ntime = {0:17.13f}\n".format(tsim))
+            self.file.write(f"####\ntime = {tsim:17.13f}\n")
             to_print = (
                 "E    = {:31.26e} ; E_k    = {:11.6e} ; "
                 "E_l    = {:11.6e} ; E_e    = {:11.6e}\n"
diff --git a/fluidsim/solvers/plate2d/output/spectra.py b/fluidsim/solvers/plate2d/output/spectra.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL291dHB1dC9zcGVjdHJhLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL291dHB1dC9zcGVjdHJhLnB5 100644
--- a/fluidsim/solvers/plate2d/output/spectra.py
+++ b/fluidsim/solvers/plate2d/output/spectra.py
@@ -107,6 +107,6 @@
             tmax_plot = times[imax_plot]
 
             print(
-                "plot1d(tmin={0}, tmax={1}, delta_t={2:.2f},".format(
+                "plot1d(tmin={}, tmax={}, delta_t={:.2f},".format(
                     tmin, tmax, delta_t
                 )
@@ -111,6 +111,6 @@
                     tmin, tmax, delta_t
                 )
-                + " coef_compensate={0:.3f})".format(coef_compensate)
+                + f" coef_compensate={coef_compensate:.3f})"
             )
 
             print(
@@ -134,7 +134,7 @@
             ax1.set_title(
                 "1D spectra, solver "
                 + self.output.name_solver
-                + ", nh = {:5d}".format(self.nx)
+                + f", nh = {self.nx:5d}"
             )
             ax1.set_xscale("log")
             ax1.set_yscale("log")
@@ -188,6 +188,6 @@
             tmax_plot = times[imax_plot]
 
             print(
-                "plot2d(tmin={0}, tmax={1}, delta_t={2:.2f},".format(
+                "plot2d(tmin={}, tmax={}, delta_t={:.2f},".format(
                     tmin, tmax, delta_t
                 )
@@ -192,6 +192,6 @@
                     tmin, tmax, delta_t
                 )
-                + " coef_compensate={0:.3f})".format(coef_compensate)
+                + f" coef_compensate={coef_compensate:.3f})"
             )
 
             print(
@@ -215,7 +215,7 @@
             ax1.set_title(
                 "2D spectra, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
+                + f", nh = {self.nx:5d}"
             )
             ax1.set_xscale("log")
             ax1.set_yscale("log")
diff --git a/fluidsim/solvers/plate2d/solver.py b/fluidsim/solvers/plate2d/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL3NvbHZlci5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/plate2d/solver.py
+++ b/fluidsim/solvers/plate2d/solver.py
@@ -52,7 +52,7 @@
 
         """
 
-        super(InfoSolverPlate2D, self)._init_root()
+        super()._init_root()
 
         package = "fluidsim.solvers.plate2d"
         self.module_name = package + ".solver"
@@ -204,7 +204,7 @@
 
     def compute_freq_diss(self):
         """Compute the dissipation frequencies with dissipation only for w."""
-        f_d_w, f_d_hypo_w = super(Simul, self).compute_freq_diss()
+        f_d_w, f_d_hypo_w = super().compute_freq_diss()
         f_d = np.zeros_like(self.state.state_spect, dtype=np.float64)
         f_d_hypo = np.zeros_like(self.state.state_spect, dtype=np.float64)
         f_d[0] = f_d_w
diff --git a/fluidsim/solvers/plate2d/state.py b/fluidsim/solvers/plate2d/state.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL3N0YXRlLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL3N0YXRlLnB5 100644
--- a/fluidsim/solvers/plate2d/state.py
+++ b/fluidsim/solvers/plate2d/state.py
@@ -94,4 +94,4 @@
                 z_fft = np.zeros_like(w_fft)
                 self.init_state_from_wz_fft(w_fft, z_fft)
         else:
-            super(StatePlate2D, self).init_statespect_from(**kwargs)
+            super().init_statespect_from(**kwargs)
diff --git a/fluidsim/solvers/sphere/ns2d/solver.py b/fluidsim/solvers/sphere/ns2d/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zcGhlcmUvbnMyZC9zb2x2ZXIucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zcGhlcmUvbnMyZC9zb2x2ZXIucHk= 100644
--- a/fluidsim/solvers/sphere/ns2d/solver.py
+++ b/fluidsim/solvers/sphere/ns2d/solver.py
@@ -26,7 +26,7 @@
         """Init. `self` by writting the information on the solver.
         """
 
-        super(InfoSolverSphereNS2D, self)._init_root()
+        super()._init_root()
 
         here = "fluidsim.solvers.sphere.ns2d"
 
diff --git a/fluidsim/solvers/sphere/sw1l/solver.py b/fluidsim/solvers/sphere/sw1l/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zcGhlcmUvc3cxbC9zb2x2ZXIucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zcGhlcmUvc3cxbC9zb2x2ZXIucHk= 100644
--- a/fluidsim/solvers/sphere/sw1l/solver.py
+++ b/fluidsim/solvers/sphere/sw1l/solver.py
@@ -39,7 +39,7 @@
         """Init. `self` by writting the information on the solver.
         """
 
-        super(InfoSolverSphereSW1L, self)._init_root()
+        super()._init_root()
 
         here = "fluidsim.solvers.sphere.sw1l"
 
diff --git a/fluidsim/solvers/sw1l/exactlin/modified/solver.py b/fluidsim/solvers/sw1l/exactlin/modified/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2V4YWN0bGluL21vZGlmaWVkL3NvbHZlci5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2V4YWN0bGluL21vZGlmaWVkL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/sw1l/exactlin/modified/solver.py
+++ b/fluidsim/solvers/sw1l/exactlin/modified/solver.py
@@ -18,7 +18,7 @@
     """Information about the solver SW1L."""
 
     def _init_root(self):
-        super(InfoSolverSW1LExactLinModified, self)._init_root()
+        super()._init_root()
 
         sw1l = "fluidsim.solvers.sw1l"
 
diff --git a/fluidsim/solvers/sw1l/exactlin/solver.py b/fluidsim/solvers/sw1l/exactlin/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2V4YWN0bGluL3NvbHZlci5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2V4YWN0bGluL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/sw1l/exactlin/solver.py
+++ b/fluidsim/solvers/sw1l/exactlin/solver.py
@@ -24,7 +24,7 @@
     """Information about the solver SW1L."""
 
     def _init_root(self):
-        super(InfoSolverSW1LExactLin, self)._init_root()
+        super()._init_root()
 
         sw1l = "fluidsim.solvers.sw1l"
 
@@ -47,7 +47,7 @@
         if params.beta != 0:
             raise NotImplementedError("Do not use this solver for beta-plane!")
 
-        super(Simul, self).__init__(params)
+        super().__init__(params)
 
     def tendencies_nonlin(self, state_spect=None, old=None):
         oper = self.oper
diff --git a/fluidsim/solvers/sw1l/exactlin/state.py b/fluidsim/solvers/sw1l/exactlin/state.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2V4YWN0bGluL3N0YXRlLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2V4YWN0bGluL3N0YXRlLnB5 100644
--- a/fluidsim/solvers/sw1l/exactlin/state.py
+++ b/fluidsim/solvers/sw1l/exactlin/state.py
@@ -112,7 +112,7 @@
                 self.it_computed[key2] = it
 
         else:
-            result = super(StateSW1LExactLin, self).compute(
+            result = super().compute(
                 key, SAVE_IN_DICT=SAVE_IN_DICT, RAISE_ERROR=RAISE_ERROR
             )
             SAVE_IN_DICT = False
diff --git a/fluidsim/solvers/sw1l/forcing.py b/fluidsim/solvers/sw1l/forcing.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2ZvcmNpbmcucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2ZvcmNpbmcucHk= 100644
--- a/fluidsim/solvers/sw1l/forcing.py
+++ b/fluidsim/solvers/sw1l/forcing.py
@@ -101,7 +101,7 @@
     _key_forced_default = ("q_fft", "a_fft")
 
     def __init__(self, sim):
-        super(WavesVortices, self).__init__(sim)
+        super().__init__(sim)
         params = sim.params.forcing
         self.forcing_rate = 0.5 * params.forcing_rate
 
diff --git a/fluidsim/solvers/sw1l/init_fields.py b/fluidsim/solvers/sw1l/init_fields.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2luaXRfZmllbGRzLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2luaXRfZmllbGRzLnB5 100644
--- a/fluidsim/solvers/sw1l/init_fields.py
+++ b/fluidsim/solvers/sw1l/init_fields.py
@@ -45,7 +45,7 @@
 
     @classmethod
     def _complete_params_with_default(cls, params):
-        super(InitFieldsWave, cls)._complete_params_with_default(params)
+        super()._complete_params_with_default(params)
         params.init_fields._set_child(cls.tag, attribs={"eta_max": 1.0, "ikx": 2})
 
     def __call__(self):
@@ -87,7 +87,7 @@
 
     @classmethod
     def _complete_params_with_default(cls, params):
-        super(InitFieldsVortexGrid, cls)._complete_params_with_default(params)
+        super()._complete_params_with_default(params)
         params.init_fields._set_child(
             cls.tag, attribs={"omega_max": 1.0, "n_vort": 8, "sd": None}
         )
diff --git a/fluidsim/solvers/sw1l/modified/solver.py b/fluidsim/solvers/sw1l/modified/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zdzFsL21vZGlmaWVkL3NvbHZlci5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zdzFsL21vZGlmaWVkL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/sw1l/modified/solver.py
+++ b/fluidsim/solvers/sw1l/modified/solver.py
@@ -20,7 +20,7 @@
     """Information about the solver SW1L."""
 
     def _init_root(self, **kargs):
-        super(InfoSolverSW1LModified, self)._init_root()
+        super()._init_root()
 
         sw1l = "fluidsim.solvers.sw1l"
 
diff --git a/fluidsim/solvers/sw1l/onlywaves/solver.py b/fluidsim/solvers/sw1l/onlywaves/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zdzFsL29ubHl3YXZlcy9zb2x2ZXIucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zdzFsL29ubHl3YXZlcy9zb2x2ZXIucHk= 100644
--- a/fluidsim/solvers/sw1l/onlywaves/solver.py
+++ b/fluidsim/solvers/sw1l/onlywaves/solver.py
@@ -22,7 +22,7 @@
     """Information about the solver SW1L."""
 
     def _init_root(self):
-        super(InfoSolverSW1LWaves, self)._init_root()
+        super()._init_root()
 
         sw1l = "fluidsim.solvers.sw1l"
 
diff --git a/fluidsim/solvers/sw1l/onlywaves/state.py b/fluidsim/solvers/sw1l/onlywaves/state.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zdzFsL29ubHl3YXZlcy9zdGF0ZS5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zdzFsL29ubHl3YXZlcy9zdGF0ZS5weQ== 100644
--- a/fluidsim/solvers/sw1l/onlywaves/state.py
+++ b/fluidsim/solvers/sw1l/onlywaves/state.py
@@ -115,7 +115,7 @@
             result = self.oper.create_arrayX(value=0)
 
         else:
-            result = super(StateSW1LWaves, self).compute(
+            result = super().compute(
                 key, SAVE_IN_DICT=SAVE_IN_DICT, RAISE_ERROR=RAISE_ERROR
             )
             SAVE_IN_DICT = False
diff --git a/fluidsim/solvers/sw1l/output/_old_spect_energy_budg.py b/fluidsim/solvers/sw1l/output/_old_spect_energy_budg.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9fb2xkX3NwZWN0X2VuZXJneV9idWRnLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9fb2xkX3NwZWN0X2VuZXJneV9idWRnLnB5 100644
--- a/fluidsim/solvers/sw1l/output/_old_spect_energy_budg.py
+++ b/fluidsim/solvers/sw1l/output/_old_spect_energy_budg.py
@@ -17,7 +17,7 @@
         self.c2 = params.c2
         self.f = params.f
 
-        super(SpectralEnergyBudgetSW1LWaves, self).__init__(output)
+        super().__init__(output)
 
     def compute(self):
         """compute spectral energy budget the one time."""
@@ -446,9 +446,9 @@
             tmin_plot = times[imin_plot]
             tmax_plot = times[imax_plot]
 
-            to_print = "plot(tmin={0}, tmax={1}, delta_t={2:.2f})".format(
+            to_print = "plot(tmin={}, tmax={}, delta_t={:.2f})".format(
                 tmin, tmax, delta_t
             )
             print(to_print)
 
             to_print = """plot fluxes 2D
@@ -450,10 +450,10 @@
                 tmin, tmax, delta_t
             )
             print(to_print)
 
             to_print = """plot fluxes 2D
-    tmin = {0:8.6g} ; tmax = {1:8.6g} ; delta_t = {2:8.6g}
-    imin = {3:8d} ; imax = {4:8d} ; delta_i = {5:8d}""".format(
+    tmin = {:8.6g} ; tmax = {:8.6g} ; delta_t = {:8.6g}
+    imin = {:8d} ; imax = {:8d} ; delta_i = {:8d}""".format(
                 tmin_plot, tmax_plot, delta_t, imin_plot, imax_plot, delta_i_plot
             )
             print(to_print)
@@ -470,8 +470,8 @@
             title = (
                 "energy flux, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
-                + ", c = {0:.4g}, f = {1:.4g}".format(np.sqrt(self.c2), self.f)
+                + f", nh = {self.nx:5d}"
+                + ", c = {:.4g}, f = {:.4g}".format(np.sqrt(self.c2), self.f)
             )
             ax1.set_title(title)
             ax1.set_xscale("log")
diff --git a/fluidsim/solvers/sw1l/output/increments.py b/fluidsim/solvers/sw1l/output/increments.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9pbmNyZW1lbnRzLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9pbmNyZW1lbnRzLnB5 100644
--- a/fluidsim/solvers/sw1l/output/increments.py
+++ b/fluidsim/solvers/sw1l/output/increments.py
@@ -10,10 +10,10 @@
     """
 
     def __init__(self, output):
-        super(IncrementsSW1L, self).__init__(output)
+        super().__init__(output)
         params = output.sim.params
         self.c2 = params.c2
         self.f = params.f
 
     def _online_plot_saving(self, dict_results, key="eta"):
         """online plot on pdf"""
@@ -14,9 +14,9 @@
         params = output.sim.params
         self.c2 = params.c2
         self.f = params.f
 
     def _online_plot_saving(self, dict_results, key="eta"):
         """online plot on pdf"""
-        super(IncrementsSW1L, self)._online_plot_saving(dict_results, key=key)
+        super()._online_plot_saving(dict_results, key=key)
 
     def compute(self):
@@ -21,6 +21,6 @@
 
     def compute(self):
-        dict_results = super(IncrementsSW1L, self).compute()
+        dict_results = super().compute()
 
         get_var = self.sim.state.get_var
         ux = get_var("ux")
@@ -91,9 +91,9 @@
         tmin_plot = times[imin_plot]
         tmax_plot = times[imax_plot]
 
-        to_print = "plot(tmin={0}, tmax={1}, delta_t={2:.2f})".format(
+        to_print = "plot(tmin={}, tmax={}, delta_t={:.2f})".format(
             tmin, tmax, delta_t
         )
         print(to_print)
 
         to_print = """plot structure functions
@@ -95,10 +95,10 @@
             tmin, tmax, delta_t
         )
         print(to_print)
 
         to_print = """plot structure functions
-tmin = {0:8.6g} ; tmax = {1:8.6g} ; delta_t = {2:8.6g}
-imin = {3:8d} ; imax = {4:8d} ; delta_i = {5:8d}""".format(
+tmin = {:8.6g} ; tmax = {:8.6g} ; delta_t = {:8.6g}
+imin = {:8d} ; imax = {:8d} ; delta_i = {:8d}""".format(
             tmin_plot, tmax_plot, delta_t, imin_plot, imax_plot, delta_i_plot
         )
         print(to_print)
@@ -126,9 +126,9 @@
         fig, ax1 = self.output.figure_axe(size_axe=size_axe)
         ax1.set_xlabel("$r_x$")
         ax1.set_ylabel(
-            r"$\langle \delta u^{" + "{0}".format(order) + "} \\rangle$"
+            r"$\langle \delta u^{" + f"{order}" + "} \\rangle$"
         )
 
         ax1.set_title(
             "struct. functions, solver "
             + self.output.name_solver
@@ -130,9 +130,9 @@
         )
 
         ax1.set_title(
             "struct. functions, solver "
             + self.output.name_solver
-            + ", nh = {0:5d}".format(self.nx)
+            + f", nh = {self.nx:5d}"
         )
         # +', c = {0:.4g}, f = {1:.4g}'.format(np.sqrt(self.c2), self.f))
 
@@ -270,7 +270,7 @@
             tmin_plot = times[imin_plot]
             tmax_plot = times[imax_plot]
 
-            to_print = "plot(tmin={0}, tmax={1})".format(tmin, tmax)
+            to_print = f"plot(tmin={tmin}, tmax={tmax})"
             print(to_print)
 
             to_print = """plot structure functions
@@ -274,8 +274,8 @@
             print(to_print)
 
             to_print = """plot structure functions
-    tmin = {0:8.6g} ; tmax = {1:8.6g}
-    imin = {2:8d} ; imax = {3:8d}""".format(
+    tmin = {:8.6g} ; tmax = {:8.6g}
+    imin = {:8d} ; imax = {:8d}""".format(
                 tmin_plot, tmax_plot, imin_plot, imax_plot
             )
             print(to_print)
@@ -315,8 +315,8 @@
         title = (
             "struct. functions, solver "
             + self.output.name_solver
-            + ", nh = {0:5d}".format(self.nx)
-            + ", c2 = {0:.4g}, f = {1:.4g}".format(self.c2, self.f)
+            + f", nh = {self.nx:5d}"
+            + f", c2 = {self.c2:.4g}, f = {self.f:.4g}"
         )
         ax1.set_title(title)
         ax1.set_xscale("log")
diff --git a/fluidsim/solvers/sw1l/output/normal_mode.py b/fluidsim/solvers/sw1l/output/normal_mode.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9ub3JtYWxfbW9kZS5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9ub3JtYWxfbW9kZS5weQ== 100644
--- a/fluidsim/solvers/sw1l/output/normal_mode.py
+++ b/fluidsim/solvers/sw1l/output/normal_mode.py
@@ -24,7 +24,7 @@
 from fluiddyn.util import mpi
 
 
-class NormalModeBase(object):
+class NormalModeBase:
     def __init__(self, output):
         self.sim = output.sim
         self.params = output.sim.params
@@ -111,7 +111,7 @@
 
 class NormalModeDecomposition(NormalModeBase):
     def __init__(self, output):
-        super(NormalModeDecomposition, self).__init__(output)
+        super().__init__(output)
 
     @cached_property
     def qmat(self):
@@ -363,7 +363,7 @@
 
             self.bvecrot_fft = bvecrot_fft
 
-        return super(NormalModeDecompositionModified, self).compute()
+        return super().compute()
 
     def normalmodefft_from_keyfft(self, key):
         """Returns the normal mode decomposition for the state_spect key specified."""
@@ -396,6 +396,4 @@
             return key_modes, normal_mode_vec_fft
 
         else:
-            return super(
-                NormalModeDecompositionModified, self
-            ).normalmodefft_from_keyfft(key)
+            return super().normalmodefft_from_keyfft(key)
diff --git a/fluidsim/solvers/sw1l/output/print_stdout.py b/fluidsim/solvers/sw1l/output/print_stdout.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9wcmludF9zdGRvdXQucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9wcmludF9zdGRvdXQucHk= 100644
--- a/fluidsim/solvers/sw1l/output/print_stdout.py
+++ b/fluidsim/solvers/sw1l/output/print_stdout.py
@@ -11,7 +11,7 @@
     the current state of the simulation."""
 
     def _make_str_info(self):
-        to_print = super(PrintStdOutSW1L, self)._make_str_info()
+        to_print = super()._make_str_info()
 
         energyK, energyA = self.output.compute_energiesKA()
         energy = energyK + energyA
@@ -108,7 +108,7 @@
         title = (
             "info stdout, solver "
             + self.output.name_solver
-            + ", nh = {0:5d}".format(self.nx)
+            + f", nh = {self.nx:5d}"
         )
 
         try:
@@ -112,7 +112,7 @@
         )
 
         try:
-            title = title + ", c = {0:.4g}, f = {1:.4g}".format(
+            title = title + ", c = {:.4g}, f = {:.4g}".format(
                 np.sqrt(self.c2), self.f
             )
         except AttributeError:
diff --git a/fluidsim/solvers/sw1l/output/spatial_means.py b/fluidsim/solvers/sw1l/output/spatial_means.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9zcGF0aWFsX21lYW5zLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9zcGF0aWFsX21lYW5zLnB5 100644
--- a/fluidsim/solvers/sw1l/output/spatial_means.py
+++ b/fluidsim/solvers/sw1l/output/spatial_means.py
@@ -322,8 +322,8 @@
         title = (
             "mean energy, solver "
             + self.output.name_solver
-            + ", nh = {0:5d}".format(self.nx)
-            + ", c = {0:.4g}, f = {1:.4g}".format(np.sqrt(self.c2), self.f)
+            + f", nh = {self.nx:5d}"
+            + ", c = {:.4g}, f = {:.4g}".format(np.sqrt(self.c2), self.f)
         )
         ax1.set_title(title)
         norm = self.c2 / 2
@@ -351,8 +351,8 @@
         title = (
             "forcing and dissipation, solver "
             + self.output.name_solver
-            + ", nh = {0:5d}".format(self.nx)
-            + ", c = {0:.4g}, f = {1:.4g}".format(np.sqrt(self.c2), self.f)
+            + f", nh = {self.nx:5d}"
+            + ", c = {:.4g}, f = {:.4g}".format(np.sqrt(self.c2), self.f)
         )
         ax1.set_title(title)
         if "PK_tot" in dict_results:
@@ -426,7 +426,7 @@
         for k in keys:
             E = dict_results[k]
             dE_dt = abs(np.gradient(E, 1.0) / dt)
-            dE_dt_avg = "{0:11.6e}".format(dE_dt.mean())
+            dE_dt_avg = "{:11.6e}".format(dE_dt.mean())
             try:
                 axarr[i].semilogy(t, dE_dt, label=dE_dt_avg)
                 axarr[i].set_ylabel(r"$\partial_t$" + keys[i])
diff --git a/fluidsim/solvers/sw1l/output/spect_energy_budget.py b/fluidsim/solvers/sw1l/output/spect_energy_budget.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9zcGVjdF9lbmVyZ3lfYnVkZ2V0LnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9zcGVjdF9lbmVyZ3lfYnVkZ2V0LnB5 100644
--- a/fluidsim/solvers/sw1l/output/spect_energy_budget.py
+++ b/fluidsim/solvers/sw1l/output/spect_energy_budget.py
@@ -221,9 +221,9 @@
             tmin_plot = times[imin_plot]
             tmax_plot = times[imax_plot]
 
-            to_print = "plot(tmin={0}, tmax={1}, delta_t={2:.2f})".format(
+            to_print = "plot(tmin={}, tmax={}, delta_t={:.2f})".format(
                 tmin, tmax, delta_t
             )
             print(to_print)
 
             to_print = """plot fluxes 2D
@@ -225,10 +225,10 @@
                 tmin, tmax, delta_t
             )
             print(to_print)
 
             to_print = """plot fluxes 2D
-    tmin = {0:8.6g} ; tmax = {1:8.6g} ; delta_t = {2:8.6g}
-    imin = {3:8d} ; imax = {4:8d} ; delta_i = {5:8d}""".format(
+    tmin = {:8.6g} ; tmax = {:8.6g} ; delta_t = {:8.6g}
+    imin = {:8d} ; imax = {:8d} ; delta_i = {:8d}""".format(
                 tmin_plot, tmax_plot, delta_t, imin_plot, imax_plot, delta_i_plot
             )
             print(to_print)
@@ -245,8 +245,8 @@
             title = (
                 "energy flux, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
-                + ", c = {0:.4g}, f = {1:.4g}".format(np.sqrt(self.c2), self.f)
+                + f", nh = {self.nx:5d}"
+                + ", c = {:.4g}, f = {:.4g}".format(np.sqrt(self.c2), self.f)
             )
             ax1.set_title(title)
             ax1.set_xscale("log")
@@ -321,7 +321,7 @@
         else:
             self.norm_mode = norm_mode
 
-        super(SpectralEnergyBudgetSW1L, self).__init__(output)
+        super().__init__(output)
 
     def _transfer_keys_coeff(self):
         c2 = self.params.c2
@@ -551,8 +551,8 @@
             title = (
                 "Spectral Energy Budget, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
-                + ", c = {0:.4g}, f = {1:.4g}".format(np.sqrt(self.c2), self.f)
+                + f", nh = {self.nx:5d}"
+                + ", c = {:.4g}, f = {:.4g}".format(np.sqrt(self.c2), self.f)
             )
             self.axe_a.set_title(title)
             self.axe_a.legend()
@@ -580,7 +580,7 @@
             tmin_plot = times[imin_plot]
             tmax_plot = times[imax_plot]
 
-            to_print = "plot(tmin={0}, tmax={1}, delta_t={2:.2f})".format(
+            to_print = "plot(tmin={}, tmax={}, delta_t={:.2f})".format(
                 tmin, tmax, delta_t
             )
             print(to_print)
@@ -618,8 +618,8 @@
             title = (
                 "Spectral Energy Budget, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
-                + ", c = {0:.4g}, f = {1:.4g}".format(np.sqrt(self.c2), self.f)
+                + f", nh = {self.nx:5d}"
+                + ", c = {:.4g}, f = {:.4g}".format(np.sqrt(self.c2), self.f)
             )
             ax1.set_title(title)
             ax1.set_xscale("log")
@@ -727,7 +727,7 @@
 class SpectralEnergyBudgetSW1LModified(SpectralEnergyBudgetSW1L):
     def __init__(self, output):
         norm_mode = NormalModeDecompositionModified(output)
-        super(SpectralEnergyBudgetSW1LModified, self).__init__(output, norm_mode)
+        super().__init__(output, norm_mode)
 
     def _transfer_keys_coeff(self):
         c2 = self.params.c2
diff --git a/fluidsim/solvers/sw1l/output/spectra.py b/fluidsim/solvers/sw1l/output/spectra.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9zcGVjdHJhLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9zcGVjdHJhLnB5 100644
--- a/fluidsim/solvers/sw1l/output/spectra.py
+++ b/fluidsim/solvers/sw1l/output/spectra.py
@@ -17,6 +17,6 @@
         self.c2 = params.c2
         self.f = params.f
 
-        super(SpectraSW1L, self).__init__(output)
+        super().__init__(output)
 
     def _init_online_plot(self):
@@ -21,7 +21,7 @@
 
     def _init_online_plot(self):
-        super(SpectraSW1L, self)._init_online_plot()
+        super()._init_online_plot()
         if mpi.rank == 0:
             title = (
                 "spectra, solver "
                 + self.output.name_solver
@@ -24,9 +24,9 @@
         if mpi.rank == 0:
             title = (
                 "spectra, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
-                + ", c = {0:.4g}, f = {1:.4g}".format(np.sqrt(self.c2), self.f)
+                + f", nh = {self.nx:5d}"
+                + ", c = {:.4g}, f = {:.4g}".format(np.sqrt(self.c2), self.f)
             )
             self.axe.set_title(title)
 
@@ -187,9 +187,9 @@
             tmin_plot = times[imin_plot]
             tmax_plot = times[imax_plot]
 
-            to_print = "plot1d(tmin={0}, tmax={1}, delta_t={2:.2f},".format(
+            to_print = "plot1d(tmin={}, tmax={}, delta_t={:.2f},".format(
                 tmin, tmax, delta_t
             )
             print(to_print)
 
             to_print = """plot 1D spectra
@@ -191,10 +191,10 @@
                 tmin, tmax, delta_t
             )
             print(to_print)
 
             to_print = """plot 1D spectra
-    tmin = {0:8.6g} ; tmax = {1:8.6g} ; delta_t = {2:8.6g}
-    imin = {3:8d} ; imax = {4:8d} ; delta_i = {5:8d}""".format(
+    tmin = {:8.6g} ; tmax = {:8.6g} ; delta_t = {:8.6g}
+    imin = {:8d} ; imax = {:8d} ; delta_i = {:8d}""".format(
                 tmin_plot, tmax_plot, delta_t, imin_plot, imax_plot, delta_i_plot
             )
             print(to_print)
@@ -207,8 +207,8 @@
             title = (
                 "1D spectra, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
-                + ", c = {0:.4g}, f = {1:.4g}".format(np.sqrt(self.c2), self.f)
+                + f", nh = {self.nx:5d}"
+                + ", c = {:.4g}, f = {:.4g}".format(np.sqrt(self.c2), self.f)
             )
             ax.set_title(title)
             ax.set_xscale("log")
@@ -292,9 +292,9 @@
             tmin_plot = times[imin_plot]
             tmax_plot = times[imax_plot]
 
-            to_print = "plot2d(tmin={0}, tmax={1}, delta_t={2:.2f},".format(
+            to_print = "plot2d(tmin={}, tmax={}, delta_t={:.2f},".format(
                 tmin, tmax, delta_t
             )
             print(to_print)
 
             to_print = """plot 2D spectra
@@ -296,10 +296,10 @@
                 tmin, tmax, delta_t
             )
             print(to_print)
 
             to_print = """plot 2D spectra
-    tmin = {0:8.6g} ; tmax = {1:8.6g} ; delta_t = {2:8.6g}
-    imin = {3:8d} ; imax = {4:8d} ; delta_i = {5:8d}""".format(
+    tmin = {:8.6g} ; tmax = {:8.6g} ; delta_t = {:8.6g}
+    imin = {:8d} ; imax = {:8d} ; delta_i = {:8d}""".format(
                 tmin_plot, tmax_plot, delta_t, imin_plot, imax_plot, delta_i_plot
             )
             print(to_print)
@@ -312,8 +312,8 @@
             title = (
                 "2D spectra, solver "
                 + self.output.name_solver
-                + ", nh = {0:5d}".format(self.nx)
-                + ", c = {0:.4g}, f = {1:.4g}".format(np.sqrt(self.c2), self.f)
+                + f", nh = {self.nx:5d}"
+                + ", c = {:.4g}, f = {:.4g}".format(np.sqrt(self.c2), self.f)
             )
             ax.set_title(title)
             ax.set_xscale("log")
@@ -471,7 +471,7 @@
 class SpectraSW1LNormalMode(SpectraSW1L):
     def __init__(self, output):
         self.norm_mode = NormalModeBase(output)
-        super(SpectraSW1LNormalMode, self).__init__(output)
+        super().__init__(output)
 
     def compute_lin_spectra(self):
         energy_glin_fft, energy_aplin_fft, energy_amlin_fft = (
diff --git a/fluidsim/solvers/sw1l/solver.py b/fluidsim/solvers/sw1l/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zdzFsL3NvbHZlci5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zdzFsL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/sw1l/solver.py
+++ b/fluidsim/solvers/sw1l/solver.py
@@ -78,7 +78,7 @@
         - :class:`fluidsim.solvers.sw1l.forcing.ForcingSW1L`
 
         """
-        super(InfoSolverSW1L, self)._init_root()
+        super()._init_root()
 
         package = "fluidsim.solvers.sw1l"
 
@@ -152,7 +152,7 @@
                 "Equations are non-periodic in this formulation."
             )
 
-        super(Simul, self).__init__(params)
+        super().__init__(params)
 
         if mpi.rank == 0:
             self.output.print_stdout(
@@ -156,7 +156,7 @@
 
         if mpi.rank == 0:
             self.output.print_stdout(
-                "c2 = {0:6.5g} ; f = {1:6.5g} ; kd2 = {2:6.5g}".format(
+                "c2 = {:6.5g} ; f = {:6.5g} ; kd2 = {:6.5g}".format(
                     params.c2, params.f, params.kd2
                 )
             )
diff --git a/fluidsim/solvers/sw1l/state.py b/fluidsim/solvers/sw1l/state.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy9zdzFsL3N0YXRlLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy9zdzFsL3N0YXRlLnB5 100644
--- a/fluidsim/solvers/sw1l/state.py
+++ b/fluidsim/solvers/sw1l/state.py
@@ -178,8 +178,8 @@
                 )
                 init_from_keyfft(value)
             except AttributeError:
-                super(StateSW1L, self).init_statespect_from(**kwargs)
+                super().init_statespect_from(**kwargs)
         elif len(kwargs) == 2:
             if "q_fft" in kwargs and "a_fft" in kwargs:
                 self.init_from_qafft(**kwargs)
         else:
@@ -182,8 +182,8 @@
         elif len(kwargs) == 2:
             if "q_fft" in kwargs and "a_fft" in kwargs:
                 self.init_from_qafft(**kwargs)
         else:
-            super(StateSW1L, self).init_statespect_from(**kwargs)
+            super().init_statespect_from(**kwargs)
 
     def init_from_etafft(self, eta_fft):
         r"""Initialize from :math:`\hat \eta` and set velocities to zero."""
diff --git a/fluidsim/solvers/waves2d/solver.py b/fluidsim/solvers/waves2d/solver.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vc29sdmVycy93YXZlczJkL3NvbHZlci5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vc29sdmVycy93YXZlczJkL3NvbHZlci5weQ== 100644
--- a/fluidsim/solvers/waves2d/solver.py
+++ b/fluidsim/solvers/waves2d/solver.py
@@ -116,7 +116,7 @@
 
     def __init__(self, params):
         self.info_solver = info_solver
-        super(Simul, self).__init__(params)
+        super().__init__(params)
 
     def tendencies_nonlin(self, state_spect=None, old=None):
 
diff --git a/fluidsim/util/console/bench_analysis.py b/fluidsim/util/console/bench_analysis.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vdXRpbC9jb25zb2xlL2JlbmNoX2FuYWx5c2lzLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vdXRpbC9jb25zb2xlL2JlbmNoX2FuYWx5c2lzLnB5 100644
--- a/fluidsim/util/console/bench_analysis.py
+++ b/fluidsim/util/console/bench_analysis.py
@@ -24,7 +24,7 @@
         solver = solver.split(".", 1)[1]
 
     dicts = []
-    for path in glob(path_dir + "/result_bench_{}*.json".format(solver)):
+    for path in glob(path_dir + f"/result_bench_{solver}*.json"):
         with open(path) as f:
             d = json.load(f)
 
@@ -163,8 +163,8 @@
         fig, axes = plt.subplots(1, 2)
         ax0, ax1 = axes.ravel()
 
-    ax0.set_ylabel("speedup ({} scaling)".format(type_plot))
-    ax1.set_ylabel("efficiency % ({} scaling)".format(type_plot))
+    ax0.set_ylabel(f"speedup ({type_plot} scaling)")
+    ax1.set_ylabel(f"efficiency % ({type_plot} scaling)")
 
     def plot_once(ax, x, y, label, linestyle="-k"):
         """Avoid plotting the same label again."""
@@ -236,7 +236,7 @@
                 efficiency.index,
                 efficiency.values,
                 "x-",
-                label="{}, {}".format(name, name_dir),
+                label=f"{name}, {name_dir}",
             )
 
     add_hline(ax1, efficiency, 100)
diff --git a/fluidsim/util/console/profile.py b/fluidsim/util/console/profile.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vdXRpbC9jb25zb2xlL3Byb2ZpbGUucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vdXRpbC9jb25zb2xlL3Byb2ZpbGUucHk= 100644
--- a/fluidsim/util/console/profile.py
+++ b/fluidsim/util/console/profile.py
@@ -325,7 +325,7 @@
 
     long_functions = {}
 
-    key_fft = "fft{}d".format(nb_dim)
+    key_fft = f"fft{nb_dim}d"
     kinds = _kinds + (key_fft,)
     times = {k: 0.0 for k in kinds}
 
@@ -391,7 +391,7 @@
         print(name + ":")
         for k, v in d.items():
             if k == "percentage":
-                print("    " + k + " = {:.2f} %".format(v))
+                print("    " + k + f" = {v:.2f} %")
             else:
                 print("    " + k + " = " + repr(v))
 
diff --git a/fluidsim/util/console/test_profile.py b/fluidsim/util/console/test_profile.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vdXRpbC9jb25zb2xlL3Rlc3RfcHJvZmlsZS5weQ==..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vdXRpbC9jb25zb2xlL3Rlc3RfcHJvZmlsZS5weQ== 100644
--- a/fluidsim/util/console/test_profile.py
+++ b/fluidsim/util/console/test_profile.py
@@ -19,7 +19,7 @@
 if mpi.nb_proc > 1:
     pid = mpi.comm.bcast(pid, root=0)
 
-path_dir = "/tmp/tmp_test_fluidsim_profile_dir_pid{}".format(pid)
+path_dir = f"/tmp/tmp_test_fluidsim_profile_dir_pid{pid}"
 
 
 class TestsProfile(TestCase):
diff --git a/fluidsim/util/console/util.py b/fluidsim/util/console/util.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vdXRpbC9jb25zb2xlL3V0aWwucHk=..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vdXRpbC9jb25zb2xlL3V0aWwucHk= 100644
--- a/fluidsim/util/console/util.py
+++ b/fluidsim/util/console/util.py
@@ -279,7 +279,7 @@
         if n0 is None:
             n0 = 512
     else:
-        raise ConsoleError("dim should not be {}".format(dim))
+        raise ConsoleError(f"dim should not be {dim}")
 
     if n1 is None:
         n1 = n0
@@ -305,7 +305,7 @@
     t_as_str = time_as_str()
     key_solver = sim.info_solver.short_name.lower()
     pid = str(os.getpid())
-    nb_proc = "np={}".format(mpi.nb_proc)
+    nb_proc = f"np={mpi.nb_proc}"
     type_fft = sim.params.oper.type_fft.split(".")[-1].replace("_", "-")
     name_file = (
         "_".join(
diff --git a/fluidsim/util/testing.py b/fluidsim/util/testing.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vdXRpbC90ZXN0aW5nLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vdXRpbC90ZXN0aW5nLnB5 100644
--- a/fluidsim/util/testing.py
+++ b/fluidsim/util/testing.py
@@ -163,7 +163,7 @@
 
         for name, elapsed in result.getTestTimings():
             if elapsed > self.slow_test_threshold:
-                self.write_result("({:.3f}s) {}".format(elapsed, name))
+                self.write_result(f"({elapsed:.3f}s) {name}")
 
         return result
 
@@ -178,7 +178,7 @@
     """ Get the full dotted name of the test method """
 
     mod_name = obj.__class__.__module__.replace("fluidsim.", "")
-    return "%s.%s.%s" % (mod_name, obj.__class__.__name__, obj._testMethodName)
+    return f"{mod_name}.{obj.__class__.__name__}.{obj._testMethodName}"
 
 
 def deactivate_redirect_stdout(tests):
@@ -208,7 +208,7 @@
         msg = "Skipped tests"
         testRunner.write_result("\n", msg, "\n", "-" * len(msg))
         for (case, reason) in result.skipped:
-            testRunner.write_result("S  %s (%s)" % (_mname(case), reason))
+            testRunner.write_result("S  {} ({})".format(_mname(case), reason))
         for (case, reason) in result.expectedFailures:
             testRunner.write_result("X  %s" % _mname(case))
 
diff --git a/fluidsim/util/util.py b/fluidsim/util/util.py
index dcfee750993158507a19b585c48795ecd9fe3f98_Zmx1aWRzaW0vdXRpbC91dGlsLnB5..6360550698ec6db213709bb2f849713d9c584550_Zmx1aWRzaW0vdXRpbC91dGlsLnB5 100644
--- a/fluidsim/util/util.py
+++ b/fluidsim/util/util.py
@@ -408,7 +408,7 @@
 
     print(sim2.params.path_run)
 
-    sim2.output.path_run = str(path_dir) + "/State_phys_{0}x{1}".format(
+    sim2.output.path_run = str(path_dir) + "/State_phys_{}x{}".format(
         sim2.params.oper.nx, sim2.params.oper.ny
     )
     print("Save file in directory\n" + sim2.output.path_run)
diff --git a/scripts/launch/sw1l_force_vortex_grid.py b/scripts/launch/sw1l_force_vortex_grid.py
index dcfee750993158507a19b585c48795ecd9fe3f98_c2NyaXB0cy9sYXVuY2gvc3cxbF9mb3JjZV92b3J0ZXhfZ3JpZC5weQ==..6360550698ec6db213709bb2f849713d9c584550_c2NyaXB0cy9sYXVuY2gvc3cxbF9mb3JjZV92b3J0ZXhfZ3JpZC5weQ== 100755
--- a/scripts/launch/sw1l_force_vortex_grid.py
+++ b/scripts/launch/sw1l_force_vortex_grid.py
@@ -58,5 +58,5 @@
 
 sim = Simul(params)
 sim.output.print_stdout(
-    "Froude number ~ {0:3e}, Viscosity = {1:3e} + {2:3e}".format(
+    "Froude number ~ {:3e}, Viscosity = {:3e} + {:3e}".format(
           (1./params.c2), sim.params.nu_2, sim.params.nu_8))
@@ -62,5 +62,5 @@
           (1./params.c2), sim.params.nu_2, sim.params.nu_8))
-sim.output.print_stdout("Forcing rate = {0}".format(sim.params.forcing.forcing_rate))
+sim.output.print_stdout(f"Forcing rate = {sim.params.forcing.forcing_rate}")
 
 sim.time_stepping.start()
 
diff --git a/scripts/launch/sw1l_no_flow.py b/scripts/launch/sw1l_no_flow.py
index dcfee750993158507a19b585c48795ecd9fe3f98_c2NyaXB0cy9sYXVuY2gvc3cxbF9ub19mbG93LnB5..6360550698ec6db213709bb2f849713d9c584550_c2NyaXB0cy9sYXVuY2gvc3cxbF9ub19mbG93LnB5 100755
--- a/scripts/launch/sw1l_no_flow.py
+++ b/scripts/launch/sw1l_no_flow.py
@@ -74,7 +74,7 @@
 # -------Simulation commences-------------
 sim = Simul(params)
 if mpi.rank == 0:
-    print("Froude number ~ {0:3e}, Viscosity = {1:3e} + {2:3e}".format(
+    print("Froude number ~ {:3e}, Viscosity = {:3e} + {:3e}".format(
         (1./params.c2), sim.params.nu_2, sim.params.nu_8))
 
 sim.time_stepping.start()
diff --git a/scripts/ns2d.strat/coeff_diss.py b/scripts/ns2d.strat/coeff_diss.py
index dcfee750993158507a19b585c48795ecd9fe3f98_c2NyaXB0cy9uczJkLnN0cmF0L2NvZWZmX2Rpc3MucHk=..6360550698ec6db213709bb2f849713d9c584550_c2NyaXB0cy9uczJkLnN0cmF0L2NvZWZmX2Rpc3MucHk= 100644
--- a/scripts/ns2d.strat/coeff_diss.py
+++ b/scripts/ns2d.strat/coeff_diss.py
@@ -225,7 +225,7 @@
     if time_total > 1000:
         print(
             "The stationarity has not " + \
-            "reached after {} simulations.".format(it))
+            f"reached after {it} simulations.")
         should_I_stop = "non_stationarity"
 
     if ratio_x > threshold_ratio and ratio_y > threshold_ratio:
@@ -245,8 +245,8 @@
             should_I_stop = False
 
         else:
-            print("Checking stationarity... with nu8 = {}".format(nu_8_old))
+            print(f"Checking stationarity... with nu8 = {nu_8_old}")
             E, t, P_tot = compute_energy_spatial(sim)
             ratio = np.mean(np.diff(E[2:]) / np.diff(t[2:]))
             if (ratio / injection_energy_0) < 0.5 and \
                abs(nu_8_old - params.nu_8) / params.nu_8 < 0.05:
@@ -249,8 +249,8 @@
             E, t, P_tot = compute_energy_spatial(sim)
             ratio = np.mean(np.diff(E[2:]) / np.diff(t[2:]))
             if (ratio / injection_energy_0) < 0.5 and \
                abs(nu_8_old - params.nu_8) / params.nu_8 < 0.05:
-                print("Stationarity is reached.\n nu_8 = {}".format(params.nu_8))
+                print(f"Stationarity is reached.\n nu_8 = {params.nu_8}")
                 factor = 1.
                 should_I_stop = True
             else:
@@ -269,7 +269,7 @@
 
 # Create directory in path
 path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim/Coef_Diss"
-name_directory = "Coef_Diss_gamma{}".format(gamma_not)
+name_directory = f"Coef_Diss_gamma{gamma_not}"
 path = os.path.join(path_root, name_directory)
 
 if mpi.rank == 0 and not os.path.exists(path):
diff --git a/scripts/ns2d.strat/compute_anisotropy.py b/scripts/ns2d.strat/compute_anisotropy.py
index dcfee750993158507a19b585c48795ecd9fe3f98_c2NyaXB0cy9uczJkLnN0cmF0L2NvbXB1dGVfYW5pc290cm9weS5weQ==..6360550698ec6db213709bb2f849713d9c584550_c2NyaXB0cy9uczJkLnN0cmF0L2NvbXB1dGVfYW5pc290cm9weS5weQ== 100644
--- a/scripts/ns2d.strat/compute_anisotropy.py
+++ b/scripts/ns2d.strat/compute_anisotropy.py
@@ -60,7 +60,7 @@
     else:
         gamma_out = float(gamma_str)
 
-    print("Compute anisotropy nx = {} and gamma {}..".format(res_out, gamma_out))
+    print(f"Compute anisotropy nx = {res_out} and gamma {gamma_out}..")
 
     # Load data energy spectra file.
     with h5py.File(path_simulation + "/spectra1D.h5", "r") as f:
@@ -93,4 +93,4 @@
                     "sim1920_no_shear_modes/NS2D.strat_1920x480_S2pix1.571_F07_gamma05_2018-08-14_10-01-22")
 
     anisotropy = compute_anisotropy(path_simulation)
-    print("anisotropy = {}".format(anisotropy))
+    print(f"anisotropy = {anisotropy}")
diff --git a/scripts/ns2d.strat/compute_flow_features.py b/scripts/ns2d.strat/compute_flow_features.py
index dcfee750993158507a19b585c48795ecd9fe3f98_c2NyaXB0cy9uczJkLnN0cmF0L2NvbXB1dGVfZmxvd19mZWF0dXJlcy5weQ==..6360550698ec6db213709bb2f849713d9c584550_c2NyaXB0cy9uczJkLnN0cmF0L2NvbXB1dGVfZmxvd19mZWF0dXJlcy5weQ== 100644
--- a/scripts/ns2d.strat/compute_flow_features.py
+++ b/scripts/ns2d.strat/compute_flow_features.py
@@ -21,5 +21,5 @@
 path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim"
 
 if nx == 1920:
-    directory = "sim{}_no_shear_modes".format(nx)
+    directory = f"sim{nx}_no_shear_modes"
 elif nx == 3840 or nx == 7680:
@@ -25,5 +25,5 @@
 elif nx == 3840 or nx == 7680:
-    directory = "sim{}_modif_res_no_shear_modes".format(nx)
+    directory = f"sim{nx}_modif_res_no_shear_modes"
 else:
     raise ValueError(".")
 
diff --git a/scripts/ns2d.strat/compute_ratio_dissipation.py b/scripts/ns2d.strat/compute_ratio_dissipation.py
index dcfee750993158507a19b585c48795ecd9fe3f98_c2NyaXB0cy9uczJkLnN0cmF0L2NvbXB1dGVfcmF0aW9fZGlzc2lwYXRpb24ucHk=..6360550698ec6db213709bb2f849713d9c584550_c2NyaXB0cy9uczJkLnN0cmF0L2NvbXB1dGVfcmF0aW9fZGlzc2lwYXRpb24ucHk= 100644
--- a/scripts/ns2d.strat/compute_ratio_dissipation.py
+++ b/scripts/ns2d.strat/compute_ratio_dissipation.py
@@ -21,7 +21,7 @@
     else:
         gamma_out = float(gamma_str)
 
-    print("Compute dissipation nx = {} and gamma {}..".format(res_out, gamma_out))
+    print(f"Compute dissipation nx = {res_out} and gamma {gamma_out}..")
 
     # Load object parameters
     params = load_params_simul(path_simulation)
diff --git a/scripts/ns2d.strat/compute_spectrum_kykx.py b/scripts/ns2d.strat/compute_spectrum_kykx.py
index dcfee750993158507a19b585c48795ecd9fe3f98_c2NyaXB0cy9uczJkLnN0cmF0L2NvbXB1dGVfc3BlY3RydW1fa3lreC5weQ==..6360550698ec6db213709bb2f849713d9c584550_c2NyaXB0cy9uczJkLnN0cmF0L2NvbXB1dGVfc3BlY3RydW1fa3lreC5weQ== 100644
--- a/scripts/ns2d.strat/compute_spectrum_kykx.py
+++ b/scripts/ns2d.strat/compute_spectrum_kykx.py
@@ -123,6 +123,6 @@
 
 if key == "ap":
     data = np.log10(abs(ap_fft_arr)**2)
-    text = "$\hat{a}_+$"
+    text = r"$\hat{a}_+$"
 elif key == "am":
     data = np.log10(abs(am_fft_arr)**2)
@@ -127,6 +127,6 @@
 elif key == "am":
     data = np.log10(abs(am_fft_arr)**2)
-    text = "$\hat{a}_-$"
+    text = r"$\hat{a}_-$"
 else:
     raise ValueError("Not implemented")
 
@@ -266,6 +266,6 @@
 
 if SAVE:
     path_save = "/home/users/calpelin7m/Phd/docs/Manuscript/figures"
-    name = "spectrakykx_forced{}_key_{}_t_{}.png".format(index_simulation, key, tmin)
+    name = f"spectrakykx_forced{index_simulation}_key_{key}_t_{tmin}.png"
     fig1.savefig(os.path.join(path_save, name), bbox_inches="tight")
 plt.show()
diff --git a/scripts/ns2d.strat/compute_vertical_length.py b/scripts/ns2d.strat/compute_vertical_length.py
index dcfee750993158507a19b585c48795ecd9fe3f98_c2NyaXB0cy9uczJkLnN0cmF0L2NvbXB1dGVfdmVydGljYWxfbGVuZ3RoLnB5..6360550698ec6db213709bb2f849713d9c584550_c2NyaXB0cy9uczJkLnN0cmF0L2NvbXB1dGVfdmVydGljYWxfbGVuZ3RoLnB5 100644
--- a/scripts/ns2d.strat/compute_vertical_length.py
+++ b/scripts/ns2d.strat/compute_vertical_length.py
@@ -22,7 +22,7 @@
     else:
         gamma_out = float(gamma_str)
 
-    print("Compute dissipation nx = {} and gamma {}..".format(res_out, gamma_out))
+    print(f"Compute dissipation nx = {res_out} and gamma {gamma_out}..")
 
     # Load object parameters
     params = load_params_simul(path_simulation)
diff --git a/scripts/ns2d.strat/make_table_parameters.py b/scripts/ns2d.strat/make_table_parameters.py
index dcfee750993158507a19b585c48795ecd9fe3f98_c2NyaXB0cy9uczJkLnN0cmF0L21ha2VfdGFibGVfcGFyYW1ldGVycy5weQ==..6360550698ec6db213709bb2f849713d9c584550_c2NyaXB0cy9uczJkLnN0cmF0L21ha2VfdGFibGVfcGFyYW1ldGVycy5weQ== 100644
--- a/scripts/ns2d.strat/make_table_parameters.py
+++ b/scripts/ns2d.strat/make_table_parameters.py
@@ -31,9 +31,9 @@
 
 if MAKE_TABLE:
     path_table = ("/home/users/calpelin7m/" +
-                  "Phd/docs/Manuscript/buoyancy_reynolds_table_n{}.tex".format(nx))
+                  f"Phd/docs/Manuscript/buoyancy_reynolds_table_n{nx}.tex")
 
     to_print = ("\\begin{table}[h]\n"
                 "\\centering \n"
                 "\\begin{tabular}{cccc} \n"
                 "\\toprule[1.5pt] \n" + \
@@ -35,10 +35,10 @@
 
     to_print = ("\\begin{table}[h]\n"
                 "\\centering \n"
                 "\\begin{tabular}{cccc} \n"
                 "\\toprule[1.5pt] \n" + \
-                "\\bm{$\gamma$} & \\bm{$F_h$} & \\bm{$Re_8$} & \\bm{$\mathcal{R}$} \\\\ \n"
-                "\\midrule\ \n")
+                "\\bm{$\\gamma$} & \\bm{$F_h$} & \\bm{$Re_8$} & \\bm{$\\mathcal{R}$} \\\\ \n"
+                "\\midrule\\ \n")
 
 for path in path_simulations:
 
diff --git a/scripts/ns2d.strat/make_video_spectrum_kykx.py b/scripts/ns2d.strat/make_video_spectrum_kykx.py
index dcfee750993158507a19b585c48795ecd9fe3f98_c2NyaXB0cy9uczJkLnN0cmF0L21ha2VfdmlkZW9fc3BlY3RydW1fa3lreC5weQ==..6360550698ec6db213709bb2f849713d9c584550_c2NyaXB0cy9uczJkLnN0cmF0L21ha2VfdmlkZW9fc3BlY3RydW1fa3lreC5weQ== 100644
--- a/scripts/ns2d.strat/make_video_spectrum_kykx.py
+++ b/scripts/ns2d.strat/make_video_spectrum_kykx.py
@@ -168,7 +168,7 @@
 
 ax.plot(kx, 1e-0 * sim.params.N * (kx / eps)**(1/3), color="white")
 ax.plot(kx, -1e-0 * sim.params.N * (kx / eps)**(1/3), color="white")
-ax.text(kx[ikx_text], kz[ikz_text], "\hat{a}_+", color="white", fontsize=15)
+ax.text(kx[ikx_text], kz[ikz_text], r"\hat{a}_+", color="white", fontsize=15)
 
 ax2.plot(kx, 1e-0 * sim.params.N * (kx / eps)**(1/3), color="white")
 ax2.plot(kx, -1e-0 * sim.params.N * (kx / eps)**(1/3), color="white")
@@ -172,7 +172,7 @@
 
 ax2.plot(kx, 1e-0 * sim.params.N * (kx / eps)**(1/3), color="white")
 ax2.plot(kx, -1e-0 * sim.params.N * (kx / eps)**(1/3), color="white")
-ax2.text(kx[ikx_text], kz[ikz_text], "\hat{a}_-", color="white", fontsize=15)
+ax2.text(kx[ikx_text], kz[ikz_text], r"\hat{a}_-", color="white", fontsize=15)
 
 if scale == "linear":
   _im = ax.pcolormesh(KX, KZ, data_ap_plot_modified, vmin=vmin, vmax=vmax)
diff --git a/scripts/ns2d.strat/submit_legi.py b/scripts/ns2d.strat/submit_legi.py
index dcfee750993158507a19b585c48795ecd9fe3f98_c2NyaXB0cy9uczJkLnN0cmF0L3N1Ym1pdF9sZWdpLnB5..6360550698ec6db213709bb2f849713d9c584550_c2NyaXB0cy9uczJkLnN0cmF0L3N1Ym1pdF9sZWdpLnB5 100644
--- a/scripts/ns2d.strat/submit_legi.py
+++ b/scripts/ns2d.strat/submit_legi.py
@@ -12,8 +12,8 @@
 cluster.commands_setting_env.append(
     'export FLUIDSIM_PATH="/fsnet/project/meige/2015/15DELDUCA/DataSim"')
 
-name_run_root = "find_coeff_nu8_gamma={}".format(args.gamma)
+name_run_root = f"find_coeff_nu8_gamma={args.gamma}"
 
 walltime = '24:00:00'
 nb_proc = 8
 
@@ -16,8 +16,8 @@
 
 walltime = '24:00:00'
 nb_proc = 8
 
-command_to_submit = "python coeff_diss.py {}".format(args.gamma)
+command_to_submit = f"python coeff_diss.py {args.gamma}"
 
 cluster.submit_command(
     command_to_submit,
diff --git a/scripts/ns2d.strat/tests/check_freq_spectra_with_linear_mode.py b/scripts/ns2d.strat/tests/check_freq_spectra_with_linear_mode.py
index dcfee750993158507a19b585c48795ecd9fe3f98_c2NyaXB0cy9uczJkLnN0cmF0L3Rlc3RzL2NoZWNrX2ZyZXFfc3BlY3RyYV93aXRoX2xpbmVhcl9tb2RlLnB5..6360550698ec6db213709bb2f849713d9c584550_c2NyaXB0cy9uczJkLnN0cmF0L3Rlc3RzL2NoZWNrX2ZyZXFfc3BlY3RyYV93aXRoX2xpbmVhcl9tb2RlLnB5 100644
--- a/scripts/ns2d.strat/tests/check_freq_spectra_with_linear_mode.py
+++ b/scripts/ns2d.strat/tests/check_freq_spectra_with_linear_mode.py
@@ -149,7 +149,7 @@
         #     ax.loglog(omegas/omega_n, freq_spectrum[0, :, 1, i])
 
         ax.semilogy(omegas/omega_n, freq_spectrum[0, :, 1, 4])
-        ax.axvline(x=omega_n / omega_n, label="$\omega_{th}$", c="k")
+        ax.axvline(x=omega_n / omega_n, label=r"$\omega_{th}$", c="k")
 
         # Set text
         ax.text(5e-2, 1e3,
diff --git a/scripts/plot_results/plot_interm.py b/scripts/plot_results/plot_interm.py
index dcfee750993158507a19b585c48795ecd9fe3f98_c2NyaXB0cy9wbG90X3Jlc3VsdHMvcGxvdF9pbnRlcm0ucHk=..6360550698ec6db213709bb2f849713d9c584550_c2NyaXB0cy9wbG90X3Jlc3VsdHMvcGxvdF9pbnRlcm0ucHk= 100644
--- a/scripts/plot_results/plot_interm.py
+++ b/scripts/plot_results/plot_interm.py
@@ -106,7 +106,7 @@
     pol = np.polyfit(np.log(rxs[condr]), np.log(M_order[condr]), 1)
     expo = pol[0]
 
-    print('order = {0:.2f} ; expo = {1:.2f}'.format(order, expo))
+    print(f'order = {order:.2f} ; expo = {expo:.2f}')
     M_lin = np.exp((pol[1] + np.log(rxs[condr])*pol[0]))
 
 
@@ -114,10 +114,10 @@
         fig, ax1 = sim.output.figure_axe()
         title = (
             'struct. function, solver '+sim.output.name_solver+
-            ', nh = {0:5d}'.format(sim.param.nx)+
-            ', c = {0:.4g}, f = {1:.4g}'.format(np.sqrt(sim.param.c2),
+            f', nh = {sim.param.nx:5d}'+
+            ', c = {:.4g}, f = {:.4g}'.format(np.sqrt(sim.param.c2),
                                                 sim.param.f)
             )
 
         ax1.set_xlabel('$r$')
         ax1.set_ylabel(
@@ -119,9 +119,9 @@
                                                 sim.param.f)
             )
 
         ax1.set_xlabel('$r$')
         ax1.set_ylabel(
-            '$\langle |\delta v|^{'+'{0:.2f}'.format(order)+'}\\rangle/r$')
+            r'$\langle |\delta v|^{'+f'{order:.2f}'+'}\\rangle/r$')
 
         ax1.set_title(title)
         ax1.hold(True)
@@ -146,8 +146,8 @@
             rxs[cond],temp[cond]/temp[0]*M_order[0]/norm[0],
             'k:', linewidth=1)
 
-        plt.text(1, 6*10**0,'order = {0:.2f}'.format(order),fontsize=16)
-        plt.text(0.11, 4*10**-1,'expo = {0:.3f}'.format(expo),fontsize=16)
+        plt.text(1, 6*10**0,f'order = {order:.2f}',fontsize=16)
+        plt.text(0.11, 4*10**-1,f'expo = {expo:.3f}',fontsize=16)
 
 
 
@@ -166,8 +166,8 @@
     if PLOT_PDF:
         fig, ax1 = sim.output.figure_axe()
         title = ('pdf increments, solver '+sim.output.name_solver+
-', nh = {0:5d}'.format(sim.param.nx)+
-', c2 = {0:.4g}, f = {1:.4g}'.format(sim.param.c2, sim.param.f)
+f', nh = {sim.param.nx:5d}'+
+f', c2 = {sim.param.c2:.4g}, f = {sim.param.f:.4g}'
                  )
 
         ax1.set_title(title)
@@ -176,7 +176,7 @@
         ax1.set_yscale('linear')
 
         ax1.set_xlabel(key_var)
-        ax1.set_ylabel('PDF x $\delta v^'+repr(order)+'$')
+        ax1.set_ylabel(r'PDF x $\delta v^'+repr(order)+'$')
 
         colors = ['k', 'y', 'r', 'b', 'g', 'm', 'c']
 
@@ -209,8 +209,8 @@
 
 fig, ax1 = sim.output.figure_axe()
 title = ('intermittency, solver '+sim.output.name_solver+
-', nh = {0:5d}'.format(sim.param.nx)+
-', c = {0:.4g}, f = {1:.4g}'.format(np.sqrt(sim.param.c2), sim.param.f)
+f', nh = {sim.param.nx:5d}'+
+', c = {:.4g}, f = {:.4g}'.format(np.sqrt(sim.param.c2), sim.param.f)
 )
 
 ax1.set_title(title)
@@ -219,7 +219,7 @@
 ax1.set_yscale('linear')
 
 ax1.set_xlabel('order q')
-ax1.set_ylabel('$\zeta_q$')
+ax1.set_ylabel(r'$\zeta_q$')
 
 ax1.plot(orders, expos)
 
diff --git a/scripts/plot_results/plot_spatial_means.py b/scripts/plot_results/plot_spatial_means.py
index dcfee750993158507a19b585c48795ecd9fe3f98_c2NyaXB0cy9wbG90X3Jlc3VsdHMvcGxvdF9zcGF0aWFsX21lYW5zLnB5..6360550698ec6db213709bb2f849713d9c584550_c2NyaXB0cy9wbG90X3Jlc3VsdHMvcGxvdF9zcGF0aWFsX21lYW5zLnB5 100644
--- a/scripts/plot_results/plot_spatial_means.py
+++ b/scripts/plot_results/plot_spatial_means.py
@@ -65,8 +65,8 @@
 ax1.set_xlabel('$t$')
 ax1.set_ylabel('$2E(t)/c^2$')
 title = ('mean energy, solver ' + sim.output.name_solver +
-         ', nh = {0:5d}'.format(nx) +
-         ', c2 = {0:.4g}, f = {1:.4g}'.format(c2, f))
+         f', nh = {nx:5d}' +
+         f', c2 = {c2:.4g}, f = {f:.4g}')
 
 ax1.set_title(title)
 ax1.hold(True)
@@ -92,5 +92,5 @@
 size_axe[1] = z_bottom_axe
 fig, ax1 = sim.output.figure_axe(size_axe=size_axe)
 ax1.set_xlabel('t')
-ax1.set_ylabel('$P_E(t)$, $\epsilon(t)$')
+ax1.set_ylabel(r'$P_E(t)$, $\epsilon(t)$')
 title = ('forcing and dissipation, solver '+sim.output.name_solver+
@@ -96,6 +96,6 @@
 title = ('forcing and dissipation, solver '+sim.output.name_solver+
-         ', nh = {0:5d}'.format(nx)+
-         ', c2 = {0:.4g}, f = {1:.4g}'.format(c2, f))
+         f', nh = {nx:5d}'+
+         f', c2 = {c2:.4g}, f = {f:.4g}')
 ax1.set_title(title)
 ax1.hold(True)
 (l_P_tot,) = ax1.plot(t, P_tot, 'c', linewidth=2 )
@@ -103,9 +103,9 @@
 (l_epsE_hypo,) = ax1.plot(t, epsE_hypo, 'g', linewidth=2 )
 (l_epsE_tot,) = ax1.plot(t, epsE_tot, 'k', linewidth=2 )
 l_P_tot.set_label('$P_{tot}$')
-l_epsE.set_label('$\epsilon$')
-l_epsE_hypo.set_label('$\epsilon_{hypo}$')
-l_epsE_tot.set_label('$\epsilon_{tot}$')
+l_epsE.set_label(r'$\epsilon$')
+l_epsE_hypo.set_label(r'$\epsilon_{hypo}$')
+l_epsE_tot.set_label(r'$\epsilon_{tot}$')
 
 ax1.legend(loc=2)
 
@@ -113,7 +113,7 @@
 size_axe[1] = z_bottom_axe
 ax2 = fig.add_axes(size_axe)
 ax2.set_xlabel('t')
-ax2.set_ylabel('$\epsilon$ Charney PE(t)')
+ax2.set_ylabel(r'$\epsilon$ Charney PE(t)')
 title = ('dissipation Charney PE')
 ax2.set_title(title)
 ax2.hold(True)
diff --git a/scripts/plot_results/plot_time_means.py b/scripts/plot_results/plot_time_means.py
index dcfee750993158507a19b585c48795ecd9fe3f98_c2NyaXB0cy9wbG90X3Jlc3VsdHMvcGxvdF90aW1lX21lYW5zLnB5..6360550698ec6db213709bb2f849713d9c584550_c2NyaXB0cy9wbG90X3Jlc3VsdHMvcGxvdF90aW1lX21lYW5zLnB5 100644
--- a/scripts/plot_results/plot_time_means.py
+++ b/scripts/plot_results/plot_time_means.py
@@ -88,7 +88,7 @@
 fig, ax1 = sim.output.figure_axe()
 
 ax1.set_xlabel('$2E/c^2$')
-ax1.set_ylabel('$\epsilon / P$')
+ax1.set_ylabel(r'$\epsilon / P$')
 ax1.hold(True)
 
 coef_norm = Froude_numbers**-1.54
diff --git a/setup.py b/setup.py
index dcfee750993158507a19b585c48795ecd9fe3f98_c2V0dXAucHk=..6360550698ec6db213709bb2f849713d9c584550_c2V0dXAucHk= 100644
--- a/setup.py
+++ b/setup.py
@@ -144,7 +144,7 @@
     extensions = init_pythran_extensions(
         "fluidsim",
         include_dirs=np.get_include(),
-        compile_args=("-O3", "-march={}".format(compile_arch), "-DUSE_XSIMD"),
+        compile_args=("-O3", f"-march={compile_arch}", "-DUSE_XSIMD"),
     )
     return extensions
 
@@ -158,7 +158,7 @@
 
     transonize()
 
-    logger.info("Importing mpi4py: {}".format(MPI4PY))
+    logger.info(f"Importing mpi4py: {MPI4PY}")
 
     define_macros = []
     if has_cython and os.getenv("TOXENV") is not None: