📄 erweifd2d3[1].3.1.cpp
字号:
cb1=(1-ca1)/(sigmax*dx);
for(j=1;j<=je;j++)
{
caeybcr[i][j]=ca1 ;
cbeybcr[i][j]=cb1 ;
}
for(j=1;j<=jebc;j++)
{
caeybcf[i+iebc+ie][j]=ca1 ;
cbeybcf[i+iebc+ie][j]=cb1 ;
caeybcb[i+iebc+ie][j]=ca1 ;
cbeybcb[i+iebc+ie][j]=cb1 ;
}
}
sigmax=bcfactor*pow((0.5*dx),(orderbc+1));
ca1=exp(-sigmax*dt/(epsz*eps[1]));
cb1=(1-ca1)/(sigmax*dx);
for(j=1;j<=je;j++)
{
caey[ib][j]=ca1 ;
cbey[ib][j]=cb1 ;
}
for(j=1;j<=jebc;j++)
{
caeybcf[iebc+ib][j]=ca1 ;
cbeybcf[iebc+ib][j]=cb1 ;
caeybcb[iebc+ib][j]=ca1 ;
cbeybcb[iebc+ib][j]=cb1 ;
}
PRECISION**dahzxbcr=MyAllocPrecision2D(iebc,je,0.0);
PRECISION**dbhzxbcr=MyAllocPrecision2D(iebc,je,0.0);
PRECISION**dahzybcr=MyAllocPrecision2D(iebc,je,0.0);
PRECISION**dbhzybcr=MyAllocPrecision2D(iebc,je,0.0);
for(i=1;i<=iebc;i++)
{
x1=i*dx ;
x2=(i-1)*dx ;
sigmax=bcfactor*(pow(x1,(orderbc+1))-pow(x2,(orderbc+1)));
sigmaxs=sigmax*(muz/(epsz*eps[1]));
da1=exp(-sigmaxs*dt/muz);
db1=(1-da1)/(sigmaxs*dx);
for(j=1;j<=je;j++)
{
dahzxbcr[i][j]=da1 ;
dbhzxbcr[i][j]=db1 ;
dahzybcr[i][j]=da[1];
dbhzybcr[i][j]=db[1];
}
for(j=2;j<=je;j++)
{
caexbcr[i][j]=ca[1];
cbexbcr[i][j]=cb[1];
}
for(j=1;j<=jebc;j++)
{
dahzxbcf[i+ie+iebc][j]=da1 ;
dbhzxbcf[i+ie+iebc][j]=db1 ;
dahzxbcb[i+ie+iebc][j]=da1 ;
dbhzxbcb[i+ie+iebc][j]=db1 ;
}
}
/***********************************************************************
open a file "twave.txt" to record trans wave
***********************************************************************/
FILE * Fptwave;
Fptwave=fopen("twave.txt","w");
FILE * Fptest;
Fptest=fopen("test.txt","w");
/***********************************************************************
BEGIN TIME-STEPPING LOOP
***********************************************************************/
cout<<"BEGIN TIME-STEPPING LOOP"<<endl ;
for(n=1;n<=nmax;n++)
{
cout<<"Time step "<<n<<endl;
/************************************************************************/
/* Update electric fields EX in main grid */
/************************************************************************/
for(i=1;i<=ie;i++)
{
for(j=2;j<=je;j++)
{
ex[i][j]=caex[i][j]*ex[i][j]+cbex[i][j]*(hz[i][j]-hz[i][j-1]);
}
}
/************************************************************************/
/* Update electric fields EY in main grid */
/************************************************************************/
for(i=2;i<=ie;i++)
{
for(j=1;j<=je;j++)
{
ey[i][j]=caey[i][j]*ey[i][j]+cbey[i][j]*(hz[i-1][j]-hz[i][j]);
}
}
/***********************************************************************
Update EX in PML regions
***********************************************************************/
// FRONT
for(i=1;i<=iefbc;i++)
{
for(j=2;j<=jebc;j++)
{
exbcf[i][j]=caexbcf[i][j]*exbcf[i][j]
-cbexbcf[i][j]*(hzxbcf[i][j-1]+hzybcf[i][j-1]
-hzxbcf[i][j]-hzybcf[i][j]);
}
}
for(i=1;i<=ie;i++)
{
ex[i][1]=caex[i][1]*ex[i][1]
-cbex[i][1]*(hzxbcf[iebc+i][jebc]
+hzybcf[iebc+i][jebc]-hz[i][1]);
}
// BACK
for(i=1;i<=iefbc;i++)
{
for(j=2;j<=(jebc-1);j++)
{
exbcb[i][j]=caexbcb[i][j]*exbcb[i][j]
-cbexbcb[i][j]*(hzxbcb[i][j-1]+hzybcb[i][j-1]
-hzxbcb[i][j]-hzybcb[i][j]);
}
}
for(i=1;i<=ie;i++)
{
ex[i][jb]=caex[i][jb]*ex[i][jb]
-cbex[i][jb]*(hz[i][jb-1]-hzxbcb[iebc+i][1]
-hzybcb[iebc+i][1]);
}
// LEFT
for(i=1;i<=iebc;i++)
{
for(j=2;j<=je;j++)
{
exbcl[i][j]=caexbcl[i][j]*exbcl[i][j]
-cbexbcl[i][j]*(hzxbcl[i][j-1]+hzybcl[i][j-1]
-hzxbcl[i][j]-hzybcl[i][j]);
}
}
for(i=1;i<=iebc;i++)
{
exbcl[i][1]=caexbcl[i][1]*exbcl[i][1]
-cbexbcl[i][1]*(hzxbcf[i][jebc]+hzybcf[i][jebc]
-hzxbcl[i][1]-hzybcl[i][1]);
exbcl[i][jb]=caexbcl[i][jb]*exbcl[i][jb]
-cbexbcl[i][jb]*(hzxbcl[i][je]+hzybcl[i][je]
-hzxbcb[i][1]-hzybcb[i][1]);
}
// RIGHT
for(i=1;i<=iebc;i++)
{
for(j=2;j<=je;j++)
{
exbcr[i][j]=caexbcr[i][j]*exbcr[i][j]
-cbexbcr[i][j]*(hzxbcr[i][j-1]+hzybcr[i][j-1]
-hzxbcr[i][j]-hzybcr[i][j]);
}
}
for(i=1;i<=iebc;i++)
{
exbcr[i][1]=caexbcr[i][1]*exbcr[i][1]
-cbexbcr[i][1]*(hzxbcf[ie+iebc+i][jebc]
+hzybcf[ie+iebc+i][jebc]-hzxbcr[i][1]-hzybcr[i][1]);
exbcr[i][jb]=caexbcr[i][jb]*exbcr[i][jb]
-cbexbcr[i][jb]*(hzxbcr[i][je]+hzybcr[i][je]
-hzxbcb[ie+iebc+i][1]-hzybcb[ie+iebc+i][1]);
}
/***********************************************************************
Update EY in PML regions
***********************************************************************/
for(i=2;i<=iefbc;i++)
{
for(j=1;j<=jebc;j++)
{
// FRONT
eybcf[i][j]=caeybcf[i][j]*eybcf[i][j]-
cbeybcf[i][j]*(hzxbcf[i][j]+hzybcf[i][j]-
hzxbcf[i-1][j]-hzybcf[i-1][j]);
// BACK
eybcb[i][j]=caeybcb[i][j]*eybcb[i][j]-
cbeybcb[i][j]*(hzxbcb[i][j]+hzybcb[i][j]-
hzxbcb[i-1][j]-hzybcb[i-1][j]);
}
}
for(i=2;i<=iebc;i++)
{
for(j=1;j<=je;j++)
{
// LEFT
eybcl[i][j]=caeybcl[i][j]*eybcl[i][j]-
cbeybcl[i][j]*(hzxbcl[i][j]+hzybcl[i][j]-
hzxbcl[i-1][j]-hzybcl[i-1][j]);
// RIGHT
eybcr[i][j]=caeybcr[i][j]*eybcr[i][j]-
cbeybcr[i][j]*(hzxbcr[i][j]+hzybcr[i][j]-
hzxbcr[i-1][j]-hzybcr[i-1][j]);
}
}
for(j=1;j<=je;j++)
{
// LEFT
ey[1][j]=caey[1][j]*ey[1][j]-cbey[1][j]*(hz[1][j]-hzxbcl[iebc][j]-hzybcl[iebc][j]);
// RIGHT
ey[ib][j]=caey[ib][j]*ey[ib][j]-cbey[ib][j]*(hzxbcr[1][j]+hzybcr[1][j]-hz[ie][j]);
}
/************************************************************************/
/* Update magnetic fields (HZ) in main grid */
/************************************************************************/
for(i=1;i<=ie;i++)
{
for(j=1;j<=je;j++)
{
hz[i][j]=dahz[i][j]*hz[i][j]+dbhz[i][j]*(ex[i][j+1]-ex[i][j]+ey[i][j]-ey[i+1][j]);
}
}
hz[is][js]=source[n];
fprintf(Fptwave,"%d %e %e %e %e %e %e %e %e %e %e \n",n,hz[ip][jp-50],hz[ip][jp-40],hz[ip][jp-30],hz[ip][jp-20],hz[ip][jp-10],hz[ip][jp],hz[ip][jp+10],hz[ip][jp+20],hz[ip][jp+30],hz[ip][jp+40]); //out put the hz value at [ip][jp] to file
/***********************************************************************
Update HZX in PML regions
***********************************************************************/
for(i=1;i<=iefbc;i++)
{
for(j=1;j<=jebc;j++)
{
// FRONT
hzxbcf[i][j]=dahzxbcf[i][j]*hzxbcf[i][j]-
dbhzxbcf[i][j]*(eybcf[i+1][j]-eybcf[i][j]);
// BACK
hzxbcb[i][j]=dahzxbcb[i][j]*hzxbcb[i][j]-
dbhzxbcb[i][j]*(eybcb[i+1][j]-eybcb[i][j]);
}
}
// LEFT
for(i=1;i<=iebc-1;i++)
{
for(j=1;j<=je;j++)
{
hzxbcl[i][j]=dahzxbcl[i][j]*hzxbcl[i][j]-
dbhzxbcl[i][j]*(eybcl[i+1][j]-eybcl[i][j]);
}
}
for(j=1;j<=je;j++)
{
hzxbcl[iebc][j]=dahzxbcl[iebc][j]*hzxbcl[iebc][j]-
dbhzxbcl[iebc][j]*(ey[1][j]-eybcl[iebc][j]);
}
// RIGHT
for(i=2;i<=iebc;i++)
{
for(j=1;j<=je;j++)
{
hzxbcr[i][j]=dahzxbcr[i][j]*hzxbcr[i][j]-
dbhzxbcr[i][j]*(eybcr[i+1][j]-eybcr[i][j]);
}
}
for(j=1;j<=je;j++)
{
hzxbcr[1][j]=dahzxbcr[1][j]*hzxbcr[1][j]-
dbhzxbcr[1][j]*(eybcr[2][j]-ey[ib][j]);
}
/***********************************************************************
Update HZY in PML regions
***********************************************************************/
// FRONT
for(i=1;i<=iefbc;i++)
{
for(j=1;j<=jebc-1;j++)
{
hzybcf[i][j]=dahzybcf[i][j]*hzybcf[i][j]-
dbhzybcf[i][j]*(exbcf[i][j]-exbcf[i][j+1]);
}
}
for(i=1;i<=iebc;i++)
{
hzybcf[i][jebc]=dahzybcf[i][jebc]*hzybcf[i][jebc]-
dbhzybcf[i][jebc]*(exbcf[i][jebc]-exbcl[i][1]);
}
for(i=iebc+1;i<=iebc+ie;i++)
{
hzybcf[i][jebc]=
dahzybcf[i][jebc]*hzybcf[i][jebc]-
dbhzybcf[i][jebc]*(exbcf[i][jebc]-ex[i-iebc][1]);
}
for(i=iebc+ie+1;i<=iefbc;i++)
{
hzybcf[i][jebc]=
dahzybcf[i][jebc]*hzybcf[i][jebc]-
dbhzybcf[i][jebc]*(exbcf[i][jebc]-
exbcr[i-ie-iebc][1]);
}
// BACK
for(i=1;i<=iefbc;i++)
{
for(j=2;j<=jebc;j++)
{
hzybcb[i][j]=dahzybcb[i][j]*hzybcb[i][j]-
dbhzybcb[i][j]*(exbcb[i][j]-exbcb[i][j+1]);
}
}
for(i=1;i<=iebc;i++)
{
hzybcb[i][1]=dahzybcb[i][1]*hzybcb[i][1]-
dbhzybcb[i][1]*(exbcl[i][jb]-exbcb[i][2]);
}
for(i=iebc+1;i<=iebc+ie;i++)
{
hzybcb[i][1]=
dahzybcb[i][1]*hzybcb[i][1]-
dbhzybcb[i][1]*(ex[i-iebc][jb]-exbcb[i][2]);
}
for(i=iebc+ie+1;i<=iefbc;i++)
{
hzybcb[i][1]=
dahzybcb[i][1]*hzybcb[i][1]-
dbhzybcb[i][1]*(exbcr[i-ie-iebc][jb]-exbcb[i][2]);
}
for(i=1;i<=iebc;i++)
{
for(j=1;j<=je;j++)
{
// LEFT
hzybcl[i][j]=dahzybcl[i][j]*hzybcl[i][j]-
dbhzybcl[i][j]*(exbcl[i][j]-exbcl[i][j+1]);
// RIGHT
hzybcr[i][j]=dahzybcr[i][j]*hzybcr[i][j]-
dbhzybcr[i][j]*(exbcr[i][j]-exbcr[i][j+1]);
}
}
/***********************************************************************
END TIME-STEPPING LOOP
***********************************************************************/
}
FILE*FpProb ;
FpProb=fopen("Prob.txt","w");
for(j=1;j<=je;j++)
{
for(i=1;i<=ie;i++)
{
fprintf(FpProb,"%e ",hz[i][j]);
}
fprintf(FpProb,"\n");
}
fclose(FpProb);
fclose(Fptwave);
fclose(Fptest);
cout<<"dt="<<dt<<endl;
cout<<"done."<<endl;
return 0;
}
/***************************************************
* some functions used in main()
****************************************************/
PRECISION*MyAllocPrecision1D(int N,PRECISION value)
{
PRECISION*returnArray ;
N++;
int i ;
returnArray=(PRECISION*)malloc(N*sizeof(PRECISION));
if(returnArray==NULL)
{
printf("Sorry, but you don't have enough memory installed!");
exit(1);
}
for(i=0;i<N;i++)
{
returnArray[i]=value ;
}
return returnArray ;
}
PRECISION**MyAllocPrecision2D(int dimensionX,int dimensionY,PRECISION value)
{
PRECISION**returnArray ;
dimensionX++;
dimensionY++;
int i,j ;
returnArray=(PRECISION**)malloc(dimensionX*sizeof(PRECISION*));
if(returnArray==NULL)
{
printf("Sorry, but you don't have enough memory installed!");
exit(1);
}
for(i=0;i<dimensionX;i++)
{
returnArray[i]=(PRECISION*)malloc(dimensionY*sizeof(PRECISION));
if(returnArray[i]==NULL)
{
printf("Sorry, but you don't have enough memory installed!");
exit(1);
}
for(j=0;j<dimensionY;j++)
{
returnArray[i][j]=value ;
}
}
return returnArray ;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -