Finite difference modeling

Prof. dr.ir.T.N.Olsthoorn

Heemstede, Sept. 2016, 24 May 2017

Approach

In this chapter we set-up a 3D steady-state finite difference model from scratch. We do this by computing a numerical groundwater problem step by step, by hand, using finite difference, building up the pieces of the model, which we will assemble in the next chapter.

Setup of the model by specifying its dimensions

The 3D steady state FDM will be based on a regular grid consisting of rows an columns and layers. The column widths and the row heigts are constant on a per column and per row basis, but the layer thickness can vary on a cell by cell basis. The grid of a full 3D model will thus be specified in general by a vector of x cell boundary coordinates, a vector y row boundary coordinates and a full 3D array of cell top and bottom coordinates.

Notice that the arrays are interpreted as [z, y, x] or [layer row col]. This is a convenience in Python where when Phi is a 3D array of the shape of the grid [Nz, Ny, Nz] we have

Phi[k].shape is [Ny, Nx], the entire layer number i. Phi[k][j] = Nx, the entire row j of layer i. Phi[k][j][i] = the head in cell [k, j, i] which is the same as Phi[k, j, i]

In [1]:
import numpy as np

# specify a rectangular grid
x = np.arange(-1000., 1000., 25.)
y = np.arange(-1000., 1000., 25.) # backward, i.e. first row grid line has highest y
z = np.arange(-100., 0., 20.)  # backward, i.e. from top to bottom

From these coordinates we obtain the number of cells along each axis and the cell sizes and

In [3]:
# as well as the number of cells along the three axes
Nx = len(x)-1
Ny = len(y)-1
Nz = len(z)-1

sz = (Nz,Ny,Nx)  # the shape of the model
Nod = np.prod(sz) # total number of cells in the model

# from this we have the width of columns, rows and layers
dx = np.diff(x).reshape(1, 1, Nx)
dy = np.diff(y).reshape(1, Ny,1)
dz = np.abs(np.diff(z)).reshape(Nz, 1,1)

IBOUND array - telling which cells are active and which have a prescribed head

Let’s first specify which of the cells have their head prescrided and which cells are inactive. We have to tackle inactive cells early to make sure their conductance is made zero (in case there conductivities might be specified as non-zeros).

We do that my means of a so-called boundary array IBOUND (MODFLOW terminology), which is an integer array of the shape of the model grid that tells which cells have a prescribed head, which cells are inactive (i.e. which cells does not take part of the computation, such as cells that represent impermeable rock) and for which cells the head should be computed.

  • IBOUND > 0, means heads will be computed
  • IBOUND == 0, means cells are inactive
  • IBOUND <0 , means heads prescribed

In this particular example we specify that the vertical zx plane at the last row of the model will have prescribed heads equal to zero.

In [7]:
IBOUND = np.ones(sz)
IBOUND[:,-1,:] = -1  # last row of model heads are prescribed
IBOUND[:, 40:45, 20:70]=0 # these cells are inactive

This boundary array makes it easy telling which cells cells are active (head computed), inactive, and fixed-head.

In [8]:
active = (IBOUND>0).reshape(Nod)  # active is now a vector of booleans of length Nod
inact  = (IBOUND==0).reshape(Nod) # dito for inact
fxhd   = (IBOUND<0).reshape(Nod)  # dito for fxhd

Cell conductancies: defining the ease of flow between adjacent cells

The first thing to define based on the properties of the cells is the flow resistance of each cell in the 3 grid directions, x, y and z. For that we need the cell sizes from the coordinates and the hydraulic conductivities in the x, y and z direction. The latter are given as full 3D arrays kx, ky, kz whose shapes correspond to that of the model mesh.

In [9]:
k = 10.0 # m/d uniform conductivity
kx = k * np.ones(sz) # [L/T] 3D kx array
ky = k * np.ones(sz) # [L/T] 3D ky array with same values as kx
kz = k * np.ones(sz) # [L/T] 3D kz array with same values as kx

The flow resistances for each cell is the head loss across opposite cell faces due to a unit flux through the cell along the axis perperndicular to them. These resistances are cell properties that can immediately be computed for the entire grid of the model. Because we always need the resistance between the cell center and its outer faces, we use the factor 0.5 (flow over half the lenght of the cell in each direction)

