3.8. Solving the curlcurl eigenvalue problem using explicit methods

by M. Wess, 2025

This Notebook is part of the dualcellspaces documentation for the addon package implementing the Dual Cell method in NGSolve, as well as part of the td_evp documentation on the implementation of time-domain methods for resonance problems.

We show how to numerically solve the Maxwell eigenvalue problem combining dualcellspaces package with the ideas from

[1] L. Nannen and M. Wess, A Krylov Eigenvalue Solver Based on Filtered Time Domain Solutions, arXiv prepreint, 2024, https://doi.org/10.48550/arXiv.2402.08515.

Moreover we compare the time domain method to LOBPCG.

Problem setting

We want to approximate the eigenvectors and eigenfunctions of the problem to find \(\lambda\in\mathbb R, E\in H(\mathrm{curl})(\Omega)\) for some domain \(\Omega\subset\mathbb R^3\) such that

\[ \mathrm{curl}\mathrm{curl}E(x) = \tilde\lambda E(x) \]

with natural boundary conditions.

Setting up the discrete Problem

We start by doing the necessary imports

from ngsolve import *
import dualcellspaces as dcs
from ngsolve.webgui import Draw
import numpy as np
from netgen.occ import *
import matplotlib.pyplot as pl
import scipy.linalg as spl
from ngsolve.solvers import LOBPCG

Choosing \(\Omega:=[0,\pi]^3\) we generate our geometry and mesh:

maxh = 0.8
order = 2
lx = ly = lz = np.pi

geo = OCCGeometry(Box((0,0,0),(lx,ly,lz)))
mesh = Mesh(geo.GenerateMesh(maxh=maxh))
Draw(mesh)

WebGLScene

fes_E = dcs.HCurlDualCells3D(mesh, order=order)
fes_H = dcs.HCurlPrimalCells3D(mesh,order=order)

fes = fes_E*fes_H
print("total DoFs:",fes.ndof)
print("DoFs E:",fes_E.ndof)
gf = GridFunction(fes)
gfE, gfH = gf.components


#integral symbols with special integration rules
dxE = dx(intrules=fes_E.GetIntegrationRules())
dxH = dx(intrules=fes_H.GetIntegrationRules())

dxw = dx(intrules=dcs.GetIntegrationRules(2*order+4))
dSw = dx(element_boundary=True,intrules=dcs.GetIntegrationRules(2*order+4))


#mixed bilinear form
E,dE = fes_E.TnT()
H,dH = fes_H.TnT()

normal = specialcf.normal(3)
Curl = BilinearForm(E*curl(dH)*dxw+E*Cross(dH,normal)*dSw, geom_free=True).Assemble().mat

SetHeapSize(int(2e9))
#prepare mass operators

with TaskManager():
    massH_inv = fes_H.Mass(1).Inverse()
    massE = fes_E.Mass(1)
    massE_inv = fes_E.Mass(1).Inverse(freedofs=fes_E.FreeDofs())

total DoFs: 144522
DoFs E: 73614

Solving the problem with LOBPCG

To use LOBPCG we need to deal with the kernel of the \(\mathrm{curl}\)-operator. The dualcellspaces package provides the fitting \(H^1\)-space and \(\nabla\)-operator respecting the discrete exact sequence properties.

fespot = fes_E.GetPotentialSpace()
phi, dphi = fespot.TnT()
Grad = BilinearForm(grad(phi)*dE*dxE).Assemble()
gmat = Grad.mat.DeleteZeroElements(1e-12)

massPot = BilinearForm(phi*dphi*dxw, diagonal=True).Assemble()
with TaskManager():
    massPot_inv = massPot.mat.Inverse()

ccdd = Curl.T @ massH_inv @ Curl + 5000*gmat@massPot_inv@gmat.CreateTranspose()

Now we solve the stabilized problem using LOBPCG. As a preconditioner we just use the inverse of the \(E\)-mass matrix, which leads to very slow convergence. Thus we cheat a little by setting the initial vectors to projections of analytical eigenfunctions.

nevs = 5
initvecs = GridFunction(fes_E,multidim=nevs)

for i in range(nevs):
    initvecs.vecs[i].SetRandom()


initvecs.vecs[0].data = massE_inv*LinearForm(dE*CF((sin(x)*cos(y),-cos(x)*sin(y),0))*dxE).Assemble().vec
initvecs.vecs[1].data = massE_inv*LinearForm(dE*CF((0,sin(y)*cos(z),-cos(y)*sin(z)))*dxE).Assemble().vec
initvecs.vecs[2].data = massE_inv*LinearForm(dE*CF((sin(x)*cos(z),0,-cos(x)*sin(z)))*dxE).Assemble().vec
initvecs.vecs[3].data = massE_inv*LinearForm(dE*CF((2*sin(x)*cos(y)*cos(z),-cos(x)*sin(y)*cos(z),-cos(x)*cos(y)*sin(z)))*dxE).Assemble().vec

