No (Open) Source To Rely On – Finite Element Adventures, Chapter VI

Lately I have been trying to teach my kid about the importance of perserverence. It’s a good thing, because this ongoing project is the equivalent of standing in a traffic jam with your AC malfunctioning.

For various reasons, I have been trying for a while to get an eigenmode solver going. Apart for some application I want to place on my website, I think this it is the part of EM solvers I least understood, all these years. And boy, was I underestimating… In this series of posts, I will be covering only the most interesting parts of the ordeal.

It’s not over, yet. Not by a long-shot.

Step 1: Get yourself a toolbox

While starting off, I picked off using my alleged 2D finite element solver. I calculated all the matrix coefficients manually and hoped for the best. I had a lot of trouble defining the PEC boundary conditions, as well as a myriad of other things.

Soon afterwards, I noticed that I also can’t solve actual eigen-modes for the simplest of waveguides. The source code for it is in my github, should the avid reader should ever want to tackle this.

Step 2: Look Ma! A TE mode!

The next step was to look for an existing finite element library. The consideratios were pretty straight forward.

C++ or Python interface, Python is preferable.
Usable examples, specifically a waveguide eigenvalue problem. I will elaborate why this is important in a moment.
I will be able to build install it on my sparky-linux virtual machine. Sparky is light and fun, so why not.

A few libraries stood out at the time: XLife++, GetDP, FEniCS and Sparselizard. XLife++ doesn’t have a python interface and GetDP didn’t have a usable example. Sparselizard is relatively new, but it seemed the easiest to use, at the time. So let’s roll!

The first script looked something like this:

from spylizard import *
import numpy as np

import gmsh
import sys

##########
# Inputs #
##########

# Tags for various materials
skin = 1
air = 2
wholedomain = 2

# Frequency to search cutoff frequencies around
f_search = 1e9;

# Number of modes to calculate
numModesToCalc = 5

# Eigenmode to display
monToLoad = 0

mymesh = mesh("rect_WG_2D.msh")

#######################
# Initial definitions #
#######################

# Edge shape functions 'hcurl' for the electric field E.
# Fields x and y are the x and y coordinate fields.
E = field("hcurl") 
x = field("x") 
y = field("y")

# Use interpolation order 2 on the whole domain:
E.setorder(wholedomain, 2)

# Some coefficients
m0 = getmu0()
e0 = getepsilon0()
C = 1/np.sqrt(e0*m0)

k_search = 2*np.pi*f_search/C

# Electric permittivity:
er = parameter()
mr = parameter()

er.setvalue(air, 1)
mr.setvalue(air, 1)


# The waveguide is a perfect conductor. We thus force all
# tangential components of E to 0 on the waveguide skin.
E.setconstraint(skin)

##################################
# Generate weak-form formulation #
##################################

maxwell_curlE = formulation()
maxwell_EE = formulation()

# This is the weak formulation for electromagnetic waves:
maxwell_curlE += integral(wholedomain, (1/mr)*curl(dof(E))*curl(tf(E)))
maxwell_EE += integral(wholedomain, er*dof(E)*tf(E))

maxwell_curlE.generate()
maxwell_EE.generate()


# Get the stiffness and mass matrix:
A = maxwell_curlE.A()
B = maxwell_EE.A()

#############################
# Solve Eigenvalue problem  #
#############################

# Create the object to solve the generalized eigenvalue problem K*x = lambda*M*x :
eig = eigenvalue(A, B)

# Compute the 10 eigenvalues closest to the target magnitude 0.0 (i.e. the 10 first ones):
eig.compute(numModesToCalc, k_search**2.0)

lambda_r = np.array(eig.geteigenvaluerealpart())
lambda_i = np.array(eig.geteigenvalueimaginarypart())

# lam = np.array(lambda_r) + 1j*np.array(lambda_i)
lam = lambda_r
w_eig = np.sqrt(lam)*C
f_eig = w_eig/(2.0*np.pi) 


# Print the eigenfrequencies:
for fIdx in range(len(lam)):
    print("Eigenmode #{} has a cutoff frequency of {} GHz".format(fIdx + 1,f_eig[fIdx]/1e9))


# The eigenvectors are real thus we only need the real part:
Eeig_r = eig.geteigenvectorrealpart() 

fc = f_eig[monToLoad] 
kc = (2*np.pi*fc)/C

E.setdata(wholedomain, Eeig_r[monToLoad])
E.write(wholedomain, "Er.pos", 2)

Erx = compx(E)
Ery = compy(E)

# Initialize gmsh:
gmsh.initialize()

gmsh.open("Er.pos")
gmsh.merge("rect_WG_2D.geo")

# Creates graphical user interface
if 'close' not in sys.argv:
    gmsh.fltk.run()

# It finalize the Gmsh API
gmsh.finalize()

And the result was a very satisfying $TE_{10}$ mode,

Hoozah! I am the master of space and time! I can draw that mode I can probably produce faster analytically! Let’s delve into some details for a moment. If you are familiar with finite-element analysis, then you know each solved problem goes through something called a weak-form formulation. I am definitely not getting into what the that means right now, but it is sufficient, at this stage, to know that for the equation you want to solve, for example in this case

$\nabla\times \left({\frac{1}{\mu_r} \nabla\times\mathbf{E} }\right) - k_0^2\epsilon_r \mathbf{E} = -\mathrm{j}k_0 Z_0 \mathbf{J},$

exists a weak-form functional that has to be minimized. In the eigenvalue problem case, there are no sources. So it will be given by

$\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}.$

Does it always exist? Of course not. But if it doesn’t, there is no finite element solution. Moving on, the chosen library exapands the weak-form integral over the basis functions, for each finite-element. The eigen value problem is roughly phrased as

$\mathbf{A}\mathbf{E} = -k_0^2\mathbf{B}\mathbf{E}$

where the element in the matrices for the $e$ -th edge is given by

${{A}}^e = \iiint\limits_{e}{\frac{1}{\mu_r^e}\left[{\nabla\times\mathbf{N}^e}\right]\cdot\left[{\nabla\times\mathbf{N}^e}\right]^T\mathrm{d}V},$

$\mathbf{N}^e$ being the (vector) basis function for the e-th edge, and

${{B}}^e = \iiint\limits_{e}{\epsilon_r^e{\mathbf{N}^e}\cdot{\mathbf{N}^e}^T\mathrm{d}V}$ .

Notice that this formulation is generalized for 3D, but can easily be converted for 2D. In the software, all you need to do, is to load a 2D mesh… So here is a a picture of the mesh generated by gmsh:

The finite elements are in the shape of triangles, but can also have 4 edges. This script also allows plotting higher modes, of course.

This formulation, however, does not allow resolving TEM modes, nor QTEM modes. Why, you might ask? Notice that this formulation’s eigenvalues correspond to a quantity $-k_0^2$ , namely the wavenumber. Which wavenumber? The cutoff wavenumber, of course. As TEM modes have a cutoff frequency of 0, said expression is not valid for these modes.

Here are some more things that this formulation does not allow calculating explicitly:

Wave impedance
H-Field

Where these two are, of-course, tightly related, as the E\H field relation is the wave impedance.

What’s to be done, then?

Step 3: Needs More Harmonics

Luckily, wise people, not so long ago, have devised a specific formulation for eigenvalue problems. This formulation was devised as the “regular” continuous Lagrange basis functions sometimes yielded spurious solutions. There is tons of literature on this and of course, this issue is discussed thoroughly in “The Finite Element Method In Electromagnetics” by Jin. Very briefly, the curl of Lagrange basis functions is non-zero, which causes spurious solutions. The development of vector basis functions (sometimes called tangential basis functions, or Nedelec elements) removed that risk.

The weak-form formulation for this requires the definition of transverse and normal E-field components, $\mathbf{E}_t$ and $E_z$ ,

$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$

are the re-normalized transverse and normal E-field components, respectively. Notice that $\mathbf{e}_t$ is expanded over vector basis functions, while $e_z$ on the usual Lagrange basis functions.

Lucky me, there was an existing example in the Sparselizard python interface repository, called “photonic-waveguide-coupling-mode-2d“. The formulation differs a bit, but again, I’m powering through here, not going into the nitty-gritty details.

# Operators grad() and curl() in the transverse plane:
dtdtgradEl = expression(3,1,[dtdt(dx(dof(El))), dtdt(dy(dof(El))), 0])
gradtfEl = expression(3,1,[dx(tf(El)), dy(tf(El)), 0])

mode = formulation()

mode += integral(all, curl(dof(Et))*curl(tf(Et))- k0*k0*mur*epsr*(dof(Et))*tf(Et))
mode += integral(all, dtdtgradEl*tf(Et)+ dtdt(dof(Et))*tf(Et))

mode += integral(all, dtdtgradEl*gradtfEl+dtdt(dof(Et))*gradtfEl)
mode += integral(all, -k0*k0*mur*epsr*dtdt(dof(El))*tf(El))

mode.generate()

# Get the stiffness matrix K and mass matrix M:
K = mode.K()
M = mode.M()

# Create the object to solve the eigenvalue problem:
eig = eigenvalue(K, M)

Something is off, though. How come there are time derivatives here? Well, this is a very specific trait of Sparselizard. It was developed specifically to address harmonic-balance problems. Hence it allows direct formulation of time derivatives. A very nice webinar presented by Alexandre Halbach, the one who wrote Sparselizard, is available here.

This isn’t the only thing that differentiates it from other libraries, so keep up. The reason the time derivatives are defined here, is to instantly generate the two matrices needed for the eigenvalue problem. I’m going to change this, very soon, so I’m not diving too deep into this. Let’s see if we can manange solving the previous, rectangle wavguide mode, with this formulation:

All right, so far, so good. You may notice that the E-field scaling is off and that the polarization is inverse, but the mode shape is the same and that is enough for now. Let’s say I am intereseted in calculating the H-field. The analytical formulation should be something like

$-\frac{1}{\mathrm{j}\omega\mu} \nabla\times\mathbf{E}=\mathbf{H}.$

But something is missing. The values obtained for the E-field are real. No imaginary part. Let’s refer to the sparselizard manual, where it states that In-phase and Quadrature parts (and yes, I am using communiciation theory terms) of the mode are actually calulcated separately.

That means we need a different formulation in order to obtain these two harmonics. Also, for this formulation, a fundamental frequency for the harmonic balance calculation must be set.

So I attempted something like this:

setfundamentalfrequency(f0)

gradez = expression(3,1,[dx(dof(ez)), dy(dof(ez)), 0])
gradtfez = expression(3,1,[dx(tf(ez)), dy(tf(ez)), 0])

mode_1 = formulation()
mode_2 = formulation()

iC2 = 1/(C*C)

# Att
mode_1 += integral(wholedomain, (1/mr)*curl(dof(et))*curl(tf(et)) + iC2*er*dtdt(dof(et))*tf(et))

# Btt & Bzz
mode_2 += integral(wholedomain, (1/mr)*dof(et)*tf(et) + (1/mr)*gradez*gradtfez + iC2*er*dtdt(dof(ez))*tf(ez))
# Btz
mode_2 += integral(wholedomain, (1/mr)*dof(et)*gradtfez)
# Bzt
mode_2 += integral(wholedomain, (1/mr)*gradez*tf(et))

mode_1.generate()
mode_2.generate()

A = mode_1.A()
B = mode_2.A()

And got this error, in return:

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Arguments are incompatible
[0]PETSC ERROR: Dimensions of A and B do not match (2148, 2682)

Meaning it only formulated the $A_{tt}$ part of the matrix, without the zero padding. I did not find any specific information in the documentation regarding how to artificially pad the sparse matrices, so I started looking for a workaround at this point.

Interlude: Failing, Falling And Flailing

As I noted in the beginning, this is not a story of success. Every post will end in a failure. I did resolve most of the problems along the way, including the one above. but I ran into trouble in every single step. I do consider this to be a much more teaching experience and that is what I’m trying to create here.

I do take great interest in your inputs here! How would you resolve this? How would you approach these issues? What would you have done differently? I have put the codes for this chapter in my github, for you to test them out. I will post about how I solved this specific problem (before running into the next wall) in the next time.

The desired end result is a free-to-use online mode calculator for arbitrary waveguide structures. Let me know if you want to lend a hand!

RF With Care

No (Open) Source To Rely On – Finite Element Adventures, Chapter VI

Leave a Reply Cancel reply