Odd convergence properties for IPCS scheme with N-S lid driven cavity flow

General
1

Hi!

I have been toying around with Jorgen’s IPCS schemes (from his tutorial) for the Navier-Stokes equation and was evaluating their performance on the lid-driven cavity flow benchmark (see https://doi.org/10.1115/1.4052149). While I am able to obtain qualitatively accurate results that look fine, the convergence seems quite suspect:

ipcs_a_liddriven
ipcs_a_liddriven3600×2400 396 KB

The solver code I have implemented is transient, so I tested the simulation time to check if I have actually reached steady state and this was not the source of the issue. I have also compared the simulation results to other benchmarks in the literature and observed similar behavior. I also used the 2D CFL number to set the dt for the simulation with a CFL of 0.5.

Would appreciate any insight that anyone may have!

"""
Lid‐driven cavity IPCS in DOLFINx, with center‐line sampling of u and p
moved inside the solver function.  Uses Marchi et al, 2021 Tables 20 & 21
for reference.
"""

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path

import dolfinx
from dolfinx.fem import (Constant, Function, functionspace,
                         dirichletbc, locate_dofs_geometrical,
                         form, Expression)
from dolfinx.fem.petsc import (assemble_matrix, assemble_vector,
                               apply_lifting, create_vector, set_bc)
from dolfinx.mesh import create_unit_square
from basix.ufl import element
from ufl import (FacetNormal, Identity, TestFunction, TrialFunction,
                 div, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym)

from dolfinx import geometry

#formatting
plt.rcParams["text.usetex"] = True
plt.rc('font', family='serif')

# ————— Reference data from Tables 9 & 10 —————
y_ref   = np.array([0.0546875, 0.0625, 0.0703125, 0.1015625,0.125,0.171875,0.1875,
                    0.25,0.28125,0.3125,0.375,0.4375,0.453125,0.5,0.5625,0.6171875,
                    0.625,0.6875,0.734375,0.75,0.8125,0.8515625,0.875,
                    0.9375,0.953125,0.9609375,0.96875,0.9765625])

u_ref = np.array([-1.8125321e-1,-2.023291e-1,-2.229271e-1,-3.003688e-1,-3.478431e-1,-3.885676e-1,-3.844078e-1,
                  -3.1894525e-1,-2.8042785e-1,-2.4569306e-1,-1.8373161e-1,-1.2341016e-1, -1.0817529e-1,-6.205605e-2,5.6167e-4,
                  5.700436e-2,6.524839e-2,1.3357199e-1,1.8864358e-1,2.0791378e-1,2.884412e-1,3.371763e-1,3.625438e-1,
                  4.229307e-1,4.724491e-1,5.171830e-1,5.803655e-1,6.6397203e-1])



x_ref = np.array([
    0.0625,   0.0703125, 0.078125,   0.09375,
    0.125,    0.15625,   0.1875,    0.2265625,
    0.234375, 0.25,      0.3125,    0.375,
    0.4375,   0.5,       0.5625,    0.625,
    0.6875,   0.75,      0.8046875, 0.8125,
    0.859375, 0.875,     0.90625,   0.9375,
    0.9453125,0.953125,  0.9609375, 0.96875
])

v_ref = np.array([
    2.807042e-1,   2.962922e-1,   3.099493e-1,   3.329768e-1,
    3.650400e-1,   3.769157e-1,   3.678512e-1,   3.340319e-1,
    3.253866e-1,   3.0710340e-1,  2.3126780e-1,  1.6056381e-1,
    0.9296911e-1,  2.5799473e-2,  -4.184046e-2,   -1.1079783e-1,
    -1.8167907e-1,  -2.533806e-1,   -3.201954e-1,   -3.31565806e-1,
    -4.263892e-1,   -4.677741e-1,   -5.264155e-1,   -4.561505e-1,
    -4.102923e-1,   -3.551312e-1,   -2.933786e-1,   -2.283407e-1
])




def sample_function_at_points(mesh, func, phys_points):
    # build bounding‐box tree
    bb = geometry.bb_tree(mesh, mesh.topology.dim)
    # find candidate cells for each point (points as 3×N in tutorial, so transpose)
    candidates = geometry.compute_collisions_points(bb, phys_points)
    # filter to actual containing cells
    colliding  = geometry.compute_colliding_cells(mesh, candidates, phys_points)

    pts_on_proc, cells, indices = [], [], []
    for i, p in enumerate(phys_points):
        links = colliding.links(i)
        if len(links) > 0:
            pts_on_proc.append(p)
            cells.append(links[0])
            indices.append(i)

    if not pts_on_proc:
        # no sample on this rank
        return np.full(len(phys_points), np.nan)

    pts_on_proc = np.array(pts_on_proc, dtype=np.float64)
    cells       = np.array(cells, dtype=np.int32)
    # evaluate your FE function
    vals = func.eval(pts_on_proc, cells)

    # flatten scalar result
    if vals.ndim == 2 and vals.shape[1] == 1:
        vals = vals[:, 0]

    # reassemble full‐length array
    if vals.ndim == 1:
        full = np.full(len(phys_points), np.nan)
        for idx, v in zip(indices, vals):
            full[idx] = v
    else:
        # vector case
        dim = vals.shape[1]
        full = np.full((len(phys_points), dim), np.nan)
        for idx, v in zip(indices, vals):
            full[idx, :] = v

    return full

def solve_lid_cavity(Nx, T, Re):
    #Set up time stepping
    Re      = Re
    CFL     = 0.5
    T       = T
    stride  = 50
    nu_val  = 1.0 / Re
    Nx=Ny = Nx

    #Set up mesh
    domain = create_unit_square(MPI.COMM_WORLD, Nx, Ny)
    v_cg2 = element("Lagrange", domain.topology.cell_name(), 2, shape=(domain.geometry.dim, ))
    s_cg1 = element("Lagrange", domain.topology.cell_name(), 1)
    V = functionspace(domain, v_cg2)
    Q = functionspace(domain, s_cg1)

    coords = domain.geometry.x
    xs = np.unique(coords[:, 0]); ys = np.unique(coords[:, 1])
    dxs = float(np.min(np.diff(np.sort(xs))))
    dy = float(np.min(np.diff(np.sort(ys))))

    dt_value = CFL * (dxs * dy) / (dxs + dy)
    dt = PETSc.ScalarType(dt_value)

    dt_const = Constant(domain, dt)
    nu_const = Constant(domain, PETSc.ScalarType(nu_val))

    u = TrialFunction(V)
    v = TestFunction(V)
    p = TrialFunction(Q)
    q = TestFunction(Q)

    def top(x): return np.isclose(x[1], 1.0)

    def walls(x):
        return np.logical_or(
            np.isclose(x[1], 0.0),
            np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0))
        )

    wall_dofs = locate_dofs_geometrical(V, walls)
    u_noslip = np.array((0.0,)*domain.geometry.dim, dtype=PETSc.ScalarType)
    bc_walls = dirichletbc(u_noslip, wall_dofs, V)

    lid_dofs = locate_dofs_geometrical(V, top)
    u_lid   = np.array((1.0, 0.0), dtype=PETSc.ScalarType)
    bc_lid   = dirichletbc(u_lid, lid_dofs, V)

    bcs_u = [bc_lid, bc_walls]

    p_dofs   = locate_dofs_geometrical(
        Q, lambda x: np.logical_and(np.isclose(x[0], 0.0),
                                     np.isclose(x[1], 0.0))
    )
    bc_p     = [dirichletbc(Constant(domain, 0.0),p_dofs, Q)]

    u_n = Function(V)
    u_  = Function(V)
    p_n = Function(Q)
    p_  = Function(Q)
    n = FacetNormal(domain)
    U = 0.5*(u_n + u)

    f = Constant(domain, PETSc.ScalarType((0, 0)))

    def epsilon(u):
        return sym(nabla_grad(u))
    def sigma(u, p):
        return 2 * nu_const * epsilon(u) - p * Identity(len(u))

    F1 = dot((u - u_n) / dt_const, v) * dx
    F1 += dot(dot(u_n, nabla_grad(u_n)), v) * dx
    F1 += inner(sigma(U, p_n), epsilon(v)) * dx
    F1 -= dot(f, v) * dx

    a1 = form(lhs(F1))
    L1 = form(rhs(F1))

    A1 = assemble_matrix(a1, bcs=bcs_u)
    A1.assemble()
    b1 = create_vector(L1)

    F2 = (dot(nabla_grad(p - p_n), nabla_grad(q))*dx
        +(1/dt_const)*div(u_)*q*dx
        )

    a2, L2 = form(lhs(F2)), form(rhs(F2))
    A2 = assemble_matrix(a2, bcs=bc_p)
    A2.assemble()
    b2 = create_vector(L2)

    F3 = dot(u, v)*dx + dt_const*dot(nabla_grad(p_ - p_n), v)*dx
    F3 -= dot(u_, v)*dx
    a3, L3 = form(lhs(F3)), form(rhs(F3))
    A3 = assemble_matrix(a3)
    A3.assemble()
    b3 = create_vector(L3)

    solver1 = PETSc.KSP().create(domain.comm)
    solver1.setOperators(A1)
    solver1.setType(PETSc.KSP.Type.BCGS)
    pc1 = solver1.getPC()
    pc1.setType(PETSc.PC.Type.HYPRE)
    pc1.setHYPREType("boomeramg")

    solver2 = PETSc.KSP().create(domain.comm)
    solver2.setOperators(A2)
    solver2.setType(PETSc.KSP.Type.BCGS)
    pc2 = solver2.getPC()
    pc2.setType(PETSc.PC.Type.HYPRE)
    pc2.setHYPREType("boomeramg")

    solver3 = PETSc.KSP().create(domain.comm)
    solver3.setOperators(A3)
    solver3.setType(PETSc.KSP.Type.CG)
    pc3 = solver3.getPC()
    pc3.setType(PETSc.PC.Type.SOR)

    u_n.interpolate(lambda x: (np.zeros_like(x[0]), np.zeros_like(x[0])))
    p_n.interpolate(lambda x: np.zeros_like(x[0]))

    step = 0
    t = 0.0

    while t < T - 1e-6:
        t += float(dt)
        print(t)

        with b1.localForm() as loc_1:
            loc_1.set(0)
        assemble_vector(b1, L1)
        apply_lifting(b1, [a1], [bcs_u])
        b1.ghostUpdate(addv=PETSc.InsertMode.ADD,
                       mode=PETSc.ScatterMode.REVERSE)
        set_bc(b1, bcs_u)
        solver1.solve(b1, u_.x.petsc_vec)
        u_.x.scatter_forward()

        with b2.localForm() as loc_2:
            loc_2.set(0)
        assemble_vector(b2, L2)
        apply_lifting(b2, [a2], [bc_p])
        b2.ghostUpdate(addv=PETSc.InsertMode.ADD,
                       mode=PETSc.ScatterMode.REVERSE)
        set_bc(b2, bc_p)
        solver2.solve(b2, p_.x.petsc_vec)
        p_.x.scatter_forward()

        with b3.localForm() as loc_3:
            loc_3.set(0)
        assemble_vector(b3, L3)
        b3.ghostUpdate(addv=PETSc.InsertMode.ADD,
                       mode=PETSc.ScatterMode.REVERSE)
        solver3.solve(b3, u_.x.petsc_vec)
        u_.x.scatter_forward()

        u_n.x.array[:] = u_.x.array[:]
        p_n.x.array[:] = p_.x.array[:]

        step += 1

    b1.destroy()
    b2.destroy()
    b3.destroy()
    solver1.destroy()
    solver2.destroy()
    solver3.destroy()

    # ————— sample at final time —————
    mesh = u_n.function_space.mesh
    # coarse points for error
    if MPI.COMM_WORLD.rank == 0:
        phys_u     = np.column_stack([0.5*np.ones_like(y_ref), y_ref,     np.zeros_like(y_ref)])
        phys_v     = np.column_stack([x_ref,               0.5*np.ones_like(x_ref), np.zeros_like(x_ref)])
        # fine points for plotting
        y_plot     = np.linspace(0.0, 1.0, 100)
        x_plot     = np.linspace(0.0, 1.0, 100)
        phys_u_f   = np.column_stack([0.5*np.ones_like(y_plot), y_plot, np.zeros_like(y_plot)])
        phys_v_f   = np.column_stack([x_plot, 0.5*np.ones_like(x_plot), np.zeros_like(x_plot)])
    else:
        phys_u = phys_v = phys_u_f = phys_v_f = np.empty((0,3), dtype=np.float64)

    # coarse samples
    u_vals   = sample_function_at_points(mesh, u_n, phys_u)[:,0]
    p_u_vals = sample_function_at_points(mesh, p_n, phys_u)
    v_vals   = sample_function_at_points(mesh, u_n, phys_v)[:,1]
    p_v_vals = sample_function_at_points(mesh, p_n, phys_v)

    # fine samples for smooth plotting
    u_plot   = sample_function_at_points(mesh, u_n, phys_u_f)[:,0]
    v_plot   = sample_function_at_points(mesh, u_n, phys_v_f)[:,1]

    return (
        u_vals, p_u_vals,    # coarse (y_ref/x_ref) samples
        v_vals, p_v_vals,
        y_plot, x_plot,      # the fine grids
        u_plot, v_plot,dxs       # smooth samples on the fine grids
    )


if __name__ == "__main__":
    comm = MPI.COMM_WORLD
    rank = comm.rank

    Nx_list = [20,40,50]
    u_plot_list, v_plot_list = [], []
    dx_list = []
    du_list, dv_list = [], []
    err_L2_u, err_inf_u = [], []
    err_L2_v, err_inf_v = [], []

    for Nx in Nx_list:
        print(f"\nRunning Nx={Nx}")
        (u_vals, pu_vals,
            v_vals, pv_vals,
            y_plot, x_plot,
            u_plot, v_plot,dxs) = solve_lid_cavity(Nx, 40.0, 1000.0)

        du = u_vals - u_ref
        dv = v_vals - v_ref

        u_plot_list.append(u_plot)
        v_plot_list.append(v_plot)
        du_list.append(du)
        dv_list.append(dv)
        dx_list.append(dxs)

        err_L2_u.append(np.sqrt(np.mean(du**2)))
        err_inf_u.append(np.max(np.abs(du)))
        err_L2_v.append(np.sqrt(np.mean(dv**2)))
        err_inf_v.append(np.max(np.abs(dv)))

    if rank == 0:
        fig, axs = plt.subplots(2, 2, figsize=(12, 8))

        # Top‐left: Table 9 reference and simulations
        axs[0, 0].scatter(y_ref, u_ref, marker='o', color='k', label='Reference')
        for Nx, u_plot, dxs in zip(Nx_list, u_plot_list, dx_list):
            axs[0, 0].plot(y_plot, u_plot, '-', label=rf'$\Delta x={dxs:.4f}$')
        axs[0, 0].set_xlabel('y')
        axs[0, 0].set_ylabel('u')
        axs[0, 0].set_title(r'$u_{x}$ along x=0.5')
        axs[0, 0].legend()
     

        # Top‐right: Table 10 reference and simulations
        axs[0, 1].scatter(x_ref, v_ref, marker='o', color='k', label='Reference')
        for Nx, v_plot, dxs in zip(Nx_list, v_plot_list, dx_list):
            axs[0, 1].plot(x_plot, v_plot, '-', label=rf'$\Delta x={dxs:.4f}$')
        axs[0, 1].set_xlabel('x')
        axs[0, 1].set_ylabel('v')
        axs[0, 1].set_title(r'$u_{y}$ along y=0.5')
        axs[0, 1].legend()


        # Bottom‐left: Error metrics for u (RMSE & max abs vs Nx)
        axs[1, 0].plot(dx_list, err_L2_u, 'o-', label='RMSE')
        axs[1, 0].plot(dx_list, err_inf_u, 's--', label='MAE')
        axs[1, 0].set_xlabel(r'$\Delta x$')
        axs[1, 0].set_ylabel('Error')
        # axs[1, 0].set_title('Error metrics: Horizontal velocity')
        axs[1,0].set_yscale("log")
        axs[1,0].set_xscale("log")
        axs[1, 0].legend()

        # Bottom‐right: Error metrics for v (RMSE & max abs vs Nx)
        axs[1, 1].plot(dx_list, err_L2_v, 'o-', label='RMSE')
        axs[1, 1].plot(dx_list, err_inf_v, 's--', label='MAE')
        axs[1, 1].set_xlabel(r'$\Delta x$')
        axs[1, 1].set_ylabel('Error')
        axs[1,1].set_yscale("log")
        axs[1,1].set_xscale("log")
        # axs[1, 1].set_title('Error metrics: Vertical velocity')
        axs[1, 1].legend()


        plt.tight_layout()
        plt.savefig("ipcs_a_liddriven.png",dpi=300)
        plt.show()
