# Poisson problem with Neumann boundary condition

Good evening,

I am solving the singular Poisson problem, I saw that the best way to remove the singularity is to add a null space to the matrix, the theory I found behind the method for an iterative solver can be founded in the article “On the singular Neumann problem in linear elasticity”. Regarding the implementation the result I got are without any sense. I took inspiration from the post Issues Solving Pure Neumann Problems in dolfinx - #2 by dokken
Can someone tell me what I am missing?

The code is the following:

``````from mpi4py import MPI
from dolfinx import mesh, fem

import numpy as np
import pyvista

from dolfinx.fem import (Constant,  Function, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction,
div, dot, dx, grad, inner, lhs, rhs)
from dolfinx.io import XDMFFile
from dolfinx.plot import create_vtk_mesh
N = 160
domain = mesh.create_unit_square(MPI.COMM_WORLD, N, N, mesh.CellType.quadrilateral)
#domain = mesh.create_box(comm=MPI.COMM_WORLD,
#                            points=((0.0, 0.0, 0.0), (1.0, 1.0, 1.0)), n=(N, N, N))

from dolfinx.fem import (Constant, Function, FunctionSpace, dirichletbc, Expression,
extract_function_spaces, form, assemble_scalar,
locate_dofs_geometrical, locate_dofs_topological)
import numpy as np
from petsc4py import PETSc
import basix, basix.ufl_wrapper
V = FunctionSpace(domain, ("CG", 1))

def u_exact(x):
return np.cos(2*np.pi*x[0])

import ufl
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

from petsc4py.PETSc import ScalarType
x = ufl.SpatialCoordinate(domain)
f = 4*ufl.pi*ufl.pi * ufl.cos(2*ufl.pi*x[0])
g = -2*ufl.pi* ufl.sin(2*ufl.pi*x[0])

L = ufl.inner(f,v)* ufl.dx

boundaries = [(1, lambda x: np.isclose(x[0], 0)),
(2, lambda x: np.isclose(x[0], 1)),
(3, lambda x: np.isclose(x[1], 0)),
(4, lambda x: np.isclose(x[1], 1))]

facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(domain, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
with XDMFFile(domain.comm, "facet_tags.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_meshtags(facet_tag)

ds = Measure("ds", domain=domain, subdomain_data=facet_tag)

class BoundaryCondition():
def __init__(self, type, marker, values):
self._type = type
if type == "Dirichlet":
u_D = Function(V)
u_D.interpolate(values)
facets = facet_tag.find(marker)
dofs = locate_dofs_topological(V, fdim, facets)
self._bc = dirichletbc(u_D, dofs)
elif type == "Neumann":
self._bc = values*v* ds(marker)
else:
raise TypeError("Unknown boundary condition: {0:s}".format(type))
@property
def bc(self):
return self._bc

@property
def type(self):
return self._type

boundary_conditions = [BoundaryCondition("Neumann", 1, -g),
BoundaryCondition("Neumann", 2, g),
BoundaryCondition("Neumann", 3, 0),
BoundaryCondition("Neumann", 4, 0)]

bcs = []
for condition in boundary_conditions:
if condition.type == "Dirichlet":
bcs.append(condition.bc)
else:
L += condition.bc

'''def create_nullspace(rhs_form):
null_vec = fem.petsc.create_vector(rhs_form)
null_vec.set(1.0)
null_vec.normalize()
nsp = PETSc.NullSpace().create(null_vec)
return nsp

vec = fem.form(L)
nsp = create_nullspace(vec)
'''
######### ASSEMBLE THE NESTED SYSTEM
def assemble_system(lhs_form, rhs_form, bcs):
A = fem.petsc.assemble_matrix(lhs_form, bcs=bcs)
A.assemble()

b = fem.petsc.create_vector(rhs_form)

return A, b

A, b = assemble_system(fem.form(a), fem.form(L), bcs)
'''assert nsp.test(A)
A.setNullSpace(nsp)'''

# Create vector that spans the null space
nullspace = PETSc.NullSpace().create(constant=True,comm=MPI.COMM_WORLD)
assert nullspace.test(A)
A.setNullSpace(nullspace)

# orthogonalize b with respect to the nullspace ensures that
# b does not contain any component in the nullspace
nullspace.remove(b)

######### KRILOV SOLVER
ksp = PETSc.KSP().create(MPI.COMM_WORLD)
ksp.setOperators(A)
ksp.setType("cg")
ksp.getPC().setType("asm")
ksp.setTolerances(rtol=1e-10)
# Monitor the convergence of the KSP
ksp.setFromOptions()

uh = fem.Function(V)

# Not needed in this case since pure neumann
'''with b.localForm() as loc:
loc.set(0)
fem.petsc.assemble_vector(b, fem.form(L))
fem.petsc.apply_lifting(b, [fem.form(a)], [bcs])
fem.petsc.set_bc(b, bcs)'''

ksp.solve(b, uh.vector)

'''problem = fem.petsc.LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "cg", "pc_type": "asm", "ksp_monitor": None, "ksp_view":None, "ksp_rtol":1e-6, "ksp_atol":1e-10, "ksp_max_it": 1000})
uh = problem.solve()'''

V2 = fem.FunctionSpace(domain, ("CG", 3))
uex = fem.Function(V2)
uex.interpolate(u_exact)

L2_error = fem.form(ufl.inner(uh - uex, uh - uex) * ufl.dx)
error_local = fem.assemble_scalar(L2_error)
error_L2 = np.sqrt(domain.comm.allreduce(error_local, op=MPI.SUM))

if domain.comm.rank == 0:
print(f"Error_L2 : {error_L2:.2e}")

eh = uh - uex
E_H10 = np.sqrt(domain.comm.allreduce(fem.assemble_scalar(error_H10), op=MPI.SUM))
if domain.comm.rank == 0:
print(f"H01-error: {E_H10:.2e}")

from dolfinx.io import XDMFFile
with XDMFFile(domain.comm, "stokes_true/stoke.xdmf", "w") as file:
file.write_mesh(domain)
file.write_function(uh)
``````

You are using an iterative solver

Then you cannot remove the nullspace directly.

Im not sure if `asm` is suitable to solve singular problems efficiently. Why not use “boomeramg”
https://petsc.org/release/manual/ksp/#algebraic-multigrid-amg-preconditioners
Or

Thanks for the answer. Then, which are the right commands to do this? I also took inspiration from the demo of the Stoke problem Stokes equations with Taylor-Hood elements — DOLFINx 0.5.1 documentation .
But, still I am struggling on understating the correct functioning.

Boomeramg is for instance used in: https://github.com/jorgensd/dolfinx_mpc/blob/360166711fff96c8ecdd3c1c0f0033ca5af9a28d/python/demos/demo_periodic_geometrical.py#L95
and setnearnullspace is used in: https://github.com/search?q=repo%3Ajorgensd%2Fdolfinx_mpc%20setnearnullspace&type=code

Thanks for the example you have suggested. I modify the code accordingly but still doesn’t work

``````# Create PETSc nullspace
nullspace = PETSc.NullSpace().create(constant=True)
PETSc.Mat.setNearNullSpace(A, nullspace)

######### KRILOV SOLVER
ksp = PETSc.KSP().create(MPI.COMM_WORLD)
ksp.setOperators(A)

ksp.setType("cg")

ksp.getPC().setType(PETSc.PC.Type.HYPRE)
ksp.getPC().setHYPREType("boomeramg")
ksp.setTolerances(rtol=1e-10)
# Monitor the convergence of the KSP
ksp.setFromOptions()

uh = fem.Function(V)

# Not needed in this case since pure neumann
'''with b.localForm() as loc:
loc.set(0)
fem.petsc.assemble_vector(b, fem.form(L))
fem.petsc.apply_lifting(b, [fem.form(a)], [bcs])
fem.petsc.set_bc(b, bcs)'''

ksp.solve(b, uh.vector)
``````

You should take a step back and make sure that the code does what you expect.
Consider the following example: Add singular Poisson · Issue #139 · jorgensd/dolfinx-tutorial · GitHub

``````# Copyright (C) 2023 Jørgen S. Dokken
#
# Solve a singular Poisson problem with a null-space
#

from dolfinx import mesh, fem, io
from pathlib import Path
import ufl
import numpy as np
from petsc4py import PETSc
from mpi4py import MPI
import argparse
from enum import Enum

class SolverMode(Enum):
direct = 1
iterative = 2

parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)
dest="N", type=int, default=160,
help="Number of elements in each direction")
dest="iterative", default=False,
help="Direct solver if not set, else iterative")
args = parser.parse_args()

