Solve mixed formulation of incompressible flow modelled by Darcy's law
] The solver is based on a hybrid discretization scheme which yields a system of linear equations of the form
=
] where represents the interface fluxes for each interface for each cell, are the cell pressures and are the interface pressures.
More...
#include <IncompFlowSolverHybrid.hpp>
|
typedef const FlowSolution & | SolutionType |
| Type representing the solution to the problem defined by the parameters to. More...
|
|
|
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 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...
|
|
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
] The solver is based on a hybrid discretization scheme which yields a system of linear equations of the form
=
] where represents the interface fluxes for each interface for each cell, are the cell pressures and are the interface pressures.
Through a Schur complement analysis, the above system of linear equations is reduced to an equivalent system of the form
=
] where , , and . Similarly, and .
The system 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
] Specifically, the matrix is diagonal (a single scalar per grid cell) and the matrix is available as a bi-product of the Schur reduction process. Finally, in the case of mimetic discretizations using the class, the matrix is directly available through an explicit formula.
The matrix 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 and depend only on the grid topolgy. Consequently, the cost of assembling the system 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 and system right hand side on a cell-by-cell basis. This yields a number of important simplifications. In particular, the matrix 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 with the number of elements equal to the number of interfaces (i.e., neighbours) of the cell. Finally, the matrix is computed as
where denotes the static part of and 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 , and 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 –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 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 , 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 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 , computes the static (i.e., geology and geomtry-dependent) inner product matrices . Method is offered separately in order to support solve several different property models on the same underlying geometry (grid). Finally, method 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 method. At any time following a call to may the solution, represented by the type , to the most recently resolved flow problem be retrieved through the method. We note that is always a reference-to-const.
- Template Parameters
-
GridInterface | Type presenting an interface to a grid (typically a discretized geological model). The type is assumed to expose a forward iterator type,, and a pair of delimiters and to traverse the cells of a grid. The cell iterator, in turn, is expected to expose a forward iterator type, and a pair of delimiters and to traverse the faces of a cell. |
RockInterface | Type presenting an interface to reservoir properties such as permeability and porosity. The type is assumed to expose a method, through which the assigned permeability of a single grid cell, represented by a, may be recovered. |
BCInterface | Type presenting an interface to boundary conditions. The type is expected to provide a method,, from which a 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. |
InnerProduct | Type presenting a specific inner product defining a discretization of the Darcy equation. |
template<class GridInterface, class RockInterface, class BCInterface, template< class GridIF, class RockIF > class InnerProduct>
Type representing the solution to the problem defined by the parameters to.
. Always a reference-to-const. The exposes methods and from which the cell pressure in cell and outward-pointing flux across interface may be recovered.
template<class GridInterface, class RockInterface, class BCInterface, template< class GridIF, class RockIF > class InnerProduct>
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 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] | r | The reservoir properties of each grid cell. In method , we only inspect the permeability field of. |
Referenced by Opm::IncompFlowSolverHybrid< GridInterface, ReservoirProperties, BoundaryConditions, InnerProd >::init().
template<class GridInterface, class RockInterface, class BCInterface, template< class GridIF, class RockIF > class InnerProduct>
Recover the solution to the problem defined by the parameters to method.
. This solution is meaningless without a previous call to method .
- 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] | g | The grid. |
[in] | r | The reservoir properties of each grid cell. |
[in] | grav | Gravity 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] | bc | The 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. |
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] | g | The grid. |
[in] | bc | The 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. |
Referenced by Opm::IncompFlowSolverHybrid< GridInterface, ReservoirProperties, BoundaryConditions, InnerProd >::init().
template<class GridInterface, class RockInterface, class BCInterface, template< class GridIF, class RockIF > class InnerProduct>
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
-
charT | Character type of output stream. |
traits | Character traits of. |
- Parameters
-
os | Output stream into which the statistics will be printed. |
template<class GridInterface, class RockInterface, class BCInterface, template< class GridIF, class RockIF > class InnerProduct>
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 and functions makes it easy to reconstruct the system of linear equations from within MATLAB.
- Parameters
-
[in] | prefix | Prefix from which file names for the coefficient matrix and right hand side data are derived. Specifically, the matrix data is output to the file while the right hand side data is output to the file. |
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.
@encode, you may recover the flow solution from the
method.
- Template Parameters
-
FluidInterface | Type presenting an interface to fluid properties such as density, mobility &c. The type is expected to provide methods and for phase mobility and density in a single cell, respectively. |
- Parameters
-
[in] | r | The fluid properties of each grid cell. In method we query this object for the phase mobilities (i.e.,) 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] | bc | The boundary conditions describing how the current flow problem interacts with the outside world. Method inspects the actual values of the boundary conditions whilst forming the system of linear equations. |
Specifically, the method is expected to yield a valid object for which the methods , , and yield valid responses depending on the type of the object.
- Parameters
-
[in] | src | Explicit 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_tolerance | Control parameter for iterative linear solver software. The iteration process is terminated when the norm of the linear system residual is less than times the initial residual. |
[in] | linsolver_verbosity | Control 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_type | Control 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_maxit | maximum iterations allowed |
[in] | prolongate_factor | Factor to scale the prolongated coarse coarse grid correction |
[in] | smooth_steps | Number of smoothing steps to be used for pre and post smoothing in AMG |
[in] | same_matrix | Whether the matrix is the same as in the previous solve. |
The documentation for this class was generated from the following file:
|