13.3 The augmented Lagrangian method

Next: 13.4 The streamfunction-vorticity method Up: 13 SIMPLE SOLVERS FOR Previous: 13.2 The penalty method

13.3 The augmented Lagrangian method

The augmented Lagrangian method is a refinement of the penalty method described above in which the equations

$\begin{displaymath} \mbox{\rm Re }\mbox{\boldmath{$v$}}\cdot\nabla\mbox{\boldmath{$v$}}=-\nabla p+\nabla^2\mbox{\boldmath{$v$}}\end{displaymath}$

$\begin{displaymath} \nabla\cdot\mbox{\boldmath{$v$}}=0\end{displaymath}$

are solved using the iterative procedure

$\begin{displaymath} -\mbox{\rm Pen }\nabla(\nabla\cdot\mbox{\boldmath{$v$}}_n)+\... ...\mbox{\boldmath{$v$}}_{n-1}\cdot\nabla\mbox{\boldmath{$v$}}_n=0\end{displaymath}$

$\begin{displaymath} p_n=p_{n-1}-\mbox{\rm Pen }\nabla\cdot\mbox{\boldmath{$v$}}_n\end{displaymath}$

This formulation removes the requirement that the penalty parameter be large. The algorithm, if it converges, does not introduce any further error and provides an answer for the pressure. The actual algorithm that is implemented below is obtained by dividing throughout by Re, modifying

and Pen, and setting $\Delta p = p_n - p_{n-1}$ . This results in the following equations:

$\begin{displaymath} -\mbox{\rm Pen }\nabla(\nabla\cdot\mbox{\boldmath{$v$}}_n)+\... ...\mbox{\boldmath{$v$}}_{n-1}\cdot\nabla\mbox{\boldmath{$v$}}_n=0\end{displaymath}$

$\begin{displaymath} \Delta p+\mbox{\rm Pen }\nabla\cdot\mbox{\boldmath{$v$}}_n=0\end{displaymath}$

The mesh file for the problem is cav.msh as given previously. The problem file is cav_aug.prb . The first few lines, starting with P, declare values for the parameters Re, Pen, iter and nstep.

The lines starting with D declare that boundary segments tagged 1 and 2 are to be called wall and top respectively. It should be noted that different boundary tags can be given the same name. Conversely, for different problems, the same boundary tag can be assigned different names. (For example, if in another problem, the top boundary were to be used as the bottom of a different region, the tag 2 could variously be assigned both top and bottom. The interpretation of the name would be particular to the problem.)

The two problems, augment and delp, are explained below.

The macro run is generally straightforward. The algorithm will execute nstep iterations, assembling and solving the problems augment and delp at each iteration. At each iteration, there is a printout of the L1 norm of the difference between iterations and an arrow plot of the current solution is shown on the graphics screen. Various manipulations of the vector stack are necessary. This macro is explained in greater detail below.

We now investigate the Fasttalk code in more detail. We shall assume that the last estimates of the velocity $\mbox{\boldmath{$v$}}_{n-1}$ and the pressure $p_{n-1}$ are located on V100 and V201 respectively.

   A   augment
   e   D_i{Pen}D_jU1_j+D_j{1.0/Re}D_jU1_i=     \
                           {V100}_jD_jU1_i +D_i{V201}
   b   wall U1={0.0,0.0}
   b   top  U1={1.0,0.0}

The problem augment is declared by the first line. The e statement is, as explained above, a vector equation (indexed by the i suffix) representing the momentum equation. Notice that there are unknown terms on both the RHS and LHS. Although the equations are often described as if the unknowns are on the left and the known terms are on the right, as if set up for solution, there is no syntax requirement for this in Fastflo, or even a requirement to have an equals sign at all. The basic difference between this equation and that for the penalty method is that the gradient of the pressure appears as the last term.

The Dirichlet conditions in the b statements set the velocity to zero on the walls (the no-slip condition) and to the prescribed sliding velocity on the top. The use of the symbolic terms wall and top confers no particular advantage here. In more complex programming, however, such ability to declare names for boundary segments allows the programmer to code without having any particular mesh in mind; the numerical tags in the mesh can then be linked with a set of D statements at a later stage.

   A   [reduced]   delp
   e   U1={Pen}D_j{-V100}_j

These lines describe the problem delp for the equation for the pressure. It expresses the incompressibility constraint. In the Navier-Stokes equations, it is well known that including this equation at full accuracy can lead to an over-constrained problem, with the appearance of checker-board instabilities in the pressure. A common remedy is to interpolate the pressure on a sparser grid. Fastflo supports this by linking the quadratic elements with corresponding linear elements. Here the test functions from the 6-noded elements used for velocity are replaced by a smaller set of linear shape functions drawn from 3-noded elements. This is achieved by the option [reduced]. The shape functions used for interpolating data are not changed, so the velocity that is differentiated has quadratic interpolation.

