Navier-Stokes equations with a couple velocity-pressure scheme

General
1

Hello,
My advisor and I have been working on solving the Navier-Stokes equations with a couple velocity-pressure scheme. We are using FEniCSx0.8.
We implemented a manufactured solution test to check our implementation using either an explicit scheme (where nonlinear term (u.grad)u is treated explicitly) or a semi-implicit scheme. Both methods uses a Backward Euler method.
While the explicit scheme works as intended (order 1 convergence), the semi-implicit scheme does not converge in pressure (the error always stays around 0.6).
We introduced a parameter “explicit” that can be set to either 1 to the run explicit scheme or 0 to run the semi-implicit scheme so the only difference is highlighted in this part of the code when defining the bilinear and linear form of the problem.

# ############### Load/import modules
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
from dolfinx import fem, io, mesh, plot, default_real_type
from dolfinx.fem import (Constant, Function, functionspace, assemble_scalar, 
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block, LinearProblem, create_vector
from basix.ufl import element
from ufl import (FacetNormal, Identity, TestFunction, TrialFunction, cos, sin, variable, as_vector,
                 div, dot, ds, dx, inner, lhs, nabla_grad, grad, rhs, sym, diff, Measure, gt, outer,
                 SpatialCoordinate, exp, avg, conditional)
# ############### END Load/import modules


# ############### Problem Parameters
# Time interval
t = 0.0
T = 1.0
# Time discretisation
explicit=0
num_steps = 20
dt = T / num_steps
# Space discretisation
n_elem = 10
# Degree Lagrange polynomial
deg = 1 #P_deg for pressure and P_{deg+1} for velocity
# Problem parameters
nu = 0.1
# ############### END Problem Parameters


# ############### Mesh
# Create mesh
msh = mesh.create_rectangle(
    comm=MPI.COMM_WORLD,
    points=((0.0, 0.0), (1.0, 1.0)),
    n=(n_elem,n_elem),
    cell_type=mesh.CellType.triangle,
)
gdim = msh.geometry.dim
# Define spatial coordinate (x[0], x[1], etc.)
x = SpatialCoordinate(msh)
# ############### END Mesh


# ############### Problem functions
#Source term in Navier-Stokes
class f_exp:
    def __init__(self, t):
        self.t = t
    def __call__(self, x):
        values = np.zeros((gdim, x.shape[1]), dtype=PETSc.ScalarType)
        values[0] = (
            #time derivative
            -np.sin(np.pi * x[0]) * np.cos(np.pi * x[1]) * np.sin(self.t)
            #diffusion
            + nu * 2.0 * np.pi**2 * np.sin(np.pi * x[0]) * np.cos(np.pi * x[1]) * np.cos(self.t)
            #pressure gradient
            - np.pi * np.sin(np.pi * x[0]) * np.cos(np.pi * x[1]) * np.cos(self.t)
            #nonlinear term
            + np.pi * np.sin(np.pi * x[0]) * np.cos(np.pi * x[0]) * np.cos(self.t)**2
        )
        values[1] = (
            #time derivative
            np.cos(np.pi * x[0]) * np.sin(np.pi * x[1]) * np.sin(self.t)
            #diffusion
            - nu * 2.0 * np.pi**2 * np.cos(np.pi * x[0]) * np.sin(np.pi * x[1]) * np.cos(self.t)
            #pressure gradient
            - np.pi * np.cos(np.pi * x[0]) * np.sin(np.pi * x[1]) * np.cos(self.t)
            #nonlinear term
            + np.pi * np.sin(np.pi * x[1]) * np.cos(np.pi * x[1]) * np.cos(self.t)**2
        )
        return values

# Exact velocity (also used for boundary condition(
class u_exact_exp():
    def __init__(self, t):
        self.t = t
    def __call__(self, x):
        values = np.zeros((gdim, x.shape[1]), dtype=PETSc.ScalarType)
        values[0] = np.sin(np.pi * x[0]) * np.cos(np.pi * x[1])*np.cos(self.t)
        values[1] = -np.cos(np.pi * x[0]) * np.sin(np.pi * x[1])*np.cos(self.t)
        return values

# Exact pressure
class p_exact_exp():
    def __init__(self, t):
        self.t = t
    def __call__(self, x):
        values = np.cos(np.pi*x[0]) * np.cos(np.pi*x[1])*np.cos(self.t)
        return values

# Expression for location Dirichlet BC
def gamma_D_exp(x):
        return (np.isclose(x[0],0.0) | np.isclose(x[0],1.0) |
                np.isclose(x[1],0.0) | np.isclose(x[1],1.0)
                )
# ############### END Problem functions


# ############### Useful functions to compute norms, average, etc
#Compute the average of a function  over the domain
def domain_average(msh, p):
    """Compute the average of a function over the domain"""
    vol = msh.comm.allreduce(
        fem.assemble_scalar(fem.form(fem.Constant(msh, default_real_type(1.0)) * dx)), op=MPI.SUM
    )
    return (1.0 / vol) * msh.comm.allreduce(fem.assemble_scalar(fem.form(p * dx)), op=MPI.SUM)

# Compute norm L2 
def norm_L2(comm, v):
    """Compute the L2(Ω)-norm of v"""
    return np.sqrt(msh.comm.allreduce(fem.assemble_scalar(fem.form(inner(v, v) * dx)), op=MPI.SUM))

# Compute semi norm H1
def norm_sH1(comm, v):
    """Compute the semi H1(Ω)-norm of v"""
    return np.sqrt(msh.comm.allreduce(fem.assemble_scalar(
        fem.form(inner(grad(v), grad(v)) * dx)), op=MPI.SUM))
# ############### END Useful functions to compute norms, average, etc


# ############################################################################
# ################# VARIATIONAL FORMULATION
# ############################################################################

# ############### Functional Space
# Create Taylor-Hood element
Pu = element("Lagrange", msh.basix_cell(), deg+1, shape=(msh.geometry.dim,))
Pp = element("Lagrange", msh.basix_cell(), deg)

# Create variational function spaces
V = functionspace(msh,Pu)
Q = functionspace(msh,Pp)
# ############### END Functional Space


# ############### Boundary Conditions
facets = mesh.locate_entities_boundary(msh, 1, gamma_D_exp)
dofs_u = fem.locate_dofs_topological(V,msh.topology.dim - 1, facets)
# Define time dependent boundary conditions
u_bdy = Function(V)
u_exact_exp.t = t
u_exact_exp = u_exact_exp(t)
u_bdy.interpolate(u_exact_exp)
bc_u = fem.dirichletbc(u_bdy, dofs_u)
bcs=[bc_u]
# ############### END Boundary Conditions


# ############### Define variational problem
# Define trial and test functions
(u, p) = TrialFunction(V), TrialFunction(Q)
(v, q) = TestFunction(V), TestFunction(Q)

# Define approximation at previous time step
u_n = Function(V)
u_exact_exp.t = t
u_n.interpolate(u_exact_exp)

# Define/Declare source term
#f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))
f = Function(V)
f_exp.t = t
f_exp = f_exp(t)
f.interpolate(f_exp)

# Define bilinear form of variational formulation
if explicit==0:
    a_der_LHS = inner(u / dt + dot(u_n, nabla_grad(u)), v) * dx
else:
    a_der_LHS = inner(u / dt, v) * dx

a_00 = nu * inner(grad(u), grad(v)) * dx  
a_01 = -inner(p, div(v)) * dx
a_10 = -inner(div(u), q) * dx
a = form([[a_00 + a_der_LHS, a_01], [a_10, None]])

# Define linear form of variational formulation
if explicit==0:
    a_der_RHS = inner(u_n / dt, v) * dx
else:
    a_der_RHS = inner(u_n / dt - dot(u_n, nabla_grad(u_n)), v) * dx

L_0 = inner(f, v) * dx 
L_1 = inner(fem.Constant(msh, default_real_type(0.0)), q) * dx
L = form([L_0 + a_der_RHS, L_1])
# ############### END Define variational problem


# ############### Assemble matrix-vector and Define solver
# Assemble Navier-Stokes problem
A = assemble_matrix_block(a, bcs=bcs)
A.assemble()
b = assemble_vector_block(L, a, bcs=bcs)

# Create vector to store solution linear system
x = A.createVecRight()

# Create and configure solver
ksp = PETSc.KSP().create(msh.comm)  # type: ignore
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
opts = PETSc.Options()  # type: ignore
opts["mat_mumps_icntl_14"] = 80  # Increase MUMPS working memory
opts["mat_mumps_icntl_24"] = 1  # Option to support solving a singular matrix (pressure nullspace)
opts["mat_mumps_icntl_25"] = 0  # Option to support solving a singular matrix (pressure nullspace)
opts["ksp_error_if_not_converged"] = 1
ksp.setFromOptions()
# ############### END Assemble matrix-vector and Define solver


# #############################################################################
# ######################## TIME LOOP
# #############################################################################
for i in range(num_steps):
    # Update current time step
    t += dt

    # Compute velocity Dirichlet boundary condition
    u_exact_exp.t = t
    u_bdy.interpolate(u_exact_exp)
    bc_u = fem.dirichletbc(u_bdy, dofs_u)
    bcs=[bc_u]

    # Compute source term
    f_exp.t=t
    f.interpolate(f_exp)

    # Assemble/Compute LHS (bcs time dependent) and RHS
    A = assemble_matrix_block(a, bcs=bcs)
    A.assemble()
    b = assemble_vector_block(L, a, bcs=bcs)

    # Solve the variational problem
    try:
        ksp.solve(b, x)
    except PETSc.Error as e:  # type: ignore
        if e.ierr == 92:
            print("The required PETSc solver/preconditioner is not available. Exiting.")
            print(e)
            exit(0)
        else:
            raise e

    # Split the solution
    u_h = fem.Function(V)
    p_h = fem.Function(Q)
    p_h.name = "p"
    offset = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
    u_h.x.array[:offset] = x.array_r[:offset]
    u_h.x.scatter_forward()
    p_h.x.array[: (len(x.array_r) - offset)] = x.array_r[offset:]
    p_h.x.scatter_forward()

    # Subtract the average of the pressure
    p_h.x.array[:] -= domain_average(msh, p_h)

    # Update u_n
    u_n.x.array[:] = u_h.x.array


# #############################################################################
# ######################## POST-PROCESSING (FINAL ITERATION)
# #############################################################################
# ############## COMPUTE L2 and H1 ERRORS  ##########
# Function spaces for exact velocity and pressure
V_e = fem.functionspace(msh, ("Lagrange", deg + 1 + 3, (gdim,)))
Q_e = fem.functionspace(msh, ("Lagrange", deg + 3))
# Interpolate exact velocity
u_e = fem.Function(V_e)
u_exact_exp.t = t
u_e.interpolate(u_exact_exp)
# Interpolate exact pressure
p_e = fem.Function(Q_e)
p_exact_exp.t = t
p_exact_exp = p_exact_exp(t)
p_e.interpolate(p_exact_exp)

# Compute errors
e_u = norm_L2(msh.comm, u_h - u_e)
e_sH1_u = norm_sH1(msh.comm, u_h - u_e)
e_div_u = norm_L2(msh.comm, div(u_h))
p_e_avg = domain_average(msh, p_e)
e_p = norm_L2(msh.comm, p_h - (p_e - p_e_avg))

if  i == num_steps - 1 :
# Print errors
    if msh.comm.rank == 0:
        print(f"e_u = {e_u}")
        print(f"e_sH1_u = {e_sH1_u}")
        print(f"e_div_u = {e_div_u}")
        print(f"e_p = {e_p}")
# ############## END COMPUTE L2 and H1 ERRORS  ##########

I am not sure if we are doing something wrong in the above formulation when switching between explicit/semi-implicit and was wondering if you would see something wrong there.

Thank you for your help.

2

That’s unlike any semi implicit scheme I have seen for NS. Say that your solution u is stationary, ie constant in time. Here you have an extra term on the RHS which makes your scheme inconsistent.

4

The notation may be confusing but the semi-implicit scheme is used when we set “explicit=0”. So, in the part of the code you highlighted, the RHS corresponds to what is used for the explicit scheme (i.e. (u grad)u is treated explicitly and is not present in the LHS).

The semi implicit scheme is run when we set “explicit=0”. In that case the LHS has a term of the form

( u /dt + (un grad) u , v)

while the right hand-side only has the rest of the time derivative

(un / dt, v)

It seems to me that there is no extra term when we use the semi implicit scheme (i.e. “explicit=0”).

5

Well, you see why in Read before posting: How do I get my question answered? we stress that the provided example must be minimal.

Asking someone to go through 300 lines of code in their own free time (because we are all volunteers here) is no small job. I’ll try to have another look at some point this week when I have time, but you may want to try to simplify your example as much as possible to help anyone else reading your message.

6

Thank you for your reply. I do understand that the code may be too long and that it requires effort to go through, but the problem in the code arises with time dependence, so I am not sure how much I can simplify it in such a way that it becomes possible to identify the issue.
Thank you so much anyway for your time and for trying to answer my question.