# Demo - Lid driven cavity **Mikael Mortensen** (email: `mikaem@math.uio.no`), Department of Mathematics, University of Oslo. Date: **Jun 23, 2020** Copyright 2020, Mikael Mortensen. Released under CC Attribution 4.0 license **Summary.** The lid driven cavity is a classical benchmark for Navier Stokes solvers. This is a demonstration of how the Python module [shenfun](https://github.com/spectralDNS/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.

Velocity vectors for $Re=100$.

Navier Stokes equations

The nonlinear steady Navier Stokes equations are given in strong form as

$$ \begin{align*} \nu \nabla^2 \boldsymbol{u} - \nabla p &= \nabla \cdot \boldsymbol{u} \boldsymbol{u} \quad \text{in } \Omega , \\ \nabla \cdot \boldsymbol{u} &= 0 \quad \text{in } \Omega \\ \int_{\Omega} p dx &= 0 \\ \boldsymbol{u}(x, y=1) = (1, 0) \, &\text{ or }\, \boldsymbol{u}(x, y=1) = ((1-x)^2(1+x)^2, 0) \\ \boldsymbol{u}(x, y=-1) &= (0, 0) \\ \boldsymbol{u}(x=\pm 1, y) &= (0, 0) \end{align*} $$

where $\boldsymbol{u}, p$ and $\nu$ are, respectively, the fluid velocity vector, pressure and kinematic viscosity. The domain $\Omega = [-1, 1]^2$ and the nonlinear term $\boldsymbol{u} \boldsymbol{u}$ is the outer product of vector $\boldsymbol{u}$ with itself. Note that the final $\int_{\Omega} p dx = 0$ is there because there is no Dirichlet boundary condition on the pressure and the system of equations would otherwise be ill conditioned.

We want to solve these steady nonlinear Navier Stokes equations with the Galerkin method, using the shenfun Python package. The first thing we need to do then is to import all of shenfun's functionality

Note that MPI for Python (mpi4py) is a requirement for shenfun, but the current solver cannot be used with more than one processor.

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

$$ \begin{equation} v_{kl}(x, y) = \mathcal{X}_k(x) \mathcal{Y}_l(y), \label{eq:nstestfunction} \tag{1} \end{equation} $$

and solutions (trial functions) as

$$ \begin{equation} u(x, y) = \sum_{k}\sum_{l} \hat{u}_{kl} v_{kl}(x, y). \label{eq:nstrialfunction} \tag{2} \end{equation} $$

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):

$$ \begin{equation} \mathcal{X}_k(x) = L_k(x) - L_{k+2}(x), \quad \forall \, k \in \boldsymbol{k}^{N_0-2}, \label{eq:D0} \tag{3} \end{equation} $$

$$ \begin{equation} \mathcal{Y}_l(y) = L_l(y) - L_{l+2}(y), \quad \forall \, l \in \boldsymbol{l}^{N_1-2}, \label{eq:D1} \tag{4} \end{equation} $$

