[讨论]写了个UPML吸收边界的一维FDTD程序，为什么结果图形显示的振幅始终为0.5左右

图片:Graph7.jpg

1|gP :t}  
>;
W)tc,  
e('c9 Y  
我写了个UPML吸收边界的一维FDTD程序，为什么结果图形显示的振幅始终为0.5， yl'~H;su  
6!"15dPN  
源是高斯点源。软源：  xx[i1] += J0*cos(omega*time + phase);  J0=1是振幅。 !)9zH  
Zo0&<QWj  
我实际默认的振幅是1.0啊，不知道为什么？ 大家有碰到过吗？ RPiCX
pJv&  
>uuX<\cW  
i^IvT  
,Fr{i1Ky  
按照版主的提示我找了，下面是前一百个时间步源输出的结果， 看不出有什么异样。 请教大家！ z|b4w7I  
&6\rKOsn  
图片是每隔100个时间步输出的Hz的结果图形。 pb{P[-f  
AN~1E@"  
Process : 0 -- Start FDTD 6U/wFT!7$  
^5T{x>Lj  
xx[i1] = 0.960239    --PT_Sin ]owH [wvX  
xx[i1] = 0.962533    --PT_Sin Xj/X.  
xx[i1] = -325.327    --PT_Sin ,OasT!Sr  
xx[i1] = 0.2498    --PT_Sin iuHG9#n  
xx[i1] = -145.996    --PT_Sin `a6;*
ry  
xx[i1] = 0.211419    --PT_Sin F'#3wCzt  
xx[i1] = -48.8379    --PT_Sin Xj-3C[8@  
xx[i1] = 0.401967    --PT_Sin zWY6D4   
xx[i1] = -154.021    --PT_Sin |;_ yAL  
xx[i1] = -0.178958    --PT_Sin vGAPQg6*  
xx[i1] = 36.4055    --PT_Sin 2h
px%H  
xx[i1] = -0.131179    --PT_Sin UJm`GO  
xx[i1] = 77.6417    --PT_Sin 0"QE,pLe4  
xx[i1] = -0.157616    --PT_Sin Xl aNR+  
xx[i1] = 42.0433    --PT_Sin ]
Pd*w`R  
xx[i1] = -0.593708    --PT_Sin J`mp8?;%  
xx[i1] = 208.648    --PT_Sin 8%|x)  
xx[i1] = -0.422673    --PT_Sin cKfYkJ)A'  
xx[i1] = 183.003    --PT_Sin r=w%"3vb^  
xx[i1] = -0.497458    --PT_Sin <K0lS;@K  
xx[i1] = 165.057    --PT_Sin WWO jyj  
xx[i1] = -0.705738    --PT_Sin n
K|"
;  
xx[i1] = 265.784    --PT_Sin fzKKK+  
xx[i1] = -0.401242    --PT_Sin 8EE7mEmLH  
xx[i1] = 169.669    --PT_Sin $o
@?D^  
xx[i1] = -0.447727    --PT_Sin q@!:<Ra,){  
xx[i1] = 147.568    --PT_Sin rb_G0/R  
xx[i1] = -0.423352    --PT_Sin ZE
\t{s0  
xx[i1] = 171.735    --PT_Sin v[|iuOU  
xx[i1] = -0.052944    --PT_Sin n)=&=Uj`f  
xx[i1] = 31.8007    --PT_Sin =_1" d$S&  
xx[i1] = -0.0847835    --PT_Sin 3|?fGT;P  
xx[i1] = 15.2583    --PT_Sin E+2y-B)E  
xx[i1] = 0.0853313    --PT_Sin [)Ge^yI7  
xx[i1] = -12.5427    --PT_Sin r"Bf@va  
xx[i1] = 0.395852    --PT_Sin ;|^fAc~9{r  
xx[i1] = -145.757    --PT_Sin <^~F~]wnH  
xx[i1] = 0.313667    --PT_Sin RTU:J67E  
xx[i1] = -129.753    --PT_Sin UF{2Gx  
xx[i1] = 0.502903    --PT_Sin wd]Yjr#%Ii  
xx[i1] = -168.642    --PT_Sin 67g/(4&  
xx[i1] = 0.636916    --PT_Sin zR .MXr  
xx[i1] = -246.076    --PT_Sin -(iJ<  
xx[i1] = 0.456816    --PT_Sin v8X&H
 
xx[i1] = -179.002    --PT_Sin "_l[4
o[D  
xx[i1] = 0.573885    --PT_Sin ~8X'p6  
xx[i1] = -199.4    --PT_Sin MJC Yi<D  
xx[i1] = 0.492852    --PT_Sin }F.1j!71L  
xx[i1] = -200.143    --PT_Sin 5[R}MhLZ  
xx[i1] = 0.238317    --PT_Sin  A:!{+  
xx[i1] = -92.1195    --PT_Sin bCrB'&^t  
xx[i1] = 0.274161    --PT_Sin OiOL4}5(  
xx[i1] = -93.1865    --PT_Sin a@a1/3  
xx[i1] = 0.0412927    --PT_Sin i!HGM=f  
xx[i1] = -34.609    --PT_Sin 7/1S5yUr|  
xx[i1] = -0.192185    --PT_Sin B%pvk.`  
xx[i1] = 75.2772    --PT_Sin B'
;Ob  
xx[i1] = -0.175713    --PT_Sin y,x~S\>+  
xx[i1] = 69.7022    --PT_Sin J91`wA&r  
xx[i1] = -0.433669    --PT_Sin +?<jSmGW  
xx[i1] = 145.059    --PT_Sin t}tKm  
xx[i1] = -0.532358    --PT_Sin %G@aZWk Sa  
xx[i1] = 209.3    --PT_Sin gvJJ.IX]+  
xx[i1] = -0.461939    --PT_Sin X"0Q)  
xx[i1] = 172.704    --PT_Sin 4(&'V+o  
xx[i1] = -0.630279    --PT_Sin <#Lw.;(U;k  
xx[i1] = 224.702    --PT_Sin lV 9q;!/1  
xx[i1] = -0.535067    --PT_Sin /7#
&q
x8  
xx[i1] = 215.544    --PT_Sin fkG8,=  
xx[i1] = -0.390707    --PT_Sin 1=X=jPwO C  
xx[i1] = 142.803    --PT_Sin w-"&;klV  
xx[i1] = -0.43058    --PT_Sin 3q>"#+R.t  
xx[i1] = 157.506    --PT_Sin ;H=6u  
xx[i1] = -0.179445    --PT_Sin Lv4=-mWv&0  
xx[i1] = 84.0556    --PT_Sin 
:1=?/8h  
xx[i1] = -0.0186302    --PT_Sin G#V22Wca8  
xx[i1] = 0.282467    --PT_Sin 'KL(A-}!  
xx[i1] = 0.0216978    --PT_Sin e&5K]W0{  
xx[i1] = -5.42531    --PT_Sin ~V&R
eW/  
xx[i1] = 0.311431    --PT_Sin Dk-L4FS  
xx[i1] = -102.723    --PT_Sin dF,FH-  
xx[i1] = 0.38785    --PT_Sin {2x5 V#6  
xx[i1] = -155.329    --PT_Sin VJ"3G;;  
xx[i1] = 0.41939    --PT_Sin EyeLC6u  
xx[i1] = -150.535    --PT_Sin t5k&xV=~ #  
xx[i1] = 0.618776    --PT_Sin UE4#j\  
xx[i1] = -224.607    --PT_Sin 8if"U xV(  
xx[i1] = 0.536137    --PT_Sin zaZ}:N/w(z  
xx[i1] = -213.51    --PT_Sin > 95Cs`>d  
xx[i1] = 0.499646    --PT_Sin ts]7 + 6V  
xx[i1] = -179.163    --PT_Sin GCKl[<9*  
xx[i1] = 0.542162    --PT_Sin t>fB
@xHBB  
xx[i1] = -204.068    --PT_Sin Rv-o__C!  
xx[i1] = 0.309125    --PT_Sin NIfc/%  
xx[i1] = -128.745    --PT_Sin q{hq.KZ  
xx[i1] = 0.220424    --PT_Sin twtDyo(\  
xx[i1] = -74.286    --PT_Sin |-fx 0y   
xx[i1] = 0.1345    --PT_Sin K)?^
b|D  
xx[i1] = -58.0096    --PT_Sin jveRiW@  
xx[i1] = -0.15247    --PT_Sin ~roHnJ>  
xx[i1] = 47.4533    --PT_Sin JdHc'WtS!|  
xx[i1] = -0.213577    --PT_Sin ,gvX
~k  
xx[i1] = 88.5413    --PT_Sin ie!4z34  
xx[i1] = -0.333515    --PT_Sin W!k6qTz)  
xx[i1] = 114.074    --PT_Sin >DRs(~|V#  
xx[i1] = -0.543425    --PT_Sin [1C#[Vla  
xx[i1] = 200.425    --PT_Sin DR /)hAE  
xx[i1] = -0.492682    --PT_Sin @nP}q!y  
xx[i1] = 193.356    --PT_Sin Pb,^UFa=  
xx[i1] = -0.555534    --PT_Sin mSfhl(<L  
xx[i1] = 197.493    --PT_Sin 76b7-Nj"  
xx[i1] = -0.599966    --PT_Sin ^H4iHjg  
xx[i1] = 229.339    --PT_Sin arP+(1U  
xx[i1] = -0.415018    --PT_Sin c)8wO=!  
xx[i1] = 163.493    --PT_Sin )ta5y7np  
xx[i1] = -0.39503    --PT_Sin DBUwf1=qj  
xx[i1] = 139.644    --PT_Sin h+UscdUl  
xx[i1] = -0.27881    --PT_Sin o0'av+e7  
xx[i1] = 115.171    --PT_Sin 7gwZ9Fob  
xx[i1] = -0.0255249    --PT_Sin cPcV[6)5K9  
xx[i1] = 14.5886    --PT_Sin r!^\Q7  
xx[i1] = 0.0247965    --PT_Sin 2M?lgh4"  
xx[i1] = -15.1405    --PT_Sin }gW/heUE  
xx[i1] = 0.212232    --PT_Sin Z?.*.<"Sj  
xx[i1] = -66.3698    --PT_Sin m_2P{  
xx[i1] = 0.413393    --PT_Sin C]fTV{  
xx[i1] = -155.249    --PT_Sin ib_Gy77Os  
xx[i1] = 0.407178    --PT_Sin WDdi}i>2  
xx[i1] = -156.635    --PT_Sin /09=Tyy/\  
xx[i1] = 0.553133    --PT_Sin Y14R"*t~  
xx[i1] = -195.756    --PT_Sin -y?Z}5-rs  
xx[i1] = 0.599839    --PT_Sin shT[|@"C  
xx[i1] = -231.548    --PT_Sin 4E>(Y98  
xx[i1] = 0.485114    --PT_Sin *fSM'q;  
xx[i1] = -184.438    --PT_Sin Rda1X~-g  
xx[i1] = 0.526204    --PT_Sin yk<jlVF$j  
xx[i1] = -189.969    --PT_Sin nY9qYFw  
xx[i1] = 0.398243    --PT_Sin a*j <TR  
xx[i1] = -161.091    --PT_Sin w<Cmzkf  
xx[i1] = 0.203638    --PT_Sin !&O/7ywe  
xx[i1] = -76.8309    --PT_Sin u;Eu<jU1  
xx[i1] = 0.16127    --PT_Sin Wp}9%Mq~Jy  
xx[i1] = -57.8966    --PT_Sin dV{Hn {(  
xx[i1] = -0.0668099    --PT_Sin xbC8Amo;8"  
xx[i1] = 11.6205    --PT_Sin Y:%)cUxA  
xx[i1] = -0.242777    --PT_Sin ^P/D8cXa4  
xx[i1] = 94.0426    --PT_Sin Wk?|BR]O  
xx[i1] = -0.286566    --PT_Sin 7omGg~!k(  
xx[i1] = 106.608    --PT_Sin C..2y4bA}  
xx[i1] = -0.492734    --PT_Sin }+ 2"?f|]  
xx[i1] = 173.875    --PT_Sin #2jn4>  
xx[i1] = -0.543016    --PT_Sin >iH).:j  
xx[i1] = 211.046    --PT_Sin |WU`p  
xx[i1] = -0.511151    --PT_Sin xm|4\H&Bg  
xx[i1] = 189.127    --PT_Sin H?j-=Zka  
xx[i1] = -0.601816    --PT_Sin 3:joSQa  
xx[i1] = 220.344    --PT_Sin S}^s5ztm  
xx[i1] = -0.482056    --PT_Sin (Dm"e`  
xx[i1] = 191.811    --PT_Sin eCIRt/ uA  
xx[i1] = -0.363158    --PT_Sin W8$=a
  
