Waveguide number of eigenmodes

For completeness sake.
Wave equation \nabla^2\vec{E}+k^2\vec{E}=0 has an infinite number of solutions which do not satisfy divergence equation \nabla\cdot\vec{E}=0. \nabla\times\nabla\times\vec{E}-k^2\vec{E}=0 on the other hand satisfy divergence-free condition automatically (as explained here). There’s not mathematically intensive explanation in Computational electromagnetics by David Davidson, chapter “The null space of the curl operator and spurious modes”.

PS code in the question is not entirely correct because it does not properly setup spectral transformation. Here’s the corrected version:

import dolfinx
import dolfinx.io
import ufl
from mpi4py import MPI
from slepc4py import SLEPc
import numpy as np

mesh = dolfinx.RectangleMesh(
    MPI.COMM_WORLD,
    [np.array([0, 0, 0]), np.array([np.pi, np.pi, 0])], [40, 40],
    cell_type=dolfinx.cpp.mesh.CellType.triangle)

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

N1curl = ufl.FiniteElement("Nedelec 1st kind H(curl)", mesh.ufl_cell(), 1)
V = dolfinx.FunctionSpace(mesh, N1curl)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

a = ufl.inner(ufl.curl(u), ufl.curl(v)) * ufl.dx
b = ufl.inner(u, v) * ufl.dx

bc_facets = np.where(np.array(dolfinx.cpp.mesh.compute_boundary_facets(mesh.topology)) == 1)[0]
bc_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, bc_facets)

u_bc = dolfinx.Function(V)
with u_bc.vector.localForm() as loc:
    loc.set(0)
bc = dolfinx.DirichletBC(u_bc, bc_dofs)

A = dolfinx.fem.assemble_matrix(a, bcs=[bc])
A.assemble()
B = dolfinx.fem.assemble_matrix(b, bcs=[bc])
B.assemble()

eps = SLEPc.EPS().create(MPI.COMM_WORLD)
eps.setOperators(A, B)
eps.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
st = eps.getST()
st.setType(SLEPc.ST.Type.SINVERT)
eps.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)
eps.setTarget(5.5)
eps.setDimensions(20)
eps.solve()

vals = [(i, eps.getEigenvalue(i)) for i in range(eps.getConverged()) if not
        (np.isclose(eps.getEigenvalue(i), 1) or np.isclose(eps.getEigenvalue(i), 0))]
vals.sort(key=lambda x: x[1].real)

j = 0
E = dolfinx.Function(V)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "out.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    for i, _ in vals:
        print(f"eigenvalue: {eps.getEigenpair(i, E.vector).real:.12f}")
        E.name = f"E-{j:03d}-{eps.getEigenpair(i, E.vector).real:.4f}"
        xdmf.write_function(E)
        j += 1
1 Like