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.
[1]:
from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *
[2]:
from netgen.occ import *
cube1 = Box( (-1,-1,-1), (1,1,1) )
cube2 = Box( (0,0,0), (2,2,2) )
cube2.edges.hpref=1 # mark edges for geometric refinement
fichera = cube1-cube2
Draw (fichera);
[3]:
mesh = Mesh(OCCGeometry(fichera).GenerateMesh(maxh=0.4))
mesh.RefineHP(levels=2, factor=0.2)
Draw (mesh);
[4]:
# SetHeapSize(100*1000*1000)
fes = HCurl(mesh, order=3)
print ("ndof =", fes.ndof)
u,v = fes.TnT()
a = BilinearForm(curl(u)*curl(v)*dx)
m = BilinearForm(u*v*dx)
apre = BilinearForm(curl(u)*curl(v)*dx + u*v*dx)
pre = Preconditioner(apre, "direct", inverse="sparsecholesky")
ndof = 43192
[5]:
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)
/usr/lib/python3/dist-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.17.3 and <1.25.0 is required for this version of SciPy (detected version 1.26.4
warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
0 : [8.213129803313668, 17.79981324070466, 37.21982565708564, 44.01930033879474, 62.182333371265315, 65.53333262120249, 73.98716307328938, 81.08231276495654, 88.0600174829129, 91.25055695388384, 113.61980589985977, 122.21673394441416]
1 : [3.238097539135088, 5.998650630532796, 6.1613900721503585, 11.485860551493053, 12.037309680460307, 14.413520499586948, 15.234814500677576, 16.218269369624572, 17.346726952710153, 19.892283926601003, 22.638813222541838, 24.807560946532675]
2 : [3.220542688575404, 5.884876708899643, 5.888226095526421, 10.78547163214916, 10.847156783190373, 11.381986377508206, 13.343238319621873, 13.372706506711085, 14.11239328549126, 14.988639563908848, 15.344080016298419, 17.688329573541285]
3 : [3.2202954145601392, 5.880713446554482, 5.880839526591148, 10.707466526323545, 10.720069178390245, 10.944773310661384, 12.52797395497305, 13.075036550544898, 13.54069791575316, 14.030695748548794, 14.367494218084902, 15.906741581277407]
4 : [3.2202898941021925, 5.880502253291022, 5.880518097321811, 10.695114363078554, 10.699129958353117, 10.78775181154299, 12.354098370741875, 12.871053038516205, 13.445683487738647, 13.649403154487578, 14.213859959633721, 14.607665446432952]
5 : [3.2202897664212085, 5.8804933752546305, 5.8804959631163785, 10.692776367605212, 10.694999509694995, 10.721906317300862, 12.324202578887299, 12.640658945187862, 13.417157954166987, 13.50976221501871, 13.816759275477924, 14.260863114887776]
6 : [3.220289763807024, 5.880492978117827, 5.880494510444945, 10.691097181806363, 10.694161153653575, 10.698643594360192, 12.318883281766663, 12.460746604864875, 13.402394959799835, 13.454795875823102, 13.563293262896224, 14.232709516083279]
7 : [3.220289763759514, 5.880492958816468, 5.880494425037962, 10.687861627675979, 10.693985793779897, 10.6945527922947, 12.317899571105768, 12.372364426482001, 13.387450842643357, 13.43393941284152, 13.47327115795089, 14.223431358552903]
8 : [3.2202897637587014, 5.88049295783441, 5.880494420312168, 10.686464368101152, 10.693942498545187, 10.694085266048571, 12.31771105532961, 12.337392167602026, 13.373735427603892, 13.426640065252949, 13.441852510548186, 14.21926782630824]
9 : [3.220289763758684, 5.880492957782818, 5.880494420056987, 10.686108698962045, 10.693925371312586, 10.694004893999137, 12.317673740214623, 12.324680164902114, 13.365194188468868, 13.423921247805147, 13.43032072972929, 14.217191023860392]
10 : [3.2202897637586854, 5.88049295778006, 5.880494420043331, 10.68602521652086, 10.693919077435666, 10.693989757202933, 12.317666113412663, 12.320210206466525, 13.36087191675219, 13.42281682314765, 13.42565047966931, 14.21610359479629]
11 : [3.2202897637586827, 5.880492957779923, 5.880494420042617, 10.686006011017094, 10.693917414575703, 10.693986570954786, 12.31766449252742, 12.3186591999368, 13.358849289698323, 13.42234403795897, 13.423640049272196, 14.215520283663354]
12 : [3.2202897637586876, 5.880492957779889, 5.880494420042566, 10.686001612140265, 10.693917021614617, 10.693985859209521, 12.317664119922945, 12.318123649950923, 13.357930686091963, 13.422136540427969, 13.422748124706002, 14.21520352351913]
13 : [3.220289763758683, 5.880492957779889, 5.880494420042579, 10.686000605589273, 10.693916931075941, 10.693985697538217, 12.317664014900103, 12.317939026914644, 13.35751849858082, 13.422044243842372, 13.422347260146381, 14.215030378524819]
14 : [3.2202897637586823, 5.880492957779904, 5.880494420042565, 10.686000375000118, 10.693916910324218, 10.693985660615246, 12.317663977198508, 12.317875404934874, 13.357334525795956, 13.422003000662194, 13.42216585242637, 14.214935403161203]
15 : [3.2202897637586863, 5.880492957779933, 5.880494420042663, 10.686000322147715, 10.693916905559943, 10.693985652167179, 12.317663963003234, 12.3178534464673, 13.357252468936219, 13.421984373846442, 13.422083623534926, 14.214883194528081]
16 : [3.220289763758688, 5.88049295778005, 5.88049442004279, 10.686000310001575, 10.693916904465087, 10.69398565022954, 12.317663957907362, 12.317845855013298, 13.357215852581401, 13.421975918987902, 13.422046245660857, 14.21485442786466]
17 : [3.2202897637586867, 5.880492957779717, 5.880494420042691, 10.686000307205976, 10.6939169042128, 10.693985649783734, 12.317663956137796, 12.317843220978999, 13.357199456799544, 13.421972054519504, 13.422029254033875, 14.214838538198748]
18 : [3.220289763758684, 5.880492957779889, 5.880494420042589, 10.68600030655901, 10.693916904154731, 10.693985649680199, 12.317663955516236, 12.317842306796571, 13.35719207272615, 13.421970312167959, 13.422021426870389, 14.214829573551489]
19 : [3.220289763758687, 5.880492957779878, 5.880494420042545, 10.686000306412726, 10.69391690414138, 10.693985649657034, 12.317663955306264, 12.317841993507882, 13.357188798065334, 13.42196951379168, 13.422017862652558, 14.214824608612515]
[6]:
print ("Eigenvalues")
for lam in evals:
print (lam)
Eigenvalues
3.220289763758687
5.880492957779878
5.880494420042545
10.686000306412726
10.69391690414138
10.693985649657034
12.317663955306264
12.317841993507882
13.357188798065334
13.42196951379168
13.422017862652558
14.214824608612515
[7]:
gfu = GridFunction(fes, multidim=len(evecs))
for i in range(len(evecs)):
gfu.vecs[i].data = evecs[i]
Draw (Norm(gfu), mesh, order=2, min=0, max=2);
[ ]:
[ ]: