Unsteady incompressible Navier-Stokes equation for DFG3D-2 Benchmark

I’m trying to solve the DFG3D-2 Benchmark (3D laminar case (Re=20 and Re=100) - Featflow) using Crank Nicolson’s discretization. I used the same mesh from @dokken 's gmsh tutorial (Using the GMSH Python API to generate complex meshes | Jørgen S. Dokken). I modified the code from the FEniCSx tutorial to suit this particular problem to the best of my ability. The code produces no errors but the results are poor. I am using v0.6.0 on conda.
The code I used is:

gdim = 3
mesh_comm = MPI.COMM_WORLD
model_rank = 0

inlet_marker, outlet_marker, wall_marker, obstacle_marker = 1, 3, 5, 7

mesh, cell_tags, ft = read_from_msh("mesh3D.msh", MPI.COMM_WORLD, 0,
ft.name = "Facet markers"

t = 0
T = 1                       # Final time
dt = 0.01                 # Time step size
num_steps = int(T/dt)
k = Constant(mesh, PETSc.ScalarType(dt))        
mu = Constant(mesh, PETSc.ScalarType(0.001))  # Dynamic viscosity
rho = Constant(mesh, PETSc.ScalarType(1))     # Density

v_cg2 = VectorElement("CG", mesh.ufl_cell(), 2)
s_cg1 = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, v_cg2)
Q = FunctionSpace(mesh, s_cg1)

fdim = mesh.topology.dim - 1

# Define boundary conditions
class InletVelocity():
    def __init__(self, t):
        self.t = t

    def __call__(self, x):
        values = np.zeros((gdim, x.shape[1]),dtype=PETSc.ScalarType)
        values[0] = 16 * 2.25 * np.sin(self.t * np.pi/8) * x[1] * x[2] * (0.41 - x[1]) * (0.41 - x[2])/(0.41**4)
        return values

# Inlet
u_inlet = Function(V)
inlet_velocity = InletVelocity(t)
bcu_inflow = dirichletbc(u_inlet, locate_dofs_topological(V, fdim, ft.find(inlet_marker)))
# Walls
u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bcu_walls = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(wall_marker)), V)
# Obstacle
bcu_obstacle = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(obstacle_marker)), V)
bcu = [bcu_inflow, bcu_obstacle, bcu_walls]
# Outlet
bcp_outlet = dirichletbc(PETSc.ScalarType(0), locate_dofs_topological(Q, fdim, ft.find(outlet_marker)), Q)
bcp = [bcp_outlet]

u = TrialFunction(V)
v = TestFunction(V)
u_ = Function(V)
u_.name = "u"
u_s = Function(V)
u_n = Function(V)
u_n1 = Function(V)
p = TrialFunction(Q)
q = TestFunction(Q)
p_ = Function(Q)
p_.name = "p"
phi = Function(Q)

f = Constant(mesh, PETSc.ScalarType((0,0,0)))
F1 = rho / k * dot(u - u_n, v) * dx 
F1 += inner(dot(1.5 * u_n - 0.5 * u_n1, 0.5 * nabla_grad(u + u_n)), v) * dx
F1 += 0.5 * mu * inner(grad(u + u_n), grad(v))*dx - dot(p_, div(v))*dx
F1 += dot(f, v) * dx
a1 = form(lhs(F1))
L1 = form(rhs(F1))
A1 = create_matrix(a1)
b1 = create_vector(L1)

a2 = form(dot(grad(p), grad(q))*dx)
L2 = form(-1/k * dot(div(u_s), q) * dx)
A2 = assemble_matrix(a2, bcs=bcp)
b2 = create_vector(L2)

a3 = form(dot(u, v)*dx)
L3 = form(dot(u_s, v)*dx - k * dot(nabla_grad(phi), v)*dx)
A3 = assemble_matrix(a3)
b3 = create_vector(L3)

# Solver for step 1
solver1 = PETSc.KSP().create(mesh.comm)
pc1 = solver1.getPC()

# Solver for step 2
solver2 = PETSc.KSP().create(mesh.comm)
pc2 = solver2.getPC()