2

Is this working with the «first» navier stokes splitting scheme, ie Test problem 1: Channel flow (Poiseuille flow) — FEniCSx tutorial

I would recommend Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial for reasons explained in Splitting schemes — Oasisx

3

Hi @dokken Thanks for the follow up.

I have similarly tested your channel flow benchmark and it yielded similarly strange results:

channel_benchmark
channel_benchmark1500×1200 61.7 KB

I have adapted your code (similarly using CFL = 0.5 for the 2D CFL number to specify dt):

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np


from dolfinx.fem import Constant, Function, functionspace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
from dolfinx.io import VTXWriter
from dolfinx.mesh import create_unit_square
from dolfinx.plot import vtk_mesh
from basix.ufl import element
from ufl import (FacetNormal, Identity, TestFunction, TrialFunction,
                 div, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym)

import matplotlib.pyplot as plt

def evaluate_error(N):

    mesh = create_unit_square(MPI.COMM_WORLD, N, N)
    #Evaluate dx, dt base on CFL
    CFL     = 0.5
    T       = 5
    t=0
    coords = mesh.geometry.x
    xs = np.unique(coords[:, 0]); ys = np.unique(coords[:, 1])
    dxs = float(np.min(np.diff(np.sort(xs))))
    dy = float(np.min(np.diff(np.sort(ys))))

    dt = CFL * (dxs * dy) / (dxs + dy)
    num_steps = int(T/dt)+1


    v_cg2 = element("Lagrange", mesh.topology.cell_name(), 2, shape=(mesh.geometry.dim, ))
    s_cg1 = element("Lagrange", mesh.topology.cell_name(), 1)
    V = functionspace(mesh, v_cg2)
    Q = functionspace(mesh, s_cg1)

    u = TrialFunction(V)
    v = TestFunction(V)
    p = TrialFunction(Q)
    q = TestFunction(Q)

    def walls(x):
        return np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1))
    wall_dofs = locate_dofs_geometrical(V, walls)
    u_noslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
    bc_noslip = dirichletbc(u_noslip, wall_dofs, V)

    def inflow(x):
        return np.isclose(x[0], 0)
    inflow_dofs = locate_dofs_geometrical(Q, inflow)
    bc_inflow = dirichletbc(PETSc.ScalarType(8), inflow_dofs, Q)

    def outflow(x):
        return np.isclose(x[0], 1)
    outflow_dofs = locate_dofs_geometrical(Q, outflow)
    bc_outflow = dirichletbc(PETSc.ScalarType(0), outflow_dofs, Q)
    bcu = [bc_noslip]
    bcp = [bc_inflow, bc_outflow]

    u_n = Function(V)
    u_n.name = "u_n"
    U = 0.5 * (u_n + u)
    n = FacetNormal(mesh)
    f = Constant(mesh, PETSc.ScalarType((0, 0)))
    k = Constant(mesh, PETSc.ScalarType(dt))
    mu = Constant(mesh, PETSc.ScalarType(1))
    rho = Constant(mesh, PETSc.ScalarType(1))

    # Define strain-rate tensor
    def epsilon(u):
        return sym(nabla_grad(u))

    # Define stress tensor


    def sigma(u, p):
        return 2 * mu * epsilon(u) - p * Identity(len(u))


    # Define the variational problem for the first step
    p_n = Function(Q)
    p_n.name = "p_n"
    F1 = rho * dot((u - u_n) / k, v) * dx
    F1 += rho * dot(dot(u_n, nabla_grad(u_n)), v) * dx
    F1 += inner(sigma(U, p_n), epsilon(v)) * dx
    F1 += dot(p_n * n, v) * ds - dot(mu * nabla_grad(U) * n, v) * ds
    F1 -= dot(f, v) * dx
    a1 = form(lhs(F1))
    L1 = form(rhs(F1))

    A1 = assemble_matrix(a1, bcs=bcu)
    A1.assemble()
    b1 = create_vector(L1)

    # Define variational problem for step 2
    u_ = Function(V)
    a2 = form(dot(nabla_grad(p), nabla_grad(q)) * dx)
    L2 = form(dot(nabla_grad(p_n), nabla_grad(q)) * dx - (rho / k) * div(u_) * q * dx)
    A2 = assemble_matrix(a2, bcs=bcp)
    A2.assemble()
    b2 = create_vector(L2)

    # Define variational problem for step 3
    p_ = Function(Q)
    a3 = form(rho * dot(u, v) * dx)
    L3 = form(rho * dot(u_, v) * dx - k * dot(nabla_grad(p_ - p_n), 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.HYPRE)
    pc1.setHYPREType("boomeramg")

    # Solver for step 2
    solver2 = PETSc.KSP().create(mesh.comm)
    solver2.setOperators(A2)
    solver2.setType(PETSc.KSP.Type.BCGS)
    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)

    def u_exact(x):
        values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
        values[0] = 4 * x[1] * (1.0 - x[1])
        return values


    u_ex = Function(V)
    u_ex.interpolate(u_exact)

    L2_error = form(dot(u_ - u_ex, u_ - u_ex) * dx)

    for i in range(num_steps):
        # Update current time step
        t += dt

        # Step 1: Tentative veolcity step
        with b1.localForm() as loc_1:
            loc_1.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_.x.petsc_vec)
        u_.x.scatter_forward()

        # Step 2: Pressure corrrection step
        with b2.localForm() as loc_2:
            loc_2.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, p_.x.petsc_vec)
        p_.x.scatter_forward()

        # Step 3: Velocity correction step
        with b3.localForm() as loc_3:
            loc_3.set(0)
        assemble_vector(b3, L3)
        b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
        solver3.solve(b3, u_.x.petsc_vec)
        u_.x.scatter_forward()
        # Update variable with solution form this time step
        u_n.x.array[:] = u_.x.array[:]
        p_n.x.array[:] = p_.x.array[:]


        # Compute error at current time-step
        error_L2 = np.sqrt(mesh.comm.allreduce(assemble_scalar(L2_error), op=MPI.SUM))
        error_max = mesh.comm.allreduce(np.max(u_.x.petsc_vec.array - u_ex.x.petsc_vec.array), op=MPI.MAX)
        # Print error only every 20th step and at the last step
        if (i % 20 == 0):
            print(f"Time {t:.2f}, L2-error {error_L2:.2e}, Max error {error_max:.2e}")

    b1.destroy()
    b2.destroy()
    b3.destroy()
    solver1.destroy()
    solver2.destroy()
    solver3.destroy()

    return error_L2, dxs

"""
Loop
"""

Ns = [10,20,25,40,50]

results = []
dxs_ = []

for N in Ns:
    err, dxs = evaluate_error(N)
    results.append(err)
    dxs_.append(dxs)

plt.figure(figsize=(5,4))
plt.loglog(dxs_,results, marker="o")
plt.xlabel("Mesh spacing h")
plt.ylabel(r"$\|u - u_{\mathrm{exact}}\|_{L^2(\Omega)}$")
plt.title("L2 error vs mesh spacing")
plt.grid(False)
plt.tight_layout()
plt.savefig("channel_benchmark.png",dpi=300)
4

I have also tested the scheme developed in your second test problem and it performs similarly strangely.

ipcs_b_liddriven
ipcs_b_liddriven3600×2400 351 KB

"""
Lid‐driven cavity IPCS in DOLFINx, with center‐line sampling of u and p
moved inside the solver function.  Uses Marchi et al, 2021 Tables 20 & 21
for reference.
"""

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path

import dolfinx
from dolfinx.fem import (Constant, Function, functionspace,
                         dirichletbc, locate_dofs_geometrical,
                         form, Expression)
from dolfinx.fem.petsc import (assemble_matrix, assemble_vector,
                               apply_lifting, create_vector, set_bc)
from dolfinx.mesh import create_unit_square
from basix.ufl import element
from ufl import (FacetNormal, Identity, TestFunction, TrialFunction,
                 div, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym, grad)

from dolfinx import geometry

#formatting
plt.rcParams["text.usetex"] = True
plt.rc('font', family='serif')

# ————— Reference data from Tables 9 & 10 —————
y_ref   = np.array([0.0546875, 0.0625, 0.0703125, 0.1015625,0.125,0.171875,0.1875,
                    0.25,0.28125,0.3125,0.375,0.4375,0.453125,0.5,0.5625,0.6171875,
                    0.625,0.6875,0.734375,0.75,0.8125,0.8515625,0.875,
                    0.9375,0.953125,0.9609375,0.96875,0.9765625])

u_ref = np.array([-1.8125321e-1,-2.023291e-1,-2.229271e-1,-3.003688e-1,-3.478431e-1,-3.885676e-1,-3.844078e-1,
                  -3.1894525e-1,-2.8042785e-1,-2.4569306e-1,-1.8373161e-1,-1.2341016e-1, -1.0817529e-1,-6.205605e-2,5.6167e-4,
                  5.700436e-2,6.524839e-2,1.3357199e-1,1.8864358e-1,2.0791378e-1,2.884412e-1,3.371763e-1,3.625438e-1,
                  4.229307e-1,4.724491e-1,5.171830e-1,5.803655e-1,6.6397203e-1])



x_ref = np.array([
    0.0625,   0.0703125, 0.078125,   0.09375,
    0.125,    0.15625,   0.1875,    0.2265625,
    0.234375, 0.25,      0.3125,    0.375,
    0.4375,   0.5,       0.5625,    0.625,
    0.6875,   0.75,      0.8046875, 0.8125,
    0.859375, 0.875,     0.90625,   0.9375,
    0.9453125,0.953125,  0.9609375, 0.96875
])