N = args.N
iterative_solver = args.iterative

domain = mesh.create_unit_square(MPI.COMM_WORLD, N, N, mesh.CellType.quadrilateral)

V = fem.FunctionSpace(domain, ("Lagrange", 1))

def u_exact(x):
return np.cos(2*np.pi*x[0])

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

x = ufl.SpatialCoordinate(domain)
u_ex = ufl.cos(2*ufl.pi*x[0])
n = ufl.FacetNormal(domain)

L = ufl.inner(f,v)* ufl.dx + ufl.inner(g, v) * ufl.ds

A = fem.petsc.assemble_matrix(fem.form(a))
A.assemble()
A.setOption(PETSc.Mat.Option.SYMMETRIC, True)
A.setOption(PETSc.Mat.Option.SYMMETRY_ETERNAL, True)

b = fem.petsc.assemble_vector(fem.form(L))

# Create vector that spans the null space
nullspace = PETSc.NullSpace().create(constant=True,comm=MPI.COMM_WORLD)
assert nullspace.test(A)
if args.iterative:
A.setNearNullSpace(nullspace)
else:
A.setNullSpace(nullspace)
nullspace.remove(b)

ksp = PETSc.KSP().create(domain.comm)
ksp.setOperators(A)

opts = PETSc.Options()
if args.iterative:
opts["ksp_type"] = "gmres"
opts["ksp_rtol"] = 1.0e-10
opts["pc_type"] = "hypre"
opts['pc_hypre_type'] = 'boomeramg'
opts["pc_hypre_boomeramg_max_iter"] = 1
opts["pc_hypre_boomeramg_cycle_type"] = "v"
opts["mg_levels_ksp_type"] = "cg"
opts["mg_levels_pc_type"] = "jacobi"
opts["mg_levels_esteig_ksp_type"] = "cg"

else:
opts["ksp_type"] = "preonly"
opts["pc_type"] = "lu"
opts["pc_factor_solver_type"] = "mumps"
opts["mat_mumps_icntl_24"] = 1  # Option to support solving a singular matrix
opts["mat_mumps_icntl_25"] = 0  # Option to support solving a singular matrix

ksp.setFromOptions()

uh = fem.Function(V)

ksp.solve(b, uh.vector)
uh.x.scatter_forward()
iterations = ksp.getIterationNumber()
# See: https://petsc.org/release/manualpages/KSP/KSPConvergedReason/
converged_reason = ksp.getConvergedReason()
if converged_reason > 0:
print(f"PETSc solver converged in {iterations=}")
else:
raise RuntimeError(f"PETSc solver did not converge with {converged_reason=}")

ex_mean = domain.comm.allreduce(fem.assemble_scalar(fem.form(u_ex*ufl.dx)), op=MPI.SUM)
approx_mean = domain.comm.allreduce(fem.assemble_scalar(fem.form(uh*ufl.dx)), op=MPI.SUM)
uh.x.array[:] -= approx_mean - ex_mean

L2_error = fem.form(ufl.inner(uh - u_ex, uh - u_ex) * ufl.dx)
error_local = fem.assemble_scalar(L2_error)
error_L2 = np.sqrt(domain.comm.allreduce(error_local, op=MPI.SUM))

if domain.comm.rank == 0:
print(f"Error_L2 : {error_L2:.2e}")

eh = uh - u_ex
E_H10 = np.sqrt(domain.comm.allreduce(fem.assemble_scalar(error_H10), op=MPI.SUM))
if domain.comm.rank == 0:
print(f"H01-error: {E_H10:.2e}")

