Darcy flow in porous media

Solved using the (Preconditioned) Conjugate-Gradient Method

AMath 585, Winter Quarter 2020 at the University of Washington. Developed by R.J. LeVeque and distributed under the BSD license. You are free to modify and use as you please, with attribution.

These notebooks are all available on Github.

This notebook illustrates a sample application of a two-dimensional variable coefficient elliptic equation modeling the flow of a fluid through a porous medium. Realistic examples would solve the analogous system in three space dimensions.

This notebook also contains an implementation of the Conjugate Gradient algorithm to solve the system, taken mostly from ConjugateGradient2D_vcoeff.html, where it is discussed in more detail.

The Preconditioned Conjugate Gradient (PCG) algorithm is also implemented, based on diagonal preconditioner and using the corrected algorithm from PCG.html.

In [1]:
%matplotlib inline
In [2]:
from pylab import *

Animating the iteration

To use widgets, set use_widgets = True. For javascript animation, set it to False.

In [3]:
use_widgets = False

if use_widgets:
    from ipywidgets import interact
    import ipywidgets as widgets

    def animate_figs(figs):
        show_frame = lambda frameno: display(figs[frameno])
        interact(show_frame, frameno=widgets.IntSlider(min=0,max=len(figs)-1, value=0))
else:
    from jsanimate_figs import animate_figs

Flow in porous media

Consider creeping flow in a porous medium (for example groundwater or oil flowing through porous rock). In simple cases such flow is governed by Darcy's law, which says that the velocity vector is proportional to the gradient of the pressure, or in terms of the velocity components $(u(x,y), v(x,y))$ (in two dimensions),

\begin{equation*} \begin{split} u(x,y) &= -\gamma(x,y) p_x(x,y),\\ v(x,y) &= -\gamma(x,y) p_y(x,y), \end{split} \end{equation*}

where $\gamma(x,y) > 0$ is a measure of the permeability of the material at $(x,y)$, the material is more permeable and fluid flows more easily where $\gamma$ is large.

We also know that if the fluid is incompressible (as are water or oil, for practical purposes), then the divergence of the velocity must be zero everywhere,

$$ u_x(x,y) + v_y(x,y) = 0. $$

Pressure formulation

Inserting the above representation of $u, v$ into the divergence-free condition, this gives a variable-coefficient elliptic equation for the pressure:

$$ -(\gamma(x,y)p_x(x,y))_x - (\gamma(x,y)p_y(x,y))_y = 0. $$

If there were a source (or sink) of fluid distributed over the domain, then the right hand side would be equal to this source $f(x,y)$, but we'll consider the case with no source.

We could solve this elliptic equation, but we'd need to know boundary conditions for the pressure, which may not be easily available.

Stream function formulation

Instead we'll reformulate it in terms of a stream function $\psi(x,y)$ that is related to $u$ and $v$ by

\begin{equation*} \begin{split} u(x,y) &= \psi_y(x,y),\\ v(x,y) &= -\psi_x(x,y). \end{split} \end{equation*}

Note that $u_x + v_y=0$ for any function $\psi$ so the velocity field is automatically divergence free in this formulation!

To get an equation for $\psi(x,y)$, note that $p_x = (1/\gamma)\psi_y$ and $p_y = -(1/\gamma)\psi_x$ and then since $p_{xy} - p_{yx} = 0$, we obtain

$$ (\kappa(x,y)\psi_x(x,y))_x + (\kappa(x,y) \psi_y(x,y))_y = 0, $$

where $\kappa(x,y) = 1/\gamma(x,y) > 0$. This is our elliptic equation.

Note that the gradient of $\psi$ is $(-v,u)$, which is orthogonal to the velocity vector at each point. Hence contours of $\psi(x,y)$ are streamlines of the flow. This is another reason it's nice to compute $\psi$ rather than $p$: a contour plot of $\psi$ gives a visualization of the flow.

Also note that differencing $\psi$ between any two points in the plane gives the net flux across any curve connecting these two points:

\begin{equation*} \psi(x_2,y_2) - \psi(x_1,y_1) = \int_\Gamma \psi\cdot\tau \,ds = \int_\Gamma u\cdot n\,ds \end{equation*}

