IPCS solver with time-dependent pressure bc

Dear FEniCSx-community,

I am trying to adapt this flow-past-cylinder-demo to have time-dependent inflow/outflow boundary conditions. I implemented the time-dependence by InletPressure() which has a sinusoidal pressure fluctuation through time. I simplified the demo code by using the unit-square mesh instead.
The boundary conditions (see Figure 1) and the initial pressure (Figure 2) seem to be correct.

The problem is that after just a few time steps, there is no steady flow between inlet and outlet but most of the mesh cells have 0 pressure (see Figure 3). Also the velocities don’t look like expected (see Figure 4).

Figure 1:


Figure 2 (initial pressure):

Figure 3 (pressure at time frame 10):

Figure 4 (velocities at time frame 10):

Code:

import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.cpp.mesh import locate_entities_boundary, CellType
from dolfinx.fem import (Constant, Function, FunctionSpace, dirichletbc, form, locate_dofs_topological)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, create_vector, create_matrix, set_bc)
from dolfinx.io import VTXWriter
from dolfinx.mesh import create_unit_square

from ufl import (FiniteElement, TestFunction, TrialFunction, VectorElement,
                 div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs)

mesh = create_unit_square(MPI.COMM_WORLD, 50, 50, CellType.triangle)

t = 0
T = 1  # Final time
dt = 1 / 600  # Time step size
num_steps = int(T / dt)
k = Constant(mesh, PETSc.ScalarType(dt))
mu = Constant(mesh, PETSc.ScalarType(0.001))  # Dynamic viscosity
rho = Constant(mesh, PETSc.ScalarType(1))  # Density

v_cg2 = VectorElement("CG", mesh.ufl_cell(), 2)
s_cg1 = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, v_cg2)
Q = FunctionSpace(mesh, s_cg1)

gdim = 2
fdim = mesh.topology.dim - 1


# Define boundary conditions
class InletPressure():
    def __init__(self, t):
        self.t = t

    def __call__(self, x):
        return np.full(x.shape[1], 2 + np.sin(self.t))


u = TrialFunction(V)
v = TestFunction(V)
u_ = Function(V)
u_.name = "u"
u_s = Function(V)
u_n = Function(V)
u_n1 = Function(V)
p = TrialFunction(Q)
q = TestFunction(Q)
p_ = Function(Q)
p_.name = "p"
phi = Function(Q)


def inflow(x):
    return np.isclose(x[0], 0)


def outflow(x):
    return np.isclose(x[0], 1)


def walls(x):
    return np.logical_not(np.logical_or(inflow(x), outflow(x)))


# Walls
u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
facets = locate_entities_boundary(mesh, fdim, walls)
bcu_walls = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, facets), V)
bcu = [bcu_walls]

# Inlet
p_inlet = Function(Q)
inlet_pressure = InletPressure(t)
p_inlet.interpolate(inlet_pressure)
facets = locate_entities_boundary(mesh, fdim, inflow)
bcp_inflow = dirichletbc(p_inlet, locate_dofs_topological(Q, fdim, facets))

# Outlet
facets = locate_entities_boundary(mesh, fdim, outflow)
bcp_outlet = dirichletbc(PETSc.ScalarType(-2), locate_dofs_topological(Q, fdim, facets), Q)
bcp = [bcp_inflow, bcp_outlet]

u = TrialFunction(V)
v = TestFunction(V)
u_ = Function(V)
u_.name = "u"
u_s = Function(V)
u_n = Function(V)
u_n1 = Function(V)
p = TrialFunction(Q)
q = TestFunction(Q)
p_ = Function(Q)
p_.name = "p"
phi = Function(Q)

f = Constant(mesh, PETSc.ScalarType((0, 0)))
F1 = rho / k * dot(u - u_n, v) * dx
F1 += inner(dot(1.5 * u_n - 0.5 * u_n1, 0.5 * nabla_grad(u + u_n)), v) * dx
F1 += 0.5 * mu * inner(grad(u + u_n), grad(v)) * dx - dot(p_, div(v)) * dx
F1 += dot(f, v) * dx
a1 = form(lhs(F1))
L1 = form(rhs(F1))
A1 = create_matrix(a1)
b1 = create_vector(L1)

a2 = form(dot(grad(p), grad(q)) * dx)
L2 = form(-rho / k * dot(div(u_s), q) * dx)
A2 = assemble_matrix(a2, bcs=bcp)
A2.assemble()
b2 = create_vector(L2)

a3 = form(rho * dot(u, v) * dx)
L3 = form(rho * dot(u_s, v) * dx - k * dot(nabla_grad(phi), v) * dx)
A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)

# Solver for step 1
solver1 = PETSc.KSP().create(mesh.comm)
solver1.setOperators(A1)
solver1.setType(PETSc.KSP.Type.BCGS)
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.JACOBI)

# Solver for step 2
solver2 = PETSc.KSP().create(mesh.comm)
solver2.setOperators(A2)
solver2.setType(PETSc.KSP.Type.MINRES)
pc2 = solver2.getPC()
pc2.setType(PETSc.PC.Type.HYPRE)
pc2.setHYPREType("boomeramg")

# Solver for step 3
solver3 = PETSc.KSP().create(mesh.comm)
solver3.setOperators(A3)
solver3.setType(PETSc.KSP.Type.CG)
pc3 = solver3.getPC()
pc3.setType(PETSc.PC.Type.SOR)

vtx_u = VTXWriter(mesh.comm, "dfg2D-3-u.bp", [u_])
vtx_p = VTXWriter(mesh.comm, "dfg2D-3-p.bp", [p_])
vtx_u.write(t)
vtx_p.write(t)

for i in range(num_steps):
    # Update current time step
    t += dt
    # Update inlet velocity
    inlet_pressure.t = t
    p_inlet.interpolate(inlet_pressure)

    # Step 1: Tentative velocity step
    A1.zeroEntries()
    assemble_matrix(A1, a1, bcs=bcu)
    A1.assemble()
    with b1.localForm() as loc:
        loc.set(0)
    assemble_vector(b1, L1)
    apply_lifting(b1, [a1], [bcu])
    b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b1, bcu)
    solver1.solve(b1, u_s.vector)
    u_s.x.scatter_forward()

    # Step 2: Pressure correction step
    with b2.localForm() as loc:
        loc.set(0)
    assemble_vector(b2, L2)
    apply_lifting(b2, [a2], [bcp])
    b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b2, bcp)
    solver2.solve(b2, phi.vector)
    phi.x.scatter_forward()

    p_.vector.axpy(1, phi.vector)
    p_.x.scatter_forward()

    # Step 3: Velocity correction step
    with b3.localForm() as loc:
        loc.set(0)
    assemble_vector(b3, L3)
    b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    solver3.solve(b3, u_.vector)
    u_.x.scatter_forward()

    # Write solutions to file
    vtx_u.write(t)
    vtx_p.write(t)

    # Update variable with solution form this time step
    with u_.vector.localForm() as loc_, u_n.vector.localForm() as loc_n, u_n1.vector.localForm() as loc_n1:
        loc_n.copy(loc_n1)
        loc_.copy(loc_n)

vtx_u.close()
vtx_p.close()

What could be the cause of that?
And I am I correct that this flow-past-cylinder-demo uses the IPCS solver and can be used for more complex meshes?
Many thanks in advance for your help.

Best regards,
Georg

1 Like

I am facing a similar problem but also didn’t find a solution yet.

Using pressure conditions for the IPCS solver is not really well supported.

I’ve written out the mathematics of the splitting scheme at:
https://computationalphysiology.github.io/oasisx/splitting_schemes.html
highlighting issues with the boundary conditions known from literature (with references).

If you want to set DirichletBCs for the pressure condition, I would definitely not use the \phi -formulation of the pressure correction step, but solving it for p^(n+1) directly.

1 Like

@dokken Many thanks for your help. I would also be happy with the basic Chorin method in FEniCSx exactly like in this old fenics example.. Are pressure conditions also not well supported for Chorin’s method? Since IPCS is a modified version of Chorin’s method I probably don’t have to change much but I struggle to transition this old example to FEniCSx.

If you use the p^(n+1)-p^n formulation of either Chorin or IPCS, you can apply the pressure condition (although it is not natural in the sense of boundary conditions).

In your code you are applying the pressure boundary condition to the \phi. This means that
p^n+1 = p^n + bcp on any boundary with a dirichletbc, which is not correct.

What you can do is to solve for \phi, with zero Dirichlet conditions on this boundary.
Then after calling

you could assign bcp to p_.vector using dolfinx.fem.petsc.set_bc(p_.vector, bcp).

1 Like