Continue to Site

Welcome to EDAboard.com

Welcome to our site! EDAboard.com is an international Electronics Discussion Forum focused on EDA software, circuits, schematics, books, theory, papers, asic, pld, 8051, DSP, Network, RF, Analog Design, PCB, Service Manuals... and a whole lot more! To participate you need to register. Registration is free. Click here to register now.

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

Status
Not open for further replies.

zozokim

Newbie
Joined
Apr 2, 2023
Messages
1
Helped
0
Reputation
0
Reaction score
0
Trophy points
1
Activity points
41
I've been studying the UPML portion of Computational electrodynamics by ALLEN TAFLOVE and writing matlab code, but the electromagnetic field is not updating.
Please help me

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

%PML condition(polynomial grading)

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))))
    shading flat
    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: 81
Last edited by a moderator:

Status
Not open for further replies.

Part and Inventory Search

Welcome to EDABoard.com

Sponsor

Back
Top