Electromagnetic Eigenvalue problem in FEM yielding spurious solutions

Question

I have written an Electromagnetic FEM solver for waveguide Eigenvalue problems. I already wrote it in two different libraries, sparselizard and dolfinx. The relevant scripts are available here for sparselizard and here for dolfinx.

While solving the EM wave equations in two dimensions, generally this is the formulation

$$\nabla\times\left({\frac{1}{\mu_r}\nabla\times\mathbf{E}}\right)-k_0^2\epsilon_r\mathbf{E}=0$$

is used. It's weak form formulation is

$$4 \frac{1}{2}\iiint{\left[{\frac{1}{\mu_r}\left({\nabla\times\mathbf{E}}\right)\cdot\frac{1}{\mu_r}\left({\nabla\times\mathbf{E}}\right)-k_0^2\epsilon_r\mathbf{E}\cdot\mathbf{E}}\right]\mathrm{d}V}.$$

However, when the operating frequency, $k_0$ is zero, this equation becomes singular. This is very relevant for certain types of waveguides, for example, microstrips and coaxial lines. From here on forward I will refer to the direction of propagation as $\hat{z}$ and the components on the breadth of the waveguide as transversal components.

For this, it is common to separate $\mathbf{E}$ into two separate components, $\mathbf{E}_t$ for the transversal parts of the field, and $E_z$ for the part. Substituting this to the prior weak formulation obtains:

$$ F\left({\mathbf{E}}\right) = \frac{1}{2} \iint\limits_{S} \left[{ \frac{1}{\mu_r}\left({\nabla_t\times \mathbf{e}_t }\right) \cdot\left({\nabla_t\times{\mathbf{e}_t}}\right)^* - k_0^2\epsilon_r\mathbf{e}_t\cdot\mathbf{e}_t^* }\right]\mathrm{d}S + k_z^2\iint\limits{S} {\left[{ \frac{1}{\mu_r}\left({\nabla_t e_z + \mathbf{e}_t}\right) \cdot\left({\nabla_t e_z + \mathbf{e}_t}\right)^* - k_0^2\epsilon_r e_z\cdot e_z^* }\right]}\mathrm{d}S.$$

where $\mathbf{e}_t = k_z\mathbf{E}_t$ and $e_z = \mathrm{j}E_z$.

The basis functions used for $\mathbf{E}_t$ are Nédélec elements, whereas for $E_z$ simple Nodal basis functions.

I have been successful in solving the following types of waveguides (all found in my GitHub):

• Circular\rectangular\corrugated waveguide.
• Microstrip\GCPW\CPW.
• Differential lines.

Which is more than I need, most of the time, except that for some reason I cannot solve the coaxial waveguide problem. It is clear to me that the solver creates a spurious solution.

and I am pretty sure at this point that this is due to the fact that this structure does not support E-field in the $\hat{z}$ direction. An important addition is that the finite element matrix is composed of 5 sub-matrices:

$$\mathbf{B} = \left[{\matrix{\mathbf{B}_{tt} & \mathbf{B}_{tz} \\ \mathbf{B}_{zt} & \mathbf{B}_{zz}}}\right]$$

and

$$\mathbf{A} = \left[{\matrix{\mathbf{A}_{tt} & \mathbf{0} \\ \mathbf{0} & \mathbf{0}}}\right].$$

I think that this set of constraints creates a singular matrix, thus not allowing a solution.

In sparselizard, I am getting the above drawing. However, in dolfinx I am getting a solver error.

Edit #1

As @rchilton1980 suggested, I reviewed the articles about the solutions for the spurios space, generated by modes applying $Ez = 0$. The method suggested there has to do with the choice of modifying the matrices eigenvalues such that

$\tilde{A} = B$, $\tilde{B} = \mathbf{B} + \frac{\mathbf{A}}{\Theta^2}$

and $\lambda = \frac{\Theta^2}{\Theta^2 - k_z^2}$.

However, it is also suggested to modify the initial condition such that it is is $\mathbf{A}$ orthogonal and $\mathbf{B}$ orthogonal to the unwanted subspace.

@rchilton1980 suggests strongly to "deflate against" this subspace (quotes are because I don't exactly know what that means, but I have a good guess).

Fact of the matter is, both Sparselizard and DolfinX use SLEPc is their eigevalue problem solver. I know the setup for the solver in Sparselizard (still my current choice for the solver).

EPS eps;
EPSCreate( PETSC_COMM_SELF, &eps );
EPSSetOperators( eps, myA.getapetsc(), myB.getapetsc() );
EPSSetProblemType(eps, EPS_GNHEP);
EPSSetDimensions(eps, numeigenvaluestocompute, PETSC_DECIDE, PETSC_DECIDE);
EPSSetTolerances(eps, 1e-6, 100);
EPSSetType(eps, EPSKRYLOVSCHUR);
EPSSetFromOptions(eps);

ST st;
EPSGetST(eps, &st);
STSetType(st, STSINVERT);

EPSSetTarget(eps, targeteigenvaluemagnitude);
EPSSetWhichEigenpairs(eps, EPS_TARGET_MAGNITUDE);

KSP ksp;
STGetKSP(st, &ksp);
KSPSetType(ksp, "preonly");
PC pc;
KSPGetPC(ksp, &pc);
PCSetType(pc, PCLU);
PCFactorSetMatSolverType(pc, universe::solvertype);

// DO THE ACTUAL RESOLUTION:
EPSSolve( eps );

Two problems:

1. I can't manually place an initial solution set, unless I make significant changes to the library itself.
2. I don't currently know how to generate the $\mathbf{B}^{-1}\mathbf{A}$ using the existing python code. Even if I do, I'm worried about how the Dirichlet boundary conditions are precisely implemented there.

I will look into changing the solver method suggested, for now.

---

Last edit

---

Sometimes the solution is simple (and stupid). There was an issue with my generated mesh file.

The normals were not defined on most of the metals...

Sorry for the annoyance. This ones on me :)

Answer 1

You are correct, this formulation does introduce a spurious space. It's not actually due to any function space/discretization choice, but rather that change of variables step ($\mathbf e_t = k_z \mathbf E_t$), which introduces a large number of nonphysical eigenpairs that possess $k_z=0$.

For remedy, I would direct you to "Full-Wave Analysis of Dielectric Waveguides using Tangential Vector Finite Elements", by Jin-Fa Lee, Din-Kow Sun and Zoltan Cendes. You have arrived at their equation (25), but if you continue to equation (26) they describe the structure of this spurious space (and beyond that, how to suppress it during power iteration or Arnoldi iteration). In fact, give the whole paper a thorough read, because it also describes a shifting technique that makes the scheme converge extremely rapidly/robustly to propagating modes (which are typically of greatest interest). In my experience, this scheme is rather fragile without that shifting technique.