Improvements on RF waveguide design tutorial

[duarte-jfs]

@HelgeGehring Follow up from https://github.com/HelgeGehring/femwell/pull/159#issuecomment-2136159365 .

There's some code multiple times in the example, for example the mesh generation. We could either generate all meshes in the beginning or we could provide a function which returns a mesh and call it with different parameters. (Less code, easier to read, easier to maintain) there's I think also some other code more than once

Indeed, it can be made more compact. I will work on it and organize the code in a class CPW. I used that approach in another project and it made things more compact and readable.

kb and T are not used at all, right? Should we remove them?

I need to double check. There's a good chance they're not. The use of pint came as a follow up from another project and I ended up pasting things at the top hoping I'd use them and be safe.

For the other constants, we could just use the constants defined in scipy, we don't need to have the numeric values, right?

We don't need the numeric and we can use the ones from scipy. but if pint is actively used in the notebook it's best to redefine the constants as pint.Quantity to take care of unit tracking. It needs to be double checked.

We can probably generalize the calculations of I0,V0,R,L,C,G and put them in the mode class, right?

This one I would not do. At least without deeper study of multi conductor transmission line (MTL) theory and application in multimode waveguides (coupled waveguides). When I started I picked the tutorial on the repo for the CPW and the theory you follow is from MTL. However, when reading up on it, I found the math to quickly become more involved, with mathematical artifacts coming up when analyzing a waveguide as a circuit. A big issue is on the normalization of the fields, which depending on which formalism you choose, you may either conserve power and not reciprocity, not conserve power and have reciprocity or a mix of both,. The point being: I feel like it all comes down the your choice when analysing, that's why I put some effort illustrating the choices I've made to calculate i0 and v0, since those are values that made intuitive sense, and I think we should be aware of those assumptions. As for the calculation for RLGC parameters, I believe those are obtained with the assumption of single mode operation. I don't know how to scale it to multimode scenarios. One of the papers on the tutorial has a complete theory on MTL, which could help. Bottom line: I would be weary of generalizing this analysis due to the ambiguity present.

I think the images from the paper are missing :/ maybe we should also just make our own sketches there?

That one is my fault. When removing whitespaces I forgot to correct the filenames on the notebook's markdown cells.

I didn't know where you used [3] and [4], so i just put them somewhere in the beginning

Are these numbers correct? Because [4] I know it is used and [3] is to link to delft's software, which now that you mention it I'm not sure it's the correct reference or if it is linked to somewhere in the text. But if not, it should definitely be included. It's a nice software that they have.

Also, one thing I noticed @HelgeGehring : online, the tutorial has:

PI model : |Z0| = 70.88 Ohm | angle(Z0) = -0.01 pi radians PV model : |Z0| = 60.19 Ohm | angle(Z0) = -0.01 pi radians VI model : |Z0| = 65.31 Ohm | angle(Z0) = -0.00 pi radians RLGC : |Z0| = 64.97 Ohm | angle(Z0) = -0.01 pi radians

While in my laptop I get:

PI model : |Z0| = 61.77 Ohm | angle(Z0) = -0.01 pi radians PV model : |Z0| = 60.19 Ohm | angle(Z0) = -0.01 pi radians VI model : |Z0| = 60.97 Ohm | angle(Z0) = -0.00 pi radians RLGC : |Z0| = 60.66 Ohm | angle(Z0) = -0.01 pi radians

Unclear if it's something to be concerned about...

[simbilod]

The numerical mismatch is indeed a worry. There must be something different between the two versions?

[duarte-jfs]

The numerical mismatch is indeed a worry. There must be something different between the two versions?

The difference comes from the dependencies versions. I have created a new environment with a fresh install of femwell and I also get something similar to the above (not exactly the same though). I suspect that if the packages are installed the mismatch will go away. The difference in versions are below (package_name femwell_repo_version local_venv_version):

