Rectangular Cavity Code not working

[DatKnight]

Hi All,

I'd like to get started in openEMS. As a first try I wanted to make a resonant cavity model in octave. My initial thought was to use the Rectangular Waveguide tutorial as an example and try to write it as I go, referencing the wiki.
Here is my code so far, I just want to create the simulation and CSX model, then write it to a folder so I can look at it with Paraview.

close all
clear
clc

%Specify basic unit for the simulation
unit = 10^-6;

%Specify center frequency and the range over which to scan
freq0 = 900*10^6;
dfreq = 50*10^6;

%Calculate wavelength
lambda = 3*10^8/freq0/unit;

%Cavity dimensions
width = lambda/2/unit;
length = width;

%Define excitation mode, we are using RectWaveguide Ports so will have to define waveguide modes
TEmode = 'TE10';

%Mesh resolution: 30 cells per wavelength
mesh_res = lambda./[30 30 30];

%Initialise FDTD
FDTD = InitFDTD('NrTS',1e4, 'OverSampling', 5);
FDTD = SetGaussExcite(FDTD,freq0-dfreq,freq0+dfreq);

%Set our boundary conditions
BC = [0 0 0 0 0 0]; %Perfect Electrical Conductor in all dir to act as cavity walls
FDTD = SetBoundaryCond(FDTD,BC);

%Create CSXCAD mesh within cavity walls
CSX = InitCSX()
mesh.x = SmoothMeshLines([0 width], mesh_res(1));
mesh.y = SmoothMeshLines([0 width], mesh_res(2));
mesh.z = SmoothMeshLines([0 length], mesh_res(3));
CSX = DefineRectGrid(CSX, unit,mesh);

%Set up ports for measuring S params
%Create port 1 as a RectWaveGuidePort placed at 13 units from end of cavity at z origin dir
start=[mesh.x(1)   mesh.y(1)   mesh.z(11)];
stop =[mesh.x(end) mesh.y(end) mesh.z(15)];
[CSX, port{1}] = AddRectWaveGuidePort( CSX, 0, 1, start, stop, 'z', width*unit, width*unit, TE_mode, 1);

%Create port 2 as a RectWaveGuidePort placed at 13 units from end of cavity in +z dir, set excitation to 0 to make it passive
start=[mesh.x(1)   mesh.y(1)   mesh.z(end-13)];
stop =[mesh.x(end) mesh.y(end) mesh.z(end-14)];
[CSX, port{2}] = AddRectWaveGuidePort( CSX, 0, 2, start, stop, 'z', width*unit, width*unit, TE_mode, 0);

%Write model to subfolder
simpath = 'tmp_cav';
simcsx = 'rect_cav.xml';
mkdir(simpath);
WriteOpenEMS([simpath simcsx], FDTD, CSX);

But when I write it the only output I get is this:

QSocketNotifier: Can only be used with threads started with QThread

(octave-gui:255360): Gtk-WARNING **: 15:40:22.499: Theme parsing error: <data>:1:0: Expected a valid selector
CSX =

  scalar structure containing the fields:

    Properties = [](0x0)
    ATTRIBUTE =

      scalar structure containing the fields:

        CoordSystem = 0

        warning: found resolution larger than max_res
        warning: called from
            CheckMesh at line 32 column 5
            SmoothMeshLines at line 121 column 5
            cavity at line 35 column 8

And then it hangs.

Is there anything I'm doing wrong? I'd appreciate any help on this as I find it a little confusing. I'm not sure how I'm creating a larger resolution than max_res.

[VolkerMuehlhaus]

I am missing the CSXGeomPlot part in your model, for visual inspection of model and mesh. Did you ever check that?

When I check the numbers in your model, these look strange. The mesh count would be excessive as well.

lambda = 3.3333e+05
width = 1.6667e+11
length = 1.6667e+11
mesh_res =

1.1111e+04 1.1111e+04 1.1111e+04

The mistake is here, you already calculate lambda in model units and the divide width by unit factor again.

%Calculate wavelength
lambda = 3*10^8/freq0/unit;

%Cavity dimensions
width = lambda/2/unit;

[DatKnight]

Hi thank you that works and everything runs now. I've only just noticed CSXGeomPlot in the relevant tutorial, thanks!

I'm having another minor problem though, when I call up the CAD renderer with CSXGeomPlot it loads everything, but the window is very laggy and I get the following error in the console whenever I move the model about or turn features on or off:

2024-04-28 18:23:34.605 ( 179.504s) [        F1350000]vtkOpenGLRenderWindow.c:511    ERR| vtkGenericOpenGLRenderWindow (0x2331a10): GLEW could not be initialized: Unknown error
2024-04-28 18:23:34.627 ( 179.526s) [        F1350000]vtkOpenGLRenderWindow.c:511    ERR| vtkGenericOpenGLRenderWindow (0x2331a10): GLEW could not be initialized: Unknown error
2024-04-28 18:23:34.663 ( 179.562s) [        F1350000]vtkOpenGLRenderWindow.c:511    ERR| vtkGenericOpenGLRenderWindow (0x2331a10): GLEW could not be initialized: Unknown error
2024-04-28 18:23:34.686 ( 179.585s) [        F1350000]vtkOpenGLRenderWindow.c:511    ERR| vtkGenericOpenGLRenderWindow (0x2331a10): GLEW could not be initialized: Unknown error
2024-04-28 18:23:34.700 ( 179.599s) [        F1350000]vtkOpenGLRenderWindow.c:511    ERR| vtkGenericOpenGLRenderWindow (0x2331a10): GLEW could not be initialized: Unknown error
2024-04-28 18:23:34.731 ( 179.630s) [        F1350000]vtkOpenGLRenderWindow.c:511    ERR| vtkGenericOpenGLRenderWindow (0x2331a10): GLEW could not be initialized: Unknown error
2024-04-28 18:23:34.766 ( 179.664s) [        F1350000]vtkOpenGLRenderWindow.c:511    ERR| vtkGenericOpenGLRenderWindow (0x2331a10): GLEW could not be initialized: Unknown error
2024-04-28 18:23:34.792 ( 179.691s) [        F1350000]vtkOpenGLRenderWindow.c:511    ERR| vtkGenericOpenGLRenderWindow (0x2331a10): GLEW could not be initialized: Unknown error
2024-04-28 18:23:34.824 ( 179.723s) [        F1350000]vtkOpenGLRenderWindow.c:511    ERR| vtkGenericOpenGLRenderWindow (0x2331a10): GLEW could not be initialized: Unknown error
2024-04-28 18:23:34.839 ( 179.738s) [        F1350000]vtkOpenGLRenderWindow.c:511    ERR| vtkGenericOpenGLRenderWindow (0x2331a10): GLEW could not be initialized: Unknown error
2024-04-28 18:23:34.875 ( 179.774s) [        F1350000]vtkOpenGLRenderWindow.c:511    ERR| vtkGenericOpenGLRenderWindow (0x2331a10): GLEW could not be initialized: Unknown error
2024-04-28 18:23:34.908 ( 179.807s) [        F1350000]vtkOpenGLRenderWindow.c:511    ERR| vtkGenericOpenGLRenderWindow (0x2331a10): GLEW could not be initialized: Unknown error
2024-04-28 18:23:34.918 ( 179.817s) [        F1350000]vtkOpenGLRenderWindow.c:511    ERR| vtkGenericOpenGLRenderWindow (0x2331a10): GLEW could not be initialized: Unknown error
2024-04-28 18:23:34.936 ( 179.835s) [        F1350000]vtkOpenGLRenderWindow.c:511    ERR| vtkGenericOpenGLRenderWindow (0x2331a10): GLEW could not be initialized: Unknown error
2024-04-28 18:23:34.976 ( 179.875s) [        F1350000]vtkOpenGLRenderWindow.c:511    ERR| vtkGenericOpenGLRenderWindow (0x2331a10): GLEW could not be initialized: Unknown error
2024-04-28 18:23:35.001 ( 179.900s) [        F1350000]vtkOpenGLRenderWindow.c:511    ERR| vtkGenericOpenGLRenderWindow (0x2331a10): GLEW could not be initialized: Unknown error

This looks like an issue with one of the libraries so I don't think it's a openEMS issue, but do you have any idea?

[DatKnight]

I've had a go at trying to calculate the S11 and S21 values for the cavity using the rectangular waveguide tutorial as a guide, but for some reason the calcPort function seems to be returning single values for the ports' measured voltage/current amplitudes and so I get just a blank plot:

Here's my code:

close all
clear
clc

%Specify basic unit for the simulation
unit = 10^-6;

%Specify center frequency and the range over which to scan
freq0 = 900*10^6;
dfreq = 50*10^6;

%Calculate wavelength
lambda = 3*10^8/freq0/unit;

%Cavity dimensions
width = lambda/2;
length = width;

%Define excitation mode, we are using RectWaveguide Ports so will have to define waveguide modes
TE_mode = 'TE10';

%Mesh resolution: 30 cells per wavelength
mesh_res = lambda./[30 30 30];

%Initialise FDTD
FDTD = InitFDTD('NrTS',1e4, 'OverSampling', 5);
FDTD = SetGaussExcite(FDTD,freq0,freq0-dfreq);

%Set our boundary conditions
BC = [0 0 0 0 0 0]; %Perfect Electrical Conductor in all dir to act as cavity walls
FDTD = SetBoundaryCond(FDTD,BC);

%Create CSXCAD mesh within cavity walls
CSX = InitCSX();
mesh.x = SmoothMeshLines([0 width], mesh_res(1));
mesh.y = SmoothMeshLines([0 width], mesh_res(2));
mesh.z = SmoothMeshLines([0 length], mesh_res(3));
CSX = DefineRectGrid(CSX, unit,mesh);

%Set up ports for measuring S params
%Create port 1 as a RectWaveGuidePort placed at 13 units from end of cavity at z origin dir
start=[mesh.x(1)   mesh.y(1)   mesh.z(11)];
stop =[mesh.x(end) mesh.y(end) mesh.z(15)];
[CSX, port{1}] = AddRectWaveGuidePort( CSX, 0, 1, start, stop, 'z', width*unit, width*unit, TE_mode, 1);

%Create port 2 as a RectWaveGuidePort placed at 13 units from end of cavity in +z dir, set excitation to 0 to make it passive
start=[mesh.x(1)   mesh.y(1)   mesh.z(end-13)];
stop =[mesh.x(end) mesh.y(end) mesh.z(end-14)];
[CSX, port{2}] = AddRectWaveGuidePort( CSX, 0, 2, start, stop, 'z', width*unit, width*unit, TE_mode, 0);

%Write model to subfolder
simpath = 'tmp_cav';
simcsx = 'rect_cav.xml';
rmdir(simpath, 's');
mkdir(simpath);
WriteOpenEMS([simpath '/' simcsx], FDTD, CSX);

%View geometry
%CSXGeomPlot([simpath '/' simcsx]);

%Run OpenEMS
RunOpenEMS(simpath, simcsx);

%Post process
freq = linspace(freq0-dfreq,freq0+dfreq,201);
%Calculate voltages & currents for ports
port = calcPort(port, simpath, freq0);

%Calculate S11 and S21
s11 = port{1}.uf.ref./ port{1}.uf.inc;
s21 = port{2}.uf.ref./ port{1}.uf.inc;

%Draw s11 & s21 figures
figure
plot(freq*1e-6,20*log10(abs(s11)),'k-','Linewidth',2);
xlim([freq(1) freq(end)]*1e-6);
grid on;
hold on;
plot(freq*1e-6,20*log10(abs(s21)),'r--','Linewidth',2);
l = legend('S_{11}','S_{21}','Location','northeast');
set(l,'FontSize',12);
ylabel('S-Parameter (dB)','FontSize',12);
xlabel('frequency (MHz) \rightarrow','FontSize',12);

The values of port{1}.uf.ref, port{1}.uf.inc, and port{2}.uf.ref are:

octave:12> port{1}.uf.ref
ans =  9.0292e-15 + 2.8337e-15i
octave:13> port{2}.uf.ref
ans =  3.7145e-15 + 3.7804e-15i
octave:14> port{1}.uf.inc
ans = -8.9461e-15 - 3.0048e-15i

I assume this has something to do with the way I've set up the simulation but I'm not sure how I could have done it wrong? I've double checked and changed the way I setup SetGaussExcite, but other than that it all seems ok, no reason for there to be only one value returned.

