Hi everybody!
I wouldn’t think that it would be needed, but I experience that my incompressible Navier-Stokes formulation needs to be solved with smaller step sizes (or finer mesh) when I use the incremental pressure correction scheme compared to Chorin’s splitting method.
-
Can this be correct or is there something wring with my formulation?
I have attached my code below, simply type true/false under “#use chorin or IPCS?”.
Note that the program right now uses different time steps for the two methods. -
Furthermore: Does anyone know what space “P” defines? I usually use Discontinuous Lagrange (DG), but the demo in the FEniCS manual uses “P”.
-
When to use which boundary conditions? The IPCS method in the FEniCS demo has boundary parts in the weak formulation of the momemtum equation, where the Chorin demo does not. I do not completely understand this? I have below added a boundary term to the Chorin formulation, but I do not know if this is correct.
Thank you for your help!
#region: Load libraries----------------------------
from dolfin import *
import matplotlib.pyplot as plt
from ufl import nabla_div, max_value
from progress.bar import Bar
import numpy as np
from tqdm import tqdm
import sys as system
#endregion----------------------------------------
#region: Choose script conditions------------------
#get rid of solving messages
set_log_active(False)
#use chorin or IPCS?
chorin=True #if False, use IPCS
#set mesh
generate_mesh=True
#endregion-----------------------------------------
#region: Define saving files-----------------------
#get scripts name
name_of_self = system.argv[0].split('/')[-1].split('.')[0]
#for storing solutions
directory_xdmf='results/'+name_of_self+'/'
if chorin:
xdmffile_velocity = XDMFFile(directory_xdmf+'chorin_velocity.xdmf')
xdmffile_pressure = XDMFFile(directory_xdmf+'chorin_pressure.xdmf')
else:
xdmffile_velocity = XDMFFile(directory_xdmf+'IPCS_velocity.xdmf')
xdmffile_pressure = XDMFFile(directory_xdmf+'IPCS_pressure.xdmf')
xdmffile_velocity.parameters["flush_output"] = True
xdmffile_pressure.parameters["flush_output"] = True
#endregion
#region: Set parameters----------------------------
#set time
t=0.0; T = 5
if chorin:
dt = 0.0256
else:
dt = 0.001
#parameters
k = Constant(dt)
force = Constant((0., 0.))
nu = Constant(0.01)
#endregion ----------------------------------------
#region: define mesh ------------------------------
# Create mesh
mesh = UnitSquareMesh(64, 64)
n = FacetNormal(mesh)
#endregion ----------------------------------------
#region: define boundaries
class Top(SubDomain):
def inside(self, x, on_boundary):
return (near(x[1], 1.0))
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return (near(x[1], 0.0))
class leftWall(SubDomain):
def inside(self, x, on_boundary):
return (near(x[0], 0.0))
class rightWall(SubDomain):
def inside(self, x, on_boundary):
return (near(x[0], 1.0))
top=Top()
bottom = Bottom()
leftwall=leftWall()
rightwall=rightWall()
boundaries = MeshFunction("size_t", mesh,mesh.topology().dim()-1,0)
top.mark(boundaries, 1)
bottom.mark(boundaries,2)
leftwall.mark(boundaries, 3)
rightwall.mark(boundaries, 3)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
#endregion
#region: define functions -------------------------
#define vector-vaued functionspace and scalar-valued function space
V = VectorFunctionSpace(mesh, "P", 2) #for velocity
Q = FunctionSpace(mesh, "P", 1) #for pressure
#define test and trial functions
v = TestFunction(V) #for velocity
u = TrialFunction(V) #for velocity
q = TestFunction(Q) #for pressure
p = TrialFunction(Q) #for pressure
# Create functions
u0 = Function(V) #for tentative velocity
u1 = Function(V) #for storing velocity solution
p0 = Function(Q) #for previous pressure(only IPCS)
p1 = Function(Q) #for pressure
#endregion ----------------------------------------
#region: define boundary conditions ---------------
pd=Constant(1)
#pd=Expression('1/(1+exp(-8*t+5))',t=0.0,degree=2)
noslip = DirichletBC(V, (0, 0), boundaries,3)
inflow = DirichletBC(Q, pd, boundaries,1)
outflow = DirichletBC(Q, 0, boundaries,2)
bcu = [noslip]
bcp = [inflow, outflow]
#endregion ----------------------------------------
#region: define weak formulation --
if chorin:
print("-Using Chorin's method")
# Tentative velocity step
F1 = (1/k)*inner(u - u0, v)*dx \
+ inner(grad(u0)*u, v)*dx \
+ nu*inner(grad(u), grad(v))*dx - inner(force, v)*dx\
+ nu*inner(grad(u)*n,v)*ds
a1 = lhs(F1)
L1 = rhs(F1)
# Pressure update
a2 = inner(grad(p), grad(q))*dx
L2 = -(1/k)*inner(div(u1),q)*dx
# Velocity update
a3 = inner(u, v)*dx
L3 = inner(u1, v)*dx - k*inner(grad(p1), v)*dx
else:
U = 0.5*(u0 + u)
print("-Using IPCS method")
# 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((u - u0) / k, v)*dx \
+ dot(dot(u0, nabla_grad(u0)), v)*dx \
+ inner(sigma(U, p0), epsilon(v))*dx \
+ dot(p0*n, v)*ds - dot(nu*nabla_grad(U)*n, v)*ds \
- dot(force, v)*dx
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/k)*div(u1)*q*dx
# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u1, v)*dx - k*dot(nabla_grad(p1 - p0), v)*dx
# 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/k)*div(u1)*q*dx
# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u1, v)*dx - k*dot(nabla_grad(p1 - p0), v)*dx
# Assemble matrices
A1 = assemble(a1); b1 = None
A2 = assemble(a2); b2 = None
A3 = assemble(a3); b3 = None
# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]
#endregion ----------------------------------------
#region: solve problem ----------------------------
#time step and update
bar = Bar('Processing', max=T/dt)
while t < T + DOLFIN_EPS:
pd.t = t+dt
dt = min(T-t, dt)
#update
bar.next()
# Step 1: Tentative velocity step
b1 = assemble(L1)
[bc.apply(b1) for bc in bcu]
solve(A1, u1.vector(), b1, 'bicgstab', 'hypre_amg')
# Step 2: Pressure correction step
b2 = assemble(L2)
[bc.apply(b2) for bc in bcp]
solve(A2, p1.vector(), b2, 'bicgstab', 'hypre_amg')
# Step 3: Velocity correction step
b3 = assemble(L3)
solve(A3, u1.vector(), b3, 'cg', 'sor')
#save solutions
xdmffile_velocity.write(u1, t); xdmffile_pressure.write(p1, t)
t += dt
# Update previous solution
u0.assign(u1)
p0.assign(p1)
bar.finish()
#endregion ----------------------------------------