Add exporting to contact mechanics biot model

src/porepy/models/contact_mechanics_biot_model.py

@@ -77,11 +77,11 @@ def compute_aperture(self, g):
 
     def set_parameters(self):
         """
-        Set the parameters for the simulation. The stress is given in GPa.
+        Set the parameters for the simulation.
         """
-
-        self.set_mechanics_parameters()
         self.set_scalar_parameters()
+        self.set_mechanics_parameters()
+
 
     def set_mechanics_parameters(self):
         """
@@ -92,18 +92,12 @@ def set_mechanics_parameters(self):
         for g, d in gb:
             if g.dim == self.Nd:
                 # Rock parameters
-                lam = np.ones(g.num_cells)
-                mu = np.ones(g.num_cells)
+                lam = np.ones(g.num_cells) / self.scalar_scale
+                mu = np.ones(g.num_cells) / self.scalar_scale
                 C = pp.FourthOrderTensor(g.dim, mu, lam)
 
                 # Define boundary condition
                 bc = self.bc_type_mechanics(g)
-                # Default internal BC is Neumann. We change to Dirichlet for the contact
-                # problem. I.e., the mortar variable represents the displacement on the
-                # fracture faces.
-                frac_face = g.tags["fracture_faces"]
-                bc.is_neu[:, frac_face] = False
-                bc.is_dir[:, frac_face] = True
                 # BC and source values
                 bc_val = self.bc_values_mechanics(g)
                 source_val = self.source_mechanics(g)
@@ -384,6 +378,11 @@ def initial_condition(self):
                 mech_dict = {"bc_values": bc_values}
                 d[pp.STATE].update({self.mechanics_parameter_key: mech_dict})
 
+    def export_step(self):
+        pass
+
+    def export_pvd(self):
+        pass
 
 def run_biot(setup, atol=1e-10):
     """
@@ -428,25 +427,25 @@ def run_biot(setup, atol=1e-10):
     assembler = pp.Assembler(gb)
     u = d_max[pp.STATE][setup.displacement_variable]
 
+    setup.export_step()
+
     # Discretize with the biot class
     setup.discretize_biot(gb)
     # Prepare for the time loop
     errors = []
-    t = getattr(setup, "time", 0)
     dt = setup.time_step
     t_end = setup.end_time
     k = 0
-    times = [t]
-    while t < t_end:
-        t += dt
+    setup.export_step()
+    while setup.time < t_end:
+        setup.time += dt
         k += 1
         print("Time step: ", k, "/", int(np.ceil(t_end / dt)))
 
-        times.append(t)
         # Prepare for Newton
         counter_newton = 0
         converged_newton = False
-        max_newton = 10
+        max_newton = 20
         newton_errors = []
         while counter_newton <= max_newton and not converged_newton:
             print("Newton iteration number: ", counter_newton, "/", max_newton)
@@ -458,4 +457,6 @@ def run_biot(setup, atol=1e-10):
             newton_errors.append(error)
         # Prepare for next time step
         assembler.distribute_variable(sol)
+        setup.export_step()
         errors.append(newton_errors)
+    setup.export_pvd()

src/porepy/models/contact_mechanics_model.py

@@ -48,7 +48,7 @@ def create_grid(self):
             gb (pp.GridBucket): The produced grid bucket.
             Nd (int): The dimension of the matrix, i.e., the highest dimension in the
                 grid bucket.
