# Confirming Error Bounds with FENICS, pressure not converging

Hey guys,

I implemented the Taylor-Green-Vortex with time-dependent boundary conditions for the velocity and pressure to be solved with the Chorin-Projection method as used in one of the FENICS tutorials. Since we know the analytical solution of the Taylor-Green-Vortex I now tried to confirm the Error bounds given by Andreas Prohl in “Projection and Quasi-Compressibility Methods for Solving the Incompressible Navier-Stokes Equations”.

Note that by the standard norm the spatial L2 norm over the domain is meant.

The velocity error is converging well and just as expected, but the pressure error stays weirdly high independent of my spatial or temporal discretization. In fact my pressure error does not converge at all!

I used Taylor-Hood-Elements(P2-P1) so the LBB Condition should be alright and everything should be stable enough.

Why is the pressure error not converging as expected although IC, BC etc. should be sufficiently smooth? Any ideas? I dont have any more ideas . I tried setting the degree_rise of the errornorm to 0 but this was neither the reason for this weird behaviour…

Best regards!
Max

My code is below

``````from dolfin import *
import numpy as np
from math import exp

# --- Set parameter values
# - end time
T = 1
# - time discretization steps
dt = 0.01
# - spatial discretization steps
mesh_part = 65
# - Reynolds number
Re = 1

# --------------------------------------------
# --- Create MESH: Square with length 2 pi ---
p1 = Point(0,0)
p2 = Point(2*DOLFIN_PI,2*DOLFIN_PI)
mesh = RectangleMesh(p1,p2,mesh_part,mesh_part)

# ------------------------------------------
# --- initialize VARIATIONAL FORMULATION ---
# Define function spaces (P2-P1)
V = VectorFunctionSpace(mesh, "P", 2, dim=2)
Q = FunctionSpace(mesh, "P", 1)

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

# -- Setting of the Taylor-Green vortex
# - Define time-dependent pressure boundary condition
p_boundary_exp = Expression("-1/4*(cos(2*x[0])+cos(2*x[1]))*exp(-4*t/Re)", t=0.0, Re=Re,  degree=2)

# - Time-dependent velocity boundary condition
v_boundary_exp = Expression(("cos(x[0])*sin(x[1])*exp(-2*t/Re)",
"-sin(x[0])*cos(x[1])*exp(-2*t/Re)"), t=0.0, Re=Re,
degree=2)

# - Define boundary location
boundary_str = "x[0] < 100*DOLFIN_EPS | x[0] > 2*DOLFIN_PI - 100*DOLFIN_EPS | \
x[1] < 100*DOLFIN_EPS | x[1] > 2*DOLFIN_PI - 100*DOLFIN_EPS"

# - initialize boundary conditions
velocity_bc  = DirichletBC(V, v_boundary_exp, boundary_str)
pressure_bc = DirichletBC(Q, p_boundary_exp, boundary_str)

# - store boundary conditions
bcu = [velocity_bc]
bcp = [pressure_bc]

# -- Create given functions
# - initialize initial condition
u0 = interpolate(v_boundary_exp,V)
u1 = Function(V)
p0 = interpolate(p_boundary_exp,Q)
p1 = Function(Q)

# Define coefficients
# - time step size
k = Constant(dt)
# - right side
f = Constant((0, 0))

# -- Variational formulation of the CHORIN PROJECTION
# - variational Formulation of the tentative velocity step
F1 = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u, v)*dx + \
a1 = lhs(F1)
L1 = rhs(F1)

# - variational Formulation of the pressure update
L2 = -(1/k)*div(u1)*q*dx

# - variational Formulation of the velocity update
a3 = inner(u, v)*dx
L3 = inner(u1, v)*dx - k*inner(grad(p1), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Use amg preconditioner if available
prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"

# - initialise first time step
t = dt

# --- Time Iterations
while t <= T + (dt/2):
print("At time t = ","{:.5f}".format(t))

# - Update pressure boundary condition
p_boundary_exp.t = t
v_boundary_exp.t = t

# --- Computation of Chorin-Projection
# - Compute tentative velocity step
b1 = assemble(L1)
[bc.apply(A1, b1) for bc in bcu]
solve(A1, u1.vector(), b1, "gmres", "default")

# - Compute pressure correction
b2 = assemble(L2)
[bc.apply(A2, b2) for bc in bcp]
solve(A2, p1.vector(), b2 ,"cg")

# -compute velocity correction
b3 = assemble(L3)
[bc.apply(A3, b3) for bc in bcu]
solve(A3, u1.vector(), b3, "gmres", "default")

# --- Computation of Errors ----
# - computing exact solutions
u_e = interpolate(v_boundary_exp, V)
p_e = interpolate(p_boundary_exp, Q)
# -- computing errors
vel_L2 = errornorm(u_e, u1, "L2")
vel_max = np.max(np.abs(u_e.vector()-u1.vector()))
p_L2 = errornorm(p_e, p1, "L2")
p_max = np.max(np.abs(p_e.vector()-p1.vector()))
print("L2 velocity error   : ", vel_L2)
print("L2 pressure error   : ", p_L2)
print("Max velocity error  : ", vel_max)
print("Max pressure error  : ", p_max)

#-------------- NEXT TIME STEP -------------
# --- Move to next time step
u0.assign(u1)
t += dt``````

To me your choice of boundary conditions seems strange. I would not have a Dirichlet condition for both the pressure and velocity on the same boundary (see for instance: https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/navier-stokes/demo_navier-stokes.py).

If you choose to have only dirichlet conditons for the pressure, you need to either:
a) remove the pressure nullspace from the pressure update step.
b) Average the pressure after solving, as shown in: Order of convergence for a Navier-Stokes solver

Option b) I have already tried, unfortunately without success. I implemented the following code in case you want to check that.

``````p_e = interpolate(p_boundary_exp, Q)
vol_mesh = assemble(Constant(1.0)*dx(mesh))
p_exact_mean = assemble(p_e*dx(mesh))/vol_mesh
p_h_mean = assemble(p1*dx)/vol_mesh
p1.vector()[:] += p_exact_mean - p_h_mean
``````

I also eased the Dirichlet B.C. for the pressure without any major improvement.
The same goes for easing the Dirichlet B.C. for the velocity while keeping the B.C. for the pressure.

Is there any exact number of Boundary conditions and initial conditions necessary in this case?

How could I remove the pressure nullspace from the pressure update step?

