CFD simulation in 3D cylinder fails to converge

[Jack757-cmd]

Hi,

I'm trying to run a minimal CFD simulation of steady laminar flow in a straight 3D cylindrical pipe using OasisX (v1.0.0), using dolfinx0.7.2. The case is adapted from simple_aorta_cfd_mesh.py but simplified to a single cylinder domain with a parabolic inlet velocity and zero pressure outlet.

However, the simulation fails to converge — the solver.solve() call returns a large residual difference (diff > 1e-6), and the run aborts around the first few time steps. I’ve tried adjusting the time step, viscosity, and mesh size, but the issue persists.

I’m wondering if there’s something wrong in how I’m defining boundary conditions, initializing the solver, or setting up the fractional step scheme.

import argparse
import logging
import sys
from typing import List

import numpy as np
from mpi4py import MPI
import ufl
from dolfinx import fem, io
from dolfinx.io import XDMFFile
from dolfinx import default_scalar_type
import oasisx
import os

# Logging
logging.basicConfig(
    level=logging.INFO,
    format="%(message)s",
    handlers=[logging.StreamHandler(sys.stdout)],
    force=True,
)
logger = logging.getLogger("Simple-Cylinder")
logger.setLevel(logging.INFO)

# Arguments
parser = argparse.ArgumentParser(description="Simple 3D cylinder CFD")
parser.add_argument("--D", type=float, default=3.0, help="Cylinder diameter (cm)")
parser.add_argument("--L", type=float, default=10.0, help="Cylinder length (cm)")
parser.add_argument("--mesh-size", type=float, default=0.3, help="Target mesh size (cm)")
parser.add_argument("--U0", type=float, default=10.0, help="Average inlet velocity (cm/s)")
parser.add_argument("--nu", type=float, default=0.035, help="Kinematic viscosity (cm²/s)")
parser.add_argument("--dt", type=float, default=0.01, help="Timestep (s)")
parser.add_argument("--T-end", type=float, default=2.0, help="End time (s)")
parser.add_argument("--u-deg", type=int, default=2, help="Velocity element degree")
parser.add_argument("--p-deg", type=int, default=1, help="Pressure element degree")
parser.add_argument("--save-interval", type=int, default=10, help="Save interval (steps)")
parser.add_argument("--log-interval", type=int, default=5, help="Log interval (steps)")
parser.add_argument("--max-iter", type=int, default=50, help="Max inner iterations")
parser.add_argument("--outdir", type=str, default="cylinder_out", help="Output directory")

args = parser.parse_args()

comm = MPI.COMM_WORLD
rank = comm.rank
os.makedirs(args.outdir, exist_ok=True)

def create_cylinder_mesh(D=3.0, L=10.0, mesh_size=0.3):
    if rank == 0:
        logger.info(f"Creating cylinder mesh: D={D}cm, L={L}cm, mesh_size={mesh_size}cm")
        import gmsh
        gmsh.initialize()
        gmsh.option.setNumber("General.Terminal", 1)
        gmsh.model.add("cylinder")

        R = D / 2.0
        # Inlet center at origin
        inlet_center = gmsh.model.occ.addPoint(0, 0, 0, mesh_size)
        outlet_center = gmsh.model.occ.addPoint(L, 0, 0, mesh_size)
        cylinder = gmsh.model.occ.addCylinder(0, 0, 0, L, 0, 0, R)

        # Identify surfaces
        gmsh.model.occ.synchronize()
        surfaces = gmsh.model.getEntities(dim=2)
        inlet_surfs, outlet_surfs, wall_surfs = [], [], []

        for dim, tag in surfaces:
            com = gmsh.model.occ.getCenterOfMass(2, tag)
            x = com[0]
            if abs(x - 0.0) < 0.2:  # Inlet at x=0
                inlet_surfs.append(tag)
            elif abs(x - L) < 0.2:  # Outlet at x=L
                outlet_surfs.append(tag)
            else:  # Wall
                wall_surfs.append(tag)

        # Physical groups
        gmsh.model.addPhysicalGroup(2, inlet_surfs, 1, name="inlet")
        gmsh.model.addPhysicalGroup(2, outlet_surfs, 2, name="outlet")
        gmsh.model.addPhysicalGroup(2, wall_surfs, 3, name="wall")
        gmsh.model.addPhysicalGroup(3, [cylinder], 1, name="domain")

        # Mesh settings
        gmsh.option.setNumber("Mesh.MeshSizeMin", mesh_size * 0.5)
        gmsh.option.setNumber("Mesh.MeshSizeMax", mesh_size * 1.5)
        gmsh.option.setNumber("Mesh.Algorithm3D", 10)  # HXT
        gmsh.model.mesh.generate(3)
        gmsh.model.mesh.optimize("Netgen")

    comm.barrier()
    from dolfinx.io import gmshio
    domain, cell_tags, facet_tags = gmshio.model_to_mesh(
        gmsh.model if rank == 0 else None, comm, rank=0, gdim=3
    )
    if rank == 0:
        gmsh.finalize()
    comm.barrier()

    if rank == 0:
        logger.info(f"Mesh created: {domain.topology.index_map(3).size_global} cells")
    return domain, cell_tags, facet_tags

