#pragma hdrstop
//---------------------------------------------------------------------------
//#pragma argsused
#include
#include
#include
#include
#include
#include
//#include
#define PI 3.1415926535
#define e 2.718281828
float space2d(int nr, int nc);
float ***space3d(int nr, int ny, int nc);
void free_space2d(float **a, int nr);
void free_space3d(float ***b, int nr, int ny);
void wfile(char filename[], float **data, int nr, int nc);
void wfile3d(char filename[], float **data, int nr, int ny, int nc);
void wfile4d(char filename[], float **data, int nr, int ny, int nc);
void create_model(float **vp, float **vs, float **rho, float **vf, float **rhof, float **lamda, float **lamda2u, float **mu, float*por, int nr, int ny, int nc);
float **extmodel(float **init_model, int nr, int ny, int nc, int np);
int main()
{
//给定参数
int NX = 100; //x方向网格点数%%%%%
int NY = 100; //y方向网格点数%%%%%
int NZ = 300; //z方向网格点数%%%%
int NP = 20; //pml层网格点数
int i = 0;
int sx = NX / 2 + NP; //震源坐标点号
int sy = NY / 2 + NP;
int sz = NZ /6 + NP;
int NX_ext = NX + 2 * NP; //
int NY_ext = NY + 2 * NP;
int NZ_ext = NZ + 2 * NP;
int NT = 8000; //空间步长5%%%
float H = 0.01f;
float RC = 0.000001f;
float DP = NP*H;
// double DT=0.0002; //时间步长
double DT = 1.0*pow(10, -6); //时间步长%%%%
// double F0=30.0; //震源主频
double F0 = 3.0*pow(10, 3); //震源主频,主频太高会造成频散
//double T0=1.2/F0;
float T0 = 1.0 / F0;
double Vpmax = 4270.0; //模型最大纵波速度,用于稳定性计算
double Vpmin = 1500.0; //模型最小纵波速度,用于控制数值频散
double Vsmax = 2650.0;
double Vsmin = 1500.0;
float ***vs; //初始模型横波速度
float ***vp; //初始模型纵波速度
float ***rho; //初始模型密度
float ***vf; //孔隙流体速度
float ***rhof; //孔隙流体密度
float ***mu; //剪切模量
float ***lamda; //拉梅系数
float ***lamda2u;
//float **Q;
//float **R;
float ***por;
// float **D11,**D12,**D22;
// float **b;
float **sis_x; //地震记录x分量 波形记录-二维数据
float **sis_y; //地震记录y分量
float **sis_z; //地震记录z分量
float ***vx;
float ***vxx;
float ***vy;
float ***vz;
/*float **Vx_x;//流相速度
float **Vx_z;
float **Vx;
float **Vz_x;
float **Vz_z;
float **Vz;
*/
float ***txx;//应力分量
float ***tyy;//应力分量
float ***tzz;
float ***txy;
float ***txz;
float ***tyz;
float ***source;
/*float **ssx;//流体应力
float **ssz;
float **dxi,**dxi2;
float **dzj,**dzj2;
*/
// double Kb,Ks,Kf,a,D;
float tt, x, y, z, xoleft, xoright, yoleft, yoright, zoleft, zoright, d0=0.0f;
float v0, y5, y6 , y7 , y8, y9 , y10, y11, y12, y13, rho_tempx, rho_tempy, rho_tempz, muxz, muxy, muyz, lamda2u_temp, lamda_temp;
// double D11_tempx,D12_tempx,D22_tempx,R_temp,D11_tempz,D12_tempz,D22_tempz;
float xx1, xy1, xz1, yx2 , yy2, yz2 , zx3, zy3, zz3 , best_dt;
int ix, iy, iz, it;
vs = space3d(NZ, NY, NX);
vp = space3d(NZ, NY, NX);
rho = space3d(NZ, NY, NX);
vf = space3d(NZ, NY, NX);
rhof = space3d(NZ, NY, NX);
//Q=space2d(NZ,NX);
//R=space2d(NZ,NX);
por = space3d(NZ, NY, NX);
//D11=space2d(NZ,NX);D12=space2d(NZ,NX);D22=space2d(NZ,NX);
// b=space2d(NZ,NX);
mu = space3d(NZ, NY, NX);
lamda = space3d(NZ, NY, NX);
lamda2u = space3d(NZ, NY, NX);
create_model(vp, vs, rho, vf, rhof, lamda, lamda2u, mu, por, NZ, NY, NX);
char vp_name[] = "vp_ext.dat";
float ***vs_ext;
float ***vp_ext;
float ***rho_ext;
float ***mu_ext;
float ***lamda_ext;
float ***lamda2u_ext;
// float ***Q_ext;
// float ***R_ext;
float ***vf_ext;
float ***rhof_ext;
float ***por_ext;
float ***dxi;
float ***dxi2;
float ***dyj;
float ***dyj2;
float ***dzk;
float ***dzk2;
//float ***b_ext;
// float **D11_ext,**D12_ext,**D22_ext;
//数据扩充，加上吸收边界
vs_ext = extmodel(vs, NZ, NY, NX, NP);
vp_ext = extmodel(vp, NZ, NY, NX, NP);
rho_ext = extmodel(rho, NZ, NY, NX, NP);
mu_ext = extmodel(mu, NZ, NY, NX, NP);
lamda_ext = extmodel(lamda, NZ, NY, NX, NP);
lamda2u_ext = extmodel(lamda2u, NZ, NY, NX, NP);
vf_ext = extmodel(vf, NZ, NY, NX, NP);
rhof_ext = extmodel(rhof, NZ, NY, NX, NP);
//Q_ext=extmodel(Q,NZ,NX,NP);
//R_ext=extmodel(R,NZ,NX,NP);
por_ext = extmodel(por, NZ, NY, NX, NP);
// b_ext=extmodel(b,NZ,NX,NP);
// D11_ext=extmodel(D11,NZ,NX,NP);D12_ext=extmodel(D12,NZ,NX,NP);
// D22_ext=extmodel(D22,NZ,NX,NP);
wfile3d(vp_name, vp_ext, NZ_ext, NY_ext, NX_ext); //建立文件
//申请空间
vx = space3d(NZ_ext, NY_ext, NX_ext);
vy = space3d(NZ_ext, NY_ext, NX_ext);
vxx= space3d(NT/200,NZ_ext/2, NX_ext/2);
vz = space3d(NZ_ext, NY_ext, NX_ext);
/*Vx_x=space2d(NZ_ext,NX_ext);
Vx_z=space2d(NZ_ext,NX_ext);
Vx=space2d(NZ_ext,NX_ext);
Vz_x=space2d(NZ_ext,NX_ext);
Vz_z=space2d(NZ_ext,NX_ext);
Vz=space2d(NZ_ext,NX_ext);*/
txx = space3d(NZ_ext, NY_ext, NX_ext);
tyy = space3d(NZ_ext, NY_ext, NX_ext);
tzz = space3d(NZ_ext, NY_ext, NX_ext);
txy = space3d(NZ_ext, NY_ext, NX_ext);
txz = space3d(NZ_ext, NY_ext, NX_ext);
tyz = space3d(NZ_ext, NY_ext, NX_ext);
source= space3d(NZ_ext, NY_ext, NX_ext);
/* tzz_x=space2d(NZ_ext,NX_ext);
tzz_z=space2d(NZ_ext,NX_ext);
txz_x=space2d(NZ_ext,NX_ext);
txz_z=space2d(NZ_ext,NX_ext);
ssx=space2d(NZ_ext,NX_ext);
ssz=space2d(NZ_ext,NX_ext);*/
dxi = space3d(NZ_ext, NY_ext, NX_ext);
dxi2 = space3d(NZ_ext, NY_ext, NX_ext);
dyj = space3d(NZ_ext, NY_ext, NX_ext);
dyj2 = space3d(NZ_ext, NY_ext, NX_ext);
dzk = space3d(NZ_ext, NY_ext, NX_ext);
dzk2 = space3d(NZ_ext, NY_ext, NX_ext);
char dxi_name[] = "dxi.dat";
char dzk_name[] = "dzk.dat";
//(xoleft,xoright)模型的物理边界—不算PML
xoleft = DP;
xoright = (NX_ext-1)*H - DP;
yoleft = DP;
yoright = (NY_ext-1 )*H - DP;
zoleft = DP;
zoright = (NZ_ext-1 )*H - DP;
//用于对vx_x[iz][ix]等的求解，加入吸收边界，使数值衰减
for (iz = 0; iz < NZ_ext; iz++)
{
for (iy = 0; iy < NY_ext; iy++)
{
for (ix = 0; ix < NX_ext; ix++)
{
x = ix*H; y = iy*H; z = iz*H;
if ((x>=0&&x < xoleft)&&(y>=0&&y<NY_ext*H)&& (z >= 0 && z<NZ_ext*H))
{
v0 = 1500.0f;
d0 = 3.0*v0*log(1.0 / RC) / (2.0*DP);
dxi[iz][iy][ix] = d0*pow(((xoleft - x) / DP), 2);
dxi2[iz][iy][ix] = d0*pow(((xoleft - x - 0.5*H) / DP), 2);
}
if ((x >= xoright && x < NX_ext*H) && (y >= 0 && y<NY_ext*H) && (z >= 0 && z<NZ_ext*H))
{
v0 = 1500.0f;
d0 = 3.0*v0*log(1.0 / RC) / (2.0*DP);
dxi[iz][iy][ix] = d0*pow(((x - xoright) / DP), 2);
dxi2[iz][iy][ix] = d0*pow(((x + 0.5*H - xoright) / DP), 2);
}
{
if ((x >= 0 && x<NX_ext*H) && (y >= 0 && y < yoleft) && (z >= 0 && z<NZ_ext*H))
{
v0 = 1500.0f;
d0 = 3.0*v0*log(1.0 / RC) / (2.0*DP);
//d0=3.0*v0*(8.0/15.0-3.0/100.0*NP+NP*NP/1500.0)/H;
dyj[iz][iy][ix] = d0*pow(((yoleft - y) / DP), 2);
dyj2[iz][iy][ix] = d0*pow(((yoleft - y - 0.5*H) / DP), 2);
}
if ((x >= 0 && x<NX_ext*H) && (y >= yoright && y < NY_ext*H) &&(z >= 0 && z<NZ_ext*H))
{
v0 = 1500.0f;
d0 = 3.0*v0*log(1.0 / RC) / (2.0*DP);
dyj[iz][iy][ix] = d0*pow(((y - yoright) / DP), 2);
dyj2[iz][iy][ix] = d0*pow(((y + 0.5*H - yoright) / DP), 2);
}
if ((x >= 0 && x<NX_ext*H) && (y >= 0 && y<NY_ext*H) && (z >= 0 && z < zoleft) )
{
v0 = 1500.0f;
d0 = 3.0*v0*log(1.0 / RC) / (2.0*DP);
//d0=3.0*v0*(8.0/15.0-3.0/100.0*NP+NP*NP/1500.0)/H;
dzk[iz][iy][ix] = d0*pow(((zoleft - z) / DP), 2);
dzk2[iz][iy][ix] = d0*pow(((zoleft - z - 0.5*H) / DP), 2);
}
if ((x >= 0 && x<NX_ext*H) && (y >= 0 && y<NY_ext*H) && (z >= zoright && z < NZ_ext*H) )
{
v0 = 1500.0f;
d0 = 3.0*v0*log(1.0 / RC) / (2.0*DP);
dzk[iz][iy][ix] = d0*pow(((z - zoright) / DP), 2);
dzk2[iz][iy][ix] = d0*pow(((z + 0.5*H - zoright) / DP), 2);
}
}
}
}
}
/*for (iz = 0; iz<NZ_ext; iz++)
{
for (iy = 0; iy<NY_ext; iy++)
{
y = iy*H;
for (ix = 0; ix<NX_ext; ix++)
{
if (y < yoleft)
{
v0 = 1500.0;
d0 = 3.0*v0*log(1.0 / RC) / (2.0*DP);
//d0=3.0*v0*(8.0/15.0-3.0/100.0*NP+NP*NP/1500.0)/H;
dyj[iz][iy][ix] = d0*pow(((yoleft - y) / DP), 2);
dyj2[iz][iy][ix] = d0*pow(((yoleft - y - 0.5*H) / DP), 2);
}
else if (y >= 0.9999*yoright)
{
v0 = 1500.0;
d0 = 3.0*v0*log(1.0 / RC) / (2.0*DP);
dyj[iz][iy][ix] = d0*pow(((y - yoright) / DP), 2);
dyj2[iz][iy][ix] = d0*pow(((y + 0.5*H - yoright) / DP), 2);
}
else
{
dyj[iz][iy][ix] = 0.; dyj2[iz][iy][ix] = 0.;
}
}
}
}
for (iz = 0; iz<NZ_ext; iz++)
{
z = iz*H;
for (iy = 0; iy<NY_ext; iy++)
{
for (ix = 0; ix<NX_ext; ix++)
{
if (z < zoleft)
{
v0 = 1500.0;
d0 = 3.0*v0*log(1.0 / RC) / (2.0*DP);
//d0=3.0*v0*(8.0/15.0-3.0/100.0*NP+NP*NP/1500.0)/H;
dzk[iz][iy][ix] = d0*pow(((zoleft - z) / DP), 2);
dzk2[iz][iy][ix] = d0*pow(((zoleft - z - 0.5*H) / DP), 2);
}
else if (z >= 0.9999*zoright)
{
v0 = 1500.0;
d0 = 3.0*v0*log(1.0 / RC) / (2.0*DP);
dzk[iz][iy][ix] = d0*pow(((z - zoright) / DP), 2);
dzk2[iz][iy][ix] = d0*pow(((z + 0.5*H - zoright) / DP), 2);
}
else
{
dzk[iz][iy][ix] = 0.; dzk2[iz][iy][ix] = 0.;
}
}
}
}*/
wfile3d(dxi_name, dxi, NZ_ext, NY_ext, NX_ext);
wfile3d(dzk_name, dzk, NZ_ext, NY_ext, NX_ext);//建立文件
//开始时间递推
float a1 = 9.0 / 8.0;
float a2 = -1.0 / 24.0;
//pml初始参数
float xFHalfTemp1;
float xFHalfTemp2 , xFIntTemp1, xFIntTemp2, yFHalfTemp1, yFHalfTemp2, yFIntTemp1, yFIntTemp2, zFHalfTemp1, zFHalfTemp2, zFIntTemp1, zFIntTemp2 ;
float ***pmlxSxx, ***pmlySxy, ***pmlzSxz, ***pmlxSxy, ***pmlySyy, ***pmlzSyz, ***pmlxSxz, ***pmlySyz, ***pmlzSzz;
pmlxSxx= space3d(NZ_ext, NY_ext, NX_ext);
pmlySxy= space3d(NZ_ext, NY_ext, NX_ext);
pmlzSxz= space3d(NZ_ext, NY_ext, NX_ext);
pmlxSxy= space3d(NZ_ext, NY_ext, NX_ext);
pmlySyy= space3d(NZ_ext, NY_ext, NX_ext);
pmlzSyz= space3d(NZ_ext, NY_ext, NX_ext);
pmlxSxz= space3d(NZ_ext, NY_ext, NX_ext);
pmlySyz= space3d(NZ_ext, NY_ext, NX_ext);
pmlzSzz= space3d(NZ_ext, NY_ext, NX_ext);
float ***pmlxVx, ***pmlyVy, ***pmlzVz, ***pmlzVx, ***pmlxVz, ***pmlyVx, ***pmlxVy, ***pmlzVy, ***pmlyVz;
pmlxVx = space3d(NZ_ext, NY_ext, NX_ext);
pmlyVy = space3d(NZ_ext, NY_ext, NX_ext);
pmlzVz = space3d(NZ_ext, NY_ext, NX_ext);
pmlzVx = space3d(NZ_ext, NY_ext, NX_ext);
pmlxVz = space3d(NZ_ext, NY_ext, NX_ext);
pmlyVx = space3d(NZ_ext, NY_ext, NX_ext);
pmlxVy = space3d(NZ_ext, NY_ext, NX_ext);
pmlzVy = space3d(NZ_ext, NY_ext, NX_ext);
pmlyVz = space3d(NZ_ext, NY_ext, NX_ext);
float ***DxSxxPre, ***DySxyPre, ***DzSxzPre, ***DxSxyPre, ***DySyyPre, ***DzSyzPre, ***DxSxzPre, ***DySyzPre, ***DzSzzPre;
DxSxxPre = space3d(NZ_ext, NY_ext, NX_ext);
DySxyPre = space3d(NZ_ext, NY_ext, NX_ext);
DzSxzPre = space3d(NZ_ext, NY_ext, NX_ext);
DxSxyPre = space3d(NZ_ext, NY_ext, NX_ext);
DySyyPre = space3d(NZ_ext, NY_ext, NX_ext);
DzSyzPre = space3d(NZ_ext, NY_ext, NX_ext);
DxSxzPre = space3d(NZ_ext, NY_ext, NX_ext);
DySyzPre = space3d(NZ_ext, NY_ext, NX_ext);
DzSzzPre = space3d(NZ_ext, NY_ext, NX_ext);
float ***DxVxPre, ***DyVyPre, ***DzVzPre, ***DzVxPre, ***DxVzPre, ***DyVxPre, ***DxVyPre, ***DzVyPre, ***DyVzPre;
DxVxPre = space3d(NZ_ext, NY_ext, NX_ext);
DyVyPre = space3d(NZ_ext, NY_ext, NX_ext);
DzVzPre = space3d(NZ_ext, NY_ext, NX_ext);
DzVxPre = space3d(NZ_ext, NY_ext, NX_ext);
DxVzPre = space3d(NZ_ext, NY_ext, NX_ext);
DyVxPre = space3d(NZ_ext, NY_ext, NX_ext);
DxVyPre = space3d(NZ_ext, NY_ext, NX_ext);
DzVyPre = space3d(NZ_ext, NY_ext, NX_ext);
DyVzPre = space3d(NZ_ext, NY_ext, NX_ext);
//这里的代码应该修正一下，应该从数据里面提取最大最小速度
//控制网格频散
if (Vsmin / (F0*H) < 15)
printf("1\n");
//检查稳定性条件
best_dt = 6.0*H / (7.0*sqrt(2.0)*Vpmax);
if (DT >= best_dt)
printf("2 %f\n", best_dt);
clock_t start, finish;
start = clock();
sis_x = space2d(NZ, NT);
sis_y = space2d(NZ, NT);
sis_z = space2d(NZ, NT);
omp_set_num_threads(16);
#pragma omp parallel private(ix, iy,iz,rho_tempx,rho_tempy,rho_tempz,xx1 , xy1, xz1, yx2, yy2 , yz2, zx3,zy3, zz3,xFHalfTemp1,yFHalfTemp1,zFHalfTemp1,xFIntTemp1,yFIntTemp1,zFIntTemp1,xFHalfTemp2,yFHalfTemp2,zFHalfTemp2,xFIntTemp2,yFIntTemp2,zFIntTemp2,muxz,muyz,muxy,lamda2u_temp,lamda_temp,y5,y6,y7,y8,y9,y10,y11,y12,y13)
{//$OMP PARALLEL DEFAULT(SHARED) PRIVATE(x, y, z, &
//$OMP VxSource, VySource, VzSource, SxxSource, SyySource, SzzSource, SxySource, SxzSource, SyzSource, &
//$OMP DenOReciprocal, DenXReciprocal, DenYReciprocal, DenZReciprocal, LambdaO, MuO, MuXY, MuXZ, MuYZ, &
//$OMP DxSxx, DySxy, DzSxz, DxSxy, DySyy, DzSyz, DxSxz, DySyz, DzSzz, &
//$OMP xAttFactorHalfTemp1, xAttFactorIntTemp1, yAttFactorHalfTemp1, yAttFactorIntTemp1, zAttFactorHalfTemp1, zAttFactorIntTemp1, &
//$OMP xAttFactorHalfTemp2, xAttFactorIntTemp2, yAttFactorHalfTemp2, yAttFactorIntTemp2, zAttFactorHalfTemp2, zAttFactorIntTemp2, &
//$OMP BoolTemp1, BoolTemp2, BoolTemp3, BoolTemp4, BoolTemp5, BoolTemp6, &
//$OMP DxVx, DyVy, DzVz, DxVz, DzVx, DxVy, DyVx, DyVz, DzVy)
for (it = 0; it < NT; it++)
{
tt = it*DT;
if (it % 100 == 0) printf("it=%d\n", it);
//计算质点速度v
#pragma omp for
for (iz = 2; iz < NZ_ext - 2; iz++)
for (iy = 2; iy < NY_ext - 2; iy++)
for (ix = 2; ix < NX_ext - 2; ix++)
{
rho_tempx = 2 / (rho_ext[iz][iy][ix] + rho_ext[iz][iy][ix + 1]);
rho_tempy = 2 / (rho_ext[iz][iy][ix] + rho_ext[iz][iy + 1][ix]);
rho_tempz = 2 / (rho_ext[iz][iy][ix] + rho_ext[iz + 1][iy][ix]);
//4阶
/*xx1 = ((a1*(txx[iz][iy][ix + 1] - txx[iz][iy][ix]) +
a2*(txx[iz][iy][ix + 2] - txx[iz][iy][ix - 1]))) / H;
xy1 = ((a1*(txy[iz][iy][ix] - txy[iz][iy - 1][ix]) +
a2*(txy[iz][iy + 1][ix] - txy[iz][iy - 2][ix]))) / H;
xz1 = ((a1*(txz[iz][iy][ix] - txz[iz - 1][iy][ix]) +
a2*(txz[iz + 1][iy][ix] - txz[iz - 2][iy][ix]))) / H;
yx2 = ((a1*(txy[iz][iy][ix] - txy[iz][iy][ix - 1]) +
a2*(txy[iz][iy][ix + 1] - txy[iz][iy][ix - 2]))) / H;
yy2 = ((a1*(tyy[iz][iy + 1][ix] - tyy[iz][iy][ix]) +
a2*(tyy[iz][iy + 2][ix] - tyy[iz][iy - 1][ix]))) / H;
yz2 = ((a1*(tyz[iz][iy][ix] - tyz[iz - 1][iy][ix]) +
a2*(tyz[iz + 1][iy][ix] - tyz[iz - 2][iy][ix]))) / H;
zx3 = ((a1*(txz[iz][iy][ix] - txz[iz][iy][ix - 1]) +
a2*(txz[iz][iy][ix + 1] - txz[iz][iy][ix - 2]))) / H;
zy3 = ((a1*(tyz[iz][iy][ix] - tyz[iz][iy - 1][ix]) +
a2*(tyz[iz][iy + 1][ix] - tyz[iz][iy - 2][ix]))) / H;
zz3 = ((a1*(tzz[iz + 1][iy][ix] - tzz[iz][iy][ix]) +
a2*(tzz[iz + 2][iy][ix] - tzz[iz - 1][iy][ix]))) / H;*/
//二阶
xx1 = (txx[iz][iy][ix + 1] - txx[iz][iy][ix]) / H;
xy1 = (txy[iz][iy][ix] - txy[iz][iy - 1][ix]) / H;
xz1 = (txz[iz][iy][ix] - txz[iz - 1][iy][ix]) / H;
yx2 = (txy[iz][iy][ix] - txy[iz][iy][ix - 1]) / H;
yy2 = (tyy[iz][iy + 1][ix] - tyy[iz][iy][ix]) / H;
yz2 = (tyz[iz][iy][ix] - tyz[iz - 1][iy][ix]) / H;
zx3 = (txz[iz][iy][ix] - txz[iz][iy][ix - 1]) / H;
zy3 = (tyz[iz][iy][ix] - tyz[iz][iy - 1][ix]) / H;
zz3 = (tzz[iz + 1][iy][ix] - tzz[iz][iy][ix]) / H;
/*x = ix*H; y = iy*H; z = iz*H;
if(((x >= 0 && x < xoleft) && (y >= 0 && y<NY_ext*H) && (z >= 0 && z<NZ_ext*H))
||((x >= xoright && x < NX_ext*H) && (y >= 0 && y<NY_ext*H) && (z >= 0 && z<NZ_ext*H))
||((x >= 0 && x<NX_ext*H) && (y >= 0 && y < yoleft) && (z >= 0 && z<NZ_ext*H))
||((x >= 0 && x<NX_ext*H) && (y >= yoright && y < NY_ext*H) && (z >= 0 && z<NZ_ext*H))
||((x >= 0 && x<NX_ext*H) && (y >= 0 && y<NY_ext*H) && (z >= 0 && z < zoleft))
||((x >= 0 && x<NX_ext*H) && (y >= 0 && y<NY_ext*H) && (z >= zoright && z < NZ_ext*H)))*/
if((ix<NP||ix>=(NX_ext-1-NP))||(iy<NP||iy>=(NY_ext-1-NP))||(iz<NP||iz>=(NZ_ext-1-NP)))
{
xFHalfTemp1 = exp(-(dxi2[iz][iy][ix])*DT); xFIntTemp1 = exp(-(dxi[iz][iy][ix])*DT);
yFHalfTemp1 = exp(-(dyj2[iz][iy][ix])*DT); yFIntTemp1 = exp(-(dyj[iz][iy][ix])*DT);
zFHalfTemp1 = exp(-(dzk2[iz][iy][ix])*DT); zFIntTemp1 = exp(-(dzk[iz][iy][ix])*DT);
xFHalfTemp2 = -DT*dxi2[iz][iy][ix] * 0.5; xFIntTemp2 = -DT*dxi[iz][iy][ix] * 0.5;
yFHalfTemp2 = -DT*dyj2[iz][iy][ix] * 0.5; yFIntTemp2 = -DT*dyj[iz][iy][ix] * 0.5;
zFHalfTemp2 = -DT*dzk2[iz][iy][ix] * 0.5; zFIntTemp2 = -DT*dzk[iz][iy][ix] * 0.5;
pmlxSxx[iz][iy][ix] = xFHalfTemp1*pmlxSxx[iz][iy][ix] + xFHalfTemp2*(xFHalfTemp1*DxSxxPre[iz][iy][ix] + xx1);
pmlySxy[iz][iy][ix] = yFIntTemp1*pmlySxy[iz][iy][ix] + yFIntTemp2*(yFIntTemp1*DySxyPre[iz][iy][ix] + xy1);
pmlzSxz[iz][iy][ix] = zFIntTemp1*pmlzSxz[iz][iy][ix] + zFIntTemp2*(zFIntTemp1*DzSxzPre[iz][iy][ix] + xz1);
DxSxxPre[iz][iy][ix] = xx1; DySxyPre[iz][iy][ix] = xy1; DzSxzPre[iz][iy][ix] = xz1;
xx1 = xx1 + pmlxSxx[iz][iy][ix];
xy1 = xy1 + pmlySxy[iz][iy][ix];
xz1 = xz1 + pmlzSxz[iz][iy][ix];
pmlxSxy[iz][iy][ix] = xFIntTemp1*pmlxSxy[iz][iy][ix] + xFIntTemp2*(xFIntTemp1*DxSxyPre[iz][iy][ix] + yx2);
pmlySyy[iz][iy][ix] = yFHalfTemp1*pmlySyy[iz][iy][ix] + yFHalfTemp2*(yFHalfTemp1*DySyyPre[iz][iy][ix] + yy2);
pmlzSyz[iz][iy][ix] = zFIntTemp1*pmlzSyz[iz][iy][ix] + zFIntTemp2*(zFIntTemp1*DzSyzPre[iz][iy][ix] + yz2);
DxSxyPre[iz][iy][ix] = yx2; DySyyPre[iz][iy][ix] = yy2; DzSyzPre[iz][iy][ix] = yz2;
yx2 = yx2 + pmlxSxy[iz][iy][ix];
yy2 = yy2 + pmlySyy[iz][iy][ix];
yz2 = yz2 + pmlzSyz[iz][iy][ix];
pmlxSxz[iz][iy][ix] = xFIntTemp1*pmlxSxz[iz][iy][ix] + xFIntTemp2*(xFIntTemp1*DxSxzPre[iz][iy][ix] + zx3);
pmlySyz[iz][iy][ix] = yFIntTemp1*pmlySyz[iz][iy][ix] + yFIntTemp2*(yFIntTemp1*DySyzPre[iz][iy][ix] + zy3);
pmlzSzz[iz][iy][ix] = zFHalfTemp1*pmlzSzz[iz][iy][ix] + zFHalfTemp2*(zFHalfTemp1*DzSzzPre[iz][iy][ix] + zz3);
DxSxzPre[iz][iy][ix] = zx3; DySyzPre[iz][iy][ix] = zy3; DzSzzPre[iz][iy][ix] = zz3;
zx3 = zx3 + pmlxSxz[iz][iy][ix];
zy3 = zy3 + pmlySyz[iz][iy][ix];
zz3 = zz3 + pmlzSzz[iz][iy][ix];
}
vx[iz][iy][ix] = vx[iz][iy][ix] + (DT*rho_tempx)*(xx1 + xy1 + xz1);
vy[iz][iy][ix] = vy[iz][iy][ix] + (DT*rho_tempy)*(yx2 + yy2 + yz2);
vz[iz][iy][ix] = vz[iz][iy][ix] + (DT*rho_tempz)*(zx3 + zy3 + zz3);
}
#pragma omp barrier
//计算应力
#pragma omp for
for (iz = 2; iz < NZ_ext - 2; iz++)
for (iy = 2; iy < NY_ext - 2; iy++)
for (ix = 2; ix < NX_ext - 2; ix++)
{
if (mu_ext[iz][iy][ix] == 0 || mu_ext[iz][iy][ix + 1] == 0 || mu_ext[iz + 1][iy][ix] == 0 || mu_ext[iz + 1][iy][ix + 1] == 0)
muxz = 0.0;
else
muxz = 1 / (0.25*(1 / mu_ext[iz][iy][ix] + 1 / mu_ext[iz][iy][ix + 1] + 1 / mu_ext[iz + 1][iy][ix] + 1 / mu_ext[iz + 1][iy][ix + 1]));
if (mu_ext[iz][iy][ix] == 0 || mu_ext[iz][iy + 1][ix] == 0 || mu_ext[iz + 1][iy][ix] == 0 || mu_ext[iz + 1][iy + 1][ix] == 0)
muyz = 0.0;
else
muyz = 1 / (0.25*(1 / mu_ext[iz][iy][ix] + 1 / mu_ext[iz][iy + 1][ix] + 1 / mu_ext[iz + 1][iy][ix] + 1 / mu_ext[iz + 1][iy + 1][ix]));
if (mu_ext[iz][iy][ix] == 0 || mu_ext[iz][iy][ix + 1] == 0 || mu_ext[iz][iy + 1][ix] == 0 || mu_ext[iz][iy + 1][ix + 1] == 0)
muxy = 0.0;
else
muxy = 1 / (0.25*(1 / mu_ext[iz][iy][ix] + 1 / mu_ext[iz][iy][ix + 1] + 1 / mu_ext[iz][iy + 1][ix] + 1 / mu_ext[iz][iy + 1][ix + 1]));
lamda2u_temp = lamda2u_ext[iz][iy][ix];
lamda_temp = lamda_ext[iz][iy][ix];
// muxz=0.0;
if ( i==0&&ix == sx&&iy == sy&&iz == sz)
{
source[iz][iy][ix] = (1 - 2 * PI*PI*F0*F0*(tt - T0)*(tt - T0))*exp(-PI*PI*F0*F0*(tt - T0)*(tt - T0));
}
else if(i==1&&(ix == (sx+1)&&iy == sy&&iz == sz))
{
source[iz][iy][ix] = -(1 - 2 * PI*PI*F0*F0*(tt - T0)*(tt - T0))*exp(-PI*PI*F0*F0*(tt - T0)*(tt - T0));
}
else if (i == 1 && (ix == (sx - 1) && iy == sy&&iz == sz))
{
source[iz][iy][ix] = (1 - 2 * PI*PI*F0*F0*(tt - T0)*(tt - T0))*exp(-PI*PI*F0*F0*(tt - T0)*(tt - T0));
}
// Qz=Q_ext[iz][ix];R_temp=R_ext[iz][ix];lamda2u_temp=lamda2u_ext[iz][ix];lamda_temp=lamda_ext[iz][ix];
/*y5 = (a1*(vx[iz][iy][ix] - vx[iz][iy][ix - 1]) +
a2*(vx[iz][iy][ix + 1] - vx[iz][iy][ix - 2])) / H;
y6 = (a1*(vy[iz][iy][ix] - vy[iz][iy - 1][ix]) +
a2*(vy[iz][iy + 1][ix] - vy[iz][iy - 2][ix])) / H;
y7 = (a1*(vz[iz][iy][ix] - vz[iz - 1][iy][ix]) +
a2*(vz[iz + 1][iy][ix] - vz[iz - 2][iy][ix])) / H;
y8 = ((a1*(vx[iz][iy + 1][ix] - vx[iz][iy][ix]) +
a2*(vx[iz][iy + 2][ix] - vx[iz][iy - 1][ix]))) / H;
y9 = ((a1*(vy[iz][iy][ix + 1] - vy[iz][iy][ix]) +
a2*(vy[iz][iy][ix + 2] - vy[iz][iy][ix - 1]))) / H;
y10 = ((a1*(vx[iz + 1][iy][ix] - vx[iz][iy][ix]) +
a2*(vx[iz + 2][iy][ix] - vx[iz - 1][iy][ix]))) / H;
y11 = ((a1*(vz[iz][iy][ix + 1] - vz[iz][iy][ix]) +
a2*(vz[iz][iy][ix + 2] - vz[iz][iy][ix - 1]))) / H;
y12 = ((a1*(vy[iz + 1][iy][ix] - vy[iz][iy][ix]) +
a2*(vy[iz + 2][iy][ix] - vy[iz - 1][iy][ix]))) / H;
y13 = ((a1*(vz[iz][iy + 1][ix] - vz[iz][iy][ix]) +
a2*(vz[iz][iy + 2][ix] - vz[iz][iy - 1][ix]))) / H;*/
y5 = (vx[iz][iy][ix] - vx[iz][iy][ix - 1]) / H;
y6 = (vy[iz][iy][ix] - vy[iz][iy - 1][ix]) / H;
y7 = (vz[iz][iy][ix] - vz[iz - 1][iy][ix]) / H;
y8 = (vx[iz][iy + 1][ix] - vx[iz][iy][ix]) / H;
y9 = (vy[iz][iy][ix + 1] - vy[iz][iy][ix]) / H;
y10 = (vx[iz + 1][iy][ix] - vx[iz][iy][ix]) / H;
y11 = (vz[iz][iy][ix + 1] - vz[iz][iy][ix]) / H;
y12 = (vy[iz + 1][iy][ix] - vy[iz][iy][ix]) / H;
y13 = (vz[iz][iy + 1][ix] - vz[iz][iy][ix]) / H;
//x = ix*H; y = iy*H; z = iz*H;
/*if (((x >= 0 && x < xoleft) && (y >= 0 && y<NY_ext*H) && (z >= 0 && z<NZ_ext*H))
|| ((x >= xoright && x < NX_ext*H) && (y >= 0 && y<NY_ext*H) && (z >= 0 && z<NZ_ext*H))
|| ((x >= 0 && x<NX_ext*H) && (y >= 0 && y < yoleft) && (z >= 0 && z<NZ_ext*H))
|| ((x >= 0 && x<NX_ext*H) && (y >= yoright && y < NY_ext*H) && (z >= 0 && z<NZ_ext*H))
|| ((x >= 0 && x<NX_ext*H) && (y >= 0 && y<NY_ext*H) && (z >= 0 && z < zoleft))
|| ((x >= 0 && x<NX_ext*H) && (y >= 0 && y<NY_ext*H) && (z >= zoright && z < NZ_ext*H)))*/
if((ix<NP||ix>=(NX_ext-1-NP))||(iy<NP||iy>=(NY_ext-1-NP))||(iz<NP||iz>=(NZ_ext-1-NP)))
{
xFHalfTemp1 = exp(-(dxi2[iz][iy][ix])*DT); xFIntTemp1 = exp(-(dxi[iz][iy][ix])*DT);
yFHalfTemp1 = exp(-(dyj2[iz][iy][ix])*DT); yFIntTemp1 = exp(-(dyj[iz][iy][ix])*DT);
zFHalfTemp1 = exp(-(dzk2[iz][iy][ix])*DT); zFIntTemp1 = exp(-(dzk[iz][iy][ix])*DT);
xFHalfTemp2 = -DT*dxi2[iz][iy][ix] * 0.5; xFIntTemp2 = -DT*dxi[iz][iy][ix] * 0.5;
yFHalfTemp2 = -DT*dyj2[iz][iy][ix] * 0.5; yFIntTemp2 = -DT*dyj[iz][iy][ix] * 0.5;
zFHalfTemp2 = -DT*dzk2[iz][iy][ix] * 0.5; zFIntTemp2 = -DT*dzk[iz][iy][ix] * 0.5;
pmlxVx[iz][iy][ix] = xFIntTemp1*pmlxVx[iz][iy][ix] + xFIntTemp2*(xFIntTemp1*DxVxPre[iz][iy][ix] + y5);
pmlyVy[iz][iy][ix] = yFIntTemp1*pmlyVy[iz][iy][ix] + yFIntTemp2*(yFIntTemp1*DyVyPre[iz][iy][ix] + y6);
pmlzVz[iz][iy][ix] = zFIntTemp1*pmlzVz[iz][iy][ix] + zFIntTemp2*(zFIntTemp1*DzVzPre[iz][iy][ix] + y7);
DxVxPre[iz][iy][ix] = y5; DyVyPre[iz][iy][ix] = y6; DzVzPre[iz][iy][ix] = y7;
y5 = y5 + pmlxVx[iz][iy][ix];
y6 = y6 + pmlyVy[iz][iy][ix];
y7 = y7 + pmlzVz[iz][iy][ix];
pmlzVx[iz][iy][ix] = zFHalfTemp1*pmlzVx[iz][iy][ix] + zFHalfTemp2*(zFHalfTemp1*DzVxPre[iz][iy][ix] + y10);
pmlxVz[iz][iy][ix] = xFHalfTemp1*pmlxVz[iz][iy][ix] + xFHalfTemp2*(xFHalfTemp1*DxVzPre[iz][iy][ix] + y11);
DzVxPre[iz][iy][ix] = y10; DxVzPre[iz][iy][ix] = y11;
y10 = y10 + pmlzVx[iz][iy][ix];
y11 = y11 + pmlxVz[iz][iy][ix];
pmlyVx[iz][iy][ix] = yFHalfTemp1*pmlyVx[iz][iy][ix] + yFHalfTemp2*(yFHalfTemp1*DyVxPre[iz][iy][ix] + y8);
pmlxVy[iz][iy][ix] = xFHalfTemp1*pmlxVy[iz][iy][ix] + xFHalfTemp2*(xFHalfTemp1*DxVyPre[iz][iy][ix] + y9);
DyVxPre[iz][iy][ix] = y8; DxVyPre[iz][iy][ix] = y9;
y8 = y8 + pmlyVx[iz][iy][ix];
y9 = y9 + pmlxVy[iz][iy][ix];
pmlzVy[iz][iy][ix] = zFHalfTemp1*pmlzVy[iz][iy][ix] + zFHalfTemp2*(zFHalfTemp1*DzVyPre[iz][iy][ix] + y12);
pmlyVz[iz][iy][ix] = yFHalfTemp1*pmlyVz[iz][iy][ix] + yFHalfTemp2*(yFHalfTemp1*DyVzPre[iz][iy][ix] + y13);
DzVyPre[iz][iy][ix] = y12; DyVzPre[iz][iy][ix] = y13;
y12 = y12 + pmlzVy[iz][iy][ix];
y13 = y13 + pmlyVz[iz][iy][ix];
}
txx[iz][iy][ix] = txx[iz][iy][ix] + (lamda2u_temp*DT)*y5 + (lamda_temp*DT)*y6 + (lamda_temp*DT)*y7 + source[iz][iy][ix];
tyy[iz][iy][ix] = tyy[iz][iy][ix] + (lamda_temp*DT)*y5 + (lamda2u_temp*DT)*y6 + (lamda_temp*DT)*y7 + source[iz][iy][ix];
tzz[iz][iy][ix] = tzz[iz][iy][ix] + (lamda_temp*DT)*y5 + (lamda_temp*DT)*y6 + (lamda2u_temp*DT)*y7 + source[iz][iy][ix];
txz[iz][iy][ix] = txz[iz][iy][ix] + (muxz*DT)*(y10 + y11);
txy[iz][iy][ix] = txy[iz][iy][ix] + (muxy*DT)*(y8 + y9);
tyz[iz][iy][ix] = tyz[iz][iy][ix] + (muyz*DT)*(y12 + y13);
//加震源—单极子
if (i == 0)
{
tzz[sz][sy][sx] =(1 -2 * PI*PI*F0*F0*(tt - T0)*(tt - T0))*exp(-PI*PI*F0*F0*(tt - T0)*(tt - T0)); //对正应力加爆炸源
tyy[sz][sy][sx] = (1-2 * PI*PI*F0*F0*(tt - T0)*(tt - T0))*exp(-PI*PI*F0*F0*(tt - T0)*(tt - T0));
txx[sz][sy][sx] =(1 -2 * PI*PI*F0*F0*(tt - T0)*(tt - T0))*exp(-PI*PI*F0*F0*(tt - T0)*(tt - T0));
}
//加震源—偶极子
if (i == 1)
{
tzz[sz][sy][sx + 1] = tzz[sz][sy][sx + 1] - (1 - 2 * PI*PI*F0*F0*(tt - T0)*(tt - T0))*exp(-PI*PI*F0*F0*(tt - T0)*(tt - T0));; //对正应力加爆炸源
tyy[sz][sy][sx + 1] = tyy[sz][sy][sx + 1] - (1 - 2 * PI*PI*F0*F0*(tt - T0)*(tt - T0))*exp(-PI*PI*F0*F0*(tt - T0)*(tt - T0));;
txx[sz][sy][sx + 1] = txx[sz][sy][sx + 1] - (1 - 2 * PI*PI*F0*F0*(tt - T0)*(tt - T0))*exp(-PI*PI*F0*F0*(tt - T0)*(tt - T0));;
tzz[sz][sy][sx - 1] = tzz[sz][sy][sx - 1] + (1 - 2 * PI*PI*F0*F0*(tt - T0)*(tt - T0))*exp(-PI*PI*F0*F0*(tt - T0)*(tt - T0));; //对正应力加爆炸源
tyy[sz][sy][sx - 1] = tyy[sz][sy][sx - 1] + (1 - 2 * PI*PI*F0*F0*(tt - T0)*(tt - T0))*exp(-PI*PI*F0*F0*(tt - T0)*(tt - T0));;
txx[sz][sy][sx - 1] = txx[sz][sy][sx - 1] + (1 - 2 * PI*PI*F0*F0*(tt - T0)*(tt - T0))*exp(-PI*PI*F0*F0*(tt - T0)*(tt - T0));;
}
}
#pragma omp barrier
#pragma omp master
{
// tzz_z[sz][sx]=tzz_z[sz][sx]-2*PI*PI*F0*F0*(tt-T0)*exp(-PI*PI*F0*F0*(tt-T0)*(tt-T0));
//sis_y[1][it] = tzz[sz][sy][sx];
//记录地震记录
if (i == 0)
{
for (iz = NP; iz < NP + NZ; iz++)
{
sis_z[iz - NP][it] = tzz[iz][sy][sx];
sis_x[iz - NP][it] = txx[iz][sy][sx];
//sis_z[iz-NP][it]=vz_x[iz][NX_ext-iz]+vz_z[iz][NX_ext-iz];
// sis_x[iz-NP][it]=vx_x[iz][NX_ext-iz]+vx_z[iz][NX_ext-iz];
}
}
if (i == 1)
{
for (iz = NP; iz < NP + NZ; iz++)
{
sis_z[iz - NP][it] = tzz[iz][sy][sx];
sis_x[iz - NP][it] = txx[iz][sy][sx+1];
//sis_z[iz-NP][it]=vz_x[iz][NX_ext-iz]+vz_z[iz][NX_ext-iz];
// sis_x[iz-NP][it]=vx_x[iz][NX_ext-iz]+vx_z[iz][NX_ext-iz];
}
}
if (it % 200 == 0)
{
int n = it / 200;
for (iz = 2; iz < NZ_ext - 2; iz = iz + 2)
{
int a = iz / 2;
for (ix = 2; ix < NX_ext - 2; ix = ix + 2)
{
int b = ix / 2;
vxx[n][a][b] = txx[iz][sy][ix];
}
}
}
}
#pragma omp barrier
}
}
finish = clock();
printf("%f seconds\n", (double)(finish - start) / CLOCKS_PER_SEC);
//输出波场快照
char vzname[] = "vz.dat";
char vxname[] = "vx.dat";
char vxxname[] = "vxx.dat";
char vyname[] = "vy.dat";
char txx_xname[] = "txx_x.dat";
char pmlxSxxname[] = "pmlxSxx.dat";
char DxSxxPrename[] = "DxSxxPre.dat";
/* for (iz=0; iz<NZ_ext; iz++)
{
for (iy=0; iy<NY_ext; iy++)
for (ix=0; ix<NX_ext; ix++)
{
vx[iz][iy][ix]=vx_x[iz][iy][ix]+vx_y[iz][iy][ix]+vx_z[iz][iy][ix];
vy[iz][iy][ix]=vy_x[iz][iy][ix]+vy_y[iz][iy][ix]+vy_z[iz][iy][ix];
vz[iz][iy][ix]=vz_x[iz][iy][ix]+vz_y[iz][iy][ix]+vz_z[iz][iy][ix];
}
}*/
printf("NZ_ext = %d NY_ext = %d NX_ext = %d\n", NZ_ext, NY_ext, NX_ext);
wfile3d(vxname, vx, NZ_ext, NY_ext, NX_ext);
wfile4d(vxxname, vxx, NT/200, NZ_ext/2, NX_ext/2);
wfile3d(vzname, vz, NZ_ext, NY_ext, NX_ext);
wfile3d(vyname, vy, NZ_ext, NY_ext, NX_ext);
wfile3d(pmlxSxxname, pmlxSxx, NZ_ext, NY_ext, NX_ext);
wfile3d(DxSxxPrename, DxSxxPre, NZ_ext, NY_ext, NX_ext);
//wfile(ssxname, ssx, NZ_ext, NX_ext);
//输出地震记录
char sisxname[] = "sisx.dat";
char sisyname[] = "sisy.dat";
char siszname[] = "sisz.dat";
printf("NZ = %d NT = %d\n", NZ, NT);
wfile(sisxname, sis_x, NZ, NT);
wfile(sisyname, sis_y, 1, NT);
wfile(siszname, sis_z, NZ, NT);
free_space3d(tzz, NZ_ext, NY_ext);
// free_space3d(tzz_y,NZ_ext,NY_ext);
// free_space3d(tzz_z,NZ_ext,NY_ext);
free_space3d(tyy, NZ_ext, NY_ext);
// free_space3d(tyy_y,NZ_ext,NY_ext);
// free_space3d(tyy_z,NZ_ext,NY_ext);
free_space3d(txx, NZ_ext, NY_ext);
// free_space3d(txx_y,NZ_ext,NY_ext);
// free_space3d(txx_z,NZ_ext,NY_ext);
free_space3d(txy, NZ_ext, NY_ext);
// free_space3d(txy_y,NZ_ext,NY_ext);
free_space3d(txz, NZ_ext, NY_ext);
// free_space3d(txz_z,NZ_ext,NY_ext);
free_space3d(tyz, NZ_ext, NY_ext);
// free_space3d(tyz_z,NZ_ext,NY_ext);
free_space3d(vx, NZ_ext, NY_ext);
// free_space3d(vx_y,NZ_ext,NY_ext);
// free_space3d(vx_z,NZ_ext,NY_ext);
free_space3d(vy, NZ_ext, NY_ext);
// free_space3d(vy_y,NZ_ext,NY_ext);
// free_space3d(vy_z,NZ_ext,NY_ext);
free_space3d(vz, NZ_ext, NY_ext);
// free_space3d(vz_y,NZ_ext,NY_ext);
// free_space3d(vz_z,NZ_ext,NY_ext);
free_space3d(vs_ext, NZ_ext, NY_ext);
free_space3d(vp_ext, NZ_ext, NY_ext);
free_space3d(rho_ext, NZ_ext, NY_ext);
free_space3d(vs, NZ, NY);
free_space3d(vp, NZ, NY);
free_space3d(mu, NZ, NY);
free_space3d(mu_ext, NZ_ext, NY_ext);
free_space3d(lamda, NZ, NY);
free_space3d(lamda2u, NZ, NY);
free_space3d(lamda_ext, NZ_ext, NY_ext);
free_space3d(lamda2u_ext, NZ_ext, NY_ext);
free_space3d(rho, NZ, NY);
free_space3d(dxi, NZ_ext, NY_ext);
free_space3d(dxi2, NZ_ext, NY_ext);
free_space3d(dyj, NZ_ext, NY_ext);
free_space3d(dyj2, NZ_ext, NY_ext);
free_space3d(dzk, NZ_ext, NY_ext);
free_space3d(dzk2, NZ_ext, NY_ext);
free_space3d(pmlxSxx, NZ_ext, NY_ext);
free_space3d(pmlySxy, NZ_ext, NY_ext);
free_space3d(pmlzSxz, NZ_ext, NY_ext);
free_space3d(DxSxxPre, NZ_ext, NY_ext);
free_space3d(DySxyPre, NZ_ext, NY_ext);
free_space3d(DzSxzPre, NZ_ext, NY_ext);
free_space3d(pmlxSxy, NZ_ext, NY_ext);
free_space3d(pmlySyy, NZ_ext, NY_ext);
free_space3d(pmlzSyz, NZ_ext, NY_ext);
free_space3d(DxSxyPre, NZ_ext, NY_ext);
free_space3d(DySyyPre, NZ_ext, NY_ext);
free_space3d(DzSyzPre, NZ_ext, NY_ext);
free_space3d(pmlxSxz, NZ_ext, NY_ext);
free_space3d(pmlySyz, NZ_ext, NY_ext);
free_space3d(pmlzSzz, NZ_ext, NY_ext);
free_space3d(DxSxzPre, NZ_ext, NY_ext);
free_space3d(DySyzPre, NZ_ext, NY_ext);
free_space3d(DzSzzPre, NZ_ext, NY_ext);
free_space3d(pmlxVx, NZ_ext, NY_ext);
free_space3d(pmlyVy, NZ_ext, NY_ext);
free_space3d(pmlzVz, NZ_ext, NY_ext);
free_space3d(DxVxPre, NZ_ext, NY_ext);
free_space3d(DyVyPre, NZ_ext, NY_ext);
free_space3d(DzVzPre, NZ_ext, NY_ext);
free_space3d(pmlxVz, NZ_ext, NY_ext);
free_space3d(pmlzVx, NZ_ext, NY_ext);
free_space3d(DxVzPre, NZ_ext, NY_ext);
free_space3d(DzVxPre, NZ_ext, NY_ext);
free_space3d(pmlyVz, NZ_ext, NY_ext);
free_space3d(pmlzVy, NZ_ext, NY_ext);
free_space3d(DyVzPre, NZ_ext, NY_ext);
free_space3d(DzVyPre, NZ_ext, NY_ext);
free_space3d(pmlxVy, NZ_ext, NY_ext);
free_space3d(pmlyVx, NZ_ext, NY_ext);
free_space3d(DxVyPre, NZ_ext, NY_ext);
free_space3d(DyVxPre, NZ_ext, NY_ext);
free_space3d(source, NZ_ext, NY_ext);
free_space3d(vxx, NT/200, NZ_ext/2);
free_space2d(sis_x, NZ);
free_space2d(sis_z, NZ);
free_space2d(sis_y, NZ);
printf("\\Press any key to exit program...");
//getch();
return 0;
}
//申请二维动态数组
float **space2d(int nr, int nc)
{
float **a;
int i;
a = (float **)calloc(nr, sizeof(float *));
for (i = 0; i<nr; i++)
{
a[i] = (float *)calloc(nc, sizeof(float));
}
return a;
}
//申请三维动态数组
/*float **space3d(int nr,int ny,int nc)
{
float **a;
int i,j;
a=(float ***)calloc(nr,sizeof(float **));
for(i=0;i<nr;i++)
{
a[i]=(float **)calloc(ny,sizeof(float ));
}
for(i=0;i<nr;i++)
{
for(j=0;j<ny;j++)
{
a[i][j]=(float *)calloc(nc,sizeof(float));//sizeof(float) nc改为sizeof(float)*nc
}}
return a;
}/
/*float ***space3d(int nr,int ny,int nc)
{
float ***a;
a = new float **[nr];
for(int i=0;i<nr;i++)
{
a[i] = new float *[ny];
for(int j=0;j<ny;j++)
{
a[i][j]=new float [nc];
}}
return a;
}*/
float **space3d(int nr, int ny, int nc)
{
float **a;
int i, j, k;
a = (float ***)malloc(sizeof(float **)*nr);
for (i = 0; i<nr; i++)
{
a[i] = (float **)malloc(sizeof(float *)*ny);
}
for (i = 0; i<nr; i++)
{
for (j = 0; j<ny; j++)
{
a[i][j] = (float *)malloc(sizeof(float)*nc);//sizeof(float) nc改为sizeof(float)*nc
}
}
for (i = 0; i<nr; i++)
for (j = 0; j<ny; j++)
for (k = 0; k<nc; k++)
{
a[i][j][k] = 0.0f;
}
return a;
}
//释放二维动态数组
void free_space2d(float **a, int nr)
{
int i;
for (i = 0; i<nr; i++)
free(a[i]);
free(a);
}
//释放三维动态数组
/*void free_space3d(float **a,int nr,int ny)
{
for(int i=0;i<nr;i++)
{
for(int j=0;j<ny;j++)
{
delete [] a[i][j];
}
delete [] a[i];
}
delete [] a;
a=NULL;
}
*/
void free_space3d(float **a, int nr, int ny)
{
int i, j;
for (i = 0; i<nr; i++)
{
for (j = 0; j<ny; j++)
{
free(a[i][j]);
}
}
for (i = 0; i<nr; i++)
{
free(a[i]);
}
free(a);
}
//将二进制数据写入文件—二维
void wfile(char filename[], float **data, int nr, int nc)
{
int i, j;
FILE fp = fopen(filename, "wt");
/
for(int i=0;i<nr;i++)
{
fwrite(data[i],sizeof(float),nc,fp);
}
*/
for (i = 0; i<nr; i++)
{
for (j = 0; j<nc; j++)
{
fprintf(fp, "%e ", data[i][j]);
if ((j + 1) % nc == 0)
fprintf(fp, "\n");
}
fprintf(fp, "\n");
}
// fwrite(&data[i][j],1,sizeof(float),fp);
fclose(fp);
}
//将二进制数据写入文件—三维
void wfile3d(char filename[], float ***data, int nr, int ny, int nc)
{
int i, j, k;
FILE fp = fopen(filename, "wt");
/
for(int i=0;i<nr;i++)
{
fwrite(data[i],sizeof(float),nc,fp);
}
*/
{
for (i = 0; i<nr; i++)
{
for (k = 0; k<nc; k++)
{
j =70;
fprintf(fp, "%e ", data[i][j][k]);
if ((k + 1) % nc == 0)
fprintf(fp, "\n");
}
fprintf(fp, "\n");
}
}
// fwrite(&data[i][j],1,sizeof(float),fp);
fclose(fp);
}
void wfile4d(char filename[], float ***data, int nr, int ny, int nc)
{
int i, j, k;
FILE fp = fopen(filename, "wt");
/
for(int i=0;i<nr;i++)
{
fwrite(data[i],sizeof(float),nc,fp);
}
*/
{
for (i = 0; i < nr; i ++)
{
for (j = 0; j < ny; j++)
{
for (k = 0; k < nc; k++)
{
fprintf(fp, "%e ", data[i][j][k]);
if ((k + 1) % nc == 0)
fprintf(fp, "\n");
}
}
fprintf(fp, "\n");
}
}
// fwrite(&data[i][j],1,sizeof(float),fp);
fclose(fp);
}
//建立地质模型
void create_model(float **vp, float **vs, float **rho, float **vf, float **rhof, float **lamda, float **lamda2u, float **mu, float ***por, int nr, int ny, int nc)
{
//这里进行修改，最好是可视化设计或从文件读入
int ix, iy, iz;
// double por1,Ks1,Kb1,Kf1,a1,D1,tao,rou11,rou12,rou22,eta;
// eta=58.0;
for (iz = 0; iz < nr; iz++)
{
for (iy = 0; iy < ny; iy++)
{
for (ix = 0; ix < nc; ix++)
{
if (((ix - nc / 2)*(ix - nc / 2) + (iy - ny / 2)*(iy - ny / 2)) <= 100)//竖直井孔%%%%%%\
// if(iz<(nr-14-ix)||iz>=(nr+14-ix))//45度倾斜井孔
//if (iz>= 50)
//if(ix>=50)
{
vp[iz][iy][ix] = 1500.0f;
vs[iz][iy][ix] = 0.0f;
vf[iz][iy][ix] = 1500.0f;
rho[iz][iy][ix] = 1000.0f;
rhof[iz][iy][ix] = 1000.0f;
por[iz][iy][ix] = 1.0f;
// b[iz][ix]=0.0;
mu[iz][iy][ix] = rho[iz][iy][ix] * vs[iz][iy][ix] * vs[iz][iy][ix];
lamda2u[iz][iy][ix] = rho[iz][iy][ix] * vp[iz][iy][ix] * vp[iz][iy][ix];
lamda[iz][iy][ix] = lamda2u[iz][iy][ix] - 2 * mu[iz][iy][ix];
/* Kb1=rho[iz][ix]*(1-por[iz][ix])*(vp[iz][ix]*vp[iz][ix]-vs[iz][ix]*vs[iz][ix]*4.0/3.0);//骨架压缩模量，声波测井原理与应用，P39
Ks1=rho[iz][ix]*(vp[iz][ix]*vp[iz][ix]-vs[iz][ix]*vs[iz][ix]*4.0/3.0);//岩石固态颗粒的体积模量
Kf1=rhof[iz][ix]*vf[iz][ix]*vf[iz][ix];//孔隙流体的体积压缩模量
a1=1-Kb1/Ks1;
D1=a1-por[iz][ix]+por[iz][ix]*Ks1/Kf1;
Q[iz][ix]=(a1-por[iz][ix])*por[iz][ix]*Ks1/D1;
R[iz][ix]=por[iz][ix]*por[iz][ix]*Ks1/D1;
tao=0.5*(1+1/por[iz][ix]);//孔隙弯曲度
rou11=(1-por[iz][ix])*rho[iz][ix]-(1-tao)*por[iz][ix]*rhof[iz][ix];
rou12=(1-tao)*por[iz][ix]*rhof[iz][ix];
rou22=tao*por[iz][ix]*rhof[iz][ix];
D11[iz][ix]=rou11/(rou11*rou22-rou12*rou12);
D12[iz][ix]=rou12/(rou11*rou22-rou12*rou12);
D22[iz][ix]=rou22/(rou11*rou22-rou12*rou12);*/
}
else
{
vp[iz][iy][ix] = 4000.0f; //纵波
vs[iz][iy][ix] = 2340.0f; //横波
rho[iz][iy][ix] = 2500.0f; //密度
vf[iz][iy][ix] = 1500.0f;
rhof[iz][iy][ix] = 1000.0f;
por[iz][iy][ix] = 0.2f;
//b[iz][ix]=0.2*pow(10,-3);// b=eta*por*por/perm*pow(10,-3);
mu[iz][iy][ix] = rho[iz][iy][ix] * vs[iz][iy][ix] * vs[iz][iy][ix];
lamda2u[iz][iy][ix] = rho[iz][iy][ix] * vp[iz][iy][ix] * vp[iz][iy][ix];
lamda[iz][iy][ix] = lamda2u[iz][iy][ix] - 2.0*mu[iz][iy][ix];
}
}
}
}
}
/*
//将模型扩边,用于PML
//具体的操作过程是将实际模型参数放置在扩边后的数据中央，四周的数据用
//最外缘的数据填充
float **extmodel(float **inqit_model,int nz,int nx,int np)
{
float **p;
int i,j;
int nx2=nx+2*np;
int nz2=nz+2*np;
p=space2d(nz2,nx2);
for (i=np; i<nz+np; i++)
for (j=0; j<np; j++)
p[i][j]=init_model[i-np][0];
for (i=np; i<nz+np; i++)
for (j=nx+np; j<nx2; j++)
p[i][j]=init_model[i-np][nx-1];
for (i=nz; i<nz2; i++)
for (j=np; j<np+nx; j++)
p[i][j]=init_model[nz-1][j-np];
for(i=0; i<np; i++)
for(j=np; j<np+nx; j++)
p[i][j]=init_model[0][j-np];
for(i=0; i<np; i++)
for(j=0; j<np; j++)
p[i][j]=init_model[0][0];
for(i=0; i<np; i++)
for(j=nx+np; j<nx2; j++)
p[i][j]=init_model[0][nx-1];
for (i=nz+np; i<nz2; i++)
for (j=0; j<np; j++)
p[i][j]=init_model[nz-1][0];
for (i=nz+np; i<nz2; i++)
for (j=nx+np; j<nx2; j++)
p[i][j]=init_model[nz-1][nx-1];
for (i=np; i<nz+np; i++)
for (j=np; j<nx+np; j++)
p[i][j]=init_model[i-np][j-np];
return p;
}*/
//将模型扩边,用于PML
//具体的操作过程是将实际模型参数放置在扩边后的数据中央，四周的数据用
//最外缘的数据填充_3D
float **extmodel(float **init_model, int nz, int ny, int nx, int np)
{
float ***p;
int i, j, k;
int nx2 = nx + 2 * np;
int ny2 = ny + 2 * np;
int nz2 = nz + 2 * np;
p = space3d(nz2, ny2, nx2);
for (i = 0; i<np; i++)
for (k = 0; k<np; k++)
for (j = 0; j<np; j++)
p[i][k][j] = init_model[0][0][0];
for (i = 0; i<np; i++)
for (k = np; k<np + ny; k++)
for (j = 0; j<np; j++)
p[i][k][j] = init_model[0][k - np][0];
for (i = 0; i<np; i++)
for (k = np + ny; k<ny2; k++)
for (j = 0; j<np; j++)
p[i][k][j] = init_model[0][ny - 1][0];
for (i = np; i<nz + np; i++)
for (k = 0; k<np; k++)
for (j = 0; j<np; j++)
p[i][k][j] = init_model[i - np][0][0];
for (i = np; i<nz + np; i++)
for (k = np; k<np + ny; k++)
for (j = 0; j<np; j++)
p[i][k][j] = init_model[i - np][k - np][0];
for (i = np; i<nz + np; i++)
for (k = ny + np; k<ny2; k++)
for (j = 0; j<np; j++)
p[i][k][j] = init_model[i - np][ny - 1][0];
for (i = nz + np; i<nz2; i++)
for (k = 0; k<np; k++)
for (j = 0; j<np; j++)
p[i][k][j] = init_model[nz - 1][0][0];
for (i = nz + np; i<nz2; i++)
for (k = np; k<np + ny; k++)
for (j = 0; j<np; j++)
p[i][k][j] = init_model[nz - 1][k - np][0];
for (i = nz + np; i<nz2; i++)
for (k = ny + np; k<ny2; k++)
for (j = 0; j<np; j++)
p[i][k][j] = init_model[nz - 1][ny - 1][0];
for (i = 0; i<np; i++)
for (k = 0; k<np; k++)
for (j = np; j<np + nx; j++)
p[i][k][j] = init_model[0][0][j - np];
for (i = 0; i<np; i++)
for (k = np; k<np + ny; k++)
for (j = np; j<np + nx; j++)
p[i][k][j] = init_model[0][k - np][j - np];
for (i = 0; i<np; i++)
for (k = ny + np; k<ny2; k++)
for (j = np; j<np + nx; j++)
p[i][k][j] = init_model[0][ny - 1][j - np];
for (i = np; i<nz + np; i++)
for (k = 0; k<np; k++)
for (j = np; j<np + nx; j++)
p[i][k][j] = init_model[i - np][0][j - np];
for (i = np; i<nz + np; i++)
for (k = np; k<np + ny; k++)
for (j = np; j<np + nx; j++)
p[i][k][j] = init_model[i - np][k - np][j - np];
for (i = np; i<nz + np; i++)
for (k = ny + np; k<ny2; k++)
for (j = np; j<np + nx; j++)
p[i][k][j] = init_model[i - np][ny - 1][j - np];
for (i = nz + np; i<nz2; i++)
for (k = 0; k<np; k++)
for (j = np; j<np + nx; j++)
p[i][k][j] = init_model[nz - 1][0][j - np];
for (i = nz + np; i<nz2; i++)
for (k = np; k<np + ny; k++)
for (j = np; j<np + nx; j++)
p[i][k][j] = init_model[nz - 1][k - np][j - np];
for (i = nz + np; i<nz2; i++)
for (k = ny + np; k<ny2; k++)
for (j = np; j<np + nx; j++)
p[i][k][j] = init_model[nz - 1][ny - 1][j - np];
for (i = 0; i<np; i++)
for (k = 0; k<np; k++)
for (j = np + nx; j<nx2; j++)
p[i][k][j] = init_model[0][0][nx - 1];
for (i = 0; i<np; i++)
for (k = np; k<np + ny; k++)
for (j = np + nx; j<nx2; j++)
p[i][k][j] = init_model[0][k - np][nx - 1];
for (i = 0; i<np; i++)
for (k = ny + np; k<ny2; k++)
for (j = np + nx; j<nx2; j++)
p[i][k][j] = init_model[0][ny - 1][nx - 1];
for (i = np; i<nz + np; i++)
for (k = 0; k<np; k++)
for (j = np + nx; j<nx2; j++)
p[i][k][j] = init_model[i - np][0][nx - 1];
for (i = np; i<nz + np; i++)
for (k = np; k<np + ny; k++)
for (j = np + nx; j<nx2; j++)
p[i][k][j] = init_model[i - np][k - np][nx - 1];
for (i = np; i<nz + np; i++)
for (k = ny + np; k<ny2; k++)
for (j = np + nx; j<nx2; j++)
p[i][k][j] = init_model[i - np][ny - 1][nx - 1];
for (i = nz + np; i<nz2; i++)
for (k = 0; k<np; k++)
for (j = np + nx; j<nx2; j++)
p[i][k][j] = init_model[nz - 1][0][nx - 1];
for (i = nz + np; i<nz2; i++)
for (k = np; k<np + ny; k++)
for (j = np + nx; j<nx2; j++)
p[i][k][j] = init_model[nz - 1][k - np][nx - 1];
for (i = nz + np; i<nz2; i++)
for (k = ny + np; k<ny2; k++)
for (j = np + nx; j<nx2; j++)
p[i][k][j] = init_model[nz - 1][ny - 1][nx - 1];
return p;
}
