Documentation and tolerance for the contact models

src/porepy/models/contact_mechanics_biot_model.py

@@ -1,12 +1,14 @@
 """
-This is a setup class for solving the biot equations with contact between the fractures.
+This is a setup class for solving the Biot equations with contact between the fractures.
 
-The domain $[0, 2]\times[0, 1]$ with six fractures. We do not consider any fluid, and
-solve only for the linear elasticity coupled to the contact
+The class ContactMechanicsBiot inherits from ContactMechanics, which is a model for
+the purely mechanical problem with contact conditions on the fractures. Here, we add
+expand to a model where the displacement solution is coupled to a scalar variable, e.g.
+pressure (Biot equations) or temperature. Parameters, variables and discretizations are
+set in the model class, and the problem may be solved using run_biot.
 
 NOTE: This module should be considered an experimental feature, which will likely
 undergo major changes (or be deleted).
-
 """
 import numpy as np
 import porepy as pp
@@ -35,9 +37,6 @@ def __init__(self, mesh_args, folder_name):
         self.time_step = 1 * self.length_scale ** 2
         self.end_time = self.time_step * 1
 
-        # Initial scalar value
-        self.initial_scalar = 0
-
         # Whether or not to subtract the fracture pressure contribution for the contact
         # traction. This should be done if the scalar variable is pressure, but not for
         # temperature. See assign_discretizations
@@ -82,7 +81,6 @@ def set_parameters(self):
         self.set_scalar_parameters()
         self.set_mechanics_parameters()
 
-
     def set_mechanics_parameters(self):
         """
         Set the parameters for the simulation.
@@ -371,7 +369,7 @@ def initial_condition(self):
 
         for g, d in self.gb:
             # Initial value for the scalar variable.
-            initial_scalar_value = self.initial_scalar * np.ones(g.num_cells)
+            initial_scalar_value = 1.0 * np.ones(g.num_cells)
             d[pp.STATE].update({self.scalar_variable: initial_scalar_value})
             if g.dim == self.Nd:
                 bc_values = d[pp.PARAMETERS][self.mechanics_parameter_key]["bc_values"]
@@ -384,7 +382,8 @@ def export_step(self):
     def export_pvd(self):
         pass
 
-def run_biot(setup, atol=1e-10):
+
+def run_biot(setup, newton_tol=1e-10):
     """
     Function for solving the time dependent Biot equations with a non-linear Coulomb
     contact condition on the fractures.
@@ -409,6 +408,7 @@ def run_biot(setup, atol=1e-10):
             and attributes:
                 end_time: End time time of simulation.
                 time_step: Time step size
+        newton_tol: Tolerance for the Newton solver, see contact_mechanics_model.
     """
     if "gb" not in setup.__dict__:
         setup.create_grid()
@@ -436,7 +436,6 @@ def run_biot(setup, atol=1e-10):
     dt = setup.time_step
     t_end = setup.end_time
     k = 0
-    setup.export_step()
     while setup.time < t_end:
         setup.time += dt
         k += 1
@@ -445,18 +444,19 @@ def run_biot(setup, atol=1e-10):
         # Prepare for Newton
         counter_newton = 0
         converged_newton = False
-        max_newton = 20
+        max_newton = 15
         newton_errors = []
         while counter_newton <= max_newton and not converged_newton:
             print("Newton iteration number: ", counter_newton, "/", max_newton)
             # One Newton iteration:
             sol, u, error, converged_newton = pp.models.contact_mechanics_model.newton_iteration(
-                assembler, setup, u
+                assembler, setup, u, tol=newton_tol
             )
             counter_newton += 1
             newton_errors.append(error)
         # Prepare for next time step
         assembler.distribute_variable(sol)
         setup.export_step()
         errors.append(newton_errors)
+    setup.newton_errors = errors
     setup.export_pvd()

src/porepy/models/contact_mechanics_model.py

@@ -270,37 +270,45 @@ def extract_iterate(self, assembler, solution_vector):
             (np.array): displacement solution vector for the Nd grid.
 
         """
-        dof = np.cumsum(np.append(0, np.asarray(assembler.full_dof)))
+        variable_names = []
+        for pair in assembler.block_dof.keys():
+            variable_names.append(pair[1])
 
-        for pair, bi in assembler.block_dof.items():
-            g = pair[0]
-            name = pair[1]
-            # Identify edges, and update the mortar displacement iterate
-            if isinstance(g, tuple):
-                if name == self.mortar_displacement_variable:
-                    mortar_u = solution_vector[dof[bi] : dof[bi + 1]]
-                    data = self.gb.edge_props(g)
-                    data[pp.STATE]["previous_iterate"][
-                        self.mortar_displacement_variable
-                    ] = mortar_u
-                continue
-            else:
-                # g is a node (not edge)
+        dof = np.cumsum(np.append(0, np.asarray(assembler.full_dof)))
 
-                # For the fractures, update the contact force
-                if g.dim < self.gb.dim_max():
-                    if name == self.contact_traction_variable:
-                        contact = solution_vector[dof[bi] : dof[bi + 1]]
-                        data = self.gb.node_props(g)
+        for var_name in set(variable_names):
+            for pair, bi in assembler.block_dof.items():
+                g = pair[0]
+                name = pair[1]
+                if name != var_name:
+                    continue
+                if isinstance(g, tuple):
+                    # This is really an edge
+                    if name == self.mortar_displacement_variable:
+                        mortar_u = solution_vector[dof[bi] : dof[bi + 1]]
+                        data = self.gb.edge_props(g)
                         data[pp.STATE]["previous_iterate"][
-                            self.contact_traction_variable
-                        ] = contact
-
+                            self.mortar_displacement_variable
+                        ] = mortar_u.copy()
                 else:
-                    # Only need the displacements for Nd
-                    if name != self.displacement_variable:
-                        continue
-                    u = solution_vector[dof[bi] : dof[bi + 1]]
+                    data = self.gb.node_props(g)
+
+                    # g is a node (not edge)
+
+                    # For the fractures, update the contact force
+                    if g.dim < self.Nd:
+                        if name == self.contact_traction_variable:
+                            contact = solution_vector[dof[bi] : dof[bi + 1]]
+                            data = self.gb.node_props(g)
+                            data[pp.STATE]["previous_iterate"][
+                                self.contact_traction_variable
+                            ] = contact.copy()
+
+                    else:
+                        # Only need the displacements for Nd
+                        if name != self.displacement_variable:
+                            continue
+                        u = solution_vector[dof[bi] : dof[bi + 1]]
         return u
 
     def reconstruct_local_displacement_jump(self, data_edge):
@@ -410,7 +418,7 @@ def run_mechanics(setup):
     assembler.distribute_variable(sol)
 
 
-def newton_iteration(assembler, setup, u0, solver=None):
+def newton_iteration(assembler, setup, u0, tol=1e-14, solver=None):
     converged = False
     # @EK! If this is to work for both mechanics and biot, we probably need to pass the solver to this method.
     g_max = setup.gb.grids_of_dimension(setup.Nd)[0]
@@ -420,6 +428,12 @@ def newton_iteration(assembler, setup, u0, solver=None):
 
     # Assemble and solve
     A, b = assembler.assemble_matrix_rhs()
+    print("max A: {0:.2e}".format(np.max(np.abs(A))))
+    print(
+        "max: {0:.2e} and min: {1:.2e} A sum: ".format(
+            np.max(np.sum(np.abs(A), axis=1)), np.min(np.sum(np.abs(A), axis=1))
+        )
+    )
 
     if solver is None:
         sol = sps.linalg.spsolve(A, b)
@@ -435,15 +449,21 @@ def newton_iteration(assembler, setup, u0, solver=None):
     iterate_difference = l2_norm_cell(g_max, u1, u0)
 
     # The if is intended to avoid division through zero
-    if solution_norm < 1e-12 and iterate_difference < 1e-12:
+    if solution_norm < tol and iterate_difference < tol:
         converged = True
         error = np.sum((u1 - u0) ** 2)
+
     else:
-        if iterate_difference / solution_norm < 1e-10:
+        if iterate_difference / solution_norm < tol:
             converged = True
         error = np.sum((u1 - u0) ** 2) / np.sum(u1 ** 2)
 
-    print("Error: ", error)
+    print(
+        "Error, solution norm and iterate_difference/solution_norm: ",
+        error,
+        solution_norm,
+        iterate_difference / solution_norm,
+    )
 
     return sol, u1, error, converged