Hello,
I am also trying to solve the Stokes equation with pure Neumann boundary condition, but I don’t manage to define the correct nullspace.
The problem I am trying to solve is
\begin{align}
\Delta {\bf u} - \nabla p = 0 \text{ in } \Omega \\
\text{div} \, {\bf u} = 0 \text{ in } \Omega \\
(\nabla {\bf u}) \cdot {\bf n} - p \, {\bf n} = {\bf g} \text{ on } \partial \Omega
\end{align}
and g satisfies the compatibility condition \int_{\partial \Omega} {\bf g} = 0.
Here is the minimal working example with a tentative nullspace:
from mpi4py import MPI
from petsc4py import PETSc
dtype = PETSc.ScalarType
import ufl
from basix.ufl import element, mixed_element
import dolfinx
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 40, 40)
V_el = element("Lagrange", domain.basix_cell(), 2, shape=(2,))
Q_el = element("Lagrange", domain.basix_cell(), 1)
W_el = mixed_element([V_el,Q_el])
W = dolfinx.fem.functionspace(domain, W_el)
V,V_to_W = W.sub(0).collapse()
Q,Q_to_W = W.sub(1).collapse()
## Definition of the linear forms
(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
a = dolfinx.fem.form([[ ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx , ufl.inner(p, ufl.div(v)) * ufl.dx ], [ufl.inner(ufl.div(u), q) * ufl.dx, None]])
## Neumann BC
x = ufl.SpatialCoordinate(domain)
g = ufl.as_vector([2*x[0]-1, 0.2*x[1]-0.1])
L = dolfinx.fem.form([ ufl.inner(g, v) * ufl.ds, ufl.inner(dolfinx.fem.Constant(domain, dtype(0)), q) * ufl.dx ])
# Assemble the block operator and RHS vector
A = assemble_matrix_block(a)
A.assemble()
b = assemble_vector_block(L, a)
## Set the nullspace
bs = W.dofmap.index_map_bs
length0 = W.dofmap.index_map.size_local
basis = [dolfinx.la.vector(W.dofmap.index_map, bs=bs, dtype=dtype) for i in range(3)]
basis_arrays = [_b.array for _b in basis]
dofs_v = [V.sub(i).dofmap.list.flatten() for i in range(2)]
dofs_p = [Q.dofmap.list.flatten()]
basis_arrays[0][dofs_v[0]] = 1.0
basis_arrays[1][dofs_v[1]] = 1.0
basis_arrays[2][dofs_p] = 1.0
dolfinx.la.orthonormalize(basis)
basis_petsc = [
PETSc.Vec().createWithArray(x[: bs * length0], bsize=bs, comm=W.mesh.comm)
for x in basis_arrays
]
nullspace = PETSc.NullSpace().create(vectors=basis_petsc)
assert nullspace.test(A)
A.setNearNullSpace(nullspace)
# Create a block vector (x) to store the full solution, and solve
# Create a solver
ksp = PETSc.KSP().create(domain.comm)
ksp.setOperators(A)
ksp.setType("preonly")
# Set the solver type to MUMPS (LU solver)
pc = ksp.getPC()
pc.setType("lu")
pc.setFactorSolverType("superlu_dist")
pc.setFactorSetUpSolverType()
# Create a block vector (x) to store the full solution, and solve
x = A.createVecLeft()
ksp.solve(b, x)
# Create Functions and scatter x solution
u, p = dolfinx.fem.Function(V), dolfinx.fem.Function(Q)
offset = V.dofmap.index_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:]
The assert nullspace.test(A) raises an error. Is there something wrong with the way I am creating the nullspace? I am quite new to dolfinx.
Any help is appreciated!
