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/stokes_with_pressure_bc.py", line 127, in <module>
norm_u_1, norm_p_1 = block_iterative_solver()
File "/home/user/my_env/stokes_with_pressure_bc.py", line 82, in block_iterative_solver
A, P, b = block_operators()
File "/home/user/my_env/stokes_with_pressure_bc.py", line 65, in block_operators
A = fem.petsc.assemble_matrix_block(a, bcs=bcs)
File "/home/user/anaconda3/envs/my_env/lib/python3.10/functools.py", 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/petsc.py", 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/petsc.py", 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)
in_pressure.interpolate(inflow_pressure_expression)
facets = locate_entities_boundary(mesh, facetdim, inflow)
bc_in = dirichletbc(in_pressure, locate_dofs_topological(Q, facetdim, facets))
out_pressure = Function(Q)
out_pressure.interpolate(outflow_pressure_expression)
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)
A.assemble()
P = fem.petsc.assemble_matrix_block(a_p, bcs=bcs)
P.assemble()
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
null_vec.normalize()
nsp = PETSc.NullSpace().create(vectors=[null_vec])
# assert nsp.test(A)
A.setNullSpace(nsp)
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.setTolerances(rtol=1e-9)
ksp.setType("minres")
ksp.getPC().setType("fieldsplit")
ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)
ksp.getPC().setFieldSplitIS(("u", is_u), ("p", is_p))
# Configure velocity and pressure sub-solvers
ksp_u, ksp_p = ksp.getPC().getFieldSplitSubKSP()
ksp_u.setType("preonly")
ksp_u.getPC().setType("gamg")
ksp_p.setType("preonly")
ksp_p.getPC().setType("jacobi")
ksp.getPC().setUp()
Pu, _ = ksp_u.getPC().getOperators()
Pu.setBlockSize(mesh.topology.dim)
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,
Georg