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

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

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

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

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

我是用fortran语言编写的，模拟的计算模型为一个200*200网格的计算区域，中间为真空，四周各有9个网格的完全匹配吸收介质层，具体代码如下: fzJiW@-T  
      program fdtd b
8@}Jv  
c    implicit none D <iG*I  
    parameter(im=200,jm=200,k=4,np=9) Hk}P  
    real pi,c,ds,dt,miu,epson,lamda,sigma_max,tt C. .|O  
    real ex(1:im,1:jm+1),ey(1:im+1,1:jm),bz(1:im,1:im), v[$e{Dz(  
     *     bzx(1:im,1:jm),bzy(1:im,1:jm),bzxx(1:1000) Es[3Ppz  
so`\e^d  
c-----------set constent------------- D1RQkAZS  
         pi=3.1415 WL+EpNKSf  
                          c=3.0e8 7J9<
B5U  
         lamda=(3.0e-7) =c6d$  
                          ds=lamda/20 z>!./z]p  
         dt=ds/2/c ]-wyZ +a  
         sigma_max=(k+1.0)/150.0/ds/pi >)4~,-;k  
c------------initial------------------         }R*%q  
    do i=1,im b'O/u."O  
       do j=1,jm+1 7{r7  
          ex(i,j)=0 nAQ[ -NbW,  
       enddo g3
ukx$Q{>  
    enddo fHaF9o+/b  
    do i=1,im+1 !
Vl)aL  
       do j=1,jm Ku'a,\7z  
          ey(i,j)=0 #
?Ix6 {R  
       enddo ^iH[ 22b4  
    enddo -c>3|bo  
    do i=1,im Yup#aeXY/  
       do j=1,jm L 0Ckw},,  
                              bz(i,j)=0 $?ss5: S  
       enddo Y$OE[nGi%X  
    enddo Wps^wY  
LuRCkKJ  
c---------start the main loop----------     v}!lx)#  
      do tt=1,1000 O]IAIM  
           do i=1,im %qV:h#  
             if ((i.ge.np+1).and.(i.le.im-np)) then e/4C` J-  
    di=0 V dJ  
             elseif (i.le.np) then TFHYB9vV  
                  di=np-i+0.5 UR^r>  
             elseif (i.ge.im-np+1) then 
^2dQVV.  
                  di=np+i-im-0.5 EwBrOq`C  
             end if   13@emb  
            sigma_mx=sigma_max*((di/np)**k) hGKQK ^bn  
            do j=1,jm   +184|nJ<2  
                 if ((j.ge.np+1).and.(j.le.jm-np)) then &tUX(  
                    dj=0 5Y;&L!T  
                    elseif (j.le.np) then uBG!R
#T  
                    dj=np-j+0.5 o5]-Kuw`  
                    elseif (j.ge.jm-np+1) then vAP1PQX;  
                    dj=np+j-jm-0.5 odL*_<Z  
                 end if                               NP3 e^  
                 sigma_my=sigma_max*((dj/np)**k) e=##X}4zZ  
                 if((sigma_mx.eq.0).and.(sigma_my.eq.0)) then gR/?MJ(v  
         if((i.eq.100).and.(j.eq.100)) then 1TEKq#t;y  
                      t=30                                                                   dOaOWMrfdf  
                                       term=(tt-t) )W_akUL  
                                       pulse=exp(-4*pi*(term**2)/100) $mT)<N ;w  
                                       dey=ey(i+1,j)-ey(i,j) B
uvnY  
                                       dex=ex(i,j+1)-ex(i,j) ?0 cv  
                                       bz(i,j)=bz(i,j)-dt*c*(dey/ds-dex/ds)+pulse 6# bTlmcg  
                     else       0`pCgF  
                                            dey=ey(i+1,j)-ey(i,j) ,]t_9B QK  
                                            dex=ex(i,j+1)-ex(i,j) +6:
  
                                            bz(i,j)=bz(i,j)-dt*c*(dey/ds-dex/ds) -V2f.QE%  
                        end if B|o@|zF  
               else `j3 OFC{7E  
                  cpm=(1-2*pi*sigma_mx*dt)/(1+2*pi*sigma_mx*dt)
(0@b4}Z  
                                   cqm=c*dt/(1+2*pi*sigma_mx*dt)/ds 5ms]Wbh)  
                                   bzx(i,j)=cpm*bzx(i,j)-cqm*(ey(i+1,j)-ey(i,j))               P*[wB_^&UP  
                  cpm=(1-2*pi*sigma_my*dt)/(1+2*pi*sigma_my*dt) 3Z~_6P^ +N  
                                   cqm=c*dt/(1+2*pi*sigma_my*dt)/ds $J4)z&%dr  
                                   bzy(i,j)=cpm*bzy(i,j)+cqm*(ex(i,j+1)-ex(i,j)) U `lp56  
                                   bz(i,j)=bzx(i,j)+bzy(i,j) uju'Bs7  
               end if &CQ28WG X  
      enddo 1Y"9<ry  
         enddo y ~-v0/  
         do i=2,im  "O# V/(  
                if ((i.ge.np+1).and.(i.le.im-np)) then Y{~`g(~9_A  
              di=0 p! k~ufU  
                elseif (i.le.np) then `s69p'<;p  
                          di=np-i+1 ;eo}/-a_Xw  
                elseif (i.ge.im-np+1) then D\;5{,:d  
                         di=np+i-im-1 h143HXBi1+  
                 end if pSAtn  
          sigma_ex=sigma_max*((di/np)**k) ga,kKPL  
                 do j=1,jm Ze/\IBd  
                      cam=(1-2*pi*sigma_ex*dt)/(1+2*pi*sigma_ex*dt) Pc?"H!Hkn  
                      cbm=c*dt/(1+2*pi*sigma_ex*dt)/ds g3Q;]8Y&  
     ey(i,j)=cam*ey(i,j)-cbm*(bz(i,j)-bz(i-1,j)) "C3J[) qC  
                 enddo W2%@}IDm  
    enddo WNeBthq6  
         do i=1,im X!{K`~DRX  
               do j=2,jm /)RH-_63  
                   if ((j.ge.np+1).and.(j.le.jm-np)) then d %FLk=]  
                   dj=0 -k:x e:$  
                   elseif (j.le.np) then Q e/XEW  
                   dj=np-j+1 $[Ut])4 ~  
                   elseif (j.ge.jm-np+1) then .p Mwa  
                   dj=np+j-jm-1 tYgHJ~1L*  
                 end if *N: $,x
f  
               sigma_ey=sigma_max*((dj/np)**k) o/&K>]8M  
               cam=(1-2*pi*sigma_ey*dt)/(1+2*pi*sigma_ey*dt) eU)QoVt  
               cbm=c*dt/(1+2*pi*sigma_ey*dt)/ds "3
++S  
               ex(i,j)=cam*ex(i,j)+cbm*(bz(i,j)-bz(i,j-1)) P B-x_D  
            enddo d=D#cs;\  
         enddo ^?0'\
Z  
     bzxx(tt)=bz(190,100) )zy;!  
   enddo r#oJch=  
    open(7,file='bzxx-gaosi-pml-2.txt') \ C$t  
                     write(7,30) bzxx +2tFX  
                     close(7) Qe!3ae`Z  
30                format(4x,1f20.8) Jza?DhSAZ  
                    end

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

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

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

你发给我的代码，高斯制和国际单位制下的，递推公式怎么都不一样的。 Nd_fjB  
国际单位制下是： CI|lJ  
cam=exp(-sigma_ex*del_t/epson); f0Q6sVZHa  
PML层分了三种情况处理。( dj|5'<l2  
sigma_mx == 0 && sigmx_my!=0) ?0tg}0|  
sigma_my=0 && sigma_mx!=0 !,J# r  
sigma_my!=0 && sigma_mx!=0 WA{igj@\  
高斯制下是： B*7kX&Uq  
cam=(1-2*pi*sigma_ex*del_t)/(1+2*pi*sigma_ex*del_t); cw;wv+|k  
PML层只有一种情况呢 T X`X5j  
按道理来说，他们公式都应该一样的，差的只是系数而已。 7\u+%i;YZ  
zd?@xno  
是因为什么考虑才使他们不一样的呢？

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