I have a question regarding the case when solving Navier-Stokes and on some boundaries, both Neumann and Dirichlet boundary conditions are applied, on neighboring segments.
For example I considered the analytical solution:
\mathbf{u} = \begin{pmatrix}-\sin (2\pi y) \\ \sin(2\pi x) \end{pmatrix} \exp(-4\nu\pi^2 t)
p = -\cos(2\pi x) \cos(2\pi y) \exp (-8\nu \pi^2 t)
taken from https://arxiv.org/abs/1706.09252
The domain is square and on each side, half the length is an inflow and the other half an outflow. Where inflow occurs, Dirichlet boundary conditions are assigned, while Neumann boundary conditions are given on the outflow segments. I define the domains accordingly:
and solve Navier Stokes with the standard steps of the IPCS solution algorithm:
-
Momentum step \langle\frac{ 3/2 \mathbf{u} - 2 \mathbf{u}_0 + 1/2 \mathbf{u}_{00}}{\Delta t}, v \rangle_{\Omega} + \langle 2(\mathbf{u} _0 \cdot \nabla) \mathbf{u}_0 - (\mathbf{u}_{00}\cdot \nabla) \mathbf{u}_{00}, v\rangle_{\Omega} + \langle \sigma (\mathbf{u} , p_0), \epsilon(v)\rangle_{\Omega} + \\ \langle p_0, v\rangle_{\partial \Omega} - \langle \nu \nabla \mathbf{u} \cdot \mathbf{n}, v\rangle_{\partial \Omega} = 0. A Neumann boundary condition also appears, while the old pressure value p_0 is used. BDF2 time stepping.
-
Pressure step \langle \nabla p - p_0 , \nabla q \rangle_{\Omega} = -3/2 \frac{1}{\Delta t} \langle \text{div} \mathbf{u} , q\rangle_{\Omega} with Dirichlet boundary condition p = p^{n+1} on the Neumann boundaries.
-
Correction step \langle \mathbf{u}^{n+1},v \rangle = \langle \mathbf{u} + \frac{2}{3} \Delta t \nabla (p-p_0), v \rangle_{\Omega}
Simulations are run with standard P2/P1 finite elements until time T=1, with various time steps. The problem is that the pressure shows an irregular behaviour at the corners and at x=0, y=0 on the boundary, that is the points where Neumann and Dirichlet segments touch. Here is the divergence field after the beginning of the simulation:
which clearly shows at which points trouble is starting. The resulting pressure at T=1:
Can you help me figure out what is going wrong here? Things work if I assign Dirichlet over all of the boundary, but not in this case.
My code, based on the FEniCS tutorial
import matplotlib.pyplot as plt
from dolfin import *
import ufl
# Print log messages only from the root process in parallel
parameters["std_out_all_processes"] = False;
Nc = 32
#mesh = UnitSquareMesh(MPI.comm_world,Nc,Nc)
mesh = RectangleMesh(Point(-0.5, -0.5), Point(0.5, 0.5), Nc, Nc)
x,y = SpatialCoordinate(mesh)
n = FacetNormal(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
dt = 1e-2
T = 1.0
nu = Constant(0.025)
h = Constant(1.0)
# Define domains and boundary measure
class DirichletBoundary(SubDomain):
def inside(self, x, on_boundary):
return (near(x[0], -0.5, 1e-4) and x[1] > -DOLFIN_EPS) or (near(x[1], 0.5, 1e-4) and x[0] > -DOLFIN_EPS ) or (near(x[0], 0.5, 1e-4) and x[1] < DOLFIN_EPS ) or (near(x[1], -0.5, 1e-4) and x[0] < DOLFIN_EPS)
class NeumannBoundary(SubDomain):
def inside(self, x, on_boundary):
return (near(x[0], -0.5, 1e-4) and x[1] < DOLFIN_EPS) or (near(x[1], 0.5, 1e-4) and x[0] < DOLFIN_EPS ) or (near(x[0], 0.5, 1e-4) and x[1] > - DOLFIN_EPS ) or (near(x[1], -0.5, 1e-4) and x[0] > - DOLFIN_EPS)
domains = MeshFunction("size_t", mesh, 1)
domains.set_all(0)
DirichletBoundary().mark(domains, 1)
NeumannBoundary().mark(domains, 2)
ds = Measure("ds", domain=mesh, subdomain_data=domains)
with XDMFFile(MPI.comm_world, 'domains_WALL.xdmf') as f:
f.write(domains)
# Define time-dependent pressure and velocity boundary conditions
up_ex_str = ['-sin(2.0*pi*x[1])*exp(-4.0*nu*pi*pi*t)', 'sin(2.0*pi*x[0])*exp(-4.0*nu*pi*pi*t)', '-cos(2.0*pi*x[0])*cos(2.0*pi*x[1])*exp(-8.0*nu*pi*pi*t)']
u_ex = Expression((up_ex_str[0],up_ex_str[1]),element=V.ufl_element(),domain=mesh,t=0.0,nu=nu)
p_ex = Expression(up_ex_str[2],element=Q.ufl_element(),domain=mesh,t=0.0,nu=nu)
# UFL form for the source term
# Define boundary conditions
bcu = [DirichletBC(V, u_ex, domains, 2)]#, DirichletBC(V, u_ex, domains, 2)]
bcp = [DirichletBC(Q, p_ex, domains, 1)]
# Create functions
u0 = Function(V)
u00 = Function(V)
u1 = Function(V)
p1 = Function(Q)
p0 = Function(Q)
div1 = Function(Q)
fi,VV = [[u0,p0], [V,Q]]
# Initialize the solutions for old time steps
for i,ui in enumerate([('-sin(2.0*pi*x[1])*exp(-4.0*nu*pi*pi*t)', 'sin(2.0*pi*x[0])*exp(-4.0*nu*pi*pi*t)'), '-cos(2.0*pi*x[0])*cos(2.0*pi*x[1])*exp(-8.0*nu*pi*pi*t)']):
deltat = dt / 2. if i==1 else 0.
vv = interpolate(Expression(ui, element=VV[i].ufl_element(),t=0.0 + deltat, nu=nu), VV[i])
fi[i].vector()[:] = vv.vector()[:]
if not i == 1:
deltat = -dt
vv = interpolate(Expression(ui, element=VV[i].ufl_element(), t=0.0 + deltat, nu=nu), VV[i])
u00.vector()[:] = vv.vector()[:]
# Define coefficients
U = 0.5*(u0 + u)
n = FacetNormal(mesh)
k = Constant(dt)
f = Constant((0, 0))
# Define symmetric gradient
def epsilon(u):
return sym(nabla_grad(u))
# Define stress tensor
def sigma(u, p):
return 2*nu*epsilon(u) - p*Identity(len(u))
# Define variational problem for step 1
F1 = dot(( 1.5*u - 2.0*u0 +0.5*u00) / k, v)*dx \
+ dot(dot( 2.0*u0, nabla_grad(u0)), v)*dx - dot(dot( u00, nabla_grad(u00)), v)*dx\
+ inner(sigma(u, p0), epsilon(v))*dx \
+ dot(p0*n, v)*ds - dot(nu*dot(nabla_grad(u_ex),n), v)*ds\
- dot(f, v)*dx
# + dot(dot(u0, nabla_grad(u0)), v)*dx\
# + dot(grad(p0), v)*dx - dot(nu*nabla_grad(U)*n, v)*ds \
a1 = lhs(F1)
L1 = rhs(F1)
# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p0), nabla_grad(q))*dx - 1.5*(1/k)*div(u1)*q*dx
# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u1, v)*dx - 2./3.*k*dot(nabla_grad(p1 - p0), v)*dx
# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
# Attach pressure nullspace 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)
Aa.null_space = 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")
divfile = File("results/divergence.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_)
[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")
div1.assign( project(div(u1),Q) )
# Save to file
ufile << u1
pfile << p1
divfile << div1
# Move to next time step
u00.assign(u0)
u0.assign(u1)
p0.assign(p1)
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
