Out of memory while using FFT

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

[Matt J]

We don't know what the value of nt is, but isn't the cause of the problem obvious from the erorr message? Your Nfft must be absurdly big.

[Sreenidhi Yeturi]

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');

Answer 1

[Madhav Thakker]

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.