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])
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
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])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
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
error_H10 = fem.form(ufl.inner(ufl.grad(eh), ufl.grad(eh)) * ufl.dx)
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)