MATLAB Answers

Out of memory while using FFT

12 views (last 30 days)
Sreenidhi Yeturi
Sreenidhi Yeturi on 18 Apr 2021
Answered: Madhav Thakker on 18 May 2021
I'm trying to plot the field for a FDTD simulation in frequncy domain but I keep running into memory error. I've attached my code below.
Nfft = nt *8;
df = 1/(Nfft*dt);
ftEZ = abs(fft(Ez_out,Nfft));
fmax = 600e6;
numfSamps = floor(fmax/df);
fr = df*(0:numfSamps-1);
figure;
plot(fr,ftEz(1:numfSamps));
ylabel('Amplitude (V/m)');
xlabel('Frequency (Hz)');
set(gca,'FontSize',16);
hold on;
p = ftEz(1:numfSamps);
[peaks,locs] = findpeaks(p,'MinPeakDistance',width, 'MinPeakHeight',height);
f110 = 2.998e8/2/pi*sqrt(2*pi^2);
f210 = 2.998e8/2/pi*sqrt((2*pi^2)+pi^2);
nf110 = floor(f110/df);
nf210 = floor(f210/df);
fex = [f110,f210];
aex = [ftEz(nf110), ftEz(nf210)];
plot(fex,aex,'ro');
legend('FDTD sim','Exact','location','northwest');
%%%%%Error Message%%%%%
Error using fft
Requested array exceeds the maximum possible variable size.
Error in Trail_code (line 179)
ftEZ = abs(fft(Ez_out,Nfft));
More information
  2 Comments
