Taylor Green Vortex - Coupled Navier Stokes

variational formulation
1

Hi !

I’m trying to solve the Taylor-Green Vortex (TGV) problem by using directly the weak form of Navier Stokes without much success.

TGV on 0\leq x,y\leq \pi and \nu = 0.01 with this exact solution:
U(x,y,t) = -cos(x)sin(y)exp(-2\nu t)
V(x,y,t) = sin(x)cos(y)exp(-2\nu t)
P(x,y,t) =- 0.25 (cos(2x)+cos(2y))exp(-4\nu t)

On the 4 edges, is defined a DirichletBC by using the velocity component from the exact solution.
At t=0 the initial condition is the exact solution at t=0. The error between the exact and approximated solution is calculated that way:

e_V = \sqrt{\int_{\Omega}(u_{exact}-u_{approx})^2} and e_P = \sqrt{\int_{\Omega}(p_{exact}-p_{approx})^2} at each timestep.

I have at least two different questions:

  1. To make the pressure field converging I need to add 10^{-6}\times p in the weak form. I think it might come from the fact that I have no pressure condition and only \nabla p appear in the NS equation.

  2. Once this term implemented I manage to have a relatively small error between the calculated field and the exact solution with a backward differentiation formula order 1 (BDF1).

Time = 0.937500 , L2_error velocity = 0.000006 and L2_error pressure = 0.000255
Time = 1.000000 , L2_error velocity = 0.000006 and L2_error pressure = 0.000255

Problem is, when I swap the BDF1 for a BDF2 and I initialize the first two timesteps with the exact solution the error increase:

Time = 0.937500 , L2_error velocity = 0.000922 and L2_error pressure = 0.000676
Time = 1.000000 , L2_error velocity = 0.000852 and L2_error pressure = 0.000641

If you can explain me either of those questions it would be much appreciated (especially for the second one) !

My MWE:

You need to download my meshfile here:
https://drive.google.com/file/d/1JIJHlTG7-DL9ViFJZK1BEBZBOg64ntTE/view?usp=sharing

Then just need to change BDF = 1 to BDF = 2 to swap the scheme.

from dolfin import *

set_log_level(50)

BDF = 1

comm = MPI.comm_world
nu = 0.01
dt = pow(2,-4)
mesh_file = "square_valid.h5"

#Exact solution
Exact = Expression(("-cos(x[0])*sin(x[1])*exp(-2.0*nu*t)",\
                    "sin(x[0])*cos(x[1])*exp(-2.0*nu*t)",\
                    "-0.25*(cos(2*x[0])+cos(2*x[1]))*exp(-4.0*nu*t)"),\
                   nu = nu, t = 0.0, degree = 3)

#No RHS
f = Constant((0.0, 0.0))

#Exact U imposed as Dirichlet BC
ExactU = Expression(("-cos(x[0])*sin(x[1])*exp(-2.0*nu*t)",\
                     "sin(x[0])*cos(x[1])*exp(-2.0*nu*t)"),\
                    nu = nu, t = 0.0, degree = 3)

#Read the mesh and the facet, 1 = top, 2 = right, 3 = bot, 4 = left
mesh = Mesh()
try:
    with HDF5File(comm, mesh_file, "r") as h5file:
        h5file.read(mesh, "mesh", False)
        facet = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
        h5file.read(facet, "facet")
except FileNotFoundError as fnf_error:
    print(fnf_error)
    
#Taylor-Hood element [P2, P2]-P1
V = VectorElement("CG", mesh.ufl_cell(), 2)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, V * Q)

#All the test/trial functions needed to solve the problem
v, q = TestFunctions(W)
w = Function(W)
u, p = split(w)
w_prev1 = Function(W)
u_prev1, p_prev1 = split(w_prev1)
w_prev2 = Function(W)
u_prev2, p_prev2 = split(w_prev2)

#BC comming from the exact solution
topBC = DirichletBC(W.sub(0), ExactU, facet, 1)
rightBC = DirichletBC(W.sub(0), ExactU, facet, 2)
botBC = DirichletBC(W.sub(0), ExactU, facet, 3)
leftBC = DirichletBC(W.sub(0), ExactU, facet, 4)

bcs = [topBC, rightBC, botBC, leftBC]

#Weak form of Navier-Stokes with BDF1 or 2
F = (inner(u - u_prev1, v) / Constant(dt)\
     + inner(grad(u) * u_prev1, v) + inner(grad(u_prev1) * u, v) - inner(grad(u_prev1) * u_prev1, v)\
     + Constant(nu) * inner(grad(u), grad(v)) - inner(p, div(v)) - inner(q, div(u)) - inner(f, v) - Constant(1e-6) * p * q) * dx

if BDF == 2:
    F = (inner(Constant(1.5) * u - Constant(2.0) * u_prev1 + Constant(0.5) * u_prev2, v) / Constant(dt)\
         + inner(grad(u) * u_prev1, v) + inner(grad(u_prev1) * u, v) - inner(grad(u_prev1) * u_prev1, v)\
         + Constant(nu) * inner(grad(u), grad(v)) - inner(p, div(v)) - inner(q, div(u)) - inner(f, v) - Constant(1e-6) * p * q) * dx

J = derivative(F, w)
problem = NonlinearVariationalProblem(F, w, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.parameters["newton_solver"]["linear_solver"] = "mumps"

#Interpolate the exact solution on w to use it as initial condition
Time = 0.0
Exact.t = Time
w.interpolate(Exact)

if BDF == 2:
    #If BDF2 interpolate the next timestep too
    Time += dt
    Exact.t = Time
    w_prev1.interpolate(Exact)

while Time < 1.0:
    Time += dt
    if BDF == 2:
        assign(w_prev2, w_prev1)
    assign(w_prev1, w)

    ExactU.t = Time
    solver.solve()

    #Calculate the L2 error on the velocity and pressure field
    Exact.t = Time
    w_prev2.interpolate(Exact) #w_prev2 will be dumped anyway
    error_v, error_p = sqrt(assemble((w_prev2.sub(0) - w.sub(0))**2 * dx)), sqrt(assemble((w_prev2.sub(1) - w.sub(1))**2 * dx))
    if comm.rank == 0:
        print("Time = %f , L2_error velocity = %f and L2_error pressure = %f" %(Time, error_v, error_p))
2

Given the rest of the time stepping logic, should you maybe make the following replacement?

if BDF == 2:
    #If BDF2 interpolate the next timestep too
    #Time += dt
    #Exact.t = Time
    #w_prev1.interpolate(Exact)

    w_prev1.interpolate(Exact)
    Time += dt
    Exact.t = Time
    w.interpolate(Exact)

With that, the error does not increase when switching to BDF = 2.

3

Oh gosh… You’re right, I’ve swapped the first two timesteps by doing that ! Thank you a lot for pointing that :+1:

4

Hey Thibaut, just stumbled upon this as I was looking to solve a problem of mine. Is there a reason you use a nonlinear solver instead of a linear one?