登录 | 注册
开启辅助访问
设为首页加入收藏关于我们

搜索附件  
头雁微网 › 附件中心 › 技术应用 › 情报信息 › 程序.txt
板块导航
附件中心&附件聚合2.0
For Discuz! X3.5 © hgcad.com

程序.txt

 

在FDTD计算时出现反射是，应怎样调试程序:
        在用UPML吸收边界条件，出现了反射应怎样调试程序？一般的方法是什么？
学习ing！！！
从简单的复杂，先检查FDTD程序有没有问题，可以选用简单的Mur吸收条件，调试程序；FDTD程序经调试正确后，吸收边界条件改为UPML边界，此时出问题自然是UPML处了，使用debug逐步调试，结合场分布分析检测出错位置。
非常感谢版主的指导，我会逐步调试的
江海客 发表于 2012-11-28 12:33

                               
登录/注册后可看大图

非常感谢版主的指导，我会逐步调试的

不客气，有问题欢迎在此版区讨论！
在学习用UPML作为吸收边界条件计算时，在微网上下了一段经典的程序（3D的UPML程序），此程序上是（绝缘介质——UPML的情况），我将程序改写成了导电介质——UPML的情况后，发现程序的迭代出现了问题。在调试的时候我将电导率（sigma取0），理论上应该和经典程序上的一样，但实际的程序上是无法迭代的
我将我改写的程序一并传上
三维UPML的FDTD公式：导电介质——UPML情形
我觉得是程序中电场的迭代下面的迭代出现了问题：
ex(:,2:je_tot,2:ke_tot)=C3ex(:,2:je_tot,2:ke_tot).*ex(:,2:je_tot,2:ke_tot)+...
                            C41ex(:,2:je_tot,2:ke_tot).*(C5ex(:,2:je_tot,2:ke_tot).*px(:,2:je_tot,2:ke_tot)-...
                            C6ex(:,2:je_tot,2:ke_tot).*pstore(:,2:je_tot,2:ke_tot));
但这个和葛老师书上的公式一致，我就不知道该如何处理这个问题
tingo 发表于 2012-12-4 22:09

                               
登录/注册后可看大图

不客气，有问题欢迎在此版区讨论！

在学习用UPML作为吸收边界条件计算时，在微网上下了一段经典的程序（3D的UPML程序），此程序上是（绝缘介质——UPML的情况），我将程序改写成了导电介质——UPML的情况后，发现程序的迭代出现了问题。在调试的时候我将电导率（sigma取0），理论上应该和经典程序上的一样，但实际的程序上是无法迭代的
tingo 发表于 2012-12-4 22:09

                               
登录/注册后可看大图

不客气，有问题欢迎在此版区讨论！

我将我改写的程序一并传上
tingo 发表于 2012-12-4 22:09

                               
登录/注册后可看大图

不客气，有问题欢迎在此版区讨论！

三维UPML的FDTD公式：导电介质——UPML情形
tingo 发表于 2012-12-4 22:09

                               
登录/注册后可看大图

不客气，有问题欢迎在此版区讨论！
  1. %--------------------------------------------------
  2. clear;
  3. clc;
  4. close all;
  5. %------------------------------------------------
  6. % Fundamental  constants
  7. %-----------------------------------------------
  8. cc=3.0e8;                            %(cc 光波的传播速度)
  9. muz=4.0*pi*1.0e-7;                   %(muz 真空磁导率)
  10. epsz=1.0/(cc*cc*muz);                %(真空介电常数)
  11. etaz=sqrt(muz/epsz);
  12. %----------------------------------------------
  13. % Material  parameters
  14. %-----------------------------------------------

  16. mur(1)=1.0;
  17. epsr(1)=1.0;
  18. eta=etaz*sqrt(mur(1)/epsr(1));
  19. %-------------------------------------------------
  20. %   Grid parameters
  21. %-------------------------------------------------


  24. ie=137;                                  % x 方向的网格
  25. je=80;                                  % y 方向的网格
  26. ke=20;                                  % z 方向的网格

  28. ih=ie+1;                                % x 方向的边数（节点）
  29. jh=je+1;
  30. kh=ke+1;
  31. upml=10;                                %  UPML层的厚度
  32. ih_bc=upml+1;
  33. jh_bc=upml+1;
  34. kh_bc=upml+1;
  35. ie_tot=ie+2*upml;                       % 总计算区域网格的数目
  36. je_tot=je+2*upml;
  37. ke_tot=ke+2*upml;
  38. ih_tot=ie_tot+1;
  39. jh_tot=je_tot+1;
  40. kh_tot=ke_tot+1;

  42. is=round(ie_tot/2);
  43. js=round(je_tot/2);
  44. ks=round(ke_tot/2);
  45.    
  46. %---------------------------------------------
  47. %    Fundamental grid parameters
  48. %----------------------------------------------
  49. delta=0.01;
  50. % dt=delta*sqrt(epsr*mur)/(2.0*cc);
  51. dt=delta/(2.0*cc);
  52. nmax=50;

  54. %-----------------------------------------------
  55. %  Wave excitation
  56. %-----------------------------------------------
  57. rtau=50.0e-12;
  58. tau=rtau/dt;
  59. ndelay=3*tau;
  60. J0=-1.0*epsz.*10^20;
  61. %-----------------------------------------------
  62. % Initialize field arrays
  63. %------------------------------------------------
  64. ex=zeros(ie_tot, jh_tot, kh_tot);          % the field in main grid
  65. ey=zeros(ih_tot, je_tot, kh_tot);
  66. ez=zeros(ih_tot, jh_tot, ke_tot);

  68. dx=zeros(ie_tot, jh_tot, kh_tot);          % 绝缘介质——UPML中的中间变量D   
  69. dy=zeros(ih_tot, je_tot, kh_tot);
  70. dz=zeros(ih_tot, jh_tot, ke_tot);

  72. px=zeros(ie_tot, jh_tot, kh_tot);          % 导电介质——UPML中的中间变量P'和P，其中p'用表示 （这里是我加上去的）     
  73. py=zeros(ih_tot, je_tot, kh_tot);
  74. pz=zeros(ih_tot, jh_tot, ke_tot);
  75. qx=zeros(ie_tot, jh_tot, kh_tot);      
  76. qy=zeros(ih_tot, je_tot, kh_tot);
  77. qz=zeros(ih_tot, jh_tot, ke_tot);


  80. hx=zeros(ih_tot, je_tot, ke_tot);          % to magnetic field
  81. hy=zeros(ie_tot, jh_tot, ke_tot);
  82. hz=zeros(ie_tot, je_tot, kh_tot);
  83. bx=zeros(ih_tot, je_tot, ke_tot);
  84. by=zeros(ie_tot, jh_tot, ke_tot);
  85. bz=zeros(ie_tot, je_tot, kh_tot);
  86. %------------------------------------------------
  87. %  Initialize update coefficient arays
  88. %--------------------------------------------------
  89. C1ex=zeros(size(ex));
  90. C2ex=zeros(size(ex));
  91. C21ex=zeros(size(ex));                      % 这是导电介质——UPML中的迭代系数（葛德彪老师书上P105页的迭代公式（4-9-16）中的系数）
  92. C3ex=zeros(size(ex));
  93. C4ex=zeros(size(ex));
  94. C41ex=zeros(size(ex));                      % 这是导电介质——UPML中的迭代系数 （葛德彪老师书上P105页的迭代公式（4-9-16）中的系数）                  
  95. C5ex=zeros(size(ex));
  96. C6ex=zeros(size(ex));
  97. C7ex=zeros(size(ex));                       % 这是导电介质——UPML中的迭代系数 （葛德彪老师书上P105页的迭代公式（4-9-13）中的系数）
  98. C8ex=zeros(size(ex));

  100. C1ey=zeros(size(ey));
  101. C2ey=zeros(size(ey));
  102. C21ey=zeros(size(ey));
  103. C3ey=zeros(size(ey));
  104. C4ey=zeros(size(ey));
  105. C41ey=zeros(size(ey));
  106. C5ey=zeros(size(ey));
  107. C6ey=zeros(size(ey));
  108. C7ey=zeros(size(ey));
  109. C8ey=zeros(size(ey));

  111. C1ez=zeros(size(ez));
  112. C2ez=zeros(size(ez));
  113. C21ez=zeros(size(ez));
  114. C3ez=zeros(size(ez));
  115. C4ez=zeros(size(ez));
  116. C41ez=zeros(size(ez));
  117. C5ez=zeros(size(ez));
  118. C6ez=zeros(size(ez));
  119. C7ez=zeros(size(ez));
  120. C8ez=zeros(size(ez));

  122. D1hx=zeros(size(hx));
  123. D2hx=zeros(size(hx));
  124. D3hx=zeros(size(hx));
  125. D4hx=zeros(size(hx));
  126. D5hx=zeros(size(hx));
  127. D6hx=zeros(size(hx));

  129. D1hy=zeros(size(hy));
  130. D2hy=zeros(size(hy));
  131. D3hy=zeros(size(hy));
  132. D4hy=zeros(size(hy));
  133. D5hy=zeros(size(hy));
  134. D6hy=zeros(size(hy));

  136. D1hz=zeros(size(hz));
  137. D2hz=zeros(size(hz));
  138. D3hz=zeros(size(hz));
  139. D4hz=zeros(size(hz));   
  140. D5hz=zeros(size(hz));
  141. D6hz=zeros(size(hz));
  142. %--------------------------------------------------------
  143. %    Update coefficients,as described in section 7.8.2  (have some questions)
  144. %--------------------------------------------------------
  145. C1=1.0;
  146. C2=dt;
  147. C21=1.0;
  148. C3=1.0;
  149. C4=1.0/2.0/epsr(1)/epsr(1)/epsz/epsz;
  150. C41=1.0/2.0/epsr(1)/epsz;
  151. C5=2.0*epsr(1)*epsz;
  152. C6=2.0*epsr(1)*epsz;
  153. C7=1.0;
  154. C8=dt/epsr(1)/epsz;

  156. D1=1.0;
  157. D2=dt;
  158. D3=1.0;
  159. D4=1.0/2.0/epsr(1)/epsz/mur(1)/muz;
  160. D5=2.0*epsr(1)*epsz;
  161. D6=2.0*epsr(1)*epsz;

  163. %-------------------------------------------------------
  164. %   Initialize main grid update coefficients     (considering  the basic cell)
  165. %-------------------------------------------------------
  166. C1ex(:, jh_bc:jh_tot-upml, :)=C1;
  167. C2ex(:, jh_bc:jh_tot-upml, :)=C2;
  168. C21ex(:, jh_bc:jh_tot-upml, :)=C21;
  169. C3ex(:, :, kh_bc:kh_tot-upml)=C3;
  170. C4ex(:, :, kh_bc:kh_tot-upml)=C4;
  171. C41ex(:, :, kh_bc:kh_tot-upml)=C41;
  172. C5ex(ih_bc:ie_tot-upml, :, :)=C5;
  173. C6ex(ih_bc:ie_tot-upml, :, :)=C6;
  174. C7ex(1:ie_tot,1:jh_tot,1:kh_tot)=C7;
  175. C8ex(1:ie_tot,1:jh_tot,1:kh_tot)=C8;

  177. C1ey(:, :, kh_bc:kh_tot-upml)=C1;
  178. C2ey(:, :, kh_bc:kh_tot-upml)=C2;
  179. C21ey(:, :, kh_bc:kh_tot-upml)=C21;
  180. C3ey(ih_bc:ih_tot-upml, :, :)=C3;
  181. C4ey(ih_bc:ih_tot-upml, :, :)=C4;
  182. C41ey(ih_bc:ih_tot-upml, :, :)=C41;
  183. C5ey(:, jh_bc:je_tot-upml, :)=C5;
  184. C6ey(:, jh_bc:je_tot-upml, :)=C6;
  185. C7ey(1:ih_tot,1:je_tot,1:kh_tot)=C7;
  186. C8ey(1:ih_tot,1:je_tot,1:kh_tot)=C8;

  188. C1ez(ih_bc:ih_tot-upml, :, :)=C1;
  189. C2ez(ih_bc:ih_tot-upml, :, :)=C2;
  190. C21ez(ih_bc:ih_tot-upml, :, :)=C21;
  191. C3ez(:, jh_bc:jh_tot-upml, :)=C3;
  192. C4ez(:, jh_bc:jh_tot-upml, :)=C4;
  193. C41ez(:, jh_bc:jh_tot-upml, :)=C41;
  194. C5ez(:, :, kh_bc:ke_tot-upml)=C5;
  195. C6ez(:, :, kh_bc:ke_tot-upml)=C6;
  196. C7ez(1:ih_tot,1:jh_tot,1:ke_tot)=C7;
  197. C8ez(1:ih_tot,1:jh_tot,1:ke_tot)=C8;

  199. D1hx(:, jh_bc:je_tot-upml, :)=D1;
  200. D2hx(:, jh_bc:je_tot-upml, :)=D2;
  201. D3hx(:, :, kh_bc:ke_tot-upml)=D3;
  202. D4hx(:, :, kh_bc:ke_tot-upml)=D4;
  203. D5hx(ih_bc:ih_tot-upml, :, :)=D5;
  204. D6hx(ih_bc:ih_tot-upml, :, :)=D6;

  206. D1hy(:, :, kh_bc:ke_tot-upml)=D1;
  207. D2hy(:, :, kh_bc:ke_tot-upml)=D2;
  208. D3hy(ih_bc:ie_tot-upml, :, :)=D3;
  209. D4hy(ih_bc:ie_tot-upml, :, :)=D4;
  210. D5hy(:, jh_bc:jh_tot-upml, :)=D5;
  211. D6hy(:, jh_bc:jh_tot-upml, :)=D6;

  213. D1hz(ih_bc:ie_tot-upml, :, :)=D1;
  214. D2hz(ih_bc:ie_tot-upml, :, :)=D2;
  215. D3hz(:, jh_bc:je_tot-upml, :)=D3;
  216. D4hz(:, jh_bc:je_tot-upml, :)=D4;
  217. D5hz(:, :, kh_bc:kh_tot-upml)=D5;
  218. D6hz(:, :, kh_bc:kh_tot-upml)=D6;

  220. %--------------------------------------------
  221. %   Fill in PML regions
  222. %---------------------------------------------
  223. rmax=exp(-16);
  224. orderbc=4;
  225. %   x-varying material properties
  226. delbc=upml*delta;
  227. % sigmam=-log(rmax)*(orderbc+1.0)/(2.0*eta*delbc);
  228. sigmam=-log(rmax)*epsr*epsz*cc*(orderbc+1.0)/(2.0*delbc);
  229. sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1.0));
  230. kmax=1;
  231. kfactor=(kmax-1.0)/delta/(orderbc+1.0)/delbc^orderbc;

  233. for i=1:upml
  234.     % Coefficients for field components in the center of the grid cell
  235.     x1=(upml-i+1)*delta;
  236.     x2=(upml-i)*delta;
  237.     sigma=sigfactor*(x1^(orderbc+1)-x2^(orderbc+1));
  238.     ki=1+kfactor*(x1^(orderbc+1)-x2^(orderbc+1));
  239.     facm=(2*epsr(1)*epsz*ki-sigma*dt);
  240.     facp=(2*epsr(1)*epsz*ki+sigma*dt);
  241.    
  242.     C5ex(i, :, :)=facp;
  243.     C5ex(ie_tot-i+1, :, :)=facp;
  244.     C6ex(i, :, :)=facm;
  245.     C6ex(ie_tot-i+1, :, :)=facm;
  246.     D1hz(i, :, :)=facm/facp;
  247.     D1hz(ie_tot-i+1, :, :)=facm/facp;
  248.     D2hz(i, :, :)=2.0*epsr(1)*epsz*dt/facp;
  249.     D2hz(ie_tot-i+1, :, :)=2.0*epsr(1)*epsz*dt/facp;
  250.     D3hy(i, :, :)=facm/facp;
  251.     D3hy(ie_tot-i+1, :, :)=facm/facp;
  252.     D4hy(i, :, :)=1.0/facp/mur(1)/muz;
  253.     D4hy(ie_tot-i+1, :, :)=1.0/facp/mur(1)/muz;
  254.    
  255.     % Coefficients for field components on the grid cell boundary
  256.     x1=(upml-i+1.5)*delta;
  257.     x2=(upml-i+0.5)*delta;
  258.     sigma=sigfactor*(x1^(orderbc+1)-x2^(orderbc+1));
  259.     ki=1+kfactor*(x1^(orderbc+1)-x2^(orderbc+1));
  260.     facm=(2.0*epsr(1)*epsz*ki-sigma*dt);
  261.     facp=(2.0*epsr(1)*epsz*ki+sigma*dt);
  262.    
  263.     C1ez(i,:,:)=facm/facp;
  264.     C1ez(ih_tot-i+1,:,:)=facm/facp;
  265.     C2ez(i,:,:)=2.0*epsr(1)*epsz*dt/facp;
  266.     C2ez(ih_tot-i+1,:,:)=2.0*epsr(1)*epsz*dt/facp;
  267.     C21ez(i,:,:)=2.0*epsr(1)*epsz/facp;
  268.     C21ez(ih_tot-i+1,:,:)=2.0*epsr(1)*epsz/facp;
  269.     C3ey(i,:,:)=facm/facp;
  270.     C3ey(ih_tot-i+1,:,:)=facm/facp;
  271.     C4ey(i,:,:)=1.0/facp/epsr(1)/epsz;
  272.     C4ey(ih_tot-i+1,:,:)=1.0/facp/epsr(1)/epsz;
  273.     C41ey(i,:,:)=1.0/facp;
  274.     C41ey(ih_tot-i+1,:,:)=1.0/facp;
  275.     D5hx(i,:,:)=facp;
  276.     D5hx(ih_tot-i+1,:,:)=facp;
  277.     D6hx(i,:,:)=facm;
  278.     D6hx(ih_tot-i+1,:,:)=facm;
  279. end

  281. %  PEC walls
  282. C1ez(1,:,:)=-1.0;
  283. C1ez(ih_tot,:,:)=-1.0;
  284. C2ez(1,:,:)=0.0;
  285. C2ez(ih_tot,:,:)=0.0;
  286. C21ez(1,:,:)=0.0;
  287. C21ez(ih_tot,:,:)=0.0;
  288. C3ey(1,:,:)=-1.0;
  289. C3ey(ih_tot,:,:)=-1.0;
  290. C4ey(1,:,:)=0.0;
  291. C4ey(ih_tot,:,:)=0.0;
  292. C41ey(1,:,:)=0.0;
  293. C41ey(ih_tot,:,:)=0.0;
  294. %  y-varying material properties
  295. delbc=upml*delta;
  296. sigmam=-log(rmax)*epsr*epsz*cc*(orderbc+1.0)/(2.0*delbc);    %(some questions)
  297. sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1.0));
  298. kmax=1.0;
  299. kfactor=(kmax-1.0)/delta/(orderbc+1.0)/delbc^orderbc;
  300. for j=1:upml
  301.     % Coefficients for field componments in the center of the grid cell
  302.     y1=(upml-j+1)*delta;
  303.     y2=(upml-j)*delta;
  304.     sigma=sigfactor*(y1^(orderbc+1)-y2^(orderbc+1));
  305.     ki=1+kfactor*(y1^(orderbc+1)-y2^(orderbc+1));
  306.     facm=2*epsr(1)*epsz*ki-sigma*dt;
  307.     facp=2*epsr(1)*epsz*ki+sigma*dt;
  308.    
  309.     C5ey(:,j,:)=facp;
  310.     C5ey(:,je_tot-j+1,:)=facp;
  311.     C6ey(:,j,:)=facm;
  312.     C6ey(:,je_tot-j+1,:)=facm;
  313.     D1hx(:,j,:)=facm/facp;
  314.     D1hx(:,je_tot-j+1,:)=facm/facp;
  315.     D2hx(:,j,:)=2*epsr(1)*epsz*dt/facp;
  316.     D2hx(:,je_tot-j+1,:)=2*epsr(1)*epsz*dt/facp;
  317.     D3hz(:,j,:)=facm/facp;
  318.     D3hz(:,je_tot-j+1,:)=facm/facp;
  319.     D4hz(:,j,:)=1/facp/mur(1)/muz;
  320.     D4hz(:,je_tot-j+1,:)=1/facp/mur(1)/muz;
  321.    
  322.     % Coefficients for field components on the grid cell boundary
  323.     y1=(upml-j+1.5)*delta;
  324.     y2=(upml-j+0.5)*delta;
  325.     sigma=sigfactor*(y1^(orderbc+1)-y2^(orderbc+1));
  326.     ki=1+kfactor*(y1^(orderbc+1)-y2^(orderbc+1));
  327.     facm=(2*epsr(1)*epsz*ki-sigma*dt);
  328.     facp=(2*epsr(1)*epsz*ki+sigma*dt);
  329.    
  330.     C1ex(:,j,:)=facm/facp;
  331.     C1ex(:,jh_tot-j+1,:)=facm/facp;
  332.     C2ex(:,j,:)=2*epsr(1)*epsz*dt/facp;
  333.     C2ex(:,jh_tot-j+1,:)=2*epsr(1)*epsz*dt/facp;
  334.     C21ex(:,j,:)=2*epsr(1)*epsz/facp;
  335.     C21ex(:,jh_tot-j+1,:)=2*epsr(1)*epsz/facp;
  336.     C3ez(:,j,:)=facm/facp;
  337.     C3ez(:,jh_tot-j+1,:)=facm/facp;
  338.     C4ez(:,j,:)=1/facp/epsr(1)/epsz;
  339.     C4ez(:,jh_tot-j+1,:)=1/facp/epsr(1)/epsz;
  340.     C41ez(:,j,:)=1/facp;
  341.     C41ez(:,jh_tot-j+1,:)=1/facp;
  342.     D5hy(:,j,:)=facp;
  343.     D5hy(:,jh_tot-j+1,:)=facp;
  344.     D6hy(:,j,:)=facm;
  345.     D6hy(:,jh_tot-j+1,:)=facm;
  346. end
  347. % PEC walls
  348. C1ex(:,1,:)=-1;
  349. C1ex(:,jh_tot,:)=-1;
  350. C2ex(:,1,:)=0;
  351. C2ex(:,jh_tot,:)=0;
  352. C21ex(:,1,:)=0;
  353. C21ex(:,jh_tot,:)=0;
  354. C3ez(:,1,:)=-1;
  355. C3ez(:,jh_tot,:)=-1;
  356. C4ez(:,1,:)=0;
  357. C4ez(:,jh_tot,:)=0;
  358. C41ez(:,1,:)=0;
  359. C41ez(:,jh_tot,:)=0;
  360. % z-varying material properties
  361. delbc=upml*delta;
  362. sigmam=-log(rmax)*epsr(1)*epsz*cc*(orderbc+1.0)/(2.0*delbc);    %(some questions)
  363. sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1.0));
  364. kmax=1.0;
  365. kfactor=(kmax-1.0)/delta/(orderbc+1.0)/delbc^orderbc;

  367. for k=1:upml
  368.      
  369.      % Coefficients for field components in the center of the grid cell
  370.      z1=(upml-k+1)*delta;
  371.      z2=(upml-k)*delta;
  372.      sigma=sigfactor*(z1^(orderbc+1)-z2^(orderbc+1));
  373.      ki=1+kfactor*(z1^(orderbc+1)-z2^(orderbc+1));
  374.      facm=2*epsr(1)*epsz*ki-sigma*dt;
  375.      facp=2*epsr(1)*epsz*ki+sigma*dt;
  376.      C5ez(:,:,k)=facp;
  377.      C5ez(:,:,ke_tot-k+1)=facp;
  378.      C6ez(:,:,k)=facm;
  379.      C6ez(:,:,ke_tot-k+1)=facm;
  380.      D1hy(:,:,k)=facm/facp;
  381.      D1hy(:,:,ke_tot-k+1)=facm/facp;
  382.      D2hy(:,:,k)=2*epsr(1)*epsz*dt/facp;
  383.      D2hy(:,:,ke_tot-k+1)=2*epsr(1)*epsz*dt/facp;
  384.      D3hx(:,:,k)=facm/facp;
  385.      D3hx(:,:,ke_tot-k+1)=facm/facp;
  386.      D4hx(:,:,k)=1/facp/mur(1)/muz;
  387.      D4hx(:,:,ke_tot-k+1)=1/facp/mur(1)/muz;
  388.      
  389.      % Coefficients for field components on the grid cell boundary
  390.      z1=(upml-k+1.5)*delta;
  391.      z2=(upml-k+0.5)*delta;
  392.      sigma=sigfactor*(z1^(orderbc+1)-z2^(orderbc+1));
  393.      ki=1+kfactor*(z1^(orderbc+1)-z2^(orderbc+1));
  394.      facm=2*epsr(1)*epsz*ki-sigma*dt;
  395.      facp=2*epsr(1)*epsz*ki+sigma*dt;
  396.      C1ey(:,:,k)=facm/facp;
  397.      C1ey(:,:,kh_tot-k+1)=facm/facp;
  398.      C2ey(:,:,k)=2*epsr(1)*epsz*dt/facp;
  399.      C2ey(:,:,kh_tot-k+1)=2*epsr(1)*epsz*dt/facp;
  400.      C21ey(:,:,k)=2*epsr(1)*epsz/facp;
  401.      C21ey(:,:,kh_tot-k+1)=2*epsr(1)*epsz/facp;
  402.      C3ex(:,:,k)=facm/facp;
  403.      C3ex(:,:,kh_tot-k+1)=facm/facp;
  404.      C4ex(:,:,k)=1/facp/epsr(1)/epsz;
  405.      C4ex(:,:,kh_tot-k+1)=1/facp/epsr(1)/epsz;
  406.      C41ex(:,:,k)=1/facp;
  407.      C41ex(:,:,kh_tot-k+1)=1/facp;
  408.      D5hz(:,:,k)=facp;
  409.      D5hz(:,:,kh_tot-k+1)=facp;
  410.      D6hz(:,:,k)=facm;
  411.      D6hz(:,:,kh_tot-k+1)=facm;
  412. end
  413. %   PEC walls
  414. C1ey(:,:,1)=-1;
  415. C1ey(:,:,kh_tot)=-1;
  416. C2ey(:,:,1)=0;
  417. C2ey(:,:,kh_tot)=0;
  418. C21ey(:,:,1)=0;
  419. C21ey(:,:,kh_tot)=0;
  420. C3ex(:,:,1)=-1;
  421. C3ex(:,:,kh_tot)=-1;
  422. C4ex(:,:,1)=0;
  423. C4ex(:,:,kh_tot)=0;
  424. C41ex(:,:,1)=0;
  425. C41ex(:,:,kh_tot)=0;
  426. % figure
  427. % set(gcf,'Double Buffer','on')
  428. %--------------------------------------------------------
  429. %    Begin  time stepping loop
  430. %--------------------------------------------------------
  431. for n=1:nmax
  432.      % Update magnetic field
  433.      
  434.     bstore=bx;
  435.     bx(2:ie_tot,:,:)=D1hx(2:ie_tot,:,:).*bx(2:ie_tot,:,:)-...
  436.                      D2hx(2:ie_tot,:,:).*((ez(2:ie_tot,2:jh_tot,:)-ez(2:ie_tot,1:je_tot,:))-...
  437.                      (ey(2:ie_tot,:,2:kh_tot)-ey(2:ie_tot,:,1:ke_tot)))./delta;
  438.     hx(2:ie_tot,:,:)= D3hx(2:ie_tot,:,:).*hx(2:ie_tot,:,:)+...
  439.                       D4hx(2:ie_tot,:,:).*(D5hx(2:ie_tot,:,:).*bx(2:ie_tot,:,:)-...
  440.                       D6hx(2:ie_tot,:,:).*bstore(2:ie_tot,:,:));            
  441.     bstore=by;
  442.     by(:,2:je_tot,:)=D1hy(:,2:je_tot,:).*by(:,2:je_tot,:)-...
  443.                      D2hy(:,2:je_tot,:).*((ex(:,2:je_tot,2:kh_tot)-ex(:,2:je_tot,1:ke_tot))-...
  444.                      (ez(2:ih_tot,2:je_tot,:)-ez(1:ie_tot,2:je_tot,:)))./delta;
  445.     hy(:,2:je_tot,:)= D3hy(:,2:je_tot,:).*hy(:,2:je_tot,:)+...
  446.                       D4hy(:,2:je_tot,:).*(D5hy(:,2:je_tot,:).*by(:,2:je_tot,:)-...
  447.                       D6hy(:,2:je_tot,:).*bstore(:,2:je_tot,:));
  448.     bstore=bz;
  449.     bz(:,:,2:ke_tot)=D1hz(:,:,2:ke_tot).*bz(:,:,2:ke_tot)-...
  450.                      D2hz(:,:,2:ke_tot).*((ey(2:ih_tot,:,2:ke_tot)-ey(1:ie_tot,:,2:ke_tot))-...
  451.                      (ex(:,2:jh_tot,2:ke_tot)-ex(:,1:je_tot,2:ke_tot)))./delta;
  452.     hz(:,:,2:ke_tot)= D3hz(:,:,2:ke_tot).*hz(:,:,2:ke_tot)+...
  453.                       D4hz(:,:,2:ke_tot).*(D5hz(:,:,2:ke_tot).*bz(:,:,2:ke_tot)-...
  454.                       D6hz(:,:,2:ke_tot).*bstore(:,:,2:ke_tot));
  455.                   
  456.     % Update electric field
  457. %     dstore=dx;
  458. %     dx(:,2:je_tot,2:ke_tot)=C1ex(:,2:je_tot,2:ke_tot).*dx(:,2:je_tot,2:ke_tot)+...
  459. %                             C2ex(:,2:je_tot,2:ke_tot).*((hz(:,2:je_tot,2:ke_tot)-hz(:,1:je_tot-1,2:ke_tot))-...
  460. %                             (hy(:,2:je_tot,2:ke_tot)-hy(:,2:je_tot,1:ke_tot-1)))./delta;
  461. %     dx(is:is+1,js,ks)=dx(is:is+1,js,ks)+J0*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));
  462. %         
  463. %     ex(:,2:je_tot,2:ke_tot)=C3ex(:,2:je_tot,2:ke_tot).*ex(:,2:je_tot,2:ke_tot)+...
  464. %                             C4ex(:,2:je_tot,2:ke_tot).*(C5ex(:,2:je_tot,2:ke_tot).*dx(:,2:je_tot,2:ke_tot)-...
  465. %                             C6ex(:,2:je_tot,2:ke_tot).*dstore(:,2:je_tot,2:ke_tot));

  467.     qstore=qx;
  468.     qx(:,2:je_tot,2:ke_tot)=C7ex(:,2:je_tot,2:ke_tot).*qx(:,2:je_tot,2:ke_tot)+...
  469.                             C8ex(:,2:je_tot,2:ke_tot).*((hz(:,2:je_tot,2:ke_tot)-hz(:,1:je_tot-1,2:ke_tot))-...
  470.                             (hy(:,2:je_tot,2:ke_tot)-hy(:,2:je_tot,1:ke_tot-1)))./delta;
  471.                         
  472.     px(:,2:je_tot,2:ke_tot)=C1ex(:,2:je_tot,2:ke_tot).*px(:,2:je_tot,2:ke_tot)+...
  473.                             C21ex(:,2:je_tot,2:ke_tot).*(qx(:,2:je_tot,2:ke_tot)-qstore(:,2:je_tot,2:ke_tot));
  474.     pstore=px;
  475.     ex(is:is+1,js,ks)=ex(is:is+1,js,ks)+J0.*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));  
  476.    
  477.     ex(:,2:je_tot,2:ke_tot)=C3ex(:,2:je_tot,2:ke_tot).*ex(:,2:je_tot,2:ke_tot)+...
  478.                             C41ex(:,2:je_tot,2:ke_tot).*(C5ex(:,2:je_tot,2:ke_tot).*px(:,2:je_tot,2:ke_tot)-...
  479.                             C6ex(:,2:je_tot,2:ke_tot).*pstore(:,2:je_tot,2:ke_tot));      

  481. %     dstore=dy;
  482. %     dy(2:ie_tot,:,2:ke_tot)=C1ey(2:ie_tot,:,2:ke_tot).*dy(2:ie_tot,:,2:ke_tot)+...
  483. %                             C2ey(2:ie_tot,:,2:ke_tot).*((hx(2:ie_tot,:,2:ke_tot)-hx(2:ie_tot,:,1:ke_tot-1))-...
  484. %                             (hz(2:ie_tot,:,2:ke_tot)-hz(1:ie_tot-1,:,2:ke_tot)))./delta;
  485. % %     dy(is,js:js+1,ks)=dz(is,js:js+1,ks)+J0*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));
  486. %     ey(2:ie_tot,:,2:ke_tot)=C3ey(2:ie_tot,:,2:ke_tot).*ey(2:ie_tot,:,2:ke_tot)+...
  487. %                             C4ey(2:ie_tot,:,2:ke_tot).*(C5ey(2:ie_tot,:,2:ke_tot).*dy(2:ie_tot,:,2:ke_tot)-...
  488. %                             C6ey(2:ie_tot,:,2:ke_tot).*dstore(2:ie_tot,:,2:ke_tot));
  489.                         
  490. %     qstore=qy;
  491. %     qy(2:ie_tot,:,2:ke_tot)=C7ey(2:ie_tot,:,2:ke_tot).*qy(2:ie_tot,:,2:ke_tot)+...
  492. %                             C8ey(2:ie_tot,:,2:ke_tot).*((hx(2:ie_tot,:,2:ke_tot)-hx(2:ie_tot,:,1:ke_tot-1))-...
  493. %                             (hz(2:ie_tot,:,2:ke_tot)-hz(1:ie_tot-1,:,2:ke_tot)))./delta;
  494. %      
  495. %     py(2:ie_tot,:,2:ke_tot)=C1ey(2:ie_tot,:,2:ke_tot).*py(2:ie_tot,:,2:ke_tot)+...
  496. %                             C21ey(2:ie_tot,:,2:ke_tot).*(qy(2:ie_tot,:,2:ke_tot)-qstore(2:ie_tot,:,2:ke_tot));
  497. %                        
  498. %     pstore=py;
  499. %     ey(2:ie_tot,:,2:ke_tot)=C3ey(2:ie_tot,:,2:ke_tot).*ey(2:ie_tot,:,2:ke_tot)+...
  500. %                             C41ey(2:ie_tot,:,2:ke_tot).*(C5ey(2:ie_tot,:,2:ke_tot).*py(2:ie_tot,:,2:ke_tot)-...
  501. %                             C6ey(2:ie_tot,:,2:ke_tot).*pstore(2:ie_tot,:,2:ke_tot));
  502. %                 
  503. %     dstore=dz;
  504. %     dz(2:ie_tot,2:je_tot,:)=C1ez(2:ie_tot,2:je_tot,:).*dz(2:ie_tot,2:je_tot,:)+...
  505. %                             C2ez(2:ie_tot,2:je_tot,:).*((hy(2:ie_tot,2:je_tot,:)-hy(1:ie_tot-1,2:je_tot,:))-...
  506. %                             (hx(2:ie_tot,2:je_tot,:)-hx(2:ie_tot,1:je_tot-1,:)))./delta;
  507. % %     dz(is,js,ks:ks+1)=dz(is,js,ks:ks+1)+J0*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));
  508. %     ez(2:ie_tot,2:je_tot,:)=C3ez(2:ie_tot,2:je_tot,:).*ez(2:ie_tot,2:je_tot,:)+...
  509. %                             C4ez(2:ie_tot,2:je_tot,:).*(C5ez(2:ie_tot,2:je_tot,:).*dz(2:ie_tot,2:je_tot,:)-...
  510. %                             C6ez(2:ie_tot,2:je_tot,:).*dstore(2:ie_tot,2:je_tot,:));
  511.   
  512. %     qstore=qz;
  513. %     qz(2:ie_tot,2:je_tot,:)=C7ez(2:ie_tot,2:je_tot,:).*qz(2:ie_tot,2:je_tot,:)+...
  514. %                             C8ez(2:ie_tot,2:je_tot,:).*((hy(2:ie_tot,2:je_tot,:)-hy(1:ie_tot-1,2:je_tot,:))-...
  515. %                             (hx(2:ie_tot,2:je_tot,:)-hx(2:ie_tot,1:je_tot-1,:)))./delta;                  
  516. %     pz(2:ie_tot,2:je_tot,:)=C1ez(2:ie_tot,2:je_tot,:).*pz(2:ie_tot,2:je_tot,:)+...
  517. %                             C21ez(2:ie_tot,2:je_tot,:).*(qz(2:ie_tot,2:je_tot,:)-qstore(2:ie_tot,2:je_tot,:));
  518. %                        
  519. %     pz(is,js,ks:ks+1)=pz(is,js,ks:ks+1)+J0*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));  
  520. % %     c(n)=J0*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));
  521. %     pstore=pz;
  522. % %      pz(is,js,ks:ks+1)=source(n);
  523. %      ez(2:ie_tot,2:je_tot,:)=C3ez(2:ie_tot,2:je_tot,:).*ez(2:ie_tot,2:je_tot,:)+...
  524. %                              C41ez(2:ie_tot,2:je_tot,:).*(C5ez(2:ie_tot,2:je_tot,:).*pz(2:ie_tot,2:je_tot,:)-...
  525. %                              C6ez(2:ie_tot,2:je_tot,:).*pstore(2:ie_tot,2:je_tot,:));


  528.     %***********************************************************************
  529.     %     Visualize fields
  530.     %***********************************************************************

  532.     timestep=int2str(n);
  533. %     tview(:,:)=squeeze(ex(ih_bc:upml+ie,jh_bc:upml+je,ks));
  534. %     sview(:,:)=squeeze(ex(ih_bc:upml+ie,js,kh_bc:upml+ke));
  535.     tview(:,:)=squeeze(ex(is,jh_bc:upml+je,kh_bc:upml+ke));
  536.     sview(:,:)=squeeze(ex(ih_bc:upml+ie,js,kh_bc:upml+ke));
  537. %     tview(:,:)=squeeze(hz(ih_bc:upml+ie,jh_bc:upml+je,ks));
  538. %     sview(:,:)=squeeze(hz(ih_bc:upml+ie,js,kh_bc:upml+ke));
  539.     subplot('position',[0.15 0.57 0.7 0.35])
  540.     %imagesc(tview');

  542.      pcolor(tview');
  543.     shading interp;
  544.     caxis([-0.2 0.2]);
  545.     colorbar;
  546.     axis image; axis xy;
  547.     title(['E_z(i,j,k=k_s_o_u_r_c_e), time step = ',timestep]);
  548.     xlabel('i coordinate');
  549.     ylabel('j coordinate');
  550.    
  551.     subplot('position',[0.15 0.08 0.7 0.35])
  552.     %imagesc(tview');
  553.     pcolor(tview');
  554.     shading interp;
  555.     caxis([-0.2 0.2]);
  556.     colorbar;
  557.     axis image; axis xy;
  558.     title(['E_z(i,j=j_s_o_u_r_c_e,k), time step = ',timestep]);
  559.     xlabel('i coordinate');
  560.     ylabel('k coordinate');
  561.     pause(0.05)
  562.    
  563. end
复制代码
tingo 发表于 2012-12-4 22:09

                               
登录/注册后可看大图

不客气，有问题欢迎在此版区讨论！

我觉得是程序中电场的迭代下面的迭代出现了问题：
ex(:,2:je_tot,2:ke_tot)=C3ex(:,2:je_tot,2:ke_tot).*ex(:,2:je_tot,2:ke_tot)+...
                            C41ex(:,2:je_tot,2:ke_tot).*(C5ex(:,2:je_tot,2:ke_tot).*px(:,2:je_tot,2:ke_tot)-...
                            C6ex(:,2:je_tot,2:ke_tot).*pstore(:,2:je_tot,2:ke_tot));
但这个和葛老师书上的公式一致，我就不知道该如何处理这个问题
你参考这本书中的代码The Finite-Difference Time-Domain Method for Electromagntics with MATLAB Simulations
这本书现在有中文版了，名字叫MATLAB模拟的电磁学时域有限差分法。
英文电子版微网上好像有，你找找看。
这本书中对CPML讲的详细，你再看看，尽可能自己将程序调通，这是一个学习的过程。有原理不懂的可以贴出来，和大家探讨。
客服中心 搜索
关于我们
关于我们
关注我们
联系我们
帮助中心
资讯中心
企业生态
社区论坛
服务支持
资源下载
售后服务
推广服务
关注我们
官方微博
官方空间
官方微信
诚信网站 可信网站 Powered by Discuz! X3.5 © 2008-2024 头雁微网 版权所有 蜀ICP备2023016130号-1 客服QQ: 68870544  技术支持: 克米设计
返回顶部