Navier stokes - block iterative solver with pressure BC

Dear FenicsX-community,

I am trying to adapt this lid-driven-cavity-example to have an inflow and an outflow with pressure boundary conditions and a no-slip condition on the remaining walls. I can’t use the poiseuille-flow-demo because this is already a minimal example for a more complex mesh.

In the following code I mainly adapted those boundary conditions and kept the solver as it is.
I am getting the following error which is probably related to the block-diagonal preconditioner which I don’t understand fully:

Traceback (most recent call last):
  File "/home/user/my_env/", line 127, in <module>
    norm_u_1, norm_p_1 = block_iterative_solver()
  File "/home/user/my_env/", line 82, in block_iterative_solver
    A, P, b = block_operators()
  File "/home/user/my_env/", line 65, in block_operators
    A = fem.petsc.assemble_matrix_block(a, bcs=bcs)
  File "/home/user/anaconda3/envs/my_env/lib/python3.10/", line 889, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
  File "/home/user/anaconda3/envs/my_env/lib/python3.10/site-packages/dolfinx/fem/", line 431, in _assemble_matrix_block_form
    return _assemble_matrix_block_mat(A, a, bcs, diagonal, constants, coeffs)
  File "/home/user/anaconda3/envs/my_env/lib/python3.10/site-packages/dolfinx/fem/", line 461, in _assemble_matrix_block_mat
    raise RuntimeError(
RuntimeError: Diagonal sub-block (1, 1) cannot be 'None' and have DirichletBC applied. Consider assembling a zero block.

Minimal code:

from dolfinx import *
from dolfinx.mesh import locate_entities_boundary, create_rectangle, CellType
from dolfinx.fem import form, Constant, dirichletbc, FunctionSpace, locate_dofs_topological, Function
from mpi4py import MPI
import numpy as np
import ufl
from petsc4py import PETSc
from ufl import inner, grad, dx, div

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

def inflow_pressure_expression(x):
    return np.ones(x.shape[0]) * 5

def outflow_pressure_expression(x):
    return np.ones(x.shape[0]) * -5

t = 0
T = 0.2
num_steps = 20
dt = T / num_steps

mesh = create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])],
                        [32, 32], CellType.triangle)

P2 = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V, Q = FunctionSpace(mesh, P2), FunctionSpace(mesh, P1)

noslip = np.zeros(mesh.geometry.dim, dtype=PETSc.ScalarType)
facetdim = mesh.topology.dim - 1
facets = locate_entities_boundary(mesh, facetdim, walls)
bc0 = dirichletbc(noslip, locate_dofs_topological(V, facetdim, facets), V)
in_pressure = Function(Q)
facets = locate_entities_boundary(mesh, facetdim, inflow)
bc_in = dirichletbc(in_pressure, locate_dofs_topological(Q, facetdim, facets))
out_pressure = Function(Q)
facets = locate_entities_boundary(mesh, facetdim, outflow)
bc_out = dirichletbc(out_pressure, locate_dofs_topological(Q, facetdim, facets))
bcs = [bc0, bc_in, bc_out]

(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
f = Constant(mesh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))

a = form([[inner(grad(u), grad(v)) * dx, inner(p, div(v)) * dx],
          [inner(div(u), q) * dx, None]])
L = form([inner(f, v) * dx, inner(Constant(mesh, PETSc.ScalarType(0)), q) * dx])

a_p11 = form(inner(p, q) * dx)
a_p = [[a[0][0], None],
       [None, a_p11]]

def block_operators():
    A = fem.petsc.assemble_matrix_block(a, bcs=bcs)
    P = fem.petsc.assemble_matrix_block(a_p, bcs=bcs)
    b = fem.petsc.assemble_vector_block(L, a, bcs=bcs)

    null_vec = A.createVecLeft()
    offset = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
    null_vec.array[offset:] = 1.0
    nsp = PETSc.NullSpace().create(vectors=[null_vec])
    # assert nsp.test(A)

    return A, P, b

def block_iterative_solver():
    A, P, b = block_operators()

    V_map = V.dofmap.index_map
    Q_map = Q.dofmap.index_map
    offset_u = V_map.local_range[0] * V.dofmap.index_map_bs + Q_map.local_range[0]
    offset_p = offset_u + V_map.size_local * V.dofmap.index_map_bs
    is_u = PETSc.IS().createStride(V_map.size_local * V.dofmap.index_map_bs, offset_u, 1, comm=PETSc.COMM_SELF)
    is_p = PETSc.IS().createStride(Q_map.size_local, offset_p, 1, comm=PETSc.COMM_SELF)

    ksp = PETSc.KSP().create(mesh.comm)
    ksp.setOperators(A, P)
    ksp.getPC().setFieldSplitIS(("u", is_u), ("p", is_p))

    # Configure velocity and pressure sub-solvers
    ksp_u, ksp_p = ksp.getPC().getFieldSplitSubKSP()

    Pu, _ = ksp_u.getPC().getOperators()

    x = A.createVecRight()
    ksp.solve(b, x)

    u, p = Function(V), Function(Q)
    offset = V_map.size_local * V.dofmap.index_map_bs
    u.x.array[:offset] = x.array_r[:offset]
    p.x.array[:(len(x.array_r) - offset)] = x.array_r[offset:]

    norm_u, norm_p = u.x.norm(), p.x.norm()
    if MPI.COMM_WORLD.rank == 0:
        print(f"(B) Norm of velocity coefficient vector (blocked, iterative): {norm_u}")
        print(f"(B) Norm of pressure coefficient vector (blocked, iterative): {norm_p}")

    return norm_u, norm_p

for i in range(num_steps):
    t += dt
    norm_u_1, norm_p_1 = block_iterative_solver()

Many thanks in advance for your help.
Best regards,

The error is quite clear, the block (1,1) of a

which is None, cannot have a DirichletBC applied (which is does through bcs).

I.e. using something like:

a = form([[inner(grad(u), grad(v)) * dx, inner(p, div(v)) * dx],
          [inner(div(u), q) * dx, Constant(mesh, PETSc.ScalarType(0))*p*q*dx(99999)]])
L = form([inner(f, v) * dx, inner(Constant(mesh, PETSc.ScalarType(0)), q) * dx(99999)])

which assembles a zero blocks, resolves the issue.

Here I use (99999) inside dx as no cell is marked with this, improving the efficiency of the assembly code.

1 Like