# Demo - Lid driven cavity¶

- Authors
Mikael Mortensen (mikaem at math.uio.no)

- Date
Nov 14, 2019

*Summary.* The lid driven cavity is a classical benchmark for Navier Stokes solvers.
This is a demonstration of how the Python module shenfun can be used to solve the lid
driven cavity problem with full spectral accuracy using a mixed (coupled) basis
in a 2D tensor product domain. The demo also shows how to use mixed
tensor product spaces for vector valued equations. Note that the regular
lid driven cavity, where the top wall has constant velocity and the
remaining three walls are stationary, has a singularity at the two
upper corners, where the velocity is discontinuous.
Due to their global nature, spectral methods
are usually not very good at handling problems with discontinuities, and
for this reason we will also look at a regularized lid driven cavity,
where the top lid moves according to \((1-x)^2(1+x)^2\), thus removing
the corner discontinuities.

## Bases and tensor product spaces¶

With the Galerkin method we need basis functions for both velocity and pressure, as well as for the nonlinear right hand side. A Dirichlet basis will be used for velocity, whereas there is no boundary restriction on the pressure basis. For both two-dimensional bases we will use one basis function for the \(x\)-direction, \(\mathcal{X}_k(x)\), and one for the \(y\)-direction, \(\mathcal{Y}_l(y)\). And then we create two-dimensional basis functions like

and solutions (trial functions) as

For the homogeneous Dirichlet boundary condition the basis functions \(\mathcal{X}_k(x)\) and \(\mathcal{Y}_l(y)\) are chosen as composite Legendre polynomials (we could also use Chebyshev):

where \(\mathbf{k}^{N_0-2} = (0, 1, \ldots, N_0-3)\), \(\mathbf{l}^{N_1-2} = (0, 1, \ldots, N_1-3)\) and \(N = (N_0, N_1)\) is the number of quadrature points in each direction. Note that \(N_0\) and \(N_1\) do not need to be the same. The basis funciton (3) satisfies the homogeneous Dirichlet boundary conditions at \(x=\pm 1\) and (4) the same at \(y=\pm 1\). As such, the basis function \(v_{kl}(x, y)\) satisfies the homogeneous Dirichlet boundary condition for the entire domain.

With shenfun we create these homogeneous spaces, \(D_0^{N_0}(x)=\text{span}\{L_k-L_{k+2}\}_{k=0}^{N_0-2}\) and \(D_0^{N_1}(y)=\text{span}\{L_l-L_{l+2}\}_{l=0}^{N_1-2}\) as

```
N = (51, 51)
family = 'Legendre' # or use 'Chebyshev'
quad = 'LG' # for Chebyshev use 'GC' or 'GL'
D0X = Basis(N[0], family, quad=quad, bc=(0, 0))
D0Y = Basis(N[1], family, quad=quad, bc=(0, 0))
```

The spaces are here the same, but we will use `D0X`

in the \(x\)-direction and
`D0Y`

in the \(y\)-direction. But before we use these bases in
tensor product spaces, they remain identical as long as \(N_0 = N_1\).

Special attention is required by the moving lid. To get a solution with nonzero boundary condition at \(y=1\) we need to add one more basis function that satisfies that solution. In general, a nonzero boundary condition can be added on both sides of the domain using the following basis

And then the unknown component \(N_1-2\) decides the value at \(y=1\), whereas the unknown at \(N_1-1\) decides the value at \(y=-1\). Here we only need to add the \(N_1-2\) component, but for generality this is implemented in shenfun using both additional basis functions. We create the space \(D_1^{N_1}(y)=\text{span}\{\mathcal{Y}_l(y)\}_{l=0}^{N_1-1}\) as

```
D1Y = Basis(N[1], family, quad=quad, bc=(1, 0))
```

where `bc=(1, 0)`

fixes the values for \(y=1\) and \(y=-1\), respectively.
For a regularized lid driven cavity the velocity of the top lid is
\((1-x)^2(1+x)^2\) and not unity. To implement this boundary condition
instead, we can make use of sympy and
quite straight forward do

```
import sympy
x = sympy.symbols('x')
#D1Y = Basis(N[1], family, quad=quad, bc=((1-x)**2*(1+x)**2, 0))
```

Uncomment the last line to run the regularized boundary conditions. Otherwise, there is no difference at all between the regular and the regularized lid driven cavity implementations.

