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.

Thank you for your help!

#region: Load libraries----------------------------
from dolfin import *
import matplotlib.pyplot as plt
from ufl import nabla_div, max_value
from import Bar
import numpy as np
from tqdm import tqdm
import sys as system

#region: Choose script conditions------------------
#get rid of solving messages

#use chorin or IPCS?
chorin=True #if False, use IPCS

#set mesh

#region: Define saving files-----------------------
#get scripts name
name_of_self = system.argv[0].split('/')[-1].split('.')[0]

#for storing solutions

if chorin:
    xdmffile_velocity = XDMFFile(directory_xdmf+'chorin_velocity.xdmf') 
    xdmffile_pressure = XDMFFile(directory_xdmf+'chorin_pressure.xdmf')
    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

#region: Set parameters----------------------------
#set time
t=0.0; T = 5
if chorin:
    dt = 0.0256
    dt = 0.001

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))

bottom = Bottom()

boundaries = MeshFunction("size_t", mesh,mesh.topology().dim()-1,0)
top.mark(boundaries, 1)
leftwall.mark(boundaries, 3)
rightwall.mark(boundaries, 3)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

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

    # 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

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