Correct implemtation of streamfunction in lid driven cavity problem from fluid mechanics

[Jul3k]

I am attempting to solve the problem of the lid-driven cavity in fluid mechanics using the streamfunction formulation:

$$\frac{\partial}{\partial t} (\nabla^2 \psi) + \frac{\partial \psi}{\partial y} \frac{\partial}{\partial x} (\nabla^2 \psi) - \frac{\partial \psi}{\partial x} \frac{\partial}{\partial y} (\nabla^2 \psi) = \nu \nabla^4 \psi$$

Given the boundary conditions:

$$\psi = 0 \text{ on } \partial \Omega,$$ $$\frac{\partial \psi}{\partial n} = \begin{cases} 1 & \text{if } x_2 = 1, \\ 0 & \text{if } 0 \le x_2 < 1. \end{cases}$$

The weak formulation I used for the steady-state solution is based on a Fenics implementation:

$$\nu \int_{\Omega} \Delta \psi \Delta \phi + \int_{\Omega} \Delta \psi \left( \frac{\partial \psi}{\partial y} \frac{\partial \phi}{\partial x} - \frac{\partial \psi}{\partial x} \frac{\partial \phi}{\partial y} \right) = \int_{\Omega} f \phi$$

where $\nu$ where is the kinematic viscosity, $\psi$ the stream function and $\phi$ the test function.

For the first part, I used the biharmonic implementation from example 20: u*ddot(dd(u),dd(v)), which provides the solution for the Stokes flow for low Reynold numbers and matches my expectations.

My problem lies in implementing the second part of the BilinearForm, which should be div(grad(psi))*(grad(psi)[1]*grad(phi)[0]-grad(psi)[0]*grad(phi)[1]) according to the Fenics implementation. My approach was:

trace(dd(u)) * (grad(u)[1]*grad(v)[0] - grad(u)[0]*grad(v)[1])

However, grad(u)[1]*grad(v)[0] - grad(u)[0]*grad(v)[1] returns a zero matrix, rendering the second part ineffective.

What is the correct implementation of the second part of the BilinearForm?

Reference solutions for Re=10 and Re=800 can be found here for an implementation using finite differences and OpenFoam. For Re=800 I expect:

Below is a small test implementation:

from skfem import *
from skfem.helpers import *
from skfem.visuals.matplotlib import plot as fem_plot

Re = 800
nu = 1/Re

mesh = MeshTri().refined(6)
elm  = ElementTriMorley()
V    = Basis(mesh, elm)

@BilinearForm
def stream_function(u, v, w):
    # Biharmonic + 
    return nu*ddot(dd(u),dd(v)) + trace(dd(u))*(d(u)[1]*d(v)[0]-d(u)[0]*d(v)[1]) 
    
x = V.zeros()
x[V.get_dofs('top').facet['u_n']] = 1 # Derivative on boundary defines the velocity

A = asm(stream_function, V)
b = V.zeros()

u = solve(*condense(A, b, x=x, D=V.get_dofs()))

fem_plot(V,u, shading='gouraud', colorbar=True)

[Jul3k]

Sorry, I just realized I should have asked this as a discussion.

[kinnala]

Before looking into the question, I made a sanity check and read the code. It seems that your bilinear form is actually not bilinear, it is nonlinear in $u$.

Probably that is intended since you are studying high Re case, I did not look into the formulation in detail. However, BilinearForm does only support bilinear forms.

To solve nonlinear problems you usually implement Newton’s method in a for loop and, with the help of the derivative of the nonlinear weak form, you write a bilinear form (so-called Jacobian) which describes a single Newton iteration.

[Jul3k]

Thank you. This makes sense and I see that now. The bilinear form is actually not linear. I will try to adapt ex10 to my problem.

[Jul3k]

As the solution of the of the nonlinear streamfunction requires a iterative approach anyway , I was looking into the following pure streamfunction "one-level time-stepping schemes of IMplicit-EXplicit (IMEX) type. For IMEX schemes the convective term is treated explicitly, while the diffusive term is diagonally implicit". [A High Order Compact Scheme for the Pure-Streamfunction Formulation
of the Navier-Stokes Equations, Ben-Artzi et al., Chapter 4]. Do you think that such an approach could be implemented in SKFEM?

Step 1: Computation of $\psi^{n+1/2}$

$$ \left( \tilde{\Delta} - \frac{1}{4} k \nu \tilde{\Delta}^2 \right) \psi^{n+1/2} = \left( \tilde{\Delta} \psi^n + \frac{1}{4} k \nu \tilde{\Delta}^2 \right) \psi^n + \frac{1}{2} k \tilde{C} (\psi^n) \tag{5.25} $$

Step 2: Computation of $\psi^{n+1}$

$$ \left( \tilde{\Delta} - \frac{1}{2} k \nu \tilde{\Delta}^2 \right) \psi^{n+1} = \left( \tilde{\Delta} \psi^n + \frac{1}{2} k \nu \tilde{\Delta}^2 \right) \psi^n + k \tilde{C} (\psi^{n+1/2}) \tag{5.26} $$

where $k$ is the time step size and $\tilde{C}$ represents the convective term.

Below is a first attempt at the implementation. The convective term is commented out at the moment as I struggle how to compute the non-linear explicit form. The time stepped solution does not converge to the result that I obatin by only solving the strokes flow as described in the above example, so there must be a mistake in the implementation. I highly appreciate any assistance as I am lacking the mathematical background.

from skfem import *
from skfem.models.poisson import laplace

# Constants
dt = 0.01  # Time step size
Re = 100
nu = 1 / Re

# Mesh and basis functions setup
mesh = MeshTri().refined(5)
element = ElementTriMorley()
V = Basis(mesh, element)

# Define bilinear forms for U, D, C
@BilinearForm
def biharmonic(u, v, w):
    return nu * ddot(dd(u), dd(v))

@BilinearForm
def convective(u, v, w):
    # TODO: It is non-linear so it cannot be implemented as a LinearForm
    return dot(d(u),d(v)) *(d(u)[1]-d(u)[0]) 
    #return trace(dd(u))*(d(u)[1]*d(v)[0]-d(u)[0]*d(v)[1])

# Assemble matrices
U = asm(laplace, V)     # Laplacian matrix for diffusive term
D = asm(biharmonic, V)  # Biharmonic matrix for second diffusive contribution
C = asm(convective, V)  # Convective term (nonlinear)

# Initial condition
phi = V.zeros()
x = V.zeros()
x[V.get_dofs('top').facet['u_n']] = 1 # Derivative on boundary defines the velocity

# Time-stepping
for step in range(10):  # Example: 10 time steps
    # Prediction Step TODO: Plus convective term
    #b2 = (U @ phi) - dt/2 * (C @ phi) + dt/4 * (D @ phi)
    b2 = (U @ phi) + dt/4 * (D @ phi)
    b2[V.get_dofs().nodal['u']] = 0

    phi_p = solve(*condense(U - dt/4 * D, b2, x=x, D=V.get_dofs()))

    # Correction Step TODO: Plus convective term
    b3 = (U @ phi) + dt/2 *(D @ phi)

    phi = solve(*condense(U - dt/2 * D, b3, x=x, D=V.get_dofs()))

    grad_phi = V.interpolate(phi).grad
    u = V.project(grad_phi[1])
    v = V.project(-grad_phi[0])

    du = V.project(V.interpolate(u).grad[1])
    dv = V.project(V.interpolate(v).grad[0])

    w = dv - du

    figure(figsize=(9,7))
    ax = subplot(221)
    title('Ψ'); axis('equal')
    fem_plot(V,phi, ax=ax, shading='gouraud', colorbar=True)
    ax = subplot(222)
    title('w'); axis('equal')
    fem_plot(V,w, ax=ax, shading='gouraud', colorbar=True)
    ax = subplot(223)
    title('u'); axis('equal')
    fem_plot(V,u, ax=ax, shading='gouraud', colorbar=True)
    ax = subplot(224); title('v'); axis('equal')
    fem_plot(V,v, ax=ax, shading='gouraud', colorbar=True)

[Jul3k]

Just by experimenting I discovered that by changeing the sign in :

@BilinearForm
def biharmonic(u, v, w):
    return -nu * ddot(dd(u), dd(v))

The solution for the stokes flow converges to the expected value. So the only problem is how to implement the explicit nonlinear convective part

[Jul3k]

I was able to solve the Navier-Stokes equations in pure streamfunction formulation using the following scheme to time step:

$$\frac{\Delta \psi}{dt} + \frac{\psi}{dy} \cdot \frac{\Delta \psi}{dx} - \frac{\psi}{dx} \cdot \frac{\Delta \psi}{dx} - \nu \cdot \Delta^2 \cdot \psi = 0$$ $$\frac{\Delta \psi^-}{k} + \frac{\psi^*}{dy} \cdot \frac{\Delta \psi^+}{dx} - \frac{\psi^*}{dx} \cdot \frac{\Delta \psi^+}{dx} - \nu \cdot \Delta^2 \cdot \psi^+ = 0$$

with $\psi^+=\frac{\psi^{n+1} + \psi^{n}}{2}$, $\psi^-=\psi^{n+1} - \psi^{n}$ and $\psi^*=\psi^n$. Rearanged to $\psi^{n+1}$ this gives:

$$\psi^{n+1} \left( \underbrace{\Delta}_{\text{Laplace [L]}} + \underbrace{\frac{k\nu}{2} \Delta^2}_{\text{Biharmonic [H]}} + \frac{k}{2} \underbrace{\left( \frac{\psi^*}{dy} \frac{\Delta}{dx} - \frac{\psi^*}{dx} \frac{\Delta}{dy} \right)}_{\text{Convective [C]}} \right) = \psi^{n} \left( \underbrace{\Delta}_{\text{Laplace [L]}} - \underbrace{\frac{k\nu}{2} \Delta^2}_{\text{Biharmonic [H]}} - \frac{k}{2} \underbrace{\left( \frac{\psi^*}{dy} \frac{\Delta}{dx} - \frac{\psi^*}{dx} \frac{\Delta}{dy} \right)}_{\text{Convective [C]}} \right)$$

To linearize the convective term I ended up with the following weak formulation:

$$C := \frac{\psi^*}{dy} \left(\frac{\phi}{dx}\frac{\psi}{dx^2} + \frac{\phi}{dy}\frac{\psi}{dxdy} \right) - \frac{\psi^*}{dx} \left(\frac{\phi}{dy}\frac{\psi}{dy^2} + \frac{\phi}{dx}\frac{\psi}{dxdy} \right)$$

The essential parts of the implementation are:

@BilinearForm
def conv(u, v, w):
    return d(w.ppsi)[1]*(d(v)[0]*dd(u)[0,0] + d(v)[1]*dd(u)[0,1]) - d(w.ppsi)[0] * (d(v)[1]*dd(u)[1,1] + d(v)[0]*dd(u)[0,1])

L = asm(laplace, V)
H = asm(biharmonic, V)

for i in range(N):
    C = asm(conv,V, ppsi=psi)

    A = L + H * dt/2 + C*dt/2
    B = L - H * dt/2 - C*dt/2

    b = B @ psi
    psi = solve(*condense(A,b, x=x, D=V.get_dofs()))

@kinnala Do you see a way how this computation can be speed up? Maybe the convective term can be split into two bilinear forms and the matrix prefactorized similar as in example 19? Would this be a nice example to add to the repository?

The simulation runs fine an yiels the solutions I would expect for the lid driven cavity problem and also for a turbulent flow forming a Kármán vortex street:

Below is again a small example code for the lid driven cavity:

from skfem import *
from skfem.helpers import *
from skfem.visuals.matplotlib import plot as fem_plot
from skfem.models.poisson import laplace

# Constants
Re = 400
nu = 1 / Re

# Mesh and basis functions setup
mesh = MeshTri().refined(6)
element = ElementTriMorley()
V = Basis(mesh, element)

# Define bilinear forms for U, D, C
@BilinearForm
def biharmonic(u, v, w):
    return nu * ddot(dd(u), dd(v))

@BilinearForm
def conv(u, v, w):
    return d(w.ppsi)[1]*(d(v)[0]*dd(u)[0,0] + d(v)[1]*dd(u)[0,1]) - d(w.ppsi)[0] * (d(v)[1]*dd(u)[1,1] + d(v)[0]*dd(u)[0,1])


# Assemble matrices
L = asm(laplace, V)
H = asm(biharmonic, V)

# Initial condition
psi = V.zeros()
ppsi = V.zeros()
x = V.zeros()
x[V.get_dofs('top').facet['u_n']] = 1 # Derivative on boundary defines the velocity

dt = 0.1  # Time step size
T = 400
t = 0
i = 0
# Time-stepping


while (t<T):  # Example: 10 time steps  
    i += 1; t += dt

    C = asm(conv,V, ppsi=psi)# + dt/2*(psi-ppsi))
    A = L + H * dt/2 + C*dt/2
    B = L - H * dt/2 - C*dt/2

    b = B @ psi
    ppsi = psi
    psi = solve(*condense(A,b, x=x, D=V.get_dofs()))
    
    
    psii = V.interpolate(psi)
    u = V.project(psii.grad[1])
    v = V.project(-psii.grad[0])
    
    du = V.project(V.interpolate(u).grad[1])
    dv = V.project(V.interpolate(v).grad[0])

    w = dv - du
    
    res = np.linalg.norm(psi - ppsi)
    print(f"{i:3d} | {dt:4.3f} | {t:4.1f} | {res:4.3f}")
    
    dt = dt/res
    
    if i%4 != 1:
        continue
    
    figure(figsize=(9,7))
    ax = subplot(221); title('Ψ'); axis('equal')
    fem_plot(V,psi, ax=ax, shading='gouraud', levels=13, colorbar=True)
    ax = subplot(222); title('velo'); axis('equal')
    fem_plot(V,sqrt(u**2+v**2) , ax=ax, shading='gouraud', colorbar=True)
    ax = subplot(223); title('u'); axis('equal')
    fem_plot(V, u, ax=ax, shading='gouraud', colorbar=True)
    ax = subplot(224); title('v'); axis('equal')
    fem_plot(V, v, ax=ax, shading='gouraud', colorbar=True)
    show()

[kinnala]

I cannot see an obvious way to prefactorize the matrix A because it is different for each time step. Using iterative matrix solvers may sometimes lead to improvements in run time because it allows using the previous timestep result as an initial guess. However, designing proper preconditioners might end up being difficult for this class of problems.

I am happy to merge any examples that can be run in few seconds. All examples are run by the CI test suite to avoid breaking older examples and we want it to run quickly.

[Jul3k]

I thought about solving the convective term explicit. Then the left hand side does not change and it could be preconditioned.

A = L + H * dt/2          #LHS
B = L - H * dt/2 - C*dt   #RHS

I am however currently struggeling to implement the preconditioning based on the provided examples as they are already quite complex. I think my problem is how to apply the boundary conditions to the system matrix prior to the preconditioning.