主网格与UPML区域分离主迭代程序

////////////////////////////////////////////////////////////////////////// // 时域循环主迭代 ////////////////////////////////////////////////////////////////////////// void begin_time_stepping_loop() { double*** bstore = zeros(ih_tot,jh_tot,kh_tot); double*** dstore = zeros(ih_tot,jh_tot,kh_tot); int ii=0,jj=0,kk=0; FILE* EMFile = fopen("D://电磁场空间分布值.TXT","w"); if (!EMFile) { cerr<<"创建文件失败!"; return ; } FILE* RefE = fopen("D://参考点电场数值.TXT","w"); double sig = 0 ; double ca=(1.0-(dt*sig)/(2.0*epsz*epsr))/(1.0+(dt*sig)/(2.0*epsz*epsr)); double cb=(dt/epsz/epsr/delta)/(1.0+(dt*sig)/(2.0*epsz*epsr)); double cp=1.0; double cq=dt/muz/delta; ////////////////////////////////////////////////////////////////////////// // 时域主循环 ////////////////////////////////////////////////////////////////////////// for(int n=1;n<=nmax;n++){ ez[is][js][ks] +=J0*(n-ndelay)*exp(-pow(n-ndelay,2)/pow(tau,2));//激励源当前时刻的数值; ////////////////////////////////////////////////////////////////////////// // update magnetic field in main grid ////////////////////////////////////////////////////////////////////////// for(ii=upml;ii<=ie_tot-upml;ii++) for(jj=upml;jj<je_tot-upml;jj++) for(kk=upml;kk<ke_tot-upml;kk++) hx[ii][jj][kk] = cp*hx[ii][jj][kk]+cq*((ey[ii][jj][kk+1]-ey[ii][jj][kk])-(ez[ii][jj+1][kk]-ez[ii][jj][kk])); for(ii=upml;ii<ie_tot-upml;ii++) for(jj=upml;jj<=je_tot-upml;jj++) for(kk=upml;kk<ke_tot-upml;kk++) hy[ii][jj][kk] = cp*hy[ii][jj][kk]+cq*((ex[ii][jj][kk]-ex[ii][jj][kk+1])+(ez[ii+1][jj][kk]-ez[ii][jj][kk])); for(ii=upml;ii<ie_tot-upml;ii++) for(jj=upml;jj<je_tot-upml;jj++) for(kk=upml;kk<=ke_tot-upml;kk++) hz[ii][jj][kk] = cp*hz[ii][jj][kk]+cq*((ex[ii][jj+1][kk]-ex[ii][jj][kk])+(ey[ii][jj][kk]-ey[ii+1][jj][kk])); ////////////////////////////////////////////////////////////////////////// // update electric field in main grid ////////////////////////////////////////////////////////////////////////// for(ii=upml;ii<ie_tot-upml;ii++) for(jj=upml;jj<=je_tot-upml;jj++) for(kk=upml;kk<=ke_tot-upml;kk++) ex[ii][jj][kk] = ca*ex[ii][jj][kk]+cb*((hz[ii][jj][kk]-hz[ii][jj-1][kk])+(hy[ii][jj][kk-1]-hy[ii][jj][kk])); for(ii=upml;ii<=ie_tot-upml;ii++) for(jj=upml;jj<je_tot-upml;jj++) for(kk=upml;kk<=ke_tot-upml;kk++) ey[ii][jj][kk] = ca*ey[ii][jj][kk]+cb*((hx[ii][jj][kk]-hx[ii][jj][kk-1])+(hz[ii-1][jj][kk]-hz[ii][jj][kk])); for(ii=upml;ii<=ie_tot-upml;ii++) for(jj=upml;jj<=je_tot-upml;jj++) for(kk=upml;kk<ke_tot-upml;kk++) ez[ii][jj][kk] =ca*ez[ii][jj][kk]+cb*((hx[ii][jj-1][kk]-hx[ii][jj][kk])+(hy[ii][jj][kk]-hy[ii-1][jj][kk])); ////////////////////////////////////////////////////////////////////////// // update magnetic field in upml region ////////////////////////////////////////////////////////////////////////// store_em(bstore,bx,ih_tot,je_tot,ke_tot);//bstore = bx ; size of bx = ih_tot,je_tot,ke_tot for(ii=1;ii<ie_tot;ii++)// bx(2:ie_tot,:,:) ; size of bx = ih_tot,je_tot,ke_tot ;debug 2010-11-10 for(jj=0;jj<je_tot;jj++) for(kk=0;kk<ke_tot;kk++){ if (!(ii>=upml&&ii<=ie_tot-upml&&jj>=upml&&jj<je_tot-upml&&kk>=upml&&kk<ke_tot-upml)) { bx[ii][jj][kk] = D1hx[ii][jj][kk]*bx[ii][jj][kk]- D2hx[ii][jj][kk]*((ez[ii][jj+1][kk]-ez[ii][jj][kk])-(ey[ii][jj][kk+1]-ey[ii][jj][kk]))/delta; hx[ii][jj][kk] = D3hx[ii][jj][kk]*hx[ii][jj][kk]+ D4hx[ii][jj][kk]*(D5hx[ii][jj][kk]*bx[ii][jj][kk]-D6hx[ii][jj][kk]*bstore[ii][jj][kk]); } } ////////////////////////////////////////////////////////////////////////// store_em(bstore,by,ie_tot,jh_tot,ke_tot);//bstore = by ; size of by = ie_tot,jh_tot,ke_tot; for(ii=0;ii<ie_tot;ii++)// by(:,2:je_tot,:) size of by = ie_tot,jh_tot,ke_tot;--------debug 2010-11-10 for(jj=1;jj<je_tot;jj++) for(kk=0;kk<ke_tot;kk++){ if (!(ii>=upml&&ii<ie_tot-upml&&jj>=upml&&jj<=je_tot-upml&&kk>=upml&&kk<ke_tot-upml)) { by[ii][jj][kk] = D1hy[ii][jj][kk]*by[ii][jj][kk]- D2hy[ii][jj][kk]*((ex[ii][jj][kk+1]-ex[ii][jj][kk])-(ez[ii+1][jj][kk]-ez[ii][jj][kk]))/delta; hy[ii][jj][kk] = D3hy[ii][jj][kk]*hy[ii][jj][kk]+ D4hy[ii][jj][kk]*(D5hy[ii][jj][kk]*by[ii][jj][kk]-D6hy[ii][jj][kk]*bstore[ii][jj][kk]); } } ////////////////////////////////////////////////////////////////////////// store_em(bstore,bz,ie_tot,je_tot,kh_tot);//bstore = bz ; size of bz = ie_tot,je_tot,kh_tot for(ii=0;ii<ie_tot;ii++)//bz(:,:,2:ke_tot) size of bz = ie_tot,je_tot,kh_tot for(jj=0;jj<je_tot;jj++) for(kk=1;kk<ke_tot;kk++){ if (!(ii>=upml&&ii<ie_tot-upml&&jj>=upml&&jj<je_tot-upml&&kk>=upml&&kk<=ke_tot-upml)) { bz[ii][jj][kk] = D1hz[ii][jj][kk]*bz[ii][jj][kk]- D2hz[ii][jj][kk]*((ey[ii+1][jj][kk]-ey[ii][jj][kk])-(ex[ii][jj+1][kk]-ex[ii][jj][kk]))/delta; hz[ii][jj][kk] = D3hz[ii][jj][kk]*hz[ii][jj][kk]+ D4hz[ii][jj][kk]*(D5hz[ii][jj][kk]*bz[ii][jj][kk]-D6hz[ii][jj][kk]*bstore[ii][jj][kk]); } } ////////////////////////////////////////////////////////////////////////// // update electric field in upml region ////////////////////////////////////////////////////////////////////////// store_em(dstore,dx,ie_tot,jh_tot,kh_tot);//dstore = dx ;size of dx = ie_tot,jh_tot,kh_tot for(ii=0;ii<ie_tot;ii++)//dx(:,2:je_tot,2:ke_tot) for(jj=1;jj<je_tot;jj++) for(kk=1;kk<ke_tot;kk++){ if (!(ii>=upml&&ii<ie_tot-upml&&jj>=upml&&jj<=je_tot-upml&&kk>=upml&&kk<=ke_tot-upml)) { dx[ii][jj][kk]=C1ex[ii][jj][kk]*dx[ii][jj][kk] + C2ex[ii][jj][kk]*((hz[ii][jj][kk]-hz[ii][jj-1][kk])-(hy[ii][jj][kk]-hy[ii][jj][kk-1]))/delta; ex[ii][jj][kk]=C3ex[ii][jj][kk]*ex[ii][jj][kk]+ C4ex[ii][jj][kk]*(C5ex[ii][jj][kk]*dx[ii][jj][kk]-C6ex[ii][jj][kk]*dstore[ii][jj][kk]); } } ////////////////////////////////////////////////////////////////////////// store_em(dstore,dy,ih_tot,je_tot,kh_tot);//dstore = dy size of dy = ih_tot,je_tot,kh_tot for(ii=1;ii<ie_tot;ii++)//dy(2:ie_tot,:,2:ke_tot) for(jj=0;jj<je_tot;jj++) for(kk=1;kk<ke_tot;kk++){ if (!(ii>=upml&&ii<=ie_tot-upml&&jj>=upml&&jj<je_tot-upml&&kk>=upml&&kk<=ke_tot-upml)) { dy[ii][jj][kk]=C1ey[ii][jj][kk]*dy[ii][jj][kk]+ C2ey[ii][jj][kk]*((hx[ii][jj][kk]-hx[ii][jj][kk-1])-(hz[ii][jj][kk]-hz[ii-1][jj][kk]))/delta; ey[ii][jj][kk]=C3ey[ii][jj][kk]*ey[ii][jj][kk]+ C4ey[ii][jj][kk]*(C5ey[ii][jj][kk]*dy[ii][jj][kk]-C6ey[ii][jj][kk]*dstore[ii][jj][kk]); } } ////////////////////////////////////////////////////////////////////////// store_em(dstore,dz,ih_tot,jh_tot,ke_tot);//dstore = dz ; size of dz = ih_tot,jh_tot,ke_tot for(ii=1;ii<ie_tot;ii++)//dz(2:ie_tot,2:je_tot,:)--------debug 2010-11-10 for(jj=1;jj<je_tot;jj++) for(kk=0;kk<ke_tot;kk++){ if (!(ii>=upml&&ii<=ie_tot-upml&&jj>=upml&&jj<=je_tot-upml&&kk>=upml&&kk<ke_tot-upml)) { dz[ii][jj][kk]=C1ez[ii][jj][kk]*dz[ii][jj][kk]+ C2ez[ii][jj][kk]*((hy[ii][jj][kk]-hy[ii-1][jj][kk])-(hx[ii][jj][kk]-hx[ii][jj-1][kk]))/delta; ez[ii][jj][kk]=C3ez[ii][jj][kk]*ez[ii][jj][kk]+ C4ez[ii][jj][kk]*(C5ez[ii][jj][kk]*dz[ii][jj][kk]-C6ez[ii][jj][kk]*dstore[ii][jj][kk]); } } /**/ ////////////////////////////////////////////////////////////////////////// if (n==nmax-1) { for(int idx=0;idx<ie;idx++){ for(int idy=0;idy<je;idy++){ fprintf(EMFile,"%le ",ez[idx+upml][idy+upml][ks+upml]); } fprintf(EMFile,"/n"); } } fprintf(RefE,"%le /n",ez[is/2][js][ks]); } fclose(EMFile); fclose(RefE); }

 

代码中

主网格为空气

吸收层分开单独做循环

 

该版本实现了功能，但还比较浪费内存。