My understanding is that the condition number (for stiffness matrix of poisson) should scale like
h^-2. Here is a MWE, more or less from the demos, where I run this for a few levels of refinement using nx = [10,20,40,80,…]. I do this for a 1D, 2D, and 3D unit mesh. The 1D and 2D scale no problem (by a factor of 4 for each halving of h), however in 3D it scales by a factor of 8… or h^-3. The smallest eigenvalue is shrinking faster than h^2 in 3D.
What am I missing? (I hope it is very obvious)
from fenics import *
from slepc4py import SLEPc
nx=10
#mesh = UnitIntervalMesh(nx)
#mesh = UnitSquareMesh(nx,nx)
mesh = UnitCubeMesh(nx,nx,nx)
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(v), grad(u)) * dx
bc = DirichletBC(V, Constant(1), "on_boundary")
A = PETScMatrix()
assemble(a, tensor=A)
bc.apply(A)
eigenSolver = SLEPcEigenSolver(A)
eigenSolver.parameters["spectrum"] = "smallest magnitude"
eigenSolver.solve(0) # compute one eigenvalue
min_eig = eigenSolver.get_eigenvalue(0)[0]
eigenSolver.parameters["spectrum"] = "largest magnitude"
eigenSolver.solve(0) # compute one eigenvalue
max_eig = eigenSolver.get_eigenvalue(0)[0]
print(max_eig/min_eig)
This is a consequence of strong boundary enforcement. For some reason it is apparent in 3D, your largest eigenvalue doesn’t scale like it should because of how your system is modified with bc.apply(A). Either perform the eigenvalue computation on a reduced system (remove row/columns of strong Dirichlet dofs) or apply Dirichlet conditions weakly.
It would likely benefit you to use assemble_system which preserves symmetry of the underlying operator when applying boundary conditions. Particularly for the efficient computation of eigenvalue approximations.
This ambiguity (unless you read the documentation in detail) is a “feature” of legacy DOLFIN, and one of many reasons I highly recommend using DOLFINx instead of legacy DOLFIN.
The issue lies here; (strong) BC enforcement in FEM zeros-out rows and columns in the stiffness matrix and sets a 1 at that spot in the diagonal. Consequently, there will be repeating eigenvalues of 1 in your spectrum.
While the condition number of the stiffness matrix does not depend on spatial dimension, the value of the smallest and largest eigenvalues do. (Not entirely sure at which rate, but they decrease with d). By the time you moved to 3D, the largest eigenvalue becomes so small, that the 1 dominates ..
A ‘fix’ in dolfinx is to alter the value set on the diagonal through A = assemble_matrix(fem.form(a), bcs=[bc], diag=some_value). To do that properly you’d have to scale that similarly to how the largest and smallers eigenvalues scale. Or you set it to zero, and use look for the smallest nonzero eigenvalue (see Exercise 7: Use of Deflation Subspaces — SLEPc 3.24.1 documentation)
Here is a dolfinx translated code snippet to play with:
from mpi4py import MPI
from dolfinx import fem, mesh
from dolfinx.fem.petsc import assemble_matrix
import ufl
from petsc4py import PETSc
from slepc4py import SLEPc
d = 2
nx = 10
diag = 0.1
match d:
case 1:
domain = mesh.create_unit_interval(MPI.COMM_WORLD, nx)
case 2:
domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, nx)
case 3:
domain = mesh.create_unit_cube(MPI.COMM_WORLD, nx, nx, nx)
V = fem.functionspace(domain, ("Lagrange", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(ufl.grad(v), ufl.grad(u)) * ufl.dx
# Define boundary condition
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, domain.topology.dim - 1, boundary_facets)
bc = fem.dirichletbc(PETSc.ScalarType(0), boundary_dofs, V)
# Assemble matrix
A = assemble_matrix(fem.form(a), bcs=[bc], diag=diag)
A.assemble()
# Solve for eigenvalues
E = SLEPc.EPS().create(MPI.COMM_WORLD)
E.setOperators(A)
# Find smallest magnitude eigenvalue
E.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_MAGNITUDE)
E.solve()
min_eig = E.getEigenpair(0).real
# Find largest magnitude eigenvalue
E.setWhichEigenpairs(SLEPc.EPS.Which.LARGEST_MAGNITUDE)
E.solve()
max_eig = E.getEigenpair(0).real
print(max_eig/min_eig, max_eig, min_eig)
To further supply on top of @Stein’s code and comments, you can also extract a sub-matrix from your PETSc matrix, and use that as input to your eigenvalue problem.
I here illustrate this with standard PETSc KSP, but that is just to show the procedure:
# # Extract a sub matrix without boundary condition dofs and use it with PETSc
# Author: Jørgen S. Dokken
# SPDX License identifier: MIT
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx.fem.petsc
import ufl
import numpy as np
M = 30
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, M, M)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
u, v= ufl.TrialFunction(V), ufl.TestFunction(V)
a = dolfinx.fem.form(ufl.inner(u, v) * ufl.dx)
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
exterior_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, exterior_facets)
bcs = [dolfinx.fem.dirichletbc(0., exterior_dofs, V)]
A = dolfinx.fem.petsc.assemble_matrix(a, bcs, diag=0.0)
A.assemble()
num_dofs_local = V.dofmap.index_map.size_local*V.dofmap.index_map_bs
num_ghosts = V.dofmap.index_map.num_ghosts*V.dofmap.index_map_bs
not_dirichlet = np.full(num_dofs_local+num_ghosts, 1, dtype=np.int32)
not_dirichlet[exterior_dofs] = 0
idx_set = np.flatnonzero(not_dirichlet).astype(np.int32)
idx_set_row = idx_set[idx_set<num_dofs_local]
idx_set_col = idx_set[idx_set<num_dofs_local]
isx_row = PETSc.IS(A.getComm()).createGeneral(idx_set_row)
isx_col = PETSc.IS(A.getComm()).createGeneral(idx_set_col)
lgm_row = A.getLGMap()[0].applyIS(isx_row)
lgm_col = A.getLGMap()[1].applyIS(isx_col)
A_sub = A.createSubMatrix(lgm_row, lgm_col)
A_sub.assemble()
assert np.isclose(A.norm(0), A_sub.norm(0))
assert np.isclose(A.norm(2), A_sub.norm(2))
assert np.isclose(A.norm(3), A_sub.norm(3))
x = ufl.SpatialCoordinate(mesh)
L = dolfinx.fem.form(ufl.inner(x[0]+2*x[1], v)*ufl.dx)
b = dolfinx.fem.petsc.assemble_vector(L)
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
b_sub = b.getSubVector(lgm_row)
ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A_sub)
ksp.setType("preonly")
ksp.setErrorIfNotConverged(True)
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
ksp.getPC().getFactorMatrix().setMumpsIcntl(14, 1000)
x= b.duplicate()
x_sub = x.getSubVector(lgm_row)
ksp.solve(b_sub, x_sub)
x.restoreSubVector(lgm_row, x_sub)
uh = dolfinx.fem.Function(V)
x.copy(uh.x.petsc_vec)
uh.x.scatter_forward()
with dolfinx.io.VTXWriter(uh.function_space.mesh.comm, "u_sub.bp", [uh]) as bp:
bp.write(0.0)
As I am using the PETSc interface to extract the matrix, you could always achieve this in legacy fenics as well.
The symmetric approach of applying bcs in legacy happens through assemble_system.
This would give you an equivalent matrix to A in my post above (except for the diagonal values).
Then, you would have to create the not_dirichlet marker for the index set, i.e along the lines of:
from dolfin import *
import numpy as np
mesh = UnitSquareMesh(10, 10)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx
L = inner(Constant((0,0)), v)*dx
bc = DirichletBC(V, (0.0, 0.0), "on_boundary")
exterior_dofs = np.array(list(bc.get_boundary_values().keys()))
im = V.dofmap().index_map()
num_dofs_local = im.size(im.MapSize.OWNED)
num_ghosts = im.size(im.MapSize.UNOWNED)
not_dirichlet = np.full(num_dofs_local+num_ghosts, 1, dtype=np.int32)
not_dirichlet[exterior_dofs] = 0
idx_set = np.flatnonzero(not_dirichlet).astype(np.int32)
print(mesh.mpi_comm().rank, idx_set)
The rest should follow by assembling A, and getting the PETSc matrix with: