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个回答

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

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

qq_32787919
qq_32787919 您好,我这个没有办法,就是说上面的程序可以正常运行,但是加在openmp就有问题,所以我把所有的程序贴上来,希望懂这个的,帮我看看,这个基本的openmp语句都没有问题。
12 个月之前 回复

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

qq_32787919
qq_32787919 全放上来,才能说明问题。
12 个月之前 回复
Csdn user default icon
上传中...
上传图片
插入图片
抄袭、复制答案,以达到刷声望分或其他目的的行为,在CSDN问答是严格禁止的,一经发现立刻封号。是时候展现真正的技术了!
其他相关推荐
c++openmp问题 程序正常运行可以,放在openmp中不运行!!急急急
#pragma hdrstoprnrn//---------------------------------------------------------------------------rnrn//#pragma argsusedrn#include rn#include rn#include rn#include rn#include rn#includern//#include rn#define PI 3.1415926535rn#define e 2.718281828rnrnrnrnfloat **space2d(int nr, int nc);rnfloat ***space3d(int nr, int ny, int nc);rnvoid free_space2d(float **a, int nr);rnvoid free_space3d(float ***b, int nr, int ny);rnvoid wfile(char filename[], float **data, int nr, int nc);rnvoid wfile3d(char filename[], float ***data, int nr, int ny, int nc);rnvoid wfile4d(char filename[], float ***data, int nr, int ny, int nc);rnvoid 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);rnrnfloat ***extmodel(float ***init_model, int nr, int ny, int nc, int np);rnrnrnint main()rnrn //给定参数rnrn int NX = 100; //x方向网格点数%%%%%rn int NY = 100; //y方向网格点数%%%%%rn int NZ = 300; //z方向网格点数%%%%rn int NP = 20; //pml层网格点数rn int i = 0;rnrn int sx = NX / 2 + NP; //震源坐标点号rn int sy = NY / 2 + NP;rn int sz = NZ /6 + NP;rnrn int NX_ext = NX + 2 * NP; //rn int NY_ext = NY + 2 * NP;rn int NZ_ext = NZ + 2 * NP;rnrn int NT = 8000; //空间步长5%%%rn float H = 0.01f;rn float RC = 0.000001f;rn float DP = NP*H;rnrn // double DT=0.0002; //时间步长rn double DT = 1.0*pow(10, -6); //时间步长%%%%rn // double F0=30.0; //震源主频rn double F0 = 3.0*pow(10, 3); //震源主频,主频太高会造成频散rn //double T0=1.2/F0; rn float T0 = 1.0 / F0;rn double Vpmax = 4270.0; //模型最大纵波速度,用于稳定性计算rn double Vpmin = 1500.0; //模型最小纵波速度,用于控制数值频散rn double Vsmax = 2650.0;rn double Vsmin = 1500.0;rnrn float ***vs; //初始模型横波速度rn float ***vp; //初始模型纵波速度rn float ***rho; //初始模型密度rn float ***vf; //孔隙流体速度rn float ***rhof; //孔隙流体密度rn float ***mu; //剪切模量rn float ***lamda; //拉梅系数rn float ***lamda2u;rn //float **Q;rn //float **R;rn float ***por;rn // float **D11,**D12,**D22;rn // float **b;rn float **sis_x; //地震记录x分量 波形记录-二维数据rn float **sis_y; //地震记录y分量rn float **sis_z; //地震记录z分量rnrnrnrn float ***vx;rn float ***vxx;rnrn float ***vy;rnrn float ***vz;rn /*float **Vx_x;//流相速度rn float **Vx_z;rn float **Vx;rn float **Vz_x;rn float **Vz_z;rn float **Vz;rn */rnrn float ***txx;//应力分量rnrn float ***tyy;//应力分量rnrnrn float ***tzz;rnrn float ***txy;rn float ***txz;rn float ***tyz;rn float ***source;rn /*float **ssx;//流体应力rn float **ssz;rn float **dxi,**dxi2;rn float **dzj,**dzj2;rn */rn // double Kb,Ks,Kf,a,D;rn float tt, x, y, z, xoleft, xoright, yoleft, yoright, zoleft, zoright, d0=0.0f;rn float v0, y5, y6 , y7 , y8, y9 , y10, y11, y12, y13, rho_tempx, rho_tempy, rho_tempz, muxz, muxy, muyz, lamda2u_temp, lamda_temp;rn // double D11_tempx,D12_tempx,D22_tempx,R_temp,D11_tempz,D12_tempz,D22_tempz;rn float xx1, xy1, xz1, yx2 , yy2, yz2 , zx3, zy3, zz3 , best_dt;rn int ix, iy, iz, it;rn rnrn vs = space3d(NZ, NY, NX);rn vp = space3d(NZ, NY, NX);rn rho = space3d(NZ, NY, NX);rn vf = space3d(NZ, NY, NX);rn rhof = space3d(NZ, NY, NX);rn //Q=space2d(NZ,NX);rn //R=space2d(NZ,NX);rn por = space3d(NZ, NY, NX);rn //D11=space2d(NZ,NX);D12=space2d(NZ,NX);D22=space2d(NZ,NX);rn // b=space2d(NZ,NX);rn mu = space3d(NZ, NY, NX);rn lamda = space3d(NZ, NY, NX);rn lamda2u = space3d(NZ, NY, NX);rn create_model(vp, vs, rho, vf, rhof, lamda, lamda2u, mu, por, NZ, NY, NX);rnrn char vp_name[] = "vp_ext.dat";rn float ***vs_ext;rn float ***vp_ext;rn float ***rho_ext;rn float ***mu_ext;rn float ***lamda_ext;rn float ***lamda2u_ext;rn // float ***Q_ext;rn // float ***R_ext;rn float ***vf_ext;rn float ***rhof_ext;rn float ***por_ext;rn float ***dxi;rn float ***dxi2;rn float ***dyj;rn float ***dyj2;rn float ***dzk;rn float ***dzk2;rn //float ***b_ext;rn // float **D11_ext,**D12_ext,**D22_ext;rnrnrn //数据扩充,加上吸收边界rn vs_ext = extmodel(vs, NZ, NY, NX, NP);rn vp_ext = extmodel(vp, NZ, NY, NX, NP);rn rho_ext = extmodel(rho, NZ, NY, NX, NP);rn mu_ext = extmodel(mu, NZ, NY, NX, NP);rn lamda_ext = extmodel(lamda, NZ, NY, NX, NP);rn lamda2u_ext = extmodel(lamda2u, NZ, NY, NX, NP);rn vf_ext = extmodel(vf, NZ, NY, NX, NP);rn rhof_ext = extmodel(rhof, NZ, NY, NX, NP);rn //Q_ext=extmodel(Q,NZ,NX,NP);rn //R_ext=extmodel(R,NZ,NX,NP);rn por_ext = extmodel(por, NZ, NY, NX, NP);rn // b_ext=extmodel(b,NZ,NX,NP);rn // D11_ext=extmodel(D11,NZ,NX,NP);D12_ext=extmodel(D12,NZ,NX,NP);rn // D22_ext=extmodel(D22,NZ,NX,NP);rnrn wfile3d(vp_name, vp_ext, NZ_ext, NY_ext, NX_ext); //建立文件rnrn //申请空间rnrn vx = space3d(NZ_ext, NY_ext, NX_ext);rnrnrn vy = space3d(NZ_ext, NY_ext, NX_ext);rn vxx= space3d(NT/200,NZ_ext/2, NX_ext/2);rnrn vz = space3d(NZ_ext, NY_ext, NX_ext);rn /*Vx_x=space2d(NZ_ext,NX_ext);rn Vx_z=space2d(NZ_ext,NX_ext);rn Vx=space2d(NZ_ext,NX_ext);rn Vz_x=space2d(NZ_ext,NX_ext);rn Vz_z=space2d(NZ_ext,NX_ext);rn Vz=space2d(NZ_ext,NX_ext);*/rn txx = space3d(NZ_ext, NY_ext, NX_ext);rnrnrn tyy = space3d(NZ_ext, NY_ext, NX_ext);rnrn tzz = space3d(NZ_ext, NY_ext, NX_ext);rnrn txy = space3d(NZ_ext, NY_ext, NX_ext);rnrn txz = space3d(NZ_ext, NY_ext, NX_ext);rnrn tyz = space3d(NZ_ext, NY_ext, NX_ext);rn source= space3d(NZ_ext, NY_ext, NX_ext);rnrn /* tzz_x=space2d(NZ_ext,NX_ext);rn tzz_z=space2d(NZ_ext,NX_ext);rn txz_x=space2d(NZ_ext,NX_ext);rn txz_z=space2d(NZ_ext,NX_ext);rn ssx=space2d(NZ_ext,NX_ext);rn ssz=space2d(NZ_ext,NX_ext);*/rnrn dxi = space3d(NZ_ext, NY_ext, NX_ext);rn dxi2 = space3d(NZ_ext, NY_ext, NX_ext);rn dyj = space3d(NZ_ext, NY_ext, NX_ext);rn dyj2 = space3d(NZ_ext, NY_ext, NX_ext);rn dzk = space3d(NZ_ext, NY_ext, NX_ext);rn dzk2 = space3d(NZ_ext, NY_ext, NX_ext);rn char dxi_name[] = "dxi.dat";rn char dzk_name[] = "dzk.dat";rn //(xoleft,xoright)模型的物理边界—不算PMLrn xoleft = DP;rn xoright = (NX_ext-1)*H - DP;rn yoleft = DP;rn yoright = (NY_ext-1 )*H - DP;rn zoleft = DP;rn zoright = (NZ_ext-1 )*H - DP;rnrnrn //用于对vx_x[iz][ix]等的求解,加入吸收边界,使数值衰减rn for (iz = 0; iz < NZ_ext; iz++)rn rn for (iy = 0; iy < NY_ext; iy++)rn rn for (ix = 0; ix < NX_ext; ix++)rn rn x = ix*H; y = iy*H; z = iz*H;rn if ((x>=0&&x < xoleft)&&(y>=0&&y= 0 && z= xoright && x < NX_ext*H) && (y >= 0 && y= 0 && z= 0 && x= 0 && y < yoleft) && (z >= 0 && z= 0 && x= yoright && y < NY_ext*H) &&(z >= 0 && z= 0 && x= 0 && y= 0 && z < zoleft) )rn rn v0 = 1500.0f;rn d0 = 3.0*v0*log(1.0 / RC) / (2.0*DP);rn //d0=3.0*v0*(8.0/15.0-3.0/100.0*NP+NP*NP/1500.0)/H; rn dzk[iz][iy][ix] = d0*pow(((zoleft - z) / DP), 2);rn dzk2[iz][iy][ix] = d0*pow(((zoleft - z - 0.5*H) / DP), 2);rn rn if ((x >= 0 && x= 0 && y= zoright && z < NZ_ext*H) )rn rn v0 = 1500.0f;rn d0 = 3.0*v0*log(1.0 / RC) / (2.0*DP);rn dzk[iz][iy][ix] = d0*pow(((z - zoright) / DP), 2);rn dzk2[iz][iy][ix] = d0*pow(((z + 0.5*H - zoright) / DP), 2);rn rn rn rn rn rn rnrnrn /*for (iz = 0; iz= 0.9999*yoright)rn rn v0 = 1500.0;rn d0 = 3.0*v0*log(1.0 / RC) / (2.0*DP);rn dyj[iz][iy][ix] = d0*pow(((y - yoright) / DP), 2);rn dyj2[iz][iy][ix] = d0*pow(((y + 0.5*H - yoright) / DP), 2);rn rn elsern rn dyj[iz][iy][ix] = 0.; dyj2[iz][iy][ix] = 0.;rn rn rn rn rnrn for (iz = 0; iz= 0.9999*zoright)rn rn v0 = 1500.0;rn d0 = 3.0*v0*log(1.0 / RC) / (2.0*DP);rn dzk[iz][iy][ix] = d0*pow(((z - zoright) / DP), 2);rn dzk2[iz][iy][ix] = d0*pow(((z + 0.5*H - zoright) / DP), 2);rn rn elsern rn dzk[iz][iy][ix] = 0.; dzk2[iz][iy][ix] = 0.;rn rn rn rn */rn wfile3d(dxi_name, dxi, NZ_ext, NY_ext, NX_ext);rn wfile3d(dzk_name, dzk, NZ_ext, NY_ext, NX_ext);//建立文件rn //开始时间递推rn float a1 = 9.0 / 8.0;rn float a2 = -1.0 / 24.0;rn //pml初始参数rn rn float xFHalfTemp1;rn float xFHalfTemp2 , xFIntTemp1, xFIntTemp2, yFHalfTemp1, yFHalfTemp2, yFIntTemp1, yFIntTemp2, zFHalfTemp1, zFHalfTemp2, zFIntTemp1, zFIntTemp2 ;rn float ***pmlxSxx, ***pmlySxy, ***pmlzSxz, ***pmlxSxy, ***pmlySyy, ***pmlzSyz, ***pmlxSxz, ***pmlySyz, ***pmlzSzz;rn pmlxSxx= space3d(NZ_ext, NY_ext, NX_ext);rn pmlySxy= space3d(NZ_ext, NY_ext, NX_ext);rn pmlzSxz= space3d(NZ_ext, NY_ext, NX_ext);rn pmlxSxy= space3d(NZ_ext, NY_ext, NX_ext);rn pmlySyy= space3d(NZ_ext, NY_ext, NX_ext);rn pmlzSyz= space3d(NZ_ext, NY_ext, NX_ext);rn pmlxSxz= space3d(NZ_ext, NY_ext, NX_ext);rn pmlySyz= space3d(NZ_ext, NY_ext, NX_ext);rn pmlzSzz= space3d(NZ_ext, NY_ext, NX_ext);rn float ***pmlxVx, ***pmlyVy, ***pmlzVz, ***pmlzVx, ***pmlxVz, ***pmlyVx, ***pmlxVy, ***pmlzVy, ***pmlyVz;rn pmlxVx = space3d(NZ_ext, NY_ext, NX_ext);rn pmlyVy = space3d(NZ_ext, NY_ext, NX_ext);rn pmlzVz = space3d(NZ_ext, NY_ext, NX_ext);rn pmlzVx = space3d(NZ_ext, NY_ext, NX_ext);rn pmlxVz = space3d(NZ_ext, NY_ext, NX_ext);rn pmlyVx = space3d(NZ_ext, NY_ext, NX_ext);rn pmlxVy = space3d(NZ_ext, NY_ext, NX_ext);rn pmlzVy = space3d(NZ_ext, NY_ext, NX_ext);rn pmlyVz = space3d(NZ_ext, NY_ext, NX_ext);rn float ***DxSxxPre, ***DySxyPre, ***DzSxzPre, ***DxSxyPre, ***DySyyPre, ***DzSyzPre, ***DxSxzPre, ***DySyzPre, ***DzSzzPre;rn DxSxxPre = space3d(NZ_ext, NY_ext, NX_ext);rn DySxyPre = space3d(NZ_ext, NY_ext, NX_ext);rn DzSxzPre = space3d(NZ_ext, NY_ext, NX_ext);rn DxSxyPre = space3d(NZ_ext, NY_ext, NX_ext);rn DySyyPre = space3d(NZ_ext, NY_ext, NX_ext);rn DzSyzPre = space3d(NZ_ext, NY_ext, NX_ext);rn DxSxzPre = space3d(NZ_ext, NY_ext, NX_ext);rn DySyzPre = space3d(NZ_ext, NY_ext, NX_ext);rn DzSzzPre = space3d(NZ_ext, NY_ext, NX_ext);rn float ***DxVxPre, ***DyVyPre, ***DzVzPre, ***DzVxPre, ***DxVzPre, ***DyVxPre, ***DxVyPre, ***DzVyPre, ***DyVzPre;rn DxVxPre = space3d(NZ_ext, NY_ext, NX_ext);rn DyVyPre = space3d(NZ_ext, NY_ext, NX_ext);rn DzVzPre = space3d(NZ_ext, NY_ext, NX_ext);rn DzVxPre = space3d(NZ_ext, NY_ext, NX_ext);rn DxVzPre = space3d(NZ_ext, NY_ext, NX_ext);rn DyVxPre = space3d(NZ_ext, NY_ext, NX_ext);rn DxVyPre = space3d(NZ_ext, NY_ext, NX_ext);rn DzVyPre = space3d(NZ_ext, NY_ext, NX_ext);rn DyVzPre = space3d(NZ_ext, NY_ext, NX_ext);rnrn //这里的代码应该修正一下,应该从数据里面提取最大最小速度rnrn //控制网格频散rn if (Vsmin / (F0*H) < 15)rn printf("1\n");rnrn //检查稳定性条件rn best_dt = 6.0*H / (7.0*sqrt(2.0)*Vpmax);rn if (DT >= best_dt)rn printf("2 %f\n", best_dt);rnrn clock_t start, finish;rn start = clock();rnrn sis_x = space2d(NZ, NT);rn sis_y = space2d(NZ, NT);rn sis_z = space2d(NZ, NT);rn omp_set_num_threads(16);rn #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) rn //$OMP PARALLEL DEFAULT(SHARED) PRIVATE(x, y, z, &rn //$OMP VxSource, VySource, VzSource, SxxSource, SyySource, SzzSource, SxySource, SxzSource, SyzSource, &rn //$OMP DenOReciprocal, DenXReciprocal, DenYReciprocal, DenZReciprocal, LambdaO, MuO, MuXY, MuXZ, MuYZ, &rn //$OMP DxSxx, DySxy, DzSxz, DxSxy, DySyy, DzSyz, DxSxz, DySyz, DzSzz, &rn //$OMP xAttFactorHalfTemp1, xAttFactorIntTemp1, yAttFactorHalfTemp1, yAttFactorIntTemp1, zAttFactorHalfTemp1, zAttFactorIntTemp1, &rn //$OMP xAttFactorHalfTemp2, xAttFactorIntTemp2, yAttFactorHalfTemp2, yAttFactorIntTemp2, zAttFactorHalfTemp2, zAttFactorIntTemp2, &rn //$OMP BoolTemp1, BoolTemp2, BoolTemp3, BoolTemp4, BoolTemp5, BoolTemp6, &rn //$OMP DxVx, DyVy, DzVz, DxVz, DzVx, DxVy, DyVx, DyVz, DzVy)rnrnrn for (it = 0; it < NT; it++)rn rn tt = it*DT;rn if (it % 100 == 0) printf("it=%d\n", it);rnrn //计算质点速度vrn #pragma omp forrn for (iz = 2; iz < NZ_ext - 2; iz++)rn for (iy = 2; iy < NY_ext - 2; iy++)rn for (ix = 2; ix < NX_ext - 2; ix++)rn rn rho_tempx = 2 / (rho_ext[iz][iy][ix] + rho_ext[iz][iy][ix + 1]);rn rho_tempy = 2 / (rho_ext[iz][iy][ix] + rho_ext[iz][iy + 1][ix]);rn rho_tempz = 2 / (rho_ext[iz][iy][ix] + rho_ext[iz + 1][iy][ix]);rnrn //4阶rn /*xx1 = ((a1*(txx[iz][iy][ix + 1] - txx[iz][iy][ix]) +rn a2*(txx[iz][iy][ix + 2] - txx[iz][iy][ix - 1]))) / H;rn xy1 = ((a1*(txy[iz][iy][ix] - txy[iz][iy - 1][ix]) +rn a2*(txy[iz][iy + 1][ix] - txy[iz][iy - 2][ix]))) / H;rn xz1 = ((a1*(txz[iz][iy][ix] - txz[iz - 1][iy][ix]) +rn a2*(txz[iz + 1][iy][ix] - txz[iz - 2][iy][ix]))) / H;rnrn yx2 = ((a1*(txy[iz][iy][ix] - txy[iz][iy][ix - 1]) +rn a2*(txy[iz][iy][ix + 1] - txy[iz][iy][ix - 2]))) / H;rn yy2 = ((a1*(tyy[iz][iy + 1][ix] - tyy[iz][iy][ix]) +rn a2*(tyy[iz][iy + 2][ix] - tyy[iz][iy - 1][ix]))) / H;rn yz2 = ((a1*(tyz[iz][iy][ix] - tyz[iz - 1][iy][ix]) +rn a2*(tyz[iz + 1][iy][ix] - tyz[iz - 2][iy][ix]))) / H;rnrn zx3 = ((a1*(txz[iz][iy][ix] - txz[iz][iy][ix - 1]) +rn a2*(txz[iz][iy][ix + 1] - txz[iz][iy][ix - 2]))) / H;rn zy3 = ((a1*(tyz[iz][iy][ix] - tyz[iz][iy - 1][ix]) +rn a2*(tyz[iz][iy + 1][ix] - tyz[iz][iy - 2][ix]))) / H;rn zz3 = ((a1*(tzz[iz + 1][iy][ix] - tzz[iz][iy][ix]) +rn a2*(tzz[iz + 2][iy][ix] - tzz[iz - 1][iy][ix]))) / H;*/rn //二阶rn xx1 = (txx[iz][iy][ix + 1] - txx[iz][iy][ix]) / H;rn xy1 = (txy[iz][iy][ix] - txy[iz][iy - 1][ix]) / H;rn xz1 = (txz[iz][iy][ix] - txz[iz - 1][iy][ix]) / H;rn yx2 = (txy[iz][iy][ix] - txy[iz][iy][ix - 1]) / H;rn yy2 = (tyy[iz][iy + 1][ix] - tyy[iz][iy][ix]) / H;rn yz2 = (tyz[iz][iy][ix] - tyz[iz - 1][iy][ix]) / H;rn zx3 = (txz[iz][iy][ix] - txz[iz][iy][ix - 1]) / H;rn zy3 = (tyz[iz][iy][ix] - tyz[iz][iy - 1][ix]) / H;rn zz3 = (tzz[iz + 1][iy][ix] - tzz[iz][iy][ix]) / H;rnrn /*x = ix*H; y = iy*H; z = iz*H;rn if(((x >= 0 && x < xoleft) && (y >= 0 && y= 0 && z= xoright && x < NX_ext*H) && (y >= 0 && y= 0 && z= 0 && x= 0 && y < yoleft) && (z >= 0 && z= 0 && x= yoright && y < NY_ext*H) && (z >= 0 && z= 0 && x= 0 && y= 0 && z < zoleft))rn rn ||((x >= 0 && x= 0 && y= zoright && z < NZ_ext*H)))*/rn if((ix=(NX_ext-1-NP))||(iy=(NY_ext-1-NP))||(iz=(NZ_ext-1-NP)))rn rn xFHalfTemp1 = exp(-(dxi2[iz][iy][ix])*DT); xFIntTemp1 = exp(-(dxi[iz][iy][ix])*DT);rn yFHalfTemp1 = exp(-(dyj2[iz][iy][ix])*DT); yFIntTemp1 = exp(-(dyj[iz][iy][ix])*DT);rn zFHalfTemp1 = exp(-(dzk2[iz][iy][ix])*DT); zFIntTemp1 = exp(-(dzk[iz][iy][ix])*DT);rn xFHalfTemp2 = -DT*dxi2[iz][iy][ix] * 0.5; xFIntTemp2 = -DT*dxi[iz][iy][ix] * 0.5;rn yFHalfTemp2 = -DT*dyj2[iz][iy][ix] * 0.5; yFIntTemp2 = -DT*dyj[iz][iy][ix] * 0.5;rn zFHalfTemp2 = -DT*dzk2[iz][iy][ix] * 0.5; zFIntTemp2 = -DT*dzk[iz][iy][ix] * 0.5;rn pmlxSxx[iz][iy][ix] = xFHalfTemp1*pmlxSxx[iz][iy][ix] + xFHalfTemp2*(xFHalfTemp1*DxSxxPre[iz][iy][ix] + xx1);rn pmlySxy[iz][iy][ix] = yFIntTemp1*pmlySxy[iz][iy][ix] + yFIntTemp2*(yFIntTemp1*DySxyPre[iz][iy][ix] + xy1);rn pmlzSxz[iz][iy][ix] = zFIntTemp1*pmlzSxz[iz][iy][ix] + zFIntTemp2*(zFIntTemp1*DzSxzPre[iz][iy][ix] + xz1);rn DxSxxPre[iz][iy][ix] = xx1; DySxyPre[iz][iy][ix] = xy1; DzSxzPre[iz][iy][ix] = xz1;rn xx1 = xx1 + pmlxSxx[iz][iy][ix];rn xy1 = xy1 + pmlySxy[iz][iy][ix];rn xz1 = xz1 + pmlzSxz[iz][iy][ix];rnrn pmlxSxy[iz][iy][ix] = xFIntTemp1*pmlxSxy[iz][iy][ix] + xFIntTemp2*(xFIntTemp1*DxSxyPre[iz][iy][ix] + yx2);rn pmlySyy[iz][iy][ix] = yFHalfTemp1*pmlySyy[iz][iy][ix] + yFHalfTemp2*(yFHalfTemp1*DySyyPre[iz][iy][ix] + yy2);rn pmlzSyz[iz][iy][ix] = zFIntTemp1*pmlzSyz[iz][iy][ix] + zFIntTemp2*(zFIntTemp1*DzSyzPre[iz][iy][ix] + yz2);rn DxSxyPre[iz][iy][ix] = yx2; DySyyPre[iz][iy][ix] = yy2; DzSyzPre[iz][iy][ix] = yz2;rn yx2 = yx2 + pmlxSxy[iz][iy][ix];rn yy2 = yy2 + pmlySyy[iz][iy][ix];rn yz2 = yz2 + pmlzSyz[iz][iy][ix];rnrn pmlxSxz[iz][iy][ix] = xFIntTemp1*pmlxSxz[iz][iy][ix] + xFIntTemp2*(xFIntTemp1*DxSxzPre[iz][iy][ix] + zx3);rn pmlySyz[iz][iy][ix] = yFIntTemp1*pmlySyz[iz][iy][ix] + yFIntTemp2*(yFIntTemp1*DySyzPre[iz][iy][ix] + zy3);rn pmlzSzz[iz][iy][ix] = zFHalfTemp1*pmlzSzz[iz][iy][ix] + zFHalfTemp2*(zFHalfTemp1*DzSzzPre[iz][iy][ix] + zz3);rn DxSxzPre[iz][iy][ix] = zx3; DySyzPre[iz][iy][ix] = zy3; DzSzzPre[iz][iy][ix] = zz3;rn zx3 = zx3 + pmlxSxz[iz][iy][ix];rn zy3 = zy3 + pmlySyz[iz][iy][ix];rn zz3 = zz3 + pmlzSzz[iz][iy][ix];rn rn rn vx[iz][iy][ix] = vx[iz][iy][ix] + (DT*rho_tempx)*(xx1 + xy1 + xz1);rn vy[iz][iy][ix] = vy[iz][iy][ix] + (DT*rho_tempy)*(yx2 + yy2 + yz2);rn vz[iz][iy][ix] = vz[iz][iy][ix] + (DT*rho_tempz)*(zx3 + zy3 + zz3);rnrnrnrn rnrn #pragma omp barrierrnrn //计算应力rn #pragma omp forrn for (iz = 2; iz < NZ_ext - 2; iz++)rn for (iy = 2; iy < NY_ext - 2; iy++)rn for (ix = 2; ix < NX_ext - 2; ix++)rn rnrn 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)rn muxz = 0.0;rn elsern 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]));rn 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)rn muyz = 0.0;rn elsern 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]));rn 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)rn muxy = 0.0;rn elsern 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]));rnrnrn lamda2u_temp = lamda2u_ext[iz][iy][ix];rn lamda_temp = lamda_ext[iz][iy][ix];rn // muxz=0.0;rn if ( i==0&&ix == sx&&iy == sy&&iz == sz)rn rn source[iz][iy][ix] = (1 - 2 * PI*PI*F0*F0*(tt - T0)*(tt - T0))*exp(-PI*PI*F0*F0*(tt - T0)*(tt - T0));rn rn else if(i==1&&(ix == (sx+1)&&iy == sy&&iz == sz))rn rn source[iz][iy][ix] = -(1 - 2 * PI*PI*F0*F0*(tt - T0)*(tt - T0))*exp(-PI*PI*F0*F0*(tt - T0)*(tt - T0));rn rn else if (i == 1 && (ix == (sx - 1) && iy == sy&&iz == sz))rn rn source[iz][iy][ix] = (1 - 2 * PI*PI*F0*F0*(tt - T0)*(tt - T0))*exp(-PI*PI*F0*F0*(tt - T0)*(tt - T0));rn rnrn // Qz=Q_ext[iz][ix];R_temp=R_ext[iz][ix];lamda2u_temp=lamda2u_ext[iz][ix];lamda_temp=lamda_ext[iz][ix];rnrn /*y5 = (a1*(vx[iz][iy][ix] - vx[iz][iy][ix - 1]) +rn a2*(vx[iz][iy][ix + 1] - vx[iz][iy][ix - 2])) / H;rnrn y6 = (a1*(vy[iz][iy][ix] - vy[iz][iy - 1][ix]) +rn a2*(vy[iz][iy + 1][ix] - vy[iz][iy - 2][ix])) / H;rnrn y7 = (a1*(vz[iz][iy][ix] - vz[iz - 1][iy][ix]) +rn a2*(vz[iz + 1][iy][ix] - vz[iz - 2][iy][ix])) / H;rnrn y8 = ((a1*(vx[iz][iy + 1][ix] - vx[iz][iy][ix]) +rn a2*(vx[iz][iy + 2][ix] - vx[iz][iy - 1][ix]))) / H;rnrn y9 = ((a1*(vy[iz][iy][ix + 1] - vy[iz][iy][ix]) +rn a2*(vy[iz][iy][ix + 2] - vy[iz][iy][ix - 1]))) / H;rnrn y10 = ((a1*(vx[iz + 1][iy][ix] - vx[iz][iy][ix]) +rn a2*(vx[iz + 2][iy][ix] - vx[iz - 1][iy][ix]))) / H;rnrn y11 = ((a1*(vz[iz][iy][ix + 1] - vz[iz][iy][ix]) +rn a2*(vz[iz][iy][ix + 2] - vz[iz][iy][ix - 1]))) / H;rnrn y12 = ((a1*(vy[iz + 1][iy][ix] - vy[iz][iy][ix]) +rn a2*(vy[iz + 2][iy][ix] - vy[iz - 1][iy][ix]))) / H;rnrn y13 = ((a1*(vz[iz][iy + 1][ix] - vz[iz][iy][ix]) +rn a2*(vz[iz][iy + 2][ix] - vz[iz][iy - 1][ix]))) / H;*/rn y5 = (vx[iz][iy][ix] - vx[iz][iy][ix - 1]) / H;rn y6 = (vy[iz][iy][ix] - vy[iz][iy - 1][ix]) / H;rn y7 = (vz[iz][iy][ix] - vz[iz - 1][iy][ix]) / H;rn y8 = (vx[iz][iy + 1][ix] - vx[iz][iy][ix]) / H;rn y9 = (vy[iz][iy][ix + 1] - vy[iz][iy][ix]) / H;rn y10 = (vx[iz + 1][iy][ix] - vx[iz][iy][ix]) / H;rn y11 = (vz[iz][iy][ix + 1] - vz[iz][iy][ix]) / H;rn y12 = (vy[iz + 1][iy][ix] - vy[iz][iy][ix]) / H;rn y13 = (vz[iz][iy + 1][ix] - vz[iz][iy][ix]) / H;rn //x = ix*H; y = iy*H; z = iz*H;rn /*if (((x >= 0 && x < xoleft) && (y >= 0 && y= 0 && z= xoright && x < NX_ext*H) && (y >= 0 && y= 0 && z= 0 && x= 0 && y < yoleft) && (z >= 0 && z= 0 && x= yoright && y < NY_ext*H) && (z >= 0 && z= 0 && x= 0 && y= 0 && z < zoleft))rnrn || ((x >= 0 && x= 0 && y= zoright && z < NZ_ext*H)))*/rn if((ix=(NX_ext-1-NP))||(iy=(NY_ext-1-NP))||(iz=(NZ_ext-1-NP)))rn rn xFHalfTemp1 = exp(-(dxi2[iz][iy][ix])*DT); xFIntTemp1 = exp(-(dxi[iz][iy][ix])*DT);rn yFHalfTemp1 = exp(-(dyj2[iz][iy][ix])*DT); yFIntTemp1 = exp(-(dyj[iz][iy][ix])*DT);rn zFHalfTemp1 = exp(-(dzk2[iz][iy][ix])*DT); zFIntTemp1 = exp(-(dzk[iz][iy][ix])*DT);rn xFHalfTemp2 = -DT*dxi2[iz][iy][ix] * 0.5; xFIntTemp2 = -DT*dxi[iz][iy][ix] * 0.5;rn yFHalfTemp2 = -DT*dyj2[iz][iy][ix] * 0.5; yFIntTemp2 = -DT*dyj[iz][iy][ix] * 0.5;rn zFHalfTemp2 = -DT*dzk2[iz][iy][ix] * 0.5; zFIntTemp2 = -DT*dzk[iz][iy][ix] * 0.5;rn pmlxVx[iz][iy][ix] = xFIntTemp1*pmlxVx[iz][iy][ix] + xFIntTemp2*(xFIntTemp1*DxVxPre[iz][iy][ix] + y5);rn pmlyVy[iz][iy][ix] = yFIntTemp1*pmlyVy[iz][iy][ix] + yFIntTemp2*(yFIntTemp1*DyVyPre[iz][iy][ix] + y6);rn pmlzVz[iz][iy][ix] = zFIntTemp1*pmlzVz[iz][iy][ix] + zFIntTemp2*(zFIntTemp1*DzVzPre[iz][iy][ix] + y7);rn DxVxPre[iz][iy][ix] = y5; DyVyPre[iz][iy][ix] = y6; DzVzPre[iz][iy][ix] = y7;rn y5 = y5 + pmlxVx[iz][iy][ix];rn y6 = y6 + pmlyVy[iz][iy][ix];rn y7 = y7 + pmlzVz[iz][iy][ix];rnrn pmlzVx[iz][iy][ix] = zFHalfTemp1*pmlzVx[iz][iy][ix] + zFHalfTemp2*(zFHalfTemp1*DzVxPre[iz][iy][ix] + y10);rn pmlxVz[iz][iy][ix] = xFHalfTemp1*pmlxVz[iz][iy][ix] + xFHalfTemp2*(xFHalfTemp1*DxVzPre[iz][iy][ix] + y11);rn DzVxPre[iz][iy][ix] = y10; DxVzPre[iz][iy][ix] = y11;rn y10 = y10 + pmlzVx[iz][iy][ix];rn y11 = y11 + pmlxVz[iz][iy][ix];rnrn pmlyVx[iz][iy][ix] = yFHalfTemp1*pmlyVx[iz][iy][ix] + yFHalfTemp2*(yFHalfTemp1*DyVxPre[iz][iy][ix] + y8);rn pmlxVy[iz][iy][ix] = xFHalfTemp1*pmlxVy[iz][iy][ix] + xFHalfTemp2*(xFHalfTemp1*DxVyPre[iz][iy][ix] + y9);rn DyVxPre[iz][iy][ix] = y8; DxVyPre[iz][iy][ix] = y9;rn y8 = y8 + pmlyVx[iz][iy][ix];rn y9 = y9 + pmlxVy[iz][iy][ix];rnrn pmlzVy[iz][iy][ix] = zFHalfTemp1*pmlzVy[iz][iy][ix] + zFHalfTemp2*(zFHalfTemp1*DzVyPre[iz][iy][ix] + y12);rn pmlyVz[iz][iy][ix] = yFHalfTemp1*pmlyVz[iz][iy][ix] + yFHalfTemp2*(yFHalfTemp1*DyVzPre[iz][iy][ix] + y13);rn DzVyPre[iz][iy][ix] = y12; DyVzPre[iz][iy][ix] = y13;rn y12 = y12 + pmlzVy[iz][iy][ix];rn y13 = y13 + pmlyVz[iz][iy][ix];rn rn rn txx[iz][iy][ix] = txx[iz][iy][ix] + (lamda2u_temp*DT)*y5 + (lamda_temp*DT)*y6 + (lamda_temp*DT)*y7 + source[iz][iy][ix];rn tyy[iz][iy][ix] = tyy[iz][iy][ix] + (lamda_temp*DT)*y5 + (lamda2u_temp*DT)*y6 + (lamda_temp*DT)*y7 + source[iz][iy][ix];rn tzz[iz][iy][ix] = tzz[iz][iy][ix] + (lamda_temp*DT)*y5 + (lamda_temp*DT)*y6 + (lamda2u_temp*DT)*y7 + source[iz][iy][ix];rn txz[iz][iy][ix] = txz[iz][iy][ix] + (muxz*DT)*(y10 + y11);rn txy[iz][iy][ix] = txy[iz][iy][ix] + (muxy*DT)*(y8 + y9);rn tyz[iz][iy][ix] = tyz[iz][iy][ix] + (muyz*DT)*(y12 + y13);rnrn //加震源—单极子rn if (i == 0)rn rn tzz[sz][sy][sx] =(1 -2 * PI*PI*F0*F0*(tt - T0)*(tt - T0))*exp(-PI*PI*F0*F0*(tt - T0)*(tt - T0)); //对正应力加爆炸源rn tyy[sz][sy][sx] = (1-2 * PI*PI*F0*F0*(tt - T0)*(tt - T0))*exp(-PI*PI*F0*F0*(tt - T0)*(tt - T0));rn txx[sz][sy][sx] =(1 -2 * PI*PI*F0*F0*(tt - T0)*(tt - T0))*exp(-PI*PI*F0*F0*(tt - T0)*(tt - T0));rn rn //加震源—偶极子 rnrn if (i == 1)rn rn 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));; //对正应力加爆炸源rn 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));;rn 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));;rn 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));; //对正应力加爆炸源rn 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));;rn 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));;rn rnrn rnrn #pragma omp barrierrn #pragma omp masterrn rn rnrn // tzz_z[sz][sx]=tzz_z[sz][sx]-2*PI*PI*F0*F0*(tt-T0)*exp(-PI*PI*F0*F0*(tt-T0)*(tt-T0));rn //sis_y[1][it] = tzz[sz][sy][sx];rnrn //记录地震记录rn if (i == 0)rn rn for (iz = NP; iz < NP + NZ; iz++)rn rn sis_z[iz - NP][it] = tzz[iz][sy][sx];rnrn sis_x[iz - NP][it] = txx[iz][sy][sx];rn //sis_z[iz-NP][it]=vz_x[iz][NX_ext-iz]+vz_z[iz][NX_ext-iz]; rn // sis_x[iz-NP][it]=vx_x[iz][NX_ext-iz]+vx_z[iz][NX_ext-iz];rn rn rn if (i == 1)rn rn for (iz = NP; iz < NP + NZ; iz++)rn rn sis_z[iz - NP][it] = tzz[iz][sy][sx];rnrn sis_x[iz - NP][it] = txx[iz][sy][sx+1];rn //sis_z[iz-NP][it]=vz_x[iz][NX_ext-iz]+vz_z[iz][NX_ext-iz]; rn // sis_x[iz-NP][it]=vx_x[iz][NX_ext-iz]+vx_z[iz][NX_ext-iz];rn rn rn if (it % 200 == 0)rn rn int n = it / 200;rn for (iz = 2; iz < NZ_ext - 2; iz = iz + 2)rn rn int a = iz / 2;rn for (ix = 2; ix < NX_ext - 2; ix = ix + 2)rn rn int b = ix / 2;rnrn vxx[n][a][b] = txx[iz][sy][ix];rn rn rn rn rn #pragma omp barrierrn rn rn finish = clock();rn printf("%f seconds\n", (double)(finish - start) / CLOCKS_PER_SEC);rnrn //输出波场快照rn char vzname[] = "vz.dat";rn char vxname[] = "vx.dat";rn char vxxname[] = "vxx.dat";rn char vyname[] = "vy.dat";rn char txx_xname[] = "txx_x.dat";rn char pmlxSxxname[] = "pmlxSxx.dat";rn char DxSxxPrename[] = "DxSxxPre.dat";rn /* for (iz=0; iz=(nr+14-ix))//45度倾斜井孔rn//if (iz>= 50)rn //if(ix>=50)rn rn vp[iz][iy][ix] = 1500.0f;rn vs[iz][iy][ix] = 0.0f;rn vf[iz][iy][ix] = 1500.0f;rn rho[iz][iy][ix] = 1000.0f;rn rhof[iz][iy][ix] = 1000.0f;rn por[iz][iy][ix] = 1.0f;rn // b[iz][ix]=0.0;rnrn mu[iz][iy][ix] = rho[iz][iy][ix] * vs[iz][iy][ix] * vs[iz][iy][ix];rn lamda2u[iz][iy][ix] = rho[iz][iy][ix] * vp[iz][iy][ix] * vp[iz][iy][ix];rn lamda[iz][iy][ix] = lamda2u[iz][iy][ix] - 2 * mu[iz][iy][ix];rn /* Kb1=rho[iz][ix]*(1-por[iz][ix])*(vp[iz][ix]*vp[iz][ix]-vs[iz][ix]*vs[iz][ix]*4.0/3.0);//骨架压缩模量,声波测井原理与应用,P39rn Ks1=rho[iz][ix]*(vp[iz][ix]*vp[iz][ix]-vs[iz][ix]*vs[iz][ix]*4.0/3.0);//岩石固态颗粒的体积模量rn Kf1=rhof[iz][ix]*vf[iz][ix]*vf[iz][ix];//孔隙流体的体积压缩模量rn a1=1-Kb1/Ks1;rn D1=a1-por[iz][ix]+por[iz][ix]*Ks1/Kf1;rn Q[iz][ix]=(a1-por[iz][ix])*por[iz][ix]*Ks1/D1;rn R[iz][ix]=por[iz][ix]*por[iz][ix]*Ks1/D1;rnrn tao=0.5*(1+1/por[iz][ix]);//孔隙弯曲度rn rou11=(1-por[iz][ix])*rho[iz][ix]-(1-tao)*por[iz][ix]*rhof[iz][ix];rn rou12=(1-tao)*por[iz][ix]*rhof[iz][ix];rn rou22=tao*por[iz][ix]*rhof[iz][ix];rnrn D11[iz][ix]=rou11/(rou11*rou22-rou12*rou12);rn D12[iz][ix]=rou12/(rou11*rou22-rou12*rou12);rn D22[iz][ix]=rou22/(rou11*rou22-rou12*rou12);*/rn rn elsern rnrn vp[iz][iy][ix] = 4000.0f; //纵波rn vs[iz][iy][ix] = 2340.0f; //横波rn rho[iz][iy][ix] = 2500.0f; //密度rn vf[iz][iy][ix] = 1500.0f;rn rhof[iz][iy][ix] = 1000.0f;rn por[iz][iy][ix] = 0.2f;rn //b[iz][ix]=0.2*pow(10,-3);// b=eta*por*por/perm*pow(10,-3);rnrn mu[iz][iy][ix] = rho[iz][iy][ix] * vs[iz][iy][ix] * vs[iz][iy][ix];rn lamda2u[iz][iy][ix] = rho[iz][iy][ix] * vp[iz][iy][ix] * vp[iz][iy][ix];rn lamda[iz][iy][ix] = lamda2u[iz][iy][ix] - 2.0*mu[iz][iy][ix];rnrn rnrn rn rnrn rnrnrnrn/*rn//将模型扩边,用于PMLrn//具体的操作过程是将实际模型参数放置在扩边后的数据中央,四周的数据用rn//最外缘的数据填充rnfloat **extmodel(float **inqit_model,int nz,int nx,int np)rnrnfloat **p;rnint i,j;rnint nx2=nx+2*np;rnint nz2=nz+2*np;rnrnp=space2d(nz2,nx2);rnfor (i=np; i
openmp openmp
openmpopenmpopenmpopenmpopenmpopenmpopenmpopenmpopenmpopenmpopenmpopenmp
openMP程序
openmp程序,需要的快下载阿,多核并行算法很好用的!mpi更适合集群计算,openmp是多核并行计算的首选
openMp
你好,我用的intel c9.0+vc6,我用OpenMp并行我的迭代,为什么效果不明显啊,和串行差不多
openmp
用vc+2005和openmp开发的多线程优化了某一算法生成了dll,请问该dll在运行时依赖于哪些东西,比如frame work 2.0?还有Common Language Runtime?谢谢!
OpenMP
我想请问下,在多核计算机中,用什么样的命令调试运行openmp程序呢?
并行15puzll openMP C程序
15puzzle openMP C西安电子科技大学 霍老师并行计算程序作业第5题 开4个进程才能运行无错.
OPENMP问题
用服务器做测试。rn4颗CPU共32核。rnrn用omp_set_num_threads()设置线程为32,64,128,256rn然后后面rn omp_set_num_threads(32);//设置线程为32,64,128,256rn #pragma omp parallel for rn for()rn rn rn无论设置为哪个,看性能管理器,CPU利用率只有5%。无法充分复利用性能,请用需要怎么才能达到利用率接近100%的效果?rn[img=http://img.my.csdn.net/uploads/201301/30/1359511920_4387.jpg][/img]rnrnrn又自己看了下,每次运行调用的不同的8个线程。感觉是只掉用成功了一个核的8个线程。剩下3个CPU没调用成功。
openmp的问题
#pragma omp parallel for rnfor(int i = 0 ; i < 10; i++)rnrn如果是双核的机器上跑,怎么让第一个线程跑i = 0~3, 第二个线程跑 4~9?
OpenMP 运行库函数
运行库函数就是运行时会用到的库函数,库函数是直接调用的,是被封装起来的,我们看不到。 因为属于完成老师布置的翻译,所以,写在这里分享一下。              OMP_SET_NUM_THREADS    目的:      设置下一个并行区域使用的线程数。必须是正整数。 说明:    1. 动态线程机制的修改影响此函数。  2.可行的:指定线程的最大数量,可
OpenMP并行程序设计之OpenMP使用入门
OpenMP并行程序设计之OpenMP使用入门
OpenMP和manifest,了解openmp的请进
rn在控制台程序中练习openMP的例子,编译可以通过,但是运行时出错。(VS2005 + XP SP2)rnrnint TestOpenMP()rnrn #pragma omp parallel forrn for (int i = 0; i < 10; i++ )rn rn printf("i = %d\n", i);rn rn return 0;rnrnrnrn除了弹出错误信息如下:应用程序正常初始化(0xc0000142)失败.请单击确定终止应用程序。rnrnrn在程序的控制台中还有如下错误信息,rnError: r6034rnAn application has made an attempt to load vcompd.dll without using a manifest.rnthis is an unsupported way to load Visual C++ DLLs. You need to modify your applicationrnto build with a manifest.rnFor more information , see the "Visual C++ Libraries as shared Side-by-Side Assembles" rntopic in the product documentation.
还可以使用OpenMP
还可以使用OpenMP来实现程序的并行运算
C语言多线程 openmp
C语言程序的openmp多线程处理,需要对每个线程都要创建一套自己私有的全局指针变量(改指针变量是跨文件的,而且里面包括各种指针、结构体),即使用threadprivate声明各个线程中的变量也是同一块内存,怎样才能各自独立的使用这些全局指针?
两个openmp简单问题,急!
看到一篇博客上介绍omp举例,请问这样递归中用pragma不会递归创建无数的线程吗:rn[code=C/C++]void QuickSort (int numList[], int nLower, int nUpper) rn rn if (nLower < nUpper) rn rn // create partitionsrn int nSplit = Partition (numList, nLower, nUpper);rn #pragma omp parallel sectionsrn rn #pragma omp sectionrn QuickSort (numList, nLower, nSplit - 1);rn rn #pragma omp sectionrn QuickSort (numList, nSplit + 1, nUpper);rn rn rn[/code]
OPENMP与GPU并行的问题!!急
[img=https://img-bbs.csdn.net/upload/201506/30/1435647578_653088.png][/img]rn其中,每个图像卷积操作都是一样的函数,只是图像不同。使用的是同一个GPU即GTX960.rn图像大小为640*480.rn现在的问题是,我只对一路的图像进行卷积操作时,时间大概在10ms。rn但是当四路并行处理时,每路的时间就会增加到35ms左右,而且还会有较大的波动,从20ms到45ms不等的波动。rn我现在怀疑是GPU的计算单元即SM SP同时对四路进行处理,造成了等待延迟造成的。不知道有没有大神能帮我看看,问题到底出在哪里?急rn
OpenMp并行编程 程序
该工程包括三个文件:stdafa.h;OpenMp.cpp;stdafx.cpp。主要内容在OpenMp.cpp文件中。实现的功能是:在一个元素个数为10000的随机数组中(元素值为0—9之间),在不同的线程数时,查找某一个数字出现的次数和查找所用到的时间。
openMP 小程序改错
写了一个计算3在数组里面出现的次数,现在的问题是不能多线程运算,而且统计出来的3的次数是始终不变的,求高手帮忙看一下rnrn#includern#includern#includern#define N 10000rn rn rnint main(int argc, char *argv[])rn int i; int array[N];rnint count,nthreads, tid,chunk;rnchunk = 100;rncount = 0;rndouble start, end;rn rnrnstart = omp_get_wtime();rn#pragma opm parallel shared(nthreads,chunk)rnrn tid=omp_get_thread_num();rnif(tid==0)rn rn nthreads=omp_get_num_threads();rn printf("Starting count with %d threads\n",nthreads);rnrnrn#pragma omp for schedule (static,chunk)rnfor( i = 0; irn#includern#includern#define N 10000rn rn rnint main(int argc, char *argv[])rn int i; int array[N];rnint count,nthreads, tid,chunk;rnchunk = 100;rncount = 0;rndouble start, end;rn rnrnstart = omp_get_wtime();rn#pragma opm parallel shared(nthreads,chunk)rnrn tid=omp_get_thread_num();rnif(tid==0)rn rn nthreads=omp_get_num_threads();rn printf("Starting count with %d threads\n",nthreads);rnrnrn#pragma omp for schedule (static,chunk)rnfor( i = 0; i
OpenMP入门程序
1.源文件:omp.cc#include<iostream> using namespace std; //#include<omp.h>int main() { #pragma omp parallel for(int i=0;i<4;i++) cout<<i<<endl; return 0; }2.编译文件:test.shg++ -fopenmp -o run
第一个OpenMP程序
[size=medium]今天开始学习并行程序的开发,刚刚知道Visual Studio 2008里的Visual C++就可以实现,于是在里面写了个程序。 并设置了C/C++-&gt;Language-&gt;将OpenMP支持改为“是/(OpenMP)”。 于是欣喜若狂的Run,结果…… 出了一个错误的对话框: 没有找到VCOMP90D.DLL,因此这个应用程序未能启动。重新安装应...
OpenMP例子程序
这是OpenMP的例子程序,只要仔细研读,很快就能入门,不需要看网上所谓的教程。
最简单的OpenMP程序
最简单的OpenMP程序 在Win7平台VS2010下编译成功
第一个OPENMP程序
环境:mac OSX gcc 版本: 4.9gcc -v :gcc version 4.9.2 20141029 (prerelease) (GCC) 如果gcc版本是4.2,直接编译含omp.h的文件时,会报错(No such file or directory)安装4.9的具体步骤如下:http://stackoverflow.com/questions/20340117/omp-h-lib
关于openmp程序计时的问题?
在一个无操作系统的片上执行一个openmp程序,在计时的时候出现了错误:rn 例如:rn double start = clock();rn rn /* 中间计算时间的程序*/rn .......rn #pragram omp parallel forrn for( ...) rn ........rn double end = clock();rn 在最后得出的结果time = end - start 是一个溢出的结果,也就是说end 比 start要小,程序执行的时间很短,不可能clock时间会轮回一次的,rn请教大家为什么????
win7 codeblock下openmp运行出错
代码如下rn#include "stdio.h"rn#include "omp.h"rn#include "tchar.h"rn//#include rnrnrnint _tmain(int argc,_TCHAR*argv[])rnrn #pragma omp parallelrn for(int i=0;i<5;i++)rn printf("hello world i=%d\n",i);rn return 0;rnrn编译时rnrnmingw32-g++.exe -Wall -fexceptions -g -fopenmp -c G:\11Download\10\main.cpp -o obj\Debug\main.ornmingw32-g++.exe -o bin\Debug\10.exe obj\Debug\main.o -lgomp -lpthread rnOutput size is 25.83 KBrnProcess terminated with status 0 (0 minutes, 0 seconds)rn0 errors, 0 warningsrnrn运行结果rnhello world i=0rnhello world i=1rnhello world i=2rnhello world i=3rnhello world i=4rnhello world i=0rnhello world i=1rnhello world i=2rnhello world i=3rnhello world i=4rn这是哪出了问题啊,要如何设置
VS2005运行简单OpenMP出错。
已在project的property页面上的C/C++ -> Language 选项卡中将OpenMP Support这项设置为Yes,也添加了#include ,编译通过,但是在执行时出错,说是找不到vcompd.dll。拷贝dll到debug目录下也不行。
Openmp中的数据处理子句
1. privateprivate 子句是将一个或者多个变量声明为线程私有的变量, 声明后指定的线程有他的私有副本,其他线程无法访问, 即使在并行程序外由同名的变量也不会起任何作用, 而且并行程序中的变量不会共享到并行程序之外。code: (c++) #include <iostream> #include <omp.h>using namespace std;int main() { in
OpenMp入门
1.介绍 OpenMp提供了对并行算法的高层的抽象描述,程序员通过在源代码中加入专用的pragma来指明自己的意图,由此编译器可以自动将程序进行并行化,并在必要之处加入同步互斥以及通信。当选择忽略这些pragma,或者编译器不支持OpenMp时,程序又可退化为通常的程序(一般为串行),代码仍然可以正常运作,只是不能利用多线程来加速程序执行。 2.配置 下面是两种配置方式 (1)CodeBlo
openmp讲义
本文档为OpenMP讲义 详细介绍了OpenMP编程相关知识 这也是我之前学OpenMP所选择的资源
C++ and OpenMP
关于如何在C++环境中使用OpenMP的一个小资料
OpenMP资料收集
OpenMP资料收集
OpenMP并行程序设计
详细地讲解了共享存储编程技术。对于学习OpenMP编程的同学来说,有很大的帮助。对关键部分已经做了提示。
OpenMP求素数
用OpenMP求一定区间内素数个数,代码较简单,适合初学者学习,欢迎下载
OPENMP程序设计
OPENMP 程序设计 OpenMP概述 OpenMP编程简介 OpenMP计算实例
OpenMP 程序设计
OpenMP 程序设计OpenMP 程序设计
openMP的ppt
上课时用的关于openmp的ppt简单通俗易懂
openmp命令
我在linux下安得gcc,其中包含openmp编译器。我用vi编写了一个hello.c程序rn不知道用什么命令编译和运行?
openMP疑问
openMP是什么?和C有关系吗?
openmp测试程序
#include &amp;lt;stdio.h&amp;gt;#include &amp;lt;unistd.h&amp;gt;#include &amp;lt;omp.h&amp;gt;void main() {  int j;  //#pragma omp parallel for  //#pragma omp taskwait  //#pragma omp barrier  //#pragma omp for ordered   //#...
OpenMP初试
体验了一把OpenMP,从运行结果上看使用OpenMP在多核机器上,性能还是有大幅提升的。下面上方法。   在Visual Studio Community 2013中创建Win32命令行程序。代码如下: #include #include #include #include using namespace std; int main() { clock_t startTime
相关热词 c# 去空格去转义符 c#用户登录窗体代码 c# 流 c# linux 可视化 c# mvc 返回图片 c# 像素空间 c# 日期 最后一天 c#字典序排序 c# 截屏取色 c#中的哪些属于托管机制