Stokes problem in 3D -DOLFINx

General
1

Good day to the community, I have a problem reproducing the code of [1]. I have already updated the code to DOLFINX version 0.8.0 but I get the following error:

Traceback (most recent call last):
  File "/home/oscar/viskex/nuevo1.py", line 79, in <module>
    A = fem.petsc.assemble_matrix(a, bcs=bcs)
  File "/usr/lib/python3.10/functools.py", line 889, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/petsc.py", line 424, in assemble_matrix
    A = _cpp.fem.petsc.create_matrix(a._cpp_object)
AttributeError: 'list' object has no attribute '_cpp_object'

I attach the code that I have updated:

import numpy as np
import dolfinx 

import ufl
from basix.ufl import element, mixed_element
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem
from dolfinx.fem import form
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
from ufl import div, grad, inner


N = 16
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, N, N, N,
                                     cell_type=dolfinx.mesh.CellType.hexahedron)
P2 = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))
P1 = element("Lagrange", mesh.basix_cell(), 1)
TH = mixed_element([P2, P1])
W = dolfinx.fem.functionspace(mesh, TH)


def right(x, tol=1e-14): return x[0] > 1 - tol


def top_bottom(x, tol=1e-14):
    return np.logical_or(
        np.isclose(x[1], 0, atol=tol), np.isclose(x[1], 1, atol=tol))


V, V_to_W = W.sub(0).collapse()

# No-slip boundary condition for velocity
nonslip = dolfinx.fem.Function(V)
nonslip.x.array[:]=1
nonslip_facets = dolfinx.mesh.locate_entities_boundary(
    mesh, mesh.topology.dim-1, top_bottom)
nonslip_dofs = dolfinx.fem.locate_dofs_topological(
    (W.sub(0), V), mesh.topology.dim-1, nonslip_facets)
bc0 = dolfinx.fem.dirichletbc(nonslip, nonslip_dofs, W.sub(0))

# Inflow boundary condition for velocity


def inflow(x):
    zero_comp = np.zeros(x.shape[1], dtype=np.float64)
    return (-np.sin(np.pi*x[1]), zero_comp, zero_comp)


inflow_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim-1,
                                                      right)
inflow_dofs = dolfinx.fem.locate_dofs_topological((W.sub(0), V), mesh.topology.dim-1,
                                                  inflow_facets)
u_bc = dolfinx.fem.Function(V)
u_bc.interpolate(inflow)
bc1 = dolfinx.fem.dirichletbc(u_bc, inflow_dofs, W.sub(0))

bcs = [bc0, bc1]

(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
f = dolfinx.fem.Constant(mesh, (0., 0., 0.))

dx = ufl.Measure("dx", domain=mesh)
# a = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx + \
    # ufl.div(v)*p*dx + q*ufl.div(u)*dx
# L = ufl.inner(f, v) * dx


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(dolfinx.fem.Constant(mesh, PETSc.ScalarType(0)), q) * dx])  # type: ignore


A = fem.petsc.assemble_matrix(a, bcs=bcs)
A.assemble()
b = dolfinx.fem.petsc.assemble_vector(L)
dolfinx.fem.petsc.apply_lifting(b, [a], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(b, bcs)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

# Preconditioner matrix form
Pf = dolfinx.fem.form(ufl.inner(ufl.grad(u), ufl.grad(v)) * dx + p * q * dx)
P = dolfinx.fem.petsc.assemble_matrix(Pf, bcs=bcs)
P.assemble()

ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A, P)
ksp.setType(PETSc.KSP.Type.MINRES)
pc = ksp.getPC()
pc.setType("hypre")
pc.setHYPREType("boomeramg")
pc.setUp()

# To view Preconditioner info
# pc.view()


wh = dolfinx.fem.Function(W)
ksp.solve(b, wh.vector)

# To view solver info
# ksp.view()

uh = wh.sub(0).collapse()
p = wh.sub(1).collapse()

print(f"Converged Reason {ksp.getConvergedReason()}" +
      f"\nNumber of iterations {ksp.getIterationNumber()}")

with dolfinx.io.VTXWriter(mesh.comm, "u.bp", [uh]) as vtx:
    vtx.write(0.)

with dolfinx.io.VTXWriter(mesh.comm, "p.bp", [uh]) as vtx:
    vtx.write(0.)

Thank you in advance for your help.

2

You seem to be mixing up two different codes.
You are using a mixed_element, hence your forms shouldn’t be block forms.
If you are willing to use block forms, then define the two spaces for velocity and pressure separately, and use block assembly from dolfinx.fem.petsc (e.g., assemble_matrix_block.