# Out of memory while using FFT

10 views (last 30 days)
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));
##### 2 CommentsShowHide 1 older comment
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');

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.