Solution Method

It is useful to write equation 6.33 in operator form. For any pressure $p$ we set

$$(Gp)\hackscore{i} = p\hackscore{,j}$$ (136)

and a velocity $v$ we set

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

We $K=(\kappa\hackscore{ij})$ we can write the Darcy problem 6.33 as

$$\begin{array}{rcl} u + K \, Gp & = & g \\ Du & = & f \end{array}$$ (138)

We solve this equation by minimising the functional

$$J(u,p):=\Vert K^{-\frac{1}{2}}(u + K \, G p - g)\Vert^2\hackscore{0} + \Vert\lambda (Du-f) \Vert\hackscore{0}^2$$ (139)

over all suitable $u$ and $p$ . In this equation we set $\Vert p\Vert^2\hackscore{0}=(p,p)\hackscore{0}$ with

$$(p,q)\hackscore{0} = \int\hackscore{\Omega } p\cdot q \, dx$$ (140)

The factor $\lambda>0$ is a weighting factor. A simple calculation shows that one has to solve

$$( K^{-1} (v + K \, Gq) , u +K \, G p - g)\hackscore{0} + (\lambda Dv,\lambda (Du-f) )\hackscore{0} =0$$ (141)

for all velocities $v$ and pressure $q$ which fulfill the homogeneous boundary conditions 6.34. This so-called variational equation can be translated back into operator notation

$$\begin{array}{rcl} (K^{-1}+ D^*\lambda^2 D)u + Gp & = & D^*\lambda f + K^{-1} g \\ G^*u + G^*K \, G p & = & G^*g \\ \end{array}$$ (142)

where $D^*$ and $Q^*$ denote the adjoint operators with respect to $(.,.)\hackscore{0}$ . In [28] it has been shown that this problem is continuous and coercive and therefore has a unique solution. Also standard FEM methods can be used for discretization. It is also possible to solve the problem in 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 permeability ( $\alpha\hackscore{1} \ll 1$ or $\alpha\hackscore{0} \gg 1$ )

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

$$v= (K^{-1}+ D^*\lambda^2 D)^{-1} ( D^*\lambda f + K^{-1} g - Gp)$$ (143)

(notice that $K^{-1}+\lambda D^*D$ is coercive) which is inserted into the second equation

$$G^* (K^{-1}+\lambda D^*D)^{-1} (\lambda D^*f + K^{-1} g - Gp) + G^* KG p = G^*g$$ (144)

which is

$$G^* ( K - (K^{-1}+ D^*\lambda^2 D)^{-1} ) G p = G^* (g-(K^{-1}+D^*\lambda^2 D)^{-1} ( D^*\lambda f + K^{-1} g) )$$ (145)

We use the PCG method to solve this. The residual $r$ is given as

$$\begin{array}{rcl} r & =& G^* \left( g - K\, G p - v \right) \end{array}$$ (146)

for the current pressure approximation $p$ and current velocity $v$ defined by equation 6.43. So in a particular implementation we use $\hat{r}=g-K\, Gp-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

$$(K^{-1}+ D^*\lambda^2 D)v = Gp$$ (147)

To derive a preconditioner we use the identity

$$\begin{array}{rcl} \par G^* ( K - (K^{-1}+ D^*\lambda^2 D)^{-... ...K G \\ & \approx & G^* \frac{\lambda^2}{dx^2} K^2 G \end{array}$$ (148)

where $dx$ is the local mesh size and we use the approximation

$$D^* \lambda^2 D \approx \frac{\lambda^2}{dx^2}$$ (149)

Therefore we use $G^* \frac{\lambda^2}{dx^2} K^2 G$ as a preconditioner.To evalute the preconditioner we need to solve the equation

$$G^* \frac{\lambda^2}{dx^2} K^2 G p = G^* \hat{r}$$ (150)

It remains to answer the question how to choose $\lambda$ . We need to balance the first and second term in $J(u,p)$ in equation 6.39. We inspect $J$ for $(\hat{u}, \hat{p})$ which is a perturbed exact solution $(u,p)$ . Assuming $\hat{u}=u+u\hackscore{0}e^{ik^tx}$ and $\hat{p}=p+p\hackscore{0}e^{ik^tx}$ and constant $K$ we get

$$J(\hat{u},\hat{p}) = C \left[ ( \Vert K^{-1}\Vert\hackscore{2} \v... ...t ) + \lambda^2 \Vert k\Vert\hackscore{2}^2 \vert u\hackscore{0}\vert^2 \right]$$ (151)

with some constant $C>0$ . The first two terms and the third term correspond to the first term and second term in the definition of $J(u,p)$ in equation 6.39. For small $\Vert k\Vert\hackscore{2}$ (i.e. for a smooth perturbation) $J(\hat{u},\hat{p})$ is dominated by $\Vert K^{-1}\Vert\hackscore{2} \vert u\hackscore{0}\vert^2$ . To scale the second term which is corresponds to the incompressibility condition for the velocity we need to meet the condition $\Vert K^{-1}\Vert\hackscore{2} = \lambda^2 \Vert k\Vert\hackscore{2}^2$ . Taking the boundary conditions into consideration the smallest possible value for $\Vert k\Vert\hackscore{2} = \frac{\pi}{l}$ where $l$ is the longest edge of the domain. This leads to the

$$\lambda = \Vert K^{-1}\Vert\hackscore{2}^{\frac{1}{2}} \frac{l}{\pi}$$ (152)

Notice that with this setting the preconditioner $G^* \frac{\lambda^2}{dx^2} K^2 G$ becomes equivalent to $G^* K G$ if $K$ is a diagonal matrix and the mesh has a constant mesh size.

The residual norm used in the PCG is given as

$$\Vert r\Vert\hackscore{PCG}^2 = \int r (G^* \frac{\lambda^2}{dx^2... ...}{dx^2} K^2 G)^{-1} G^* \hat{r} \; dx \approx \int \hat{r} K^{-1} \hat{r} \; dx$$ (153)

The iteration is terminated if

$$\Vert r\Vert\hackscore{PCG} \le$$ (154)

where we set

$$= + \cdot \left(\frac{1}{\Vert K^{-\frac{1}{2}}v\Vert\hackscore{0}} + \frac{1}{\Vert K^{\frac{1}{2}} G p\Vert\hackscore{0}} \right)^{-1}$$ (155)

where rtol is a given relative tolerance and atol is a given absolute tolerance (typically $=0$ ). Notice that if $Gp$ 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.43. In practice one use the initial guess for $p$ to get a first value for ATOL. If the stopping criterion is met in the PCG iteration, a new $v$ is calculated from the current pressure approximation and ATOL is recalculated. If 6.54 is still fulfilled the calculation is terminated and $(v,p)$ is returned. Otherwise PCG is restarted with a new ATOL.