diff --git a/src/porepy/__init__.py b/src/porepy/__init__.py index 8281316d2d..1cc6dea781 100644 --- a/src/porepy/__init__.py +++ b/src/porepy/__init__.py @@ -60,9 +60,6 @@ import porepy.numerics -# Models -import porepy.models - # Transport related from porepy.numerics.fv.upwind import Upwind from porepy.numerics.interface_laws.hyperbolic_interface_laws import UpwindCoupling diff --git a/src/porepy/models/__init__.py b/src/porepy/models/__init__.py index 549f16303d..e69de29bb2 100644 --- a/src/porepy/models/__init__.py +++ b/src/porepy/models/__init__.py @@ -1 +0,0 @@ -from . import contact_mechanics # , contact_mechanics_biot diff --git a/src/porepy/models/contact_mechanics_biot_model.py b/src/porepy/models/contact_mechanics_biot_model.py new file mode 100644 index 0000000000..64b3cefe0e --- /dev/null +++ b/src/porepy/models/contact_mechanics_biot_model.py @@ -0,0 +1,303 @@ +""" +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 +""" +import numpy as np +import scipy.sparse as sps +from scipy.spatial.distance import cdist +import porepy as pp + +from porepy.utils import assign_discretizations +import porepy.models.contact_mechanics_model as contact_model + + +class ContactMechanicsBiot(contact_model.ContactMechanics): + def __init__(self, mesh_args, folder_name): + super().__init__(mesh_args, folder_name) + + # Temperature + self.scalar_variable = "T" + self.mortar_scalar_variable = "mortar_" + self.scalar_variable + self.scalar_coupling_term = "robin_" + self.scalar_variable + self.scalar_parameter_key = "temperature" + + # Scaling coefficients + self.scalar_scale = 1 + self.length_scale = 1 + + # Time + self.time_step = 1e0 * self.length_scale ** 2 + self.end_time = self.time_step * 1 + + self.T_0 = 0 + self.s_0 = 1 + + def bc_type(self, g, key, t=0): + if key == self.mechanics_parameter_key: + # Use parent class method for mechanics + bc = super().bc_type(g) + elif key == self.scalar_parameter_key: + # Define boundary regions + all_bf, *_ = self.domain_boundary_sides(g) + # Define boundary condition on faces + bc = pp.BoundaryCondition(g, all_bf, "dir") + else: + raise ValueError("No BCs implemented for keyword " + str(key)) + return bc + + def bc_values(self, g, key, t=0): + # Set the boundary values + if key == self.mechanics_parameter_key: + bc_values = super().bc_values(g) + elif key == self.scalar_parameter_key: + bc_values = np.zeros(g.num_faces) + else: + raise ValueError("No BC values implemented for keyword " + str(key)) + return bc_values + + def source(self, g, key, t=0): + if key == self.mechanics_parameter_key: + values = super().source(g) + elif key == self.scalar_parameter_key: + values = np.zeros(g.num_cells) + else: + raise ValueError("No source values implemented for keyword " + str(key)) + return values + + def biot_alpha(self): + return 1 + + def compute_aperture(self, g): + apertures = np.ones(g.num_cells) + if g.dim < self.Nd: + apertures *= 0.1 + return apertures + + def set_parameters(self): + """ + Set the parameters for the simulation. The stress is given in GPa. + """ + + self.set_mechanics_parameters() + self.set_scalar_parameters() + + def set_mechanics_parameters(self): + """ + Set the parameters for the simulation. + """ + gb = self.gb + + for g, d in gb: + if g.dim == self.Nd: + # Rock parameters + lam = np.ones(g.num_cells) + mu = np.ones(g.num_cells) + C = pp.FourthOrderTensor(g.dim, mu, lam) + + # Define boundary condition + bc = self.bc_type(g, self.mechanics_parameter_key) + # 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, self.mechanics_parameter_key) + source_val = self.source(g, self.mechanics_parameter_key) + + pp.initialize_data( + g, + d, + self.mechanics_parameter_key, + { + "bc": bc, + "bc_values": bc_val, + "source": source_val, + "fourth_order_tensor": C, + "time_step": self.time_step, + "biot_alpha": self.biot_alpha(), + }, + ) + + elif g.dim == self.Nd - 1: + friction = self._set_friction_coefficient(g) + pp.initialize_data( + g, + d, + self.mechanics_parameter_key, + {"friction_coefficient": friction, "time_step": self.time_step}, + ) + # Should we keep this, @EK? + for e, d in gb.edges(): + mg = d["mortar_grid"] + + # Parameters for the surface diffusion. + mu = 1 + lmbda = 1 + + pp.initialize_data( + mg, d, self.mechanics_parameter_key, {"mu": mu, "lambda": lmbda} + ) + + def set_scalar_parameters(self): + gb = self.gb + self.Nd = gb.dim_max() + + tensor_scale = self.scalar_scale / self.length_scale ** 2 + kappa = 1 * tensor_scale + mass_weight = 1 + alpha = self.biot_alpha() + for g, d in gb: + bc = self.bc_type(g, self.scalar_parameter_key) + bc_values = self.bc_values(g, self.scalar_parameter_key) + source_values = self.source(g, self.scalar_parameter_key, 0) + + a = self.compute_aperture(g) + cross_sectional_area = np.power(a, self.gb.dim_max() - g.dim) * np.ones( + g.num_cells + ) + diffusivity = pp.SecondOrderTensor(self.Nd, kappa * np.ones(g.num_cells)) + + pp.initialize_data( + g, + d, + self.scalar_parameter_key, + { + "bc": bc, + "bc_values": bc_values, + "mass_weight": mass_weight, + "biot_alpha": alpha, + "source": source_values, + "second_order_tensor": diffusivity, + "aperture": cross_sectional_area, + "time_step": self.time_step, + }, + ) + + # 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) + 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, + self.scalar_parameter_key, + {"normal_diffusivity": normal_diffusivity}, + ) + + def assign_discretisations_and_variables(self): + assign_discretizations.contact_mechanics_and_biot_discretizations(self) + + def discretize_biot(self, gb): + """ + To save computational time, the full Biot equation (without contact mechanics) + is discretized once. This computing the same terms multiple times. + """ + g = gb.grids_of_dimension(gb.dim_max())[0] + d = gb.node_props(g) + biot = pp.Biot( + mechanics_keyword=self.mechanics_parameter_key, + flow_keyword=self.scalar_parameter_key, + vector_variable=self.displacement_variable, + scalar_variable=self.scalar_variable, + ) + biot.discretize(g, d) + + def initial_condition(self): + """ + Initial guess for Newton iteration, scalar variable and bc_values (for time + discretization). + """ + super().initial_condition() + + for g, d in self.gb: + # Initial value for the scalar variable. + initial_scalar_value = self.T_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"] + mech_dict = {"bc_values": bc_values} + d[pp.STATE].update({self.mechanics_parameter_key: mech_dict}) + + +def run_biot(setup, atol=1e-10): + """ + Function for solving the time dependent Biot equations with a non-linear Coulomb + contact condition on the fractures. + + The parameter keyword from the elasticity is assumed the same as the + parameter keyword from the contact condition. + + In addition to the standard parameters for Biot we also require the following + under the mechanics keyword (returned from setup.set_parameters): + 'friction_coeff' : The coefficient of friction + 'c' : The numerical parameter in the non-linear complementary function. + + Arguments: + setup: A setup class with methods: + set_parameters(): assigns data to grid bucket. + assign_discretizations_and_variables(): assign the appropriate + discretizations and variables to each node and edge of the grid + bucket. + create_grid(): Create grid bucket and set rotations for all fractures. + initial_condition(): Set initial guesses for the iterates (contact + traction and mortar displacement) and the scalar variable. + and attributes: + end_time: End time time of simulation. + time_step: Time step size + """ + if "gb" not in setup.__dict__: + setup.create_grid() + gb = setup.gb + # Extract the grids we use + ambient_dim = gb.dim_max() + g_max = gb.grids_of_dimension(ambient_dim)[0] + d_max = gb.node_props(g_max) + + # set parameters + setup.set_parameters() + setup.initial_condition() + + setup.assign_discretisations_and_variables() + # Define rotations + # Set up assembler and get initial condition + assembler = pp.Assembler(gb) + + u = d_max[pp.STATE][setup.displacement_variable] + + # Discretize with the biot class + setup.discretize_biot(gb) + errors = [] + + t = 0.0 + dt = setup.time_step + T = setup.end_time + k = 0 + times = [t] + while t < T: + t += dt + k += 1 + print("Time step: ", k, "/", int(np.ceil(T / dt))) + + times.append(t) + # Prepare for Newton + counter_newton = 0 + converged_newton = False + max_newton = 10 + 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 + ) + counter_newton += 1 + newton_errors.append(error) + # Prepare for next time step + assembler.distribute_variable(sol) + errors.append(newton_errors) diff --git a/src/porepy/models/contact_mechanics.py b/src/porepy/models/contact_mechanics_model.py similarity index 91% rename from src/porepy/models/contact_mechanics.py rename to src/porepy/models/contact_mechanics_model.py index 97432d780c..ca6b9487b2 100644 --- a/src/porepy/models/contact_mechanics.py +++ b/src/porepy/models/contact_mechanics_model.py @@ -240,9 +240,10 @@ def initial_condition(self): if mg.dim == self.Nd - 1: size = mg.num_cells * self.Nd state = { + self.mortar_displacement_variable: np.zeros(mg.num_cells * self.Nd), "previous_iterate": { self.mortar_displacement_variable: np.zeros(size) - } + }, } pp.set_state(d, state) @@ -391,44 +392,52 @@ def run_mechanics(setup): while counter_newton <= max_newton and not converged_newton: print("Newton iteration number: ", counter_newton, "/", max_newton) + sol, u0, error, converged_newton = newton_iteration(assembler, setup, u0) counter_newton += 1 - # Re-discretize the nonlinear term - assembler.discretize(term_filter=setup.friction_coupling_term) - - # Assemble and solve - A, b = assembler.assemble_matrix_rhs() - # if gb.num_cells() > 6e4: #4 - # sol = solvers.amg(gb, A, b) - # else: - sol = sps.linalg.spsolve(A, b) + viz.write_vtk({"ux": u0[::2], "uy": u0[1::2]}) + errors.append(error) - # Obtain the current iterate for the displacement, and distribute the current - # iterates for mortar displacements and contact traction. - u1 = setup.extract_iterate(assembler, sol) + if counter_newton > max_newton and not converged_newton: + raise ValueError("Newton iterations did not converge") + assembler.distribute_variable(sol) - viz.write_vtk({"ux": u1[::2], "uy": u1[1::2]}) - # Calculate the error - solution_norm = l2_norm_cell(g_max, u1) - iterate_difference = l2_norm_cell(g_max, u1, u0) +def newton_iteration(assembler, setup, u0, 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] - # The if is intended to avoid division through zero - if solution_norm < 1e-12 and iterate_difference < 1e-12: - converged_newton = True - error = np.sum((u1 - u0) ** 2) - else: - if iterate_difference / solution_norm < 1e-10: - converged_newton = True - error = np.sum((u1 - u0) ** 2) / np.sum(u1 ** 2) + # Re-discretize the nonlinear term + assembler.discretize(term_filter=setup.friction_coupling_term) - print("Error: ", error) - errors.append(error) - # Prepare for next iteration - u0 = u1 + # Assemble and solve + A, b = assembler.assemble_matrix_rhs() - if counter_newton > max_newton and not converged_newton: - raise ValueError("Newton iterations did not converge") - assembler.distribute_variable(sol) + if solver is None: + sol = sps.linalg.spsolve(A, b) + # else: + # # sol = solvers.amg(gb, A, b) + + # Obtain the current iterate for the displacement, and distribute the current + # iterates for mortar displacements and contact traction. + u1 = setup.extract_iterate(assembler, sol) + + # Calculate the error + solution_norm = l2_norm_cell(g_max, u1) + 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: + converged = True + error = np.sum((u1 - u0) ** 2) + else: + if iterate_difference / solution_norm < 1e-10: + converged = True + error = np.sum((u1 - u0) ** 2) / np.sum(u1 ** 2) + + print("Error: ", error) + + return sol, u1, error, converged def l2_norm_cell(g, u, uref=None): diff --git a/src/porepy/numerics/contact_mechanics/contact_conditions.py b/src/porepy/numerics/contact_mechanics/contact_conditions.py index d0c647e343..85d1b08e88 100644 --- a/src/porepy/numerics/contact_mechanics/contact_conditions.py +++ b/src/porepy/numerics/contact_mechanics/contact_conditions.py @@ -22,11 +22,8 @@ def __init__(self, keyword, ambient_dimension): self.dim = ambient_dimension - self.surface_variable = "mortar_u" - self.contact_variable = "contact_force" - - self.friction_parameter_key = "friction" - self.surface_parameter_key = "surface" + self.mortar_displacement_variable = "mortar_u" + self.contact_variable = "contact_traction" self.traction_discretization = "traction_discretization" self.displacement_discretization = "displacement_discretization" @@ -80,9 +77,7 @@ def discretize(self, g_h, g_l, data_h, data_l, data_edge): # Process input parameters_l = data_l[pp.PARAMETERS] - friction_coefficient = parameters_l[self.friction_parameter_key][ - "friction_coefficient" - ] + friction_coefficient = parameters_l[self.keyword]["friction_coefficient"] if np.asarray(friction_coefficient).size == 1: friction_coefficient = friction_coefficient * np.ones(g_l.num_cells) @@ -124,7 +119,7 @@ def discretize(self, g_h, g_l, data_h, data_l, data_edge): displacement_jump_global_coord = ( mg.mortar_to_slave_avg(nd=self.dim) * mg.sign_of_mortar_sides(nd=self.dim) - * data_edge[pp.STATE]["previous_iterate"][self.surface_variable] + * data_edge[pp.STATE]["previous_iterate"][self.mortar_displacement_variable] ) # Rotated displacement jumps. These are in the local coordinates, on # the lower-dimensional grid 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 169166bf28..c0ee5a0345 100644 --- a/src/porepy/numerics/interface_laws/contact_mechanics_interface_laws.py +++ b/src/porepy/numerics/interface_laws/contact_mechanics_interface_laws.py @@ -688,9 +688,8 @@ class DivUCoupling: """ def __init__(self, variable, discr_master, discr_slave): - # self.mechanics_keyword = discr_master.mechanics_keyword - # Set variable names for the vector variable (displacement), used to access - # solutions from previous time steps + # Set variable names for the vector variable on the nodes (displacement), used + # to access solutions from previous time steps. self.variable = variable # The terms are added by calls to assemble methods of DivU discretizations, # namely assemble_int_bound_displacement_trace for the master and diff --git a/src/porepy/utils/assign_discretizations.py b/src/porepy/utils/assign_discretizations.py index f59ca587b4..e45bfb0cfa 100644 --- a/src/porepy/utils/assign_discretizations.py +++ b/src/porepy/utils/assign_discretizations.py @@ -65,7 +65,7 @@ def contact_mechanics_discretizations(setup): d[pp.PRIMARY_VARIABLES] = {} -def contact_mechanics_and_biot_discretizations(setup, subtract_fracture_pressure): +def contact_mechanics_and_biot_discretizations(setup, subtract_fracture_pressure=True): """ Assign the discretizations for fracture deformation with a coupled scalar (pressure) in both dimensions. No fracture intersections are allowed (for now). @@ -73,12 +73,11 @@ def contact_mechanics_and_biot_discretizations(setup, subtract_fracture_pressure Setup should have a gb field, and the following names specified: Parameter keys: mechanics_parameter_key - friction_parameter_key scalar_parameter_key Variables: displacement_variable - higher-dimensional displacements mortar_displacement_variable - displacement on the internal boundary - contact_variable - represents traction on the fracture + contact_traction_variable - represents traction on the fracture scalar_variable - scalar (pressure) in both dimensions mortar_scalar_variable - darcy flux subtract_fracture_pressure (bool): Whether or not to subtract the fracture pressure @@ -94,12 +93,10 @@ def contact_mechanics_and_biot_discretizations(setup, subtract_fracture_pressure # Define discretization # For the Nd domain we solve linear elasticity with mpsa. mpsa = pp.Mpsa(key_m) - empty_discr = pp.VoidDiscretization( - setup.friction_parameter_key, ndof_cell=ambient_dim - ) + empty_discr = pp.VoidDiscretization(key_m, ndof_cell=ambient_dim) # Scalar discretizations (all dimensions) diff_disc_s = IE_discretizations.ImplicitMpfa(key_s) - mass_disc_s = IE_discretizations.ImplicitMassMatrix(key_s, variable=var_s) + mass_disc_s = IE_discretizations.ImplicitMassMatrix(key_s, var_s) source_disc_s = pp.ScalarSource(key_s) # Coupling discretizations # All dimensions @@ -108,7 +105,7 @@ def contact_mechanics_and_biot_discretizations(setup, subtract_fracture_pressure ) # Nd grad_p_disc = pp.GradP(key_m) - stabilization_disc_s = pp.BiotStabilization(key_s, variable=var_s) + stabilization_disc_s = pp.BiotStabilization(key_s, var_s) # Assign node discretizations for g, d in gb: @@ -130,11 +127,11 @@ def contact_mechanics_and_biot_discretizations(setup, subtract_fracture_pressure } elif g.dim == ambient_dim - 1: d[pp.PRIMARY_VARIABLES] = { - setup.contact_variable: {"cells": ambient_dim}, + setup.contact_traction_variable: {"cells": ambient_dim}, var_s: {"cells": 1}, } d[pp.DISCRETIZATION] = { - setup.contact_variable: {"empty": empty_discr}, + setup.contact_traction_variable: {"empty": empty_discr}, var_s: { "diffusion": diff_disc_s, "mass": mass_disc_s, @@ -145,9 +142,9 @@ def contact_mechanics_and_biot_discretizations(setup, subtract_fracture_pressure d[pp.PRIMARY_VARIABLES] = {} # Define edge discretizations for the mortar grid - contact_law = pp.ColoumbContact(setup.friction_parameter_key, ambient_dim) + contact_law = pp.ColoumbContact(setup.mechanics_parameter_key, ambient_dim) contact_discr = pp.PrimalContactCoupling( - setup.friction_parameter_key, mpsa, contact_law + setup.mechanics_parameter_key, mpsa, contact_law ) # Account for the mortar displacements effect on scalar balance in the # matrix, as an internal boundary contribution, @@ -178,7 +175,7 @@ def contact_mechanics_and_biot_discretizations(setup, subtract_fracture_pressure d[pp.COUPLING_DISCRETIZATION] = { setup.friction_coupling_term: { g_h: (var_d, "mpsa"), - g_l: (setup.contact_variable, "empty"), + g_l: (setup.contact_traction_variable, "empty"), (g_h, g_l): (setup.mortar_displacement_variable, contact_discr), }, setup.scalar_coupling_term: { diff --git a/test/integration/test_contact_coupling_biot.py b/test/integration/test_contact_coupling_biot.py deleted file mode 100644 index bdce512c94..0000000000 --- a/test/integration/test_contact_coupling_biot.py +++ /dev/null @@ -1,603 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Integration tests for contact mechanics with pressure coupling. - -We have the full Biot equations in the matrix, and mass conservation and contact -conditions in the non-intersecting fracture(s). For the contact mechanical part of this -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 -from test_contact_mechanics import SetupContactMechanics -from porepy.utils.derived_discretizations import implicit_euler as discretizations - - -class TestContactMechanicsBiot(unittest.TestCase): - def _solve(self, model): - _ = solve_biot(model) - gb = model.gb - - nd = gb.dim_max() - - g2 = gb.grids_of_dimension(2)[0] - g1 = gb.grids_of_dimension(1)[0] - - d_m = gb.edge_props((g1, g2)) - d_1 = gb.node_props(g1) - - mg = d_m["mortar_grid"] - - u_mortar = d_m[pp.STATE][model.surface_variable] - contact_force = d_1[pp.STATE][model.contact_variable] - fracture_pressure = d_1[pp.STATE][model.scalar_variable] - - displacement_jump_global_coord = ( - mg.mortar_to_slave_avg(nd=nd) * mg.sign_of_mortar_sides(nd=nd) * u_mortar - ) - projection = d_m["tangential_normal_projection"] - - project_to_local = projection.project_tangential_normal(int(mg.num_cells / 2)) - u_mortar_local = project_to_local * displacement_jump_global_coord - u_mortar_local_decomposed = u_mortar_local.reshape((2, -1), order="F") - - contact_force = contact_force.reshape((2, -1), order="F") - - return u_mortar_local_decomposed, contact_force, fracture_pressure - - def test_pull_top_positive_opening(self): - - model = SetupContactMechanicsBiot( - ux_bottom=0, uy_bottom=0, ux_top=0, uy_top=0.001 - ) - - u_mortar, contact_force, fracture_pressure = self._solve(model) - - # All components should be open in the normal direction - self.assertTrue(np.all(u_mortar[1] < 0)) - - # By symmetry (reasonable to expect from this grid), the jump in tangential - # deformation should be zero. - self.assertTrue(np.abs(np.sum(u_mortar[0])) < 1e-5) - - # The contact force in normal direction should be zero - - # NB: This assumes the contact force is expressed in local coordinates - self.assertTrue(np.all(np.abs(contact_force) < 1e-7)) - - # Check that the dilation of the fracture yields a negative fracture pressure - self.assertTrue(np.all(fracture_pressure < -1e-7)) - - def test_pull_bottom_positive_opening(self): - - model = SetupContactMechanicsBiot( - ux_bottom=0, uy_bottom=-0.001, ux_top=0, uy_top=0 - ) - - u_mortar, contact_force, fracture_pressure = self._solve(model) - - # All components should be open in the normal direction - self.assertTrue(np.all(u_mortar[1] < 0)) - - # By symmetry (reasonable to expect from this grid), the jump in tangential - # deformation should be zero. - self.assertTrue(np.abs(np.sum(u_mortar[0])) < 1e-5) - - # The contact force in normal direction should be zero - - # NB: This assumes the contact force is expressed in local coordinates - self.assertTrue(np.all(np.abs(contact_force) < 1e-7)) - - # Check that the dilation of the fracture yields a negative fracture pressure - self.assertTrue(np.all(fracture_pressure < -1e-7)) - - def test_push_top_zero_opening(self): - - model = SetupContactMechanicsBiot( - ux_bottom=0, uy_bottom=0, ux_top=0, uy_top=-0.001 - ) - - u_mortar, contact_force, fracture_pressure = self._solve(model) - - # All components should be closed in the normal direction - self.assertTrue(np.abs(np.sum(u_mortar[1])) < 1e-5) - - # Contact force in normal direction should be negative - self.assertTrue(np.all(contact_force[1] < 0)) - - # Compression of the domain yields a (slightly) positive fracture pressure - self.assertTrue(np.all(fracture_pressure > 1e-10)) - - def test_push_bottom_zero_opening(self): - - model = SetupContactMechanicsBiot( - ux_bottom=0, uy_bottom=0.001, ux_top=0, uy_top=0 - ) - - u_mortar, contact_force, fracture_pressure = self._solve(model) - - # All components should be closed in the normal direction - self.assertTrue(np.abs(np.sum(u_mortar[1])) < 1e-5) - - # Contact force in normal direction should be negative - self.assertTrue(np.all(contact_force[1] < 0)) - - # Compression of the domain yields a (slightly) positive fracture pressure - self.assertTrue(np.all(fracture_pressure > 1e-10)) - - def test_positive_fracture_pressure_positive_opening(self): - - model = SetupContactMechanicsBiot( - ux_bottom=0, uy_bottom=0, ux_top=0, uy_top=0, source_value=0.001 - ) - - u_mortar, contact_force, fracture_pressure = self._solve(model) - - # All components should be open in the normal direction - self.assertTrue(np.all(u_mortar[1] < 0)) - - # By symmetry (reasonable to expect from this grid), the jump in tangential - # deformation should be zero. - self.assertTrue(np.abs(np.sum(u_mortar[0])) < 1e-5) - - # The contact force in normal direction should be zero - - # NB: This assumes the contact force is expressed in local coordinates - self.assertTrue(np.all(np.abs(contact_force) < 1e-7)) - - # Fracture pressure is positive - self.assertTrue(np.all(fracture_pressure > 1e-7)) - - -def solve_biot(setup): - """ - Function for solving the time dependent Biot equations with a non-linear Coulomb - contact condition on the fractures. - - See solve_contact_mechanics in test_contact_mechanics for assumptions on the - mechanical parameters set. In addition, we require parameters and discretizations - for the pressure and coupling terms, see SetupContactMechanicsBiot. - - Arguments: - setup: A setup class with methods: - set_parameters(): assigns data for the contact mechanics problem. See - test_contact_mechanics. - set_scalar_parameters(): assign data for the scalar parameter, here - pressure - create_grid(): Create and return the grid bucket - initial_condition(): Returns initial conditions. - and attributes: - end_time: End time time of simulation. - """ - gb = setup.create_grid() - # Extract the grids we use - dim = gb.dim_max() - g_max = gb.grids_of_dimension(dim)[0] - d_max = gb.node_props(g_max) - - # set parameters - setup.set_parameters(gb) - setup.set_scalar_parameters() - setup.initial_condition() - - # Shorthand for some parameters - dt = d_max[pp.PARAMETERS][setup.scalar_parameter_key]["time_step"] - setup.assign_discretisations() - # Define rotations - pp.contact_conditions.set_projections(gb) - # Set up assembler and get initial condition - assembler = pp.Assembler(gb) - - u_k_minus_one = d_max[pp.STATE][setup.displacement_variable].reshape( - (dim, -1), order="F" - ) - - # Discretize with the Biot class - setup.discretize_biot(gb) - - def l2_error_cell(g, u, uref=None): - if uref is None: - norm = np.reshape(u ** 2, (g.dim, g.num_cells), order="F") * g.cell_volumes - else: - norm = ( - np.reshape((u - uref) ** 2, (g.dim, g.num_cells), order="F") - * g.cell_volumes - ) - return np.sum(norm) - - t = 0.0 - T = setup.end_time - k = 0 - times = [] - newton_it = 0 - - while t < T: - t += dt - k += 1 - times.append(t) - # Prepare for Newton iteration - counter_newton = 0 - converged_newton = False - max_newton = 12 - while counter_newton <= max_newton and not converged_newton: - counter_newton += 1 - # Rediscretize the contact conditions (remaining discretizations assumed - # constant in time). - assembler.discretize(term_filter=setup.friction_coupling_term) - - # Reassemble and solve - A, b = assembler.assemble_matrix_rhs() - sol = sps.linalg.spsolve(A, b) - - # Update the previous iterate of the mortar displacement and contact - # traction, and obtain current matrix displacement iterate. - u_k = distribute_iterate( - assembler, setup, sol, setup.surface_variable, setup.contact_variable - ) - u_k = u_k.reshape((dim, -1), order="F") - # Calculate the error - solution_norm = l2_error_cell(g_max, u_k) - iterate_difference = l2_error_cell(g_max, u_k, u_k_minus_one) - if iterate_difference / solution_norm < 1e-10: - converged_newton = True - - # Prepare for next Newton iteration - u_k_minus_one = u_k - newton_it += 1 - assembler.distribute_variable(sol) - - return sol - - -def distribute_iterate( - assembler, setup, values, mortar_displacement_variable, contact_traction_variable -): - """ Update the previous iterate of the mortar displacement and contact traction, - and obtain current matrix displacement iterate. - - Method is a tailored copy from assembler.distribute_variable. - """ - dof = np.cumsum(np.append(0, np.asarray(assembler.full_dof))) - var_name = setup.displacement_variable - - for pair, bi in assembler.block_dof.items(): - g = pair[0] - name = pair[1] - # Avoid edges - if not isinstance(g, pp.Grid): - if name == mortar_displacement_variable: - mortar_u = values[dof[bi] : dof[bi + 1]] - data = setup.gb.edge_props(g) - data[pp.STATE]["previous_iterate"][ - mortar_displacement_variable - ] = mortar_u - continue - # Only interested in highest dimension - if g.dim < setup.gb.dim_max(): - if name == contact_traction_variable: - contact = values[dof[bi] : dof[bi + 1]] - data = setup.gb.node_props(g) - data[pp.STATE]["previous_iterate"][contact_traction_variable] = contact - - continue - # Only need the displacement - if name != var_name: - continue - - u = values[dof[bi] : dof[bi + 1]] - return u - - -class SetupContactMechanicsBiot(SetupContactMechanics): - def __init__(self, ux_bottom, uy_bottom, ux_top, uy_top, source_value=0): - super().__init__(ux_bottom, uy_bottom, ux_top, uy_top) - self.scalar_variable = "p" - self.scalar_coupling_term = "robin_" + self.scalar_variable - self.scalar_parameter_key = "flow" - self.pressure_source_value = source_value - - self.scalar_scale = 1 - self.length_scale = 1 - self.time_step = 1 - self.end_time = 1 - - def biot_alpha(self): - return 1 - - def sources(self, g, key, t=0): - if key == self.mechanics_parameter_key: - values = np.zeros((g.dim, g.num_cells)) - values = values.ravel("F") - elif key == self.scalar_parameter_key: - if g.dim == 2: - values = np.zeros(g.num_cells) - else: - values = ( - self.pressure_source_value * self.time_step * np.ones(g.num_cells) - ) - else: - raise ValueError("No BC values implemented for keyword " + str(key)) - return values - - def set_scalar_parameters(self): - gb = self.gb - ambient_dim = gb.dim_max() - tensor_scale = self.scalar_scale / self.length_scale ** 2 - k_frac = 100 - a = 1e-3 - for g, d in gb: - # Define boundary regions - top = g.face_centers[ambient_dim - 1] > self.box["ymax"] - 1e-9 - bot = g.face_centers[ambient_dim - 1] < self.box["ymin"] + 1e-9 - - # Define boundary condition - bc = pp.BoundaryCondition(g, top + bot, "dir") - bc_values = np.zeros(g.num_faces) - - alpha = self.biot_alpha() - - if g.dim == ambient_dim: - kxx = 1 * tensor_scale * np.ones(g.num_cells) - K = pp.SecondOrderTensor(ambient_dim, kxx) - - mass_weight = 1e-1 - - pp.initialize_data( - g, - d, - self.scalar_parameter_key, - { - "bc": bc, - "bc_values": bc_values, - "mass_weight": mass_weight, - "aperture": np.ones(g.num_cells), - "biot_alpha": alpha, - "time_step": self.time_step, - "source": self.sources(g, self.scalar_parameter_key, 0), - "second_order_tensor": K, - }, - ) - # Add Biot alpha and time step to the mechanical parameters - d[pp.PARAMETERS].update_dictionaries( - self.mechanics_parameter_key, - {"biot_alpha": alpha, "time_step": self.time_step}, - ) - - elif g.dim == ambient_dim - 1: - kxx = k_frac * tensor_scale * np.ones(g.num_cells) - K = pp.SecondOrderTensor(ambient_dim, kxx) - mass_weight = 1e-1 # compressibility - cross_sectional_area = a * np.ones(g.num_cells) - pp.initialize_data( - g, - d, - self.scalar_parameter_key, - { - "bc": bc, - "bc_values": bc_values, - "mass_weight": mass_weight, - "source": self.sources(g, self.scalar_parameter_key, 0), - "second_order_tensor": K, - "aperture": cross_sectional_area, - "time_step": self.time_step, - "biot_alpha": alpha, - }, - ) - pp.initialize_data( - g, - d, - self.mechanics_parameter_key, - {"biot_alpha": 1}, # Enters in the div d term for the fracture - ) - else: - # No intersections yet - raise NotImplementedError - - # 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) - mg = data_edge["mortar_grid"] - k = k_frac * tensor_scale * np.ones(mg.num_cells) - k_n = 2 / a * k - data_edge = pp.initialize_data( - e, data_edge, self.scalar_parameter_key, {"normal_diffusivity": k_n} - ) - - def assign_discretisations(self): - gb = self.gb - ambient_dim = gb.dim_max() - key_s, key_m = self.scalar_parameter_key, self.mechanics_parameter_key - var_s, var_d = self.scalar_variable, self.displacement_variable - # Define discretization - # For the 2D domain we solve linear elasticity with mpsa. - mpsa = pp.Mpsa(key_m) - - empty_discr = pp.VoidDiscretization( - self.friction_parameter_key, ndof_cell=ambient_dim - ) - diff_disc_s = discretizations.ImplicitMpfa(key_s) - mass_disc_s = discretizations.ImplicitMassMatrix(key_s, variable=var_s) - source_disc_s = pp.ScalarSource(key_s) - div_u_disc = pp.DivU( - key_m, variable=var_d, mortar_variable=self.surface_variable - ) - grad_p_disc = pp.GradP(key_m) - div_u_disc_frac = pp.DivU( - key_m, variable=var_d, mortar_variable=self.surface_variable - ) - stabilisiation_disc_s = pp.BiotStabilization(key_s, variable=var_s) - coloumb = pp.ColoumbContact(self.friction_parameter_key, ambient_dim) - - # Define discretization parameters - for g, d in gb: - if g.dim == ambient_dim: - d[pp.PRIMARY_VARIABLES] = { - var_d: {"cells": ambient_dim}, - var_s: {"cells": 1}, - } - d[pp.DISCRETIZATION] = { - var_d: {"mpsa": mpsa}, - var_s: { - "diffusion": diff_disc_s, - "mass": mass_disc_s, - "stabilisation": stabilisiation_disc_s, - "source": source_disc_s, - }, - var_d + "_" + var_s: {"grad_p": grad_p_disc}, - var_s + "_" + var_d: {"div_u": div_u_disc}, - } - elif g.dim == ambient_dim - 1: - d[pp.PRIMARY_VARIABLES] = { - self.contact_variable: {"cells": ambient_dim}, - var_s: {"cells": 1}, - } - d[pp.DISCRETIZATION] = { - self.contact_variable: {"empty": empty_discr}, - var_s: { - "diffusion": diff_disc_s, - "mass": mass_disc_s, - "source": source_disc_s, - }, - } - else: - d[pp.PRIMARY_VARIABLES] = {} - - coloumb = pp.ColoumbContact(self.friction_parameter_key, ambient_dim) - contact = pp.PrimalContactCoupling(self.friction_parameter_key, mpsa, coloumb) - # Add the div_u contribution from the mortar displacements to the mass balance - # equations in matrix and fractures. - div_u_coupling = pp.DivUCoupling( - self.displacement_variable, div_u_disc, div_u_disc_frac - ) - # Add the matrix pressure contribution to the reconstructed force on the - # fracture faces, which enters into the force balance equations on the fractures. - # This discretization needs the keyword used to store the grad p discretization: - matrix_pressure_to_contact = pp.MatrixScalarToForceBalance( - key_m, mass_disc_s, mass_disc_s - ) - # Subtract the fracture pressure, to ensure that the contact conditions are - # formulated on the _contact_ forces only. This is needed because the total - # forces on a fracture surface is the sum of contact and pressure forces. - fracture_pressure_to_contact = pp.FractureScalarToForceBalance( - mass_disc_s, mass_disc_s - ) - for e, d in gb.edges(): - g_l, g_h = gb.nodes_of_edge(e) - - if g_h.dim == ambient_dim: - d[pp.PRIMARY_VARIABLES] = { - self.surface_variable: {"cells": ambient_dim}, - self.scalar_variable: {"cells": 1}, - } - - d[pp.COUPLING_DISCRETIZATION] = { - self.friction_coupling_term: { - g_h: (var_d, "mpsa"), - g_l: (self.contact_variable, "empty"), - (g_h, g_l): (self.surface_variable, contact), - }, - self.scalar_coupling_term: { - g_h: (var_s, "diffusion"), - g_l: (var_s, "diffusion"), - e: (self.scalar_variable, pp.RobinCoupling(key_s, diff_disc_s)), - }, - "div_u_coupling": { - g_h: ( - var_s, - "mass", - ), # This is really the div_u, but this is not implemented - g_l: (var_s, "mass"), - e: (self.surface_variable, div_u_coupling), - }, - "matrix_pressure_to_force_balance": { - g_h: (var_s, "mass"), - g_l: (var_s, "mass"), - e: (self.surface_variable, matrix_pressure_to_contact), - }, - "fracture_pressure_to_force_balance": { - g_h: (var_s, "mass"), - g_l: (var_s, "mass"), - e: (self.surface_variable, fracture_pressure_to_contact), - }, - } - else: - d[pp.PRIMARY_VARIABLES] = {} - - def discretize_biot(self, gb): - """ - Discretization of Biot equations is done once, instead of separate - discretization methods for each of the classes DivU, GradP and BiotStabilization. - """ - g = gb.grids_of_dimension(gb.dim_max())[0] - d = gb.node_props(g) - biot = pp.Biot( - mechanics_keyword=self.mechanics_parameter_key, - flow_keyword=self.scalar_parameter_key, - vector_variable=self.displacement_variable, - scalar_variable=self.scalar_variable, - ) - biot.discretize(g, d) - - def initial_condition(self): - """ - Initial guess for Newton iteration. - """ - gb = self.gb - - ambient_dimension = gb.dim_max() - - for g, d in gb: - nc_nd = g.num_cells * ambient_dimension - # Initial value for the scalar variable - initial_scalar_value = 0 * np.ones(g.num_cells) - if g.dim == ambient_dimension: - # Initialize displacement variable - key_m = self.mechanics_parameter_key - bc_dict = {"bc_values": d[pp.PARAMETERS][key_m]["bc_values"]} - pp.set_state( - d, - { - self.displacement_variable: np.zeros(nc_nd), - self.scalar_variable: initial_scalar_value, - key_m: bc_dict, - }, - ) - elif g.dim == ambient_dimension - 1: - # Initialize contact variable - traction = np.vstack( - (np.zeros(g.num_cells), -100 * np.ones(g.num_cells)) - ).ravel(order="F") - pp.set_state( - d, - { - self.contact_variable: traction, - self.scalar_variable: initial_scalar_value, - "previous_iterate": {self.contact_variable: traction}, - }, - ) - - for e, d in gb.edges(): - mg = d["mortar_grid"] - - if mg.dim == 1: - nc_nd = mg.num_cells * ambient_dimension - pp.set_state( - d, - { - self.surface_variable: np.zeros(nc_nd), - self.scalar_variable: np.zeros(mg.num_cells), - "previous_iterate": {self.surface_variable: np.zeros(nc_nd)}, - }, - ) - - def _set_friction_coefficient(self, g): - friction_coefficient = 0.5 * np.ones(g.num_cells) - return friction_coefficient - - -if __name__ == "__main__": - unittest.main() diff --git a/test/integration/test_contact_mechanics.py b/test/integration/test_contact_mechanics.py index 436fcd2a91..5eef833954 100644 --- a/test/integration/test_contact_mechanics.py +++ b/test/integration/test_contact_mechanics.py @@ -6,12 +6,13 @@ import scipy.sparse.linalg as spla import porepy as pp +import porepy.models.contact_mechanics_model as model class TestContactMechanics(unittest.TestCase): - def _solve(self, model): - pp.models.contact_mechanics.run_mechanics(model) - gb = model.gb + def _solve(self, setup): + model.run_mechanics(setup) + gb = setup.gb nd = gb.dim_max() @@ -23,8 +24,8 @@ def _solve(self, model): mg = d_m["mortar_grid"] - u_mortar = d_m[pp.STATE][model.mortar_displacement_variable] - contact_force = d_1[pp.STATE][model.contact_traction_variable] + u_mortar = d_m[pp.STATE][setup.mortar_displacement_variable] + contact_force = d_1[pp.STATE][setup.contact_traction_variable] displacement_jump_global_coord = ( mg.mortar_to_slave_avg(nd=nd) * mg.sign_of_mortar_sides(nd=nd) * u_mortar @@ -41,9 +42,9 @@ def _solve(self, model): def test_pull_top_positive_opening(self): - model = SetupContactMechanics(ux_bottom=0, uy_bottom=0, ux_top=0, uy_top=0.001) + setup = SetupContactMechanics(ux_south=0, uy_bottom=0, ux_north=0, uy_top=0.001) - u_mortar, contact_force = self._solve(model) + u_mortar, contact_force = self._solve(setup) # All components should be open in the normal direction self.assertTrue(np.all(u_mortar[1] < 0)) @@ -59,9 +60,11 @@ def test_pull_top_positive_opening(self): def test_pull_bottom_positive_opening(self): - model = SetupContactMechanics(ux_bottom=0, uy_bottom=-0.001, ux_top=0, uy_top=0) + setup = SetupContactMechanics( + ux_south=0, uy_bottom=-0.001, ux_north=0, uy_top=0 + ) - u_mortar, contact_force = self._solve(model) + u_mortar, contact_force = self._solve(setup) # All components should be open in the normal direction self.assertTrue(np.all(u_mortar[1] < 0)) @@ -77,9 +80,11 @@ def test_pull_bottom_positive_opening(self): def test_push_top_zero_opening(self): - model = SetupContactMechanics(ux_bottom=0, uy_bottom=0, ux_top=0, uy_top=-0.001) + setup = SetupContactMechanics( + ux_south=0, uy_bottom=0, ux_north=0, uy_top=-0.001 + ) - u_mortar, contact_force = self._solve(model) + u_mortar, contact_force = self._solve(setup) # All components should be closed in the normal direction self.assertTrue(np.abs(np.sum(u_mortar[1])) < 1e-5) @@ -89,9 +94,9 @@ def test_push_top_zero_opening(self): def test_push_bottom_zero_opening(self): - model = SetupContactMechanics(ux_bottom=0, uy_bottom=0.001, ux_top=0, uy_top=0) + setup = SetupContactMechanics(ux_south=0, uy_bottom=0.001, ux_north=0, uy_top=0) - u_mortar, contact_force = self._solve(model) + u_mortar, contact_force = self._solve(setup) # All components should be closed in the normal direction self.assertTrue(np.abs(np.sum(u_mortar[1])) < 1e-5) @@ -100,17 +105,17 @@ def test_push_bottom_zero_opening(self): self.assertTrue(np.all(contact_force[1] < 0)) -class SetupContactMechanics(pp.models.contact_mechanics.ContactMechanics): - def __init__(self, ux_bottom, uy_bottom, ux_top, uy_top): +class SetupContactMechanics(model.ContactMechanics): + def __init__(self, ux_south, uy_bottom, ux_north, uy_top): mesh_args = { "mesh_size_frac": 0.5, "mesh_size_min": 0.023, "mesh_size_bound": 0.5, } super().__init__(mesh_args, folder_name="dummy") - self.ux_bottom = ux_bottom + self.ux_south = ux_south self.uy_bottom = uy_bottom - self.ux_top = ux_top + self.ux_north = ux_north self.uy_top = uy_top def create_grid(self, rotate_fracture=False): @@ -137,64 +142,18 @@ def create_grid(self, rotate_fracture=False): self.gb = gb self.Nd = gb.dim_max() - def set_parameters(self): - """ - Set the parameters for the simulation. The stress is given in GPa. - """ - gb = self.gb - - for g, d in gb: - if g.dim == self.Nd: - # Rock parameters - rock = pp.Granite() - lam = rock.LAMBDA * np.ones(g.num_cells) / pp.GIGA - mu = rock.MU * np.ones(g.num_cells) / pp.GIGA - - k = pp.FourthOrderTensor(g.dim, mu, lam) - - # Define boundary regions - top = g.face_centers[g.dim - 1] > np.max(g.nodes[1]) - 1e-9 - bot = g.face_centers[g.dim - 1] < np.min(g.nodes[1]) + 1e-9 - - # Define boundary condition on sub_faces - bc = pp.BoundaryConditionVectorial(g, top + bot, "dir") - frac_face = g.tags["fracture_faces"] - bc.is_neu[:, frac_face] = False - bc.is_dir[:, frac_face] = True - - # Set the boundary values - u_bc = np.zeros((g.dim, g.num_faces)) - - u_bc[0, bot] = self.ux_bottom - u_bc[1, bot] = self.uy_bottom - u_bc[0, top] = self.ux_top - u_bc[1, top] = self.uy_top - - pp.initialize_data( - g, - d, - self.mechanics_parameter_key, - { - "bc": bc, - "bc_values": u_bc.ravel("F"), - "source": 0, - "fourth_order_tensor": k, - }, - ) - - elif g.dim == 1: - friction = self._set_friction_coefficient(g) - pp.initialize_data( - g, - d, - self.mechanics_parameter_key, - {"friction_coefficient": friction}, - ) - - for e, d in gb.edges(): - mg = d["mortar_grid"] - - pp.initialize_data(mg, d, self.mechanics_parameter_key, {}) + def bc_values(self, g): + _, _, _, north, south, _, _ = self.domain_boundary_sides(g) + values = np.zeros((g.dim, g.num_faces)) + values[0, south] = self.ux_south + values[1, south] = self.uy_bottom + values[0, north] = self.ux_north + values[1, north] = self.uy_top + return values.ravel("F") + + def bc_type(self, g): + _, _, _, north, south, _, _ = self.domain_boundary_sides(g) + return pp.BoundaryConditionVectorial(g, north + south, "dir") if __name__ == "__main__": diff --git a/test/integration/test_contact_mechanics_biot.py b/test/integration/test_contact_mechanics_biot.py new file mode 100644 index 0000000000..e787ae2ad9 --- /dev/null +++ b/test/integration/test_contact_mechanics_biot.py @@ -0,0 +1,273 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Integration tests for contact mechanics with pressure coupling. + +We have the full Biot equations in the matrix, and mass conservation and contact +conditions in the non-intersecting fracture(s). For the contact mechanical part of this +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 + + +class TestContactMechanicsBiot(unittest.TestCase): + def _solve(self, setup): + model.run_biot(setup) + gb = setup.gb + + nd = gb.dim_max() + + g2 = gb.grids_of_dimension(2)[0] + g1 = gb.grids_of_dimension(1)[0] + + d_m = gb.edge_props((g1, g2)) + d_1 = gb.node_props(g1) + + mg = d_m["mortar_grid"] + + u_mortar = d_m[pp.STATE][setup.mortar_displacement_variable] + contact_force = d_1[pp.STATE][setup.contact_traction_variable] + fracture_pressure = d_1[pp.STATE][setup.scalar_variable] + + displacement_jump_global_coord = ( + mg.mortar_to_slave_avg(nd=nd) * mg.sign_of_mortar_sides(nd=nd) * u_mortar + ) + projection = d_m["tangential_normal_projection"] + + project_to_local = projection.project_tangential_normal(int(mg.num_cells / 2)) + u_mortar_local = project_to_local * displacement_jump_global_coord + u_mortar_local_decomposed = u_mortar_local.reshape((2, -1), order="F") + + contact_force = contact_force.reshape((2, -1), order="F") + + return u_mortar_local_decomposed, contact_force, fracture_pressure + + def test_pull_north_positive_opening(self): + + setup = SetupContactMechanicsBiot( + ux_south=0, uy_south=0, ux_north=0, uy_north=0.001 + ) + + u_mortar, contact_force, fracture_pressure = self._solve(setup) + + # All components should be open in the normal direction + self.assertTrue(np.all(u_mortar[1] < 0)) + + # By symmetry (reasonable to expect from this grid), the jump in tangential + # deformation should be zero. + self.assertTrue(np.abs(np.sum(u_mortar[0])) < 1e-5) + + # The contact force in normal direction should be zero + + # NB: This assumes the contact force is expressed in local coordinates + self.assertTrue(np.all(np.abs(contact_force) < 1e-7)) + + # Check that the dilation of the fracture yields a negative fracture pressure + self.assertTrue(np.all(fracture_pressure < -1e-7)) + + def test_pull_south_positive_opening(self): + + setup = SetupContactMechanicsBiot( + ux_south=0, uy_south=-0.001, ux_north=0, uy_north=0 + ) + + u_mortar, contact_force, fracture_pressure = self._solve(setup) + + # All components should be open in the normal direction + self.assertTrue(np.all(u_mortar[1] < 0)) + + # By symmetry (reasonable to expect from this grid), the jump in tangential + # deformation should be zero. + self.assertTrue(np.abs(np.sum(u_mortar[0])) < 1e-5) + + # The contact force in normal direction should be zero + + # NB: This assumes the contact force is expressed in local coordinates + self.assertTrue(np.all(np.abs(contact_force) < 1e-7)) + + # Check that the dilation of the fracture yields a negative fracture pressure + self.assertTrue(np.all(fracture_pressure < -1e-7)) + + def test_push_north_zero_opening(self): + + setup = SetupContactMechanicsBiot( + ux_south=0, uy_south=0, ux_north=0, uy_north=-0.001 + ) + + u_mortar, contact_force, fracture_pressure = self._solve(setup) + + # All components should be closed in the normal direction + self.assertTrue(np.abs(np.sum(u_mortar[1])) < 1e-5) + + # Contact force in normal direction should be negative + self.assertTrue(np.all(contact_force[1] < 0)) + + # Compression of the domain yields a (slightly) positive fracture pressure + self.assertTrue(np.all(fracture_pressure > 1e-10)) + + def test_push_south_zero_opening(self): + + setup = SetupContactMechanicsBiot( + ux_south=0, uy_south=0.001, ux_north=0, uy_north=0 + ) + + u_mortar, contact_force, fracture_pressure = self._solve(setup) + + # All components should be closed in the normal direction + self.assertTrue(np.abs(np.sum(u_mortar[1])) < 1e-5) + + # Contact force in normal direction should be negative + self.assertTrue(np.all(contact_force[1] < 0)) + + # Compression of the domain yields a (slightly) positive fracture pressure + self.assertTrue(np.all(fracture_pressure > 1e-10)) + + def test_positive_fracture_pressure_positive_opening(self): + + setup = SetupContactMechanicsBiot( + ux_south=0, uy_south=0, ux_north=0, uy_north=0, source_value=0.001 + ) + + u_mortar, contact_force, fracture_pressure = self._solve(setup) + + # All components should be open in the normal direction + self.assertTrue(np.all(u_mortar[1] < 0)) + + # By symmetry (reasonable to expect from this grid), the jump in tangential + # deformation should be zero. + self.assertTrue(np.abs(np.sum(u_mortar[0])) < 1e-5) + + # The contact force in normal direction should be zero + + # NB: This assumes the contact force is expressed in local coordinates + self.assertTrue(np.all(np.abs(contact_force) < 1e-7)) + + # Fracture pressure is positive + self.assertTrue(np.all(fracture_pressure > 1e-7)) + + +def distribute_iterate( + assembler, setup, values, mortar_displacement_variable, contact_traction_variable +): + """ Update the previous iterate of the mortar displacement and contact traction, + and obtain current matrix displacement iterate. + + Method is a tailored copy from assembler.distribute_variable. + """ + dof = np.cumsum(np.append(0, np.asarray(assembler.full_dof))) + var_name = setup.displacement_variable + + for pair, bi in assembler.block_dof.items(): + g = pair[0] + name = pair[1] + # Avoid edges + if not isinstance(g, pp.Grid): + if name == mortar_displacement_variable: + mortar_u = values[dof[bi] : dof[bi + 1]] + data = setup.gb.edge_props(g) + data[pp.STATE]["previous_iterate"][ + mortar_displacement_variable + ] = mortar_u + continue + # Only interested in highest dimension + if g.dim < setup.gb.dim_max(): + if name == contact_traction_variable: + contact = values[dof[bi] : dof[bi + 1]] + data = setup.gb.node_props(g) + data[pp.STATE]["previous_iterate"][contact_traction_variable] = contact + + continue + # Only need the displacement + if name != var_name: + continue + + u = values[dof[bi] : dof[bi + 1]] + return u + + +class SetupContactMechanicsBiot(model.ContactMechanicsBiot): + def __init__(self, ux_south, uy_south, ux_north, uy_north, source_value=0): + mesh_args = { + "mesh_size_frac": 0.5, + "mesh_size_min": 0.023, + "mesh_size_bound": 0.5, + } + super().__init__(mesh_args, folder_name="dummy") + self.ux_south = ux_south + self.uy_south = uy_south + self.ux_north = ux_north + self.uy_north = uy_north + self.scalar_source_value = source_value + + def create_grid(self, rotate_fracture=False): + """ + Method that creates and returns the GridBucket of a 2D domain with six + fractures. The two sides of the fractures are coupled together with a + mortar grid. + """ + if rotate_fracture: + self.frac_pts = np.array([[0.7, 0.3], [0.3, 0.7]]) + else: + self.frac_pts = np.array([[0.3, 0.7], [0.5, 0.5]]) + frac_edges = np.array([[0], [1]]) + + self.box = {"xmin": 0, "ymin": 0, "xmax": 1, "ymax": 1} + + network = pp.FractureNetwork2d(self.frac_pts, frac_edges, domain=self.box) + # Generate the mixed-dimensional mesh + gb = network.mesh(self.mesh_args) + + # Set projections to local coordinates for all fractures + pp.contact_conditions.set_projections(gb) + + self.gb = gb + self.Nd = gb.dim_max() + + def source(self, g, key, t=0): + if key == self.mechanics_parameter_key: + values = super().source(g, key, t) + elif key == self.scalar_parameter_key: + if g.dim == self.Nd: + values = np.zeros(g.num_cells) + else: + values = self.scalar_source_value * np.ones(g.num_cells) + else: + raise ValueError("No BC values implemented for keyword " + str(key)) + return values + + def bc_type(self, g, key, t=0): + _, _, _, north, south, _, _ = self.domain_boundary_sides(g) + if key == self.mechanics_parameter_key: + bc = pp.BoundaryConditionVectorial(g, north + south, "dir") + elif key == self.scalar_parameter_key: + # Define boundary condition on faces + bc = pp.BoundaryCondition(g, north + south, "dir") + else: + raise ValueError("No BCs implemented for keyword " + str(key)) + return bc + + def bc_values(self, g, key, t=0): + # Set the boundary values + if key == self.mechanics_parameter_key: + _, _, _, north, south, _, _ = self.domain_boundary_sides(g) + values = np.zeros((g.dim, g.num_faces)) + values[0, south] = self.ux_south + values[1, south] = self.uy_south + values[0, north] = self.ux_north + values[1, north] = self.uy_north + values = values.ravel("F") + elif key == self.scalar_parameter_key: + values = np.zeros(g.num_faces) + else: + raise ValueError("No BC values implemented for keyword " + str(key)) + return values + + +if __name__ == "__main__": + unittest.main()