diff --git a/src/porepy/models/contact_mechanics_biot_model.py b/src/porepy/models/contact_mechanics_biot_model.py
index 19ac384d2d..8293281ec1 100644
--- a/src/porepy/models/contact_mechanics_biot_model.py
+++ b/src/porepy/models/contact_mechanics_biot_model.py
@@ -133,7 +133,7 @@ def set_mechanics_parameters(self):
                     {"friction_coefficient": friction, "time_step": self.time_step},
                 )
         # Should we keep this, @EK?
-        for e, d in gb.edges():
+        for _, d in gb.edges():
             mg = d["mortar_grid"]
 
             # Parameters for the surface diffusion.
@@ -181,10 +181,13 @@ def set_scalar_parameters(self):
 
         # Assign diffusivity in the normal direction of the fractures.
         for e, data_edge in self.gb.edges():
-            g1, g2 = self.gb.nodes_of_edge(e)
+            g1, _ = self.gb.nodes_of_edge(e)
+            
             a = self.compute_aperture(g1)
             mg = data_edge["mortar_grid"]
+            
             normal_diffusivity = 2 / kappa * mg.slave_to_mortar_int() * a
+            
             data_edge = pp.initialize_data(
                 e,
                 data_edge,
diff --git a/src/porepy/models/contact_mechanics_model.py b/src/porepy/models/contact_mechanics_model.py
index fffe39b97a..c03c391d11 100644
--- a/src/porepy/models/contact_mechanics_model.py
+++ b/src/porepy/models/contact_mechanics_model.py
@@ -147,7 +147,7 @@ def set_parameters(self):
                     {"friction_coefficient": friction},
                 )
         # Should we keep this, @EK?
-        for e, d in gb.edges():
+        for _, d in gb.edges():
             mg = d["mortar_grid"]
 
             # Parameters for the surface diffusion.
@@ -240,7 +240,7 @@ def initial_condition(self):
                 state = {}
             pp.set_state(d, state)
 
-        for e, d in self.gb.edges():
+        for _, d in self.gb.edges():
             mg = d["mortar_grid"]
 
             if mg.dim == self.Nd - 1:
@@ -267,7 +267,7 @@ def extract_iterate(self, assembler, solution_vector):
 
         Returns:
             (np.array): displacement solution vector for the Nd grid.
-            
+
         """
         dof = np.cumsum(np.append(0, np.asarray(assembler.full_dof)))
 
diff --git a/src/porepy/numerics/fv/biot.py b/src/porepy/numerics/fv/biot.py
index faee14700f..cd0bb0c138 100644
--- a/src/porepy/numerics/fv/biot.py
+++ b/src/porepy/numerics/fv/biot.py
@@ -795,7 +795,7 @@ def compute_stress(self, g, u, data):
 class GradP():
     """ Class for the pressure gradient term of the Biot equation.
     """
-    
+
     def __init__(self, keyword):
         """ Set the discretization, with the keyword used for storing various
         information associated with the discretization.
@@ -814,7 +814,7 @@ def _key(self):
             String, on the form self.keyword + '_'.
 
         """
-        return self.keyword + "_"    
+        return self.keyword + "_"
 
     def ndof(self, g):
         """ Return the number of degrees of freedom associated to the method.
diff --git a/src/porepy/numerics/fv/mpsa.py b/src/porepy/numerics/fv/mpsa.py
index 6bc69a1e63..9e1db02415 100644
--- a/src/porepy/numerics/fv/mpsa.py
+++ b/src/porepy/numerics/fv/mpsa.py
@@ -20,8 +20,8 @@
 
 
 class Mpsa():
-    
-    
+
+
     def __init__(self, keyword):
         """ Set the discretization, with the keyword used for storing various
         information associated with the discretization.
@@ -1132,7 +1132,7 @@ def _mpsa_local(
     if g.dim == 2:
         g = g.copy()
 
-        cell_centers, face_normals, face_centers, R, _, nodes = pp.map_geometry.map_grid(
+        cell_centers, face_normals, face_centers, _, _, nodes = pp.map_geometry.map_grid(
             g
         )
         g.cell_centers = cell_centers
@@ -1762,7 +1762,7 @@ def _inverse_gradient(
         * pp.fvutils.invert_diagonal_blocks(grad, size_of_blocks, method=inverter)
         * rows2blk_diag
     )
-    logger.debug("max igrad: ", np.max(np.abs(igrad)))
+    logger.debug("max igrad: " + str(np.max(np.abs(igrad))))
     return igrad
 
 
@@ -2277,7 +2277,7 @@ def _eliminate_ncasym_neumann(
     dof_elim = subfno_nd.ravel("C")[remove_singular]
     # and eliminate the rows corresponding to these subfaces
     pp.utils.sparse_mat.zero_rows(ncasym, dof_elim)
-    logger.debug("number of ncasym eliminated: ", np.sum(dof_elim.size))
+    logger.debug("number of ncasym eliminated: " + str(np.sum(dof_elim.size)))
     ## the following is some code to enforce symmetric G. Comment for now
     # # Find the equations for the x-values
     # x_row = np.arange(0, round(ncasym.shape[0]/nd))
diff --git a/src/porepy/numerics/fv/upwind.py b/src/porepy/numerics/fv/upwind.py
index 94dc72631b..2075435db4 100644
--- a/src/porepy/numerics/fv/upwind.py
+++ b/src/porepy/numerics/fv/upwind.py
@@ -78,10 +78,6 @@ def assemble_matrix_rhs(self, g, data):
                 conditions. The size of the vector will depend on the discretization.
 
         """
-        matrix_dictionary = data[pp.DISCRETIZATION_MATRICES][self.keyword]
-        if not self.matrix_keyword in matrix_dictionary.keys():
-            self.discretize(g, data)
-
         return self.assemble_matrix(g, data), self.assemble_rhs(g, data)
 
     def assemble_matrix(self, g, data):
diff --git a/src/porepy/numerics/interface_laws/contact_mechanics_interface_laws.py b/src/porepy/numerics/interface_laws/contact_mechanics_interface_laws.py
index 0dd02a2065..95bfbd532e 100644
--- a/src/porepy/numerics/interface_laws/contact_mechanics_interface_laws.py
+++ b/src/porepy/numerics/interface_laws/contact_mechanics_interface_laws.py
@@ -94,7 +94,7 @@ def discretize(self, g_h, g_l, data_h, data_l, data_edge):
             # List of surface diffusion discretizations - one per side.
             A_list = []
 
-            for proj, side_grid in mg.project_to_side_grids():
+            for _, side_grid in mg.project_to_side_grids():
 
                 unity = np.ones(side_grid.num_cells)
 
@@ -547,11 +547,12 @@ class FractureScalarToForceBalance:
 
     For the contact mechanics, we only want to consider the _contact_ traction. Thus, we
     have to subtract the pressure contribution, i.e.
-    
+
         \lambda_contact - p_check I \dot n = boundary_traction_hat,
 
     since the full tractions experienced by a fracture surface are the sum of the
     contact forces and the fracture pressure force.
+
     """
 
     def __init__(self, discr_master, discr_slave):
diff --git a/src/porepy/params/bc.py b/src/porepy/params/bc.py
index 39737112c1..c372fc68a1 100644
--- a/src/porepy/params/bc.py
+++ b/src/porepy/params/bc.py
@@ -54,7 +54,7 @@ class BoundaryCondition(AbstractBoundaryCondition):
             well as Dirichlet faces.
         is_dir (np.ndarary, boolean, size g.num_faces): Element i is true if
             face i has been assigned a Neumann condition.
-        is_rob (np.ndarray, boolean, size g.num_faces): Element i is true if 
+        is_rob (np.ndarray, boolean, size g.num_faces): Element i is true if
             face i has been assigned a Robin condition.
 
     """
@@ -275,7 +275,7 @@ class BoundaryConditionVectorial(AbstractBoundaryCondition):
             well as Dirichlet faces.
         is_dir (np.ndarary, boolean, size g.dim x g.num_faces): Element i is true if
             face i has been assigned a Neumann condition.
-        is_rob (np.ndarray, boolean, size g.dim x g.num_faces): Element i is true if 
+        is_rob (np.ndarray, boolean, size g.dim x g.num_faces): Element i is true if
             face i has been assigned a Robin condition.
 
     """
diff --git a/src/porepy/utils/derived_discretizations/implicit_euler.py b/src/porepy/utils/derived_discretizations/implicit_euler.py
index a487379b74..d872f3364b 100644
--- a/src/porepy/utils/derived_discretizations/implicit_euler.py
+++ b/src/porepy/utils/derived_discretizations/implicit_euler.py
@@ -132,7 +132,7 @@ class ImplicitUpwind(pp.Upwind):
     Multiply all contributions by the time step and advection weight.
     """
 
-    def assemble_matrix_rhs(self, g, data, d_name="darcy_flux"):
+    def assemble_matrix_rhs(self, g, data):
         if g.dim == 0:
             data["flow_faces"] = sps.csr_matrix([0.0])
             return sps.csr_matrix([0.0]), np.array([0.0])
@@ -140,7 +140,7 @@ def assemble_matrix_rhs(self, g, data, d_name="darcy_flux"):
         parameter_dictionary = data[pp.PARAMETERS][self.keyword]
         dt = parameter_dictionary["time_step"]
         w = parameter_dictionary["advection_weight"] * dt
-        a, b = super().assemble_matrix_rhs(g, data, d_name)
+        a, b = super().assemble_matrix_rhs(g, data)
         a = a * w
         b = b * w
         return a, b
diff --git a/test/integration/test_contact_mechanics.py b/test/integration/test_contact_mechanics.py
index 5eef833954..f0ccb8f428 100644
--- a/test/integration/test_contact_mechanics.py
+++ b/test/integration/test_contact_mechanics.py
@@ -3,7 +3,6 @@
 """
 import numpy as np
 import unittest
-import scipy.sparse.linalg as spla
 
 import porepy as pp
 import porepy.models.contact_mechanics_model as model
diff --git a/test/integration/test_contact_mechanics_biot.py b/test/integration/test_contact_mechanics_biot.py
index e787ae2ad9..336d7487c6 100644
--- a/test/integration/test_contact_mechanics_biot.py
+++ b/test/integration/test_contact_mechanics_biot.py
@@ -8,9 +8,7 @@
 test, please refer to test_contact_mechanics.
 """
 import numpy as np
-import scipy.sparse as sps
 import unittest
-import scipy.sparse.linalg as spla
 
 import porepy as pp
 import porepy.models.contact_mechanics_biot_model as model