# Solver for step 3
solver3 = PETSc.KSP().create(mesh.comm)
pc3 = solver3.getPC()

n = -FacetNormal(mesh) # Normal pointing out of obstacle
dObs = Measure("ds", domain=mesh, subdomain_data=ft, subdomain_id=obstacle_marker)
u_t = inner(as_vector((n[1], -n[0], 0)), u_)
drag = form(2 / 0.225 * (mu / rho * inner(grad(u_t), n) * n[1] - p_ * n[0]) * dObs)
lift = form(-2 / 0.225 * (mu / rho * inner(grad(u_t), n) * n[0] + p_ * n[1]) * dObs)
if mesh.comm.rank == 0:
    C_D = np.zeros(num_steps, dtype=PETSc.ScalarType)
    C_L = np.zeros(num_steps, dtype=PETSc.ScalarType)
    t_u = np.zeros(num_steps, dtype=np.float64)
#    t_p = np.zeros(num_steps, dtype=np.float64)

progress = tqdm.autonotebook.tqdm(desc="Solving PDE", total=num_steps)
for i in range(num_steps):
    # Update current time step
    t += dt
    # Update inlet velocity
    inlet_velocity.t = t

    # Step 1: Tentative velocity step
    assemble_matrix(A1, a1, bcs=bcu)
    with b1.localForm() as loc:
    assemble_vector(b1, L1)
    apply_lifting(b1, [a1], [bcu])
    b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b1, bcu)
    solver1.solve(b1, u_s.vector)

    # Step 2: Pressure corrrection step
    with b2.localForm() as loc:
    assemble_vector(b2, L2)
    apply_lifting(b2, [a2], [bcp])
    b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b2, bcp)
    solver2.solve(b2, phi.vector)

    p_.vector.axpy(1, phi.vector)

    # Step 3: Velocity correction step
    with b3.localForm() as loc:
    assemble_vector(b3, L3)
    b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    solver3.solve(b3, u_.vector)

    # Update variable with solution form this time step
    with u_.vector.localForm() as loc_, u_n.vector.localForm() as loc_n, u_n1.vector.localForm() as loc_n1:

    drag_coeff = mesh.comm.gather(assemble_scalar(drag), root=0)
    lift_coeff = mesh.comm.gather(assemble_scalar(lift), root=0)
    if mesh.comm.rank == 0:
        t_u[i] = t
        C_D[i] = sum(drag_coeff)
        C_L[i] = sum(lift_coeff)

The main changes are the inlet boundary conditions and the CD and CL formulae. I ran the program upto T = 1 with dt = 0.01.
The CD and CL plots I got are:

I can’t seem to figure out what is wrong with the code. I got the same results using the incremental pressure correction scheme(IPCS). Any suggestions on what’s wrong will be beneficial. Also, could anyone suggest any other benchmark other than FeatFlow to compare CD and CL results from Crank Nicolson, IPCS or Chorin’s projection with.

Have you looked at the actual flow profile? Does it look sensible?

Have you check that the solvers converge?

The flow profile results had a velocity of 1.2e-38 when I checked the xdmf in ParaView.
I don’t know how to check if the solvers converged.

You are setting a 0 bc on the outer boundaries. You would need to use the Glyphs filter or Clip to look at the interior flow profile. See the figure below.

I am working from a laptop, and do not have the time to run your whole problem.

Have you also checked what you get as drag and lift when you use the original drag and lift coefficients from my tutorial?

CD increased a little bit because the denominator is 0.1(mine was 0.225). CL remained 0.

Note that the meshes in the DFG suite are tailored to be accurate in certain areas. You might get better results using a second order mesh, or a finer mesh, or a smaller time step.

If you are getting the same results with both a traditional IPCS and a more complicated CN/Adam-Bashforth formulation, it indicates that you at least have a consistent formulation.

I ran the program with the CN/Adam-Bashforth formulation for Taylor-Hood elements of P3 for Velocity and P2 for Pressure up to T=0.1. The results were just as poor despite having 3.6 million dofs.

I am using an intel i5 MacBook pro. It took 53 sec for each iteration despite running it across 4 processors. Is there any way of speeding it up, perhaps by changing the solver configuration?