From 98a6761d43ded555c2659dbec0badd4447acbc95 Mon Sep 17 00:00:00 2001
From: Jan Thorbecke <janth@xs4all.nl>
Date: Thu, 2 Nov 2017 17:04:44 +0100
Subject: [PATCH] added OpenMP + more cleaning
---
README | 6 ++-
raytime/Makefile | 4 +-
raytime/getParameters.c | 15 +-----
raytime/raytime.c | 24 +++------
raytime/raytime.h | 103 ---------------------------------------
raytime/readModel.c | 2 +-
raytime/writeSrcRecPos.c | 7 +--
7 files changed, 17 insertions(+), 144 deletions(-)
diff --git a/README b/README
index 904c6d4..cbeda68 100644
--- a/README
+++ b/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:
diff --git a/raytime/Makefile b/raytime/Makefile
index 0d70c1b..fdf35dd 100644
--- a/raytime/Makefile
+++ b/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
diff --git a/raytime/getParameters.c b/raytime/getParameters.c
index 1651e1a..4970036 100644
--- a/raytime/getParameters.c
+++ b/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];
diff --git a/raytime/raytime.c b/raytime/raytime.c
index 0626201..d38eaf4 100644
--- a/raytime/raytime.c
+++ b/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];
diff --git a/raytime/raytime.h b/raytime/raytime.h
index 4645539..045e99e 100644
--- a/raytime/raytime.h
+++ b/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;
diff --git a/raytime/readModel.c b/raytime/readModel.c
index e54c6d1..bed640a 100644
--- a/raytime/readModel.c
+++ b/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;
diff --git a/raytime/writeSrcRecPos.c b/raytime/writeSrcRecPos.c
index 1fb1530..faf2978 100644
--- a/raytime/writeSrcRecPos.c
+++ b/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);
--
GitLab