基于MPICH的VTI介质有限差分正演与逆时偏移

//#############################################################a //##s //##s 2D Acoustic VTI Medium RTM & Pick Up ADCIG //##s //##s---------------------------------------------------------- //##s PML , P&SV , VTI , Acoustic , RTM , -SV , Adcig, //##s Migration , Both sides receive , SU or dat(v,e,d) , //##s Adcig( poynting , s-r , offset-domain ), //##s //##s---------------------------------------------------------- //##s Rong Tao //##s 2016 //############################################################a #include<stdio.h> #include<malloc.h> #include<math.h> #include<stdlib.h> #include "mpi.h" #include "/home/Toa/hc/alloc.c" #include "/home/Toa/hc/alloc.h" #include "/home/Toa/hc/complex.c" #include "/home/Toa/hc/complex.h" #include "/home/Toa/hc/bhdr.h" #include "/home/Toa/hc/hdr.h" #include "/home/Toa/hc/fft.h" #include "/home/Toa/hc/fft.c" #include "/home/Toa/hc/mute_direct.h" #include "/home/Toa/hc/cjbsegy.h" /********* SEG-Y header *********/ typedef cjbsegy segy; /********** SU & SEG-Y **********/ #define SU_NKEYS 80 /* Number of key header words */ #define HDRBYTES 240 /* Bytes in the trace header */ #define EBCBYTES 3200 /* Bytes in the card image EBCDIC block */ #define BNYBYTES 400 /* Bytes in the binary coded block */ #define SEGY_HDRBYTES 240 /* Bytes in the tape trace header */ #define SEGY_NKEYS 71 /* Number of mandated header fields */ #define BHED_NKEYS 27 /* Number of mandated binary fields */ #define pi 3.1415926535898 /********* SEG-Y header *********/ Y_3200 y3200; bhed head_400; cjbsegy tr, vtr; /********* SEG-Y header *********/ /********** SU function *********/ void swap_short_2(short *tni2); void swap_u_short_2(unsigned short *tni2); void swap_int_4(int *tni4); void swap_u_int_4(unsigned int *tni4); void swap_long_4(long *tni4); void swap_u_long_4(unsigned long *tni4); void swap_float_4(float *tnf4); void swap_double_8(double *tndd8); void swaphval(segy *tr, int index); /********* SU function *********/ /* MAIN */ main(int argc,char *argv[]) { /*************************************** function ********************************************/ void model_vti_get_boundry(int nx,int nz,int vnx,int vnz,int nt,int npd,float dx,float dz, float vdx,float vdz,float favg,float tmax,float dt,float dtout,float pfac, float **vp,float **epsilu,float **deta,float **rho, int ns_sxd,int ds_sxd,int fs_sxd,int zs_sxd,int is, float **p_cal_x,float **p_cal_z, float **p_top_x,float **p_bottom_x,float **p_left_x,float **p_right_x, float **p_top_z,float **p_bottom_z,float **p_left_z,float **p_right_z,int _Circle_, int mm,int wtype,int hsx,int myid,float *mu_v,int flag_snap,int seismic); void mute_directwave(int flag_mu,int nx,int nt,float dt,float favg, float dx,float dz,int fs_sxd,int ds_sxd,int zs_sxd,int is, float mu_v,float **p_cal,int tt); void RTM_corr_adcig(int nx,int nz,int vnx,int vnz,int nt,int npd,float dx,float dz, float vdx,float vdz,float favg,float tmax,float dt,float dtout, float pfac,float **vp,float **epsilu,float **deta,float **rho,char FN5[], int ns_sxd,int ds_sxd,int fs_sxd,int ds_initial,int fs_initial,int zs_sxd,int is,int flag_cdp, float **p_top_x,float **p_bottom_x,float **p_left_x,float **p_right_x,//////////////////////////////// float **p_top_z,float **p_bottom_z,float **p_left_z,float **p_right_z,float **p_obs_x,float **p_obs_z, int mm,int wtype,int hsx,int myid,float **mig_is,float **mig_ns0, float ***adcig_is,float ***adcig_ns0,float ***Ixhz_is,float ***Ixhz_ns0,int nh, int flag_snap,int seismic,int flag_adcig); void smooth1float(float *v,int r,int n); void smooth2float(int nx,int rx,int nz,int rz,float **v); void pad_vv(int nx,int nz,int npd,float **ee); void read_v_e_d_r(char FN1[],char FN2[],char FN3[],int nx,int nz,float **vv,float **epsilu,float **deta, float **rho0,int npd,int seismic); /*************************** parameter statement ***********************/ int i,j,k,l,m,ih,is,nx,nz,nt,vnx,vnz,i_start,i_end,mm,wtype,hsx,ia,nxs,flag_cdp,flag_adcig,nh; int ns_sxd,ds_sxd,fs_sxd,zs_sxd,fs,ds,npd; float dx,dz,vdx,vdz,tmax,dt,dtout,pfac,favg; int myid,numprocs,_Circle_,flag_mu,flag_snap,seismic,nangle,dangle,fangle; float mu_v; /***************** wave float **************/ float **p_cal_x,**p_cal_z,**p_top_x,**p_bottom_x,**p_left_x,**p_right_x; float **p_top_z,**p_bottom_z,**p_left_z,**p_right_z; /***************** image float *************/ float **mig_is,**mig_ns,**mig_ns0,***adcig_is,***adcig_ns,***adcig_ns0,**adcig2d; float ***Ixhz_ns0,***Ixhz_ns,***Ixhz_is; float **vp,**rho,**deta,**epsilu; float **vps,**rhos,**detas,**epsilus; //a########################################################################### //#### input the parameter and confirm #### //a########################################################################### /************************ dat document **********************/ /* Input velocity */char FN1[250]={"vel_1801_862.dat"}; /* Input epsilu */char FN2[250]={"eps_1801_862.dat"}; /* Input deta */char FN3[250]={"del_1801_862.dat"}; /* Cal data */char FN4[250]={"shot_cal.dat"}; /* Obs data */char FN5[250]={"shot_obs.dat"}; /* Migration */char FN6[250]={"mig_ns.dat"}; /* IS tempor migration */char FN7[250]={"mig_is_tempor.dat"}; /* Adcig initial */char FN8[250]={"adcig.dat"}; /* Adcig smooth */char FN9[250]={"adcig_smooth.dat"}; /* Migration adcig stack */char FN10[250]={"mig_stack_adcig.dat"}; /* Half offset I(x,h,z) */char FN11[250]={"Image_x_h_z.dat"}; /*************************** type ***************************/ /* Wavelet type (1-3) */ wtype=1; /* 1-mute , 0->don't mute */ flag_mu=1; /* 1-snap , 0->nosnap */ flag_snap=0; /* 1.dat, 2.su (v,e,d) */ seismic=1; /*************************** cdp ****************************/ /* Activate "nxs"(=1) */ flag_cdp=1; /* 1-activate ;0-invaliable */ /* CDP of each shot */ nxs=501; /* Must be odd number */ /*************************** cdp ****************************/ /* choice adcig type */ flag_adcig=1; /* 1-poynting */ /* 2-source-receivers domain */ /* 3-half-offset-domain */ /*************************** wave ***************************/ hsx=1;mm=4;npd=20;_Circle_=15; /******************** observation system ********************/ favg=20; pfac=1000.0; nx=1801; dx=10.0; nz=862; dz=5.0; tmax=6.5; dt=0.5;//ms nt=11001; ns_sxd=250; fs_sxd=255; ds_sxd=5; zs_sxd=1; nangle=90; fangle=0; dangle=1; /*************************v****************************/ vdz=dz;vdx=dx;vnx=nx;vnz=nz;dtout=dt; /************************FILE**************************/ FILE *fp4,*fp6,*fp7,*fp8,*fp9,*fp10,*fp11; fp4=fopen(FN4,"wb"); /* Cal data */ fp6=fopen(FN6,"wb"); /* Migration */ fp7=fopen(FN7,"wb"); /* IS tempor migration */ fp8=fopen(FN8,"wb"); /* Adcig initial */ fp9=fopen(FN9,"wb"); /* Adcig smooth */ fp10=fopen(FN10,"wb"); /* Migration adcig stack */ if(flag_adcig==3) fp11=fopen(FN11,"wb"); /* Image half-offset */ /************************** read_file **************