Waveguide number of eigenmodes

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)

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

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

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.

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

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?`

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