porepy.numerics.fv.biot module

Modules contains discretization of poro-elasticity by the multi-point stress approximation.

The discretization scheme is described in

J.M. Nordbotten (2016): STABLE CELL-CENTERED FINITE VOLUME DISCRETIZATION FOR BIOT

EQUATIONS, SINUM.

The module contains four classes: Biot is the main class, responsible for discretization of the poro-elastic system (relying heavily on its parent class pp.Mpsa). This is a proper discretization class, in that it has discretize() and assemble_matrix_rhs() methods.

The other classes, GradP, DivU and BiotStabilization, are containers for the other terms in the poro-elasticity system. These have empty discretize() methods, and should not be used directly.

Because of the number of variables and equations, and their somewhat difficult relation, the most convenient way to set up a discretization for poro-elasticity is to use pp.BiotContactMechanicsModel (designed for fractures and contact mechanics, but will turn into a standard poro-elasticity equation for non-fractured domains).

class Biot(mechanics_keyword='mechanics', flow_keyword='flow', vector_variable='displacement', scalar_variable='pressure')[source]

Bases: Mpsa

Discretization class for poro-elasticity, based on MPSA.

This is a subclass of pp.Mpsa()

Parameters
  • mechanics_keyword (str) –

  • flow_keyword (str) –

  • vector_variable (str) –

  • scalar_variable (str) –

mechanics_keyword

Keyword used to identify the parameter dictionary associated with the mechanics subproblem. Defaults to “mechanics”.

Type

str

flow_keyword

Keyword used to identify the parameter dictionary associated with the flow subproblem. Defaults to “flow”.

Type

str

vector_variable

Name for the vector variable, used throughout the discretization and solution. Defaults to “displacement”.

Type

str

scalar_variable

Name for the vector variable, used throughout the discretization and solution. Defaults to “pressure”.

Type

str

div_u_matrix_key

Keyword used to identify the discretization matrix of the term div(u). Defaults to “div_u”.

Type

str

bound_div_u_matrix_key

Keyword used to identify the discretization matrix of the boundary condition for the term div(u). Defaults to “bound_div_u”.

Type

str

grad_p_matrix_key

Keyword used to identify the discretization matrix of the term grad(p). Defaults to “grad_p”.

Type

str

stabilization_matrix_key

Keyword used to identify the discretization matrix for the stabilization term. Defaults to “biot_stabilization”.

Type

str

bound_pressure_matrix_key

Keyword used to identify the discretization matrix for the pressure contribution to boundary displacement reconstruction. Defaults to “bound_displacement_pressure”.

Type

str

See also pp.Mpsa for further attributes.

assemble_matrix(sd, sd_data)[source]

Assemble the poro-elastic system matrix.

The discretization is presumed stored in the data dictionary.

Parameters
  • sd (pp.Grid) – Grid for discretization

  • sd_data (dictionary) – Data for discretization, as well as matrices with discretization of the sub-parts of the system.

Returns

Block matrix with the combined MPSA/MPFA

discretization.

Return type

scipy.sparse.bmat

Raises

NotImplementedError – If invoked. The method is included to be compatible with the wider discretization class, but assembly should be handled by the ContactMechanicsBiot class.

assemble_matrix_rhs(sd, sd_data)[source]

Method inherited from superclass should not be used.

Parameters
  • sd (pp.Grid) – Grid to be discretized.

  • sd_data (dict) – Data dictionary for this grid.

Returns

System matrix of this discretization. The size of

the matrix will depend on the specific discretization.

np.ndarray: Right-hand side vector with representation of boundary

conditions. The size of the vector will depend on the discretization.

Return type

scipy.sparse.csr_matrix

Raises

NotImplementedError – If invoked. The method is included to be compatible with the wider discretization class, but assembly should be handled by the ContactMechanicsBiot class.

assemble_rhs(sd, sd_data)[source]

Return the right-hand side for a poro-elastic system.

Parameters
  • sd (pp.Grid) – Grid to be discretized.

  • sd_data (dict) – Data dictionary for this grid.