The pressure basis that comes with no restrictions for the boundary is a little trickier. The reason for this has to do with inf-sup stability. The obvious choice of basis functions are the regular Legendre polynomials \(L_k(x)\) in \(x\) and \(L_l(y)\) in the \(y\)-directions. The problem is that for the natural choice of \((k, l) \in \mathbf{k}^{N_0} \times \mathbf{l}^{N_1}\) there are nullspaces and the problem is not well-defined. It turns out that the proper choice for the pressure basis is simply the regular Legendre basis functions, but for \((k, l) \in \mathbf{k}^{N_0-2} \times \mathbf{l}^{N_1-2}\). The bases \(P^{N_0}(x)=\text{span}\{L_k(x)\}_{k=0}^{N_0-3}\) and \(P^{N_1}(y)=\text{span}\{L_l(y)\}_{l=0}^{N_1-3}\) are created as

```
PX = Basis(N[0], family, quad=quad)
PY = Basis(N[1], family, quad=quad)
PX.slice = lambda: slice(0, N[0]-2)
PY.slice = lambda: slice(0, N[1]-2)
```

Note that we still use these spaces with the same \(N_0 \cdot N_1\) quadrature points in real space, but the two highest frequencies have been set to zero.

We have now created all relevant function spaces for the problem at hand. It remains to combine these spaces into tensor product spaces, and to combine tensor product spaces into mixed (coupled) tensor product spaces. From the Dirichlet bases we create two different tensor product spaces, whereas one is enough for the pressure

With shenfun the tensor product spaces are created as

```
V1 = TensorProductSpace(comm, (D0X, D1Y))
V0 = TensorProductSpace(comm, (D0X, D0Y))
P = TensorProductSpace(comm, (PX, PY))
```

These tensor product spaces are all scalar valued. The velocity is a vector, and a vector requires a mixed basis like \(W_1^{\mathbf{N}} = V_1^{\mathbf{N}} \times V_0^{\mathbf{N}}\). The mixed basis is created in shenfun as

```
W1 = MixedTensorProductSpace([V1, V0])
W0 = MixedTensorProductSpace([V0, V0])
```

Note that the second mixed basis, \(W_0^{\mathbf{N}} = V_0^{\mathbf{N}} \times V_0^{\mathbf{N}}\), uses homogeneous boundary conditions throughout.

## Mixed variational form¶

We now formulate a variational problem using the Galerkin method: Find \(\mathbf{u} \in W_1^{\mathbf{N}}\) and \(p \in P^{\mathbf{N}}\) such that

Note that we are using test functions \(\mathbf{v}\) with homogeneous boundary conditions.

The first obvious issue with Eq (11) is the nonlinearity. In other words we will need to linearize and iterate to be able to solve these equations with the Galerkin method. To this end we will introduce the solution on iteration \(k \in [0, 1, \ldots]\) as \(\mathbf{u}^k\) and compute the nonlinearity using only known solutions \(\int_{\Omega} (\nabla \cdot \mathbf{u}^k\mathbf{u}^k) \cdot \mathbf{v}\, dxdy\). Using further integration by parts we end up with the equations to solve for iteration number \(k+1\) (using \(\mathbf{u} = \mathbf{u}^{k+1}\) and \(p=p^{k+1}\) for simplicity)

Note that the nonlinear term may also be integrated by parts and evaluated as \(\int_{\Omega}-\mathbf{u}^k\mathbf{u}^k \, \colon \nabla \mathbf{v} \, dxdy\). All boundary integrals disappear since we are using test functions with homogeneous boundary conditions.

Since we are to solve for \(\mathbf{u}\) and \(p\) at the same time, we formulate a mixed (coupled) problem: find \((\mathbf{u}, p) \in W_1^{\mathbf{N}} \times P^{\mathbf{N}}\) such that

where bilinear (\(a\)) and linear (\(L\)) forms are given as

Note that the bilinear form will assemble to a block matrix, whereas the right hand side linear form will assemble to a block vector. The bilinear form does not change with the solution and as such it does not need to be reassembled inside an iteration loop.

The algorithm used to solve the equations are:

Set \(k = 0\)

Guess \(\mathbf{u}^0 = (0, 0)\)

while not converged:

assemble \(L((\mathbf{v}, q); \mathbf{u}^{k})\)

solve \(a((\mathbf{u}, p), (\mathbf{v}, q)) = L((\mathbf{v}, q); \mathbf{u}^{k})\) for \(\mathbf{u}^{k+1}, p^{k+1}\)

compute error = \(\int_{\Omega} (\mathbf{u}^{k+1}-\mathbf{u}^{k})^2 \, dxdy\)

if error \(<\) some tolerance then converged = True

\(k\) += \(1\)

## Implementation of solver¶

We will now implement the coupled variational problem described in previous sections. First of all, since we want to solve for the velocity and pressure in a coupled solver, we have to create a mixed tensor product space \(VQ = W_1^{\mathbf{N}} \times P^{\mathbf{N}}\) that couples velocity and pressure

```
VQ = MixedTensorProductSpace([W1, P]) # Coupling velocity and pressure
```

We can now create test- and trialfunctions for the coupled space \(VQ\), and then split them up into components afterwards:

```
up = TrialFunction(VQ)
vq = TestFunction(VQ)
u, p = up
v, q = vq
```

Note

The test function `v`

is using homogeneous Dirichlet boundary conditions even
though it is derived from `VQ`

, which contains `W1`

. It is currently not (and will
probably never be) possible to use test functions with inhomogeneous
boundary conditions.

With the basisfunctions in place we may assemble the different blocks of the final coefficient matrix. For this we also need to specify the kinematic viscosity, which is given here in terms of the Reynolds number:

```
Re = 100.
nu = 2./Re
A = inner(grad(v), -nu*grad(u))
G = inner(div(v), p)
D = inner(q, div(u))
```

Note

The inner products may also be assembled with one single line, as

```
AA = inner(grad(v), -nu*grad(u)) + inner(div(v), p) + inner(q, div(u))
```

But note that this requires addition, not subtraction, of inner products,
and it is not possible to move the negation to `-inner(grad(v), nu*grad(u))`

.
This is because the `inner()`

function returns a list of
tensor product matrices of type `TPMatrix`

, and you cannot
negate a list.

The assembled subsystems `A, G`

and `D`

are lists containg the different blocks of
the complete, coupled, coefficient matrix. `A`

actually contains 4
tensor product matrices of type `TPMatrix`

. The first two
matrices are for vector component zero of the test function `v[0]`

and
trial function `u[0]`

, the
matrices 2 and 3 are for components 1. The first two matrices are as such for

```
A[0:2] = inner(grad(v[0]), -nu*grad(u[0]))
```

Breaking it down the inner product is mathematically

We can now insert for test function \(\mathbf{v}[0]\)

and trialfunction

where \(\hat{\mathbf{u}}\) are the unknown degrees of freedom for the velocity vector. Notice that the sum over the second index runs all the way to \(N_1-1\), whereas the other indices runs to either \(N_0-3\) or \(N_1-3\). This is because of the additional basis functions required for the inhomogeneous boundary condition.

Inserting for these basis functions into (17), we obtain after a few trivial manipulations

We see that each tensor product matrix (both A[0] and A[1]) is composed as outer products of two smaller matrices, one for each dimension. The first tensor product matrix, A[0], is

where \(C\in \mathbb{R}^{N_0-2 \times N_1-2}\) and \(F \in \mathbb{R}^{N_0-2 \times N_1}\). Note that due to the inhomogeneous boundary conditions this last matrix \(F\) is actually not square. However, remember that all contributions from the two highest degrees of freedom (\(\hat{\mathbf{u}}[0]_{m,N_1-2}\) and \(\hat{\mathbf{u}}[0]_{m,N_1-1}\)) are already known and they can, as such, be moved directly over to the right hand side of the linear algebra system that is to be solved. More precisely, we can split the tensor product matrix into two contributions and obtain

where the first term on the right hand side is square and the second term is known and can be moved to the right hand side of the linear algebra equation system.

All the parts of the matrices that are to be moved to the right hand side can be extracted from A, G and D as follows

```
# Extract the boundary matrices
bc_mats = extract_bc_matrices([A, G, D])
```

These matrices are applied to the solution below (see `BlockMatrix BM`

).
Furthermore, this leaves us with square submatrices (A, G, D), which make up a
symmetric block matrix

This matrix, and the matrix responsible for the boundary degrees of freedom, can be assembled from the pieces we already have as

```
M = BlockMatrix(A+G+D)
BM = BlockMatrix(bc_mats)
```

We now have all the matrices we need in order to solve the Navier Stokes equations. However, we also need some work arrays for iterations and we need to assemble the constant boundary contribution to the right hand side

```
# Create Function to hold solution
uh_hat = Function(VQ)
ui_hat = uh_hat[0]
D1Y.bc.apply_after(ui_hat[0], True) # Fixes the values of the boundary dofs
# New solution (iterative)
uh_new = Function(VQ)
ui_new = uh_new[0]
D1Y.bc.apply_after(ui_new[0], True)
# Compute the constant contribution to rhs due to nonhomogeneous boundary conditions
bh_hat0 = Function(VQ)
bh_hat0 = BM.matvec(-uh_hat, bh_hat0) # Negative because moved to right hand side
bi_hat0 = bh_hat0[0]
```

Note that `bh_hat0`

now contains the part of the right hand side that is
due to the non-symmetric part of assembled matrices. The line with
`D1Y.bc.apply_after(ui_hat[0], True)`

ensures the known boundary values of
the solution are fixed for `ui_hat`

.

The nonlinear right hand side also requires some additional attention. Nonlinear terms are usually computed in physical space before transforming to spectral. For this we need to evaluate the velocity vector on the quadrature mesh. We also need a rank 2 Array to hold the outer product \(\mathbf{u}\mathbf{u}\). The required arrays and spaces are created as

```
bh_hat = Function(VQ)
# Create arrays to hold velocity vector solution
ui = Array(W1)
# Create work arrays for nonlinear part
QT = MixedTensorProductSpace([W1, W0]) # for uiuj
uiuj = Array(QT)
uiuj_hat = Function(QT)
```

The right hand side \(L((\mathbf{v}, q);\mathbf{u}^{k});\) is computed in its
own function `compute_rhs`

as

```
def compute_rhs(ui_hat, bh_hat):
global ui, uiuj, uiuj_hat, V1, bh_hat0
bh_hat.fill(0)
ui = W1.backward(ui_hat, ui)
uiuj = outer(ui, ui, uiuj)
uiuj_hat = uiuj.forward(uiuj_hat)
bi_hat = bh_hat[0]
#bi_hat = inner(v, div(uiuj_hat), output_array=bi_hat)
bi_hat = inner(grad(v), -uiuj_hat, output_array=bi_hat)
bh_hat += bh_hat0
return bh_hat
```

Here `outer()`

is a shenfun function that computes the
outer product of two vectors and returns the product in a rank two
array (here `uiuj`

). With `uiuj`

forward transformed to `uiuj_hat`

we can assemble the linear form either as `inner(v, div(uiuj_hat)`

or
`inner(grad(v), -uiuj_hat)`

. Also notice that the constant contribution
from the inhomogeneous boundary condition, `bh_hat0`

,
is added to the right hand side vector.

Now all that remains is to guess an initial solution and solve iteratively until convergence. For initial solution we simply set the velocity and pressure to zero and solve the Stokes equations:

```
from scipy.sparse.linalg import splu
uh_hat, Ai = M.solve(bh_hat0, u=uh_hat, constraints=((2, 0, 0),), return_system=True) # Constraint for component 2 of mixed space
Alu = splu(Ai)
uh_new[:] = uh_hat
```

Note that the `BlockMatrix`

given by `M`

has a solve method that sets up
a sparse coefficient matrix `Ai`

of size \(\mathbb{R}^{3(N_0-2)(N_1-2) \times 3(N_0-2)(N_1-2)}\),
and then solves using scipy.sparse.linalg.spsolve.
The matrix `Ai`

is then pre-factored for reuse with splu.
Also note that the `constraints=((2, 0, 0),)`

keyword argument
ensures that the pressure integrates to zero, i.e., \(\int_{\Omega} pdxdy=0\).
Here the number 2 tells us that block component 2 in the mixed space
(the pressure) should be integrated, dof 0 should be fixed, and it
should be fixed to 0.

With an initial solution from the Stokes equations we are ready to start iterating. However, for convergence it is necessary to add some underrelaxation \(\alpha\), and update the solution each time step as

where \(\hat{\mathbf{u}}^*\) and \(\hat{p}^*\) are the newly computed velocity
and pressure returned from `M.solve`

. Without underrelaxation the solution
will quickly blow up. The iteration loop goes as follows

```
converged = False
count = 0
alfa = 0.5
while not converged:
count += 1
bh_hat = compute_rhs(ui_hat, bh_hat)
uh_new = M.solve(bh_hat, u=uh_new, constraints=((2, 0, 0),), Alu=Alu) # Constraint for component 2 of mixed space
error = np.linalg.norm(ui_hat-ui_new)
uh_hat[:] = alfa*uh_new + (1-alfa)*uh_hat
converged = abs(error) < 1e-10 or count >= 10000
print('Iteration %d Error %2.4e' %(count, error))
up = uh_hat.backward()
u, p = up
X = V0.local_mesh(True)
plt.figure()
plt.quiver(X[0], X[1], u[0], u[1])
```

The last three lines plots the velocity vectors that are shown in Figure Velocity vectors for Re=100. The solution is apparently nice and smooth, but hidden underneath are Gibbs oscillations from the corner discontinuities. This is painfully obvious when switching from Legendre to Chebyshev polynomials. With Chebyshev the same plot looks like Figure Velocity vectors for Re=100 using Chebyshev. However, choosing instead the regularized lid, the solutions will be nice and smooth, both for Legendre and Chebyshev polynomials.

## Complete solver¶

A complete solver can be found in demo NavierStokesDrivenCavity.py.