ccdd = Curl.T @ massH_inv @ Curl + 5000*gmat@massPot_inv@gmat.CreateTranspose()   

with TaskManager():
    lams,vecs = LOBPCG(ccdd, massE, pre=massE_inv, maxit=100, initial = initvecs.vecs,
                       printrates = False)
print("approximated evs: ",lams)

approximated evs:   2.00015
 2.00017
 2.00021
  3.0004
 858.118

The filtered problem

The idea from [1] in short is to use the solution of the time-domain Maxwell problem

\[\begin{split} \begin{aligned} \partial_t^2 e(t,x) + \mathrm{curl}\mathrm{curl}e(t,x)&=0\\ e(0,x)&=E_0(x) \end{aligned} \end{split}\]

To define the filtered operator \(\tilde \Pi_\alpha\) by $\( (\tilde \Pi_{\alpha} E^0)(x) = \int_0^\infty e(t,x)\alpha(t)dt \)\( for some weight function \)\alpha\( with final support in \)\mathbb R$.

It can be shown that for an eigenpair \((\tilde\omega, E)\) of the Maxwell eigenvalue problem \((\beta(\tilde\omega),E)\) is an eigenpair of \(\Pi_\alpha\) for

\[ \tilde \beta(s):=\int_0^\infty \cos(ts)\alpha(t)dt. \]

Discretizations

The same arguments hold true if we replace the continuous problem by the discrete eigenvalue problem

\[ \mathbf S \mathbf E = \lambda\mathbf M \mathbf E \]

where \(\mathbf S,\mathbf M\) are the discretization of the \(\mathrm {curl}\mathrm {curl}\) and the mass matrix and \(\mathbf E\) is the coefficient vector of the discrete eigenfunctions.

Discretizing the semi-discrete time domain problem by a time-stepping and the integral by a trapezoidal rule leads to the discrete filter

\[ (\Pi_\alpha)(\mathbf{e^0})=\tau \sum_{j=0}^M \alpha(\tau j) \mathbf e^j \]

where \(\mathbf e^j\) is the time-domain solution at time step number \(j\).

Then if \(\lambda,\mathbf E\) is an eigenpair of the discrete problem \((\beta(\lambda),\mathbf E)\) is an eigenpair of the operator \(\Pi_\alpha\) with

\[ \beta(s):= \tau \sum_{j=0}^M q_j(s)\alpha(\tau j) \]

and \(q_j\) is the approximation of the scalar problem

\[ \partial_t^2 l(t)+s l(t)=0,\quad l(0) = 1 \]

by the chosen time-stepping.

Choosing the filter and solving the problem

The idea is now to choose the weight function \(\alpha\) such that the resulting filter \(\beta\) is maximal in a region where the sought after eigenvalues are and close to zero otherwise. Then applying e.g., an Arnoldi algorithm to the filtered operator \(\Pi_\alpha\) will yield approximations to the eigenfunctions where \(\beta\) is maximal, which we know to be eigenvectors of the original problems corresponding to eigenvalues in the interesting region. Lastly we solve the original problem projected onto the space spanned by the Krylov space of \(\Pi_\alpha\).

We define the filtered operator as a class. Since we use a mixed method for the spacial discretization we need to supply the (inverse) mass matrices for the two fieds, as well as the discrete curl operator and the parameters for the time-discretization and filtering.

class FilteredC(BaseMatrix):
    def __init__(self, mata, matm_inv, tau, weights, freedofs = None):
        super().__init__()
        self.dt = tau
        self.weights = weights
        self.nsteps = len(weights)

        self.mata = mata
        self.matm_inv = matm_inv

        self.freedofs=freedofs


        self.vecu = self.mata.CreateColVector()
        self.tmpvec1 = self.mata.CreateColVector()
        self.tmpvec2 = self.mata.CreateColVector()


    def CreateColVector(self):
        return self.mata.CreateColVector()

    def Shape(self):
        return self.mata.shape

    def CreateVector(self,col):
        return self.mata.CreateVector(col)

    def Mult(self,rhs,out):
        with TaskManager():
            self.vecu.data = rhs
            tau = self.dt
            out.data = tau*self.weights[0]*self.vecu
            t = 0

            unew = self.tmpvec1
            uold = self.tmpvec2
            uold.data = self.vecu

            for i in range(1,self.nsteps):
                t += tau
                #print("\r time = {}, step = {}".format(t,i),end="")

                unew.data = 2*self.vecu - uold
                unew.data -= tau**2 * self.matm_inv@self.mata * self.vecu
                if self.freedofs:
                    unew.data[~self.freedofs] = 0.
                uold.data = self.vecu
                self.vecu.data = unew.data

                out.data += tau*self.weights[i]*self.vecu

To obtain a stable method we need to respect the CFL condition of the time-domain problem. Thus we apply a simple power iteration to determine the largest magnitude eigenvalue of the discrete problem:

massinv = massE_inv
mass = massE
stiffness = Curl.T@massH_inv@Curl
tol = 1e-4
tmpv = stiffness.CreateVector()
tmpv2 = stiffness.CreateVector()
tmpv.SetRandom()
tmpv /= tmpv.Norm()
mus = [1]
for i in range(500):
    #print("i = {}".format(i))
    tmpv2.data = tmpv
    tmpv.data = massinv@stiffness*tmpv2
    mu = InnerProduct(tmpv2,tmpv)/InnerProduct(tmpv2,tmpv2)
    tmpv /= tmpv.Norm()
    mus.append(mu)
    if mus[-1]/mus[-2]<1+tol:
        break

lammax = mus[-1]*(1+tol*200) #some safety factors
print("estimated largest eigenvalue: ",lammax)
tau = np.sqrt(4/lammax)
print("by power iteration stable for tau = {}".format(tau))

estimated largest eigenvalue:  10978.779644127553
by power iteration stable for tau = 0.019087671908812852

Next we need to pick suitable parameters for the time-domain filtering. We look for eigenvalues \(\lambda\in[1,4]\) and plot the resulting discrete filter.

# parameters for the filter function
lam_min = 1
lam_max = 4
endT = 10
L = int(endT/tau)


w_min = np.sqrt(lam_min)
w_max = np.sqrt(lam_max)

weightf = lambda t: 4/np.pi*np.cos((w_max+w_min)/2*t)*(w_max-w_min)/2*np.sinc((w_max-w_min)/2*t/np.pi)

weights = weightf(tau*np.arange(L))


def beta(lam):
    if np.isscalar(lam):
        q = 1
    else:
        q = np.ones(lam.shape)

    q_old = q
    out = tau*weights[0]*q
    for alpha in weights[1:]:
        q_new = 2*q-tau**2*lam*q-q_old
        q_old = q
        q = q_new
        out += tau*alpha*q
    return out

pl.figure()
pl.plot(np.arange(0,10,0.01),beta(np.arange(0,10,0.01)));
pl.xlim((0,10))

(0.0, 10.0)

maxsteps = 15

C = FilteredC(stiffness,massinv,tau,weights)

solveevery = 1

tmpvec = massinv.CreateVector()
tmpvec.SetRandom()
tmpvec.data = 1/tmpvec.Norm()*tmpvec
tmpvec2 = massinv.CreateVector()

vecs = MultiVector(tmpvec,0)
vecs.Append(tmpvec)
for i in range(maxsteps):
    print("\n Krylowstep = {}".format(i))

    tmpvec.data = C*vecs[-1]
    vecs.AppendOrthogonalize(tmpvec)

    if i%solveevery == 0:
        tvecs = MultiVector(tmpvec,len(vecs))
        tvecs.data = stiffness*vecs
        tvecs2 = MultiVector(tmpvec,len(vecs))
        tvecs2.data = mass*vecs

        matm_proj = InnerProduct(tvecs2,vecs)
        mats_proj = InnerProduct(tvecs,vecs)
        lam,v = spl.eigh(mats_proj.NumPy(),matm_proj.NumPy())

        #compute residuals
        eigf =(vecs*Matrix(v.real)).Evaluate()
        res = []
        tmpvec2[:] = 0
        for i in range(len(lam)):
            tmpvec2.data = massinv@stiffness*eigf[i]
            tmpvec2.data -= lam[i]*eigf[i]
            res.append(tmpvec2.Norm()/eigf[i].Norm())

print("approximated eigenvalues: ", lam)
print("residuals: ",res)

 Krylowstep = 0

 Krylowstep = 1

 Krylowstep = 2

 Krylowstep = 3

 Krylowstep = 4

 Krylowstep = 5

 Krylowstep = 6

 Krylowstep = 7

 Krylowstep = 8

 Krylowstep = 9

 Krylowstep = 10

 Krylowstep = 11

 Krylowstep = 12

 Krylowstep = 13

 Krylowstep = 14

approximated eigenvalues:  [5.16333662e-08 2.00009114e+00 2.00011540e+00 3.00024317e+00
 3.00024740e+00 5.00101566e+00 5.99834381e+00 8.89790511e+00
 1.10063602e+01 1.18263189e+01 2.23031012e+01 3.01821577e+01
 6.45826541e+01 2.06235882e+02 6.39703430e+02 1.12440071e+03]
residuals:  [0.0021394255329461064, 8.282502214398499e-06, 1.628315939636587e-05, 0.0021202456949313167, 0.0009724849514352617, 0.01912760191472676, 24.628820887078806, 1.3624106879472084, 0.4146277628135319, 17.760874441900715, 25.26322229292085, 49.280350189986564, 121.5324435480984, 333.16306979196423, 658.4167603303343, 1052.4785767586175]