added OpenMP + more cleaning

README

@@ -19,11 +19,13 @@ REFERENCES
 
 Finite-difference modeling experiments for seismic interferometry
 Jan Thorbecke and Deyan Draganov
-2011, Geophysics, Vol. 76, no 6 (November-December); p H1--H18, doi: 10.1190/GEO2010-0039.1
+2011, Geophysics, Vol. 76, no. 6 (November-December); p H1--H18, doi: 10.1190/GEO2010-0039.1
 
 -2- If the Machenko code has helped you in your research please refer to our paper in your publications:
 
-Hopefully, a published reference to the paper can be put here.
+Implementation of the Marchenko method
+Jan Thorbecke, Evert Slob, Joeri Brackenhoff, Joost van der Neut, and Kees Wapenaar
+2017, Geophysics, Vol. 82, no. 6 (November-December); p. WB29--WB45, doi: 10.1190/GEO2017-0108.1
 
 These papers can be downloaded from:
 

raytime/Makefile

@@ -7,8 +7,8 @@ include ../Make_include
 ALLINC  = -I.
 LIBS    += -L$L -lgenfft -lm $(LIBSM)
 #LIBS    += -L$L -lgenfft -lm -lc
-OPTC = -g -Wall -fsignaling-nans -O0
-#OPTC += -fopenmp -Waddress
+#OPTC = -g -Wall -fsignaling-nans -O0
+OPTC += -fopenmp -Waddress
 #OPTC := $(subst -O3 -ffast-math, -O1 -g ,$(OPTC))
 #PGI options for compiler feedback
 #OPTC += -Mprof=lines

raytime/getParameters.c

@@ -20,17 +20,11 @@
 *           The Netherlands
 **/
 
-float gaussGen();
-
-int optncr(int n);
-
 int getModelInfo(char *file_name, int *n1, int *n2, float *d1, float *d2, float *f1, float *f2, float *min, float *max, int *axis, int zeroch, int verbose);
 
 int recvPar(recPar *rec, float sub_x0, float sub_z0, float dx, float dz, int nx, int nz);
 
