Classes | Public Types | Public Member Functions | List of all members
Opm::IncompFlowSolverHybrid< GridInterface, RockInterface, BCInterface, InnerProduct > Class Template Reference

Solve mixed formulation of incompressible flow modelled by Darcy's law

\[@f{aligned}{ v &= -K\lamda(\nabla p + \rho\vec{g}), \\ \nabla\cdot v &= q. \end

] The solver is based on a hybrid discretization scheme which yields a system of linear equations of the form

\[@f{bmatrix}{ B & C & D \\ C^{T} & 0 & 0 \\ D^{T} & 0 & 0 }\, @f{bmatrix}{v \\ -p \\ \pi\end

=

\begin{bmatrix}f \\ g \\ h\end{bmatrix}

] where $v$ represents the interface fluxes for each interface for each cell, $p$ are the cell pressures and $\pi$ are the interface pressures. More...

#include <IncompFlowSolverHybrid.hpp>

Inheritance diagram for Opm::IncompFlowSolverHybrid< GridInterface, RockInterface, BCInterface, InnerProduct >:
Inheritance graph

Public Types

typedef const FlowSolution & SolutionType
 Type representing the solution to the problem defined by the parameters to. More...
 

Public Member Functions

template<class Point >
void init (const GridInterface &g, const RockInterface &r, const Point &grav, const BCInterface &bc)
 All-in-one initialization routine. Enumerates all grid connections, allocates sufficient space, defines the structure of the global system of linear equations for the contact pressures, and computes the permeability dependent inner products for all of the grid's cells. More...
 
void clear ()
 Clear all topologic, geometric and rock-dependent information currently held in internal data structures. Release all memory resources in preparation of constructing a solver for a new problem. Method. More...
 
void initSystemStructure (const GridInterface &g, const BCInterface &bc)
 Compute structure of coefficient matrix in final system of linear equations for this flow problem. Allocates all storage needed by the flow solver itself. More...
 
template<class Point >
void computeInnerProducts (const RockInterface &r, const Point &grav)
 Compute static (i.e., independent of saturation) parts of the spatially varying inner products $ (v,K^{-1}w) $ for each grid cell. This is often a fairly expensive process, so in order to increase the speed of the system assembly we do not wish to perform this computation for each time step unless we have to. Moreover, if the permeability field changes but the grid connections remain static, we do not have to build a new matrix structure but need only compute new values for these static inner products. More...
 
template<class FluidInterface >
void solve (const FluidInterface &r, const std::vector< double > &sat, const BCInterface &bc, const std::vector< double > &src, double residual_tolerance=1e-8, int linsolver_verbosity=1, int linsolver_type=1, bool same_matrix=false, int linsolver_maxit=0, double prolongate_factor=1.6, int smooth_steps=1)
 Construct and solve system of linear equations for the pressure values on each interface/contact between neighbouring grid cells. Recover cell pressure and interface fluxes. Following a call to. More...
 
double postProcessFluxes ()
 Postprocess the solution fluxes. This method modifies the solution object so that out-fluxes of twin faces (that is, the two faces on a cell-cell intersection) will be made antisymmetric. More...
 
SolutionType getSolution ()
 Recover the solution to the problem defined by the parameters to method. More...
 
template<typename charT , class traits >
void printStats (std::basic_ostream< charT, traits > &os)
 Print statistics about the connections in the current model. This is mostly for debugging purposes and should only rarely be used in client code. More...
 
void printSystem (const std::string &prefix)
 Output current system of linear equations to permanent storage in files. One file for the coefficient matrix and one file for the right hand side. This is mostly useful whilst debugging the solver and should rarely be used from client code. More...
 

Detailed Description

template<class GridInterface, class RockInterface, class BCInterface, template< class GridIF, class RockIF > class InnerProduct>
class Opm::IncompFlowSolverHybrid< GridInterface, RockInterface, BCInterface, InnerProduct >

Solve mixed formulation of incompressible flow modelled by Darcy's law

\[@f{aligned}{ v &= -K\lamda(\nabla p + \rho\vec{g}), \\ \nabla\cdot v &= q. \end

] The solver is based on a hybrid discretization scheme which yields a system of linear equations of the form

\[@f{bmatrix}{ B & C & D \\ C^{T} & 0 & 0 \\ D^{T} & 0 & 0 }\, @f{bmatrix}{v \\ -p \\ \pi\end

=

\begin{bmatrix}f \\ g \\ h\end{bmatrix}

] where $v$ represents the interface fluxes for each interface for each cell, $p$ are the cell pressures and $\pi$ are the interface pressures.

Through a Schur complement analysis, the above system of linear equations is reduced to an equivalent system of the form

\[@f{bmatrix}{ B & C & D \\ 0 & -L & -F \\ 0 & 0 & S }\, @f{bamtrix}{v \\ -p \\ \pi\end{bmatrix}

=

\begin{bmatrix}f \\ \hat{g} \\ r\end{bmatrix}

] where $L = C^{T}B^{-1}C$, $F=C^{T}B^{-1}D$, and $S = D^{T}B^{-1}D - F^{T}L^{-1}F$. Similarly, $\hat{g} = g - C^{T}B^{-1}f$ and $r = D^{T}B^{-1}f + F^{T}L^{-1}\hat{g} - h$.

The system $S\pi = r$ is the system of linear equations which is actually solved using separte linear system solver software. Then, a back substitution process yields the cell pressures and interface fluxes by solving the simpler systems

\[@f{aligned}{ Lp &= \hat{g} + F\pi \\ Bv &= f + Cp - D\pi \end{bmatrix}

] Specifically, the matrix $L$ is diagonal (a single scalar per grid cell) and the matrix $B^{-1}$ is available as a bi-product of the Schur reduction process. Finally, in the case of mimetic discretizations using the

MimeticIPEvaluator

class, the matrix $B^{-1}$ is directly available through an explicit formula.

The matrix $B$ is a product of two parts–one which depends only on the grid (i.e., the topology and geomtry) and the geological properties (i.e., the permeability field) and one part which depends only on the temporally varying fluid data (i.e., viscosities and saturation dependent relative permeabilities). Furthermore, the matrices $C$ and $D$ depend only on the grid topolgy. Consequently, the cost of assembling the system $S$ may be reduced, at least for simulations involving multiple pressure solves, by performing all of the static (i.e., the parts dependent only on the grid and geological properties) assembly once, at system initalization. Then the system assembly process consists of calculating the saturation dependent contributions and adding these into the system of linear equations.

The code is structured in a manner similar to traditional FEM software implementations. In particular, we assemble the coefficient matrix $S$ and system right hand side $r$ on a cell-by-cell basis. This yields a number of important simplifications. In particular, the matrix $D$ reduces to the identity on a single cell. Its effects are produced only through a local-to-global map of degrees of freedom. Similarly, on a single grid cell, the matrix $C=[1,1,\dots,1]^T$ with the number of elements equal to the number of interfaces (i.e., neighbours) of the cell. Finally, the matrix $B^{-1}$ is computed as

\[B^{-1} = \hat{B}^{-1}/\lambda_T\]

where $\hat{B}^{-1}$ denotes the static part of $B^{-1}$ and $\lambda_T$ denotes the total fluid mobility in the grid cell. For ease of implementation of the back substitution process, we compute and store the values of $L$, $F$ and $\hat{g}$ during each system assembly process.

A final quirk of the implementation is that we always solve for every interface pressure, even if a given pressure value is known through, e.g., a prescribed pressure boundary condition. This feature enables changing the type and value of a set of boundary conditions between pressure solves without having to reconstruct the sparsity structure of the coefficient matrix $S$–at least while no connections are introduced or removed. Periodic boundary condtions introduce new connections (i.e., new non-zero coefficient matrix entries) so in order to gain maximum efficiency, you should initialize the flow solver system only after you know the existence of all periodic boundary conditions which will be employed in a given simulation.

The use of non-periodic boundary conditions after periodic connections have been introduced is fully supported. This mode introduces a performance caveat. Specifically, a few entries which would otherwise not have been entered into the matrix due to being automatically zero will become explicit zeros instead. Consequently, each iteration of an iterative solver for the system $S\pi = r$ becomes slightly more expensive as explicit zero-multiplications occur. Moreover, introducing explicit non-zero representation of matrix entries which are always zero increases the demand for memory resources.

The flow solver is initialized in a three-step process. The first step, represented by method

clear()

, releases any previously held data in the internal data structures and prepares the solver system for defining a new problem. The second step, represented by the

initSystemStructure()

enumerates the primary degrees of freedom (i.e., the interface pressure values) and determines, and allocates, the coefficient matrix non-zero structure. The third and final step, represented by method

computeInnerProducts()

, computes the static (i.e., geology and geomtry-dependent) inner product matrices $\hat{B}^{-1}$. Method

computeInnerProducts()

is offered separately in order to support solve several different property models on the same underlying geometry (grid). Finally, method

init()

is a convenience method which calls the other initialization methods in sequence.

Following solver intialization, a sequence of flow problems differing only by boundary condition type/value and/or differing fluid saturation values may be resolved by separate calls to the

solve()

method. At any time following a call to

solve()

may the solution, represented by the type

SolutionType

, to the most recently resolved flow problem be retrieved through the

getSolution()

method. We note that

SolutionType

is always a reference-to-const.

Template Parameters
GridInterfaceType presenting an interface to a grid (typically a discretized geological model). The type is assumed to expose a forward iterator type,
CellIter
, and a pair of delimiters
cellbegin()
and
cellend()
to traverse the cells of a grid. The cell iterator, in turn, is expected to expose a forward iterator type
FaceIter
, and a pair of delimiters
facebegin()
and
faceend()
to traverse the faces of a cell.
RockInterfaceType presenting an interface to reservoir properties such as permeability and porosity. The type is assumed to expose a method,
permeability()
through which the assigned permeability of a single grid cell, represented by a
GridInterface::CellIter
, may be recovered.
BCInterfaceType presenting an interface to boundary conditions. The type is expected to provide a method,
flowCond()
, from which a
BCInterface::FlowBC
boundary condtion type of any given boundary face may be recovered. The flow boundary condition type is expected to provide methods for quering the type of boundary condition (i.e., prescribed pressure values, prescribed flux values, or prescribed pressure drops in the case of periodic boundary conditions) as well as the numerical values of these conditions.
InnerProductType presenting a specific inner product defining a discretization of the Darcy equation.

Member Typedef Documentation

template<class GridInterface, class RockInterface, class BCInterface, template< class GridIF, class RockIF > class InnerProduct>
typedef const FlowSolution& Opm::IncompFlowSolverHybrid< GridInterface, RockInterface, BCInterface, InnerProduct >::SolutionType

Type representing the solution to the problem defined by the parameters to.

solve()

. Always a reference-to-const. The

SolutionType

exposes methods

pressure(c)

and

outflux(f)

from which the cell pressure in cell

*c

and outward-pointing flux across interface

*f

may be recovered.

Member Function Documentation

template<class GridInterface, class RockInterface, class BCInterface, template< class GridIF, class RockIF > class InnerProduct>
void Opm::IncompFlowSolverHybrid< GridInterface, RockInterface, BCInterface, InnerProduct >::clear ( )
inline

Clear all topologic, geometric and rock-dependent information currently held in internal data structures. Release all memory resources in preparation of constructing a solver for a new problem. Method.

clear()

must be called prior to any other method of the class.

Referenced by Opm::IncompFlowSolverHybrid< GridInterface, ReservoirProperties, BoundaryConditions, InnerProd >::init().

template<class GridInterface, class RockInterface, class BCInterface, template< class GridIF, class RockIF > class InnerProduct>
template<class Point >
void Opm::IncompFlowSolverHybrid< GridInterface, RockInterface, BCInterface, InnerProduct >::computeInnerProducts ( const RockInterface &  r,
const Point &  grav 
)
inline

Compute static (i.e., independent of saturation) parts of the spatially varying inner products $ (v,K^{-1}w) $ for each grid cell. This is often a fairly expensive process, so in order to increase the speed of the system assembly we do not wish to perform this computation for each time step unless we have to. Moreover, if the permeability field changes but the grid connections remain static, we do not have to build a new matrix structure but need only compute new values for these static inner products.

Parameters
[in]rThe reservoir properties of each grid cell. In method
computeInnerProducts()
, we only inspect the permeability field of
r
.

Referenced by Opm::IncompFlowSolverHybrid< GridInterface, ReservoirProperties, BoundaryConditions, InnerProd >::init().

template<class GridInterface, class RockInterface, class BCInterface, template< class GridIF, class RockIF > class InnerProduct>
SolutionType Opm::IncompFlowSolverHybrid< GridInterface, RockInterface, BCInterface, InnerProduct >::getSolution ( )
inline

Recover the solution to the problem defined by the parameters to method.

solve()

. This solution is meaningless without a previous call to method

solve()

.

Returns
The current solution.
template<class GridInterface, class RockInterface, class BCInterface, template< class GridIF, class RockIF > class InnerProduct>
template<class Point >
void Opm::IncompFlowSolverHybrid< GridInterface, RockInterface, BCInterface, InnerProduct >::init ( const GridInterface &  g,
const RockInterface &  r,
const Point &  grav,
const BCInterface &  bc 
)
inline

All-in-one initialization routine. Enumerates all grid connections, allocates sufficient space, defines the structure of the global system of linear equations for the contact pressures, and computes the permeability dependent inner products for all of the grid's cells.

Parameters
[in]gThe grid.
[in]rThe reservoir properties of each grid cell.
[in]gravGravity vector. Its Euclidian two-norm value represents the strength of the gravity field (in units of m/s^2) while its direction is the direction of gravity in the current model.
[in]bcThe boundary conditions describing how the current flow problem interacts with the outside world. This is used only for the purpose of introducing additional couplings in the case of periodic boundary conditions. The specific values of the boundary conditions are not inspected in
init()
.
template<class GridInterface, class RockInterface, class BCInterface, template< class GridIF, class RockIF > class InnerProduct>
void Opm::IncompFlowSolverHybrid< GridInterface, RockInterface, BCInterface, InnerProduct >::initSystemStructure ( const GridInterface &  g,
const BCInterface &  bc 
)
inline

Compute structure of coefficient matrix in final system of linear equations for this flow problem. Allocates all storage needed by the flow solver itself.

Parameters
[in]gThe grid.
[in]bcThe boundary conditions describing how the current flow problem interacts with the outside world. This is used only for the purpose of introducing additional couplings in the case of periodic boundary conditions. The specific values of the boundary conditions are not inspected in
init()
.

Referenced by Opm::IncompFlowSolverHybrid< GridInterface, ReservoirProperties, BoundaryConditions, InnerProd >::init().

template<class GridInterface, class RockInterface, class BCInterface, template< class GridIF, class RockIF > class InnerProduct>
double Opm::IncompFlowSolverHybrid< GridInterface, RockInterface, BCInterface, InnerProduct >::postProcessFluxes ( )
inline

Postprocess the solution fluxes. This method modifies the solution object so that out-fluxes of twin faces (that is, the two faces on a cell-cell intersection) will be made antisymmetric.

Returns
The maximum modification made to the fluxes.
template<class GridInterface, class RockInterface, class BCInterface, template< class GridIF, class RockIF > class InnerProduct>
template<typename charT , class traits >
void Opm::IncompFlowSolverHybrid< GridInterface, RockInterface, BCInterface, InnerProduct >::printStats ( std::basic_ostream< charT, traits > &  os)
inline

Print statistics about the connections in the current model. This is mostly for debugging purposes and should only rarely be used in client code.

Template Parameters
charTCharacter type of output stream.
traitsCharacter traits of
charT
.
Parameters
osOutput stream into which the statistics will be printed.
template<class GridInterface, class RockInterface, class BCInterface, template< class GridIF, class RockIF > class InnerProduct>
void Opm::IncompFlowSolverHybrid< GridInterface, RockInterface, BCInterface, InnerProduct >::printSystem ( const std::string &  prefix)
inline

Output current system of linear equations to permanent storage in files. One file for the coefficient matrix and one file for the right hand side. This is mostly useful whilst debugging the solver and should rarely be used from client code.

The system is stored in a format which is suitable for importing into MATLAB or MATLAB-compatible software. In particular, using the MATLAB

LOAD

and

SPCONVERT

functions makes it easy to reconstruct the system of linear equations from within MATLAB.

Parameters
[in]prefixPrefix from which file names for the coefficient matrix and right hand side data are derived. Specifically, the matrix data is output to the file
prefix +
"-mat.dat"
while the right hand side data is output to the file
prefix + "-rhs.dat"
.
template<class GridInterface, class RockInterface, class BCInterface, template< class GridIF, class RockIF > class InnerProduct>
template<class FluidInterface >
void Opm::IncompFlowSolverHybrid< GridInterface, RockInterface, BCInterface, InnerProduct >::solve ( const FluidInterface &  r,
const std::vector< double > &  sat,
const BCInterface &  bc,
const std::vector< double > &  src,
double  residual_tolerance = 1e-8,
int  linsolver_verbosity = 1,
int  linsolver_type = 1,
bool  same_matrix = false,
int  linsolver_maxit = 0,
double  prolongate_factor = 1.6,
int  smooth_steps = 1 
)
inline

Construct and solve system of linear equations for the pressure values on each interface/contact between neighbouring grid cells. Recover cell pressure and interface fluxes. Following a call to.

solve()
@encode, you may recover the flow solution from the
@code getSolution()

method.

Template Parameters
FluidInterfaceType presenting an interface to fluid properties such as density, mobility &c. The type is expected to provide methods
phaseMobilities()
and
phaseDensities()
for phase mobility and density in a single cell, respectively.
Parameters
[in]rThe fluid properties of each grid cell. In method
solve()
we query this object for the phase mobilities (i.e.,
r.phaseMobilities()
) and the phase densities (i.e.,
phaseDensities() @encode) of each phase.
@param [in] sat
Saturation of primary phase. One scalar value for each
grid cell. This parameter currently limits @code
IncompFlowSolverHybrid
to two-phase flow problems.
[in]bcThe boundary conditions describing how the current flow problem interacts with the outside world. Method
solve()
inspects the actual values of the boundary conditions whilst forming the system of linear equations.

Specifically, the

bc.flowCond(bid)

method is expected to yield a valid

FlowBC

object for which the methods

pressure()

,

pressureDifference()

, and

outflux()

yield valid responses depending on the type of the object.

Parameters
[in]srcExplicit source terms. One scalar value for each grid cell representing the rate (in units of m^3/s) of fluid being injected into (>0) or extracted from (<0) a given grid cell.
[in]residual_toleranceControl parameter for iterative linear solver software. The iteration process is terminated when the norm of the linear system residual is less than
residual_tolerance
times the initial residual.
[in]linsolver_verbosityControl parameter for iterative linear solver software. Verbosity level 0 prints nothing, level 1 prints summary information, level 2 prints data for each iteration.
[in]linsolver_typeControl parameter for iterative linear solver software. Type 0 selects a ILU0/CG solver, type 1 selects AMG/CG, type 2 selects KAMG/CG, type 3 selects AMG/CG with fast Gauss-Seidel smoothing
[in]linsolver_maxitmaximum iterations allowed
[in]prolongate_factorFactor to scale the prolongated coarse coarse grid correction
[in]smooth_stepsNumber of smoothing steps to be used for pre and post smoothing in AMG
[in]same_matrixWhether the matrix is the same as in the previous solve.
The documentation for this class was generated from the following file: