[已解决]关于高斯制下完全匹配层的差分格式问题

公式对的啊，没问题啊。应该是code里面的问题了。 `*|LI  
你有没有试过国际单位制下面的FDTD+PML啊。那种情况出问题么？

要不你把代码发给我看看，或者我发一个国际单位制的2D FDTD+PML给你参考一下？

国际单位制下的FDTD+PML我已经编过了，没有问题，我用刚才我发的那个我自己构造的高斯制完全匹配层控制方程加上FDTD模拟真空场加边界失败了，那我把我编的这个高斯制程序代码发给你看看，你帮我瞧瞧

好的，我的email: gwei_zhao@126.com

我是用fortran语言编写的，模拟的计算模型为一个200*200网格的计算区域，中间为真空，四周各有9个网格的完全匹配吸收介质层，具体代码如下: 7n3x19T  
      program fdtd =;-ju@d  
c    implicit none :."n@sA@  
    parameter(im=200,jm=200,k=4,np=9) Dg2#Gv0B  
    real pi,c,ds,dt,miu,epson,lamda,sigma_max,tt (? #U&  
    real ex(1:im,1:jm+1),ey(1:im+1,1:jm),bz(1:im,1:im), <`VJU2  
     *     bzx(1:im,1:jm),bzy(1:im,1:jm),bzxx(1:1000) }5c'ui!3H  
d)J] Y=j
 
c-----------set constent------------- &>^Ympr  
         pi=3.1415 OuJy$e  
                          c=3.0e8 g he=mQ-  
         lamda=(3.0e-7) ;@;
ie8H  
                          ds=lamda/20 ObHz+qRG  
         dt=ds/2/c RSVN(-wIi)  
         sigma_max=(k+1.0)/150.0/ds/pi :!*;0~#  
c------------initial------------------         sNan"  
    do i=1,im 0kr& c;~  
       do j=1,jm+1 zdQu%q  
          ex(i,j)=0 W*WH .1&  
       enddo CIs1*:Q9  
    enddo otnY{r
*  
    do i=1,im+1 E/Gs',Y  
       do j=1,jm /3:IE%o  
          ey(i,j)=0 JbD)}(G;  
       enddo U\!LZ?gC  
    enddo U3R`mHr0  
    do i=1,im sMpC4E  
       do j=1,jm e+[J[<8  
                              bz(i,j)=0 ]KV8u1H>  
       enddo ]Pl6:FB8%@  
    enddo 6"c!tJc7j  
j$i8@]  
c---------start the main loop----------     ^Y*.Ktp,o  
      do tt=1,1000 !/q&0a  
           do i=1,im b}{9 :n/SC  
             if ((i.ge.np+1).and.(i.le.im-np)) then 4O,a`:d1$6  
    di=0
L08
;z  
             elseif (i.le.np) then Svp
Ts  
                  di=np-i+0.5 VDI S`E  
             elseif (i.ge.im-np+1) then oDiv9jm  
                  di=np+i-im-0.5 w"Y55EURB  
             end if   ofhZ@3  
            sigma_mx=sigma_max*((di/np)**k) `uJ l<kHI  
            do j=1,jm   L\'qAfRZ  
                 if ((j.ge.np+1).and.(j.le.jm-np)) then %!/liS  
                    dj=0 ve/6-J!5Y.  
                    elseif (j.le.np) then aRb:.\ \zc  
                    dj=np-j+0.5 LCm}v&~%A  
                    elseif (j.ge.jm-np+1) then ra#)*fG,~  
                    dj=np+j-jm-0.5 lSsFI30  
                 end if                               sf&K<C](  
                 sigma_my=sigma_max*((dj/np)**k) \YF!< 2|[  
                 if((sigma_mx.eq.0).and.(sigma_my.eq.0)) then KKTfxNxJn  
         if((i.eq.100).and.(j.eq.100)) then :s5<AT Q  
                      t=30                                                                   T%vbD*nt.  
                                       term=(tt-t) )RKhEm%Vr2  
                                       pulse=exp(-4*pi*(term**2)/100) vWU4ZBT8G  
                                       dey=ey(i+1,j)-ey(i,j) Q X5#$-H@  
                                       dex=ex(i,j+1)-ex(i,j) #3ro?w  
                                       bz(i,j)=bz(i,j)-dt*c*(dey/ds-dex/ds)+pulse &=t(N
I$  
                     else       g,]5&C T3v  
                                            dey=ey(i+1,j)-ey(i,j) 4]Nr$FY  
                                            dex=ex(i,j+1)-ex(i,j) gJfL$S'w  
                                            bz(i,j)=bz(i,j)-dt*c*(dey/ds-dex/ds) G Z[5m[  
                        end if |4*2xDcl  
               else *:\9T#h  
                  cpm=(1-2*pi*sigma_mx*dt)/(1+2*pi*sigma_mx*dt) S{|)9EKw  
                                   cqm=c*dt/(1+2*pi*sigma_mx*dt)/ds {mr)n3  
                                   bzx(i,j)=cpm*bzx(i,j)-cqm*(ey(i+1,j)-ey(i,j))               1+ARV&bc  
                  cpm=(1-2*pi*sigma_my*dt)/(1+2*pi*sigma_my*dt) OL+40J  
                                   cqm=c*dt/(1+2*pi*sigma_my*dt)/ds z_0lMX`  
                                   bzy(i,j)=cpm*bzy(i,j)+cqm*(ex(i,j+1)-ex(i,j)) xB]v  
                                   bz(i,j)=bzx(i,j)+bzy(i,j) wI]>0geb*  
               end if @x>2|`65Y  
      enddo g.VIe  
         enddo [ :(M<u`y>  
         do i=2,im XA^:n+Yo  
                if ((i.ge.np+1).and.(i.le.im-np)) then X#C7r@H  
              di=0 >knR>96  
                elseif (i.le.np) then &!aLOx*3`  
                          di=np-i+1 s~ZC!-[;  
                elseif (i.ge.im-np+1) then }^bL'
 
                         di=np+i-im-1 +s 0
Bt '  
                 end if _BA_lkN+D  
          sigma_ex=sigma_max*((di/np)**k) g\h7`-#t  
                 do j=1,jm ?mUu(D:7D  
                      cam=(1-2*pi*sigma_ex*dt)/(1+2*pi*sigma_ex*dt) J0z0%p   
                      cbm=c*dt/(1+2*pi*sigma_ex*dt)/ds `r bqYU0  
     ey(i,j)=cam*ey(i,j)-cbm*(bz(i,j)-bz(i-1,j)) !oRm.cO  
                 enddo q2pao?aa  
    enddo y:Ab5/bHy  
         do i=1,im 31\^9w__8  
               do j=2,jm ,TA[el%#  
                   if ((j.ge.np+1).and.(j.le.jm-np)) then Qy,qQA/  
                   dj=0 xE[tD? M{  
                   elseif (j.le.np) then %QZ!Tb  
                   dj=np-j+1 1:Gd{z  
                   elseif (j.ge.jm-np+1) then l\y*wr`  
                   dj=np+j-jm-1 xR%a
yT.  
                 end if r:Tb{cA  
               sigma_ey=sigma_max*((dj/np)**k) $?Mz[X  
               cam=(1-2*pi*sigma_ey*dt)/(1+2*pi*sigma_ey*dt) L#N.pd  
               cbm=c*dt/(1+2*pi*sigma_ey*dt)/ds u1Yp5jp^K  
               ex(i,j)=cam*ex(i,j)+cbm*(bz(i,j)-bz(i,j-1)) 0cU^ue%  
            enddo _NW OSt  
         enddo f__WnW5h  
     bzxx(tt)=bz(190,100) 7VBw@Rh  
   enddo 7anpz%  
    open(7,file='bzxx-gaosi-pml-2.txt') lR3^&d72?  
                     write(7,30) bzxx wx'Tv  
                     close(7) uzA'D~)P  
30                format(4x,1f20.8) Q]:%Jj2  
                    end

对了，你能不能把对应的国际单位制下的程序也发一下，这样我对比一下，看起来容易点。 [{ pc1U-  
没办法，我家机器上现在除了游戏，没其他的了。

好，明天到教研室发给你哈，寝室网速很慢

关于程序代码我已经发到你邮箱里面了，请查收一下，我的编程语言可能不是很规范，你先看看，要是有必要的话，你就说声，我重新整理一下搞个注释得比较详尽的版本给你，时间匆忙，没来得及整理好就发给你了，实在抱歉！！ iT%UfN/q=I  
或者你可以直接用ｑｑ与我及时沟通，我的ｑｑ号码是：２４８２８５４９８，臜们坛外交流！！衷心感谢！！

你发给我的代码，高斯制和国际单位制下的，递推公式怎么都不一样的。 {R#nGsrt;  
国际单位制下是： z}&<DYD  
cam=exp(-sigma_ex*del_t/epson); f`cz@  
PML层分了三种情况处理。( N>~*Jp2;  
sigma_mx == 0 && sigmx_my!=0) LDNp
EX~  
sigma_my=0 && sigma_mx!=0 J+TYm%A;-  
sigma_my!=0 && sigma_mx!=0 iZ:-V8{  
高斯制下是： v!x[1[  
cam=(1-2*pi*sigma_ex*del_t)/(1+2*pi*sigma_ex*del_t); !nu['6I%  
PML层只有一种情况呢 =_XcG!"  
按道理来说，他们公式都应该一样的，差的只是系数而已。 l}wBthwCc  
 }QI*Ns  
是因为什么考虑才使他们不一样的呢？

非常抱歉，今天我们有集体活动到晚上才来到教研室，关于你说的这个情况，国际制下的这个PML分三种情况是出于应用指数差分的需要，因为当其中某一参数为零时就不能出现在分母里面。我最初的版本时不用指数差分的，和高斯制的分类方法就完全一样，等会我把我的那个国际制版本发给你 ，前面是我搞错了，sorry!