Dear Fenics community,
I am solving for an FSI problem setting in FEniCS wherein I am trying to implement SUPG stabilization scheme for the navier stokes equation. As of now I am using the projection-method (Chorins method) or more precisely the incremental pressure correction method (IPCS) suggested in the fenics tutorial and this forces me to use Taylor hood elements. Also I am new to stabilization schemes and haven’t really done them before.
I have already gone through the following tutorials or notes and have got a decent idea of how these stabilisation schemes work:
- https://david-kamensky.eng.ucsd.edu/teaching/mae-207-fea-for-coupled-problems
- https://fenicsproject.org/docs/dolfin/1.3.0/python/demo/documented/stokes-stabilized/python/documentation.html
This brings me to three fundamental questions:
Question1 : Most of these tutorials add the NS and continuity together - add SUPG stabilization in this setting - and then solve using the newtons solver. However I am using the pressure correction scheme, so does SUPG stabilization also work for the pressure correction scheme? If so how exactly shall I proceed doing that? (I have attached a copy of my working non-dimensional code below - including just the fluid equation as a reference)
from dolfin import *
from mshr import *
import time
T = 3 # Total time
dt = 0.01 # Time step
# Non-dimensional numbers
Re = 75 # Reynolds number
Fr = 10000 # Froude number
mesh1 = RectangleMesh(Point(0.0, 0.0), Point(2.5, 1.0), 100, 40)
pv1 = File("mesh1.pvd")
pv1 << mesh1
# MPI-rank
mpi_comm = mpi_comm_world()
MPI.barrier(mpi_comm)
my_rank = MPI.rank(mpi_comm)
# Define function spaces
V = VectorFunctionSpace(mesh1, 'P', 2) # Fluid velocity on mesh1
Q = FunctionSpace(mesh1, 'P', 1) # Fluid pressure and temperature on mesh1
# Define boundaries for combined domain
bottom = 'near(x[1], 0)'
top = 'near(x[1], 1)'
left = 'near(x[0], 0)'
right = 'near(x[0], 2.5)'
# Define boundary conditions for combined domain
bcu_left = DirichletBC(V, Constant((1, 0)), left)
bcp_right = DirichletBC(Q, Constant(0), right)
bcu_bottom = DirichletBC(V, Constant((0, 0)), bottom)
bcu_top = DirichletBC(V, Constant((0, 0)), top)
bcu = [bcu_bottom, bcu_top, bcu_left]
bcp = [bcp_right]
# Define trial and test functions
u1 = TrialFunction(V) # Trial function for fluid velocity on mesh1
v = TestFunction(V) # Test function for fluid veloctiy on mesh1
p = TrialFunction(Q) # Trial function for fluid pressure on mesh1
q = TestFunction(Q) # Test function for fluid pressure on mesh1
# Define functions
u_n = Function(V)
us_ = Function(V)
p_n = Function(Q)
p_ = Function(Q)
# Defining integration domains
dx1 = dolfin.dx(mesh1)
ds1 = dolfin.ds(mesh1)
# Define expressions
U = 0.5*(u_n + u1)
n = FacetNormal(mesh1)
tsp = dt
dt = Constant(dt)
Re = Constant(Re)
# Let us treat this a lagrange muliplier term .............
# I have kept this so that I can relate it to my actual local lagrange multiplier term
# which arises due to solid motion in FSI. I think I have to apply SUPG stabilisation to this term as well, right?
f = Constant((0, -1/(Fr*Fr)))
# Define symmetric gradient
def epsilon(u):
return sym(nabla_grad(u))
# Define stress tensor
def sigma(u, p):
return (2/Re)*epsilon(u) - p*Identity(len(u))
# Define variational problem for step 1
F1 = dot((u1 - u_n) / dt, v)*dx1 + dot(dot(u_n, nabla_grad(u_n)), v)*dx1 + inner(sigma(U, p_n), epsilon(v))*dx1 + dot(p_n*n, v)*ds1 - dot((1/Re)*nabla_grad(U)*n, v)*ds1 - dot(f, v)*dx1
a1 = lhs(F1)
L1 = rhs(F1)
# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx1
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx1 - (1/dt)*div(us_)*q*dx1
# Define variational problem for step 3
a3 = dot(u1, v)*dx1
L3 = dot(us_, v)*dx1 - dt*dot(nabla_grad(p_ - p_n), v)*dx1
# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
# Time
t = 0
c = 0
# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]
file_u = File("us.pvd")
file_p = File("p.pvd")
# Create progress bar
progress = Progress('Time-stepping')
set_log_level(PROGRESS)
start_time = time.time()
while t < T:
# Update current time
t += tsp
c += 1
if my_rank == 0:
print "Step 1: Tentative velocity step"
b1 = assemble(L1)
[bc.apply(b1) for bc in bcu]
solve(A1, us_.vector(), b1, 'bicgstab', 'hypre_amg')
if my_rank == 0:
print "Step 2: Pressure correction step"
b2 = assemble(L2)
[bc.apply(b2) for bc in bcp]
solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')
if my_rank == 0:
print "Step 3: Velocity correction step"
b3 = assemble(L3)
solve(A3, us_.vector(), b3, 'bicgstab', 'sor')
# Print output files
if c > 20:
c = 0
file_u << us_
file_p << p_
# Update progress bar
if my_rank == 0:
progress.update(t / T)
# Update previous solution
u_n.assign(us_)
p_n.assign(p_)
end_time = time.time()
print(" total time", end_time - start_time)
Question2 : I have noticed that Oasis also uses the pressure correction scheme and all the elements for velocity and pressure are linear elements? This paper by Mortensen explains the basics of the fast solver and they do not discuss the stabilization schemes in the paper. So how does the optimized solver function well with P1P1 pair of elements?
Questions3 : Does using such a SUPG stabilisation (for NS eq…) for my FSI setting using the DLM/FD method come at a cost of reduced accuracy and loss in convergence?
I know these might be trivial questions, but can someone provide more insight and help me regarding these? Thanks in advance.