-                
+
         """
         # List the fracture points
         self.frac_pts = np.array([[0.2, 0.8], [0.5, 0.5]])
@@ -89,6 +89,12 @@ def domain_boundary_sides(self, g):
     def bc_type(self, g):
         all_bf, *_ = self.domain_boundary_sides(g)
         bc = pp.BoundaryConditionVectorial(g, all_bf, "dir")
+        # Default internal BC is Neumann. We change to Dirichlet for the contact
+        # problem. I.e., the mortar variable represents the displacement on the
+        # fracture faces.
+        frac_face = g.tags["fracture_faces"]
+        bc.is_neu[:, frac_face] = False
+        bc.is_dir[:, frac_face] = True
         return bc
 
     def bc_values(self, g):
@@ -116,12 +122,7 @@ def set_parameters(self):
 
                 # Define boundary condition
                 bc = self.bc_type(g)
-                # Default internal BC is Neumann. We change to Dirichlet for the contact
-                # problem. I.e., the mortar variable represents the displacement on the
-                # fracture faces.
-                frac_face = g.tags["fracture_faces"]
-                bc.is_neu[:, frac_face] = False
-                bc.is_dir[:, frac_face] = True
+
                 # BC and source values
                 bc_val = self.bc_values(g)
                 source_val = self.source(g)

src/porepy/numerics/fv/biot.py

@@ -315,11 +315,9 @@ def _discretize_mech(self, g, data):
                 term.
        """
         parameters_m = data[pp.PARAMETERS][self.mechanics_keyword]
-        parameters_f = data[pp.PARAMETERS][self.flow_keyword]
         matrices_m = data[pp.DISCRETIZATION_MATRICES][self.mechanics_keyword]
         matrices_f = data[pp.DISCRETIZATION_MATRICES][self.flow_keyword]
         bound_mech = parameters_m["bc"]
-        bound_flow = parameters_f["bc"]
         constit = parameters_m["fourth_order_tensor"]
 
         eta = parameters_m.get("mpsa_eta", fvutils.determine_eta(g))
@@ -354,7 +352,7 @@ def _discretize_mech(self, g, data):
         if bound_mech.num_faces == subcell_topology.num_subfno_unique:
             subface_rhs = True
         else:
-            #  If they har given on the faces, expand the boundary conditions them
+            # If they are given on the faces, expand the boundary conditions
             bound_mech = pp.fvutils.boundary_to_sub_boundary(
                 bound_mech, subcell_topology
             )
@@ -1257,6 +1255,7 @@ def assemble_int_bound_displacement_source(
         normal_component = rotation.project_normal(g.num_cells)
 
         biot_alpha = data[pp.PARAMETERS][self.flow_keyword]["biot_alpha"]
+#        aperture = data[pp.PARAMETERS][self.flow_keyword]["aperture"]
         if biot_alpha != 1:
             warnings.warn(
                 "Are you sure you want a non-unitary biot alpha for the fracture?"

test/integration/test_contact_mechanics.py

@@ -153,7 +153,14 @@ def bc_values(self, g):
 
     def bc_type(self, g):
         _, _, _, north, south, _, _ = self.domain_boundary_sides(g)
-        return pp.BoundaryConditionVectorial(g, north + south, "dir")
+        bc = pp.BoundaryConditionVectorial(g, north + south, "dir")
+        # Default internal BC is Neumann. We change to Dirichlet for the contact
+        # problem. I.e., the mortar variable represents the displacement on the
+        # fracture faces.
+        frac_face = g.tags["fracture_faces"]
+        bc.is_neu[:, frac_face] = False
+        bc.is_dir[:, frac_face] = True
+        return bc
 
 
 if __name__ == "__main__":

test/integration/test_contact_mechanics_biot.py

@@ -237,7 +237,14 @@ def source_scalar(self, g):
 
     def bc_type_mechanics(self, g):
         _, _, _, north, south, _, _ = self.domain_boundary_sides(g)
-        return pp.BoundaryConditionVectorial(g, north + south, "dir")
+        bc = pp.BoundaryConditionVectorial(g, north + south, "dir")
+        # Default internal BC is Neumann. We change to Dirichlet for the contact
+        # problem. I.e., the mortar variable represents the displacement on the
+        # fracture faces.
+        frac_face = g.tags["fracture_faces"]
+        bc.is_neu[:, frac_face] = False
+        bc.is_dir[:, frac_face] = True
+        return bc
 
     def bc_type_scalar(self, g):
         _, _, _, north, south, _, _ = self.domain_boundary_sides(g)