I've tried adding a dump box as in the tutorial, but that doesn't change anything.

[VolkerMuehlhaus]

The reason for the single frequency result is here:

port = calcPort(port, simpath, freq0);

Instead of freq0, you need to use the freq vector

[DatKnight]

Thanks so much! It's working now:

Evidently I need to figure out where the best port placement is now, is there any particular reason the ports are made 4 units wide in the example or can I make them any size? I imagine in reality they would need to be less than 1/10th of the wavelength but idk how that applies in the simulation.

[DatKnight]

Hi All,

I've gone back to the drawing board and redesigned the cavity taking into account the guide wavelength. I've placed the input and output ports at a quarter wavelength from the ends of the cavity, which should be around the optimum, and plotted S11, S21, and the Et field dump.

My issue is that first off, during the simulation I get this warning:

RunFDTD: Warning: Max. number of timesteps was reached before the end-criteria of -50dB was reached... 
	You may want to choose a higher number of max. timesteps... 
Time for 50000 iterations with 140864.00 cells : 253.66 sec
Speed: 27.77 MCells/s

Then after that the S11/S21 plot looks like this:

And the Et field dump doesn't seem to actually resonate, but I think that might be because the excitation doesn't last long enough, and I'm not entirely sure that it is the right thing to look at in the first place.

Am I barking up the wrong tree? Should I be looking to do something else? I'm not really certain where to go from here.

My code is:

close all
clear
clc

%Specify basic unit for the simulation in m
unit = 10^-6;

%Specify center frequency and the range over which to scan
freq0 = 900*10^6;
dfreq = 50*10^6;

%Calculate wavelength
lambda = 3*10^8/freq0;

%Cavity dimensions in m
width = 0.5;
height = 0.25;
length = 0.5755;

%Convert dimensions and lambda to simulation units
width = width/unit;
length = length/unit;
height = height/unit;
lambda = lambda/unit;

%Define excitation mode, we are using RectWaveguide Ports so will have to define waveguide modes
TE_mode = 'TE10';

%Create mesh wavelength to define our mesh resolution
mesh_lambda = 40;

%Mesh resolution: 40 cells per wavelength
mesh_res = lambda./[mesh_lambda mesh_lambda mesh_lambda];

%Initialise FDTD
FDTD = InitFDTD('NrTS',5e4, 'OverSampling', 5);
FDTD = SetGaussExcite(FDTD,freq0,freq0-dfreq);

%Set our boundary conditions
BC = [0 0 0 0 0 0]; %Perfect Electrical Conductor in all dir to act as cavity walls
FDTD = SetBoundaryCond(FDTD,BC);

%Create CSXCAD mesh within cavity walls
CSX = InitCSX();
mesh.x = SmoothMeshLines([0 width], mesh_res(1));
mesh.y = SmoothMeshLines([0 height], mesh_res(2));
mesh.z = SmoothMeshLines([0 length], mesh_res(3));
CSX = DefineRectGrid(CSX, unit,mesh);

%Set up ports for measuring S params

%Create port 1 as a RectWaveGuidePort placed a quarter wavelength from end of cavity at z origin dir, set excitation to 1 to make it active
start=[mesh.x(1)   mesh.y(1)   mesh.z(mesh_lambda/4)];
stop =[mesh.x(end) mesh.y(end) mesh.z(mesh_lambda/4 + 1)];
[CSX, port{1}] = AddRectWaveGuidePort( CSX, 0, 1, start, stop, 'z', width*unit, width*unit, TE_mode, 1);

%Create port 2 as a RectWaveGuidePort placed a quarter wavelength units from end of cavity in +z dir, set excitation to 0 to make it passive
start=[mesh.x(1)   mesh.y(1)   mesh.z(end - mesh_lambda/4)];
stop =[mesh.x(end) mesh.y(end) mesh.z(end - mesh_lambda/4 - 1)];
[CSX, port{2}] = AddRectWaveGuidePort( CSX, 0, 2, start, stop, 'z', width*unit, width*unit, TE_mode, 0);

%Create dump box to record time-domain E field results
CSX = AddDump(CSX,'Et','FileType',1,'SubSampling','2,2,2');
start = [mesh.x(1)   mesh.y(1)   mesh.z(1)];
stop  = [mesh.x(end) mesh.y(end) mesh.z(end)];
CSX = AddBox(CSX,'Et',0 , start,stop);

