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

图片:Graph7.jpg

rwUKg[ 1N  
Q>}*l|Ci  
X}$uvB}+>  
我写了个UPML吸收边界的一维FDTD程序，为什么结果图形显示的振幅始终为0.5， -h<Rby  
bl;C=n  
源是高斯点源。软源：  xx[i1] += J0*cos(omega*time + phase);  J0=1是振幅。 vo_m$/O  
a,sU-w!X'  
我实际默认的振幅是1.0啊，不知道为什么？ 大家有碰到过吗？ {qHf%
y&[  
_A%8oYS  
%XieKL  
.oFkx*Ln  
按照版主的提示我找了，下面是前一百个时间步源输出的结果， 看不出有什么异样。 请教大家！ qm1;^j&y  
2nGQD{  
图片是每隔100个时间步输出的Hz的结果图形。 H%:~&_D  
[EETx-  
Process : 0 -- Start FDTD $pm5G} .  
5ZkMd!$y  
xx[i1] = 0.960239    --PT_Sin _n;V iQMu  
xx[i1] = 0.962533    --PT_Sin J9 iQW  
xx[i1] = -325.327    --PT_Sin Vg)]F+E  
xx[i1] = 0.2498    --PT_Sin `kdP)lI `  
xx[i1] = -145.996    --PT_Sin ,!?&LdPt>  
xx[i1] = 0.211419    --PT_Sin ?lYi![.o  
xx[i1] = -48.8379    --PT_Sin 3,cZ*4('d  
xx[i1] = 0.401967    --PT_Sin 0oFRcU  
xx[i1] = -154.021    --PT_Sin K2glkGK  
xx[i1] = -0.178958    --PT_Sin )0I;+9:D=  
xx[i1] = 36.4055    --PT_Sin `|i[*+WC  
xx[i1] = -0.131179    --PT_Sin Htf|VpzMb  
xx[i1] = 77.6417    --PT_Sin b^Xq(q>5  
xx[i1] = -0.157616    --PT_Sin L-lDvc?5c  
xx[i1] = 42.0433    --PT_Sin sinnHQ  
xx[i1] = -0.593708    --PT_Sin b4$-?f?V  
xx[i1] = 208.648    --PT_Sin hC[MYAaF  
xx[i1] = -0.422673    --PT_Sin Qr4c':
8  
xx[i1] = 183.003    --PT_Sin nRmZu\(Ow|  
xx[i1] = -0.497458    --PT_Sin Q!<b"8V]  
xx[i1] = 165.057    --PT_Sin )J"Lne*"  
xx[i1] = -0.705738    --PT_Sin 4M&6q(389  
xx[i1] = 265.784    --PT_Sin p Rn vd|  
xx[i1] = -0.401242    --PT_Sin M#Vl{ b  
xx[i1] = 169.669    --PT_Sin jHj*S9:`  
xx[i1] = -0.447727    --PT_Sin Nn],s
Es  
xx[i1] = 147.568    --PT_Sin <-:gaA`KM  
xx[i1] = -0.423352    --PT_Sin PKhH0O\_U  
xx[i1] = 171.735    --PT_Sin @,RrAL}|  
xx[i1] = -0.052944    --PT_Sin =sk]/64h``  
xx[i1] = 31.8007    --PT_Sin `4'['x  
xx[i1] = -0.0847835    --PT_Sin \)?[1b&[_  
xx[i1] = 15.2583    --PT_Sin [C]u!\(IF  
xx[i1] = 0.0853313    --PT_Sin  d
!%:Ok  
xx[i1] = -12.5427    --PT_Sin :N<ZO`l?  
xx[i1] = 0.395852    --PT_Sin W^Jh'^E  
xx[i1] = -145.757    --PT_Sin >d&_e[j  
xx[i1] = 0.313667    --PT_Sin pbe" w=<  
xx[i1] = -129.753    --PT_Sin 4Y[uqn[  
xx[i1] = 0.502903    --PT_Sin ~97T0{E3  
xx[i1] = -168.642    --PT_Sin rQ^$)%uP  
xx[i1] = 0.636916    --PT_Sin .O
Hjn|  
xx[i1] = -246.076    --PT_Sin Vb`m3  
xx[i1] = 0.456816    --PT_Sin i-'rS/R  
xx[i1] = -179.002    --PT_Sin 0# D4;v  
xx[i1] = 0.573885    --PT_Sin JR1/\F<}  
xx[i1] = -199.4    --PT_Sin 85<zl|ZD  
xx[i1] = 0.492852    --PT_Sin b NB
pt}$  
xx[i1] = -200.143    --PT_Sin D<zgs2Ex  
xx[i1] = 0.238317    --PT_Sin (q!tI*}  
xx[i1] = -92.1195    --PT_Sin # 4AyA$t  
xx[i1] = 0.274161    --PT_Sin ugtb`d{ Sl  
xx[i1] = -93.1865    --PT_Sin $4\,a^  
xx[i1] = 0.0412927    --PT_Sin P%v7(bqL4+  
xx[i1] = -34.609    --PT_Sin *t'qn  
xx[i1] = -0.192185    --PT_Sin $~<);dYu0  
xx[i1] = 75.2772    --PT_Sin p xrd D7  
xx[i1] = -0.175713    --PT_Sin =8?gx$r2  
xx[i1] = 69.7022    --PT_Sin 8lo /BGxS>  
xx[i1] = -0.433669    --PT_Sin T=|oZ  
xx[i1] = 145.059    --PT_Sin |5 sI=?p&t  
xx[i1] = -0.532358    --PT_Sin ei>8{v&g  
xx[i1] = 209.3    --PT_Sin piE9qXn  
xx[i1] = -0.461939    --PT_Sin n"'1.  
xx[i1] = 172.704    --PT_Sin *5bKJgwJ  
xx[i1] = -0.630279    --PT_Sin _N5$>2  
xx[i1] = 224.702    --PT_Sin DbdxHuKa>  
xx[i1] = -0.535067    --PT_Sin k{3:$, b  
xx[i1] = 215.544    --PT_Sin =X?jId{  
xx[i1] = -0.390707    --PT_Sin iD)R*vnAi  
xx[i1] = 142.803    --PT_Sin C<GS._V&  
xx[i1] = -0.43058    --PT_Sin 5fK#*(x  
xx[i1] = 157.506    --PT_Sin Q]7Rqslz  
xx[i1] = -0.179445    --PT_Sin x(nWyVB  
xx[i1] = 84.0556    --PT_Sin  xA DjQ%B  
xx[i1] = -0.0186302    --PT_Sin M~Yho".  
xx[i1] = 0.282467    --PT_Sin 2-9'zN0u  
xx[i1] = 0.0216978    --PT_Sin aO9a G*9T  
xx[i1] = -5.42531    --PT_Sin }%B^Vl%ZZ  
xx[i1] = 0.311431    --PT_Sin t'bzhPQO)f  
xx[i1] = -102.723    --PT_Sin [.O3z*[9#  
xx[i1] = 0.38785    --PT_Sin $\xS~w  
xx[i1] = -155.329    --PT_Sin X0Q};,  
xx[i1] = 0.41939    --PT_Sin 72qbxPY13h  
xx[i1] = -150.535    --PT_Sin SbmakNWJ}  
xx[i1] = 0.618776    --PT_Sin !A(*?0`  
xx[i1] = -224.607    --PT_Sin DS,"^K  
xx[i1] = 0.536137    --PT_Sin 0^VA,QkQ\  
xx[i1] = -213.51    --PT_Sin ]g jhrD  
xx[i1] = 0.499646    --PT_Sin RzG<&a3B3s  
xx[i1] = -179.163    --PT_Sin )E<<
  
xx[i1] = 0.542162    --PT_Sin p.9VyM  
xx[i1] = -204.068    --PT_Sin $okGqu8z.O  
xx[i1] = 0.309125    --PT_Sin ^[{\ZX  
xx[i1] = -128.745    --PT_Sin &<Bx1\ ~V  
xx[i1] = 0.220424    --PT_Sin _VFxzM9f  
xx[i1] = -74.286    --PT_Sin 8v7;{4^  
xx[i1] = 0.1345    --PT_Sin ad).X:Qs  
xx[i1] = -58.0096    --PT_Sin w\U fq  
xx[i1] = -0.15247    --PT_Sin 2 w2JFdm  
xx[i1] = 47.4533    --PT_Sin "pb,|U
 
xx[i1] = -0.213577    --PT_Sin IG?044Y  
xx[i1] = 88.5413    --PT_Sin `Z*k M VN  
xx[i1] = -0.333515    --PT_Sin  hfpSxL  
xx[i1] = 114.074    --PT_Sin D}1Z TX_  
xx[i1] = -0.543425    --PT_Sin !JtVp&?  
xx[i1] = 200.425    --PT_Sin __\Tv>Y  
xx[i1] = -0.492682    --PT_Sin >E&mNp  
xx[i1] = 193.356    --PT_Sin \vVGfG?6  
xx[i1] = -0.555534    --PT_Sin zmH8#  
xx[i1] = 197.493    --PT_Sin \K+LKa)  
xx[i1] = -0.599966    --PT_Sin x j6-~<  
xx[i1] = 229.339    --PT_Sin i?uJ<BdU[  
xx[i1] = -0.415018    --PT_Sin Z}(,OZh  
xx[i1] = 163.493    --PT_Sin wV- kB4^4
 
xx[i1] = -0.39503    --PT_Sin 3+gp_7L  
xx[i1] = 139.644    --PT_Sin ^lCys  
xx[i1] = -0.27881    --PT_Sin
n22k<@y  
xx[i1] = 115.171    --PT_Sin P8jXruZr  
xx[i1] = -0.0255249    --PT_Sin SqL8MKN)  
xx[i1] = 14.5886    --PT_Sin u?[dy n  
xx[i1] = 0.0247965    --PT_Sin +5Yf9  
xx[i1] = -15.1405    --PT_Sin a $pxt!6  
xx[i1] = 0.212232    --PT_Sin %/17K2g  
xx[i1] = -66.3698    --PT_Sin |cwGc\ES  
xx[i1] = 0.413393    --PT_Sin :jB8Q$s  
xx[i1] = -155.249    --PT_Sin E9_aNYD
 
xx[i1] = 0.407178    --PT_Sin 9ZjSM,+  
xx[i1] = -156.635    --PT_Sin \h
pD  
xx[i1] = 0.553133    --PT_Sin U$R+&@;  
xx[i1] = -195.756    --PT_Sin ` V}e$  
xx[i1] = 0.599839    --PT_Sin : [o0Va2 d  
xx[i1] = -231.548    --PT_Sin !ydJ{\;  
xx[i1] = 0.485114    --PT_Sin lg_X|yhL  
xx[i1] = -184.438    --PT_Sin R8a4F^{*  
xx[i1] = 0.526204    --PT_Sin b&BSigrvou  
xx[i1] = -189.969    --PT_Sin kH Y  
xx[i1] = 0.398243    --PT_Sin Np>[mNmga  
xx[i1] = -161.091    --PT_Sin *Z"9QX  
xx[i1] = 0.203638    --PT_Sin }Wf\\  
xx[i1] = -76.8309    --PT_Sin &D,gKT~  
xx[i1] = 0.16127    --PT_Sin &K@2kq,  
xx[i1] = -57.8966    --PT_Sin <^?64  
xx[i1] = -0.0668099    --PT_Sin H"8fnN=xB  
xx[i1] = 11.6205    --PT_Sin N>?R,XM V  
xx[i1] = -0.242777    --PT_Sin zG|}| //}  
xx[i1] = 94.0426    --PT_Sin ngM>Tzirt  
xx[i1] = -0.286566    --PT_Sin PYDf|S7  
xx[i1] = 106.608    --PT_Sin nd(O;XBI  
xx[i1] = -0.492734    --PT_Sin uc@f#(-  
xx[i1] = 173.875    --PT_Sin 4R5+"h:  
xx[i1] = -0.543016    --PT_Sin (;2J}XQvO~  
xx[i1] = 211.046    --PT_Sin &KC^Vn3Nj  
xx[i1] = -0.511151    --PT_Sin gz6BfHQG  
xx[i1] = 189.127    --PT_Sin 9V*h:[6a(  
xx[i1] = -0.601816    --PT_Sin hC@oyC(4  
xx[i1] = 220.344    --PT_Sin WyH2` xxX  
xx[i1] = -0.482056    --PT_Sin z}v6!u|iZu  
xx[i1] = 191.811    --PT_Sin FCv3ZF?K  
xx[i1] = -0.363158    --PT_Sin ^v.~FFK  
xx[i1] = 132.74    --PT_Sin 5#+G7 'k  
xx[i1] = -0.327607    --PT_Sin _ F&BSu  
xx[i1] = 123.647    --PT_Sin `z<k7ig  
xx[i1] = -0.0890027    --PT_Sin 0\f3La  
xx[i1] = 45.2235    --PT_Sin p!<Y 'G  
xx[i1] = 0.0493984    --PT_Sin (Hqy^EOZ  
xx[i1] = -23.2674    --PT_Sin )En*5-1  
xx[i1] = 0.141235    --PT_Sin $ BV4i$  
xx[i1] = -47.8013    --PT_Sin E"7 iU  
xx[i1] = 0.380387    --PT_Sin |,c\R"8xS  
xx[i1] = -133.801    --PT_Sin FBpf_=(_1  
xx[i1] = 0.435991    --PT_Sin y=)Cid  
xx[i1] = -170.321    --PT_Sin 1*aw~nY
0  
xx[i1] = 0.489207    --PT_Sin ^<fN  
xx[i1] = -176.668    --PT_Sin w8M,35b  
xx[i1] = 0.615248    --PT_Sin .b*%c
?e  
xx[i1] = -227.859    --PT_Sin T0xU}  
xx[i1] = 0.522702    --PT_Sin 8KKI.i8`  
xx[i1] = -204.76    --PT_Sin swv1>52{  
xx[i1] = 0.487312    --PT_Sin {e4`D1B  
xx[i1] = -176.404    --PT_Sin :4]^PB@dl  
xx[i1] = 0.459206    --PT_Sin +~lZ]a7k  
xx[i1] = -175.971    --PT_Sin
%k(V 2]WF  
epa)~/sA  

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

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

谢谢版主！

TA x
9<'  
y([""z3<w  
顶起！

你把你的代码贴出来看看

Kn#3^>D  
/////////////////////////////////////////////////////////////////////////////////////// vYYLn9}5  
//Set the PML matrices 1pUIZ$@?`  
/////////////////////////////////////////////////////////////////////////////////////// k 5<[N2D|!  
void FDTD_1D_EyHz::Set_PML_Param_1D() Is87 9_Z  
{ ]^>#?yEA3  
    for (i = 0; i < nx; i++) ";)SA,Z  
    { Ik92='Z  
        K_E1_a = 1.0;                // j3!]wolY  
        K_E1_b = dt/(mu_0*mu_r*dx);            // l c<&f  
F$Cf\#{3  
        if (i < nx-1) =,*4:
TU  
        { jslfq@5v  
            K_E2_a = 1.0;            // ?pT\Ft V  
            K_E2_b = dt/(eps_0*eps_r*dx);    // dt/(eps_r*dx);     5y}kI  
        }             64#6L.Q-c  
    } DwQaj"1<%  
[n53eC  
    //PML parameters SqiLp!Y`  
    double ka_max = 1; g{`rWKj  
    int exponent = 4; Ly0U')D:  
    double R_err = 1e-16; v{R:F  
    double eta = sqrt(mu_0*mu_r/eps_0/eps_r); )u. ut8![T  
    double sigma_x, sigma_max, ka_x; x8]9Xe:_>O  
         `=]I-5#.W  
    sigma_max= -(exponent+1)*log(R_err)/(2*eta*n_PML*dx); D]I]I!2c  
    //sigma_max= -(exponent+1)*log(R_err)*eps_0/(2*n_PML*dx*sqrt(mu_0*mu_r*eps_0*eps_r)); P5:X7[  
\2rC
T~x  
    for (i = 0; i < n_PML; i++) ,W'?F9Y\  
    { vFKt=o$ g  
        sigma_x         = sigma_max*pow( (n_PML - i)/((double) n_PML) ,exponent); . PzlhTL7  
        ka_x            = 1 + (ka_max - 1)*pow( (n_PML-i)/((double) n_PML) ,exponent); // ka_x = 1 Vd8BQB,Q  
&DqeO8?Q  
        K_E1_a[nx-i-1]  = (2*eps_0*ka_x - sigma_x*dt)/(2*eps_0*ka_x + sigma_x*dt);// }ph;~og}y  
~}SOd<n)|  
        K_E1_a       = K_E1_a[nx-i-1];   // nUAs:Q  
I+) Acy;  
        K_E1_b[nx-i-1]  = 2*eps_0*dt/(2*eps_0*ka_x+sigma_x*dt)/(mu_0*mu_r*dx);// N"i'[!H%  
>3S^9{d  
        K_E1_b       = K_E1_b[nx-i-1];    // (z[|\6O  
bS0z\!1  
Jy,Dcl  
        //i+0.5 4_`ss+gk  
        sigma_x         = sigma_max*pow( (n_PML - i - 0.5)/n_PML ,exponent); y=L9E?  
        ka_x            = 1 + (ka_max - 1)*pow( (n_PML - i - 0.5)/n_PML ,exponent); ([-xM%BI6  
(IbT5  
        K_E2_a[nx-i-2]  = (2*eps_0*ka_x - sigma_x*dt)/(2*eps_0*ka_x + sigma_x*dt); `.0QY<;  
]FJpe^ ua  
        K_E2_a       = K_E2_a[nx-i-2]; t^|+|>S  
=uZOpeviQ  
        K_E2_b[nx-i-2]  = 2*dt/(2*eps_0*ka_x + sigma_x*dt)/(eps_r*dx);//?eps_0*  c`'2  
zuWfR&U|W  
        K_E2_b       = K_E2_b[nx-i-2]; lgxG:zAC  
    } HWr")%EhD  
} %N"9'g>  
}y-b<J?H  
////////////////////////////////////////////////////////////////////// N6GvzmG#g  
//Compute the Hz component ]LvpYRU$P  
////////////////////////////////////////////////////////////////////// hL\gI(B  
void FDTD_1D_EyHz::calc_Hz_1D(double t) k5>K/;*9  
{ ^JDiI7  
    switch (source_type) 4eSV(u)4  
    { +["t@Q4IQ  
        case 1: //Gaussian pulse Wr\rruH6  
            Hz_1D[nx/2] = H0*exp( -pow( (t-t0)/tw , 2) ); x)X=sX.  
            break; 0n<t/74  
        case 2: //Sinusoidal plane wave %"R|tlG  
            Hz_1D[nx/2] = H0*cos( omega*t + phi); x<  Td  
            break; e3+'m  
        case 3: //Wave packet(sinusoidal modulated Gaussian pulse) MDfE(cn2q  
            Hz_1D[nx/2] = H0*cos( omega*t + phi)*exp( -pow( (t-t0)/tw , 2) ); &^HqbLz  
    } \78w1Rkl  
//#pragma omp parallel for    这里加入并行就出现错误。可能因为系数 w$JG:y#  
    for (i = 1; i < nx - 1; i++) 5WR(jl+M  
    { @wWro?s'p  
        Hz_1D = K_E1_a*Hz_1D - K_E1_b*(Ey_1D - Ey_1D[i-1]);//H振幅为0.5 ^9})@,(D  
        //Hz_1D = Hz_1D + (0.5)*(Ey_1D - Ey_1D[i+1]);//这样设置数据还是H振幅为0.5 _084GK9{W  
    } FW)G5^Tf  
=6mnXpM.  
} zc QFIP  
rb+&]  
////////////////////////////////////////////////////////////////////// vW\#2[j[  
//Compute the Ey component E7eOKNVC#  
////////////////////////////////////////////////////////////////////// ''9FB5
 
void FDTD_1D_EyHz::calc_Ey_1D() B}= WxG|)  
{ (6C%w)8'  
//#pragma omp parallel for shared() ")l_>y?  
    for (i = 0; i < nx - 1; i++) S`LS/)  
    { ]jUxL=]r  
        Ey_1D = K_E2_a*Ey_1D - K_E2_b*(Hz_1D[i+1] - Hz_1D);//E振幅为0.5 p;S<WJv k  
        //Ey_1D = Ey_1D + (0.5)*(Hz_1D[i-1] - Hz_1D);//这样设置数据还是E振幅为0.5 t83n`LC  
        if(i%100 == 1) Rt$Q*`u   
            cout << "Ey_1D[ i ] = " << Ey_1D <<endl; G22=8V  
    } 3{mu77  
} N'L3Oa\%  
>Q0HqOq  
{DPobyvwFk  
J%?
'Q{
  
D`;Q?fC  
V 0<>Xo%  
g-c\;  
============================================================================= ~E2xIhV  
H];QDix?  
EJ*  
这段代码是核心部分，包括UPML设置和EH迭代，一维波沿x方向传输。 8
3*"58  
bg\~"  
代码中的有些数组的[ i ]显示不出来 8ODrW!o  
Qrt[MJ+#  
=============================================================================

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