diff --git a/src/porepy/models/contact_mechanics_biot_model.py b/src/porepy/models/contact_mechanics_biot_model.py index 813453feb2..f5318809e0 100644 --- a/src/porepy/models/contact_mechanics_biot_model.py +++ b/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,13 +444,13 @@ 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) @@ -459,4 +458,5 @@ def run_biot(setup, atol=1e-10): assembler.distribute_variable(sol) setup.export_step() errors.append(newton_errors) + setup.newton_errors = errors setup.export_pvd() diff --git a/src/porepy/models/contact_mechanics_model.py b/src/porepy/models/contact_mechanics_model.py index 256dc84c00..508f8ace5c 100644 --- a/src/porepy/models/contact_mechanics_model.py +++ b/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