Numeric errors lead to amplification/solution instability

How to use openEMS. Discussion on examples, tutorials etc

Moderator: thorsten

Hale_812
Posts: 173
Joined: Fri 13 May 2016, 02:54

Numeric errors lead to amplification/solution instability

Post by Hale_812 » Thu 23 Jun 2016, 10:34

Looks like dense mesh (used for oblique surfaces and coax) causes heavy amplification of the signal and solution instability.

This often happens when the antenna is not fed perfectly. The power is reflected back to boundary, but the error amplification is stronger than PML absorption. Then, even far beyound the pulse time, the feeding coax becomes red, power leaks from coax at the PML boundary to nearby space... and all the space is slowly filled with E-field of incredible strength. The solution can not stop until INF, or NaN vaues are reached. It runs for days.

Similar thing can happen when simulating simple horn antenna fed with perfect rect.port. The pulse has already passed the feeding waveguide, but the power has started amplifying at the antenna edges with enormous gain. At the moment I am observing the situation and will try to upload vector animation of the problem. Seems like 5 PML_8 boundaries are still stronger than the generation, especially when 6th back boundary is PML_16... but still it does not look normal.

see the partial simulation here https://sendvid.com/77ucxmn3


And here is the self-amplifying and leaking at PML coax. http://sendvid.com/voscm93r
Last edited by Hale_812 on Fri 24 Jun 2016, 03:52, edited 1 time in total.

thorsten
Posts: 1460
Joined: Mon 27 Jun 2011, 12:26

Re: Numeric errors lead to amplification/solution instabilit

Post by thorsten » Thu 23 Jun 2016, 19:33

Hi,

the FDTD method itself is unconditionally stable, regardless of the mesh density.
But what is not stable is the PML and there is no stable PML possible, since this concept does not exist in reality.
What the PML cannot handle are evanescence fields. Propagating waves are not an issue. Thus it is usually enough to place a PML far enough away from any structure that might cause evanescence fields.

regards
Thorsten

thorsten
Posts: 1460
Joined: Mon 27 Jun 2011, 12:26

Re: Numeric errors lead to amplification/solution instabilit

Post by thorsten » Thu 23 Jun 2016, 19:35

One more note. Use only! pml_8. Every other PML is not tested and I think they are all unstable... One day I need to have a closer look, but really, you never need anything else than pml_8...

Hale_812
Posts: 173
Joined: Fri 13 May 2016, 02:54

Re: Numeric errors lead to amplification/solution instabilit

Post by Hale_812 » Fri 24 Jun 2016, 03:34

Thank you for the notes. I will try to make the problem larger... the problem is that it is already taking 7GB of Ram and twice of that during operator preparation.

BTW, In the leaking coax video, PML8 was used. In horn antenna video, all walls are PML8, only the back wall is PML16 (it serves as the feeder termination, so I tried to improve its absorption properties). Looks, like PML16 works better...I will try to move boundaries farther away.

I have a question. Do you know, what is the fundamental difference between the radiation boundary of HFSS and classic PML. I mean, I understand that PML is simulated as an antireflective coating with a high loss tangent, or kind of that. But what are BC for a radiation boundary? Why is there no such in OpenEMS FEM?

P.S. for me, it looks like the problem with FEM: with a mesh of a certain density, it will converge, but with denser mesh, the weak field rounding error will accumulate giving slightly worse solution, worse dynamic range etc. It was a big problem when I was calculating coupled chiral lattices... the 0xPi and 2xPi polarization solution were just not fitting together... small difference in mesh, end the noise was driving resonances away. For me, the situation with growing field looks similar.
Last edited by Hale_812 on Fri 24 Jun 2016, 10:07, edited 1 time in total.

Hale_812
Posts: 173
Joined: Fri 13 May 2016, 02:54

Re: Numeric errors lead to amplification/solution instabilit

Post by Hale_812 » Fri 24 Jun 2016, 04:24

OK. And now I see another thing that confused me.
During simulation I got "Excitation signal length is: 180 timesteps"
BUT! While simulating, I've checked the profile, and the port was still oscillating at high power at the "Timestep: 4298"...

So, definitely, I should make the excitation pulse band wider and the airbox 8 times larger.

thorsten
Posts: 1460
Joined: Mon 27 Jun 2011, 12:26

Re: Numeric errors lead to amplification/solution instabilit

