This page was generated from unit-2.4-Maxwell/Maxwellevp.ipynb.

2.4.1 Maxwell eigenvalue problem

We solve the Maxwell eigenvalue problem

\[\int \operatorname{curl} u \, \operatorname{curl} v = \lambda \int u v\]

for \(u, v \; \bot \; \nabla H^1\) using a PINVIT solver from the ngsolve solvers module.

The orthogonality to gradient fields is important to eliminate the huge number of zero eigenvalues. The orthogonal sub-space is implemented using a Poisson projection:

\[P u = u - \nabla \Delta^{-1} \operatorname{div} u\]

The algorithm and example is take form the Phd thesis of Sabine Zaglmayr, p 145-150.

import netgen.gui
from netgen.csg import *
from ngsolve import *

geom = CSGeometry()

cube1 = OrthoBrick(Pnt(-1,-1,-1), Pnt(1,1,1))
cube2 = OrthoBrick(Pnt(0,0,0), Pnt(2,2,2))
fichera = cube1-cube2
geom.Add (fichera)

# all edges defined by intersection of faces of
# cube2 and cube2 are marked for geometric edge refinement
geom.SingularEdge (cube2,cube2, 1)
geom.SingularPoint (cube2,cube2,cube2, 1)

mesh = Mesh(geom.GenerateMesh(maxh=0.5))
# two levels of hp-refinement
mesh.RefineHP(2, factor=0.2)
Draw (mesh)

SetHeapSize(100*1000*1000)

fes = HCurl(mesh, order=3)
print ("ndof =", fes.ndof)
u,v = fes.TnT()

a = BilinearForm(fes)
a += curl(u)*curl(v)*dx

m = BilinearForm(fes)
m += u*v*dx

apre = BilinearForm(fes)
apre += curl(u)*curl(v)*dx + u*v*dx
pre = Preconditioner(apre, "direct", inverse="sparsecholesky")

ndof = 30460

with TaskManager():
    a.Assemble()
    m.Assemble()
    apre.Assemble()

    # build gradient matrix as sparse matrix (and corresponding scalar FESpace)
    gradmat, fesh1 = fes.CreateGradient()


    gradmattrans = gradmat.CreateTranspose() # transpose sparse matrix
    math1 = gradmattrans @ m.mat @ gradmat   # multiply matrices
    math1[0,0] += 1     # fix the 1-dim kernel
    invh1 = math1.Inverse(inverse="sparsecholesky")

    # build the Poisson projector with operator Algebra:
    proj = IdentityMatrix() - gradmat @ invh1 @ gradmattrans @ m.mat

    projpre = proj @ pre.mat

    evals, evecs = solvers.PINVIT(a.mat, m.mat, pre=projpre, num=12, maxit=20)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-3-b9c9cf11f0eb> in <module>
     18     projpre = proj @ pre.mat
     19
---> 20     evals, evecs = solvers.PINVIT(a.mat, m.mat, pre=projpre, num=12, maxit=20)

/usr/lib/python3/dist-packages/ngsolve/eigenvalues.py in PINVIT(mata, matm, pre, num, maxit, printrates)
     47                 msmall[j,k] = InnerProduct(Mv, vecs[k])
     48
---> 49         ev,evec = scipy.linalg.eigh(a=asmall, b=msmall)
     50         lams[:] = ev[0:num]
     51         if printrates:

/usr/local/lib/python3.6/dist-packages/scipy/linalg/decomp.py in eigh(a, b, lower, eigvals_only, overwrite_a, overwrite_b, turbo, eigvals, type, check_finite)
    372
    373     """
--> 374     a1 = _asarray_validated(a, check_finite=check_finite)
    375     if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
    376         raise ValueError('expected square matrix')

/usr/local/lib/python3.6/dist-packages/scipy/_lib/_util.py in _asarray_validated(a, check_finite, sparse_ok, objects_ok, mask_ok, as_inexact)
    237             raise ValueError('masked arrays are not supported')
    238     toarray = np.asarray_chkfinite if check_finite else np.asarray
--> 239     a = toarray(a)
    240     if not objects_ok:
    241         if a.dtype is np.dtype('O'):

/usr/lib/python3/dist-packages/numpy/lib/function_base.py in asarray_chkfinite(a, dtype, order)
   1213     if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all():
   1214         raise ValueError(
-> 1215             "array must not contain infs or NaNs")
   1216     return a
   1217

ValueError: array must not contain infs or NaNs

print ("Eigenvalues")
for lam in evals:
    print (lam)

Eigenvalues

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-4-22ead1bc97a9> in <module>
      1 print ("Eigenvalues")
----> 2 for lam in evals:
      3     print (lam)

NameError: name 'evals' is not defined

Draw (mesh)
gfu = GridFunction(fes, multidim=len(evecs))
for i in range(len(evecs)):
    gfu.vecs[i].data = evecs[i]
Draw (Norm(gfu), mesh, "u", sd=4)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-5-77fa5933e2a2> in <module>
      1 Draw (mesh)
----> 2 gfu = GridFunction(fes, multidim=len(evecs))
      3 for i in range(len(evecs)):
      4     gfu.vecs[i].data = evecs[i]
      5 Draw (Norm(gfu), mesh, "u", sd=4)

NameError: name 'evecs' is not defined