Continue to Site

# My 3D FDTD (UPML) MATLAB code is not updating and I need help.

#### zozokim

##### Newbie
I've been studying the UPML portion of Computational electrodynamics by ALLEN TAFLOVE and writing matlab code, but the electromagnetic field is not updating.

I can't remove things like emoticon :/), so I'm attaching it as a notepad.
Thank you.

I apologize for my poor English.
This is my MATLAB code.
Code:
--------------------------------------------------------------------------------------------------------

clear all; close all; clc;

%%
%Defining constants

e0 = 8.854*10^-12;                      % Permittivity of vacuum [F/m]
u0 = 4*pi*10^-7;                        % Permeability of vacuum [H/m]
c0 = 1/(e0*u0)^0.5;                     % Speed of light         [m/s]
f0 = 10e6;                              % Frequency of Source    [Hz]
w0 = 2*pi*f0;                           % Angular frequency      [rad/s]
l0 = c0/f0;                             % lambda                 [m]

X = 100;                                % X size [m]
Y = 100;                                % Y size [m]
Z = 100;                                % Z size [m]
Mx = 100;                               % X cell
My = 100;                               % Y cell
Mz = 100;                               % Z cell
N = 2000;                                % Time step
dx = X/Mx;                              % Increment of X
dy = Y/My;                              % Increment of Y
dz = Z/Mz;                              % Increment of Z
dt = (dx^-2+dy^-2+dz^-2)^(-0.5)/c0;    % Courant number (Magic time step)
x_axis = (1:Mx)*dx;                     % Physical length of X
y_axis = (1:My)*dy;                     % Physical length of X

sigx = zeros(Mx+1,My+1,Mz+1);                 % Conductivity x-axis
sigy = zeros(Mx+1,My+1,Mz+1);                 % Conductivity y-axis
sigz = zeros(Mx+1,My+1,Mz+1);                 % Conductivity z-axis (In 2D, sigz = 0)
e = e0*ones(Mx+1,My+1,Mz+1);                  % Permittivity matrix of PML region
u = u0*ones(Mx+1,My+1,Mz+1);                  % Permeability matrix of PML region
eta = sqrt(u0/e0);                      % Intrinsic impedance of PML region

p = 25;                                 % PML grid
kx = 1;                                 % Kappa x
ky = 1;                                 % Kappa y
kz = 1;                                 % Kappa z (In 2D, Kappa z = 1)
m = 4;                                  % Grading order
R0 = exp(-16);                          % Reflection errors
sig_max = -((m+1)*log(R0))/(2*eta*p);   % Sigma max

for n = 1:p
sigx(p+1-n,:,:) = (n/p)^m*sig_max;
sigy(:,p+1-n,:) = (n/p)^m*sig_max;
sigz(:,:,p+1-n) = (n/p)^m*sig_max;
sigx(end-p+n,:,:) = (n/p)^m*sig_max;
sigy(:,end-p+n,:) = (n/p)^m*sig_max;
sigz(:,:,end-p+n) = (n/p)^m*sig_max;
end

%Initializing dimensions
Dx = zeros(Mx,My+1,Mz+1);
Ex = zeros(Mx,My+1,Mz+1);
Dy = zeros(Mx+1,My,Mz+1);
Ey = zeros(Mx+1,My,Mz+1);
Dz = zeros(Mx+1,My+1,Mz);
Ez = zeros(Mx+1,My+1,Mz);
Bx = zeros(Mx+1,My,Mz);
Hx = zeros(Mx+1,My,Mz);
By = zeros(Mx,My+1,Mz);
Hy = zeros(Mx,My+1,Mz);
Bz = zeros(Mx,My,Mz+1);
Hz = zeros(Mx,My,Mz+1);

% monitor = zeros(Mx,My,N);

Wx = 2*e0*kx-sigx*dt;                   %Coefficient of update equations
Wy = 2*e0*ky-sigy*dt;                   %Coefficient of update equations
Wz = 2*e0*kz-sigz*dt;                   %Coefficient of update equations
Qx = 2*e0*kx+sigx*dt;                   %Coefficient of update equations
Qy = 2*e0*ky+sigy*dt;                   %Coefficient of update equations
Qz = 2*e0*kz+sigz*dt;                   %Coefficient of update equations

clims = [-70 -30];

% v = VideoWriter('2D FDTD PML1');
% open(v);

%%
%UPML Update equation by auxiliary differential equation

% Source
t = (1:N)*dt;
E_in = exp(-2 * pi^2 * (f0)^2 * (t - 1/(f0)).^2 );
figure(2);
plot(t, E_in);

i1 = 1:Mx;
j1 = 1:My;
k1 = 1:Mz;
i2 = 2:Mx;
j2 = 2:My;
k2 = 2:Mz;

% Define coefficients of update equations
T1 = (2*e0*dt);

Dx_C1 = Wy(i1,j2,k2)./Qy(i1,j2,k2);
Dx_C2 = T1./Qy(i1,j2,k2)/dy;
Dx_C3 = T1./Qy(i1,j2,k2)/dz;
Dy_C1 = Wz(i2,j1,k2)./Qz(i2,j1,k2);
Dy_C2 = T1./Qz(i2,j1,k2)/dz;
Dy_C3 = T1./Qz(i2,j1,k2)/dx;
Dz_C1 = Wx(i2,j2,k1)./Qx(i2,j2,k1);
Dz_C2 = T1./Qx(i2,j2,k1)/dx;
Dz_C3 = T1./Qx(i2,j2,k1)/dy;

Ex_C1 = Wz(i1,j2,k2)./Qz(i1,j2,k2);
Ex_C2 = 1./(Qz(i1,j2,k2).*e(i1,j2,k2)).*Qx(i1,j2,k2);
Ex_C3 = 1./(Qz(i1,j2,k2).*e(i1,j2,k2)).*Wx(i1,j2,k2);
Ey_C1 = Wx(i2,j1,k2)./Qx(i2,j1,k2);
Ey_C2 = 1./(Qx(i2,j1,k2).*e(i2,j1,k2)).*Qy(i2,j1,k2);
Ey_C3 = 1./(Qx(i2,j1,k2).*e(i2,j1,k2)).*Wy(i2,j1,k2);
Ez_C1 = Wy(i2,j2,k1)./Qy(i2,j2,k1);
Ez_C2 = 1./(Qy(i2,j2,k1).*e(i2,j2,k1)).*Qz(i2,j2,k1);
Ez_C3 = 1./(Qy(i2,j2,k1).*e(i2,j2,k1)).*Wz(i2,j2,k1);

Bx_C1 = Wy(:,j1,k1)./Qy(:,j1,k1);
Bx_C2 = T1./Qy(:,j1,k1)/dy;
Bx_C3 = T1./Qy(:,j1,k1)/dz;
By_C1 = Wz(i1,:,k1)./Qz(i1,:,k1);
By_C2 = T1./Qz(i1,:,k1)/dz;
By_C3 = T1./Qz(i1,:,k1)/dx;
Bz_C1 = Wx(i1,j1,:)./Qx(i1,j1,:);
Bz_C2 = T1./Qx(i1,j1,:)/dx;
Bz_C3 = T1./Qx(i1,j1,:)/dy;

Hx_C1 = Wz(:,j1,k1)./Qz(:,j1,k1);
Hx_C2 = 1./(Qz(:,j1,k1).*u(:,j1,k1)).*Qx(:,j1,k1);
Hx_C3 = 1./(Qz(:,j1,k1).*u(:,j1,k1)).*Wx(:,j1,k1);
Hy_C1 = Wx(i1,:,k1)./Qx(i1,:,k1);
Hy_C2 = 1./(Qx(i1,:,k1).*u(i1,:,k1)).*Qy(i1,:,k1);
Hy_C3 = 1./(Qx(i1,:,k1).*u(i1,:,k1)).*Wy(i1,:,k1);
Hz_C1 = Wy(i1,j1,:)./Qy(i1,j1,:);
Hz_C2 = 1./(Qy(i1,j1,:).*u(i1,j1,:)).*Qz(i1,j1,:);
Hz_C3 = 1./(Qy(i1,j1,:).*u(i1,j1,:)).*Wz(i1,j1,:);

figure(1);
for n = 1:N
t(n)

Dx_temp = Dx(:,j2,k2);
Dy_temp = Dy(i2,:,k2);
Dz_temp = Dz(i2,j2,:);

Ex_temp = Ex(:,j2,k2);
Ey_temp = Ey(i2,:,k2);
Ez_temp = Ez(i2,j2,:);

Bx_temp = Bx(:,j1,k1);
By_temp = By(i1,:,k1);
Bz_temp = Bz(i1,j1,:);

Hx_temp = Hx(:,j1,k1);
Hy_temp = Hy(i1,:,k1);
Hz_temp = Hz(i1,j1,:);

Dx(:,j2,k2) = Dx_C1.*Dx_temp + Dx_C2.*(Hz(:,j2,k2) - Hz(:,j2-1,k2)) - Dx_C3.*(Hy(:,j2,k2) - Hy(:,j2,k2-1));
Dy(i2,:,k2) = Dy_C1.*Dy_temp + Dy_C2.*(Hx(i2,:,k2) - Hx(i2,:,k2-1)) - Dy_C3.*(Hz(i2,:,k2) - Hz(i2-1,:,k2));
Dz(i2,j2,:) = Dz_C1.*Dz_temp + Dz_C2.*(Hy(i2,j2,:) - Hy(i2-1,j2,:)) - Dz_C3.*(Hx(i2,j2,:) - Hx(i2,j2-1,:));

Ex(:,j2,k2) = Ex_C1.*Ex_temp + Ex_C2.*Dx(:,j2,k2) - Ex_C3.*Dx_temp;
Ey(i2,:,k2) = Ey_C1.*Ey_temp + Ey_C2.*Dy(i2,:,k2) - Ey_C3.*Dy_temp;
Ez(i2,j2,:) = Ez_C1.*Ez_temp + Ez_C2.*Dz(i2,j2,:) - Ez_C3.*Dz_temp;

Bx(:,j1,k1) = Bx_C1.*Bx_temp - Bx_C2.*(Ez(:,j1+1,k1) - Ez(:,j1,k1)) + Bx_C3.*(Ey(:,j1,k1+1) - Ey(:,j1,k1));
By(i1,:,k1) = By_C1.*By_temp - By_C2.*(Ex(i1,:,k1+1) - Ex(i1,:,k1)) + By_C3.*(Ez(i1+1,:,k1) - Ez(i1,:,k1));
Bz(i1,j1,:) = Bz_C1.*Bz_temp - Bz_C2.*(Ey(i1+1,j1,:) - Ey(i1,j1,:)) + Bz_C3.*(Ex(i1,j1+1,:) - Ex(i1,j1,:));

Hx(:,j1,k1) = Hx_C1.*Hx_temp + Hx_C2.*Bx(:,j1,k1) - Hx_C3.*Bx_temp;
Hy(i1,:,k1) = Hy_C1.*Hy_temp + Hy_C2.*By(i1,:,k1) - Hy_C3.*By_temp;
Hz(i1,j1,:) = Hz_C1.*Hz_temp + Hz_C2.*Bz(i1,j1,:) - Hz_C3.*Bz_temp;

Ez(round(X/2),round(Y/2),round(Z/2)) = Ez(round(X/2),round(Y/2),round(Z/2)) - exp(-2 * pi^2 * (f0)^2 * (n*dt - 1/(f0)).^2 );

E_total = sqrt(abs(Ex(i1,j1,k1)).^2+abs(Ey(i1,j1,k1)).^2+abs(Ez(i1,j1,k1)).^2);
%     subplot(1,2,1)
%pcolor(db(E_total(:,:,round(Mz/2))))
figure(1);
pcolor(y_axis, x_axis, (E_total(:,:,round(Z/2))))
caxis([0 5]);
colorbar;

drawnow;
%     frame = getframe(gcf);
%     writeVideo(v,frame);
end

% close(v);

---------------------------------------------------------------------------------------------------------------------------

#### Attachments

• clear all; close all; clc;.txt
6.9 KB · Views: 27
Last edited by a moderator: