 # Does IPCS require smaller steps than Chorin?

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.

1. 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.

2. Furthermore: Does anyone know what space “P” defines? I usually use Discontinuous Lagrange (DG), but the demo in the FEniCS manual uses “P”.

3. 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.

``````#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.split('/')[-1].split('.')

#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.0))

class Bottom(SubDomain):
def inside(self, x, on_boundary):
return (near(x, 0.0))

class leftWall(SubDomain):
def inside(self, x, on_boundary):
return (near(x, 0.0))

class rightWall(SubDomain):
def inside(self, x, on_boundary):
return (near(x, 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 \
a1 = lhs(F1)
L1 = rhs(F1)

# Pressure update
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")
def epsilon(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 \
+ 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

# 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

# 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 ----------------------------------------``````

Please note that there are many IPCS methods, depending on how you do temporal discretization and the linearization of the non linear term. See for instance: https://www.sciencedirect.com/science/article/pii/S0045782520303145