From e15fd4233dacab324036f6033b8e2d089077b9b1 Mon Sep 17 00:00:00 2001 From: IvarStefansson Date: Wed, 3 Jul 2019 10:34:09 +0200 Subject: [PATCH] Move discretization assignment to the models --- .../models/contact_mechanics_biot_model.py | 201 ++++++++++++++-- src/porepy/numerics/fv/mpsa.py | 1 + src/porepy/utils/assign_discretizations.py | 222 ------------------ 3 files changed, 180 insertions(+), 244 deletions(-) delete mode 100644 src/porepy/utils/assign_discretizations.py diff --git a/src/porepy/models/contact_mechanics_biot_model.py b/src/porepy/models/contact_mechanics_biot_model.py index 64b3cefe0e..0e569f81e9 100644 --- a/src/porepy/models/contact_mechanics_biot_model.py +++ b/src/porepy/models/contact_mechanics_biot_model.py @@ -9,8 +9,8 @@ from scipy.spatial.distance import cdist import porepy as pp -from porepy.utils import assign_discretizations import porepy.models.contact_mechanics_model as contact_model +from porepy.utils.derived_discretizations import implicit_euler as IE_discretizations class ContactMechanicsBiot(contact_model.ContactMechanics): @@ -18,21 +18,26 @@ def __init__(self, mesh_args, folder_name): super().__init__(mesh_args, folder_name) # Temperature - self.scalar_variable = "T" + self.scalar_variable = "p" self.mortar_scalar_variable = "mortar_" + self.scalar_variable self.scalar_coupling_term = "robin_" + self.scalar_variable - self.scalar_parameter_key = "temperature" + self.scalar_parameter_key = "flow" # Scaling coefficients self.scalar_scale = 1 self.length_scale = 1 # Time - self.time_step = 1e0 * self.length_scale ** 2 + self.time_step = 1 * self.length_scale ** 2 self.end_time = self.time_step * 1 - self.T_0 = 0 - self.s_0 = 1 + # Initial scalar value + self.initial_scalar = 0 + + # Whether or not to subtract the fracture pressure contribution for the contact + # traction. This should be done if the scalar variable is pressure, but not for + # temperature. See assign_discretizations + self.subtract_fracture_pressure = True def bc_type(self, g, key, t=0): if key == self.mechanics_parameter_key: @@ -190,8 +195,162 @@ def set_scalar_parameters(self): {"normal_diffusivity": normal_diffusivity}, ) - def assign_discretisations_and_variables(self): - assign_discretizations.contact_mechanics_and_biot_discretizations(self) + 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): """ @@ -217,7 +376,7 @@ def initial_condition(self): for g, d in self.gb: # Initial value for the scalar variable. - initial_scalar_value = self.T_0 * np.ones(g.num_cells) + initial_scalar_value = self.initial_scalar * np.ones(g.num_cells) d[pp.STATE].update({self.scalar_variable: initial_scalar_value}) if g.dim == self.Nd: bc_values = d[pp.PARAMETERS][self.mechanics_parameter_key]["bc_values"] @@ -254,35 +413,33 @@ def run_biot(setup, atol=1e-10): if "gb" not in setup.__dict__: setup.create_grid() gb = setup.gb + # Extract the grids we use - ambient_dim = gb.dim_max() - g_max = gb.grids_of_dimension(ambient_dim)[0] + g_max = gb.grids_of_dimension(setup.Nd)[0] d_max = gb.node_props(g_max) - # set parameters + # Assign parameters, variables and discretizations setup.set_parameters() setup.initial_condition() - - setup.assign_discretisations_and_variables() - # Define rotations - # Set up assembler and get 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] # Discretize with the biot class setup.discretize_biot(gb) + # Prepare for the time loop errors = [] - - t = 0.0 + t = getattr(setup, "time", 0) dt = setup.time_step - T = setup.end_time + t_end = setup.end_time k = 0 times = [t] - while t < T: + while t < t_end: t += dt k += 1 - print("Time step: ", k, "/", int(np.ceil(T / dt))) + print("Time step: ", k, "/", int(np.ceil(t_end / dt))) times.append(t) # Prepare for Newton diff --git a/src/porepy/numerics/fv/mpsa.py b/src/porepy/numerics/fv/mpsa.py index cabce8af2c..3ec1c02e8c 100644 --- a/src/porepy/numerics/fv/mpsa.py +++ b/src/porepy/numerics/fv/mpsa.py @@ -122,6 +122,7 @@ def discretize(self, g, data): ) 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 diff --git a/src/porepy/utils/assign_discretizations.py b/src/porepy/utils/assign_discretizations.py deleted file mode 100644 index e45bfb0cfa..0000000000 --- a/src/porepy/utils/assign_discretizations.py +++ /dev/null @@ -1,222 +0,0 @@ -""" -Convenience methods for assigning discretization methods for some common coupled -problems. - -""" -import porepy as pp -from porepy.utils.derived_discretizations import implicit_euler as IE_discretizations - - -def contact_mechanics_discretizations(setup): - """ - Setup should have a gb field, and the following names specified: - Parameter keys: - mechanics_parameter_key - friction_parameter_key - Variables: - displacement_variable - higher-dimensional displacements - mortar_displacement_variable - displacement on the internal boundary - contact_variable - represents traction on the fracture - """ - gb = setup.gb - ambient_dim = gb.dim_max() - # Define discretization - # For the Nd domain we solve linear elasticity with mpsa. - mpsa = pp.Mpsa(setup.mechanics_parameter_key) - empty_discr = pp.VoidDiscretization( - setup.friction_parameter_key, ndof_cell=ambient_dim - ) - - # Define discretization parameters - for g, d in gb: - if g.dim == ambient_dim: - d[pp.PRIMARY_VARIABLES] = { - setup.displacement_variable: {"cells": ambient_dim} - } - d[pp.DISCRETIZATION] = {setup.displacement_variable: {"mpsa": mpsa}} - elif g.dim == ambient_dim - 1: - d[pp.PRIMARY_VARIABLES] = {setup.contact_variable: {"cells": ambient_dim}} - d[pp.DISCRETIZATION] = {setup.contact_variable: {"empty": empty_discr}} - else: - d[pp.PRIMARY_VARIABLES] = {} - - # And define a Robin condition on the mortar grid - coloumb = pp.ColoumbContact(setup.friction_parameter_key, ambient_dim) - contact_discr = pp.PrimalContactCoupling( - setup.friction_parameter_key, mpsa, coloumb - ) - - for e, d in gb.edges(): - g_l, g_h = gb.nodes_of_edge(e) - - if g_h.dim == ambient_dim: - d[pp.PRIMARY_VARIABLES] = { - setup.mortar_displacement_variable: {"cells": ambient_dim} - } - - d[pp.COUPLING_DISCRETIZATION] = { - setup.friction_coupling_term: { - g_h: (setup.displacement_variable, "mpsa"), - g_l: (setup.contact_variable, "empty"), - (g_h, g_l): (setup.mortar_displacement_variable, contact_discr), - } - } - else: - d[pp.PRIMARY_VARIABLES] = {} - - -def contact_mechanics_and_biot_discretizations(setup, subtract_fracture_pressure=True): - """ - Assign the discretizations for fracture deformation with a coupled scalar (pressure) - in both dimensions. No fracture intersections are allowed (for now). - - Setup should have a gb field, and the following names specified: - Parameter keys: - mechanics_parameter_key - scalar_parameter_key - Variables: - displacement_variable - higher-dimensional displacements - mortar_displacement_variable - displacement on the internal boundary - contact_traction_variable - represents traction on the fracture - scalar_variable - scalar (pressure) in both dimensions - mortar_scalar_variable - darcy flux - subtract_fracture_pressure (bool): Whether or not to subtract the fracture pressure - contribution to the fracture force balance equation. This is needed for the - pressure case, where the forces on the fracture surfaces are the sum of the - contact force and the pressure force. It is not, however, needed for TM - simulations, where there is no force from the fracture temperature. - """ - gb = setup.gb - ambient_dim = gb.dim_max() - key_s, key_m = setup.scalar_parameter_key, setup.mechanics_parameter_key - var_s, var_d = setup.scalar_variable, setup.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=ambient_dim) - # 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=setup.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 gb: - if g.dim == ambient_dim: - d[pp.PRIMARY_VARIABLES] = { - var_d: {"cells": ambient_dim}, - var_s: {"cells": 1}, - } - d[pp.DISCRETIZATION] = { - var_d: {"mpsa": mpsa}, - var_s: { - "diffusion": diff_disc_s, - "mass": mass_disc_s, - "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 == ambient_dim - 1: - d[pp.PRIMARY_VARIABLES] = { - setup.contact_traction_variable: {"cells": ambient_dim}, - var_s: {"cells": 1}, - } - d[pp.DISCRETIZATION] = { - setup.contact_traction_variable: {"empty": empty_discr}, - var_s: { - "diffusion": diff_disc_s, - "mass": mass_disc_s, - "source": source_disc_s, - }, - } - else: - d[pp.PRIMARY_VARIABLES] = {} - - # Define edge discretizations for the mortar grid - contact_law = pp.ColoumbContact(setup.mechanics_parameter_key, ambient_dim) - contact_discr = pp.PrimalContactCoupling( - setup.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( - setup.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 subtract_fracture_pressure: - fracture_scalar_to_force_balance = pp.FractureScalarToForceBalance( - mass_disc_s, mass_disc_s - ) - for e, d in gb.edges(): - g_l, g_h = gb.nodes_of_edge(e) - - if g_h.dim == ambient_dim: - d[pp.PRIMARY_VARIABLES] = { - setup.mortar_displacement_variable: {"cells": ambient_dim}, - setup.mortar_scalar_variable: {"cells": 1}, - } - - d[pp.COUPLING_DISCRETIZATION] = { - setup.friction_coupling_term: { - g_h: (var_d, "mpsa"), - g_l: (setup.contact_traction_variable, "empty"), - (g_h, g_l): (setup.mortar_displacement_variable, contact_discr), - }, - setup.scalar_coupling_term: { - g_h: (var_s, "diffusion"), - g_l: (var_s, "diffusion"), - e: ( - setup.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: (setup.mortar_displacement_variable, div_u_coupling), - }, - "matrix_scalar_to_force_balance": { - g_h: (var_s, "mass"), - g_l: (var_s, "mass"), - e: ( - setup.mortar_displacement_variable, - matrix_scalar_to_force_balance, - ), - }, - } - if subtract_fracture_pressure: - d[pp.COUPLING_DISCRETIZATION].update( - { - "matrix_scalar_to_force_balance": { - g_h: (var_s, "mass"), - g_l: (var_s, "mass"), - e: ( - setup.mortar_displacement_variable, - fracture_scalar_to_force_balance, - ), - } - } - ) - else: - raise ValueError( - "assign_discretizations assumes no fracture intersections." - )