results = Path("results")
results.mkdir(exist_ok=True)
with io.XDMFFile(domain.comm, results / "solution.xdmf", "w") as file:
file.write_mesh(domain)
file.write_function(uh)

out_ex = fem.Expression(u_ex, V.element.interpolation_points())
exact = fem.Function(V)
exact.interpolate(out_ex)
with io.XDMFFile(domain.comm, results / "exact.xdmf", "w") as file:
file.write_mesh(domain)
file.write_function(exact)
``````

which gives you the correct solution for both iterative and direct solvers (direct solver on very coarse grids might fail).

2 Likes

Thanks a lot is really helpful, for what I saw the problem is my code was related to how I deal with the right hand side, true? Also, why this formulation is not working for the cg method?

I suggest you spend some time comparing the resulting vectors on the RHS, Where I use automatic differentiation for the RHS source term an Neumann condition.

You should also visualize the solution and the exact solution in Paraview, and see if the problem is determined up to a constant.

I’ve had some issues with CG on coarse grids for this problem.

Sorry to interrupt, I referenced your example to solve the thermal expansion of a cylinder at uniform temperature (boundary conditions: fixed displacement in the z-direction on the top and bottom faces); I used the null space and the direct solver settings you mentioned, and got the following output: Some settings were not used.

``````Iteration: 0, rel. residual: 58477.55547242168
Iteration: 1, rel. residual: 7.1058772642900544e-09
WARNING! There are options you set that were not used!
WARNING! could be spelling mistake, etc!
There are 3 unused database options. They are:
Option left: name:-mat_mumps_icntl_24 value: 1 source: code
Option left: name:-mat_mumps_icntl_25 value: 0 source: code
Option left: name:-pc_factor_solver_type value: mumps source: code
``````

Could you please explain what is going on? Below are my code and geo file:

``````import dolfinx
from dolfinx import la
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx.io import gmshio
from dolfinx import io, default_scalar_type
from ufl import (tr, inner, lhs, rhs, Identity, TestFunction, TrialFunction, dx, grad)
from dolfinx.fem import (functionspace, FunctionSpace, Function, form, Constant, dirichletbc, locate_dofs_topological)
from dolfinx.fem.petsc import (assemble_matrix, assemble_vector, apply_lifting, set_bc)

def null_space(V: FunctionSpace):
"""Build PETSc nullspace for 3D elasticity"""
# Create vectors that will span the nullspace
bs = V.dofmap.index_map_bs
length0 = V.dofmap.index_map.size_local
basis = [la.vector(V.dofmap.index_map, bs=bs, dtype=PETSc.ScalarType) for i in range(6)]
b = [b.array for b in basis]
# Get dof indices for each subspace (x, y and z dofs)
dofs = [V.sub(i).dofmap.list.flatten() for i in range(3)]
# Set the three translational rigid body modes
for i in range(3):
b[i][dofs[i]] = 1.0
# Set the three rotational rigid body modes
x = V.tabulate_dof_coordinates()
dofs_block = V.dofmap.list.flatten()
x0, x1, x2 = x[dofs_block, 0], x[dofs_block, 1], x[dofs_block, 2]
b[3][dofs[0]] = -x1
b[3][dofs[1]] = x0
b[4][dofs[0]] = x2
b[4][dofs[2]] = -x0
b[5][dofs[2]] = x1
b[5][dofs[1]] = -x2

