📄 air_pml_lsj.cpp
字号:
# include<math.h>
# include<iostream>
# include<fstream>
# include<iomanip>
# include<stdlib.h>
#include "str.h"
using namespace std;
#define C 299792458
#define Pi 3.1415926535897932384626433832795
#define epsilon0 8.85418e-12
#define mu0 1.256637e-6
double ***Alloc(int nx,int ny,int nz,double val)
{
double ***p;
int i,j,k;
p=new double **[nx];
for(i=0;i<nx;i++){
p[i]=new double *[ny];
for(j=0;j<ny;j++){
p[i][j]=new double[nz];
for(k=0;k<nz;k++){
p[i][j][k]=val;
}
}
}
return p;
}
double *Alloc(int nx,double val)
{
if (nx<=0) {
return NULL;
}
double *p;
int i;
p=new double [nx];
for(i=0;i<nx;i++){
p[i]=val;
}
return p;
}
int main ()
{
int I, J, K,It,Jt,Kt, N;
int i, j, k, is,js,ks,n;
double max_frequency;
double dx,delta_t;
//高斯脉冲的变量
double puls_gauss;
//PML边界条件的变量
double refmax; //最大反射系数
double temp,depth_pml;
double ca,cb,da,db,re,rm;
double sigmam_air,rhomam_air,eta0;
int n_pml=8; //完全匹配层网格数
eta0=sqrt(mu0/epsilon0);
refmax=1e-5;
max_frequency=20e9;
double srcT=1/(2*max_frequency);
double srcT0=3*srcT;
//-------空气网格数目----------
I=50;
J=20;
K=100;
//-------加入PML层后在空间所画网格数目----------
It=I+2*n_pml;
Jt=J+2*n_pml;
Kt=K+2*n_pml;
//激励源的位置
is=It/2;
js=Jt/2;
ks=Kt/2;
//
dx=(C/max_frequency)/15;//0.2e-3; //空间步长
delta_t=0.8*dx/(sqrt(3)*C); //时间步长
//------空气对应的电场场计算相关系数----------------
ca=1;
cb=delta_t/(epsilon0*dx);
da=1;
db=delta_t/(mu0*dx);
//allocate PML coefficients
double *sig=Alloc(n_pml,0);
double *rho=Alloc(n_pml,0);
double *ca0=Alloc(n_pml,0);
double *cb0=Alloc(n_pml,0);
double *da0=Alloc(n_pml,0);
double *db0=Alloc(n_pml,0);
double *sigf=Alloc(n_pml,0);
double *rhof=Alloc(n_pml,0);
double *sigh=Alloc(n_pml,0);
double *rhoh=Alloc(n_pml,0);
double *ca0f=Alloc(n_pml,0);
double *ca0h=Alloc(n_pml,0);
double *cb0f=Alloc(n_pml,0);
double *cb0h=Alloc(n_pml,0);
double *da0f=Alloc(n_pml,0);
double *da0h=Alloc(n_pml,0);
double *db0f=Alloc(n_pml,0);
double *db0h=Alloc(n_pml,0);
//------申请电磁场及其相关计算系数的动态数组-------------
double ***Ex= Alloc(It,Jt+1,Kt+1,0);
double ***Exy= Alloc(It,Jt+1,Kt+1,0);
double ***Exz= Alloc(It,Jt+1,Kt+1,0);
double ***caExy=Alloc(It,Jt+1,Kt+1,ca);
double ***cbExy=Alloc(It,Jt+1,Kt+1,cb);
double ***caExz=Alloc(It,Jt+1,Kt+1,ca);
double ***cbExz=Alloc(It,Jt+1,Kt+1,cb);
//
double ***Ey= Alloc(It+1,Jt,Kt+1,0);
double ***Eyx= Alloc(It+1,Jt,Kt+1,0);
double ***Eyz= Alloc(It+1,Jt,Kt+1,0);
double ***caEyx=Alloc(It+1,Jt,Kt+1,ca);
double ***cbEyx=Alloc(It+1,Jt,Kt+1,cb);
double ***caEyz=Alloc(It+1,Jt,Kt+1,ca);
double ***cbEyz=Alloc(It+1,Jt,Kt+1,cb);
//
double ***Ez= Alloc(It+1,Jt+1,Kt,0);
double ***Ezx= Alloc(It+1,Jt+1,Kt,0);
double ***Ezy= Alloc(It+1,Jt+1,Kt,0);
double ***caEzx=Alloc(It+1,Jt+1,Kt,ca);
double ***cbEzx=Alloc(It+1,Jt+1,Kt,cb);
double ***caEzy=Alloc(It+1,Jt+1,Kt,ca);
double ***cbEzy=Alloc(It+1,Jt+1,Kt,cb);
//
double ***Hx= Alloc(It+1,Jt,Kt,0);
double ***Hxy= Alloc(It+1,Jt,Kt,0);
double ***Hxz= Alloc(It+1,Jt,Kt,0);
double ***daHxy=Alloc(It+1,Jt,Kt,da);
double ***daHxz=Alloc(It+1,Jt,Kt,da);
double ***dbHxy=Alloc(It+1,Jt,Kt,db);
double ***dbHxz=Alloc(It+1,Jt,Kt,db);
//
double ***Hy= Alloc(It,Jt+1,Kt,0);
double ***Hyx= Alloc(It,Jt+1,Kt,0);
double ***Hyz= Alloc(It,Jt+1,Kt,0);
double ***daHyx=Alloc(It,Jt+1,Kt,da);
double ***daHyz=Alloc(It,Jt+1,Kt,da);
double ***dbHyx=Alloc(It,Jt+1,Kt,db);
double ***dbHyz=Alloc(It,Jt+1,Kt,db);
//
double ***Hz= Alloc(It,Jt,Kt+1,0);
double ***Hzx= Alloc(It,Jt,Kt+1,0);
double ***Hzy= Alloc(It,Jt,Kt+1,0);
double ***daHzx=Alloc(It,Jt,Kt+1,da);
double ***daHzy=Alloc(It,Jt,Kt+1,da);
double ***dbHzx=Alloc(It,Jt,Kt+1,db);
double ***dbHzy=Alloc(It,Jt,Kt+1,db);
//
//------------修正PML区域电磁场计算参数---------------
int m;
double retmp,rmtmp;
m=3;
depth_pml=dx*n_pml;
sigmam_air=-(m+1)*log(refmax)/(2*eta0*depth_pml); //-(m+1)ln[R(0)]/(2*d*eta) 3<=m<=4 //空气对应PML的最大电导
rhomam_air=sigmam_air*eta0*eta0; //空气对应PML的最大磁阻
retmp=delta_t/epsilon0;
rmtmp=delta_t/mu0;
for(n=0;n<n_pml;n++)
{
sig[n]=sigmam_air*pow((n+0.5)/(n_pml+0.5),2);
rho[n]=rhomam_air*pow((n+1)/(n_pml+0.5),2);
re=sig[n]*delta_t/epsilon0;
rm=rho[n]*delta_t/mu0;
ca0[n]=exp(-re); //空气对应的PML电磁场计算参数
cb0[n]=(1-ca0[n])/(sig[n]*dx);
da0[n]=exp(-rm);
db0[n]=(1-da0[n])/(rho[n]*dx);
}
//----------------PML区域电场相关参数------------------------
//-------------left and right PML---------------
//Eyx
for(j=0;j<Jt;j++)
{
for(k=0;k<=Kt;k++)
{
n=n_pml-1;
for(i=1;i<=n_pml;i++)
{
caEyx[i][j][k]=ca0[n];
cbEyx[i][j][k]=cb0[n];
n--;
}
n=0;
for(i=n_pml+I;i<It;i++)
{
caEyx[i][j][k]=ca0[n];
cbEyx[i][j][k]=cb0[n];
n++;
}
}
}
//Ezx
for(j=0;j<=Jt;j++)
{
for(k=0;k<Kt;k++)
{
n=n_pml-1;
for(i=1;i<=n_pml;i++)
{
caEzx[i][j][k]=ca0[n];
cbEzx[i][j][k]=cb0[n];
n--;
}
n=0;
for(i=n_pml+I;i<It;i++)
{
caEzx[i][j][k]=ca0[n];
cbEzx[i][j][k]=cb0[n];
n++;
}
}
}
//------------front and back PML----------------
//Exy
for(i=0;i<It;i++)
{
for(k=0;k<=Kt;k++)
{
n=n_pml-1;
for(j=1;j<=n_pml;j++)
{
caExy[i][j][k]=ca0[n];
cbExy[i][j][k]=cb0[n];
n--;
}
n=0;
for(j=n_pml+J;j<Jt;j++)
{
caExy[i][j][k]=ca0[n];
cbExy[i][j][k]=cb0[n];
n++;
}
}
}
//Ezy
for(i=0;i<=It;i++)
{
for(k=0;k<Kt;k++)
{
n=n_pml-1;
for(j=1;j<=n_pml;j++)
{
caEzy[i][j][k]=ca0[n];
cbEzy[i][j][k]=cb0[n];
n--;
}
n=0;
for(j=n_pml+J;j<Jt;j++)
{
caEzy[i][j][k]=ca0[n];
cbEzy[i][j][k]=cb0[n];
n++;
}
}
}
//------------down and up PML----------------
//Exz
for(i=0;i<It;i++)
{
for(j=0;j<=Jt;j++)
{
n=n_pml-1;
for(k=1;k<=n_pml;k++)
{
caExz[i][j][k]=ca0[n];
cbExz[i][j][k]=cb0[n];
n--;
}
n=0;
for(k=n_pml+K;k<Kt;k++)
{
caExz[i][j][k]=ca0[n];
cbExz[i][j][k]=cb0[n];
n++;
}
}
}
//Eyz
for(i=0;i<=It;i++)
{
for(j=0;j<Jt;j++)
{
n=n_pml-1;
for(k=1;k<=n_pml;k++)
{
caEyz[i][j][k]=ca0[n];
cbEyz[i][j][k]=cb0[n];
n--;
}
n=0;
for(k=n_pml+K;k<Kt;k++)
{
caEyz[i][j][k]=ca0[n];
cbEyz[i][j][k]=cb0[n];
n++;
}
}
}
//----------------PML区域磁场相关参数设置------------------------
//-------------left and right PML---------------
//Hyx
for(j=0;j<=Jt;j++)
{
for(k=0;k<Kt;k++)
{
n=n_pml-1;
for(i=0;i<n_pml;i++)
{
daHyx[i][j][k]=da0[n];
dbHyx[i][j][k]=db0[n];
n--;
}
n=0;
for(i=n_pml+I;i<It;i++)
{
daHyx[i][j][k]=da0[n];
dbHyx[i][j][k]=db0[n];
n++;
}
}
}
//Hzx
for(j=0;j<Jt;j++)
{
for(k=0;k<=Kt;k++)
{
n=n_pml-1;
for(i=0;i<n_pml;i++)
{
daHzx[i][j][k]=da0[n];
dbHzx[i][j][k]=db0[n];
n--;
}
n=0;
for(i=n_pml+I;i<It;i++)
{
daHzx[i][j][k]=da0[n];
dbHzx[i][j][k]=db0[n];
n++;
}
}
}
//------------Front and back PML----------------
//Hxy
for(i=0;i<=It;i++)
{
for(k=0;k<Kt;k++)
{
n=n_pml-1;
for(j=0;j<n_pml;j++)
{
daHxy[i][j][k]=da0[n];
dbHxy[i][j][k]=db0[n];
n--;
}
n=0;
for(j=n_pml+J;j<Jt;j++)
{
daHxy[i][j][k]=da0[n];
dbHxy[i][j][k]=db0[n];
n++;
}
}
}
//Hzy
for(i=0;i<It;i++)
{
for(k=0;k<=Kt;k++)
{
n=n_pml-1;
for(j=0;j<n_pml;j++)
{
daHzy[i][j][k]=da0[n];
dbHzy[i][j][k]=db0[n];
n--;
}
n=0;
for(j=n_pml+J;j<Jt;j++)
{
daHzy[i][j][k]=da0[n];
dbHzy[i][j][k]=db0[n];
n++;
}
}
}
//------------down and up PML-----------------------
//Hxz
for(i=0;i<=It;i++)
{
for(j=0;j<Jt;j++)
{
n=n_pml-1;
for(k=0;k<n_pml;k++)
{
daHxz[i][j][k]=da0[n];
dbHxz[i][j][k]=db0[n];
n--;
}
n=0;
for(k=n_pml+K;k<Kt;k++)
{
daHxz[i][j][k]=da0[n];
dbHxz[i][j][k]=db0[n];
n++;
}
}
}
//Hyz
for(i=0;i<It;i++)
{
for(j=0;j<=Jt;j++)
{
n=n_pml-1;
for(k=0;k<n_pml;k++)
{
daHyz[i][j][k]=da0[n];
dbHyz[i][j][k]=db0[n];
n--;
}
n=0;
for(k=n_pml+K;k<Kt;k++)
{
daHyz[i][j][k]=da0[n];
dbHyz[i][j][k]=db0[n];
n++;
}
}
}
//--------------------计算传播时间和高斯脉冲生成时间-------------------
N=600;
//----------------主循环--------------------
cout<<"主循环"<<endl;
cout<<"完成进度:"<<setw(4)<<setw(4)<<N<<endl;
putchar(8);putchar(8);putchar(8);putchar(8);putchar(8);
putchar(8);putchar(8);putchar(8);putchar(8);
//save source value
FILE *pSrcfile;
CStr str;
str.Format("srcData.txt",n);
pSrcfile=fopen(str,"w");
if (pSrcfile==NULL) { return 0;
}
double t;
for(n=1;n<=N;n++)
{
//----------------计算高斯脉冲-----------------------
t=n*delta_t;
temp= (t-srcT0)/srcT;
puls_gauss=exp(-temp*temp);
cout<<"step="<<n<<endl;
//------------------电场的循环迭代---------------------
for(i=0;i<It;i++)
{
for(j=1;j<Jt;j++)
{
for(k=1;k<Kt;k++)
{
Exy[i][j][k]=caExy[i][j][k]*Exy[i][j][k]+cbExy[i][j][k]*(Hz[i][j][k]-Hz[i][j-1][k]);
Exz[i][j][k]=caExz[i][j][k]*Exz[i][j][k]+cbExz[i][j][k]*(-Hy[i][j][k]+Hy[i][j][k-1]);
Ex[i][j][k]=Exy[i][j][k]+Exz[i][j][k];
}
}
}
//boundary
for(k=0;k<=Kt;k+=Kt){
for(i=0;i<It;i++){
for(j=0;j<=Jt;j++){
Ex[i][j][k]=0;
}
}
}
for(j=0;j<=Jt;j+=Jt){
for(i=0;i<It;i++){
for(k=0;k<=Kt;k++){
Ex[i][j][k]=0;
}
}
}
////
for(i=1;i<It;i++)
{
for(j=0;j<Jt;j++)
{
for(k=1;k<Kt;k++)
{
Eyx[i][j][k]=caEyx[i][j][k]*Eyx[i][j][k]+cbEyx[i][j][k]*(-Hz[i][j][k]+Hz[i-1][j][k]);
Eyz[i][j][k]=caEyz[i][j][k]*Eyz[i][j][k]+cbEyz[i][j][k]*(Hx[i][j][k]-Hx[i][j][k-1]);
Ey[i][j][k]=Eyx[i][j][k]+Eyz[i][j][k];
}
}
}
//boundary
for(k=0;k<=Kt;k+=Kt){
for(i=0;i<=It;i++){
for(j=0;j<Jt;j++){
Ey[i][j][k]=0;
}
}
}
for(i=0;i<=It;i+=It){
for(j=0;j<Jt;j+=Jt){
for(k=0;k<=Kt;k++){
Ey[i][j][k]=0;
}
}
}
////
for(i=1;i<It;i++)
{
for(j=1;j<Jt;j++)
{
for(k=0;k<Kt;k++)
{
Ezx[i][j][k]=caEzx[i][j][k]*Ezx[i][j][k]+cbEzx[i][j][k]*(Hy[i][j][k]-Hy[i-1][j][k]);
Ezy[i][j][k]=caEzy[i][j][k]*Ezy[i][j][k]+cbEzy[i][j][k]*(-Hx[i][j][k]+Hx[i][j-1][k]);
Ez[i][j][k]=Ezx[i][j][k]+Ezy[i][j][k];
}
}
}
//boundary
for(j=0;j<=Jt;j+=Jt){
for(i=0;i<=It;i++){
for(k=0;k<Kt;k++){
Ez[i][j][k]=0;
}
}
}
for(i=0;i<=It;i+=It){
for(j=0;j<=Jt;j+=Jt){
for(k=0;k<Kt;k++){
Ez[i][j][k]=0;
}
}
}
////
//---------------------------加入高斯脉冲-----------------
Ex[is][js][ks]=delta_t*puls_gauss/epsilon0+Ex[is][js][ks];
//-----------------------------循环迭代H-----------------------------
//hx
for(i=0;i<=It;i++){
for(j=0;j<Jt;j++){
for(k=0;k<Kt;k++){
Hxy[i][j][k]=daHxy[i][j][k]*Hxy[i][j][k]-dbHxy[i][j][k]*(Ez[i][j+1][k]-Ez[i][j][k]);
Hxz[i][j][k]=daHxz[i][j][k]*Hxz[i][j][k]+dbHxz[i][j][k]*(Ey[i][j][k+1]-Ey[i][j][k]);
Hx[i][j][k]=Hxy[i][j][k]+Hxz[i][j][k];
}
}
}
//hy
for(i=0;i<It;i++){
for(j=0;j<=Jt;j++){
for(k=0;k<Kt;k++){
Hyx[i][j][k]=daHyx[i][j][k]*Hyx[i][j][k]+dbHyx[i][j][k]*(Ez[i+1][j][k]-Ez[i][j][k]);
Hyz[i][j][k]=dbHyz[i][j][k]*Hyz[i][j][k]-dbHyz[i][j][k]*(Ex[i][j][k+1]-Ex[i][j][k]);
Hy[i][j][k]=Hyx[i][j][k]+Hyz[i][j][k];
}
}
}
//hz
for(i=0;i<It;i++){
for(j=0;j<Jt;j++){
for(k=0;k<=Kt;k++){
Hzx[i][j][k]=daHzx[i][j][k]*Hzx[i][j][k]-dbHzx[i][j][k]*(Ey[i+1][j][k]-Ey[i][j][k]);
Hzy[i][j][k]=daHzy[i][j][k]*Hzy[i][j][k]+dbHzy[i][j][k]*(Ex[i][j+1][k]-Ex[i][j][k]);
Hz[i][j][k]=Hzx[i][j][k]+Hzy[i][j][k];
}
}
}
fprintf(pSrcfile,"%10.5E\t%10.5E\n",t,Ex[It/2+10][Jt/2][Kt/2]);
/*
if (fmod(n,20)==0) {
FILE *pfile;
CStr str;
str.Format("data%d.txt",n);
pfile=fopen(str,"w");
for(j=0;j<=Jt;j++){
for(k=0;k<=Kt;k++){
fprintf(pfile,"%10.5E\t",Ex[is+5][j][k]);
}
fprintf(pfile,"\n");
}
fclose(pfile);
}
*/
// cout<<setw(4)<<n;
// putchar(8);putchar(8);putchar(8);putchar(8);
}
fclose(pSrcfile);
//-----------------------------------结束-----------------------------------
return 0;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -