pml2d.m

二维FDTD程序

function pml2Dexample()
%%% 仅分裂PML中的场。
% PML test  
%% the problem:在100*50的区域中，在点（50，25）处加H激励H,运行100 time steps 用pml和mur分别计算j=1处的Hz(i,1).
               %计算空间扩大到400*400，不加边界，运行相同的时间步，在相同的地方取Hr，
               %局部误差：(Hz-Hr)/Hr(50,1)
nx=100;
Hz_mur=zeros(nx,1);R_mur=zeros(nx,1);
Hz_pml0=zeros(nx,1);R_pml0=zeros(nx,1);
Hz_pml1=zeros(nx,1);R_pml1=zeros(nx,1);
Hz_pml2=zeros(nx,1);R_pml2=zeros(nx,1);
Hz_non=zeros(nx,1);
%%calculate
nx=400;ny=400;
Hz_non=local_maincal(nx,ny);% calculate the hz without reflection bigger calculation domain
nx=100;ny=50;
npml=0;
Hz_mur=local_maincal(nx,ny,npml);
R_mur=(Hz_mur-Hz_non)/Hz_non(50,1);
R0=1e-3;%the theoretical reflection coefficient at normal incidence for the PML over PEC
npml=4;
order=0;%
Hz_pml40=local_maincal(nx,ny,npml,order,R0);%nargin==0 mur//nargin==2 pml
R_pml40=(Hz_pml40-Hz_non)/ Hz_non(50,1);
R0=1e-3;%the theoretical reflection coefficient at normal incidence for the PML over PEC
npml=4;
order=1;%
Hz_pml41=local_maincal(nx,ny,npml,order,R0);%nargin==0 mur//nargin==2 pml
R_pml41=(Hz_pml41-Hz_non)/ Hz_non(50,1);
R0=1e-3;%the theoretical reflection coefficient at normal incidence for the PML over PEC
npml=4;
order=2;%
Hz_pml42=local_maincal(nx,ny,npml,order,R0);%nargin==0 mur//nargin==2 pml
R_pml42=(Hz_pml42-Hz_non)/ Hz_non(50,1);
R0=1e-3;%the theoretical reflection coefficient at normal incidence for the PML over PEC
npml=2;
order=1;%
Hz_pml32=local_maincal(nx,ny,npml,order,R0);%nargin==0 mur//nargin==2 pml
R_pml32=(Hz_pml32-Hz_non)/ Hz_non(50,1);
%%plot the result
figure;plot(1:100,R_mur,'k',1:100,R_pml40,'b',1:100,R_pml41,'r',1:100,R_pml42,'g',1:100,R_pml32,'b--')
xlabel('the mesh number')
ylabel('the local error')
legend('1 mur-1','2 PML(4,C,0.001)','3 PML(4,L,0.001)','4 PML(4,P,0.001)','5,PML(3,P,0.001)')

% nx=400;ny=400;
% Hz_non=local_maincal(nx,ny);% calculate the hz without reflection bigger calculation domain
% nx=100;ny=50;
% R0=1e-3;%the theoretical reflection coefficient at normal incidence for the PML over PEC
% npml=4;
% order=2;%
% Hz_pml41=local_maincal(nx,ny,npml,order,R0);%nargin==0 mur//nargin==2 pml
% R_pml41=(Hz_pml41-Hz_non)/ Hz_non(50,1);
% R0=1e-3;%the theoretical reflection coefficient at normal incidence for the PML over PEC
% npml=6;
% order=2;%
% Hz_pml61=local_maincal(nx,ny,npml,order,R0);%nargin==0 mur//nargin==2 pml
% R_pml61=(Hz_pml61-Hz_non)/ Hz_non(50,1);
% R0=1e-3;%the theoretical reflection coefficient at normal incidence for the PML over PEC
% npml=8;
% order=2;%
% Hz_pml81=local_maincal(nx,ny,npml,order,R0);%nargin==0 mur//nargin==2 pml
% R_pml81=(Hz_pml81-Hz_non)/ Hz_non(50,1);
% R0=1e-3;%the theoretical reflection coefficient at normal incidence for the PML over PEC
% npml=10;
% order=2;%
% Hz_pml101=local_maincal(nx,ny,npml,order,R0);%nargin==0 mur//nargin==2 pml
% R_pml101=(Hz_pml101-Hz_non)/ Hz_non(50,1);
% %%%plot the result
% figure;plot(1:100,R_pml41,'k',1:100,R_pml61,'b',1:100,R_pml81,'r',1:100,R_pml101,'g')
% xlabel('the mesh number')
% ylabel('the local error')
% legend('1 PML(4,L,0.001)','2 PML(6,L,0.001)','3 PML(8,L,0.001)','4 PML(10,L,0.001)')


function Hz=local_maincal(nx,ny,npml,order,R0)
%%% premeters for PML
max_iteration=100;%总时间步
plot_modulus=30;%每时间步存储场量
%%%物理常数
light_speed=299792458;%真空光速，单位m
eta_0=120*pi;
mu_0=1.2566e-6;%真空磁导率
epsilon_0=8.854e-12;%真空介电常数
dx=1.5e-2;dy=1.5e-2;
dt=25e-12;
dtmudx=dt/(dx*mu_0);dtmudy=dt/(dy*mu_0);%factor in the main region/vacuum
dtepsdx=dt/(dx*epsilon_0);dtepsdy=dt/(dy*epsilon_0);
if nargin==2
    npml=0;
end
nx=nx+2*npml;ny=ny+2*npml;
ex=zeros(nx,ny+1);ey=zeros(nx+1,ny);hz=zeros(nx,ny);
if nargin==5%pml
    hzx=zeros(nx,ny);hzy=zeros(nx,ny);
    sigmax=zeros(nx-1,ny);sigmay=zeros(nx,ny-1);sigmamx=zeros(nx,ny);sigmamy=zeros(nx,ny);
    
%     maxsigmax=-log(R0)*(order+1)*epsilon_0*light_speed/(2*dx); 
%     maxsigmay=-log(R0)*(order+1)*epsilon_0*light_speed/(2*dy);
%     maxsigmamx=mu_0*maxsigmax/epsilon_0;
%     maxsigmamy=mu_0*maxsigmay/epsilon_0;
%     index=1:npml;
%     sigmax(npml:-1:1,:)=(maxsigmax*((index'-1)/(npml)).^order)*ones(1,ny);
%     sigmax(nx-npml:nx-1,:)=(maxsigmax*((index'-1)/(npml)).^order)*ones(1,ny);
%     sigmamx(npml:-1:1,:)=(maxsigmamx*((index'-0.5)/(npml)).^order)*ones(1,ny);
%     sigmamx(nx-npml+1:nx,:)=(maxsigmamx*((index'-0.5)/(npml)).^order)*ones(1,ny);
%     sigmay(:,npml:-1:1)=ones(nx,1)*(maxsigmay*((index-1)/(npml)).^order);
%     sigmay(:,ny-npml:ny-1)=ones(nx,1)*(maxsigmay*((index-1)/(npml)).^order);
%     sigmamy(:,npml:-1:1)=ones(nx,1)*(maxsigmamy*((index-0.5)/(npml)).^order);
%     sigmamy(:,ny-npml+1:ny)=ones(nx,1)*(maxsigmamy*((index-0.5)/(npml)).^order);
    maxsigmax=-log(R0)*(order+1)*epsilon_0*light_speed/(2*dx);
    maxsigmay=-log(R0)*(order+1)*epsilon_0*light_speed/(2*dy);
    maxsigmamx=mu_0*maxsigmax/epsilon_0;
    maxsigmamy=mu_0*maxsigmay/epsilon_0;
    index=0.5:npml-0.5;
    sigmamx(npml:-1:1,:)=maxsigmamx/(dx*(order+1)*(dx*npml)^order)*((index'*dx+0.5*dx).^(order+1)-(index'*dx-0.5*dx).^(order+1))*ones(1,ny);
    sigmamx(nx-npml+1:nx,:)=maxsigmamx/(dx*(order+1)*(dx*npml)^order)*((index'*dx+0.5*dx).^(order+1)-(index'*dx'-0.5*dx).^(order+1))*ones(1,ny);
    sigmamy(:,npml:-1:1)=ones(nx,1)*maxsigmamy/(dx*(order+1)*(dy*npml)^order)*((index*dy+0.5*dy).^(order+1)-(index*dy-0.5*dy).^(order+1));
    sigmamy(:,ny-npml+1:ny)=ones(nx,1)*maxsigmamy/(dx*(order+1)*(dy*npml)^order)*((index*dy+0.5*dy).^(order+1)-(index*dy-0.5*dy).^(order+1));
    index=1:npml-1;
    sigmax(npml-1:-1:1,:)=maxsigmax/(dx*(order+1)*(dx*npml)^order)*((index'*dx+0.5*dx).^(order+1)-(index'*dx-0.5*dx).^(order+1))*ones(1,ny);
    sigmax(nx-npml+1:nx-1,:)=maxsigmax/(dx*(order+1)*(dx*npml)^order)*((index'*dx+0.5*dx).^(order+1)-(index'*dx-0.5*dx).^(order+1))*ones(1,ny);
    sigmay(:,npml-1:-1:1)=ones(nx,1)*maxsigmay/(dx*(order+1)*(dy*npml)^order)*((index*dy+0.5*dy).^(order+1)-(index*dy-0.5*dy).^(order+1));
    sigmay(:,ny-npml+1:ny-1)=ones(nx,1)*maxsigmay/(dx*(order+1)*(dy*npml)^order)*((index*dy+0.5*dy).^(order+1)-(index*dy-0.5*dy).^(order+1));
    index=0;
    sigmax(npml,:)=maxsigmax/(dx*(order+1)*(dx*npml)^order)*(index*dx+0.5*dx)^(order+1)*ones(1,ny);
    sigmax(nx-npml,:)=maxsigmax/(dx*(order+1)*(dx*npml)^order)*(index*dx+0.5*dx)^(order+1)*ones(1,ny);
    sigmay(:,npml)=ones(nx,1)*maxsigmay/(dx*(order+1)*(dy*npml)^order)*(index*dy+0.5*dy)^(order+1);
    sigmay(:,ny-npml)=ones(nx,1)*maxsigmay/(dx*(order+1)*(dy*npml)^order)*(index*dy+0.5*dy)^(order+1);
end 
tempexy=zeros(nx,1);tempexy_=zeros(nx,1);%存储前一个时间步的临近边界磁场值tempexy(:,:)=ex(:,:,No-1)
tempeyx=zeros(1,ny);tempeyx_=zeros(1,ny);

for iteration=1:max_iteration
    iteration
    current_simulatedtime=(iteration-1)*dt;
    if current_simulatedtime<1e-9;
        stimulus=1/320*(10-15*cos(2*pi*current_simulatedtime*1e9)+6*cos(4*pi*current_simulatedtime*1e9)-cos(6*pi*current_simulatedtime*1e9));
    else
        stimulus=0;
    end
    
    %%%%%%%%%%%%%%%%%%%%update the hz
    %主体区域
    m=1+npml:nx-npml;n=1+npml:ny-npml;
    hz(m,n)=hz(m,n)-dtmudx*(ey(m+1,n)-ey(m,n))+dtmudy*(ex(m,n+1)-ex(m,n));
    m=nx/2;n=ny/2;%the excitation
    hz(m,n)=stimulus;
    % the PML boundary
    if nargin==5
        m=[1:npml,nx-npml+1:nx];n=1:ny;% right and left with the corner
        hzx(m,n)=exp(-sigmamx(m,n)*dt/mu_0).*hzx(m,n)-((1-exp(-sigmamx(m,n)*dt/mu_0))./(dx*sigmamx(m,n))).*(ey(m+1,n)-ey(m,n));
        m=npml+1:nx-npml;n=[1:npml,ny-npml+1:ny];%up and down without the corner
        hzx(m,n)=hzx(m,n)-dtmudx*(ey(m+1,n)-ey(m,n));
        
        m=1:nx;n=[1:npml,ny-npml+1:ny];%up and down with the corner
        hzy(m,n)=exp(-sigmamy(m,n)*dt/mu_0).*hzy(m,n)+((1-exp(-sigmamy(m,n)*dt/mu_0))./(dy*sigmamy(m,n))).*(ex(m,n+1)-ex(m,n));
        hz(m,n)=hzx(m,n)+hzy(m,n);
        m=[1:npml,nx-npml+1:nx];n=npml+1:ny-npml;%right and left without the corner
        hzy(m,n)=hzy(m,n)+dtmudy*(ex(m,n+1)-ex(m,n));
        hz(m,n)=hzx(m,n)+hzy(m,n);
    end
    %%%%%%%%%%%%%%%%%%%%%%update the ex
    ny=ny+1;
    m=1+npml:nx-npml;n=1+(npml+1):ny-(npml+1);%the main area
    ex(m,n)=ex(m,n)+dtepsdy*(hz(m,n)-hz(m,n-1));
    if nargin==3
        %处理边界 mur
        n=1;%y=ymin
            m=1:nx;
            ex(m,n)=tempexy_(m,1)+(light_speed*dt-dy)/(light_speed*dt+dy)*(ex(m,n+1)-ex(m,n));
        n=ny;%y=ymax
            ex(m,n)=tempexy(m,1)+(light_speed*dt-dy)/(light_speed*dt+dy)*(ex(m,n-1)-ex(m,n));%空气
        %记录临时变量
        tempexy_=ex(:,2);
        tempexy=ex(:,ny-1);
    end
    if nargin==5
        m=1:nx;n=[2:npml,ny-npml+1:ny-1];%down with the corner
        ex(m,n)=exp(-sigmay(m,n-1)*dt/epsilon_0).*ex(m,n)+((1-exp(-sigmay(m,n-1)*dt/epsilon_0))./(dy*sigmay(m,n-1))).*(hz(m,n)-hz(m,n-1));
        m=1:nx;n=[npml+1,ny-npml];%up and down the cross section
%         ex(m,n)=ex(m,n)+dtepsdy*(hz(m,n)-hz(m,n-1));
        ex(m,n)=exp(-sigmay(m,n-1)*dt/epsilon_0).*ex(m,n)+((1-exp(-sigmay(m,n-1)*dt/epsilon_0))./(dy*sigmay(m,n-1))).*(hz(m,n)-hz(m,n-1));
        m=1:nx;n=[1,ny];%up and down the PEC boundary
        ex(m,n)=0;
%         m=1:nx;n=1;%mur-1
%             ex(m,n)=tempexy_(m,1)+(light_speed*dt-dy)/(light_speed*dt+dy)*(ex(m,n+1)-ex(m,n));
%         n=ny;
%             ex(m,n)=tempexy(m,1)+(light_speed*dt-dy)/(light_speed*dt+dy)*(ex(m,n-1)-ex(m,n));
%         %记录临时变量
%         tempexy_=ex(:,2);
%         tempexy=ex(:,ny-1);
        
        m=[1:npml,nx-npml+1:nx];n=1+npml+1:ny-npml-1;%right and left without the corner 
        ex(m,n)=ex(m,n)+dtepsdy*(hz(m,n)-hz(m,n-1));
    end
    ny=ny-1;
    
    %%%%%%%%%%%%%%%%%%%%%update the ey
    nx=nx+1;
    m=1+(npml+1):nx-(npml+1);n=1+npml:ny-npml;%the main area
    ey(m,n)=ey(m,n)-dtepsdx*(hz(m,n)-hz(m-1,n));
    if nargin==3
        %处理边界mur
        m=1;%x=xmin
            n=1:ny;
            ey(m,n)=tempeyx_(1,n)+(light_speed*dt-dx)/(light_speed*dt+dx)*(ey(m+1,n)-ey(m,n));
        m=nx;%x=xmax
            ey(m,n)=tempeyx(1,n)+(light_speed*dt-dx)/(light_speed*dt+dx)*(ey(m-1,n)-ey(m,n));
        %记录临时变量
        tempeyx_=ey(2,:);
        tempeyx=ey(nx-1,:);
    end
    if nargin==5
        m=[2:npml,nx-npml+1:nx-1];n=1:ny;%left with the corner
        ey(m,n)=exp(-sigmax(m-1,n)*dt/epsilon_0).*ey(m,n)-((1-exp(-sigmax(m-1,n)*dt/epsilon_0))./(dx*sigmax(m-1,n))).*(hz(m,n)-hz(m-1,n));
        m=[npml+1,nx-npml];n=1:ny;%right and left the cross section
        ey(m,n)=exp(-sigmax(m-1,n)*dt/epsilon_0).*ey(m,n)-((1-exp(-sigmax(m-1,n)*dt/epsilon_0))./(dx*sigmax(m-1,n))).*(hz(m,n)-hz(m-1,n));
%         ey(m,n)=ey(m,n)-dtepsdx*(hz(m,n)-hz(m-1,n));
        m=[1,nx];n=1:ny;%right and left the PEC boundary
        ey(m,n)=0;
%         m=1;n=1:ny;%mur-1
%             ey(m,n)=tempeyx_(1,n)+(light_speed*dt-dx)/(light_speed*dt+dx)*(ey(m+1,n)-ey(m,n));
%         m=nx;
%             ey(m,n)=tempeyx(1,n)+(light_speed*dt-dx)/(light_speed*dt+dx)*(ey(m-1,n)-ey(m,n));
%         %记录临时变量
%         tempeyx_=ey(2,:);
%         tempeyx=ey(nx-1,:);
        
        m=1+npml+1:nx-npml-1;n=[1:npml,ny-npml+1:ny];%up and down without the corner
        ey(m,n)=ey(m,n)-dtepsdx*(hz(m,n)-hz(m-1,n));
    end
    nx=nx-1;
    %%%%%%%%%%%%%save the magnitude of E every plot_modulus step
%     if iteration==90|iteration==120
%         figure;mesh(1:ny,1:nx,hz(1:nx,1:ny));
%     end
end
if nargin==2
    Hz=hz(151:250,176);
else
    Hz=hz(npml+1:nx-npml,npml+1);
%     Hz=hz(1:nx,1);
end