diff --git a/test/h_transport_problem/test_initialisation.py b/test/h_transport_problem/test_initialisation.py
index e0cfa77c0..347f57a00 100644
--- a/test/h_transport_problem/test_initialisation.py
+++ b/test/h_transport_problem/test_initialisation.py
@@ -2,9 +2,18 @@
 import pytest
 import fenics
 from pathlib import Path
+import os
 
 
 def test_initialisation_from_xdmf(tmpdir):
+    """
+    Test that initialise_solutions interpolates correctly
+    from an XDMF-file
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
+    """
+
     mesh = fenics.UnitSquareMesh(5, 5)
     V = fenics.VectorFunctionSpace(mesh, "P", 1, 2)
     u = fenics.Function(V)
@@ -19,7 +28,6 @@ def test_initialisation_from_xdmf(tmpdir):
     d = tmpdir.mkdir("Initial solutions")
     file1 = d.join("u_1out.xdmf")
     file2 = d.join("u_2out.xdmf")
-    print(Path(file1))
     with fenics.XDMFFile(str(Path(file1))) as f:
         f.write_checkpoint(
             u.sub(0), "1", 2, fenics.XDMFFile.Encoding.HDF5, append=False
diff --git a/test/h_transport_problem/test_solving.py b/test/h_transport_problem/test_solving.py
index 8aa7d5dae..a292b7374 100644
--- a/test/h_transport_problem/test_solving.py
+++ b/test/h_transport_problem/test_solving.py
@@ -14,7 +14,9 @@ def test_default_dt_min_value():
 
     t = 0
     dt = festim.Stepsize(
-        initial_value=0.5, stepsize_change_ratio=1.1, t_stop=430, stepsize_stop_max=0.5
+        initial_value=0.5,
+        stepsize_change_ratio=1.1,
+        max_stepsize=None if t < 430 else 0.5,
     )
 
     my_settings = festim.Settings(
@@ -32,7 +34,8 @@ def test_default_dt_min_value():
     my_problem.F = f.dot(f.grad(my_problem.u), f.grad(my_problem.v)) * f.dx
 
     du = f.TrialFunction(my_problem.u.function_space())
-    my_problem.J = f.derivative(my_problem.F, my_problem.u, du)  # Define the Jacobian
+    # Define the Jacobian
+    my_problem.J = f.derivative(my_problem.F, my_problem.u, du)
     # run & test
     my_problem.update(t, dt)
 
diff --git a/test/simulation/test_construction.py b/test/simulation/test_construction.py
index 767652a66..4020d63ed 100644
--- a/test/simulation/test_construction.py
+++ b/test/simulation/test_construction.py
@@ -2,85 +2,66 @@
 import pytest
 
 
-def test_setting_traps():
-    """Checks traps can be set with the expected types (F.Trap, list, or F.Traps)"""
+class TestSettingTrapsMaterialsExprots:
     my_sim = F.Simulation()
-    my_mat = F.Material(1, 1, 0)
-    trap1 = F.Trap(1, 1, 1, 1, [my_mat], density=1)
-    trap2 = F.Trap(2, 2, 2, 2, [my_mat], density=1)
-
-    combinations = [trap1, [trap1], [trap1, trap2], F.Traps([trap1, trap2])]
-
-    for combination in combinations:
-        my_sim.traps = combination
+    mat1 = F.Material(1, 1, 1)
+    mat2 = F.Material(2, 2, 2)
+    trap1 = F.Trap(1, 1, 1, 1, [mat1], density=1)
+    trap2 = F.Trap(2, 2, 2, 2, [mat1], density=1)
+    export1 = F.XDMFExport("solute")
+    export2 = F.XDMFExport("solute")
+    wrong_types = ["coucou", True]
 
+    @pytest.mark.parametrize(
+        "trap", [trap1, [trap1], [trap1, trap2], F.Traps([trap1, trap2])]
+    )
+    def test_setting_traps(self, trap):
+        """Checks traps can be set with the expected types (F.Trap, list, or F.Traps)"""
 
-def test_setting_traps_wrong_type():
-    """Checks an error is raised when traps is set with the wrong type"""
-    my_sim = F.Simulation()
+        self.my_sim.traps = trap
 
-    combinations = ["coucou", True]
+    @pytest.mark.parametrize("trap", wrong_types)
+    def test_setting_traps_wrong_type(self, trap):
+        """Checks an error is raised when traps is set with the wrong type"""
 
-    for combination in combinations:
         with pytest.raises(
             TypeError,
             match="Accepted types for traps are list, festim.Traps or festim.Trap",
         ):
-            my_sim.traps = combination
-
-
-def test_setting_materials():
-    """Checks materials can be set with the expected types (F.Material, list, or F.Materials)"""
-    my_sim = F.Simulation()
-    mat1 = F.Material(1, 1, 1)
-    mat2 = F.Material(2, 2, 2)
+            self.my_sim.traps = trap
 
-    combinations = [mat1, [mat1], [mat1, mat2], F.Materials([mat1, mat2])]
+    @pytest.mark.parametrize(
+        "mat", [mat1, [mat1], [mat1, mat2], F.Materials([mat1, mat2])]
+    )
+    def test_setting_materials(self, mat):
+        """Checks materials can be set with the expected types (F.Material, list, or F.Materials)"""
 
-    for combination in combinations:
-        my_sim.materials = combination
+        self.my_sim.materials = mat
 
+    @pytest.mark.parametrize("mat", wrong_types)
+    def test_setting_materials_wrong_type(self, mat):
+        """Checks an error is raised when materials is set with the wrong type"""
 
-def test_setting_materials_wrong_type():
-    """Checks an error is raised when materials is set with the wrong type"""
-    my_sim = F.Simulation()
-
-    combinations = ["coucou", True]
-
-    for combination in combinations:
         with pytest.raises(
             TypeError,
             match="accepted types for materials are list, festim.Material or festim.Materials",
         ):
-            my_sim.materials = combination
+            self.my_sim.materials = mat
 
+    @pytest.mark.parametrize(
+        "exp", [export1, [export1], [export1, export2], F.Exports([export1, export2])]
+    )
+    def test_setting_exports(self, exp):
+        """Checks exports can be set with the expected types (F.Material, list, or F.Exports)"""
 
-def test_setting_exports():
-    """Checks exports can be set with the expected types (F.Material, list, or F.Exports)"""
-    my_sim = F.Simulation()
-    export1 = F.XDMFExport("solute")
-    export2 = F.XDMFExport("solute")
-
-    combinations = [
-        export1,
-        [export1],
-        [export1, export2],
-        F.Exports([export1, export2]),
-    ]
-
-    for trap_combination in combinations:
-        my_sim.exports = trap_combination
-
-
-def test_setting_exports_wrong_type():
-    """Checks an error is raised when exports is set with the wrong type"""
-    my_sim = F.Simulation()
+        self.my_sim.exports = exp
 
-    combinations = ["coucou", True]
+    @pytest.mark.parametrize("exp", wrong_types)
+    def test_setting_exports_wrong_type(self, exp):
+        """Checks an error is raised when exports is set with the wrong type"""
 
-    for trap_combination in combinations:
         with pytest.raises(
             TypeError,
             match="accepted types for exports are list, festim.Export or festim.Exports",
         ):
-            my_sim.exports = trap_combination
+            self.my_sim.exports = exp
diff --git a/test/simulation/test_initialise.py b/test/simulation/test_initialise.py
index da78696f0..83653d93d 100644
--- a/test/simulation/test_initialise.py
+++ b/test/simulation/test_initialise.py
@@ -79,6 +79,9 @@ def test_TXTExport_times_added_to_milestones(tmpdir):
     """Creates a Simulation object and checks that, if no dt.milestones
      are given and TXTExport.times are given, TXTExport.times are
     are added to dt.milestones by .initialise()
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     # tmpdir
     d = tmpdir.mkdir("test_folder")
@@ -120,8 +123,13 @@ def test_TXTExport_times_added_to_milestones(tmpdir):
 )
 @pytest.mark.parametrize("sys", ["cylindrical", "spherical"])
 def test_cartesian_and_surface_flux_warning(quantity, sys):
-    """Creates a Simulation object and checks that, if either a cylindrical
+    """
+    Creates a Simulation object and checks that, if either a cylindrical
     or spherical meshes are given with a SurfaceFlux, a warning is raised.
+
+    Args:
+        quantity (festim.DerivedQuantity): a festim.DerivedQuantity object
+        sys (str): type of the coordinate system
     """
     # build
     my_model = F.Simulation()
diff --git a/test/simulation/test_meshing.py b/test/simulation/test_meshing.py
index d5b7081df..b8a717393 100644
--- a/test/simulation/test_meshing.py
+++ b/test/simulation/test_meshing.py
@@ -6,6 +6,9 @@
 def test_define_markers(tmpdir):
     """Checks that markers can be defined from XDMF files and that the mesh
     functions have the correct values.
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     # build
     mesh = fenics.UnitSquareMesh(16, 16)
diff --git a/test/simulation/test_postprocessing_integration.py b/test/simulation/test_postprocessing_integration.py
index b5f1deb27..6681dc398 100644
--- a/test/simulation/test_postprocessing_integration.py
+++ b/test/simulation/test_postprocessing_integration.py
@@ -9,6 +9,7 @@
 class TestPostProcessing:
     @pytest.fixture
     def my_sim(self):
+        """A pytest fixture defining the festim.Simulation object"""
         my_sim = festim.Simulation({})
         my_sim.t = 0
         my_sim.mesh = festim.MeshFromRefinements(10, 1)
@@ -42,6 +43,13 @@ def my_sim(self):
         return my_sim
 
     def test_derived_quantities_size(self, my_sim):
+        """
+        Checks that the length of data produced by F.DerivedQuantities
+        equals to the number of post-processing calls
+
+        Args:
+            my_sim (festim.Simulation): the simulation object
+        """
         derived_quantities = festim.DerivedQuantities(
             [
                 festim.SurfaceFlux("solute", 1),
@@ -64,6 +72,12 @@ def test_derived_quantities_size(self, my_sim):
         assert my_sim.exports[0].data[i][0] == t
 
     def test_pure_diffusion(self, my_sim):
+        """
+        Checks that run_post_processing() assigns data correctly
+
+        Args:
+            my_sim (festim.Simulation): the simulation object
+        """
         my_sim.materials = festim.Materials(
             [
                 festim.Material(1, D_0=1, E_D=1, borders=[0, 0.5]),
@@ -164,6 +178,10 @@ def test_performance_xdmf_export_every_N_iterations(self, my_sim, tmpdir):
         """Runs run_post_processing several times with different export.mode
         values and checks that the xdmf
         files have the correct timesteps
+
+        Args:
+            my_sim (festim.Simulation): a festim.Simulation object
+            tmpdir (os.PathLike): path to the pytest temporary folder
         """
         # build
         d = tmpdir.mkdir("test_folder")
@@ -205,8 +223,8 @@ def test_xdmf_export_only_last_timestep(self, my_sim, tmpdir):
         correct timesteps
 
         Args:
-            my_sim (_type_): _description_
-            tmpdir (_type_): _description_
+            my_sim (festim.Simulation): the simulation object
+            tmpdir (os.PathLike): path to the pytest temporary folder
         """
         d = tmpdir.mkdir("test_folder")
 
diff --git a/test/system/test_chemical_potential.py b/test/system/test_chemical_potential.py
index 60b090e70..7649095e4 100644
--- a/test/system/test_chemical_potential.py
+++ b/test/system/test_chemical_potential.py
@@ -7,8 +7,16 @@
 
 
 def compute_error(exact, computed, t, norm):
-    exact_sol = fenics.Expression(sp.printing.ccode(exact), degree=4, t=t)
+    """
+    An auxiliary method to compute the error
 
+    Args:
+        exact (sympy.Expr): exact solution
+        computed (fenics.Function): computed solution
+        t (float): simulation time
+        norm (str): type of norm (maximum absolute error or L2 norm)
+    """
+    exact_sol = fenics.Expression(sp.printing.ccode(exact), degree=4, t=t)
     if norm == "error_max":
         mesh = computed.function_space().mesh()
         vertex_values_u = computed.compute_vertex_values(mesh)
@@ -24,6 +32,9 @@ def test_run_MMS_chemical_pot(tmpdir):
     """
     Test function run() with conservation of chemical potential henry
     (1 material)
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     d = tmpdir.mkdir("Solution_Test")
     u = 1 + sp.sin(2 * fenics.pi * festim.x) * festim.t + festim.t
diff --git a/test/system/test_misc.py b/test/system/test_misc.py
index 1dbc3fdbb..ec3e3804e 100644
--- a/test/system/test_misc.py
+++ b/test/system/test_misc.py
@@ -4,7 +4,13 @@
 import os
 
 
-def test_convective_flux():
+def test_convective_flux(tmpdir):
+    """
+    Tests that convective boundary condition works correctly
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
+    """
     sim = F.Simulation()
 
     sim.mesh = F.MeshFromVertices(np.linspace(0, 1, num=50))
@@ -19,7 +25,13 @@ def test_convective_flux():
 
     sim.materials = F.Materials([F.Material(1, D_0=1, E_D=0, thermal_cond=2)])
 
-    sim.exports = F.Exports([F.XDMFExport("T", checkpoint=False)])
+    sim.exports = F.Exports(
+        [
+            F.XDMFExport(
+                "T", checkpoint=False, filename="{}/temperature.xdmf".format(tmpdir)
+            )
+        ]
+    )
 
     sim.settings = F.Settings(1e-10, 1e-10, transient=False)
 
@@ -186,6 +198,9 @@ def test_txt_export_desired_times(tmp_path, final_time, stepsize, export_times):
 def test_txt_export_all_times(tmp_path):
     """
     Tests that TXTExport can be exported at all timesteps
+
+    Args:
+        tmp_path (os.PathLike): path to a temporary folder
     """
     my_model = F.Simulation()
 
@@ -214,6 +229,9 @@ def test_txt_export_all_times(tmp_path):
 def test_txt_export_steady_state(tmp_path):
     """
     Tests that TXTExport can be exported in steady state
+
+    Args:
+        tmp_path (os.PathLike): path to a temporary folder
     """
     my_model = F.Simulation()
 
diff --git a/test/system/test_system.py b/test/system/test_system.py
index 41576d03a..671068b94 100644
--- a/test/system/test_system.py
+++ b/test/system/test_system.py
@@ -7,6 +7,15 @@
 
 
 def compute_error(exact, computed, t, norm):
+    """
+    An auxiliary method to compute the error
+
+    Args:
+        exact (sympy.Expr): exact solution
+        computed (fenics.Function): computed solution
+        t (float): simulation time
+        norm (str): type of norm (maximum absolute error or L2 norm)
+    """
     exact_sol = fenics.Expression(sp.printing.ccode(exact), degree=4, t=t)
 
     if norm == "error_max":
@@ -26,6 +35,9 @@ def compute_error(exact, computed, t, norm):
 def test_run_temperature_stationary(tmpdir):
     """
     Check that the temperature module works well in 1D stationary
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     d = tmpdir.mkdir("Solution_Test")
     u = 1 + 2 * festim.x**2
@@ -82,6 +94,9 @@ def test_run_temperature_stationary(tmpdir):
 def test_run_temperature_transient(tmpdir):
     """
     Check that the temperature module works well in 1D transient
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     d = tmpdir.mkdir("Solution_Test")
     u = 1 + 2 * festim.x**2 + festim.t
@@ -128,8 +143,7 @@ def test_run_temperature_transient(tmpdir):
     my_dt = festim.Stepsize(
         initial_value=0.5,
         stepsize_change_ratio=1,
-        t_stop=40,
-        stepsize_stop_max=0.5,
+        max_stepsize=lambda t: None if t < 40 else 0.5,
         dt_min=1e-5,
     )
 
@@ -161,6 +175,9 @@ def test_run_temperature_transient(tmpdir):
 def test_run_MMS(tmpdir):
     """
     Test function run() for several refinements
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     d = tmpdir.mkdir("Solution_Test")
     u = 1 + sp.sin(2 * fenics.pi * festim.x) * festim.t
@@ -283,6 +300,9 @@ def run(h):
 def test_run_MMS_chemical_pot(tmpdir):
     """
     Test function run() with conservation of chemical potential (1 material)
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     d = tmpdir.mkdir("Solution_Test")
     u = 1 + sp.sin(2 * fenics.pi * festim.x) * festim.t + festim.t
@@ -409,6 +429,9 @@ def test_run_chemical_pot_mass_balance(tmpdir):
     increases.
     Creates a model with a constant concentration of mobile (c_m(t=0)=1,
     non-flux conditions at surfaces) with a varying temperature
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     d = tmpdir.mkdir("Solution_Test")
     my_materials = festim.Materials(
@@ -465,6 +488,9 @@ def test_run_chemical_pot_mass_balance(tmpdir):
 def test_run_MMS_soret(tmpdir):
     """
     MMS test with soret effect
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     d = tmpdir.mkdir("Solution_Test")
     u = 1 + festim.x**2 + festim.t
@@ -557,6 +583,9 @@ def run(h):
 def test_run_MMS_steady_state(tmpdir):
     """
     MMS test with one trap at steady state
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     d = tmpdir.mkdir("Solution_Test")
     u = 1 + festim.x
@@ -681,6 +710,9 @@ def test_chemical_pot_T_solve_stationary(tmpdir):
     type solve_stationary for temperature
 
     adapted to catch bug described in issue #310
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     d = tmpdir.mkdir("Solution_Test")
     my_materials = festim.Materials(
@@ -732,6 +764,9 @@ def test_chemical_pot_T_solve_stationary(tmpdir):
 def test_export_particle_flux_with_chemical_pot(tmpdir):
     """Checks that surface particle fluxes can be computed with conservation
     of chemical potential
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     d = tmpdir.mkdir("Solution_Test")
     my_materials = festim.Materials(
diff --git a/test/unit/test_concentration.py b/test/unit/test_concentration.py
index 215629040..f7cdeb39f 100644
--- a/test/unit/test_concentration.py
+++ b/test/unit/test_concentration.py
@@ -5,6 +5,10 @@
 
 
 class TestGetComp:
+    """
+    General test for the get_comp method of the festim.Concentration class
+    """
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
 
@@ -44,6 +48,10 @@ def test_get_comp_from_xdmf(self, tmpdir):
 
 
 class TestInitialise:
+    """
+    General test for the initialise method of the festim.Concentration class
+    """
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
     u = f.Function(V)
diff --git a/test/unit/test_exports/test_derived_quantities/test_average_surface.py b/test/unit/test_exports/test_derived_quantities/test_average_surface.py
index 864535a2f..556dd7865 100644
--- a/test/unit/test_exports/test_derived_quantities/test_average_surface.py
+++ b/test/unit/test_exports/test_derived_quantities/test_average_surface.py
@@ -1,22 +1,26 @@
 from festim import AverageSurface
 import fenics as f
+import pytest
 
 
-def test_title_H():
-    surface = 1
-    field = "solute"
-    my_average = AverageSurface(field, surface)
-    assert my_average.title == "Average {} surface {}".format(field, surface)
+@pytest.mark.parametrize("field, surface", [("solute", 1), ("T", 2)])
+def test_title(field, surface):
+    """
+    A simple test to check that the title is set
+    correctly in festim.AverageSurface
 
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+    """
 
-def test_title_T():
-    surface = 2
-    field = "T"
     my_average = AverageSurface(field, surface)
     assert my_average.title == "Average {} surface {}".format(field, surface)
 
 
 class TestCompute:
+    """Test that the average surface export computes the correct value"""
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
 
diff --git a/test/unit/test_exports/test_derived_quantities/test_average_volume.py b/test/unit/test_exports/test_derived_quantities/test_average_volume.py
index 1b0a5d901..172f9b5e6 100644
--- a/test/unit/test_exports/test_derived_quantities/test_average_volume.py
+++ b/test/unit/test_exports/test_derived_quantities/test_average_volume.py
@@ -1,22 +1,26 @@
 from festim import AverageVolume
 import fenics as f
+import pytest
 
 
-def test_title_H():
-    volume = 1
-    field = "solute"
-    my_average = AverageVolume(field, volume)
-    assert my_average.title == "Average {} volume {}".format(field, volume)
+@pytest.mark.parametrize("field,volume", [("solute", 1), ("T", 2)])
+def test_title(field, volume):
+    """
+    A simple test to check that the title is set
+    correctly in festim.AverageVolume
 
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        volume (int): the volume id
+    """
 
-def test_title_T():
-    volume = 2
-    field = "T"
     my_average = AverageVolume(field, volume)
     assert my_average.title == "Average {} volume {}".format(field, volume)
 
 
 class TestCompute:
+    """Test that the average volume export computes the correct value"""
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
 
diff --git a/test/unit/test_exports/test_derived_quantities/test_derived_quantities.py b/test/unit/test_exports/test_derived_quantities/test_derived_quantities.py
index f61cbb585..d7eb90124 100644
--- a/test/unit/test_exports/test_derived_quantities/test_derived_quantities.py
+++ b/test/unit/test_exports/test_derived_quantities/test_derived_quantities.py
@@ -26,12 +26,20 @@ class TestMakeHeader:
     max_vol_1 = MaximumVolume("T", 2)
 
     def test_simple(self):
+        """
+        Tests that a proper header is made for the output .csv file
+        when there is only one festim.DerivedQuantity object
+        """
         my_derv_quant = DerivedQuantities([self.surface_flux_1])
         header = my_derv_quant.make_header()
         expected_header = ["t(s)", self.surface_flux_1.title]
         assert header == expected_header
 
     def test_two_quantities(self):
+        """
+        Tests that a proper header is made for the output .csv file
+        when there are two festim.DerivedQuantity objects
+        """
         my_derv_quant = DerivedQuantities(
             [
                 self.surface_flux_1,
@@ -43,6 +51,10 @@ def test_two_quantities(self):
         assert header == expected_header
 
     def test_all_quantities(self):
+        """
+        Tests that a proper header is made for the output .csv file
+        when there are many festim.DerivedQuantity objects
+        """
         my_derv_quant = DerivedQuantities(
             [
                 self.surface_flux_1,
@@ -66,6 +78,11 @@ def test_all_quantities(self):
 
 
 class TestAssignMeasuresToQuantities:
+    """
+    Tests that measure attributes are properly assigned to all
+    festim.DerivedQuantity objects
+    """
+
     my_quantities = DerivedQuantities(
         [
             SurfaceFlux("solute", 2),
@@ -81,19 +98,27 @@ class TestAssignMeasuresToQuantities:
     my_quantities.assign_measures_to_quantities(dx, ds)
 
     def test_quantities_have_dx(self):
+        """Check for the volume measure"""
         for quantity in self.my_quantities:
             assert quantity.dx == self.dx
 
     def test_quantities_have_ds(self):
+        """Check for the surface measure"""
         for quantity in self.my_quantities:
             assert quantity.ds == self.ds
 
     def test_quantities_have_n(self):
+        """Check for the normal vector of the surface"""
         for quantity in self.my_quantities:
             assert quantity.n == self.n
 
 
 class TestAssignPropertiesToQuantities:
+    """
+    Tests that property attributes are properly assigned to all
+    festim.DerivedQuantity objects
+    """
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
     my_quantities = DerivedQuantities(
@@ -113,23 +138,29 @@ class TestAssignPropertiesToQuantities:
     my_quantities.assign_properties_to_quantities(my_mats)
 
     def test_quantities_have_D(self):
+        """Check for diffusivity"""
         for quantity in self.my_quantities:
             assert quantity.D == self.my_mats.D
 
     def test_quantities_have_S(self):
+        """Check for solubility"""
         for quantity in self.my_quantities:
             assert quantity.S == self.my_mats.S
 
     def test_quantities_have_Q(self):
+        """Check for Soret"""
         for quantity in self.my_quantities:
             assert quantity.Q == self.my_mats.Q
 
     def test_quantities_have_thermal_cond(self):
+        """Check for thermal conductivity"""
         for quantity in self.my_quantities:
             assert quantity.thermal_cond == self.my_mats.thermal_cond
 
 
 class TestCompute:
+    """Test that the derived qunatities compute the correct value"""
+
     surface_flux_1 = SurfaceFlux("solute", 2)
     surface_flux_2 = SurfaceFlux("T", 3)
     average_vol_1 = AverageVolume("solute", 1)
@@ -168,6 +199,7 @@ class TestCompute:
     my_mats.thermal_cond = f.interpolate(f.Constant(2), V)
 
     def test_simple(self):
+        """Check for the case of one festim.DerivedQuantity object"""
         my_derv_quant = DerivedQuantities([self.surface_flux_1])
         for quantity in my_derv_quant:
             quantity.function = self.label_to_function[quantity.field]
@@ -182,6 +214,7 @@ def test_simple(self):
         assert my_derv_quant.data[0] == expected_data
 
     def test_two_quantities(self):
+        """Check for the case of two festim.DerivedQuantity objects"""
         my_derv_quant = DerivedQuantities(
             [
                 self.surface_flux_1,
@@ -202,6 +235,7 @@ def test_two_quantities(self):
         assert my_derv_quant.data[0] == expected_data
 
     def test_all_quantities(self):
+        """Check for the case of many festim.DerivedQuantity objects"""
         my_derv_quant = DerivedQuantities(
             [
                 self.surface_flux_1,
diff --git a/test/unit/test_exports/test_derived_quantities/test_hydrogen_flux.py b/test/unit/test_exports/test_derived_quantities/test_hydrogen_flux.py
index 384fc77b6..1b1b71181 100644
--- a/test/unit/test_exports/test_derived_quantities/test_hydrogen_flux.py
+++ b/test/unit/test_exports/test_derived_quantities/test_hydrogen_flux.py
@@ -2,5 +2,10 @@
 
 
 def test_field_is_solute():
+    """
+    Tests that the festim.SurfaceQuantity field is set to solute
+    when festim.HydrogenFlux is used
+    """
+
     my_flux = HydrogenFlux(2)
     assert my_flux.field == "solute"
diff --git a/test/unit/test_exports/test_derived_quantities/test_maximum_surface.py b/test/unit/test_exports/test_derived_quantities/test_maximum_surface.py
index e65dc175f..b97a707d6 100644
--- a/test/unit/test_exports/test_derived_quantities/test_maximum_surface.py
+++ b/test/unit/test_exports/test_derived_quantities/test_maximum_surface.py
@@ -1,23 +1,27 @@
 from festim import MaximumSurface
 import fenics as f
 import numpy as np
+import pytest
 
 
-def test_title_H():
-    surface = 1
-    field = "solute"
-    my_max = MaximumSurface(field, surface)
-    assert my_max.title == "Maximum {} surface {}".format(field, surface)
+@pytest.mark.parametrize("field,surface", [("solute", 1), ("T", 2)])
+def test_title(field, surface):
+    """
+    A simple test to check that the title is set
+    correctly in festim.MaximumSurface
 
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+    """
 
-def test_title_T():
-    surface = 2
-    field = "T"
     my_max = MaximumSurface(field, surface)
     assert my_max.title == "Maximum {} surface {}".format(field, surface)
 
 
 class TestCompute:
+    """Test that the maximum surface export computes the correct value"""
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
 
diff --git a/test/unit/test_exports/test_derived_quantities/test_maximum_volume.py b/test/unit/test_exports/test_derived_quantities/test_maximum_volume.py
index 4f3ff73d8..2dffe5d97 100644
--- a/test/unit/test_exports/test_derived_quantities/test_maximum_volume.py
+++ b/test/unit/test_exports/test_derived_quantities/test_maximum_volume.py
@@ -1,23 +1,27 @@
 from festim import MaximumVolume
 import fenics as f
 import numpy as np
+import pytest
 
 
-def test_title_H():
-    volume = 1
-    field = "solute"
-    my_max = MaximumVolume(field, volume)
-    assert my_max.title == "Maximum {} volume {}".format(field, volume)
+@pytest.mark.parametrize("field,volume", [("solute", 1), ("T", 2)])
+def test_title(field, volume):
+    """
+    A simple test to check that the title is set
+    correctly in festim.MaximumVolume
 
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        volume (int): the volume id
+    """
 
-def test_title_T():
-    volume = 2
-    field = "T"
     my_max = MaximumVolume(field, volume)
     assert my_max.title == "Maximum {} volume {}".format(field, volume)
 
 
 class TestCompute:
+    """Test that the maximum volume export computes the correct value"""
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
 
diff --git a/test/unit/test_exports/test_derived_quantities/test_minimum_surface.py b/test/unit/test_exports/test_derived_quantities/test_minimum_surface.py
index bc816619f..c2b7801a2 100644
--- a/test/unit/test_exports/test_derived_quantities/test_minimum_surface.py
+++ b/test/unit/test_exports/test_derived_quantities/test_minimum_surface.py
@@ -1,23 +1,27 @@
 from festim import MinimumSurface
 import fenics as f
 import numpy as np
+import pytest
 
 
-def test_title_H():
-    surface = 1
-    field = "solute"
-    my_min = MinimumSurface(field, surface)
-    assert my_min.title == "Minimum {} surface {}".format(field, surface)
+@pytest.mark.parametrize("field,surface", [("solute", 1), ("T", 2)])
+def test_title(field, surface):
+    """
+    A simple test to check that the title is set
+    correctly in festim.MinimumSurface
 
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+    """
 
-def test_title_T():
-    surface = 2
-    field = "T"
     my_min = MinimumSurface(field, surface)
     assert my_min.title == "Minimum {} surface {}".format(field, surface)
 
 
 class TestCompute:
+    """Test that the minimum surface export computes the correct value"""
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
 
diff --git a/test/unit/test_exports/test_derived_quantities/test_minimum_volume.py b/test/unit/test_exports/test_derived_quantities/test_minimum_volume.py
index e1285c13c..61fd983bb 100644
--- a/test/unit/test_exports/test_derived_quantities/test_minimum_volume.py
+++ b/test/unit/test_exports/test_derived_quantities/test_minimum_volume.py
@@ -1,23 +1,27 @@
 from festim import MinimumVolume
 import fenics as f
 import numpy as np
+import pytest
 
 
-def test_title_H():
-    volume = 1
-    field = "solute"
-    my_min = MinimumVolume(field, volume)
-    assert my_min.title == "Minimum {} volume {}".format(field, volume)
+@pytest.mark.parametrize("field,volume", [("solute", 1), ("T", 2)])
+def test_title(field, volume):
+    """
+    A simple test to check that the title is set
+    correctly in festim.MinimumVolume
 
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+    """
 
-def test_title_T():
-    volume = 2
-    field = "T"
     my_min = MinimumVolume(field, volume)
     assert my_min.title == "Minimum {} volume {}".format(field, volume)
 
 
 class TestCompute:
+    """Test that the minimum volume export computes the correct value"""
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
 
diff --git a/test/unit/test_exports/test_derived_quantities/test_point_value.py b/test/unit/test_exports/test_derived_quantities/test_point_value.py
index 7c0fa4f31..fc49133f4 100644
--- a/test/unit/test_exports/test_derived_quantities/test_point_value.py
+++ b/test/unit/test_exports/test_derived_quantities/test_point_value.py
@@ -5,6 +5,13 @@
 
 @pytest.mark.parametrize("field", ["solute", "T"])
 def test_title(field):
+    """
+    A simple test to check that the title is set
+    correctly in festim.PointValue
+
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+    """
     x = 1
     my_value = PointValue(field, x)
     assert my_value.title == "{} value at [{}]".format(field, x)
@@ -14,6 +21,7 @@ def test_title(field):
     "mesh,x", [(f.UnitIntervalMesh(10), 1), (f.UnitCubeMesh(10, 10, 10), (1, 0, 1))]
 )
 def test_point_compute(mesh, x):
+    """Test that the point value export computes the correct value"""
     V = f.FunctionSpace(mesh, "P", 1)
     c = f.interpolate(f.Expression("x[0]", degree=1), V)
 
diff --git a/test/unit/test_exports/test_derived_quantities/test_surface_flux.py b/test/unit/test_exports/test_derived_quantities/test_surface_flux.py
index 668d805e6..55446254e 100644
--- a/test/unit/test_exports/test_derived_quantities/test_surface_flux.py
+++ b/test/unit/test_exports/test_derived_quantities/test_surface_flux.py
@@ -5,21 +5,24 @@
 import pytest
 
 
-def test_title_H():
-    surface = 1
-    field = "solute"
-    my_h_flux = SurfaceFlux(field, surface)
-    assert my_h_flux.title == "Flux surface {}: {}".format(surface, field)
+@pytest.mark.parametrize("field,surface", [("solute", 1), ("T", 2)])
+def test_title(field, surface):
+    """
+    A simple test to check that the title is set
+    correctly in festim.SurfaceFlux
 
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+    """
 
-def test_title_heat():
-    surface = 2
-    field = "T"
     my_h_flux = SurfaceFlux(field, surface)
     assert my_h_flux.title == "Flux surface {}: {}".format(surface, field)
 
 
 class TestCompute:
+    """Test that the surface flux export computes the correct value"""
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
 
diff --git a/test/unit/test_exports/test_derived_quantities/test_surface_quantity.py b/test/unit/test_exports/test_derived_quantities/test_surface_quantity.py
index 45e70aa68..e790ef9d2 100644
--- a/test/unit/test_exports/test_derived_quantities/test_surface_quantity.py
+++ b/test/unit/test_exports/test_derived_quantities/test_surface_quantity.py
@@ -3,8 +3,13 @@
 import pytest
 
 
-def test_wrong_type_for_surface():
-    for surface in [2.1, [0], [1, 1], "coucou", ["coucou"], True]:
-        with pytest.raises(TypeError, match="surface should be an int"):
-            print(surface)
-            F.SurfaceQuantity(field=0, surface=surface)
+@pytest.mark.parametrize("surface", [2.1, [0], [1, 1], "coucou", ["coucou"], True])
+def test_wrong_type_for_surface(surface):
+    """
+    Tests that error is raised when the surface attribute is set with wrong type
+
+    Args:
+        surface (): wrong type for the surface attribute
+    """
+    with pytest.raises(TypeError, match="surface should be an int"):
+        F.SurfaceQuantity(field=0, surface=surface)
diff --git a/test/unit/test_exports/test_derived_quantities/test_thermal_flux.py b/test/unit/test_exports/test_derived_quantities/test_thermal_flux.py
index 935d58484..fac0e5f46 100644
--- a/test/unit/test_exports/test_derived_quantities/test_thermal_flux.py
+++ b/test/unit/test_exports/test_derived_quantities/test_thermal_flux.py
@@ -2,5 +2,9 @@
 
 
 def test_field_is_T():
+    """
+    Tests that the festim.SurfaceQuantity field is set to T
+    when festim.ThermalFlux is used
+    """
     my_flux = ThermalFlux(2)
     assert my_flux.field == "T"
diff --git a/test/unit/test_exports/test_derived_quantities/test_total_surface.py b/test/unit/test_exports/test_derived_quantities/test_total_surface.py
index 647da8166..f5e152595 100644
--- a/test/unit/test_exports/test_derived_quantities/test_total_surface.py
+++ b/test/unit/test_exports/test_derived_quantities/test_total_surface.py
@@ -1,23 +1,27 @@
 from this import s
 from festim import TotalSurface
 import fenics as f
+import pytest
 
 
-def test_title_H():
-    surface = 1
-    field = "solute"
-    my_total = TotalSurface(field, surface)
-    assert my_total.title == "Total {} surface {}".format(field, surface)
+@pytest.mark.parametrize("field,surface", [("solute", 1), ("T", 2)])
+def test_title(field, surface):
+    """
+    A simple test to check that the title is set
+    correctly in festim.TotalSurface
 
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+    """
 
-def test_title_T():
-    surface = 2
-    field = "T"
     my_total = TotalSurface(field, surface)
     assert my_total.title == "Total {} surface {}".format(field, surface)
 
 
 class TestCompute:
+    """Test that the total surface export computes the correct value"""
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
 
diff --git a/test/unit/test_exports/test_derived_quantities/test_total_volume.py b/test/unit/test_exports/test_derived_quantities/test_total_volume.py
index 5a12ed43b..8f2e3cfa5 100644
--- a/test/unit/test_exports/test_derived_quantities/test_total_volume.py
+++ b/test/unit/test_exports/test_derived_quantities/test_total_volume.py
@@ -1,22 +1,26 @@
 from festim import TotalVolume
 import fenics as f
+import pytest
 
 
-def test_title_H():
-    volume = 1
-    field = "solute"
-    my_total = TotalVolume(field, volume)
-    assert my_total.title == "Total {} volume {}".format(field, volume)
+@pytest.mark.parametrize("field,volume", [("solute", 1), ("T", 2)])
+def test_title(field, volume):
+    """
+    A simple test to check that the title is set
+    correctly in festim.TotalVolume
 
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+    """
 
-def test_title_T():
-    volume = 2
-    field = "T"
     my_total = TotalVolume(field, volume)
     assert my_total.title == "Total {} volume {}".format(field, volume)
 
 
 class TestCompute:
+    """Test that the total volume export computes the correct value"""
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
 
diff --git a/test/unit/test_exports/test_derived_quantities/test_volume_quantity.py b/test/unit/test_exports/test_derived_quantities/test_volume_quantity.py
index c6c1d212a..4a2cb56b2 100644
--- a/test/unit/test_exports/test_derived_quantities/test_volume_quantity.py
+++ b/test/unit/test_exports/test_derived_quantities/test_volume_quantity.py
@@ -3,8 +3,13 @@
 import pytest
 
 
-def test_wrong_type_for_volume():
-    for volume in [2.1, [0], [1, 1], "coucou", ["coucou"], True]:
-        with pytest.raises(TypeError, match="volume should be an int"):
-            print(volume)
-            F.VolumeQuantity(field=0, volume=volume)
+@pytest.mark.parametrize("volume", [2.1, [0], [1, 1], "coucou", ["coucou"], True])
+def test_wrong_type_for_volume(volume):
+    """
+    Tests that error is raised when the surface volume is set with wrong type
+
+    Args:
+        volume (): wrong type for the surface attribute
+    """
+    with pytest.raises(TypeError, match="volume should be an int"):
+        F.VolumeQuantity(field=0, volume=volume)
diff --git a/test/unit/test_exports/test_trap_density_xdmf_export.py b/test/unit/test_exports/test_trap_density_xdmf_export.py
index 3a0d880c7..a0f7a41d1 100644
--- a/test/unit/test_exports/test_trap_density_xdmf_export.py
+++ b/test/unit/test_exports/test_trap_density_xdmf_export.py
@@ -10,6 +10,9 @@ def test_trap_density_xdmf_export_intergration_with_simultion(tmpdir):
     Creates a festim simulation and exports the trap density as an .xmdf file.
     An equivalent fenics function is created and is compared to that that read
     from the .xdmf file created. Ensures compatability with festim.Simulation()
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     density_expr = 2 + festim.x**2 + 2 * festim.y
 
@@ -58,6 +61,9 @@ def test_trap_density_xdmf_export_write(tmpdir):
     Creates a festim density function and exports as an .xmdf file.
     An equivalent fenics function is created and is compared to that that read
     from the .xdmf file created.
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     # build
     mesh = UnitSquareMesh(30, 30)
@@ -96,6 +102,9 @@ def test_trap_density_xdmf_export_traps_materials_mixed(tmpdir):
     .xmdf file. An equivalent fenics function is created and is compared to that
     that read from the .xdmf file created. Ensures compatability with
     festim.Simulation()
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     density_expr = 2e06 + festim.x**2
 
diff --git a/test/unit/test_helpers.py b/test/unit/test_helpers.py
index a6b4f4352..f1931df44 100644
--- a/test/unit/test_helpers.py
+++ b/test/unit/test_helpers.py
@@ -8,6 +8,7 @@
     t,
 )
 from fenics import Constant, Expression, UserExpression
+import pytest
 
 
 def test_energy_converter():
@@ -19,45 +20,44 @@ def test_energy_converter():
         assert energy_in_eV == expected_value
 
 
-def test_as_constant():
-    assert isinstance(as_constant(3), Constant)
-    assert isinstance(as_constant(3.0), Constant)
-    assert isinstance(as_constant(-2.0), Constant)
-    assert isinstance(as_constant(Constant(2.0)), Constant)
+@pytest.mark.parametrize("constant", [3, 3.0, -2.0, Constant(2.0)])
+def test_as_constant(constant):
+    assert isinstance(as_constant(constant), Constant)
 
 
-def test_as_expression():
-    assert isinstance(as_expression(3 * t), Expression)
-    assert isinstance(as_expression(Expression("2 + x[0]", degree=2)), Expression)
+class CustomExpr(UserExpression):
+    def __init__(self):
+        super().__init__()
 
-    class CustomExpr(UserExpression):
-        def __init__(self):
-            super().__init__()
+    def eval(self, x, values):
+        values[0] = x
 
-        def eval(self, x, values):
-            values[0] = x
 
-    assert isinstance(as_expression(CustomExpr()), UserExpression)
-
-
-def test_as_constant_or_expression():
-    # constants
-    assert isinstance(as_constant_or_expression(3), Constant)
-    assert isinstance(as_constant_or_expression(3.0), Constant)
-    assert isinstance(as_constant_or_expression(-2.0), Constant)
-    assert isinstance(as_constant_or_expression(Constant(2.0)), Constant)
-
-    # expressions
-    assert isinstance(as_constant_or_expression(3 * t), Expression)
-    assert isinstance(
-        as_constant_or_expression(Expression("2 + x[0]", degree=2)), Expression
-    )
-
-    class CustomExpr(UserExpression):
-        def __init__(self):
-            super().__init__()
-
-        def eval(self, x, values):
-            values[0] = x
-
-    assert isinstance(as_constant_or_expression(CustomExpr()), UserExpression)
+@pytest.mark.parametrize(
+    "expression,type",
+    [
+        (3 * t, Expression),
+        (Expression("2 + x[0]", degree=2), Expression),
+        (CustomExpr(), UserExpression),
+    ],
+)
+def test_as_expression(expression, type):
+    assert isinstance(as_expression(expression), type)
+
+
+@pytest.mark.parametrize(
+    "expression,type",
+    [
+        # constants
+        (3, Constant),
+        (3.0, Constant),
+        (-2.0, Constant),
+        (Constant(2.0), Constant),
+        # expressions
+        (3 * t, Expression),
+        (Expression("2 + x[0]", degree=2), Expression),
+        (CustomExpr(), UserExpression),
+    ],
+)
+def test_as_constant_or_expression(expression, type):
+    assert isinstance(as_constant_or_expression(expression), type)
diff --git a/test/unit/test_materials.py b/test/unit/test_materials.py
index 9abe14c1a..90fd0bdcb 100644
--- a/test/unit/test_materials.py
+++ b/test/unit/test_materials.py
@@ -154,6 +154,8 @@ def test_non_matching_properties():
 
 
 class TestCheckBorders:
+    """General test for the check_borders method of the festim.Materials class"""
+
     def test_works(self):
         materials = [
             F.Material(id=1, D_0=None, E_D=None, borders=[0.5, 0.7]),
diff --git a/test/unit/test_mobile.py b/test/unit/test_mobile.py
index c6415842b..adbb0d709 100644
--- a/test/unit/test_mobile.py
+++ b/test/unit/test_mobile.py
@@ -5,6 +5,7 @@
 
 
 def test_mobile_create_diffusion_form():
+    """Tests the create_diffusion_form method of the festim.Mobile class"""
     # build
     Index._globalcount = 8
     mesh = f.UnitIntervalMesh(10)
@@ -37,6 +38,10 @@ def test_mobile_create_diffusion_form():
 
 
 def test_mobile_create_source_form_one_dict():
+    """
+    Tests the create_source_form method of the festim.Mobile class
+    for the case of one festim.Source object
+    """
     # build
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
@@ -58,6 +63,10 @@ def test_mobile_create_source_form_one_dict():
 
 
 def test_mobile_create_source_form_several_sources():
+    """
+    Tests the create_source_form method of the festim.Mobile class
+    for the case of several festim.Source objects
+    """
     # build
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
@@ -84,6 +93,7 @@ def test_mobile_create_source_form_several_sources():
 
 
 def test_mobile_create_form():
+    """Tests the create_form method of the festim.Mobile class"""
     # build
     Index._globalcount = 8
     mesh = f.UnitIntervalMesh(10)
@@ -119,6 +129,8 @@ def add_functions(trap, V, id=1):
 
 
 class TestCreateDiffusionForm:
+    """General test for the create_diffusion_form method of the festim.Mobile class"""
+
     mesh = f.UnitIntervalMesh(10)
     my_mesh = festim.Mesh(mesh)
     my_temp = festim.Temperature(value=100)
@@ -130,6 +142,7 @@ class TestCreateDiffusionForm:
     mat2 = festim.Material(2, D_0=2, E_D=2, S_0=3, E_S=4)
 
     def test_with_traps_transient(self):
+        """Check for the case of transient simulation with traps"""
         # build
         Index._globalcount = 8
         my_mobile = festim.Mobile()
@@ -186,6 +199,7 @@ def test_with_traps_transient(self):
         assert my_mobile.F.equals(expected_form)
 
     def test_with_trap_conglo_transient(self):
+        """Check for the case of transient simulation with a trap conglomerate"""
         # build
         Index._globalcount = 8
         my_mobile = festim.Mobile()
@@ -260,6 +274,10 @@ def test_error_soret_cylindrical_spherical(self):
 
 
 class TestInitialise:
+    """
+    Tests the initialise method of the festim.Mobile class
+    """
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
     u = f.Function(V)
@@ -277,6 +295,9 @@ def test_from_expresion(self):
 
 
 def test_fluxes():
+    """
+    Tests the create_fluxes_form method of the festim.Mobile class
+    """
     Kr_0 = 2
     E_Kr = 3
     order = 2
diff --git a/test/unit/test_radioactive_decay.py b/test/unit/test_radioactive_decay.py
index 374a9a69c..23dd5b443 100644
--- a/test/unit/test_radioactive_decay.py
+++ b/test/unit/test_radioactive_decay.py
@@ -2,30 +2,29 @@
 import pytest
 
 
-def test_init():
-    rd = RadioactiveDecay(0.5, 100)
-    assert rd.decay_constant == 0.5
-    assert rd.volume == 100
-
+class TestRadioactiveDecay:
+    """
+    General test for the festim.RadioactiveDecay class
+    """
 
-def test_decay_constant_setter():
     rd = RadioactiveDecay(0.5, 100)
-    rd.decay_constant = 0.7
-    assert rd.decay_constant == 0.7
-
 
-def test_decay_constant_setter_invalid_type():
-    rd = RadioactiveDecay(0.5, 100)
-    with pytest.raises(TypeError):
-        rd.decay_constant = "invalid"
+    def test_init(self):
+        assert self.rd.decay_constant == 0.5
+        assert self.rd.volume == 100
 
+    def test_decay_constant_setter(self):
+        self.rd.decay_constant = 0.7
+        assert self.rd.decay_constant == 0.7
 
-def test_decay_constant_setter_negative_value():
-    rd = RadioactiveDecay(0.5, 100)
-    with pytest.raises(ValueError):
-        rd.decay_constant = -0.5
+    def test_decay_constant_setter_invalid_type(self):
+        with pytest.raises(TypeError):
+            self.rd.decay_constant = "invalid"
 
+    def test_decay_constant_setter_negative_value(self):
+        with pytest.raises(ValueError):
+            self.rd.decay_constant = -0.5
 
-def test_form():
-    rd = RadioactiveDecay(0.5, 100)
-    assert rd.form(200) == -100
+    def test_form(self):
+        self.rd.decay_constant = 0.5
+        assert self.rd.form(200) == -100
diff --git a/test/unit/test_temperature.py b/test/unit/test_temperature.py
index c1e975858..ed36020ef 100644
--- a/test/unit/test_temperature.py
+++ b/test/unit/test_temperature.py
@@ -100,6 +100,9 @@ def test_heat_transfer_create_functions_transient(tmpdir):
     Creates a function, writes it to an XDMF file, then a HeatTransferProblem
     class is created from this file and the error norm between the written and
     read fuctions is computed to ensure they are the same.
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     # create function to be comapared
     mesh = fenics.UnitIntervalMesh(10)
@@ -135,6 +138,9 @@ def test_temperature_from_xdmf_create_functions(tmpdir):
     Creates a function, writes it to an XDMF file, then a TemperatureFromXDMF
     class is created from this file and the error norm between the written and
     read fuctions is computed to ensure they are the same.
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     # create function to be comapared
     mesh = fenics.UnitSquareMesh(10, 10)
@@ -160,6 +166,9 @@ def test_temperature_from_xdmf_label_checker(tmpdir):
     """Test for the label check test within the TemperatureFromXDMF class,
     ensures that a ValueError is raised when reading a file with an
     incorrect label.
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     # create function to be written
     mesh = fenics.UnitSquareMesh(10, 10)
@@ -177,8 +186,13 @@ def test_temperature_from_xdmf_label_checker(tmpdir):
 
 
 def test_temperature_from_xdmf_transient_case(tmpdir):
-    """Test that the TemperatureFromXdmf class works in a transient
-    h transport case"""
+    """
+    Test that the TemperatureFromXdmf class works in a transient
+    h transport case
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
+    """
     # create temperature field xdmf
     my_model = festim.Simulation(log_level=20)
     my_model.mesh = festim.MeshFromVertices(vertices=np.linspace(0, 1, num=100))
@@ -209,6 +223,9 @@ def test_temperature_from_xdmf(tmpdir):
     """
     Tests that .is_steady_state() can be run for
     a TemperatureFromXDMF object
+
+    Args:
+        tmpdir (os.PathLike): path to the pytest temporary folder
     """
     mesh = fenics.UnitSquareMesh(10, 10)
     V = fenics.FunctionSpace(mesh, "CG", 1)
diff --git a/test/unit/test_theta.py b/test/unit/test_theta.py
index bcb434a55..a00123793 100644
--- a/test/unit/test_theta.py
+++ b/test/unit/test_theta.py
@@ -5,6 +5,10 @@
 
 
 class TestInitialise:
+    """
+    General test for the initialise method of the festim.Theta class
+    """
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
     u = f.Function(V)
@@ -13,6 +17,7 @@ class TestInitialise:
     T.create_functions(festim.Mesh(mesh))
 
     def test_from_expresion_chemical_pot(self):
+        """test initialisation from an expression"""
         my_theta = festim.Theta()
         S = f.interpolate(f.Constant(2), self.V)
         my_theta.materials = festim.Materials([festim.Material(1, 1, 0, S_0=2, E_S=0)])
@@ -32,6 +37,10 @@ def test_from_expresion_chemical_pot(self):
 
 
 class TestCreateDiffusionForm:
+    """
+    General test for the create_diffusion_form method of the festim.Theta class
+    """
+
     mesh = f.UnitIntervalMesh(10)
     my_mesh = festim.Mesh(mesh)
     my_temp = festim.Temperature(value=100)
@@ -41,6 +50,7 @@ class TestCreateDiffusionForm:
     V = f.FunctionSpace(my_mesh.mesh, "CG", 1)
 
     def test_sieverts(self):
+        """test with the Sievert solubility law"""
         # build
         mat1 = festim.Material(1, D_0=1, E_D=1, S_0=2, E_S=3, solubility_law="sievert")
         Index._globalcount = 8
@@ -72,6 +82,7 @@ def test_sieverts(self):
         assert my_theta.F.equals(expected_form)
 
     def test_henry(self):
+        """test with the Henry solubility law"""
         # build
         mat2 = festim.Material(2, D_0=2, E_D=2, S_0=3, E_S=4, solubility_law="henry")
 
@@ -105,6 +116,10 @@ def test_henry(self):
 
 
 def test_get_concentration_for_a_given_material():
+    """
+    Test that the get_concentration_for_a_given_material
+    method works correctly
+    """
     # build
     S_0 = 2
     E_S = 0.5
diff --git a/test/unit/test_traps/test_extrinsic_trap.py b/test/unit/test_traps/test_extrinsic_trap.py
index 59bdb9ff4..e36ba4010 100644
--- a/test/unit/test_traps/test_extrinsic_trap.py
+++ b/test/unit/test_traps/test_extrinsic_trap.py
@@ -3,6 +3,10 @@
 
 
 class TestExtrinsicTrap:
+    """
+    General test for the ExtrinsicTrap class
+    """
+
     my_trap = festim.ExtrinsicTrap(
         1,
         1,
@@ -28,6 +32,9 @@ class TestExtrinsicTrap:
     dt = festim.Stepsize(initial_value=1)
 
     def test_that_form_parameters_are_expressions(self):
+        """
+        Checks that the attributes are of correct types
+        """
         prms = [
             self.my_trap.phi_0,
             self.my_trap.n_amax,
@@ -41,6 +48,10 @@ def test_that_form_parameters_are_expressions(self):
             assert isinstance(prm, (f.Expression, f.Constant))
 
     def test_create_form_density(self):
+        """
+        Checks that the forumlation produced by the create_form_density
+        function produces the expected formulation
+        """
         density = self.my_trap.density[0]
         T = self.my_temp
         expected_form = (
diff --git a/test/unit/test_traps/test_trap.py b/test/unit/test_traps/test_trap.py
index b5ed81860..8711ed080 100644
--- a/test/unit/test_traps/test_trap.py
+++ b/test/unit/test_traps/test_trap.py
@@ -9,49 +9,54 @@ def add_functions(trap, V, id=1):
     trap.test_function = f.TestFunction(V)
 
 
-def test_error_wrong_type_material():
-    """Checks that an error is raised when the wrong type is given to
+@pytest.mark.parametrize("mats", [True, [True, "mat_name"], 1, [1, 2]])
+def test_error_wrong_type_material(mats):
+    """
+    Checks that an error is raised when the wrong type is given to
     materials
+
+    Args:
+        mats (): wrong type objects for the Trap.materials attribute
     """
     msg = "Accepted types for materials are str or festim.Material"
     with pytest.raises(TypeError, match=msg):
-        festim.Trap(1, 1, 1, 1, materials=True, density=1)
+        festim.Trap(1, 1, 1, 1, materials=mats, density=1)
 
-    with pytest.raises(TypeError, match=msg):
-        festim.Trap(1, 1, 1, 1, materials=[True, "mat_name"], density=1)
-
-    with pytest.raises(TypeError, match=msg):
-        festim.Trap(1, 1, 1, 1, materials=1, density=1)
 
-    with pytest.raises(TypeError, match=msg):
-        festim.Trap(1, 1, 1, 1, materials=[1, 2], density=1)
-
-
-def test_error_if_duplicate_material():
+class TestDuplicateMaterial:
     mat1 = festim.Material(1, D_0=1, E_D=0, name="name1")
     mat2 = festim.Material(2, D_0=1, E_D=0, name="name2")
 
     materials = festim.Materials([mat1, mat2])
-    with pytest.raises(ValueError, match="Duplicate materials in trap"):
-        festim.Trap(1, 1, 1, 1, [mat1, mat1], 1).make_materials(materials)
-
-    with pytest.raises(ValueError, match="Duplicate materials in trap"):
-        festim.Trap(1, 1, 1, 1, ["name1", "name1", mat2], 1).make_materials(materials)
 
-    with pytest.raises(ValueError, match="Duplicate materials in trap"):
-        festim.Trap(1, 1, 1, 1, ["name1", mat1, mat1], 1).make_materials(materials)
-
-    with pytest.raises(ValueError, match="Duplicate materials in trap"):
-        festim.Trap(1, 1, 1, 1, ["name2", mat2], 1).make_materials(materials)
-
-    with pytest.raises(ValueError, match="Duplicate materials in trap"):
-        festim.Trap(1, 1, 1, 1, [mat2, mat2], 1).make_materials(materials)
+    @pytest.mark.parametrize(
+        "mats_list",
+        [
+            [mat1, mat1],
+            ["name1", "name1", mat2],
+            ["name1", mat1, mat1],
+            ["name2", mat2],
+            [mat2, mat2],
+            ["name1", mat1],
+        ],
+    )
+    def test_error_if_duplicate_material(self, mats_list):
+        """
+        Checks that an error is raised when there are duplicates in
+        the materials attribute of the festim.Trap class
 
-    with pytest.raises(ValueError, match="Duplicate materials in trap"):
-        festim.Trap(1, 1, 1, 1, ["name1", mat1], 1).make_materials(materials)
+        Args:
+            mats_list (): list containing objects with duplicate names
+        """
+        with pytest.raises(ValueError, match="Duplicate materials in trap"):
+            festim.Trap(1, 1, 1, 1, mats_list, 1).make_materials(self.materials)
 
 
 class TestCreateTrappingForm:
+    """
+    General test for the create_trapping_form method of the festim.Trap class
+    """
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
     my_mobile = festim.Mobile()
@@ -68,6 +73,10 @@ class TestCreateTrappingForm:
     mat2 = festim.Material(2, D_0=2, E_D=2, S_0=3, E_S=4, name="mat2")
 
     def test_steady_state(self):
+        """
+        Test that create_trapping_form creates the correct formulation in
+        the steady-state case
+        """
         # build
         my_trap = festim.Trap(
             k_0=1, E_k=2, p_0=3, E_p=4, materials=self.mat1, density=1 + festim.x
@@ -105,6 +114,10 @@ def test_steady_state(self):
         assert my_trap.F_trapping.equals(expected_form)
 
     def test_transient(self):
+        """
+        Test that create_trapping_form creates the correct formulation in
+        the transient case
+        """
         # build
         my_trap = festim.Trap(
             k_0=1, E_k=2, p_0=3, E_p=4, materials=self.mat1, density=1 + festim.x
@@ -150,6 +163,10 @@ def test_transient(self):
         assert my_trap.F_trapping.equals(expected_form)
 
     def test_chemical_potential(self):
+        """
+        Test that create_trapping_form creates the correct formulation
+        with chemical potential conservation
+        """
         # build
         my_trap = festim.Trap(
             k_0=1, E_k=2, p_0=3, E_p=4, materials=self.mat1, density=1 + festim.x
@@ -197,6 +214,10 @@ def test_chemical_potential(self):
         assert my_trap.F_trapping.equals(expected_form)
 
     def test_2_materials(self):
+        """
+        Test that create_trapping_form creates the correct formulation
+        with two materials
+        """
         # build
         my_trap = festim.Trap(
             k_0=1,
@@ -242,6 +263,10 @@ def test_2_materials(self):
         assert my_trap.F_trapping.equals(expected_form)
 
     def test_multi_parameters_trap(self):
+        """
+        Test that create_trapping_form creates the correct formulation
+        with a trap conglomerate
+        """
         # build
         my_trap = festim.Trap(
             k_0=[1, 2],
@@ -286,6 +311,10 @@ def test_multi_parameters_trap(self):
         assert my_trap.F_trapping.equals(expected_form)
 
     def test_steady_state_trap_not_defined_everywhere(self):
+        """
+        Test that create_trapping_form creates the correct formulation
+        in steady-state with a trap distribution
+        """
         # build
         my_trap = festim.Trap(
             k_0=1, E_k=2, p_0=3, E_p=4, materials=self.mat1, density=1 + festim.x
@@ -358,7 +387,8 @@ def test_expression_as_density(self):
         assert my_trap.F_trapping.equals(expected_form)
 
     def test_user_expression_as_density(self):
-        """Test that create_trapping_form creates the correct formulation when
+        """
+        Test that create_trapping_form creates the correct formulation when
         a fenics.UserExpression is given as density
         """
 
@@ -481,6 +511,10 @@ def test_2_materials_names_and_object(self):
 
 
 class TestCreateSourceForm:
+    """
+    General test for the create_source_form method of the festim.Trap class
+    """
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
     my_mobile = festim.Mobile()
@@ -516,6 +550,10 @@ def test(self):
 
 
 class TestCreateForm:
+    """
+    General test for the create_form method of the festim.Trap class
+    """
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
     my_mobile = festim.Mobile()
@@ -540,6 +578,7 @@ def test_form_is_zero_by_default(self):
         assert my_trap.F == 0
 
     def test_1_mat_steady(self):
+        """Tests the case of one material in steady-state"""
         # build
         my_trap = festim.Trap(1, 1, 1, 1, materials=self.mat1, density=1)
         add_functions(my_trap, self.V, id=1)
@@ -558,6 +597,7 @@ def test_1_mat_steady(self):
         assert my_trap.F.equals(expected_form)
 
     def test_1_mat_transient(self):
+        """Tests the case of one material in transient"""
         # build
         my_trap = festim.Trap(1, 1, 1, 1, materials=self.mat1, density=1)
         add_functions(my_trap, self.V, id=1)
@@ -578,6 +618,7 @@ def test_1_mat_transient(self):
         assert my_trap.F.equals(expected_form)
 
     def test_2_mats_transient(self):
+        """Tests the case of two materials in transient"""
         # build
         my_trap = festim.Trap(1, 1, 1, 1, materials=[self.mat1, self.mat2], density=1)
         add_functions(my_trap, self.V, id=1)
@@ -598,6 +639,7 @@ def test_2_mats_transient(self):
         assert my_trap.F.equals(expected_form)
 
     def test_1_mat_and_source(self):
+        """Tests the case of one material with a source term"""
         # build
         my_trap = festim.Trap(1, 1, 1, 1, materials=self.mat2, density=1)
         my_trap.sources = [festim.Source(1 + festim.x + festim.y, volume=1, field="1")]
diff --git a/test/unit/test_traps/test_traps.py b/test/unit/test_traps/test_traps.py
index 267fd8f78..1c1cec14a 100644
--- a/test/unit/test_traps/test_traps.py
+++ b/test/unit/test_traps/test_traps.py
@@ -21,6 +21,10 @@ def add_functions(trap, V, id=1):
 
 
 class TestCreateTrappingForms:
+    """
+    General test for the create_forms method of the festim.Traps class
+    """
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
     my_mobile = festim.Mobile()
@@ -45,6 +49,7 @@ class TestCreateTrappingForms:
     add_functions(trap2, V, id=2)
 
     def test_one_trap_steady_state(self):
+        """Tests the case of one trap in steady-state"""
         my_traps = festim.Traps([self.trap1])
 
         my_traps.create_forms(self.my_mobile, self.my_mats, self.my_temp, self.dx)
@@ -53,6 +58,7 @@ def test_one_trap_steady_state(self):
             assert trap.F is not None
 
     def test_one_trap_transient(self):
+        """Tests the case of one trap in transient"""
         my_traps = festim.Traps([self.trap1])
 
         my_traps.create_forms(
@@ -63,6 +69,7 @@ def test_one_trap_transient(self):
             assert trap.F is not None
 
     def test_two_traps_transient(self):
+        """Tests the case of two traps in transient"""
         my_traps = festim.Traps([self.trap1, self.trap2])
 
         my_traps.create_forms(
@@ -74,6 +81,10 @@ def test_two_traps_transient(self):
 
 
 class TestGetTrap:
+    """
+    General test for the get_trap method of the festim.Traps class
+    """
+
     mesh = f.UnitIntervalMesh(10)
     V = f.FunctionSpace(mesh, "P", 1)
     my_mobile = festim.Mobile()