# Create mesh
domain, cell_tags, facet_tags = create_cylinder_mesh(
    D=args.D, L=args.L, mesh_size=args.mesh_size
)
domain.name = "domain"
facet_tags.name = "facet_tags"

fdim = domain.topology.dim - 1
tdim = domain.topology.dim
domain.topology.create_entities(fdim)
domain.topology.create_connectivity(fdim, tdim)
gdim = domain.geometry.dim

# Boundary tags
INLET_TAG = 1
OUTLET_TAG = 2
WALL_TAG = 3

if rank == 0:
    logger.info(f"Boundary tags: Inlet={INLET_TAG}, Outlet={OUTLET_TAG}, Wall={WALL_TAG}")

# Compute inlet area and center
def compute_inlet_area_and_center(domain, facet_tags, inlet_tag, comm):
    fdim = domain.topology.dim - 1
    tdim = domain.topology.dim
    domain.topology.create_connectivity(fdim, 0)
    domain.topology.create_connectivity(fdim, tdim)
    inlet_facets = facet_tags.find(inlet_tag)
    
    local_area = 0.0
    local_center = np.zeros(gdim)
    if len(inlet_facets) > 0:
        conn = domain.topology.connectivity(fdim, 0)
        for f in inlet_facets:
            verts = conn.links(f)
            coords = domain.geometry.x[verts, :gdim]
            if len(verts) == 3:  # Triangle
                v0, v1, v2 = coords[0], coords[1], coords[2]
                area = 0.5 * np.linalg.norm(np.cross(v1 - v0, v2 - v0))
                centroid = (v0 + v1 + v2) / 3.0
            elif len(verts) == 4:  # Quad
                v0, v1, v2, v3 = coords[0], coords[1], coords[2], coords[3]
                area = 0.5 * (np.linalg.norm(np.cross(v2 - v0, v3 - v1)))
                centroid = (v0 + v1 + v2 + v3) / 4.0
            else:
                continue
            local_area += area
            local_center += centroid * area
    
    global_area = comm.allreduce(local_area, op=MPI.SUM)
    global_center = comm.allreduce(local_center, op=MPI.SUM)
    if global_area > 0:
        global_center /= global_area
    return global_area, global_center

inlet_area, inlet_center = compute_inlet_area_and_center(domain, facet_tags, INLET_TAG, comm)
inlet_radius = np.sqrt(inlet_area / np.pi)

if rank == 0:
    logger.info(f"Inlet area: {inlet_area:.4f} cm², radius: {inlet_radius:.4f} cm")
    logger.info(f"Inlet center: {inlet_center}")

# Parabolic inlet velocity
class ParabolicInlet:
    def __init__(self, U0: float, R: float, center: np.ndarray):
        self.U0 = U0  # Average velocity
        self.R = R    # Inlet radius
        self.center = center  # [y_center, z_center]

    def eval_component(self, x, comp: int):
        if comp != 0:  # Only non-zero in x-direction
            return np.zeros(x.shape[1], dtype=default_scalar_type)
        y = x[1] - self.center[0]
        z = x[2] - self.center[1]
        r2 = y**2 + z**2
        return 2 * self.U0 * (1.0 - r2 / (self.R ** 2))