-int writesufile(char *filename, float *data, int n1, int n2, float f1, float f2, float d1, float d2);
-
-int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *src, shotPar *shot, bndPar *bnd, rayPar *ray, int verbose)
+int getParameters(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPar *ray, int verbose)
 {
 	int nx, nz, nsrc, ix, axis, is0;
 	int idzshot, idxshot;
@@ -42,7 +36,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
 	float dxrcv,dzrcv,dxspread,dzspread;
 	float xmax, zmax;
 	float xsrc1, xsrc2, zsrc1, zsrc2;
-	float *xsrca, *zsrca, rrcv;
+	float *xsrca, *zsrca;
 	float rsrc, oxsrc, ozsrc, dphisrc, ncsrc;
 	size_t nsamp;
 	int nxsrc, nzsrc;
@@ -277,11 +271,6 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
 	
 	recvPar(rec, sub_x0, sub_z0, dx, dz, nx, nz);
 
-	/* receivers are on a circle, use default interpolation to real (not on a grid-point) receiver position */
-	if (getparfloat("rrcv", &rrcv)) { 
-		if (!getparint("rec_int_p",&rec->int_p)) rec->int_p=3;
-	}
-
 	if (verbose) {
 		if (rec->n) {
 			dxrcv = rec->xr[MIN(1,rec->n-1)]-rec->xr[0];

raytime/raytime.c

@@ -29,11 +29,11 @@ void name_ext(char *filename, char *extension);
 
 void threadAffinity(void);
 
-int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *src, shotPar *shot, bndPar *bnd, rayPar *ray, int verbose);
+int getParameters(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPar *ray, int verbose);
 
 int getWaveParameter(float *slowness, icoord size, float dgrid, fcoord s, fcoord r, rayPar ray, fcoord *T, float *Jr);
 
-int readModel(modPar mod, bndPar bnd, float *velocity, float *slowness);
+int readModel(modPar mod, float *velocity, float *slowness);
 
 int defineSource(wavPar wav, srcPar src, modPar mod, float **src_nwav, int reverse, int verbose);
 
@@ -72,15 +72,11 @@ char *sdoc[] = {
 "   xsrca= ............ defines source array x-positions",
 "   zsrca= ............ defines source array z-positions",
 "   wav_random=1 ...... 1 generates (band limited by fmax) noise signatures ",
-"   src_multiwav=0 .... use traces in file_src as areal source",
 "   src_at_rcv=1 ...... inject wavefield at receiver coordinates (1), inject at source (0)",
 "" ,
 " PLANE WAVE SOURCE DEFINITION:",
 "   plane_wave=0 ...... model plane wave with nsrc= sources",
 "   nsrc=1 ............ number of sources per (plane-wave) shot ",
-"   src_angle=0 ....... angle of plane source array",
-"   src_velo=1500 ..... velocity to use in src_angle definition",
-"   src_window=0 ...... length of taper at edges of source array",
 "",
 " RANDOM SOURCE DEFINITION FOR SEISMIC INTERFEROMTERY:",
 "   src_random=0 ...... 1 enables nsrc random sources positions in one modeling",
@@ -89,11 +85,6 @@ char *sdoc[] = {
 "   xsrc2=0 ........... right bound for x-position of sources",
 "   zsrc1=0 ........... left bound for z-position of sources",
 "   zsrc2=0 ........... right bound for z-position of sources",
-"   tsrc1=0.0 ......... begin time interval for random sources being triggered",
-"   tsrc2=tmod ........ end time interval for random sources being triggered",
-"   tactive=tsrc2 ..... end time for random sources being active",
-"   tlength=tsrc2-tsrc1 average duration of random source signal",
-"   length_random=1 ... duration of source is rand*tlength",
 "   amplitude=0 ....... distribution of source amplitudes",
 "   distribution=0 .... random function for amplitude and tlength 0=flat 1=Gaussian ",
 "   seed=10 ........... seed for start of random sequence ",
@@ -126,10 +117,7 @@ int main(int argc, char **argv)
 {
 	modPar mod;
 	recPar rec;
-	snaPar sna;
-	wavPar wav;
 	srcPar src;
-	bndPar bnd;
 	shotPar shot;
 	rayPar ray;
     float *velocity, *slowness;
@@ -151,7 +139,7 @@ int main(int argc, char **argv)
 	requestdoc(0);
 
 	if (!getparint("verbose",&verbose)) verbose=0;
-	getParameters(&mod, &rec, &sna, &wav, &src, &shot, &bnd, &ray, verbose);
+	getParameters(&mod, &rec, &src, &shot, &ray, verbose);
 
 	/* allocate arrays for model parameters: the different schemes use different arrays */
 
@@ -162,7 +150,7 @@ int main(int argc, char **argv)
 
 	/* read velocity and density files */
 
-	readModel(mod, bnd, velocity, slowness);
+	readModel(mod, velocity, slowness);
 
 	/* read and/or define source wavelet(s) */
 
@@ -186,7 +174,6 @@ int main(int argc, char **argv)
 	/* Sinking source and receiver arrays: 
 	   If P-velocity==0 the source and receiver 
 	   postions are placed deeper until the P-velocity changes. 
-	   The free-surface position is stored in bnd.surface[ix].
 	   Setting the option rec.sinkvel only sinks the receiver position 
        (not the source) and uses the velocity 
 	   of the first receiver to sink through to the next layer. */
@@ -256,6 +243,9 @@ int main(int argc, char **argv)
         	grid.z = mod.nz;
         	grid.y = 1;
 
+#pragma omp parallel for default(shared) \
+private (coordgx) \
+private (irec,Time,Jr) 
         	for (irec=0; irec<rec.n; irec++) {
             	coordgx.x=mod.x0+rec.xr[irec];
             	coordgx.z=mod.z0+rec.zr[irec];

raytime/raytime.h

@@ -2,33 +2,15 @@
 #include<stdio.h>
 #include<math.h>
 
-typedef struct _compType { /* Receiver Type */
-	int vz;
-	int vx;
-	int p;
-	int txx;
-	int tzz;
-	int txz;
-	int pp;
-	int ss;
-	int ud;
-} compType;
-
 typedef struct _receiverPar { /* Receiver Parameters */
 	char *file_rcv;
-	compType type;
 	int n;
 	int nt;
-	int delay;
-	int skipdt;
 	int max_nrec;
 	int *z;
 	int *x;
 	float *zr;
 	float *xr;
-	int int_p;
-	int int_vx;
-	int int_vz;
 	int scale;
 	int sinkdepth;
 	int sinkvel;
@@ -36,79 +18,19 @@ typedef struct _receiverPar { /* Receiver Parameters */
 	float rho;
 } recPar;
 
-typedef struct _snapshotPar { /* Snapshot Parameters */
-	char *file_snap;
-	char *file_beam;
-	compType type;
-	int nsnap;
-	int delay;
-	int skipdt;
-	int skipdz;
-	int skipdx;
-	int nz;
-	int nx;
-	int z1;
-	int z2;
-	int x1;
-	int x2;
-	int vxvztime;
-	int beam;
-	int withbnd;
-} snaPar;
-
 typedef struct _modelPar { /* Model Parameters */
-	int iorder;
-	int ischeme;
-	int grid_dir;
 	int sh;
 	char *file_cp;
-	char *file_ro;
-	char *file_cs;
-	char *file_qp;
-	char *file_qs;
 	float dz;
 	float dx;
 	float dt;
-	float tmod;
-	int nt;
 	float z0;
 	float x0;
 	/* medium max/min values */
 	float cp_min;
 	float cp_max;
-	float cs_min;
-	float cs_max;
-	float ro_min;
-	float ro_max;
 	int nz;
 	int nx;
-	int naz;
-	int nax;
-	/* Vx: rox */
-	int ioXx;
-	int ioXz;
-	int ieXx;
-	int ieXz;
-	/* Vz: roz */
-	int ioZx;
-	int ioZz;
-	int ieZx;
-	int ieZz;
-	/* P, Txx, Tzz: lam, l2m */
-	int ioPx;
-	int ioPz;
-	int iePx;
-	int iePz;
-	/* Txz: muu */
-	int ioTx;
-	int ioTz;
-	int ieTx;
-	int ieTz;
-	/* attenuation / dissipative medium */
-	float Qp;
-	float Qs;
-	float fw;
-	float qr;
 } modPar;
 
 typedef struct _waveletPar { /* Wavelet Parameters */
@@ -137,8 +59,6 @@ typedef struct _sourcePar { /* Source Array Parameters */
 	int circle;
 	int array;
 	int random;
-	float *tbeg;
-	float *tend;
 	int multiwav;
 	float angle;
 	float velo;
@@ -158,29 +78,6 @@ typedef struct _shotPar { /* Shot Parameters */
 	int *x;
 } shotPar;
 
-typedef struct _boundPar { /* Boundary Parameters */
-	int top;
-	int bot;
-	int lef;
-	int rig;
-	float *tapz;
-	float *tapx;
-	float *tapxz;
-	int cfree;
-	int ntap;
-	int *surface;
-    int npml;
-    float R; /* reflection at side of model */
-    float m; /* scaling order */
-    float *pml_Vx;
-    float *pml_nzVx;
-    float *pml_nxVz;
-    float *pml_nzVz;
-    float *pml_nxP;
-    float *pml_nzP;
-
-} bndPar;
-
 typedef struct _raypar { /* ray-tracing parameters */
     int smoothwindow;
     int useT2;

raytime/readModel.c

@@ -26,7 +26,7 @@
 **/
 
 
-int readModel(modPar mod, bndPar bnd, float *velocity, float *slowness)
+int readModel(modPar mod, float *velocity, float *slowness)
 {
     FILE    *fpcp;
     size_t  nread;

raytime/writeSrcRecPos.c

@@ -126,12 +126,7 @@ int writeSrcRecPos(modPar *mod, recPar *rec, srcPar *src, shotPar *shot)
 		dum[rec->x[is]*nz+MIN(nz-1,rec->z[is]+1)] = -1.0;
 
 //		vmess("receiver position %d at grid[ix=%d, iz=%d] = (x=%f z=%f)", ir, ix+ioPx, rec.z[ir]+ioPz, rec.xr[ir]+mod.x0, rec.zr[ir]+mod.z0);
-		if (rec->int_vx==3) {
-			fprintf(fp, "(%f, %f)\n", rec->xr[is]*dx+sub_x0, rec->zr[is]*dz+sub_z0);
-		}
-		else {
-			fprintf(fp, "(%f, %f)\n", rec->x[is]*dx+sub_x0, rec->z[is]*dz+sub_z0);
-		}
+		fprintf(fp, "(%f, %f)\n", rec->x[is]*dx+sub_x0, rec->z[is]*dz+sub_z0);
 	}
 	fclose(fp);
 	writesufile("SrcRecPositions.su", dum, nz, nx, sub_z0, sub_x0, dz, dx);