diff --git a/docs/source/conf.py b/docs/source/conf.py
index 5765a588f..2016a27ea 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -37,7 +37,6 @@
 # Add the directory containing your Python script to the Python path
 sys.path.insert(0, os.path.abspath("."))
 
-import map
 
 m = map.generate_map()
 current_dir = os.path.dirname(__file__)
diff --git a/examples/advection_diffusion_multi_material.py b/examples/advection_diffusion_multi_material.py
index a1869335a..96850d2b8 100644
--- a/examples/advection_diffusion_multi_material.py
+++ b/examples/advection_diffusion_multi_material.py
@@ -27,7 +27,6 @@ def half(x):
     )
 
     # Split domain in half and set an interface tag of 5
-    gdim = mesh.geometry.dim
     tdim = mesh.topology.dim
     fdim = tdim - 1
     top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary)
@@ -80,10 +79,10 @@ def velocity_func(x):
         values[1] = 0  # Second component remains zero
         return values
 
-    mesh2, mt2, ct2 = generate_mesh(n=10)
-    submesh, cell_map, v_map = dolfinx.mesh.create_submesh(mesh2, ct2.dim, ct2.find(4))[
-        0:3
-    ]
+    mesh2, _mt2, ct2 = generate_mesh(n=10)
+    submesh, _cell_map, _v_map = dolfinx.mesh.create_submesh(
+        mesh2, ct2.dim, ct2.find(4)
+    )[0:3]
     v_cg = basix.ufl.element(
         "Lagrange", submesh.topology.cell_name(), 2, shape=(submesh.geometry.dim,)
     )
diff --git a/examples/convergence_rates.py b/examples/convergence_rates.py
index 6286c8e5e..a6b24d702 100644
--- a/examples/convergence_rates.py
+++ b/examples/convergence_rates.py
@@ -112,7 +112,9 @@ def mms_source(x, t):
     my_problem.settings = F.Settings(
         atol=1e-8,
         rtol=1e-10,
-        final_time=1,  # final time shouldn't be too long so that a potential error at the initial timestep is not negligible
+        # final time shouldn't be too long so that a potential error at the
+        # initial timestep is not negligible
+        final_time=1,
     )
 
     # Forward euler isn't great so dt should be small
@@ -127,7 +129,8 @@ def mms_source(x, t):
     my_problem.run()
 
     computed_solution = my_problem.u
-    # we use the exact final time of the simulation which may differ from the one specified in the settings
+    # we use the exact final time of the simulation which may differ from
+    # the one specified in the settings
     final_time_sim = my_problem.t.value
 
     def exact_solution_end(x):
diff --git a/examples/multi_material_1d.py b/examples/multi_material_1d.py
index 2fb2a815c..4810f27bb 100644
--- a/examples/multi_material_1d.py
+++ b/examples/multi_material_1d.py
@@ -43,7 +43,8 @@
 left_surface = F.SurfaceSubdomain1D(id=1, x=vertices[0])
 right_surface = F.SurfaceSubdomain1D(id=2, x=vertices[-1])
 
-# the ids here are arbitrary in 1D, you can put anything as long as it's not the same as the surfaces
+# the ids here are arbitrary in 1D, you can put anything as long as it's
+# not the same as the surfaces
 my_model.interfaces = [
     F.Interface(6, (left_domain, middle_domain)),
     F.Interface(7, (middle_domain, right_domain)),
diff --git a/examples/multi_material_2d.py b/examples/multi_material_2d.py
index 7a64cbe8c..201961511 100644
--- a/examples/multi_material_2d.py
+++ b/examples/multi_material_2d.py
@@ -24,7 +24,6 @@ def half(x):
     )
 
     # Split domain in half and set an interface tag of 5
-    gdim = mesh.geometry.dim
     tdim = mesh.topology.dim
     fdim = tdim - 1
     top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary)
@@ -125,11 +124,17 @@ def half(x):
     from dolfinx.mesh import EntityMap  # noqa: F401
 
     legacy_entity_map = False
-    entity_map_wrapper = lambda e_map: list(e_map.values())
+
+    def entity_map_wrapper(e_map):
+        return list(e_map.values())
+
     entity_maps = [sd.cell_map for sd in my_model.volume_subdomains]
 except ImportError:
     legacy_entity_map = True
-    entity_map_wrapper = lambda e_map: e_map
+
+    def entity_map_wrapper(e_map):
+        return e_map
+
     entity_maps = {
         sd.submesh: sd.parent_to_submesh for sd in my_model.volume_subdomains
     }
diff --git a/examples/multi_material_2d_bis.py b/examples/multi_material_2d_bis.py
index cc6ad33da..933747c13 100644
--- a/examples/multi_material_2d_bis.py
+++ b/examples/multi_material_2d_bis.py
@@ -24,7 +24,6 @@ def half(x):
     )
 
     # Split domain in half and set an interface tag of 5
-    gdim = mesh.geometry.dim
     tdim = mesh.topology.dim
     fdim = tdim - 1
     top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary)
@@ -168,10 +167,14 @@ def half(x):
     from dolfinx.mesh import EntityMap  # noqa: F401
 
     legacy_entity_map = False
-    entity_map_wrapper = lambda e_map: list(e_map.values())
+
+    def entity_map_wrapper(e_map):
+        return list(e_map.values())
 except ImportError:
     legacy_entity_map = True
-    entity_map_wrapper = lambda e_map: e_map
+
+    def entity_map_wrapper(e_map):
+        return e_map
 
 
 form = dolfinx.fem.form(
diff --git a/examples/multi_material_transient.py b/examples/multi_material_transient.py
index 8495336fa..54b30cfc0 100644
--- a/examples/multi_material_transient.py
+++ b/examples/multi_material_transient.py
@@ -36,8 +36,8 @@
 left_surface = F.SurfaceSubdomain1D(id=1, x=vertices[0])
 right_surface = F.SurfaceSubdomain1D(id=2, x=vertices[-1])
 
-# the ids here are arbitrary in 1D, you can put anything as long as it's not the same as the surfaces
-# TODO remove mesh and meshtags from these arguments
+# the ids here are arbitrary in 1D, you can put anything as long as it's not the same
+# as the surfaces
 my_model.interfaces = [
     F.Interface(6, (left_domain, middle_domain)),
     F.Interface(7, (middle_domain, right_domain)),
diff --git a/examples/surface_reaction_example.py b/examples/surface_reaction_example.py
index 2f022e0ba..14d41ea80 100644
--- a/examples/surface_reaction_example.py
+++ b/examples/surface_reaction_example.py
@@ -1,4 +1,5 @@
 import dolfinx.fem as fem
+import matplotlib.pyplot as plt
 import numpy as np
 import ufl
 
@@ -98,7 +99,6 @@ def compute(self, ds):
 my_model.initialise()
 my_model.run()
 
-import matplotlib.pyplot as plt
 
 plt.stackplot(
     H_flux_left.t,
diff --git a/src/festim/__init__.py b/src/festim/__init__.py
index 823d519d9..d0a824587 100644
--- a/src/festim/__init__.py
+++ b/src/festim/__init__.py
@@ -1,3 +1,4 @@
+# ruff: noqa: F401, E402
 from importlib import metadata
 
 try:
diff --git a/src/festim/advection.py b/src/festim/advection.py
index 2ecf944f3..10081f823 100644
--- a/src/festim/advection.py
+++ b/src/festim/advection.py
@@ -10,8 +10,7 @@
 
 
 class AdvectionTerm:
-    """
-    Advection term class
+    """Advection term class.
 
     args:
         velocity: the velocity field or function
@@ -22,7 +21,6 @@ class AdvectionTerm:
         velocity: the velocity field or function
         subdomain: the volume subdomain where the velocity is to be applied
         species: the species to which the velocity field is acting on
-
     """
 
     velocity: fem.Function
@@ -94,9 +92,8 @@ def species(self, value):
 
 
 class VelocityField(Value):
-    """
-    A class to handle input values of velocity fields from users and convert them to a
-    relevent fenics object
+    """A class to handle input values of velocity fields from users and convert them to
+    a relevent fenics object.
 
     Args:
         input_value: The value of the user input
@@ -124,7 +121,7 @@ def convert_input_value(
         function_space: fem.FunctionSpace,
         t: fem.Constant | None = None,
     ):
-        """Converts a user given value to a relevent fenics object
+        """Converts a user given value to a relevent fenics object.
 
         Args:
             function_space: the function space of the fenics object
@@ -165,7 +162,7 @@ def convert_input_value(
         nmm_interpolate(self.fenics_object, vel)
 
     def update(self, t: fem.Constant):
-        """Updates the velocity field
+        """Updates the velocity field.
 
         Args:
             t: the time
diff --git a/src/festim/boundary_conditions/dirichlet_bc.py b/src/festim/boundary_conditions/dirichlet_bc.py
index f947ab98c..4ae5604cd 100644
--- a/src/festim/boundary_conditions/dirichlet_bc.py
+++ b/src/festim/boundary_conditions/dirichlet_bc.py
@@ -14,9 +14,7 @@
 
 
 class DirichletBCBase:
-    """
-    Dirichlet boundary condition class
-    u = value
+    """Dirichlet boundary condition class u = value.
 
     Args:
         subdomain: The surface subdomain where the boundary condition is applied
@@ -28,7 +26,6 @@ class DirichletBCBase:
         value_fenics: The value of the boundary condition in fenics format
         bc_expr: The expression of the boundary condition that is used to
             update the `value_fenics`
-
     """
 
     subdomain: _subdomain.SurfaceSubdomain
@@ -67,14 +64,14 @@ def value_fenics(self, value: None | fem.Function | fem.Constant | np.ndarray):
         if not isinstance(value, (fem.Function, fem.Constant, np.ndarray)):
             # FIXME: Should we allow sending in a callable here?
             raise TypeError(
-                "Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, or a np.ndarray not"
+                "Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, or a np.ndarray not"  # noqa: E501
                 + f"{type(value)}"
             )
         self._value_fenics = value
 
     @property
     def time_dependent(self) -> bool:
-        """Returns true if the value of the boundary condition is time dependent"""
+        """Returns true if the value of the boundary condition is time dependent."""
         if self.value is None:
             return False
         if isinstance(self.value, fem.Constant):
@@ -92,17 +89,20 @@ def define_surface_subdomain_dofs(
     ) -> npt.NDArray[np.int32] | tuple[npt.NDArray[np.int32], npt.NDArray[np.int32]]:
         """Defines the facets and the degrees of freedom of the boundary condition.
 
-        Given the input meshtags, find all facets matching the boundary condition subdomain ID,
+        Given the input meshtags, find all facets matching the boundary
+        condition subdomain ID,
         and locate all DOFs associated with the input function space(s).
 
         Note:
-            For sub-spaces, a tuple of sub-spaces are expected as input, and a tuple of arrays
+            For sub-spaces, a tuple of sub-spaces are expected as input, and a
+            tuple of arrays
             associated to each of the function spaces are returned.
 
         Args:
             facet_meshtags: MeshTags describing some facets in the domain
             mesh:
-            function_space: The function space or a tuple of function spaces: (sub, collapsed)
+            function_space: The function space or a tuple of function spaces:
+            (sub, collapsed)
         """
         mesh = (
             function_space[0].mesh
@@ -111,11 +111,11 @@ def define_surface_subdomain_dofs(
         )
         if facet_meshtags.topology != mesh.topology._cpp_object:
             raise ValueError(
-                "Mesh of function-space is not the same as the one used for the meshtags"
+                "Mesh of function-space is not the same as the one used for the meshtags"  # noqa: E501
             )
         if mesh.topology.dim - 1 != facet_meshtags.dim:
             raise ValueError(
-                f"Meshtags of dimension {facet_meshtags.dim}, expected {mesh.topology.dim - 1}"
+                f"Meshtags of dimension {facet_meshtags.dim}, expected {mesh.topology.dim - 1}"  # noqa: E501
             )
         bc_dofs = fem.locate_dofs_topological(
             function_space, facet_meshtags.dim, facet_meshtags.find(self.subdomain.id)
@@ -124,7 +124,7 @@ def define_surface_subdomain_dofs(
         return bc_dofs
 
     def update(self, t: float):
-        """Updates the boundary condition value
+        """Updates the boundary condition value.
 
         Args:
             t: the time
@@ -144,11 +144,13 @@ class FixedConcentrationBC(DirichletBCBase):
     Args:
         subdomain (festim.Subdomain): the surface subdomain where the boundary
             condition is applied
-        value: The value of the boundary condition. It can be a function of space and/or time
+        value: The value of the boundary condition. It can be a function of
+        space and/or time
         species: The name of the species
 
     Attributes:
-        temperature_dependent (bool): True if the value of the bc is temperature dependent
+        temperature_dependent (bool): True if the value of the bc is
+        temperature dependent
 
     Examples:
 
@@ -202,17 +204,18 @@ def create_value(
         K_S: fem.Function = None,
     ):
         """Creates the value of the boundary condition as a fenics object and sets it to
-        self.value_fenics.
-        If the value is a constant, it is converted to a `dolfinx.fem.Constant`.
-        If the value is a function of t, it is converted to  `dolfinx.fem.Constant`.
-        Otherwise, it is converted to a `dolfinx.fem.Function`.Function and the
-        expression of the function is stored in `bc_expr`.
+        self.value_fenics. If the value is a constant, it is converted to a
+        `dolfinx.fem.Constant`. If the value is a function of t, it is converted to
+        `dolfinx.fem.Constant`. Otherwise, it is converted to a
+        `dolfinx.fem.Function`.Function and the expression of the function is stored in
+        `bc_expr`.
 
         Args:
             function_space: the function space
             temperature: The temperature
             t: the time
-            K_S: The solubility of the species. If provided, the value of the boundary condition
+            K_S: The solubility of the species. If provided, the value of the
+            boundary condition
                 is divided by K_S (change of variable method).
         """
         mesh = function_space.mesh
@@ -313,11 +316,10 @@ def create_value(
 class FixedTemperatureBC(DirichletBCBase):
     def create_value(self, function_space: fem.FunctionSpace, t: fem.Constant):
         """Creates the value of the boundary condition as a fenics object and sets it to
-        self.value_fenics.
-        If the value is a constant, it is converted to a `dolfinx.fem.Constant`.
-        If the value is a function of t, it is converted to a `dolfinx.fem.Constant`.
-        Otherwise, it is converted to a` dolfinx.fem.Function` and the
-        expression of the function is stored in `bc_expr`.
+        self.value_fenics. If the value is a constant, it is converted to a
+        `dolfinx.fem.Constant`. If the value is a function of t, it is converted to a
+        `dolfinx.fem.Constant`. Otherwise, it is converted to a` dolfinx.fem.Function`
+        and the expression of the function is stored in `bc_expr`.
 
         Args:
             function_space: the function space
diff --git a/src/festim/boundary_conditions/flux_bc.py b/src/festim/boundary_conditions/flux_bc.py
index 5ae48c777..09eebe6a8 100644
--- a/src/festim/boundary_conditions/flux_bc.py
+++ b/src/festim/boundary_conditions/flux_bc.py
@@ -7,8 +7,7 @@
 
 
 class FluxBCBase:
-    """
-    Flux boundary condition class
+    """Flux boundary condition class.
 
     Ensuring the gradient of the solution u at a boundary:
 
@@ -30,11 +29,12 @@ class FluxBCBase:
         subdomain (festim.SurfaceSubdomain): the surface subdomain where the boundary
             condition is applied
         value (float, fem.Constant, callable): the value of the boundary condition
-        value_fenics (fem.Function or fem.Constant): the value of the boundary condition in
+        value_fenics (fem.Function or fem.Constant): the value of the boundary
+        condition in
             fenics format
-        bc_expr (fem.Expression): the expression of the boundary condition that is used to
+        bc_expr (fem.Expression): the expression of the boundary condition that
+        is used to
             update the value_fenics
-
     """
 
     def __init__(self, subdomain: SurfaceSubdomain, value):
@@ -57,7 +57,8 @@ def value_fenics(self, value):
             value, fem.Function | fem.Constant | np.ndarray | ufl.core.expr.Expr
         ):
             raise TypeError(
-                f"Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, np.ndarray or ufl.core.expr.Expr not {type(value)}"  # noqa: E501
+                f"Value must be a dolfinx.fem.Function, dolfinx.fem.Constant,"
+                f" np.ndarray or ufl.core.expr.Expr not {type(value)}"
             )
         self._value_fenics = value
 
@@ -87,10 +88,9 @@ def temperature_dependent(self):
 
     def create_value_fenics(self, mesh, temperature, t: fem.Constant):
         """Creates the value of the boundary condition as a fenics object and sets it to
-        self.value_fenics.
-        If the value is a constant, it is converted to a fenics.Constant.
-        If the value is a function of t, it is converted to a fenics.Constant.
-        Otherwise, it is converted to a ufl Expression
+        self.value_fenics. If the value is a constant, it is converted to a
+        fenics.Constant. If the value is a function of t, it is converted to a
+        fenics.Constant. Otherwise, it is converted to a ufl Expression.
 
         Args:
             mesh (dolfinx.mesh.Mesh) : the mesh
@@ -109,7 +109,7 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant):
                 # only t is an argument
                 if not isinstance(self.value(t=float(t)), (float, int)):
                     raise ValueError(
-                        f"self.value should return a float or an int, not {type(self.value(t=float(t)))} "
+                        f"self.value should return a float or an int, not {type(self.value(t=float(t)))} "  # noqa: E501
                     )
                 self.value_fenics = F.as_fenics_constant(
                     mesh=mesh, value=self.value(t=float(t))
@@ -126,7 +126,7 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant):
                 self.value_fenics = self.value(**kwargs)
 
     def update(self, t):
-        """Updates the flux bc value
+        """Updates the flux bc value.
 
         Args:
             t (float): the time
@@ -138,22 +138,27 @@ def update(self, t):
 
 
 class ParticleFluxBC(FluxBCBase):
-    """
-    Particle flux boundary condition class
+    """Particle flux boundary condition class.
+
     Ensuring the gradient of the solution c at a boundary:
     -D * grad(c) * n = f
-    where D is the material diffusivity, n is the outwards normal vector of the boundary, f is a function of space and time.
+    where D is the material diffusivity, n is the outwards normal vector of the
+    boundary, f is a function of space and time.
 
     Args:
-        subdomain (festim.SurfaceSubdomain): the surface subdomain where the particle flux
+        subdomain (festim.SurfaceSubdomain): the surface subdomain where the
+        particle flux
             is applied
         value (float, fem.Constant, callable): the value of the particle flux
         species (festim.Species): the species to which the flux is applied
-        species_dependent_value (dict): a dictionary containing the species that the value. Example: {"name": species}
-            where "name" is the variable name in the callable value and species is a festim.Species object.
+        species_dependent_value (dict): a dictionary containing the species
+        that the value. Example: {"name": species}
+            where "name" is the variable name in the callable value and species
+            is a festim.Species object.
 
     Attributes:
-        subdomain (festim.SurfaceSubdomain): the surface subdomain where the particle flux
+        subdomain (festim.SurfaceSubdomain): the surface subdomain where the
+        particle flux
             is applied
         value (float or fem.Constant): the value of the particle flux
         species (festim.Species): the species to which the flux is applied
@@ -161,8 +166,10 @@ class ParticleFluxBC(FluxBCBase):
             fenics format
         bc_expr (fem.Expression): the expression of the particle flux that is used to
             update the value_fenics
-        species_dependent_value (dict): a dictionary containing the species that the value. Example: {"name": species}
-            where "name" is the variable name in the callable value and species is a festim.Species object.
+        species_dependent_value (dict): a dictionary containing the species
+        that the value. Example: {"name": species}
+            where "name" is the variable name in the callable value and species
+            is a festim.Species object.
 
 
     Examples:
@@ -176,11 +183,14 @@ class ParticleFluxBC(FluxBCBase):
         .. testcode:: ParticleFluxBC
 
             ParticleFluxBC(subdomain=my_subdomain, value=1, species="H")
-            ParticleFluxBC(subdomain=my_subdomain, value=lambda x: 1 + x[0], species="H")
+            ParticleFluxBC(subdomain=my_subdomain, value=lambda x: 1 + x[0],
+            species="H")
             ParticleFluxBC(subdomain=my_subdomain, value=lambda t: 1 + t, species="H")
             ParticleFluxBC(subdomain=my_subdomain, value=lambda T: 1 + T, species="H")
-            ParticleFluxBC(subdomain=my_subdomain, value=lambda x, t: 1 + x[0] + t, species="H")
-            ParticleFluxBC(subdomain=my_subdomain, value=lambda c1: 2 * c1**2, species="H", species_dependent_value={"c1": species1})
+            ParticleFluxBC(subdomain=my_subdomain, value=lambda x, t: 1 + x[0]
+            + t, species="H")
+            ParticleFluxBC(subdomain=my_subdomain, value=lambda c1: 2 * c1**2,
+            species="H", species_dependent_value={"c1": species1})
     """
 
     def __init__(self, subdomain, value, species, species_dependent_value={}):
@@ -191,10 +201,9 @@ def __init__(self, subdomain, value, species, species_dependent_value={}):
 
     def create_value_fenics(self, mesh, temperature, t: fem.Constant):
         """Creates the value of the boundary condition as a fenics object and sets it to
-        self.value_fenics.
-        If the value is a constant, it is converted to a fenics.Constant.
-        If the value is a function of t, it is converted to a fenics.Constant.
-        Otherwise, it is converted to a ufl Expression
+        self.value_fenics. If the value is a constant, it is converted to a
+        fenics.Constant. If the value is a function of t, it is converted to a
+        fenics.Constant. Otherwise, it is converted to a ufl Expression.
 
         Args:
             mesh (dolfinx.mesh.Mesh) : the mesh
@@ -213,7 +222,7 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant):
                 # only t is an argument
                 if not isinstance(self.value(t=float(t)), (float, int)):
                     raise ValueError(
-                        f"self.value should return a float or an int, not {type(self.value(t=float(t)))} "
+                        f"self.value should return a float or an int, not {type(self.value(t=float(t)))} "  # noqa: E501
                     )
                 self.value_fenics = F.as_fenics_constant(
                     mesh=mesh, value=self.value(t=float(t))
@@ -239,11 +248,12 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant):
 
 
 class HeatFluxBC(FluxBCBase):
-    """
-    Heat flux boundary condition class
+    """Heat flux boundary condition class.
+
     Ensuring the gradient of the solution T at a boundary:
     -lambda * grad(T) * n = f
-    where lambda is the thermal conductivity , n is the outwards normal vector of the boundary, f is a function of space and time.
+    where lambda is the thermal conductivity , n is the outwards normal vector
+    of the boundary, f is a function of space and time.
 
     Args:
         subdomain (festim.SurfaceSubdomain): the surface subdomain where the heat flux
diff --git a/src/festim/boundary_conditions/henrys_bc.py b/src/festim/boundary_conditions/henrys_bc.py
index c92083f0f..071ca08e8 100644
--- a/src/festim/boundary_conditions/henrys_bc.py
+++ b/src/festim/boundary_conditions/henrys_bc.py
@@ -5,14 +5,13 @@
 
 
 def henrys_law(T, H_0, E_H, pressure):
-    """Applies the Henry's law to compute the concentration at the boundary"""
+    """Applies the Henry's law to compute the concentration at the boundary."""
     H = H_0 * ufl.exp(-E_H / k_B / T)
     return H * pressure
 
 
 class HenrysBC(FixedConcentrationBC):
-    """
-    Henrys boundary condition class
+    """Henrys boundary condition class.
 
     c = H * pressure
     H = H_0 * exp(-E_H / k_B / T)
@@ -21,7 +20,8 @@ class HenrysBC(FixedConcentrationBC):
         subdomain (festim.Subdomain): the subdomain where the boundary
             condition is applied
         species (str): the name of the species
-        H_0 (float or fem.Constant): the Henrys constant pre-exponential factor (H/m3/Pa)
+        H_0 (float or fem.Constant): the Henrys constant pre-exponential factor
+        (H/m3/Pa)
         E_H (float or fem.Constant): the Henrys constant activation energy (eV)
         pressure (float or callable): the pressure at the boundary (Pa)
 
@@ -30,7 +30,8 @@ class HenrysBC(FixedConcentrationBC):
             condition is applied
         value (float or fem.Constant): the value of the boundary condition
         species (festim.Species or str): the name of the species
-        H_0 (float or fem.Constant): the Henrys constant pre-exponential factor (H/m3/Pa)
+        H_0 (float or fem.Constant): the Henrys constant pre-exponential factor
+        (H/m3/Pa)
         E_H (float or fem.Constant): the Henrys constant activation energy (eV)
         pressure (float or callable): the pressure at the boundary (Pa)
 
@@ -43,11 +44,16 @@ class HenrysBC(FixedConcentrationBC):
 
         .. testcode:: HenrysBC
 
-            HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=1e5, species="H")
-            HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda x: 1e5 + x[0], species="H")
-            HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda t: 1e5 + t, species="H")
-            HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda T: 1e5 + T, species="H")
-            HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda x, t: 1e5 + x[0] + t, species="H")
+            HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=1e5,
+            species="H")
+            HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda
+            x: 1e5 + x[0], species="H")
+            HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda
+            t: 1e5 + t, species="H")
+            HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda
+            T: 1e5 + T, species="H")
+            HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda
+            x, t: 1e5 + x[0] + t, species="H")
     """
 
     def __init__(self, subdomain, H_0, E_H, pressure, species) -> None:
@@ -62,7 +68,7 @@ def __init__(self, subdomain, H_0, E_H, pressure, species) -> None:
         super().__init__(value=value, species=species, subdomain=subdomain)
 
     def create_new_value_function(self):
-        """Creates a new value function based on the pressure attribute
+        """Creates a new value function based on the pressure attribute.
 
         Raises:
             ValueError: if the pressure function is not supported
diff --git a/src/festim/boundary_conditions/sieverts_bc.py b/src/festim/boundary_conditions/sieverts_bc.py
index 911d5dcf9..d15b317b7 100644
--- a/src/festim/boundary_conditions/sieverts_bc.py
+++ b/src/festim/boundary_conditions/sieverts_bc.py
@@ -5,14 +5,13 @@
 
 
 def sieverts_law(T, S_0, E_S, pressure):
-    """Applies the Sieverts law to compute the concentration at the boundary"""
+    """Applies the Sieverts law to compute the concentration at the boundary."""
     S = S_0 * ufl.exp(-E_S / k_B / T)
     return S * pressure**0.5
 
 
 class SievertsBC(FixedConcentrationBC):
-    """
-    Sieverts boundary condition class
+    """Sieverts boundary condition class.
 
     c = S * sqrt(pressure)
     S = S_0 * exp(-E_S / k_B / T)
@@ -21,7 +20,8 @@ class SievertsBC(FixedConcentrationBC):
         subdomain (festim.Subdomain): the subdomain where the boundary
             condition is applied
         species (str): the name of the species
-        S_0 (float or fem.Constant): the Sieverts constant pre-exponential factor (H/m3/Pa0.5)
+        S_0 (float or fem.Constant): the Sieverts constant pre-exponential
+        factor (H/m3/Pa0.5)
         E_S (float or fem.Constant): the Sieverts constant activation energy (eV)
         pressure (float or callable): the pressure at the boundary (Pa)
 
@@ -30,7 +30,8 @@ class SievertsBC(FixedConcentrationBC):
             condition is applied
         value (float or fem.Constant): the value of the boundary condition
         species (festim.Species or str): the name of the species
-        S_0 (float or fem.Constant): the Sieverts constant pre-exponential factor (H/m3/Pa0.5)
+        S_0 (float or fem.Constant): the Sieverts constant pre-exponential
+        factor (H/m3/Pa0.5)
         E_S (float or fem.Constant): the Sieverts constant activation energy (eV)
         pressure (float or callable): the pressure at the boundary (Pa)
 
@@ -43,11 +44,16 @@ class SievertsBC(FixedConcentrationBC):
 
         .. testcode:: SievertsBC
 
-            SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, pressure=1e5, species="H")
-            SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, pressure=lambda x: 1e5 + x[0], species="H")
-            SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, pressure=lambda t: 1e5 + t, species="H")
-            SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, pressure=lambda T: 1e5 + T, species="H")
-            SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, pressure=lambda x, t: 1e5 + x[0] + t, species="H")
+            SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, pressure=1e5,
+            species="H")
+            SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2,
+            pressure=lambda x: 1e5 + x[0], species="H")
+            SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2,
+            pressure=lambda t: 1e5 + t, species="H")
+            SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2,
+            pressure=lambda T: 1e5 + T, species="H")
+            SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2,
+            pressure=lambda x, t: 1e5 + x[0] + t, species="H")
     """
 
     def __init__(self, subdomain, S_0, E_S, pressure, species) -> None:
@@ -62,7 +68,7 @@ def __init__(self, subdomain, S_0, E_S, pressure, species) -> None:
         super().__init__(value=value, species=species, subdomain=subdomain)
 
     def create_new_value_function(self):
-        """Creates a new value function based on the pressure attribute
+        """Creates a new value function based on the pressure attribute.
 
         Raises:
             ValueError: if the pressure function is not supported
diff --git a/src/festim/boundary_conditions/surface_reaction.py b/src/festim/boundary_conditions/surface_reaction.py
index ca8efe92e..f56e4c56e 100644
--- a/src/festim/boundary_conditions/surface_reaction.py
+++ b/src/festim/boundary_conditions/surface_reaction.py
@@ -6,8 +6,8 @@
 
 
 class SurfaceReactionBCpartial(ParticleFluxBC):
-    """Boundary condition representing a surface reaction
-    A + B <-> C
+    """Boundary condition representing a surface reaction A + B <-> C.
+
     where A, B are the reactants and C is the product
     the forward reaction rate is K_r = k_r0 * exp(-E_kr / (k_B * T))
     and the backward reaction rate is K_d = k_d0 * exp(-E_kd / (k_B * T))
@@ -79,8 +79,8 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant):
 
 
 class SurfaceReactionBC:
-    """Boundary condition representing a surface reaction
-    A + B <-> C
+    """Boundary condition representing a surface reaction A + B <-> C.
+
     where A, B are the reactants and C is the product
     the forward reaction rate is K_r = k_r0 * exp(-E_kr / (k_B * T))
     and the backward reaction rate is K_d = k_d0 * exp(-E_kd / (k_B * T))
diff --git a/src/festim/coupled_heat_hydrogen_problem.py b/src/festim/coupled_heat_hydrogen_problem.py
index 64ce4cfec..d7dab6f10 100644
--- a/src/festim/coupled_heat_hydrogen_problem.py
+++ b/src/festim/coupled_heat_hydrogen_problem.py
@@ -11,8 +11,7 @@
 
 
 class CoupledTransientHeatTransferHydrogenTransport:
-    """
-    Coupled heat transfer and hydrogen transport transient problem
+    """Coupled heat transfer and hydrogen transport transient problem.
 
     Args:
         heat_problem: the heat transfer problem
@@ -37,8 +36,6 @@ class CoupledTransientHeatTransferHydrogenTransport:
                 heat_problem=my_heat_transfer_model,
                 hydrogen_problem=my_h_transport_model,
             )
-
-
     """
 
     heat_problem: HeatTransferProblem
diff --git a/src/festim/exports/average_surface.py b/src/festim/exports/average_surface.py
index ac3354d06..da6b46b5e 100644
--- a/src/festim/exports/average_surface.py
+++ b/src/festim/exports/average_surface.py
@@ -6,7 +6,7 @@
 
 
 class AverageSurface(SurfaceQuantity):
-    """Computes the average value of a field on a given surface
+    """Computes the average value of a field on a given surface.
 
     Args:
         field (festim.Species): species for which the average surface is computed
@@ -25,9 +25,8 @@ def title(self):
     def compute(
         self, u: fem.Function | ufl.indexed.Indexed, ds: ufl.Measure, entity_maps=None
     ):
-        """
-        Computes the average value of the field on the defined surface
-        subdomain, and appends it to the data list
+        """Computes the average value of the field on the defined surface subdomain, and
+        appends it to the data list.
 
         Args:
             u: field for which the average value is computed
diff --git a/src/festim/exports/average_volume.py b/src/festim/exports/average_volume.py
index f2ff01a8d..01feeae9e 100644
--- a/src/festim/exports/average_volume.py
+++ b/src/festim/exports/average_volume.py
@@ -4,12 +4,13 @@
 
 
 class AverageVolume(VolumeQuantity):
-    """Computes the average value of a field in a given volume
+    """Computes the average value of a field in a given volume.
 
     Args:
         field (festim.Species): species for which the average volume is computed
         volume (festim.VolumeSubdomain): volume subdomain
-        filename (str, optional): name of the file to which the average volume is exported
+        filename (str, optional): name of the file to which the average volume
+        is exported
 
     Attributes:
         see `festim.VolumeQuantity`
@@ -20,10 +21,8 @@ def title(self):
         return f"Average {self.field.name} volume {self.volume.id}"
 
     def compute(self, u, dx, entity_maps=None):
-        """
-        Computes the average value of solution function within the defined volume
-        subdomain, and appends it to the data list
-        """
+        """Computes the average value of solution function within the defined volume
+        subdomain, and appends it to the data list."""
         self.value = assemble_scalar(
             u * dx(self.volume.id), entity_maps=entity_maps
         ) / assemble_scalar(1 * dx(self.volume.id))
diff --git a/src/festim/exports/maximum_surface.py b/src/festim/exports/maximum_surface.py
index 25f45d5f7..1a29f5d7e 100644
--- a/src/festim/exports/maximum_surface.py
+++ b/src/festim/exports/maximum_surface.py
@@ -7,12 +7,13 @@
 
 
 class MaximumSurface(SurfaceQuantity):
-    """Computes the maximum value of a field on a given surface
+    """Computes the maximum value of a field on a given surface.
 
     Args:
         field (festim.Species): species for which the maximum surface is computed
         surface (festim.SurfaceSubdomain): surface subdomain
-        filename (str, optional): name of the file to which the maximum surface is exported
+        filename (str, optional): name of the file to which the maximum surface
+        is exported
 
     Attributes:
         see `festim.SurfaceQuantity`
@@ -23,10 +24,8 @@ def title(self):
         return f"Maximum {self.field.name} surface {self.surface.id}"
 
     def compute(self):
-        """
-        Computes the maximum value of the field on the defined surface
-        subdomain, and appends it to the data list
-        """
+        """Computes the maximum value of the field on the defined surface subdomain, and
+        appends it to the data list."""
         solution = self.field.post_processing_solution
         entities = self.facet_meshtags.find(self.surface.id)
         if isinstance(solution, dolfinx.fem.Function):
diff --git a/src/festim/exports/maximum_volume.py b/src/festim/exports/maximum_volume.py
index 043b4e284..d3bb257e3 100644
--- a/src/festim/exports/maximum_volume.py
+++ b/src/festim/exports/maximum_volume.py
@@ -7,16 +7,16 @@
 
 
 class MaximumVolume(VolumeQuantity):
-    """Computes the maximum value of a field in a given volume
+    """Computes the maximum value of a field in a given volume.
 
     Args:
         field (festim.Species): species for which the maximum volume is computed
         volume (festim.VolumeSubdomain): volume subdomain
-        filename (str, optional): name of the file to which the maximum volume is exported
+        filename (str, optional): name of the file to which the maximum volume
+        is exported
 
     Attributes:
         see `festim.VolumeQuantity`
-
     """
 
     @property
@@ -24,10 +24,8 @@ def title(self):
         return f"Maximum {self.field.name} volume {self.volume.id}"
 
     def compute(self):
-        """
-        Computes the maximum value of solution function within the defined volume
-        subdomain, and appends it to the data list
-        """
+        """Computes the maximum value of solution function within the defined volume
+        subdomain, and appends it to the data list."""
         solution = self.field.post_processing_solution
         entities = self.volume_meshtags.find(self.volume.id)
 
diff --git a/src/festim/exports/minimum_surface.py b/src/festim/exports/minimum_surface.py
index 454f6b763..1d2e934d5 100644
--- a/src/festim/exports/minimum_surface.py
+++ b/src/festim/exports/minimum_surface.py
@@ -7,12 +7,13 @@
 
 
 class MinimumSurface(SurfaceQuantity):
-    """Computes the minimum value of a field on a given surface
+    """Computes the minimum value of a field on a given surface.
 
     Args:
         field (festim.Species): species for which the minimum surface is computed
         surface (festim.SurfaceSubdomain): surface subdomain
-        filename (str, optional): name of the file to which the minimum surface is exported
+        filename (str, optional): name of the file to which the minimum surface
+        is exported
 
     Attributes:
         see `festim.SurfaceQuantity`
@@ -23,10 +24,8 @@ def title(self):
         return f"Minimum {self.field.name} surface {self.surface.id}"
 
     def compute(self):
-        """
-        Computes the minimum value of the field on the defined surface
-        subdomain, and appends it to the data list
-        """
+        """Computes the minimum value of the field on the defined surface subdomain, and
+        appends it to the data list."""
         solution = self.field.post_processing_solution
         entities = self.facet_meshtags.find(self.surface.id)
         if isinstance(solution, dolfinx.fem.Function):
diff --git a/src/festim/exports/minimum_volume.py b/src/festim/exports/minimum_volume.py
index 730b3bf29..e17ce5aff 100644
--- a/src/festim/exports/minimum_volume.py
+++ b/src/festim/exports/minimum_volume.py
@@ -7,12 +7,13 @@
 
 
 class MinimumVolume(VolumeQuantity):
-    """Computes the minimum value of a field in a given volume
+    """Computes the minimum value of a field in a given volume.
 
     Args:
         field (festim.Species): species for which the minimum volume is computed
         volume (festim.VolumeSubdomain): volume subdomain
-        filename (str, optional): name of the file to which the minimum volume is exported
+        filename (str, optional): name of the file to which the minimum volume
+        is exported
 
     Attributes:
         see `festim.VolumeQuantity`
@@ -23,10 +24,8 @@ def title(self):
         return f"Minimum {self.field.name} volume {self.volume.id}"
 
     def compute(self):
-        """
-        Computes the minimum value of solution function within the defined volume
-        subdomain, and appends it to the data list
-        """
+        """Computes the minimum value of solution function within the defined volume
+        subdomain, and appends it to the data list."""
         solution = self.field.post_processing_solution
         entities = self.volume_meshtags.find(self.volume.id)
 
diff --git a/src/festim/exports/surface_flux.py b/src/festim/exports/surface_flux.py
index fe70e8f8c..b2ff9de21 100644
--- a/src/festim/exports/surface_flux.py
+++ b/src/festim/exports/surface_flux.py
@@ -8,7 +8,7 @@
 
 
 class SurfaceFlux(SurfaceQuantity):
-    """Computes the flux of a field on a given surface
+    """Computes the flux of a field on a given surface.
 
     Args:
         field: species for which the surface flux is computed
@@ -39,7 +39,7 @@ def __init__(
     def compute(
         self, u: fem.Function | ufl.indexed.Indexed, ds: ufl.Measure, entity_maps=None
     ):
-        """Computes the value of the flux at the surface
+        """Computes the value of the flux at the surface.
 
         Args:
             u: field for which the flux is computed
diff --git a/src/festim/exports/surface_quantity.py b/src/festim/exports/surface_quantity.py
index b6f9d2726..515ab368e 100644
--- a/src/festim/exports/surface_quantity.py
+++ b/src/festim/exports/surface_quantity.py
@@ -5,7 +5,7 @@
 
 
 class SurfaceQuantity:
-    """Export SurfaceQuantity
+    """Export SurfaceQuantity.
 
     Args:
         field: species for which the surface flux is computed
@@ -75,8 +75,8 @@ def field(self, value):
         self._field = value
 
     def write(self, t):
-        """If the filename doesnt exist yet, create it and write the header,
-        then append the time and value to the file"""
+        """If the filename doesnt exist yet, create it and write the header, then append
+        the time and value to the file."""
 
         if self.filename is not None:
             if self._first_time_export:
diff --git a/src/festim/exports/total_surface.py b/src/festim/exports/total_surface.py
index d3f51327a..21719575c 100644
--- a/src/festim/exports/total_surface.py
+++ b/src/festim/exports/total_surface.py
@@ -6,7 +6,7 @@
 
 
 class TotalSurface(SurfaceQuantity):
-    """Computes the total value of a field on a given surface
+    """Computes the total value of a field on a given surface.
 
     Args:
         field (`festim.Species`): species for which the total volume is computed
@@ -24,9 +24,8 @@ def title(self):
     def compute(
         self, u: fem.Function | ufl.indexed.Indexed, ds: ufl.Measure, entity_maps=None
     ):
-        """
-        Computes the total value of the field on the defined surface
-        subdomain, and appends it to the data list
+        """Computes the total value of the field on the defined surface subdomain, and
+        appends it to the data list.
 
         Args:
             u: field for which the total value is computed
diff --git a/src/festim/exports/total_volume.py b/src/festim/exports/total_volume.py
index 312499421..194ce4ec6 100644
--- a/src/festim/exports/total_volume.py
+++ b/src/festim/exports/total_volume.py
@@ -5,7 +5,7 @@
 
 
 class TotalVolume(VolumeQuantity):
-    """Computes the total value of a field in a given volume
+    """Computes the total value of a field in a given volume.
 
     Args:
         field (festim.Species): species for which the total volume is computed
@@ -21,9 +21,8 @@ def title(self):
         return f"Total {self.field.name} volume {self.volume.id}"
 
     def compute(self, u, dx: ufl.Measure, entity_maps=None):
-        """
-        Computes the value of the total volume of the field in the volume subdomain
-        and appends it to the data list
+        """Computes the value of the total volume of the field in the volume subdomain
+        and appends it to the data list.
 
         Args:
             u: field for which the total volume is computed
diff --git a/src/festim/exports/volume_quantity.py b/src/festim/exports/volume_quantity.py
index aed892a91..1daaaac6a 100644
--- a/src/festim/exports/volume_quantity.py
+++ b/src/festim/exports/volume_quantity.py
@@ -5,7 +5,7 @@
 
 
 class VolumeQuantity:
-    """Export VolumeQuantity
+    """Export VolumeQuantity.
 
     Args:
         field: species for which the volume quantity is computed
@@ -63,8 +63,8 @@ def field(self, value):
         self._field = value
 
     def write(self, t):
-        """If the filename doesnt exist yet, create it and write the header,
-        then append the time and value to the file"""
+        """If the filename doesnt exist yet, create it and write the header, then append
+        the time and value to the file."""
 
         if self.filename is not None:
             if self._first_time_export:
diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py
index 7159360fd..87ad6a631 100644
--- a/src/festim/exports/vtx.py
+++ b/src/festim/exports/vtx.py
@@ -8,7 +8,7 @@
 
 
 class ExportBaseClass:
-    """Export functions to VTX file
+    """Export functions to VTX file.
 
     Args:
         filename: The name of the output file
@@ -49,7 +49,7 @@ def filename(self):
 
 
 class VTXTemperatureExport(ExportBaseClass):
-    """Export temperature field functions to VTX file
+    """Export temperature field functions to VTX file.
 
     Args:
         filename: The name of the output file
@@ -74,7 +74,7 @@ def __init__(
 
 
 class VTXSpeciesExport(ExportBaseClass):
-    """Export species field functions to VTX file
+    """Export species field functions to VTX file.
 
     Args:
         filename: The name of the output file
@@ -117,8 +117,7 @@ def field(self) -> list[Species]:
 
     @field.setter
     def field(self, value: Species | list[Species]):
-        """
-        Update the field to export.
+        """Update the field to export.
 
         Args:
             value: The species to export
@@ -149,9 +148,9 @@ def field(self, value: Species | list[Species]):
         self._field = val
 
     def get_functions(self) -> list[fem.Function]:
-        """
-        Returns list of species for a given subdomain. If using legacy mode, return the
-        whole species.
+        """Returns list of species for a given subdomain.
+
+        If using legacy mode, return the whole species.
         """
 
         legacy_output: bool = False
diff --git a/src/festim/exports/xdmf.py b/src/festim/exports/xdmf.py
index 727bf0a10..bd86171bc 100644
--- a/src/festim/exports/xdmf.py
+++ b/src/festim/exports/xdmf.py
@@ -10,7 +10,7 @@
 
 
 class XDMFExport(ExportBaseClass):
-    """Export functions to XDMFfile
+    """Export functions to XDMFfile.
 
     Args:
         filename: The name of the output file
@@ -43,20 +43,20 @@ def field(self, value: Species | list[Species]):
             for element in value:
                 if not isinstance(element, Species):
                     raise TypeError(
-                        f"Each element in the list must be a species, got {type(element)}."
+                        f"Each element in the list must be a species, got {type(element)}."  # noqa: E501
                     )
             val = value
         elif isinstance(value, Species):
             val = [value]
         else:
             raise TypeError(
-                f"field must be of type festim.Species or a list of festim.Species, got "
+                f"field must be of type festim.Species or a list of festim.Species, got "  # noqa: E501
                 f"{type(value)}."
             )
         self._field = val
 
     def define_writer(self, comm: mpi4py.MPI.Intracomm) -> None:
-        """Define the writer
+        """Define the writer.
 
         Args:
             comm (mpi4py.MPI.Intracomm): the MPI communicator
@@ -64,7 +64,7 @@ def define_writer(self, comm: mpi4py.MPI.Intracomm) -> None:
         self._writer = XDMFFile(comm, self.filename, "w")
 
     def write(self, t: float):
-        """Write functions to VTX file
+        """Write functions to VTX file.
 
         Args:
             t (float): the time of export
diff --git a/src/festim/heat_transfer_problem.py b/src/festim/heat_transfer_problem.py
index 600b95e6b..710bdcc8b 100644
--- a/src/festim/heat_transfer_problem.py
+++ b/src/festim/heat_transfer_problem.py
@@ -79,10 +79,12 @@ def initialise(self):
         self.initialise_exports()
 
     def define_function_space(self):
-        """Creates the function space of the model, creates a mixed element if
-        model is multispecies. Creates the main solution and previous solution
-        function u and u_n. Create global DG function spaces of degree 0 and 1
-        for the global diffusion coefficient"""
+        """Creates the function space of the model, creates a mixed element if model is
+        multispecies.
+
+        Creates the main solution and previous solution function u and u_n. Create
+        global DG function spaces of degree 0 and 1 for the global diffusion coefficient
+        """
 
         degree = 1
         element_CG = basix.ufl.element(
@@ -99,7 +101,7 @@ def define_function_space(self):
         self.test_function = ufl.TestFunction(self.function_space)
 
     def create_dirichletbc_form(self, bc):
-        """Creates a dirichlet boundary condition form
+        """Creates a dirichlet boundary condition form.
 
         Args:
             bc (festim.FixedTemperatureBC): the boundary condition
@@ -128,7 +130,7 @@ def create_dirichletbc_form(self, bc):
         return form
 
     def create_source_values_fenics(self):
-        """For each source create the value_fenics"""
+        """For each source create the value_fenics."""
         for source in self.sources:
             # create value_fenics for all source objects
             source.value.convert_input_value(
@@ -138,7 +140,7 @@ def create_source_values_fenics(self):
             )
 
     def create_flux_values_fenics(self):
-        """For each heat flux create the value_fenics"""
+        """For each heat flux create the value_fenics."""
         for bc in self.boundary_conditions:
             # create value_fenics for all F.HeatFluxBC objects
             if isinstance(bc, boundary_conditions.HeatFluxBC):
@@ -149,8 +151,8 @@ def create_flux_values_fenics(self):
                 )
 
     def create_initial_conditions(self):
-        """For each initial condition, create the value_fenics and assign it to
-        the previous solution of the condition's species"""
+        """For each initial condition, create the value_fenics and assign it to the
+        previous solution of the condition's species."""
 
         if not self.initial_condition:
             return
@@ -171,7 +173,7 @@ def create_initial_conditions(self):
         self.u_n.interpolate(self.initial_condition.expr_fenics)
 
     def create_formulation(self):
-        """Creates the formulation of the model"""
+        """Creates the formulation of the model."""
 
         self.formulation = 0
 
@@ -216,8 +218,8 @@ def create_formulation(self):
                 )
 
     def initialise_exports(self):
-        """Defines the export writers of the model, if field is given as
-        a string, find species object in self.species"""
+        """Defines the export writers of the model, if field is given as a string, find
+        species object in self.species."""
 
         for export in self.exports:
             if isinstance(export, exports.XDMFExport):
@@ -233,7 +235,7 @@ def initialise_exports(self):
                 )
 
     def post_processing(self):
-        """Post processes the model"""
+        """Post processes the model."""
 
         for export in self.exports:
             # skip if it isn't time to export
diff --git a/src/festim/helpers.py b/src/festim/helpers.py
index 0a1ba5290..104a6888c 100644
--- a/src/festim/helpers.py
+++ b/src/festim/helpers.py
@@ -12,7 +12,7 @@
 def as_fenics_constant(
     value: float | int | fem.Constant, mesh: dolfinx.mesh.Mesh
 ) -> fem.Constant:
-    """Converts a value to a dolfinx.Constant
+    """Converts a value to a dolfinx.Constant.
 
     Args:
         value: the value to convert
@@ -41,7 +41,7 @@ def as_mapped_function(
     temperature: fem.Function | fem.Constant | ufl.core.expr.Expr | None = None,
 ) -> ufl.core.expr.Expr:
     """Maps a user given callable function to the mesh, time or temperature within
-    festim as needed
+    festim as needed.
 
     Args:
         value: the callable to convert
@@ -76,7 +76,7 @@ def as_fenics_interp_expr_and_function(
 ) -> tuple[fem.Expression, fem.Function]:
     """Takes a user given callable function, maps the function to the mesh, time or
     temperature within festim as needed. Then creates the fenics interpolation
-    expression and function objects
+    expression and function objects.
 
     Args:
         value: the callable to convert
@@ -104,9 +104,8 @@ def as_fenics_interp_expr_and_function(
 
 
 class Value:
-    """
-    A class to handle input values from users and convert them to a relevent fenics
-    object
+    """A class to handle input values from users and convert them to a relevent fenics
+    object.
 
     Args:
         input_value: The value of the user input
@@ -119,7 +118,6 @@ class Value:
         explicit_time_dependent : True if the user input value is explicitly time
             dependent
         temperature_dependent : True if the user input value is temperature dependent
-
     """
 
     input_value: (
@@ -177,7 +175,7 @@ def input_value(self, value):
 
     @property
     def explicit_time_dependent(self) -> bool:
-        """Returns true if the value given is time dependent"""
+        """Returns true if the value given is time dependent."""
         if self.input_value is None:
             return False
         if isinstance(self.input_value, fem.Constant | ufl.core.expr.Expr):
@@ -190,7 +188,7 @@ def explicit_time_dependent(self) -> bool:
 
     @property
     def temperature_dependent(self) -> bool:
-        """Returns true if the value given is temperature dependent"""
+        """Returns true if the value given is temperature dependent."""
         if self.input_value is None:
             return False
         if isinstance(self.input_value, fem.Constant | ufl.core.expr.Expr):
@@ -208,8 +206,8 @@ def convert_input_value(
         temperature: fem.Function | fem.Constant | ufl.core.expr.Expr | None = None,
         up_to_ufl_expr: bool | None = False,
     ):
-        """Converts a user given value to a relevent fenics object depending
-        on the type of the value provided
+        """Converts a user given value to a relevent fenics object depending on the type
+        of the value provided.
 
         Args:
             function_space: the function space of the fenics object, optional
@@ -264,7 +262,7 @@ def convert_input_value(
                 )
 
     def update(self, t: float):
-        """Updates the value
+        """Updates the value.
 
         Args:
             t: the time
@@ -285,9 +283,13 @@ def update(self, t: float):
 
 # Define the appropriate method based on the version
 if version.parse(dolfinx_version) > version.parse("0.9.0"):
-    get_interpolation_points = lambda element: element.interpolation_points
+
+    def get_interpolation_points(element):
+        return element.interpolation_points
 else:
-    get_interpolation_points = lambda element: element.interpolation_points()
+
+    def get_interpolation_points(element):
+        return element.interpolation_points()
 
 
 def nmm_interpolate(
@@ -322,8 +324,7 @@ def nmm_interpolate(
 def is_it_time_to_export(
     times: list | None, current_time: float, atol=0, rtol=1.0e-5
 ) -> bool:
-    """
-    Checks if the exported field should be written to a file or not based on the
+    """Checks if the exported field should be written to a file or not based on the
     current time and the times in `export.times`
 
     After a successful match, the corresponding time is removed from the list to
@@ -356,7 +357,7 @@ def is_it_time_to_export(
 
 def convergenceTest(snes, it, norms):
     global _residual0
-    xnorm, gnorm, f = norms  # ||x_k||, ||x_k-x_k-1||, ||F(x_k)||
+    _xnorm, gnorm, f = norms  # ||x_k||, ||x_k-x_k-1||, ||F(x_k)||
 
     rtol, atol, stol, max_its = snes.getTolerances()
 
@@ -378,7 +379,7 @@ def convergenceTest(snes, it, norms):
 def SnesMonitor(snes, iter, rnorm):
     global _prev_xnorm
     if MPI.COMM_WORLD.rank == 0:
-        rtol, atol, stol, max_its = snes.getTolerances()
+        rtol, atol, stol, _max_its = snes.getTolerances()
         x = snes.getSolution()
         xnorm = x.norm()
 
@@ -390,7 +391,7 @@ def SnesMonitor(snes, iter, rnorm):
 
         dolfinx.log.log(
             dolfinx.log.LogLevel.INFO,
-            f"SNES {iter=} ; {rnorm=:.5e} ({atol=}) ; {relative_residual=:.5e} ({rtol=}) ; {stepsize_rel=:.5e} ({stol=:.5e})",
+            f"SNES {iter=} ; {rnorm=:.5e} ({atol=}) ; {relative_residual=:.5e} ({rtol=}) ; {stepsize_rel=:.5e} ({stol=:.5e})",  # noqa: E501
         )
 
         # Update previous xnorm
diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py
index 14d9c35b4..51ceee8c3 100644
--- a/src/festim/hydrogen_transport_problem.py
+++ b/src/festim/hydrogen_transport_problem.py
@@ -41,7 +41,6 @@
     is_it_time_to_export,
     nmm_interpolate,
 )
-from festim.material import SolubilityLaw
 
 from .mesh import CoordinateSystem, Mesh
 
@@ -52,8 +51,7 @@
 
 
 class HydrogenTransportProblem(problem.ProblemBase):
-    """
-    Hydrogen Transport Problem.
+    """Hydrogen Transport Problem.
 
     Args:
         mesh: The mesh
@@ -134,7 +132,6 @@ class HydrogenTransportProblem(problem.ProblemBase):
                 species=[F.Species(name="H"), F.Species(name="Trap")],
             )
             my_model.initialise()
-
     """
 
     _temperature_as_function: fem.Function
@@ -282,7 +279,7 @@ def volume_meshtags(self, value):
 
     @property
     def _unpacked_bcs(self):
-        """Returns all boundary conditions, including fluxes from surface reactions"""
+        """Returns all boundary conditions, including fluxes from surface reactions."""
         all_boundary_conditions = []
         for bc in self.boundary_conditions:
             if isinstance(bc, boundary_conditions.SurfaceReactionBC):
@@ -331,7 +328,7 @@ def initialise(self):
         self.initialise_exports()
 
     def create_implicit_species_value_fenics(self):
-        """For each implicit species, create the value_fenics"""
+        """For each implicit species, create the value_fenics."""
         for reaction in self.reactions:
             for reactant in reaction.reactant:
                 if isinstance(reactant, _species.ImplicitSpecies):
@@ -341,7 +338,7 @@ def create_implicit_species_value_fenics(self):
                     )
 
     def create_species_from_traps(self):
-        """Generate a species and reaction per trap defined in self.traps"""
+        """Generate a species and reaction per trap defined in self.traps."""
 
         for trap in self.traps:
             trap.create_species_and_reaction()
@@ -349,12 +346,13 @@ def create_species_from_traps(self):
             self.reactions.append(trap.reaction)
 
     def define_temperature(self):
-        """Sets the value of temperature_fenics_value. The type depends on
-        self.temperature. If self.temperature is a function on t only, create
-        a fem.Constant. Else, create an dolfinx.fem.Expression (stored in
-        self.temperature_expr) to be updated, a dolfinx.fem.Function object
-        is created from the Expression (stored in self.temperature_fenics_value).
-        Raise a ValueError if temperature is None.
+        """Sets the value of temperature_fenics_value.
+
+        The type depends on self.temperature. If self.temperature is a function on t
+        only, create a fem.Constant. Else, create an dolfinx.fem.Expression (stored in
+        self.temperature_expr) to be updated, a dolfinx.fem.Function object is created
+        from the Expression (stored in self.temperature_fenics_value). Raise a
+        ValueError if temperature is None.
         """
         # check if temperature is None
         if self.temperature is None:
@@ -410,8 +408,8 @@ def define_temperature(self):
                 self.temperature_fenics.interpolate(self.temperature_expr)
 
     def initialise_exports(self):
-        """Defines the export writers of the model, if field is given as
-        a string, find species object in self.species"""
+        """Defines the export writers of the model, if field is given as a string, find
+        species object in self.species."""
 
         for export in self.exports:
             if isinstance(export, exports.ExportBaseClass):
@@ -513,9 +511,8 @@ def initialise_exports(self):
                 export.data = []
 
     def _get_temperature_field_as_function(self) -> dolfinx.fem.Function:
-        """
-        Based on the type of the temperature_fenics attribute, converts
-        it as a Function to be used in VTX export
+        """Based on the type of the temperature_fenics attribute, converts it as a
+        Function to be used in VTX export.
 
         Returns:
             the temperature field of the simulation
@@ -537,7 +534,7 @@ def _get_temperature_field_as_function(self) -> dolfinx.fem.Function:
             return temperature_field
 
     def define_D_global(self, species):
-        """Defines the global diffusion coefficient for a given species
+        """Defines the global diffusion coefficient for a given species.
 
         Args:
             species (F.Species): the species
@@ -636,7 +633,7 @@ def define_function_spaces(self, element_degree: int = 1):
     def assign_functions_to_species(self):
         """Creates the solution, prev solution, test function and post-processing
         solution for each species, as well as a collapsed function space for each
-        species"""
+        species."""
 
         sub_solutions = list(ufl.split(self.u))
         sub_prev_solution = list(ufl.split(self.u_n))
@@ -656,7 +653,7 @@ def assign_functions_to_species(self):
             spe.test_function = sub_test_functions[idx]
 
     def define_boundary_conditions(self):
-        """Defines the boundary conditions of the model"""
+        """Defines the boundary conditions of the model."""
 
         for bc in self._unpacked_bcs:
             if isinstance(bc.species, str):
@@ -672,7 +669,7 @@ def define_boundary_conditions(self):
         super().define_boundary_conditions()
 
     def create_dirichletbc_form(self, bc):
-        """Creates a dirichlet boundary condition form
+        """Creates a dirichlet boundary condition form.
 
         Args:
             bc (festim.DirichletBC): the boundary condition
@@ -712,7 +709,7 @@ def create_dirichletbc_form(self, bc):
         return form
 
     def convert_source_input_values_to_fenics_objects(self):
-        """For each source create the value_fenics"""
+        """For each source create the value_fenics."""
         for source in self.sources:
             # create value_fenics for all F.ParticleSource objects
             if isinstance(source, _source.ParticleSource):
@@ -724,7 +721,7 @@ def convert_source_input_values_to_fenics_objects(self):
                 )
 
     def convert_advection_term_to_fenics_objects(self):
-        """For each advection term convert the input value"""
+        """For each advection term convert the input value."""
 
         for advec_term in self.advection_terms:
             advec_term.velocity.convert_input_value(
@@ -732,7 +729,7 @@ def convert_advection_term_to_fenics_objects(self):
             )
 
     def create_flux_values_fenics(self):
-        """For each particle flux create the value_fenics"""
+        """For each particle flux create the value_fenics."""
         for bc in self.boundary_conditions:
             # create value_fenics for all F.ParticleFluxBC objects
             if isinstance(bc, boundary_conditions.ParticleFluxBC):
@@ -743,8 +740,8 @@ def create_flux_values_fenics(self):
                 )
 
     def create_initial_conditions(self):
-        """For each initial condition, create the value_fenics and assign it to
-        the previous solution of the condition's species"""
+        """For each initial condition, create the value_fenics and assign it to the
+        previous solution of the condition's species."""
 
         if len(self.initial_conditions) > 0 and not self.settings.transient:
             raise ValueError(
@@ -773,7 +770,7 @@ def create_initial_conditions(self):
             )
 
     def create_formulation(self):
-        """Creates the formulation of the model"""
+        """Creates the formulation of the model."""
 
         self.formulation = 0
 
@@ -809,7 +806,7 @@ def create_formulation(self):
                             )
                         case _:
                             raise NotImplementedError(
-                                f"Unknown coordinate system {self.mesh.coordinate_system!s}"
+                                f"Unknown coordinate system {self.mesh.coordinate_system!s}"  # noqa: E501
                             )
 
                 if self.settings.transient:
@@ -929,7 +926,7 @@ def update_time_dependent_values(self):
                 advec_term.velocity.update(t=t)
 
     def update_post_processing_solutions(self):
-        """Updates the post-processing solutions of each species"""
+        """Updates the post-processing solutions of each species."""
 
         for spe in self.species:
             spe.post_processing_solution.x.array[:] = self.u.x.array[
@@ -937,7 +934,7 @@ def update_post_processing_solutions(self):
             ]
 
     def post_processing(self):
-        """Post processes the model"""
+        """Post processes the model."""
 
         self.update_post_processing_solutions()
 
@@ -1058,8 +1055,8 @@ def __init__(
         surface_to_volume: dict | None = None,
         petsc_options: dict | None = None,
     ):
-        """Class for a multi-material hydrogen transport problem
-        For other arguments see ``festim.HydrogenTransportProblem``.
+        """Class for a multi-material hydrogen transport problem For other arguments see
+        ``festim.HydrogenTransportProblem``.
 
         Args:
             interfaces (list, optional): list of interfaces (``festim.Interface``
@@ -1218,8 +1215,7 @@ def define_temperature(self):
                 subdomain.sub_T = sub_T
 
     def create_dirichletbc_form(self, bc: boundary_conditions.FixedConcentrationBC):
-        """
-        Creates the ``value_fenics`` attribute for a given
+        """Creates the ``value_fenics`` attribute for a given
         ``festim.FixedConcentrationBC`` and returns the appropriate
         ``dolfinx.fem.DirichletBC`` object.
 
@@ -1270,8 +1266,8 @@ def create_dirichletbc_form(self, bc: boundary_conditions.FixedConcentrationBC):
         return form
 
     def create_initial_conditions(self):
-        """For each intial condition, create the value_fenics and assign it to
-        the previous solution of the condition's species"""
+        """For each intial condition, create the value_fenics and assign it to the
+        previous solution of the condition's species."""
 
         for condition in self.initial_conditions:
             idx = self.species.index(condition.species)
@@ -1299,13 +1295,13 @@ def create_initial_conditions(self):
     def define_function_spaces(
         self, subdomain: _subdomain.VolumeSubdomain, element_degree=1
     ):
-        """
-        Creates appropriate function space and functions for a given subdomain (submesh)
-        based on the number of species existing in this subdomain. Then stores the
-        functionspace, the current solution (``u``) and the previous solution (``u_n``)
-        functions. It also populates the correspondance dicts attributes of the species
-        (eg. ``species.subdomain_to_solution``, ``species.subdomain_to_test_function``,
-        etc) for easy access to the right subfunctions, sub-testfunctions etc.
+        """Creates appropriate function space and functions for a given subdomain
+        (submesh) based on the number of species existing in this subdomain. Then stores
+        the functionspace, the current solution (``u``) and the previous solution
+        (``u_n``) functions. It also populates the correspondance dicts attributes of
+        the species (eg. ``species.subdomain_to_solution``,
+        ``species.subdomain_to_test_function``, etc) for easy access to the right
+        subfunctions, sub-testfunctions etc.
 
         Args:
             subdomain (F.VolumeSubdomain): a subdomain of the geometry
@@ -1358,7 +1354,7 @@ def define_function_spaces(
             species.subdomain_to_post_processing_solution[subdomain].name = name
 
     def convert_source_input_values_to_fenics_objects(self):
-        """For each source create the value_fenics"""
+        """For each source create the value_fenics."""
         for source in self.sources:
             # create value_fenics for all F.ParticleSource objects
             if isinstance(source, _source.ParticleSource):
@@ -1373,7 +1369,7 @@ def convert_source_input_values_to_fenics_objects(self):
                     )
 
     def convert_advection_term_to_fenics_objects(self):
-        """For each advection term convert the input value"""
+        """For each advection term convert the input value."""
 
         for advec_term in self.advection_terms:
             if isinstance(advec_term, AdvectionTerm):
@@ -1389,8 +1385,7 @@ def define_boundary_conditions(self):
         super().define_boundary_conditions()
 
     def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain):
-        """
-        Creates the variational formulation for each subdomain and stores it in
+        """Creates the variational formulation for each subdomain and stores it in
         ``subdomain.F``
 
         Args:
@@ -1433,7 +1428,7 @@ def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain):
                         )
                     case _:
                         raise ValueError(
-                            f"Unsupported coordinate system {self.mesh.coordinate_system}"
+                            f"Unsupported coordinate system {self.mesh.coordinate_system}"  # noqa: E501
                         )
 
         # add reaction terms
@@ -1499,8 +1494,8 @@ def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain):
         subdomain.F = form
 
     def create_formulation(self):
-        """
-        Takes all the formulations for each subdomain and adds the interface conditions.
+        """Takes all the formulations for each subdomain and adds the interface
+        conditions.
 
         Finally compute the jacobian matrix and store it in the ``J`` attribute,
         adds the ``entity_maps`` to the forms and store them in the ``forms`` attribute
@@ -1567,11 +1562,12 @@ def create_formulation(self):
         )
 
     def override_solution_attributes(self, reaction: _reaction.Reaction):
-        """
-        Reaction.reaction_term() relies on the .solution attribute of the species
+        """Reaction.reaction_term() relies on the .solution attribute of the species
         however, in the discontinuous class, this attribute doesn't really make sense
         since there is one solution per subdomain.
-        Therefore we temporarily override the .solution attribute based on the reactants,
+
+        Therefore we temporarily override the .solution attribute based on the
+        reactants,
         products, and `others` if there are implicit species
         """
         list_of_species_to_override = reaction.reactant + reaction.product
@@ -1636,7 +1632,7 @@ def create_solver(self):
                 del opts[f"{prefix}{k}"]
 
     def create_flux_values_fenics(self):
-        """For each particle flux create the ``value_fenics`` attribute"""
+        """For each particle flux create the ``value_fenics`` attribute."""
         for bc in self._unpacked_bcs:
             if isinstance(bc, boundary_conditions.ParticleFluxBC):
                 volume_subdomain = self.surface_to_volume[bc.subdomain]
@@ -1808,7 +1804,7 @@ def post_processing(self):
                 export.t.append(float(self.t))
 
     def iterate(self):
-        """Iterates the model for a given time step"""
+        """Iterates the model for a given time step."""
         if self.show_progress_bar:
             self.progress_bar.update(
                 min(self.dt.value, abs(self.settings.final_time - self.t.value))
@@ -1819,12 +1815,12 @@ def iterate(self):
 
         # Solve main problem
         if Version(dolfinx.__version__) == Version("0.9.0"):
-            nb_its, converged = self.solver.solve()
+            nb_its, _converged = self.solver.solve()
         elif Version(dolfinx.__version__) > Version("0.9.0"):
             _ = self.solver.solve()
             converged_reason = self.solver.solver.getConvergedReason()
             assert converged_reason > 0, (
-                f"Non-linear solver did not converge. Reason code: {converged_reason}. \n See https://petsc.org/release/manualpages/SNES/SNESConvergedReason/ for more information."
+                f"Non-linear solver did not converge. Reason code: {converged_reason}. \n See https://petsc.org/release/manualpages/SNES/SNESConvergedReason/ for more information."  # noqa: E501
             )
             nb_its = self.solver.solver.getIterationNumber()
 
@@ -1888,7 +1884,7 @@ def initialise(self):
         super().initialise()
 
     def create_formulation(self):
-        """Creates the formulation of the model"""
+        """Creates the formulation of the model."""
 
         self.formulation = 0
 
@@ -1963,7 +1959,7 @@ def create_formulation(self):
                         )
 
     def add_reaction_term(self, reaction: _reaction.Reaction):
-        """Adds the reaction term to the formulation"""
+        """Adds the reaction term to the formulation."""
 
         products = (
             reaction.product
@@ -2044,7 +2040,7 @@ def override_post_processing_solution(self):
             )  # NOTE: do we need this line since it's in initialise?
 
     def update_post_processing_solutions(self):
-        """Updates the post-processing solutions after each time step"""
+        """Updates the post-processing solutions after each time step."""
         # need to compute c = theta * K_S
         # this expression is stored in species.dg_expr
 
@@ -2054,7 +2050,7 @@ def update_post_processing_solutions(self):
             spe.post_processing_solution.interpolate(spe.dg_expr)
 
     def create_dirichletbc_form(self, bc: boundary_conditions.FixedConcentrationBC):
-        """Creates a dirichlet boundary condition form
+        """Creates a dirichlet boundary condition form.
 
         Args:
             bc (festim.DirichletBC): the boundary condition
diff --git a/src/festim/initial_condition.py b/src/festim/initial_condition.py
index 7153dffc4..e10768a80 100644
--- a/src/festim/initial_condition.py
+++ b/src/festim/initial_condition.py
@@ -15,8 +15,7 @@
 
 
 class InitialConditionBase:
-    """
-    Base initial condition class
+    """Base initial condition class.
 
     Args:
         value: the value of the initial condition.
@@ -67,8 +66,7 @@ def volume(self, value):
 
 
 class InitialConcentration(InitialConditionBase):
-    """
-    Initial concentration class
+    """Initial concentration class.
 
     Args:
         value: the value of the initial concentration of a given species.
@@ -172,8 +170,7 @@ def create_expr_fenics(
 
 
 class InitialTemperature(InitialConditionBase):
-    """
-    Initial temperature class
+    """Initial temperature class.
 
     Args:
         value: the value of the initial temperature
@@ -247,8 +244,7 @@ def read_function_from_file(
     order: int = 1,
     mesh: dolfinx.mesh.Mesh | None = None,
 ) -> fem.Function:
-    """
-    Read a function from a file
+    """Read a function from a file.
 
     note::
         The function is read from a file using adios4dolfinx. For more information
diff --git a/src/festim/material.py b/src/festim/material.py
index c66a9cb97..1ffff3dc0 100644
--- a/src/festim/material.py
+++ b/src/festim/material.py
@@ -29,8 +29,7 @@ def from_string(cls, s: str):
 
 
 class Material:
-    """
-    Material class
+    """Material class.
 
     Args:
         D_0: the pre-exponential factor of the
@@ -147,7 +146,7 @@ def D(self, value):
             raise TypeError("D must be of type fem.Function")
 
     def get_D_0(self, species=None):
-        """Returns the pre-exponential factor of the diffusion coefficient
+        """Returns the pre-exponential factor of the diffusion coefficient.
 
         Args:
             species (festim.Species or str, optional): the species we want the
@@ -176,7 +175,7 @@ def get_D_0(self, species=None):
             raise TypeError("D_0 must be either a float, int or a dict")
 
     def get_E_D(self, species=None):
-        """Returns the activation energy of the diffusion coefficient
+        """Returns the activation energy of the diffusion coefficient.
 
         Args:
             species (festim.Species or str, optional): the species we want the
@@ -205,7 +204,7 @@ def get_E_D(self, species=None):
             raise TypeError("E_D must be either a float, int or a dict")
 
     def get_K_S_0(self, species=None) -> float:
-        """Returns the pre-exponential factor of the solubility coefficient
+        """Returns the pre-exponential factor of the solubility coefficient.
 
         Args:
             species: the species we want the pre-exponential
@@ -233,7 +232,7 @@ def get_K_S_0(self, species=None) -> float:
             raise TypeError("K_S_0 must be either a float, int or a dict")
 
     def get_E_K_S(self, species=None) -> float:
-        """Returns the activation energy of the solubility coefficient
+        """Returns the activation energy of the solubility coefficient.
 
         Args:
             species: the species we want the activation
@@ -261,7 +260,7 @@ def get_E_K_S(self, species=None) -> float:
             raise TypeError("E_K_S must be either a float, int or a dict")
 
     def get_diffusion_coefficient(self, mesh=None, temperature=None, species=None):
-        """Defines the diffusion coefficient
+        """Defines the diffusion coefficient.
 
         Args:
 
@@ -319,7 +318,7 @@ def get_diffusion_coefficient(self, mesh=None, temperature=None, species=None):
             raise ValueError("D_0 and E_D must be either floats or dicts")
 
     def get_solubility_coefficient(self, mesh, temperature, species=None):
-        """Defines the solubility coefficient
+        """Defines the solubility coefficient.
 
         Args:
 
diff --git a/src/festim/mesh/mesh.py b/src/festim/mesh/mesh.py
index 032897c7a..76fa9003f 100644
--- a/src/festim/mesh/mesh.py
+++ b/src/festim/mesh/mesh.py
@@ -26,7 +26,7 @@ def from_string(cls, s: str):
             return cls.SPHERICAL
         else:
             raise ValueError(
-                "coordinate_system must be one of 'cartesian', 'cylindrical', or 'spherical'"
+                "coordinate_system must be one of 'cartesian', 'cylindrical', or 'spherical'"  # noqa: E501
             )
 
     def __str__(self):
@@ -41,8 +41,7 @@ def __str__(self):
 
 
 class Mesh:
-    """
-    Mesh class
+    """Mesh class.
 
     Args:
         mesh: The mesh. Defaults to None.
@@ -129,11 +128,13 @@ def n(self):
         return ufl.FacetNormal(self._mesh)
 
     def define_meshtags(self, surface_subdomains, volume_subdomains, interfaces=None):
-        """Defines the facet and volume meshtags of the mesh
+        """Defines the facet and volume meshtags of the mesh.
 
         Args:
-            surface_subdomains (list of festim.SufaceSubdomains): the surface subdomains of the model
-            volume_subdomains (list of festim.VolumeSubdomains): the volume subdomains of the model
+            surface_subdomains (list of festim.SufaceSubdomains): the surface
+            subdomains of the model
+            volume_subdomains (list of festim.VolumeSubdomains): the volume
+            subdomains of the model
             interfaces (dict, optional): the interfaces between volume
                 subdomains {int: [VolumeSubdomain, VolumeSubdomain]}. Defaults to None.
 
@@ -180,6 +181,14 @@ def define_meshtags(self, surface_subdomains, volume_subdomains, interfaces=None
                     )
                 tags_volumes[:] = vol.id
 
+        # make sure there are no zeros in the tags_volumes
+        assert np.all(tags_volumes != 0), (
+            "All cells must be tagged with a non-zero value"
+            "This error can be caused by a volume subdomain not being defined"
+            "correctly, or by a mesh that is too coarse to capture the geometry"
+            "of the volume subdomains"
+        )
+
         volume_meshtags = meshtags(
             self._mesh, self.vdim, mesh_cell_indices, tags_volumes
         )
@@ -211,7 +220,7 @@ def define_meshtags(self, surface_subdomains, volume_subdomains, interfaces=None
 
     def check_mesh_dim_coords(self):
         """Checks if the used coordinates can be applied for geometry with the specified
-        dimensions"""
+        dimensions."""
 
         if self.coordinate_system == CoordinateSystem.SPHERICAL and self.vdim != 1:
             raise AttributeError(
diff --git a/src/festim/mesh/mesh_1d.py b/src/festim/mesh/mesh_1d.py
index d07d2c58e..45b4b66a5 100644
--- a/src/festim/mesh/mesh_1d.py
+++ b/src/festim/mesh/mesh_1d.py
@@ -9,8 +9,7 @@
 
 
 class Mesh1D(Mesh):
-    """
-    1D Mesh
+    """1D Mesh.
 
     Args:
         vertices (list or np.ndarray): the mesh x-coordinates (m)
@@ -34,7 +33,7 @@ def vertices(self, value):
         self._vertices = np.sort(np.unique(value)).astype(np.float64)
 
     def generate_mesh(self):
-        """Generates a 1D mesh"""
+        """Generates a 1D mesh."""
 
         if MPI.COMM_WORLD.rank == 0:
             mesh_points = np.reshape(self.vertices, (len(self.vertices), 1))
@@ -54,7 +53,7 @@ def generate_mesh(self):
         )
 
     def check_borders(self, volume_subdomains):
-        """Checks that the borders of the subdomain are within the domain
+        """Checks that the borders of the subdomain are within the domain.
 
         Args:
             volume_subdomains (list of festim.VolumeSubdomain1D): the volume subdomains
diff --git a/src/festim/mesh/mesh_from_xdmf.py b/src/festim/mesh/mesh_from_xdmf.py
index 9d0d20a1d..a318f06bc 100644
--- a/src/festim/mesh/mesh_from_xdmf.py
+++ b/src/festim/mesh/mesh_from_xdmf.py
@@ -6,8 +6,7 @@
 
 
 class MeshFromXDMF(Mesh):
-    """
-    Mesh read from the XDMF files
+    """Mesh read from the XDMF files.
 
     Args:
         volume_file (str): path to the volume file
@@ -23,7 +22,8 @@ class MeshFromXDMF(Mesh):
         volume_file (str): path to the volume file
         facet_file (str): path to the facet file
         mesh_name (str): name of the mesh in the volume XDMF file.
-        surface_meshtags_name (str): name of the surface meshtags in the facet XDMF file.
+        surface_meshtags_name (str): name of the surface meshtags in the facet
+        XDMF file.
         volume_meshtags_name (str): name of the volume meshtags in the volume XDMF file
         mesh (fenics.mesh.Mesh): the fenics mesh
     """
@@ -48,7 +48,7 @@ def __init__(
         super().__init__(mesh=mesh)
 
     def define_surface_meshtags(self):
-        """Creates the facet meshtags
+        """Creates the facet meshtags.
 
         Returns:
             dolfinx.MeshTags: the facet meshtags
@@ -61,7 +61,7 @@ def define_surface_meshtags(self):
         return facet_meshtags
 
     def define_volume_meshtags(self):
-        """Creates the volume meshtags
+        """Creates the volume meshtags.
 
         Returns:
             dolfinx.MeshTags: the volume meshtags
diff --git a/src/festim/problem.py b/src/festim/problem.py
index 15901c797..675ef7474 100644
--- a/src/festim/problem.py
+++ b/src/festim/problem.py
@@ -1,3 +1,4 @@
+import warnings
 from typing import Any
 
 from mpi4py import MPI
@@ -10,7 +11,6 @@
 from dolfinx import fem
 from dolfinx.nls.petsc import NewtonSolver
 from packaging.version import Version
-import warnings
 
 import festim as F
 from festim.mesh.mesh import Mesh as _Mesh
@@ -21,9 +21,10 @@
 
 
 class ProblemBase:
-    """
-    Base class for :py:class:`HeatTransferProblem <festim.heat_transfer_problem.HeatTransferProblem>` and
-    :py:class:`HydrogenTransportProblem <festim.hydrogen_transport_problem.HydrogenTransportProblem>`.
+    """Base class for :py:class:`HeatTransferProblem
+    <festim.heat_transfer_problem.HeatTransferProblem>` and
+    :py:class:`HydrogenTransportProblem
+    <festim.hydrogen_transport_problem.HydrogenTransportProblem>`.
 
     Attributes:
         show_progress_bar: If `True` a progress bar is displayed during the simulation
@@ -85,8 +86,8 @@ def timesteps(self):
         return self._timesteps
 
     def define_meshtags_and_measures(self):
-        """Defines the facet and volume meshtags of the model which are used
-        to define the measures fo the model, dx and ds"""
+        """Defines the facet and volume meshtags of the model which are used to define
+        the measures fo the model, dx and ds."""
 
         if isinstance(self.mesh, F.MeshFromXDMF):
             # TODO: fix naming inconsistency between facet and surface meshtags
@@ -119,16 +120,15 @@ def define_meshtags_and_measures(self):
         )
 
     def define_boundary_conditions(self):
-        """Defines the dirichlet boundary conditions of the model"""
+        """Defines the dirichlet boundary conditions of the model."""
         for bc in self.boundary_conditions:
             if isinstance(bc, F.DirichletBCBase):
                 form = self.create_dirichletbc_form(bc)
                 self.bc_forms.append(form)
 
     def get_petsc_options(self) -> dict[str, Any]:
-        """
-        Gets the PETSc options to pass to the NewtonProblem solver. Default
-        options are updated with user-provided options, if any.
+        """Gets the PETSc options to pass to the NewtonProblem solver. Default options
+        are updated with user-provided options, if any.
 
         Returns:
             the petsc options to pass to the NewtonProblem solver.
@@ -165,7 +165,7 @@ def get_petsc_options(self) -> dict[str, Any]:
         return petsc_options
 
     def create_solver(self):
-        """Creates the solver of the model"""
+        """Creates the solver of the model."""
 
         if Version(dolfinx.__version__) == Version("0.9.0"):
             problem = fem.petsc.NonlinearProblem(
@@ -226,7 +226,7 @@ def create_solver(self):
                 del opts[f"{prefix}{k}"]
 
     def run(self):
-        """Runs the model"""
+        """Runs the model."""
 
         if self.settings.transient:
             # Solve transient
@@ -250,7 +250,7 @@ def run(self):
             self.post_processing()
 
     def iterate(self):
-        """Iterates the model for a given time step"""
+        """Iterates the model for a given time step."""
         self._timesteps.append(float(self.t))
 
         if self.show_progress_bar:
@@ -271,12 +271,12 @@ def iterate(self):
 
         # solve main problem
         if Version(dolfinx.__version__) == Version("0.9.0"):
-            nb_its, converged = self.solver.solve(self.u)
+            nb_its, _converged = self.solver.solve(self.u)
         elif Version(dolfinx.__version__) > Version("0.9.0"):
             _ = self.solver.solve()
             converged_reason = self.solver.solver.getConvergedReason()
             assert converged_reason > 0, (
-                f"Non-linear solver did not converge. Reason code: {converged_reason}. \n See https://petsc.org/release/manualpages/SNES/SNESConvergedReason/ for more information."
+                f"Non-linear solver did not converge. Reason code: {converged_reason}. \n See https://petsc.org/release/manualpages/SNES/SNESConvergedReason/ for more information."  # noqa: E501
             )
             nb_its = self.solver.solver.getIterationNumber()
 
diff --git a/src/festim/reaction.py b/src/festim/reaction.py
index 336b18f80..946bfc7e2 100644
--- a/src/festim/reaction.py
+++ b/src/festim/reaction.py
@@ -13,22 +13,26 @@ class Reaction:
     """A reaction between two species, with a forward and backward rate.
 
     Arguments:
-        reactant (Union[F.Species, F.ImplicitSpecies], List[Union[F.Species, F.ImplicitSpecies]]): The reactant.
+        reactant (Union[F.Species, F.ImplicitSpecies], List[Union[F.Species,
+        F.ImplicitSpecies]]): The reactant.
         product (Optional[Union[F.Species, List[F.Species]]]): The product.
         k_0 (float): The forward rate constant pre-exponential factor.
         E_k (float): The forward rate constant activation energy.
         p_0 (float): The backward rate constant pre-exponential factor.
         E_p (float): The backward rate constant activation energy.
-        volume (F.VolumeSubdomain1D): The volume subdomain where the reaction takes place.
+        volume (F.VolumeSubdomain1D): The volume subdomain where the reaction
+        takes place.
 
     Attributes:
-        reactant (Union[F.Species, F.ImplicitSpecies], List[Union[F.Species, F.ImplicitSpecies]]): The reactant.
+        reactant (Union[F.Species, F.ImplicitSpecies], List[Union[F.Species,
+        F.ImplicitSpecies]]): The reactant.
         product (Optional[Union[F.Species, List[F.Species]]]): The product.
         k_0 (float): The forward rate constant pre-exponential factor.
         E_k (float): The forward rate constant activation energy.
         p_0 (float): The backward rate constant pre-exponential factor.
         E_p (float): The backward rate constant activation energy.
-        volume (F.VolumeSubdomain1D): The volume subdomain where the reaction takes place.
+        volume (F.VolumeSubdomain1D): The volume subdomain where the reaction
+        takes place.
 
     Examples:
 
@@ -53,7 +57,6 @@ class Reaction:
             # compute the reaction term at a given temperature
             temperature = 300.0
             reaction_term = reaction.reaction_term(temperature)
-
     """
 
     def __init__(
@@ -63,8 +66,8 @@ def __init__(
         E_k: float,
         volume: VS1D,
         product: Union[_Species, list[_Species]] | None = [],
-        p_0: float = None,
-        E_p: float = None,
+        p_0: float | None = None,
+        E_p: float | None = None,
     ) -> None:
         self.k_0 = k_0
         self.E_k = E_k
@@ -84,7 +87,7 @@ def reactant(self, value):
             value = [value]
         if len(value) == 0:
             raise ValueError(
-                "reactant must be an entry of one or more species objects, not an empty list."
+                "reactant must be an entry of one or more species objects, not an empty list."  # noqa: E501
             )
         for i in value:
             if not isinstance(i, (_Species, _ImplicitSpecies)):
@@ -101,7 +104,7 @@ def __repr__(self) -> str:
             products = " + ".join([str(product) for product in self.product])
         else:
             products = self.product
-        return f"Reaction({reactants} <--> {products}, {self.k_0}, {self.E_k}, {self.p_0}, {self.E_p})"
+        return f"Reaction({reactants} <--> {products}, {self.k_0}, {self.E_k}, {self.p_0}, {self.E_p})"  # noqa: E501
 
     def __str__(self) -> str:
         reactants = " + ".join([str(reactant) for reactant in self.reactant])
@@ -114,8 +117,8 @@ def __str__(self) -> str:
     def reaction_term(
         self,
         temperature,
-        reactant_concentrations: list = None,
-        product_concentrations: list = None,
+        reactant_concentrations: list | None = None,
+        product_concentrations: list | None = None,
     ) -> Expr:
         """Compute the reaction term at a given temperature.
 
@@ -146,11 +149,11 @@ def reaction_term(
                     + " when no products are present."
                 )
         else:
-            if self.p_0 == None:
+            if self.p_0 is None:
                 raise ValueError(
                     "p_0 cannot be None when reaction products are present."
                 )
-            elif self.E_p == None:
+            elif self.E_p is None:
                 raise ValueError(
                     "E_p cannot be None when reaction products are present."
                 )
diff --git a/src/festim/settings.py b/src/festim/settings.py
index 16fdbdad6..c3652211b 100644
--- a/src/festim/settings.py
+++ b/src/festim/settings.py
@@ -30,7 +30,6 @@ class Settings:
         stepsize (festim.Stepsize): stepsize for a transient
             simulation.
         convergence_criterion: resiudal or incremental (for Newton solver)
-
     """
 
     def __init__(
diff --git a/src/festim/source.py b/src/festim/source.py
index 054bbde1e..9d8dc23f6 100644
--- a/src/festim/source.py
+++ b/src/festim/source.py
@@ -8,8 +8,7 @@
 
 
 class SourceBase:
-    """
-    Source base class
+    """Source base class.
 
     Args:
         value: the value of the source
@@ -72,8 +71,7 @@ def volume(self, value):
 
 
 class ParticleSource(SourceBase):
-    """
-    Particle source class
+    """Particle source class.
 
     Args:
         value: the value of the source
@@ -119,8 +117,7 @@ def species(self, value):
 
 
 class HeatSource(SourceBase):
-    """
-    Heat source class
+    """Heat source class.
 
     Args:
         value: the value of the source
diff --git a/src/festim/species.py b/src/festim/species.py
index 9b9aa412b..8f8d8babb 100644
--- a/src/festim/species.py
+++ b/src/festim/species.py
@@ -10,8 +10,7 @@
 
 
 class Species:
-    """
-    Hydrogen species class for H transport simulation.
+    """Hydrogen species class for H transport simulation.
 
     Args:
         name: a name given to the species. Defaults to None.
@@ -53,8 +52,6 @@ class Species:
 
             Species(name="H")
             Species(name="Trap", mobile=False)
-
-
     """
 
     name: str | None
@@ -113,9 +110,7 @@ def concentration(self):
 
     @property
     def legacy(self) -> bool:
-        """
-        Check if we are using FESTIM 1.0 implementation or FESTIM 2.0
-        """
+        """Check if we are using FESTIM 1.0 implementation or FESTIM 2.0."""
         if not self.subdomain_to_solution:
             return True
         else:
@@ -176,10 +171,9 @@ def concentration(self):
 
     def create_value_fenics(self, mesh, t: fem.Constant):
         """Creates the value of the density as a fenics object and sets it to
-        self.value_fenics.
-        If the value is a constant, it is converted to a fenics.Constant.
-        If the value is a function of t, it is converted to a fenics.Constant.
-        Otherwise, it is converted to a ufl Expression
+        self.value_fenics. If the value is a constant, it is converted to a
+        fenics.Constant. If the value is a function of t, it is converted to a
+        fenics.Constant. Otherwise, it is converted to a ufl Expression.
 
         Args:
             mesh (dolfinx.mesh.Mesh) : the mesh
@@ -231,8 +225,7 @@ def update_density(self, t):
 
 
 def find_species_from_name(name: str, species: list):
-    """Returns the correct species object from a list of species
-    based on a string
+    """Returns the correct species object from a list of species based on a string.
 
     Args:
         name (str): the name of the species
@@ -243,7 +236,6 @@ def find_species_from_name(name: str, species: list):
 
     Raises:
         ValueError: if the species name is not found in the list of species
-
     """
     for spe in species:
         if spe.name == name:
diff --git a/src/festim/stepsize.py b/src/festim/stepsize.py
index e555e6594..e8f92512c 100644
--- a/src/festim/stepsize.py
+++ b/src/festim/stepsize.py
@@ -2,8 +2,7 @@
 
 
 class Stepsize:
-    """
-    A class for evaluating the stepsize of transient simulations.
+    """A class for evaluating the stepsize of transient simulations.
 
     Args:
         initial_value (float, int): initial stepsize.
@@ -113,8 +112,7 @@ def max_stepsize(self, value):
         self._max_stepsize = value
 
     def get_max_stepsize(self, t):
-        """
-        Returns the maximum stepsize at time t.
+        """Returns the maximum stepsize at time t.
 
         Args:
             t (float): the current time
@@ -152,9 +150,7 @@ def modify_value(self, value, nb_iterations, t=None):
         return updated_value
 
     def is_adapt(self, t):
-        """
-        Methods that defines if the stepsize should be
-        adapted or not
+        """Methods that defines if the stepsize should be adapted or not.
 
         Args:
             t (float): the current time
@@ -165,8 +161,8 @@ def is_adapt(self, t):
         return True
 
     def next_milestone(self, current_time: float):
-        """Returns the next milestone that the simulation must pass.
-        Returns None if there are no more milestones.
+        """Returns the next milestone that the simulation must pass. Returns None if
+        there are no more milestones.
 
         Args:
             current_time (float): current time.
diff --git a/src/festim/subdomain/interface.py b/src/festim/subdomain/interface.py
index 02060d31c..8b568d40e 100644
--- a/src/festim/subdomain/interface.py
+++ b/src/festim/subdomain/interface.py
@@ -54,9 +54,8 @@ def from_string(cls, s: str) -> "InterfaceMethod":
 class InterfaceBase(ABC):
     """Abstract base class for interfaces between subdomains.
 
-    Provides common functionality for handling interfaces in discontinuous
-    finite element problems, including integration data computation and
-    restriction handling.
+    Provides common functionality for handling interfaces in discontinuous finite
+    element problems, including integration data computation and restriction handling.
     """
 
     def __init__(
@@ -76,9 +75,9 @@ def __init__(
     def pad_parent_maps(self):
         """Pad parent-to-submesh maps for correct sparsity pattern.
 
-        This is a workaround to ensure the sparsity pattern is correct when
-        assembling forms with interface integrals. It pads the mapping between
-        parent mesh cells and submesh cells for DOLFINx versions that require it.
+        This is a workaround to ensure the sparsity pattern is correct when assembling
+        forms with interface integrals. It pads the mapping between parent mesh cells
+        and submesh cells for DOLFINx versions that require it.
         """
         try:
             # No padding needed for latest version of DOLFINx
diff --git a/src/festim/subdomain/surface_subdomain.py b/src/festim/subdomain/surface_subdomain.py
index bdf673374..5a7e2803d 100644
--- a/src/festim/subdomain/surface_subdomain.py
+++ b/src/festim/subdomain/surface_subdomain.py
@@ -5,8 +5,7 @@
 
 
 class SurfaceSubdomain:
-    """
-    Surface subdomain class
+    """Surface subdomain class.
 
     Args:
         id: the id of the surface subdomain
@@ -21,9 +20,12 @@ class SurfaceSubdomain:
         .. testcode:: SurfaceSubdomain
 
             SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[0], 1.0))
-            SurfaceSubdomain(id=1, locator=lambda x: np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1.0)))
-            SurfaceSubdomain(id=1, locator=lambda x: np.logical_and(np.isclose(x[0], 0.0), np.isclose(x[1], 1.0)))
-            SurfaceSubdomain(id=1, locator=lambda x: np.logical_and(np.isclose(x[0], 0.0), x[1] <= 0.5))
+            SurfaceSubdomain(id=1, locator=lambda x:
+            np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1.0)))
+            SurfaceSubdomain(id=1, locator=lambda x:
+            np.logical_and(np.isclose(x[0], 0.0), np.isclose(x[1], 1.0)))
+            SurfaceSubdomain(id=1, locator=lambda x:
+            np.logical_and(np.isclose(x[0], 0.0), x[1] <= 0.5))
     """
 
     id: int
@@ -56,8 +58,7 @@ def locate_boundary_facet_indices(self, mesh: dolfinx.mesh.Mesh) -> np.ndarray:
 
 
 class SurfaceSubdomain1D(SurfaceSubdomain):
-    """
-    Surface subdomain class for 1D cases
+    """Surface subdomain class for 1D cases.
 
     Args:
         id: the id of the surface subdomain
@@ -88,8 +89,8 @@ def __init__(self, id: int, x: float) -> None:
 
 
 def find_surface_from_id(id: int, surfaces: list):
-    """Returns the correct surface subdomain object from a list of surface ids
-    based on an int
+    """Returns the correct surface subdomain object from a list of surface ids based on
+    an int.
 
     Args:
         id (int): the id of the surface subdomain
@@ -100,7 +101,6 @@ def find_surface_from_id(id: int, surfaces: list):
 
     Raises:
         ValueError: if the surface name is not found in the list of surfaces
-
     """
     for surf in surfaces:
         if surf.id == id:
diff --git a/src/festim/subdomain/volume_subdomain.py b/src/festim/subdomain/volume_subdomain.py
index 56bdbd806..3815b5ebc 100644
--- a/src/festim/subdomain/volume_subdomain.py
+++ b/src/festim/subdomain/volume_subdomain.py
@@ -20,11 +20,10 @@
 
 
 class VolumeSubdomain:
-    """
-    Volume subdomain class
+    """Volume subdomain class.
 
     Args:
-        id: the id of the volume subdomain
+        id: the id of the volume subdomain (> 0)
         submesh: the submesh of the volume subdomain
         cell_map: the cell map of the volume subdomain
         parent_mesh: the parent mesh of the volume subdomain
@@ -58,6 +57,7 @@ class VolumeSubdomain:
     def __init__(
         self, id, material, locator: Callable | None = None, name: str | None = None
     ):
+        assert id != 0, "Volume subdomain id cannot be 0"
         self.id = id
         self.material = material
         self.locator = locator
@@ -115,7 +115,7 @@ def transfer_meshtag(self, mesh: dolfinx.mesh.Mesh, tag: dolfinx.mesh.MeshTags):
         )
 
     def locate_subdomain_entities(self, mesh: Mesh) -> npt.NDArray[np.int32]:
-        """Locates all cells in subdomain borders within domain
+        """Locates all cells in subdomain borders within domain.
 
         Args:
             mesh: the mesh of the model
@@ -131,8 +131,7 @@ def locate_subdomain_entities(self, mesh: Mesh) -> npt.NDArray[np.int32]:
 
 
 class VolumeSubdomain1D(VolumeSubdomain):
-    """
-    Volume subdomain class for 1D cases
+    """Volume subdomain class for 1D cases.
 
     Args:
         id (int): the id of the volume subdomain
@@ -166,8 +165,8 @@ def __init__(self, id, borders, material) -> None:
 
 
 def find_volume_from_id(id: int, volumes: list):
-    """Returns the correct volume subdomain object from a list of volume ids
-    based on an int
+    """Returns the correct volume subdomain object from a list of volume ids based on an
+    int.
 
     Args:
         id (int): the id of the volume subdomain
@@ -178,7 +177,6 @@ def find_volume_from_id(id: int, volumes: list):
 
     Raises:
         ValueError: if the volume name is not found in the list of volumes
-
     """
     for vol in volumes:
         if vol.id == id:
diff --git a/src/festim/trap.py b/src/festim/trap.py
index 0e5553a32..653c06fcc 100644
--- a/src/festim/trap.py
+++ b/src/festim/trap.py
@@ -28,13 +28,15 @@ class Trap(_Species):
         volume (F.VolumeSubdomain1D): The volume subdomain where the trap is.
         trapped_concentration (_Species): The immobile trapped concentration
         trap_reaction (_Reaction): The reaction for trapping the mobile conc.
-        empty_trap_sites (F.ImplicitSpecies): The implicit species for the empty trap sites
+        empty_trap_sites (F.ImplicitSpecies): The implicit species for the
+        empty trap sites
 
     Examples:
 
         .. testsetup:: Trap
 
-            from festim import Trap, Species, VolumeSubdomain, Material, HydrogenTransportProblem
+            from festim import Trap, Species, VolumeSubdomain, Material
+            from fsetim import HydrogenTransportProblem
 
             my_mat = Material(D_0=1, E_D=1, name="test_mat")
             my_vol = VolumeSubdomain(id=1, material=my_mat)
@@ -42,7 +44,8 @@ class Trap(_Species):
 
         .. testcode:: Trap
 
-            trap = Trap(name="Trap", mobile_species=H, k_0=1.0, E_k=0.2, p_0=0.1, E_p=0.3, n=100, volume=my_vol)
+            trap = Trap(name="Trap", mobile_species=H, k_0=1.0, E_k=0.2,
+            p_0=0.1, E_p=0.3, n=100, volume=my_vol)
 
             my_model = HydrogenTransportProblem()
             my_model.traps = [trap]
@@ -85,8 +88,6 @@ class Trap(_Species):
             )
             my_model.species = [cm, ct]
             my_model.reactions = [trap_reaction]
-
-
     """
 
     def __init__(
@@ -105,7 +106,7 @@ def __init__(
         self.reaction = None
 
     def create_species_and_reaction(self):
-        """create the immobile trapped species object and the reaction for trapping"""
+        """Create the immobile trapped species object and the reaction for trapping."""
         self.trapped_concentration = _Species(name=self.name, mobile=False)
         self.empty_trap_sites = _ImplicitSpecies(
             n=self.n, others=[self.trapped_concentration]
diff --git a/test/benchmark.py b/test/benchmark.py
index 64dfd99f2..387907484 100644
--- a/test/benchmark.py
+++ b/test/benchmark.py
@@ -187,8 +187,8 @@ def siverts_law(T, S_0, E_S, pressure):
 
 
 def test_festim_vs_fenics_permeation_benchmark():
-    """Runs a problem with pure fenicsx and the same problem with FESTIM and
-    raise ValueError if difference is too high"""
+    """Runs a problem with pure fenicsx and the same problem with FESTIM and raise
+    ValueError if difference is too high."""
 
     start = time.time()
     fenics_test_permeation_problem(mesh_size=20001)
diff --git a/test/system_tests/test_1_mobile_1_trap.py b/test/system_tests/test_1_mobile_1_trap.py
index 39d0e43ec..f08c0152d 100644
--- a/test/system_tests/test_1_mobile_1_trap.py
+++ b/test/system_tests/test_1_mobile_1_trap.py
@@ -18,9 +18,7 @@
 
 
 def test_1_mobile_1_trap_MMS_steady_state():
-    """
-    MMS test with 1 mobile species and 1 trap at steady state
-    """
+    """MMS test with 1 mobile species and 1 trap at steady state."""
 
     def u_exact(mod):
         return lambda x: 1.5 + mod.sin(3 * mod.pi * x[0])
@@ -42,7 +40,10 @@ def v_exact(mod):
     E_k = 1.5
     p_0 = 0.2
     E_p = 0.1
-    T_expr = lambda x: 500 + 100 * x[0]
+
+    def T_expr(x):
+        return 500 + 100 * x[0]
+
     T.interpolate(T_expr)
     n_trap = 3
     E_D = 0.1
@@ -117,10 +118,8 @@ def v_exact(mod):
 
 
 def test_1_mobile_1_trap_MMS_transient():
-    """
-    MMS test with 1 mobile species and 1 trap in 0.1s transient, the value at the last time step is
-    compared to an analytical solution
-    """
+    """MMS test with 1 mobile species and 1 trap in 0.1s transient, the value at the
+    last time step is compared to an analytical solution."""
 
     final_time = 0.1
 
@@ -138,7 +137,10 @@ def u_exact_alt(mod):
 
     D_0 = 1
     E_D = 0.1
-    T_expr = lambda x: 600 + 50 * x[0]
+
+    def T_expr(x):
+        return 600 + 50 * x[0]
+
     T.interpolate(T_expr)
     D = D_0 * ufl.exp(-E_D / (F.k_B * T))
 
@@ -163,12 +165,16 @@ def u_exact_alt(mod):
         F.DirichletBC(subdomain=right, value=H_analytical_ufl, species=H),
     ]
 
-    init_value = lambda x: 1 + ufl.sin(2 * ufl.pi * x[0])
+    def init_value(x):
+        return 1 + ufl.sin(2 * ufl.pi * x[0])
+
     my_model.initial_conditions = [
         F.InitialConcentration(value=init_value, species=H, volume=vol)
     ]
 
-    f = lambda x, t: 4 * t - ufl.div(D * ufl.grad(H_analytical_ufl(x_1d, t)))
+    def f(x, t):
+        return 4 * t - ufl.div(D * ufl.grad(H_analytical_ufl(x_1d, t)))
+
     my_model.sources = [F.ParticleSource(value=f, volume=vol, species=H)]
 
     my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=final_time)
@@ -185,8 +191,8 @@ def u_exact_alt(mod):
 
 
 def test_1_mobile_1_trap_MMS_3D():
-    """Tests that a steady simulation can be run in a 3D domain with
-    1 mobile and 1 trapped species"""
+    """Tests that a steady simulation can be run in a 3D domain with 1 mobile and 1
+    trapped species."""
 
     def u_exact(mod):
         return lambda x: 1.5 + mod.sin(3 * mod.pi * x[0]) + mod.cos(3 * mod.pi * x[1])
@@ -208,7 +214,10 @@ def v_exact(mod):
     E_k = 1.5
     p_0 = 0.2
     E_p = 0.1
-    T_expr = lambda x: 500 + 100 * x[0]
+
+    def T_expr(x):
+        return 500 + 100 * x[0]
+
     T.interpolate(T_expr)
     n_trap = 3
     E_D = 0.1
@@ -298,8 +307,8 @@ def locate_boundary_facet_indices(self, mesh):
 
 
 def test_1_mobile_1_trap_MMS_2D():
-    """Tests that a steady simulation can be run in a 2D domain with
-    1 mobile and 1 trapped species"""
+    """Tests that a steady simulation can be run in a 2D domain with 1 mobile and 1
+    trapped species."""
 
     def u_exact(mod):
         return lambda x: 1.5 + mod.sin(3 * mod.pi * x[0]) + mod.cos(3 * mod.pi * x[1])
@@ -321,7 +330,10 @@ def v_exact(mod):
     E_k = 1.5
     p_0 = 0.2
     E_p = 0.1
-    T_expr = lambda x: 500 + 100 * x[0]
+
+    def T_expr(x):
+        return 500 + 100 * x[0]
+
     T.interpolate(T_expr)
     n_trap = 3
     E_D = 0.1
diff --git a/test/system_tests/test_1_mobile_species.py b/test/system_tests/test_1_mobile_species.py
index 9168a0e3b..b75f3334e 100644
--- a/test/system_tests/test_1_mobile_species.py
+++ b/test/system_tests/test_1_mobile_species.py
@@ -18,9 +18,7 @@
 
 
 def test_1_mobile_MMS_steady_state():
-    """
-    MMS test with one mobile species at steady state
-    """
+    """MMS test with one mobile species at steady state."""
 
     def u_exact(mod):
         return lambda x: 1 + mod.sin(2 * mod.pi * x[0])
@@ -33,7 +31,10 @@ def u_exact(mod):
 
     D_0 = 1
     E_D = 0.1
-    T_expr = lambda x: 500 + 100 * x[0]
+
+    def T_expr(x):
+        return 500 + 100 * x[0]
+
     T.interpolate(T_expr)
     D = D_0 * ufl.exp(-E_D / (F.k_B * T))
 
@@ -72,10 +73,8 @@ def u_exact(mod):
 
 
 def test_1_mobile_MMS_transient():
-    """
-    MMS test with 1 mobile species in 0.1s transient, the value at the last time step is
-    compared to an analytical solution
-    """
+    """MMS test with 1 mobile species in 0.1s transient, the value at the last time step
+    is compared to an analytical solution."""
 
     final_time = 0.1
 
@@ -93,7 +92,10 @@ def u_exact_alt(mod):
 
     D_0 = 1
     E_D = 0.1
-    T_expr = lambda x: 600 + 50 * x[0]
+
+    def T_expr(x):
+        return 600 + 50 * x[0]
+
     T.interpolate(T_expr)
     D = D_0 * ufl.exp(-E_D / (F.k_B * T))
 
@@ -118,12 +120,16 @@ def u_exact_alt(mod):
         F.DirichletBC(subdomain=right, value=H_analytical_ufl, species=H),
     ]
 
-    init_value = lambda x: 1 + ufl.sin(2 * ufl.pi * x[0])
+    def init_value(x):
+        return 1 + ufl.sin(2 * ufl.pi * x[0])
+
     my_model.initial_conditions = [
         F.InitialConcentration(value=init_value, species=H, volume=vol)
     ]
 
-    f = lambda x, t: 4 * t - ufl.div(D * ufl.grad(H_analytical_ufl(x, t)))
+    def f(x, t):
+        return 4 * t - ufl.div(D * ufl.grad(H_analytical_ufl(x, t)))
+
     my_model.sources = [F.ParticleSource(value=f, volume=vol, species=H)]
 
     my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=final_time)
@@ -140,8 +146,8 @@ def u_exact_alt(mod):
 
 
 def test_1_mobile_MMS_2D():
-    """Tests that a steady simulation can be run in a 2D domain with
-    1 mobile species"""
+    """Tests that a steady simulation can be run in a 2D domain with 1 mobile
+    species."""
 
     def u_exact(mod):
         return lambda x: 1 + mod.sin(2 * mod.pi * x[0]) + mod.cos(2 * mod.pi * x[1])
@@ -154,7 +160,10 @@ def u_exact(mod):
 
     D_0 = 1
     E_D = 0.1
-    T_expr = lambda x: 500 + 100 * x[0]
+
+    def T_expr(x):
+        return 500 + 100 * x[0]
+
     T.interpolate(T_expr)
     D = D_0 * ufl.exp(-E_D / (F.k_B * T))
 
@@ -205,8 +214,8 @@ def locate_boundary_facet_indices(self, mesh):
 
 
 def test_1_mobile_MMS_3D():
-    """Tests that a steady simulation can be run in a 3D domain with
-    1 mobile species"""
+    """Tests that a steady simulation can be run in a 3D domain with 1 mobile
+    species."""
 
     def u_exact(mod):
         return lambda x: 1 + mod.sin(2 * mod.pi * x[0]) + mod.cos(2 * mod.pi * x[1])
@@ -219,7 +228,10 @@ def u_exact(mod):
 
     D_0 = 1
     E_D = 0.1
-    T_expr = lambda x: 500 + 100 * x[0]
+
+    def T_expr(x):
+        return 500 + 100 * x[0]
+
     T.interpolate(T_expr)
     D = D_0 * ufl.exp(-E_D / (F.k_B * T))
 
diff --git a/test/system_tests/test_advection_term_integration.py b/test/system_tests/test_advection_term_integration.py
index 3d3e728a0..0eec01693 100644
--- a/test/system_tests/test_advection_term_integration.py
+++ b/test/system_tests/test_advection_term_integration.py
@@ -15,7 +15,7 @@
 def test_MMS_1_species_1_trap_with_advection():
     """MMS coupled heat and hydrogen test with 1 mobile species and 1 trap in a 1s
     transient, the values of the temperature, mobile and trapped solutions at the last
-    time step is compared to an analytical solution"""
+    time step is compared to an analytical solution."""
 
     test_mesh_2d = create_unit_square(MPI.COMM_WORLD, 200, 200)
     x_2d = ufl.SpatialCoordinate(test_mesh_2d)
@@ -38,7 +38,10 @@ def test_MMS_1_species_1_trap_with_advection():
 
     V = fem.functionspace(test_mesh_2d, ("Lagrange", 1))
     T = fem.Function(V)
-    T_expr = lambda x: 100 + 200 * x[0] + 100 * x[1]
+
+    def T_expr(x):
+        return 100 + 200 * x[0] + 100 * x[1]
+
     T.interpolate(T_expr)
 
     # create velocity field
@@ -63,8 +66,11 @@ def velocity_func(x):
     u.interpolate(velocity_func)
 
     # define hydrogen problem
-    exact_mobile_solution = lambda x: 200 * x[0] ** 2 + 300 * x[1] ** 2
-    exact_trapped_solution = lambda x: 10 * x[0] ** 2 + 10 * x[1] ** 2
+    def exact_mobile_solution(x):
+        return 200 * x[0] ** 2 + 300 * x[1] ** 2
+
+    def exact_trapped_solution(x):
+        return 10 * x[0] ** 2 + 10 * x[1] ** 2
 
     D = D_0 * ufl.exp(-E_D / (k_B * T))
     k = k_0 * ufl.exp(-E_k / (k_B * T))
@@ -207,8 +213,8 @@ def velocity_func(x):
             values[1] = 0  # Second component remains zero
             return values
 
-        mesh2, mt2, ct2 = generate_mesh(n=10)
-        submesh, cell_map, v_map = dolfinx.mesh.create_submesh(
+        mesh2, _mt2, ct2 = generate_mesh(n=10)
+        submesh, _cell_map, _v_map = dolfinx.mesh.create_submesh(
             mesh2, ct2.dim, ct2.find(4)
         )[0:3]
         v_cg = basix.ufl.element(
diff --git a/test/system_tests/test_coupled_heat_hydrogen.py b/test/system_tests/test_coupled_heat_hydrogen.py
index df6b657df..f9865e692 100644
--- a/test/system_tests/test_coupled_heat_hydrogen.py
+++ b/test/system_tests/test_coupled_heat_hydrogen.py
@@ -9,7 +9,7 @@
 def test_MMS_coupled_problem():
     """MMS coupled heat and hydrogen test with 1 mobile species and 1 trap in a 1s
     transient, the values of the temperature, mobile and trapped solutions at the last
-    time step is compared to an analytical solution"""
+    time step is compared to an analytical solution."""
 
     # coupled simulation properties
     density, heat_capacity = 1.2, 2.6
@@ -38,7 +38,8 @@ def test_MMS_coupled_problem():
     test_traps = F.ImplicitSpecies(n=n_trap, others=[test_mobile, test_trapped])
 
     # define temperature sim
-    exact_T_solution = lambda x, t: 3 * x[0] ** 2 + 10 * t
+    def exact_T_solution(x, t):
+        return 3 * x[0] ** 2 + 10 * t
 
     dTdt = 10
     mms_T_source = (
@@ -67,18 +68,29 @@ def test_MMS_coupled_problem():
     )
 
     # define hydrogen problem
-    exact_mobile_solution = lambda x, t: 2 * x[0] ** 2 + 15 * t
-    exact_trapped_solution = lambda x, t: 4 * x[0] ** 2 + 12 * t
+    def exact_mobile_solution(x, t):
+        return 2 * x[0] ** 2 + 15 * t
 
-    exact_mobile_intial_cond = lambda x: exact_mobile_solution(x, t=0)
-    exact_trapped_intial_cond = lambda x: exact_trapped_solution(x, t=0)
+    def exact_trapped_solution(x, t):
+        return 4 * x[0] ** 2 + 12 * t
+
+    def exact_mobile_intial_cond(x):
+        return exact_mobile_solution(x, t=0)
+
+    def exact_trapped_intial_cond(x):
+        return exact_trapped_solution(x, t=0)
 
     dmobiledt = 15
     dtrappeddt = 12
 
-    D = lambda x, t: D_0 * ufl.exp(-E_D / (k_B * exact_T_solution(x, t)))
-    k = lambda x, t: k_0 * ufl.exp(-E_k / (k_B * exact_T_solution(x, t)))
-    p = lambda x, t: p_0 * ufl.exp(-E_p / (k_B * exact_T_solution(x, t)))
+    def D(x, t):
+        return D_0 * ufl.exp(-E_D / (k_B * exact_T_solution(x, t)))
+
+    def k(x, t):
+        return k_0 * ufl.exp(-E_k / (k_B * exact_T_solution(x, t)))
+
+    def p(x, t):
+        return p_0 * ufl.exp(-E_p / (k_B * exact_T_solution(x, t)))
 
     def mms_mobile_source(x, t):
         return (
@@ -176,9 +188,14 @@ def mms_trapped_source(x, t):
     mobile_computed = test_mobile.post_processing_solution
     trapped_computed = test_trapped.post_processing_solution
 
-    exact_final_T = lambda x: exact_T_solution(x, t=final_time)
-    exact_final_mobile = lambda x: exact_mobile_solution(x, t=final_time)
-    exact_final_trapped = lambda x: exact_trapped_solution(x, t=final_time)
+    def exact_final_T(x):
+        return exact_T_solution(x, t=final_time)
+
+    def exact_final_mobile(x):
+        return exact_mobile_solution(x, t=final_time)
+
+    def exact_final_trapped(x):
+        return exact_trapped_solution(x, t=final_time)
 
     L2_error_T = error_L2(T_computed, exact_final_T)
     L2_error_mobile = error_L2(mobile_computed, exact_final_mobile)
@@ -193,7 +210,7 @@ def mms_trapped_source(x, t):
 def test_coupled_problem_non_matching_mesh():
     """MMS coupled heat and hydrogen test with 1 mobile species in a 1s transient with
     mismatched meshes, the value of the mobile solution at the last time step is
-    compared to an analytical solution"""
+    compared to an analytical solution."""
 
     # coupled simulation properties
     density, heat_capacity = 1.3, 2.3
@@ -216,7 +233,8 @@ def test_coupled_problem_non_matching_mesh():
     test_mobile = F.Species("mobile", mobile=True)
 
     # define temperature sim
-    exact_T_solution = lambda x, t: 2 * x[0] ** 2 + 5 * t
+    def exact_T_solution(x, t):
+        return 2 * x[0] ** 2 + 5 * t
 
     dTdt = 5
 
@@ -249,11 +267,13 @@ def test_coupled_problem_non_matching_mesh():
     )
 
     # define hydrogen problem
-    exact_mobile_solution = lambda x, t: 4 * x[0] ** 2 + 10 * t
+    def exact_mobile_solution(x, t):
+        return 4 * x[0] ** 2 + 10 * t
 
     dmobiledt = 10
 
-    D = lambda x, t: D_0 * ufl.exp(-E_D / (k_B * exact_T_solution(x, t)))
+    def D(x, t):
+        return D_0 * ufl.exp(-E_D / (k_B * exact_T_solution(x, t)))
 
     def mms_mobile_source(x, t):
         return dmobiledt - ufl.div(D(x, t) * ufl.grad(exact_mobile_solution(x, t)))
@@ -302,7 +322,8 @@ def mms_mobile_source(x, t):
     test_coupled_problem.initialise()
     test_coupled_problem.run()
 
-    exact_solution = lambda x: exact_mobile_solution(x, t=final_time)
+    def exact_solution(x):
+        return exact_mobile_solution(x, t=final_time)
 
     computed_solution = test_mobile.post_processing_solution
 
diff --git a/test/system_tests/test_cylindrical_coordinates.py b/test/system_tests/test_cylindrical_coordinates.py
index eaa313e98..0335c8715 100644
--- a/test/system_tests/test_cylindrical_coordinates.py
+++ b/test/system_tests/test_cylindrical_coordinates.py
@@ -6,14 +6,13 @@
 
 
 def test_run_MMS_cylindrical():
-    """
-    Tests that festim produces the correct concentration field in cylindrical
-    coordinates
-    """
+    """Tests that festim produces the correct concentration field in cylindrical
+    coordinates."""
 
     my_mesh = F.Mesh1D(vertices=np.linspace(1, 2, 500), coordinate_system="cylindrical")
 
-    u_exact = lambda x: 1 + x[0] ** 2
+    def u_exact(x):
+        return 1 + x[0] ** 2
 
     f = -4
 
@@ -67,10 +66,8 @@ def test_run_MMS_cylindrical():
 
 
 def test_run_MMS_cylindrical_mixed_domain():
-    """
-    Tests that festim produces the correct concentration field in cylindrical
-    coordinates in a discontinuous domain with two materials
-    """
+    """Tests that festim produces the correct concentration field in cylindrical
+    coordinates in a discontinuous domain with two materials."""
 
     my_model = F.HydrogenTransportProblemDiscontinuous()
 
@@ -98,7 +95,8 @@ def c_exact_left(x):
     def c_exact_right(x):
         return K_S_right / K_S_left * c_exact_left(x)
 
-    lap_c = lambda r: 4 - 2 * r_interface / r
+    def lap_c(r):
+        return 4 - 2 * r_interface / r
 
     mat_1 = F.Material(D_0=D, E_D=0, K_S_0=K_S_left, E_K_S=0, solubility_law="sievert")
     mat_2 = F.Material(D_0=D, E_D=0, K_S_0=K_S_right, E_K_S=0, solubility_law="sievert")
@@ -130,8 +128,11 @@ def c_exact_right(x):
 
     my_model.temperature = 500
 
-    f_left = lambda x: -D * lap_c(x[0])
-    f_right = lambda x: -D * K_S_right / K_S_left * lap_c(x[0])
+    def f_left(x):
+        return -D * lap_c(x[0])
+
+    def f_right(x):
+        return -D * K_S_right / K_S_left * lap_c(x[0])
 
     my_model.sources = [
         F.ParticleSource(value=f_left, volume=vol_1, species=H),
diff --git a/test/system_tests/test_fluxes.py b/test/system_tests/test_fluxes.py
index 1d3237a86..70304488b 100644
--- a/test/system_tests/test_fluxes.py
+++ b/test/system_tests/test_fluxes.py
@@ -11,11 +11,10 @@
 
 
 def test_flux_bc_1_mobile_MMS_steady_state():
-    """
-    MMS test with a flux BC considering one mobile species at steady state
-    """
+    """MMS test with a flux BC considering one mobile species at steady state."""
 
-    u_exact = lambda x: 1 + 2 * x[0] ** 2
+    def u_exact(x):
+        return 1 + 2 * x[0] ** 2
 
     V = fem.functionspace(test_mesh_1d.mesh, ("Lagrange", 1))
     T = fem.Function(V)
@@ -60,11 +59,10 @@ def test_flux_bc_1_mobile_MMS_steady_state():
 
 
 def test_flux_bc_heat_transfer_steady_state():
-    """
-    MMS test with a flux BC in a heat transfer problem at steady state
-    """
+    """MMS test with a flux BC in a heat transfer problem at steady state."""
 
-    u_exact = lambda x: 1 + 2 * x[0] ** 2
+    def u_exact(x):
+        return 1 + 2 * x[0] ** 2
 
     thermal_cond = 2.3
 
diff --git a/test/system_tests/test_heat_transfer.py b/test/system_tests/test_heat_transfer.py
index ded5284f1..96e41fcd6 100644
--- a/test/system_tests/test_heat_transfer.py
+++ b/test/system_tests/test_heat_transfer.py
@@ -11,9 +11,7 @@
 
 
 def test_heat_transfer_steady_state():
-    """
-    MMS test with heat transfer simulation
-    """
+    """MMS test with heat transfer simulation."""
 
     def u_exact(mod):
         return lambda x: 1 + mod.sin(2 * mod.pi * x[0])
diff --git a/test/system_tests/test_io.py b/test/system_tests/test_io.py
index ffae401a8..4039431cc 100644
--- a/test/system_tests/test_io.py
+++ b/test/system_tests/test_io.py
@@ -9,9 +9,7 @@
 
 
 def test_writing_and_reading_of_species_function_using_checkpoints(tmpdir):
-    """
-    Tests that a model can write a checkpoint file and another model can read it.
-    """
+    """Tests that a model can write a checkpoint file and another model can read it."""
     mesh = dolfinx.mesh.create_unit_square(
         MPI.COMM_WORLD, nx=10, ny=10, cell_type=dolfinx.cpp.mesh.CellType.quadrilateral
     )
@@ -19,7 +17,7 @@ def test_writing_and_reading_of_species_function_using_checkpoints(tmpdir):
     my_model.mesh = F.Mesh(mesh)
 
     my_mat = F.Material(name="mat", D_0=1, E_D=0)
-    vol = F.VolumeSubdomain(id=0, material=my_mat)
+    vol = F.VolumeSubdomain(id=1, material=my_mat)
     surf = F.SurfaceSubdomain(id=1)
     my_model.subdomains = [vol, surf]
 
@@ -117,8 +115,8 @@ def test_writing_and_reading_of_species_function_using_checkpoints(tmpdir):
 
 def test_VTXExport_times_added_to_milestones(tmpdir):
     """Creates a HydrogenTransportProblem object and checks that, if no
-    stepsize.milestones are given and VTXExport.times are given, VTXExport.times are
-    are added to stepsize.milestones by .initialise()
+    stepsize.milestones are given and VTXExport.times are given, VTXExport.times are are
+    added to stepsize.milestones by .initialise()
 
     Args:
         tmpdir (os.PathLike): path to the pytest temporary folder
@@ -156,8 +154,8 @@ def test_VTXExport_times_added_to_milestones(tmpdir):
 
 
 def test_vtx_writer_called_only_at_specified_times(tmpdir):
-    """test that the VTXWriter.write function is called the number of times specified in
-    the export.times"""
+    """Test that the VTXWriter.write function is called the number of times specified in
+    the export.times."""
 
     filename = str(tmpdir.join("my_export.bp"))
 
diff --git a/test/system_tests/test_misc.py b/test/system_tests/test_misc.py
index 47d855942..f565b2d9a 100644
--- a/test/system_tests/test_misc.py
+++ b/test/system_tests/test_misc.py
@@ -44,10 +44,11 @@ def test_petsc_options():
 
 
 def test_D_global_on_2d_mesh():
-    """
-    Test that the D_global is defined correctly on a 2D mesh with two different
-    materials. The D_global should be defined as a piecewise constant function
-    with two different values, one for each material.
+    """Test that the D_global is defined correctly on a 2D mesh with two different
+    materials.
+
+    The D_global should be defined as a piecewise constant function with two different
+    values, one for each material.
     """
     mesh, mt, ct = generate_mesh(20)
 
@@ -98,7 +99,7 @@ def test_D_global_on_2d_mesh():
     ],
 )
 def test_min_max_vol_on_2d_mesh(species):
-    """Added test that catches bug #908"""
+    """Added test that catches bug #908."""
     mesh, mt, ct = generate_mesh(10)
 
     my_model = F.HydrogenTransportProblem()
@@ -183,7 +184,7 @@ def test_min_max_vol_on_2d_mesh(species):
 
 
 def test_temp_dependent_bc_mixed_domain_temperature_as_function():
-    """Test to catch bug 986"""
+    """Test to catch bug 986."""
     mesh, mt, ct = generate_mesh(8)
 
     V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
@@ -248,8 +249,8 @@ def test_del():
 
 
 def test_multispecies_with_immobile():
-    """test to catch bug 1036, that the diffusion coefficent for a non-mobile species
-    is not requested in the dict given to material"""
+    """Test to catch bug 1036, that the diffusion coefficent for a non-mobile species is
+    not requested in the dict given to material."""
 
     my_model = F.HydrogenTransportProblem()
     my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=50))
@@ -287,9 +288,7 @@ def test_multispecies_with_immobile():
 
 
 def test_implicit_species_bug_reaction():
-    """
-    This test catches the bug described in issue #1084
-    """
+    """This test catches the bug described in issue #1084."""
     tungsten = F.Material(
         D_0=1,
         E_D=0,
@@ -336,7 +335,7 @@ def test_implicit_species_bug_reaction():
 
 
 def test_timesteps():
-    """Test that the timesteps are correctly saved in the model"""
+    """Test that the timesteps are correctly saved in the model."""
     my_model = F.HydrogenTransportProblem()
     my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=50))
 
diff --git a/test/system_tests/test_multi_mat_penalty.py b/test/system_tests/test_multi_mat_penalty.py
index 116da435e..db693e58c 100644
--- a/test/system_tests/test_multi_mat_penalty.py
+++ b/test/system_tests/test_multi_mat_penalty.py
@@ -25,7 +25,6 @@ def half(x):
     )
 
     # Split domain in half and set an interface tag of 5
-    gdim = mesh.geometry.dim
     tdim = mesh.topology.dim
     fdim = tdim - 1
     top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary)
@@ -71,12 +70,12 @@ def test_2_materials_2d_mms(tmpdir):
     K_S_bot = 6.0
     D_top = 2.0
     D_bot = 5.0
-    c_exact_top_ufl = lambda x: (
-        3 + ufl.sin(ufl.pi * (2 * x[1] + 0.5)) + 0.1 * ufl.cos(2 * ufl.pi * x[0])
-    )
-    c_exact_top_np = lambda x: (
-        3 + np.sin(np.pi * (2 * x[1] + 0.5)) + 0.1 * np.cos(2 * np.pi * x[0])
-    )
+
+    def c_exact_top_ufl(x):
+        return 3 + ufl.sin(ufl.pi * (2 * x[1] + 0.5)) + 0.1 * ufl.cos(2 * ufl.pi * x[0])
+
+    def c_exact_top_np(x):
+        return 3 + np.sin(np.pi * (2 * x[1] + 0.5)) + 0.1 * np.cos(2 * np.pi * x[0])
 
     def c_exact_bot_ufl(x):
         return K_S_bot / K_S_top**2 * c_exact_top_ufl(x) ** 2
diff --git a/test/system_tests/test_multi_material.py b/test/system_tests/test_multi_material.py
index 86526d7b4..700c8cebf 100644
--- a/test/system_tests/test_multi_material.py
+++ b/test/system_tests/test_multi_material.py
@@ -3,6 +3,7 @@
 import dolfinx
 import dolfinx.fem.petsc
 import numpy as np
+import pytest
 import ufl
 
 import festim as F
@@ -25,7 +26,6 @@ def half(x):
     )
 
     # Split domain in half and set an interface tag of 5
-    gdim = mesh.geometry.dim
     tdim = mesh.topology.dim
     fdim = tdim - 1
     top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary)
@@ -71,16 +71,15 @@ def test_2_materials_2d_mms(tmpdir):
     K_S_bot = 6.0
     D_top = 2.0
     D_bot = 5.0
-    c_exact_top_ufl = lambda x: (
-        1 + ufl.sin(ufl.pi * (2 * x[0] + 0.5)) + ufl.cos(2 * ufl.pi * x[1])
-    )
+
+    def c_exact_top_ufl(x):
+        return 1 + ufl.sin(ufl.pi * (2 * x[0] + 0.5)) + ufl.cos(2 * ufl.pi * x[1])
 
     def c_exact_bot_ufl(x):
         return K_S_bot / K_S_top * c_exact_top_ufl(x)
 
-    c_exact_top_np = lambda x: (
-        1 + np.sin(np.pi * (2 * x[0] + 0.5)) + np.cos(2 * np.pi * x[1])
-    )
+    def c_exact_top_np(x):
+        return 1 + np.sin(np.pi * (2 * x[0] + 0.5)) + np.cos(2 * np.pi * x[1])
 
     def c_exact_bot_np(x):
         return K_S_bot / K_S_top * c_exact_top_np(x)
@@ -126,12 +125,14 @@ def c_exact_bot_np(x):
         F.FixedConcentrationBC(bottom_surface, value=c_exact_bot_ufl, species=H),
     ]
 
-    source_top_val = lambda x: (
-        8 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1]))
-    )
-    source_bottom_val = lambda x: (
-        40 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1]))
-    )
+    def source_top_val(x):
+        return 8 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1]))
+
+    def source_bottom_val(x):
+        return (
+            40 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1]))
+        )
+
     my_model.sources = [
         F.ParticleSource(volume=top_domain, species=H, value=source_top_val),
         F.ParticleSource(volume=bottom_domain, species=H, value=source_bottom_val),
@@ -263,8 +264,8 @@ def test_3_materials_transient(tmpdir):
     left_surface = F.SurfaceSubdomain1D(id=1, x=vertices[0])
     right_surface = F.SurfaceSubdomain1D(id=2, x=vertices[-1])
 
-    # the ids here are arbitrary in 1D, you can put anything as long as it's not the same as the surfaces
-    # TODO remove mesh and meshtags from these arguments
+    # the ids here are arbitrary in 1D, you can put anything as long as it's not the
+    # same as the surfaces
     my_model.interfaces = [
         F.Interface(6, (left_domain, middle_domain)),
         F.Interface(7, (middle_domain, right_domain)),
@@ -396,3 +397,52 @@ def test_2_mats_particle_flux_bc(tmpdir):
 
     my_model.initialise()
     my_model.run()
+
+
+def test_all_cells_are_not_tagged():
+    """Checks that an error is raised when not all cells are tagged with a non-zero
+    value.
+
+    This can be caused by a volume subdomain not being defined correctly, or by a mesh
+    that is too coarse to capture the geometry of the volume subdomains.
+    """
+
+    my_model = F.HydrogenTransportProblemDiscontinuous()
+
+    mat = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0)
+
+    vol1 = F.VolumeSubdomain1D(id=1, borders=[0, 0.5], material=mat)
+    vol2 = F.VolumeSubdomain1D(id=2, borders=[0.5, 1], material=mat)
+
+    surf1 = F.SurfaceSubdomain1D(id=1, x=0)
+    surf2 = F.SurfaceSubdomain1D(id=2, x=1)
+
+    my_model.subdomains = [vol1, vol2, surf1, surf2]
+
+    my_model.surface_to_volume = {
+        surf1: vol1,
+        surf2: vol2,
+    }
+    my_model.interfaces = [F.Interface(id=3, subdomains=[vol1, vol2])]
+
+    my_model.mesh = F.Mesh1D(np.linspace(0, 1, 10))
+
+    my_model.species = [F.Species("H", subdomains=[vol1, vol2])]
+
+    my_model.boundary_conditions = [
+        F.FixedConcentrationBC(species=my_model.species[0], subdomain=surf1, value=1),
+    ]
+
+    my_model.temperature = 300
+
+    my_model.settings = F.Settings(transient=False, atol=1e-9, rtol=1e-9)
+
+    my_model.exports = [
+        F.AverageVolume(field=my_model.species[0], volume=vol1),
+        F.AverageVolume(field=my_model.species[0], volume=vol2),
+    ]
+
+    with pytest.raises(
+        AssertionError, match="All cells must be tagged with a non-zero value"
+    ):
+        my_model.initialise()
diff --git a/test/system_tests/test_multi_material_change_of_var.py b/test/system_tests/test_multi_material_change_of_var.py
index 8091a2663..14e68af7b 100644
--- a/test/system_tests/test_multi_material_change_of_var.py
+++ b/test/system_tests/test_multi_material_change_of_var.py
@@ -24,7 +24,6 @@ def half(x):
     )
 
     # Split domain in half and set an interface tag of 5
-    gdim = mesh.geometry.dim
     tdim = mesh.topology.dim
     fdim = tdim - 1
     top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary)
@@ -156,16 +155,15 @@ def test_2_materials_2d_mms(tmpdir):
     K_S_bot = 6.0
     D_top = 2.0
     D_bot = 5.0
-    c_exact_top_ufl = lambda x: (
-        1 + ufl.sin(ufl.pi * (2 * x[0] + 0.5)) + ufl.cos(2 * ufl.pi * x[1])
-    )
+
+    def c_exact_top_ufl(x):
+        return 1 + ufl.sin(ufl.pi * (2 * x[0] + 0.5)) + ufl.cos(2 * ufl.pi * x[1])
 
     def c_exact_bot_ufl(x):
         return K_S_bot / K_S_top * c_exact_top_ufl(x)
 
-    c_exact_top_np = lambda x: (
-        1 + np.sin(np.pi * (2 * x[0] + 0.5)) + np.cos(2 * np.pi * x[1])
-    )
+    def c_exact_top_np(x):
+        return 1 + np.sin(np.pi * (2 * x[0] + 0.5)) + np.cos(2 * np.pi * x[1])
 
     def c_exact_bot_np(x):
         return K_S_bot / K_S_top * c_exact_top_np(x)
@@ -201,12 +199,14 @@ def c_exact_bot_np(x):
         F.FixedConcentrationBC(bottom_surface, value=c_exact_bot_ufl, species=H),
     ]
 
-    source_top_val = lambda x: (
-        8 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1]))
-    )
-    source_bottom_val = lambda x: (
-        40 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1]))
-    )
+    def source_top_val(x):
+        return 8 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1]))
+
+    def source_bottom_val(x):
+        return (
+            40 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1]))
+        )
+
     my_model.sources = [
         F.ParticleSource(volume=top_domain, species=H, value=source_top_val),
         F.ParticleSource(volume=bottom_domain, species=H, value=source_bottom_val),
diff --git a/test/system_tests/test_non_homogeneous_density.py b/test/system_tests/test_non_homogeneous_density.py
index c021641f8..5b2fa4176 100644
--- a/test/system_tests/test_non_homogeneous_density.py
+++ b/test/system_tests/test_non_homogeneous_density.py
@@ -6,12 +6,21 @@
 
 import festim as F
 
-step_function_space = lambda x: ufl.conditional(ufl.gt(x[0], 0.5), 10, 0.0)
-step_function_time_non_homogeneous = lambda x, t: ufl.conditional(
-    ufl.gt(t, 50), 10 - 10 * x[0], 0.0
-)
-step_function_time_homogeneous = lambda t: 10.0 if t > 50 else 2.0
-simple_space = lambda x: 10 - 10 * x[0]
+
+def step_function_space(x):
+    return ufl.conditional(ufl.gt(x[0], 0.5), 10, 0.0)
+
+
+def step_function_time_non_homogeneous(x, t):
+    return ufl.conditional(ufl.gt(t, 50), 10 - 10 * x[0], 0.0)
+
+
+def step_function_time_homogeneous(t):
+    return 10.0 if t > 50 else 2.0
+
+
+def simple_space(x):
+    return 10 - 10 * x[0]
 
 
 @pytest.mark.parametrize(
diff --git a/test/system_tests/test_reaction_not_in_every_volume.py b/test/system_tests/test_reaction_not_in_every_volume.py
index b1cdea67b..5b9fd8e22 100644
--- a/test/system_tests/test_reaction_not_in_every_volume.py
+++ b/test/system_tests/test_reaction_not_in_every_volume.py
@@ -4,8 +4,8 @@
 
 
 def test_sim_reaction_not_in_every_volume():
-    """Tests that a steady simulation can be run if a reaction is not defined
-    in every volume"""
+    """Tests that a steady simulation can be run if a reaction is not defined in every
+    volume."""
     my_model = F.HydrogenTransportProblem()
     my_model.mesh = F.Mesh1D(np.linspace(0, 2e-04, num=1001))
 
diff --git a/test/system_tests/test_spherical_coordinates.py b/test/system_tests/test_spherical_coordinates.py
index d477f484d..5a59b493b 100644
--- a/test/system_tests/test_spherical_coordinates.py
+++ b/test/system_tests/test_spherical_coordinates.py
@@ -7,15 +7,14 @@
 
 
 def test_run_MMS_spherical():
-    """
-    Tests that festim produces the correct concentration field in spherical
-    coordinates
-    """
+    """Tests that festim produces the correct concentration field in spherical
+    coordinates."""
 
     my_mesh = F.Mesh1D(vertices=np.linspace(1, 2, 1000), coordinate_system="spherical")
-    V = fem.functionspace(my_mesh.mesh, ("Lagrange", 1))
+    fem.functionspace(my_mesh.mesh, ("Lagrange", 1))
 
-    u_exact = lambda x: 3 + 2 * x[0] ** 2
+    def u_exact(x):
+        return 3 + 2 * x[0] ** 2
 
     f = -12
 
@@ -68,10 +67,8 @@ def test_run_MMS_spherical():
 
 
 def test_run_MMS_spherical_mixed_domain():
-    """
-    Tests that festim produces the correct concentration field in spherical
-    coordinates in a discontinuous domain with two materials
-    """
+    """Tests that festim produces the correct concentration field in spherical
+    coordinates in a discontinuous domain with two materials."""
 
     my_model = F.HydrogenTransportProblemDiscontinuous()
 
@@ -99,7 +96,8 @@ def c_exact_left(x):
     def c_exact_right(x):
         return K_S_right / K_S_left * c_exact_left(x)
 
-    lap_c = lambda r: -(4 * r_interface / r) + 6
+    def lap_c(r):
+        return -(4 * r_interface / r) + 6
 
     mat_1 = F.Material(D_0=D, E_D=0, K_S_0=K_S_left, E_K_S=0, solubility_law="sievert")
     mat_2 = F.Material(D_0=D, E_D=0, K_S_0=K_S_right, E_K_S=0, solubility_law="sievert")
@@ -131,8 +129,11 @@ def c_exact_right(x):
 
     my_model.temperature = 500
 
-    f_left = lambda x: -D * lap_c(x[0])
-    f_right = lambda x: -D * K_S_right / K_S_left * lap_c(x[0])
+    def f_left(x):
+        return -D * lap_c(x[0])
+
+    def f_right(x):
+        return -D * K_S_right / K_S_left * lap_c(x[0])
 
     my_model.sources = [
         F.ParticleSource(value=f_left, volume=vol_1, species=H),
diff --git a/test/system_tests/test_surface_reactions.py b/test/system_tests/test_surface_reactions.py
index 6754c10f1..7efaa2bad 100644
--- a/test/system_tests/test_surface_reactions.py
+++ b/test/system_tests/test_surface_reactions.py
@@ -22,13 +22,12 @@ def compute(self, u, ds, entity_maps=None):
 
 
 def test_2_isotopes_no_pressure():
-    """
-    Runs a simple 1D hydrogen transport problem with 2 isotopes
-    and 3 surface reactions on the right boundary.
+    """Runs a simple 1D hydrogen transport problem with 2 isotopes and 3 surface
+    reactions on the right boundary.
 
-    Then checks that the fluxes of the isotopes are consistent with the surface reactions
-    by computing the flux of H and D (from the gradient) and comparing it to the fluxes
-    computed from the surface reactions.
+    Then checks that the fluxes of the isotopes are consistent with the
+    surface reactions by computing the flux of H and D (from the gradient)
+    and comparing it to the fluxes computed from the surface reactions.
 
     Example: H + D <-> HD, H + H <-> HH, D + D <-> DD
     -D grad(c_H) n = 2*kr*c_H*c_H + kr*c_H*c_D
@@ -149,13 +148,12 @@ def test_2_isotopes_no_pressure():
 
 
 def test_2_isotopes_with_pressure():
-    """
-    Runs a simple 1D hydrogen transport problem with 2 isotopes
-    and 3 surface reactions on the right boundary.
+    """Runs a simple 1D hydrogen transport problem with 2 isotopes and 3 surface
+    reactions on the right boundary.
 
-    Then checks that the fluxes of the isotopes are consistent with the surface reactions
-    by computing the flux of H and D (from the gradient) and comparing it to the fluxes
-    computed from the surface reactions.
+    Then checks that the fluxes of the isotopes are consistent with the
+    surface reactions by computing the flux of H and D (from the gradient) and
+    comparing it to the fluxes computed from the surface reactions.
 
     Example: H + D <-> HD, H + H <-> HH, D + D <-> DD
     -D grad(c_H) n = 2*kr*c_H*c_H + kr*c_H*c_D - kd*P_H2
@@ -276,9 +274,8 @@ def test_2_isotopes_with_pressure():
 
 
 def test_pressure_varies_in_time():
-    """
-    Runs a problem with a surface reaction and a time-dependent pressure
-    on the right boundary.
+    """Runs a problem with a surface reaction and a time-dependent pressure on the right
+    boundary.
 
     Then checks that the flux is consistent with the surface reaction
     """
diff --git a/test/test_advection.py b/test/test_advection.py
index 88ed2f857..17c00e2fc 100644
--- a/test/test_advection.py
+++ b/test/test_advection.py
@@ -91,7 +91,7 @@ def test_subdomain_accepts_None_value():
 
 
 def test_velocity_field_update():
-    """Test a explicit time dependent input value is updated correctly"""
+    """Test a explicit time dependent input value is updated correctly."""
 
     test_mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
     V = fem.functionspace(test_mesh, ("Lagrange", 1))
diff --git a/test/test_average_surface.py b/test/test_average_surface.py
index 5ab11170e..74f875736 100644
--- a/test/test_average_surface.py
+++ b/test/test_average_surface.py
@@ -6,7 +6,7 @@
 
 
 def test_average_surface_compute_1D():
-    """Test that the average surface export computes the correct value"""
+    """Test that the average surface export computes the correct value."""
 
     # BUILD
     L = 6.0
@@ -16,7 +16,7 @@ def test_average_surface_compute_1D():
     dummy_volume = F.VolumeSubdomain1D(
         id=1, borders=[0, L], material=F.Material(D_0=1, E_D=1, name="dummy")
     )
-    facet_meshtags, temp = my_mesh.define_meshtags(
+    facet_meshtags, _temp = my_mesh.define_meshtags(
         surface_subdomains=[dummy_surface], volume_subdomains=[dummy_volume]
     )
     ds = ufl.Measure("ds", domain=my_mesh.mesh, subdomain_data=facet_meshtags)
diff --git a/test/test_average_volume.py b/test/test_average_volume.py
index 7bb93d42e..de1fddd8b 100644
--- a/test/test_average_volume.py
+++ b/test/test_average_volume.py
@@ -8,13 +8,13 @@
 
 
 def test_average_volume_compute_1D():
-    """Test that the average volume export computes the correct value"""
+    """Test that the average volume export computes the correct value."""
 
     # BUILD
     L = 6.0
     my_mesh = F.Mesh1D(np.linspace(0, L, 10000))
     dummy_volume = F.VolumeSubdomain1D(id=1, borders=[0, L], material=dummy_mat)
-    temp, cell_meshtags = my_mesh.define_meshtags(
+    _temp, cell_meshtags = my_mesh.define_meshtags(
         surface_subdomains=[], volume_subdomains=[dummy_volume]
     )
     dx = ufl.Measure("dx", domain=my_mesh.mesh, subdomain_data=cell_meshtags)
diff --git a/test/test_complex_reaction.py b/test/test_complex_reaction.py
index bcbbb7923..170e68c98 100644
--- a/test/test_complex_reaction.py
+++ b/test/test_complex_reaction.py
@@ -5,10 +5,10 @@
 
 
 def concentration_A_exact(t, c_A_0, k, p):
-    """Analytical solution for the concentration of species A in a reaction A + B <-> C + D
-    assuming [A]_0 = [B]_0 and [C]_0 = [D]_0 = 0
-    can be obtained by solving
-    d[A]/dt = -k[A][B] + p[C][D]
+    """Analytical solution for the concentration of species A in a reaction A + B <-> C
+    + D assuming [A]_0 = [B]_0 and [C]_0 = [D]_0 = 0 can be obtained by solving d[A]/dt
+    = -k[A][B] + p[C][D]
+
     where [A] = [B] and [C] = [D] = 0 at t=0
     Moreover, [C] = [D] = [A]_0 - [A]
     The equation then becomes
@@ -39,8 +39,8 @@ def concentration_A_exact(t, c_A_0, k, p):
 
 
 def model_test_reaction(stepsize=1, k=350e-4, p=120e-4, c_A_0=1):
-    """Creates a festim model with a single reaction and runs it.
-    The reaction is A + B <-> C + D
+    """Creates a festim model with a single reaction and runs it. The reaction is A + B
+    <-> C + D.
 
     Args:
         stepsize (float): the stepsize
@@ -142,7 +142,8 @@ def model_test_reaction(stepsize=1, k=350e-4, p=120e-4, c_A_0=1):
 
 
 def compute_error(model):
-    """Computes the relative L2 error norm between the analytical solution and the numerical solution
+    """Computes the relative L2 error norm between the analytical solution and the
+    numerical solution.
 
     Args:
         model (F.HydrogenTransportModel): the numerical festim model
@@ -170,14 +171,15 @@ def compute_error(model):
 
 @pytest.mark.parametrize("k, p, c_A_0", [(350e-4, 120e-4, 3), (200e-4, 100e-4, 2)])
 def test_reaction(k, p, c_A_0):
-    """Test the reaction A + B <-> C + D with a festim model and compare the results with the analytical solution"""
+    """Test the reaction A + B <-> C + D with a festim model and compare the results
+    with the analytical solution."""
     model = model_test_reaction(stepsize=0.5, k=k, p=p, c_A_0=c_A_0)
     error = compute_error(model)
     assert error < 1e-2
 
 
 def test_surface_reaction_in_discontinuous_case():
-    """Test the SurfaceReactionBC in a discontinuous problem setup"""
+    """Test the SurfaceReactionBC in a discontinuous problem setup."""
     my_model = F.HydrogenTransportProblemDiscontinuous()
 
     my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=1000))
diff --git a/test/test_coupled_heat_hydrogen_problem.py b/test/test_coupled_heat_hydrogen_problem.py
index c905d74ac..b7e51f0fb 100644
--- a/test/test_coupled_heat_hydrogen_problem.py
+++ b/test/test_coupled_heat_hydrogen_problem.py
@@ -15,13 +15,13 @@
 )
 def test_error_raised_when_wrong_hydrogen_problem_given(object):
     """Test TypeError is raised when an object that isnt a
-    festim.HydrogenTransportProblem object is given to hydrogen_problem"""
+    festim.HydrogenTransportProblem object is given to hydrogen_problem."""
 
     test_heat_problem = F.HeatTransferProblem()
 
     with pytest.raises(
         TypeError,
-        match="hydrogen_problem must be a festim.HydrogenTransportProblem object",
+        match=r"hydrogen_problem must be a festim.HydrogenTransportProblem object",
     ):
         F.CoupledTransientHeatTransferHydrogenTransport(
             heat_problem=test_heat_problem, hydrogen_problem=object
@@ -37,7 +37,7 @@ def test_error_raised_when_wrong_hydrogen_problem_given(object):
 )
 def test_error_raised_when_wrong_type_hydrogen_problem_given(object):
     """Test TypeError is raised when an object that isnt a
-    festim.HydrogenTransportProblem object is given to hydrogen_problem"""
+    festim.HydrogenTransportProblem object is given to hydrogen_problem."""
 
     test_heat_problem = F.HeatTransferProblem()
 
@@ -58,14 +58,14 @@ def test_error_raised_when_wrong_type_hydrogen_problem_given(object):
     ["coucou", 1, 1.0],
 )
 def test_error_raised_when_wrong_heat_problem_given(object):
-    """Test TypeError is raised when an object that isnt a
-    festim.HeatTransferProblem object is given to heat_problem"""
+    """Test TypeError is raised when an object that isnt a festim.HeatTransferProblem
+    object is given to heat_problem."""
 
     test_hydrogen_problem = F.HydrogenTransportProblem()
 
     with pytest.raises(
         TypeError,
-        match="heat_problem must be a festim.HeatTransferProblem object",
+        match=r"heat_problem must be a festim.HeatTransferProblem object",
     ):
         F.CoupledTransientHeatTransferHydrogenTransport(
             heat_problem=object, hydrogen_problem=test_hydrogen_problem
@@ -74,7 +74,7 @@ def test_error_raised_when_wrong_heat_problem_given(object):
 
 def test_initial_dt_values_are_the_same():
     """Test that the smallest of the stepsize intial_value values given is used in both
-    the heat_problem and the hydrogen_problem"""
+    the heat_problem and the hydrogen_problem."""
 
     test_heat_problem = F.HeatTransferProblem(
         mesh=test_mesh,
@@ -107,7 +107,7 @@ def test_initial_dt_values_are_the_same():
 
 
 def test_dts_always_the_same():
-    """Test that the correct dt value is modified"""
+    """Test that the correct dt value is modified."""
 
     dt = F.Stepsize(
         1000, growth_factor=1000, cutback_factor=0.9, target_nb_iterations=4
@@ -150,7 +150,7 @@ def test_dts_always_the_same():
 
 def test_error_raised_when_final_times_not_the_same():
     """Test that an error is raised when the final time values given in the heat_problem
-    and the hydrogen_problem are not the same"""
+    and the hydrogen_problem are not the same."""
 
     test_heat_problem = F.HeatTransferProblem(
         mesh=test_mesh,
@@ -184,7 +184,7 @@ def test_error_raised_when_final_times_not_the_same():
 
 def test_error_raised_when_both_problems_not_transient():
     """Test that an error is raised when the heat_problem and hydorgen_problem are not
-    set to transient"""
+    set to transient."""
 
     test_heat_problem = F.HeatTransferProblem(
         settings=F.Settings(atol=1, rtol=0.9999, transient=False),
diff --git a/test/test_dirichlet_bc.py b/test/test_dirichlet_bc.py
index 0dd366644..bb632552c 100644
--- a/test/test_dirichlet_bc.py
+++ b/test/test_dirichlet_bc.py
@@ -15,7 +15,7 @@
 
 
 def test_init():
-    """Test that the attributes are set correctly"""
+    """Test that the attributes are set correctly."""
     # create a DirichletBC object
     subdomain = F.SurfaceSubdomain1D(1, x=0)
     value = 1.0
@@ -31,9 +31,8 @@ def test_init():
 
 
 def test_value_fenics():
-    """Test that the value_fenics attribute can be set to a valid value
-    and that an invalid type throws an error
-    """
+    """Test that the value_fenics attribute can be set to a valid value and that an
+    invalid type throws an error."""
     # create a DirichletBC object
     subdomain = F.SurfaceSubdomain1D(1, x=0)
     value = 1.0
@@ -53,7 +52,7 @@ def test_value_fenics():
 
 
 def test_callable_for_value():
-    """Test that the value attribute can be a callable function of x and t"""
+    """Test that the value attribute can be a callable function of x and t."""
 
     subdomain = F.SurfaceSubdomain1D(1, x=1)
     vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
@@ -97,7 +96,7 @@ def value(x, t):
 
 
 def test_value_callable_x_t_T():
-    """Test that the value attribute can be a callable function of x, t, and T"""
+    """Test that the value attribute can be a callable function of x, t, and T."""
 
     subdomain = F.SurfaceSubdomain1D(1, x=1)
     vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
@@ -145,7 +144,7 @@ def value(x, t, T):
 
 @pytest.mark.parametrize("value", [lambda t: t, lambda t: 1.0 + t])
 def test_callable_t_only(value):
-    """Test that the value attribute can be a callable function of t only"""
+    """Test that the value attribute can be a callable function of t only."""
 
     subdomain = F.SurfaceSubdomain1D(1, x=1)
     vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
@@ -186,7 +185,7 @@ def test_callable_t_only(value):
 
 
 def test_callable_x_only():
-    """Test that the value attribute can be a callable function of x only"""
+    """Test that the value attribute can be a callable function of x only."""
 
     # BUILD
     subdomain = F.SurfaceSubdomain1D(1, x=1)
@@ -247,8 +246,8 @@ def value(x):
     ],
 )
 def test_integration_with_HTransportProblem(value):
-    """test that different callable functions can be applied to a dirichlet
-    boundary condition, asserting in each case they match an expected value"""
+    """Test that different callable functions can be applied to a dirichlet boundary
+    condition, asserting in each case they match an expected value."""
     subdomain = F.SurfaceSubdomain1D(1, x=1)
     vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
 
@@ -311,8 +310,8 @@ def test_integration_with_HTransportProblem(value):
     ],
 )
 def test_define_value_error_if_ufl_conditional_t_only(value):
-    """Test that a ValueError is raised when the value attribute is a callable
-    of t only and contains a ufl conditional"""
+    """Test that a ValueError is raised when the value attribute is a callable of t only
+    and contains a ufl conditional."""
 
     subdomain = F.SurfaceSubdomain1D(1, x=1)
     species = F.Species("test")
@@ -322,14 +321,14 @@ def test_define_value_error_if_ufl_conditional_t_only(value):
     t = fem.Constant(mesh, 0.0)
     V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
     with pytest.raises(
-        ValueError, match="self.value should return a float or an int, not "
+        ValueError, match=r"self.value should return a float or an int, not "
     ):
         bc.create_value(V, temperature=None, t=t)
 
 
 def test_species_predefined():
-    """Test a ValueError is raised when the species defined in the boundary
-    condition is not predefined in the model"""
+    """Test a ValueError is raised when the species defined in the boundary condition is
+    not predefined in the model."""
 
     subdomain = F.SurfaceSubdomain1D(1, x=1)
     vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
@@ -362,9 +361,9 @@ def test_species_predefined():
     ],
 )
 def test_integration_with_a_multispecies_HTransportProblem(value_A, value_B):
-    """test that a mixture of callable functions can be applied to dirichlet
-    boundary conditions in a multispecies case, asserting in each case they
-    match an expected value"""
+    """Test that a mixture of callable functions can be applied to dirichlet boundary
+    conditions in a multispecies case, asserting in each case they match an expected
+    value."""
     subdomain_A = F.SurfaceSubdomain1D(1, x=0)
     subdomain_B = F.SurfaceSubdomain1D(2, x=1)
     vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
@@ -438,7 +437,7 @@ def test_integration_with_a_multispecies_HTransportProblem(value_A, value_B):
     ],
 )
 def test_bc_time_dependent_attribute(input, expected_value):
-    """Test that the time_dependent attribute is correctly set"""
+    """Test that the time_dependent attribute is correctly set."""
     subdomain = F.SurfaceSubdomain1D(1, x=0)
     species = F.Species("test")
     my_bc = F.DirichletBC(subdomain, input, species)
@@ -463,7 +462,7 @@ def test_bc_time_dependent_attribute(input, expected_value):
     ],
 )
 def test_bc_temperature_dependent_attribute(input, expected_value):
-    """Test that the temperature_dependent attribute is correctly set"""
+    """Test that the temperature_dependent attribute is correctly set."""
     subdomain = F.SurfaceSubdomain1D(1, x=0)
     species = F.Species("test")
     my_bc = F.DirichletBC(subdomain, input, species)
diff --git a/test/test_flux_bc.py b/test/test_flux_bc.py
index 04ff92183..d048d3277 100644
--- a/test/test_flux_bc.py
+++ b/test/test_flux_bc.py
@@ -14,7 +14,7 @@
 
 
 def test_init():
-    """Test that the attributes are set correctly"""
+    """Test that the attributes are set correctly."""
     # create a DirichletBC object
     subdomain = F.SurfaceSubdomain1D(1, x=0)
     value = 1.0
@@ -46,7 +46,7 @@ def test_init():
     ],
 )
 def test_create_value_fenics_type(value, expected_type):
-    """Test that"""
+    """Test that."""
     # BUILD
     left = F.SurfaceSubdomain1D(1, x=0)
     my_species = F.Species("test")
@@ -72,7 +72,7 @@ def test_create_value_fenics_type(value, expected_type):
     ],
 )
 def test_create_value_fenics_value(value, expected_value):
-    """Test that"""
+    """Test that."""
     # BUILD
     left = F.SurfaceSubdomain1D(1, x=0)
     my_species = F.Species("test")
@@ -90,7 +90,8 @@ def test_create_value_fenics_value(value, expected_value):
 
 
 def test_create_value_fenics_dependent_conc():
-    """Test that the value_fenics of ParticleFluxBC is set correctly when the value is dependent on the concentration"""
+    """Test that the value_fenics of ParticleFluxBC is set correctly when the value is
+    dependent on the concentration."""
     # BUILD
     left = F.SurfaceSubdomain1D(1, x=0)
     my_species = F.Species("test")
@@ -119,14 +120,17 @@ def test_value_fenics_setter_error():
 
     with pytest.raises(
         TypeError,
-        match="Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, np.ndarray or ufl.core.expr.Expr not <class 'str'>",
+        match=(
+            r"Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, np.ndarray or "  # noqa: E501
+            r"ufl.core.expr.Expr not <class 'str'>"
+        ),
     ):
         bc.value_fenics = "coucou"
 
 
 def test_ValueError_raised_when_callable_returns_wrong_type():
     """The create_value_fenics method should raise a ValueError when the callable
-    returns an object which is not a float or int"""
+    returns an object which is not a float or int."""
 
     surface = F.SurfaceSubdomain(id=1)
     species = F.Species("test")
@@ -141,7 +145,10 @@ def my_value(t):
 
     with pytest.raises(
         ValueError,
-        match="self.value should return a float or an int, not <class 'ufl.conditional.Conditional'",
+        match=(
+            r"self.value should return a float or an int, not <class "
+            r"'ufl.conditional.Conditional'"
+        ),
     ):
         bc.create_value_fenics(mesh, T, t)
 
@@ -160,7 +167,7 @@ def my_value(t):
     ],
 )
 def test_bc_time_dependent_attribute(input, expected_value):
-    """Test that the time_dependent attribute is correctly set"""
+    """Test that the time_dependent attribute is correctly set."""
     surface = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
     my_species = F.Species("test")
     my_bc = F.ParticleFluxBC(subdomain=surface, value=input, species=my_species)
@@ -169,7 +176,8 @@ def test_bc_time_dependent_attribute(input, expected_value):
 
 
 def test_bc_time_dependent_attribute_raises_error_when_value_none():
-    """Test that the time_dependent attribute raises a TypeError when the value is None"""
+    """Test that the time_dependent attribute raises a TypeError when the value is
+    None."""
     surface = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
     my_flux_bc = F.FluxBCBase(subdomain=surface, value=None)
 
@@ -197,7 +205,7 @@ def test_bc_time_dependent_attribute_raises_error_when_value_none():
     ],
 )
 def test_bc_temperature_dependent_attribute(input, expected_value):
-    """Test that the temperature_dependent attribute is correctly set"""
+    """Test that the temperature_dependent attribute is correctly set."""
     surface = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
     my_species = F.Species("test")
     my_bc = F.ParticleFluxBC(subdomain=surface, value=input, species=my_species)
@@ -206,7 +214,7 @@ def test_bc_temperature_dependent_attribute(input, expected_value):
 
 
 def test_HeatFluxBC_init():
-    """Test that the attributes are set correctly"""
+    """Test that the attributes are set correctly."""
     # create a DirichletBC object
     subdomain = F.SurfaceSubdomain1D(1, x=0)
     value = 1.0
@@ -235,7 +243,7 @@ def test_HeatFluxBC_init():
     ],
 )
 def test_create_value_fenics_type_HeatFluxBC(value, expected_type):
-    """Test that"""
+    """Test that."""
     # BUILD
     left = F.SurfaceSubdomain1D(1, x=0)
     t = F.as_fenics_constant(0, mesh)
@@ -266,7 +274,7 @@ def test_create_value_fenics_type_HeatFluxBC(value, expected_type):
     ],
 )
 def test_create_value_fenics_value_HeatFluxBC(value, expected_value):
-    """Test that"""
+    """Test that."""
     # BUILD
     left = F.SurfaceSubdomain1D(1, x=0)
     t = F.as_fenics_constant(0, mesh)
diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py
index d33e002a2..b1ddfa097 100644
--- a/test/test_h_transport_problem.py
+++ b/test/test_h_transport_problem.py
@@ -27,7 +27,7 @@
     ],
 )
 def test_temperature_setter_type(value):
-    """Test that the temperature type is correctly set"""
+    """Test that the temperature type is correctly set."""
     my_model = F.HydrogenTransportProblem(
         mesh=test_mesh,
     )
@@ -53,7 +53,7 @@ def test_temperature_setter_type(value):
     ],
 )
 def test_time_dependent_temperature_attribute(input, expected_value):
-    """Test that the temperature_time_dependent attribute is correctly set"""
+    """Test that the temperature_time_dependent attribute is correctly set."""
 
     my_model = F.HydrogenTransportProblem()
     my_model.temperature = input
@@ -62,7 +62,7 @@ def test_time_dependent_temperature_attribute(input, expected_value):
 
 
 def test_define_temperature_value_error_raised():
-    """Test that a ValueError is raised when the temperature is None"""
+    """Test that a ValueError is raised when the temperature is None."""
 
     # BUILD
     my_model = F.HydrogenTransportProblem(mesh=test_mesh)
@@ -94,9 +94,9 @@ def test_define_temperature_value_error_raised():
     ],
 )
 def test_define_temperature(input, expected_type):
-    """Test that the define_temperature method correctly sets the
-    temperature_fenics attribute to either a fem.Constant or a
-    fem.Function depending on the type of input"""
+    """Test that the define_temperature method correctly sets the temperature_fenics
+    attribute to either a fem.Constant or a fem.Function depending on the type of
+    input."""
 
     # BUILD
     my_model = F.HydrogenTransportProblem(mesh=test_mesh)
@@ -121,8 +121,8 @@ def test_define_temperature(input, expected_type):
     ],
 )
 def test_define_temperature_error_if_ufl_conditional_t_only(input):
-    """Test that a ValueError is raised when the temperature attribute is a callable
-    of t only and contains a ufl conditional"""
+    """Test that a ValueError is raised when the temperature attribute is a callable of
+    t only and contains a ufl conditional."""
     my_model = F.HydrogenTransportProblem(mesh=test_mesh)
     my_model.t = fem.Constant(test_mesh.mesh, 0.0)
 
@@ -130,13 +130,13 @@ def test_define_temperature_error_if_ufl_conditional_t_only(input):
 
     with pytest.raises(
         ValueError,
-        match="self.temperature should return a float or an int, not",
+        match=r"self.temperature should return a float or an int, not",
     ):
         my_model.define_temperature()
 
 
 def test_iterate():
-    """Test that the iterate method updates the solution and time correctly"""
+    """Test that the iterate method updates the solution and time correctly."""
     # BUILD
     my_model = F.HydrogenTransportProblem()
 
@@ -201,8 +201,8 @@ def test_iterate():
     ],
 )
 def test_update_time_dependent_values_temperature(T_function, expected_values):
-    """Test that different time-dependent callable functions for the
-    temperature are updated at each time step and match an expected value"""
+    """Test that different time-dependent callable functions for the temperature are
+    updated at each time step and match an expected value."""
 
     # BUILD
     my_model = F.HydrogenTransportProblem(
@@ -227,7 +227,7 @@ def test_update_time_dependent_values_temperature(T_function, expected_values):
 
 
 def test_initialise_exports_find_species_with_one_field():
-    """Test that a species can be found from the model species if given as a string"""
+    """Test that a species can be found from the model species if given as a string."""
 
     # BUILD
     my_model = F.HydrogenTransportProblem(
@@ -246,11 +246,11 @@ def test_initialise_exports_find_species_with_one_field():
 
 
 def test_define_D_global_different_temperatures():
-    """Test that the D_global object is correctly defined when the temperature
-    is different in the volume subdomains"""
+    """Test that the D_global object is correctly defined when the temperature is
+    different in the volume subdomains."""
     D_0, E_D = 1.5, 0.1
     my_mat = F.Material(D_0=D_0, E_D=E_D, name="my_mat")
-    surf = F.SurfaceSubdomain1D(id=1, x=0)
+    F.SurfaceSubdomain1D(id=1, x=0)
     H = F.Species("H")
 
     my_model = F.HydrogenTransportProblem(
@@ -268,7 +268,7 @@ def test_define_D_global_different_temperatures():
     my_model.t = 1
     my_model.define_temperature()
 
-    D_computed, D_expr = my_model.define_D_global(H)
+    D_computed, _D_expr = my_model.define_D_global(H)
 
     computed_values = [D_computed.x.array[0], D_computed.x.array[-1]]
 
@@ -281,8 +281,8 @@ def test_define_D_global_different_temperatures():
 
 
 def test_define_D_global_different_materials():
-    """Test that the D_global object is correctly defined when the material
-    is different in the volume subdomains"""
+    """Test that the D_global object is correctly defined when the material is different
+    in the volume subdomains."""
     D_0_left, E_D_left = 1.0, 0.1
     D_0_right, E_D_right = 2.0, 0.2
     my_mat_L = F.Material(D_0=D_0_left, E_D=E_D_left, name="my_mat_L")
@@ -304,7 +304,7 @@ def test_define_D_global_different_materials():
     my_model.t = 0
     my_model.define_temperature()
 
-    D_computed, D_expr = my_model.define_D_global(H)
+    D_computed, _D_expr = my_model.define_D_global(H)
 
     computed_values = [D_computed.x.array[0], D_computed.x.array[-1]]
 
@@ -317,9 +317,9 @@ def test_define_D_global_different_materials():
 
 
 def test_initialise_exports_multiple_exports_same_species():
-    """Test that the diffusion coefficient within the D_global object function is the same
-    for multiple exports of the same species, and that D_global object is only
-    created once per species"""
+    """Test that the diffusion coefficient within the D_global object function is the
+    same for multiple exports of the same species, and that D_global object is only
+    created once per species."""
 
     D_0, E_D = 1.5, 0.1
     my_mat = F.Material(D_0=D_0, E_D=E_D, name="my_mat")
@@ -359,7 +359,8 @@ def test_initialise_exports_multiple_exports_same_species():
 
 
 def test_export_resets_quantities():
-    """Test that the export.data and export.t are correctly reset every time a simulation is initiated."""
+    """Test that the export.data and export.t are correctly reset every time a
+    simulation is initiated."""
     my_mat = F.Material(D_0=1, E_D=0)
     H = F.Species("H")
     surf = F.SurfaceSubdomain1D(id=1, x=4)
@@ -389,7 +390,7 @@ def test_export_resets_quantities():
 
 def test_define_D_global_multispecies():
     """Test that the D_global object is correctly defined when there are multiple
-    species in one subdomain"""
+    species in one subdomain."""
     A = F.Species("A")
     B = F.Species("B")
 
@@ -399,7 +400,7 @@ def test_define_D_global_multispecies():
     my_mat = F.Material(
         D_0={A: D_0_A, B: D_0_B}, E_D={A: E_D_A, B: E_D_B}, name="my_mat"
     )
-    surf = F.SurfaceSubdomain1D(id=1, x=1)
+    F.SurfaceSubdomain1D(id=1, x=1)
 
     my_model = F.HydrogenTransportProblem(
         mesh=F.Mesh1D(np.linspace(0, 1, num=101)),
@@ -415,8 +416,8 @@ def test_define_D_global_multispecies():
     my_model.t = 0
     my_model.define_temperature()
 
-    D_A_computed, D_A_expr = my_model.define_D_global(A)
-    D_B_computed, D_B_expr = my_model.define_D_global(B)
+    D_A_computed, _D_A_expr = my_model.define_D_global(A)
+    D_B_computed, _D_B_expr = my_model.define_D_global(B)
 
     computed_values = [D_A_computed.x.array[-1], D_B_computed.x.array[-1]]
 
@@ -429,8 +430,8 @@ def test_define_D_global_multispecies():
 
 
 def test_post_processing_update_D_global():
-    """Test that the D_global object is updated at each time
-    step when temperture is time dependent"""
+    """Test that the D_global object is updated at each time step when temperture is
+    time dependent."""
     my_mesh = F.Mesh1D(np.linspace(0, 1, num=11))
     my_mat = F.Material(D_0=1.5, E_D=0.1, name="my_mat")
     surf = F.SurfaceSubdomain1D(id=1, x=1)
@@ -480,8 +481,8 @@ def test_post_processing_update_D_global():
 
 
 def test_post_processing_update_D_global_2():
-    """Test that the D_global object is updated at each time
-    step when temperture is time dependent"""
+    """Test that the D_global object is updated at each time step when temperture is
+    time dependent."""
     my_mesh = F.Mesh1D(np.linspace(0, 1, num=11))
     my_mat = F.Material(D_0=1.5, E_D=0.1, name="my_mat")
     surf = F.SurfaceSubdomain1D(id=1, x=1)
@@ -524,8 +525,8 @@ def test_post_processing_update_D_global_2():
 
 
 def test_post_processing_update_D_global_volume_1():
-    """Test that the D_global object is updated at each time
-    step when temperture is time dependent"""
+    """Test that the D_global object is updated at each time step when temperture is
+    time dependent."""
     my_mesh = F.Mesh1D(np.linspace(0, 1, num=11))
     my_mat = F.Material(D_0=1.5, E_D=0.1, name="my_mat")
     my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat)
@@ -562,8 +563,8 @@ def test_post_processing_update_D_global_volume_1():
 
 
 def test_post_processing_update_D_global_volume_2():
-    """Test that the D_global object is updated at each time
-    step when temperture is time dependent"""
+    """Test that the D_global object is updated at each time step when temperture is
+    time dependent."""
     my_mesh = F.Mesh1D(np.linspace(0, 1, num=11))
     my_mat = F.Material(D_0=1.5, E_D=0.1, name="my_mat")
     my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat)
@@ -617,7 +618,7 @@ def test_update_time_dependent_bcs_with_time_dependent_temperature(
     temperature_value, bc_value, expected_values
 ):
     """Test that temperature dependent bcs are updated at each time step when the
-    temperature is time dependent, and match an expected value"""
+    temperature is time dependent, and match an expected value."""
 
     # BUILD
     H = F.Species("H")
@@ -666,8 +667,8 @@ def test_update_time_dependent_bcs_with_time_dependent_temperature(
     ],
 )
 def test_update_time_dependent_values_source(source_value, expected_values):
-    """Test that time dependent sources are updated at each time step,
-    and match an expected value"""
+    """Test that time dependent sources are updated at each time step, and match an
+    expected value."""
     # BUILD
     my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat)
     H = F.Species("H")
@@ -714,8 +715,8 @@ def test_update_time_dependent_values_source(source_value, expected_values):
 def test_update_sources_with_time_dependent_temperature(
     temperature_value, source_value, expected_values
 ):
-    """Test that temperature dependent source terms are updated at each time step
-    when the temperature is time dependent, and match an expected value"""
+    """Test that temperature dependent source terms are updated at each time step when
+    the temperature is time dependent, and match an expected value."""
 
     # BUILD
     H = F.Species("H")
@@ -753,7 +754,7 @@ def test_update_sources_with_time_dependent_temperature(
 
 def test_convert_source_input_values_to_fenics_objects_multispecies():
     """Test that the define_sources method correctly sets the value_fenics attribute in
-    a multispecies case"""
+    a multispecies case."""
     # BUILD
     my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat)
     H, D = F.Species("H"), F.Species("D")
@@ -784,7 +785,7 @@ def test_convert_source_input_values_to_fenics_objects_multispecies():
 
 # TODO replace this by a proper MMS test
 def test_run_in_steady_state():
-    """Test that the run method works in steady state"""
+    """Test that the run method works in steady state."""
     # BUILD
     my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat)
     my_model = F.HydrogenTransportProblem(
@@ -806,20 +807,20 @@ def test_run_in_steady_state():
 
 def test_species_setter():
     """Test that a TypeError is rasied when a species of type other than F.Species is
-    given"""
+    given."""
 
     my_model = F.HydrogenTransportProblem()
 
     with pytest.raises(
         TypeError,
-        match="elements of species must be of type festim.Species not <class 'int'>",
+        match=r"elements of species must be of type festim.Species not <class 'int'>",
     ):
         my_model.species = [1, 2, 3]
 
 
 def test_create_initial_conditions_ValueError_raised_when_not_transient():
-    """Test that ValueError is raised if initial conditions are defined in
-    a steady state simulation"""
+    """Test that ValueError is raised if initial conditions are defined in a steady
+    state simulation."""
 
     my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat)
     H = F.Species("H")
@@ -852,8 +853,8 @@ def test_create_initial_conditions_ValueError_raised_when_not_transient():
     ],
 )
 def test_create_initial_conditions_expr_fenics(input_value, expected_value):
-    """Test that after calling create_initial_conditions, the prev_solution
-    attribute of the species has the correct value at x=4.0."""
+    """Test that after calling create_initial_conditions, the prev_solution attribute of
+    the species has the correct value at x=4.0."""
 
     # BUILD
     vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 4], material=dummy_mat)
@@ -931,8 +932,8 @@ def test_create_species_from_trap():
 def test_create_initial_conditions_value_fenics_multispecies(
     input_value_1, input_value_2, expected_value_1, expected_value_2
 ):
-    """Test that after calling create_initial_conditions, the prev_solution
-    attribute of each species has the correct value at x=4.0 in a multispecies case"""
+    """Test that after calling create_initial_conditions, the prev_solution attribute of
+    each species has the correct value at x=4.0 in a multispecies case."""
 
     # BUILD
     test_mesh = F.Mesh1D(vertices=np.linspace(0, 4, num=101))
@@ -970,7 +971,7 @@ def test_create_initial_conditions_value_fenics_multispecies(
 
 
 def test_adaptive_timestepping_grows():
-    """Tests that the stepsize grows"""
+    """Tests that the stepsize grows."""
     # BUILD
     my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat)
     my_model = F.HydrogenTransportProblem(
@@ -1006,7 +1007,7 @@ def test_adaptive_timestepping_grows():
 
 
 def test_adaptive_timestepping_shrinks():
-    """Tests that the stepsize shrinks"""
+    """Tests that the stepsize shrinks."""
     # BUILD
     my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat)
     my_model = F.HydrogenTransportProblem(
@@ -1053,8 +1054,8 @@ def test_adaptive_timestepping_shrinks():
     ],
 )
 def test_reinstantiation_of_class(attribute, value):
-    """Test that when an attribute defaults to empty list, when the class
-    is reinstantiated the list is not passed to the new class object"""
+    """Test that when an attribute defaults to empty list, when the class is
+    reinstantiated the list is not passed to the new class object."""
 
     model_1 = F.HydrogenTransportProblem()
     getattr(model_1, attribute).append(value)
@@ -1064,8 +1065,8 @@ def test_reinstantiation_of_class(attribute, value):
 
 
 def test_define_meshtags_and_measures_with_custom_fenics_mesh():
-    """Test that the define_meshtags_and_measures method works when the mesh is
-    a custom fenics mesh"""
+    """Test that the define_meshtags_and_measures method works when the mesh is a custom
+    fenics mesh."""
 
     # BUILD
     mesh_1D = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10)
@@ -1096,7 +1097,7 @@ def test_define_meshtags_and_measures_with_custom_fenics_mesh():
 
 
 def test_error_raised_when_custom_fenics_mesh_wrong_facet_meshtags_type():
-    """Test the facet_meshtags type hinting raises error when given as wrong type"""
+    """Test the facet_meshtags type hinting raises error when given as wrong type."""
 
     # BUILD
     mesh_1D = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10)
@@ -1104,12 +1105,12 @@ def test_error_raised_when_custom_fenics_mesh_wrong_facet_meshtags_type():
     my_model = F.HydrogenTransportProblem(mesh=my_mesh)
 
     # TEST
-    with pytest.raises(TypeError, match="value must be of type dolfinx.mesh.MeshTags"):
+    with pytest.raises(TypeError, match=r"value must be of type dolfinx.mesh.MeshTags"):
         my_model.facet_meshtags = [0, 1]
 
 
 def test_error_raised_when_custom_fenics_mesh_wrong_volume_meshtags_type():
-    """Test the volume_meshtags type hinting raises error when given as wrong type"""
+    """Test the volume_meshtags type hinting raises error when given as wrong type."""
 
     # BUILD
     mesh_1D = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10)
@@ -1117,7 +1118,7 @@ def test_error_raised_when_custom_fenics_mesh_wrong_volume_meshtags_type():
     my_model = F.HydrogenTransportProblem(mesh=my_mesh)
 
     # TEST
-    with pytest.raises(TypeError, match="value must be of type dolfinx.mesh.MeshTags"):
+    with pytest.raises(TypeError, match=r"value must be of type dolfinx.mesh.MeshTags"):
         my_model.volume_meshtags = [0, 1]
 
 
@@ -1135,8 +1136,8 @@ def test_error_raised_when_custom_fenics_mesh_wrong_volume_meshtags_type():
     ],
 )
 def test_update_time_dependent_values_flux(bc_value, expected_values):
-    """Test that time dependent fluxes are updated at each time step,
-    and match an expected value"""
+    """Test that time dependent fluxes are updated at each time step, and match an
+    expected value."""
     # BUILD
     my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat)
     surface = F.SurfaceSubdomain1D(id=2, x=0)
@@ -1187,8 +1188,8 @@ def test_update_time_dependent_values_flux(bc_value, expected_values):
 def test_update_fluxes_with_time_dependent_temperature(
     temperature_value, bc_value, expected_values
 ):
-    """Test that temperature dependent flux terms are updated at each time step
-    when the temperature is time dependent, and match an expected value"""
+    """Test that temperature dependent flux terms are updated at each time step when the
+    temperature is time dependent, and match an expected value."""
 
     # BUILD
     H = F.Species("H")
@@ -1226,7 +1227,7 @@ def test_update_fluxes_with_time_dependent_temperature(
 
 def test_create_flux_values_fenics_multispecies():
     """Test that the create_flux_values_fenics method correctly sets the value_fenics
-    attribute in a multispecies case"""
+    attribute in a multispecies case."""
     # BUILD
     my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat)
     surface = F.SurfaceSubdomain1D(id=2, x=0)
@@ -1258,7 +1259,7 @@ def test_create_flux_values_fenics_multispecies():
 
 def test_not_implemented_error_raised_with_D_as_function():
     """Test that NotImplementedError is raised when D is a function and more than one
-    volume subdomain has been defined"""
+    volume subdomain has been defined."""
 
     # BUILD
     test_mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=101))
@@ -1293,7 +1294,7 @@ def test_not_implemented_error_raised_with_D_as_function():
 
 def test_define_D_global_with_D_as_function():
     """Test that if a function is given as diffusion coeff of a material, that it is
-    passed to D_global"""
+    passed to D_global."""
 
     # BUILD
     test_mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=101))
@@ -1323,7 +1324,7 @@ def test_define_D_global_with_D_as_function():
 @pytest.mark.parametrize("coord_sys", ["cylindrical", "spherical"])
 def test_exports_cyl_sph_coord_system_raise_not_implemented(coord_sys):
     """Test that NotImplementedError is raised when trying to use exports in a
-    cylindrical or spherical coordinate system"""
+    cylindrical or spherical coordinate system."""
 
     # BUILD
     my_mesh = F.Mesh1D(np.linspace(0, 1, num=11), coordinate_system=coord_sys)
@@ -1352,7 +1353,7 @@ def test_exports_cyl_sph_coord_system_raise_not_implemented(coord_sys):
 
 def test_traps_with_CG_elements():
     """Test that when creating species from traps, if the element for traps is set to
-    CG, the created species has a CG element as well"""
+    CG, the created species has a CG element as well."""
 
     # BUILD
     my_mesh = F.Mesh1D(np.linspace(0, 1, num=11))
@@ -1376,7 +1377,7 @@ def test_traps_with_CG_elements():
 
 
 def test_surface_reaction_BC_discontinuous():
-    """Test that surface reaction BC conserves flux"""
+    """Test that surface reaction BC conserves flux."""
 
     # BUILD
     my_model = F.HydrogenTransportProblemDiscontinuous()
@@ -1461,7 +1462,7 @@ def test_surface_reaction_BC_discontinuous():
 
 def test_temperature_as_function_in_discontinuous():
     """Test that a discontinuous problem can be initialised with a temperature field
-    given as a function"""
+    given as a function."""
 
     # BUILD
     my_model = F.HydrogenTransportProblemDiscontinuous()
@@ -1503,7 +1504,7 @@ def test_temperature_as_function_in_discontinuous():
 
 def test_wrong_element_for_immobile_species():
     """Test that an error is raised when an immobile species is defined on a subdomain
-    with an element that is incorrect"""
+    with an element that is incorrect."""
     my_model = F.HydrogenTransportProblem()
 
     with pytest.raises(ValueError, match="element_immobile should be in"):
@@ -1514,8 +1515,8 @@ def test_wrong_element_for_immobile_species():
     "element_type, expected_element", [("CG", "CG"), ("DG", "DG"), ("P", "CG")]
 )
 def test_element_for_immobile_species(element_type, expected_element):
-    """Test that the element for immobile species is correctly set when the user sets
-    it to CG or DG"""
+    """Test that the element for immobile species is correctly set when the user sets it
+    to CG or DG."""
     my_model = F.HydrogenTransportProblem()
 
     my_model.element_immobile = element_type
@@ -1524,7 +1525,7 @@ def test_element_for_immobile_species(element_type, expected_element):
 
 def test_nb_dofs_dg_or_cg():
     """Test that the number of dofs is different when using CG or DG elements for
-    immobile species"""
+    immobile species."""
 
     # BUILD
     my_model = F.HydrogenTransportProblem()
diff --git a/test/test_heat_transfer_problem.py b/test/test_heat_transfer_problem.py
index 39210a320..9cfcb7e73 100644
--- a/test/test_heat_transfer_problem.py
+++ b/test/test_heat_transfer_problem.py
@@ -110,7 +110,7 @@ def exact_solution(x):
 
 
 def test_MMS_T_dependent_thermal_cond():
-    """MMS test with space T dependent thermal cond"""
+    """MMS test with space T dependent thermal cond."""
 
     def thermal_conductivity(T):
         return 3 * T + 2
@@ -159,10 +159,8 @@ def mms_source(x):
 
 
 def test_heat_transfer_transient(tmpdir):
-    """
-    MMS test for transient heat transfer
-    constant thermal conductivity density and heat capacity
-    """
+    """MMS test for transient heat transfer constant thermal conductivity density and
+    heat capacity."""
     density = 2
     heat_capacity = 3
     thermal_conductivity = 4
@@ -207,7 +205,8 @@ def exact_solution(x, t):
     my_problem.settings = F.Settings(
         atol=1e-8,
         rtol=1e-10,
-        final_time=1,  # final time shouldn't be too long so that a potential error at the initial timestep is not negligible
+        # final time shouldn't be too long so that a potential error at the initial timestep is not negligible  # noqa: E501
+        final_time=1,
     )
 
     # Forward euler isn't great so dt should be small
@@ -222,7 +221,7 @@ def exact_solution(x, t):
     my_problem.run()
 
     computed_solution = my_problem.u
-    # we use the exact final time of the simulation which may differ from the one specified in the settings
+    # we use the exact final time of the simulation which may differ from the one specified in the settings  # noqa: E501
     final_time_sim = my_problem.t.value
 
     def exact_solution_end(x):
@@ -233,8 +232,8 @@ def exact_solution_end(x):
 
 
 def test_MES():
-    """Method of Exact Solution test for transient heat transfer
-    with thermal cond. k = 2, T = 0 on surfaces, and source term q = 8 k
+    """Method of Exact Solution test for transient heat transfer with thermal cond. k =
+    2, T = 0 on surfaces, and source term q = 8 k.
 
     Analytical solution: T(x) = 4 x (1 - x)
     """
@@ -331,7 +330,8 @@ def mms_source(x, t):
     my_problem.settings = F.Settings(
         atol=1e-8,
         rtol=1e-10,
-        final_time=1,  # final time shouldn't be too long so that a potential error at the initial timestep is not negligible
+        # final time shouldn't be too long so that a potential error at the initial timestep is not negligible  # noqa: E501
+        final_time=1,
     )
 
     # Forward euler isn't great so dt should be small
@@ -346,7 +346,7 @@ def mms_source(x, t):
     my_problem.run()
 
     computed_solution = my_problem.u
-    # we use the exact final time of the simulation which may differ from the one specified in the settings
+    # we use the exact final time of the simulation which may differ from the one specified in the settings  # noqa: E501
     final_time_sim = my_problem.t.value
 
     def exact_solution_end(x):
@@ -357,7 +357,7 @@ def exact_solution_end(x):
 
 
 def test_sources():
-    """Tests the sources setter of the HeatTransferProblem class"""
+    """Tests the sources setter of the HeatTransferProblem class."""
     htp = F.HeatTransferProblem()
     vol = F.VolumeSubdomain1D(1, borders=[0, 1], material=None)
     # Test that setting valid sources works
@@ -369,7 +369,7 @@ def test_sources():
     assert htp.sources == valid_sources
 
     # Test that setting invalid sources raises a TypeError
-    with pytest.raises(TypeError, match="festim.HeatSource objects"):
+    with pytest.raises(TypeError, match=r"festim.HeatSource objects"):
         spe = F.Species("H")
         htp.sources = [
             F.ParticleSource(1, vol, spe),
@@ -378,7 +378,7 @@ def test_sources():
 
 
 def test_boundary_conditions():
-    """Tests the boundary_conditions setter of the HeatTransferProblem class"""
+    """Tests the boundary_conditions setter of the HeatTransferProblem class."""
     htp = F.HeatTransferProblem()
     left = F.SurfaceSubdomain(1)
 
@@ -399,7 +399,8 @@ def test_boundary_conditions():
 
 @pytest.mark.parametrize("mesh", [mesh_1D, mesh_2D, mesh_3D])
 def test_meshtags_from_xdmf(tmp_path, mesh):
-    """Test that the facet and volume meshtags are read correctly from the mesh XDMF files"""
+    """Test that the facet and volume meshtags are read correctly from the mesh XDMF
+    files."""
     # create mesh functions
     fdim = mesh.topology.dim - 1
     vdim = mesh.topology.dim
@@ -482,9 +483,9 @@ def test_meshtags_from_xdmf(tmp_path, mesh):
 
 
 def test_raise_error_non_unique_vol_ids():
-    """Test that an error is raised if the volume ids are not unique"""
+    """Test that an error is raised if the volume ids are not unique."""
     my_problem = F.HeatTransferProblem()
-    my_problem.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 2000))
+    my_problem.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 11))
     left = F.SurfaceSubdomain1D(id=1, x=0)
     right = F.SurfaceSubdomain1D(id=2, x=1)
     mat = F.Material(D_0=1, E_D=None)
@@ -503,7 +504,8 @@ def test_raise_error_non_unique_vol_ids():
 
 
 def test_initial_condition():
-    """Test that an error is raised when initial conditions are defined for a transient simulation"""
+    """Test that an error is raised when initial conditions are defined for a transient
+    simulation."""
     my_problem = F.HeatTransferProblem()
 
     my_problem.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 3000))
@@ -536,7 +538,7 @@ def test_initial_condition():
 
 
 def test_adaptive_timestepping_grows():
-    """Tests that the stepsize grows"""
+    """Tests that the stepsize grows."""
     # BUILD
     my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat)
     my_model = F.HeatTransferProblem(
@@ -570,7 +572,7 @@ def test_adaptive_timestepping_grows():
 
 
 def test_adaptive_timestepping_shrinks():
-    """Tests that the stepsize shrinks"""
+    """Tests that the stepsize shrinks."""
     # BUILD
     my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat)
     my_model = F.HeatTransferProblem(
@@ -616,8 +618,8 @@ def test_adaptive_timestepping_shrinks():
     ],
 )
 def test_update_time_dependent_values_HeatFluxBC(bc_value, expected_values):
-    """Test that time dependent fluxes are updated at each time step,
-    and match an expected value"""
+    """Test that time dependent fluxes are updated at each time step, and match an
+    expected value."""
     # BUILD
     my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat)
     surface = F.SurfaceSubdomain1D(id=2, x=0)
diff --git a/test/test_helpers.py b/test/test_helpers.py
index 971069e3f..a4fff0129 100644
--- a/test/test_helpers.py
+++ b/test/test_helpers.py
@@ -26,7 +26,7 @@
     ],
 )
 def test_temperature_type_and_processing(value):
-    """Test that the temperature type is correctly set"""
+    """Test that the temperature type is correctly set."""
 
     if not isinstance(value, fem.Constant | int | float):
         with pytest.raises(TypeError):
@@ -39,7 +39,7 @@ def test_temperature_type_and_processing(value):
     "input_value, expected_output_type", [(1.0, fem.Constant), (3, fem.Constant)]
 )
 def test_value_convert_float_int_inputs(input_value, expected_output_type):
-    """Test that float and  value is correctly converted"""
+    """Test that float and  value is correctly converted."""
 
     test_value = F.Value(input_value)
 
@@ -63,7 +63,7 @@ def test_value_convert_float_int_inputs(input_value, expected_output_type):
     ],
 )
 def test_value_convert_up_to_ufl_inputs(input_value, expected_output_type):
-    """Test that float and  value is correctly converted"""
+    """Test that float and  value is correctly converted."""
 
     my_mesh = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10)
     V = fem.functionspace(my_mesh, ("Lagrange", 1))
@@ -97,7 +97,7 @@ def test_value_convert_up_to_ufl_inputs(input_value, expected_output_type):
     ],
 )
 def test_value_convert_callable_inputs(input_value, expected_output_type):
-    """Test that float and  value is correctly converted"""
+    """Test that float and  value is correctly converted."""
 
     my_mesh = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 12)
     my_t = fem.Constant(my_mesh, default_scalar_type(8))
@@ -117,13 +117,13 @@ def test_value_convert_callable_inputs(input_value, expected_output_type):
 
 
 def test_error_raised_wehn_input_value_is_not_accepted():
-    """Test that an error is raised when the input value is not accepted"""
+    """Test that an error is raised when the input value is not accepted."""
 
     with pytest.raises(
         TypeError,
         match=(
-            "Value must be a float, int, fem.Constant, np.ndarray, fem.Expression, "
-            "ufl.core.expr.Expr, fem.Function, or callable not coucou"
+            r"Value must be a float, int, fem.Constant, np.ndarray, fem.Expression, "
+            r"ufl.core.expr.Expr, fem.Function, or callable not coucou"
         ),
     ):
         F.Value("coucou")
@@ -147,7 +147,7 @@ def test_error_raised_wehn_input_value_is_not_accepted():
     ],
 )
 def test_time_dependent_values(input_value, expected_output):
-    """Test that the time_dependent attribute is correctly set"""
+    """Test that the time_dependent attribute is correctly set."""
 
     test_value = F.Value(input_value)
 
@@ -172,7 +172,7 @@ def test_time_dependent_values(input_value, expected_output):
     ],
 )
 def test_temperature_dependent_values(input_value, expected_output):
-    """Test that the time_dependent attribute is correctly set"""
+    """Test that the time_dependent attribute is correctly set."""
 
     test_value = F.Value(input_value)
 
@@ -187,7 +187,7 @@ def test_temperature_dependent_values(input_value, expected_output):
     ],
 )
 def test_input_values_of_constants_and_functions_are_accepted(value):
-    """Test that the input values of constants and functions are accepted"""
+    """Test that the input values of constants and functions are accepted."""
 
     test_value = F.Value(value)
 
@@ -197,9 +197,11 @@ def test_input_values_of_constants_and_functions_are_accepted(value):
 
 
 def test_input_values_of_expressions_are_accepted():
-    """Test that the input values of constants and functions are accepted"""
+    """Test that the input values of constants and functions are accepted."""
+
+    def my_func(x):
+        return 1.0 + x[0]
 
-    my_func = lambda x: 1.0 + x[0]
     kwargs = {}
     kwargs["x"] = x
     mapped_func = my_func(**kwargs)
@@ -217,7 +219,7 @@ def test_input_values_of_expressions_are_accepted():
 
 def test_ValueError_raised_when_callable_returns_wrong_type():
     """The create_value_fenics method should raise a ValueError when the callable
-    returns an object which is not a float or int"""
+    returns an object which is not a float or int."""
 
     def my_value(t):
         return ufl.conditional(ufl.lt(t, 0.5), 100, 0)
@@ -229,7 +231,10 @@ def my_value(t):
 
     with pytest.raises(
         ValueError,
-        match="self.value should return a float or an int, not <class 'ufl.conditional.Conditional'",
+        match=(
+            r"self.value should return a float or an int, not <class "
+            r"'ufl.conditional.Conditional'"
+        ),
     ):
         test_value.convert_input_value(
             function_space=test_function_space, temperature=T, t=t
@@ -251,7 +256,7 @@ def my_value(t):
     ],
 )
 def test_value_representation(value):
-    """Test that the representation of the value is correct"""
+    """Test that the representation of the value is correct."""
 
     test_value = F.Value(value)
 
@@ -270,7 +275,7 @@ def test_value_representation(value):
 )
 def test_velocity_field_convert_input_error_when_t_not_only_arg(value):
     """Test when an input value of type callable is converted that a Type error is
-    rasied when t is not the only arg"""
+    rasied when t is not the only arg."""
 
     test_value = F.VelocityField(value)
     t = F.as_fenics_constant(value=1.0, mesh=test_mesh.mesh)
@@ -283,7 +288,7 @@ def test_velocity_field_convert_input_error_when_t_not_only_arg(value):
 
 def test_velocity_field_convert_input_error_when_callable_doesnt_return_fem_func():
     """Test when an input value of type callable is converted that a Type error is
-    rasied when t is not the only arg"""
+    rasied when t is not the only arg."""
 
     def example_func(t):
         return 2 * t
@@ -293,7 +298,10 @@ def example_func(t):
 
     with pytest.raises(
         ValueError,
-        match="A time dependent advection field should return an fem.Function, not a <class 'ufl.algebra.Product'>",
+        match=(
+            r"A time dependent advection field should return an fem.Function, not a "
+            r"<class 'ufl.algebra.Product'>"
+        ),
     ):
         test_value.convert_input_value(function_space=test_function_space, t=t)
 
diff --git a/test/test_henrysbc.py b/test/test_henrysbc.py
index 998ea78e5..fae1a6e9f 100644
--- a/test/test_henrysbc.py
+++ b/test/test_henrysbc.py
@@ -10,13 +10,14 @@
 
 
 def henrys_law(T, H_0, E_H, pressure):
-    """Applies the Henrys law to compute the concentration at the boundary"""
+    """Applies the Henrys law to compute the concentration at the boundary."""
     H = H_0 * ufl.exp(-E_H / F.k_B / T)
     return H * pressure
 
 
 def test_raise_error():
-    """Test that a value error is raised if the pressure function is not supported in HenrysBC"""
+    """Test that a value error is raised if the pressure function is not supported in
+    HenrysBC."""
     with pytest.raises(ValueError, match="pressure function not supported"):
         F.HenrysBC(subdomain=None, H_0=1.0, E_H=1.0, pressure=lambda c: c, species="H")
 
diff --git a/test/test_initial_condition.py b/test/test_initial_condition.py
index 6ced0b703..d99defe38 100644
--- a/test/test_initial_condition.py
+++ b/test/test_initial_condition.py
@@ -15,7 +15,7 @@
 
 
 def test_init():
-    """Test that the attributes are set correctly"""
+    """Test that the attributes are set correctly."""
     # create an InitialConcentration object
     value = 1.0
     species = F.Species("test")
@@ -38,8 +38,8 @@ def test_init():
     ],
 )
 def test_create_value_fenics(input_value, expected_type):
-    """Test that after calling .create_expr_fenics, the prev_solution
-    attribute of the species has the correct value at x=1.0."""
+    """Test that after calling .create_expr_fenics, the prev_solution attribute of the
+    species has the correct value at x=1.0."""
 
     # BUILD
 
@@ -66,7 +66,7 @@ def test_create_value_fenics(input_value, expected_type):
 
 
 def test_warning_raised_when_giving_time_as_arg():
-    """Test that a warning is raised if the value is given with t in its arguments"""
+    """Test that a warning is raised if the value is given with t in its arguments."""
 
     vol = F.VolumeSubdomain(id=1, material=dummy_mat)
 
@@ -75,20 +75,21 @@ def test_warning_raised_when_giving_time_as_arg():
     my_species = F.Species("test")
     my_species.prev_solution = fem.Function(V)
 
-    my_value = lambda t: 1.0 + t
+    def my_value(t):
+        return 1.0 + t
 
     init_cond = F.InitialConcentration(value=my_value, species=my_species, volume=vol)
 
     T = fem.Constant(test_mesh.mesh, 10.0)
 
     with pytest.raises(
-        ValueError, match="Initial condition cannot be a function of time."
+        ValueError, match=r"Initial condition cannot be a function of time."
     ):
         init_cond.create_expr_fenics(test_mesh.mesh, T, V)
 
 
 def test_warning_raised_when_giving_time_as_arg_initial_temperature():
-    """Test that a warning is raised if the value is given with t in its arguments"""
+    """Test that a warning is raised if the value is given with t in its arguments."""
 
     # give function to species
     V = fem.functionspace(test_mesh.mesh, ("Lagrange", 1))
@@ -97,12 +98,13 @@ def test_warning_raised_when_giving_time_as_arg_initial_temperature():
 
     vol = F.VolumeSubdomain(id=1, material=dummy_mat)
 
-    my_value = lambda t: 1.0 + t
+    def my_value(t):
+        return 1.0 + t
 
     init_cond = F.InitialTemperature(value=my_value, volume=vol)
 
     with pytest.raises(
-        ValueError, match="Initial condition cannot be a function of time."
+        ValueError, match=r"Initial condition cannot be a function of time."
     ):
         init_cond.create_expr_fenics(test_mesh.mesh, V)
 
@@ -116,8 +118,8 @@ def test_warning_raised_when_giving_time_as_arg_initial_temperature():
     ],
 )
 def test_create_value_fenics_initial_temperature(input_value, expected_type):
-    """Test that after calling .create_expr_fenics, the prev_solution
-    attribute of the species has the correct value at x=1.0."""
+    """Test that after calling .create_expr_fenics, the prev_solution attribute of the
+    species has the correct value at x=1.0."""
 
     # BUILD
 
@@ -140,10 +142,8 @@ def test_create_value_fenics_initial_temperature(input_value, expected_type):
 
 
 def test_checkpointing_single_species(tmpdir):
-    """
-    Writes a P1 function to a file and reads it back in as initial condition
-    for one species.
-    """
+    """Writes a P1 function to a file and reads it back in as initial condition for one
+    species."""
     # build initial condition
     mesh = dolfinx.mesh.create_unit_square(
         MPI.COMM_WORLD, nx=6, ny=6, cell_type=dolfinx.cpp.mesh.CellType.quadrilateral
@@ -169,7 +169,7 @@ def f(x):
     my_problem.mesh = F.Mesh(mesh)
 
     mat = F.Material(D_0=1, E_D=0.1, name="dummy_mat")
-    vol = F.VolumeSubdomain(id=0, material=mat)
+    vol = F.VolumeSubdomain(id=1, material=mat)
     my_problem.subdomains = [vol]
 
     function_initial_value = F.read_function_from_file(
@@ -191,10 +191,8 @@ def f(x):
 
 
 def test_checkpointing_multiple_species(tmpdir):
-    """
-    Writes two P1 functions to a file and reads them back in as initial conditions
-    for two species.
-    """
+    """Writes two P1 functions to a file and reads them back in as initial conditions
+    for two species."""
     # build initial condition
     mesh = dolfinx.mesh.create_unit_square(
         MPI.COMM_WORLD, nx=6, ny=6, cell_type=dolfinx.cpp.mesh.CellType.quadrilateral
@@ -227,7 +225,7 @@ def f2(x):
     my_problem.mesh = F.Mesh(mesh)
 
     mat = F.Material(D_0=1, E_D=0.1, name="dummy_mat")
-    vol = F.VolumeSubdomain(id=0, material=mat)
+    vol = F.VolumeSubdomain(id=1, material=mat)
     my_problem.subdomains = [vol]
 
     my_problem.initial_conditions = [
@@ -272,11 +270,11 @@ def f2(x):
     ],
 )
 def test_error_raised_with_volume_setter(vol):
-    """Test that the volume setter works correctly"""
+    """Test that the volume setter works correctly."""
 
     spe = F.Species("test")
     with pytest.raises(
-        TypeError, match="volume must be of type festim.VolumeSubdomain"
+        TypeError, match=r"volume must be of type festim.VolumeSubdomain"
     ):
         F.InitialConcentration(value=1.0, species=spe, volume=vol)
 
@@ -290,21 +288,23 @@ def test_error_raised_with_volume_setter(vol):
     ],
 )
 def test_error_raised_with_species_setter(species):
-    """Test that the volume setter works correctly"""
+    """Test that the volume setter works correctly."""
 
     vol = F.VolumeSubdomain(id=1, material=dummy_mat)
-    with pytest.raises(TypeError, match="species must be of type festim.Species"):
+    with pytest.raises(TypeError, match=r"species must be of type festim.Species"):
         F.InitialConcentration(value=1.0, species=species, volume=vol)
 
 
 def test_create_initial_temperature_from_function():
-    """Test that the initial temperature can be created from a function"""
+    """Test that the initial temperature can be created from a function."""
 
     # create a volume subdomain
     vol = F.VolumeSubdomain(id=1, material=dummy_mat)
 
     # create an initial temperature from a function
-    T = lambda x: 300 + 10 * x[0]
+    def T(x):
+        return 300 + 10 * x[0]
+
     V = fem.functionspace(test_mesh.mesh, ("Lagrange", 1))
     T_func = fem.Function(V)
     T_func.interpolate(T)
@@ -320,7 +320,7 @@ def test_create_initial_temperature_from_function():
 
 def test_initial_condition_discontinuous():
     """Test the initial condition in a multispecies case with a discontinuous volume
-    subdomain"""
+    subdomain."""
 
     my_model = F.HydrogenTransportProblemDiscontinuous()
 
@@ -368,10 +368,10 @@ def test_initial_condition_discontinuous():
 
     my_model.initialise()
 
-    spe1_left, spe1_to_vol1 = vol1.u_n.function_space.sub(0).collapse()
-    spe2_left, spe2_to_vol1 = vol1.u_n.function_space.sub(1).collapse()
-    spe1_right, spe1_to_vol2 = vol2.u_n.function_space.sub(0).collapse()
-    spe2_right, spe2_to_vol2 = vol2.u_n.function_space.sub(1).collapse()
+    _spe1_left, spe1_to_vol1 = vol1.u_n.function_space.sub(0).collapse()
+    _spe2_left, spe2_to_vol1 = vol1.u_n.function_space.sub(1).collapse()
+    _spe1_right, spe1_to_vol2 = vol2.u_n.function_space.sub(0).collapse()
+    _spe2_right, spe2_to_vol2 = vol2.u_n.function_space.sub(1).collapse()
 
     prev_solution_spe1_left = vol1.u_n.x.array[spe1_to_vol1]
     prev_solution_spe2_left = vol1.u_n.x.array[spe2_to_vol1]
@@ -389,7 +389,7 @@ def test_initial_condition_discontinuous():
 
 def test_initial_condition_continuous_multimaterial():
     """Test the initial condition in multi-material continous case that the condition is
-    only appilied in the correct volume subdomain"""
+    only appilied in the correct volume subdomain."""
 
     my_model = F.HydrogenTransportProblem()
 
@@ -433,7 +433,7 @@ def test_initial_condition_continuous_multimaterial():
 
 
 def test_initial_condition_mixed_domain():
-    """Test the initial condition in a multi-material discontinous case"""
+    """Test the initial condition in a multi-material discontinous case."""
 
     my_model = F.HydrogenTransportProblemDiscontinuous()
 
@@ -478,7 +478,7 @@ def test_initial_condition_mixed_domain():
 
 
 def test_initial_condition_mixed_domain_multispecies():
-    """Test the initial condition in a multispecies multi-material discontinous case"""
+    """Test the initial condition in a multispecies multi-material discontinous case."""
 
     my_model = F.HydrogenTransportProblemDiscontinuous()
 
@@ -523,8 +523,8 @@ def test_initial_condition_mixed_domain_multispecies():
 
     my_model.initialise()
 
-    spe1_left, spe1_to_vol1 = vol1.u_n.function_space.sub(0).collapse()
-    spe2_left, spe2_to_vol1 = vol1.u_n.function_space.sub(1).collapse()
+    _spe1_left, spe1_to_vol1 = vol1.u_n.function_space.sub(0).collapse()
+    _spe2_left, spe2_to_vol1 = vol1.u_n.function_space.sub(1).collapse()
 
     prev_solution_spe1_left = vol1.u_n.x.array[spe1_to_vol1]
     prev_solution_spe2_left = vol1.u_n.x.array[spe2_to_vol1]
diff --git a/test/test_material.py b/test/test_material.py
index 7c425b558..b6d166344 100644
--- a/test/test_material.py
+++ b/test/test_material.py
@@ -8,7 +8,7 @@
 
 
 def test_define_diffusion_coefficient():
-    """Test that the diffusion coefficient is correctly defined"""
+    """Test that the diffusion coefficient is correctly defined."""
     T, D_0, E_D = 10, 1.2, 0.5
     dum_spe = F.Species("dummy")
     my_mat = F.Material(D_0=D_0, E_D=E_D)
@@ -21,7 +21,7 @@ def test_define_diffusion_coefficient():
 
 def test_multispecies_dict_strings():
     """Test that the diffusion coefficient is correctly defined when keys are
-    strings"""
+    strings."""
     T = 500
     D_0_A, D_0_B = 1, 2
     E_D_A, E_D_B = 0.1, 0.2
@@ -43,7 +43,7 @@ def test_multispecies_dict_strings():
 
 def test_multispecies_dict_objects():
     """Test that the diffusion coefficient is correctly defined when keys are
-    festim.Species objects"""
+    festim.Species objects."""
     T = 500
     D_0_A, D_0_B = 1, 2
     E_D_A, E_D_B = 0.1, 0.2
@@ -65,8 +65,8 @@ def test_multispecies_dict_objects():
 
 
 def test_multispecies_dict_objects_and_strings():
-    """Test that the diffusion coefficient is correctly defined when keys
-    are a mix of festim.Species objects and strings"""
+    """Test that the diffusion coefficient is correctly defined when keys are a mix of
+    festim.Species objects and strings."""
     T = 500
     D_0_A, D_0_B = 1, 2
     E_D_A, E_D_B = 0.1, 0.2
@@ -88,8 +88,8 @@ def test_multispecies_dict_objects_and_strings():
 
 
 def test_multispecies_dict_different_keys():
-    """Test that a value error is raised when the keys of the D_0 and E_D
-    are not the same"""
+    """Test that a value error is raised when the keys of the D_0 and E_D are not the
+    same."""
     A = F.Species("A")
     my_mat = F.Material(D_0={"A": 1, "B": 2}, E_D={"A": 0.1, "B": 0.2, "C": 0.3})
 
@@ -98,8 +98,7 @@ def test_multispecies_dict_different_keys():
 
 
 def test_D_0_type_raises_error():
-    """Test that a value error is raised in the get_diffusion_coefficient
-    function"""
+    """Test that a value error is raised in the get_diffusion_coefficient function."""
     # TODO remove this when material class is updated
     A = F.Species("A")
     my_mat = F.Material(D_0=[1, 1], E_D=0.1)
@@ -109,8 +108,8 @@ def test_D_0_type_raises_error():
 
 
 def test_error_raised_when_species_not_given_with_dict():
-    """Test that a value error is raised when a species has not been given in
-    the get_diffusion_coefficient function when using a dict for properties"""
+    """Test that a value error is raised when a species has not been given in the
+    get_diffusion_coefficient function when using a dict for properties."""
     A = F.Species("A")
     B = F.Species("B")
     my_mat = F.Material(D_0={A: 1, B: 1}, E_D={A: 0.1, B: 0.1})
@@ -122,8 +121,8 @@ def test_error_raised_when_species_not_given_with_dict():
 
 
 def test_error_raised_when_species_not_not_in_D_0_dict():
-    """Test that a value error is raised when a species has not been given but
-    has no value in the dict"""
+    """Test that a value error is raised when a species has not been given but has no
+    value in the dict."""
     A = F.Species("A")
     B = F.Species("B")
     J = F.Species("J")
@@ -134,8 +133,7 @@ def test_error_raised_when_species_not_not_in_D_0_dict():
 
 
 def test_D_0_raises_ValueError_if_species_not_provided_in_dict():
-    """Test that a value error is raised in the get_diffusion_coefficient
-    function"""
+    """Test that a value error is raised in the get_diffusion_coefficient function."""
     # TODO remove this when material class is updated
     A = F.Species("A")
     B = F.Species("B")
@@ -146,8 +144,7 @@ def test_D_0_raises_ValueError_if_species_not_provided_in_dict():
 
 
 def test_D_0_raises_ValueError_if_species_given_not_in_dict_keys():
-    """Test that a value error is raised in the get_diffusion_coefficient
-    function"""
+    """Test that a value error is raised in the get_diffusion_coefficient function."""
     # TODO remove this when material class is updated
     A = F.Species("A")
     B = F.Species("B")
@@ -159,7 +156,7 @@ def test_D_0_raises_ValueError_if_species_given_not_in_dict_keys():
 
 
 def test_raises_TypeError_when_D_0_is_not_correct_type():
-    """Test that a TypeError is raised when D_0 is not a float or a dict"""
+    """Test that a TypeError is raised when D_0 is not a float or a dict."""
 
     my_mat = F.Material(D_0=[1, 2], E_D=1)
 
@@ -168,8 +165,7 @@ def test_raises_TypeError_when_D_0_is_not_correct_type():
 
 
 def test_E_D_raises_ValueError_if_species_not_provided_in_dict():
-    """Test that a value error is raised in the get_diffusion_coefficient
-    function"""
+    """Test that a value error is raised in the get_diffusion_coefficient function."""
     # TODO remove this when material class is updated
     A = F.Species("A")
     B = F.Species("B")
@@ -180,8 +176,7 @@ def test_E_D_raises_ValueError_if_species_not_provided_in_dict():
 
 
 def test_E_D_raises_ValueError_if_species_given_not_in_dict_keys():
-    """Test that a value error is raised in the get_diffusion_coefficient
-    function"""
+    """Test that a value error is raised in the get_diffusion_coefficient function."""
     # TODO remove this when material class is updated
     A = F.Species("A")
     B = F.Species("B")
@@ -193,7 +188,7 @@ def test_E_D_raises_ValueError_if_species_given_not_in_dict_keys():
 
 
 def test_raises_TypeError_when_E_D_is_not_correct_type():
-    """Test that a TypeError is raised when E_D is not a float or a dict"""
+    """Test that a TypeError is raised when E_D is not a float or a dict."""
 
     my_mat = F.Material(D_0=1, E_D=[1, 2])
 
@@ -211,15 +206,14 @@ def test_raises_TypeError_when_E_D_is_not_correct_type():
     ],
 )
 def test_raises_TypeError_when_D_is_not_correct_type(input_value):
-    """Test that a TypeError is raised when D is not an fem.Function"""
+    """Test that a TypeError is raised when D is not an fem.Function."""
 
-    with pytest.raises(TypeError, match="D must be of type fem.Function"):
+    with pytest.raises(TypeError, match=r"D must be of type fem.Function"):
         F.Material(D=input_value)
 
 
 def test_get_diffusion_coefficient_returns_function_when_given_to_D():
-    """Test that the diffusion coefficient is correctly defined when D is a
-    function"""
+    """Test that the diffusion coefficient is correctly defined when D is a function."""
 
     V = fem.functionspace(test_mesh.mesh, ("Lagrange", 1))
     D = fem.Function(V)
@@ -232,23 +226,23 @@ def test_get_diffusion_coefficient_returns_function_when_given_to_D():
 
 
 def test_error_raised_when_D_and_D_0_given():
-    """Test that a value error is raised when both D and D_0 are given"""
+    """Test that a value error is raised when both D and D_0 are given."""
 
     V = fem.functionspace(test_mesh.mesh, ("Lagrange", 1))
     D_func = fem.Function(V)
 
     with pytest.raises(
         ValueError,
-        match="D_0 and D cannot be set at the same time. Please set only one of them.",
+        match=r"D_0 and D cannot be set at the same time. Please set only one of them.",
     ):
         F.Material(D=D_func, D_0=1)
 
 
 def test_error_raised_when_D_and_D_0_both_None():
-    """Test that a value error is raised when both D and D_0 are None"""
+    """Test that a value error is raised when both D and D_0 are None."""
 
     with pytest.raises(
         ValueError,
-        match="D_0 and D cannot both be None. Please set one of them.",
+        match=r"D_0 and D cannot both be None. Please set one of them.",
     ):
         F.Material()
diff --git a/test/test_maximum_surface.py b/test/test_maximum_surface.py
index 54eb62de2..096c06b70 100644
--- a/test/test_maximum_surface.py
+++ b/test/test_maximum_surface.py
@@ -5,7 +5,7 @@
 
 
 def test_maximum_surface_compute_1D():
-    """Test that the maximum surface export computes the correct value"""
+    """Test that the maximum surface export computes the correct value."""
 
     # BUILD
     L = 4.0
@@ -14,7 +14,7 @@ def test_maximum_surface_compute_1D():
     dummy_surface = F.SurfaceSubdomain1D(id=1, x=4)
 
     # create mesh tags
-    ft, ct = my_mesh.define_meshtags(
+    ft, _ct = my_mesh.define_meshtags(
         surface_subdomains=[dummy_surface],
         volume_subdomains=[
             F.VolumeSubdomain1D(id=1, material=F.Material(D_0=1, E_D=0), borders=[0, L])
diff --git a/test/test_maximum_volume.py b/test/test_maximum_volume.py
index ea3719d82..ef47c1d68 100644
--- a/test/test_maximum_volume.py
+++ b/test/test_maximum_volume.py
@@ -5,7 +5,7 @@
 
 
 def test_maximum_volume_compute_1D():
-    """Test that the maximum volume export computes the right value"""
+    """Test that the maximum volume export computes the right value."""
 
     # BUILD
     L = 6
@@ -14,7 +14,7 @@ def test_maximum_volume_compute_1D():
     dummy_volume = F.VolumeSubdomain1D(id=1, borders=[0, L], material=dummy_material)
 
     # create mesh tags
-    ft, ct = my_mesh.define_meshtags(
+    _ft, ct = my_mesh.define_meshtags(
         surface_subdomains=[F.SurfaceSubdomain1D(id=1, x=0)],
         volume_subdomains=[dummy_volume],
         interfaces=None,
diff --git a/test/test_mesh.py b/test/test_mesh.py
index f00872ef6..a087149f2 100644
--- a/test/test_mesh.py
+++ b/test/test_mesh.py
@@ -1,3 +1,4 @@
+# ruff: noqa: E402
 import logging
 import os
 
@@ -76,7 +77,8 @@ def test_vdim_changes_when_mesh_changes():
 
 @pytest.mark.parametrize("mesh", [mesh_1D, mesh_2D, mesh_3D])
 def test_meshtags_from_xdmf(tmp_path, mesh):
-    """Test that the facet and volume meshtags are read correctly from the mesh XDMF files"""
+    """Test that the facet and volume meshtags are read correctly from the mesh XDMF
+    files."""
     # create mesh functions
     fdim = mesh.topology.dim - 1
     vdim = mesh.topology.dim
@@ -160,23 +162,24 @@ def test_meshtags_from_xdmf(tmp_path, mesh):
 
 @pytest.mark.parametrize("vertices", [[1, 2, 3, 4], [0, 0.1, 0.2, 0.3, 0.4, 0.5]])
 def test_mesh_vertices_from_list(vertices):
-    """Check that giving vertices as a list is correctly processed and ends up as a np.ndarray for the mesh"""
+    """Check that giving vertices as a list is correctly processed and ends up as a
+    np.ndarray for the mesh."""
     my_mesh = F.Mesh1D(vertices=vertices)
 
     assert isinstance(my_mesh.vertices, np.ndarray)
 
 
 def test_error_raised_when_mesh_is_wrong_type():
-    """Test that an TypeError is raised when the mesh is not a dolfinx mesh"""
+    """Test that an TypeError is raised when the mesh is not a dolfinx mesh."""
 
-    with pytest.raises(TypeError, match="Mesh must be of type dolfinx.mesh.Mesh"):
+    with pytest.raises(TypeError, match=r"Mesh must be of type dolfinx.mesh.Mesh"):
         F.Mesh(
             mesh="mesh",
         )
 
 
 def test_create_1D_mesh_parallel(cluster):
-    """Test creating a 1D mesh in parallel using ipyparallel"""
+    """Test creating a 1D mesh in parallel using ipyparallel."""
 
     def create_mesh():
         import numpy as np
@@ -193,7 +196,7 @@ def create_mesh():
 @pytest.mark.parametrize("mesh", [mesh_2D, mesh_3D])
 def test_attr_error_with_incompitable_mesh_in_spherical(mesh):
     """Test that an AttributeError is raised when trying to use an incompatible mesh
-    with a spherical coordinate system"""
+    with a spherical coordinate system."""
 
     with pytest.raises(
         AttributeError,
@@ -204,7 +207,7 @@ def test_attr_error_with_incompitable_mesh_in_spherical(mesh):
 
 def test_attr_error_with_3D_mesh_in_cylindrical():
     """Test that an AttributeError is raised when trying to use an incompatible mesh
-    with a cylindrical coordinate system"""
+    with a cylindrical coordinate system."""
 
     with pytest.raises(
         AttributeError,
@@ -216,7 +219,7 @@ def test_attr_error_with_3D_mesh_in_cylindrical():
 @pytest.mark.parametrize("system", ["cyl", "cart", 1.0, "coucou", mesh_1D])
 def test_coordinate_system_setter(system):
     if isinstance(system, str):
-        err_msg = "coordinate_system must be one of 'cartesian', 'cylindrical', or 'spherical'"
+        err_msg = "coordinate_system must be one of 'cartesian', 'cylindrical', or 'spherical'"  # noqa: E501
     else:
         err_msg = "coordinate_system must be of type str or CoordinateSystem"
     with pytest.raises(
diff --git a/test/test_minimum_surface.py b/test/test_minimum_surface.py
index b6a64e9aa..416481819 100644
--- a/test/test_minimum_surface.py
+++ b/test/test_minimum_surface.py
@@ -5,7 +5,7 @@
 
 
 def test_minimum_surface_export_compute_1D():
-    """Test that the minimum surface export computes the correct value"""
+    """Test that the minimum surface export computes the correct value."""
 
     # BUILD
     L = 4.0
@@ -14,7 +14,7 @@ def test_minimum_surface_export_compute_1D():
     dummy_surface = F.SurfaceSubdomain1D(id=1, x=4)
 
     # create mesh tags
-    ft, ct = my_mesh.define_meshtags(
+    ft, _ct = my_mesh.define_meshtags(
         surface_subdomains=[dummy_surface],
         volume_subdomains=[
             F.VolumeSubdomain1D(id=1, material=F.Material(D_0=1, E_D=0), borders=[0, L])
diff --git a/test/test_minimum_volume.py b/test/test_minimum_volume.py
index 820cab40f..a6e42ccbd 100644
--- a/test/test_minimum_volume.py
+++ b/test/test_minimum_volume.py
@@ -5,7 +5,7 @@
 
 
 def test_minimum_volume_compute_1D():
-    """Test that the minimum volume export computes the right value"""
+    """Test that the minimum volume export computes the right value."""
 
     # BUILD
     L = 6
@@ -14,7 +14,7 @@ def test_minimum_volume_compute_1D():
     dummy_volume = F.VolumeSubdomain1D(id=1, borders=[0, L], material=dummy_material)
 
     # create mesh tags
-    ft, ct = my_mesh.define_meshtags(
+    _ft, ct = my_mesh.define_meshtags(
         surface_subdomains=[F.SurfaceSubdomain1D(id=1, x=0)],
         volume_subdomains=[dummy_volume],
         interfaces=None,
diff --git a/test/test_multispecies_problem.py b/test/test_multispecies_problem.py
index 4defae394..6529e81f8 100644
--- a/test/test_multispecies_problem.py
+++ b/test/test_multispecies_problem.py
@@ -33,9 +33,9 @@ def relative_error_computed_to_analytical(
 
 
 def test_multispecies_permeation_problem():
-    """Test running a problem with 2 mobile species permeating through a 1D
-    domain, with different diffusion coefficients, checks that the computed
-    permeation flux matches the analytical solution"""
+    """Test running a problem with 2 mobile species permeating through a 1D domain, with
+    different diffusion coefficients, checks that the computed permeation flux matches
+    the analytical solution."""
 
     # festim model
     L = 3e-04
diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py
index d5da1e391..011821539 100644
--- a/test/test_permeation_problem.py
+++ b/test/test_permeation_problem.py
@@ -36,9 +36,8 @@ def relative_error_computed_to_analytical(
 
 
 def test_permeation_problem(mesh_size=1001):
-    """Test running a problem with a mobile species permeating through a 1D
-    domain, checks that the computed permeation flux matches the analytical
-    solution"""
+    """Test running a problem with a mobile species permeating through a 1D domain,
+    checks that the computed permeation flux matches the analytical solution."""
 
     # festim model
     L = 3e-04
@@ -107,7 +106,7 @@ def test_permeation_problem(mesh_size=1001):
     times = outgassing_flux.t
     flux_values = outgassing_flux.data
 
-    # -------------------------- analytical solution -------------------------------------
+    # -------------------------- analytical solution -------------------------------------  # noqa: E501
 
     D = my_mat.get_diffusion_coefficient(my_mesh.mesh, my_model.temperature)
 
@@ -127,8 +126,7 @@ def test_permeation_problem(mesh_size=1001):
 
 
 def test_permeation_problem_multi_volume(tmp_path):
-    """Same permeation problem as above but with 4 volume subdomains instead
-    of 1"""
+    """Same permeation problem as above but with 4 volume subdomains instead of 1."""
 
     L = 3e-04
     vertices = np.linspace(0, L, num=1001)
diff --git a/test/test_reaction.py b/test/test_reaction.py
index a34861763..850214a58 100644
--- a/test/test_reaction.py
+++ b/test/test_reaction.py
@@ -1,10 +1,10 @@
 from mpi4py import MPI
 
+import numpy as np
 import pytest
 from dolfinx.fem import Function, functionspace
 from dolfinx.mesh import create_unit_cube
 from ufl import exp
-import numpy as np
 
 import festim as F
 
@@ -12,7 +12,7 @@
 
 
 def test_reaction_init():
-    """Test that the Reaction class initialises correctly"""
+    """Test that the Reaction class initialises correctly."""
     # create two species
     species1 = F.Species("A")
     species2 = F.Species("B")
@@ -41,7 +41,7 @@ def test_reaction_init():
 
 
 def test_reaction_repr():
-    """Test that the Reaction __repr__ method returns the expected string"""
+    """Test that the Reaction __repr__ method returns the expected string."""
 
     # create two species
     species1 = F.Species("A")
@@ -67,7 +67,7 @@ def test_reaction_repr():
 
 
 def test_reaction_repr_2_products():
-    """Test that the Reaction __repr__ method returns the expected string"""
+    """Test that the Reaction __repr__ method returns the expected string."""
 
     # create two species
     species1 = F.Species("A")
@@ -94,7 +94,7 @@ def test_reaction_repr_2_products():
 
 
 def test_reaction_repr_0_products():
-    """Test that the Reaction __repr__ method returns the expected string"""
+    """Test that the Reaction __repr__ method returns the expected string."""
 
     # create two species
     species1 = F.Species("A")
@@ -113,7 +113,7 @@ def test_reaction_repr_0_products():
 
 
 def test_reaction_str():
-    """Test that the Reaction __str__ method returns the expected string"""
+    """Test that the Reaction __str__ method returns the expected string."""
 
     # create two species
     species1 = F.Species("A")
@@ -139,7 +139,8 @@ def test_reaction_str():
 
 
 def test_reaction_str_2_products():
-    """Test that the Reaction __str__ method returns the expected string when there are 2 products"""
+    """Test that the Reaction __str__ method returns the expected string when there are
+    2 products."""
 
     # create two species
     species1 = F.Species("A")
@@ -166,7 +167,8 @@ def test_reaction_str_2_products():
 
 
 def test_reaction_str_no_products():
-    """Test that the Reaction __str__ method returns the expected string when there are 2 products"""
+    """Test that the Reaction __str__ method returns the expected string when there are
+    2 products."""
 
     # create two species
     species1 = F.Species("A")
@@ -186,7 +188,8 @@ def test_reaction_str_no_products():
 
 @pytest.mark.parametrize("temperature", [300.0, 350, 370, 500.0])
 def test_reaction_reaction_term(temperature):
-    """Test that the Reaction.reaction_term method returns the expected reaction term"""
+    """Test that the Reaction.reaction_term method returns the expected reaction
+    term."""
 
     mesh = create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)
     V = functionspace(mesh, ("Lagrange", 1))
@@ -258,7 +261,8 @@ def arrhenius(pre, act, T):
 
 @pytest.mark.parametrize("temperature", [300.0, 350, 370, 500.0])
 def test_reaction_reaction_term_2_products(temperature):
-    """Test that the Reaction.reaction_term method returns the expected reaction term with two products"""
+    """Test that the Reaction.reaction_term method returns the expected reaction term
+    with two products."""
 
     mesh = create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)
     V = functionspace(mesh, ("Lagrange", 1))
@@ -304,7 +308,10 @@ def test_reactant_setter_raises_error_with_zero_length_list():
     """Test a value error is raised when the first reactant is given a wrong type."""
     with pytest.raises(
         ValueError,
-        match="reactant must be an entry of one or more species objects, not an empty list.",
+        match=(
+            r"reactant must be an entry of one or more species objects, not an empty "
+            r"list."
+        ),
     ):
         F.Reaction(
             reactant=[],
@@ -320,7 +327,7 @@ def test_reactant_setter_raises_error_with_wrong_type():
     """Test a type error is raised when the first reactant is given a wrong type."""
     with pytest.raises(
         TypeError,
-        match="reactant must be an F.Species or F.ImplicitSpecies, not <class 'str'>",
+        match=r"reactant must be an F.Species or F.ImplicitSpecies, not <class 'str'>",
     ):
         F.Reaction(
             reactant=["A", F.Species("B")],
@@ -336,7 +343,7 @@ def test_reactant_setter_raises_error_with_wrong_type():
 def test_product_setter_raise_error_p_0_no_product():
     with pytest.raises(
         ValueError,
-        match="p_0 must be None, not 2 when no products are present.",
+        match=r"p_0 must be None, not 2 when no products are present.",
     ):
         reaction = F.Reaction(
             reactant=[F.Species("A")],
@@ -351,7 +358,7 @@ def test_product_setter_raise_error_p_0_no_product():
 def test_no_E_p_with_product():
     with pytest.raises(
         ValueError,
-        match="E_p cannot be None when reaction products are present.",
+        match=r"E_p cannot be None when reaction products are present.",
     ):
         reaction = F.Reaction(
             reactant=[F.Species("A")],
@@ -367,7 +374,7 @@ def test_no_E_p_with_product():
 def test_no_p_0_with_product():
     with pytest.raises(
         ValueError,
-        match="p_0 cannot be None when reaction products are present.",
+        match=r"p_0 cannot be None when reaction products are present.",
     ):
         reaction = F.Reaction(
             reactant=[F.Species("A")],
@@ -383,7 +390,7 @@ def test_no_p_0_with_product():
 def test_product_setter_raise_error_E_p_no_product():
     with pytest.raises(
         ValueError,
-        match="E_p must be None, not 2 when no products are present.",
+        match=r"E_p must be None, not 2 when no products are present.",
     ):
         reaction = F.Reaction(
             reactant=[F.Species("A")],
@@ -401,7 +408,7 @@ def test_product_setter_raise_error_E_p_no_product():
 vol1 = F.VolumeSubdomain1D(id=1, borders=[0, 0.5], material=mat)
 vol2 = F.VolumeSubdomain1D(id=2, borders=[0.5, 1], material=mat)
 my_model = F.HydrogenTransportProblemDiscontinuous()
-my_model.mesh = F.Mesh1D(np.linspace(0, 1, 100))
+my_model.mesh = F.Mesh1D(np.linspace(0, 1, 101))
 
 my_model.subdomains = [vol1, vol2]
 
@@ -430,11 +437,9 @@ def test_product_setter_raise_error_E_p_no_product():
 
 @pytest.mark.parametrize("reaction", [reac1, reac2])
 def test_override_solution_attributes(reaction):
-    """
-    Tests the HydrogenTransportProblemDiscontinuous.override_solution_attributes method
-    Checks that the .solution attribute is the expected one based on the volume of the
-    reaction
-    """
+    """Tests the HydrogenTransportProblemDiscontinuous.override_solution_attributes
+    method Checks that the .solution attribute is the expected one based on the volume
+    of the reaction."""
 
     # RUN
     my_model.override_solution_attributes(reaction)
diff --git a/test/test_settings.py b/test/test_settings.py
index 1a2275420..4d6a0bbeb 100644
--- a/test/test_settings.py
+++ b/test/test_settings.py
@@ -6,7 +6,7 @@
 
 @pytest.mark.parametrize("test_type", [int, F.Stepsize, float])
 def test_stepsize_value(test_type):
-    """Test that the stepsize is correctly set"""
+    """Test that the stepsize is correctly set."""
     test_value = 23.0
     my_settings = F.Settings(atol=1, rtol=0.1)
     my_settings.stepsize = test_type(test_value)
@@ -16,7 +16,7 @@ def test_stepsize_value(test_type):
 
 
 def test_stepsize_value_wrong_type():
-    """Checks that an error is raised when the wrong type is given"""
+    """Checks that an error is raised when the wrong type is given."""
     my_settings = F.Settings(atol=1, rtol=0.1)
 
     with pytest.raises(TypeError):
@@ -49,7 +49,7 @@ def test_stepsize_value_wrong_type():
 #     ],
 # )
 # def test_tolerances_value(rtol, atol):
-#     """Tests that callable tolerances are called & return correct float before passed to fenics"""
+#     """Tests that callable tolerances are called & return correct float before passed to fenics"""  # noqa: E501
 
 #     # BUILD
 #     test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0]))
diff --git a/test/test_sievertsbc.py b/test/test_sievertsbc.py
index a70d697a5..992104887 100644
--- a/test/test_sievertsbc.py
+++ b/test/test_sievertsbc.py
@@ -10,13 +10,14 @@
 
 
 def sieverts_law(T, S_0, E_S, pressure):
-    """Applies the Sieverts law to compute the concentration at the boundary"""
+    """Applies the Sieverts law to compute the concentration at the boundary."""
     S = S_0 * ufl.exp(-E_S / F.k_B / T)
     return S * pressure**0.5
 
 
 def test_raise_error():
-    """Test that a value error is raised if the pressure function is not supported in SievertsBC"""
+    """Test that a value error is raised if the pressure function is not supported in
+    SievertsBC."""
     with pytest.raises(ValueError, match="pressure function not supported"):
         F.SievertsBC(
             subdomain=None, S_0=1.0, E_S=1.0, pressure=lambda c: c, species="H"
diff --git a/test/test_source.py b/test/test_source.py
index 676d346a9..7290e9998 100644
--- a/test/test_source.py
+++ b/test/test_source.py
@@ -13,7 +13,7 @@
 
 
 def test_init():
-    """Test that the attributes are set correctly"""
+    """Test that the attributes are set correctly."""
     # create a Source object
     volume = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=dummy_mat)
     value = 1.0
@@ -46,7 +46,7 @@ def test_init():
     ],
 )
 def test_create_fenics_object(value, expected_type):
-    """Test that the correct fenics object is created depending on the value input"""
+    """Test that the correct fenics object is created depending on the value input."""
 
     # BUILD
     vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
@@ -83,7 +83,7 @@ def test_create_fenics_object(value, expected_type):
     ],
 )
 def test_source_explicit_time_dependent_attribute(input, expected_value):
-    """Test that the time_dependent attribute is correctly set"""
+    """Test that the time_dependent attribute is correctly set."""
     volume = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
     species = F.Species("test")
     my_source = F.ParticleSource(input, volume, species)
@@ -108,7 +108,7 @@ def test_source_explicit_time_dependent_attribute(input, expected_value):
     ],
 )
 def test_source_temperature_dependent_attribute(input, expected_value):
-    """Test that the temperature_dependent attribute is correctly set"""
+    """Test that the temperature_dependent attribute is correctly set."""
     volume = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
     species = F.Species("test")
     my_source = F.ParticleSource(input, volume, species)
@@ -117,8 +117,8 @@ def test_source_temperature_dependent_attribute(input, expected_value):
 
 
 def test_ValueError_raised_when_callable_returns_wrong_type():
-    """The create_value method should raise a ValueError when the callable
-    returns an object which is not a float or int"""
+    """The create_value method should raise a ValueError when the callable returns an
+    object which is not a float or int."""
 
     vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
     species = F.Species("test")
@@ -135,7 +135,10 @@ def my_value(t):
 
     with pytest.raises(
         ValueError,
-        match="self.value should return a float or an int, not <class 'ufl.conditional.Conditional'",
+        match=(
+            r"self.value should return a float or an int, not <class "
+            r"'ufl.conditional.Conditional'"
+        ),
     ):
         source.value.convert_input_value(
             function_space=V, temperature=T, t=t, up_to_ufl_expr=True
@@ -143,8 +146,8 @@ def my_value(t):
 
 
 def test_ValueError_raised_when_callable_returns_wrong_type_heat_source():
-    """The create_value method should raise a ValueError when the callable
-    returns an object which is not a float or int"""
+    """The create_value method should raise a ValueError when the callable returns an
+    object which is not a float or int."""
 
     vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
 
@@ -159,7 +162,10 @@ def my_value(t):
 
     with pytest.raises(
         ValueError,
-        match="self.value should return a float or an int, not <class 'ufl.conditional.Conditional'",
+        match=(
+            r"self.value should return a float or an int, not <class "
+            r"'ufl.conditional.Conditional'"
+        ),
     ):
         source.value.convert_input_value(function_space=V, t=t, up_to_ufl_expr=True)
 
@@ -179,12 +185,12 @@ def my_value(t):
 )
 def test_TypeError_is_raised_when_volume_wrong_type(volume_input):
     """Test that a TypeError is raised when the volume is not of type
-    festim.VolumeSubdomain"""
+    festim.VolumeSubdomain."""
 
     my_spe = F.Species("test")
     with pytest.raises(
         TypeError,
-        match="volume must be of type festim.VolumeSubdomain",
+        match=r"volume must be of type festim.VolumeSubdomain",
     ):
         F.ParticleSource(volume=volume_input, value=1.0, species=my_spe)
 
@@ -203,12 +209,12 @@ def test_TypeError_is_raised_when_volume_wrong_type(volume_input):
 )
 def test_TypeError_is_raised_when_species_wrong_type(species_input):
     """Test that a TypeError is raised when the species is not of type
-    festim.Species"""
+    festim.Species."""
 
     my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=dummy_mat)
 
     with pytest.raises(
         TypeError,
-        match="species must be of type festim.Species",
+        match=r"species must be of type festim.Species",
     ):
         F.ParticleSource(volume=my_vol, value=1.0, species=species_input)
diff --git a/test/test_species.py b/test/test_species.py
index e5d2af306..dc2f4d7ea 100644
--- a/test/test_species.py
+++ b/test/test_species.py
@@ -10,9 +10,8 @@
 
 
 def test_assign_functions_to_species():
-    """Test that checks if the function assign_functions_to_species
-    creates the correct attributes for each species
-    """
+    """Test that checks if the function assign_functions_to_species creates the correct
+    attributes for each species."""
 
     mesh = F.Mesh1D(
         vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0])
@@ -38,8 +37,7 @@ def test_assign_functions_to_species():
 
 def test_species_repr_and_str():
     """Test that the __repr__ and __str__ methods of the Species class returns the
-    expected string.
-    """
+    expected string."""
     # create a species
     species = F.Species("A")
 
@@ -53,9 +51,8 @@ def test_species_repr_and_str():
 
 
 def test_implicit_species_repr_and_str():
-    """Test that the __repr__ and __str__ methods of the ImplicitSpecies class
-    returns the expected string.
-    """
+    """Test that the __repr__ and __str__ methods of the ImplicitSpecies class returns
+    the expected string."""
     # create two species
     species1 = F.Species("A")
     species2 = F.Species("B")
@@ -73,9 +70,7 @@ def test_implicit_species_repr_and_str():
 
 
 def test_implicit_species_concentration():
-    """Test that the concentration of an implicit species is computed
-    correctly.
-    """
+    """Test that the concentration of an implicit species is computed correctly."""
     # create two species
     species1 = F.Species("A")
     species2 = F.Species("B")
@@ -99,10 +94,8 @@ def test_implicit_species_concentration():
 
 
 def test_implicit_species_concentration_with_no_solution():
-    """Test that a ValueError is raised when on of the 'others' species
-    has no solution and the concentration of the implicit species is
-    requested.
-    """
+    """Test that a ValueError is raised when on of the 'others' species has no solution
+    and the concentration of the implicit species is requested."""
     # create two species
     species1 = F.Species("A")
     species2 = F.Species("B")
@@ -121,8 +114,8 @@ def test_implicit_species_concentration_with_no_solution():
 
 
 def test_create_species_and_reaction():
-    """test that the trapped_concentration and trap_reaction attributes
-    are correctly set"""
+    """Test that the trapped_concentration and trap_reaction attributes are correctly
+    set."""
 
     # BUILD
     my_mobile_species = F.Species("test_mobile")
@@ -147,15 +140,16 @@ def test_create_species_and_reaction():
 
 
 def test_implicit_species_wrong_type():
-    """Test that a ValueError is raised when the value of n in an ImplicitSpecies is
-    a function of time only but doesn't return a float or an int.
-    """
+    """Test that a ValueError is raised when the value of n in an ImplicitSpecies is a
+    function of time only but doesn't return a float or an int."""
     mesh = create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)
 
-    density = lambda t: "coucou"
+    def density(t):
+        return "coucou"
+
     species = F.ImplicitSpecies(n=density)
 
     with pytest.raises(
-        ValueError, match="self.value should return a float or an int, not "
+        ValueError, match=r"self.value should return a float or an int, not "
     ):
         species.create_value_fenics(mesh=mesh, t=0.0)
diff --git a/test/test_stepsize.py b/test/test_stepsize.py
index 835c1fd79..a8e258ff9 100644
--- a/test/test_stepsize.py
+++ b/test/test_stepsize.py
@@ -6,7 +6,7 @@
 
 @pytest.mark.parametrize("growth_factor, target", [(10, 5), (1.2, 2), (1, 1)])
 def test_adaptive_stepsize_grows(growth_factor, target):
-    """Checks that the stepsize is increased correctly
+    """Checks that the stepsize is increased correctly.
 
     Args:
         growth_factor (float): the growth factor
@@ -28,7 +28,7 @@ def test_adaptive_stepsize_grows(growth_factor, target):
 
 @pytest.mark.parametrize("cutback_factor, target", [(0.8, 5), (0.5, 2), (1, 1)])
 def test_adaptive_stepsize_shrinks(cutback_factor, target):
-    """Checks that the stepsize is shrinks correctly
+    """Checks that the stepsize is shrinks correctly.
 
     Args:
         cutback_factor (float): the cutback factor
@@ -50,8 +50,7 @@ def test_adaptive_stepsize_shrinks(cutback_factor, target):
 
 @pytest.mark.parametrize("nb_its, target", [(1, 4), (5, 4), (4, 4)])
 def test_max_stepsize(nb_its, target):
-    """Checks that the stepsize is capped at max
-    stepsize.
+    """Checks that the stepsize is capped at max stepsize.
 
     Args:
         nb_its (int): the current number of iterations
@@ -75,10 +74,7 @@ def test_max_stepsize(nb_its, target):
 
 
 def test_stepsize_is_unchanged():
-    """
-    Checks that the stepsize is unchanged when reaches
-    the target nb iterations
-    """
+    """Checks that the stepsize is unchanged when reaches the target nb iterations."""
     my_stepsize = F.Stepsize(initial_value=2)
     my_stepsize.target_nb_iterations = 5
 
@@ -89,7 +85,7 @@ def test_stepsize_is_unchanged():
 
 
 def test_custom_stepsize_not_adaptive():
-    """Checks that a custom stepsize that isn't adaptive is unchanged"""
+    """Checks that a custom stepsize that isn't adaptive is unchanged."""
 
     class CustomStepsize(F.Stepsize):
         def is_adapt(self, t):
@@ -103,7 +99,7 @@ def is_adapt(self, t):
 
 
 def test_growth_factor_setter():
-    """Checks that the growth factor setter works correctly"""
+    """Checks that the growth factor setter works correctly."""
     stepsize = F.Stepsize(1)
 
     # Test that setting a growth factor less than 1 raises a ValueError
@@ -116,7 +112,7 @@ def test_growth_factor_setter():
 
 
 def test_cutback_factor_setter():
-    """Checks that the cutback factor setter works correctly"""
+    """Checks that the cutback factor setter works correctly."""
     stepsize = F.Stepsize(1)
 
     # Test that setting a cutback factor greater than 1 raises a ValueError
@@ -129,10 +125,11 @@ def test_cutback_factor_setter():
 
 
 def test_max_stepsize_setter():
-    """Checks that the maximum stepsize setter works correctly"""
+    """Checks that the maximum stepsize setter works correctly."""
     stepsize = F.Stepsize(initial_value=1)
 
-    # Test that setting a maximum stepsize less than initial stepsize raises a ValueError
+    # Test that setting a maximum stepsize less than initial
+    # stepsize raises a ValueError
     with pytest.raises(
         ValueError, match="maximum stepsize cannot be less than initial stepsize"
     ):
@@ -154,8 +151,7 @@ def test_max_stepsize_setter():
     ],
 )
 def test_next_milestone(milestones, current_time, expected_value):
-    """Checks that the next milestone is
-     identified and set correctly.
+    """Checks that the next milestone is identified and set correctly.
 
     Args:
         milestone (float): next milestone
@@ -170,9 +166,7 @@ def test_next_milestone(milestones, current_time, expected_value):
 
 
 def test_overshoot_milestone():
-    """Test that stepsize is modified
-    when going to overshoot a milestone.
-    """
+    """Test that stepsize is modified when going to overshoot a milestone."""
 
     my_stepsize = F.Stepsize(initial_value=0.1)
     my_stepsize.growth_factor = 1
@@ -193,9 +187,7 @@ def test_overshoot_milestone():
     [(0, 10), (10, 10), (100, 10), (1000, None), (1001, None)],
 )
 def test_get_max_stepsize(t, expected_value):
-    """Tests get_max_stepsize when
-    max_stepsize is a callable.
-    """
+    """Tests get_max_stepsize when max_stepsize is a callable."""
     my_stepsize = F.Stepsize(initial_value=2)
 
     my_stepsize.max_stepsize = lambda t: 10 if t < 1000 else None
@@ -204,9 +196,9 @@ def test_get_max_stepsize(t, expected_value):
 
 
 def test_milestones_without_adaptivity_raises_error():
-    """Checks that setting milestones without adaptivity raises an error"""
+    """Checks that setting milestones without adaptivity raises an error."""
 
     with pytest.raises(
-        ValueError, match="Milestones are only relevant if the stepsize is adaptive."
+        ValueError, match=r"Milestones are only relevant if the stepsize is adaptive."
     ):
         F.Stepsize(initial_value=2, milestones=[1, 2, 3])
diff --git a/test/test_subdomains.py b/test/test_subdomains.py
index 9f87ab145..7b8927662 100644
--- a/test/test_subdomains.py
+++ b/test/test_subdomains.py
@@ -8,7 +8,7 @@
 
 
 def test_different_surface_ids():
-    """Checks that different surface ids are correctly set"""
+    """Checks that different surface ids are correctly set."""
     my_test_model = F.HydrogenTransportProblem()
     my_test_model.species = [F.Species("H")]
     L = 1e-04
@@ -54,7 +54,7 @@ def test_different_volume_ids():
 
 
 def test_non_matching_volume_borders():
-    """Checks that non-matching borders raise an error"""
+    """Checks that non-matching borders raise an error."""
     mesh = F.Mesh1D(vertices=np.linspace(0, 5, 6))
     vol_subdomain_1 = F.VolumeSubdomain1D(id=1, borders=[0, 2], material=None)
     vol_subdomain_2 = F.VolumeSubdomain1D(id=1, borders=[3, 5], material=None)
@@ -65,7 +65,7 @@ def test_non_matching_volume_borders():
 
 
 def test_matching_volume_borders_non_ascending_order():
-    """Checks that subdomain placed in non ascending order still passes"""
+    """Checks that subdomain placed in non ascending order still passes."""
     mesh = F.Mesh1D(vertices=np.linspace(0, 8, 9))
     vol_subdomain_1 = F.VolumeSubdomain1D(id=1, borders=[0, 2], material=None)
     vol_subdomain_2 = F.VolumeSubdomain1D(id=2, borders=[4, 8], material=None)
@@ -76,7 +76,7 @@ def test_matching_volume_borders_non_ascending_order():
 
 
 def test_borders_out_of_domain():
-    """Checks that borders outside of the domain raise an error"""
+    """Checks that borders outside of the domain raise an error."""
     mesh = F.Mesh1D(vertices=np.linspace(0, 2))
     subdomains = [F.VolumeSubdomain1D(id=1, borders=[1, 15], material=None)]
     with pytest.raises(ValueError, match="borders dont match domain borders"):
@@ -84,7 +84,7 @@ def test_borders_out_of_domain():
 
 
 def test_borders_inside_domain():
-    """Checks that borders inside of the domain raise an error"""
+    """Checks that borders inside of the domain raise an error."""
     mesh = F.Mesh1D(vertices=np.linspace(0, 20))
     subdomains = [F.VolumeSubdomain1D(id=1, borders=[1, 6], material=None)]
     with pytest.raises(ValueError, match="borders dont match domain borders"):
@@ -92,7 +92,7 @@ def test_borders_inside_domain():
 
 
 def test_raise_error_with_no_volume_subdomain():
-    """Checks that error is rasied when no volume subdomain is defined"""
+    """Checks that error is rasied when no volume subdomain is defined."""
     mesh = F.Mesh1D(vertices=np.linspace(0, 20))
 
     with pytest.raises(ValueError, match="No volume subdomains defined"):
@@ -101,7 +101,7 @@ def test_raise_error_with_no_volume_subdomain():
 
 @pytest.mark.parametrize("input", [1, 2, 3, 4])
 def test_find_volume_from_int(input):
-    """test that the correct volume is returned when input is an int"""
+    """Test that the correct volume is returned when input is an int."""
 
     vol_1 = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=None)
     vol_2 = F.VolumeSubdomain1D(id=2, borders=[1, 2], material=None)
@@ -114,7 +114,8 @@ def test_find_volume_from_int(input):
 
 
 def test_ValueError_raised_when_id_not_found_in_volumes_subdomains():
-    """test that a ValueError is raised when an id is not found in the list of volume subdomains"""
+    """Test that a ValueError is raised when an id is not found in the list of volume
+    subdomains."""
 
     volumes = [F.VolumeSubdomain1D(id=1, borders=[0, 1], material=None)]
 
@@ -123,9 +124,9 @@ def test_ValueError_raised_when_id_not_found_in_volumes_subdomains():
 
 
 def test_ValueError_rasied_when_volume_ids_are_not_unique():
-    """Checks"""
+    """Checks."""
     my_test_model = F.HydrogenTransportProblem(
-        mesh=F.Mesh1D(np.linspace(0, 2, num=10)), species=[F.Species("H")]
+        mesh=F.Mesh1D(np.linspace(0, 2, num=11)), species=[F.Species("H")]
     )
 
     vol_1 = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=None)
@@ -142,7 +143,8 @@ def test_ValueError_rasied_when_volume_ids_are_not_unique():
 
 
 def test_ValueError_raised_when_id_not_found_in_surface_subdomains():
-    """test that a ValueError is raised when an id is not found in the list of surface subdomains"""
+    """Test that a ValueError is raised when an id is not found in the list of surface
+    subdomains."""
 
     surfaces = [F.SurfaceSubdomain(id=1)]
 
@@ -152,7 +154,7 @@ def test_ValueError_raised_when_id_not_found_in_surface_subdomains():
 
 @pytest.mark.parametrize("input_id, output_index", [(1, 2), (4, 0), (7, 1), (9, 3)])
 def test_find_surface_from_id(input_id, output_index):
-    """test that the correct surface is returned when input is an int"""
+    """Test that the correct surface is returned when input is an int."""
 
     surf_1 = F.SurfaceSubdomain(id=7)
     surf_2 = F.SurfaceSubdomain1D(id=4, x=0)
@@ -165,8 +167,8 @@ def test_find_surface_from_id(input_id, output_index):
 
 
 def test_volume_subdomain_properties():
-    """Tests that the volume subdomain property obtains the correct
-    subdomains from the the model subdomains list"""
+    """Tests that the volume subdomain property obtains the correct subdomains from the
+    the model subdomains list."""
 
     my_model = F.HydrogenTransportProblem()
     my_model.subdomains = [
@@ -209,8 +211,29 @@ def test_name_attribute():
 
 
 def test_name_setter():
-    """Makes sure that a type error is raised when the name type variable is not
-    a string"""
+    """Makes sure that a type error is raised when the name type variable is not a
+    string."""
 
     with pytest.raises(TypeError, match="Name must be a string"):
         F.VolumeSubdomain(id=1, material=None, name=1)
+
+
+def test_all_cells_are_not_tagged():
+    """Checks that an error is raised when not all cells are tagged with a non-zero
+    value.
+
+    This can be caused by a volume subdomain not being defined correctly, or by a mesh
+    that is too coarse to capture the geometry of the volume subdomains.
+    """
+
+    mat = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0)
+
+    vol1 = F.VolumeSubdomain1D(id=1, borders=[0, 0.5], material=mat)
+    vol2 = F.VolumeSubdomain1D(id=2, borders=[0.5, 1], material=mat)
+
+    mesh = F.Mesh1D(np.linspace(0, 1, 10))
+
+    with pytest.raises(
+        AssertionError, match="All cells must be tagged with a non-zero value"
+    ):
+        mesh.define_meshtags(surface_subdomains=[], volume_subdomains=[vol1, vol2])
diff --git a/test/test_surface_quantity.py b/test/test_surface_quantity.py
index 5d4b529b7..c2943bcaf 100644
--- a/test/test_surface_quantity.py
+++ b/test/test_surface_quantity.py
@@ -10,7 +10,7 @@
 
 
 def test_surface_flux_export_compute():
-    """Test that the surface flux export computes the correct value"""
+    """Test that the surface flux export computes the correct value."""
 
     # BUILD
     L = 4.0
@@ -57,7 +57,7 @@ def test_surface_flux_export_compute():
 
 @pytest.mark.parametrize("value", ["my_export.csv", "my_export.txt"])
 def test_title_generation(tmp_path, value):
-    """Test that the title is made to be written to the header in a csv or txt file"""
+    """Test that the title is made to be written to the header in a csv or txt file."""
     my_export = F.SurfaceFlux(
         filename=os.path.join(tmp_path, f"{value}"),
         field=F.Species("TEST"),
@@ -73,7 +73,7 @@ def test_title_generation(tmp_path, value):
 
 
 def test_write_overwrite(tmp_path):
-    """Test that the write method overwrites the file if it already exists"""
+    """Test that the write method overwrites the file if it already exists."""
     filename = os.path.join(tmp_path, "my_export.csv")
     my_export = F.SurfaceFlux(
         filename=filename,
@@ -102,7 +102,7 @@ def test_write_overwrite(tmp_path):
 
 
 def test_filename_setter_raises_TypeError():
-    """Test that a TypeError is raised when the filename is not a string"""
+    """Test that a TypeError is raised when the filename is not a string."""
 
     with pytest.raises(TypeError, match="filename must be of type str"):
         F.SurfaceQuantity(
@@ -113,7 +113,8 @@ def test_filename_setter_raises_TypeError():
 
 
 def test_filename_setter_raises_ValueError(tmp_path):
-    """Test that a ValueError is raised when the filename does not end with .csv or .txt"""
+    """Test that a ValueError is raised when the filename does not end with .csv or
+    .txt."""
 
     with pytest.raises(ValueError):
         F.SurfaceQuantity(
@@ -124,7 +125,7 @@ def test_filename_setter_raises_ValueError(tmp_path):
 
 
 def test_field_setter_raises_TypeError():
-    """Test that a TypeError is raised when the field is not a F.Species"""
+    """Test that a TypeError is raised when the field is not a F.Species."""
 
     with pytest.raises(TypeError):
         F.SurfaceQuantity(
@@ -135,7 +136,7 @@ def test_field_setter_raises_TypeError():
 
 @pytest.mark.parametrize("value", ["my_export.csv", "my_export.txt"])
 def test_writer(tmp_path, value):
-    """Test that the writes values at each timestep to either a csv or txt file"""
+    """Test that the writes values at each timestep to either a csv or txt file."""
     my_export = F.SurfaceFlux(
         filename=os.path.join(tmp_path, f"{value}"),
         field=F.Species("test"),
@@ -153,11 +154,10 @@ def test_writer(tmp_path, value):
 
 
 def test_surface_setter_raises_TypeError():
-    """Test that a TypeError is raised when the surface is not a
-    F.SurfaceSubdomain"""
+    """Test that a TypeError is raised when the surface is not a F.SurfaceSubdomain."""
 
     with pytest.raises(
-        TypeError, match="surface should be an int or F.SurfaceSubdomain"
+        TypeError, match=r"surface should be an int or F.SurfaceSubdomain"
     ):
         F.SurfaceQuantity(
             field=F.Species("H"),
diff --git a/test/test_total_surface.py b/test/test_total_surface.py
index cce43aa7a..7c9bcb89a 100644
--- a/test/test_total_surface.py
+++ b/test/test_total_surface.py
@@ -6,7 +6,7 @@
 
 
 def test_total_surface_compute_1D():
-    """Test that the total surface export computes the correct value"""
+    """Test that the total surface export computes the correct value."""
 
     # BUILD
     L = 6
@@ -16,7 +16,7 @@ def test_total_surface_compute_1D():
     dummy_volume = F.VolumeSubdomain1D(
         id=1, borders=[0, L], material=F.Material(D_0=1, E_D=1, name="dummy")
     )
-    facet_meshtags, temp = my_mesh.define_meshtags(
+    facet_meshtags, _temp = my_mesh.define_meshtags(
         surface_subdomains=[dummy_surface], volume_subdomains=[dummy_volume]
     )
     ds = ufl.Measure("ds", domain=my_mesh.mesh, subdomain_data=facet_meshtags)
diff --git a/test/test_total_volume.py b/test/test_total_volume.py
index 2b2b43e39..226d0c513 100644
--- a/test/test_total_volume.py
+++ b/test/test_total_volume.py
@@ -11,13 +11,13 @@
 
 
 def test_total_volume_export_compute():
-    """Test that the total volume export computes the correct value"""
+    """Test that the total volume export computes the correct value."""
 
     # BUILD
     L = 4.0
     my_mesh = F.Mesh1D(np.linspace(0, L, 10000))
     dummy_volume = F.VolumeSubdomain1D(id=1, borders=[0, L], material=dummy_mat)
-    temp, cell_meshtags = my_mesh.define_meshtags(
+    _temp, cell_meshtags = my_mesh.define_meshtags(
         surface_subdomains=[], volume_subdomains=[dummy_volume]
     )
     dx = ufl.Measure("dx", domain=my_mesh.mesh, subdomain_data=cell_meshtags)
@@ -48,7 +48,7 @@ def test_total_volume_export_compute():
 
 @pytest.mark.parametrize("value", ["my_export.csv", "my_export.txt"])
 def test_title_generation(tmp_path, value):
-    """Test that the title is made to be written to the header in a csv or txt file"""
+    """Test that the title is made to be written to the header in a csv or txt file."""
     my_export = F.TotalVolume(
         filename=os.path.join(tmp_path, f"{value}"),
         field=F.Species("TEST"),
diff --git a/test/test_volume_quantity.py b/test/test_volume_quantity.py
index d93e07c66..7f530ef58 100644
--- a/test/test_volume_quantity.py
+++ b/test/test_volume_quantity.py
@@ -7,7 +7,7 @@
 
 
 def test_filename_setter_raises_TypeError():
-    """Test that a TypeError is raised when the filename is not a string"""
+    """Test that a TypeError is raised when the filename is not a string."""
 
     with pytest.raises(TypeError, match="filename must be of type str"):
         F.VolumeQuantity(
@@ -18,7 +18,8 @@ def test_filename_setter_raises_TypeError():
 
 
 def test_filename_setter_raises_ValueError():
-    """Test that a ValueError is raised when the filename does not end with .csv or .txt"""
+    """Test that a ValueError is raised when the filename does not end with .csv or
+    .txt."""
 
     with pytest.raises(ValueError):
         F.VolumeQuantity(
@@ -29,7 +30,7 @@ def test_filename_setter_raises_ValueError():
 
 
 def test_field_setter_raises_TypeError():
-    """Test that a TypeError is raised when the field is not a F.Species"""
+    """Test that a TypeError is raised when the field is not a F.Species."""
 
     with pytest.raises(TypeError):
         F.VolumeQuantity(
diff --git a/test/test_vtx.py b/test/test_vtx.py
index edbb2e531..bb51bb9d8 100644
--- a/test/test_vtx.py
+++ b/test/test_vtx.py
@@ -11,7 +11,7 @@
 
 
 def test_vtx_export_one_function(tmpdir):
-    """Test can add one function to a vtx export"""
+    """Test can add one function to a vtx export."""
     u = dolfinx.fem.Function(V)
     sp = F.Species("H")
     sp.post_processing_solution = u
@@ -24,7 +24,7 @@ def test_vtx_export_one_function(tmpdir):
 
 
 def test_vtx_export_two_functions(tmpdir):
-    """Test can add two functions to a vtx export"""
+    """Test can add two functions to a vtx export."""
     u = dolfinx.fem.Function(V)
     v = dolfinx.fem.Function(V)
 
@@ -43,8 +43,8 @@ def test_vtx_export_two_functions(tmpdir):
 
 @pytest.mark.skip(reason="Not implemented")
 def test_vtx_export_subdomain():
-    """Test that given multiple subdomains in problem,
-    only correct functions are extracted from species"""
+    """Test that given multiple subdomains in problem, only correct functions are
+    extracted from species."""
     pass
 
 
@@ -55,7 +55,7 @@ def test_vtx_suffix_converter(tmpdir):
 
 
 def test_vtx_DG(tmpdir):
-    """Test VTX export setup for DG formulation"""
+    """Test VTX export setup for DG formulation."""
     my_model = F.HydrogenTransportProblemDiscontinuous()
     my_model.mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0]))
     my_mat = F.Material(D_0=3, E_D=2, K_S_0=1, E_K_S=0, name="mat")
@@ -129,7 +129,7 @@ def test_vtx_integration_with_h_transport_problem(tmpdir, checkpoint):
     ],
 )
 def test_vtx_temperature(T, tmpdir):
-    """Tests that VTX temperature exports work with HydrogenTransportProblem"""
+    """Tests that VTX temperature exports work with HydrogenTransportProblem."""
     my_model = F.HydrogenTransportProblem()
     my_model.mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0]))
     my_mat = F.Material(D_0=1, E_D=0, name="mat")
@@ -154,7 +154,7 @@ def test_vtx_temperature(T, tmpdir):
 
 
 def test_field_attribute_is_always_list():
-    """Test that the field attribute is always a list"""
+    """Test that the field attribute is always a list."""
     my_export = F.VTXSpeciesExport("my_export.bp", field=F.Species("H"))
     assert isinstance(my_export.field, list)
 
@@ -164,18 +164,20 @@ def test_field_attribute_is_always_list():
 
 @pytest.mark.parametrize("field", [["H", 1], 1, [F.Species("H"), 1]])
 def test_field_attribute_raises_error_when_invalid_type(field):
-    """Test that the field attribute raises an error if the type is not festim.Species or list"""
+    """Test that the field attribute raises an error if the type is not festim.Species
+    or list."""
     with pytest.raises(TypeError):
         F.VTXSpeciesExport("my_export.bp", field=field)
 
 
 def test_filename_raises_error_when_wrong_type():
-    """Test that the filename attribute raises an error if the extension is not .bp"""
+    """Test that the filename attribute raises an error if the extension is not .bp."""
     with pytest.raises(TypeError):
         F.VTXSpeciesExport(1, field=[F.Species("H")])
 
 
 def test_filename_temp_raises_error_when_wrong_type():
-    """Test that the filename attribute for VTXTemperature export raises an error if the extension is not .bp"""
+    """Test that the filename attribute for VTXTemperature export raises an error if the
+    extension is not .bp."""
     with pytest.raises(TypeError):
         F.VTXTemperatureExport(1)
diff --git a/test/test_xdmf.py b/test/test_xdmf.py
index fa3387b0f..b990fa8d3 100644
--- a/test/test_xdmf.py
+++ b/test/test_xdmf.py
@@ -11,7 +11,7 @@
 
 
 def test_init():
-    """Tests the initialisation of XDMFExport"""
+    """Tests the initialisation of XDMFExport."""
     species = F.Species("H")
     my_export = F.XDMFExport(filename="my_export.xdmf", field=species)
 
@@ -20,7 +20,7 @@ def test_init():
 
 
 def test_write(tmp_path):
-    """Tests the write method of XDMFExport creates a file"""
+    """Tests the write method of XDMFExport creates a file."""
     species = F.Species("H")
     filename = os.path.join(tmp_path, "test.xdmf")
     my_export = F.XDMFExport(filename=filename, field=species)
@@ -38,7 +38,7 @@ def test_write(tmp_path):
 
 
 def test_integration_with_HTransportProblem(tmp_path):
-    """Tests that XDMFExport can be used in conjunction with HTransportProblem"""
+    """Tests that XDMFExport can be used in conjunction with HTransportProblem."""
     my_model = F.HydrogenTransportProblem()
     my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1))
     my_mat = F.Material(D_0=1.9e-7, E_D=0.2, name="my_mat")
@@ -59,7 +59,7 @@ def test_integration_with_HTransportProblem(tmp_path):
 
 
 def test_field_attribute_is_always_list():
-    """Test that the field attribute is always a list"""
+    """Test that the field attribute is always a list."""
     my_export = F.XDMFExport("my_export.xdmf", field=F.Species("H"))
     assert isinstance(my_export.field, list)
 
@@ -75,15 +75,13 @@ def test_vtx_suffix_converter(tmpdir):
 
 @pytest.mark.parametrize("field", [["H", 2], 1, [F.Species("H"), 1]])
 def test_field_attribute_raises_error_when_invalid_type(field):
-    """
-    Test that the field attribute raises an error if the type
-    is not festim.Species or list.
-    """
+    """Test that the field attribute raises an error if the type is not festim.Species
+    or list."""
     with pytest.raises(TypeError):
         F.XDMFExport("my_export.xdmf", field=field)
 
 
 def test_filename_raises_error_when_wrong_type():
-    """Test that the filename attribute raises an error if the file is not str"""
+    """Test that the filename attribute raises an error if the file is not str."""
     with pytest.raises(TypeError):
         F.XDMFExport(1, field=[F.Species("H")])