inlet_profile = ParabolicInlet(args.U0, inlet_radius, inlet_center[1:])

def zero_func(x):
    return np.zeros(x.shape[1], dtype=default_scalar_type)

# Function spaces
el_u = ("Lagrange", args.u_deg)
el_p = ("Lagrange", args.p_deg)

# Boundary conditions
bcs_u_list = []
for comp in range(gdim):
    bc_inlet = oasisx.DirichletBC(
        lambda x, c=comp: inlet_profile.eval_component(x, c),
        oasisx.LocatorMethod.TOPOLOGICAL,
        (facet_tags, INLET_TAG)
    )
    bc_wall = oasisx.DirichletBC(
        zero_func,
        oasisx.LocatorMethod.TOPOLOGICAL,
        (facet_tags, WALL_TAG)
    )
    bcs_u_list.append([bc_inlet, bc_wall])

bcs_p = [oasisx.PressureBC(0.0, (facet_tags, OUTLET_TAG))]

# Solver options
options = {
    "low_memory_version": False,
}
solver_options = {
    "tentative": {"ksp_type": "cg", "pc_type": "hypre", "ksp_monitor": True},
    "pressure": {"ksp_type": "cg", "pc_type": "hypre", "ksp_monitor": True},
    "scalar": {"ksp_type": "cg", "pc_type": "hypre", "ksp_monitor": True},
}

if rank == 0:
    logger.info("Creating solver...")

solver = oasisx.FractionalStep_AB_CN(
    domain, el_u, el_p, bcs_u=bcs_u_list, bcs_p=bcs_p,
    solver_options=solver_options, options=options, body_force=None,
)

# Initialize from rest
for i in range(gdim):
    solver._u1[i].x.array[:] = 0.0
    solver._u2[i].x.array[:] = 0.0
solver._p.x.array[:] = 0.0

# Output writers
try:
    geom_deg = domain.geometry.cmap.degree
except AttributeError:
    geom_deg = domain.geometry.cmaps[0].degree

Ve = ufl.VectorElement("Lagrange", domain.ufl_cell(), geom_deg)
Pe = ufl.FiniteElement("Lagrange", domain.ufl_cell(), geom_deg)
V_out = fem.FunctionSpace(domain, Ve)
Q_out = fem.FunctionSpace(domain, Pe)

u_out = fem.Function(V_out); u_out.name = "u"
p_out = fem.Function(Q_out); p_out.name = "p"

u_file = XDMFFile(comm, os.path.join(args.outdir, "u_velocity.xdmf"), "w")
u_file.write_mesh(domain)
p_file = XDMFFile(comm, os.path.join(args.outdir, "p_pressure.xdmf"), "w")
p_file.write_mesh(domain)

# Time stepping
dt = args.dt
nu = args.nu
num_steps = int(np.ceil(args.T_end / dt))
Re = args.U0 * args.D / nu

if rank == 0:
    logger.info(f"Starting simulation: steps={num_steps}, dt={dt}s, Re={Re:.1f}")

# Initial output
u_out.interpolate(solver.u)
p_out.interpolate(solver._p)
u_file.write_function(u_out, 0.0)
p_file.write_function(p_out, 0.0)

# Time loop
t = 0.0
for step in range(1, num_steps + 1):
    t += dt
    diff = solver.solve(dt, nu, max_iter=args.max_iter)
    
    if diff > 1e-6:
        if rank == 0:
            logger.error(f"Failed to converge at step {step}, t={t:.4f}, diff={diff:.3e}")
        sys.exit(1)
    
    if step % args.save_interval == 0 or step == num_steps:
        u_out.interpolate(solver.u)
        p_out.interpolate(solver._p)
        u_file.write_function(u_out, t)
        p_file.write_function(p_out, t)
        if rank == 0:
            logger.info(f"Saved at step {step}, t={t:.4f}")
    
    if step % args.log_interval == 0 or step == num_steps:
        ux_local = solver._u1[0].x.array