Post by thorsten » Fri 24 Jun 2016, 18:29

But what are BC for a radiation boundary? Why is there no such in OpenEMS FEM?
Both the MUR and PML boundary conditions are radiation boundary conditions. I think in FEM they use the same kind of BC's. It's just a question of naming...

regards
Thorsten

Hale_812
Posts: 173
Joined: Fri 13 May 2016, 02:54

Re: Numeric errors lead to amplification/solution instabilit

Post by Hale_812 » Mon 27 Jun 2016, 06:25

Yes, I looked at it precisely, I guess that is MUR.
It was said to me many times that HFSS Radiation boundary is better than PML, but works only far in far zone. If you do not wish to increase mesh density at the problem boundaries, it is the best choice. PML can be brought much closer, but it is frequency dependent (I guess, because its mesh is frequency dependent), always give some amount of reflection and hell slow to simulate.

Hale_812
Posts: 173
Joined: Fri 13 May 2016, 02:54

Re: Numeric errors lead to amplification/solution instabilit

Post by Hale_812 » Fri 15 Jul 2016, 06:16

So, here is another project. Not quite an amplification of errors, I got rid of that in this simple setup (no coaxials, no port junctions, no band-gap structures, short feed)

But still, can't get an antenna pattern. Can't even guess what's wrong.
%note: if you set a half density, the pattern will be wrong, with a gain of 6dB below the real thing. I.e. the gain should be around 26dB. And the H-plane should be quite narrow, just a bit wider than E-plane below -10dBc level

Code: Select all

confirm_recursive_rmdir (0,"local")
close all
clear
clc

%% prepare simulation folder
Sim_Path = 'C:/Users/474/Documents/OpenEMS/';
Sim_Dir = 'tmp_SuperBigHorn';
Sim_CSX = 'SuperBigHorn.xml';
[status, message, messageid] = rmdir([Sim_Path Sim_Dir], 's' ); % clear previous directory
[status, message, messageid] = mkdir([Sim_Path Sim_Dir] ); % create empty simulation folder

%% setup the simulation
physical_constants;
unit = 1e-3; % all length in mm


% setup simulation quality
wave_subsample = 40; %resolution denominator - part of lambda. recommended in the range of 10-20!!!! FIX improved with 3-pass smoothing
stop_dB = -40; % stop calculation when reached wave decay -xx dB 
freespace = 1*4+1 ; %added free space X lambda/4


% frequency range of interest
f_start =  6.5e9;
%f_stop  =  15.1e9;
%f0 = 0.5*(f_start+f_stop); % frequency of interest
f0=8.5e9;
f_stop  =  2*f0-f_start;

lambda = c0 / (f0)/ unit; %wavelength

%Rect waveguide TE-mode definition
TE_mode = 'TE10';


%visual geometry parameters
a  = 22.86; % Waveguide width in x-direction
b = 10.16; % Waveguide height in y-direction
horn.width = 249;
horn.height = 184;
horn.edge = 567;

%Calculate important geometry parameters
horn.side_length = sqrt((horn.edge)^2-(horn.height-b)^2 /4); %side length
horn.top_length = sqrt((horn.edge)^2-(horn.width-a)^2 /4); %top length


% horn length in z-direction
horn.length = sqrt( (horn.edge)^2-(horn.height-b)^2 /4-(horn.width-a)^2 /4 );

%Feeder virtual parameters (WG is used instaed the probe)
horn.feed_length = 43.05; %feeder cavity
probe.indent = 7.05; %lambda/4 not used
probe.height = 0.73*lambda/4; %lambda/4 not used
horn.thickness = 4; %make it thicker to catch grid points for sure


%Horn opening angle in x, y
horn.angle = [asin(sqrt(((horn.width-a)^2)/(4*(horn.edge)^2-(horn.height-b)^2 ) ))  , asin(sqrt(((horn.height-b)^2)/(4*(horn.edge)^2-(horn.width-a)^2 )))]; % a-expansion b-expansion angles

%Initial size of the simulation box
SimBox = [freespace*lambda/4+a+2.*horn.length*tan(horn.angle(1))+horn.thickness.*2, freespace*lambda/4+b+2.*horn.length*tan(horn.angle(2))+horn.thickness.*2, freespace*lambda/6+horn.feed_length+horn.length+horn.thickness*2];




%% setup FDTD parameter & excitation function
FDTD = InitFDTD('EndCriteria', 10 .^(stop_dB./10)); % decay stop level 1e-4 is recommended
FDTD = SetGaussExcite(FDTD,0.5*(f_start+f_stop),0.5*(f_stop-f_start));
BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_16' 'PML_8'}; % boundary conditions
FDTD = SetBoundaryCond( FDTD, BC );

%% setup CSXCAD geometry & mesh
% currently, openEMS cannot automatically generate a mesh
max_res = lambda/ wave_subsample; % cell size mm: lambda/parts !!! not very actual because of 3-pass smooth refinement
CSX = InitCSX();

%create fixed lines for the simulation box, substrate and port
mesh.x = [-SimBox(1)/2, -a/2, a/2, SimBox(1)/2];
mesh.y = [-SimBox(2)/2,-b/2, b/2, SimBox(2)/2];
%create fixed lines for the simulation box and given number of lines inside the substrate
mesh.z = [-horn.feed_length, -horn.feed_length+probe.indent, 0, horn.length, SimBox(3)-horn.feed_length-horn.thickness ];


%% create horn
% horn feed rect waveguide
CSX = AddMetal(CSX, 'metal');
start = [-a/2-horn.thickness -b/2 -horn.feed_length];
stop  = [-a/2                 b/2 0];
CSX = AddBox(CSX,'metal',10,start,stop);
start = [a/2+horn.thickness -b/2 -horn.feed_length];
stop  = [a/2                 b/2 0];
CSX = AddBox(CSX,'metal',10,start,stop);
start = [-a/2-horn.thickness b/2+horn.thickness -horn.feed_length];
stop  = [ a/2+horn.thickness b/2                0];
CSX = AddBox(CSX,'metal',10,start,stop);
start = [-a/2-horn.thickness -b/2-horn.thickness -horn.feed_length];
stop  = [ a/2+horn.thickness -b/2                0];
CSX = AddBox(CSX,'metal',10,start,stop);

p=[]; clear p;
% horn opening
%XZ plane: x                                     z
p(2,1) = -a/2                               ;p(1,1) = 0;           %point 1
p(2,2) = -a/2-horn.length*tan(horn.angle(1));p(1,2) = horn.length/cos(horn.angle(2)); %point 2
p(2,3) = a/2+horn.length*tan(horn.angle(1)) ;p(1,3) = horn.length/cos(horn.angle(2)); %point 3
p(2,4) = a/2                                ;p(1,4) = 0;           %point 4
CSX = AddLinPoly( CSX, 'metal', 10, 'y', 0, p, horn.thickness, 'Transform', {'Rotate_X',-horn.angle(2),'Translate',['0,' num2str(b/2) ',0']});
CSX = AddLinPoly( CSX, 'metal', 10, 'y', 0, p, -horn.thickness, 'Transform', {'Rotate_X',horn.angle(2),'Translate',['0,' num2str(-b/2) ',0']});

p=[]; clear p;
%ZY plane: Z                                     y
p(1,1) = b/2                                  ;p(2,1) = 0;           %point 1
p(1,2) = b/2+horn.thickness/cos(horn.angle(2));p(2,2) = 0; %point 2
p(1,3) = b/2+horn.length*tan(horn.angle(2))+horn.thickness*cos(horn.angle(2)); p(2,3) = horn.length/cos(horn.angle(1))-horn.thickness*sin(horn.angle(2))/cos(horn.angle(1)); %point 3
p(1,4) = b/2+horn.length*tan(horn.angle(2));p(2,4) = horn.length/cos(horn.angle(1)); %point 4
p(1,5) = -p(1,4) ;p(2,5) = p(2,4); %point 5
p(1,6) = -p(1,3); p(2,6) = p(2,3); %point 6
p(1,7) = -p(1,2);p(2,7) = p(2,2); %point 7
p(1,8) = -p(1,1);p(2,8) = p(2,1); %point 8
CSX = AddLinPoly( CSX, 'metal', 10, 'x', 0, p, -horn.thickness, 'Transform', {'Rotate_Y',-horn.angle(1),'Translate',[num2str(-a/2) ', 0, 0']});
CSX = AddLinPoly( CSX, 'metal', 10, 'x', 0, p, horn.thickness, 'Transform', {'Rotate_Y',horn.angle(1),'Translate',[num2str(a/2) ', 0, 0']});


%Refining the grid
mesh.x=sort(mesh.x); mex=mesh.x(1);
for ix=2:size(mesh.x,2)
  if mesh.x(ix)>(mex(end)+max_res/1000); %removes duplicates
    mex(end+1)=mesh.x(ix);
  endif
endfor
mesh.x=mex;
mex=[]; clear mex;
mesh.x=sort(mesh.x);
mesh.x = SmoothMeshLines( mesh.x, max_res, 1.4); % create a smooth mesh between specified fixed mesh lines
mex=mesh.x(1);
for ix=2:size(mesh.x,2)
  if mesh.x(ix)>(mex(end)+max_res/15); %removes extra-dense lines
    mex(end+1)=mesh.x(ix);
  endif
endfor
mesh.x=mex; mex=[]; clear mex;

mesh.y=sort(mesh.y); mey=mesh.y(1);
for iy=2:size(mesh.y,2)
  if mesh.y(iy)>(mey(end)+max_res/1000); %removes duplicates
    mey(end+1)=mesh.y(iy);
  endif
endfor
mesh.y=mey;
mey=[]; clear mey;
mesh.y=sort(mesh.y);
mesh.y = SmoothMeshLines( mesh.y, max_res, 1.4 );
mey=mesh.y(1);
for iy=2:size(mesh.y,2)
  if mesh.y(iy)>(mey(end)+max_res/15); %removes extra-dense lines
    mey(end+1)=mesh.y(iy);
  endif
endfor
mesh.y=mey; mey=[]; clear mey;

mesh.z=sort(mesh.z);
mesh.z = SmoothMeshLines( mesh.z, max_res, 1.4 );
mez=mesh.z(1);
for iz=2:size(mesh.z,2)
  if mesh.z(iz)>(mez(end)+max_res/15); %removes extra-dense lines
    mez(end+1)=mesh.z(iz);
  endif
endfor
mesh.z=mez; mez=[]; clear mez;

%2x times refinement
mesh.x = SmoothMeshLines( mesh.x, max_res, 1.4); % create a smooth mesh between specified fixed mesh lines
mesh.y = SmoothMeshLines( mesh.y, max_res, 1.4 );
mesh.z = SmoothMeshLines( mesh.z, max_res, 1.4 );
%3x times refinement
mesh.x = SmoothMeshLines( mesh.x, max_res, 1.4); % create a smooth mesh between specified fixed mesh lines
mesh.y = SmoothMeshLines( mesh.y, max_res, 1.4 );
mesh.z = SmoothMeshLines( mesh.z, max_res, 1.4 );




CSX = DefineRectGrid( CSX, unit, mesh );


% horn aperture
A = (a + 2*tan(horn.angle(1))*horn.length)*unit * (b + 2*tan(horn.angle(2))*horn.length)*unit;

%% apply the excitation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%WG port
start=[-a/2, -b/2, -horn.feed_length ];
stop =[ a/2,  b/2, -horn.feed_length+probe.indent ];
[CSX, port] = AddRectWaveGuidePort( CSX, 0, 1, start, stop, 2, a*unit, b*unit, TE_mode, 1);

%% nf2ff calc
start = [mesh.x(11) mesh.y(11) mesh.z(11)];
stop  = [mesh.x(end-10) mesh.y(end-10) mesh.z(end-10)];
[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop, 'Directions', [1 1 1 1 0 1]);

CSX = AddDump(CSX,'Et','FileType',0,'SubSampling','4,4,4');
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 openEMS compatible xml-file
WriteOpenEMS([Sim_Path Sim_Dir '/' Sim_CSX], FDTD, CSX);

%% show the structure
%CSXGeomPlot([Sim_Path Sim_Dir '/' Sim_CSX],"--RenderDiscMaterial");
if isOctave()
    fflush(stdout);
end
%!!!!!  FIX THE PATH. Very rough workaround for blocking csxCAD call
system(['start /B /I c:\OpenEMS\AppCSXCAD.exe --disableEdit --RenderDiscMaterial ' [Sim_Path Sim_Dir '/' Sim_CSX]]);

%% run openEMS
RunOpenEMS([Sim_Path Sim_Dir], Sim_CSX);

%% postprocessing & do the plots
freq = linspace(f_start,f_stop,201);

port = calcPort(port, [Sim_Path Sim_Dir], freq);

Zin = port.uf.tot ./ port.if.tot;
s11 = port.uf.ref ./ port.uf.inc;
P_in = 0.5 * port.uf.inc .* conj( port.if.inc ); % antenna feed power
disp( ['Port impedance: Zin = ' num2str(sum(real(Zin))/length(Zin)) ' Ohm'] );
disp( ['Port best reflection: S11 = ' num2str(min(20*log10(abs(s11)))) ' Ohm'] );


plot( freq/1e9, 20*log10(abs(s11)), 'k-', 'Linewidth', 2 );
ylim([-60 0]);
grid on
title( 'reflection coefficient S_{11}' );
xlabel( 'frequency f / GHz' );
ylabel( 'reflection coefficient |S_{11}|' );

drawnow

%break %!!!!!!

%% NFFF contour plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% calculate the far field at phi=0 degrees and at phi=90 degrees
thetaRange = (0:2:359) - 180;
disp( 'calculating far field at phi=[0 90] deg...' );
nf2ff = CalcNF2FF(nf2ff, [Sim_Path Sim_Dir], f0, thetaRange*pi/180, [0 90]*pi/180);

Dlog=10*log10(nf2ff.Dmax);
G_a = 4*pi*A/(c0/f0)^2;
e_a = nf2ff.Dmax/G_a;

% display some antenna parameter
disp( ['radiated power: Prad = ' num2str(nf2ff.Prad) ' Watt']);
disp( ['directivity: Dmax = ' num2str(Dlog) ' dBi'] );
disp( ['aperture efficiency: e_a = ' num2str(e_a*100) '%'] );
 

%%
% normalized directivity
figure
plotFFdB(nf2ff,'xaxis','theta','param',[1 2]);
drawnow

% polar plot
figure
polarFF(nf2ff,'xaxis','theta','param',[1 2],'logscale',[-19 21], 'xtics', 10);
drawnow
%   polar( nf2ff.theta, nf2ff.E_norm{1}(:,1) )

%% calculate 3D pattern
phiRange = sort( unique( [-180:5:-100 -100:2.5:-50 -50:1:50 50:2.5:100 100:5:180] ) );
thetaRange = sort( unique([ 0:1:50 50:2.:100 100:5:180 ]));

disp( 'calculating 3D far field...' );
nf2ff = CalcNF2FF(nf2ff, [Sim_Path Sim_Dir], f0, thetaRange*pi/180, phiRange*pi/180, 'Verbose',2,'Outfile','nf2ff_3D.h5');

figure
plotFF3D(nf2ff);

%%
E_far_normalized = nf2ff.E_norm{1}/max(nf2ff.E_norm{1}(:));
DumpFF2VTK([Sim_Path Sim_Dir '/Horn_Pattern.vtk'],E_far_normalized,thetaRange,phiRange,'scale',1e-3);

thorsten
Posts: 1460
Joined: Mon 27 Jun 2011, 12:26

Re: Numeric errors lead to amplification/solution instabilit

Post by thorsten » Fri 15 Jul 2016, 20:49

You used PML_16 at z_max, but put you nf2ff to "mesh.z(end-10)" that is inside a 16 cell wide PML. Thus you killed your field record for the nf2ff...

Please do not use any other than PML_8, it is of no use, fully untested, potentially unstable and even if everything works you gain virtually nothing from it...

And you seem to have the start of your wavguideport at mesh.z(1) too. That means your excitation is inside the PML... That may end nasty...
Maybe better like this?:

Code: Select all

start=[-a/2, -b/2, mesh.z(9) ];
stop =[ a/2,  b/2, mesh.z(9)+probe.indent ];
[CSX, port] = AddRectWaveGuidePort( CSX, 0, 1, start, stop, 2, a*unit, b*unit, TE_mode, 1);
regards
Thorsten

thorsten
Posts: 1460
Joined: Mon 27 Jun 2011, 12:26

Re: Numeric errors lead to amplification/solution instabilit

Post by thorsten » Fri 15 Jul 2016, 21:26

I have been working on a "real" feed for a rect. waveguide.
Please have a look, maybe you can use it to feed your horn. This way you can model air all around.

regards
Thorsten
Attachments
Rect_Waveguide_Feed.m
(3.04 KiB) Downloaded 476 times

Post Reply