Preconditioned Inverse Iteration for Laplace Eigenvalue Problems
We will use the Preconditioned INVerse ITeration (PINVIT) scheme, developed by Knyazef and Neymeyr, to compute the lowest eigenvalues of the Dirichlet Laplacian. We seek \((u,\lambda) \in H^1_0(\Omega)\times \mathbb{R}\) such that for any \(v\in H^1_0(\Omega)\)
\[\int_\Omega \nabla u \cdot \nabla v \; d\vec{x} = \lambda \int_\Omega uv\;d\vec{x}.\]
In the process of coding the PINVIT scheme, we will show how to use the VectorMapping
class to
map PETSc vectors to NGSolve vectors and vice versa. We will also show how to use the Matrix
class
to create a PETSc matrix from an NGSolve BilinearForm
.
First, we need to construct the distributed mesh for the finite element space discretising the PDE.
from ngsolve import Mesh, unit_square
import netgen.meshing as ngm
import numpy as np
from scipy.linalg import eigh
from mpi4py.MPI import COMM_WORLD
mesh = Mesh(unit_square.GenerateMesh(maxh=0.2, comm=COMM_WORLD))
We now proceed to construct a linear polynomial finite element space, with \(H^1\) conformity, and discretise the mass matrix that represents the \(L^2\) inner product. We also fetch a compatible vector from the mass matrix.
from ngsolve import H1, BilinearForm, dx
fes = H1(mesh, order=1, dirichlet="left|right|top|bottom")
u,v = fes.TnT()
m = BilinearForm(u*v*dx).Assemble()
M = m.mat
ngsVec = M.CreateColVector()
We are now ready to create a VectorMapping
, which will first be used to construct
a PETSc vector from the ngsVec
just initialized.
The only information that the VectorMapping
class needs is the finite element space
associated to the vector GridFunction
. This is because the NGSolve FESpace
class
contains information about the way the degrees of freedom are distributed and which degrees of
freedom are not constrained by the boundary conditions.
from ngsPETSc import VectorMapping
Map = VectorMapping(fes)
petscVec = Map.petscVec(ngsVec)
print("Vector type is {} and it has size {}.".format(petscVec.type,petscVec.size))
We now use the Matrix
class to create a PETSc matrix from an NGSolve BilinearForm
.
Once the Matrix
class has been set up, it is possible to access the corresponding PETSc matrix
object as Matrix().mat
. By default, if COMM_WORLD.GetSize()
is larger than one,
mat
is initialized as a PETSc mpiaij
which is the default sparse parallel matrix in PETSc,
otherwise mat
is initialized as a PETSc seqaij
which is the default serial matrix in PETSc.
We can also view the matrix using the Matrix().view()
method.
from ngsPETSc import Matrix
M = Matrix(m.mat, fes)
print("Matrix type is {} and it has size {}.".format(M.mat.type,M.mat.size))
M.view()
There are other matrix formats available. To mention a few:
dense
, a dense format,cusparse
, CUDA sparse format, for NVIDIA GPUs.aijmkl
, Intel MKL format.
We solve the discretised problem by looking for the eigenvalues of the generalised eigenproblem \(A\vec{u}_h = \lambda M\vec{u}_h\) where \(A\) and \(M\) are the finite element discretisations respectively of the stiffness matrix for the Laplacian and the mass matrix. Since we have already constructed the mass matrix, it remains to construct the stiffness matrix:
from ngsolve import grad, Preconditioner, GridFunction
a = BilinearForm(fes)
a += grad(u)*grad(v)*dx
a.Assemble()
A = Matrix(a.mat, fes)
At the heart of the PINVIT scheme there is an iteration similar to the Rayleigh quotient iteration for a generalised eigenvalue problem. More details on the latter can be found in Nick Trefethen’s Numerical Linear Algebra, Lecture 27:
\[\vec{u}_h^{(n+1)} = \omega_1^{(n)}\vec{u}_{h}^{(n)}+\omega_2^{(n)} \vec{\omega}_h^{(n)}, \qquad \vec{\omega}_h^{(n)}= P^{-1}(A\vec{u}_h^{(n)}-\rho_n M\vec{u}_h^{(n)}),\]
where \(P^{-1}\) is an approximate inverse of the stiffness matrix \(A\), \(\omega_i\) are step sizes, and \(\rho_n\) is the Rayleigh quotient corresponding to \(\vec{u}_h^{(n)}\), i.e.
\[\rho_{n} = \frac{(\vec{u}_h^{(n)}, A \vec{u}_h^{(n)})}{(\vec{u}_h^{(n)}, M\vec{u}_h^{(n)})}.\]
The choice of \(\omega_{1,2}^{(n)}\) is important to obtain a convergent PINVIT scheme. We do this by solving the optimization problem,
\[\vec{u}_h^{(n+1)} = \underset{\vec{v}\in <\vec{u}_h^{(n)},\, \vec{\omega}_h^{(n)}>}{arg\;min} \frac{(\vec{u}_h^{(n+1)}, A \vec{u}_h^{(n+1)})}{(\vec{u}_h^{(n+1)}, M\vec{u}_h^{(n+1)})}\]
which is equivalent to a small generalised eigenvalue problem, i.e.
\[\begin{split}\begin{bmatrix} \vec{u}_h^{(n)}\cdot A \vec{u}_h^{(n)} & \vec{u_h}^{(n)}\cdot A \vec{\omega}_h^{(n)}\\ \vec{\omega}_h^{(n)}\cdot A \vec{u}_h^{(n)} & \vec{\omega}_h^{(n)}\cdot A \vec{\omega}_h^{(n)} \end{bmatrix} = \omega \begin{bmatrix} \vec{u}_h^{(n)}\cdot M \vec{u}_h^{(n)} & \vec{u_h}^{(n)}\cdot M \vec{\omega}_h^{(n)}\\ \vec{\omega}_h^{(n)}\cdot M \vec{u}_h^{(n)} & \vec{\omega}_h^{(n)}\cdot M \vec{\omega}_h^{(n)} \end{bmatrix}.\end{split}\]
def stepChoice(Asc,Msc,w,u0):
Au0 = u0.duplicate(); Asc.mult(u0,Au0)
Mu0 = u0.duplicate(); Msc.mult(u0,Mu0)
Aw = w.duplicate(); Asc.mult(w,Aw)
Mw = w.duplicate(); Msc.mult(w,Mw)
smallA = np.array([[u0.dot(Au0),u0.dot(Aw)],[w.dot(Au0),w.dot(Aw)]])
smallM = np.array([[u0.dot(Mu0),u0.dot(Mw)],[w.dot(Mu0),w.dot(Mw)]])
_, evec = eigh(a=smallA, b=smallM)
return (float(evec[0,0]),float(evec[1,0]))
We then construct a PETSc preconditioner acting as an approximate inverse of \(A\). For this example, we use an algebraic multigrid preconditioner built with HYPRE.
from petsc4py import PETSc
pc = PETSc.PC()
pc.create(PETSc.COMM_WORLD)
pc.setOperators(A.mat)
pc.setType(PETSc.PC.Type.HYPRE)
pc.setUp()
We now implement the iteration itself:
from math import pi
itMax = 10
u0 = A.mat.createVecLeft()
w = A.mat.createVecLeft()
u0.setRandom()
for it in range(itMax):
Au0 = u0.duplicate(); A.mat.mult(u0,Au0)
Mu0 = u0.duplicate(); M.mat.mult(u0,Mu0)
rho = Au0.dot(u0)/Mu0.dot(u0)
print("[{}] Eigenvalue estimate: {}".format(it,rho/(pi**2)))
u = Au0+rho*Mu0
pc.apply(u,w)
alpha = stepChoice(A.mat,M.mat,w,u0)
u0 = alpha[0]*u0+alpha[1]*w
Iteration |
1 |
2 |
3 |
4 |
5 |
6 |
---|---|---|---|---|---|---|
Eigenvalue |
7.187387 |
2.172190 |
2.133502 |
2.131257 |
2.130477 |
2.130475 |