Solution Method

If we assume that $\eta$ is independent from the velocity and pressure equations 6.1 and 6.2 can be written in the block form

$$\left[ \begin{array}{cc} A & B^{*} \\ B & 0 \\ \end{array} \right... ...p \\ \end{array} \right] =\left[ \begin{array}{c} G \\ 0 \\ \end{array} \right]$$ (105)

where $A$ is coercive, self-adjoint linear operator in a suitable Hilbert space, $B$ is the $(-1) \cdot$ divergence operator and $B^{*}$ is it adjoint operator (=gradient operator). For more details on the mathematics see references [11,12].

If $v\hackscore{0}$ and $p\hackscore{0}$ are given initial guesses for velocity and pressure we calculate a correction $dv$ for the velocity by solving the first equation of equation 6.5

$$A dv\hackscore{1} = G - A v\hackscore{0} - B^{*} p\hackscore{0}$$ (106)

We then insert the new approximation $v\hackscore{1}=v\hackscore{0}+dv\hackscore{1}$ to calculate a correction $dp\hackscore{2}$ for the pressure and an additional correction $dv\hackscore{2}$ for the velocity by solving

$$\begin{array}{rcl} B A^{-1} B^{*} dp\hackscore{2} & = & Bv\ha... ...1} \\ A dv\hackscore{2} & = & B^{*} dp\hackscore{2} \end{array}$$ (107)

The new velocity and pressure are then given by $v\hackscore{2}=v\hackscore{1}-dv\hackscore{2}$ and $p\hackscore{2}=p\hackscore{0}+dp\hackscore{2}$ which will fulfill the block system 6.5. This solution strategy is called the Uzawa scheme .

There is a problem with this scheme: In practice we will use an iterative scheme to solve any problem for operator $A$ . So we will be unable to calculate the operator $B A^{-1} B^{*}$ required for $dp\hackscore{2}$ explicitly. In fact, we need to use another iterative scheme to solve the first equation in 6.7 where in each iteration step an iterative solver for $A$ is applied. Another issue is the fact that the viscosity $\eta$ may depend on velocity or pressure and so we need to iterate over the three equations 6.6 and 6.7.

In the following we will use the two norms

$$\Vert v\Vert\hackscore{1}^2 = \int\hackscore{\Omega} v\hackscore{j,k}v\hackscore{j,k} \; dx \Vert p\Vert\hackscore{0}^2= \int\hackscore{\Omega} p^2 \; dx.$$ (108)

for velocity $v$ and pressure $p$ . The iteration is terminated if the stopping criteria

$$\max(\Vert Bv\hackscore{1}\Vert\hackscore{0},\Vert v\hackscore{2}... ...ore{0}\Vert\hackscore{1}) \le \tau \cdot \Vert v\hackscore{2}\Vert\hackscore{1}$$ (109)

for a given given tolerance $0<\tau<1$ is meet. Notice that because of the first equation of 6.7 we have that $\Vert Bv\hackscore{1}\Vert\hackscore{0}$ equals the norm of $B A^{-1} B^{*} dp\hackscore{2}$ and consequently provides a norm for the pressure correction.

We want to optimize the tolerance choice for solving 6.6 and 6.7. To do this we write the iteration scheme as a fixed point problem. Here we consider the errors produced by the iterative solvers being used. From Equation 6.6 we have

$$v\hackscore{1} = e\hackscore{1} + v\hackscore{0} + A^{-1} ( G - Av\hackscore{0} - B^{*} p\hackscore{0} )$$ (110)

where $e\hackscore{1}$ is the error when solving 6.6. We will use a sparse matrix solver so we have not full control on the norm $\Vert.\Vert\hackscore{s}$ used in the stopping criteria for this equation. In fact we will have a stopping criteria of the form

$$\Vert A e\hackscore{1} \Vert\hackscore{s} = \Vert G - A v\hacksco... ...ckscore{1} \Vert G - A v\hackscore{0} - B^{*} p\hackscore{0} \Vert\hackscore{s}$$ (111)

where $\tau\hackscore{1}$ is the tolerance which we need to choose. This translates into the condition

$$\Vert e\hackscore{1} \Vert\hackscore{1} \le K \tau\hackscore{1} \Vert dv\hackscore{1} - e\hackscore{1} \Vert\hackscore{1}$$ (112)

The constant $K$ represents some uncertainty combining a variety of unknown factors such as the norm being used and the condition number of the stiffness matrix. From the first equation of 6.7 we have

$$p\hackscore{2} = p\hackscore{0} + (B A^{-1} B^{*})^{-1} (e\hackscore{2} + Bv\hackscore{1} )$$ (113)

where $e\hackscore{2}$ represents the error when solving 6.7. We use an iterative preconditioned conjugate gradient method (PCG) with iteration operator $B A^{-1} B^{*}$ using the $\Vert.\Vert\hackscore{0}$ norm. As suitable preconditioner for the iteration operator we use $\frac{1}{\eta}$ [16], ie the evaluation of the preconditioner $P$ for a given pressure increment $q$ is the solution of

$$\frac{1}{\eta} (Pq) = q \; .$$ (114)

Note that in each evaluation of the iteration operator $q=B A^{-1} B^{*} s$ one needs to solve the problem

$$A w = B^{*} s$$ (115)

with sufficient accuracy to return $q=Bw$ . We assume that the desired tolerance is sufficiently small, for instance one can take $\tau\hackscore{2}^2$ where $\tau\hackscore{2}$ is the tolerance for 6.7.

In an implementation we use the fact that the residual $r$ is given as

$$r= B (v\hackscore{1} - A^{-1} B^{*} dp) = B (v\hackscore{1} - A^{-1} B^{*} dp) = B (v\hackscore{1}-dv\hackscore{2}) = B v\hackscore{2}$$ (116)

In particular we have $e\hackscore{2} = B v\hackscore{2}$ So the residual $r$ is represented by the updated velocity $v\hackscore{2}=v\hackscore{1}-dv\hackscore{2}$ . In practice, if one uses the velocity to represent the residual $r$ there is no need to recover $dv\hackscore{2}$ in 6.7 after $dp\hackscore{2}$ has been calculated. In PCG the iteration is terminated if

$$\Vert P^{\frac{1}{2}}B v\hackscore{2} \Vert\hackscore{0} \le \tau\hackscore{2} \Vert P^{\frac{1}{2}}B v\hackscore{1} \Vert\hackscore{0}$$ (117)

where $\tau\hackscore{2}$ is the given tolerance. This translates into

$$\Vert e\hackscore{2}\Vert\hackscore{0} = \Vert B v\hackscore{2} \... ...\hackscore{0} \le M \tau\hackscore{2} \Vert B v\hackscore{1} \Vert\hackscore{0}$$ (118)

where $M$ is taking care of the fact that $P^{\frac{1}{2}}$ is dropped.

As we assume that there is no significant error from solving with the operator $A$ we have

$$v\hackscore{2} = v\hackscore{1} - dv\hackscore{2} = v\hackscore{1} - A^{-1} B^{*}dp$$ (119)

Combining the equations 6.10, 6.13 and 6.19 and setting the errors to zero we can write the solution process as a fix point problem

$$v = \Phi(v,p) p = \Psi(u,p)$$ (120)

with suitable functions $\Phi(v,p)$ and $\Psi(v,p)$ representing the iteration operator without errors. In fact for a linear problem, $\Phi$ and $\Psi$ are constant. With this notation we can write the update step in the form $p\hackscore{2}= \delta p + \Psi(v\hackscore{0},p\hackscore{0})$ and $v\hackscore{2}= \delta v + \Phi(v\hackscore{0},p\hackscore{0})$ where the total error $\delta p$ and $\delta v$ are given as

$$\begin{array}{rcl} \delta p & = & (B A^{-1} B^{*})^{-1} ( e\h... ...a v & = & e\hackscore{1} - A^{-1} B^{*}\delta p \;. \end{array}$$ (121)

Notice that $B\delta v = - e\hackscore{2}=-Bv\hackscore{2}$ . Our task is now to choose the tolerances $\tau\hackscore{1}$ and $\tau\hackscore{2}$ such that the global errors $\delta p$ and $\delta v$ do not stop the convergence of the iteration process.

To measure convergence we use

$$\epsilon = \max(\Vert v\hackscore{2}-v\Vert\hackscore{1}, \Vert B A^{-1} B^{*} (p\hackscore{2}-p)\Vert\hackscore{0})$$ (122)

In practice using the fact that $B A^{-1} B^{*} (p\hackscore{2}-p\hackscore{0}) = B v\hackscore{1}$ and assuming that $v\hackscore{2}$ gives a better approximation to the true $v$ than $v\hackscore{0}$ we will use the estimate

$$\epsilon = \max(\Vert v\hackscore{2}-v\hackscore{0}\Vert\hackscore{1}, \Vert B v\hackscore{1}\Vert\hackscore{0})$$ (123)

to estimate the progress of the iteration step after the step is completed. Note that the estimate of $\epsilon$ used in the stopping criteria 6.9. If $\chi^{-}$ is the convergence rate assuming exact calculations, i.e. $e\hackscore{1}=0$ and $e\hackscore{2}=0$ , we are expecting to maintain $\epsilon \le \chi^{-} \cdot \epsilon^{-}$ . For the pressure increment we get:

$$\begin{array}{rcl} \Vert B A^{-1} B^{*} (p\hackscore{2}-p)\Ve... ...core{2} \Vert B v\hackscore{1}\Vert\hackscore{0} \\ \end{array}$$ (124)

So we choose the value for $\tau\hackscore{2}$ from

$$M \tau\hackscore{2} \Vert B v\hackscore{1}\Vert\hackscore{0} \le (\chi^{-})^2 \epsilon^{-}$$ (125)

in order to make the perturbation for the termination of the pressure iteration a second order effect. We use a similar argument for the velocity:

$$\begin{array}{rcl} \Vert v\hackscore{2}-v\Vert\hackscore{1} &... ... + K \tau\hackscore{1}) \chi^{-} \cdot \epsilon^{-} \end{array}$$ (126)

So we choose the value for $\tau\hackscore{1}$ from

$$K \tau\hackscore{1} \le \chi^{-}$$ (127)

Assuming we have estimates for $M$ and $K$ (if none is available, we use the value $1$ .) we can use 6.27 and 6.25 to get appropriate values for the tolerances. After the step has been completed we can calculate a new convergence rate $\chi =\frac{\epsilon}{\epsilon^{-}}$ . For partial reasons we restrict $\chi$ to be less or equal a given maximum value $\chi\hackscore{max}\le 1$ . If we see $\chi \le \chi^{-} (1+\chi^{-})$ our choices for the tolerances was suitable. Otherwise, we need to adjust the values for $K$ and $M$ . From the estimates 6.24 and 6.26 we establish

$$\chi \le ( 1 + \max(M \frac{\tau\hackscore{2} \Vert B v\hackscore... ...ert\hackscore{0}}{\chi^{-} \epsilon^{-}},K \tau\hackscore{1} ) ) \cdot \chi^{-}$$ (128)

If we assume that this inequality would be an equation if we would have chosen the right values $M^{+}$ and $K^{+}$ then we get

$$\chi = ( 1 + \max(M^{+} \frac{\chi^{-}}{M},K^{+} \frac{\chi^{-}}{K}) ) \cdot \chi^{-}$$ (129)

From this equation we set for than our choice for $K$ was not good enough. In this case we can calculate a new value

$$K^{+} = \frac{\chi-\chi^{-}}{(\chi^{-})^2} K$$ (130)

In practice we will use

$$K^{+} = \max(\frac{\chi-\chi^{-}}{(\chi^{-})^2} K,\frac{1}{2}K,1)$$ (131)

where the second term is used to reduce a potential overestimate of $K$ . The same identity is used for to update $M$ . The updated $M^{+}$ and $K^{+}$ are then use in the next iteration step to control the tolerances.

In some cases one can observe that there is a significant change in the velocity but the new velocity $v\hackscore{1}$ has still a small divergence, i.e. we have $\Vert Bv\hackscore{1}\Vert\hackscore{0} \ll \Vert v\hackscore{1}-v\hackscore{0}\Vert\hackscore{1}$ . In this case we will get a small pressure increment and consequently only very small changes to the velocity as a result of the second update step which therefore can be skipped and we can directly repeat the first update step until the increment in velocity becomes significant relative to its divergence. In practice we will ignore the second half of the iteration step as long as

$$\Vert Bv\hackscore{1}\Vert\hackscore{0} \le \theta \cdot \Vert v\hackscore{1}-v\hackscore{0}\Vert$$ (132)

where $0<\theta<1$ is a given factor. In this case we will also check the stopping criteria with $v\hackscore{1}\rightarrow v\hackscore{2}$ but we will not correct $M$ in this case.

Starting from an initial guess $v\hackscore{0}$ and $p\hackscore{0}$ for velocity and pressure the solution procedure is implemented as follows.

1. calculate viscosity $\eta(v\hackscore{0},p)\hackscore{0}$ and assemble operator $A$ from $\eta$ .
2. calculate the tolerance $\tau\hackscore{1}$ from 6.27.
3. Solve 6.6 for $dv\hackscore{1}$ with tolerance $\tau\hackscore{1}$ .
4. update $v\hackscore{1}=v\hackscore{0}+dv\hackscore{1}$
5. if $Bv\hackscore{1}$ is large (see 6.32):
   1. calculate the tolerance $\tau\hackscore{2}$ from 6.25.
   2. Solve 6.7 for $dp\hackscore{2}$ and $v\hackscore{2}$ with tolerance $\tau\hackscore{2}$ .
   3. update $p\hackscore{2}\leftarrow p\hackscore{0}+ dp\hackscore{2}$
6. else:
   1. update $p\hackscore{2}\leftarrow p$ and $v\hackscore{2}\leftarrow v\hackscore{1}$
7. calculate convergence measure $\epsilon$ and convergence rate $\chi$
8. if stopping criteria 6.9 holds:
   1. return $v\hackscore{2}$ and $p\hackscore{2}$
9. else:
   1. update $M$ and $K$
   2. goto step 1 with $v\hackscore{0}\leftarrow v\hackscore{2}$ and $p\hackscore{0}\leftarrow p\hackscore{2}$ .