Hello again, engineers of culture. Welcome back to the death-defying maze which is finite-element eigen-value problems.
In the last post, I was stuck with the fact that the matrix for the eigenvalue problem was actually comprised only of the
sub-matrix. I didn’t find in the documentation a way to artificially pad the matrix, so I looked for a workaround.
Initially I tried something like this:
mode_1 += integral(wholedomain, (1/mr)*curl(dof(et))*curl(tf(et)) + iC2*er*dtdt(dof(et))*tf(et))
mode_1 += integral(wholedomain, 0*((1/mr)*gradez*tf(et) + (1/mr)*dof(et)*tf(et) + iC2*er*dtdt(dof(ez))*tf(ez) + (1/mr)*gradez*gradtfez + (1/mr)*dof(et)*gradtfez))
mode_2 += integral(wholedomain, (1/mr)*(gradez*gradtfez + dof(et)*gradtfez))
mode_2 += integral(wholedomain, iC2*er*dtdt(dof(ez))*tf(ez))
mode_2 += integral(wholedomain, (1/mr)*(gradez*tf(et) + dof(et)*tf(et)))
mode_2 += integral(wholedomain, 0*((1/mr)*curl(dof(et))*curl(tf(et)) + iC2*er*dtdt(dof(et))*tf(et)))
This also didn’t work. It is a good time to remind ourselves that Sparselizard uses sparse matrices. Namely zeros do not enter the matrix.
After a lot of trial and error, I ended up with this formulation
mode_1 += integral(wholedomain, (1/mr)*curl(dof(et))*curl(tf(et)) + iC2*er*dtdt(dof(et))*tf(et))
mode_1 += integral(wholedomain, 1*((1/mr)*gradez*tf(et) + (1/mr)*dof(et)*tf(et) + iC2*er*dtdt(dof(ez))*tf(ez) + (1/mr)*gradez*gradtfez + (1/mr)*dof(et)*gradtfez))
mode_1 += integral(wholedomain, -1*((1/mr)*gradez*tf(et) + (1/mr)*dof(et)*tf(et) + iC2*er*dtdt(dof(ez))*tf(ez) + (1/mr)*gradez*gradtfez + (1/mr)*dof(et)*gradtfez))
mode_2 += integral(wholedomain, 1*((1/mr)*curl(dof(et))*curl(tf(et)) + iC2*er*dtdt(dof(et))*tf(et)))
mode_2 += integral(wholedomain, -1*((1/mr)*curl(dof(et))*curl(tf(et)) + iC2*er*dtdt(dof(et))*tf(et)))
mode_2 += integral(wholedomain, (1/mr)*(gradez*gradtfez + dof(et)*gradtfez))
mode_2 += integral(wholedomain, iC2*er*dtdt(dof(ez))*tf(ez))
mode_2 += integral(wholedomain, (1/mr)*(gradez*tf(et) + dof(et)*tf(et)))
A bit odd, and if Alex is reading this he is probably starting to jitter in anger right now, but it worked. Now the matrices are the same size and we are ready to solve eigenvalues again. I also noticed that the ordering of the adding and subtractions is critical. Both for the size of the matrices, and for the stability of the solution.
Here’s a riddle for you. Notice that the fundamental mode is actually calculated twice, this time.
Eigenmode #1 has a cutoff frequency of 1.3626946878782629 GHz
Eigenmode #2 has a cutoff frequency of 1.3626946879111876 GHz
Eigenmode #3 has a cutoff frequency of 3.536676520063204e-06 GHz
Eigenmode #4 has a cutoff frequency of 1.0316264790935775e-06 GHz
Eigenmode #5 has a cutoff frequency of 1.0316264790935775e-06 GHz
Why is that, you might ask? That is commonly referred to as mode degeneracy. For some reason, when allowed to calculated both real and imaginary parts of the mode (in-phase and quadrature), then it also calculates then as separate eigenvalues. This could have to do with the internal definitions of the eigenmode solver, but it doesn’t really matter. I’ll leave it to the reader to find a way to filter the degenerate modes.
In order to calculate the H-field, I calculated the curl operator manually. Curl operations are only allowed on field classes and not on array3x1.
The real and imaginary parts can be used to calculated the H-field, as such:
pol1 = (whichPol == 0)*2 + (whichPol == 1)*3
pol2 = (whichPol == 0)*3 + (whichPol == 1)*2
E_r = array3x1(compx(et.harmonic(pol1))/kz,compy(et.harmonic(pol1))/kz,el.harmonic(pol2))
Et_r = array3x1(compx(et.harmonic(pol1))/kz,compy(et.harmonic(pol1))/kz,0)
Ez_r = array3x1(0,0,ez.harmonic(pol2))
Ht_r = (1/(w0*m0*mr))*array3x1(
-(compy(et.harmonic(pol1)) + dy(ez.harmonic(pol1))),
(compx(et.harmonic(pol1)) + dx(ez.harmonic(pol1))),
0)
Hz_r = array3x1(0,0,compz(-(1/(w0*m0*mr*kz))*curl(et.harmonic(pol2))))
Notice how I renormalized the fields, and when I used harmonic #1 and #2. These are correspondent to the real and imaginary parts. Don’t worry, I thoroughly verified this script. Let’s plot the real and imaginary parts of the H-field:


Alrighty then, time to move forward
Step 4: A Feeble ATEMpt
So now I can finally get around to the main reason I started this process. T-Line impedance calculator! Let’s start with a simple Microstrip Line, shall we? First, a geometry and mesh are required.

And a fancy new script, also available in my github
from spylizard import *
import numpy as np
import gmsh
import sys
##########
# Inputs #
##########
# Tags for various materials
gnd = 1
line = 2
air = 3
fr4 = 4
bounds = 5
sub_epsr = 4.5
# Fundamental frequency for harmonics calculation
f0 = 2.0e9
# Number of modes to calculate
numModesToCalc = 5
# Eigenmode to display
monToLoad = 0
# Choose to display in-phase or quadrature
whichPol = 0
mymesh = mesh("MSL_XY.msh")
# Edge shape functions 'hcurl' for the electric field E.
# Fields x and y are the x and y coordinate fields.
x = field('x')
y = field('y')
et = field("hcurl",[2,3])
ez = field("h1",[2,3])
wholedomain = selectunion([air,fr4])
metals = selectunion([gnd,line,bounds])
# Use interpolation order 2 on the whole domain:
et.setorder(wholedomain, 2)
ez.setorder(wholedomain, 2)
# The cutoff frequency for a 0.2 m width is freq = 0.75 GHz in theory.
# With this code and a fine enough mesh you will get the same value.
m0 = getmu0()
e0 = getepsilon0()
C = 1/np.sqrt(e0*m0)
k0 = 2.0*getpi()*f0/C
w0 = 2*getpi()*f0
# Electric permittivity:
er = parameter()
mr = parameter()
# Set Dirichlet Et = 0 boundary conditions on metalic shell
er.setvalue(air, 1.0)
mr.setvalue(air, 1.0)
er.setvalue(fr4, sub_epsr)
mr.setvalue(fr4, 1.0)
# The waveguide is a perfect conductor. We thus force all
# tangential components of E to 0 on the waveguide skin.
et.setconstraint(metals)
ez.setconstraint(metals)
#########################
# Weak form formulation #
#########################
# Set the fundamental frequency for the multi-harmonic calculation
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)
mode_1 += integral(wholedomain, (1/mr)*curl(dof(et))*curl(tf(et)) + iC2*er*dtdt(dof(et))*tf(et))
mode_1 += integral(wholedomain, 1*((1/mr)*gradez*tf(et) + (1/mr)*dof(et)*tf(et) + iC2*er*dtdt(dof(ez))*tf(ez) + (1/mr)*gradez*gradtfez + (1/mr)*dof(et)*gradtfez))
mode_1 += integral(wholedomain, -1*((1/mr)*gradez*tf(et) + (1/mr)*dof(et)*tf(et) + iC2*er*dtdt(dof(ez))*tf(ez) + (1/mr)*gradez*gradtfez + (1/mr)*dof(et)*gradtfez))
mode_2 += integral(wholedomain, 1*((1/mr)*curl(dof(et))*curl(tf(et)) + iC2*er*dtdt(dof(et))*tf(et)))
mode_2 += integral(wholedomain, -1*((1/mr)*curl(dof(et))*curl(tf(et)) + iC2*er*dtdt(dof(et))*tf(et)))
mode_2 += integral(wholedomain, (1/mr)*(gradez*gradtfez + dof(et)*gradtfez))
mode_2 += integral(wholedomain, iC2*er*dtdt(dof(ez))*tf(ez))
mode_2 += integral(wholedomain, (1/mr)*(gradez*tf(et) + dof(et)*tf(et)))
mode_1.generate()
mode_2.generate()
A = mode_1.A()
B = mode_2.A()
# Create the object to solve the generalized eigenvalue problem K*x = lambda*M*x :
eig = eigenvalue(A, B)
kz2_search = (k0**2.0)*0.5*(sub_epsr + 1)
# Compute the 10 eigenvalues closest to the target magnitude 0.0 (i.e. the 10 first ones):
eig.compute(numModesToCalc, -kz2_search)
kz2_found = -np.array(eig.geteigenvaluerealpart())
kc_found = np.sqrt(k0**2.0 - kz2_found)
fc_found = kc_found*C/(2.0*np.pi)
# Print the eigenfrequencies:
for fIdx in range(len(fc_found)):
print("Eigenmode #{} has a cutoff frequency of {} GHz".format(fIdx + 1,fc_found[fIdx]/1e9))
kz = np.sqrt(kz2_found[monToLoad])
########################
# Calculate eigenmodes #
########################
Evec = eig.geteigenvectorrealpart()
et.setdata(wholedomain,Evec[monToLoad])
ez.setdata(wholedomain,Evec[monToLoad])
print("Drawing mode with fc = {} GHz\n".format(fc_found[monToLoad]/1e9))
pol1 = (whichPol == 0)*2 + (whichPol == 1)*3
pol2 = (whichPol == 0)*3 + (whichPol == 1)*2
# One way to phrase all three components of the E-field
E_r = array3x1(compx(et.harmonic(pol1))/kz,compy(et.harmonic(pol1))/kz,ez.harmonic(pol2))
# Can also disect it to it's transverse and normal components.
Et_r = array3x1(compx(et.harmonic(pol1))/kz,compy(et.harmonic(pol1))/kz,0)
Ez_r = array3x1(0,0,ez.harmonic(pol2))
# And now calculate the H-field, finally
Ht_r = (1/(w0*m0*mr))*array3x1(
-(compy(et.harmonic(pol1)) + dy(ez.harmonic(pol1))),
(compx(et.harmonic(pol1)) + dx(ez.harmonic(pol1))),
0)
Hz_r = array3x1(0,0,compz(-(1/(w0*m0*mr*kz))*curl(et.harmonic(pol2))))
H_r = array3x1(compx(Ht_r),compy(Ht_r),compz(Hz_r))
E_r.write(wholedomain, "Er.pos", 2)
H_r.write(wholedomain, "Hr.pos", 2)
# Initialize gmsh:
gmsh.initialize()
gmsh.open("Er.pos")
gmsh.merge("MSL_XY.msh")
# Creates graphical user interface
if 'close' not in sys.argv:
gmsh.fltk.run()
# It finalize the Gmsh API
gmsh.finalize()
Even while running, something isn’t right.
Eigenmode #1 has a cutoff frequency of nan GHz
Eigenmode #2 has a cutoff frequency of nan GHz/home/gadi/eclipse-workspace/sparselizard_python_tests/WG_Eigenmode_Multiharmonic_MSL.py:115: RuntimeWarning: invalid value encountered in sqrt
kc_found = np.sqrt(k0**2.0 - kz2_found)
Eigenmode #3 has a cutoff frequency of nan GHz
Eigenmode #4 has a cutoff frequency of nan GHz
Eigenmode #5 has a cutoff frequency of nan GHz
The NaNs (not-a-number) are happening as the script is trying to calculate the square root of a negative number. Well, this could be because of a numerical error, due to the fact that the cutoff frequency is very close to zero. I can live with that! Let’s plot the mode then!

Where is the mode. Where is it!! What is this singular piece of dung in the middle of my structure!
This is the part that I am sparing you how much time and effort it took me to arrive to certain conclusions. I have tried looking for higher modes, lower modes, change the search values for the eigenmode solver, everything. I did understand, however that this is a numerical problem.
Interlude: Eigenbrick Wall
We meet again, in the realm of the unknown. Have you encountered this problem? How would you approach this? Let me know!
I can share with you that I solved this problem 60% by accident. I did take a logical step, but the result wasn’t what I exactly intended it to be.
See ya’ll next week!