In [10]:
# half cell flow resistances
Rx = 0.5 * dx / (dy * dz) / kx  # [T/L2], flow resistance half cell in x-direction
Ry = 0.5 * dy / (dz * dx) / ky  # same in y-direction
Rz = 0.5 * dz / (dx * dy) / kz  # same in z-direction

Make inactive cells inactive by setting their resistance to np.Inf (infinite):

In [11]:
Rx = Rx.reshape(Nod,); Rx[inact] = np.Inf; Rx=Rx.reshape(sz)
Ry = Ry.reshape(Nod,); Ry[inact] = np.Inf; Ry=Ry.reshape(sz)
Rz = Rz.reshape(Nod,); Rz[inact] = np.Inf; Rz=Rz.reshape(sz)

From this we compute the conductance between each pair of adjacent cells across their connecting cell face. The conductance is just the reverse of the resistance of the two connected half cells. This resistance is the sum of the resistances of the two connected half cells because these resistances are placed in series with respect to the flow.

In [12]:
# conductances between adjacent cells
Cx = 1 / (Rx[:, :, :-1] + Rx[:, :,1:]) # [L2/T] in x-direction
Cy = 1 / (Ry[:, :-1,:] + Ry[:, 1:,:]) # idem in y-direction
Cz = 1 / (Rz[:-1,:,:] + Rz[1:,:,:]) # idem in z-direction

Setting up the system matrix - set of water balance equations

The system matrix has size of (Nod, Nod) allowing a connection between each pair of cells. Of course only cells that share their cell face are connected in reality. In a 3D model this means that each cell is connected to its 6 neighbors instead of to all other cells in the model. This means that most of the matrix entries will be zero.

To be able to indentify adjacent cells we generate cell numbers in an array that has the size of the model grid:

In [14]:
NOD = np.arange(Nod).reshape(sz) # this is a full 3D array of node numbers (cell numbers)

With this array it’s easy to identify adjacent cells by their cell number. Thus we generate arrays with the cel numbers of right hand neigbor of the cells (east neighbor), the left hand neighbor (the west neigbor), the north neighbor, south neighbor, the top neighbor and the bottom neighbor as follows

In [15]:
IE = NOD[:, :, 1:] # numbers of the eastern neighbors of each cell
IW = NOD[:, :, :-1] # same western neighbors
IN = NOD[:, :-1,:] # same northern neighbors
IS = NOD[:, 1:,:] # southern neighbors
IT = NOD[:-1,:,:] # top neighbors
IB = NOD[1:,:,:] # bottom neighbors

Notice that the shape of the IE and IW is the same as that of Cx, the size of IN and IS is the same as that of Cy and the size of IT and IB is the same as that of Cx.

To put the conductances into the system matrix we need their row and column indices together with their value, so that we can say a[j,i] = value. Because we have the numbers of adjacent cells in the arrays IE, IW etc, we can immediately place all the system matrix coefficiencts at the place into a sparse matrix.

In [16]:
import scipy.sparse as sp

R = lambda x : x.ravel() # define short hand for x.ravel()

# notice the call signature:
#      csc_matrix( (data, (row_index, col_index) ), (M,N)); This is a tuple within tuple.
A = sp.csc_matrix(( np.concatenate(( R(Cx), R(Cx), R(Cy), R(Cy), R(Cz), R(Cz)) ),\
                   (np.concatenate(( R(IE), R(IW), R(IN), R(IS), R(IB), R(IT)) ),\
                    np.concatenate(( R(IW), R(IE), R(IS), R(IN), R(IT), R(IB)) ),\
                  )),(Nod,Nod))

We now have to define the diagonal elements of the system matrix A, i.e. the values a[i,i] for i=[0:Nod].

These are just the negative sum of the row coefficients. Hence we sum A over the second axis (axis=1) to get them in a [Nod,1] sized vector. (Notice stat sparace matrix derived vectors keep their orientation, contrary to vectors obained from numpy arrays, which produce dimensionless vectors).

Generate the diagonal values:

In [17]:
# to use the vector of diagonal values int a call of sp.diags() we need to have it aa a
# standard nondimensional numpy vector.
# To get this:
# - first turn the matrix obtained by A.sum(axis=1) into a np.array by np.array( .. )
# - then take the whole column to loose the array orientation (to get a dimensionless numpy vector)
adiag = np.array(-A.sum(axis=1))[:,0]

Then generate a diagonal array from these values, that we can add it to A, to complete A.

In [18]:
Adiag = sp.diags(adiag)

More complex alternative: Generate diagonal array by calling the csr_matrix constructor:

Adiag = sp.csr_matrix((adiag,(np.arange(Nod),np.arange(Nod))),(Nod,Nod))

Boundary conditions

For this chapter we only use fixed flow and fixed head boundary conditions.

Fixed flows

Fixed flow boundary conditions are specified by an 3D array of the size of the grid. Each values specifies the inflow for the corresponding cell (injections are positive). Cells without a specified flow are, in fact, cells where the specified flow is zero. Hence the fixed-flows array is a full 3D array with flow values that are zero where no flow enters or leaves the cells and have non-zero values elsewhere.

For this example, we specify a single extraction of Q=-1200 m3/d in cell [30,25,2]:

In [19]:
FQ = np.zeros(sz)    # all flows zero. Note sz is the shape of the model grid
FQ[2, 30, 25] = -1200  # [m3/d] extraction in this cell

The righ-hand size of the matrix equation to be solved, the vector RHS, contains the flows. So we can generate it by assignment of FQ and converting it to a numpy vector

In [20]:
RHS = FQ.reshape(Nod)

See further down how we use RHS for only the active and non-fixed head rows.

The next step is to add fixed head boundary conditions.

Fixed heads

Fixed heads are known heads. This implies that in the set of equations that represent the model, i.e

\(A \times Phi = RHS\)

Some of the Phis are prescribed and should not be computed as defined by IBOUND and contained in the boolean vectors active, fxhd and inact specified and computed above.

Now that we know which cell have fixed heads, we can multiply out these heads with the corresponding columns of the system matrix, which yields a vector of constant values with dimension flow [m3/d] that can be added to the fixed flow vector in the RHS vector. The RHS vector is now the sum of the FQ and the contribution from the fixed heads.

Notice that the fixed heads will be obtained from the given array HI of the initial heads, where the head in the cells where IBOUND>0 correspond with the fixed heads.

In [21]:
HI = np.zeros(sz)

We reshape FQ and HI to a column vector to allow matrix multiplication

In [22]:
RHS = FQ.reshape(Nod,1) - A[:,fxhd].dot(HI.reshape(Nod,1)[fxhd])

We have now the complete RHS of the matrix equation to solve:

\(A \times Phi = RHS\)

Solving the matrix equation for the unknown heads

We use the sparse matrix solver in module scipy.sparse.linalg to compute the unknown heads.

In [23]:
from scipy.sparse.linalg import spsolve # import with from to use its short name

Of course we only need the active rows and columns of A and the active rows from RHS.

But first allocate a full-fledged vector of heads to store the result.

In [24]:
Phi = HI.flatten()

Then compute the unknown heads (i.e. the active cels only).

Remark: If we want to select a submatrix from A defiend by a given vectors of row and column indices, we can do so in sequence: Rows (I) first, columns (J) next, like so:

\(A[I][:,J]\)

which we apply in the next line

In [25]:
 Phi[active] = spsolve( (A+Adiag)[active][:,active] ,RHS[active] )

At this point we solved the problem and now have the heads for all cells in the vector Phi.

We didn’t touch the rows and columns that are inactive. So the heads of these inactive cells whatever they are in HI are know still in Phi. Just to make sure we detect them and won’t use them, set them to NaN (Not a Number).

In [26]:
Phi[inact] = np.NaN

Finally we reshape the head vector to that of the model grid.

In [27]:
Phi=Phi.reshape(sz) # reshape vector Phi to 3D shape of the grid
In [28]:
Phi # show Phi
Out[28]:
array([[[ 1.77107246,  1.77132861,  1.77183761, ...,  1.56399507,
          1.56320139,  1.56280385],
        [ 1.77081631,  1.77107575,  1.77159133, ...,  1.56360172,
          1.56280525,  1.56240631],
        [ 1.77030071,  1.77056677,  1.77109553, ...,  1.56281367,
          1.56201158,  1.56160982],
        ...,
        [ 0.05529699,  0.05524583,  0.0551436 , ...,  0.03419641,
          0.03426788,  0.03430386],
        [ 0.02763571,  0.02761016,  0.02755909, ...,  0.01708962,
          0.01712512,  0.01714299],
        [ 0.        ,  0.        ,  0.        , ...,  0.        ,
          0.        ,  0.        ]],

       [[ 1.77107246,  1.77132861,  1.77183761, ...,  1.56399507,
          1.56320139,  1.56280385],
        [ 1.77081631,  1.77107575,  1.77159133, ...,  1.56360172,
          1.56280525,  1.56240631],
        [ 1.77030071,  1.77056677,  1.77109553, ...,  1.56281367,
          1.56201158,  1.56160982],
        ...,
        [ 0.05529699,  0.05524583,  0.0551436 , ...,  0.03419641,
          0.03426788,  0.03430386],
        [ 0.02763571,  0.02761016,  0.02755909, ...,  0.01708962,
          0.01712512,  0.01714299],
        [ 0.        ,  0.        ,  0.        , ...,  0.        ,
          0.        ,  0.        ]],

       [[ 1.77107246,  1.77132861,  1.77183761, ...,  1.56399507,
          1.56320139,  1.56280385],
        [ 1.77081631,  1.77107575,  1.77159133, ...,  1.56360172,
          1.56280525,  1.56240631],
        [ 1.77030071,  1.77056677,  1.77109553, ...,  1.56281367,
          1.56201158,  1.56160982],
        ...,
        [ 0.05529699,  0.05524583,  0.0551436 , ...,  0.03419641,
          0.03426788,  0.03430386],
        [ 0.02763571,  0.02761016,  0.02755909, ...,  0.01708962,
          0.01712512,  0.01714299],
        [ 0.        ,  0.        ,  0.        , ...,  0.        ,
          0.        ,  0.        ]],

       [[ 1.77107246,  1.77132861,  1.77183761, ...,  1.56399507,
          1.56320139,  1.56280385],
        [ 1.77081631,  1.77107575,  1.77159133, ...,  1.56360172,
          1.56280525,  1.56240631],
        [ 1.77030071,  1.77056677,  1.77109553, ...,  1.56281367,
          1.56201158,  1.56160982],
        ...,
        [ 0.05529699,  0.05524583,  0.0551436 , ...,  0.03419641,
          0.03426788,  0.03430386],
        [ 0.02763571,  0.02761016,  0.02755909, ...,  0.01708962,
          0.01712512,  0.01714299],
        [ 0.        ,  0.        ,  0.        , ...,  0.        ,
          0.        ,  0.        ]]])

Plotting the heads as contours

Import the required plotting module and setup the plot.

In [29]:
%matplotlib notebook
import matplotlib.pyplot as plt # combines namespace of numpy and pyplot

For coordinates of the cells use their centers.

In [30]:
xm = 0.5 * (x[:-1] + x[1:]) # [L] coordinates of column centers
ym = 0.5 * (y[:-1] + y[1:]) # [L] coordinates of row centers
layer = 2 # contours for this layer
nc = 50   # number of contours in total

Plot the results using plt functions like in Matlab

In [31]:
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.title("Contours (%d in total) of the head in layer %d with inactive section" % (nc, layer))
plt.contour(xm, ym, Phi[layer], nc)
Out[31]:
<matplotlib.contour.QuadContourSet at 0x10e3410f0>
  • The white area is the aquifer part that was defined as inactive (impervious).
  • The red trough is the well location.
  • The two flanks and the front sides are closed (no FQ and no fixed head).
  • The head at the back side is prescribed and maintained at zero.

Conclusion

In this chapter we have developed, from scratch, a full 3D finite difference model for which fixed heads, fixed flows and inactive subareas are be prescribed. While developing we also computed a concrete example of which the head contours were finally shown.

To turn this model into a general Python function that can compute arbitrary 3D steady-state models of this kind, we only have to gather the developed lines and put them under a function definition, add some help strings (a doc string) and add error checking for convenience, while the data for this specific case, like the grid dimensions (x,y,z) the conductivities (kx,ky,kz), the prescribed flows FQ and initial heads HI plus the boundary array IBOUND that tells which cells have fixed heads and which are inactive, are to be passed as in user-given arguments at the function call like so:

Phi = fdm3(x,y,z,kx,ky,kz,FQ,IH,IBOUND) # function that solves arbitrary 3D steady steate finite difference model

We do this in the next section.