v_ref = np.array([
    2.807042e-1,   2.962922e-1,   3.099493e-1,   3.329768e-1,
    3.650400e-1,   3.769157e-1,   3.678512e-1,   3.340319e-1,
    3.253866e-1,   3.0710340e-1,  2.3126780e-1,  1.6056381e-1,
    0.9296911e-1,  2.5799473e-2,  -4.184046e-2,   -1.1079783e-1,
    -1.8167907e-1,  -2.533806e-1,   -3.201954e-1,   -3.31565806e-1,
    -4.263892e-1,   -4.677741e-1,   -5.264155e-1,   -4.561505e-1,
    -4.102923e-1,   -3.551312e-1,   -2.933786e-1,   -2.283407e-1
])




def sample_function_at_points(mesh, func, phys_points):
    # build bounding‐box tree
    bb = geometry.bb_tree(mesh, mesh.topology.dim)
    # find candidate cells for each point (points as 3×N in tutorial, so transpose)
    candidates = geometry.compute_collisions_points(bb, phys_points)
    # filter to actual containing cells
    colliding  = geometry.compute_colliding_cells(mesh, candidates, phys_points)

    pts_on_proc, cells, indices = [], [], []
    for i, p in enumerate(phys_points):
        links = colliding.links(i)
        if len(links) > 0:
            pts_on_proc.append(p)
            cells.append(links[0])
            indices.append(i)

    if not pts_on_proc:
        # no sample on this rank
        return np.full(len(phys_points), np.nan)

    pts_on_proc = np.array(pts_on_proc, dtype=np.float64)
    cells       = np.array(cells, dtype=np.int32)
    # evaluate your FE function
    vals = func.eval(pts_on_proc, cells)

    # flatten scalar result
    if vals.ndim == 2 and vals.shape[1] == 1:
        vals = vals[:, 0]

    # reassemble full‐length array
    if vals.ndim == 1:
        full = np.full(len(phys_points), np.nan)
        for idx, v in zip(indices, vals):
            full[idx] = v
    else:
        # vector case
        dim = vals.shape[1]
        full = np.full((len(phys_points), dim), np.nan)
        for idx, v in zip(indices, vals):
            full[idx, :] = v

    return full

