I am about to generalise the geothermal subsurface flow simulator I'm
working on so that it can simulate dual-porosity systems, and need some
advice on the best way to to this within the PETSc framework.
1) Dual porosity:
For anyone not familiar with the idea, dual porosity models are a way of
simulating flow in fractured media. The main flow in the system takes
place along the fractures, but fluid is also stored in the 'matrix' rock
between the fractures. The rocks in the fracture and matrix cells
usually have different properties (e.g. permeability, porosity).
Starting from a conventional single-porosity finite volume model, dual
porosity effectively adds a one-dimensional sub-model inside each cell
(or possibly only in some cells), representing storage in the 'matrix'
rock. The original single-porosity cells now represent the fractures in
the rock. Their effective volumes are reduced but their connectivity to
other fracture cells is unchanged. Each one gets an additional
connection to a matrix cell. However, the matrix cells are usually not
connected to each other. In the simplest case there is just one matrix
cell for every fracture cell, though it is often desirable to use two or
more, to get better representation of the flow in and out of the matrix.
In that case the higher-level matrix cells form a one-dimensional
sub-model for each cell.
2) Simple solution method:
The simplest way to implement this approach is just to add in all the
matrix cells and treat them the same way as the fracture cells. If the
original number of single-porosity cells was n, and the number of matrix
cells per fracture cell is m, then the Jacobian for the system is of
size (m+1) * n. (In practice we are almost always solving for more than
one variable per cell, e.g. pressure and temperature, so if there are p
variables per cell, the total Jacobian size is really (m+1) * n * p, but
if we think of the Jacobian as a block matrix with blocksize p then the
total number of blocks is still (m+1) * n.)
2) More efficient solution method:
The straight-forward approach is not very efficient, because it doesn't
take any advantage of the particular sparsity pattern of the Jacobian.
Because all the matrix rock sub-models are one-dimensional, the Jacobian
has a block-tridiagonal high-level structure, in which all of the the
block matrices (of size n) are themselves block-diagonal (with blocksize
p), except for the upper left one which has the same sparsity pattern as
the Jacobian for the original single-porosity system.
An efficient way to solve this linear system is described by, among
others, Zyvoloski et al (2008):
http://www.sciencedirect.com/science/article/pii/S0309170807001741
The method is detailed on p. 537, equations 1 - 6. What it amounts to is
taking successive 2-block-by-2-block sub-systems, starting from the
bottom right, and replacing the upper left sub-matrix in each case by
its Schur complement (with a similar modification to the right hand side
vector). Because all the sub-matrices are block diagonal, these can be
computed quite cheaply (and in parallel, if the sub-matrices are
distributed appropriately across the processors). When you get back up
to the top left of the whole matrix, the remaining modified sub-matrix
of size n can be solved on its own to get the solution in the fracture
cells. You can then back-substitute to get the solution in the matrix
rock cells.
In other words, after a relatively inexpensive reduction process, you
wind up solving a linear system of size n (same as the original
single-porosity system, and with the same sparsity pattern- only the
diagonal and RHS are modified) instead of (m+1) * n. This whole process
is usually a lot faster than solving the whole (m+1) * n sized linear
system.
3) PETSc implementation
So, how would I implement this approach in the PETSc framework? Because
the flow equations are non-linear, solving these linear systems happens
multiple times each time-step as part of a SNES solve.
One way might be to form the whole Jacobian but somehow use a modified
KSP solve which would implement the reduction process, do a KSP solve on
the reduced system of size n, and finally back-substitute to find the
unknowns in the matrix rock cells.
Another way might be to form only the reduced-size Jacobian and the
other block-diagonal matrices separately, use KSP to solve the reduced
system but first incorporate the reduction process into the Jacobian
calculation routine, and somewhere a post-solve step to back-substitute
for the unknowns in the matrix cells. However currently we are using
finite differences to compute these Jacobians and it seems to me it
would be messy to try to do that separately for each of the
sub-matrices. Doing it the first way above would avoid all that.
Any suggestions for what might be a good approach? or any other ideas
that could be easier to implement with PETSc but have similar
efficiency? I didn't see anything currently in PETSc specifically for
solving block-tridiagonal systems.
Thanks, Adrian
--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611