When solved, this equation produces a result on the subset of corner nodes, and this vector can only be used in expressions containing other similar vectors; furthermore, it cannot be used in graphics or as data to other problems (even reduced ones). Consequently, a command expand is provided which linearly interpolates onto the other nodes, creating a full vector which can be used for all purposes. The overall construction is equivalent to finding $\Delta p$ from the equation $\Delta p = -\mbox{\rm Pen }\nabla\cdot\mbox{\boldmath{$v$}}$ .

   < run
       while iter < nstep
         augment
         solve
         V500=V200-V100
         test=L1 V500
         show test
         delp
         solve
         expand 1
         V401=V401+V101
         popp
         V200=V100
         black
         arrow
         popp
         iter=iter+1
       endwhile
   >

The macro run assembles and solves the problems augment and delp. This takes place nstep times, as controlled by the while ... endwhile construction. Within each iteration, the parameter iter is incremented by 1. We assume now and verify later that the last estimates of $\mbox{\boldmath{$v$}}_{n-1}$ and $p_{n-1}$ are located on V100 and V201 respectively. The commands augment and solve assemble and solve the problem for the velocity components. On solving augment, the vector stack is pushed and an initial value for $\mbox{\boldmath{$v$}}_n$ is calculated and stored on V100. The current vector stack thus reads

$\begin{displaymath} \begin{array} {ccc} {\tt V100} = \mbox{\boldmath{$v$}}_n & {... ...box{\boldmath{$v$}}_{n-1} & {\tt V301} = p_{n-1} \\ \end{array}\end{displaymath}$

The three commands V500=V200-V100, test=L1 V500, show test calculate the difference between the velocity estimates, calculate the L1 norm of the difference, and print out the L1 norm.

Now it is time to update the pressure. The commands delp and solve assemble and solve the pressure update equation. After this process, the vector stack has the form

$\begin{displaymath} \begin{array} {cccc} {\tt V101} = p_n-p_{n-1} , & {\tt V200}... ...}=\mbox{\boldmath{$v$}}_{n-1} , & {\tt V401}=p_{n-1}\end{array}\end{displaymath}$

At this stage, the pressure has been calculated on 3-noded triangles, hence the use of the expand 1 command to interpolate the pressure on 6-noded triangles. Next, the command V401=V401+V101 calculates

and assigns it to V401. All that remains is to adjust the vector stack so that V100 and V201 contain the latest values for $\mbox{\boldmath{$v$}}$ and

respectively. This is accomplished by the commands popp, V200=V100, popp. The remaining commands show an arrow plot and update the parameter iter.

In fluid dynamics problems, it is often useful to calculate the streamfunction $\psi$ for the flow. This is achieved by the problem statement

   A   strf
   e   D_jD_jU1=[curlv] {-V100}
   b   top  U1={0.0}
   b   wall U1={0.0}

which is incorporated in the file cav_aug.prb . When assembled and solved, this problem finds $\psi$ from the equation $-\nabla^2\psi=\mbox{\boldmath{$k$}}\cdot\mbox{curl }\mbox{\boldmath{$v$}}$ which results from taking the curl of the defining equation $\mbox{\boldmath{$v$}}=\mbox{curl}(\mbox{\boldmath{$k$}}\psi)$ ,where $\mbox{\boldmath{$k$}}$ is the out-of-plane unit normal. Each value of $\psi$ corresponds to a streamline with the property of being everywhere tangential to the flow. As a result, fluid does not cross a streamline at any point. For steady flow, streamlines are equivalent to streak lines made visible when dye is injected into the fluid flow.

The special option [curlv] used in the e statement is exactly what is needed to express the 2D operator on the RHS of the equation. The boundary conditions are especially simple in this example, because the entire boundary is a streamline, and the streamfunction can be set there to any constant value. A more general condition uses the defining equation $\mbox{\boldmath{$v$}}=\mbox{curl}(\mbox{\boldmath{$k$}}\psi)$ in the form

$\begin{displaymath} \mbox{\boldmath{$t$}}\cdot\mbox{\boldmath{$v$}}=\mbox{\boldm... ...dmath{$t$}}\cdot\nabla\psi=\mbox{\boldmath{$n$}}\cdot\nabla\psi\end{displaymath}$

where $\mbox{\boldmath{$t$}}=-\mbox{\boldmath{$k$}}\times\mbox{\boldmath{$n$}}$ . This is a Neumann condition on $\psi$ . However, there does have to be at least one point at which $\psi$ is specified, otherwise the problem is degenerate, because any constant could be added to the solution. The simplest solution is to specify $\psi$ by a Dirichlet condition on just one contiguous wall-type boundary.

To complete this example, the following macro stream assembles the problem strf, solves it, and makes a contour plot of $\psi$ on the graphics screen. The vector stack is then popped to return the current velocity estimate to V100. The macro stream is also included in the file cav_aug.prb .

   < stream
       strf
       solve
       black
       contour 101
       popp
   >

At this stage, it is time to load the files for the augmented Lagrangian method and to run the code. Again, you might wish to experiment with different sorts of control over the iteration, different meshes, and different values for Re and Pen. [

Another interesting exercise is to try to obtain a plot of the pressure

, and to observe how sensitive

is to the accuracy of the calculation.] The converged results obtained by the augmented Lagrangian method look similar to those obtained in the previous section by the penalty method, although they differ in some details and because the Reynolds number is higher. As a concluding comment, the augmented Lagrangian method should be a more reliable algorithm than the penalty method.