xx[i1] = 132.74    --PT_Sin :{:?D\%6  
xx[i1] = -0.327607    --PT_Sin B"m:<@ "  
xx[i1] = 123.647    --PT_Sin !X%!7wsc  
xx[i1] = -0.0890027    --PT_Sin 5 ?~-Vv31s  
xx[i1] = 45.2235    --PT_Sin 2ZbY|8X$r  
xx[i1] = 0.0493984    --PT_Sin x}<G
!*3  
xx[i1] = -23.2674    --PT_Sin T[h}A"yK;  
xx[i1] = 0.141235    --PT_Sin s
Nfb %r  
xx[i1] = -47.8013    --PT_Sin V-;nj,.mY  
xx[i1] = 0.380387    --PT_Sin X/-KkC  
xx[i1] = -133.801    --PT_Sin d Zz^9:C+  
xx[i1] = 0.435991    --PT_Sin 0ITA3v8{  
xx[i1] = -170.321    --PT_Sin CS5jJi"pD3  
xx[i1] = 0.489207    --PT_Sin NzAtdcwR  
xx[i1] = -176.668    --PT_Sin (1SO;8k\  
xx[i1] = 0.615248    --PT_Sin o+-Ge J  
xx[i1] = -227.859    --PT_Sin ./nYXREO|  
xx[i1] = 0.522702    --PT_Sin udD*E~1q  
xx[i1] = -204.76    --PT_Sin 7G[ GHc>  
xx[i1] = 0.487312    --PT_Sin 'Z2N{65  
xx[i1] = -176.404    --PT_Sin M.:@<S  
xx[i1] = 0.459206    --PT_Sin A Ok7G?Y  
xx[i1] = -175.971    --PT_Sin
5Vnr"d  
n5k^v$'  

未注册仅能浏览部分内容，查看全部内容及附件请先 登录 或 注册

那应该是你程序问题吧，你检查一下就好了。 c?p0#3%L#  
比如你源是加在Z_source这一点的，你看这一点的幅度是多少，Z_source +/- 1 位置的幅度是多少.

谢谢版主！

Qww^P/vm  
P2t_T'R}  
顶起！

你把你的代码贴出来看看

c3#q0Ma  
/////////////////////////////////////////////////////////////////////////////////////// e9:P9Di(b  
//Set the PML matrices ;UpJ=?W  
/////////////////////////////////////////////////////////////////////////////////////// K}K)`bifw  
void FDTD_1D_EyHz::Set_PML_Param_1D() 8lb-}=  
{ -TS? fne)  
    for (i = 0; i < nx; i++) ESv:1o`?n  
    { hfv%,,e  
        K_E1_a = 1.0;                // >0T Za  
        K_E1_b = dt/(mu_0*mu_r*dx);            // v)+@XU2wZ  
Q;w
B{vr$  
        if (i < nx-1) az2Xch]  
        { 8(Fu  
            K_E2_a = 1.0;            // [O1|75  
            K_E2_b = dt/(eps_0*eps_r*dx);    // dt/(eps_r*dx);     !^L-T?y.2  
        }             >)3VbO  
    } vYdlSe=6G  
pX_b6%yX(  
    //PML parameters !.-.#<<_a  
    double ka_max = 1; DEtf(lW_  
    int exponent = 4; GO~k '  
    double R_err = 1e-16; U0IE1_R  
    double eta = sqrt(mu_0*mu_r/eps_0/eps_r); Q1T@oxV  
    double sigma_x, sigma_max, ka_x; PUCx]5  
         A?,A(-0C  
    sigma_max= -(exponent+1)*log(R_err)/(2*eta*n_PML*dx); y<*-tZV[  
    //sigma_max= -(exponent+1)*log(R_err)*eps_0/(2*n_PML*dx*sqrt(mu_0*mu_r*eps_0*eps_r)); O,irpQ  
Bm}iU~(Z`  
    for (i = 0; i < n_PML; i++) IT&i,`cJ~F  
    { )5GQJiY  
        sigma_x         = sigma_max*pow( (n_PML - i)/((double) n_PML) ,exponent); 6pm~sD  
        ka_x            = 1 + (ka_max - 1)*pow( (n_PML-i)/((double) n_PML) ,exponent); // ka_x = 1 Q7(eq0na  
+NR n0 z(  
        K_E1_a[nx-i-1]  = (2*eps_0*ka_x - sigma_x*dt)/(2*eps_0*ka_x + sigma_x*dt);// Y&GuDLUF  
=<.F3lo\s  
        K_E1_a       = K_E1_a[nx-i-1];   // J3IRP/*z  
]EN&S
Wh  
        K_E1_b[nx-i-1]  = 2*eps_0*dt/(2*eps_0*ka_x+sigma_x*dt)/(mu_0*mu_r*dx);// FUqt)YHi  
uA%Ts*aN  
        K_E1_b       = K_E1_b[nx-i-1];    // FWY[=S  
u E.^w;~2=  
>3P9 i ;W  
        //i+0.5 ox4W$YdMG  
        sigma_x         = sigma_max*pow( (n_PML - i - 0.5)/n_PML ,exponent); 9I kUZW  
        ka_x            = 1 + (ka_max - 1)*pow( (n_PML - i - 0.5)/n_PML ,exponent); L[]BzsIv  
"@)lH  
        K_E2_a[nx-i-2]  = (2*eps_0*ka_x - sigma_x*dt)/(2*eps_0*ka_x + sigma_x*dt); VYigxhP7  
y\z > /q  
        K_E2_a       = K_E2_a[nx-i-2]; |X3">U +-  
O^NP0E  
        K_E2_b[nx-i-2]  = 2*dt/(2*eps_0*ka_x + sigma_x*dt)/(eps_r*dx);//?eps_0* mPG7Zy$z  
s.rT]  
        K_E2_b       = K_E2_b[nx-i-2]; YxyG\J\|,  
    } >UWStzH<  
} I4~^TrznRa  
]/44Ygz/  
////////////////////////////////////////////////////////////////////// (Oq
Hfv  
//Compute the Hz component a.,i.2  
////////////////////////////////////////////////////////////////////// ^1VbH3M  
void FDTD_1D_EyHz::calc_Hz_1D(double t) W>$BF[x!{  
{ GA@ Ue9  
    switch (source_type) c[:Wf<%|  
    { M=[th  
        case 1: //Gaussian pulse 
X+at%L=  
            Hz_1D[nx/2] = H0*exp( -pow( (t-t0)/tw , 2) ); [%~^kq=|  
            break; Hg whe=P  
        case 2: //Sinusoidal plane wave =YHt9fb$c  
            Hz_1D[nx/2] = H0*cos( omega*t + phi); h.^o)T  
            break; i|4_m  
        case 3: //Wave packet(sinusoidal modulated Gaussian pulse) OL9]*G?F  
            Hz_1D[nx/2] = H0*cos( omega*t + phi)*exp( -pow( (t-t0)/tw , 2) ); 5~U:@Tp  
    } U`1l8'W}:#  
//#pragma omp parallel for    这里加入并行就出现错误。可能因为系数 ]1}h8/  
    for (i = 1; i < nx - 1; i++) JY@X2'>v/  
    { t'aSF{%  
        Hz_1D = K_E1_a*Hz_1D - K_E1_b*(Ey_1D - Ey_1D[i-1]);//H振幅为0.5 N&x:K+Zm.  
        //Hz_1D = Hz_1D + (0.5)*(Ey_1D - Ey_1D[i+1]);//这样设置数据还是H振幅为0.5 ,R~eY?{a  
    } =G>.-Qfs  
<jFSj=cIL  
} PG"@A  
BSDk9Oc  
////////////////////////////////////////////////////////////////////// 1i+FL''  
//Compute the Ey component bpp*  
////////////////////////////////////////////////////////////////////// ytz8=\p_b  
void FDTD_1D_EyHz::calc_Ey_1D() ugxw!cj  
{ $T/#1w P  
//#pragma omp parallel for shared() Qi:j)uDW  
    for (i = 0; i < nx - 1; i++) Mj'lASI  
    { ttj2b$M,  
        Ey_1D = K_E2_a*Ey_1D - K_E2_b*(Hz_1D[i+1] - Hz_1D);//E振幅为0.5 #>bT<  
        //Ey_1D = Ey_1D + (0.5)*(Hz_1D[i-1] - Hz_1D);//这样设置数据还是E振幅为0.5 pL)xqKj  
        if(i%100 == 1) T.2ZBG~|[  
            cout << "Ey_1D[ i ] = " << Ey_1D <<endl; G_+Ph^  
    } y\Dn^  
} 5p )IV>G  
@'gl~J7  
}{mG/(LX8  
k/bque  
&57U? oY  
:l`i4
kx  
o8tS  
============================================================================= !7c'<[+Hm  

$}4K`Iu  
R9vY:oN%  
这段代码是核心部分，包括UPML设置和EH迭代，一维波沿x方向传输。 }P.K2ku  
PGY9*0n  
代码中的有些数组的[ i ]显示不出来 N8 M'0i?  
cAibB&`~  
=============================================================================

用户被禁言,该主题自动屏蔽!