qq_32787919
qq_32787919
采纳率100%
2018-09-21 07:01 阅读 911

c++openmp问题 程序正常运行可以,放在openmp中不运行!!急急急

#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;

}

  • 点赞
  • 写回答
  • 关注问题
  • 收藏
  • 复制链接分享

2条回答 默认 最新

  • 已采纳
    shen_wei shen_wei 2018-09-21 07:26

    此类问题要么单步自己调试,要么做成可下载工程。。

    问答已经把代码给完全破坏了,代码也太长了。。。码农都习惯编译 --调试---解决BUG!!!!

    点赞 评论 复制链接分享
  • caozhy 从今以后生命中的每一秒都属于我爱的人 2018-09-22 02:56

    那么你凭什么认为你把代码一丢,就可以等到回答呢?

    点赞 评论 复制链接分享

相关推荐