[petsc-users] Jacobian matrix for dual porosity model

[Adrian Croucher]

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.croucher at auckland.ac.nz
tel: +64 (0)9 923 4611