Sreenidhi Yeturi
Sreenidhi Yeturi on 19 Apr 2021
Thanks for commenting Matt, here's the complete code.
scalefactor=11; % This scales the resolution,
c = 2.99792458e+8;
no = 376.7303134617706554679;
mu = no/c;
e0 = 1/(no*c);
epsr = 4.2;
cfln = 0.995; % Stability criterion
tw = 1e-9; % half width of pulse 10 had low resonance
to = 5*tw; % Time delay for sources
Lx = 1;
Ly = 1;
Lz = 1;
nx = 11*scalefactor;
ny = 11*scalefactor;
nz = 11*scalefactor;
dx = Lx/nx;
dy = Ly/ny;
dz = Lz/nz;
S=5.8e1;
eps = e0;
dt = cfln/(c*sqrt((1/dx)^2+(1/dy)^2+(1/dz)^2));
T = 4.5e-9*2;
simtime = T/dt;
nt = floor(simtime / dt);
nMax = 50; % Time steps
% Injection and Output point ranges
% Source Injection points
i1_src = 4;
j1_src = 4;
k1_src = 4;
i2_src = 4;
j2_src = 4;
k2_src = 6;
% Field ouput points
i1_fld = 8;
j1_fld = 8;
k1_fld = 8;
i2_fld = 8;
j2_fld = 8;
k2_fld = 9;
% Electric & Magnetic Field Coefficients
% Electric Field
cExy = dt/(dy*eps);
cExz = dt/(dz*eps);
cEyz = dt/(dz*eps);
cEyx = dt/(dx*eps);
cEzx = dt/(dx*eps);
cEzy = dt/(dy*eps);
% Magnetic Field
cHxy = dt/(mu*dy);
cHxz = dt/(mu*dz);
cHyz = dt/(mu*dz);
cHyx = dt/(mu*dx);
cHzx = dt/(mu*dx);
cHzy = dt/(mu*dy);
% PEC Boundary conditoins
% Set all tangential E fields to zero and never update them. This will keep
% them zero. Update loops only modify inner cells
Hx=zeros(nx,ny-1,nz-1);
Hy=zeros(nx-1,ny,nz-1);
Hz=zeros(nx-1,ny-1,nz);
Ex=zeros(nx-1,ny,nz);
Ey=zeros(nx,ny-1,nz);
Ez=zeros(nx,ny,nz-1);
% Ez_out = zeros(nt,1);
for n=1:nMax
%Main update loops for the Electric fields:
%Ex
for k = 2:nz-1
for j = 2:ny-1
for i = 1:nx-1 %1->2
Ex(i,j,k) = Ex(i,j,k)+cExy*((Hz(i,j,k)-Hz(i,j-1,k))-cExz*(Hy(i,j,k)-Hy(i,j,k-1)));
end
end
end
%Ey
for k = 2:nz-1
for j = 1:ny-1
for i = 2:nx-1
Ey(i,j,k) = Ey(i,j,k)+cEyz*((Hx(i,j,k)-Hx(i,j,k-1))-cEyx*(Hz(i,j,k)-Hz(i-1,j,k)));
end
end
end
%Ez
for k = 1:nz-1 %1->2
for j = 2:ny-1
for i = 2:nx-1
Ez(i,j,k) = Ez(i,j,k)+cEzx*((Hy(i,j,k)-Hy(i-1,j,k))-cEzy*(Hx(i,j,k)-Hx(i,j-1,k)));
end
end
end
%Source Injection
t = (n-0.5)*dt;
Jz =(-2/tw)*(t-to)*exp(-(t-to)^2/tw^2);
for k = k1_src:k2_src-1
for j = j1_src:j2_src
for i = i1_src:i2_src
Ez(i,j,k) = Ez(i,j,k)+dt*Jz;
end
end
end
%Main update loops for the Magnetic fields:
%Hx
for k = 1:nz-1 %1->2 all 3
for j = 1:ny-1
for i = 1:nx %-1
Hx(i,j,k) = Hx(i,j,k)-cHxy*(Ez(i,j+1,k)-Ez(i,j,k))+cHxz*(Ey(i,j,k+1)-Ey(i,j,k));
end
end
end
%Hy
for k = 1:nz-1 %1->2 all 3
for j = 1:ny %-1
for i = 1:nx-1
Hy(i,j,k) = Hy(i,j,k)-cHyz*(Ex(i,j,k+1)-Ex(i,j,k))+cHyx*(Ez(i+1,j,k)-Ez(i,j,k));
end
end
end
%Hz
for k = 1:nz %-1 %1->2 all 3
for j = 1:ny-1
for i = 1:nx-1
Hz(i,j,k) = Hz(i,j,k)-cHzx*(Ey(i+1,j,k)-Ey(i,j,k))+cHzy*(Ex(i,j+1,k)-Ex(i,j,k));
end
end
end
%%%%%%%%%% Output %%%%%%%%%%%%
Ez_out = zeros(n);
for k = k1_fld:k2_fld-1
for j = j1_fld:j2_fld
for i = i1_fld:i2_fld
Ez_out(n) = Ez_out(n) + Ez(i,j,k);
end
end
end
clc
round(n/nMax*100) %gives progress update percentage on home screen
%surf(squeeze(Ez(:,:,0.1*scalefactor))); %Z-slice with resonator in center
surf(squeeze(Ex(1*scalefactor,:,:)));
%surf(squeeze(Ey(:,0.1*scalefactor,:)));
set(gcf,'renderer','opengl')
%hold on
%whitebg('black');
grid off
set(gcf,'Position',[50 50 1500 800]);
axis([-1*scalefactor 5*scalefactor -2*scalefactor 7*scalefactor -1/scalefactor 1/scalefactor])
view([3,4,5])
Mov(n)=getframe; % comment this out for nonanimated analysis
end
movie(Mov,30,30); % comment this out for nonanimated analysis
%Post Process
Nfft = nt *8;
df = 1/(Nfft*dt);
ftEZ = abs(fft(Ez_out,Nfft));
fmax = 600e6;
numfSamps = floor(fmax/df);
fr = df*(0:numfSamps-1);
figure;
plot(fr,ftEz(1:numfSamps));
ylabel('Amplitude (V/m)');
xlabel('Frequency (Hz)');
set(gca,'FontSize',16);
hold on;
p = ftEz(1:numfSamps);
[peaks,locs] = findpeaks(p,'MinPeakDistance',width, 'MinPeakHeight',height);
f110 = 2.998e8/2/pi*sqrt(2*pi^2);
f210 = 2.998e8/2/pi*sqrt((2*pi^2)+pi^2);
nf110 = floor(f110/df);
nf210 = floor(f210/df);
fex = [f110,f210];
aex = [ftEz(nf110), ftEz(nf210)];
plot(fex,aex,'ro');
legend('FDTD sim','Exact','location','northwest');

Sign in to comment.

Answers (1)

Madhav Thakker
Madhav Thakker on 18 May 2021
Hi Sreenidhi,
The variable nt has value 3.5886e+13 which is due to which the resulting Fourier transform array exceeds the maxiumum variable size.
Hope this helps.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by