diff --git a/examples/example5/test_upwind_examples.py b/examples/example5/test_upwind_examples.py index 50741b89d1..4d0f974440 100644 --- a/examples/example5/test_upwind_examples.py +++ b/examples/example5/test_upwind_examples.py @@ -38,6 +38,8 @@ def test_upwind_example0(self, if_export=False): time_step = advect.cfl(g, data) data[pp.PARAMETERS]["transport"]["time_step"] = time_step + advect.discretize(g, data) + U, rhs = advect.assemble_matrix_rhs(g, data) rhs = time_step * rhs U = time_step * U @@ -92,6 +94,7 @@ def test_upwind_example1(self, if_export=False): time_step = advect.cfl(g, data) data[pp.PARAMETERS]["transport"]["time_step"] = time_step + advect.discretize(g, data) U, rhs = advect.assemble_matrix_rhs(g, data) rhs = time_step * rhs U = time_step * U @@ -168,6 +171,7 @@ def funp_ex(pt): } data = pp.initialize_default_data(g, {}, "flow", specified_parameters) solver = pp.MVEM("flow") + solver.discretize(g, data) D_flow, b_flow = solver.assemble_matrix_rhs(g, data) solver_source = pp.DualScalarSource("flow") @@ -194,6 +198,7 @@ def funp_ex(pt): # Advect solver advect = pp.Upwind("transport") + advect.discretize(g, data) U, rhs = advect.assemble_matrix_rhs(g, data) time_step = advect.cfl(g, data) diff --git a/src/porepy/__init__.py b/src/porepy/__init__.py index 30c1c49891..fb3c7e6f6f 100644 --- a/src/porepy/__init__.py +++ b/src/porepy/__init__.py @@ -26,12 +26,14 @@ __all__ = [] # Numerics +from porepy.numerics.discretization import VoidDiscretization + # Control volume, elliptic -from porepy.numerics.fv.mpsa import Mpsa, FracturedMpsa +from porepy.numerics.fv.mpsa import Mpsa from porepy.numerics.fv.fv_elliptic import FVElliptic from porepy.numerics.fv.tpfa import Tpfa from porepy.numerics.fv.mpfa import Mpfa -from porepy.numerics.fv.biot import Biot, GradP, DivD, BiotStabilization +from porepy.numerics.fv.biot import Biot, GradP, DivU, BiotStabilization from porepy.numerics.fv.source import ScalarSource # Virtual elements, elliptic @@ -50,9 +52,8 @@ from porepy.numerics.interface_laws.elliptic_interface_laws import ( RobinCoupling, FluxPressureContinuity, - RobinContact, - StressDisplacementContinuity, ) + from porepy.numerics.interface_laws.cell_dof_face_dof_map import CellDofFaceDofMap from porepy.numerics.mixed_dim.assembler import Assembler @@ -64,6 +65,16 @@ from porepy.numerics.fv.mass_matrix import MassMatrix from porepy.numerics.fv.mass_matrix import InvMassMatrix +# Contact mechanics +from porepy.numerics.interface_laws.contact_mechanics_interface_laws import ( + PrimalContactCoupling, + DivUCoupling, + MatrixScalarToForceBalance, + FractureScalarToForceBalance, +) +from porepy.numerics.contact_mechanics.contact_conditions import ColoumbContact +from porepy.numerics.contact_mechanics import contact_conditions + # Grids from porepy.grids.grid import Grid from porepy.grids.fv_sub_grid import FvSubGrid @@ -114,7 +125,8 @@ from porepy.fracs import meshing, fracture_importer, mortars from porepy.grids import structured, simplex, coarsening, partition, refinement from porepy.numerics.fv import fvutils -from porepy.utils import error +from porepy.utils import error, grid_utils +from porepy.utils.tangential_normal_projection import TangentialNormalProjection # Constants, units and keywords from porepy.utils.common_constants import * diff --git a/src/porepy/geometry/map_geometry.py b/src/porepy/geometry/map_geometry.py index 05d04c0803..8ef75c70d0 100644 --- a/src/porepy/geometry/map_geometry.py +++ b/src/porepy/geometry/map_geometry.py @@ -31,13 +31,18 @@ def force_point_collinearity(pts): return pts[:, 0, np.newaxis] * (1 - dist) + pts[:, end, np.newaxis] * dist -def map_grid(g, tol=1e-5): +def map_grid(g, tol=1e-5, R=None): """ If a 2d or a 1d grid is passed, the function return the cell_centers, face_normals, and face_centers using local coordinates. If a 3d grid is passed nothing is applied. The return vectors have a reduced number of rows. Parameters: g (grid): the grid. + tol (double, optional): Tolerance used to check that the grid is linear or planar. + Defaults to 1e-5. + R (np.array size 3x3, optional ): Rotation matrix. The first dim rows should map + vectors onto the tangential space of the grid. If not provided, a rotation + matrix will be computed. Returns: cell_centers: (g.dim x g.num_cells) the mapped centers of the cells. @@ -52,9 +57,11 @@ def map_grid(g, tol=1e-5): face_normals = g.face_normals face_centers = g.face_centers nodes = g.nodes - R = np.eye(3) if g.dim == 0 or g.dim == 3: + if R is None: + R = np.eye(3) + return ( cell_centers, face_normals, @@ -64,12 +71,12 @@ def map_grid(g, tol=1e-5): nodes, ) - if g.dim == 1 or g.dim == 2: - - if g.dim == 2: - R = project_plane_matrix(g.nodes, tol=tol) - else: - R = project_line_matrix(g.nodes, tol=tol) + else: # g.dim == 1 or g.dim == 2: + if R is None: + if g.dim == 2: + R = project_plane_matrix(g.nodes, tol=tol) + else: + R = project_line_matrix(g.nodes, tol=tol) face_centers = np.dot(R, face_centers) diff --git a/src/porepy/grids/grid.py b/src/porepy/grids/grid.py index 239ce8119b..37834d77e8 100644 --- a/src/porepy/grids/grid.py +++ b/src/porepy/grids/grid.py @@ -743,6 +743,34 @@ def cell_connection_map(self): return c2c + def sign_of_faces(self, faces): + """ Get the direction of the normal vector (inward or outwards from a cell) + of faces. Only boundary faces are permissible. + + Parameters: + faces: (ndarray) indices of faces that you want to know the sign for. The + faces must be boundary faces. + + Returns: + (ndarray) the sign of the faces + + Raises: + ValueError if a target face is internal. + + """ + + IA = np.argsort(faces) + IC = np.argsort(IA) + + fi, _, sgn = sps.find(self.cell_faces[faces[IA], :]) + if fi.size != faces.size: + raise ValueError("sign of internal faces does not make sense") + + I = np.argsort(fi) + sgn = sgn[I] + sgn = sgn[IC] + return sgn + def bounding_box(self): """ Return the bounding box of the grid. diff --git a/src/porepy/grids/mortar_grid.py b/src/porepy/grids/mortar_grid.py index 81bcae896d..78d495d533 100644 --- a/src/porepy/grids/mortar_grid.py +++ b/src/porepy/grids/mortar_grid.py @@ -1,9 +1,11 @@ """ Module containing the class for the mortar grid. """ from __future__ import division +import warnings import numpy as np from scipy import sparse as sps + # Module level constants, used to define sides of a mortar grid. # This is in essence an Enum, but that led to trouble in pickling a GridBucket. NONE_SIDE = 0 @@ -96,7 +98,7 @@ def __init__(self, dim, side_grids, face_cells, name=""): # to the co-dimensional grid. If this assumption is not satisfied we # need to change the following lines - # Creation of the high_to_mortar_int, besically we start from the face_cells + # Creation of the high_to_mortar_int, basically we start from the face_cells # map and we split the relation # low_dimensional_cell -> 2 high_dimensional_face # as @@ -199,7 +201,7 @@ def compute_geometry(self): [g.cell_centers for g in self.side_grids.values()] ) - # ------------------------------------------------------------------------------# + ### Methods to update the mortar grid, or the neighboring grids. def update_mortar(self, side_matrix, side_grids): """ @@ -237,8 +239,6 @@ def update_mortar(self, side_matrix, side_grids): self._check_mappings() - # ------------------------------------------------------------------------------# - def update_slave(self, side_matrix): """ Update the low_to_mortar_int map when the lower dimensional grid is changed. @@ -262,15 +262,11 @@ def update_slave(self, side_matrix): self._slave_to_mortar_int = sps.bmat(matrix, format="csc") self._check_mappings() - # ------------------------------------------------------------------------------# - def update_master(self, matrix): # Make a comment here self._master_to_mortar_int = self._master_to_mortar_int * matrix self._check_mappings() - # ------------------------------------------------------------------------------# - def num_sides(self): """ Shortcut to compute the number of sides, it has to be 2 or 1. @@ -280,9 +276,35 @@ def num_sides(self): """ return len(self.side_grids) - # ------------------------------------------------------------------------------# + ### + def project_to_side_grids(self): + """ Generator for the side grids (pp.Grid) representation of the mortar + cells, and projection operators from the mortar cells, combining cells on all + the sides, to the specific side grids. + + Yields: + grid (pp.Grid): PorePy grid representing one of the sides of the + mortar grid. Can be used for standard discretizations. + proj (sps.csc_matrix): Projection from the mortar cells to this + side grid. - def master_to_mortar_int(self): + """ + counter = 0 + for grid in self.side_grids.values(): + nc = grid.num_cells + rows = np.arange(nc) + cols = rows + counter + data = np.ones(nc) + proj = sps.coo_matrix( + (data, (rows, cols)), shape=(nc, self.num_cells) + ).tocsc() + + counter += nc + yield proj, grid + + ## Methods to construct projection matrices + + def master_to_mortar_int(self, nd=1): """ Project values from faces of master to the mortar, by summing quantities from the master side. @@ -292,14 +314,18 @@ def master_to_mortar_int(self): This mapping is intended for extensive properties, e.g. fluxes. + Parameters: + nd (int, optional): Spatial dimension of the projected quantity. + Defaults to 1 (mapping for scalar quantities). + Returns: sps.matrix: Projection matrix with column sum unity. Size: g_master.num_faces x mortar_grid.num_cells. """ - return self._master_to_mortar_int + return self._convert_to_vector_variable(self._master_to_mortar_int, nd) - def slave_to_mortar_int(self): + def slave_to_mortar_int(self, nd=1): """ Project values from cells on the slave side to the mortar, by summing quantities from the slave side. @@ -309,14 +335,18 @@ def slave_to_mortar_int(self): This mapping is intended for extensive properties, e.g. sources. + Parameters: + nd (int, optional): Spatial dimension of the projected quantity. + Defaults to 1 (mapping for scalar quantities). + Returns: sps.matrix: Projection matrix with column sum unity. Size: g_slave.num_cells x mortar_grid.num_cells. """ - return self._slave_to_mortar_int + return self._convert_to_vector_variable(self._slave_to_mortar_int, nd) - def master_to_mortar_avg(self): + def master_to_mortar_avg(self, nd=1): """ Project values from faces of master to the mortar, by averaging quantities from the master side. @@ -327,15 +357,21 @@ def master_to_mortar_avg(self): This mapping is intended for intensive properties, e.g. pressures. + Parameters: + nd (int, optional): Spatial dimension of the projected quantity. + Defaults to 1 (mapping for scalar quantities). + Returns: sps.matrix: Projection matrix with row sum unity. Size: g_master.num_faces x mortar_grid.num_cells. """ row_sum = self._master_to_mortar_int.sum(axis=1).A.ravel() - return sps.diags(1.0 / row_sum) * self._master_to_mortar_int + return self._convert_to_vector_variable( + sps.diags(1.0 / row_sum) * self._master_to_mortar_int, nd + ) - def slave_to_mortar_avg(self): + def slave_to_mortar_avg(self, nd=1): """ Project values from cells at the slave to the mortar, by averaging quantities from the slave side. @@ -346,19 +382,25 @@ def slave_to_mortar_avg(self): This mapping is intended for intensive properties, e.g. pressures. + Parameters: + nd (int, optional): Spatial dimension of the projected quantity. + Defaults to 1 (mapping for scalar quantities). + Returns: sps.matrix: Projection matrix with row sum unity. Size: g_slave.num_cells x mortar_grid.num_cells. """ row_sum = self._slave_to_mortar_int.sum(axis=1).A.ravel() - return sps.diags(1.0 / row_sum) * self._slave_to_mortar_int + return self._convert_to_vector_variable( + sps.diags(1.0 / row_sum) * self._slave_to_mortar_int, nd + ) # IMPLEMENTATION NOTE: The reverse projections, from mortar to master/slave are # found by taking transposes, and switching average and integration (since we are # changing which side we are taking the area relative to. - def mortar_to_master_int(self): + def mortar_to_master_int(self, nd=1): """ Project values from the mortar to faces of master, by summing quantities from the mortar side. @@ -368,14 +410,18 @@ def mortar_to_master_int(self): This mapping is intended for extensive properties, e.g. fluxes. + Parameters: + nd (int, optional): Spatial dimension of the projected quantity. + Defaults to 1 (mapping for scalar quantities). + Returns: sps.matrix: Projection matrix with column sum unity. Size: mortar_grid.num_cells x g_master.num_faces. """ - return self.master_to_mortar_avg().T + return self._convert_to_vector_variable(self.master_to_mortar_avg().T, nd) - def mortar_to_slave_int(self): + def mortar_to_slave_int(self, nd=1): """ Project values from the mortar to cells at the slave, by summing quantities from the mortar side. @@ -385,14 +431,18 @@ def mortar_to_slave_int(self): This mapping is intended for extensive properties, e.g. fluxes. + Parameters: + nd (int, optional): Spatial dimension of the projected quantity. + Defaults to 1 (mapping for scalar quantities). + Returns: sps.matrix: Projection matrix with column sum unity. Size: mortar_grid.num_cells x g_slave_num_faces. """ - return self.slave_to_mortar_avg().T + return self._convert_to_vector_variable(self.slave_to_mortar_avg().T, nd) - def mortar_to_master_avg(self): + def mortar_to_master_avg(self, nd=1): """ Project values from the mortar to faces of master, by averaging quantities from the mortar side. @@ -403,14 +453,18 @@ def mortar_to_master_avg(self): This mapping is intended for intensive properties, e.g. pressures. + Parameters: + nd (int, optional): Spatial dimension of the projected quantity. + Defaults to 1 (mapping for scalar quantities). + Returns: sps.matrix: Projection matrix with row sum unity. Size: mortar_grid.num_cells x g_master.num_faces. """ - return self.master_to_mortar_int().T + return self._convert_to_vector_variable(self.master_to_mortar_int().T, nd) - def mortar_to_slave_avg(self): + def mortar_to_slave_avg(self, nd=1): """ Project values from the mortar to slave, by averaging quantities from the mortar side. @@ -421,14 +475,72 @@ def mortar_to_slave_avg(self): This mapping is intended for intensive properties, e.g. pressures. + Parameters: + nd (int, optional): Spatial dimension of the projected quantity. + Defaults to 1 (mapping for scalar quantities). + Returns: sps.matrix: Projection matrix with row sum unity. Size: mortar_grid.num_cells x g_slave.num_faces. """ - return self.slave_to_mortar_int().T + return self._convert_to_vector_variable(self.slave_to_mortar_int().T, nd) - # ------------------------------------------------------------------------------# + def _convert_to_vector_variable(self, matrix, nd): + """ Convert the scalar projection to a vector quantity. If the prescribed + dimension is 1 (default for all the above methods), the projection matrix + will in effect not be altered. + """ + return sps.kron(matrix, sps.eye(nd)).tocsc() + + def sign_of_mortar_sides(self, nd=1): + """ Assign positive or negative weight to the two sides of a mortar grid. + + This is needed e.g. to make projection operators into signed projections, + for variables that have no particular defined sign conventions. + + This function defines a convention for what is a positive jump between + the mortar sides. + + Example: Take the difference between right and left variables, and + project to the slave grid by + + mortar_to_slave_avg() * sign_of_mortar_sides() + + NOTE: The flux variables in flow and transport equations are defined as + positive from master to slave. Hence the two sides have different + conventions, and there is no need to adjust the signs further. + + IMPLEMENTATION NOTE: This method will probably not be meaningful if + applied to mortar grids where the two side grids are non-matching. + + Parameters: + nd (int, optional): Spatial dimension of the projected quantity. + Defaults to 1 (mapping for scalar quantities). + + Returns: + sps.diag_matrix: Diagonal matrix with positive signs on variables + belonging to the first of the side_grids. + Size: mortar_grid.num_cells x mortar_grid.num_cells + + """ + nc = self.num_cells + if self.num_sides() == 1: + warnings.warn( + "Is it really meaningful to ask for signs of a one sided mortar grid?" + ) + return sps.dia_matrix((np.ones(nc * nd), 0), shape=(nd * nc, nd * nc)) + elif self.num_sides() == 2: + # From the numbering of the mortar cells (see __init__, the case + # num_sides() == 2)), we know that the cells are numbered first + # on one side, then on the other. + data = np.hstack( + ( + np.ones(self.side_grids[1].num_cells * nd), + -np.ones(self.side_grids[2].num_cells * nd), + ) + ) + return sps.dia_matrix((data, 0), shape=(nd * nc, nd * nc)) def cell_diameters(self): diams = np.empty(self.num_sides(), dtype=np.object) @@ -436,8 +548,6 @@ def cell_diameters(self): diams[pos] = g.cell_diameters() return np.concatenate(diams).ravel() - # ------------------------------------------------------------------------------# - def _check_mappings(self, tol=1e-4): row_sum = self._master_to_mortar_int.sum(axis=1) if not (row_sum.min() > tol): diff --git a/src/porepy/grids/partition.py b/src/porepy/grids/partition.py index ed69945148..055ea97848 100644 --- a/src/porepy/grids/partition.py +++ b/src/porepy/grids/partition.py @@ -49,9 +49,10 @@ def partition_metis(g, num_part): c2c = g.cell_connection_map() # Convert the cells into the format required by pymetis - adjacency_list = [c2c.getrow(i).indices for i in range(c2c.shape[0])] + adjacency_list = [list(c2c.getrow(i).indices) for i in range(c2c.shape[0])] # Call pymetis - part = pymetis.part_graph(num_part, adjacency=adjacency_list) + # It seems it is important that num_part is an int, not an np.int. + part = pymetis.part_graph(int(num_part), adjacency=adjacency_list) # The meaning of the first number returned by pymetis is not clear (poor # documentation), only return the partitioning. diff --git a/src/porepy/models/__init__.py b/src/porepy/models/__init__.py new file mode 100644 index 0000000000..e69de29bb2 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..1f30198e0e --- /dev/null +++ b/src/porepy/models/contact_mechanics_biot_model.py @@ -0,0 +1,462 @@ +""" +This is a setup class for solving the Biot equations with contact between the fractures. + +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 + +import porepy.models.contact_mechanics_model as contact_model +from porepy.utils.derived_discretizations import implicit_euler as IE_discretizations + + +class ContactMechanicsBiot(contact_model.ContactMechanics): + def __init__(self, mesh_args, folder_name): + super().__init__(mesh_args, folder_name) + + # Temperature + self.scalar_variable = "p" + self.mortar_scalar_variable = "mortar_" + self.scalar_variable + self.scalar_coupling_term = "robin_" + self.scalar_variable + self.scalar_parameter_key = "flow" + + # Scaling coefficients + self.scalar_scale = 1 + self.length_scale = 1 + + # Time + # The time attribute may be used e.g. to update BCs. + self.time = 0 + self.time_step = 1 * self.length_scale ** 2 + self.end_time = self.time_step * 1 + + # 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 + self.subtract_fracture_pressure = True + + def bc_type_mechanics(self, g): + # Use parent class method for mechanics + return super().bc_type(g) + + def bc_type_scalar(self, g): + # Define boundary regions + all_bf, *_ = self.domain_boundary_sides(g) + # Define boundary condition on faces + return pp.BoundaryCondition(g, all_bf, "dir") + + def bc_values_mechanics(self, g): + # Set the boundary values + return super().bc_values(g) + + def bc_values_scalar(self, g): + return np.zeros(g.num_faces) + + def source_mechanics(self, g): + return super().source(g) + + def source_scalar(self, g): + return np.zeros(g.num_cells) + + 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. + """ + self.set_scalar_parameters() + self.set_mechanics_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) / 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) + # BC and source values + bc_val = self.bc_values_mechanics(g) + source_val = self.source_mechanics(g) + + 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 _, 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_scalar(g) + bc_values = self.bc_values_scalar(g) + source_values = self.source_scalar(g) + + 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, _ = 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(self): + """ + Assign discretizations to the nodes and edges of the grid bucket. + + Note the attribute subtract_fracture_pressure: Indicates whether or not to + subtract the fracture pressure contribution for the contact traction. This + should not be done if the scalar variable is temperature. + """ + # Shorthand + 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 Nd domain we solve linear elasticity with mpsa. + mpsa = pp.Mpsa(key_m) + empty_discr = pp.VoidDiscretization(key_m, ndof_cell=self.Nd) + # Scalar discretizations (all dimensions) + diff_disc_s = IE_discretizations.ImplicitMpfa(key_s) + mass_disc_s = IE_discretizations.ImplicitMassMatrix(key_s, var_s) + source_disc_s = pp.ScalarSource(key_s) + # Coupling discretizations + # All dimensions + div_u_disc = pp.DivU( + key_m, + key_s, + variable=var_d, + mortar_variable=self.mortar_displacement_variable, + ) + # Nd + grad_p_disc = pp.GradP(key_m) + stabilization_disc_s = pp.BiotStabilization(key_s, var_s) + + # Assign node discretizations + for g, d in self.gb: + if g.dim == self.Nd: + d[pp.DISCRETIZATION] = { + var_d: {"mpsa": mpsa}, + var_s: { + "diffusion": diff_disc_s, + "mass": mass_disc_s, + "stabilization": stabilization_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 == self.Nd - 1: + d[pp.DISCRETIZATION] = { + self.contact_traction_variable: {"empty": empty_discr}, + var_s: { + "diffusion": diff_disc_s, + "mass": mass_disc_s, + "source": source_disc_s, + }, + } + + # Define edge discretizations for the mortar grid + contact_law = pp.ColoumbContact(self.mechanics_parameter_key, self.Nd) + contact_discr = pp.PrimalContactCoupling( + self.mechanics_parameter_key, mpsa, contact_law + ) + # Account for the mortar displacements effect on scalar balance in the matrix, + # as an internal boundary contribution, fracture, aperture changes appear as a + # source contribution. + div_u_coupling = pp.DivUCoupling( + self.displacement_variable, div_u_disc, div_u_disc + ) + # Account for the pressure contributions to the force balance on the fracture + # (see contact_discr). + # This discretization needs the keyword used to store the grad p discretization: + grad_p_key = key_m + matrix_scalar_to_force_balance = pp.MatrixScalarToForceBalance( + grad_p_key, mass_disc_s, mass_disc_s + ) + if self.subtract_fracture_pressure: + fracture_scalar_to_force_balance = pp.FractureScalarToForceBalance( + mass_disc_s, mass_disc_s + ) + + for e, d in self.gb.edges(): + g_l, g_h = self.gb.nodes_of_edge(e) + + if g_h.dim == self.Nd: + d[pp.COUPLING_DISCRETIZATION] = { + self.friction_coupling_term: { + g_h: (var_d, "mpsa"), + g_l: (self.contact_traction_variable, "empty"), + (g_h, g_l): (self.mortar_displacement_variable, contact_discr), + }, + self.scalar_coupling_term: { + g_h: (var_s, "diffusion"), + g_l: (var_s, "diffusion"), + e: ( + self.mortar_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.mortar_displacement_variable, div_u_coupling), + }, + "matrix_scalar_to_force_balance": { + g_h: (var_s, "mass"), + g_l: (var_s, "mass"), + e: ( + self.mortar_displacement_variable, + matrix_scalar_to_force_balance, + ), + }, + } + if self.subtract_fracture_pressure: + d[pp.COUPLING_DISCRETIZATION].update( + { + "matrix_scalar_to_force_balance": { + g_h: (var_s, "mass"), + g_l: (var_s, "mass"), + e: ( + self.mortar_displacement_variable, + fracture_scalar_to_force_balance, + ), + } + } + ) + else: + raise ValueError( + "assign_discretizations assumes no fracture intersections." + ) + + def assign_variables(self): + """ + Assign primary variables to the nodes and edges of the grid bucket. + """ + # First for the nodes + for g, d in self.gb: + if g.dim == self.Nd: + d[pp.PRIMARY_VARIABLES] = { + self.displacement_variable: {"cells": self.Nd}, + self.scalar_variable: {"cells": 1}, + } + elif g.dim == self.Nd - 1: + d[pp.PRIMARY_VARIABLES] = { + self.contact_traction_variable: {"cells": self.Nd}, + self.scalar_variable: {"cells": 1}, + } + else: + d[pp.PRIMARY_VARIABLES] = {} + # Then for the edges + for e, d in self.gb.edges(): + _, g_h = self.gb.nodes_of_edge(e) + + if g_h.dim == self.Nd: + d[pp.PRIMARY_VARIABLES] = { + self.mortar_displacement_variable: {"cells": self.Nd}, + self.mortar_scalar_variable: {"cells": 1}, + } + + def discretize_biot(self, gb): + """ + To save computational time, the full Biot equation (without contact mechanics) + is discretized once. This is to avoid 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 = np.zeros(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 export_step(self): + pass + + def export_pvd(self): + pass + + +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. + + 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 + newton_tol: Tolerance for the Newton solver, see contact_mechanics_model. + """ + if "gb" not in setup.__dict__: + setup.create_grid() + gb = setup.gb + + # Extract the grids we use + g_max = gb.grids_of_dimension(setup.Nd)[0] + d_max = gb.node_props(g_max) + + # Assign parameters, variables and discretizations + setup.set_parameters() + setup.initial_condition() + setup.assign_variables() + setup.assign_discretisations() + # Set up assembler and get initial condition for the displacements + 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 = [] + dt = setup.time_step + t_end = setup.end_time + k = 0 + while setup.time < t_end: + setup.time += dt + k += 1 + print("Time step: ", k, "/", int(np.ceil(t_end / dt))) + + # Prepare for Newton + counter_newton = 0 + converged_newton = False + 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, 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() diff --git a/src/porepy/models/contact_mechanics_model.py b/src/porepy/models/contact_mechanics_model.py new file mode 100644 index 0000000000..508f8ace5c --- /dev/null +++ b/src/porepy/models/contact_mechanics_model.py @@ -0,0 +1,486 @@ +""" +This is a setup class for solving linear elasticity with contact between the fractures. + +The setup handles parameters, variables and discretizations. Default (unitary-like) +parameters are set. A "run script" function for setting up the class and solving the +nonlinear contact mechanics problem is also provided. + +NOTE: This module should be considered an experimental feature, which will likely +undergo major changes (or be deleted). + +""" +import numpy as np +import scipy.sparse as sps +from scipy.spatial.distance import cdist + +import porepy as pp + + +class ContactMechanics: + def __init__(self, mesh_args, folder_name): + self.mesh_args = mesh_args + self.folder_name = folder_name + + # Variables + self.displacement_variable = "u" + self.mortar_displacement_variable = "mortar_u" + self.contact_traction_variable = "contact_traction" + + # Keyword + self.mechanics_parameter_key = "mechanics" + + # Terms of the equations + self.friction_coupling_term = "fracture_force_balance" + + def create_grid(self): + """ + Method that creates a GridBucket of a 2D domain with one fracture and sets + projections to local coordinates for all fractures. + + The method requires the following attribute: + mesh_args (dict): Containing the mesh sizes. + + The method assigns the following attributes to self: + frac_pts (np.array): Nd x (number of fracture points), the coordinates of + the fracture endpoints. + box (dict): The bounding box of the domain, defined through minimum and + maximum values in each dimension. + 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]]) + # Each column defines one fracture + 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 = self.gb.dim_max() + + def domain_boundary_sides(self, g): + """ + Obtain indices of the faces of a grid that lie on each side of the domain + boundaries. + """ + tol = 1e-10 + box = self.box + east = g.face_centers[0] > box["xmax"] - tol + west = g.face_centers[0] < box["xmin"] + tol + north = g.face_centers[1] > box["ymax"] - tol + south = g.face_centers[1] < box["ymin"] + tol + if self.Nd == 2: + top = np.zeros(g.num_faces, dtype=bool) + bottom = top.copy() + else: + top = g.face_centers[2] > box["zmax"] - tol + bottom = g.face_centers[2] < box["zmin"] + tol + all_bf = g.get_boundary_faces() + return all_bf, east, west, north, south, top, bottom + + 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): + # Values for all Nd components, facewise + values = np.zeros((self.Nd, g.num_faces)) + # Reshape according to PorePy convention + values = values.ravel("F") + return values + + def source(self, g): + return 0 + + def set_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) + + # BC and source values + bc_val = self.bc_values(g) + source_val = self.source(g) + + pp.initialize_data( + g, + d, + self.mechanics_parameter_key, + { + "bc": bc, + "bc_values": bc_val, + "source": source_val, + "fourth_order_tensor": C, + }, + ) + + elif g.dim == self.Nd - 1: + friction = self._set_friction_coefficient(g) + pp.initialize_data( + g, + d, + self.mechanics_parameter_key, + {"friction_coefficient": friction}, + ) + # Should we keep this, @EK? + for _, 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 assign_variables(self): + """ + Assign variables to the nodes and edges of the grid bucket. + """ + gb = self.gb + for g, d in gb: + if g.dim == self.Nd: + d[pp.PRIMARY_VARIABLES] = { + self.displacement_variable: {"cells": self.Nd} + } + elif g.dim == self.Nd - 1: + d[pp.PRIMARY_VARIABLES] = { + self.contact_traction_variable: {"cells": self.Nd} + } + else: + d[pp.PRIMARY_VARIABLES] = {} + + for e, d in gb.edges(): + + if e[0].dim == self.Nd: + d[pp.PRIMARY_VARIABLES] = { + self.mortar_displacement_variable: {"cells": self.Nd} + } + + else: + d[pp.PRIMARY_VARIABLES] = {} + + def assign_discretizations(self): + """ + Assign discretizations to the nodes and edges of the grid bucket. + """ + # For the Nd domain we solve linear elasticity with mpsa. + Nd = self.Nd + gb = self.gb + mpsa = pp.Mpsa(self.mechanics_parameter_key) + # We need a void discretization for the contact traction variable defined on + # the fractures. + empty_discr = pp.VoidDiscretization(self.mechanics_parameter_key, ndof_cell=Nd) + + for g, d in gb: + if g.dim == Nd: + d[pp.DISCRETIZATION] = {self.displacement_variable: {"mpsa": mpsa}} + elif g.dim == Nd - 1: + d[pp.DISCRETIZATION] = { + self.contact_traction_variable: {"empty": empty_discr} + } + + # Define the contact condition on the mortar grid + coloumb = pp.ColoumbContact(self.mechanics_parameter_key, Nd) + contact = pp.PrimalContactCoupling(self.mechanics_parameter_key, mpsa, coloumb) + + for e, d in gb.edges(): + g_l, g_h = gb.nodes_of_edge(e) + if g_h.dim == Nd: + d[pp.COUPLING_DISCRETIZATION] = { + self.friction_coupling_term: { + g_h: (self.displacement_variable, "mpsa"), + g_l: (self.contact_traction_variable, "empty"), + (g_h, g_l): (self.mortar_displacement_variable, contact), + } + } + + def initial_condition(self): + """ + Initial guess for Newton iteration. + """ + + for g, d in self.gb: + if g.dim == self.Nd: + # Initialize displacement variable + state = {self.displacement_variable: np.zeros(g.num_cells * self.Nd)} + + elif g.dim == self.Nd - 1: + # Initialize contact variable + traction = np.vstack( + (np.zeros((g.dim, g.num_cells)), -1 * np.ones(g.num_cells)) + ).ravel(order="F") + state = {"previous_iterate": {self.contact_traction_variable: traction}} + else: + state = {} + pp.set_state(d, state) + + for _, d in self.gb.edges(): + mg = d["mortar_grid"] + + 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) + + def extract_iterate(self, assembler, solution_vector): + """ + Extract parts of the solution for current iterate. + + The iterate solutions in d[pp.STATE]["previous_iterate"] are updated for the + mortar displacements and contact traction are updated. + Method is a tailored copy from assembler.distribute_variable. + + Parameters: + assembler (pp.Assembler): assembler for self.gb. + solution_vector (np.array): solution vector for the current iterate. + + Returns: + (np.array): displacement solution vector for the Nd grid. + + """ + variable_names = [] + for pair in assembler.block_dof.keys(): + variable_names.append(pair[1]) + + dof = np.cumsum(np.append(0, np.asarray(assembler.full_dof))) + + 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.mortar_displacement_variable + ] = mortar_u.copy() + else: + 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): + """ + Reconstruct the displacement jump in local coordinates. + + Args: + data_edge (dictionary): The dictionary on the gb edge. Should contain + - a mortar grid + - a projection, obtained by calling + pp.contact_conditions.set_projections(self.gb) + Returns: + (np.array): ambient_dim x g_l.num_cells. First 1-2 dimensions are in the + tangential direction of the fracture, last dimension is normal. + """ + mg = data_edge["mortar_grid"] + mortar_u = data_edge[pp.STATE][self.mortar_displacement_variable] + displacement_jump_global_coord = ( + mg.mortar_to_slave_avg(nd=self.Nd) + * mg.sign_of_mortar_sides(nd=self.Nd) + * mortar_u + ) + projection = data_edge["tangential_normal_projection"] + # Rotated displacement jumps. these are in the local coordinates, on + project_to_local = projection.project_tangential_normal(int(mg.num_cells / 2)) + u_mortar_local = project_to_local * displacement_jump_global_coord + return u_mortar_local.reshape((self.Nd, -1), order="F") + + def _set_friction_coefficient(self, g): + + nodes = g.nodes + + tips = nodes[:, [0, -1]] + + fc = g.cell_centers + D = cdist(fc.T, tips.T) + D = np.min(D, axis=1) + R = 200 + beta = 10 + friction_coefficient = 0.5 * (1 + beta * np.exp(-R * D ** 2)) + # friction_coefficient = 0.5 * np.ones(g.num_cells) + return friction_coefficient + + +def run_mechanics(setup): + """ + Function for solving linear elasticity with a non-linear Coulomb contact. + + In addition to the standard parameters for mpsa 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: + create_grid(): Create and return the grid bucket + set_parameters(): assigns data to grid bucket. + assign_variables(): assigns variables on grid bucket nodes and edges. + assign_discretizations(): assigns discretizations on grid bucket nodes + and edges. + initial_condition(): Returns initial guess for 'u' and 'lam'. + and attributes: + folder_name: returns a string. The data from the simulation will be + written to the file 'folder_name/' + setup.out_name and the vtk files to + 'res_plot/' + setup.out_name + """ + # Define mixed-dimensional grid. Avoid overwriting existing gb. + if "gb" in setup.__dict__: + gb = setup.gb + else: + gb = setup.create_grid() + gb = setup.gb + + # Pick up grid of highest dimension - there should be a single one of these + g_max = gb.grids_of_dimension(setup.Nd)[0] + # Set simulation parameters and assign variables and discretizations + setup.set_parameters() + setup.initial_condition() + setup.assign_variables() + setup.assign_discretizations() + + # Set up assembler and discretize + assembler = pp.Assembler(gb) + assembler.discretize() + + # Prepare for iteration + + u0 = gb.node_props(g_max)[pp.STATE][setup.displacement_variable] + errors = [] + + counter_newton = 0 + converged_newton = False + max_newton = 15 + + viz = pp.Exporter(g_max, name="mechanics", folder=setup.folder_name) + + 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 + viz.write_vtk({"ux": u0[::2], "uy": u0[1::2]}) + errors.append(error) + + if counter_newton > max_newton and not converged_newton: + raise ValueError("Newton iterations did not converge") + assembler.distribute_variable(sol) + + +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] + + # Re-discretize the nonlinear term + assembler.discretize(term_filter=setup.friction_coupling_term) + + # 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) + # 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 < tol and iterate_difference < tol: + converged = True + error = np.sum((u1 - u0) ** 2) + + else: + if iterate_difference / solution_norm < tol: + converged = True + error = np.sum((u1 - u0) ** 2) / np.sum(u1 ** 2) + + print( + "Error, solution norm and iterate_difference/solution_norm: ", + error, + solution_norm, + iterate_difference / solution_norm, + ) + + return sol, u1, error, converged + + +def l2_norm_cell(g, u, uref=None): + """ + Compute the cell volume weighted norm of a vector-valued cellwise quantity. + + Args: + g (pp.Grid) + u (np.array): Vector-valued function. + """ + 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) diff --git a/src/porepy/numerics/__init__.py b/src/porepy/numerics/__init__.py index 934034839f..050c592e26 100644 --- a/src/porepy/numerics/__init__.py +++ b/src/porepy/numerics/__init__.py @@ -1,2 +1 @@ from . import interface_laws, mixed_dim -from . import fracture_deformation diff --git a/src/porepy/numerics/contact_mechanics/__init__.py b/src/porepy/numerics/contact_mechanics/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/src/porepy/numerics/contact_mechanics/contact_conditions.py b/src/porepy/numerics/contact_mechanics/contact_conditions.py new file mode 100644 index 0000000000..dc9e89fc93 --- /dev/null +++ b/src/porepy/numerics/contact_mechanics/contact_conditions.py @@ -0,0 +1,464 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon May 13 08:53:05 2019 + +@author: eke001 + +For details on the conditions discretized herein, see + +Berge et al., 2019: Finite volume discretization for poroelastic media with fractures +modeled by contact mechanics. +""" +import numpy as np +import scipy.sparse as sps + +import porepy as pp + + +class ColoumbContact: + def __init__(self, keyword, ambient_dimension): + self.keyword = keyword + + self.dim = ambient_dimension + + self.mortar_displacement_variable = "mortar_u" + self.contact_variable = "contact_traction" + + self.traction_discretization = "traction_discretization" + self.displacement_discretization = "displacement_discretization" + self.rhs_discretization = "contact_rhs" + + def _key(self): + return self.keyword + "_" + + def _discretization_key(self): + return self._key() + pp.keywords.DISCRETIZATION + + def ndof(self, g): + return g.num_cells * self.dim + + def discretize(self, g_h, g_l, data_h, data_l, data_edge): + """ Discretize the contact conditions using a semi-smooth Newton + approach. + + The function relates the contact forces, represented on the + lower-dimensional grid, to the jump in displacement between the two + adjacent mortar grids. The function provides a (linearized) + disrcetizaiton of the contact conditions, as described in Berge et al. + + The discertization is stated in the coordinate system defined by the + projection operator associated with the surface. The contact forces + should be interpreted as tangential and normal to this plane. + + NOTES: + Quantities stated in the global coordinate system (e.g. + displacements on the adjacent mortar grids) must be projected to the + local system, using the same projection operator, when paired with the + produced discretization (that is, in the global assembly). + There is a numerical parameter c_num. The sensitivity is currently + unknown. + + Assumptions and other noteworthy aspects: TODO: Rewrite this when the + implementation is ready. + * The contact surface is planar, so that all cells on the surface can + be described by a single normal vector. + * The contact forces are represented directly in the local + coordinate system of the surface. The first self.dim - 1 elements + of the contact vector are the tangential components of the first + cell, then the normal component, then tangential of the second cell + etc. + + """ + + # CLARIFICATIONS NEEDED: + # 1) Do projection and rotation commute on non-matching grids? The + # gut feel says yes, but I'm not sure. + + # Process input + parameters_l = data_l[pp.PARAMETERS] + 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) + + # Numerical parameter, value and sensitivity is currently unknown. + # The thesis of Huber is probably a good place to look for information. + c_num = 100 + + mg = data_edge["mortar_grid"] + + # TODO: Implement a single method to get the normal vector with right sign + # thus the right local coordinate system. + + # Pick the projection operator (defined elsewhere) for this surface. + # IMPLEMENATION NOTE: It is paramount that this projection is used for all + # operations relating to this surface, or else directions of normal vectors + # will get confused. + projection = data_edge["tangential_normal_projection"] + + # The contact force is already computed in local coordinates + contact_force = data_l[pp.STATE]["previous_iterate"][self.contact_variable] + + # Pick out the tangential and normal direction of the contact force. + # The contact force of the first cell is in the first self.dim elements + # of the vector, second cell has the next self.dim etc. + # By design the tangential force is the first self.dim-1 components of + # each cell, while the normal force is the last component. + normal_indices = np.arange(self.dim - 1, contact_force.size, self.dim) + tangential_indices = np.setdiff1d(np.arange(contact_force.size), normal_indices) + contact_force_normal = contact_force[normal_indices] + contact_force_tangential = contact_force[tangential_indices].reshape( + (self.dim - 1, g_l.num_cells), order="F" + ) + + # The displacement jump (in global coordinates) is found by switching the + # sign of the second mortar grid, and then sum the displacements on the + # two sides (which is really a difference since one of the sides have + # its sign switched). + 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.mortar_displacement_variable] + ) + # Rotated displacement jumps. These are in the local coordinates, on + # the lower-dimensional grid + displacement_jump_normal = ( + projection.project_normal(g_l.num_cells) * displacement_jump_global_coord + ) + # The jump in the tangential direction is in g_l.dim columns, one per + # dimension in the tangential direction. + displacement_jump_tangential = ( + projection.project_tangential(g_l.num_cells) + * displacement_jump_global_coord + ).reshape((self.dim - 1, g_l.num_cells), order="F") + + # The friction bound is computed from the previous state of the contact + # force and normal component of the displacement jump. + # Note that the displacement jump is rotated before adding to the contact force + friction_bound = friction_coefficient * np.clip( + -contact_force_normal + c_num * displacement_jump_normal, 0, np.inf + ) + + num_cells = friction_coefficient.size + + # Find contact and sliding region + + # Contact region is determined from the normal direction. + penetration_bc = self._penetration( + contact_force_normal, displacement_jump_normal, c_num + ) + sliding_bc = self._sliding( + contact_force_tangential, + displacement_jump_tangential, + friction_bound, + c_num, + ) + + # Structures for storing the computed coefficients. + displacement_weight = [] # Multiplies displacement jump + traction_weight = [] # Multiplies the normal forces + rhs = np.array([]) # Goes to the right hand side. + + # Zero vectors of the size of the tangential space and the full space, + # respectively. These are needed to complement the discretization + # coefficients to be determined below. + zer = np.array([0] * (self.dim - 1)) + zer1 = np.array([0] * (self.dim)) + zer1[-1] = 1 + + # Loop over all mortar cells, discretize according to the current state of + # the contact + # The loop computes three parameters: + # L will eventually multiply the displacement jump, and be associated with + # the coefficient in a Robin boundary condition (using the terminology of + # the mpsa implementation) + # r is the right hand side term + # IS: Comment about the traction weight? + + for i in range(num_cells): + if sliding_bc[i] & penetration_bc[i]: # in contact and sliding + # The equation for the normal direction is computed from equation + # (24)-(25) in Berge et al. + # Compute coeffecients L, r, v + loc_displacement_tangential, r, v = self._L_r( + contact_force_tangential[:, i], + displacement_jump_tangential[:, i], + friction_bound[i], + c_num, + ) + + # There is no interaction between displacement jumps in normal and + # tangential direction + L = np.hstack((loc_displacement_tangential, np.atleast_2d(zer).T)) + loc_displacement_weight = np.vstack((L, zer1)) + # Right hand side is computed from (24-25). In the normal + # direction, zero displacement is enforced. + # This assumes that the original distance, g, between the fracture + # walls is zero. + r = np.vstack((r + friction_bound[i] * v, 0)) + # Unit contribution from tangential force + loc_traction_weight = np.eye(self.dim) + # Zero weight on normal force + loc_traction_weight[-1, -1] = 0 + # Contribution from normal force + loc_traction_weight[:-1, -1] = -friction_coefficient[i] * v.ravel() + + elif ~sliding_bc[i] & penetration_bc[i]: # In contact and sticking + # Weight for contact force computed according to (23) + loc_traction_tangential = ( + -friction_coefficient[i] + * displacement_jump_tangential[:, i].ravel("F") + / friction_bound[i] + ) + # Unit coefficient for all displacement jumps + loc_displacement_weight = np.eye(self.dim) + + # Tangential traction dependent on normal one + loc_traction_weight = np.zeros((self.dim, self.dim)) + loc_traction_weight[:-1, -1] = loc_traction_tangential + + # The right hand side is the previous tangential jump, and zero + # in the normal direction. + r = np.hstack((displacement_jump_tangential[:, i], 0)).T + + elif ~penetration_bc[i]: # not in contact + # This is a free boundary, no conditions on displacement + loc_displacement_weight = np.zeros((self.dim, self.dim)) + + # Free boundary conditions on the forces. + loc_traction_weight = np.eye(self.dim) + r = np.zeros(self.dim) + + else: # should never happen + raise AssertionError("Should not get here") + + # Depending on the state of the system, the weights in the tangential direction may + # become huge or tiny compared to the other equations. This will + # impede convergence of an iterative solver for the linearized + # system. As a partial remedy, rescale the condition to become + # closer to unity. + w_diag = np.diag(loc_displacement_weight) + np.diag(loc_traction_weight) + W_inv = np.diag(1 / w_diag) + loc_displacement_weight = W_inv.dot(loc_displacement_weight) + loc_traction_weight = W_inv.dot(loc_traction_weight) + r = r.ravel() / w_diag + + # Append to the list of global coefficients. + displacement_weight.append(loc_displacement_weight) + traction_weight.append(loc_traction_weight) + rhs = np.hstack((rhs, r)) + + traction_discretization_coefficients = sps.block_diag(traction_weight) + displacement_discretization_coefficients = sps.block_diag(displacement_weight) + + data_l[pp.DISCRETIZATION_MATRICES][self.keyword][ + self.traction_discretization + ] = traction_discretization_coefficients + data_l[pp.DISCRETIZATION_MATRICES][self.keyword][ + self.displacement_discretization + ] = displacement_discretization_coefficients + data_l[pp.DISCRETIZATION_MATRICES][self.keyword][self.rhs_discretization] = rhs + + def assemble_matrix_rhs(self, g, data): + # Generate matrix for the coupling. This can probably be generalized + # once we have decided on a format for the general variables + traction_coefficient = data[pp.DISCRETIZATION_MATRICES][self.keyword][ + self.traction_discretization + ] + displacement_coefficient = data[pp.DISCRETIZATION_MATRICES][self.keyword][ + self.displacement_discretization + ] + + rhs = data[pp.DISCRETIZATION_MATRICES][self.keyword][self.rhs_discretization] + + return traction_coefficient, displacement_coefficient, rhs + + # Active and inactive boundary faces + def _sliding(self, Tt, ut, bf, ct): + """ Find faces where the frictional bound is exceeded, that is, the face is + sliding. + + Arguments: + Tt (np.array, nd-1 x num_faces): Tangential forces. + u_hat (np.array, nd-1 x num_faces): Displacements in tangential + direction. + bf (np.array, num_faces): Friction bound. + ct (double): Numerical parameter that relates displacement jump to + tangential forces. See Huber et al for explanation. + + Returns: + boolean, size num_faces: True if |-Tt + ct*ut| > bf for a face + + """ + # Use thresholding to not pick up faces that are just about sticking + # Not sure about the sensitivity to the tolerance parameter here. + return self._l2(-Tt + ct * ut) - bf > 1e-10 + + def _penetration(self, Tn, un, cn): + """ Find faces that are in contact. + + Arguments: + Tn (np.array, num_faces): Normal forces. + un (np.array, num_faces): Displament in normal direction. + ct (double): Numerical parameter that relates displacement jump to + normal forces. See Huber et al for explanation. + + Returns: + boolean, size num_faces: True if |-Tt + ct*ut| > bf for a face + + """ + # Not sure about the sensitivity to the tolerance parameter here. + tol = 1e-8 * cn + return (-Tn + cn * un) > tol + + # Below here are different help function for calculating the Newton step + def _ef(self, Tt, cut, bf): + # Compute part of (25) in Berge et al. + return bf / self._l2(-Tt + cut) + + def _Ff(self, Tt, cut, bf): + # Implementation of the term Q involved in the calculation of (25) in Berge + # et al. + numerator = -Tt.dot((-Tt + cut).T) + + # Regularization to avoid issues during the iterations to avoid dividing by + # zero if the faces are not in contact durign iterations. + denominator = max(bf, self._l2(-Tt)) * self._l2(-Tt + cut) + + return numerator / denominator + + def _M(self, Tt, cut, bf): + """ Compute the coefficient M used in Eq. (25) in Berge et al. + """ + Id = np.eye(Tt.shape[0]) + return self._ef(Tt, cut, bf) * (Id - self._Ff(Tt, cut, bf)) + + def _hf(self, Tt, cut, bf): + return self._ef(Tt, cut, bf) * self._Ff(Tt, cut, bf).dot(-Tt + cut) + + def _L_r(self, Tt, ut, bf, c): + """ + Compute the coefficient L, defined in Eq. (25) in Berge et al. + + Arguments: + Tt: Tangential forces. np array, two or three elements + ut: Tangential displacement. Same size as Tt + bf: Friction bound for this mortar cell. + c: Numerical parameter + + + """ + if Tt.ndim <= 1: + Tt = np.atleast_2d(Tt).T + ut = np.atleast_2d(ut).T + + cut = c * ut + # Identity matrix + Id = np.eye(Tt.shape[0]) + + # Shortcut if the friction coefficient is effectively zero. + # Numerical tolerance here is likely somewhat arbitrary. + if bf <= 1e-10: + return ( + 0 * Id, + bf * np.ones((Id.shape[0], 1)), + (-Tt + cut) / self._l2(-Tt + cut), + ) + + # Compute the coefficient M + coeff_M = self._M(Tt, cut, bf) + + # Regularization during the iterations requires computations of parameters + # alpha, beta, delta + alpha = -Tt.T.dot(-Tt + cut) / (self._l2(-Tt) * self._l2(-Tt + cut)) + delta = min(self._l2(-Tt) / bf, 1) + + if alpha < 0: + beta = 1 / (1 - alpha * delta) + else: + beta = 1 + + # The expression (I - beta * M)^-1 + IdM_inv = np.linalg.inv(Id - beta * coeff_M) + + v = IdM_inv.dot(-Tt + cut) / self._l2(-Tt + cut) + + return c * (IdM_inv - Id), -IdM_inv.dot(self._hf(Tt, cut, bf)), v + + def _l2(self, x): + x = np.atleast_2d(x) + return np.sqrt(np.sum(x ** 2, axis=0)) + + +def set_projections(gb): + """ Define a local coordinate system, and projection matrices, for all + grids of co-dimension 1. + + The function adds one item to the data dictionary of all GridBucket edges + that neighbors a co-dimension 1 grid, defined as: + key: tangential_normal_projection, value: pp.TangentialNormalProjection + provides projection to the surface of the lower-dimensional grid + + Note that grids of co-dimension 2 and higher are ignored in this construction, + as we do not plan to do contact mechanics on these objects. + + It is assumed that the surface is planar. + + """ + # Information on the vector normal to the surface is not available directly + # from the surface grid (it could be constructed from the surface geometry, + # which spans the tangential plane). We instead get the normal vector from + # the adjacent higher dimensional grid. + # We therefore access the grids via the edges of the mixed-dimensional grid. + for e, d_m in gb.edges(): + + mg = d_m["mortar_grid"] + # Only consider edges where the lower-dimensional neighbor is of co-dimension 1 + if not mg.dim == (gb.dim_max() - 1): + continue + + # Neigboring grids + _, g_h = gb.nodes_of_edge(e) + + # Find faces of the higher dimensional grid that coincide with the mortar + # grid. Go via the master to mortar projection + # Convert matrix to csr, then the relevant face indices are found from + # the row indices + faces_on_surface = mg.master_to_mortar_int().tocsr().indices + + # Find out whether the boundary faces have outwards pointing normal vectors + # Negative sign implies that the normal vector points inwards. + sgn = g_h.sign_of_faces(faces_on_surface) + + # Unit normal vector + unit_normal = g_h.face_normals[: g_h.dim] / g_h.face_areas + # Ensure all normal vectors on the relevant surface points outwards + unit_normal[:, faces_on_surface] *= sgn + + # Now we need to pick out *one* normal vector of the higher dimensional grid + # which coincides with this mortar grid. This could probably have been + # done with face tags, but we instead project the normal vectors onto the + # mortar grid to kill off all irrelevant faces. Restriction to a single + # normal vector is done in the construction of the projection object + # (below). + # NOTE: Use a single normal vector to span the tangential and normal space, + # thus assuming the surface is planar. + outwards_unit_vector_mortar = mg.master_to_mortar_int().dot(unit_normal.T).T + + # NOTE: The normal vector is based on the first cell in the mortar grid, + # and will be pointing from that cell towards the other side of the + # mortar grid. This defines the positive direction in the normal direction. + # Although a simpler implementation seems to be possible, going via the + # first element in faces_on_surface, there is no guarantee that this will + # give us a face on the positive (or negative) side, hence the more general + # approach is preferred. + # + # NOTE: The basis for the tangential direction is determined by the + # construction internally in TangentialNormalProjection. + projection = pp.TangentialNormalProjection( + outwards_unit_vector_mortar[:, 0].reshape((-1, 1)) + ) + + # Store the projection operator in the mortar data + d_m["tangential_normal_projection"] = projection diff --git a/src/porepy/numerics/discretization.py b/src/porepy/numerics/discretization.py new file mode 100644 index 0000000000..211470f6e0 --- /dev/null +++ b/src/porepy/numerics/discretization.py @@ -0,0 +1,96 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" Module with a do-nothing discretization class. + + +""" + +import numpy as np +import scipy.sparse as sps + + +class VoidDiscretization: + """ Do-nothing discretization object. Used if a discretizaiton object + is needed for technical reasons, but not really necessary. + + Attributes: + keyword (str): Keyword used to identify parameters and discretization + matrices for this object. + ndof_cell (int): Number of degrees of freedom per cell in a grid. + ndof_face (int): Number of degrees of freedom per face in a grid. + ndof_node (int): Number of degrees of freedom per node in a grid. + + """ + + def __init__(self, keyword, ndof_cell=0, ndof_face=0, ndof_node=0): + """ Set the discretization, with the keyword used for storing various + information associated with the discretization. + + Paramemeters: + keyword (str): Identifier of all information used for this + discretization. + ndof_cell (int, optional): Number of degrees of freedom per cell + in a grid. Defaults to 0. + ndof_face (int, optional): Number of degrees of freedom per face + in a grid. Defaults to 0. + ndof_node (int, optional): Number of degrees of freedom per node + in a grid. Defaults to 0. + + """ + self.keyword = keyword + self.ndof_cell = ndof_cell + self.ndof_face = ndof_face + self.ndof_node = ndof_node + + def _key(self): + """ Get the keyword of this object, on a format friendly to access relevant + fields in the data dictionary + + Returns: + String, on the form self.keyword + '_'. + + """ + return self.keyword + "_" + + def ndof(self, g): + """ Abstract method. Return the number of degrees of freedom associated to the + method. + + Parameters + g (grid): Computational grid + + Returns: + int: the number of degrees of freedom. + + """ + return ( + g.num_cells * self.ndof_cell + + g.num_faces * self.ndof_face + + g.num_nodes * self.ndof_node + ) + + def discretize(self, g, data): + """ Construct discretization matrices. + + Parameters: + g (pp.Grid): Grid to be discretized. + data (dictionary): With discretization parameters. + + """ + pass + + def assemble_matrix_rhs(self, g, data): + """ Assemble discretization matrix and rhs vector, both empty. + + Parameters: + g (pp.Grid): Grid to be discretized. + data (dictionary): With discretization parameters. + + Returns: + sps.csc_matrix: Of specified dimensions relative to the grid. Empty. + np.array: Of specified dimensions relative to the grid. All zeros. + + """ + ndof = self.ndof(g) + + return sps.csc_matrix((ndof, ndof)), np.zeros(ndof) diff --git a/src/porepy/numerics/fracture_deformation.py b/src/porepy/numerics/fracture_deformation.py deleted file mode 100644 index 10dd5636fb..0000000000 --- a/src/porepy/numerics/fracture_deformation.py +++ /dev/null @@ -1,30 +0,0 @@ -""" This is a resurrection of a small part of the old fracture_deformation -model. It contains a single function that should really be placed somewhere els. - - -""" - -import numpy as np -import scipy.sparse as sps - - -def sign_of_faces(g, faces): - """ - returns the sign of faces as defined by g.cell_faces. - Parameters: - g: (Grid Object) - faces: (ndarray) indices of faces that you want to know the sign for. The - faces must be boundary faces. - Returns: - sgn: (ndarray) the sign of the faces - """ - - IA = np.argsort(faces) - IC = np.argsort(IA) - - fi, _, sgn = sps.find(g.cell_faces[faces[IA], :]) - assert fi.size == faces.size, "sign of internal faces does not make sense" - I = np.argsort(fi) - sgn = sgn[I] - sgn = sgn[IC] - return sgn diff --git a/src/porepy/numerics/fv/__init__.py b/src/porepy/numerics/fv/__init__.py index 16e3323bc4..e7fec0eb52 100644 --- a/src/porepy/numerics/fv/__init__.py +++ b/src/porepy/numerics/fv/__init__.py @@ -3,7 +3,7 @@ from .mpfa import Mpfa -from .mpsa import FracturedMpsa, Mpsa +from .mpsa import Mpsa from .source import ScalarSource diff --git a/src/porepy/numerics/fv/biot.py b/src/porepy/numerics/fv/biot.py index 85f9fb9048..556cefb514 100644 --- a/src/porepy/numerics/fv/biot.py +++ b/src/porepy/numerics/fv/biot.py @@ -1,6 +1,7 @@ import scipy.sparse as sps import scipy.sparse.linalg as la import numpy as np +import warnings import porepy as pp @@ -22,7 +23,8 @@ def __init__( """ self.mechanics_keyword = mechanics_keyword self.flow_keyword = flow_keyword - # Set variable names for the vector and scalar variable + # Set variable names for the vector and scalar variable, used to access + # solutions from previous time steps self.vector_variable = vector_variable self.scalar_variable = scalar_variable @@ -62,7 +64,7 @@ def rhs_bound(self, g, data): TODO: Boundary effects of coupling terms. - There is an assumption on constant mechanics BCs, see DivD.assemble_matrix(). + There is an assumption on constant mechanics BCs, see DivU.assemble_matrix(). Parameters: g: grid, or subclass, with geometry fields computed. @@ -96,10 +98,10 @@ def rhs_bound(self, g, data): p_bound = -div_flow * bound_flux * p * dt s_bound = -div_mech * bound_stress * d # Note that the following is zero only if the previous time step is zero. - # See comment in the DivD class + # See comment in the DivU class biot_alpha = data[pp.PARAMETERS][self.flow_keyword]["biot_alpha"] - div_d_rhs = -0 * biot_alpha * matrices_m["bound_div_d"] * d - return np.hstack((s_bound, p_bound + div_d_rhs)) + div_u_rhs = -0 * biot_alpha * matrices_f["bound_div_u"] * d + return np.hstack((s_bound, p_bound + div_u_rhs)) def rhs_time(self, g, data): """ Time component of the right hand side (dependency on previous time @@ -126,18 +128,15 @@ def rhs_time(self, g, data): self.scalar_variable: np.zeros(g.num_cells), } - d = self.extractD(g, state[self.vector_variable], as_vector=True) + d = self.extract_vector(g, state[self.vector_variable], as_vector=True) p = state[self.scalar_variable] parameter_dictionary = data[pp.PARAMETERS][self.mechanics_keyword] matrix_dictionaries = data[pp.DISCRETIZATION_MATRICES] - d_scaling = parameter_dictionary.get("displacement_scaling", 1) - div_d = matrix_dictionaries[self.mechanics_keyword]["div_d"] + div_u = matrix_dictionaries[self.flow_keyword]["div_u"] - div_d_rhs = np.squeeze( - parameter_dictionary["biot_alpha"] * div_d * d * d_scaling - ) + div_u_rhs = np.squeeze(parameter_dictionary["biot_alpha"] * div_u * d) p_cmpr = matrix_dictionaries[self.flow_keyword]["mass"] * p mech_rhs = np.zeros(g.dim * g.num_cells) @@ -147,7 +146,7 @@ def rhs_time(self, g, data): # discretization. stab_time = matrix_dictionaries[self.flow_keyword]["biot_stabilization"] * p - return np.hstack((mech_rhs, div_d_rhs + p_cmpr + stab_time)) + return np.hstack((mech_rhs, div_u_rhs + p_cmpr + stab_time)) def discretize(self, g, data): """ Discretize flow and mechanics equations using FV methods. @@ -239,13 +238,12 @@ def assemble_matrix(self, g, data): # Time step size dt = param[self.flow_keyword]["time_step"] - d_scaling = param[self.mechanics_keyword].get("displacement_scaling", 1) # Matrix for left hand side A_biot = sps.bmat( [ [A_mech, grad_p], [ - matrices_m["div_d"] * biot_alpha * d_scaling, + matrices_f["div_u"] * biot_alpha, matrices_f["mass"] + dt * A_flow + stabilization, ], ] @@ -317,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)) @@ -356,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 ) @@ -394,11 +390,11 @@ def _discretize_mech(self, g, data): # trace of strain matrix div = self._subcell_gradient_to_cell_scalar(g, cell_node_blocks) - div_d = div * igrad * rhs_cells + div_u = div * igrad * rhs_cells - # The boundary discretization of the div_d term is represented directly + # The boundary discretization of the div_u term is represented directly # on the cells, instead of going via the faces. - bound_div_d = div * igrad * rhs_bound + bound_div_u = div * igrad * rhs_bound # Call discretization of grad_p-term rhs_jumps, grad_p_face = self.discretize_biot_grad_p( @@ -428,8 +424,8 @@ def _discretize_mech(self, g, data): # Add discretizations to data matrices_m["stress"] = stress matrices_m["bound_stress"] = bound_stress - matrices_m["div_d"] = div_d - matrices_m["bound_div_d"] = bound_div_d + matrices_f["div_u"] = div_u + matrices_f["bound_div_u"] = bound_div_u matrices_m["grad_p"] = grad_p matrices_f["biot_stabilization"] = stabilization matrices_m["bound_displacement_cell"] = disp_cell @@ -464,19 +460,22 @@ def discretize_biot_grad_p(self, g, subcell_topology, alpha, bound_exclusion): imbalance, and thus induce additional displacement gradients in the sub-cells. An additional system is set up, which applies non-zero conditions to the traction continuity equation. This can be expressed as a linear system on the form + (i) A * grad_u = I (ii) B * grad_u + C * u_cc = 0 (iii) 0 D * u_cc = 0 + Thus (i)-(iii) can be inverted to express the additional displacement gradients due to imbalance in pressure forces as in terms of the cell center variables. Thus we can compute the basis functions 'grad_p_jumps' on the sub-cells. To ensure traction continuity, as soon as a convention is chosen for what side the force evaluation should be considered on, an additional term, called - 'grad_p_face', is added to the full force. This latter term represnts the force + 'grad_p_face', is added to the full force. This latter term represents the force due to cell-center pressure acting on the face from the chosen side. The pair subfno_unique-unique_subfno gives the side convention. The full force on the face is therefore given by - t = stress * u + bound_stress * u_b + (grad_p_jumps + grad_p_face) * p + + t = stress * u + bound_stress * u_b + alpha * (grad_p_jumps + grad_p_face) * p The strategy is as follows. 1. compute product normal_vector * alpha and get a map for vector problems @@ -489,12 +488,11 @@ def discretize_biot_grad_p(self, g, subcell_topology, alpha, bound_exclusion): num_subhfno = subcell_topology.subhfno.size num_subfno_unique = subcell_topology.num_subfno_unique num_subfno = subcell_topology.num_subfno - num_cno = subcell_topology.num_cno - - num_nodes = np.diff(g.face_nodes.indptr) # Step 1 + # The implementation is valid for tensor Biot coefficients, but for the + # moment, we only allow for scalar inputs. # Take Biot's alpha as a tensor alpha_tensor = pp.SecondOrderTensor(nd, alpha * np.ones(g.num_cells)) @@ -507,7 +505,8 @@ def discretize_biot_grad_p(self, g, subcell_topology, alpha, bound_exclusion): nAlpha_grad, cell_node_blocks, sub_cell_index = fvutils.scalar_tensor_vector_prod( g, alpha_tensor, subcell_topology ) - # transfer nAlpha to a face-based + # transfer nAlpha to a subface-based quantity by pairing expressions on the + # two sides of the subface unique_nAlpha_grad = subcell_topology.pair_over_subfaces(nAlpha_grad) # convenience method for reshaping nAlpha from face-based @@ -532,6 +531,7 @@ def map_tensor(mat, nd, ind): # as a force on the faces. The right hand side is thus formed of the # unit vector. def build_rhs_units_single_dimension(dim): + # EK: Can we skip argument dim? vals = np.ones(num_subfno_unique) ind = subcell_topology.subfno_unique mat = sps.coo_matrix( @@ -553,7 +553,7 @@ def build_rhs_units_single_dimension(dim): ].A.ravel("F") # NOTE: For some reason one should not multiply with the sign, but I don't # understand why. It should not matter much for the Biot alpha term since - # by construction the biot_alpha_jumps and biot_alpha_force will cancell for + # by construction the biot_alpha_jumps and biot_alpha_force will cancel for # Neumann boundaries. We keep the sign matrix as an Identity matrix to remember # where it should be multiplied: sgn_nd = np.tile(np.abs(sgn), (g.dim, 1)) @@ -567,7 +567,7 @@ def build_rhs_units_single_dimension(dim): sgn_diag_F = sps.diags(sgn_nd.ravel("F")) sgn_diag_C = sps.diags(sgn_nd.ravel("C")) - # Remembering the ordering of the local equations: + # Recall the ordering of the local equations: # First stress equilibrium for the internal subfaces. # Then the stress equilibrium for the Neumann subfaces. # Then the Robin subfaces. @@ -697,7 +697,7 @@ def slv(b): return slv # ----------------------- Methods for post processing ------------------------- - def extractD(self, g, u, dims=None, as_vector=False): + def extract_vector(self, g, u, dims=None, as_vector=False): """ Extract displacement field from solution. Parameters: @@ -725,7 +725,7 @@ def extractD(self, g, u, dims=None, as_vector=False): else: return vals - def extractP(self, g, u): + def extract_scalar(self, g, u): """ Extract pressure field from solution. Parameters: @@ -757,7 +757,7 @@ def compute_flux(self, g, u, data): flux_discr = data[pp.DISCRETIZATION_MATRICES][self.flow_keyword]["flux"] bound_flux = data[pp.DISCRETIZATION_MATRICES][self.flow_keyword]["bound_flux"] bound_val = data[pp.PARAMETERS][self.flow_keyword]["bc_values"] - p = self.extractP(g, u) + p = self.extract_scalar(g, u) flux = flux_discr * p + bound_flux * bound_val return flux @@ -781,17 +781,35 @@ def compute_stress(self, g, u, data): stress_discr = matrix_dictionary["stress"] bound_stress = matrix_dictionary["bound_stress"] bound_val = data[pp.PARAMETERS][self.mechanics_keyword]["bc_values"] - d = self.extractD(g, u, as_vector=True) + d = self.extract_vector(g, u, as_vector=True) stress = np.squeeze(stress_discr * d) + (bound_stress * bound_val) return stress -class GradP( - pp.numerics.interface_laws.elliptic_discretization.VectorEllipticDiscretization -): +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. + + Paramemeters: + keyword (str): Identifier of all information used for this + discretization. + """ + self.keyword = keyword + + def _key(self): + """ Get the keyword of this object, on a format friendly to access relevant + fields in the data dictionary + + Returns: + String, on the form self.keyword + '_'. + + """ + return self.keyword + "_" + def ndof(self, g): """ Return the number of degrees of freedom associated to the method. @@ -910,7 +928,7 @@ def assemble_int_bound_displacement_trace( displacement on internal boundaries. The intended use is when the internal boundary is coupled to another - node in the GridBucket sence. Specific usage depends on the + node in the GridBucket sense. Specific usage depends on the interface condition between the nodes; this method will typically be used to impose displacement continuity on an interface. @@ -958,7 +976,7 @@ def assemble_int_bound_displacement_trace( hf2f = pp.fvutils.map_hf_2_f(g=g) num_nodes = np.diff(g.face_nodes.indptr) weight = sps.kron(sps.eye(g.dim), sps.diags(1 / num_nodes)) - # hf2f adds all subface values to one face values. For the displacement we want + # hf2f adds all subface values to one face value. For the displacement we want # to take the average, therefore we divide each face by the number of subfaces. cc[2, self_ind] += proj_avg * weight * hf2f * bp else: @@ -968,22 +986,30 @@ def enforce_neumann_int_bound(self, *_): pass -class DivD( - pp.numerics.interface_laws.elliptic_discretization.VectorEllipticDiscretization -): +class DivU: """ Class for the displacement divergence term of the Biot equation. """ def __init__( - self, keyword="mechanics", variable="displacement", mortar_variable="traction" + self, + mechanics_keyword="mechanics", + flow_keyword="flow", + variable="displacement", + mortar_variable="mortar_displacement", ): - """ Set the two keywords. + """ Set the mechanics keyword and specify the variables. The keywords are used to access and store parameters and discretization matrices. + The variable names are used to obtain the previous solution for the time + discretization. Consequently, they are those of the unknowns contributing to + the DivU term (displacements), not the scalar variable. """ - super().__init__(keyword) - # Set variable name for the vector variable (displacement) + self.flow_keyword = flow_keyword + self.mechanics_keyword = mechanics_keyword + # We also need to specify the names of the displacement variables on the node + # and adjacent edges. T + # Set variable name for the vector variable (displacement). self.variable = variable # The following is only used for mixed-dimensional problems. # Set the variable used for contact mechanics. @@ -1071,14 +1097,14 @@ def assemble_matrix(self, g, data): ValueError if the displacement divergence term has not already been discretized. """ - matrix_dictionary = data[pp.DISCRETIZATION_MATRICES][self.keyword] - if not "div_d" in matrix_dictionary: + matrix_dictionary = data[pp.DISCRETIZATION_MATRICES][self.flow_keyword] + if not "div_u" in matrix_dictionary: raise ValueError( - """DivD class requires a pre-computed discretization to be + """DivU class requires a pre-computed discretization to be stored in the matrix dictionary.""" ) - biot_alpha = data[pp.PARAMETERS][self.keyword]["biot_alpha"] - return matrix_dictionary["div_d"] * biot_alpha + biot_alpha = data[pp.PARAMETERS][self.flow_keyword]["biot_alpha"] + return matrix_dictionary["div_u"] * biot_alpha def assemble_rhs(self, g, data): """ Return the right-hand side for a discretization of the displacement @@ -1095,39 +1121,42 @@ def assemble_rhs(self, g, data): np.ndarray: Zero right hand side vector with representation of boundary conditions. """ - parameter_dictionary = data[pp.PARAMETERS][self.keyword] - matrix_dictionary = data[pp.DISCRETIZATION_MATRICES][self.keyword] + parameter_dictionary_mech = data[pp.PARAMETERS][self.mechanics_keyword] + parameter_dictionary_flow = data[pp.PARAMETERS][self.flow_keyword] + matrix_dictionary = data[pp.DISCRETIZATION_MATRICES][self.flow_keyword] # For IE and constant BCs, the boundary part cancels, as the contribution from # successive timesteps (n and n+1) appear on the rhs with opposite signs. For # transient BCs, use the below with the appropriate version of d_bound_i. - d_bound_1 = parameter_dictionary["bc_values"] + # Get bc values from mechanics + d_bound_1 = parameter_dictionary_mech["bc_values"] - d_bound_0 = data[pp.STATE][self.keyword]["bc_values"] - biot_alpha = parameter_dictionary["biot_alpha"] + d_bound_0 = data[pp.STATE][self.mechanics_keyword]["bc_values"] + # and coupling parameter from flow + biot_alpha = parameter_dictionary_flow["biot_alpha"] rhs_bound = ( - -matrix_dictionary["bound_div_d"] * (d_bound_1 - d_bound_0) * biot_alpha + -matrix_dictionary["bound_div_u"] * (d_bound_1 - d_bound_0) * biot_alpha ) # Time part d_cell = data[pp.STATE][self.variable] - d_scaling = parameter_dictionary.get("displacement_scaling", 1) - div_d = matrix_dictionary["div_d"] - rhs_time = np.squeeze(biot_alpha * div_d * d_cell * d_scaling) + div_u = matrix_dictionary["div_u"] + rhs_time = np.squeeze(biot_alpha * div_u * d_cell) return rhs_bound + rhs_time - def assemble_int_bound_stress( + def assemble_int_bound_displacement_trace( self, g, data, data_edge, grid_swap, cc, matrix, rhs, self_ind ): - """Assemble the contribution the stress mortar on an internal boundary, - manifested as a stress boundary condition. + """Assemble the contribution from the displacement mortar on an internal boundary, + manifested as a displacement boundary condition. The intended use is when the internal boundary is coupled to another node by an interface law. Specific usage depends on the interface condition between the nodes; this method will typically be - used to impose the effect of the stress mortar on the divergence term. + used to impose the effect of the displacement mortar on the divergence term on + the higher-dimensional grid. Implementations of this method will use an interplay between the grid on the node and the mortar grid on the relevant edge. @@ -1154,29 +1183,101 @@ def assemble_int_bound_stress( mg = data_edge["mortar_grid"] if grid_swap: - proj = mg.slave_to_mortar_avg() + proj = mg.mortar_to_slave_avg(nd=g.dim) else: - proj = mg.master_to_mortar_avg() - # Expand indices as Fortran. - proj_int = sps.kron(proj, sps.eye(g.dim)).tocsr() + proj = mg.mortar_to_master_avg(nd=g.dim) - matrix_dictionary = data[pp.DISCRETIZATION_MATRICES][self.keyword] - biot_alpha = data[pp.PARAMETERS][self.keyword]["biot_alpha"] - bound_div_d = matrix_dictionary["bound_div_d"] + matrix_dictionary = data[pp.DISCRETIZATION_MATRICES][self.flow_keyword] + biot_alpha = data[pp.PARAMETERS][self.flow_keyword]["biot_alpha"] + bound_div_u = matrix_dictionary["bound_div_u"] - lam_k = data_edge[pp.STATE][self.mortar_variable] + u_bound_previous = data_edge[pp.STATE][self.mortar_variable] - if bound_div_d.shape[1] != proj_int.shape[1]: + if bound_div_u.shape[1] != proj.shape[0]: raise ValueError( """Inconsistent shapes. Did you define a sub-face boundary condition but only a face-wise mortar?""" ) - # The mortar will act as a boundary condition for the div_d term. - # We assume implicit Euler in Biot, thus the div_d term appares - # on the rhs as div_d^{k-1}. This results in a contribution to the + # The mortar will act as a boundary condition for the div_u term. + # We assume implicit Euler in Biot, thus the div_u term appears + # on the rhs as div_u^{k-1}. This results in a contribution to the + # rhs for the coupling variable also. + cc[self_ind, 2] += biot_alpha * bound_div_u * proj + rhs[self_ind] += biot_alpha * bound_div_u * proj * u_bound_previous + + def assemble_int_bound_displacement_source( + self, g, data, data_edge, cc, matrix, rhs, self_ind + ): + """Assemble the contribution from the displacement mortar on an internal boundary, + manifested as a source term. Only the normal component of the mortar displacement + is considered. + + The intended use is when the internal boundary is coupled to another + node by an interface law. Specific usage depends on the + interface condition between the nodes; this method will typically be + used to impose the effect of the displacement mortar on the divergence term on + the lower-dimensional grid. + + Implementations of this method will use an interplay between the grid + on the node and the mortar grid on the relevant edge. + + Parameters: + g (Grid): Grid which the condition should be imposed on. + data (dictionary): Data dictionary for the node in the + mixed-dimensional grid. + data_edge (dictionary): Data dictionary for the edge in the + mixed-dimensional grid. + grid_swap (boolean): If True, the grid g is identified with the @ + slave side of the mortar grid in data_adge. + cc (block matrix, 3x3): Block matrix for the coupling condition. + The first and second rows and columns are identified with the + master and slave side; the third belongs to the edge variable. + The discretization of the relevant term is done in-place in cc. + matrix (block matrix 3x3): Discretization matrix for the edge and + the two adjacent nodes. + self_ind (int): Index in cc and matrix associated with this node. + Should be either 1 or 2. + + """ + + mg = data_edge["mortar_grid"] + + # From the mortar displacements, we want to + # 1) Take the jump between the two mortar sides, + # 2) Project to the slave grid and + # 3) Extract the normal component. + + # Define projections and rotations + nd = g.dim + 1 + proj = mg.mortar_to_slave_avg(nd=nd) + jump_on_slave = proj * mg.sign_of_mortar_sides(nd=nd) + rotation = data_edge["tangential_normal_projection"] + normal_component = rotation.project_normal(g.num_cells) + + biot_alpha = data[pp.PARAMETERS][self.flow_keyword]["biot_alpha"] + + # Project the previous solution to the slave grid + previous_displacement_jump_global_coord = ( + jump_on_slave * data_edge[pp.STATE][self.mortar_variable] + ) + # Rotated displacement jumps. These are in the local coordinates, on + # the lower-dimensional grid + previous_displacement_jump_normal = ( + normal_component * previous_displacement_jump_global_coord + ) + # The same procedure is applied to the unknown displacements, by assembling the + # jump operator, projection and normal component extraction in the coupling matrix. + # Finally, we integrate over the cell volume. + # The jump on the slave is defined to be negative for an open fracture (!), + # hence the negative sign. + vol = sps.dia_matrix((g.cell_volumes, 0), shape=(g.num_cells, g.num_cells)) + cc[self_ind, 2] -= biot_alpha * vol * normal_component * jump_on_slave + + # We assume implicit Euler in Biot, thus the div_u term appears + # on the rhs as div_u^{k-1}. This results in a contribution to the # rhs for the coupling variable also. - cc[self_ind, 2] += biot_alpha * bound_div_d * proj_int.T - rhs[self_ind] += biot_alpha * bound_div_d * proj_int.T * lam_k + # See note above on sign. This term is on the rhs, yielding the opposite sign. + rhs[self_ind] += biot_alpha * vol * previous_displacement_jump_normal def enforce_neumann_int_bound(self, *_): pass @@ -1244,7 +1345,7 @@ def discretize(self, g, data): discretize method of the Biot class. """ raise NotImplementedError( - """No discretize method implemented for the DivD + """No discretize method implemented for the DivU class. See the Biot class.""" ) diff --git a/src/porepy/numerics/fv/mpsa.py b/src/porepy/numerics/fv/mpsa.py index b6638e6ba9..9188a6182c 100644 --- a/src/porepy/numerics/fv/mpsa.py +++ b/src/porepy/numerics/fv/mpsa.py @@ -19,9 +19,27 @@ logger = logging.getLogger(__name__) -class Mpsa( - pp.numerics.interface_laws.elliptic_discretization.VectorEllipticDiscretization -): +class Mpsa: + def __init__(self, keyword): + """ Set the discretization, with the keyword used for storing various + information associated with the discretization. + + Paramemeters: + keyword (str): Identifier of all information used for this + discretization. + """ + self.keyword = keyword + + def _key(self): + """ Get the keyword of this object, on a format friendly to access relevant + fields in the data dictionary + + Returns: + String, on the form self.keyword + '_'. + + """ + return self.keyword + "_" + def ndof(self, g): """ Return the number of degrees of freedom associated to the method. @@ -113,15 +131,32 @@ def discretize(self, g, data): partial = parameter_dictionary.get("partial_update", False) inverter = parameter_dictionary.get("inverter", None) + max_memory = parameter_dictionary.get("max_memory", None) if not partial: - stress, bound_stress, bound_displacement_cell, bound_displacement_face = mpsa( - g, c, bnd, eta=eta, hf_eta=hf_eta, inverter=inverter - ) - matrix_dictionary["stress"] = stress - matrix_dictionary["bound_stress"] = bound_stress - matrix_dictionary["bound_displacement_cell"] = bound_displacement_cell - matrix_dictionary["bound_displacement_face"] = bound_displacement_face + if max_memory is None: + stress, bound_stress, bound_displacement_cell, bound_displacement_face = mpsa( + g, c, bnd, eta=eta, hf_eta=hf_eta, inverter=inverter + ) + matrix_dictionary["stress"] = stress + matrix_dictionary["bound_stress"] = bound_stress + # Should be face_displacement_cell and _face + matrix_dictionary["bound_displacement_cell"] = bound_displacement_cell + matrix_dictionary["bound_displacement_face"] = bound_displacement_face + + else: + stress, bound_stress = mpsa( + g, + c, + bnd, + eta=eta, + hf_eta=hf_eta, + inverter=inverter, + max_memory=max_memory, + ) + matrix_dictionary["stress"] = stress + matrix_dictionary["bound_stress"] = bound_stress + # This option else: raise NotImplementedError( """Partial discretiation for the Mpsa class is not @@ -487,416 +522,6 @@ def enforce_neumann_int_bound( pass -class FracturedMpsa(Mpsa): - """ - Subclass of MPSA for discretizing a fractured domain. Adds DOFs on each - fracture face which describe the fracture deformation. - """ - - def __init__(self, keyword, given_traction=False, **kwargs): - Mpsa.__init__(self, keyword, **kwargs) - if not hasattr(self, "keyword"): - raise AttributeError("Mpsa must assign keyword") - self.given_traction_flag = given_traction - - def ndof(self, g): - """ - Return the number of degrees of freedom associated to the method. - In this case number of cells times dimension (stress dof). - - Parameter - --------- - g: grid, or a subclass. - - Return - ------ - dof: the number of degrees of freedom. - - """ - num_fracs = np.sum(g.tags["fracture_faces"]) - return g.dim * (g.num_cells + num_fracs) - - def assemble_matrix_rhs(self, g, data, discretize=True, **kwargs): - """ - Return the matrix and right-hand side for a discretization of a second - order elliptic equation using a FV method with a multi-point stress - approximation with dofs added on the fracture interfaces. - - Parameters - ---------- - g : grid, or a subclass, with geometry fields computed. - data: dictionary to store the data. For details on necessary keywords, - see method discretize() - discretize (boolean, optional): default True. Whether to discetize - prior to matrix assembly. If False, data should already contain - discretization. - - Return - ------ - matrix: sparse csr (g.dim * g_num_cells + 2 * {#of fracture faces}, - 2 * {#of fracture faces}) - Discretization matrix. - rhs: array (g.dim * g_num_cells + g.dim * num_frac_faces) - Right-hand side which contains the boundary conditions and the scalar - source term. - """ - if discretize: - self.discretize_fractures(g, data, **kwargs) - - parameter_dictionary = data[pp.PARAMETERS][self.keyword] - matrix_dictionary = data[pp.DISCRETIZATION_MATRICES][self.keyword] - stress = matrix_dictionary["stress"] - bound_stress = matrix_dictionary["bound_stress"] - b_e = matrix_dictionary["b_e"] - A_e = matrix_dictionary["A_e"] - - if self.given_traction_flag: - L, b_l = self.given_traction(g, stress, bound_stress) - else: - L, b_l = self.given_slip_distance(g, stress, bound_stress) - - bc_val = parameter_dictionary["bc_values"] - - frac_faces = np_matlib.repmat(g.tags["fracture_faces"], g.dim, 1) - if parameter_dictionary["bc"].bc_type == "scalar": - frac_faces = frac_faces.ravel("F") - elif parameter_dictionary["bc"].bc_type == "vectorial": - bc_val = bc_val.ravel("F") - else: - raise ValueError("Unknown boundary type") - - slip_distance = parameter_dictionary["slip_distance"] - - A = sps.vstack((A_e, L), format="csr") - rhs = np.hstack((b_e * bc_val, b_l * (slip_distance + bc_val))) - - return A, rhs - - def rhs(self, g, data): - """ - Return the matrix and right-hand side for a discretization of a second - order elliptic equation using a FV method with a multi-point stress - approximation with dofs added on the fracture interfaces. - - Parameters - ---------- - g : grid, or a subclass, with geometry fields computed. - data: dictionary to store the data. For details on necessary keywords, - see method discretize() - discretize (boolean, optional): default True. Whether to discetize - prior to matrix assembly. If False, data should already contain - discretization. - - Return - ------ - matrix: sparse csr (g.dim * g_num_cells + 2 * {#of fracture faces}, - 2 * {#of fracture faces}) - Discretization matrix. - rhs: array (g.dim * g_num_cells + g.dim * num_frac_faces) - Right-hand side which contains the boundary conditions and the scalar - source term. - """ - - parameter_dictionary = data[pp.PARAMETERS][self.keyword] - matrix_dictionary = data[pp.DISCRETIZATION_MATRICES][self.keyword] - stress = matrix_dictionary["stress"] - bound_stress = matrix_dictionary["bound_stress"] - b_e = matrix_dictionary["b_e"] - - if self.given_traction_flag: - _, b_l = self.given_traction(g, stress, bound_stress) - else: - _, b_l = self.given_slip_distance(g, stress, bound_stress) - - bc_val = parameter_dictionary["bc_values"] - - frac_faces = np_matlib.repmat(g.tags["fracture_faces"], 3, 1) - if parameter_dictionary["bc"].bc_type == "scalar": - frac_faces = frac_faces.ravel("F") - - elif parameter_dictionary["bc"].bc_type == "vectorial": - bc_val = bc_val.ravel("F") - else: - raise ValueError("Unknown boundary type") - - slip_distance = parameter_dictionary["slip_distance"] - - rhs = np.hstack((b_e * bc_val, b_l * (slip_distance + bc_val))) - - return rhs - - def traction(self, g, data, sol): - """ - Extract the traction on the faces from fractured fv solution. - - Parameters - ---------- - g : grid, or a subclass, with geometry fields computed. - sol : array (g.dim * (g.num_cells + {#of fracture faces})) - Solution, stored as [cell_disp, fracture_disp] - - Return - ------ - T : array (g.dim * g.num_faces) - traction on each face - - """ - parameter_dictionary = data[pp.PARAMETERS][self.keyword] - matrix_dictionary = data[pp.DISCRETIZATION_MATRICES][self.keyword] - bc_val = parameter_dictionary["bc_values"].copy() - frac_disp = self.extract_frac_u(g, sol) - cell_disp = self.extract_u(g, sol) - - frac_faces = (g.frac_pairs).ravel("C") - - if parameter_dictionary["bc"].bc_type == "vectorial": - bc_val = bc_val.ravel("F") - - frac_ind = pp.utils.mcolon.mcolon( - g.dim * frac_faces, g.dim * frac_faces + g.dim - ) - bc_val[frac_ind] = frac_disp - - T = ( - matrix_dictionary["stress"] * cell_disp - + matrix_dictionary["bound_stress"] * bc_val - ) - - return T - - def extract_u(self, g, sol): - """ Extract the cell displacement from fractured fv solution. - - Parameters - ---------- - g : grid, or a subclass, with geometry fields computed. - sol : array (g.dim * (g.num_cells + {#of fracture faces})) - Solution, stored as [cell_disp, fracture_disp] - - Return - ------ - u : array (g.dim * g.num_cells) - displacement at each cell - - """ - # pylint: disable=invalid-name - return sol[: g.dim * g.num_cells] - - def extract_frac_u(self, g, sol): - """ Extract the fracture displacement from fractured fv solution. - - Parameters - ---------- - g : grid, or a subclass, with geometry fields computed. - sol : array (g.dim * (g.num_cells + {#of fracture faces})) - Solution, stored as [cell_disp, fracture_disp] - - Return - ------ - u : array (g.dim *{#of fracture faces}) - displacement at each fracture face - - """ - # pylint: disable=invalid-name - return sol[g.dim * g.num_cells :] - - def discretize_fractures(self, g, data, faces=None, **kwargs): - """ - Discretize the vector elliptic equation by the multi-point stress and added - degrees of freedom on the fracture faces - - The method computes fluxes over faces in terms of displacements in - adjacent cells (defined as the two cells sharing the face). - - The name of data in the input dictionary (data) are: - param : Parameter(Class). Contains the following parameters: - tensor : fourth_order_tensor - Permeability defined cell-wise. If not given a identity permeability - is assumed and a warning arised. - bc : boundary conditions (optional) - bc_val : dictionary (optional) - Values of the boundary conditions. The dictionary has at most the - following keys: 'dir' and 'neu', for Dirichlet and Neumann boundary - conditions, respectively. - apertures : (np.ndarray) (optional) apertures of the cells for scaling of - the face normals. - - Parameters - ---------- - g : grid, or a subclass, with geometry fields computed. - data: dictionary to store the data. - """ - - # dir_bound = g.get_all_boundary_faces() - # bound = bc.BoundaryCondition(g, dir_bound, ['dir'] * dir_bound.size) - parameter_dictionary = data[pp.PARAMETERS][self.keyword] - matrix_dictionary = data[pp.DISCRETIZATION_MATRICES][self.keyword] - - frac_faces = g.tags["fracture_faces"] - - bound = parameter_dictionary["bc"] - - if bound.bc_type == "scalar": - bound.is_dir[frac_faces] = True - bound.is_neu[frac_faces] = False - elif bound.bc_type == "vectorial": - bound.is_dir[:, frac_faces] = True - bound.is_neu[:, frac_faces] = False - else: - raise ValueError("Unknow boundary condition type: " + bound.bc_type) - if np.sum(bound.is_dir * bound.is_neu) != 0: - raise AssertionError("Found faces that are both dirichlet and neuman") - # Discretize with normal mpsa - self.discretize(g, data, **kwargs) - - stress, bound_stress = ( - matrix_dictionary["stress"], - matrix_dictionary["bound_stress"], - ) - - # Create A and rhs - div = pp.fvutils.vector_divergence(g) - a = div * stress - b = div * bound_stress - - # we find the matrix indices of the fracture - if faces is None: - frac_faces = g.frac_pairs - frac_faces_left = frac_faces[0] - frac_faces_right = frac_faces[1] - else: - raise NotImplementedError("not implemented given faces") - - int_b_left = pp.utils.mcolon.mcolon( - g.dim * frac_faces_left, g.dim * frac_faces_left + g.dim - ) - int_b_right = pp.utils.mcolon.mcolon( - g.dim * frac_faces_right, g.dim * frac_faces_right + g.dim - ) - int_b_ind = np.ravel((int_b_left, int_b_right), "C") - - # We find the sign of the left and right faces. - sgn_left = _sign_matrix(g, frac_faces_left) - sgn_right = _sign_matrix(g, frac_faces_right) - # The displacement on the internal boundary face are considered unknowns, - # so we move them over to the lhs. The rhs now only consists of the - # external boundary faces - b_internal = b[:, int_b_ind] - b_external = b.copy() - pp.utils.sparse_mat.zero_columns(b_external, int_b_ind) - - bound_stress_external = bound_stress.copy().tocsc() - pp.utils.sparse_mat.zero_columns(bound_stress_external, int_b_ind) - # We assume that the traction on the left hand side is equal but - # opisite - - frac_stress_diff = ( - sgn_left * bound_stress[int_b_left, :] - + sgn_right * bound_stress[int_b_right, :] - )[:, int_b_ind] - internal_stress = sps.hstack( - ( - sgn_left * stress[int_b_left, :] + sgn_right * stress[int_b_right, :], - frac_stress_diff, - ) - ) - - A = sps.vstack((sps.hstack((a, b_internal)), internal_stress), format="csr") - # negative sign since we have moved b_external from lhs to rhs - d_b = -b_external - # sps.csr_matrix((int_b_left.size, g.num_faces * g.dim)) - d_t = ( - -sgn_left * bound_stress_external[int_b_left] - - sgn_right * bound_stress_external[int_b_right] - ) - - b_matrix = sps.vstack((d_b, d_t), format="csr") - - matrix_dictionary["b_e"] = b_matrix - matrix_dictionary["A_e"] = A - - def given_traction(self, g, stress, bound_stress, faces=None, **kwargs): - # we find the matrix indices of the fracture - if faces is None: - frac_faces = g.frac_pairs - frac_faces_left = frac_faces[0] - frac_faces_right = frac_faces[1] - else: - raise NotImplementedError("not implemented given faces") - - int_b_left = pp.utils.mcolon.mcolon( - g.dim * frac_faces_left, g.dim * frac_faces_left + g.dim - ) - int_b_right = pp.utils.mcolon.mcolon( - g.dim * frac_faces_right, g.dim * frac_faces_right + g.dim - ) - int_b_ind = np.ravel((int_b_left, int_b_right), "C") - - # We find the sign of the left and right faces. - sgn_left = _sign_matrix(g, frac_faces_left) - sgn_right = _sign_matrix(g, frac_faces_right) - - # We obtain the stress from boundary conditions on the domain boundary - bound_stress_external = bound_stress.copy().tocsc() - pp.utils.sparse_mat.zero_columns(bound_stress_external, int_b_ind) - bound_stress_external = bound_stress_external.tocsc() - - # We construct the L matrix, i.e., we set the traction on the left - # fracture side - frac_stress = (sgn_left * bound_stress[int_b_left, :])[:, int_b_ind] - - L = sps.hstack((sgn_left * stress[int_b_left, :], frac_stress)) - - # negative sign since we have moved b_external from lhs to rhs - d_t = ( - sps.csr_matrix( - (np.ones(int_b_left.size), (np.arange(int_b_left.size), int_b_left)), - (int_b_left.size, g.num_faces * g.dim), - ) - - sgn_left * bound_stress_external[int_b_left] - ) # \ - # + sgn_right * bound_stress_external[int_b_right] - - return L, d_t - - def given_slip_distance(self, g, stress, bound_stress, faces=None): - # we find the matrix indices of the fracture - if faces is None: - frac_faces = g.frac_pairs - frac_faces_left = frac_faces[0] - frac_faces_right = frac_faces[1] - else: - raise NotImplementedError("not implemented given faces") - - int_b_left = pp.utils.mcolon.mcolon( - g.dim * frac_faces_left, g.dim * frac_faces_left + g.dim - ) - int_b_right = pp.utils.mcolon.mcolon( - g.dim * frac_faces_right, g.dim * frac_faces_right + g.dim - ) - int_b_ind = np.ravel((int_b_left, int_b_right), "C") - - # We construct the L matrix, by assuming that the relative displacement - # is given - L = sps.hstack( - ( - sps.csr_matrix((int_b_left.size, g.dim * g.num_cells)), - sps.identity(int_b_left.size), - -sps.identity(int_b_right.size), - ) - ) - - d_f = sps.csr_matrix( - (np.ones(int_b_left.size), (np.arange(int_b_left.size), int_b_left)), - (int_b_left.size, g.num_faces * g.dim), - ) - - return L, d_f - - -# ------------------------------------------------------------------------------# - - def mpsa( g, constit, @@ -1024,7 +649,7 @@ def mpsa( # entire grid. # TODO: We may want to estimate the memory need, and give a warning if # this seems excessive - stress, bound_stress, hf_cell, hf_bound = _mpsa_local( + return _mpsa_local( g, constit, bound, @@ -1033,16 +658,20 @@ def mpsa( hf_disp=hf_disp, hf_eta=hf_eta, ) + else: + if hf_disp: + raise ValueError("Mpsa options max_memory and hf_disp are incompatible") + # Estimate number of partitions necessary based on prescribed memory # usage peak_mem = _estimate_peak_memory_mpsa(g) - num_part = np.ceil(peak_mem / max_memory) + num_part = np.ceil(peak_mem / max_memory).astype(np.int) logger.info("Split MPSA discretization into " + str(num_part) + " parts") # Let partitioning module apply the best available method - part = pp.grid.partition.partition_metis(g, num_part) + part = pp.partition.partition_metis(g, num_part) # Empty fields for stress and bound_stress. Will be expanded as we go. # Implementation note: It should be relatively straightforward to @@ -1067,7 +696,13 @@ def mpsa( # Perform local discretization. loc_stress, loc_bound_stress, loc_faces = mpsa_partial( - g, constit, bound, eta=eta, inverter=inverter, nodes=active_nodes + g, + constit, + bound, + eta=eta, + inverter=inverter, + nodes=active_nodes, + hf_disp=False, ) # Eliminate contribution from faces already covered @@ -1080,7 +715,7 @@ def mpsa( stress += loc_stress bound_stress += loc_bound_stress - return stress, bound_stress, hf_cell, hf_bound + return stress, bound_stress def mpsa_update_partial( @@ -1321,42 +956,45 @@ def mpsa_partial( pp.fvutils.zero_out_sparse_rows(stress_glob, eliminate_ind) pp.fvutils.zero_out_sparse_rows(bound_stress_glob, eliminate_ind) - # If we are returning the subface displacement reconstruction matrices we have - # to do some more work. The following is equivalent to what is done for the stresses, - # but as they are working on faces, the displacement reconstruction has to work on - # subfaces. - # First, we find the mappings from local subfaces to global subfaces - subcell_topology = pp.fvutils.SubcellTopology(g) - l2g_sub_faces = np.where(np.in1d(subcell_topology.fno_unique, l2g_faces))[0] - # We now create a fake grid, just to be able to use the function map_subgrid_to_grid. - subgrid = pp.CartGrid([1] * g.dim) - subgrid.num_faces = subcell_topology.fno_unique.size - subgrid.num_cells = g.num_cells - sub_face_map, _ = pp.fvutils.map_subgrid_to_grid( - subgrid, l2g_sub_faces, l2g_cells, is_vector=True - ) - # The sub_face_map is now a map from local sub_faces to global subfaces. - # Next we need to mat the the local sub face reconstruction "hf_cell_loc" - # onto the global grid. The cells are ordered the same, so we can use the - # cell_map from the stress computation. Similarly for the faces. - hf_cell_glob = sub_face_map * hf_cell_loc * cell_map - hf_bound_glob = sub_face_map * hf_bound_loc * face_map.T - # Next we need to eliminate the subfaces outside the active faces. - # We map from outside faces to outside subfaces - sub_outside = np.where(np.in1d(subcell_topology.fno_unique, outside))[0] - # Then expand the indices. - # The indices are ordered as first all variables of subface 1 then all variables - # of subface 2, etc. Duplicate indices for each dimension and multipy by g.dim to - # obtain correct x-index. - sub_eliminate_ind = g.dim * np.tile(sub_outside, (g.dim, 1)) - # Next add an increment to the y (and possible z) dimension to obtain correct index - # For them - sub_eliminate_ind += np.atleast_2d(np.arange(0, g.dim)).T - sub_eliminate_ind = sub_eliminate_ind.ravel("F") - # now kill the contribution of these faces - pp.fvutils.zero_out_sparse_rows(hf_cell_glob, sub_eliminate_ind) - pp.fvutils.zero_out_sparse_rows(hf_bound_glob, sub_eliminate_ind) - return stress_glob, bound_stress_glob, hf_cell_glob, hf_bound_glob, active_faces + if hf_disp: + # If we are returning the subface displacement reconstruction matrices we have + # to do some more work. The following is equivalent to what is done for the stresses, + # but as they are working on faces, the displacement reconstruction has to work on + # subfaces. + # First, we find the mappings from local subfaces to global subfaces + subcell_topology = pp.fvutils.SubcellTopology(g) + l2g_sub_faces = np.where(np.in1d(subcell_topology.fno_unique, l2g_faces))[0] + # We now create a fake grid, just to be able to use the function map_subgrid_to_grid. + subgrid = pp.CartGrid([1] * g.dim) + subgrid.num_faces = subcell_topology.fno_unique.size + subgrid.num_cells = g.num_cells + sub_face_map, _ = pp.fvutils.map_subgrid_to_grid( + subgrid, l2g_sub_faces, l2g_cells, is_vector=True + ) + # The sub_face_map is now a map from local sub_faces to global subfaces. + # Next we need to mat the the local sub face reconstruction "hf_cell_loc" + # onto the global grid. The cells are ordered the same, so we can use the + # cell_map from the stress computation. Similarly for the faces. + hf_cell_glob = sub_face_map * hf_cell_loc * cell_map + hf_bound_glob = sub_face_map * hf_bound_loc * face_map.T + # Next we need to eliminate the subfaces outside the active faces. + # We map from outside faces to outside subfaces + sub_outside = np.where(np.in1d(subcell_topology.fno_unique, outside))[0] + # Then expand the indices. + # The indices are ordered as first all variables of subface 1 then all variables + # of subface 2, etc. Duplicate indices for each dimension and multipy by g.dim to + # obtain correct x-index. + sub_eliminate_ind = g.dim * np.tile(sub_outside, (g.dim, 1)) + # Next add an increment to the y (and possible z) dimension to obtain correct index + # For them + sub_eliminate_ind += np.atleast_2d(np.arange(0, g.dim)).T + sub_eliminate_ind = sub_eliminate_ind.ravel("F") + # now kill the contribution of these faces + pp.fvutils.zero_out_sparse_rows(hf_cell_glob, sub_eliminate_ind) + pp.fvutils.zero_out_sparse_rows(hf_bound_glob, sub_eliminate_ind) + return stress_glob, bound_stress_glob, hf_cell_glob, hf_bound_glob, active_faces + else: + return stress_glob, bound_stress_glob, active_faces def _mpsa_local( @@ -1446,6 +1084,41 @@ def _mpsa_local( if bound.bc_type != "vectorial": raise AttributeError("MPSA must be given a vectorial boundary condition") + + if g.dim == 1: + tpfa_key = "tpfa_elasticity" + discr = pp.Tpfa(tpfa_key) + params = pp.Parameters(g) + + # Implicitly set Neumann boundary conditions on the whole domain. + # More general values should be permissible, but it will require handling + # of rotated boundary conditions. + if np.any(bound.is_dir): + # T + raise ValueError("have not considered Dirichlet boundary values here") + + bnd = pp.BoundaryCondition(g) + params["bc"] = bnd + + # The elasticity tensor here is set to 2*mu + lmbda, that is, the standard + # diagonal term in the stiffness matrix + k = pp.SecondOrderTensor(3, 2 * constit.mu + constit.lmbda) + params["second_order_tensor"] = k + params["aperture"] = np.ones(g.num_cells) + + d = { + pp.PARAMETERS: {tpfa_key: params}, + pp.DISCRETIZATION_MATRICES: {tpfa_key: {}}, + } + discr.discretize(g, d) + matrix_dictionary = d[pp.DISCRETIZATION_MATRICES][tpfa_key] + return ( + matrix_dictionary["flux"], + matrix_dictionary["bound_flux"], + matrix_dictionary["bound_pressure_cell"], + matrix_dictionary["bound_pressure_face"], + ) + # The grid coordinates are always three-dimensional, even if the grid is # really 2D. This means that there is not a 1-1 relation between the number # of coordinates of a point / vector and the real dimension. This again @@ -1457,11 +1130,21 @@ def _mpsa_local( # proper 2D. if g.dim == 2: g = g.copy() - g.cell_centers = np.delete(g.cell_centers, (2), axis=0) - g.face_centers = np.delete(g.face_centers, (2), axis=0) - g.face_normals = np.delete(g.face_normals, (2), axis=0) - g.nodes = np.delete(g.nodes, (2), axis=0) + cell_centers, face_normals, face_centers, _, _, nodes = pp.map_geometry.map_grid( + g + ) + g.cell_centers = cell_centers + g.face_normals = face_normals + g.face_centers = face_centers + g.nodes = nodes + + # The stiffness matrix should also be rotated before deleting rows and + # columns. However, for isotropic media, the standard __init__ for the + # FourthOrderTensor, followed by the below deletions will in effect generate + # just what we wanted (assuming we are happy with the Lame parameters, + # and do not worry about plane-strain / plane-stress consistency). + # That is all to say, this is a bit inconsistent, but it may just end up okay. constit = constit.copy() constit.values = np.delete(constit.values, (2, 5, 6, 7, 8), axis=0) constit.values = np.delete(constit.values, (2, 5, 6, 7, 8), axis=1) @@ -2078,7 +1761,7 @@ def _inverse_gradient( * pp.fvutils.invert_diagonal_blocks(grad, size_of_blocks, method=inverter) * rows2blk_diag ) - print("max igrad: ", np.max(np.abs(igrad))) + logger.debug("max igrad: " + str(np.max(np.abs(igrad)))) return igrad @@ -2593,7 +2276,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) - print("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/__init__.py b/src/porepy/numerics/interface_laws/__init__.py index fb9c7cf727..49cb98a041 100644 --- a/src/porepy/numerics/interface_laws/__init__.py +++ b/src/porepy/numerics/interface_laws/__init__.py @@ -1,3 +1,3 @@ """ Discretization of coupling terms for mixed-dimensional problems. """ -from . import elliptic_discretization +from . import elliptic_discretization, contact_mechanics_interface_laws diff --git a/src/porepy/numerics/interface_laws/contact_mechanics_interface_laws.py b/src/porepy/numerics/interface_laws/contact_mechanics_interface_laws.py new file mode 100644 index 0000000000..12d61bcb35 --- /dev/null +++ b/src/porepy/numerics/interface_laws/contact_mechanics_interface_laws.py @@ -0,0 +1,780 @@ +""" +Implementation of contact conditions for fracture mechanics, using a primal formulation. + +We provide a class for coupling the higher-dimensional mechanical discretization to the +tractions on the fractures. Also, in the case of coupled physics (Biot and the like), +classes handling the arising coupling terms are provided. +""" + +import numpy as np +import scipy.sparse as sps + +import porepy as pp + + +class PrimalContactCoupling(object): + """ Implement the coupling conditions for the pure mechanics problem. + + The primary variables for this formulation are displacement in the ambient dimension, + displacements at the boundary of the highest dimensional grid (represented as mortar + variables), and contact forces on grids of co-dimension 1. + + The conditions represented here are + 1) KKT condition for the traction / displacement in the normal direction. + 2) Conditions for the tangential traction / displacement, according + to whether the fracture is sliding, sticking or free. + 3) Linear elasticity on the surface displacements, with the tangential contact + force as a driving force. + 4) The mortar displacements act as Dirichlet boundary conditions for the + higher-dimensional domain. + + """ + + def __init__(self, keyword, discr_master, discr_slave, use_surface_discr=False): + self.keyword = keyword + self.mortar_displacement_variable = "mortar_u" + self.discr_master = discr_master + self.discr_slave = discr_slave + + self.SURFACE_DISCRETIZATION_KEY = "surface_smoother" + + self.use_surface_discr = use_surface_discr + + def _key(self): + return self.keyword + "_" + + def _discretization_key(self): + return self._key() + pp.keywords.DISCRETIZATION + + def ndof(self, mg): + """ Get the number of dof for this coupling. + + It is assumed that this method will only be called for mortar grids of + co-dimension 1. If the assumption is broken, this will not work. + """ + return (mg.dim + 1) * mg.num_cells + + def discretize(self, g_h, g_l, data_h, data_l, data_edge): + + # Discretize the surface PDE + parameter_dictionary_edge = data_edge[pp.PARAMETERS][self.keyword] + matrix_dictionary_edge = data_edge[pp.DISCRETIZATION_MATRICES][self.keyword] + + mg = data_edge["mortar_grid"] + + # Projection onto the tangential space of the mortar grid + + # Tangential_normal projection + tangential_normal_projection = data_edge["tangential_normal_projection"] + + normal_projection = tangential_normal_projection.project_normal() + + # The right hand side of the normal diffusion considers only the tangential part + # of the normal forces. + matrix_dictionary_edge["contact_force_map"] = normal_projection + + # Keyword to control if the surface discretization should be rediscretized + # It is a linear term, so we may save time during Newton iterations here + discretize_surface = parameter_dictionary_edge.get("discretize_surface", True) + + if self.use_surface_discr and discretize_surface: + # Discretize the surface pde if asked for. + + # Lame parameters to be used for discretizing the surface elliptic equation. + mu = parameter_dictionary_edge["mu"] + lmbda = parameter_dictionary_edge["lambda"] + + # Parameter used when mapping surface grids to their lower-dimensional planes. + # This is necessary for the mapping function, but at this point in the + # simulation workflow, it should not really be an issue. + deviation_from_plane_tol = 1e-5 + + # List of surface diffusion discretizations - one per side. + A_list = [] + + for _, side_grid in mg.project_to_side_grids(): + + unity = np.ones(side_grid.num_cells) + + # Create an finite volume discretization for elasticity. + # Define parameters for the surface diffusion in an appropriate form. + mpsa = pp.Mpsa(self.keyword) + + # The stiffness matrix is istropic, thus we need not care about the + # basis used for mapping grid coordinates into the tangential space. + # Simply define the parameters directly in 2d space. + stiffness = pp.FourthOrderTensor( + side_grid.dim, mu * unity, lmbda * unity + ) + + bc = pp.BoundaryConditionVectorial(side_grid) + + mpsa_parameters = pp.initialize_data( + side_grid, + {}, + self.keyword, + {"fourth_order_tensor": stiffness, "bc": bc}, + ) + + # Project the side grid into its natural dimension. + g = side_grid.copy() + # Use the same projection matrix as in the projections used on the + # variables. + rot = tangential_normal_projection.projection[:, :, 0] + if rot.shape == (2, 2): + rot = np.vstack((np.hstack((rot, np.zeros((2, 1)))), np.zeros((3)))) + cell_centers, face_normals, face_centers, _, _, nodes = pp.map_geometry.map_grid( + g, deviation_from_plane_tol, R=rot + ) + g.cell_centers = cell_centers + g.face_normals = face_normals + g.face_centers = face_centers + g.nodes = nodes + + mpsa.discretize(g, mpsa_parameters) + + # We are only interested in the elasticity discretization as a smoother. + # Construct the discretiation matrix, and disregard all other output. + A_loc = ( + pp.fvutils.vector_divergence(side_grid) + * mpsa_parameters[pp.DISCRETIZATION_MATRICES][self.keyword][ + "stress" + ] + ) + + # The local discretization must be mapped to the full mortar degrees of freedom. + # This entails a projection onto the normal plane, followed by a restriction to this + # side grid + + # Projection to remove degrees of freedom in the normal direction to the grid + # This should be used after the projection to the tangent space, + # when we know which rows are + tangential_projection = tangential_normal_projection.project_tangential( + side_grid.num_cells + ) + A_list.append(A_loc * tangential_projection) + + # Concatenate discretization matrices + A = sps.block_diag([mat for mat in A_list]) + + # The discretization is still a non-square matrix, it needs to be expanded to + # be compatible with the block assembler. + # The final equations should relate to continuity of the normal froces + matrix_dictionary_edge[self.SURFACE_DISCRETIZATION_KEY] = A + + # Discretization of the contact mechanics is done by a ColumbContact + # object. + # The resulting equations are located at the lower-dimensional grid, + # however, the discretization is inherently linked to the mortar grid. + # It is therefore constructed here. + + self.discr_slave.discretize(g_h, g_l, data_h, data_l, data_edge) + + def assemble_matrix_rhs( + self, g_master, g_slave, data_master, data_slave, data_edge, matrix + ): + + """ Assemble the dicretization of the interface law, and its impact on + the neighboring domains. + Parameters: + g_master: Grid on one neighboring subdomain. + g_slave: Grid on the other neighboring subdomain. + data_master: Data dictionary for the master suddomain + data_slave: Data dictionary for the slave subdomain. + data_edge: Data dictionary for the edge between the subdomains + matrix: original discretization matrix, to which the coupling terms will be + added. + + """ + matrix_dictionary_edge = data_edge[pp.DISCRETIZATION_MATRICES][self.keyword] + + ambient_dimension = g_master.dim + + master_ind = 0 + slave_ind = 1 + mortar_ind = 2 + + # Generate matrix for the coupling. This can probably be generalized + # once we have decided on a format for the general variables + mg = data_edge["mortar_grid"] + projection = data_edge["tangential_normal_projection"] + + dof_master = self.discr_master.ndof(g_master) + dof_slave = self.discr_slave.ndof(g_slave) + + if not dof_master == matrix[master_ind, master_ind].shape[1]: + raise ValueError( + """The number of dofs of the master discretization given + in RobinCoupling must match the number of dofs given by the matrix + """ + ) + elif not dof_slave == matrix[master_ind, slave_ind].shape[1]: + raise ValueError( + """The number of dofs of the slave discretization given + in RobinCoupling must match the number of dofs given by the matrix + """ + ) + elif not mg.num_cells * ambient_dimension == matrix[master_ind, 2].shape[1]: + raise ValueError( + """The number of dofs of the edge discretization given + in the PrimalContactCoupling must match the number of dofs given by the matrix + """ + ) + + # We know the number of dofs from the master and slave side from their + # discretizations + # dof = np.array([dof_master, dof_slave, mg.num_cells]) + dof = np.array( + [ + matrix[master_ind, master_ind].shape[1], + matrix[slave_ind, slave_ind].shape[1], + mg.num_cells, + ] + ) + cc = np.array([sps.coo_matrix((i, j)) for i in dof for j in dof]) + cc = cc.reshape((3, 3)) + + rhs = np.empty(3, dtype=np.object) + rhs[master_ind] = np.zeros(dof_master) + rhs[slave_ind] = np.zeros(dof_slave) + rhs[mortar_ind] = np.zeros(mg.num_cells * ambient_dimension) + + # IMPLEMENTATION NOTE: The current implementation is geared towards + # using mpsa for the mechanics problem. A more general approach would + # be possible - for an example see the flow problem with the RobinCoupling + # and EllipticDiscretization and its subclasses. However, at present such a general + # framework currently seems over the top, hence this more mundane approach. + + ### Equation for the master side + # The mortar variable acts as a Dirichlet boundary condition for the master. + master_bound_stress = data_master[pp.DISCRETIZATION_MATRICES][ + self.discr_master.keyword + ]["bound_stress"] + master_stress = data_master[pp.DISCRETIZATION_MATRICES][ + self.discr_master.keyword + ]["stress"] + master_bc_values = data_master[pp.PARAMETERS][self.discr_master.keyword][ + "bc_values" + ] + master_divergence = pp.fvutils.vector_divergence(g_master) + + # The mortar variable (boundary displacement) takes the form of a Dirichlet + # condition for the master side. The MPSA convention is to have + # - div * bound_stress * bc_values + # on the rhs. Accordingly, the contribution from the mortar variable (boundary + # displacement) on the left hand side is positive: + # div * bound_stress * u_mortar + cc[master_ind, mortar_ind] = ( + master_divergence + * master_bound_stress + * mg.mortar_to_master_avg(nd=ambient_dimension) + ) + + ### Equation for the slave side + # + # These are the contact conditions, which dictate relations between + # the contact forces on the slave, and the displacement jumps. + # + # NOTE: Both the contact conditions and the contact stresses are defined in the + # local coordinate system of the surface. The displacements must therefore + # be rotated to this local coordinate system during assembly. + traction_discr, displacement_jump_discr, rhs_slave = self.discr_slave.assemble_matrix_rhs( + g_slave, data_slave + ) + # The contact forces. Can be applied directly, these are in their own + # local coordinate systems. + cc[slave_ind, slave_ind] = traction_discr + + # The contact condition discretization gives coefficients for the mortar + # variables. To finalize the relation with the contact conditions, we + # (from the right) 1) assign +- signs to the two sides of the mortar, so that + # summation in reality is a difference, 2) project to the mortar grid + # 3) project to the local coordinates of the fracture, 4) assign the + # coefficients of the displacement jump. + cc[slave_ind, mortar_ind] = ( + displacement_jump_discr + * projection.project_tangential_normal(g_slave.num_cells) + * mg.mortar_to_slave_avg(nd=ambient_dimension) + * mg.sign_of_mortar_sides(nd=ambient_dimension) + ) + + # Right hand side system. In the local (surface) coordinate system. + rhs[slave_ind] = rhs_slave + + ### Equation for the mortar rows + + # This is first a stress balance: stress from the higher dimensional + # domain (both interior and bound_stress) should match with the contact stress: + # + # traction_slave + traction_master = 0 + # + # Optionally, a diffusion term can be added in the tangential direction + # of the stresses, this is currently under implementation. + + # A diagonal operator is needed to switch the sign of vectors on + # higher-dimensional faces that point into the fracture surface. The effect is to + # switch direction of the stress on boundary for the higher dimensional domain: The + # contact forces are defined as negative in contact, whereas the sign of the higher + # dimensional stresses are defined according to the direction of the normal vector. + faces_on_fracture_surface = mg.master_to_mortar_int().tocsr().indices + sign_switcher = pp.grid_utils.switch_sign_if_inwards_normal( + g_master, ambient_dimension, faces_on_fracture_surface + ) + + ## First, we obtain T_master = stress * u_master + bound_stress * u_mortar + # Stress contribution from the higher dimensional domain, projected onto + # the mortar grid + # Switch the direction of the vectors to obtain the traction as defined + # by the outwards pointing normal vector. + traction_from_master = ( + mg.master_to_mortar_int(nd=ambient_dimension) + * sign_switcher + * master_stress + ) + cc[mortar_ind, master_ind] = traction_from_master + # Stress contribution from boundary conditions. + rhs[mortar_ind] = -( + mg.master_to_mortar_int(nd=ambient_dimension) + * sign_switcher + * master_bound_stress + * master_bc_values + ) + # The stress contribution from the mortar variables, mapped to the higher + # dimensional domain via a boundary condition, and back again by a + # projection operator. + # Switch the direction of the vectors, so that for all faces, a positive + # force points into the fracture surface. + traction_from_mortar = ( + mg.master_to_mortar_int(nd=ambient_dimension) + * sign_switcher + * master_bound_stress + * mg.mortar_to_master_avg(nd=ambient_dimension) + ) + cc[mortar_ind, mortar_ind] = traction_from_mortar + + ## Second, the contact stress is mapped to the mortar grid. + # We have for the positive (first) and negative (second) side of the mortar that + # T_slave = T_master_pos = -T_master_neg, + # so we need to map the slave traction with the corresponding signs to match the + # mortar tractions. + + # The contact forces are defined in the surface coordinate system. + # Map to the mortar grid, and rotate back again to the global coordinates + # (note the inverse rotation is given by a transpose). + # Finally, the contact stresses will be felt in different directions by + # the two sides of the mortar grids (Newton's third law), hence + # adjust the signs + contact_traction_to_mortar = ( + mg.sign_of_mortar_sides(nd=ambient_dimension) + * projection.project_tangential_normal(mg.num_cells).T + * mg.slave_to_mortar_int(nd=ambient_dimension) + ) + # Minus to obtain -T_slave + T_master = 0. + cc[mortar_ind, slave_ind] = -contact_traction_to_mortar + + if self.use_surface_discr: + restrict_to_tangential_direction = projection.project_tangential( + mg.num_cells + ) + + # The first block contains the surface diffusion component. This has + # the surface diffusion operator for the mortar variables, and a + # mapping of contact forces on the slave variables. + # The second block gives continuity of forces in the normal direction. + surface_discr = matrix_dictionary_edge[self.SURFACE_DISCRETIZATION_KEY] + + cc[mortar_ind, mortar_ind] += ( + restrict_to_tangential_direction.T * surface_discr + ) + + matrix += cc + + return matrix, rhs + + +class MatrixScalarToForceBalance: + """ + This class adds the matrix scalar (pressure) contribution to the force balance posed + on the mortar grid by PrimalContactCoupling. + + We account for the scalar variable contribution to the forces on the higher-dimensional + internal boundary, i.e. the last term of: + + boundary_traction_hat = stress * u_hat + bound_stress * u_mortar + gradP * p_hat + + Note that with this approach to discretization of the boundary pressure force, it + will only be included for nonzero values of the biot_alpha coefficient. + + If the scalar is e.g. pressure, subtraction of the pressure contribution is needed: + + T_contact - p_check I \dot n = boundary_traction_hat + + This is taken care of by FracturePressureToForceBalance. + + """ + + def __init__(self, keyword, discr_master, discr_slave): + """ + Parameters: + keyword used for storage of the gradP discretization. If the GradP class is + used, this is the keyword associated with the mechanical parameters. + discr_master and + discr_slave are the discretization objects operating on the master and slave + pressure, respectively. Used for #DOFs. In FV, one cell variable is + expected. + """ + # Set node discretizations + self.discr_master = discr_master + self.discr_slave = discr_slave + # Keyword used to retrieve gradP discretization. + self.keyword = keyword + + def discretize(self, g_h, g_l, data_h, data_l, data_edge): + """ + Nothing to do + """ + pass + + def assemble_matrix_rhs( + self, g_master, g_slave, data_master, data_slave, data_edge, matrix + ): + """ + Assemble the pressure contributions of the interface force balance law. + + Parameters: + g_master: Grid on one neighboring subdomain. + g_slave: Grid on the other neighboring subdomain. + data_master: Data dictionary for the master suddomain + data_slave: Data dictionary for the slave subdomain. + data_edge: Data dictionary for the edge between the subdomains + matrix: original discretization matrix, to which the coupling terms will be + added. + """ + + ambient_dimension = g_master.dim + + master_ind = 0 + slave_ind = 1 + mortar_ind = 2 + + # Generate matrix for the coupling. This can probably be generalized + # once we have decided on a format for the general variables + mg = data_edge["mortar_grid"] + + dof_master = self.discr_master.ndof(g_master) + dof_slave = self.discr_slave.ndof(g_slave) + + if not dof_master == matrix[master_ind, master_ind].shape[1]: + raise ValueError( + """The number of dofs of the master discretization given + in RobinCoupling must match the number of dofs given by the matrix + """ + ) + elif not dof_slave == matrix[master_ind, slave_ind].shape[1]: + raise ValueError( + """The number of dofs of the slave discretization given + in RobinCoupling must match the number of dofs given by the matrix + """ + ) + elif not mg.num_cells * ambient_dimension == matrix[master_ind, 2].shape[1]: + raise ValueError( + """The number of dofs of the edge discretization given + in the PrimalContactCoupling must match the number of dofs given by the matrix + """ + ) + + # We know the number of dofs from the master and slave side from their + # discretizations + dof = np.array( + [ + matrix[master_ind, master_ind].shape[1], + matrix[slave_ind, slave_ind].shape[1], + mg.num_cells * ambient_dimension, + ] + ) + cc = np.array([sps.coo_matrix((i, j)) for i in dof for j in dof]) + cc = cc.reshape((3, 3)) + + rhs = np.empty(3, dtype=np.object) + rhs[master_ind] = np.zeros(dof_master) + rhs[slave_ind] = np.zeros(dof_slave) + rhs[mortar_ind] = np.zeros(mg.num_cells * ambient_dimension) + + master_scalar_gradient = data_master[pp.DISCRETIZATION_MATRICES][self.keyword][ + "grad_p" + ] + + # We want to modify the stress balance posed on the edge to account for the + # scalar (usually pressure) contribution. + # In the purely mechanical case, stress from the higher dimensional + # domain (both interior and bound_stress) should match the contact stress: + # -T_slave + T_master = 0, + # see PrimalContactCoupling. + # The following modification is needed: + # Add the scalar gradient contribution to the traction on the master + # boundary. + + # A diagonal operator is needed to switch the sign of vectors on + # higher-dimensional faces that point into the fracture surface, see + # PrimalContactCoupling. + faces_on_fracture_surface = mg.master_to_mortar_int().tocsr().indices + sign_switcher = pp.grid_utils.switch_sign_if_inwards_normal( + g_master, ambient_dimension, faces_on_fracture_surface + ) + + # i) Obtain pressure stress contribution from the higher dimensional domain. + # ii) Switch the direction of the vectors, so that for all faces, a positive + # force points into the fracture surface (along the outwards normal on the + # boundary). + # iii) Map to the mortar grid. + # iv) Minus according to - alpha grad p already in the discretization matrix + master_scalar_to_master_traction = ( + mg.master_to_mortar_int(nd=ambient_dimension) + * sign_switcher + * master_scalar_gradient + ) + cc[mortar_ind, master_ind] = master_scalar_to_master_traction + + matrix += cc + + return matrix, rhs + + +class FractureScalarToForceBalance: + """ + This class adds the fracture pressure contribution to the force balance posed on the + mortar grid by PrimalContactCoupling and modified to account for matrix pressure by + MatrixPressureToForceBalance. + + For the contact mechanics, we only want to consider the _contact_ traction. Thus, we + have to subtract the pressure contribution, i.e. + + T_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): + """ + Parameters: + keyword used for storage of the gradP discretization. If the GradP class is + used, this is the keyword associated with the mechanical parameters. + discr_master and + discr_slave are the discretization objects operating on the master and slave + pressure, respectively. Used for #DOFs. In FV, one cell variable is + expected. + """ + # Set node discretizations + self.discr_master = discr_master + self.discr_slave = discr_slave + + def discretize(self, g_h, g_l, data_h, data_l, data_edge): + """ + Nothing to do + """ + pass + + def assemble_matrix_rhs( + self, g_master, g_slave, data_master, data_slave, data_edge, matrix + ): + """ + Assemble the pressure contributions of the interface force balance law. + + Parameters: + g_master: Grid on one neighboring subdomain. + g_slave: Grid on the other neighboring subdomain. + data_master: Data dictionary for the master suddomain + data_slave: Data dictionary for the slave subdomain. + data_edge: Data dictionary for the edge between the subdomains + matrix: original discretization matrix, to which the coupling terms will be + added. + """ + + ambient_dimension = g_master.dim + + master_ind = 0 + slave_ind = 1 + mortar_ind = 2 + + # Generate matrix for the coupling. This can probably be generalized + # once we have decided on a format for the general variables + mg = data_edge["mortar_grid"] + + dof_master = self.discr_master.ndof(g_master) + dof_slave = self.discr_slave.ndof(g_slave) + + if not dof_master == matrix[master_ind, master_ind].shape[1]: + raise ValueError( + """The number of dofs of the master discretization given + in RobinCoupling must match the number of dofs given by the matrix + """ + ) + elif not dof_slave == matrix[master_ind, slave_ind].shape[1]: + raise ValueError( + """The number of dofs of the slave discretization given + in RobinCoupling must match the number of dofs given by the matrix + """ + ) + elif not mg.num_cells * ambient_dimension == matrix[master_ind, 2].shape[1]: + raise ValueError( + """The number of dofs of the edge discretization given + in the PrimalContactCoupling must match the number of dofs given by the matrix + """ + ) + + # We know the number of dofs from the master and slave side from their + # discretizations + dof = np.array( + [ + matrix[master_ind, master_ind].shape[1], + matrix[slave_ind, slave_ind].shape[1], + mg.num_cells * ambient_dimension, + ] + ) + cc = np.array([sps.coo_matrix((i, j)) for i in dof for j in dof]) + cc = cc.reshape((3, 3)) + + rhs = np.empty(3, dtype=np.object) + rhs[master_ind] = np.zeros(dof_master) + rhs[slave_ind] = np.zeros(dof_slave) + rhs[mortar_ind] = np.zeros(mg.num_cells * ambient_dimension) + + ## Ensure that the contact variable is only the force from the contact of the + # two sides of the fracture. This requires subtraction of the pressure force. + + # Construct the dot product between normals on fracture faces and the identity + # matrix. Similar sign switching as above is needed (this one operating on + # fracture faces only). + faces_on_fracture_surface = mg.master_to_mortar_int().tocsr().indices + sgn = g_master.sign_of_faces(faces_on_fracture_surface) + fracture_normals = g_master.face_normals[ + :ambient_dimension, faces_on_fracture_surface + ] + outwards_fracture_normals = sgn * fracture_normals + + data = outwards_fracture_normals.ravel("F") + row = np.arange(g_master.dim * mg.num_cells) + col = np.tile(np.arange(mg.num_cells), (g_master.dim, 1)).ravel("F") + n_dot_I = sps.csc_matrix((data, (row, col))) + # i) The scalar contribution to the contact stress is mapped to the mortar grid + # and multiplied by n \dot I, with n being the outwards normals on the two sides. + # Note that by using different normals for the two sides, we do not need to + # adjust the slave pressure with the corresponding signs by applying + # sign_of_mortar_sides as done in PrimalContactCoupling. + # iii) The contribution should be subtracted so that we balance the master + # forces by + # T_contact - n dot I p, + # hence the minus. + slave_pressure_to_contact_traction = -(n_dot_I * mg.slave_to_mortar_int(nd=1)) + # Minus to obtain -T_slave + T_master = 0, i.e. from placing the two + # terms on the same side of the equation, as also done in PrimalContactCoupling. + cc[mortar_ind, slave_ind] = -slave_pressure_to_contact_traction + + matrix += cc + + return matrix, rhs + + +class DivUCoupling: + """ + Coupling conditions for DivU term. + + For mixed-dimensional flow in coupled to matrix mechanics, i.e. Biot in the matrix + and conservation of a scalar quantity (usually fluid mass) in matrix and fractures. + We have assumed a primal displacement mortar variable, which will contribute + to the div u term in fracture ("div aperture") and matrix. + """ + + def __init__(self, variable, discr_master, discr_slave): + # 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 + self.discr_master = discr_master + # assemble_int_bound_displacement_source for the slave. + self.discr_slave = discr_slave + + def discretize(self, g_h, g_l, data_h, data_l, data_edge): + """ + Nothing to do + """ + pass + + def assemble_matrix_rhs( + self, g_master, g_slave, data_master, data_slave, data_edge, matrix + ): + """ + Assemble the mortar displacement's contribution as a internal Dirichlet + contribution for the higher dimension, and source term for the lower dimension. + Parameters: + g_master: Grid on one neighboring subdomain. + g_slave: Grid on the other neighboring subdomain. + data_master: Data dictionary for the master suddomain + data_slave: Data dictionary for the slave subdomain. + data_edge: Data dictionary for the edge between the subdomains + matrix: original discretization matrix, to which the coupling terms will be + added. + """ + ambient_dimension = g_master.dim + + master_ind = 0 + slave_ind = 1 + mortar_ind = 2 + + # Generate matrix for the coupling. This can probably be generalized + # once we have decided on a format for the general variables + mg = data_edge["mortar_grid"] + + dof_master = self.discr_master.ndof(g_master) + dof_slave = self.discr_slave.ndof(g_slave) + + if not dof_master == matrix[master_ind, master_ind].shape[1]: + raise ValueError( + """The number of dofs of the master discretization given + in RobinCoupling must match the number of dofs given by the matrix + """ + ) + elif not dof_slave == matrix[master_ind, slave_ind].shape[1]: + raise ValueError( + """The number of dofs of the slave discretization given + in RobinCoupling must match the number of dofs given by the matrix + """ + ) + elif not mg.num_cells * ambient_dimension == matrix[master_ind, 2].shape[1]: + raise ValueError( + """The number of dofs of the edge discretization given + in the PrimalContactCoupling must match the number of dofs given by the matrix + """ + ) + + # We know the number of dofs from the master and slave side from their + # discretizations + dof = np.array( + [ + matrix[master_ind, master_ind].shape[1], + matrix[slave_ind, slave_ind].shape[1], + mg.num_cells * ambient_dimension, + ] + ) + cc = np.array([sps.coo_matrix((i, j)) for i in dof for j in dof]) + cc = cc.reshape((3, 3)) + rhs = np.empty(3, dtype=np.object) + rhs[master_ind] = np.zeros(dof_master) + rhs[slave_ind] = np.zeros(dof_slave) + rhs[mortar_ind] = np.zeros(mg.num_cells * ambient_dimension) + + grid_swap = False + # Let the DivU class assemble the contribution from mortar to master + self.discr_master.assemble_int_bound_displacement_trace( + g_master, data_master, data_edge, grid_swap, cc, matrix, rhs, master_ind + ) + # and from mortar to slave. + self.discr_slave.assemble_int_bound_displacement_source( + g_slave, data_slave, data_edge, cc, matrix, rhs, slave_ind + ) + matrix += cc + + return matrix, rhs diff --git a/src/porepy/numerics/interface_laws/elliptic_discretization.py b/src/porepy/numerics/interface_laws/elliptic_discretization.py index eabc74b9f8..d942091628 100644 --- a/src/porepy/numerics/interface_laws/elliptic_discretization.py +++ b/src/porepy/numerics/interface_laws/elliptic_discretization.py @@ -332,296 +332,3 @@ def enforce_neumann_int_bound( matrix (scipy.sparse.matrix): Discretization matrix to be modified. """ raise NotImplementedError("Method not implemented") - - -class VectorEllipticDiscretization: - """ Superclass for finite volume discretizations of the vector elliptic equation. - - Should not be used by itself, instead use a subclass that implements an - actual discretization method. Known subclass is Mpsa. - """ - - def __init__(self, keyword): - """ Set the discretization, with the keyword used for storing various - information associated with the discretization. - - Paramemeters: - keyword (str): Identifier of all information used for this - discretization. - """ - self.keyword = keyword - - def _key(self): - """ Get the keyword of this object, on a format friendly to access relevant - fields in the data dictionary - - Returns: - String, on the form self.keyword + '_'. - - """ - return self.keyword + "_" - - def ndof(self, g): - """ Abstract method. Return the number of degrees of freedom associated to the - method. - - Parameter: - g (grid): Computational grid - - Returns: - int: the number of degrees of freedom. - - """ - raise NotImplementedError("Method not implemented") - - def extract_displacement(self, g, solution_array, d): - """ Abstract method. Extract the displacement part of a solution. - - The implementation will depend on what the primary variables of the specific - implementation are. - - Parameters: - g (grid): To which the solution array belongs. - solution_array (np.array): Solution for this grid obtained from - either a mono-dimensional or a mixed-dimensional problem. - d (dictionary): Data dictionary associated with the grid. Not used, - but included for consistency reasons. - - Returns: - np.array (g.num_cells): Displacement solution vector. - - Raises: - NotImplementedError if the method - """ - raise NotImplementedError("Method not implemented") - - def extract_traction(self, g, solution_array, d): - """ Abstract method. Extract the traction part of a solution. - - The implementation will depend on what the primary variables of the specific - implementation are. - - Parameters: - g (grid): To which the solution array belongs. - solution_array (np.array): Solution for this grid obtained from - either a mono-dimensional or a mixed-dimensional problem. Will - correspond to the displacement solution. - d (dictionary): Data dictionary associated with the grid. - - Returns: - np.array (g.num_faces): Traction vector. - - """ - raise NotImplementedError() - - # ------------------------------------------------------------------------------# - - def assemble_matrix_rhs(self, g, data): - """ Return the matrix and right-hand side for a discretization of a second - order elliptic vector equation. - - - - Parameters: - g : grid, or a subclass, with geometry fields computed. - data: dictionary to store the data. For details on necessary keywords, - see method discretize() - discretize (boolean, optional): default True. Whether to discetize prior to - matrix assembly. If False, data should already contain discretization. - - Returns: - matrix: sparse csr (self.ndof, self.ndof) discretization matrix. - rhs: np.ndarray (self.ndof) Right-hand side which contains the boundary - conditions and the vector source term. - """ - raise NotImplementedError("Method not implemented") - - def assemble_matrix(self, g, data): - """Return the matrix for a discretization of a second order elliptic vector - equation. - - Parameters: - g (Grid): Computational grid, with geometry fields computed. - data (dictionary): With data stored. - - Returns: - scipy.sparse.csr_matrix: System matrix of this discretization. The - size of the matrix will depend on the specific discretization. - """ - raise NotImplementedError("Method not implemented") - - # ------------------------------------------------------------------------------# - - def assemble_rhs(self, g, data): - """ Return the right-hand side for a discretization of a second order elliptic - vector equation. - - Also discretize the necessary operators if the data dictionary does not - contain a discretization of the boundary equation. - - Parameters: - g (Grid): Computational grid, with geometry fields computed. - data (dictionary): With data stored. - - Returns: - np.ndarray: Right hand side vector with representation of boundary - conditions. The size of the vector will depend on the discretization. - """ - raise NotImplementedError("Method not implemented") - - def assemble_int_bound_stress( - self, g, data, data_edge, grid_swap, cc, matrix, rhs, self_ind - ): - """Assemble the contribution from an internal boundary, manifested as a - stress boundary condition. - - The intended use is when the internal boundary is coupled to another - node in a mixed-dimensional method. Specific usage depends on the - interface condition between the nodes; this method will typically be - used to impose stress continuity on a higher-dimensional domain. - - Implementations of this method will use an interplay between the grid - on the node and the mortar grid on the relevant edge. - - Parameters: - g (Grid): Grid which the condition should be imposed on. - data (dictionary): Data dictionary for the node in the - mixed-dimensional grid. - data_edge (dictionary): Data dictionary for the edge in the - mixed-dimensional grid. - grid_swap (boolean): If True, the grid g is identified with the @ - slave side of the mortar grid in data_adge. - cc (block matrix, 3x3): Block matrix for the coupling condition. - The first and second rows and columns are identified with the - master and slave side; the third belongs to the edge variable. - The discretization of the relevant term is done in-place in cc. - matrix (block matrix 3x3): Discretization matrix for the edge and - the two adjacent nodes. - rhs (block_array 3x1): Right hand side contribution for the edge and - the two adjacent nodes. - self_ind (int): Index in cc and matrix associated with this node. - Should be either 1 or 2. - - """ - raise NotImplementedError("Method not implemented") - - def assemble_int_bound_source( - self, g, data, data_edge, grid_swap, cc, matrix, rhs, self_ind - ): - """ Abstract method. Assemble the contribution from an internal - boundary, manifested as a body force term. - - The intended use is when the internal boundary is coupled to another - node in a mixed-dimensional method. Specific usage depends on the - interface condition between the nodes; this method will typically be - used to impose stress continuity on a lower-dimensional domain. - - Implementations of this method will use an interplay between the grid on - the node and the mortar grid on the relevant edge. - - Parameters: - g (Grid): Grid which the condition should be imposed on. - data (dictionary): Data dictionary for the node in the - mixed-dimensional grid. - data_edge (dictionary): Data dictionary for the edge in the - mixed-dimensional grid. - grid_swap (boolean): If True, the grid g is identified with the @ - slave side of the mortar grid in data_adge. - cc (block matrix, 3x3): Block matrix for the coupling condition. - The first and second rows and columns are identified with the - master and slave side; the third belongs to the edge variable. - The discretization of the relevant term is done in-place in cc. - matrix (block matrix 3x3): Discretization matrix for the edge and - the two adjacent nodes. - rhs (block_array 3x1): Right hand side contribution for the edge and - the two adjacent nodes. - self_ind (int): Index in cc and matrix associated with this node. - Should be either 1 or 2. - """ - raise NotImplementedError("Method not implemented") - - def assemble_int_bound_displacement_trace( - self, g, data, data_edge, grid_swap, cc, matrix, rhs, self_ind - ): - """ Abstract method. Assemble the contribution from an internal - boundary, manifested as a condition on the boundary displacement. - - The intended use is when the internal boundary is coupled to another - node in a mixed-dimensional method. Specific usage depends on the - interface condition between the nodes; this method will typically be - used to impose stress continuity on a higher-dimensional domain. - - Implementations of this method will use an interplay between the grid on - the node and the mortar grid on the relevant edge. - - Parameters: - g (Grid): Grid which the condition should be imposed on. - data (dictionary): Data dictionary for the node in the - mixed-dimensional grid. - data_edge (dictionary): Data dictionary for the edge in the - mixed-dimensional grid. - grid_swap (boolean): If True, the grid g is identified with the @ - slave side of the mortar grid in data_adge. - cc (block matrix, 3x3): Block matrix for the coupling condition. - The first and second rows and columns are identified with the - master and slave side; the third belongs to the edge variable. - The discretization of the relevant term is done in-place in cc. - matrix (block matrix 3x3): Discretization matrix for the edge and - the two adjacent nodes. - rhs (block_array 3x1): Right hand side contribution for the edge and - the two adjacent nodes. - self_ind (int): Index in cc and matrix associated with this node. - Should be either 1 or 2. - - """ - raise NotImplementedError("Method not implemented") - - def assemble_int_bound_displacement_cell( - self, g, data, data_edge, grid_swap, cc, matrix, rhs, self_ind - ): - """ Abstract method. Assemble the contribution from an internal - boundary, manifested as a condition on the cell displacement. - - The intended use is when the internal boundary is coupled to another - node in a mixed-dimensional method. Specific usage depends on the - interface condition between the nodes; this method will typically be - used to impose stress continuity on a lower-dimensional domain. - - Implementations of this method will use an interplay between the grid on - the node and the mortar grid on the relevant edge. - - Parameters: - g (Grid): Grid which the condition should be imposed on. - data (dictionary): Data dictionary for the node in the - mixed-dimensional grid. - data_edge (dictionary): Data dictionary for the edge in the - mixed-dimensional grid. - grid_swap (boolean): If True, the grid g is identified with the @ - slave side of the mortar grid in data_adge. - cc (block matrix, 3x3): Block matrix for the coupling condition. - The first and second rows and columns are identified with the - master and slave side; the third belongs to the edge variable. - The discretization of the relevant term is done in-place in cc. - matrix (block matrix 3x3): Discretization matrix for the edge and - the two adjacent nodes. - rhs (block_array 3x1): Right hand side contribution for the edge and - the two adjacent nodes. - self_ind (int): Index in cc and matrix associated with this node. - Should be either 1 or 2. - """ - raise NotImplementedError("Method not implemented") - - def enforce_neumann_int_bound( - self, g_master, data_edge, matrix, swap_grid, self_ind - ): - """ Enforce Neumann boundary conditions on a given system matrix. - - The method is void for finite volume approaches, but is implemented - to be compatible with the general framework. - - Parameters: - g (Grid): On which the equation is discretized - data (dictionary): Of data related to the discretization. - matrix (scipy.sparse.matrix): Discretization matrix to be modified. - """ - raise NotImplementedError("Method not implemented") diff --git a/src/porepy/numerics/interface_laws/elliptic_interface_laws.py b/src/porepy/numerics/interface_laws/elliptic_interface_laws.py index 5847e42bb1..d9fde73923 100644 --- a/src/porepy/numerics/interface_laws/elliptic_interface_laws.py +++ b/src/porepy/numerics/interface_laws/elliptic_interface_laws.py @@ -23,7 +23,6 @@ class RobinCoupling(object): """ def __init__(self, keyword, discr_master, discr_slave=None): - # @ALL should the node discretization default to Tpfa? self.keyword = keyword if discr_slave is None: discr_slave = discr_master @@ -75,8 +74,6 @@ def discretize(self, g_h, g_l, data_h, data_l, data_edge): Eta = sps.diags(np.divide(inv_k, proj * aperture_h[cells_h])) - # @ALESSIO, @EIRIK: the tpfa and vem couplers use different sign - # conventions here. We should be very careful. matrix_dictionary_edge["Robin_discr"] = -inv_M * Eta def assemble_matrix_rhs( @@ -398,678 +395,3 @@ def assemble_matrix_rhs( ) return matrix, rhs - - -# ------------------------------------------------------------------------------ - - -class RobinContact(object): - """ - Contact condition for elastic problem. This condition defines a Robin condition - for the stress and displacement jump between slave and master boundaries. g_slave - and g_master must have the same dimension. - - The contact condition is Newton's third law - \sigma \cdot n_slave = -\sigma \cdot n_master, - i.e., traction on the two sides must be equal and oposite, and a Robin-type condition - on the displacement jump - MW * \lambda + RW [u] = robin_rhs - where MW and RW are matrices of size (g_slave.dim, g.slave.dim), and - \labmda = \sigma \cdot \n_slave. - The jump operator [\cdot] is given by - [v] = v_slave - v_master, - and robin_rhs is a given rhs. - """ - - def __init__(self, keyword, discr_master, discr_slave=None): - self.keyword = keyword - if discr_slave is None: - discr_slave = discr_master - self.discr_master = discr_master - self.discr_slave = discr_slave - - def _key(self): - return self.keyword + "_" - - def _discretization_key(self): - return self._key() + pp.keywords.DISCRETIZATION - - def ndof(self, mg): - return (mg.dim + 1) * mg.num_cells - - def discretize(self, g_h, g_l, data_h, data_l, data_edge): - """ - Discretize the Mortar coupling. - We assume the following two sub-dictionaries to be present in the data_edge - dictionary: - parameter_dictionary, storing all parameters. - Stored in data[pp.PARAMETERS][self.keyword]. - matrix_dictionary, for storage of discretization matrices. - Stored in data[pp.DISCRETIZATION_MATRICES][self.keyword] - - parameter_dictionary contains the entries: - robin_weigth (list): a list of mortar_grid.num_cells np.ndarrays of - shape (mortar_grid.dim + 1, mortar_grid.dim + 1) giving the displacement - jump weight. - mortar_weigth (list): a list of mortar_grid.num_cells np.ndarrays of - shape (mortar_grid.dim + 1, mortar_grid.dim + 1) giving the mortar - - matrix_dictionary will be updated with the following entries: - mortar_weigth: sps.csc_matrix (mg.num_cells * mg.dim, mg.num_cells * mg.dim) - The weight matrix applied to the mortar variables. - robin_weigth: sps.csc_matrix (mg.num_cells * mg.dim, mg.num_cells * mg.dim) - The weight matrix applied to the displacement jump. - - Parameters: - g_h: Grid of the master domanin. - g_l: Grid of the slave domain. - data_h: Data dictionary for the master domain. - data_l: Data dictionary for the slave domain. - data_edge: Data dictionary for the edge between the domains. - - """ - matrix_dictionary_edge = data_edge[pp.DISCRETIZATION_MATRICES][self.keyword] - parameter_dictionary_edge = data_edge[pp.PARAMETERS][self.keyword] - - mortar_weight = sps.block_diag(parameter_dictionary_edge["mortar_weight"]) - robin_weight = sps.block_diag(parameter_dictionary_edge["robin_weight"]) - robin_rhs = parameter_dictionary_edge["robin_rhs"] - matrix_dictionary_edge["mortar_weight"] = mortar_weight - matrix_dictionary_edge["robin_weight"] = robin_weight - matrix_dictionary_edge["robin_rhs"] = robin_rhs - - def assemble_matrix_rhs( - self, g_master, g_slave, data_master, data_slave, data_edge, matrix - ): - """ Assemble the dicretization of the interface law, and its impact on - the neighboring domains. - Parameters: - g_master: Grid on one neighboring subdomain. - g_slave: Grid on the other neighboring subdomain. - data_master: Data dictionary for the master suddomain - data_slave: Data dictionary for the slave subdomain. - data_edge: Data dictionary for the edge between the subdomains - matrix_master: original discretization for the master subdomain - matrix_slave: original discretization for the slave subdomain - - """ - matrix_dictionary_edge = data_edge[pp.DISCRETIZATION_MATRICES][self.keyword] - - self.discretize(g_master, g_slave, data_master, data_slave, data_edge) - - if not g_master.dim == g_slave.dim: - raise AssertionError("Slave and master must have same dimension") - - master_ind = 0 - slave_ind = 1 - - # Generate matrix for the coupling. This can probably be generalized - # once we have decided on a format for the general variables - mg = data_edge["mortar_grid"] - - dof_master = self.discr_master.ndof(g_master) - dof_slave = self.discr_slave.ndof(g_slave) - - if not dof_master == matrix[master_ind, master_ind].shape[1]: - raise ValueError( - """The number of dofs of the master discretization given - in RobinContact must match the number of dofs given by the matrix - """ - ) - elif not dof_slave == matrix[master_ind, slave_ind].shape[1]: - raise ValueError( - """The number of dofs of the slave discretization given - in RobinContact must match the number of dofs given by the matrix - """ - ) - elif not self.ndof(mg) == matrix[master_ind, 2].shape[1]: - raise ValueError( - """The number of dofs of the edge discretization given - in RobinContact must match the number of dofs given by the matrix - """ - ) - # We know the number of dofs from the master and slave side from their - # discretizations - dof = np.array( - [ - matrix[master_ind, master_ind].shape[1], - matrix[slave_ind, slave_ind].shape[1], - self.ndof(mg), - ] - ) - cc = np.array([sps.coo_matrix((i, j)) for i in dof for j in dof]) - cc_master = cc.reshape((3, 3)) - cc_slave = cc_master.copy() - cc_mortar = cc_master.copy() - - # EK: For some reason, the following lines were necessary to apease python - rhs_slave = np.empty(3, dtype=np.object) - rhs_slave[master_ind] = np.zeros(dof_master) - rhs_slave[slave_ind] = np.zeros(dof_slave) - rhs_slave[2] = np.zeros(self.ndof(mg)) - # I got some problems with pointers when doing rhs_master = rhs_slave.copy() - # so just reconstruct everything. - rhs_master = np.empty(3, dtype=np.object) - rhs_master[master_ind] = np.zeros(dof_master) - rhs_master[slave_ind] = np.zeros(dof_slave) - rhs_master[2] = np.zeros(self.ndof(mg)) - - # The convention, for now, is to put the master grid information - # in the first column and row in matrix, slave grid in the second - # and mortar variables in the third - # If master and slave is the same grid, they should contribute to the same - # row and coloumn. When the assembler assigns matrix[idx] it will only add - # the slave information because of duplicate indices (master and slave is the same). - # We therefore write the both master and slave info to the slave index. - if g_master == g_slave: - master_ind = 1 - else: - master_ind = 0 - - # Obtain the displacement trace u_master - self.discr_master.assemble_int_bound_displacement_trace( - g_master, - data_master, - data_edge, - False, - cc_master, - matrix, - rhs_master, - master_ind, - ) - # set \sigma_master = -\lamba - self.discr_master.assemble_int_bound_stress( - g_master, - data_master, - data_edge, - False, - cc_master, - matrix, - rhs_master, - master_ind, - ) - # Obtain the displacement trace u_slave - self.discr_slave.assemble_int_bound_displacement_trace( - g_slave, data_slave, data_edge, True, cc_slave, matrix, rhs_slave, slave_ind - ) - # set \sigma_slave = \lamba - self.discr_slave.assemble_int_bound_stress( - g_slave, data_slave, data_edge, True, cc_slave, matrix, rhs_slave, slave_ind - ) - # We now have to flip the sign of some of the matrices - # First we flip the sign of the master stress because the mortar stress - # is defined from the slave stress. Then, stress_master = -\lambda - cc_master[master_ind, 2] = -cc_master[master_ind, 2] - # Then we flip the sign for the master displacement since the displacement - # jump is defined as u_slave - u_master - cc_master[2, master_ind] = -cc_master[2, master_ind] - rhs_master[2] = -rhs_master[2] - # Note that cc_master[2, 2] is fliped twice, first in Newton's third law, - # then for the displacement jump. - - # now, the matrix cc = cc_slave + cc_master expresses the stress and displacement - # continuities over the mortar grid. - # cc[0] -> stress_master = mortar_stress - # cc[1] -> stress_slave = -mortar_stress - # cc[2] -> mortar_weight * lambda + robin_weight * (u_slave - u_master) = robin_rhs - - # We don't want to enforce the displacement jump, but a Robin condition. - # We therefore add the mortar variable to the last equation. - cc_mortar[2, 2] = matrix_dictionary_edge["mortar_weight"] - - # The displacement jump is scaled by a matrix in the Robin condition: - robin_weight = matrix_dictionary_edge["robin_weight"] - cc_sm = cc_master + cc_slave - rhs = rhs_slave + rhs_master - rhs[2] = robin_weight * rhs[2] - for i in range(3): - cc_sm[2, i] = robin_weight * cc_sm[2, i] - - # Now define the complete Robin condition: - # mortar_weight * \lambda + "robin_weight" * [u] = robin_rhs - matrix += cc_sm + cc_mortar - rhs[2] += matrix_dictionary_edge["robin_rhs"] - - # The following two functions might or might not be needed when using - # a finite element discretization (see RobinCoupling for flow). - self.discr_master.enforce_neumann_int_bound( - g_master, data_edge, matrix, False, master_ind - ) - self.discr_slave.enforce_neumann_int_bound( - g_slave, data_edge, matrix, True, slave_ind - ) - - return matrix, rhs - - -class StressDisplacementContinuity(RobinContact): - """ - Contact condition for elastic problem. This condition defines continuity for - the stress and displacement jump between slave and master boundaries. g_slave - and g_master must have the same dimension. - - This contact condition is equivalent as if the slave and master domain was - a single connected domain (the discrete solution will be different as the - discretization will be slightly different). - """ - - def discretize(self, g_h, g_l, data_h, data_l, data_edge): - """ Nothing really to do here - - Parameters: - g_h: Grid of the master domanin. - g_l: Grid of the slave domain. - data_h: Data dictionary for the master domain. - data_l: Data dictionary for the slave domain. - data_edge: Data dictionary for the edge between the domains. - - """ - pass - - def assemble_matrix_rhs( - self, g_master, g_slave, data_master, data_slave, data_edge, matrix - ): - """ Assemble the dicretization of the interface law, and its impact on - the neighboring domains. - Parameters: - ---------- - g_master: Grid on one neighboring subdomain. - g_slave: Grid on the other neighboring subdomain. - data_master: Data dictionary for the master suddomain - data_slave: Data dictionary for the slave subdomain. - data_edge: Data dictionary for the edge between the subdomains - matrix_master: original discretization for the master subdomain - matrix_slave: original discretization for the slave subdomain - - """ - - self.discretize(g_master, g_slave, data_master, data_slave, data_edge) - - if not g_master.dim == g_slave.dim: - raise AssertionError("Slave and master must have same dimension") - - master_ind = 0 - slave_ind = 1 - - # Generate matrix for the coupling. This can probably be generalized - # once we have decided on a format for the general variables - mg = data_edge["mortar_grid"] - - dof_master = self.discr_master.ndof(g_master) - dof_slave = self.discr_slave.ndof(g_slave) - - if not dof_master == matrix[master_ind, master_ind].shape[1]: - raise ValueError( - """The number of dofs of the master discretization given - in FluxPressureContinuity must match the number of dofs given by the matrix - """ - ) - elif not dof_slave == matrix[master_ind, slave_ind].shape[1]: - raise ValueError( - """The number of dofs of the slave discretization given - in FluxPressureContinuity must match the number of dofs given by the matrix - """ - ) - elif not self.ndof(mg) == matrix[master_ind, 2].shape[1]: - raise ValueError( - """The number of dofs of the edge discretization given - in FluxPressureContinuity must match the number of dofs given by the matrix - """ - ) - # We know the number of dofs from the master and slave side from their - # discretizations - dof = np.array( - [ - matrix[master_ind, master_ind].shape[1], - matrix[slave_ind, slave_ind].shape[1], - self.ndof(mg), - ] - ) - cc = np.array([sps.coo_matrix((i, j)) for i in dof for j in dof]) - cc_master = cc.reshape((3, 3)) - cc_slave = cc_master.copy() - # The rhs is just zeros - # EK: For some reason, the following lines were necessary to apease python - rhs_slave = np.empty(3, dtype=np.object) - rhs_slave[master_ind] = np.zeros(dof_master) - rhs_slave[slave_ind] = np.zeros(dof_slave) - rhs_slave[2] = np.zeros(self.ndof(mg)) - # I got some problems with pointers when doing rhs_master = rhs_slave.copy() - # so just reconstruct everything. - rhs_master = np.empty(3, dtype=np.object) - rhs_master[master_ind] = np.zeros(dof_master) - rhs_master[slave_ind] = np.zeros(dof_slave) - rhs_master[2] = np.zeros(self.ndof(mg)) - - # The convention, for now, is to put the master grid information - # in the first column and row in matrix, slave grid in the second - # and mortar variables in the third - # If master and slave is the same grid, they should contribute to the same - # row and coloumn. When the assembler assigns matrix[idx] it will only add - # the slave information because of duplicate indices (master and slave is the same). - # We therefore write the both master and slave info to the slave index. - if g_master == g_slave: - master_ind = 1 - else: - master_ind = 0 - - # Obtain the displacement trace u_master - self.discr_master.assemble_int_bound_displacement_trace( - g_master, - data_master, - data_edge, - False, - cc_master, - matrix, - rhs_master, - master_ind, - ) - # set \sigma_master = -\lamba - self.discr_master.assemble_int_bound_stress( - g_master, - data_master, - data_edge, - False, - cc_master, - matrix, - rhs_master, - master_ind, - ) - # Obtain the displacement trace u_slave - self.discr_slave.assemble_int_bound_displacement_trace( - g_slave, data_slave, data_edge, True, cc_slave, matrix, rhs_slave, slave_ind - ) - # set \sigma_slave = \lamba - self.discr_slave.assemble_int_bound_stress( - g_slave, data_slave, data_edge, True, cc_slave, matrix, rhs_slave, slave_ind - ) - # We now have to flip the sign of some of the matrices - # First we flip the sign of the master stress because the mortar stress - # is defined from the slave stress. Then, stress_master = -\lambda - cc_master[master_ind, 2] = -cc_master[master_ind, 2] - # Then we flip the sign for the master displacement since the displacement - # jump is defined as u_slave - u_master - cc_master[2, master_ind] = -cc_master[2, master_ind] - rhs_master[2] = -rhs_master[2] - # Note that cc_master[2, 2] is fliped twice, first in Newton's third law, - # then for the displacement jump. - - # now, the matrix cc = cc_slave + cc_master expresses the stress and displacement - # continuities over the mortar grid. - # cc[0] -> stress_master = mortar_stress - # cc[1] -> stress_slave = -mortar_stress - # cc[2] -> u_slave - u_master = 0 - - matrix += cc_master + cc_slave - rhs = rhs_slave + rhs_master - # The following two functions might or might not be needed when using - # a finite element discretization (see RobinCoupling for flow). - self.discr_master.enforce_neumann_int_bound( - g_master, data_edge, matrix, False, master_ind - ) - - # Consider this terms only if the grids are of the same dimension - if g_master.dim == g_slave.dim: - self.discr_slave.enforce_neumann_int_bound( - g_slave, data_edge, matrix, True, slave_ind - ) - - return matrix, rhs - - -class RobinContactBiotPressure(RobinContact): - """ - This condition adds the fluid pressure contribution to the Robin contact condition. - The Robin condition says: - MW * lambda + RW * [u] = robin_rhs, - where MW (mortar_weight) and RW (robin_weight) are two matrices, and - [u] = u_slave - u_master is the displacement jump from the slave to the master. - In Biot the displacement on the contact boundary (u_slave and u_master) will be a - linear function of cell center displacement (u), mortar stress (lambda) and cell - centere fluid pressure (p): - A * u + B * p + C * lam = u_slave/u_master - This class adds the contribution B. - - To enforce the full continuity this interface law must be used in combination with - the RobinContact conditions which adds the contributions A and C - """ - - def discretize(self, g_h, g_l, data_h, data_l, data_edge): - """ Discretize the robin weight (RW) - - Parameters: - g_h: Grid of the master domanin. - g_l: Grid of the slave domain. - data_h: Data dictionary for the master domain. - data_l: Data dictionary for the slave domain. - data_edge: Data dictionary for the edge between the domains. - - """ - matrix_dictionary_edge = data_edge[pp.DISCRETIZATION_MATRICES][self.keyword] - parameter_dictionary_edge = data_edge[pp.PARAMETERS][self.keyword] - - robin_weight = sps.block_diag(parameter_dictionary_edge["robin_weight"]) - matrix_dictionary_edge["robin_weight"] = robin_weight - - def assemble_matrix_rhs( - self, g_master, g_slave, data_master, data_slave, data_edge, matrix - ): - """ Assemble the dicretization of the interface law, and its impact on - the neighboring domains. - Parameters: - ---------- - g_master: Grid on one neighboring subdomain. - g_slave: Grid on the other neighboring subdomain. - data_master: Data dictionary for the master suddomain - data_slave: Data dictionary for the slave subdomain. - data_edge: Data dictionary for the edge between the subdomains - matrix_master: original discretization for the master subdomain - matrix_slave: original discretization for the slave subdomain - - """ - matrix_dictionary_edge = data_edge[pp.DISCRETIZATION_MATRICES][self.keyword] - - self.discretize(g_master, g_slave, data_master, data_slave, data_edge) - - if not g_master.dim == g_slave.dim: - raise AssertionError("Slave and master must have same dimension") - - master_ind = 0 - slave_ind = 1 - - # Generate matrix for the coupling. This can probably be generalized - # once we have decided on a format for the general variables - mg = data_edge["mortar_grid"] - - # We know the number of dofs from the master and slave side from their - # discretizations - dof = np.array( - [ - matrix[master_ind, master_ind].shape[1], - matrix[slave_ind, slave_ind].shape[1], - self.ndof(mg), - ] - ) - cc = np.array([sps.coo_matrix((i, j)) for i in dof for j in dof]) - cc_master = cc.reshape((3, 3)) - cc_slave = cc_master.copy() - - # The rhs is just zeros - # EK: For some reason, the following lines were necessary to apease python - rhs = np.empty(3, dtype=np.object) - rhs[master_ind] = np.zeros(matrix[master_ind, master_ind].shape[1]) - rhs[slave_ind] = np.zeros(matrix[slave_ind, slave_ind].shape[1]) - rhs[2] = np.zeros(self.ndof(mg)) - - # The convention, for now, is to put the master grid information - # in the first column and row in matrix, slave grid in the second - # and mortar variables in the third - # If master and slave is the same grid, they should contribute to the same - # row and coloumn. When the assembler assigns matrix[idx] it will only add - # the slave information because of duplicate indices (master and slave is the same). - # We therefore write the both master and slave info to the slave index. - if g_master == g_slave: - master_ind = 1 - else: - master_ind = 0 - - # Obtain the contribution of the cell centered pressure on the displacement - # trace u_master - self.discr_master.assemble_int_bound_displacement_trace( - g_master, data_master, data_edge, False, cc_master, matrix, rhs, master_ind - ) - # Equivalent for u_slave - self.discr_slave.assemble_int_bound_displacement_trace( - g_slave, data_slave, data_edge, True, cc_slave, matrix, rhs, slave_ind - ) - # We now have to flip the sign of some of the matrices - # First we flip the sign of the master stress because the mortar stress - # is defined from the slave stress. Then, stress_master = -\lambda - cc_master[master_ind, 2] = -cc_master[master_ind, 2] - # Then we flip the sign for the master displacement since the displacement - # jump is defined as u_slave - u_master - cc_master[2, master_ind] = -cc_master[2, master_ind] - - matrix += cc_master + cc_slave - - # The displacement jump is scaled by a matrix in the Robin condition: - robin_weight = matrix_dictionary_edge["robin_weight"] - - for i in range(3): - matrix[2, i] = robin_weight * matrix[2, i] - - # The following two functions might or might not be needed when using - # a finite element discretization (see RobinCoupling for flow). - self.discr_master.enforce_neumann_int_bound( - g_master, data_edge, matrix, False, master_ind - ) - self.discr_slave.enforce_neumann_int_bound( - g_slave, data_edge, matrix, True, slave_ind - ) - - return matrix, rhs - - -class DivU_StressMortar(RobinContactBiotPressure): - """ - This condition adds the stress mortar contribution to the div u term in the - fluid mass conservation equation of the Biot equations. When fractures are - present the divergence of u (div_u) will be a function of cell centere displacement, - boundary conditions and the stress mortar (lambda): - div_u = A * u + B * u_bc_val + C * lambda - The class adds the contribution C, while the DivD discretization adds A and B. - """ - - def discretize(self, g_h, g_l, data_h, data_l, data_edge): - """ Nothing really to do here - - Parameters: - g_h: Grid of the master domanin. - g_l: Grid of the slave domain. - data_h: Data dictionary for the master domain. - data_l: Data dictionary for the slave domain. - data_edge: Data dictionary for the edge between the domains. - - """ - pass - - def assemble_matrix_rhs( - self, g_master, g_slave, data_master, data_slave, data_edge, matrix - ): - """ Assemble the dicretization of the interface law, and its impact on - the neighboring domains. - Parameters: - ---------- - g_master: Grid on one neighboring subdomain. - g_slave: Grid on the other neighboring subdomain. - data_master: Data dictionary for the master suddomain - data_slave: Data dictionary for the slave subdomain. - data_edge: Data dictionary for the edge between the subdomains - matrix_master: original discretization for the master subdomain - matrix_slave: original discretization for the slave subdomain - - """ - - if not g_master.dim == g_slave.dim: - raise AssertionError("Slave and master must have same dimension") - - master_ind = 0 - slave_ind = 1 - - # Generate matrix for the coupling. This can probably be generalized - # once we have decided on a format for the general variables - mg = data_edge["mortar_grid"] - - # We know the number of dofs from the master and slave side from their - # discretizations - dof = np.array( - [ - matrix[master_ind, master_ind].shape[1], - matrix[slave_ind, slave_ind].shape[1], - self.ndof(mg), - ] - ) - cc = np.array([sps.coo_matrix((i, j)) for i in dof for j in dof]) - cc_master = cc.reshape((3, 3)) - cc_slave = cc_master.copy() - - # When we do time stepping in Biot the mortar variable from the previous - # time step will add a contribution to the rhs due to Backward Euler: - # \partial div_u / \partial_t = (\div_u^k - \div_u^{k-1})/dt. - rhs_slave = np.empty(3, dtype=np.object) - rhs_slave[master_ind] = np.zeros(matrix[master_ind, master_ind].shape[1]) - rhs_slave[slave_ind] = np.zeros(matrix[slave_ind, slave_ind].shape[1]) - rhs_slave[2] = np.zeros(self.ndof(mg)) - # I got some problems with pointers when doing rhs_master = rhs_slave.copy() - # so just reconstruct everything. - rhs_master = np.empty(3, dtype=np.object) - rhs_master[master_ind] = np.zeros(matrix[master_ind, master_ind].shape[1]) - rhs_master[slave_ind] = np.zeros(matrix[slave_ind, slave_ind].shape[1]) - rhs_master[2] = np.zeros(self.ndof(mg)) - - # The convention, for now, is to put the master grid information - # in the first column and row in matrix, slave grid in the second - # and mortar variables in the third - # If master and slave is the same grid, they should contribute to the same - # row and coloumn. When the assembler assigns matrix[idx] it will only add - # the slave information because of duplicate indices (master and slave is the same). - # We therefore write the both master and slave info to the slave index. - if g_master == g_slave: - master_ind = 1 - else: - master_ind = 0 - - # lambda acts as a boundary condition on the div_u term. Assemble it for the master. - self.discr_master.assemble_int_bound_stress( - g_master, - data_master, - data_edge, - False, - cc_master, - matrix, - rhs_master, - master_ind, - ) - # Equivalent for the slave - self.discr_slave.assemble_int_bound_stress( - g_slave, data_slave, data_edge, True, cc_slave, matrix, rhs_slave, slave_ind - ) - # We now have to flip the sign of some of the matrices - # First we flip the sign of the master stress because the mortar stress - # is defined from the slave stress. Then, stress_master = -\lambda - cc_master[master_ind, 2] = -cc_master[master_ind, 2] - rhs_master[master_ind] = -rhs_master[master_ind] - - matrix += cc_master + cc_slave - rhs = [s + m for s, m in zip(rhs_slave, rhs_master)] - - # The following two functions might or might not be needed when using - # a finite element discretization (see RobinCoupling for flow). - self.discr_master.enforce_neumann_int_bound( - g_master, data_edge, matrix, False, master_ind - ) - self.discr_slave.enforce_neumann_int_bound( - g_slave, data_edge, matrix, True, slave_ind - ) - - return matrix, rhs diff --git a/src/porepy/numerics/interface_laws/hyperbolic_interface_laws.py b/src/porepy/numerics/interface_laws/hyperbolic_interface_laws.py index 6960a67259..63c09f20df 100644 --- a/src/porepy/numerics/interface_laws/hyperbolic_interface_laws.py +++ b/src/porepy/numerics/interface_laws/hyperbolic_interface_laws.py @@ -80,7 +80,7 @@ def assemble_matrix_rhs( not_flag = 1 - flag # assemble matrices - # Transport out off upper equals lambda + # Transport out of upper equals lambda cc[0, 2] = div * hat_P_avg.T # transport out of lower is -lambda @@ -108,7 +108,8 @@ def assemble_matrix_rhs( # rhs is zero rhs = np.squeeze([np.zeros(dof[0]), np.zeros(dof[1]), np.zeros(dof[2])]) - return matrix + cc, rhs + matrix += cc + return matrix, rhs def cfl( self, diff --git a/src/porepy/numerics/mixed_dim/assembler.py b/src/porepy/numerics/mixed_dim/assembler.py index e47abe29c0..a3476f814c 100644 --- a/src/porepy/numerics/mixed_dim/assembler.py +++ b/src/porepy/numerics/mixed_dim/assembler.py @@ -6,7 +6,6 @@ """ import numpy as np import scipy.sparse as sps - import porepy as pp @@ -150,14 +149,6 @@ def assemble_matrix_rhs(self, matrix_format="csr", add_matrices=True): else: sps_matrix = sps.csr_matrix - # Initialize the global matrix. - # This gives us a set of matrices (essentially one per term per variable) - # and a simial set of rhs vectors. Furthermore, we get block indices - # of variables on individual nodes and edges, and count the number of - # dofs per local variable. - # For details, and some nuances, see documentation of the funciton - # _initialize_matrix_rhs. - matrix, rhs = self._initialize_matrix_rhs(sps_matrix) # If there are no variables - most likely if the active_variables do not # match any of the decleared variables, we can return now. if len(self.full_dof) == 0: @@ -167,7 +158,126 @@ def assemble_matrix_rhs(self, matrix_format="csr", add_matrices=True): mat, vec = self._assign_matrix_vector(self.full_dof, sps_matrix) return mat, vec else: - return matrix, rhs + return self._initialize_matrix_rhs(sps_matrix) + + # Assemble + matrix, rhs = self._operate_on_gb("assemble", matrix_format=matrix_format) + + # At this stage, all assembly is done. The remaining step is optionally to + # add the matrices associated with different terms, and anyhow convert + # the matrix to a sps. block matrix. + if add_matrices: + size = np.sum(self.full_dof) + full_matrix = sps_matrix((size, size)) + full_rhs = np.zeros(size) + + for mat in matrix.values(): + full_matrix += sps.bmat(mat, matrix_format) + + for vec in rhs.values(): + full_rhs += np.concatenate(tuple(vec)) + + return full_matrix, full_rhs + + else: + for k, v in matrix.items(): + matrix[k] = sps.bmat(v, matrix_format) + for k, v in rhs.items(): + rhs[k] = np.concatenate(tuple(v)) + + return matrix, rhs + + def discretize(self, variable_filter=None, term_filter=None): + """ Run the discretization operation on discretizations specified in + the mixed-dimensional grid. + + Only active variables will be considered. Moreover, the discretization + operation can be filtered to only consider specified variables, or terms. + If the variable filter is active, only discretizations where all variables + survive the filter will be discretized (for diagonal terms, the variable + must survive, for off-diagonal terms, both terms must survive). + + Filtering on terms works on the more detailed levels of indivdiual terms + in a multi-physics discretization (say, zoom-in on the advection term + in a advection-diffusion system). + + The filters can be combined to select specified terms for specified equations. + + Example (discretization internal to a node or edge: + For a discretizaiton of the form + + data[pp.DISCRETIZATION] = {'temp': {'advection': Foo(), 'diffusion': Bar()}, + 'pressure' : {'diffusion': FlowFoo()}} + + variable_filter = ['temp'] will discretize all temp terms + + term_filter = ['diffusion'] will discretize duffusion for both the temp and + pressure variable + + variable_filter = ['temp'], term_filter = ['diffusion'] will only discretize + the diffusion term for variable temp + + Example (coupling terms): + Variable filter works as intenal to nodes / edges. + The term filter acts on the identifier of a coupling, so + + dd[[pp.COUPLING_DISCRETIZATION]] = {'coupling_id' : {g1: {'temp': 'diffusion'}, + g2: {'pressure': diffusion'}, + (g1, g2): {'coupling_variable': FooBar()}}} + + will survive term_filter = ['coupling_id'] + + + Parameters: + variable_filter (optional): List of variables to be discretized. If + None (default), all active variables are discretized. + term_filter (optional): List of terms to be discretized. If None + (default), all terms for all active variables are discretized. + + """ + self._operate_on_gb( + "discretize", variable_filter=variable_filter, term_filter=term_filter + ) + + def _operate_on_gb(self, operation, **kwargs): + """ Helper method, loop over the GridBucket, identify nodes / edges + variables and discretizations, and perform an operation on these. + + Implemented actions are discretizaiton and assembly. + + """ + + if operation == "discretize": + variable_keys = kwargs.get("variable_filter", None) + if variable_keys is None: + variable_filter = lambda x: True + else: + variable_filter = lambda x: x in variable_keys + term_keys = kwargs.get("term_filter", None) + if term_keys is None: + term_filter = lambda x: True + else: + term_filter = lambda x: x in term_keys + elif operation == "assemble": + # Initialize the global matrix. + # This gives us a set of matrices (essentially one per term per variable) + # and a simial set of rhs vectors. Furthermore, we get block indices + # of variables on individual nodes and edges, and count the number of + # dofs per local variable. + # For details, and some nuances, see documentation of the funciton + # _initialize_matrix_rhs. + matrix_format = kwargs.get("matrix_format", "csc") + if matrix_format == "csc": + sps_matrix = sps.csc_matrix + else: + sps_matrix = sps.csr_matrix + + matrix, rhs = self._initialize_matrix_rhs(sps_matrix) + + else: + # We will only reach this if someone has invoked this private method + # from the outside. + raise ValueError("Unknown gb operation " + str(operation)) # Loop over all grids, discretize (if necessary) and assemble. This # will populate the main diagonal of the equation. @@ -192,23 +302,32 @@ def assemble_matrix_rhs(self, matrix_format="csr", add_matrices=True): # we should do something. # Loop over all discretizations for term, d in discr.items(): - # Assemble the matrix and right hand side. This will also - # discretize if not done before. - loc_A, loc_b = d.assemble_matrix_rhs(g, data) - - # Assign values in global matrix: Create the same key used - # defined when initializing matrices (see that function) - var_key_name = self._variable_term_key(term, row, col) - - # Check if the current block is None or not, it could - # happend based on the problem setting. Better to stay - # on the safe side. - if matrix[var_key_name][ri, ci] is None: - matrix[var_key_name][ri, ci] = loc_A - else: - matrix[var_key_name][ri, ci] += loc_A - # The right hand side vector is always initialized. - rhs[var_key_name][ri] += loc_b + + if operation == "discretize": + if ( + variable_filter(row) + and variable_filter(col) + and term_filter(term) + ): + d.discretize(g, data) + elif operation == "assemble": + # Assemble the matrix and right hand side. This will also + # discretize if not done before. + loc_A, loc_b = d.assemble_matrix_rhs(g, data) + + # Assign values in global matrix: Create the same key used + # defined when initializing matrices (see that function) + var_key_name = self._variable_term_key(term, row, col) + + # Check if the current block is None or not, it could + # happend based on the problem setting. Better to stay + # on the safe side. + if matrix[var_key_name][ri, ci] is None: + matrix[var_key_name][ri, ci] = loc_A + else: + matrix[var_key_name][ri, ci] += loc_A + # The right hand side vector is always initialized. + rhs[var_key_name][ri] += loc_b # Loop over all edges for e, data_edge in self.gb.edges(): @@ -239,20 +358,29 @@ def assemble_matrix_rhs(self, matrix_format="csr", add_matrices=True): else: # Loop over all discretizations for term, d in discr.items(): - # Assemble the matrix and right hand side. This will also - # discretize if not done before. - loc_A, loc_b = d.assemble_matrix_rhs(g, data_edge) - - # Assign values in global matrix - var_key_name = self._variable_term_key(term, row, col) - # Check if the current block is None or not, it could - # happend based on the problem setting. Better to stay - # on the safe side. - if matrix[var_key_name][ri, ci] is None: - matrix[var_key_name][ri, ci] = loc_A - else: - matrix[var_key_name][ri, ci] += loc_A - rhs[var_key_name][ri] += loc_b + if operation == "discretize": + if ( + variable_filter(row) + and variable_filter(col) + and term_filter(term) + ): + d.discretize(g, data) + elif operation == "assemble": + # Assemble the matrix and right hand side. This will also + # discretize if not done before. + + loc_A, loc_b = d.assemble_matrix_rhs(g, data_edge) + + # Assign values in global matrix + var_key_name = self._variable_term_key(term, row, col) + # Check if the current block is None or not, it could + # happend based on the problem setting. Better to stay + # on the safe side. + if matrix[var_key_name][ri, ci] is None: + matrix[var_key_name][ri, ci] = loc_A + else: + matrix[var_key_name][ri, ci] += loc_A + rhs[var_key_name][ri] += loc_b # Then, discretize the interaction between the edge variables of # this edge, and the adjacent node variables. @@ -260,7 +388,7 @@ def assemble_matrix_rhs(self, matrix_format="csr", add_matrices=True): if discr is None: continue - for key, terms in discr.items(): + for coupling_key, terms in discr.items(): edge_vals = terms.get(e) edge_key = edge_vals[0] @@ -320,7 +448,9 @@ def assemble_matrix_rhs(self, matrix_format="csr", add_matrices=True): # Key to the matrix dictionary used to access this coupling # discretization. - mat_key = self._variable_term_key(key, edge_key, slave_key, master_key) + mat_key = self._variable_term_key( + coupling_key, edge_key, slave_key, master_key + ) # Edge discretization object e_discr = edge_vals[1] @@ -331,107 +461,122 @@ def assemble_matrix_rhs(self, matrix_format="csr", add_matrices=True): # used. The fourth alternative, none of them are active, is not # considered valid, and raises an error message. if mi is not None and si is not None: + if operation == "discretize": + if ( + variable_filter(master_key) + and variable_filter(slave_key) + and variable_filter(edge_key) + ): + e_discr.discretize( + g_master, g_slave, data_master, data_slave, data_edge + ) + + elif operation == "assemble": + + # Assign a local matrix, which will be populated with the + # current state of the local system. + # Local here refers to the variable and term on the two + # nodes, together with the relavant mortar variable and term + # Associate the first variable with master, the second with + # slave, and the final with edge. + loc_mat, _ = self._assign_matrix_vector( + self.full_dof[[mi, si, ei]], sps_matrix + ) - # Assign a local matrix, which will be populated with the - # current state of the local system. - # Local here refers to the variable and term on the two - # nodes, together with the relavant mortar variable and term - # Associate the first variable with master, the second with - # slave, and the final with edge. - loc_mat, _ = self._assign_matrix_vector( - self.full_dof[[mi, si, ei]], sps_matrix - ) - - # Pick out the discretizations on the master and slave node - # for the relevant variables. - # There should be no contribution or modification of the - # [0, 1] and [1, 0] terms, since the variables are only - # allowed to communicate via the edges. - loc_mat[0, 0] = matrix[mat_key_master][mi, mi] - loc_mat[1, 1] = matrix[mat_key_slave][si, si] - - # Run the discretization, and assign the resulting matrix - # to a temporary construct - tmp_mat, loc_rhs = e_discr.assemble_matrix_rhs( - g_master, g_slave, data_master, data_slave, data_edge, loc_mat - ) - # The edge column and row should be assigned to mat_key - matrix[mat_key][(ei), (mi, si, ei)] = tmp_mat[(2), (0, 1, 2)] - matrix[mat_key][(mi, si), (ei)] = tmp_mat[(0, 1), (2)] - # Also update the discretization on the master and slave - # nodes - matrix[mat_key_master][mi, mi] = tmp_mat[0, 0] - matrix[mat_key_slave][si, si] = tmp_mat[1, 1] + # Pick out the discretizations on the master and slave node + # for the relevant variables. + # There should be no contribution or modification of the + # [0, 1] and [1, 0] terms, since the variables are only + # allowed to communicate via the edges. + loc_mat[0, 0] = matrix[mat_key_master][mi, mi] + loc_mat[1, 1] = matrix[mat_key_slave][si, si] + + # Run the discretization, and assign the resulting matrix + # to a temporary construct + tmp_mat, loc_rhs = e_discr.assemble_matrix_rhs( + g_master, + g_slave, + data_master, + data_slave, + data_edge, + loc_mat, + ) + # The edge column and row should be assigned to mat_key + matrix[mat_key][(ei), (mi, si, ei)] = tmp_mat[(2), (0, 1, 2)] + matrix[mat_key][(mi, si), (ei)] = tmp_mat[(0, 1), (2)] + # Also update the discretization on the master and slave + # nodes + matrix[mat_key_master][mi, mi] = tmp_mat[0, 0] + matrix[mat_key_slave][si, si] = tmp_mat[1, 1] - # Finally take care of the right hand side - rhs[mat_key][[mi, si, ei]] += loc_rhs + # Finally take care of the right hand side + rhs[mat_key][[mi, si, ei]] += loc_rhs elif mi is not None: # si is None # The operation is a simplified version of the full option above. - loc_mat, _ = self._assign_matrix_vector( - self.full_dof[[mi, ei]], sps_matrix - ) - loc_mat[0, 0] = matrix[mat_key_master][mi, mi] - tmp_mat, loc_rhs = e_discr.assemble_matrix_rhs( - g_master, data_master, data_edge, loc_mat - ) - matrix[mat_key][(ei), (mi, ei)] = tmp_mat[(1), (0, 1)] - matrix[mat_key][mi, ei] = tmp_mat[0, 1] + if operation == "discretize": + if ( + variable_filter(master_key) + and variable_filter(edge_key) + and term_filter(term) + ): + e_discr.discretize(g_master, data_master, data_edge) + elif operation == "assemble": + + loc_mat, _ = self._assign_matrix_vector( + self.full_dof[[mi, ei]], sps_matrix + ) + loc_mat[0, 0] = matrix[mat_key_master][mi, mi] + tmp_mat, loc_rhs = e_discr.assemble_matrix_rhs( + g_master, data_master, data_edge, loc_mat + ) + matrix[mat_key][(ei), (mi, ei)] = tmp_mat[(1), (0, 1)] + matrix[mat_key][mi, ei] = tmp_mat[0, 1] - # Also update the discretization on the master and slave - # nodes - matrix[mat_key_master][mi, mi] = tmp_mat[0, 0] + # Also update the discretization on the master and slave + # nodes + matrix[mat_key_master][mi, mi] = tmp_mat[0, 0] - rhs[mat_key][[mi, ei]] += loc_rhs + rhs[mat_key][[mi, ei]] += loc_rhs elif si is not None: # mi is None # The operation is a simplified version of the full option above. - loc_mat, _ = self._assign_matrix_vector( - self.full_dof[[si, ei]], sps_matrix - ) - loc_mat[0, 0] = matrix[mat_key_slave][si, si] - tmp_mat, loc_rhs = e_discr.assemble_matrix_rhs( - g_slave, data_slave, data_edge, loc_mat - ) - matrix[mat_key][ei, (si, ei)] = tmp_mat[1, (0, 1)] - matrix[mat_key][si, ei] = tmp_mat[0, 1] + if operation == "discretize": + if ( + variable_filter(slave_key) + and variable_filter(edge_key) + and term_filter(term) + ): + e_discr.discretize(g_slave, data_slave, data_edge) + elif operation == "assemble": + + loc_mat, _ = self._assign_matrix_vector( + self.full_dof[[si, ei]], sps_matrix + ) + loc_mat[0, 0] = matrix[mat_key_slave][si, si] + tmp_mat, loc_rhs = e_discr.assemble_matrix_rhs( + g_slave, data_slave, data_edge, loc_mat + ) + matrix[mat_key][ei, (si, ei)] = tmp_mat[1, (0, 1)] + matrix[mat_key][si, ei] = tmp_mat[0, 1] - # Also update the discretization on the master and slave - # nodes - matrix[mat_key_slave][si, si] = tmp_mat[0, 0] + # Also update the discretization on the master and slave + # nodes + matrix[mat_key_slave][si, si] = tmp_mat[0, 0] - rhs[mat_key][[si, ei]] += loc_rhs + rhs[mat_key][[si, ei]] += loc_rhs else: raise ValueError( "Invalid combination of variables on node-edge relation" ) - # At this stage, all assembly is done. The remaining step is optionally to - # add the matrices associated with different terms, and anyhow convert - # the matrix to a sps. block matrix. - if add_matrices: - size = np.sum(self.full_dof) - full_matrix = sps_matrix((size, size)) - full_rhs = np.zeros(size) - - for mat in matrix.values(): - full_matrix += sps.bmat(mat, matrix_format) - - for vec in rhs.values(): - full_rhs += np.concatenate(tuple(vec)) - - return full_matrix, full_rhs - - else: - for k, v in matrix.items(): - matrix[k] = sps.bmat(v, matrix_format) - for k, v in rhs.items(): - rhs[k] = np.concatenate(tuple(v)) - + if operation == "assemble": return matrix, rhs + else: + return None def _identify_dofs(self): """ @@ -865,10 +1010,11 @@ def distribute_variable(self, values, variable_names=None, use_state=True): name = pair[1] if name != var_name: continue - if isinstance(g, pp.Grid): - data = self.gb.node_props(g) - else: # This is really an edge + if isinstance(g, tuple): + # This is really an edge data = self.gb.edge_props(g) + else: + data = self.gb.node_props(g) if pp.STATE in data.keys(): data[pp.STATE][var_name] = values[dof[bi] : dof[bi + 1]] @@ -892,13 +1038,39 @@ def merge_variable(self, var): for pair, bi in self.block_dof.items(): g = pair[0] var_name = pair[1] - if isinstance(g, pp.Grid): - data = self.gb.node_props(g) - else: # This is really an edge + if isinstance(g, tuple): + # This is really an edge data = self.gb.edge_props(g) + else: + data = self.gb.node_props(g) if var_name == var: loc_value = data[pp.STATE][var_name] else: loc_value = 0 values[dof[bi] : dof[bi + 1]] = loc_value return values + + def dof_ind(self, g, name): + """ Get the indices in the global system of variables associated with a + given node / edge (in the GridBucket sense) and a given variable. + + Parameters: + g (pp.Grid or pp.GridBucket edge): Either a grid, or an edge in the + GridBucket. + name (str): Name of a variable. Should be an active variable. + + Returns: + np.array (int): Index of degrees of freedom for this variable. + + """ + block_ind = self.block_dof[(g, name)] + dof_start = np.hstack((0, np.cumsum(self.full_dof))) + return np.arange(dof_start[block_ind], dof_start[block_ind + 1]) + + def num_dof(self): + """ Get total number of unknowns of the identified variables. + + Returns: + int: Number of unknowns. Size of solution vector. + """ + return self.full_dof.sum() diff --git a/src/porepy/params/bc.py b/src/porepy/params/bc.py index 8a82fd26b5..c372fc68a1 100644 --- a/src/porepy/params/bc.py +++ b/src/porepy/params/bc.py @@ -41,9 +41,8 @@ class BoundaryCondition(AbstractBoundaryCondition): """ Class to store information on boundary conditions. - The BCs are specified by face number, and can have type Dirichlet or - Neumann (Robin may be included later). For details on default values etc., - see constructor. + The BCs are specified by face number, and can have type Dirichlet, Neumann + or Robin. For details on default values etc. see constructor. Attributes: num_faces (int): Number of faces in the grid @@ -55,11 +54,13 @@ 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 + face i has been assigned a Robin condition. """ def __init__(self, g, faces=None, cond=None): - """Constructor for BoundaryConditions. + """Constructor for BoundaryCondition. The conditions are specified by face numbers. Faces that do not get an explicit condition will have Neumann conditions assigned. @@ -68,7 +69,8 @@ def __init__(self, g, faces=None, cond=None): g (grid): For which boundary conditions are set. faces (np.ndarray): Faces for which conditions are assigned. cond (list of str): Conditions on the faces, in the same order as - used in faces. Should be as long as faces. + used in faces. Should be as long as faces. The list elements + should be one of "dir", "neu", "rob". Example: # Assign Dirichlet condititons on the left side of a grid; implicit @@ -262,9 +264,58 @@ class BoundaryConditionVectorial(AbstractBoundaryCondition): refer to the above class BoundaryCondition. NOTE: g.dim > 1 for the procedure to make sense + + Attributes: + num_faces (int): Number of faces in the grid + dim (int): Dimension of the boundary. One less than the dimension of + the grid. + is_neu (np.ndarray boolean, size g.dim x g.num_faces): Element i is true if + face i has been assigned a Neumann condition. Tacitly assumes that + the face is on the boundary. Should be false for internal faces, as + 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 + face i has been assigned a Robin condition. + """ def __init__(self, g, faces=None, cond=None): + """Constructor for BoundaryConditionVectorial. + + The conditions are specified by face numbers. Faces that do not get an + explicit condition will have Neumann conditions assigned. + + Parameters: + g (grid): For which boundary conditions are set. + faces (np.ndarray): Faces for which conditions are assigned. + cond (list of str): Conditions on the faces, in the same order as + used in faces. Should be as long as faces. To set uniform condition + in all spatial directions for a face, use 'dir', 'neu', or 'rob'. + + NOTE: For more general combinations of boundary conditions, it is + recommended to first construct a BoundaryConditionVectorial object, + and then access the attributes is_dir, is_neu, is_rob to set the + conditions. + + Example: + # Assign Dirichlet condititons on the left side of a grid; implicit + # Neumann conditions on the rest + g = pp.CartGrid([2, 2]) + west_face = pp.bc.face_on_side(g, 'west') + bound_cond = pp.BoundaryConditionVectorial(g, faces=west_face, cond=['dir', + 'dir']) + + Example: + Assign Dirichlet condition in the x-direction, Robin in the z-direction. + g = pp.CartGrid([2, 2, 2]) + bc = pp.BoundaryConditionVectorial(g) + target_face = 0 + bc.is_neu[[0, 2], target_face] = False + bc.is_dir[0, target_face] = True + bc.is_rob[2, target_face] = True + + """ self.num_faces = g.num_faces self.dim = g.dim @@ -351,7 +402,7 @@ def set_bc(self, faces, cond): self.is_neu[1, faces[j]] = True self.is_neu[2, faces[j]] = False else: - raise ValueError("Boundary should be Dirichlet or Neumann") + raise ValueError(f"Unknown boundary condition {s}") def face_on_side(g, side, tol=1e-8): diff --git a/src/porepy/params/data.py b/src/porepy/params/data.py index e745d50376..0ad05da331 100644 --- a/src/porepy/params/data.py +++ b/src/porepy/params/data.py @@ -246,8 +246,8 @@ def initialize_data(g, data, keyword, specified_parameters=None): in data, the new keyword is added using the update_dictionaries method. Args: + g: The grid. Can be either standard grid, or mortar grid. data: Outer data dictionary, to which the parameters will be added. - g: The grid. keyword: String identifying the parameters. specified_parameters: A dictionary with specified parameters, defaults to empty dictionary. diff --git a/src/porepy/params/rock.py b/src/porepy/params/rock.py index c3c5b87c9e..d9e1fdb854 100644 --- a/src/porepy/params/rock.py +++ b/src/porepy/params/rock.py @@ -3,6 +3,8 @@ Contains standard values (e.g. found in Wikipedia) for permeability, elastic moduli etc. +Note that thermal expansion coefficients are linear (m/mK) for rocks, but +volumetric (m^3/m^3) for fluids. """ import porepy as pp @@ -157,7 +159,7 @@ def __init__(self, theta_ref=None): self.PERMEABILITY = 1e-8 * pp.DARCY self.POROSITY = 0.01 # Reported range for Young's modulus by jsg is 10-70GPa - self.YOUNG_MODULUS = 4.0 * pp.GIGA * pp.PASCAL + self.YOUNG_MODULUS = 40.0 * pp.GIGA * pp.PASCAL # Reported range for Poisson's ratio is 0.125-0.25 self.POISSON_RATIO = 0.2 @@ -166,9 +168,8 @@ def __init__(self, theta_ref=None): self.LAMBDA, self.MU = lame_from_young_poisson( self.YOUNG_MODULUS, self.POISSON_RATIO ) - self.THERMAL_EXPANSION = ( - 8e-6 * pp.METER / (pp.METER * pp.CELSIUS) - ) # from engineeringtoolbox.com + # Units of thermal expansion: m^3 / m^3 K, i.e. volumetric. From engineeringtoolbox.com + self.THERMAL_EXPANSION = 8e-6 * pp.METER / (pp.METER * pp.CELSIUS) if theta_ref is None: self.theta_ref = 20.0 * pp.CELSIUS else: diff --git a/src/porepy/params/tensor.py b/src/porepy/params/tensor.py index 219193af4d..2726c8fbc4 100644 --- a/src/porepy/params/tensor.py +++ b/src/porepy/params/tensor.py @@ -183,7 +183,7 @@ def __init__(self, dim, mu, lmbda, phi=None): Parameters ---------- - dim (int) dimension, should be 2 or 3 + dim (int) dimension, should be 1, 2 or 3 mu (numpy.ndarray), First lame parameter, 1-D, one value per cell lmbda (numpy.ndarray), Second lame parameter, 1-D, one value per cell phi (Optional numpy.ndarray), 1-D one value per cell, never been used. @@ -191,7 +191,7 @@ def __init__(self, dim, mu, lmbda, phi=None): """ # Check arguments - if dim > 3 or dim < 2: + if dim > 3 or dim < 1: raise ValueError("Dimension should be between 1 and 3") if not isinstance(mu, np.ndarray): @@ -220,6 +220,49 @@ def __init__(self, dim, mu, lmbda, phi=None): self.dim = dim # Basis for the contributions of mu, lmbda and phi is hard-coded + if dim == 1: + mu_mat = np.array( + [ + [2, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1], + ] + ) + + lmbda_mat = np.array( + [ + [1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1], + ] + ) + + phi_mat = np.array( + [ + [0, 1, 1, 1, 0, 1, 1, 1, 0], + [1, 0, 0, 0, 1, 0, 0, 0, 1], + [1, 0, 0, 0, 1, 0, 0, 0, 1], + [1, 0, 0, 0, 1, 0, 0, 0, 1], + [0, 1, 1, 1, 0, 1, 1, 1, 0], + [1, 0, 0, 0, 1, 0, 0, 0, 1], + [1, 0, 0, 0, 1, 0, 0, 0, 1], + [1, 0, 0, 0, 1, 0, 0, 0, 1], + [0, 1, 1, 1, 0, 1, 1, 1, 0], + ] + ) + if dim == 2: mu_mat = np.array( [ diff --git a/src/porepy/params/water.py b/src/porepy/params/water.py index b648830a0a..4e5e786da1 100644 --- a/src/porepy/params/water.py +++ b/src/porepy/params/water.py @@ -1,3 +1,10 @@ +""" Hard coded typical parameters that may be of use in simulations. + +Contains standard values (e.g. found in Wikipedia) for density, thermal properties etc. + +Note that thermal expansion coefficients are linear (m/mK) for rocks, but +volumetric (m^3/m^3) for fluids. +""" import numpy as np import porepy as pp @@ -9,7 +16,11 @@ def __init__(self, theta_ref=None): else: self.theta_ref = theta_ref + self.COMPRESSIBILITY = 4e-10 / pp.PASCAL # Moderate dependency on theta + self.BULK = 1 / self.COMPRESSIBILITY + def thermal_expansion(self, delta_theta): + """ Units: m^3 / m^3 K, i.e. volumetric """ return ( 0.0002115 + 1.32 * 1e-6 * delta_theta @@ -17,6 +28,7 @@ def thermal_expansion(self, delta_theta): ) def density(self, theta=None): # theta in CELSIUS + """ Units: kg / m^3 """ if theta is None: theta = self.theta_ref theta_0 = 10 * (pp.CELSIUS) @@ -24,6 +36,7 @@ def density(self, theta=None): # theta in CELSIUS return rho_0 / (1.0 + self.thermal_expansion(theta - theta_0)) def thermal_conductivity(self, theta=None): # theta in CELSIUS + """ Units: W / m K """ if theta is None: theta = self.theta_ref return ( @@ -34,11 +47,13 @@ def thermal_conductivity(self, theta=None): # theta in CELSIUS ) def specific_heat_capacity(self, theta=None): # theta in CELSIUS + """ Units: J / kg K """ if theta is None: theta = self.theta_ref - return (4245 - 1.841 * theta) / self.density(theta) + return 4245 - 1.841 * theta def dynamic_viscosity(self, theta=None): # theta in CELSIUS + """Units: Pa s""" if theta is None: theta = self.theta_ref theta = pp.CELSIUS_to_KELVIN(theta) diff --git a/src/porepy/utils/derived_discretizations/implicit_euler.py b/src/porepy/utils/derived_discretizations/implicit_euler.py new file mode 100644 index 0000000000..d872f3364b --- /dev/null +++ b/src/porepy/utils/derived_discretizations/implicit_euler.py @@ -0,0 +1,243 @@ +""" +Module for extending the Upwind, MPFA and MassMatrix discretizations in Porepy to handle +implicit Euler time-stepping. Flux terms are multiplied by time step and the mass term +has a rhs contribution from the previous time step. +See the parent discretizations for further documentation. +""" +import porepy as pp +import numpy as np +import scipy.sparse as sps + + +class ImplicitMassMatrix(pp.MassMatrix): + """ + Return rhs contribution based on the previous solution, which is stored in the + pp.STATE field of the data dictionary. + """ + + def __init__(self, keyword="flow", variable="pressure"): + """ Set the discretization, with the keyword used for storing various + information associated with the discretization. The time discretisation also + requires the previous solution, thus the variable needs to be specified. + + Paramemeters: + keyword (str): Identifier of all information used for this + discretization. + """ + super().__init__(keyword) + self.variable = variable + + def assemble_rhs(self, g, data): + """ Overwrite MassMatrix method to return the correct rhs for an IE time + discretization, e.g. of the Biot problem. + """ + matrix_dictionary = data[pp.DISCRETIZATION_MATRICES][self.keyword] + previous_solution = data[pp.STATE][self.variable] + + return matrix_dictionary["mass"] * previous_solution + + +class ImplicitMpfa(pp.Mpfa): + """ + Multiply all contributions by the time step. + """ + + def assemble_matrix_rhs(self, g, data): + """ Overwrite MPFA method to be consistent with the Biot dt convention. + """ + a, b = super().assemble_matrix_rhs(g, data) + dt = data[pp.PARAMETERS][self.keyword]["time_step"] + a = a * dt + b = b * dt + return a, b + + def assemble_int_bound_flux( + self, g, data, data_edge, grid_swap, cc, matrix, rhs, self_ind + ): + """ + Overwrite the MPFA method to be consistent with the Biot dt convention + """ + dt = data[pp.PARAMETERS][self.keyword]["time_step"] + + div = g.cell_faces.T + + bound_flux = data[pp.DISCRETIZATION_MATRICES][self.keyword]["bound_flux"] + # Projection operators to grid + mg = data_edge["mortar_grid"] + + if grid_swap: + proj = mg.mortar_to_slave_int() + else: + proj = mg.mortar_to_master_int() + + if g.dim > 0 and bound_flux.shape[0] != g.num_faces: + # If bound flux is gven as sub-faces we have to map it from sub-faces + # to faces + hf2f = pp.fvutils.map_hf_2_f(nd=1, g=g) + bound_flux = hf2f * bound_flux + if g.dim > 0 and bound_flux.shape[1] != proj.shape[0]: + raise ValueError( + """Inconsistent shapes. Did you define a + sub-face boundary condition but only a face-wise mortar?""" + ) + + cc[self_ind, 2] += dt * div * bound_flux * proj + + def assemble_int_bound_source( + self, g, data, data_edge, grid_swap, cc, matrix, rhs, self_ind + ): + """ Abstract method. Assemble the contribution from an internal + boundary, manifested as a source term. + + The intended use is when the internal boundary is coupled to another + node in a mixed-dimensional method. Specific usage depends on the + interface condition between the nodes; this method will typically be + used to impose flux continuity on a lower-dimensional domain. + + Implementations of this method will use an interplay between the grid on + the node and the mortar grid on the relevant edge. + + Parameters: + g (Grid): Grid which the condition should be imposed on. + data (dictionary): Data dictionary for the node in the + mixed-dimensional grid. + data_edge (dictionary): Data dictionary for the edge in the + mixed-dimensional grid. + grid_swap (boolean): If True, the grid g is identified with the @ + slave side of the mortar grid in data_adge. + cc (block matrix, 3x3): Block matrix for the coupling condition. + The first and second rows and columns are identified with the + master and slave side; the third belongs to the edge variable. + The discretization of the relevant term is done in-place in cc. + matrix (block matrix 3x3): Discretization matrix for the edge and + the two adjacent nodes. + rhs (block_array 3x1): Right hand side contribution for the edge and + the two adjacent nodes. + self_ind (int): Index in cc and matrix associated with this node. + Should be either 1 or 2. + + """ + mg = data_edge["mortar_grid"] + + if grid_swap: + proj = mg.mortar_to_master_int() + else: + proj = mg.mortar_to_slave_int() + dt = data[pp.PARAMETERS][self.keyword]["time_step"] + cc[self_ind, 2] -= proj * dt + + +class ImplicitUpwind(pp.Upwind): + """ + Multiply all contributions by the time step and advection weight. + """ + + 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]) + + 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) + a = a * w + b = b * w + return a, b + + +class ImplicitUpwindCoupling(pp.UpwindCoupling): + """ + Multiply the advective mortar fluxes by the time step and advection weight. + """ + + def assemble_matrix_rhs( + self, g_master, g_slave, data_master, data_slave, data_edge, matrix + ): + """ + Construct the matrix (and right-hand side) for the coupling conditions. + Note: the right-hand side is not implemented now. + + Parameters: + matrix: Uncoupled discretization matrix. + g_master: grid of higher dimension + g_slave: grid of lower dimension + data_master: dictionary which stores the data for the higher dimensional + grid + data_slave: dictionary which stores the data for the lower dimensional + grid + data: dictionary which stores the data for the edges of the grid + bucket + + Returns: + cc: block matrix which store the contribution of the coupling + condition. See the abstract coupling class for a more detailed + description. + """ + + # Normal component of the velocity from the higher dimensional grid + + # @ALL: This should perhaps be defined by a globalized keyword + parameter_dictionary_master = data_master[pp.PARAMETERS][self.keyword] + lam_flux = data_edge[pp.PARAMETERS][self.keyword]["darcy_flux"] + dt = parameter_dictionary_master["time_step"] + w = parameter_dictionary_master["advection_weight"] * dt + # Retrieve the number of degrees of both grids + # Create the block matrix for the contributions + g_m = data_edge["mortar_grid"] + + # We know the number of dofs from the master and slave side from their + # discretizations + dof = np.array([matrix[0, 0].shape[1], matrix[1, 1].shape[1], g_m.num_cells]) + cc = np.array([sps.coo_matrix((i, j)) for i in dof for j in dof]) + cc = cc.reshape((3, 3)) + + # Projection from mortar to upper dimenional faces + hat_P_avg = g_m.master_to_mortar_avg() + # Projection from mortar to lower dimensional cells + check_P_avg = g_m.slave_to_mortar_avg() + + # mapping from upper dim cellls to faces + # The mortars always points from upper to lower, so we don't flip any + # signs + div = np.abs(pp.numerics.fv.fvutils.scalar_divergence(g_master)) + + # Find upwind weighting. if flag is True we use the upper weights + # if flag is False we use the lower weighs + flag = (lam_flux > 0).astype(np.float) + not_flag = 1 - flag + + # assemble matrices + # Transport out of upper equals lambda + cc[0, 2] = div * hat_P_avg.T + + # transport out of lower is -lambda + cc[1, 2] = -check_P_avg.T + + # Discretisation of mortars + # CHANGE from UpwindCoupling: multiply the discretization of the advective + # mortar fluxes by dt and advection weight (e.g. heat capacity) + + # If fluid flux(lam_flux) is positive we use the upper value as weight, + # i.e., T_masterat * fluid_flux = lambda. + # We set cc[2, 0] = T_masterat * fluid_flux + cc[2, 0] = sps.diags(w * lam_flux * flag) * hat_P_avg * div.T + + # If fluid flux is negative we use the lower value as weight, + # i.e., T_check * fluid_flux = lambda. + # we set cc[2, 1] = T_check * fluid_flux + cc[2, 1] = sps.diags(w * lam_flux * not_flag) * check_P_avg + + # The rhs of T * fluid_flux = lambda + # Recover the information for the grid-grid mapping + cc[2, 2] = -sps.eye(g_m.num_cells) + + if data_master["node_number"] == data_slave["node_number"]: + # All contributions to be returned to the same block of the + # global matrix in this case + cc = np.array([np.sum(cc, axis=(0, 1))]) + + # rhs is zero + rhs = np.squeeze([np.zeros(dof[0]), np.zeros(dof[1]), np.zeros(dof[2])]) + matrix += cc + return matrix, rhs diff --git a/src/porepy/utils/grid_util.py b/src/porepy/utils/grid_util.py deleted file mode 100644 index fbcba6e99b..0000000000 --- a/src/porepy/utils/grid_util.py +++ /dev/null @@ -1,26 +0,0 @@ -import numpy as np -import scipy.sparse as sps - - -def sign_of_boundary_faces(g): - """ - returns the sign of boundary faces as defined by g.cell_faces. - Parameters: - g: (Grid Object) - faces: (ndarray) indices of faces that you want to know the sign for. The - faces must be boundary faces. - - Returns: - sgn: (ndarray) the sign of the faces - """ - faces = g.get_all_boundary_faces() - - IA = np.argsort(faces) - IC = np.argsort(IA) - - fi, _, sgn = sps.find(g.cell_faces[faces[IA], :]) - assert fi.size == faces.size, "sign of internal faces does not make sense" - I = np.argsort(fi) - sgn = sgn[I] - sgn = sgn[IC] - return sgn diff --git a/src/porepy/utils/grid_utils.py b/src/porepy/utils/grid_utils.py new file mode 100644 index 0000000000..c7f4d87515 --- /dev/null +++ b/src/porepy/utils/grid_utils.py @@ -0,0 +1,38 @@ +import numpy as np +import scipy.sparse as sps + + +def switch_sign_if_inwards_normal(g, nd, faces): + """ Construct a matrix that changes sign of quantities on faces with a + normal that points into the grid. + + Parameters: + g (pp.Grid): Grid. + nd (int): Number of quantities per face; this will for instance be the + number of components in a face-vector. + faces (np.array-like of ints): Index for which faces to be considered. Should only + contain boundary faces. + + Returns: + sps.dia_matrix: Diagonal matrix which switches the sign of faces if the + normal vector of the face points into the grid g. Faces not considered + will have a 0 diagonal term. If nd > 1, the first nd rows are associated + with the first face, then nd elements of the second face etc. + + """ + + faces = np.asarray(faces) + + # Find out whether the boundary faces have outwards pointing normal vectors + # Negative sign implies that the normal vector points inwards. + sgn = g.sign_of_faces(faces) + + # Create vector with the sign in the places of faces under consideration, + # zeros otherwise + sgn_mat = np.zeros(g.num_faces) + sgn_mat[faces] = sgn + # Duplicate the numbers, the operator is intended for vector quantities + sgn_mat = np.tile(sgn_mat, (nd, 1)).ravel(order="F") + + # Create the diagonal matrix. + return sps.dia_matrix((sgn_mat, 0), shape=(sgn_mat.size, sgn_mat.size)) diff --git a/src/porepy/utils/tangential_normal_projection.py b/src/porepy/utils/tangential_normal_projection.py new file mode 100644 index 0000000000..da351d78f6 --- /dev/null +++ b/src/porepy/utils/tangential_normal_projection.py @@ -0,0 +1,261 @@ +""" +Geometric projections related to the tangential and normal spaces of a set of +vectors. +""" + +import numpy as np +import scipy.sparse as sps + + +class TangentialNormalProjection: + """ Represent a set of projections into tangent and normal vectors. + + The spaces are defined by the normal vector (see __init__ documentation). + The basis for the tangential space is arbitrary (arbitrary direction in 2d, + rotation angle in 3d). The basis for the tangential is stored in the attribute + tangential_basis. + + Attributes: + num_vecs (int): Number of tangent/normal spaces represented by this object. + dim (int): Dimension of the ambient space. + tangential_basis (np.array, size: dim x dim-1 x num_vec): Basis vectors for the + tangential space. + projection (np.array, size dim x dim x num_vecs): Projection matrices onto the + tangential and normal space. The first dim-1 rows represent projection to the + tangential space, the final row is the normal component. + normal (np.array, size dim x num_vecs): Unit normal vectors. + + """ + + def __init__(self, normals, dim=None): + if dim is None: + dim = normals.shape[0] + + # Normalize vectors + normals = normals / np.linalg.norm(normals, axis=0) + + self.num_vecs = normals.shape[1] + self.dim = dim + + # Compute normal and tangential basis + basis, normal = self._decompose_vector(normals) + + basis = basis.reshape((dim, dim, self.num_vecs)) + self.tangential_basis = basis[:, :-1, :] + + # The projection is found by inverting the basis vectors + self.projection = self._invert_3d_matrix(basis) + self.normals = normal + + ## Methods for genertation of projection matrices + + def project_tangential_normal(self, num=None): + """ Define a projection matrix to decompose a matrix into tangential + and normal components. + + The intended usage is to decompose a grid-based vector variable into the + tangent and normal spaces of the grid, with the tacit understanding that there is + a single normal vector shared for all the cells (or faces) in the grid. + + The method can also create projection matrix based on unequal normal vectors. + One projection will be generated per column in self.normal. To activate + this behavior, set num=None. + + Parameters: + num (int, optional): Number of (equal) projections to be generated. + Will correspond to the number of cells / faces in the grid. + The projection matrix will have num * self.dim columns. If not + specified (default), one projection will be generated per vector in + self.normals. + NOTE: If self.num_vecs > 1, but num is not None, only the first + given normal vector will be used to generate the tangential space. + + Returns: + scipy.sparse.csc_matrix: Projection matrix, structure as a block + diagonal matrix, with block size dim x dim. + For each block, the first dim-1 rows projects onto the tangent + space, the final row projects onto the normal space. + size: ((self.dim * num) x (self.dim * num). If num is not None, + size: ((self.dim * num_vecs) x (self.dim * num_vecs) + + """ + if num is None: + return sps.block_diag( + [self.projection[:, :, i] for i in range(self.projection.shape[-1])] + ).tocsc() + else: + return sps.block_diag( + [self.projection[:, :, 0] for i in range(num)] + ).tocsc() + + def project_tangential(self, num=None): + """ Define a projection matrix of a specific size onto the tangent space. + + The intended usage is to project a grid-based vector variable onto the + tangent space of the grid, with the tacit understanding that there is + a single normal vector shared for all the cells (or faces) in the grid. + + The method can also create projection matrix based on unequal normal vectors. + One projection will be generated per column in self.normal. To activate + this behavior, set num=None. + + Parameters: + num (int, optional): Number of (equal) projections to be generated. + Will correspond to the number of cells / faces in the grid. + The projection matrix will have num * self.dim columns. If not + specified (default), one projection will be generated per vector in + self.normals. + NOTE: If self.num_vecs > 1, but num is not None, only the first + given normal vector will be used to generate the tangential space. + + Returns: + scipy.sparse.csc_matrix: Tangential projection matrix, structure as a block + diagonal matrix. The first (dim-1) x dim block projects onto the first + tangent space, etc. + size: ((self.dim - 1) * num) x (self.dim * num). If num is not None, + size: ((self.dim - 1) * num_vecs) x (self.dim * num_vecs) + + """ + # Find type and size of projection. + if num is None: + num = self.num_vecs + + size_proj = self.dim * num + + # Construct the full projection matrix - tangential and normal + full_projection = self.project_tangential_normal(num) + + # Generate restriction matrix to the tangential space only + rows = np.arange(num * (self.dim - 1)) + cols = np.setdiff1d( + np.arange(size_proj), np.arange(self.dim - 1, size_proj, self.dim) + ) + data = np.ones_like(rows) + remove_normal_components = sps.csc_matrix( + (data, (rows, cols)), shape=(rows.size, size_proj) + ) + + # Return the restricted matrix. + return remove_normal_components * full_projection + + def project_normal(self, num=None): + """ Define a projection matrix of a specific size onto the normal space. + + The intended usage is to project a grid-based vector variable onto the + normal space of the grid, with the tacit understanding that there is + a single normal vector shared for all the cells (or faces) in the grid. + + The method can also create projection matrix based on unequal normal vectors. + One projection will be generated per column in self.normal. To activate + this behavior, set num=None. + + Parameters: + num (int, optional): Number of (equal) projections to be generated. + Will correspond to the number of cells / faces in the grid. + The projection matrix will have num * self.dim columns. If not + specified (default), one projection will be generated per vector in + self.normals. + NOTE: If self.num_vecs > 1, but num is not None, only the first + given normal vector will be used to generate the normal space. + + Returns: + scipy.sparse.csc_matrix: Tangential projection matrix, structure as a block + diagonal matrix. The first 1 x dim block projects onto the first + tangent space, etc. + size: num x (self.dim * num). If num is not None. + size: num_vecs x (self.dim * num_vecs) els. + + """ + # Find mode and size of projection + if num is None: + num = self.num_vecs + + size_proj = self.dim * num + + # Generate full projection matrix + full_projection = self.project_tangential_normal(num) + + # Construct restriction matrix to normal space. + rows = np.arange(num) + cols = np.arange(self.dim - 1, size_proj, self.dim) + data = np.ones_like(rows) + remove_tangential_components = sps.csc_matrix( + (data, (rows, cols)), shape=(rows.size, size_proj) + ) + + # Return the restricted matrix + return remove_tangential_components * full_projection + + def local_projection(self, ind=None): + """ Get the local projection matrix (refe) + + Paremeters: + ind (int, optional): Index (referring to the order of the normal vectors + provided to __init__) of the basis to return. Defaults to the first one. + + Returns: + np.array (self.dim x self.dim): Local projection matrix. Multiplication + gives projection to the tangential space (first self.dim - 1 rows) + and normal space (last) + + """ + if ind is None: + ind = 0 + return self.projection[:, :, ind] + + ### Helper functions below + + def _decompose_vector(self, nc): + if self.dim == 3: + t1 = np.random.rand(self.dim, 1) * np.ones(self.num_vecs) + t2 = np.random.rand(self.dim, 1) * np.ones(self.num_vecs) + normal, tc1, tc2 = self._gram_schmidt(nc, t1, t2) + basis = np.hstack([tc1, tc2, normal]) + else: + t1 = np.random.rand(self.dim, 1) * np.ones(self.num_vecs) + normal, tc1 = self._gram_schmidt(nc, t1) + basis = np.hstack([tc1, normal]) + return basis, normal + + def _gram_schmidt(self, u1, u2, u3=None): + """ + Perform a Gram Schmidt procedure for the vectors u1, u2 and u3 to obtain a set of + orhtogonal vectors. + + Parameters: + u1: ndArray + u2: ndArray + u3: ndArray + + Returns: + u1': ndArray u1 / ||u1|| + u2': ndarray (u2 - u2*u1 * u1) / ||u2|| + u3': (optional) ndArray (u3 - u3*u2' - u3*u1')/||u3|| + """ + u1 = u1 / np.sqrt(np.sum(u1 ** 2, axis=0)) + + u2 = u2 - np.sum(u2 * u1, axis=0) * u1 + u2 = u2 / np.sqrt(np.sum(u2 ** 2, axis=0)) + + if u3 is None: + return u1, u2 + u3 = u3 - np.sum(u3 * u1, axis=0) * u1 - np.sum(u3 * u2, axis=0) * u2 + u3 = u3 / np.sqrt(np.sum(u3 ** 2, axis=0)) + + return u1, u2, u3 + + def _invert_3d_matrix(self, M): + """ + Find the inverse of the (m,m,k) 3D ndArray M. The inverse is intrepreted as the + 2d inverse of M[:, :, i] for i = 0...k + + Parameters: + M: (m, m, k) ndArray + + Returns: + M_inv: Inverse of M + """ + M_inv = np.zeros(M.shape) + for i in range(M.shape[-1]): + M_inv[:, :, i] = np.linalg.inv(M[:, :, i]) + return M_inv diff --git a/test/integration/test_biot.py b/test/integration/test_biot.py index 8564d493ea..a02006b4ea 100644 --- a/test/integration/test_biot.py +++ b/test/integration/test_biot.py @@ -65,38 +65,6 @@ def test_no_dynamics_2d(self): self.assertTrue(np.isclose(sol, np.zeros(g.num_cells * (g.dim + 1))).all()) - # def test_uniform_displacement(self): - # # Uniform displacement in mechanics (enforced by boundary conditions). - # # Constant pressure boundary conditions. - # g_list = setup_grids.setup_2d() - # for g in g_list: - # bound_faces = g.get_all_boundary_faces() - # bound = bc.BoundaryCondition(g, bound_faces.ravel('F'), - # ['dir'] * bound_faces.size) - # flux, bound_flux, div_flow = self.mpfa_discr(g, bound) - # - # a_flow = div_flow * flux - # - # stress, bound_stress, grad_p, div_d, \ - # stabilization, bound_div_d, div_mech = self.mpsa_discr(g, bound) - # - # a_mech = div_mech * stress - # - # a_biot = sps.bmat([[a_mech, grad_p], - # [div_d, a_flow + stabilization]]) - # - # const_bound_val_mech = 1 - # bval_mech = const_bound_val_mech * np.ones(g.num_faces * g.dim) - # bval_flow = np.ones(g.num_faces) - # rhs = np.hstack((-div_mech * bound_stress * bval_mech, - # div_flow * bound_flux * bval_flow\ - # + div_flow * bound_div_d * bval_mech)) - # sol = np.linalg.solve(a_biot.todense(), rhs) - # - # sz_mech = g.num_cells * g.dim - # self.assertTrue(np.isclose(sol[:sz_mech],) - # const_bound_val_mech * np.ones(sz_mech)).all() - def test_face_vector_to_scalar(self): # Test of function face_vector_to_scalar nf = 3 @@ -159,7 +127,7 @@ def test_assemble_biot(self): term_11_2: pp.BiotStabilization(kw_f), }, v_0 + "_" + v_1: {term_01: pp.GradP(kw_m)}, - v_1 + "_" + v_0: {term_10: pp.DivD(kw_m)}, + v_1 + "_" + v_0: {term_10: pp.DivU(kw_m)}, } # Assemble. Also discretizes the flow terms (fluid_mass and fluid_flux) general_assembler = pp.Assembler(gb) @@ -256,7 +224,7 @@ def test_assemble_biot_rhs_transient(self): term_11_2: pp.BiotStabilization(kw_f), }, v_0 + "_" + v_1: {term_01: pp.GradP(kw_m)}, - v_1 + "_" + v_0: {term_10: pp.DivD(kw_m)}, + v_1 + "_" + v_0: {term_10: pp.DivU(kw_m)}, } times = np.arange(5) diff --git a/test/integration/test_contact_mechanics.py b/test/integration/test_contact_mechanics.py new file mode 100644 index 0000000000..048b8fc4ce --- /dev/null +++ b/test/integration/test_contact_mechanics.py @@ -0,0 +1,167 @@ +""" +Various integration tests for contact mechanics. +""" +import numpy as np +import unittest + +import porepy as pp +import porepy.models.contact_mechanics_model as model + + +class TestContactMechanics(unittest.TestCase): + def _solve(self, setup): + model.run_mechanics(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] + + 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 + + def test_pull_top_positive_opening(self): + + setup = SetupContactMechanics(ux_south=0, uy_bottom=0, ux_north=0, uy_top=0.001) + + u_mortar, contact_force = 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)) + + def test_pull_bottom_positive_opening(self): + + setup = SetupContactMechanics( + ux_south=0, uy_bottom=-0.001, ux_north=0, uy_top=0 + ) + + u_mortar, contact_force = 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)) + + def test_push_top_zero_opening(self): + + setup = SetupContactMechanics( + ux_south=0, uy_bottom=0, ux_north=0, uy_top=-0.001 + ) + + 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) + + # Contact force in normal direction should be negative + self.assertTrue(np.all(contact_force[1] < 0)) + + def test_push_bottom_zero_opening(self): + + setup = SetupContactMechanics(ux_south=0, uy_bottom=0.001, ux_north=0, uy_top=0) + + 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) + + # Contact force in normal direction should be negative + self.assertTrue(np.all(contact_force[1] < 0)) + + +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_south = ux_south + self.uy_bottom = uy_bottom + self.ux_north = ux_north + self.uy_top = uy_top + + def create_grid(self): + """ + 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. + """ + rotate_fracture = getattr(self, "rotate_fracture", False) + 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 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) + 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__": + unittest.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..c4570f470d --- /dev/null +++ b/test/integration/test_contact_mechanics_biot.py @@ -0,0 +1,266 @@ +#!/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 unittest + +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): + """ + 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. + """ + rotate_fracture = getattr(self, "rotate_fracture", False) + 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_scalar(self, g): + if g.dim == self.Nd: + values = np.zeros(g.num_cells) + else: + values = self.scalar_source_value * np.ones(g.num_cells) + return values + + def bc_type_mechanics(self, g): + _, _, _, north, south, _, _ = self.domain_boundary_sides(g) + 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) + # Define boundary condition on faces + return pp.BoundaryCondition(g, north + south, "dir") + + def bc_values_mechanics(self, g): + # Set the boundary values + _, _, _, 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 + return values.ravel("F") + + +if __name__ == "__main__": + unittest.main() diff --git a/test/integration/test_elliptic_interface_laws.py b/test/integration/test_elliptic_interface_laws.py deleted file mode 100644 index 27be77d6e3..0000000000 --- a/test/integration/test_elliptic_interface_laws.py +++ /dev/null @@ -1,505 +0,0 @@ -""" -Module for testing the vector elliptic couplings in the interface_laws. -""" - - -import numpy as np -import scipy.sparse as sps -import unittest - -import porepy as pp - - -class TestTwoGridCoupling(unittest.TestCase): - """ - In this test we set up a coupling between two grids of dimension 2: - g_slave g_master - |-----| | |------| - | | x | | - |-----| | |------| - g_mortar - There is one cell per grid and they are coupled together by a single mortar - variable. - - We define a random Robin and Mortar weights and test if we recover the - condition on the interface. - """ - - def test_robin_coupling(self): - """ - Test a Robin condition on the interface - """ - self.kw = "mech" - gb, _ = define_gb() - mortar_weight = np.random.rand(gb.dim_max()) - robin_weight = np.random.rand(gb.dim_max()) - rhs = np.random.rand(gb.dim_max()) - self.assign_parameters(gb, mortar_weight, robin_weight, rhs) - self.assign_discretization(gb) - assembler = pp.Assembler(gb) - matrix, rhs = assembler.assemble_matrix_rhs() - u = sps.linalg.spsolve(matrix, rhs) - assembler.distribute_variable(u) - self.check_solution(gb) - - def test_continuity_coupling(self): - """ - Test a continuity condition on the interface. This is equivalent to - zero mortar weight and identity matrix for the robin weight. These - matrices are only used to check the solution. We check the solution - against a reference computed on a single grid including both g_s and g_m. - """ - self.kw = "mech" - gb, gb_full = define_gb() - # We assign weighs according to the condition. - mortar_weight = np.zeros(gb.dim_max()) - robin_weight = np.ones(gb.dim_max()) - rhs = np.zeros(gb.dim_max()) - # Assign data to coupling gb - self.assign_parameters(gb, mortar_weight, robin_weight, rhs) - self.assign_discretization(gb, robin=False) - # Assign data to mono gb - self.assign_parameters(gb_full, mortar_weight, robin_weight, rhs) - self.assign_discretization(gb_full, robin=False) - - assembler = pp.Assembler(gb) - matrix, rhs = assembler.assemble_matrix_rhs() - u = sps.linalg.spsolve(matrix, rhs) - assembler.distribute_variable(u) - self.check_solution(gb) - - assembler_full = pp.Assembler(gb_full) - matrix, rhs = assembler_full.assemble_matrix_rhs() - u_full = sps.linalg.spsolve(matrix, rhs) - # compare solutions - # We need to rearange the solutions because the ordering of the dofs are not the same - # Also, we don't have equality because the weak symmetry breaks when a cell has to many - # Neumann conditions (see comments in mpsa) - us = [] - ID = [] # to appease porpy 3.5 - for g, d in gb: - us.append(d[pp.STATE][self.kw]) - ID.append(g.grid_num - 1) - us = np.hstack([np.array(us)[ID].ravel()]) - IA = np.array([0, 1, 4, 5, 2, 3, 6, 7]) - sol = us[IA] - self.assertTrue(np.all(np.abs(sol - u_full) < 1e-4)) - - def check_solution(self, gb): - for e, d in gb.edges(): - mg = d["mortar_grid"] - gs, gm = gb.nodes_of_edge(e) - - ds = gb.node_props(gs) - dm = gb.node_props(gm) - - us = ds[pp.STATE][self.kw] - um = dm[pp.STATE][self.kw] - lam = d[pp.STATE][self.kw] - - bdcs = ds[pp.DISCRETIZATION_MATRICES][self.kw]["bound_displacement_cell"] - bdcm = dm[pp.DISCRETIZATION_MATRICES][self.kw]["bound_displacement_cell"] - bdfs = ds[pp.DISCRETIZATION_MATRICES][self.kw]["bound_displacement_face"] - bdfm = dm[pp.DISCRETIZATION_MATRICES][self.kw]["bound_displacement_face"] - - bc_val_s = ds[pp.PARAMETERS][self.kw]["bc_values"] - bc_val_m = dm[pp.PARAMETERS][self.kw]["bc_values"] - - RW = d[pp.PARAMETERS][self.kw]["robin_weight"] - MW = d[pp.PARAMETERS][self.kw]["mortar_weight"] - rhs = d[pp.PARAMETERS][self.kw]["robin_rhs"].reshape( - (gs.dim, -1), order="F" - ) - slv2mrt_nd = sps.kron(mg.slave_to_mortar_int(), sps.eye(gs.dim)).tocsr() - mstr2mrt_nd = sps.kron(mg.master_to_mortar_int(), sps.eye(gs.dim)).tocsr() - - hf2fs = pp.fvutils.map_hf_2_f(g=gs) / 2 - hf2fm = pp.fvutils.map_hf_2_f(g=gm) / 2 - jump_u = ( - slv2mrt_nd - * hf2fs - * (bdcs * us + bdfs * (bc_val_s + slv2mrt_nd.T * lam)) - - mstr2mrt_nd - * hf2fm - * (bdcm * um - bdfm * (bc_val_m + mstr2mrt_nd.T * lam)) - ).reshape((gs.dim, -1), order="F") - lam_nd = lam.reshape((gs.dim, -1), order="F") - for i in range(len(RW)): - rhs_robin = MW[i].dot(lam_nd[:, i]) + RW[i].dot(jump_u[:, i]) - self.assertTrue(np.allclose(rhs_robin, rhs[:, i])) - - def assign_discretization(self, gb, robin=True): - for _, d in gb: - d[pp.PRIMARY_VARIABLES] = {self.kw: {"cells": gb.dim_max()}} - d[pp.DISCRETIZATION] = {self.kw: {"mortar": pp.Mpsa(self.kw)}} - - if robin: - contact = pp.RobinContact(self.kw, pp.Mpsa(self.kw)) - else: - contact = pp.StressDisplacementContinuity(self.kw, pp.Mpsa(self.kw)) - for e, d in gb.edges(): - g1, g2 = gb.nodes_of_edge(e) - d[pp.PRIMARY_VARIABLES] = {self.kw: {"cells": gb.dim_max()}} - d[pp.COUPLING_DISCRETIZATION] = { - self.kw: { - g1: (self.kw, "mortar"), - g2: (self.kw, "mortar"), - e: (self.kw, contact), - } - } - - def assign_parameters(self, gb, mortar_weight, robin_weight, rhs): - for g, d in gb: - if g.grid_num == 1: - dir_faces = g.face_centers[0] < 1e-10 - elif g.grid_num == 2: - dir_faces = g.face_centers[0] > 2 - 1e-10 - elif g.grid_num == 3: - dir_faces = (g.face_centers[0] < 1e-10) + ( - g.face_centers[0] > 2 - 1e-10 - ) - - bc_val = np.zeros((g.dim, g.num_faces)) - bc_val[0, dir_faces] = 0.1 * g.face_centers[0, dir_faces] - bc = pp.BoundaryConditionVectorial(g, dir_faces, "dir") - C = pp.FourthOrderTensor( - gb.dim_max(), np.ones(g.num_cells), np.ones(g.num_cells) - ) - data = { - "bc": bc, - "bc_values": bc_val.ravel("F"), - "fourth_order_tensor": C, - "source": np.zeros(g.num_cells * g.dim), - "inverter": "python", - } - pp.initialize_data(g, d, self.kw, data) - - for _, d in gb.edges(): - mg = d["mortar_grid"] - MW = sps.diags(mortar_weight) - RW = sps.diags(robin_weight) - data = { - "mortar_weight": [MW] * mg.num_cells, - "robin_weight": [RW] * mg.num_cells, - "robin_rhs": np.tile(rhs, (mg.num_cells)), - } - pp.initialize_data(mg, d, self.kw, data) - - -class TestBiotTwoGridCoupling(unittest.TestCase): - """ - In this test we set up a coupling between two grids of dimension 2: - g_slave g_master - |-----| | |------| - | | x | | - |-----| | |------| - g_mortar - There is one cell per grid and they are coupled together by a single mortar - variable. - - We define a random Robin and Mortar weights and test if we recover the - condition on the interface. - """ - - def test_robin_coupling(self): - """ - Test a Robin condition on the interface - """ - self.kw = "mech" - self.kw_f = "flow" - gb, _ = define_gb() - mortar_weight = np.random.rand(gb.dim_max()) - robin_weight = np.random.rand(gb.dim_max()) - rhs = np.random.rand(gb.dim_max()) - self.assign_parameters(gb, mortar_weight, robin_weight, rhs) - self.assign_discretization(gb) - - # Discretize - for g, d in gb: - pp.Biot( - self.kw, self.kw_f, vector_variable="u", scalar_variable="p" - ).discretize(g, d) - - assembler = pp.Assembler(gb) - matrix, rhs = assembler.assemble_matrix_rhs() - u = sps.linalg.spsolve(matrix, rhs) - assembler.distribute_variable(u) - self.check_solution(gb) - - def test_continuity_coupling(self): - """ - Test a continuity condition on the interface. This is equivalent to - zero mortar weight and identity matrix for the robin weight. These - matrices are only used to check the solution. - """ - self.kw = "mech" - self.kw_f = "flow" - gb, gb_full = define_gb() - # We assign weighs according to the condition. - mortar_weight = np.zeros(gb.dim_max()) - robin_weight = np.ones(gb.dim_max()) - rhs = np.zeros(gb.dim_max()) - # Assign data to coupling gb - self.assign_parameters(gb, mortar_weight, robin_weight, rhs) - self.assign_discretization(gb) - # Assign data to mono gb - self.assign_parameters(gb_full, mortar_weight, robin_weight, rhs) - self.assign_discretization(gb_full) - - # Discretize - for g, d in gb: - pp.Biot(self.kw, self.kw_f).discretize(g, d) - for g, d in gb_full: - pp.Biot(self.kw, self.kw_f).discretize(g, d) - - # Assemble and solve - assembler = pp.Assembler(gb) - matrix, rhs = assembler.assemble_matrix_rhs() - u = sps.linalg.spsolve(matrix, rhs) - assembler.distribute_variable(u) - self.check_solution(gb) - - assembler_full = pp.Assembler(gb_full) - matrix, rhs = assembler_full.assemble_matrix_rhs() - u_full = sps.linalg.spsolve(matrix, rhs) - assembler_full.distribute_variable(u_full) - # Compare solutions: - # We need to rearange the solutions because the ordering of the dofs are not the - # same for the two grids. - us = [] - ps = [] - ID = [] - for g, d in gb: - us.append(d[pp.STATE]["u"]) - ps.append(d[pp.STATE]["p"]) - ID.append(g.grid_num - 1) - us = np.hstack([np.array(us)[ID].ravel()]) - ps = np.hstack([np.array(ps)[ID].ravel()]) - IA = np.array([0, 1, 4, 5, 2, 3, 6, 7]) - IAp = np.array([0, 2, 1, 3]) - sol = np.hstack([us[IA], ps[IAp]]) - - us = [] - ps = [] - for g, d in gb_full: - us.append(d[pp.STATE]["u"]) - ps.append(d[pp.STATE]["p"]) - sol_full = np.hstack([np.array(us), np.array(ps)]) - # Note, we don't have equality because the weak symmetry breaks when a cell has to many - # Neumann conditions (see comments in mpsa) - self.assertTrue(np.all(np.abs(sol - sol_full) < 5e-4)) - - def check_solution(self, gb): - for e, d in gb.edges(): - mg = d["mortar_grid"] - gs, gm = gb.nodes_of_edge(e) - - ds = gb.node_props(gs) - dm = gb.node_props(gm) - - us = ds[pp.STATE]["u"] - um = dm[pp.STATE]["u"] - lam = d[pp.STATE]["lam_u"] - ps = ds[pp.STATE]["p"] - pm = dm[pp.STATE]["p"] - - bdcs = ds[pp.DISCRETIZATION_MATRICES][self.kw]["bound_displacement_cell"] - bdcm = dm[pp.DISCRETIZATION_MATRICES][self.kw]["bound_displacement_cell"] - bdfs = ds[pp.DISCRETIZATION_MATRICES][self.kw]["bound_displacement_face"] - bdfm = dm[pp.DISCRETIZATION_MATRICES][self.kw]["bound_displacement_face"] - bdps = ds[pp.DISCRETIZATION_MATRICES][self.kw][ - "bound_displacement_pressure" - ] - bdpm = dm[pp.DISCRETIZATION_MATRICES][self.kw][ - "bound_displacement_pressure" - ] - - bc_val_s = ds[pp.PARAMETERS][self.kw]["bc_values"] - bc_val_m = dm[pp.PARAMETERS][self.kw]["bc_values"] - - RW = d[pp.PARAMETERS][self.kw]["robin_weight"] - MW = d[pp.PARAMETERS][self.kw]["mortar_weight"] - rhs = d[pp.PARAMETERS][self.kw]["robin_rhs"].reshape( - (gs.dim, -1), order="F" - ) - slv2mrt_nd = sps.kron(mg.slave_to_mortar_int(), sps.eye(gs.dim)).tocsr() - mstr2mrt_nd = sps.kron(mg.master_to_mortar_int(), sps.eye(gs.dim)).tocsr() - - hf2fs = pp.fvutils.map_hf_2_f(g=gs) / 2 - hf2fm = pp.fvutils.map_hf_2_f(g=gm) / 2 - jump_u = ( - slv2mrt_nd - * hf2fs - * (bdcs * us + bdfs * (bc_val_s + slv2mrt_nd.T * lam) + bdps * ps) - - mstr2mrt_nd - * hf2fm - * (bdcm * um - bdfm * (bc_val_m + mstr2mrt_nd.T * lam) + bdpm * pm) - ).reshape((gs.dim, -1), order="F") - lam_nd = lam.reshape((gs.dim, -1), order="F") - for i in range(len(RW)): - rhs_robin = MW[i].dot(lam_nd[:, i]) + RW[i].dot(jump_u[:, i]) - self.assertTrue(np.allclose(rhs_robin, rhs[:, i])) - - def assign_discretization(self, gb): - gradP_disc = pp.GradP(self.kw) - divU_disc = pp.DivD(self.kw, variable="u", mortar_variable="lam_u") - - for g, d in gb: - d[pp.DISCRETIZATION] = { - "u": {"div_sigma": pp.Mpsa(self.kw)}, - "p": { - "flux": pp.Mpfa(self.kw_f), - "mass": pp.MassMatrix(self.kw_f), - "stab": pp.BiotStabilization(self.kw_f, variable="p"), - }, - "u_p": {"grad_p": gradP_disc}, - "p_u": {"div_u": divU_disc}, - } - - d[pp.PRIMARY_VARIABLES] = {"u": {"cells": g.dim}, "p": {"cells": 1}} - - gradP_disp = pp.numerics.interface_laws.elliptic_interface_laws.RobinContactBiotPressure( - self.kw, gradP_disc - ) - div_u_lam = pp.numerics.interface_laws.elliptic_interface_laws.DivU_StressMortar( - self.kw, divU_disc - ) - for e, d in gb.edges(): - g1, g2 = gb.nodes_of_edge(e) - d[pp.PRIMARY_VARIABLES] = {"lam_u": {"cells": 2}, "lam_p": {"cells": 1}} - d[pp.COUPLING_DISCRETIZATION] = { - "u_contribution": { - g1: ("u", "div_sigma"), - g2: ("u", "div_sigma"), - (g1, g2): ("lam_u", pp.RobinContact(self.kw, pp.Mpsa(self.kw))), - }, - "p_contribution_to_displacement": { - g1: ("p", "flux"), - g2: ("p", "flux"), - (g1, g2): ("lam_u", gradP_disp), - }, - "lam_u_contr_2_div_u": { - g1: ("p", "flux"), - g2: ("p", "flux"), - (g1, g2): ("lam_u", div_u_lam), - }, - "lam_p": { - g1: ("p", "flux"), - g2: ("p", "flux"), - (g1, g2): ( - "lam_p", - pp.FluxPressureContinuity(self.kw_f, pp.Mpfa(self.kw_f)), - ), - }, - } - - def assign_parameters(self, gb, mortar_weight, robin_weight, rhs): - for g, d in gb: - if g.grid_num == 1: - dir_faces = g.face_centers[0] < 1e-10 - elif g.grid_num == 2: - dir_faces = g.face_centers[0] > 2 - 1e-10 - elif g.grid_num == 3: - dir_faces = (g.face_centers[0] < 1e-10) + ( - g.face_centers[0] > 2 - 1e-10 - ) - bc_val = np.zeros((g.dim, g.num_faces)) - bc_val[0, dir_faces] = 0.1 * g.face_centers[0, dir_faces] - bc = pp.BoundaryConditionVectorial(g, dir_faces, "dir") - C = pp.FourthOrderTensor( - gb.dim_max(), np.ones(g.num_cells), np.ones(g.num_cells) - ) - alpha = 1 / np.pi - data = { - "bc": bc, - "bc_values": bc_val.ravel("F"), - "fourth_order_tensor": C, - "source": np.zeros(g.num_cells * g.dim), - "inverter": "python", - "biot_alpha": alpha, - } - data_f = { - "bc": pp.BoundaryCondition(g, g.get_boundary_faces(), "dir"), - "bc_values": np.zeros(g.num_faces), - "second_order_tensor": pp.SecondOrderTensor( - g.dim, np.ones(g.num_cells) - ), - "inverter": "python", - "aperture": np.ones(g.num_cells), - "biot_alpha": alpha, - "mass_weight": 1e-1, - } - pp.initialize_data(g, d, self.kw, data) - pp.initialize_data(g, d, self.kw_f, data_f) - state = { - "u": np.zeros(g.dim * g.num_cells), - "p": np.zeros(g.num_cells), - self.kw: {"bc_values": bc_val.ravel("F")}, - } - pp.set_state(d, state) - - for _, d in gb.edges(): - mg = d["mortar_grid"] - MW = sps.diags(mortar_weight) - RW = sps.diags(robin_weight) - data = { - "mortar_weight": [MW] * mg.num_cells, - "robin_weight": [RW] * mg.num_cells, - "robin_rhs": np.tile(rhs, (mg.num_cells)), - } - pp.initialize_data(mg, d, self.kw, data) - pp.set_state(d, {"lam_u": np.zeros((mg.dim + 1) * mg.num_cells)}) - - -def define_gb(): - """ - Construct grids - """ - g_s = pp.CartGrid([1, 2], [1, 2]) - g_m = pp.CartGrid([1, 2], [1, 2]) - g_full = pp.CartGrid([2, 2], [2, 2]) - g_m.nodes[0] += 1 - g_s.compute_geometry() - g_m.compute_geometry() - g_full.compute_geometry() - - g_s.grid_num = 1 - g_m.grid_num = 2 - g_full.grid_num = 3 - - gb = pp.GridBucket() - gb_full = pp.GridBucket() - gb.add_nodes([g_s, g_m]) - gb_full.add_nodes([g_full]) - - contact_s = np.where(g_s.face_centers[0] > 1 - 1e-10)[0] - contact_m = np.where(g_m.face_centers[0] < 1 + 1e-10)[0] - data = np.ones(contact_s.size, dtype=np.bool) - - shape = (g_s.num_faces, g_m.num_faces) - slave_master = sps.csc_matrix((data, (contact_m, contact_s)), shape=shape) - - mortar_grid, _, _ = pp.grids.partition.extract_subgrid(g_s, contact_s, faces=True) - - gb.add_edge([g_s, g_m], slave_master) - - gb.assign_node_ordering() - gb_full.assign_node_ordering() - - # Slave and master is defined by the node number. - # In python 3.5 the node-nombering does not follow the one given in gb.add_edge - # I guess also the face_face mapping given on the edge also should change, - # but this is not used - g_1, _ = gb.nodes_of_edge([g_s, g_m]) - if g_1.grid_num == 2: - g_m = g_s - g_s = g_1 - slave_master = slave_master.T - - mg = pp.BoundaryMortar(mortar_grid.dim, mortar_grid, slave_master.T) - gb.set_edge_prop([g_s, g_m], "mortar_grid", mg) - return gb, gb_full - - -if __name__ == "__main__": - unittest.main() diff --git a/test/integration/test_fractured_mpsa.py b/test/integration/test_fractured_mpsa.py deleted file mode 100644 index 13ee8883cb..0000000000 --- a/test/integration/test_fractured_mpsa.py +++ /dev/null @@ -1,359 +0,0 @@ -""" -Tests of class FracturedMpsa in module porepy.numerics.fv.mpsa. -""" -import numpy as np -import unittest - -import porepy as pp - - -class BasicsTest(unittest.TestCase): - def __init__(self, *args, **kwargs): - f = pp.Fracture(np.array([[1, 1, 4], [3, 4, 1], [2, 2, 4]])) - box = {"xmin": 0, "ymin": 0, "zmin": 0, "xmax": 5, "ymax": 5, "zmax": 5} - network = pp.FractureNetwork3d([f], domain=box) - mesh_args = {"mesh_size_min": 5, "mesh_size_frac": 5} - self.gb3d = network.mesh(mesh_args) - unittest.TestCase.__init__(self, *args, **kwargs) - - # ------------------------------------------------------------------------------# - - def test_zero_force(self): - """ - test that nothing moves if nothing is touched - """ - g = self.gb3d.grids_of_dimension(3)[0] - - bound = pp.BoundaryConditionVectorial(g, g.get_all_boundary_faces(), "dir") - - specified_parameters = {"bc": bound, "inverter": "python"} - data = pp.initialize_default_data(g, {}, "mechanics", specified_parameters) - solver = pp.FracturedMpsa("mechanics") - - A, b = solver.assemble_matrix_rhs(g, data) - - u = np.linalg.solve(A.A, b) - T = solver.traction(g, data, u) - - self.assertTrue(np.all(np.abs(u) < 1e-10)) - self.assertTrue(np.all(np.abs(T) < 1e-10)) - - def test_unit_slip(self): - """ - test unit slip of fractures - """ - frac = np.array([[1, 1, 1], [1, 2, 1], [2, 2, 1], [2, 1, 1]]).T - physdims = np.array([3, 3, 2]) - g = pp.meshing.cart_grid( - [frac], [3, 3, 2], physdims=physdims - ).grids_of_dimension(3)[0] - - bound = pp.BoundaryConditionVectorial(g, g.get_all_boundary_faces(), "dir") - - frac_slip = np.zeros((g.dim, g.num_faces)) - frac_bnd = g.tags["fracture_faces"] - frac_slip[:, frac_bnd] = np.ones((g.dim, np.sum(frac_bnd))) - - specified_parameters = { - "bc": bound, - "slip_distance": frac_slip.ravel("F"), - "inverter": "python", - } - data = pp.initialize_default_data(g, {}, "mechanics", specified_parameters) - solver = pp.FracturedMpsa("mechanics") - - A, b = solver.assemble_matrix_rhs(g, data) - - u = np.linalg.solve(A.A, b) - - u_f = solver.extract_frac_u(g, u) - u_c = solver.extract_u(g, u) - u_c = u_c.reshape((3, -1), order="F") - - # obtain fracture faces and cells - frac_faces = g.frac_pairs - frac_left = frac_faces[0] - frac_right = frac_faces[1] - - cell_left = np.ravel(np.argwhere(g.cell_faces[frac_left, :])[:, 1]) - cell_right = np.ravel(np.argwhere(g.cell_faces[frac_right, :])[:, 1]) - - # Test traction - T = solver.traction(g, data, u) - T = T.reshape((3, -1), order="F") - T_left = T[:, frac_left] - T_right = T[:, frac_right] - - self.assertTrue(np.allclose(T_left, T_right)) - - # we have u_lhs - u_rhs = 1 so u_lhs should be positive - self.assertTrue(np.all(u_c[:, cell_left] > 0)) - self.assertTrue(np.all(u_c[:, cell_right] < 0)) - mid_ind = int(round(u_f.size / 2)) - u_left = u_f[:mid_ind] - u_right = u_f[mid_ind:] - self.assertTrue(np.all(np.abs(u_left - u_right - 1) < 1e-10)) - - # fracture displacement should be symetric since everything else is - # symetric - self.assertTrue(np.allclose(u_left, 0.5)) - self.assertTrue(np.allclose(u_right, -0.5)) - - def test_non_zero_bc_val(self): - """ - We mixed bc_val on domain boundary and fracture displacement in - x-direction. - """ - frac = np.array([[1, 1, 1], [1, 2, 1], [2, 2, 1], [2, 1, 1]]).T - physdims = np.array([3, 3, 2]) - - g = pp.meshing.cart_grid( - [frac], [3, 3, 2], physdims=physdims - ).grids_of_dimension(3)[0] - - # Define boundary conditions - bc_val = np.zeros((g.dim, g.num_faces)) - frac_slip = np.zeros((g.dim, g.num_faces)) - - frac_bnd = g.tags["fracture_faces"] - dom_bnd = g.tags["domain_boundary_faces"] - - frac_slip[0, frac_bnd] = np.ones(np.sum(frac_bnd)) - bc_val[:, dom_bnd] = g.face_centers[:, dom_bnd] - - bound = pp.BoundaryConditionVectorial(g, g.get_all_boundary_faces(), "dir") - - specified_parameters = { - "bc": bound, - "bc_values": bc_val.ravel("F"), - "slip_distance": frac_slip.ravel("F"), - "inverter": "python", - } - data = pp.initialize_default_data(g, {}, "mechanics", specified_parameters) - solver = pp.FracturedMpsa("mechanics") - - A, b = solver.assemble_matrix_rhs(g, data) - - u = np.linalg.solve(A.A, b) - - u_f = solver.extract_frac_u(g, u) - u_c = solver.extract_u(g, u) - u_c = u_c.reshape((3, -1), order="F") - - # Test traction - frac_faces = g.frac_pairs - frac_left = frac_faces[0] - frac_right = frac_faces[1] - - T = solver.traction(g, data, u) - T = T.reshape((3, -1), order="F") - T_left = T[:, frac_left] - T_right = T[:, frac_right] - - self.assertTrue(np.allclose(T_left, T_right)) - - # we have u_lhs - u_rhs = 1 so u_lhs should be positive - mid_ind = int(round(u_f.size / 2)) - u_left = u_f[:mid_ind] - u_right = u_f[mid_ind:] - - true_diff = np.atleast_2d(np.array([1, 0, 0])).T - u_left = u_left.reshape((3, -1), order="F") - u_right = u_right.reshape((3, -1), order="F") - self.assertTrue(np.all(np.abs(u_left - u_right - true_diff) < 1e-10)) - - # should have a positive displacement for all cells - self.assertTrue(np.all(u_c > 0)) - - def test_given_traction_on_fracture(self): - """ - We specify the traction on the fracture faces. - """ - frac = np.array([[1, 1, 1], [1, 2, 1], [2, 2, 1], [2, 1, 1]]).T - normal_ind = 2 - physdims = np.array([3, 3, 2]) - - g = pp.meshing.cart_grid( - [frac], [3, 3, 2], physdims=physdims - ).grids_of_dimension(3)[0] - - # Define boundary conditions - bc_val = np.zeros((g.dim, g.num_faces)) - frac_traction = np.zeros((g.dim, g.num_faces)) - - frac_bnd = g.tags["fracture_faces"] - # Positive values in the normal direction correspond to a normal force - # pointing from the fracture to the matrix. - frac_traction[normal_ind, frac_bnd] = np.ones(np.sum(frac_bnd)) - - bound = pp.BoundaryConditionVectorial(g, g.get_all_boundary_faces(), "dir") - - # Even though we now prescribe the traction, the discretisation uses - # the same parameter function "get_slip_distance" - specified_parameters = { - "bc": bound, - "bc_values": bc_val.ravel("F"), - "slip_distance": frac_traction.ravel("F"), - "inverter": "python", - } - data = pp.initialize_default_data(g, {}, "mechanics", specified_parameters) - solver = pp.FracturedMpsa("mechanics", given_traction=True) - - A, b = solver.assemble_matrix_rhs(g, data) - - u = np.linalg.solve(A.A, b) - - u_f = solver.extract_frac_u(g, u) - u_c = solver.extract_u(g, u) - u_c = u_c.reshape((3, -1), order="F") - - # Test traction - frac_faces = g.frac_pairs - frac_left = frac_faces[0] - frac_right = frac_faces[1] - - T = solver.traction(g, data, u) - T = T.reshape((3, -1), order="F") - T_left = T[:, frac_left] - T_right = T[:, frac_right] - self.assertTrue(np.all(np.isclose(T_left - T_right, 0))) - - # we have u_lhs - u_rhs = 1 so u_lhs should be positive - mid_ind = int(round(u_f.size / 2)) - u_left = u_f[:mid_ind] - u_right = u_f[mid_ind:] - - u_left = u_left.reshape((3, -1), order="F") - u_right = u_right.reshape((3, -1), order="F") - # The normal displacements should be equal and of opposite direction. - self.assertTrue(np.all(np.isclose(u_left + u_right, 0))) - # They should also correspond to an opening of the fracture - self.assertTrue(np.all((u_left - u_right)[normal_ind] > 0.2)) - - # The maximum displacement magnitude should be observed at the fracture - self.assertTrue(np.all(np.abs(u_c) < np.max(u_left[normal_ind]))) - - def test_domain_cut_in_two(self): - """ - test domain cut in two. We place 1 dirichlet on top. zero dirichlet on - bottom and 0 neumann on sides. Further we place 1 displacement on - fracture. this should give us displacement 1 on top cells and 0 on - bottom cells and zero traction on all faces - """ - - frac = np.array([[0, 0, 1], [0, 3, 1], [3, 3, 1], [3, 0, 1]]).T - g = pp.meshing.cart_grid([frac], [3, 3, 2]).grids_of_dimension(3)[0] - - tol = 1e-6 - frac_bnd = g.tags["fracture_faces"] - top = g.face_centers[2] > 2 - tol - bot = g.face_centers[2] < tol - - dir_bound = top | bot | frac_bnd - - bound = pp.BoundaryConditionVectorial(g, dir_bound, "dir") - - bc_val = np.zeros((g.dim, g.num_faces)) - bc_val[:, top] = np.ones((g.dim, np.sum(top))) - frac_slip = np.zeros((g.dim, g.num_faces)) - frac_slip[:, frac_bnd] = np.ones(np.sum(frac_bnd)) - - specified_parameters = { - "bc": bound, - "bc_values": bc_val.ravel("F"), - "slip_distance": frac_slip.ravel("F"), - "inverter": "python", - } - data = pp.initialize_default_data(g, {}, "mechanics", specified_parameters) - solver = pp.FracturedMpsa("mechanics") - - A, b = solver.assemble_matrix_rhs(g, data) - - u = np.linalg.solve(A.A, b) - - u_f = solver.extract_frac_u(g, u) - u_c = solver.extract_u(g, u) - u_c = u_c.reshape((3, -1), order="F") - T = solver.traction(g, data, u) - - top_cells = g.cell_centers[2] > 1 - - mid_ind = int(round(u_f.size / 2)) - u_left = u_f[:mid_ind] - u_right = u_f[mid_ind:] - - self.assertTrue(np.allclose(u_left, 1)) - self.assertTrue(np.allclose(u_right, 0)) - self.assertTrue(np.allclose(u_c[:, top_cells], 1)) - self.assertTrue(np.allclose(u_c[:, ~top_cells], 0)) - self.assertTrue(np.allclose(T, 0)) - - def test_vectorial_bc(self): - """ - We mixed bc_val on domain boundary and fracture displacement in - x-direction. - """ - frac = np.array([[1, 1, 1], [1, 2, 1], [2, 2, 1], [2, 1, 1]]).T - physdims = np.array([3, 3, 2]) - - g = pp.meshing.cart_grid( - [frac], [3, 3, 2], physdims=physdims - ).grids_of_dimension(3)[0] - data = {"param": pp.Parameters(g)} - - # Define boundary conditions - bc_val = np.zeros((g.dim, g.num_faces)) - frac_slip = np.zeros((g.dim, g.num_faces)) - - frac_bnd = g.tags["fracture_faces"] - dom_bnd = g.tags["domain_boundary_faces"] - - frac_slip[0, frac_bnd] = np.ones(np.sum(frac_bnd)) - bc_val[:, dom_bnd] = g.face_centers[:, dom_bnd] - - bound = pp.BoundaryConditionVectorial(g, g.get_all_boundary_faces(), "dir") - specified_parameters = { - "bc": bound, - "bc_values": bc_val.ravel("F"), - "slip_distance": frac_slip.ravel("F"), - "inverter": "python", - } - data = pp.initialize_default_data(g, {}, "mechanics", specified_parameters) - solver = pp.FracturedMpsa("mechanics") - - A, b = solver.assemble_matrix_rhs(g, data) - - u = np.linalg.solve(A.A, b) - - u_f = solver.extract_frac_u(g, u) - u_c = solver.extract_u(g, u) - u_c = u_c.reshape((3, -1), order="F") - - # Test traction - frac_faces = g.frac_pairs - frac_left = frac_faces[0] - frac_right = frac_faces[1] - - T = solver.traction(g, data, u) - T = T.reshape((3, -1), order="F") - T_left = T[:, frac_left] - T_right = T[:, frac_right] - - self.assertTrue(np.allclose(T_left, T_right)) - - # we have u_lhs - u_rhs = 1 so u_lhs should be positive - mid_ind = int(round(u_f.size / 2)) - u_left = u_f[:mid_ind] - u_right = u_f[mid_ind:] - - true_diff = np.atleast_2d(np.array([1, 0, 0])).T - u_left = u_left.reshape((3, -1), order="F") - u_right = u_right.reshape((3, -1), order="F") - self.assertTrue(np.all(np.abs(u_left - u_right - true_diff) < 1e-10)) - - # should have a positive displacement for all cells - self.assertTrue(np.all(u_c > 0)) - - -if __name__ == "__main__": - unittest.main() diff --git a/test/integration/test_mpfa_mpsa_partial_update.py b/test/integration/test_mpfa_mpsa_partial_update.py index 0ed4d002ed..107122a9f9 100644 --- a/test/integration/test_mpfa_mpsa_partial_update.py +++ b/test/integration/test_mpfa_mpsa_partial_update.py @@ -146,7 +146,7 @@ def test_inner_cell_node_keyword(self): nodes_of_cell = np.array([14, 15, 20, 21]) faces_of_cell = np.array([14, 15, 42, 47]) - partial_stress, partial_bound, _, _, active_faces = mpsa.mpsa_partial( + partial_stress, partial_bound, active_faces = mpsa.mpsa_partial( g, stiffness, bnd, nodes=nodes_of_cell, inverter="python" ) @@ -173,7 +173,7 @@ def test_bound_cell_node_keyword(self): inner_cell = 10 nodes_of_cell = np.array([12, 13, 18, 19]) faces_of_cell = np.array([12, 13, 40, 45]) - partial_stress, partial_bound, _, _, active_faces = mpsa.mpsa_partial( + partial_stress, partial_bound, active_faces = mpsa.mpsa_partial( g, perm, bnd, nodes=nodes_of_cell, inverter="python" ) @@ -221,7 +221,7 @@ def test_one_cell_a_time_node_keyword(self): ind = np.zeros(g.num_cells) ind[ci] = 1 nodes = np.squeeze(np.where(cn * ind > 0)) - partial_stress, partial_bound, _, _, active_faces = mpsa.mpsa_partial( + partial_stress, partial_bound, active_faces = mpsa.mpsa_partial( g, stiffness, bnd, nodes=nodes, inverter="python" ) diff --git a/test/integration/test_mpsa_robin.py b/test/integration/test_mpsa_robin.py index fcdafd5272..26f7294b3d 100644 --- a/test/integration/test_mpsa_robin.py +++ b/test/integration/test_mpsa_robin.py @@ -39,8 +39,8 @@ def T_ex(faces): u_bound = np.zeros((2, g.num_faces)) - sgn_n = pp.numerics.fracture_deformation.sign_of_faces(g, neu_ind) - sgn_r = pp.numerics.fracture_deformation.sign_of_faces(g, rob_ind) + sgn_n = g.sign_of_faces(neu_ind) + sgn_r = g.sign_of_faces(rob_ind) u_bound[:, dir_ind] = u_ex(g.face_centers[:, dir_ind]) u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n @@ -86,8 +86,8 @@ def T_ex(faces): u_bound = np.zeros((2, g.num_faces)) - sgn_n = pp.numerics.fracture_deformation.sign_of_faces(g, neu_ind) - sgn_r = pp.numerics.fracture_deformation.sign_of_faces(g, rob_ind) + sgn_n = g.sign_of_faces(neu_ind) + sgn_r = g.sign_of_faces(rob_ind) u_bound[:, dir_ind] = u_ex(g.face_centers[:, dir_ind]) u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n @@ -133,8 +133,8 @@ def T_ex(faces): u_bound = np.zeros((2, g.num_faces)) - sgn_n = pp.numerics.fracture_deformation.sign_of_faces(g, neu_ind) - sgn_r = pp.numerics.fracture_deformation.sign_of_faces(g, rob_ind) + sgn_n = g.sign_of_faces(neu_ind) + sgn_r = g.sign_of_faces(rob_ind) u_bound[:, dir_ind] = u_ex(g.face_centers[:, dir_ind]) u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n @@ -181,8 +181,8 @@ def T_ex(faces): u_bound = np.zeros((2, g.num_faces)) - sgn_n = pp.numerics.fracture_deformation.sign_of_faces(g, neu_ind) - sgn_r = pp.numerics.fracture_deformation.sign_of_faces(g, rob_ind) + sgn_n = g.sign_of_faces(neu_ind) + sgn_r = g.sign_of_faces(rob_ind) u_bound[:, dir_ind] = u_ex(g.face_centers[:, dir_ind]) u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n @@ -231,8 +231,8 @@ def T_ex(faces): u_bound = np.zeros((3, g.num_faces)) - sgn_n = pp.numerics.fracture_deformation.sign_of_faces(g, neu_ind) - sgn_r = pp.numerics.fracture_deformation.sign_of_faces(g, rob_ind) + sgn_n = g.sign_of_faces(neu_ind) + sgn_r = g.sign_of_faces(rob_ind) u_bound[:, dir_ind] = u_ex(g.face_centers[:, dir_ind]) u_bound[:, neu_ind] = T_ex(neu_ind) * sgn_n diff --git a/test/integration/test_upwind_coupling.py b/test/integration/test_upwind_coupling.py index b222fdcf55..a62d92a155 100644 --- a/test/integration/test_upwind_coupling.py +++ b/test/integration/test_upwind_coupling.py @@ -3,7 +3,6 @@ import unittest import porepy as pp -from porepy.utils.grid_util import sign_of_boundary_faces from test.integration import _helper_test_upwind_coupling from test.test_utils import permute_matrix_vector @@ -55,7 +54,7 @@ def test_upwind_coupling_2d_1d_bottom_top(self): add_constant_darcy_flux(gb, upwind, [0, 1, 0], a) assembler = pp.Assembler(gb) - + assembler.discretize() U_tmp, rhs = assembler.assemble_matrix_rhs() grids = np.empty(gb.num_graph_nodes() + gb.num_graph_edges(), dtype=np.object) @@ -142,6 +141,7 @@ def test_upwind_coupling_2d_1d_left_right(self): add_constant_darcy_flux(gb, upwind, [1, 0, 0], a) assembler = pp.Assembler(gb) + assembler.discretize() U_tmp, rhs = assembler.assemble_matrix_rhs() @@ -278,6 +278,7 @@ def test_upwind_coupling_2d_1d_left_right_cross(self): add_constant_darcy_flux(gb, upwind, [1, 0, 0], a) assembler = pp.Assembler(gb) + assembler.discretize() U_tmp, rhs = assembler.assemble_matrix_rhs() @@ -867,6 +868,7 @@ def test_upwind_coupling_3d_2d_bottom_top(self): add_constant_darcy_flux(gb, upwind, [0, 0, 1], a) assembler = pp.Assembler(gb) + assembler.discretize() U_tmp, rhs = assembler.assemble_matrix_rhs() @@ -950,6 +952,7 @@ def test_upwind_coupling_3d_2d_left_right(self): add_constant_darcy_flux(gb, upwind, [1, 0, 0], a) assembler = pp.Assembler(gb) + assembler.discretize() U_tmp, rhs = assembler.assemble_matrix_rhs() @@ -1126,6 +1129,7 @@ def test_upwind_coupling_3d_2d_1d_0d(self): add_constant_darcy_flux(gb, upwind, [1, 0, 0], a) assembler = pp.Assembler(gb) + assembler.discretize() U_tmp, rhs = assembler.assemble_matrix_rhs() grids = np.empty(gb.num_graph_nodes() + gb.num_graph_edges(), dtype=np.object) @@ -1267,6 +1271,7 @@ def test_upwind_2d_beta_positive(self): add_constant_darcy_flux(gb, upwind, [2, 0, 0], a) assembler = pp.Assembler(gb) + assembler.discretize() U_tmp, rhs = assembler.assemble_matrix_rhs() @@ -1426,6 +1431,7 @@ def test_upwind_2d_full_beta_bc_dir(self): add_constant_darcy_flux(gb, upwind, [1, 1, 0], a) assembler = pp.Assembler(gb) + assembler.discretize() U_tmp, rhs = assembler.assemble_matrix_rhs() @@ -1607,7 +1613,9 @@ def add_constant_darcy_flux(gb, upwind, flux, a): p_h = gb.node_props(g_h, pp.PARAMETERS) darcy_flux = p_h["transport"]["darcy_flux"] sign = np.zeros(g_h.num_faces) - sign[g_h.get_all_boundary_faces()] = sign_of_boundary_faces(g_h) + sign[g_h.get_all_boundary_faces()] = g_h.sign_of_faces( + g_h.get_all_boundary_faces() + ) mg = d["mortar_grid"] sign = mg.master_to_mortar_avg() * sign # d["param"] = pp.Parameters(g_h) diff --git a/test/unit/test_elliptic_interface_laws.py b/test/unit/test_elliptic_interface_laws.py deleted file mode 100644 index e2d3aac33d..0000000000 --- a/test/unit/test_elliptic_interface_laws.py +++ /dev/null @@ -1,346 +0,0 @@ -""" -Module for testing the vector elliptic couplings in the interface_laws. -""" - - -import numpy as np -import scipy.sparse as sps -import unittest - -import porepy as pp - - -class TestTwoGridCoupling(unittest.TestCase): - """ - In this test we set up a coupling between two grids of dimension 2: - g_slave g_master - |-----| | |------| - | | x | | - |-----| | |------| - g_mortar - There is one cell per grid and they are coupled together by a single mortar - variable. - """ - - def test_robin_assembly_master(self): - """ - We test the RobinContact interface law. This gives a Robin condition - on the mortar grid. - - Test the assembly of the master terms. - """ - self.kw = "mech" - gb = define_gb() - mortar_weight = np.zeros(gb.dim_max()) - robin_weight = np.ones(gb.dim_max()) - rhs = np.ones(gb.dim_max()) - self.assign_parameters(gb, mortar_weight, robin_weight, rhs) - varargs = get_variables(gb) - - robin_contact = pp.RobinContact(self.kw, MockId(), MockZero()) - matrix, rhs = robin_contact.assemble_matrix_rhs(*varargs) - - # known coupling - A = np.array( - [ - [2, 0, 2, 0, -2, 0], - [0, 2, 0, 2, 0, -2], - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [-2, 0, 0, 0, 2, 0], - [0, -2, 0, 0, 0, 2], - ] - ) - - matrix = sps.bmat(matrix) - self.assertTrue(np.allclose(A, matrix.A)) - self.assertTrue( - np.allclose(np.hstack(rhs.ravel()), np.array([0, 0, 0, 0, 1, 1])) - ) - - def test_robin_assembly_slave(self): - """ - We test the RobinContact interface law. This gives a Robin condition - on the mortar grid. - - Test the assembly of the slave terms. - """ - self.kw = "mech" - gb = define_gb() - mortar_weight = np.zeros(gb.dim_max()) - robin_weight = np.ones(gb.dim_max()) - rhs = np.ones(gb.dim_max()) - self.assign_parameters(gb, mortar_weight, robin_weight, rhs) - varargs = get_variables(gb) - - robin_contact = pp.RobinContact(self.kw, MockZero(), MockId()) - matrix, rhs = robin_contact.assemble_matrix_rhs(*varargs) - - # known coupling - A = np.array( - [ - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [2, 0, 2, 0, 2, 0], - [0, 2, 0, 2, 0, 2], - [0, 0, 2, 0, 2, 0], - [0, 0, 0, 2, 0, 2], - ] - ) - - matrix = sps.bmat(matrix) - self.assertTrue(np.allclose(A, matrix.A)) - self.assertTrue( - np.allclose(np.hstack(rhs.ravel()), np.array([0, 0, 0, 0, 1, 1])) - ) - - def test_robin_assembly_mortar(self): - """ - We test the RobinContact interface law. This gives a Robin condition - on the mortar grid. - - Test the assembly of the mortar terms. - """ - self.kw = "mech" - gb = define_gb() - mortar_weight = np.ones(gb.dim_max()) - robin_weight = np.zeros(gb.dim_max()) - rhs = np.ones(gb.dim_max()) - self.assign_parameters(gb, mortar_weight, robin_weight, rhs) - varargs = get_variables(gb) - - robin_contact = pp.RobinContact(self.kw, MockZero(), MockZero()) - matrix, rhs = robin_contact.assemble_matrix_rhs(*varargs) - - # known coupling - A = np.array( - [ - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 1, 0], - [0, 0, 0, 0, 0, 1], - ] - ) - - matrix = sps.bmat(matrix) - self.assertTrue(np.allclose(A, matrix.A)) - self.assertTrue( - np.allclose(np.hstack(rhs.ravel()), np.array([0, 0, 0, 0, 1, 1])) - ) - - def test_continuity_assembly_master(self): - """ - We test the RobinContact interface law. This gives a Robin condition - on the mortar grid. - - Test the assembly of the master terms. - """ - self.kw = "mech" - gb = define_gb() - varargs = get_variables(gb) - robin_contact = pp.StressDisplacementContinuity(self.kw, MockId(), MockZero()) - matrix, rhs = robin_contact.assemble_matrix_rhs(*varargs) - - # known coupling - A = np.array( - [ - [2, 0, 2, 0, -2, 0], - [0, 2, 0, 2, 0, -2], - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [-2, 0, 0, 0, 2, 0], - [0, -2, 0, 0, 0, 2], - ] - ) - - matrix = sps.bmat(matrix) - - self.assertTrue(np.allclose(A, matrix.A)) - self.assertTrue( - np.allclose(np.hstack(rhs.ravel()), np.array([0, 0, 0, 0, 0, 0])) - ) - - def test_continuity_assembly_slave(self): - """ - We test the StressDisplacementContinuity interface law. This gives a continuity - of stress and displacement on the mortar grid. - - Test the assembly of the slave terms. - """ - self.kw = "mech" - gb = define_gb() - varargs = get_variables(gb) - - robin_contact = pp.StressDisplacementContinuity(self.kw, MockZero(), MockId()) - matrix, rhs = robin_contact.assemble_matrix_rhs(*varargs) - - # known coupling - A = np.array( - [ - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [2, 0, 2, 0, 2, 0], - [0, 2, 0, 2, 0, 2], - [0, 0, 2, 0, 2, 0], - [0, 0, 0, 2, 0, 2], - ] - ) - - matrix = sps.bmat(matrix) - self.assertTrue(np.allclose(A, matrix.A)) - self.assertTrue( - np.allclose(np.hstack(rhs.ravel()), np.array([0, 0, 0, 0, 0, 0])) - ) - - def test_continuity_assembly_mortar(self): - """ - We test the StressDisplacementContinuity interface law. This gives a continuity - of stress and displacement on the mortar grid. - - Test the assembly of the mortar terms. - """ - self.kw = "mech" - gb = define_gb() - varargs = get_variables(gb) - - robin_contact = pp.StressDisplacementContinuity(self.kw, MockZero(), MockZero()) - matrix, rhs = robin_contact.assemble_matrix_rhs(*varargs) - - # known coupling - A = np.array( - [ - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0], - ] - ) - matrix = sps.bmat(matrix) - self.assertTrue(np.allclose(A, matrix.A)) - self.assertTrue( - np.allclose(np.hstack(rhs.ravel()), np.array([0, 0, 0, 0, 0, 0])) - ) - - def assign_parameters(self, gb, mortar_weight, robin_weight, rhs): - for e, d in gb.edges(): - mg = d["mortar_grid"] - MW = sps.diags(mortar_weight) - RW = sps.diags(robin_weight) - data = { - "mortar_weight": [MW] * mg.num_cells, - "robin_weight": [RW] * mg.num_cells, - "robin_rhs": np.tile(rhs, (mg.num_cells)), - } - d[pp.PARAMETERS] = pp.Parameters(e, self.kw, data) - d[pp.DISCRETIZATION_MATRICES] = {self.kw: {}} - - -class MockId(object): - """ - returns an identity mapping - """ - - def ndof(self, g): - return g.dim * g.num_cells - - def assemble_int_bound_displacement_trace(self, *vargs): - identity_mapping(self, *vargs) - - def assemble_int_bound_stress(self, *vargs): - identity_mapping(self, *vargs) - - def enforce_neumann_int_bound(self, *varargs): - pass - - -class MockZero(object): - """ - Do nothing - """ - - def ndof(self, g): - return g.dim * g.num_cells - - def assemble_int_bound_displacement_trace(self, *vargs): - pass - - def assemble_int_bound_stress(self, *vargs): - pass - - def enforce_neumann_int_bound(self, *varargs): - pass - - -def identity_mapping(RC, g, data, data_edge, swap, cc, matrix, rhs, ind): - dof_master = g.dim * g.num_cells - dof_slave = g.dim * g.num_cells - dof_mortar = g.dim * data_edge["mortar_grid"].num_cells - - cc[ind, 0] += sps.diags(np.ones(dof_slave), shape=(dof_slave, dof_master)) - cc[ind, 1] += sps.diags(np.ones(dof_slave), shape=(dof_slave, dof_slave)) - cc[ind, 2] += sps.diags(np.ones(dof_slave), shape=(dof_slave, dof_mortar)) - - cc[2, ind] += sps.diags(np.ones(dof_mortar), shape=(dof_mortar, dof_master)) - cc[2, 2] += sps.diags(np.ones(dof_mortar), shape=(dof_mortar, dof_mortar)) - - -def get_variables(gb): - for e, data_edge in gb.edges(): - g_slave, g_master = gb.nodes_of_edge(e) - data_slave = gb.node_props(g_slave) - data_master = gb.node_props(g_master) - break - dof = np.array( - [ - g_master.dim * g_master.num_cells, - g_slave.dim * g_slave.num_cells, - g_master.dim * data_edge["mortar_grid"].num_cells, - ] - ) - matrix = np.array([sps.coo_matrix((i, j)) for i in dof for j in dof]) - matrix = matrix.reshape((3, 3)) - - rhs = np.empty(3, dtype=np.object) - rhs[0] = np.zeros(g_master.dim * g_master.num_cells) - rhs[1] = np.zeros(g_slave.dim * g_slave.num_cells) - rhs[2] = np.zeros(data_edge["mortar_grid"].num_cells * g_slave.dim) - - return g_master, g_slave, data_master, data_slave, data_edge, matrix - - -def define_gb(): - """ - Construct grids - """ - g1 = pp.CartGrid([1, 1]) - g2 = pp.CartGrid([1, 1]) - g1.compute_geometry() - g2.compute_geometry() - - g1.grid_num = 1 - g2.grid_num = 2 - - gb = pp.GridBucket() - gb.add_nodes([g1, g2]) - gb.add_edge([g1, g2], None) - mortar_grid = pp.Grid( - 1, - np.array([[0, 0, 0], [0, 1, 0]]).T, - sps.csc_matrix(([True, True], [0, 1], [0, 1, 2])), - sps.csc_matrix(([1, -1], [0, 1], [0, 2])), - "mortar_grid", - ) - mortar_grid.compute_geometry() - face_faces = sps.csc_matrix(([True], [0], [0, 0, 1, 1, 1]), shape=(4, 4)) - mg = pp.BoundaryMortar(1, mortar_grid, face_faces) - mg.num_cells = 1 - gb.set_edge_prop([g1, g2], "mortar_grid", mg) - return gb - - -if __name__ == "__main__": - unittest.main() diff --git a/test/unit/test_rock.py b/test/unit/test_rock.py index a2cd97d430..5080762b1a 100644 --- a/test/unit/test_rock.py +++ b/test/unit/test_rock.py @@ -11,8 +11,8 @@ def test_lame_from_young_poisson(self): e = 1 nu = 0.1 lam, mu = pp.params.rock.lame_from_young_poisson(e, nu) - self.assertEqual(lam, 0.11363636363636363) - self.assertEqual(mu, 0.45454545454545453) + self.assertTrue(np.allclose(lam, 0.11363636363636363)) + self.assertTrue(np.allclose(mu, 0.45454545454545453)) def test_poisson_from_lame(self): lam = 1 @@ -31,18 +31,18 @@ def test_sand_stone(self): self.assertEqual(R.POROSITY, 0.2) self.assertEqual(R.YOUNG_MODULUS, 5 * pp.KILOGRAM / pp.CENTI ** 2 * 1e5) self.assertEqual(R.POISSON_RATIO, 0.1) - self.assertEqual(R.LAMBDA, 568181818.1818181) - self.assertEqual(R.MU, 2272727272.7272725) + self.assertTrue(np.allclose(R.LAMBDA, 568181818.1818181)) + self.assertTrue(np.allclose(R.MU, 2272727272.7272725)) def test_granite(self): R = pp.Granite() self.assertEqual(R.PERMEABILITY, 1e-8 * pp.DARCY) self.assertEqual(R.POROSITY, 0.01) - self.assertEqual(R.YOUNG_MODULUS, 4 * pp.GIGA * pp.PASCAL) + self.assertEqual(R.YOUNG_MODULUS, 40 * pp.GIGA * pp.PASCAL) self.assertEqual(R.POISSON_RATIO, 0.2) self.assertEqual(R.DENSITY, 2700 * pp.KILOGRAM / pp.METER ** 3) - self.assertEqual(R.LAMBDA, 1111111111.1111112) - self.assertEqual(R.MU, 1666666666.6666667) + self.assertTrue(np.allclose(R.LAMBDA, 11111111111.1111112)) + self.assertTrue(np.allclose(R.MU, 16666666666.6666667)) def test_shale(self): R = pp.Shale() @@ -50,8 +50,8 @@ def test_shale(self): self.assertEqual(R.POROSITY, 0.01) self.assertEqual(R.YOUNG_MODULUS, 1.5 * pp.KILOGRAM / pp.CENTI ** 2 * 1e5) self.assertEqual(R.POISSON_RATIO, 0.3) - self.assertEqual(R.LAMBDA, 865384615.3846153) - self.assertEqual(R.MU, 576923076.9230769) + self.assertTrue(np.allclose(R.LAMBDA, 865384615.3846153)) + self.assertTrue(np.allclose(R.MU, 576923076.9230769)) if __name__ == "__main__": diff --git a/test/unit/test_tangential_normal_projection.py b/test/unit/test_tangential_normal_projection.py new file mode 100644 index 0000000000..8ca6e08705 --- /dev/null +++ b/test/unit/test_tangential_normal_projection.py @@ -0,0 +1,116 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Unit tests for class pp.TangentialNormalProjection +""" + +import numpy as np +import unittest + +import porepy as pp + + +class TestTangentialNormalProjection(unittest.TestCase): + def setUp(self): + # 2d vectors + self.n2 = np.array([[0, 1, -2], [1, 1, 0]]) + + self.n3 = np.vstack((self.n2, np.array([0, 1, 0]))) + + s2 = np.sqrt(2) + self.n2_normalized = np.array([[0, 1.0 / s2, -1.0], [1, 1.0 / s2, 0]]) + s3 = np.sqrt(3) + self.n3_normalized = np.array( + [[0, 1.0 / s3, -1.0], [1, 1.0 / s3, 0], [0, 1.0 / s3, 0]] + ) + + def test_optional_arguments(self): + # Test automatic check of ambient dimension + proj = pp.TangentialNormalProjection(self.n2) + self.assertTrue(proj.dim == 2) + + proj = pp.TangentialNormalProjection(self.n3) + self.assertTrue(proj.dim == 3) + + def test_normal_vectors(self): + # Test that the normal vectors have + proj = pp.TangentialNormalProjection(self.n2) + for i in range(self.n2.shape[1]): + self.assertTrue( + np.allclose(np.sum(proj.normals[:, i] * self.n2_normalized[:, i]), 1) + ) + + proj = pp.TangentialNormalProjection(self.n3) + + for i in range(self.n3.shape[1]): + self.assertTrue( + np.allclose(np.sum(proj.normals[:, i] * self.n3_normalized[:, i]), 1) + ) + + def _verify_orthonormal(self, proj): + # Check that the basis is an orthonormal set + for i in range(proj.num_vecs): + for j in range(proj.dim): + for k in range(proj.dim): + if j == k: + truth = 1 + else: + truth = 0 + self.assertTrue( + np.allclose( + proj.projection[:, j, i].dot(proj.projection[:, k, i]), + truth, + ) + ) + + def test_computed_basis_2d(self): + # Test that the computed basis functions are orthonormal, and that the + # correct normal vector is constructed + proj = pp.TangentialNormalProjection(self.n2) + self._verify_orthonormal(proj) + + known_projection_of_normal = np.array([0, 1]) + for i in range(self.n2.shape[1]): + + # Check that the projection of the normal vector only has a component in the normal direction + projected_normal = proj.projection[:, :, i].dot(self.n2_normalized[:, i]) + self.assertTrue(np.allclose(projected_normal, known_projection_of_normal)) + + def test_computed_basis_3d(self): + # Test that the computed basis functions are orthonormal, and that the + # correct normal vector is constructed + proj = pp.TangentialNormalProjection(self.n3) + self._verify_orthonormal(proj) + + known_projection_of_normal = np.array([0, 0, 1]) + for i in range(self.n3.shape[1]): + + # Check that the projection of the normal vector only has a component in the normal direction + projected_normal = proj.projection[:, :, i].dot(self.n3_normalized[:, i]) + self.assertTrue(np.allclose(projected_normal, known_projection_of_normal)) + + def test_projections_num_keyword(self): + # Tests of the generated projection operators, using a single tangential/ + # normal space, but generating several (equal) projection matrices. + + dim = 3 + + # Random normal and tangential space + proj = pp.TangentialNormalProjection(np.random.rand(dim, 1)) + + num_reps = 4 + + # Random vector to be generated + vector = np.random.rand(dim, 1) + + projection = proj.project_tangential_normal(num=num_reps) + + proj_vector = projection * np.tile(vector, (num_reps, 1)) + + for i in range(dim): + for j in range(num_reps): + self.assertTrue(proj_vector[i + j * dim], proj_vector[i]) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/unit/test_vtk.py b/test/unit/test_vtk.py index 42f38c44f3..58e6110447 100644 --- a/test/unit/test_vtk.py +++ b/test/unit/test_vtk.py @@ -4,6 +4,7 @@ import porepy as pp + if_vtk = "vtk" in sys.modules if not if_vtk: import warnings @@ -182,10 +183,13 @@ def test_gb_1(self): gb.add_node_props(["scalar_dummy", "dummy_vector"]) for g, d in gb: - d[pp.STATE] = { - "dummy_scalar": np.ones(g.num_cells) * g.dim, - "dummy_vector": np.ones((3, g.num_cells)) * g.dim, - } + pp.set_state( + d, + { + "dummy_scalar": np.ones(g.num_cells) * g.dim, + "dummy_vector": np.ones((3, g.num_cells)) * g.dim, + }, + ) folder = "./test_vtk/" file_name = "grid" @@ -222,10 +226,13 @@ def test_gb_2(self): gb.add_node_props(["dummy_scalar", "dummy_vector"]) for g, d in gb: - d[pp.STATE] = { - "dummy_scalar": np.ones(g.num_cells) * g.dim, - "dummy_vector": np.ones((3, g.num_cells)) * g.dim, - } + pp.set_state( + d, + { + "dummy_scalar": np.ones(g.num_cells) * g.dim, + "dummy_vector": np.ones((3, g.num_cells)) * g.dim, + }, + ) folder = "./test_vtk/" file_name = "grid"