Returns

Right hand side vector with representation of boundary

conditions. The size of the vector will depend on the discretization.

Return type

np.ndarray

Raises

NotImplementedError – If invoked. The method is included to be compatible with the wider discretization class, but assembly should be handled by the ContactMechanicsBiot class.

compute_flux(sd, u, sd_data)[source]

Compute flux field corresponding to a solution.

Parameters
  • sd (Grid) – grid, or a subclass.

  • u (np.ndarray) – Solution variable, representing displacements and pressure.

  • sd_data (dictionary) – dictionary related to grid and problem. Should contain boundary discretization.

Returns

Flux over all faces

Return type

np.ndarray

compute_stress(sd, u, sd_data)[source]

Compute stress field corresponding to a solution.

Parameters
  • sd (Grid) – grid, or a subclass.

  • u (np.ndarray) – Solution variable, representing displacements and pressure.

  • sd_data (dictionary) – dictionary related to grid and problem. Should contain boundary discretization.

Returns

Stress over all faces. Stored as

all stress values on the first face, then the second etc.

Return type

np.ndarray, sd.dim * sd.num_faces

discretize(sd, sd_data)[source]

Discretize the mechanics terms in a poromechanical system.

NOTE: This function does not discretize purely flow-related terms (Darcy flow and compressibility).

The parameters needed for the discretization are stored in the dictionary sd_data, which should contain the following mandatory keywords:

Related to mechanics equation (in sd_data[pp.PARAMETERS][self.mechanics_keyword]):

fourth_order_tensor: Fourth order tensor representing elastic moduli. bc: BoundaryCondition object for mechanics equation.

Used in mpsa.

In addition, the following parameters are optional:
Related to coupling terms:
biot_alpha (double between 0 and 1): Biot’s coefficient.

Defaults to 1.

Related to numerics:
inverter (str): Which method to use for block inversion. See

pp.fvutils.invert_diagonal_blocks for detail, and for default options.

mpfa_eta (double): Location of continuity point in MPSA.

Defaults to 1/3 for simplex grids, 0 otherwise.

The discretization is stored in the data dictionary, in the form of several matrices representing different coupling terms. For details, and how to combine these, see self.assemble_matrix()

Parameters
  • sd (pp.Grid) – Grid to be discretized.

  • sd_data (dictionary) – Containing data for discretization. See above for specification.

Return type

None

extract_scalar(sd, u)[source]

Extract pressure field from solution.

Parameters
  • sd (Grid) – grid, or a subclass.

  • u (np.ndarray) – Solution variable, representing displacements and pressure.

Returns

Pressure part of solution vector.

Return type

np.ndarray

extract_vector(sd, u, dims=None, as_vector=False)[source]

Extract displacement field from solution.

Parameters
  • sd (Grid) – grid, or a subclass.

  • u (np.ndarray) – Solution variable, representing displacements and pressure.

  • dims (list of int, optional) – Which dimension to extract. If None, all dimensions are returned.

Returns

Displacement variables in the specified

dimensions.

Return type

list of np.ndarray

ndof(sd)[source]

Return the number of degrees of freedom associated to the method.

Parameters

sd (pp.Grid) – A grid.

Returns

The number of degrees of freedom.

Return type

int

update_discretization(sd, sd_data)[source]

Update discretization.

The updates can generally come as a combination of two forms:
  1. The discretization on part of the grid should be recomputed.

  2. The old discretization can be used (in parts of the grid), but the numbering of unknowns has changed, and the discretization should be reordered accordingly.

Information on the basis for the update should be stored in a field

sd_data[‘update_discretization’]

This should be a dictionary which could contain keys:

modified_cells, modified_faces

define cells, faces and nodes that have been modified (either parameters, geometry or topology), and should be rediscretized. It is up to the discretization method to implement the change necessary by this modification. Note that depending on the computational stencil of the discretization method, a grid quantity may be rediscretized even if it is not marked as modified.

The dictionary sd_data[‘update_discretization’] should further have keys:

cell_index_map, face_index_map

these should specify sparse matrices that map old to new indices. If not provided, the cell and face bookkeeping will be assumed constant.

Parameters
  • sd (pp.Grid) – Grid to be rediscretized.

  • sd_data (dictionary) – With discretization parameters.

class BiotStabilization(keyword='mechanics', variable='pressure')[source]

Bases: Discretization

Class for the stabilization term of the Biot equation.

assemble_matrix(sd, sd_data)[source]

Return the matrix and right-hand side for a discretization of the stabilization term of the Biot equation.

Parameters
  • sd (pp.Grid) – Computational grid, with geometry fields computed.

  • sd_data (dictionary) – With data stored.

Returns

System matrix of this discretization. The

size of the matrix will depend on the specific discretization.

Return type

scipy.sparse.csr_matrix

Raises
  • ValueError if the stabilization term has not already been

  • discretized.

assemble_matrix_rhs(sd, sd_data)[source]

Return the matrix and right-hand side for a discretization of the stabilization term of the Biot equation.

Parameters
  • sd (pp.Grid) – grid, or a subclass, with geometry fields computed.

  • sd_data (dict) – dictionary to store the data. For details on necessary keywords, see method discretize()

Returns

sparse csr (sd.dim * g.num_cells, sd.dim * sd.num_cells)

Discretization matrix.

rhs: array (sd.dim * sd.num_cells) Right-hand side.

Return type

matrix

assemble_rhs(sd, sd_data)[source]

Return the right-hand side for the stabilization part of the displacement divergence term.

For the time being, we assume an IE temporal discretization.

Parameters
  • sd (pp.Grid) – Computational grid.

  • sd_data (dictionary) – With data stored.

Returns

Zero right-hand side vector with representation of boundary

conditions.

Return type

np.ndarray

discretize(sd, sd_data)[source]

Discretize the stabilization term of the Biot equation.

Parameters
  • g (pp.Grid) – grid, or a subclass, with geometry fields computed.

  • sd_data (dict) – Data dictionary.

  • sd (Grid) –

Raises
  • NotImplementedError, the discretization should be performed using the

  • discretize method of the Biot class.

ndof(sd)[source]

Return the number of degrees of freedom associated to the method.

Parameters

sd (pp.Grid) – A grid.

Returns

The number of degrees of freedom.

Return type

int

class DivU(mechanics_keyword='mechanics', flow_keyword='flow', variable='displacement', mortar_variable='mortar_displacement')[source]

Bases: Discretization

Class for the displacement divergence term of the Biot equation.

Parameters
  • mechanics_keyword (str) –

  • flow_keyword (str) –

  • variable (str) –

  • mortar_variable (str) –

assemble_int_bound_displacement_source(sd, sd_data, intf, intf_data, cc, matrix, rhs, self_ind)[source]

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
  • sd (pp.Grid) – Grid which the condition should be imposed on.

  • sd_data (dictionary) – Data dictionary for the node in the mixed-dimensional grid.

  • intf_data (dictionary) – Data dictionary for the edge in the mixed-dimensional grid.

  • cc (block matrix, 3x3) – Block matrix for the coupling condition. The first and second rows and columns are identified with the primary and secondary 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 0 or 1.

  • intf (MortarGrid) –

assemble_int_bound_displacement_trace(sd, sd_data, intf, intf_data, grid_swap, cc, matrix, rhs, self_ind)[source]

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 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.

Parameters
  • sd (pp.Grid) – Grid which the condition should be imposed on.

  • sd_data (dictionary) – Data dictionary for the node in the mixed-dimensional grid.

  • intf_data (dictionary) – Data dictionary for the edge in the mixed-dimensional grid.

  • grid_swap (boolean) – If True, the grid g is identified with the @ secondary side of the mortar grid in intf_data.

  • cc (block matrix, 3x3) – Block matrix for the coupling condition. The first and second rows and columns are identified with the primary and secondary 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.

  • intf (MortarGrid) –

  • rhs (ndarray) –

assemble_matrix(sd, sd_data)[source]

Return the matrix and right-hand side for a discretization of the displacement divergence term of the Biot equation.

Parameters
  • sd (pp.Grid) – Computational grid, with geometry fields computed.

  • sd_data (dictionary) – With data stored.

Returns

System matrix of this discretization. The

size of the matrix will depend on the specific discretization.

Return type

scipy.sparse.csr_matrix

Raises
  • ValueError if the displacement divergence term has not already been

  • discretized.

assemble_matrix_rhs(sd, sd_data)[source]

Return the matrix and right-hand side for a discretization of the displacement divergence term of the Biot equation.

Parameters
  • sd (pp.Grid) – grid, or a subclass, with geometry fields computed.

  • sd_data (dict) – dictionary to store the data. For details on necessary keywords, see method discretize()

Returns

sparse csr (sd.dim * g_num_cells, sd.dim * g_num_cells) Discretization matrix. rhs: array (sd.dim * g_num_cells) Right-hand side.

Return type

matrix

assemble_rhs(sd, sd_data)[source]

Return the right-hand side for a discretization of the displacement divergence term.

For the time being, we assume an IE temporal discretization.

Parameters
  • sd (pp.Grid) – Computational grid.

  • sd_data (dictionary) – With data stored.

Returns

Zero right-hand side vector with representation of boundary

conditions.

Return type

np.ndarray

discretize(sd, sd_data)[source]

Discretize the displacement divergence term of the Biot equation.

Parameters
  • g (pp.Grid) – grid, or a subclass, with geometry fields computed.

  • sd_data (dict) – For entries, see above.

  • sd (Grid) –

Raises
  • NotImplementedError, the discretization should be performed using the

  • discretize method of the Biot class.

ndof(sd)[source]

Return the number of degrees of freedom associated to the method.

Parameters

sd (pp.Grid) – A grid.

Returns

The number of degrees of freedom.

Return type

int

class GradP(keyword)[source]

Bases: Discretization

Class for the pressure gradient term of the Biot equation.

Parameters

keyword (str) –

assemble_matrix(sd, sd_data)[source]

Return the matrix and right-hand side for a discretization of the pressure gradient term of the Biot equation.

Parameters
  • sd (pp.Grid) – Computational grid, with geometry fields computed.

  • sd_data (dictionary) – With data stored.

Returns

System matrix of this discretization. The

size of the matrix will depend on the specific discretization.

Return type

scipy.sparse.csr_matrix

Raises

ValueError if the pressure gradient term has not already been discretized.

assemble_matrix_rhs(sd, sd_data)[source]

Return the matrix and right-hand side for a discretization of the pressure gradient term of the Biot equation.

Parameters
  • sd (pp.Grid) – grid, or a subclass, with geometry fields computed.

  • sd_data (dict) – dictionary to store the data. For details on necessary keywords, see method discretize().

Returns

sparse csr (sd.dim * g_num_cells, sd.dim * sd.num_cells).

Discretization matrix.

rhs: array (sd.dim * sd.num_cells) Right-hand side.

Return type

matrix

assemble_rhs(sd, sd_data)[source]

Return the right-hand side for a discretization of the pressure gradient term.

The contribution is nonzero only for nonzero p_reference, which mainly is of relevance for thermal stress on the form

thermal_expansion * (T-T_reference) * I.

Parameters
  • sd (pp.Grid) – Computational grid.

  • sd_data (dictionary) – With data stored.

Returns

Right-hand side vector with representation of reference pressure

contribution.

Return type

np.ndarray

discretize(sd, sd_data)[source]

Discretize the pressure gradient term of the Biot equation.

Parameters
  • sd (pp.Grid) – grid, or a subclass, with geometry fields computed.

  • sd_data (dict) – For entries, see above.

Raises
  • NotImplementedError, the discretization should be performed using the

  • discretize method of the Biot class.

ndof(sd)[source]

Return the number of degrees of freedom associated to the method.

Parameters

sd (pp.Grid) – A grid.

Returns

The number of degrees of freedom.

Return type

int