Waveguide number of eigenmodes

dolfinx
1

I’m trying to implement this Dolfin waveguide example in DolfinX.
The implementation below calculates only one eigenpair. If I set solver to SLEPc.EPS.Type.LAPACK, it calculats about a 2500 eigenpairs.

What settings should I use to get a number of eigenpairs requested by eps.setDimensions?

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

c0 = 299792458

mesh = dolfinx.RectangleMesh(
    MPI.COMM_WORLD,
    [np.array([0, 0, 0]), np.array([1, 0.5, 0])], [30, 30],
    cell_type=dolfinx.cpp.mesh.CellType.triangle)
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

N1curl = ufl.FiniteElement("N1curl", mesh.ufl_cell(), 2)
V = dolfinx.FunctionSpace(mesh, N1curl)

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


a = ufl.inner(ufl.curl(v), ufl.curl(u)) * ufl.dx
b = ufl.inner(v, u) * 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)
eps.setProblemType(SLEPc.EPS.ProblemType.GHEP)
eps.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_REAL)
# eps.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_REAL)
# eps.setTarget(10)
eps.setDimensions(42)
eps.solve()

solution = dolfinx.Function(V)

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "out.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    for i in range(eps.getConverged()):
        print(eps.getEigenpair(i, solution.vector))
        solution.name = f"f-{i}"
        xdmf.write_function(solution)
2

And when I substitute \nabla\times\nabla\times\vec{E}=\nabla\underbrace{\left(\nabla\cdot\vec{E}\right)}_{0}-\nabla^2\vec{E}, I get 42 funny looking eigenvectors.
image

3

The analysis and why you’re seeing “funny” eigenvectors is well known, e.g. this book. See also this demo.

4

But I do use Nedelec elements. My program setup is similar to the one, shown in the demo except for the target magnitude option. When I use target_magnitude or target_real, I only get zero vectors.

5

Ok, problem turned out to be a selection of a target value. It should be somewhat above the expected eigenvalue. SLEPc configurations are

eps = SLEPc.EPS().create(MPI.COMM_WORLD)
eps.setOperators(A, B)
eps.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
eps.setProblemType(SLEPc.EPS.ProblemType.GHEP)
eps.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)
eps.setPowerShiftType(SLEPc.EPS.PowerShiftType.CONSTANT)
eps.setTarget(100)
eps.setDimensions(10)
eps.solve()
1 Like
6

Now formulation \int_\Omega{\nabla\times v \nabla\times E\;dx}=\int_\Omega{v E\;dx} produces the same results as the demo. But when applying the aforementioned substitution \nabla\times\nabla\times E = -\nabla^2 E I only get negative eigenvalues and those “funny” eigenvectors (if eps.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_REAL)) or no eigenpairs at all (if eps.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_MAGNITUDE)). Why formulation -\int_\Omega{\nabla v \nabla E\;dx}=\int_\Omega{v E\;dx} cannot be used to solve this problem?`

7

There are many insidious issues with the discretisation of the Maxwell operator. Like I said above, I recommend reading Peter Monk’s book or perhaps Jin’s book.

1 Like
8

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