Manufactured Solution for Navier Stokes

[gcorsi1]

I’m struggling to get decent results with a Manufactured Solution for the Navier-Stokes equations, which I solve with a standard IPCS-like method, which gives good results for the Taylor-Green 2D (periodic) solution of NS, and is based on Oasis.

The fields I am solving for are taken from an article and are defined on a 2D square of side length h:

\begin{align}\mathbf{u}_{MS,0}& = -\sin(2\pi t)( y(y-h) +x)/(2\nu), \\ \mathbf{u}_{MS,1} &= -\sin(2\pi t)( x(x-h) -y)/(2\nu),\\ p &= -\sin(2\pi t) (x+y). \end{align}.

These fields are such that the viscous terms balance the pressure terms, i.e. \nu \Delta \mathbf{u} = \nabla p. There remain two terms that have to be added to the RHS of the momentum equation, since these fields do not solve NS:

\mathbf{f} = \left( \partial \mathbf{u}_{MS}/ \partial t + ( \mathbf{u}_{MS}\cdot \nabla) \mathbf{u}_{MS} \right),

where MS of course stands for manufactured solution. The boundary conditions are Dirichlet on two of the sides of the square domain, and Neumann on the others.

However, I first start from the simpler case resulting from assigning DirichletBC conditions on all of the outer boundary.

At every time step, I calculate the extra contribution to the RHS and add it to the system of equations

    x = SpatialCoordinate(mesh)
    u_ms = as_vector( (-sin(2.0*pi*t)/(2.0*nu)*(x[1]*(x[1]-1.0) + x[0]), -sin(2.0*pi*t)/(2.0*nu)*(x[0]*(x[0]-1.0) - x[1])) )
    u_ddt_ = as_vector( (-2.0*pi*cos(2.0*pi*t)/(2.0*nu)*(x[1]*(x[1]-1.0) + x[0]),-2.0*pi*cos(2.0*pi*t)/(2.0*nu)*(x[0]*(x[0]-1.0) - x[1])) )
    f_ = u_ddt_ + dot(u_ms, nabla_grad(u_ms))

then the assembly is (v is test function of the velocity discrete space):

    srcv = assemble( v*f_[i]*dx )
    f_tmp[ui].axpy(1.0, srcv)

where f_tmp is the RHS vector of the discrete system of the momentum equation, that is, when the tentative velocity solution is solved for. I initialize the velocity and pressure fields with the expressions above, and I calculate the L^2 error as follows:

    ue = interpolate(ue, V)
    uen = norm(ue.vector())
    ue.vector().axpy(-1, q_[ui].vector())
    error = norm(ue.vector()) / uen

where ue is one of the expressions above and V is the discrete space where the solutions are defined, the vectors of the solutions being q\_ [ui] for each component of the velocity, or pressure. If I have Neumann BCs on some boundaries, I also add an additional term to the RHS:

   grdu = assemble( inner(v, dot(n,nabla_grad(u_ms[i])) )*ds )
   f_tmp[ui].axpy(1.0, grdu)

that is the normal gradient of the known solution, \mathbf{u}_{MS}. Velocity and pressure spaces are P2/P1 as a starting point.

The problem is that results for this case don’t work as expected: if I assign DirichletBC everywhere on the boundary, the error stays of the order of 10^{-1} for the velocity components. This does not improve if time step is decreased, or if velocity and pressure degrees are increased to 4 and 3, respectively. So no convergence is observed.

If, on part of the boundary Neumann BCs are assigned, things get much worse, and the solution looks like garbage. I cannot understand what it is that ruins the results, since the base code works fine for the Taylor Green 2D case.

Do you have any hints on how to implement correctly Manufactured Solutions for Navier-Stokes in FEniCS (with Chorin like splitting)??

[dokken]

As you have not supplied your code, it is tricky to know where things go wrong in your case.

Note that by assigning dirichlet conditions on the whole outer boundary, you have to normalize the pressure as it is only determined up to a constant.

The first thing I would do is to verify that the Taylor-Green example gives correct convergence rates when you use Dirichlet conditions instead of Periodic conditions, as IPCS-like methods are prone to boundary conditions.

The second thing I would do is to remove the manual implementation of the source function with something like this:

from dolfin import *
import ufl
mesh = UnitSquareMesh(10,10)
t = Constant(0.8)
h = Constant(1)
nu = Constant(1)
x,y = SpatialCoordinate(mesh)

u_man = ufl.as_vector((-sin(2*pi*t)*(y*(y-h)+x)/2/nu,
                       -sin(2*pi*t)*(x*(x-h)-y)/2/nu))
p_man = -sin(2*pi*t)*(x+y)
def E(u,p,t):
    return ufl.diff(u, t) + grad(u)*u -nu*div(grad(u))+grad(p)

Second thing I would do is to verify your scheme against the Turek DFG2D benchmark.
I have implemented the IPCS splitting scheme for this case (with Multimesh, but you can just remove every dI and dO term and replace dX with dx).

[gcorsi1]

Thanks @dokken, I realize I have not provided much code, but since I use Oasis I don’t know exactly how to provide it in full.

I tried the steps you suggested:

• when I change from periodic to Dirichlet BCs for the Taylor-Green 2D case, the convergence rates deteriorate slightly but are still OK, the pressure convergence order with the time step goes from around 2 to 1.66 which is expected I guess because of artificial boundary layers at Dirichlet boundaries. Convergence rate for velocity with respect to time step is 1.95. I only checked for two different time step values, but the results are still acceptable I think.
• I implemented your code for the Turek benchmark (mesh and post-processing) into Oasis, and I think it works well, I messed up in the calculation of lift and drag but the pressure difference seems good, the orange and blue data overlap, and the absolute error looks similar to the figure in your article on arxiv:
  
  

I’m a bit confused regarding why the manufactured solution case still doesn’t work, and gives errors in order 10^{-1} regardless of what I do. The only difference with the other cases should be in the addition of a source term, and maybe that the fields are more unsteady in time.

At every time step, now I assemble the RHS contribution for the momentum equations as follows

    x,y = SpatialCoordinate(mesh)
    h = Constant(1.0)
    t_ = Constant(t)
    u_man = ufl.as_vector((-sin(2*pi*t_)*(y*(y-h)+x)/2/nu,
                        -sin(2*pi*t_)*(x*(x-h)-y)/2/nu))
    p_man = -sin(2*pi*t_)*(x+y)
    if inner_iter == 1:
        src = [ E(u_man,p_man,t_)[0], E(u_man,p_man,t_)[1] ]
        for i,ui in enumerate(u_components):
            srcv = assemble( v*src[i]*dx )
            f_tmp[ui].axpy(1.0, srcv)

where the same is done for every component of velocity ui, and E was defined as in your reply (and v is a test function of the velocity space)

    def E(u,p,t):
        return ufl.diff(u, t) + grad(u)*u -nu*div(grad(u))+grad(p)

Is there any part of code in particular I should provide? I will try to make it fit into a discourse answer.

[dokken]

What happens if you used this latter approach on the Taylor-Green problem? As one knows the analytical solutions there, it should produce the same results (as you are using a square domain in both cases).

If it works for Taylor-green by adding the source term in this way (you should only change u_man and p_man), it suggests that the temporal discretization of the new problem is not fine enough.

If it doesn’t work, it suggests that Oasis has some underlying assumptions.

[gcorsi1]

Dear @dokken,

@dokken

What happens if you used this latter approach on the Taylor-Green problem? As one knows the analytical solutions there, it should produce the same results (as you are using a square domain in both cases).

that sounds like a nice idea, however the forcing term would be zero for the Taylor Green flow as it’s a solution for Navier-Stokes, therefore I should add only some terms of the function E as defined above, e.g. I could add only the viscous and convective parts, but which errors should I check for then? I am unsure since the NS equations are nonlinear so maybe these operations are not as simple as they look.

[dokken]

You could also check the stability of your problem with respect to temporal discretization by using the Navier Stokes demo in FEniCS

[gcorsi1]

Thanks for your suggestion @dokken, using the Navier Stokes demo improved things, as the velocity error now shows some kind of convergence with the value of the time step. Errors for the pressure are lower than before, but still very high though, and they don’t converge. Since it’s simpler I can now provide the full code, even though I still don’t understand why this would work differently that Oasis.

import matplotlib.pyplot as plt
from dolfin import *
import ufl

Nc = 20
mesh = UnitSquareMesh(MPI.comm_world,Nc,Nc)
x,y = SpatialCoordinate(mesh)
n = FacetNormal(mesh)

# Define function spaces (P2-P1)
V = VectorFunctionSpace(mesh, "Lagrange", 2)
Q = FunctionSpace(mesh, "Lagrange", 1)
VV = [V,Q]
# Define trial and test functions
u = TrialFunction(V)
p = TrialFunction(Q)
v = TestFunction(V)
q = TestFunction(Q)

# Set parameter values
dt = 0.01
T = 1.1
nu = Constant(0.1)
h = Constant(1.0)

# Define domains and boundary measure
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary 

domains = MeshFunction("size_t", mesh, 1)
domains.set_all(0)
DirichletBoundary().mark(domains, 1)
ds = Measure("ds", domain=mesh, subdomain_data=domains)

# Define time-dependent pressure and velocity boundary conditions
up_ex_str = ['-sin(2.0*pi*t)/(2.0*nu)*(x[1]*(x[1]-1.0) + x[0])', '-sin(2.0*pi*t)/(2.0*nu)*(x[0]*(x[0]-1.0) - x[1])', '-sin(2.0*pi*t)*(x[0] + x[1])']
u_ex = Expression((up_ex_str[0],up_ex_str[1]),t=0.0,nu=nu,degree=4)
p_ex = Expression(up_ex_str[2],t=0.0,nu=nu,degree=4)
# UFL form for the source term
def E(u,p,t_):
    return ufl.diff(u, t_) + grad(u)*u -nu*div(grad(u))+grad(p)
def u_man(x,y,t_): 
    return ufl.as_vector((-sin(2*pi*t_)*(y*(y-h)+x)/2/nu, -sin(2*pi*t_)*(x*(x-h)-y)/2/nu))
def p_man(x,y,t_):
    return -sin(2*pi*t_)*(x+y)
# Define boundary conditions
noslip  = DirichletBC(V, u_ex, domains, 1)
bcu = [noslip]
bcp = []

# Create functions
u0 = Function(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 + \
     nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Pressure update
a2 = inner(grad(p), grad(q))*dx
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)

#---------------------------------------------------
# ATTACH PRESSURE NULL-SPACE IF NEEDED
if bcp == []:
    null_vec = Vector(p1.vector())
    Q.dofmap().set(null_vec, 1.0)
    null_vec *= 1.0 / null_vec.norm('l2')
    Aa = as_backend_type(A2)
    null_space = VectorSpaceBasis([null_vec])
    Aa.set_nullspace(null_space)

# Use amg preconditioner if available -- define pressure solver
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

# Create files for storing solution
ufile = File("results/velocity.pvd")
pfile = File("results/pressure.pvd")

# Time-stepping
t = dt
while t < T + DOLFIN_EPS:

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

    # Compute tentative velocity step
    b1 = assemble(L1)
    t_ = Constant(t)
    um,pm = u_man(x,y,t_), p_man(x,y,t_)
    ff = assemble(inner(E(um,pm,t_),v)*dx)
    b1.axpy(1.0,ff)
    [bc.apply(A1, b1) for bc in bcu]
    solve(A1, u1.vector(), b1, "bicgstab", "default")

    # Pressure correction
    b2 = assemble(L2)
    [bc.apply(A2, b2) for bc in bcp]
    [bc.apply(p1.vector()) for bc in bcp]
    if bcp == []:
        null_space.orthogonalize(b2)
    solve(A2, p1.vector(), b2, "bicgstab", prec)
    if bcp == []:
            normalize(p1.vector())

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

    # Save to file
    ufile << u1
    pfile << p1

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

    #compute errors
    err = {}
    SS = [u1,p1]
    for i, ui in enumerate([u_ex,p_ex]):
        ue = interpolate(ui, VV[i])
        uen = norm(ue.vector())
        ue.vector().axpy(-1, SS[i].vector())
        error = norm(ue.vector()) / uen
        err[i] = "{0:2.6e}".format(norm(ue.vector()) / uen)
    if MPI.rank(MPI.comm_world) == 0:
        print("Error is ", err, " at time = ", t)

# Plot solution
plt.figure()
plot(p1, title="Pressure")

plt.figure()
plot(u1, title="Velocity")
try:
    plt.show()
except:
    pass

The velocity looks OK

but the pressure is just wrong, from the fields I have specified it should be linear

[dokken]

I would guess one of the issues is that since you use a Dirichlet condition on all boundaries, the null space normalization assumes that integral(p*dx) =0. This is not the case for your manufactured solution.