def solve_lid_cavity(Nx, T, Re):
    #Set up time stepping
    Re      = Re
    CFL     = 0.5
    T       = T
    stride  = 50
    nu_val  = 1.0 / Re
    Nx=Ny = Nx

    #Set up mesh
    domain = create_unit_square(MPI.COMM_WORLD, Nx, Ny)
    v_cg2 = element("Lagrange", domain.topology.cell_name(), 2, shape=(domain.geometry.dim, ))
    s_cg1 = element("Lagrange", domain.topology.cell_name(), 1)
    V = functionspace(domain, v_cg2)
    Q = functionspace(domain, s_cg1)

    coords = domain.geometry.x
    xs = np.unique(coords[:, 0]); ys = np.unique(coords[:, 1])
    dxs = float(np.min(np.diff(np.sort(xs))))
    dy = float(np.min(np.diff(np.sort(ys))))

    dt_value = CFL * (dxs * dy) / (dxs + dy)
    dt = PETSc.ScalarType(dt_value)

    dt_const = Constant(domain, dt)
    nu_const = Constant(domain, PETSc.ScalarType(nu_val))


    def top(x): return np.isclose(x[1], 1.0)

    def walls(x):
        return np.logical_or(
            np.isclose(x[1], 0.0),
            np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0))
        )

    wall_dofs = locate_dofs_geometrical(V, walls)
    u_noslip = np.array((0.0,)*domain.geometry.dim, dtype=PETSc.ScalarType)
    bc_walls = dirichletbc(u_noslip, wall_dofs, V)

    lid_dofs = locate_dofs_geometrical(V, top)
    u_lid   = np.array((1.0, 0.0), dtype=PETSc.ScalarType)
    bc_lid   = dirichletbc(u_lid, lid_dofs, V)

    bcs_u = [bc_lid, bc_walls]

    p_dofs   = locate_dofs_geometrical(
        Q, lambda x: np.logical_and(np.isclose(x[0], 0.0),
                                     np.isclose(x[1], 0.0))
    )
    bc_p     = [dirichletbc(Constant(domain, 0.0),p_dofs, Q)]
    

    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(domain, PETSc.ScalarType((0, 0)))

    F1 = dot((u - u_n) / dt_const, 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 * nu_const * 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 = assemble_matrix(a1, bcs=bcs_u)
    A1.assemble()
    b1 = create_vector(L1)

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

    F3 = dot(u, v)*dx + dt_const*dot(nabla_grad(phi), v)*dx
    F3 -= dot(u_s, v)*dx
    a3, L3 = form(lhs(F3)), form(rhs(F3))
    A3 = assemble_matrix(a3)
    A3.assemble()
    b3 = create_vector(L3)

    solver1 = PETSc.KSP().create(domain.comm)
    solver1.setOperators(A1)
    solver1.setType(PETSc.KSP.Type.BCGS)
    pc1 = solver1.getPC()
    pc1.setType(PETSc.PC.Type.JACOBI)

    solver2 = PETSc.KSP().create(domain.comm)
    solver2.setOperators(A2)
    solver2.setType(PETSc.KSP.Type.BCGS)
    pc2 = solver2.getPC()
    pc2.setType(PETSc.PC.Type.HYPRE)
    pc2.setHYPREType("boomeramg")

    solver3 = PETSc.KSP().create(domain.comm)
    solver3.setOperators(A3)
    solver3.setType(PETSc.KSP.Type.CG)
    pc3 = solver3.getPC()
    pc3.setType(PETSc.PC.Type.SOR)


    step = 0
    t = 0.0

    while t < T - 1e-6:
        t += float(dt)
        print(t)

        with b1.localForm() as loc_1:
            loc_1.set(0)
        assemble_vector(b1, L1)
        apply_lifting(b1, [a1], [bcs_u])
        b1.ghostUpdate(addv=PETSc.InsertMode.ADD,
                       mode=PETSc.ScatterMode.REVERSE)
        set_bc(b1, bcs_u)
        solver1.solve(b1, u_s.x.petsc_vec)
        u_s.x.scatter_forward()

        with b2.localForm() as loc_2:
            loc_2.set(0)
        assemble_vector(b2, L2)
        apply_lifting(b2, [a2], [bc_p])
        b2.ghostUpdate(addv=PETSc.InsertMode.ADD,
                       mode=PETSc.ScatterMode.REVERSE)
        set_bc(b2, bc_p)
        solver2.solve(b2, phi.x.petsc_vec)
        phi.x.scatter_forward()
        p_.x.petsc_vec.axpy(1, phi.x.petsc_vec)
        p_.x.scatter_forward()

        with b3.localForm() as loc_3:
            loc_3.set(0)
        assemble_vector(b3, L3)
        b3.ghostUpdate(addv=PETSc.InsertMode.ADD,
                       mode=PETSc.ScatterMode.REVERSE)
        solver3.solve(b3, u_.x.petsc_vec)
        u_.x.scatter_forward()

        u_n.x.array[:] = u_.x.array[:]
        # p_n.x.array[:] = p_.x.array[:]

        step += 1
        
    b2.destroy()
    b3.destroy()
    solver1.destroy()
    solver2.destroy()
    solver3.destroy()

    # ————— sample at final time —————
    mesh = u_n.function_space.mesh
    # coarse points for error
    if MPI.COMM_WORLD.rank == 0:
        phys_u     = np.column_stack([0.5*np.ones_like(y_ref), y_ref,     np.zeros_like(y_ref)])
        phys_v     = np.column_stack([x_ref,               0.5*np.ones_like(x_ref), np.zeros_like(x_ref)])
        # fine points for plotting
        y_plot     = np.linspace(0.0, 1.0, 100)
        x_plot     = np.linspace(0.0, 1.0, 100)
        phys_u_f   = np.column_stack([0.5*np.ones_like(y_plot), y_plot, np.zeros_like(y_plot)])
        phys_v_f   = np.column_stack([x_plot, 0.5*np.ones_like(x_plot), np.zeros_like(x_plot)])
    else:
        phys_u = phys_v = phys_u_f = phys_v_f = np.empty((0,3), dtype=np.float64)

    # coarse samples
    u_vals   = sample_function_at_points(mesh, u_n, phys_u)[:,0]
    # p_u_vals = sample_function_at_points(mesh, p_n, phys_u)
    v_vals   = sample_function_at_points(mesh, u_n, phys_v)[:,1]
    # p_v_vals = sample_function_at_points(mesh, p_n, phys_v)

    # fine samples for smooth plotting
    u_plot   = sample_function_at_points(mesh, u_n, phys_u_f)[:,0]
    v_plot   = sample_function_at_points(mesh, u_n, phys_v_f)[:,1]

    return (
        u_vals,    # coarse (y_ref/x_ref) samples
        v_vals,
        y_plot, x_plot,      # the fine grids
        u_plot, v_plot,dxs       # smooth samples on the fine grids
    )


if __name__ == "__main__":
    comm = MPI.COMM_WORLD
    rank = comm.rank

    Nx_list = [20,40,50,80]
    u_plot_list, v_plot_list = [], []
    dx_list = []
    du_list, dv_list = [], []
    err_L2_u, err_inf_u = [], []
    err_L2_v, err_inf_v = [], []

    for Nx in Nx_list:
        print(f"\nRunning Nx={Nx}")
        (u_vals,
            v_vals,
            y_plot, x_plot,
            u_plot, v_plot,dxs) = solve_lid_cavity(Nx, 30.0, 1000.0)

        du = u_vals - u_ref
        dv = v_vals - v_ref

        u_plot_list.append(u_plot)
        v_plot_list.append(v_plot)
        du_list.append(du)
        dv_list.append(dv)
        dx_list.append(dxs)

        err_L2_u.append(np.sqrt(np.mean(du**2)))
        err_inf_u.append(np.max(np.abs(du)))
        err_L2_v.append(np.sqrt(np.mean(dv**2)))
        err_inf_v.append(np.max(np.abs(dv)))

    if rank == 0:
        fig, axs = plt.subplots(2, 2, figsize=(12, 8))

        # Top‐left: Table 9 reference and simulations
        axs[0, 0].scatter(y_ref, u_ref, marker='o', color='k', label='Reference')
        for Nx, u_plot, dxs in zip(Nx_list, u_plot_list, dx_list):
            axs[0, 0].plot(y_plot, u_plot, '-', label=rf'$\Delta x={dxs:.4f}$')
        axs[0, 0].set_xlabel('y')
        axs[0, 0].set_ylabel('u')
        axs[0, 0].set_title(r'$u_{x}$ along x=0.5')
        axs[0, 0].legend()
     

        # Top‐right: Table 10 reference and simulations
        axs[0, 1].scatter(x_ref, v_ref, marker='o', color='k', label='Reference')
        for Nx, v_plot, dxs in zip(Nx_list, v_plot_list, dx_list):
            axs[0, 1].plot(x_plot, v_plot, '-', label=rf'$\Delta x={dxs:.4f}$')
        axs[0, 1].set_xlabel('x')
        axs[0, 1].set_ylabel('v')
        axs[0, 1].set_title(r'$u_{y}$ along y=0.5')
        axs[0, 1].legend()


        # Bottom‐left: Error metrics for u (RMSE & max abs vs Nx)
        axs[1, 0].plot(dx_list, err_L2_u, 'o-', label='RMSE')
        axs[1, 0].plot(dx_list, err_inf_u, 's--', label='MAE')
        axs[1, 0].set_xlabel(r'$\Delta x$')
        axs[1, 0].set_ylabel('Error')
        # axs[1, 0].set_title('Error metrics: Horizontal velocity')
        axs[1,0].set_yscale("log")
        axs[1,0].set_xscale("log")
        axs[1, 0].legend()

        # Bottom‐right: Error metrics for v (RMSE & max abs vs Nx)
        axs[1, 1].plot(dx_list, err_L2_v, 'o-', label='RMSE')
        axs[1, 1].plot(dx_list, err_inf_v, 's--', label='MAE')
        axs[1, 1].set_xlabel(r'$\Delta x$')
        axs[1, 1].set_ylabel('Error')
        axs[1,1].set_yscale("log")
        axs[1,1].set_xscale("log")
        # axs[1, 1].set_title('Error metrics: Vertical velocity')
        axs[1, 1].legend()


        plt.tight_layout()
        plt.savefig("ipcs_b_liddriven.png",dpi=300)
        plt.show()
5

Hi Pavan,

Could you try running direct solvers (ie LU factorization for all steps)?

As iterative solvers has a stopping tolerance (I think around 1e-6 by default that will definitely disrupt converge rates at the error size you are seeing. Alternatively, try setting atol and rtol small (like 1e-14).

6

Hi @dokken

That worked perfectly! Setting solver.atol / solver.rtol to 1e-14 gives sensible results!

channel_benchmark
channel_benchmark1500×1200 52.1 KB

7

For posterity’s sake:

Turns out playing with the linear solver tolerances is just a partial solution to the original lid driven cavity problem. It turns out it takes longer to reach steady state. In the original simulation, i ran it to t=30. I have run it to t=100 and we get sensible results.

ipcs_a2_liddriven
ipcs_a2_liddriven3600×2400 353 KB