Finite Difference Discretization of Darcy's law and solving with Picard method

Question

I am trying to discretize Darcy's Law using finite differences and then solving the resulting linear system of equations with the Picard method. So far only in 1D and the steady-state (no time derivatives).

Darcy's law reads: $$C(h)\frac{\partial h}{\partial t} = \nabla\cdot\left(K(h)\nabla h\right) + \frac{\partial K(h)}{\partial z}$$ the left hand side is zero, since I am (so far) only interested in the steady state.

Using finite differences with $K$ and $h$ values located at the cell centers, I get: $$0=\frac{K_{i+1/2}\left(\frac{h_{i+1}-h_{i}}{dz}\right)-K_{i-1/2}\left(\frac{h_{i}-h_{i-1}}{dz}\right)}{dz} + \frac{K_{i+1}-K_{i-1}}{2dz}$$ now I replace $K_{i+1/2}$ and $K_{i-1/2}$ by interpolation terms: $$0=\frac{\frac{K_{i}+K_{i+1}}{2}\left(\frac{h_{i+1}-h_{i}}{dz}\right)-\frac{K_{i}+K_{i-1}}{2}\left(\frac{h_{i}-h_{i-1}}{dz}\right)}{dz} + \frac{K_{i+1}-K_{i-1}}{2dz}$$ and simplify it to: $$0=\frac{1}{2dz^2}\left[h_{i-1}\left(K_{i-1}+K_i\right)-h_i\left(K_{i-1}+K_i+K_{i+1}\right)+h_{i+1}\left(K_i+K_{i+1}\right)+K_{i+1}dz-K_{i-1}dz\right]$$

I believe that until here things are okay, correct me if I'm wrong.

Now I transform the above equation into a linear system of equations. At the top (between $h_N$ and $h_{N+1}$) there are Dirichlet boundary conditions with $h_{N+1/2}=h_B$ and at the bottom (between $h_0$ and $h_1$) there are neumann boundary conditions $\nabla h_{1/2}=0$: $$\frac{1}{2dz^2} \begin{bmatrix} -1 & 1 & 0 & 0 & \dots & 0 & 0 & 0 & 0 \\ c_{i-1} & c_i & c_{i+1} & 0 & \dots & 0 & 0 & 0 & 0 \\ 0 & c_{i-1} & c_i & c_{i+1} & \ddots & 0 & 0 & 0 & 0 \\ 0 & 0 & c_{i-1} & c_i & \ddots & 0 & 0 & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \ddots & c_i & c_{i+1} & 0 & 0 \\ 0 & 0 & 0 & 0 & \ddots & c_{i-1} & c_i & c_{i+1} & 0 \\ 0 & 0 & 0 & 0 & \dots & 0 & c_{i-1} & c_i & c_{i+1} \\ 0 & 0 & 0 & 0 & \dots & 0 & 0 & 1 & 1 \\ \end{bmatrix} \begin{bmatrix} h_0\\ h_1\\ h_2\\ h_3\\ \vdots\\ h_{N-2}\\ h_{N-1}\\ h_{N}\\ h_{N+1} \end{bmatrix} +\frac{1}{2dz} \begin{bmatrix} 0\\ K_2-K_0\\ K_3-K_1\\ K_4-K_2\\ \vdots\\ K_{N-2}-K_{N-4}\\ K_{N-1}-K_{N-3}\\ K_{N}-K_{N-2}\\ 0 \end{bmatrix} = \begin{bmatrix} 0\\ 0\\ 0\\ 0\\ \vdots\\ 0\\ 0\\ 0\\ h_Bdz^2 \end{bmatrix} $$ with $c_i = -K_{i-1}-K_i-K_{i+1}$, $c_{i+1}=K_i+K_{i+1}$ and $c_{i-1} = K_{i-1}+K_i$. $h_N$ referes to top most cell and $h_0$ to bottom most cell. The boundary of the domain is located between $h_0$ and $h_1$ and $h_N$ and $h_{N+1}$.

This is there I'm not so sure anymore. Did I get the matrix right? And then how to get from here to solving it using the Picard method, which as far as I understand is only good for first-order ODEs and not second-order as Darcy's Eq. is.

Answer 1

So here is what I came up with:

The boundary conditions were wrong: the bottom boundary was supposed to read $h_{N}=h_{N-1}$ and the top one had to be $h_{0}=h$.

This in turn lead to a matrix which differed in the top and bottom rows: $$\frac{1}{2dz^2} \begin{bmatrix} 1 & 0 & 0 & 0 & \dots & 0 & 0 & 0 & 0 \\ c_{i-1} & c_i & c_{i+1} & 0 & \dots & 0 & 0 & 0 & 0 \\ 0 & c_{i-1} & c_i & c_{i+1} & \ddots & 0 & 0 & 0 & 0 \\ 0 & 0 & c_{i-1} & c_i & \ddots & 0 & 0 & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \ddots & c_i & c_{i+1} & 0 & 0 \\ 0 & 0 & 0 & 0 & \ddots & c_{i-1} & c_i & c_{i+1} & 0 \\ 0 & 0 & 0 & 0 & \dots & 0 & c_{i-1} & c_i & c_{i+1} \\ 0 & 0 & 0 & 0 & \dots & 0 & 0 & -1/2 & 1 \\ \end{bmatrix} \begin{bmatrix} h_0\\ h_1\\ h_2\\ h_3\\ \vdots\\ h_{N-2}\\ h_{N-1}\\ h_{N}\\ h_{N+1} \end{bmatrix} +\frac{1}{2dz} \begin{bmatrix} 0\\ K_2-K_0\\ K_3-K_1\\ K_4-K_2\\ \vdots\\ K_{N-2}-K_{N-4}\\ K_{N-1}-K_{N-3}\\ K_{N}-K_{N-2}\\ 1 \end{bmatrix} = \begin{bmatrix} h\\ 0\\ 0\\ 0\\ \vdots\\ 0\\ 0\\ 0\\ 0 \end{bmatrix} $$

The Picard method uses the above matrix and right hand side (rhs) to iteratively solve the ODE by updating the matrix and rhs after each iteration using the previously gained approximation of h.

Here is a Python implementation of the Picard method:

  dz = 1
  h = array([-20.7]+([-61.5] * (40/dz)))
  r = []
  A, b = build_soe(h, dz)
  r.append(sum(abs(A.dot(h)-b)))
  for i in range(45):
      A, b = build_soe(h, dz)
      h_ = solve(A,b)
      r.append(sum(abs(abs(h_)-abs(h))))
      h = h_

build_soe generates the matrix A and right-hand-side b as described above. The code reaches machine precision after less than 45 iterations.

A more detailed description can be found here.