Add singular Poisson

[jorgensd]

Minimal example

# Copyright (C) 2023 Jørgen S. Dokken
#
# Solve a singular Poisson problem with a null-space
#
# SPDX-License-Identifier:   MIT


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)
parser.add_argument("--N",
                    dest="N", type=int, default=160,
                    help="Number of elements in each direction")
parser.add_argument("--iterative", action="store_true",
                    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])
f = -ufl.div(ufl.grad(u_ex))
n = ufl.FacetNormal(domain)
g = -ufl.dot(ufl.grad(u_ex), n)

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

# 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
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}")

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)