where $\Gamma$ is any curve connecting the two points, $\tau$ is the tangent to the curve, and $n$ is the normal.

In particular, this means that if we plot contour lines that are equally spaced in $\psi$, then the distance between two particular lines is equal to the total flux through the tube between the lines. Hence the flow must be faster where the lines are close together and slower where the lines are farther apart.

Setting $\psi=$ constant along a boundary forces that boundary to be a streamline of the flow.

We will use a rectangular domain in two space dimensions, $0\leq x \leq 2$ and $0\leq y \leq 1$, and set $\psi(x,0) = 0$ along the bottom boundary and $\psi(x,1) = 1$ along the top boundary so these are streamlines, i.e. solid walls that the flow cannot go through $v(x,0) = v(x,1) = 0$ all along these boundaries.

At the left and right boundaries we will impose $\psi(0,y) = y$ and $\psi(1,y) = y$, which corresponds to flow that has $u = \psi_y = 1$ everywhere at both left and right boundaries, so we have specified the horizontal component of the velocity at each side. We can think of this as imposing flow going through the domain from left to right, for example, with prescribed inflow/outflow velocities at the ends.

Note: Much of the code below is copied from earlier notebooks, so the solution of the elliptic equation is called $u$ instead of $\psi$, but don't confuse this with the horizontal velocity!

Define the domain and grid

In [4]:
ax = 0.
bx = 2.
mx = 199
dx = (bx-ax)/(mx+1.)
x_full = linspace(0,2,mx+2)

ay = 0.
by = 1.
my = 99
dy = (by-ay)/(my+1.)
y_full = linspace(0,1,my+2)

X_full, Y_full = meshgrid(x_full, y_full, indexing='ij')
X = X_full[1:-1, 1:-1]
Y = Y_full[1:-1, 1:-1]

Define gamma and kappa

For illustration we will use a piecewise constant permeability function $\gamma(x,y)$ that has a small value inside a rectangle and larger value outside. Then we define $\kappa(x,y) = 1/\gamma(x,y)$:

In [5]:
in_rect = lambda X,Y: logical_and(logical_and(X>0.9, X<1.1), \
                                  logical_and(Y>0.3, Y<0.9))

gamma = lambda X,Y: where(in_rect(X,Y), 0.1, 1.)
kappa = lambda X,Y: 1/gamma(X,Y)

#kappa = lambda X,Y: ones(X.shape)  # constant kappa for testing

kappa_full = kappa(X_full, Y_full)

f_fcn = lambda X,Y: zeros(X.shape)  # no source or sink

Define the boundary conditions

In this case we don't know the true solution, but we want to impose $\psi(x,y) = y$ along the boundaries. So we set Ubc_full to be $y$ everywhere in the domain and will only use these values along the boundary later.

In [6]:
ubc_fcn = lambda X,Y: Y
Ubc_full = ubc_fcn(X_full, Y_full)

A function to plot the approximate solution

In [7]:
def make_plot(U,k,rnorm):
    
    # colormap for kappa:
    kappa_min = kappa_full.min()
    kappa_max = 2*kappa_full.max()
    levels = linspace(kappa_min,kappa_max,21)
    cmap = get_cmap('Purples')  # white to purple
    vmin = kappa_min; vmax = kappa_max
    
    fig = figure(figsize=(10,6))
    U_full = Ubc_full.copy()  # for boundary values
    U_full[1:-1, 1:-1] = U
    
    contourf(X_full,Y_full,kappa_full,levels,
             cmap=cmap,vmin=vmin,vmax=vmax)
    
    axis('scaled')
    colorbar(shrink=0.8, label='kappa')
    
    # plot u, which corresponds to the stream function psi for this problem!
    # so contours of u are approximate streamlines of the flow
    
    Ulevels = linspace(0,1,21)
    contour(X_full,Y_full,U_full,Ulevels,
             colors='k')
    
    title('iteration %3i:   2-norm(r) = %.2e\n' % (k,rnorm) \
          + 'Streamlines of approximate solution')
    
    close(fig)
    return fig
In [8]:
make_plot(Y, 0, nan)
Out[8]:

Define kappa function and evaluate at midpoints

The difference stencil is chosen so that the resulting matrix $A$ is symmetric, as discussed in section 2.15 of the text in 1D. This can be extended to 2D by using it separately for the $(\kappa u_x)_x$ and $(\kappa u_y)_y$ terms, and so the $(i,j)$ element of the grid function representing $Au$ is:

\begin{align*} \frac{1}{\Delta x^2} \left(\kappa_{i+1/2,j}(U_{i+1,j} - U_{i,j}) - \kappa_{i-1/2,j}(U_{ij} - U_{i-1,j})\right) \\ + \frac{1}{\Delta y^2} \left(\kappa_{i,j+1/2}(U_{i,j+1} - U_{i,j}) - \kappa_{i,j-1/2}(U_{ij} - U_{i,j-1})\right) \end{align*}

To clarify the implementation we define 4 different arrays below with, for example, kappa_imh_j holding values $\kappa_{i-1/2,j} = \kappa(x_i - \Delta x/2, y_j)$ and kappa_iph_j holding values $\kappa_{i+1/2,j} = \kappa(x_i + \Delta x/2, y_j)$. These two arrays hold mostly the same values, except shifted and with values included near one $x$ boundary and not the other. So for a large problem when trying to minimize storage and computation time, a single slightly larger array could hold all these values. But here we opt for clarity.

In [9]:
X_imh_j = X - dx/2.
Y_imh_j = Y
kappa_imh_j = kappa(X_imh_j, Y_imh_j)

X_iph_j = X + dx/2.
Y_iph_j = Y
kappa_iph_j = kappa(X_iph_j, Y_iph_j)

X_i_jmh = X
Y_i_jmh = Y - dy/2.
kappa_i_jmh = kappa(X_i_jmh, Y_i_jmh)

X_i_jph = X
Y_i_jph = Y + dy/2.
kappa_i_jph = kappa(X_i_jph, Y_i_jph)
In [10]:
def matvec(V):
    """
    Given a grid function v on an mx by my grid. 
    Return b = A*v
    """
   
    # pad V with zeros around border, needed for computing centered differences:
    V_full = zeros((mx+2, my+2))
    V_full[1:-1, 1:-1] = V
    
    # set array B of shape (mx,my) at interior points (same shape as V) 
    # Note that slicing with [1:-1] corresponds to interior points,
    # [:-2] correspond to points to the left (or below) and 
    # [2:] are points to the right (or above)
    
    B =   (kappa_iph_j*(V_full[2:, 1:-1] - V_full[1:-1, 1:-1]) \
         - kappa_imh_j*(V_full[1:-1, 1:-1] - V_full[:-2, 1:-1])) / dx**2 \
        + (kappa_i_jph*(V_full[1:-1, 2:] - V_full[1:-1, 1:-1]) \
         - kappa_i_jmh*(V_full[1:-1, 1:-1] - V_full[1:-1, :-2])) / dy**2
                    
    return B
In [11]:
def inner_product(V,W):
    """
    Given two grid functions V and W on an mx by my grid. 
    Return their inner product when viewed as vectors of unknowns
    """
    return sum(V*W)

The Conjugate-Gradient algorithm

The next cell implements the C-G algorithm. This is identical to the corresponding cell in ConjugateGradient2D_vcoeff.html.

In [12]:
maxiter = 1000
kplot = 100
verbose = False
tol = 1e-14  # stop if the residual falls below tol

print('Solving on %i by %i grid with %i unknowns' \
      % (mx,my,mx*my))
print('    using tol = %.2e with maxiter = %i' % (tol,maxiter))

# right-hand side:
F_full = f_fcn(X_full, Y_full)
F = F_full[1:-1, 1:-1]  # at interior points

# adjust for Dirichlet BCs:
F[0,:] = F[0,:] - kappa_imh_j[0,:]*Ubc_full[0, 1:-1] / dx**2
F[-1,:] = F[-1,:] - kappa_iph_j[-1,:]*Ubc_full[-1, 1:-1] / dx**2
F[:,0] = F[:,0] - kappa_i_jmh[:,0]*Ubc_full[1:-1, 0] / dy**2
F[:,-1] = F[:,-1] - kappa_i_jph[:,-1]*Ubc_full[1:-1, -1] / dy**2


# initial guess:
U0_full = Y_full
U0 = U0_full[1:-1, 1:-1]  # interior points

U = U0.copy() # current iterate
r = F - matvec(U)  # initial residual
p = r.copy()  # initial direction

k = 0
rTr_k = inner_product(r,r)
rnorm = sqrt(rTr_k)
rnorms = [rnorm]  # build up list of residuals

if verbose:
    print('iteration %3i:   2-norm(r) = %.2e,   max-norm(E) = %.2e' \
          % (k,rnorm))
    
figs = []  # for the list of figures we generate

fig = make_plot(U,k,rnorm) # plot initial guess
figs.append(fig)

rTr_km = rTr_k  # in general will hold r^T * r at iteration k-1

for k in range(1,maxiter+1):
    w = matvec(p)   # the only matrix-vector multiply
    a = rTr_km / inner_product(p,w)  # alpha_{k-1} in CG algorithm
    U = U + a*p
    r = r - a*w

    rTr_k = inner_product(r,r)
    rnorm = sqrt(rTr_k)
    rnorms.append(rnorm)

    if mod(k,kplot)==0 or k==maxiter:
        # every kplot iterations create a plot:
        fig = make_plot(U,k,rnorm)
        figs.append(fig)
        
    if verbose:
        print('iteration %3i:   2-norm(r) = %.2e' \
              % (k,rnorm))
    
    # check for convergence:
    if rnorm < tol:
        print('Satisfied rnorm < tol after %i iterations' % k)
        break
        
    # determine next search direction:
    b = rTr_k / rTr_km   # beta_{k-1} in CG algorithm
    rTr_km = rTr_k       # for next iteration
    p = r + b*p           # next search direction

print('Stopped after %i iterations with 2-norm(r) = %.2e' \
   % (k,rnorm))

Solving on 199 by 99 grid with 19701 unknowns
    using tol = 1.00e-14 with maxiter = 1000
Stopped after 1000 iterations with 2-norm(r) = 1.87e-03
In [13]:
animate_figs(figs)
Out[13]:
In [14]:
semilogy(range(1,len(rnorms)+1), rnorms, 'b-')
grid(True)

Save these values for comparison below:

In [15]:
rnorms_noprecond = rnorms.copy()

Preconditioned conjugate gradient

Next we modify this to implement the preconditioned conjugate gradient method on p. 95 of the text, but with some corrections as discussed in the notebook PCG.html.

In particular, the lines defining $\alpha_{k-1}$ and $\beta_k-1$ are incorrect and should read:

$$ \alpha_{k-1} = (z_{k-1}^T r_{k-1}) / (p_{k-1}^T w_{k-1}) \quad\text{($z$ instead of $r$ in the numerator)} $$

and $$ \beta_{k-1} = (z_k^T r_k) / (z_{k-1}^T r_{k-1}) \quad\text{($z$ instead of $r$ in two places)} $$

Diagonal preconditioning

Here we only consider the simplest (but often effective) preconditioner in which $M$ is the diagonal part of $A$. For the variable coefficient problem this can help a lot, particularly if the values of $\kappa$ vary greatly.

See PCG.html for discussion of the fact that for this problem, where we have set it up with a matrix $A$ that is symmetric negative definite, we want to define $M$ based on the absolution values of the diagonal elements.

For this 2D problem the diagonal element of $A$ corresponding to the central coefficient of the 5-point stencil at each grid point $(i,j)$, with absolute value

$$ (\kappa_{i-1/2,j} + \kappa_{i+1/2,j})/\Delta x^2 + (\kappa_{i,j-1/2} + \kappa_{i,j+1/2})/\Delta y^2. $$

We can define a grid function defined on the 2D grid that has these values:

In [16]:
Mdiag = (kappa_imh_j + kappa_iph_j) / dx**2 \
        + (kappa_i_jmh + kappa_i_jph) / dy**2

Then if $r$ is a grid function on the same grid representing the residual at some point, we can solve the linear system $Mz = r$ by simply computing z = r/Mdiag in Python (componentwise division of the values at each grid point).

In [17]:
maxiter = 1000
kplot = 100
verbose = False
tol = 1e-14  # stop if the residual falls below tol

print('Solving on %i by %i grid with %i unknowns' \
      % (mx,my,mx*my))
print('    using tol = %.2e with maxiter = %i' % (tol,maxiter))

# right-hand side:
F_full = f_fcn(X_full, Y_full)
F = F_full[1:-1, 1:-1]  # at interior points

# adjust for Dirichlet BCs:
F[0,:] = F[0,:] - kappa_imh_j[0,:]*Ubc_full[0, 1:-1] / dx**2
F[-1,:] = F[-1,:] - kappa_iph_j[-1,:]*Ubc_full[-1, 1:-1] / dx**2
F[:,0] = F[:,0] - kappa_i_jmh[:,0]*Ubc_full[1:-1, 0] / dy**2
F[:,-1] = F[:,-1] - kappa_i_jph[:,-1]*Ubc_full[1:-1, -1] / dy**2


# initial guess:
U0_full = Y_full
U0 = U0_full[1:-1, 1:-1]  # interior points

U = U0.copy() # current iterate
r = F - matvec(U)  # initial residual

z = r / Mdiag  # apply preconditioner

k = 0
zTr_k = inner_product(z,r)  # modified for PCG
p = z.copy()

rTr_k = inner_product(r,r)  # for checking convergence
rnorm = sqrt(rTr_k)
rnorms = [rnorm]  # build up list of residuals

if verbose:
    print('iteration %3i:   2-norm(r) = %.2e,   max-norm(E) = %.2e' \
          % (k,rnorm))
    
figs = []  # for the list of figures we generate

fig = make_plot(U,k,rnorm) # plot initial guess
figs.append(fig)

zTr_km = zTr_k  # in general will hold z^T * r at iteration k-1

for k in range(1,maxiter+1):
    w = matvec(p)   # the only matrix-vector multiply
    a = zTr_km / inner_product(p,w)  # alpha_{k-1} in PCG algorithm
    U = U + a*p
    r = r - a*w
    z = r / Mdiag  # apply preconditioner

    zTr_k = inner_product(z,r)  # modified for PCG
    
    rTr_k = inner_product(r,r)  # for checking convergence
    rnorm = sqrt(rTr_k)
    rnorms.append(rnorm)

    if mod(k,kplot)==0 or k==maxiter:
        # every kplot iterations create a plot:
        fig = make_plot(U,k,rnorm)
        figs.append(fig)
        
    if verbose:
        print('iteration %3i:   2-norm(r) = %.2e' \
              % (k,rnorm))
    
    # check for convergence:
    if rnorm < tol:
        print('Satisfied rnorm < tol after %i iterations' % k)
        break
        
    # determine next search direction:
    b = zTr_k / zTr_km   # beta_{k-1} in PCG algorithm
    zTr_km = zTr_k       # for next iteration
    p = z + b*p           # next search direction

print('Stopped after %i iterations with 2-norm(r) = %.2e' \
   % (k,rnorm))

Solving on 199 by 99 grid with 19701 unknowns
    using tol = 1.00e-14 with maxiter = 1000
Satisfied rnorm < tol after 709 iterations
Stopped after 709 iterations with 2-norm(r) = 8.86e-15
In [18]:
animate_figs(figs)
Out[18]:
In [19]:
semilogy(range(1,len(rnorms_noprecond)+1), rnorms_noprecond, 'b-', 
         label='no preconditioner')
semilogy(range(1,len(rnorms)+1), rnorms, 'r-', 
         label='with preconditioner')
grid(True)
legend()
Out[19]:
<matplotlib.legend.Legend at 0x171c5d240>

You might want to try changing gamma so that the value in the rectangle is $0.01$ rather than $0.1$, in which case this region is less permeable and the flow has greater tendency to go around it. In this case the condition number is even larger and PCG makes an even more dramatic improvement in convergence.