_basis = [x._cpp_object for x in basis]
dolfinx.cpp.la.orthonormalize(_basis)
assert dolfinx.cpp.la.is_orthonormal(_basis)
basis_petsc = [
PETSc.Vec().createWithArray(x[: bs * length0], bsize=3, comm=V.mesh.comm)
for x in b
]
return PETSc.NullSpace().create(vectors=basis_petsc)
# Mesh
mesh, cell_markers, facet_markers = gmshio.read_from_msh("unitf.msh", MPI.COMM_WORLD, gdim=3)
# Parameters
E = 2e5
nu = 0.3
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
alpha = 1e-5
def eps(v):
def sigma(v, dT):
return (lmbda*tr(eps(v))-alpha*(3*lmbda + 2*mu)*dT)*Identity(3)+2.0*mu*eps(v)
#Function space
VU = functionspace(mesh, ("CG", 1, (mesh.geometry.dim, )))
#Boundary conditions:top and bottom, z = 0
U_D = Constant(mesh, default_scalar_type(0))
fdim = mesh.topology.dim - 1
facets_top = facet_markers.find(5)
dofs_top_z = locate_dofs_topological(VU.sub(2), fdim, facets_top)
bcU1 = dirichletbc(U_D, dofs_top_z, VU.sub(2))

facets_bottom = facet_markers.find(6)
dofs_bottom_z = locate_dofs_topological(VU.sub(2), fdim, facets_bottom)
bcU2 = dirichletbc(U_D, dofs_bottom_z, VU.sub(2))
bcU = [bcU1, bcU2]
# Variational formula
dU = TrialFunction(VU)
U_ = TestFunction(VU)
Wint = inner(sigma(dU, 1000), eps(U_)) * dx
aM = lhs(Wint)
LM = rhs(Wint)
# Assemble matrix
A = assemble_matrix(form(aM), bcs = bcU)
A.assemble()
b = assemble_vector(form(LM))
apply_lifting(b, [form(aM)], bcs=[bcU])
set_bc(b, bcU)
# Set solver options
opts = PETSc.Options()
opts["ksp_type"] = "preonly"
opts["pc_type"] = "lu"
opts["pc_factor_solver_type"] = "mumps"
opts["mat_mumps_icntl_24"] = 1  # Option to support solving a singular matrix
opts["mat_mumps_icntl_25"] = 0  # Option to support solving a singular matrix
#Create PETSc Krylov solver and turn convergence monitoring on
solver = PETSc.KSP().create(mesh.comm)
solver.setFromOptions()
solver.setOperators(A)

null_space = null_space(VU)
A.setNullSpace(null_space)
null_space.remove(b)
# Solve
U = Function(VU)
solver.setMonitor(lambda _, its, rnorm: print(f"Iteration: {its}, rel. residual: {rnorm}"))
solver.solve(b, U.vector)
# Save results
vtkFile = io.VTKFile(MPI.COMM_WORLD, "results/U.pvd", "w")
vtkFile.write_function(U)
vtkFile.close()
``````
``````// Gmsh project created on Tue May 21 09:56:22 2024
r1 = 3;
N = 31;
Circle(1) = {0, 0, 0, r1, 0, 2*Pi};
Transfinite Curve {1} = N Using Progression 1;
Curve Loop(1) = {1};
Plane Surface(1) = {1};
Extrude {0, 0, 10} {Surface{1};Layers {8};}
Physical Surface("cylinder", 4) = {2};
Physical Surface("top", 5) = {3};
Physical Surface("bottom", 6) = {1};
Physical Volume("v", 7) = {1};
``````

Dear @French_fries

There is a typo in my options, it should be:

``````opts = PETSc.Options()
opts["ksp_type"] = "preonly"
opts["pc_type"] = "lu"
opts["pc_factor_mat_solver_type"] = "mumps"
opts["mat_mumps_icntl_24"] = 1  # Option to support solving a singular matrix
opts["mat_mumps_icntl_25"] = 0  # Option to support solving a singular matrix

solver.setFromOptions()
``````

However, as an experienced user of FEniCS I would expect you to:

• Use a built in mesh for your posts, as the current example requires extra work from the user to reproduce the error.
• Tell us what version of DOLFINx you are using
• If an external mesh is necessary, tell us what version of GMSH you used, how you converted the geometry etc.

Thank you for your response and suggestions. I would like to know from which folder in dolfix0.8.0 I can check the correct spelling of these settings so that I can resolve this issue myself.

I am not sure what you mean by this reply.
These are PETSc options, thus one should check the PETSc documentation for the spelling: MATSOLVERMUMPS — PETSc 3.21.5 documentation

Thank you, this is exactly what I wanted.