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

图片:Graph7.jpg

6G(k{S  
nK:39D$(  
GJ*AyYG  
我写了个UPML吸收边界的一维FDTD程序，为什么结果图形显示的振幅始终为0.5， hP7nt  
A.y$.(  
源是高斯点源。软源：  xx[i1] += J0*cos(omega*time + phase);  J0=1是振幅。 ZQyT$l~b  
Y`M.hYBXk  
我实际默认的振幅是1.0啊，不知道为什么？ 大家有碰到过吗？ BjB2YO& /  
{_ #   

9`b*Y*d  
9oaq%Sf  
按照版主的提示我找了，下面是前一百个时间步源输出的结果， 看不出有什么异样。 请教大家！ [X-Q{c4  
N%T-Q9k  
图片是每隔100个时间步输出的Hz的结果图形。 >V]>h&`  
bDr'W  
Process : 0 -- Start FDTD >Mn"k\j4  
r2
Q"NVw  
xx[i1] = 0.960239    --PT_Sin ALKhZFuz  
xx[i1] = 0.962533    --PT_Sin F`-? 3]\3  
xx[i1] = -325.327    --PT_Sin I($u L@$  
xx[i1] = 0.2498    --PT_Sin C6Kz6_DQZ  
xx[i1] = -145.996    --PT_Sin
7Fq|Zc`P  
xx[i1] = 0.211419    --PT_Sin 0]"j,  
xx[i1] = -48.8379    --PT_Sin {!-w|&bF  
xx[i1] = 0.401967    --PT_Sin 9)=as/o  
xx[i1] = -154.021    --PT_Sin >6aCBS?2  
xx[i1] = -0.178958    --PT_Sin _ p?q/-[4  
xx[i1] = 36.4055    --PT_Sin E&97;VH  
xx[i1] = -0.131179    --PT_Sin ] $%{nj<  
xx[i1] = 77.6417    --PT_Sin ?56Zw"89  
xx[i1] = -0.157616    --PT_Sin vbSz&+52;  
xx[i1] = 42.0433    --PT_Sin xd>2TW l#  
xx[i1] = -0.593708    --PT_Sin H(0d(c1s  
xx[i1] = 208.648    --PT_Sin [*1c.&%(  
xx[i1] = -0.422673    --PT_Sin '- Z4GcL  
xx[i1] = 183.003    --PT_Sin ~:JKXa?  
xx[i1] = -0.497458    --PT_Sin Y&#<{j':  
xx[i1] = 165.057    --PT_Sin  1 ft.ZJ  
xx[i1] = -0.705738    --PT_Sin B aXzz  
xx[i1] = 265.784    --PT_Sin Y(
&phv&  
xx[i1] = -0.401242    --PT_Sin x.d9mjLN8m  
xx[i1] = 169.669    --PT_Sin }$b/g  
xx[i1] = -0.447727    --PT_Sin N%^mR>.`  
xx[i1] = 147.568    --PT_Sin :?60pu=  
xx[i1] = -0.423352    --PT_Sin 87*R#((  
xx[i1] = 171.735    --PT_Sin @]cpPW-b  
xx[i1] = -0.052944    --PT_Sin au GN~"n^  
xx[i1] = 31.8007    --PT_Sin 
B[k"xs  
xx[i1] = -0.0847835    --PT_Sin {2V=BDS|?K  
xx[i1] = 15.2583    --PT_Sin AKS(WNGEp  
xx[i1] = 0.0853313    --PT_Sin "Uyw
7  
xx[i1] = -12.5427    --PT_Sin 4ba1c  
xx[i1] = 0.395852    --PT_Sin K[ylyQ1  
xx[i1] = -145.757    --PT_Sin !{SEm"J^  
xx[i1] = 0.313667    --PT_Sin ']qC,;2  
xx[i1] = -129.753    --PT_Sin _/KW5  
xx[i1] = 0.502903    --PT_Sin {xOu*8J  
xx[i1] = -168.642    --PT_Sin $+?6U  
xx[i1] = 0.636916    --PT_Sin C#gQJ=!B  
xx[i1] = -246.076    --PT_Sin ag] nVE/  
xx[i1] = 0.456816    --PT_Sin ntjUnd&v\  
xx[i1] = -179.002    --PT_Sin #M_QSD}&  
xx[i1] = 0.573885    --PT_Sin /1O6;'8He  
xx[i1] = -199.4    --PT_Sin hwexv 9""  
xx[i1] = 0.492852    --PT_Sin K<V(h#(.@  
xx[i1] = -200.143    --PT_Sin u52@{@Ad  
xx[i1] = 0.238317    --PT_Sin %';n9M  
xx[i1] = -92.1195    --PT_Sin s cn!,  
xx[i1] = 0.274161    --PT_Sin & ??)gMM[  
xx[i1] = -93.1865    --PT_Sin z`TI<B
 
xx[i1] = 0.0412927    --PT_Sin yLI=&7/e@  
xx[i1] = -34.609    --PT_Sin {8t;nsdm!  
xx[i1] = -0.192185    --PT_Sin %.Mtn%:I*  
xx[i1] = 75.2772    --PT_Sin >f_D|;EV  
xx[i1] = -0.175713    --PT_Sin A^g81s.5  
xx[i1] = 69.7022    --PT_Sin 9%)'QDVGLf  
xx[i1] = -0.433669    --PT_Sin >(\[$  
xx[i1] = 145.059    --PT_Sin )`g[k"yB3  
xx[i1] = -0.532358    --PT_Sin S46[2-v1  
xx[i1] = 209.3    --PT_Sin X-t4irZ)  
xx[i1] = -0.461939    --PT_Sin Ir]b.6B  
xx[i1] = 172.704    --PT_Sin hR. EZ|.  
xx[i1] = -0.630279    --PT_Sin /0(4wZe~?  
xx[i1] = 224.702    --PT_Sin C@KYg/nYw  
xx[i1] = -0.535067    --PT_Sin PY`V]|J  
xx[i1] = 215.544    --PT_Sin S?D2`b  
xx[i1] = -0.390707    --PT_Sin E6n;_{Se/S  
xx[i1] = 142.803    --PT_Sin 0V1kZ.  
xx[i1] = -0.43058    --PT_Sin $bMeL7CN  
xx[i1] = 157.506    --PT_Sin v}iJ:'  
xx[i1] = -0.179445    --PT_Sin xNjA>S\]W5  
xx[i1] = 84.0556    --PT_Sin *aTM3k)Zs  
xx[i1] = -0.0186302    --PT_Sin +[*UC"  
xx[i1] = 0.282467    --PT_Sin (L~3nN;rr  
xx[i1] = 0.0216978    --PT_Sin $-o39A#  
xx[i1] = -5.42531    --PT_Sin ]H.+=V;1  
xx[i1] = 0.311431    --PT_Sin 8_KXli}7=  
xx[i1] = -102.723    --PT_Sin (spX3n%p  
xx[i1] = 0.38785    --PT_Sin Jq.26I=  
xx[i1] = -155.329    --PT_Sin GifD>c |z  
xx[i1] = 0.41939    --PT_Sin Ju:=-5r"'  
xx[i1] = -150.535    --PT_Sin 0`OqD d  
xx[i1] = 0.618776    --PT_Sin bG\1<:6B  
xx[i1] = -224.607    --PT_Sin jkfI,T  
xx[i1] = 0.536137    --PT_Sin =lV
frna  
xx[i1] = -213.51    --PT_Sin DrW]`%Ql  
xx[i1] = 0.499646    --PT_Sin V @8X.R>  
xx[i1] = -179.163    --PT_Sin 1P6~IZVN  
xx[i1] = 0.542162    --PT_Sin Th"7p:SE?  
xx[i1] = -204.068    --PT_Sin @ cv`}k  
xx[i1] = 0.309125    --PT_Sin Wmp\J3  
xx[i1] = -128.745    --PT_Sin
SLBKXj|  
xx[i1] = 0.220424    --PT_Sin FmnA+fA  
xx[i1] = -74.286    --PT_Sin J\2F%kBej?  
xx[i1] = 0.1345    --PT_Sin 4,)=r3;&!  
xx[i1] = -58.0096    --PT_Sin W}(dhgf  
xx[i1] = -0.15247    --PT_Sin b"x:IDW qG  
xx[i1] = 47.4533    --PT_Sin K<N0%c~  
xx[i1] = -0.213577    --PT_Sin M`"2;  
xx[i1] = 88.5413    --PT_Sin :Z&ipd!yY  
xx[i1] = -0.333515    --PT_Sin +LrW#K;  
xx[i1] = 114.074    --PT_Sin w$2-t  
xx[i1] = -0.543425    --PT_Sin {\
.2h  
xx[i1] = 200.425    --PT_Sin hf%W grO.  
xx[i1] = -0.492682    --PT_Sin fV[xv4D.  
xx[i1] = 193.356    --PT_Sin p-ry{"XA  
xx[i1] = -0.555534    --PT_Sin /ZD/!YD&R  
xx[i1] = 197.493    --PT_Sin 4N*^%  
xx[i1] = -0.599966    --PT_Sin 1f~_# EIC  
xx[i1] = 229.339    --PT_Sin  f0:)  
xx[i1] = -0.415018    --PT_Sin C?b Mj[$  
xx[i1] = 163.493    --PT_Sin hI/p9 `w  
xx[i1] = -0.39503    --PT_Sin dU+1@_  
xx[i1] = 139.644    --PT_Sin {x-
g?HB
 
xx[i1] = -0.27881    --PT_Sin L9G=+T9  
xx[i1] = 115.171    --PT_Sin 1tg   
xx[i1] = -0.0255249    --PT_Sin n NAJ8z}Nt  
xx[i1] = 14.5886    --PT_Sin .p0;y3so4  
xx[i1] = 0.0247965    --PT_Sin `M\L6o  
xx[i1] = -15.1405    --PT_Sin \*_qP*vq@  
xx[i1] = 0.212232    --PT_Sin f!1KGP  
xx[i1] = -66.3698    --PT_Sin 4;%=ohD:!  
xx[i1] = 0.413393    --PT_Sin v^KJU +  
xx[i1] = -155.249    --PT_Sin 64zO%F*  
xx[i1] = 0.407178    --PT_Sin ?t<wp3bZ  
xx[i1] = -156.635    --PT_Sin k ^+h>B-;  
xx[i1] = 0.553133    --PT_Sin d[ {=/~0  
xx[i1] = -195.756    --PT_Sin $|AvT;4  
xx[i1] = 0.599839    --PT_Sin I|BLAm6j  
xx[i1] = -231.548    --PT_Sin P^&+ehp  
xx[i1] = 0.485114    --PT_Sin ^gv)[  
xx[i1] = -184.438    --PT_Sin *r(iegO$  
xx[i1] = 0.526204    --PT_Sin E4 JS   
xx[i1] = -189.969    --PT_Sin SR8[ 7MU  
xx[i1] = 0.398243    --PT_Sin ~~h9yvW7&  
xx[i1] = -161.091    --PT_Sin igz&7U8gg  
xx[i1] = 0.203638    --PT_Sin {'{ssCL  
xx[i1] = -76.8309    --PT_Sin *6k (xL  
xx[i1] = 0.16127    --PT_Sin "zm.jNn  
xx[i1] = -57.8966    --PT_Sin 6`EyzB%.$  
xx[i1] = -0.0668099    --PT_Sin .llAiv  
xx[i1] = 11.6205    --PT_Sin [;};qQ-C2  
xx[i1] = -0.242777    --PT_Sin s;$ eq);  
xx[i1] = 94.0426    --PT_Sin Hjlx,:'M  
xx[i1] = -0.286566    --PT_Sin mB_ba1r  
xx[i1] = 106.608    --PT_Sin LG51e7_gFi  
xx[i1] = -0.492734    --PT_Sin $z` jR*  
xx[i1] = 173.875    --PT_Sin %f?#) 01>  
xx[i1] = -0.543016    --PT_Sin @ /c{gD  
xx[i1] = 211.046    --PT_Sin {K:/(\  
xx[i1] = -0.511151    --PT_Sin c*LnLK/m  
xx[i1] = 189.127    --PT_Sin 7rsrC  
xx[i1] = -0.601816    --PT_Sin $k}+,tHtJO  
xx[i1] = 220.344    --PT_Sin
I /RvU,  
xx[i1] = -0.482056    --PT_Sin i"_JF-IbN  
xx[i1] = 191.811    --PT_Sin G.cs-f  
xx[i1] = -0.363158    --PT_Sin MJ>(HJY6?%  
xx[i1] = 132.74    --PT_Sin f61
~%@fE  
xx[i1] = -0.327607    --PT_Sin 4?8GK  
xx[i1] = 123.647    --PT_Sin 40+E#z)  
xx[i1] = -0.0890027    --PT_Sin k%uRG_  
xx[i1] = 45.2235    --PT_Sin Vd|/]Zj  
xx[i1] = 0.0493984    --PT_Sin ,![C8il,  
xx[i1] = -23.2674    --PT_Sin 1eKJ46W  
xx[i1] = 0.141235    --PT_Sin VRMlr.T+  
xx[i1] = -47.8013    --PT_Sin sd>#Hn
 
xx[i1] = 0.380387    --PT_Sin !$Mv)c/_u  
xx[i1] = -133.801    --PT_Sin hydn" 9;  
xx[i1] = 0.435991    --PT_Sin $YL}rM  
xx[i1] = -170.321    --PT_Sin .8g&V|  
xx[i1] = 0.489207    --PT_Sin @-Gf+*GZys  
xx[i1] = -176.668    --PT_Sin XMuZ
'I  
xx[i1] = 0.615248    --PT_Sin nj)M$'  
xx[i1] = -227.859    --PT_Sin P`RM"'Om  
xx[i1] = 0.522702    --PT_Sin \#~~,k 6f  
xx[i1] = -204.76    --PT_Sin {jr>Z"/q  
xx[i1] = 0.487312    --PT_Sin P]IN
YH  
xx[i1] = -176.404    --PT_Sin H2|w  
xx[i1] = 0.459206    --PT_Sin R-Uj\M>  
xx[i1] = -175.971    --PT_Sin
n j1 cqh  
7dxY07yu  

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

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

谢谢版主！

%oJ_,m_(  
)+'FTz` c  
顶起！

你把你的代码贴出来看看

 ntK#7(U'  
/////////////////////////////////////////////////////////////////////////////////////// 8s^CE[TA  
//Set the PML matrices Awy-kou[C  
/////////////////////////////////////////////////////////////////////////////////////// ]U?)_P@}  
void FDTD_1D_EyHz::Set_PML_Param_1D() Jse;@K5y  
{ X <QSi
 
    for (i = 0; i < nx; i++) z !2-
U  
    { |&!04~s;E  
        K_E1_a = 1.0;                // 2IDN?Mw  
        K_E1_b = dt/(mu_0*mu_r*dx);            // ]
'+PJdA  
H?A&P4nZ  
        if (i < nx-1) r:.3P  
        { [HNWM/ff7+  
            K_E2_a = 1.0;            // m?wPZ^u  
            K_E2_b = dt/(eps_0*eps_r*dx);    // dt/(eps_r*dx);     kYMKVR  
        }             8NU<lV`  
    } l`"i'P   
eb%`ox@&  
    //PML parameters MtWzGE=?  
    double ka_max = 1; EMK>7 aks  
    int exponent = 4; www#.D%'U  
    double R_err = 1e-16; vw(X9xa  
    double eta = sqrt(mu_0*mu_r/eps_0/eps_r); ffDh0mDN  
    double sigma_x, sigma_max, ka_x; YhQ;>Ko  
         nM]Sb|1:  
    sigma_max= -(exponent+1)*log(R_err)/(2*eta*n_PML*dx); +$_.${uwV  
    //sigma_max= -(exponent+1)*log(R_err)*eps_0/(2*n_PML*dx*sqrt(mu_0*mu_r*eps_0*eps_r)); Y.FqWJP=p  
J" :R,w`  
    for (i = 0; i < n_PML; i++) g>].m8DZ'  
    { jFAnhbbCE  
        sigma_x         = sigma_max*pow( (n_PML - i)/((double) n_PML) ,exponent); B'sgCU  
        ka_x            = 1 + (ka_max - 1)*pow( (n_PML-i)/((double) n_PML) ,exponent); // ka_x = 1 Am>^{qh9  
e\o>(is  
        K_E1_a[nx-i-1]  = (2*eps_0*ka_x - sigma_x*dt)/(2*eps_0*ka_x + sigma_x*dt);// MC=pN(l  
$646"1S  
        K_E1_a       = K_E1_a[nx-i-1];   // 4R<bfZ43  
a@+n  
        K_E1_b[nx-i-1]  = 2*eps_0*dt/(2*eps_0*ka_x+sigma_x*dt)/(mu_0*mu_r*dx);// i b$2qy  
pYXusS7S  
        K_E1_b       = K_E1_b[nx-i-1];    // qkHdr
2  
NHI(}Ea|]  
H$G`e'`OZ  
        //i+0.5 woI.
1e5  
        sigma_x         = sigma_max*pow( (n_PML - i - 0.5)/n_PML ,exponent); fSR+~Vy  
        ka_x            = 1 + (ka_max - 1)*pow( (n_PML - i - 0.5)/n_PML ,exponent); qRk<1.  
^xz*%2@  
        K_E2_a[nx-i-2]  = (2*eps_0*ka_x - sigma_x*dt)/(2*eps_0*ka_x + sigma_x*dt); FZdZGK  
!mX-g]4E  
        K_E2_a       = K_E2_a[nx-i-2]; D=ZH? d  
z''ITX)oG  
        K_E2_b[nx-i-2]  = 2*dt/(2*eps_0*ka_x + sigma_x*dt)/(eps_r*dx);//?eps_0* d#l z^Ls2  
  %4  
        K_E2_b       = K_E2_b[nx-i-2]; =]Gw9sge@  
    } VR:4|_o  
} xcf`i
:\  
Xx)PyO  
////////////////////////////////////////////////////////////////////// 4o|-v  
//Compute the Hz component 60RYw9d%0  
////////////////////////////////////////////////////////////////////// c-s A?q#|  
void FDTD_1D_EyHz::calc_Hz_1D(double t) QY)hMo=|o8  
{ O:G5n 5J  
    switch (source_type) }wY6^JF  
    { }fqz8'E9  
        case 1: //Gaussian pulse ! J7ExfEA  
            Hz_1D[nx/2] = H0*exp( -pow( (t-t0)/tw , 2) ); g$CWGB*%lm  
            break; "CH3\O\  
        case 2: //Sinusoidal plane wave hW~XE{<  
            Hz_1D[nx/2] = H0*cos( omega*t + phi); xMOq/")  
            break; TAd~#jB9  
        case 3: //Wave packet(sinusoidal modulated Gaussian pulse) o=%pR|  
            Hz_1D[nx/2] = H0*cos( omega*t + phi)*exp( -pow( (t-t0)/tw , 2) ); 3kU4?D]  
    } ^"  
//#pragma omp parallel for    这里加入并行就出现错误。可能因为系数 S!Z2aFj  
    for (i = 1; i < nx - 1; i++) 0-VC$)S  
    { vZ0K1UTEXY  
        Hz_1D = K_E1_a*Hz_1D - K_E1_b*(Ey_1D - Ey_1D[i-1]);//H振幅为0.5 LN!e_b  
        //Hz_1D = Hz_1D + (0.5)*(Ey_1D - Ey_1D[i+1]);//这样设置数据还是H振幅为0.5 = ^NTHc^*  
    } cJ^:b4j  
V$OZC;4  
} 4nvi7  
tV'>9YVdG  
////////////////////////////////////////////////////////////////////// VyF|d?b  
//Compute the Ey component -9f+O^x  
////////////////////////////////////////////////////////////////////// }={TVs^  
void FDTD_1D_EyHz::calc_Ey_1D() SO!|wag$  
{ _zuX6DO  
//#pragma omp parallel for shared() %qI.Qw$  
    for (i = 0; i < nx - 1; i++) C*C;n4AT  
    { V DN@=/  
        Ey_1D = K_E2_a*Ey_1D - K_E2_b*(Hz_1D[i+1] - Hz_1D);//E振幅为0.5 q eW{Cl~  
        //Ey_1D = Ey_1D + (0.5)*(Hz_1D[i-1] - Hz_1D);//这样设置数据还是E振幅为0.5 ;1gWz  
        if(i%100 == 1) 3 *g>kRMJ  
            cout << "Ey_1D[ i ] = " << Ey_1D <<endl; ()\=(n!J  
    } :_
0"t-  
} V{D~e0i/v  
bxXNv^  
)0\"8}!  
JA<Hm.V#  
A ,$CYLj+  
ZYMacTeJjg  
1Uy'TEk  
============================================================================= 78u9>
H  
D# Gf.c  
:"im2J  
这段代码是核心部分，包括UPML设置和EH迭代，一维波沿x方向传输。 Zb? u'Vm=u  
$F#eD0|  
代码中的有些数组的[ i ]显示不出来 YU]|N'mL2  
X`s6lV%\  
=============================================================================

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