[jorgensd]

I wouldn't use the diff as a way of identifying the convergence of the method.
Especially not on a coarse grid, where the initial flow field is 0.

I attach an adapted code (for DOLFINx v0.10) that doesn't terminate, but rather writes the output to file (and does wait for inner iterations to converge):

import argparse
import logging
import sys
from typing import List

import numpy as np
from mpi4py import MPI
import ufl
from dolfinx import fem, io
from dolfinx.io import XDMFFile
from dolfinx import default_scalar_type
import oasisx
import os
import basix.ufl

# Logging
logging.basicConfig(
    level=logging.INFO,
    format="%(message)s",
    handlers=[logging.StreamHandler(sys.stdout)],
    force=True,
)
logger = logging.getLogger("Simple-Cylinder")
logger.setLevel(logging.INFO)

# Arguments
parser = argparse.ArgumentParser(description="Simple 3D cylinder CFD")
parser.add_argument("--D", type=float, default=3.0, help="Cylinder diameter (cm)")
parser.add_argument("--L", type=float, default=10.0, help="Cylinder length (cm)")
parser.add_argument("--mesh-size", type=float, default=0.3, help="Target mesh size (cm)")
parser.add_argument("--U0", type=float, default=10.0, help="Average inlet velocity (cm/s)")
parser.add_argument("--nu", type=float, default=0.035, help="Kinematic viscosity (cm²/s)")
parser.add_argument("--dt", type=float, default=0.01, help="Timestep (s)")
parser.add_argument("--T-end", type=float, default=2.0, help="End time (s)")
parser.add_argument("--u-deg", type=int, default=2, help="Velocity element degree")
parser.add_argument("--p-deg", type=int, default=1, help="Pressure element degree")
parser.add_argument("--save-interval", type=int, default=10, help="Save interval (steps)")
parser.add_argument("--log-interval", type=int, default=5, help="Log interval (steps)")
parser.add_argument("--max-iter", type=int, default=50, help="Max inner iterations")
parser.add_argument("--outdir", type=str, default="cylinder_out", help="Output directory")

args = parser.parse_args()

comm = MPI.COMM_WORLD
rank = comm.rank
os.makedirs(args.outdir, exist_ok=True)

def create_cylinder_mesh(D=3.0, L=10.0, mesh_size=0.3):
    if rank == 0:
        logger.info(f"Creating cylinder mesh: D={D}cm, L={L}cm, mesh_size={mesh_size}cm")
        import gmsh
        gmsh.initialize()
        gmsh.option.setNumber("General.Terminal", 1)
        gmsh.model.add("cylinder")

        R = D / 2.0
        # Inlet center at origin
        cylinder = gmsh.model.occ.addCylinder(0, 0, 0, L, 0, 0, R)

        # Identify surfaces
        gmsh.model.occ.synchronize()
        surfaces = gmsh.model.getEntities(dim=2)
        inlet_surfs, outlet_surfs, wall_surfs = [], [], []

        for dim, tag in surfaces:
            com = gmsh.model.occ.getCenterOfMass(2, tag)
            x = com[0]
            if abs(x - 0.0) < 0.2:  # Inlet at x=0
                inlet_surfs.append(tag)
            elif abs(x - L) < 0.2:  # Outlet at x=L
                outlet_surfs.append(tag)
            else:  # Wall
                wall_surfs.append(tag)

        # Physical groups
        gmsh.model.addPhysicalGroup(2, inlet_surfs, 1, name="inlet")
        gmsh.model.addPhysicalGroup(2, outlet_surfs, 2, name="outlet")
        gmsh.model.addPhysicalGroup(2, wall_surfs, 3, name="wall")
        gmsh.model.addPhysicalGroup(3, [cylinder], 1, name="domain")

        # Mesh settings
        gmsh.option.setNumber("Mesh.MeshSizeMin", mesh_size * 0.5)
        gmsh.option.setNumber("Mesh.MeshSizeMax", mesh_size * 1.5)
        gmsh.option.setNumber("Mesh.Algorithm3D", 10)  # HXT
        gmsh.model.mesh.generate(3)
        gmsh.model.mesh.optimize("Netgen")

    comm.barrier()
    from dolfinx.io import gmsh as gmshio
    model = gmshio.model_to_mesh(
        gmsh.model if rank == 0 else None, comm, rank=0, gdim=3
    )
    domain, cell_tags, facet_tags  = model.mesh, model.cell_tags, model.facet_tags
    if rank == 0:
        gmsh.finalize()
    comm.barrier()

    if rank == 0:
        logger.info(f"Mesh created: {domain.topology.index_map(3).size_global} cells")
    return domain, cell_tags, facet_tags

# Create mesh
domain, cell_tags, facet_tags = create_cylinder_mesh(
    D=args.D, L=args.L, mesh_size=args.mesh_size
)
domain.name = "domain"
facet_tags.name = "facet_tags"

fdim = domain.topology.dim - 1
tdim = domain.topology.dim
domain.topology.create_entities(fdim)
domain.topology.create_connectivity(fdim, tdim)
gdim = domain.geometry.dim

# Boundary tags
INLET_TAG = 1
OUTLET_TAG = 2
WALL_TAG = 3

if rank == 0:
    logger.info(f"Boundary tags: Inlet={INLET_TAG}, Outlet={OUTLET_TAG}, Wall={WALL_TAG}")

# Compute inlet area and center
def compute_inlet_area_and_center(domain, facet_tags, inlet_tag, comm):
    fdim = domain.topology.dim - 1
    tdim = domain.topology.dim
    domain.topology.create_connectivity(fdim, 0)
    domain.topology.create_connectivity(fdim, tdim)
    inlet_facets = facet_tags.find(inlet_tag)
    
    local_area = 0.0
    local_center = np.zeros(gdim)
    if len(inlet_facets) > 0:
        conn = domain.topology.connectivity(fdim, 0)
        for f in inlet_facets:
            verts = conn.links(f)
            coords = domain.geometry.x[verts, :gdim]
            if len(verts) == 3:  # Triangle
                v0, v1, v2 = coords[0], coords[1], coords[2]
                area = 0.5 * np.linalg.norm(np.cross(v1 - v0, v2 - v0))
                centroid = (v0 + v1 + v2) / 3.0
            elif len(verts) == 4:  # Quad
                v0, v1, v2, v3 = coords[0], coords[1], coords[2], coords[3]
                area = 0.5 * (np.linalg.norm(np.cross(v2 - v0, v3 - v1)))
                centroid = (v0 + v1 + v2 + v3) / 4.0
            else:
                continue
            local_area += area
            local_center += centroid * area
    
    global_area = comm.allreduce(local_area, op=MPI.SUM)
    global_center = comm.allreduce(local_center, op=MPI.SUM)
    if global_area > 0:
        global_center /= global_area
    return global_area, global_center

inlet_area, inlet_center = compute_inlet_area_and_center(domain, facet_tags, INLET_TAG, comm)
inlet_radius = np.sqrt(inlet_area / np.pi)

if rank == 0:
    logger.info(f"Inlet area: {inlet_area:.4f} cm², radius: {inlet_radius:.4f} cm")
    logger.info(f"Inlet center: {inlet_center}")

# Parabolic inlet velocity
class ParabolicInlet:
    def __init__(self, U0: float, R: float, center: np.ndarray):
        self.U0 = U0  # Average velocity
        self.R = R    # Inlet radius
        self.center = center  # [y_center, z_center]

    def eval_component(self, x, comp: int):
        if comp != 0:  # Only non-zero in x-direction
            return np.zeros(x.shape[1], dtype=default_scalar_type)
        y = x[1] - self.center[0]
        z = x[2] - self.center[1]
        r2 = y**2 + z**2
        return 2 * self.U0 * (1.0 - r2 / (self.R ** 2))

inlet_profile = ParabolicInlet(args.U0, inlet_radius, inlet_center[1:])

def zero_func(x):
    return np.zeros(x.shape[1], dtype=default_scalar_type)

# Function spaces
el_u = ("Lagrange", args.u_deg)
el_p = ("Lagrange", args.p_deg)

# Boundary conditions
bcs_u_list = []
for comp in range(gdim):
    bc_inlet = oasisx.DirichletBC(
        lambda x, c=comp: inlet_profile.eval_component(x, c),
        oasisx.LocatorMethod.TOPOLOGICAL,
        (facet_tags, INLET_TAG)
    )
    bc_wall = oasisx.DirichletBC(
        zero_func,
        oasisx.LocatorMethod.TOPOLOGICAL,
        (facet_tags, WALL_TAG)
    )
    bcs_u_list.append([bc_inlet, bc_wall])

bcs_p = [oasisx.PressureBC(0.0, (facet_tags, OUTLET_TAG))]

# Solver options
options = {
    "low_memory_version": False,
}
solver_options = {
    "tentative": {"ksp_type": "preonly", "pc_type": "lu", "ksp_monitor": None,
               "pc_factor_mat_solver_type": "mumps", "ksp_error_if_not_converged": True},
    "pressure": {"ksp_type": "preonly", "pc_type": "lu", "ksp_monitor": None,
               "pc_factor_mat_solver_type": "mumps", "ksp_error_if_not_converged": True},
    "scalar": {"ksp_type": "preonly", "pc_type": "lu", "ksp_monitor": None,
               "pc_factor_mat_solver_type": "mumps", "ksp_error_if_not_converged": True},
}

if rank == 0:
    logger.info("Creating solver...")

solver = oasisx.FractionalStep_AB_CN(
    domain, el_u, el_p, bcs_u=bcs_u_list, bcs_p=bcs_p,
    solver_options=solver_options, options=options, body_force=None,
)

# Initialize from rest
for i in range(gdim):
    solver._u1[i].x.array[:] = 0.0
    solver._u2[i].x.array[:] = 0.0
solver._p.x.array[:] = 0.0

# Output writers
try:
    geom_deg = domain.geometry.cmap.degree
except AttributeError:
    geom_deg = domain.geometry.cmaps[0].degree

Ve = basix.ufl.element("Lagrange", domain.basix_cell(), geom_deg, shape=(domain.geometry.dim,))
Pe = basix.ufl.element("Lagrange", domain.basix_cell(), geom_deg)
V_out = fem.functionspace(domain, Ve)
Q_out = fem.functionspace(domain, Pe)

u_out = fem.Function(V_out); u_out.name = "u"
p_out = fem.Function(Q_out); p_out.name = "p"

u_file = XDMFFile(comm, os.path.join(args.outdir, "u_velocity.xdmf"), "w")
u_file.write_mesh(domain)
p_file = XDMFFile(comm, os.path.join(args.outdir, "p_pressure.xdmf"), "w")
p_file.write_mesh(domain)

# Time stepping
dt = args.dt
nu = args.nu
num_steps = int(np.ceil(args.T_end / dt))
Re = args.U0 * args.D / nu

if rank == 0:
    logger.info(f"Starting simulation: steps={num_steps}, dt={dt}s, Re={Re:.1f}")

# Initial output
u_out.interpolate(solver.u)
p_out.interpolate(solver._p)
u_file.write_function(u_out, 0.0)
p_file.write_function(p_out, 0.0)

# Time loop
t = 0.0
inner_tol = 1e-3
for step in range(1, num_steps + 1):
    t += dt
    diff = solver.solve(dt, nu, max_error=inner_tol, max_iter=args.max_iter)

    
    # if step % args.save_interval == 0 or step == num_steps:
    u_out.interpolate(solver.u)
    p_out.interpolate(solver._p)
    u_file.write_function(u_out, t)
    p_file.write_function(p_out, t)

    u_file.flush()
    p_file.flush()

    
    if rank == 0:
        logger.info(f"Saved at step {step}, t={t:.4f}")
    
    if step % args.log_interval == 0 or step == num_steps:
        ux_local = solver._u1[0].x.array

u_file.close()
p_file.close()

and executed with:

mpirun -n 3 python3 mwe.py  --dt=0.05 --mesh-size=0.08 --max-iter=2

You'll get a velocity profile that looks like