%Write model to subfolder
simpath = 'tmp_cav';
simcsx = 'rect_cav.xml';
rmdir(simpath, 's');
mkdir(simpath);
WriteOpenEMS([simpath '/' simcsx], FDTD, CSX);

%View geometry
CSXGeomPlot([simpath '/' simcsx]);

%Run OpenEMS
RunOpenEMS(simpath, simcsx);

%Post process
freq = linspace(freq0-dfreq,freq0+dfreq,201);
%Calculate voltages & currents for ports
port = calcPort(port, simpath, freq);

%Calculate S11 and S21
s11 = port{1}.uf.ref./ port{1}.uf.inc;
s21 = port{2}.uf.ref./ port{1}.uf.inc;


%Draw s11 & s21 figures
figure
plot(freq*1e-6,20*log10(abs(s11)),'k-','Linewidth',2);
xlim([freq(1) freq(end)]*1e-6);
grid on;
hold on;
plot(freq*1e-6,20*log10(abs(s21)),'r--','Linewidth',2);
l = legend('S_{11}','S_{21}','Location','northeast');
set(l,'FontSize',12);
ylabel('S-Parameter (dB)','FontSize',12);
xlabel('frequency (MHz) \rightarrow','FontSize',12);

%Plot field dumps
figure
dump_file = [simpath '/Et.h5'];
PlotArgs.slice = {width/2*unit width/2*unit 0};
PlotArgs.pauseTime=0.01;
PlotArgs.component=0;
PlotArgs.Limit = 'auto';
PlotHDF5FieldData(dump_file, PlotArgs)

[DatKnight]

I'm not sure how to check those, the energy never really went below the -50dB figure:

So I assume your comment on the excitation method is correct. I'll have to look into the documentation a bit more on how to set up a more accurate representation of real life excitation. I just went with the one presented in the waveguide tutorial as I'm not familiar with the software.

[VolkerMuehlhaus]

Yes, it seems that energy never decays.

Regarding the excitation: remember that a waveguide transmission line is very different from a cavity. For the cavity you want to inject energy, but have only weak coupling to the source and load lines, so that energy keeps circulating in the resonator for a while. The topic here is loaded Q factor vs. unloaded Q

When using waveguide ports, you can start from the waveguide example and insert two apertures to build a resonator, as shown in the example on the very top here (other software, but the picture shows the structure):
https://mwt-www.e-technik.uni-ulm.de/world/data/Web/multimedia/hf-anim/mefisto/filter.html

Another note: for the waveguide transmission line dimensions must match the operating frequency. If check your dimensions (0.5m width) with an online calculator, that would work for a frequency range 400-600 MHz in the fundamental mode. I think your simulation frequency doesn't match your dimensions in the file above.

[biergaizi]

To debug a malfunctioning simulation, enable field dump to see the raw electromagnetic field in ParaView (e.g. this example from pyems). Sometimes, something as simple a misaligned excitation port can create never-ending signal injections.

Also, a simulation speed of 6 MC/s is definitely not correct, normally it should be like 100 MC/s, not 6 MC/s, this may also be the reason of convergence failure since it was running too slowly. This may indicate compiler version, compiler settings, or simulation engine parameter problems.

[DatKnight]

I assumed that the reason the simulation speed is so slow is because I'm running it on a Dual Core i7-7600U?

I've added a field dump and dump box around the cavity:

But I don't get any .vtr files showing up in the dump folder:

[biergaizi]

I assumed that the reason the simulation speed is so slow is because I'm running it on a Dual Core i7-7600U?

Sorry for the confusion. I don't notice your simulation domain has just 4096 cells. In this case, it's completely normal, such a tiny domain size means the program spend little time in the simulation kernel, so its runtime is dominated by software overheads instead. I can reproduce this slow speed on my workstation too.

But I don't get any .vtr files showing up in the dump folder

See the documentation of Dump Box Property.

    FileType:
        0: vtk-file dump (default)
        1: hdf5-file dump

There's a file Et.h5, which is the field dump file in HDF5 format. To use the vtr format, set FileType to 0.