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

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

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

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

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

我是用fortran语言编写的，模拟的计算模型为一个200*200网格的计算区域，中间为真空，四周各有9个网格的完全匹配吸收介质层，具体代码如下: XsH(8-n0  
      program fdtd @M]uUL-ze  
c    implicit none $ 12mS  
    parameter(im=200,jm=200,k=4,np=9) y#Cp Vm#!>  
    real pi,c,ds,dt,miu,epson,lamda,sigma_max,tt jCJbmEfo9@  
    real ex(1:im,1:jm+1),ey(1:im+1,1:jm),bz(1:im,1:im), A~2U9f+\  
     *     bzx(1:im,1:jm),bzy(1:im,1:jm),bzxx(1:1000) ;f]p`!] 3  
Yg @&@S]  
c-----------set constent------------- S\\3?[!p  
         pi=3.1415 [H~Yg2O  
                          c=3.0e8 .2K4<UOAbm  
         lamda=(3.0e-7) q+j.
)e  
                          ds=lamda/20 WO}l&Q  
         dt=ds/2/c ywSV4ZtM  
         sigma_max=(k+1.0)/150.0/ds/pi 8Ce|Q8<8]  
c------------initial------------------         ~RdD6V  
    do i=1,im ,^Cl?\9"  
       do j=1,jm+1 *8QES
F9  
          ex(i,j)=0 su?{Cj6*  
       enddo cpIFjb>u{  
    enddo z XI [f  
    do i=1,im+1 g431+O0K1  
       do j=1,jm 9EPE.+ns  
          ey(i,j)=0 Z$UPLg3=;_  
       enddo 0XkLWl|k  
    enddo mYU7b8x_  
    do i=1,im ]q,5'[=~4h  
       do j=1,jm MC 8t"SB  
                              bz(i,j)=0 Fi7G S;  
       enddo zRO-oOJ  
    enddo +%O_xq
q  
o{W4@:Ib  
c---------start the main loop----------     ?&ow:OH+  
      do tt=1,1000 tQ,,krw~  
           do i=1,im .o27uB.  
             if ((i.ge.np+1).and.(i.le.im-np)) then b`W2^/D  
    di=0 uTWij4)a  
             elseif (i.le.np) then 7o+JQ&fF;  
                  di=np-i+0.5 a{4Wg:  
             elseif (i.ge.im-np+1) then gY\g+df-  
                  di=np+i-im-0.5 /.knZ_aJ!  
             end if   FW~{io]n  
            sigma_mx=sigma_max*((di/np)**k) S!LLC{  
            do j=1,jm   d9j+==S <  
                 if ((j.ge.np+1).and.(j.le.jm-np)) then ]JQ+*ZYUE  
                    dj=0 DH@]d0N  
                    elseif (j.le.np) then tK $r_*  
                    dj=np-j+0.5 #N
oY}*  
                    elseif (j.ge.jm-np+1) then O#}d!}SIp  
                    dj=np+j-jm-0.5 )aV\=a |A  
                 end if                               Y)Os]<N1  
                 sigma_my=sigma_max*((dj/np)**k) .5S< G)Ja  
                 if((sigma_mx.eq.0).and.(sigma_my.eq.0)) then KA[8NPhzZ  
         if((i.eq.100).and.(j.eq.100)) then I.4o9Z[?  
                      t=30                                                                   8!R +wy  
                                       term=(tt-t) g':/hlQ  
                                       pulse=exp(-4*pi*(term**2)/100) (f-Mm0%[  
                                       dey=ey(i+1,j)-ey(i,j) 4R c_C0O  
                                       dex=ex(i,j+1)-ex(i,j) 3?}\Hw  
                                       bz(i,j)=bz(i,j)-dt*c*(dey/ds-dex/ds)+pulse ?g~w6|U(r  
                     else      
v$WH#;(\  
                                            dey=ey(i+1,j)-ey(i,j) 9_$i.@L1  
                                            dex=ex(i,j+1)-ex(i,j) Hm>7|!  
                                            bz(i,j)=bz(i,j)-dt*c*(dey/ds-dex/ds) #hKaH -j  
                        end if mQbpv'N  
               else "k;j@  
                  cpm=(1-2*pi*sigma_mx*dt)/(1+2*pi*sigma_mx*dt) eX{:&Do  
                                   cqm=c*dt/(1+2*pi*sigma_mx*dt)/ds m'!smSx8  
                                   bzx(i,j)=cpm*bzx(i,j)-cqm*(ey(i+1,j)-ey(i,j))               xr;:gz!h  
                  cpm=(1-2*pi*sigma_my*dt)/(1+2*pi*sigma_my*dt) ;0Vyim)S]  
                                   cqm=c*dt/(1+2*pi*sigma_my*dt)/ds  L+=pEk_  
                                   bzy(i,j)=cpm*bzy(i,j)+cqm*(ex(i,j+1)-ex(i,j)) ~BUzyc%  
                                   bz(i,j)=bzx(i,j)+bzy(i,j) H3}eFl=i2  
               end if s(0S)l<  
      enddo $] gwaJ:  
         enddo 6\+ZTw  
         do i=2,im : \{>+!`w  
                if ((i.ge.np+1).and.(i.le.im-np)) then )=k8W9i8b  
              di=0 %Voq"}}N  
                elseif (i.le.np) then 7Gos-_s  
                          di=np-i+1 }L!%^siG_  
                elseif (i.ge.im-np+1) then #Epx'$9  
                         di=np+i-im-1 Wl29xY}`{!  
                 end if `<?{%ja  
          sigma_ex=sigma_max*((di/np)**k)
Q;V*M  
                 do j=1,jm N!W# N$  
                      cam=(1-2*pi*sigma_ex*dt)/(1+2*pi*sigma_ex*dt) T#o?@;  
                      cbm=c*dt/(1+2*pi*sigma_ex*dt)/ds |vGb,&3  
     ey(i,j)=cam*ey(i,j)-cbm*(bz(i,j)-bz(i-1,j)) `wMHjcUP  
                 enddo X\=m  
    enddo :2Fy`PPab  
         do i=1,im /Ezx'h3Q  
               do j=2,jm Dc1tND$X3g  
                   if ((j.ge.np+1).and.(j.le.jm-np)) then 5PcN$r"P  
                   dj=0 u|G&CV#r  
                   elseif (j.le.np) then A89n^@  
                   dj=np-j+1 ?Nbc#0pb7  
                   elseif (j.ge.jm-np+1) then 9=l6NNe)|  
                   dj=np+j-jm-1 F2N)|C<  
                 end if rfz\DvVd  
               sigma_ey=sigma_max*((dj/np)**k) v^;p]_c~2  
               cam=(1-2*pi*sigma_ey*dt)/(1+2*pi*sigma_ey*dt)
03%`ouf  
               cbm=c*dt/(1+2*pi*sigma_ey*dt)/ds 7])cu>/  
               ex(i,j)=cam*ex(i,j)+cbm*(bz(i,j)-bz(i,j-1)) oP"X-I  
            enddo [|vE*&:uO  
         enddo @)\{u$  
     bzxx(tt)=bz(190,100) `U.VfQR:  
   enddo odPdWV,&*  
    open(7,file='bzxx-gaosi-pml-2.txt') ~xp(k  
                     write(7,30) bzxx 'Nqa=_<WW  
                     close(7) O(_a6s+m  
30                format(4x,1f20.8) \ZOH3`vq  
                    end

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

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

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

你发给我的代码，高斯制和国际单位制下的，递推公式怎么都不一样的。 #<DS-^W!  
国际单位制下是： fL~@v-l#~  
cam=exp(-sigma_ex*del_t/epson); dLbSvK<(I  
PML层分了三种情况处理。( O"~CZh,:r}  
sigma_mx == 0 && sigmx_my!=0) vLIaTr gz  
sigma_my=0 && sigma_mx!=0 v2Vmcc_]9x  
sigma_my!=0 && sigma_mx!=0 x5R|,bY  
高斯制下是： 5S 4Bz  
cam=(1-2*pi*sigma_ex*del_t)/(1+2*pi*sigma_ex*del_t); 6PT"9vR`)  
PML层只有一种情况呢 |_mN:(3  
按道理来说，他们公式都应该一样的，差的只是系数而已。 +?v2MsF']  
2" u,f  
是因为什么考虑才使他们不一样的呢？

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