When I try to reduce the mesh, system reboot happens.

The same time small antenna near2far field problems compute well.

So, it seems, there's kind of memory allocation bug.

openEMS 64bit v0.0.33

Octave 4.0.2 MINGW32_NT-6.1

Windows 7 SP1 x64

i5-4590 & 8Gb RAM

Dell Optiplex 7020.

Code: Select all

```
confirm_recursive_rmdir (0,"local")
close all
clear
clc
%% setup the simulation
physical_constants;
unit = 1e-3; % all length in mm
% setup simulation quality
wave_subsample = 4; %resolution denominator - part of lambda. recommended in the range of 10-20 !!!! Here 10 crashes, 5 -reboots !!!
stop_dB = -30; % stop calculation when reached wave decay -xx dB
freespace = 1 ; %added free space X lambda/4
% frequency range of interest
f_start = 10.7e9;
f_stop = 11.3e9;
f0 = 0.5*(f_start+f_stop); % frequency of interest
lambda = c0 / (f0)/ unit; %wavelength
%waveguide TE-mode definition
%TE_mode = 'TE01'; %not needed with lumped port
% horn width in x-direction
a = 22.86; %22.86; %a
% horn height in y-direction
b = 10.16; %10.16; %b
horn.width = 118;
horn.height = 92;
horn.edge = 127+6.5;
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 = 73.1;
%horn.length = sqrt((2*a/0.7)^2-(a-b/2)^2 -(2.08*a/1.4 -a/2)^2) ;
horn.length = sqrt( (horn.edge)^2-(horn.height-b)^2 /4-(horn.width-a)^2 /4 );
horn.feed_length = 22.55+20.5; %lambda/4+4.6/((2*pi/lambda)*sqrt((lambda/a)^2-1)); %36.5;
probe.indent = 7.05; %lambda/4
probe.height = 0.73*lambda/4; %lambda/4
horn.thickness = 5;
%SMA
coax.Dout=4.53;
coax.Din=1.24;
coax.epsilon=2.02;
% horn opening angle in x, y
%horn.angle = [24.8 23.5]*pi/180; % a-expansion b-expansion angles
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
% size of the simulation box
%SimBox = [freespace*lambda/2+a+2.*horn.length*tan(horn.angle(1))+horn.thickness.*2, freespace*lambda/2+b+2.*horn.length*tan(horn.angle(2))+horn.thickness.*2, ((floor(freespace/4)*8+1)*lambda/4)+horn.feed_length+horn.length+horn.thickness*2];
SimBox = [freespace*lambda/2+a+2.*horn.length*tan(horn.angle(1))+horn.thickness.*2, freespace*lambda/2+b+2.*horn.length*tan(horn.angle(2))+horn.thickness.*2, freespace*lambda/4+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_4' '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
CSX = InitCSX();
%create fixed lines for the simulation box, substrate and port
mesh.x = [-1600, -coax.Dout/2-horn.thickness, -coax.Dout/2, -coax.Din/2, 0, coax.Din/2, coax.Dout/2, coax.Dout/2+horn.thickness, a/2, SimBox(1)/2,1800];
mesh.x = SmoothMeshLines( mesh.x, max_res, 1.4); % create a smooth mesh between specified fixed mesh lines
mesh.y = [-1300,-SimBox(2)/2,-b/2-horn.thickness-1,-b/2-horn.thickness,-b/2,-b/2+.1, -b/2+probe.height, b/2, SimBox(2)/2,800];
mesh.y = SmoothMeshLines( mesh.y, max_res, 1.4 );
%create fixed lines for the simulation box and given number of lines inside the substrate
mesh.z = [-horn.thickness-horn.feed_length,-horn.feed_length, -horn.feed_length+probe.indent-coax.Dout/2-horn.thickness, -horn.feed_length+probe.indent-coax.Dout/2, -horn.feed_length+probe.indent-coax.Din/2, -horn.feed_length+probe.indent, -horn.feed_length+probe.indent+coax.Din/2, -horn.feed_length+probe.indent+coax.Dout/2, -horn.feed_length+probe.indent+coax.Dout/2+horn.thickness, 0, horn.length, SimBox(3)-horn.feed_length-horn.thickness,horn.length, SimBox(3)-horn.feed_length-horn.thickness+4050 ];
mesh.z = SmoothMeshLines( mesh.z, max_res, 1.4 );
CSX = DefineRectGrid( CSX, unit, mesh );
%% create horn
% horn feed rect waveguide
CSX = AddMetal(CSX, 'metal');
%CSX = AddConductingSheet(CSX,'metal',56e6,70e-6); %copper
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);
start = [-a/2-horn.thickness, -b/2-horn.thickness, -horn.feed_length];
stop = [ a/2+horn.thickness, b/2+horn.thickness, -horn.feed_length-horn.thickness];
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']});
%Feeder
%% drill hole for SMA probe to enter
CSX = AddMaterial( CSX, 'teflon' );
CSX = SetMaterialProperty( CSX, 'teflon', 'Epsilon', coax.epsilon);
start = [0, -b/2-horn.thickness-1, -horn.feed_length+probe.indent];
stop = [0, -b/2+.1, -horn.feed_length+probe.indent];
CSX = AddCylinder(CSX,'teflon' ,20,start, stop,coax.Dout/2);
%Probe
start = [0, -b/2-horn.thickness-1, -horn.feed_length+probe.indent];
stop = [0, -b/2+probe.height, -horn.feed_length+probe.indent];
CSX = AddCylinder(CSX, 'metal', 30, start, stop, coax.Din/2);
%adapter
start = [0, -b/2-horn.thickness-1, -horn.feed_length+probe.indent];
stop = [0, -b/2, -horn.feed_length+probe.indent];
CSX = AddCylindricalShell(CSX, 'metal', 10, start, stop, coax.Dout/2+horn.thickness/2, horn.thickness);
%% no need for AddBox because teflon shape is set in AddCoaxialPort
start = [0, mesh.y(1), -horn.feed_length+probe.indent];
stop = [0, -b/2-horn.thickness-1, -horn.feed_length+probe.indent];
[CSX,port] = AddCoaxialPort( CSX, 40, 1, 'metal', 'teflon', start, stop, 'y',coax.Din/2 , coax.Dout/2, coax.Dout/2+horn.thickness,'ExciteAmp', 1,'FeedShift', (-b/2-horn.thickness-mesh.y(1))/4 );
%Frame
CSX = AddMaterial(CSX,'air');
CSX = SetMaterialProperty(CSX,'air','Epsilon',1.0,'Mue',1.0);
start = [-720/2, -1100, mesh.z(1)];
stop = [ 720/2, -b/2-horn.thickness, 0];
CSX = AddBox(CSX,'metal',10,start,stop);
start = [-720/2+50, -1100+50, mesh.z(1)];
stop = [ 720/2-50, -b/2-horn.thickness-50, 0];
CSX = AddBox(CSX,'air',20,start,stop);
%FLOOR
CSX = AddMaterial(CSX,'concrete');
CSX = SetMaterialProperty(CSX,'concrete','Epsilon',4.5,'Mue',1.0,'Sigma',0.7);
start = [-1600, -1100, -horn.thickness-horn.feed_length,-horn.feed_length];
stop = [ 1800, -1100-200, SimBox(3)-horn.feed_length-horn.thickness+4050 ];
CSX = AddBox(CSX,'concrete',10,start,stop);
%stand pile
start = [-600/2, -1100, 1000];
stop = [ 600/2, -1100+50, 2350];
CSX = AddBox(CSX,'metal',10,start,stop);
start = [-600/2+50, -1100, 1050];
stop = [ 600/2-50, -1100+50, 2300];
CSX = AddBox(CSX,'air',20,start,stop);
start = [-600/2, -1050, 1100];
stop = [ 600/2, -1050+50, 2220];
CSX = AddBox(CSX,'metal',10,start,stop);
start = [-600/2+50, -1050, 1150];
stop = [ 600/2-50, -1050+50, 2170];
CSX = AddBox(CSX,'air',20,start,stop);
start = [-320/2, -1000, 1200];
stop = [ 320/2, -1000+50, 2130];
CSX = AddBox(CSX,'metal',10,start,stop);
start = [-270/2+50, -1000, 1250];
stop = [ 270/2-50, -1000+50, 2080];
CSX = AddBox(CSX,'air',20,start,stop);
start = [-2.5, -950, 1270];
stop = [ 2.5, 250, 1270+600];
CSX = AddBox(CSX,'metal',10,start,stop);
%WALLS
start = [1800, -1100, -horn.thickness-horn.feed_length,-horn.feed_length];
stop = [ 1750, 150, SimBox(3)-horn.feed_length-horn.thickness+4050 ];
CSX = AddBox(CSX,'metal',10,start,stop);
start = [-1600, -1100, 700];
stop = [ -1550, 150, SimBox(3)-horn.feed_length-horn.thickness+4050 ];
CSX = AddBox(CSX,'metal',10,start,stop);
start = [-1600, -1100, SimBox(3)-horn.feed_length-horn.thickness+4000];
stop = [1800, 150, SimBox(3)-horn.feed_length-horn.thickness+4050 ];
CSX = AddBox(CSX,'metal',10,start,stop);
% 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+lambda/4 ];
%stop =[ a/2, b/2, -horn.feed_length+1.3*lambda/4 ];
%[CSX, port] = AddRectWaveGuidePort( CSX, 0, 1, start, stop, 2, a*unit, b*unit, TE_mode, 1);
%Lumped port
%start = [0, -b/2, -horn.feed_length+probe.indent];
%stop = [0, -b/2+probe.height, -horn.feed_length+probe.indent];
%[CSX, port] = AddLumpedPort(CSX, 0 ,1 , 50, start, stop, [0 1 0], true);
%% nf2ff calc
start = [mesh.x(5) mesh.y(5) mesh.z(5)];
stop = [mesh.x(end-4) mesh.y(end-4) mesh.z(end-4)];
[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop, 'Directions', [1 1 1 1 0 1]);
%% prepare simulation folder
Sim_Path = 'tmp_Environment';
Sim_CSX = 'Environment.xml';
[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory
[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder
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_CSX], FDTD, CSX);
%% show the structure
CSXGeomPlot([Sim_Path '/' Sim_CSX],"--RenderDiscMaterial");
%% run openEMS
RunOpenEMS(Sim_Path, Sim_CSX);
%% postprocessing & do the plots
freq = linspace(f_start,f_stop,201);
port = calcPort(port, Sim_Path, 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
%% 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, 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
% D_log = 20*log10(nf2ff.E_norm{1}/max(max(nf2ff.E_norm{1})));
% D_log = D_log + 10*log10(nf2ff.Dmax);
% plot( nf2ff.theta, D_log(:,1) ,'k-', nf2ff.theta, D_log(:,2) ,'r-' );
% polar plot
figure
polarFF(nf2ff,'xaxis','theta','param',[1 2],'logscale',[-15 15], '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, 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 '/Pattern.vtk'],E_far_normalized,thetaRange,phiRange,'scale',1e-3);
```