annotated-types 0.7.0 ansicon 1.89.0 blessed 1.20.0 colorama 0.4.6 contourpy 1.2.1 1.0.7 cycler 0.12.1 0.11.0 enlighten 1.12.4 femwell 0.1.11 fonttools 4.53.0 4.39.3 gmsh 4.13.1 4.12.2 jinxed 1.2.1 1.2.0 kiwisolver 1.4.5 1.4.4 markdown-it-py 3.0.0
matplotlib 3.9.0
mdurl 0.1.2 0.1.2 meshio 5.3.5
meshwell 1.0.7 1.0.7 numpy 1.26.4 1.24.2 packaging 24.0 23.0 pillow 10.3.0 9.5.0 Pint 0.23
pip 22.3.1
prefixed 0.7.1 pydantic 2.7.3 2.6.4 pydantic_core 2.18.4 2.16.3 Pygments 2.18.0 2.17.2 pygmsh 7.1.17 7.1.17 pyparsing 3.1.2 3.0.9 python-dateutil 2.9.0.post0
rich 13.7.1
scikit-fem 9.1.1 9.0.1 scipy 1.13.1 1.10.1 setuptools 65.5.0
shapely 2.0.4 2.0.3 six 1.16.0
tqdm 4.66.4 typing_extensions 4.12.1 wcwidth 0.2.13

Do any of these seem like a red flag?

The results are now:

PI model : |Z0| = 76.77 Ohm | angle(Z0) = -0.01 pi radians PV model : |Z0| = 60.19 Ohm | angle(Z0) = -0.01 pi radians VI model : |Z0| = 67.97 Ohm | angle(Z0) = -0.00 pi radians RLGC : |Z0| = 67.62 Ohm | angle(Z0) = -0.01 pi radians

However, the results when considering the symmetry plane are great:

PI model : |Z0| = 161.84 Ohm | angle(Z0) = -0.01 pi radians PV model : |Z0| = 161.84 Ohm | angle(Z0) = -0.01 pi radians VI model : |Z0| = 161.84 Ohm | angle(Z0) = -0.01 pi radians RLGC : |Z0| = 159.84 Ohm | angle(Z0) = -0.01 pi radians

Maybe it's a meshing issue?

[HelgeGehring]

Indeed, it can be made more compact. I will work on it and organize the code in a class CPW. I used that approach in another project and it made things more compact and readable.

Sounds great!

kb and T are not used at all, right? Should we remove them?

I need to double check. There's a good chance they're not. The use of pint came as a follow up from another project and I ended up pasting things at the top hoping I'd use them and be safe.

Sounds good!

For the other constants, we could just use the constants defined in scipy, we don't need to have the numeric values, right?

We don't need the numeric and we can use the ones from scipy. but if pint is actively used in the notebook it's best to redefine the constants as pint.Quantity to take care of unit tracking. It needs to be double checked.

I'm totally fine with defining them with quantities, that's a great approach to avoid problems! I just guess if we use the scipy values and define the constants based on them, it simplifies the code further :)

We can probably generalize the calculations of I0,V0,R,L,C,G and put them in the mode class, right?

This one I would not do. At least without deeper study of multi conductor transmission line (MTL) theory and application in multimode waveguides (coupled waveguides). When I started I picked the tutorial on the repo for the CPW and the theory you follow is from MTL. However, when reading up on it, I found the math to quickly become more involved, with mathematical artifacts coming up when analyzing a waveguide as a circuit. A big issue is on the normalization of the fields, which depending on which formalism you choose, you may either conserve power and not reciprocity, not conserve power and have reciprocity or a mix of both,. The point being: I feel like it all comes down the your choice when analysing, that's why I put some effort illustrating the choices I've made to calculate i0 and v0, since those are values that made intuitive sense, and I think we should be aware of those assumptions. As for the calculation for RLGC parameters, I believe those are obtained with the assumption of single mode operation. I don't know how to scale it to multimode scenarios. One of the papers on the tutorial has a complete theory on MTL, which could help. Bottom line: I would be weary of generalizing this analysis due to the ambiguity present.

Yay! Generalizing this and understanding multi conductor transmission lines sounds amazing! 🥳

I think the images from the paper are missing :/ maybe we should also just make our own sketches there?

That one is my fault. When removing whitespaces I forgot to correct the filenames on the notebook's markdown cells.

Sounds easy to fix :)

I didn't know where you used [3] and [4], so i just put them somewhere in the beginning

Are these numbers correct? Because [4] I know it is used and [3] is to link to delft's software, which now that you mention it I'm not sure it's the correct reference or if it is linked to somewhere in the text. But if not, it should definitely be included. It's a nice software that they have.

Hmm, feel free to just adjust it that the numbers are right, I went quick quickly over it, but I can give it another look. The {cite} command makes it way easier to refer to certain papers.

Also, one thing I noticed @HelgeGehring : online, the tutorial has:

PI model : |Z0| = 70.88 Ohm | angle(Z0) = -0.01 pi radians PV model : |Z0| = 60.19 Ohm | angle(Z0) = -0.01 pi radians VI model : |Z0| = 65.31 Ohm | angle(Z0) = -0.00 pi radians RLGC : |Z0| = 64.97 Ohm | angle(Z0) = -0.01 pi radians

While in my laptop I get:

PI model : |Z0| = 61.77 Ohm | angle(Z0) = -0.01 pi radians PV model : |Z0| = 60.19 Ohm | angle(Z0) = -0.01 pi radians VI model : |Z0| = 60.97 Ohm | angle(Z0) = -0.00 pi radians RLGC : |Z0| = 60.66 Ohm | angle(Z0) = -0.01 pi radians

Unclear if it's something to be concerned about...

We should for sure be worried!

My first guess would be that the integrals depend somehow on the phase of the mode.

As the modes are normalized to have a power of 1, the phase of the modes can still be arbitrary. While it often shows the same modes on the same machine, these different absolute phases can lead to different numbers in calculations. Usually, physical quantities which we calculate from such modes should not depend on the phase. Might that explain the differences?

[HelgeGehring]

Maybe try first the older gmsh version

gmsh 4.13.1 4.12.2

We had to pin it for now as the newest version created some trouble. If that doesn't fix it we can look into the other dependencies.

[duarte-jfs]

Maybe try first the older gmsh version

gmsh 4.13.1 4.12.2

We had to pin it for now as the newest version created some trouble. If that doesn't fix it we can look into the other dependencies.

I tried downgrading gmsh but the results are the same

As the modes are normalized to have a power of 1, the phase of the modes can still be arbitrary. While it often shows the same modes on the same machine, these different absolute phases can lead to different numbers in calculations. Usually, physical quantities which we calculate from such modes should not depend on the phase. Might that explain the differences?

A change in phase should not alter the result. It's the absolute value of that has changed. For the PI model to have changed the value so much, it means the current i0 calculated is now smaller than before, which is just the integral around the signal track. I assume it's either the mesh is now different or there is something in the solver

EDIT:

@HelgeGehring I have made the following successive downgrades/upgrades:

pip install gmsh==4.12.2 -> mismatch remained pip install numpy==1.24.2 -> mismatch remained pip install scikit-fem==9.0.1 -> Results match pip install numpy==1.26.4 -> Results match pip install gmsh==4.13.1 -> Results match

Something seems to be happening with scikit-fem latest version.

[duarte-jfs]

I think the images from the paper are missing :/ maybe we should also just make our own sketches there?

That one is my fault. When removing whitespaces I forgot to correct the filenames on the notebook's markdown cells.

Sounds easy to fix :)

Turns out the .py file has the correct filenames. When I launch it with jupyter lab all the images appear. Could it be something with the online website for image display? I use relative paths like this support\tuncer_results.png

[HelgeGehring]

I think hte main problem is that the pictures are just not there :D https://github.com/HelgeGehring/femwell/tree/main/docs/electronics/examples/RF_CPW_transmission_line_tutorial/support

Nice catch about the versions! Just as a sanity check: If you use all the newest versions except scikit-fem, do you also get the mismatch? Is it only the latest version? I.e. 9.0.0 gives different results compared to 9.0.1? Or is 9.0.0 also having the problem?

[duarte-jfs]

I think hte main problem is that the pictures are just not there :D

ah... that explains it XD

Nice catch about the versions! Just as a sanity check: If you use all the newest versions except scikit-fem, do you also get the mismatch? Is it only the latest version? I.e. 9.0.0 gives different results compared to 9.0.1? Or is 9.0.0 also having the problem?

Yes, the tests were with all the versions up to date (as in the list above, left column) and then I did the downgrades/upgrades. I just tried 9.0.0, but it throws an error with meshio:

Traceback (most recent call last): File "C:\Users\20230622\Downloads\RF_CPW_transmission_line_tutorial.py", line 33, in <module> from skfem.io.meshio import from_meshio File "C:\Users\20230622\AppData\Local\anaconda3\envs\phd_1\lib\site-packages\skfem\io\__init__.py", line 1, in <module> from .meshio import from_meshio # noqa File "C:\Users\20230622\AppData\Local\anaconda3\envs\phd_1\lib\site-packages\skfem\io\meshio.py", line 5, in <module> import meshio ModuleNotFoundError: No module named 'meshio'

[simbilod]

Strange, meshio is still a dependency of scikit-fem https://github.com/kinnala/scikit-fem/blob/master/requirements.txt

[HelgeGehring]

you can just install meshio additionally for now :) we just need to know where exactly it changed - from there we could also go on with git bisect to find the commit which changed it And then understand if before or after is correct.

Do you have a reference simulation which we could compare to see if the new version is correct or the old one?

[duarte-jfs]

you can just install meshio additionally for now :) we just need to know where exactly it changed - from there we could also go on with git bisect to find the commit which changed it And then understand if before or after is correct.

I just tried this, and the version 9.0.0 gives the same results as the 9.0.1

Do you have a reference simulation which we could compare to see if the new version is correct or the old one?

Not perfect ones... For the current example, the three approaches should be close together since the modes are very close to TEM modes. I have tried another example which has an analytical solution and the same numerical difference appears. The example is a shielded microstrip line surrounded by air. We can use a conformal mapping formula to get the characteristic impedance which is valid for infinitely thin metal strip and perfect conductors.

You can find the .py file for it here: microstrip.zip

With this example, scikit-fem==9.1.1 gives me:

Number of mesh elements: 20925
Modes neff: [1.00485251]
modes loss: [-0.44647038] dB/cm
Modes transversality: 0.9999998925613547
PI model : |Z0| = 124.81 Ohm  |  angle(Z0) = -0.00 pi radians
PV model : |Z0| = 115.93 Ohm  |  angle(Z0) = -0.00 pi radians
VI model : |Z0| = 120.28 Ohm  |  angle(Z0) = 0.00 pi radians
Conformal mapping: 112.84726343231252 ohm

And when downgrading to scikit-fem==9.0.1:

Number of mesh elements: 20925
Modes neff: [1.00485251]
modes loss: [-0.44647038] dB/cm
Modes transversality: 0.9999998925613547
PI model : |Z0| = 115.79 Ohm  |  angle(Z0) = -0.00 pi radians
PV model : |Z0| = 115.93 Ohm  |  angle(Z0) = -0.00 pi radians
VI model : |Z0| = 115.86 Ohm  |  angle(Z0) = -0.00 pi radians
Conformal mapping: 112.84726343231252 ohm

I also had to install and uninstall meshio for some reason. This was the cmd prompt:

pip install scikit-fem==9.1.1 & pip uninstall meshio & pip install meshio & python microstrip.py & pip install scikit-fem==9.0.1 & pip uninstall meshio & pip install meshio & python microstrip.py

I'm not sure if one version is more correct over the other, but there certainly is a difference.

You can also check the theory in this book (chapter 3.14): Collin, Robert E. 2001. Foundations for Microwave Engineering. 2. ed., [Nachdr.]. IEEE Press Series on Electromagnetic Wave Theory. New York, NY: IEEE Press.

[HelgeGehring]

ok, thanks, were narrowing it down!

i think the update of meshio was due to an update of numpy which required some changes. I/d guess you use the same numpy version for both?

It's somehow interesting that before the update the values were almost the same while afterwards they got different. Also, it seems that the PV model is not affected at all? As only i0 depends on the H-field (p should be nomalized to 1? Or is there something different because of the definition you use?), so maybe there's a change in the calculation of the H-field?

Is there a meshio version which works with both skfem versions?

It also seems that there was a change in the interface to meshio (which would make me doubt that there's a version working with both)... , see https://github.com/kinnala/scikit-fem/compare/9.0.1...9.1.1

I'd guess best would be if we could make a simplified short version of the example (i.e. only a mode solve on a simple geometry and then the integral which is different between the two), that way we can easier iterate and test. (also, the voltage / current integral is the same? The difference is just the sign in front and which field we feed in? I.e. we could there also simplify the example?)

We should then compare the exact values of the E / H fields, I'd suspect that the result of the calculation of the H-field might have changed. (easiest would be to just dump the array via numpy.save and load it with the newer version and compare) If those are the same, it must be the calculation of the line integrals.

[HelgeGehring]

The other thing I'm wondering: Is one of the line integrals you calculate along an interface at which the epsilon changes discontinuously? (I.e. the surface of the conductor?) In that case, the integral would depend on the direction in which the line integral is summed up (or more correct, on which DOFs the line integral is based, the ones with high epsilon or low epsilon) As we use the nedelec element for the in-plane components, the normal component within the plane is not continuous accross interfaces. Thus, if only part of the path is along the interface where such a change happens I'd guess there would be something wrong in the integral in the end... We should look deeper into how exactly this integral is calculated and if it's dependend on the DOFs outside or inside of the loop.

[duarte-jfs]

ok, thanks, were narrowing it down!

i think the update of meshio was due to an update of numpy which required some changes. I/d guess you use the same numpy version for both?

Yes, I have checked and both versions use numpy==1.26.4.

Is there a meshio version which works with both skfem versions?

It also seems that there was a change in the interface to meshio (which would make me doubt that there's a version working with both)... , see kinnala/scikit-fem@9.0.1...9.1.1

I have also checked that, and it successfully install meshio==5.3.5 for both versions.

I'd guess best would be if we could make a simplified short version of the example (i.e. only a mode solve on a simple geometry and then the integral which is different between the two), that way we can easier iterate and test. (also, the voltage / current integral is the same? The difference is just the sign in front and which field we feed in? I.e. we could there also simplify the example?)

I think the current example of the microstrip is good because it is the simplest geometry: just a square on a simulation region.

It's somehow interesting that before the update the values were almost the same while afterwards they got different. Also, it seems that the PV model is not affected at all? As only i0 depends on the H-field (p should be nomalized to 1? Or is there something different because of the definition you use?), so maybe there's a change in the calculation of the H-field?

We should then compare the exact values of the E / H fields, I'd suspect that the result of the calculation of the H-field might have changed. (easiest would be to just dump the array via numpy.save and load it with the newer version and compare) If those are the same, it must be the calculation of the line integrals.

Indeed! I have narrowed it down as you say to the integral calculation. The first thing I checked was to see if the modes calculated are the same: they are. Next I check the integral calculation and there is where the difference is. Here's what I did: I have the following snippet of code placed after the geometry definition, mesh generation, epsilon definition and mode calculation:

from femwell.maxwell.waveguide import Mode
mode_9_0_1 = Mode(modes[0].frequency,
                  modes[0].k,
                  modes[0].basis_epsilon_r,
                  modes[0].epsilon_r,
                  modes[0].basis,
                  np.loadtxt('E_9_0_1.txt', dtype = complex),
                  np.loadtxt('H_9_0_1.txt', dtype = complex))

mode_9_1_1 = Mode(modes[0].frequency,
                  modes[0].k,
                  modes[0].basis_epsilon_r,
                  modes[0].epsilon_r,
                  modes[0].basis,
                  np.loadtxt('E_9_1_1.txt', dtype = complex),
                  np.loadtxt('H_9_1_1.txt', dtype = complex))

print('np.all(np.isclose(mode_9_0_1.H, mode_9_1_1.H)):', np.all(np.isclose(mode_9_0_1.H, mode_9_1_1.H)))
print('np.all(np.isclose(mode_9_0_1.E, mode_9_1_1.E)):', np.all(np.isclose(mode_9_0_1.E, mode_9_1_1.E)))

@Functional(dtype=np.complex64)
def current_form(w):
    """
    What this does is it takes the normal vector to the boundary and rotates it 90deg
    Then takes the inner product with the magnetic field in it's complex form
    """
    return inner(np.array([w.n[1], -w.n[0]]), w.H)

@Functional(dtype=np.complex64)
def voltage_form(w):
    """
    What this does is it takes the normal vector to the boundary and rotates it 90deg
    Then takes the inner product with the electric field in it's complex form
    """

    return -inner(np.array([w.n[1], -w.n[0]]), w.E)

conductor = "metal_sig_interface"
line = "path_integral"

for mode in [mode_9_0_1, mode_9_1_1]:

    p0 = calculate_scalar_product(
        mode.basis, np.conjugate(mode.E), mode.basis, np.conjugate(mode.H)
    )  # The conjugate is due to femwell internal definition as conj(E) x H

    (ht, ht_basis), (hz, hz_basis) = mode.basis.split(mode.H)
    facet_basis = ht_basis.boundary(facets=mesh.boundaries[conductor])
    i0 = current_form.assemble(facet_basis, H=facet_basis.interpolate(ht))

    (et, et_basis), (ez, ez_basis) = mode.basis.split(mode.E)
    facet_basis = et_basis.boundary(facets=mesh.boundaries[line])
    v0 = voltage_form.assemble(facet_basis, E=facet_basis.interpolate(et))

    print('Power:', p0)
    print('v0:', v0)
    print('i0:', i0)

This is put into a file called post_processing.py. After I run the following cmd:

pip install scikit-fem==9.1.1 & pip uninstall meshio & pip install meshio & python post_process.py & pip install scikit-fem==9.0.1 & pip uninstall meshio & pip install meshio & python post_process.py

With this, I end up with:

np.all(np.isclose(mode_9_0_1.H, mode_9_1_1.H)): True np.all(np.isclose(mode_9_0_1.E, mode_9_1_1.E)): True Power: (1-0.0052765195478439554j) v0: (13.536407914701542+6.973158551631229j) i0: (0.056362422778659514+0.028802282930216014j) Power: (1-0.0052765195478439554j) v0: (13.536407914701542+6.973158551631229j) i0: (0.056362422778659514+0.028802282930216014j)

np.all(np.isclose(mode_9_0_1.H, mode_9_1_1.H)): True np.all(np.isclose(mode_9_0_1.E, mode_9_1_1.E)): True Power: (1-0.0052765195478439554j) v0: (13.536407914701542+6.973158551631229j) i0: (0.05832818175174991+0.030265121214258277j) Power: (1-0.0052765195478439554j) v0: (13.536407914701542+6.973158551631229j) i0: (0.05832818175174991+0.030265121214258277j)

Notice how every value stays the same except the i0 calculation.

The other thing I'm wondering: Is one of the line integrals you calculate along an interface at which the epsilon changes discontinuously? (I.e. the surface of the conductor?)

The v0 line integral goes from the middle of the conductor to the edge of the simulation region where the boundary condition is set to metal. The i0 line integral is along a line that starts at the left boundary not set to metal and loops around the conductor up to the boundary again. However, the line integral which is calculated where epsilon changes discontinuously (the v0) remains the same in both versions. It's actually the i0 which is all defined in the air that is giving an issue

[duarte-jfs]

The error reported above stems from an optimization done from version 9.0.1 to 9.1.1 in skfem on skfem.io.meshio.from_meshio(). This is being solved through issue https://github.com/kinnala/scikit-fem/issues/1150 and on the PR https://github.com/kinnala/scikit-fem/pull/1152 . The issue was that the OrientedBoundary was not properly oriented, and therefore caused errors only on the line integral.

9.0.1

PI model : |Z0| = 61.77 Ohm | angle(Z0) = -0.01 pi radians PV model : |Z0| = 60.19 Ohm | angle(Z0) = -0.01 pi radians VI model : |Z0| = 60.97 Ohm | angle(Z0) = -0.00 pi radians

9.1.1

PI model : |Z0| = 76.77 Ohm | angle(Z0) = -0.01 pi radians PV model : |Z0| = 60.19 Ohm | angle(Z0) = -0.01 pi radians VI model : |Z0| = 67.97 Ohm | angle(Z0) = -0.00 pi radians

9.1.1 with fix on oriented boundary

PI model : |Z0| = 61.77 Ohm | angle(Z0) = -0.01 pi radians PV model : |Z0| = 60.19 Ohm | angle(Z0) = -0.01 pi radians VI model : |Z0| = 60.97 Ohm | angle(Z0) = -0.00 pi radians