where $\boldsymbol{k}^{N_0-2} = (0, 1, \ldots, N_0-3)$, $\boldsymbol{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

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

$$ \begin{equation} \mathcal{Y}_l(y) = L_l(y) - L_{l+2}(y), \quad \forall \, l \in \boldsymbol{l}^{N_1-2}. \label{_auto1} \tag{5} \end{equation} $$

$$ \begin{equation} \mathcal{Y}_{N_1-2}(y) = (L_0+L_1)/2 \quad \left(=(1+y)/2\right), \label{_auto2} \tag{6} \end{equation} $$

$$ \begin{equation} \mathcal{Y}_{N_1-1}(y) = (L_0-L_1)/2 \quad \left(=(1-y)/2\right). \label{_auto3} \tag{7} \end{equation} $$

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

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

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 \boldsymbol{k}^{N_0} \times \boldsymbol{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 \boldsymbol{k}^{N_0-2} \times \boldsymbol{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

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

$$ \begin{equation} V_{1}^{\boldsymbol{N}}(\boldsymbol{x}) = D_0^{N_0}(x) \otimes D_1^{N_1}(y), \label{_auto4} \tag{8} \end{equation} $$

$$ \begin{equation} V_{0}^{\boldsymbol{N}}(\boldsymbol{x}) = D_0^{N_0}(x) \otimes D_0^{N_1}(y), \label{_auto5} \tag{9} \end{equation} $$

$$ \begin{equation} P^{\boldsymbol{N}}(\boldsymbol{x}) = P^{N_0}(x) \otimes P^{N_1}(y). \label{_auto6} \tag{10} \end{equation} $$

With shenfun the tensor product spaces are created as

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

Note that the second vector basis, $W_0^{\boldsymbol{N}} = V_0^{\boldsymbol{N}} \times V_0^{\boldsymbol{N}}$, uses homogeneous boundary conditions throughout.

Mixed variational form

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

$$ \begin{equation} \int_{\Omega} (\nu \nabla^2 \boldsymbol{u} - \nabla p ) \cdot \boldsymbol{v} \, dxdy = \int_{\Omega} (\nabla \cdot \boldsymbol{u}\boldsymbol{u}) \cdot \boldsymbol{v}\, dxdy \quad\forall \boldsymbol{v} \, \in \, W_0^{\boldsymbol{N}}, \label{eq:nsvarform} \tag{11} \end{equation} $$

$$ \begin{equation} \int_{\Omega} \nabla \cdot \boldsymbol{u} \, q \, dxdy = 0 \quad\forall q \, \in \, P^{\boldsymbol{N}}. \label{_auto7} \tag{12} \end{equation} $$

Note that we are using test functions $\boldsymbol{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 $\boldsymbol{u}^k$ and compute the nonlinearity using only known solutions $\int_{\Omega} (\nabla \cdot \boldsymbol{u}^k\boldsymbol{u}^k) \cdot \boldsymbol{v}\, dxdy$. Using further integration by parts we end up with the equations to solve for iteration number $k+1$ (using $\boldsymbol{u} = \boldsymbol{u}^{k+1}$ and $p=p^{k+1}$ for simplicity)

$$ \begin{equation} -\int_{\Omega} \nu \nabla \boldsymbol{u} \, \colon \nabla \boldsymbol{v} \, dxdy + \int_{\Omega} p \nabla \cdot \boldsymbol{v} \, dxdy = \int_{\Omega} (\nabla \cdot \boldsymbol{u}^k\boldsymbol{u}^k) \cdot \boldsymbol{v}\, dxdy \quad\forall \boldsymbol{v} \, \in \, W_0^{\boldsymbol{N}}, \label{eq:nsvarform2} \tag{13} \end{equation} $$

$$ \begin{equation} \int_{\Omega} \nabla \cdot \boldsymbol{u} \, q \, dxdy = 0 \quad\forall q \, \in \, P^{\boldsymbol{N}}. \label{_auto8} \tag{14} \end{equation} $$

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

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

$$ \begin{equation} a((\boldsymbol{u}, p), (\boldsymbol{v}, q)) = L((\boldsymbol{v}, q)) \quad \forall (\boldsymbol{v}, q) \in W_0^{\boldsymbol{N}} \times P^{\boldsymbol{N}}, \label{_auto9} \tag{15} \end{equation} $$

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

$$ \begin{equation} a((\boldsymbol{u}, p), (\boldsymbol{v}, q)) = -\int_{\Omega} \nu \nabla \boldsymbol{u} \, \colon \nabla \boldsymbol{v} \, dxdy + \int_{\Omega} p \nabla \cdot \boldsymbol{v} \, dxdy + \int_{\Omega} \nabla \cdot \boldsymbol{u} \, q \, dxdy, \label{_auto10} \tag{16} \end{equation} $$

$$ \begin{equation} L((\boldsymbol{v}, q); \boldsymbol{u}^{k}) = \int_{\Omega} (\nabla \cdot \boldsymbol{u}^{k}\boldsymbol{u}^{k}) \cdot \boldsymbol{v}\, dxdy. \label{_auto11} \tag{17} \end{equation} $$

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 $\boldsymbol{u}^0 = (0, 0)$

• while not converged:

  • assemble $L((\boldsymbol{v}, q); \boldsymbol{u}^{k})$

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

  • compute error = $\int_{\Omega} (\boldsymbol{u}^{k+1}-\boldsymbol{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^{\boldsymbol{N}} \times P^{\boldsymbol{N}}$ that couples velocity and pressure

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

Notice.

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:

Notice.

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

$$ \begin{equation} \label{eq:partialeq1} \tag{18} \int_{\Omega}-\nu \left(\frac{\partial \boldsymbol{v}[0]}{\partial x}, \frac{\partial \boldsymbol{v}[0]}{\partial y}\right) \cdot \left(\frac{\partial \boldsymbol{u}[0]}{\partial x}, \frac{\partial \boldsymbol{u}[0]}{\partial y}\right) dx dy . \end{equation} $$

We can now insert for test function $\boldsymbol{v}[0]$

$$ \begin{equation} \boldsymbol{v}[0]_{kl} = \mathcal{X}_k \mathcal{Y}_l, \quad (k, l) \in \boldsymbol{k}^{N_0-2} \times \boldsymbol{l}^{N_1-2} \label{_auto12} \tag{19} \end{equation} $$

and trialfunction

$$ \begin{equation} \boldsymbol{u}[0]_{mn} = \sum_{m=0}^{N_0-3} \sum_{n=0}^{N_1-1} \hat{\boldsymbol{u}}[0]_{mn} \mathcal{X}_m \mathcal{Y}_n, \label{_auto13} \tag{20} \end{equation} $$

where $\hat{\boldsymbol{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 (18), we obtain after a few trivial manipulations

$$ \begin{equation} -\sum_{m=0}^{N_0-3} \sum_{n=0}^{N_1-1} \nu \Big( \underbrace{\int_{-1}^{1} \frac{\partial \mathcal{X}_k(x)}{\partial x} \frac{\partial \mathcal{X}_m}{\partial x} dx \int_{-1}^{1} \mathcal{Y}_l \mathcal{Y}_n dy}_{A[0]} + \underbrace{\int_{-1}^{1} \mathcal{X}_k(x) X_m(x) dx \int_{-1}^{1} \frac{\partial \mathcal{Y}_l}{\partial y} \frac{\partial \mathcal{Y}_n}{\partial y} dy}_{A[1]} \Big) \hat{\boldsymbol{u}}[0]_{mn}. \label{_auto14} \tag{21} \end{equation} $$

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

$$ \begin{equation} \underbrace{\int_{-1}^{1} \frac{\partial \mathcal{X}_k(x)}{\partial x} \frac{\partial \mathcal{X}_m}{\partial x} dx}_{c_{km}} \underbrace{\int_{-1}^{1} \mathcal{Y}_l \mathcal{Y}_n dy}_{f_{ln}} \label{_auto15} \tag{22} \end{equation} $$

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{\boldsymbol{u}}[0]_{m,N_1-2}$ and $\hat{\boldsymbol{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

$$ \sum_{m=0}^{N_0-3}\sum_{n=0}^{N_1-1} c_{km}f_{ln} \hat{\boldsymbol{u}}[0]_{m, n} = \sum_{m=0}^{N_0-3}\sum_{n=0}^{N_1-3}c_{km}f_{ln}\hat{\boldsymbol{u}}[0]_{m, n} + \sum_{m=0}^{N_0-3}\sum_{n=N_1-2}^{N_1-1}c_{km}f_{ln}\hat{\boldsymbol{u}}[0]_{m, n}, \quad \forall (k, l) \in \boldsymbol{k}^{N_0-2} \times \boldsymbol{l}^{N_1-2}, $$

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

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

$$ M = \begin{bmatrix} A[0]+A[1] & 0 & G[0] \\ \label{eq:nsbmatrix} \tag{23} 0 & A[2]+A[3] & G[1] \\ D[0] & D[1] & 0 \end{bmatrix} $$

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

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

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 appended set_boundary_dofs() ensures the known boundary values of the solution are fixed for ui_hat and ui_new.

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 $\boldsymbol{u}\boldsymbol{u}$. The required arrays and spaces are created as

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

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:

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

$$ \begin{align*} \hat{\boldsymbol{u}}^{k+1} &= \alpha \hat{\boldsymbol{u}}^* + (1-\alpha)\hat{\boldsymbol{u}}^{k},\\ \hat{p}^{k+1} &= \alpha \hat{p}^* + (1-\alpha)\hat{p}^{k}, \end{align*} $$

where $\hat{\boldsymbol{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

The last three lines plots the velocity vectors that are shown in Figure. 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. However, choosing instead the regularized lid, the solutions will be nice and smooth, both for Legendre and Chebyshev polynomials.

Velocity vectors for Re=100 using Chebyshev.

Complete solver

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