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 + \
(1/Re)*inner(grad(u), grad(v))*dx - inner(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)
# - variational Formulation of the pressure update
a2 = inner(grad(p), grad(q))*dx
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