Solution Method

In practical applications it is an advantage to calculate the pressure $ p$ as a correction of a 'static' pressure $ p^{ref}$ which is the solution of

$\displaystyle -(\kappa\hackscore{ki}\kappa\hackscore{kj} p\hackscore{,j}^{ref})...
...} = - (\kappa\hackscore{ki} (g\hackscore{k}- u^{N}\hackscore{k}))\hackscore{,i}$    with $\displaystyle p^{ref} = p^{D}$    on $\displaystyle \Gamma\hackscore{D}$ (110)

With setting $ u \leftarrow u-u^{N}$ and $ p \leftarrow p-p^{ref}$ and

\begin{displaymath}\begin{array}{rcl} g\hackscore{i} & \leftarrow & g\hackscore{...
...re{,j }\\ f & \leftarrow & f - u^{N}\hackscore{k,k} \end{array}\end{displaymath} (111)

we can assume that $ u^{N}\hackscore{i} \; n\hackscore{i}=0$ and $ p^{D}=0$. We set

$\displaystyle V = \{ q \in H^{1}(\Omega) : q=0$    on $\displaystyle \Gamma\hackscore{D} \}$ (112)

and

$\displaystyle W = \{ v \in (L^2(\Omega))^{d} : v\hackscore{k,k} \in L^2(\Omega)$    and $\displaystyle u\hackscore{i} \; n\hackscore{i} =0$    on $\displaystyle \Gamma\hackscore{N} \}$ (113)

and define the operator $ Q: V \rightarrow (L^2(\Omega))^{d}$ defined by

$\displaystyle (Qp)\hackscore{i} = \kappa\hackscore{ij} p\hackscore{,j}$ (114)

and the operator $ D: W \rightarrow L^2(\Omega)$ defined by

$\displaystyle Dv = v\hackscore{k,k}$ (115)

In operator notation the Darcy problem 6.18 is written in the form

\begin{displaymath}\begin{array}{rcl} u + Qp & = & g \\ Du & = & f \end{array}\end{displaymath} (116)

We solve this equation by minimising the functional

$\displaystyle J(u,p):=\Vert u + Qp - g\Vert^2\hackscore{0} + \Vert Du-f\Vert\hackscore{0}^2$ (117)

over $ W \times V$ where $ \Vert.\Vert\hackscore{0}$ denotes the norm in $ L^2(\Omega)$. A simple calculation shows that one has to solve

$\displaystyle ( v + Qq , u + Qp - g) + (Dv,Du-f) =0$ (118)

for all $ v\in W$ and $ q \in V$.which translates back into operator notation

\begin{displaymath}\begin{array}{rcl} (I+D^*D)u + Qp & = & D^*f + g \\ Q^*u + Q^*Q p & = & Q^*g \\ \end{array}\end{displaymath} (119)

where $ D^*$ and $ Q^*$ denote the adjoint operators. In [19] it has been shown that this problem is continuous and coercive in $ W \times V$ and therefore has a unique solution. Also standart FEM methods can be used for discretization. It is also possible to solve the problem is coupled form, however this approach leads in some cases to a very ill-conditioned stiffness matrix in particular in the case of a very small or large permability ( $ \alpha\hackscore{1} \ll 1$ or $ \alpha\hackscore{0} \gg 1$)

The approach we are taking is to eliminate the velocity $ u$ from the problem. Assuming that $ p$ is known we have

$\displaystyle v= (I+D^*D)^{-1} (D^*f + g - Qp)$ (120)

(notice that $ (I+D^*D)$ is coercive in $ W$) which is inserted into the second equation

$\displaystyle Q^* (I+D^*D)^{-1} (D^*f + g - Qp) + Q^* Q p = Q^*g$ (121)

which is

$\displaystyle Q^* ( I - (I+D^*D)^{-1} ) Q p = Q^* (g-(I+D^*D)^{-1} (D^*f + g) )$ (122)

We use the PCG method to solve this. The residual $ r$ ($ \in V^*$) is given as

\begin{displaymath}\begin{array}{rcl} r & = & Q^* \left( g -(I+D^*D)^{-1} (D^*f ...
...- Qp) \right) \\ & =& Q^* \left( g - Qp - v \right) \end{array}\end{displaymath} (123)

So in a partical implementation we use $ \hat{r}=g-Qp-v$ to represent the residual. The evaluation of the iteration operator for a given $ p$ is then returning $ Qp+v$ where $ v$ is the solution of

$\displaystyle (I+D^*D)v = Qp$ (124)

We use $ (Q^*Q)^{-1}$ as a preconditioner for the iteration operator $ Q^* ( I - (I+D^*D)^{-1} ) Q$. So the application of the preconditioner to $ \hat{r}$ representing the residual is given by solving implemented by solving

$\displaystyle Q^*Q q = Q^*\hat{r}$ (125)

The residual norm used in the PCG is given as

$\displaystyle \Vert r\Vert\hackscore{PCG}^2 = \int r \cdot (Q^*Q)^{-1} r \; dx ...
... \cdot Q (Q^*Q)^{-1} Q^* \hat{r} \; dx \approx \Vert\hat{r}\Vert\hackscore{0}^2$ (126)

The iteration is terminated if

$\displaystyle \Vert r\Vert\hackscore{PCG} \le$   ATOL (127)

where we set

ATOL$\displaystyle =$   atol$\displaystyle +$   rtol$\displaystyle \cdot \left(\frac{1}{\Vert v\Vert\hackscore{0}} + \frac{1}{\Vert Qp\Vert\hackscore{0}} \right)^{-1}$ (128)

where rtol is a given relative tolerance and atol is a given absolute tolerance (typically $ =0$). Notice that if $ Qp$ and $ v$ both are zero, the pair $ (0,p)$ is a solution. The problem is that ATOL is depending on the solution $ p$ (and $ v$ calculated form 6.31). In partcice one use the initial guess for $ p$ to get a first value for ATOL. If the stopping crierion is met in the PCG iteration, a new $ v$ is calculated from the current pressure approximation and ATOL is recalculated. If 6.38 is still fullfilled the calculation is terminated and $ (v,p)$ is returned. Otherwise PCG is restarted with a new ATOL.

esys@esscc.uq.edu.au