Is there anything which could go wrong in the interpolation of the exact solution? (saying I am getting closer to the exact solution, the computed error is just misleading?9

Best regards!

Further I checked, my solution for the pressure is clearly non-zero since the max. nodal value as well as the L2 norm are not close to zero.

As it is quite tedious to review your code, here is my attempt of helping you.

Consider the following code, which is adapted from the Chorin splitting scheme demo in the dolfin demos. Here i compute the L2-L2 convergence rates for the splitting scheme for a wide range of time steps. Note that where I have a Dirichlet condition for the tentative velocity step, I have a neumann condition in the pressure update step.

``````# Copyright (C) 2010-2011 Anders Logg
# the Free Software Foundation, version 3
# Modified by Mikael Mortensen 2011
# This is an adapted version of the Navier-Stokes Chorin scheme
# based on the version found at:
# https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/navier-stokes/demo_navier-stokes.py
#
# Author of adapted version: Jørgen S. Dokken 2020

from dolfin import *
import numpy as np
set_log_level(LogLevel.ERROR)
def chorin(N, dt=0.1):
# Create mesh for N-th refinement lvl
mesh = UnitSquareMesh(4,4)
for i in range(N):
# Refine the original mesh to get more accurate convergence studies
mesh = refine(mesh)

# Define function spaces (P2-P1)
V = VectorFunctionSpace(mesh, "Lagrange", 2)
Q = FunctionSpace(mesh, "Lagrange", 1)

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

# Set parameter values
T = 1
nu = 0.01

# Define time-dependent pressure boundary condition
u_ex_str = ('-(cos(pi*(x[0]))*sin(pi*(x[1])))*exp(-2.0*nu*pi*pi*t)',
' (cos(pi*(x[1]))*sin(pi*(x[0])))*exp(-2.0*nu*pi*pi*t)')
p_ex_str = '-0.25*(cos(2*pi*(x[0])) + cos(2*pi*(x[1])))*exp(-4.0*nu*pi*pi*t)'
p_ex = Expression(p_ex_str,nu=nu, t=0.0, degree=3)
u_ex = Expression(u_ex_str, nu=nu, t=0.0, degree=4)

# Define boundary conditions
bc_u  = DirichletBC(V, u_ex, "on_boundary")
bcu = [bc_u]

# Create functions and initial condition for u
u0 = interpolate(u_ex, V)
u1 = Function(V)
p1 = Function(Q)

# Define coefficients
k = Constant(dt)
f = Constant((0, 0))

# Tentative velocity step
F1 = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx + \
a1 = lhs(F1)
L1 = rhs(F1)

# Pressure update
L2 = -(1/k)*div(u1)*q*dx

# Velocity update
a3 = inner(u, v)*dx
L3 = inner(u1, v)*dx - k*inner(grad(p1), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Use amg preconditioner if available
prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"

# Use nonzero guesses - essential for CG with non-symmetric BC
parameters['krylov_solver']['nonzero_initial_guess'] = True

uL2 = []
pL2 = []
# Time-stepping
t = dt
while t < T + DOLFIN_EPS:

# Update pressure boundary condition
p_ex.t = t
u_ex.t = t

# Compute tentative velocity step
b1 = assemble(L1)
[bc.apply(A1, b1) for bc in bcu]
solve(A1, u1.vector(), b1, "bicgstab", "default")

# Pressure correction
b2 = assemble(L2)
solve(A2, p1.vector(), b2, "bicgstab", prec)
# Correct for no Dirichet condition on pressure
p1.vector()[:] -= assemble(p1*dx)/assemble(1*dx(domain=mesh))

# Velocity correction
b3 = assemble(L3)
[bc.apply(A3, b3) for bc in bcu]
solve(A3, u1.vector(), b3, "bicgstab", "default")

# Move to next time step
u0.assign(u1)
t += dt

# Save error per timestep
uL2.append(errornorm(u_ex, u1, norm_type="L2"))
pL2.append(errornorm(p_ex, p1, norm_type="L2"))

# Accumulate space-time error
u_L2L2 = sqrt(sum([x**2 for x in uL2]))
p_L2L2 = sqrt(sum([x**2 for x in pL2]))
return u_L2L2, p_L2L2

for dt in [0.1*0.5**i for i in range(5)]:
e_u = []
e_p = []
for N in range(5):
ue,pe = chorin(N,dt=dt)
e_u.append(ue)
e_p.append(pe)
e_u = np.array(e_u)
e_p = np.array(e_p)
print("-"*25)
print("Error u (dt{0:.2e}): ".format(dt), e_u)
print("Error p (dt{0:.2e}): ".format(dt), e_p)

print("U rate: (dt{0:.2e})".format(dt), np.log(e_u[1:]/e_u[:-1])/np.log(0.5))
print("P rate: (dt{0:.2e})".format(dt), np.log(e_p[1:]/e_p[:-1])/np.log(0.5))
``````

This yields the following spatial errors and convergence rates:

``````-------------------------
Error u (dt1.00e-01):  [ 0.27578819  0.09456845  0.09414062  0.108523    0.11805302]
Error p (dt1.00e-01):  [ 0.14417172  0.06188845  0.05871263  0.06359949  0.06739817]
U rate: (dt1.00e-01) [ 1.54412989  0.00654156 -0.20511162 -0.12143405]
P rate: (dt1.00e-01) [ 1.22004608  0.07599932 -0.11534428 -0.08369427]
-------------------------
Error u (dt5.00e-02):  [ 0.31379833  0.07535302  0.06627579  0.0769349   0.08493375]
Error p (dt5.00e-02):  [ 0.17590836  0.05587598  0.04861583  0.05273564  0.05633692]
U rate: (dt5.00e-02) [ 2.05810048  0.18518332 -0.21515619 -0.14269986]
P rate: (dt5.00e-02) [ 1.65452377  0.20080212 -0.11735218 -0.09530229]
-------------------------
Error u (dt2.50e-02):  [ 0.40594611  0.06697529  0.05000336  0.05807345  0.0646237 ]
Error p (dt2.50e-02):  [ 0.23011328  0.05376682  0.0392755   0.0420956   0.0451641 ]
U rate: (dt2.50e-02) [ 2.59958727  0.42160399 -0.21585364 -0.1541847 ]
P rate: (dt2.50e-02) [ 2.0975562   0.45308651 -0.1000399  -0.10150706]
-------------------------
Error u (dt1.25e-02):  [ 0.55350685  0.06619998  0.03752667  0.04304778  0.04812718]
Error p (dt1.25e-02):  [ 0.31255209  0.05834376  0.03152564  0.0321475   0.03446741]
U rate: (dt1.25e-02) [ 3.06369847  0.81891467 -0.19802274 -0.16091287]
P rate: (dt1.25e-02) [ 2.42144626  0.88805264 -0.02818072 -0.10052625]
-------------------------
Error u (dt6.25e-03):  [ 0.76490413  0.07370138  0.0280446   0.03114607  0.03468403]
Error p (dt6.25e-03):  [ 0.43101697  0.07161532  0.02689323  0.02385429  0.02554243]
U rate: (dt6.25e-03) [ 3.37551549  1.39396869 -0.15132743 -0.15522138]
P rate: (dt6.25e-03) [ 2.58940451  1.41302546  0.17299432 -0.09864732]
``````

You can do a similar study for temporal convergence rates.

1 Like

Thanks for the reply and the testing code!

Actually I just figured out my problem… and to be honest I am quite embarassed by the solution, but it is only fair to share it nethertheless.
Since I only slightly modified the Expression p_boundary_exp from a previous code I did not expect this to be the cause and for some reason forgot to check this as an error source although it is quite obvious…
So since I have not had any previous experience with C++, I replaced “-rho/4” by “-1/4”. What happens is that since rho is given by the function it will be given as a float, however “-1/4” is interpreted as an integer, therefore has the integer solution 0 instead of the float .25. Now after changing “-1/4” to “-1.0/4” everything works fine.

I am really sorry for bothering you guys with my problem, I should have checked this beforehand, hope this thread still helps someone else as clueless about C++ as me in the future.

Best regards,
Max

so the source term f in your code just zero, I do not understand this.

Could you be more specific? What makes it cobfusing to you? Do you mean it should be something else? See for instance: https://en.m.wikipedia.org/wiki/Taylor–Green_vortex

I mean that using P2-P1 the convergence rate can is 3 and 2 order in space. (or come out the superconvergence 4 and 2 order) . The results is not better. When I use the pressure correction scheme, with Dirchelt bc in velocity, and the results is wrong too. So I just want to test a correct code that can have 3 and 2 order convergence with the p2-p1 element.

As you observe in the result above, the spatial convergence improves as you reduce the time step. This is natural as if the temporal discretization is too coarse, its error is dominating. Continued temporal refinement will give you your convergence rates. Note that the code I supplied is simply an adaptation of the FEniCS-demo, and there are many improvements that can be made (using splitting schemes with second order Adams-Bashforth approximations). See for instance the Oasis-paper and the corresponding software package, or my Multi-mesh splitting scheme paper and corresponding implementation or the corresponding Chorin-derivation.

Thanks a lot, I will study them carefully.

Best regards.