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