fvbiot: Consistent finite-volume discretizations for poroelasticity

This package implements a discretization of poro-mechanics by cell centered finite volume methods.

The package provides discretization of three different equations:

1. Scalar elliptic equations (Darcy flow), using Multipoint flux approximations.
2. Linear elasticity, using Multipoint stress approximations
3. Poro-mechanics, e.g. coupling terms for the combined system of 1. and 2.

Examples of usage for the three options can be found in mpfa_ex, mpsa_ex and biot_ex, respectively.

Further information on the methods implemented herein, in particular the discretization of elasticity and poro-mechanics, can be found in the following papers (all open access):

10. M. Nordbotten: Stable cell-center finite volume disrcetizaiton for Biot equations - SIAM J. Numer. Anal. 54(2) 942-968.

10. M. Nordbotten: Convergence of a cell-centered finite volume discretization for linear elasticity - SIAM J. Numer. Anal. 53(2) 2605-2625.

5. Keilegavlen, J. M. Nordbotten: Finite volume methods for elasticity with weak symmetry - Int. J. Num. Meth. Eng, 2017, doi: 10.1002/nme.5538.

The implementation of elasticity is based on the forumlation described in the latter. For information on the MPFA discretizations for the flow equation, see

1. Aavatsmark: An introduction to multipoint flux approximations for quadrilateral grids, Comput. Geosci. 6(3-4) 405-432.

The is compatible with the (freely available) Matlab Reservoir Simulation Toolbox (MRST) provided by SINTEF ICT, see http://www.sintef.no/projectweb/mrst/ In particular, the grid data structure is that of MRST, and it is assumed that MRST is in the matlab path.

The code has been tested with Matlab R2015b, and mrst-2016a.

mpfa(G, rock, activeNodes, varargin)

Discretize Darcy’s law by a multipoint flux approximation

Parameters: G - mrst grid structure, see http://www.sintef.no/projectweb/mrst/

To get an overview of the data format, type ‘help grid_structure’

rock - mrst rock structure. Should contain a field perm. To get an

overview of the data format, type ‘help makeRock’

Active nodes: Which nodes should the stencil be computed for. If empty,

all nodes in the grid is considered.

NOTE: The function computes and stores in memory the inverse of a block-diagonal matrix (which expresses gradients on the sub-cells as a function of pressures in the surrounding cells). The size of each sub-matrix will be the number of sub-cells sharing the vertex * number of dimensions, so each system will be 8x8 for Cartesian 2D, 24x24 for Cartesian 3D. For large 3D grids, in particular with simplices, this may be heavy in terms of memory needs. For this reason, the activeNodes option may be an atractive option. For a work-around that splits the discretization into several calls to mpfa, and thus reduces memory consumption, see mpfa_subgrid.

Optional parametrs, defined as keyword - value pairs: ‘eta’ - Location of continuity point on the half edges, defined

according to Aavatsmark 2002. Between 0 and 1. Default value is 0, for simplices, the value should be 1/3 (will give symmetric method, see Klausen, Radu Eigestad 2008).

bc - boundary conditions, as defined by the MRST function addBC. If none

are provided, homogeneous Neumann conditions will be assigned.

invertBlocks - method for inverting block diagonal systems. Should be

either ‘matlab’ or ‘mex’. The former is pure matlab, which will be slow for large problems. Mex is substantially faster, but requires that matlab has access to a mex / c compiler.

returnHalfEdgeValues - whether the return values for stresses etc should

be for sub edges/faces (or summed to the whole face). Default is the full faces. Can be used for a partial update of the discretization scheme.

mpfa_subgrid(G, rock, mode, varargin)

Discretize Darcy’s law by a multipoint flux approximation.

This wrapper allows for discretizaiton of parts of the domain, as opposed to the full discretization computed by mpfa.m So far, the only version implemented handles the case where the computation in memory bound (see comment in mpfa.m). A similar approach could be used if the permeability or geometry has changed in part of the domain, however, at the moment this has not been implemented.

Parameters: G - mrst grid structure, see http://www.sintef.no/projectweb/mrst/

To get an overview of the data format, type ‘help grid_structure’

rock - mrst rock structure. Should contain a field perm. To get an

overview of the data format, type ‘help makeRock’

mode - basically a parameter describing the reason for calling this

function and not mpfa.m. At the moment, this parameter defaults to 1, which is the memory bound case

Optional parametrs, defined as keyword - value pairs: ‘eta’ - Location of continuity point on the half edges, defined

according to Aavatsmark 2002. Between 0 and 1. Default value is 0, for simplices, the value should be 1/3 (will give symmetric method, see Klausen, Radu Eigestad 2008).

bc - boundary conditions, as defined by the MRST function addBC. If none

are provided, homogeneous Neumann conditions will be assigned.

invertBlocks - method for inverting block diagonal systems. Should be

either ‘matlab’ or ‘mex’. The former is pure matlab, which will be slow for large problems. Mex is substantially faster, but requires that matlab has access to a mex / c compiler.

maxMemory - maximum memory assigned to storing the inverse gradient

operator in mpfa. Note that although this gives a substantial part of the total memory consumption, in particular for large grids, maxMemory does NOT control total memory usage.

mpsa(G, C, activeNodes, varargin)

Discretize linear elasticity and Biot by a multipoint stress approximation.

The implementation is modelled on mpfa.m, but everything becomes somewhat more complex for the vector unknown.

NOTE: For poro-mechanical problems (not pure elasticity), non-homogeneous boundary conditions have not yet been implemented. This applies to the fields.

Parameters: G - mrst grid structure, see http://www.sintef.no/projectweb/mrst/

To get an overview of the data format, type ‘help grid_structure’

C - Constitutive relation. See shear_normal_stress.m for details. Active nodes: Which nodes should the stencil be computed for. If empty,

all nodes in the grid is considered.

NOTE: The function computes and stores in memory the inverse of a block-diagonal matrix (which expresses gradients on the sub-cells as a function of pressures in the surrounding cells). The size of each sub-matrix will be the number of sub-cells sharing the vertex * number of dimensions^2, so each system will be 16x16 for Cartesian 2D, 72x72 for Cartesian 3D. For large 3D grids, in particular with simplices, this may be heavy in terms of memory needs. For this reason, the activeNodes option may be an atractive option. For a work-around that splits the discretization into several calls to mpfa, and thus reduces memory consumption, see mpsa_subgrid.

Optional parametrs, defined as keyword - value pairs: ‘eta’ - Location of continuity point on the half edges, defined

according to Aavatsmark 2002. Between 0 and 1. Default value is 0, for simplices, the value should be 1/3 (will give symmetric method, see Klausen, Radu, Eigestad 2008).

bc - boundary conditions, as defined by the MRST function addBC. If none

are provided, homogeneous Neumann conditions will be assigned.

invertBlocks - method for inverting block diagonal systems. Should be

either ‘matlab’ or ‘mex’. The former is pure matlab, which will be slow for large problems. Mex is substantially faster, but requires that matlab has access to a mex / c compiler.

returnHalfEdgeValues - whether the return values for stresses etc should

be for sub edges/faces (or summed to the whole face). Default is the full faces. Can be used for a partial update of the discretization scheme.

Returns:

out –

a structure containing the following fields:

Discretization of elasticity:

A - discretization matrix div - divergence operator stress - discretized Hook’s law boundStress - Discretization of boundary conditions

Fields used for discretization of poro-mechanics

divD - Discretized volumetric stress operator gradP - Discrete pressure gradient for fluid pressure in Biot’s

equations

stabDelta - Stabilitazion term used for the flow equation.

Example: See mpsa_ex and biot_ex

mpsa_subgrid(G, stiffness, mode, varargin)

Discretize linear elasticity by a multipoint stress approximation.

This wrapper allows for discretizaiton of parts of the domain, as opposed to the full discretization computed by mpsa.m So far, the only version implemented handles the case where the computation in memory bound (see comment in mpsa.m). A similar approach could be used if the permeability or geometry has changed in part of the domain, however, at the moment this has not been implemented.

Parameters: G - mrst grid structure, see http://www.sintef.no/projectweb/mrst/

To get an overview of the data format, type ‘help grid_structure’

stiffness - stiffness matrix, as created by by shear_normal_stress. mode - basically a parameter describing the reason for calling this

function and not mpfa.m. At the moment, this parameter defaults to 1, which is the memory bound case

Optional parametrs, defined as keyword - value pairs: ‘eta’ - Location of continuity point on the half edges, defined

according to Aavatsmark 2002. Between 0 and 1. Default value is 0, for simplices, the value should be 1/3 (will give symmetric method, see Klausen, Radu Eigestad 2008).

bc - boundary conditions, as defined by the MRST function addBC. If none

are provided, homogeneous Neumann conditions will be assigned.

invertBlocks - method for inverting block diagonal systems. Should be

either ‘matlab’ or ‘mex’. The former is pure matlab, which will be slow for large problems. Mex is substantially faster, but requires that matlab has access to a mex / c compiler.

maxMemory - maximum memory assigned to storing the inverse gradient

operator in mpfa. Note that although this gives a substantial part of the total memory consumption, in particular for large grids, maxMemory does NOT control total memory usage.

shear_normal_stress(Nc, Nd, mu, lambda, phi)

Define stiffness matrix for linear elastic medium.

Parameters Nc: Number of cells in the grid Nd: Number of dimensions mu: First lame parameter, cellwise (Nc x 1 vector) lambda: Second Lame parameter, cellwise (Nc x 1 vector) phi: Extra parameter, never been used, stay away unless you know what you

are doing

Returns:Constitutive relation, one matlab cell per computational cell.

Examples

Example setting up and solving a simple 2D poro-mechanics problem.
% See also mpfa_ex and mpsa_ex for more detailed comments.
%
% The example assumes MRST is the Matlab path. For information on
% MRST-functions, confer the MRST documentation at
%   http://www.sintef.no/projectweb/mrst

Example setting up and solving a simple 2D flow problem using MPFA.
%
% The example assumes MRST is the Matlab path. For information on
% MRST-functions, confer the MRST documentation at
%   http://www.sintef.no/projectweb/mrst

Example setting up and solving a simple 2D elasticity problem using MPSA.
% Much of the syntax is similar to that used for the flow problem
% (mpfa_ex.m), see there for more information.
%
% The example assumes MRST is the Matlab path. For information on
% MRST-functions, confer the MRST documentation at
%   http://www.sintef.no/projectweb/mrst