Skip to content
Snippets Groups Projects
Commit 9881d53d authored by JanThorbecke's avatar JanThorbecke
Browse files

small changes

parent 5961b4c7
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Make_include_template
+ 4
− 4
View file @ 9881d53d
......@@ -82,8 +82,7 @@ OPTC += -fopenmp
#LDFLAGS = -static -Ofast
#############################################################################
# BLAS and LAPACK libraries
# and FFT LIBRARIES
# BLAS and LAPACK libraries and FFT LIBRARIES
#MKLROOT=/opt/intel/mkl/
MKLLIB=${MKLROOT}/lib
OPTC += -DMKL -I$(MKLROOT)/include -I$(MKLROOT)/include/fftw
......@@ -98,7 +97,8 @@ BLAS = -L${MKLROOT}/lib/ -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm
#for intel compilers
#BLAS = -mkl
#
#If you do not want to use MKL
#If you do not want or can not use MKL you can try the basic lapack and blas libaries that are default installed on most Linux
#systems
#BLAS = -llapack -lblas
#############################################################################
......@@ -108,7 +108,7 @@ BLAS = -L${MKLROOT}/lib/ -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm
#BLAS = -L$(AMDROOT)/lib -lacml -lfortran -lffio -lrt -lm
#LIBARIES
LIBS= -L$L -lgenfft $(BLAS)
LIBS= -L$L -lgenfft $(BLAS)
########################################################################
......
README
+ 36
− 7
View file @ 9881d53d
......@@ -13,19 +13,20 @@ For more information, go to http://software.seg.org/2017/00XX .
You must read and accept usage terms at:
http://software.seg.org/disclaimer.txt before use.
DOI reference of release
1.4.0 https://zenodo.org/badge/latestdoi/23060862
REFERENCES
---------
-1- If the Finite Difference code has helped you in your research please refer to our paper in your publications:
-0- DOI reference of this software release
https://zenodo.org/badge/latestdoi/23060862
-1- If the Finite Difference code has helped you in your research please refer to the following paper in your publications:
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
Download: https://janth.home.xs4all.nl/Publications/Articles/ThorbeckeDraganov2012.pdf
-2- If the Machenko code has helped you in your research please refer to our paper in your publications:
-2- If the Machenko code has helped you in your research please refer to this paper in your publications:
Implementation of the Marchenko method
Jan Thorbecke, Evert Slob, Joeri Brackenhoff, Joost van der Neut, and Kees Wapenaar
......@@ -39,6 +40,15 @@ Lele Zhang1 and Evert Slob2
2019, Geophysics, Vol. 84, no. 1 (January-February); p. A7-A11, doi: 10.1190/GEO2018-0548.1
Download: http://homepage.tudelft.nl/t4n4v/BeyondInterferometry/geo_19h.pdf
-4- If you use the fdacrtmc code of Max Holicki please refer to the following paper:
Acoustic directional snapshot wavefield decomposition
Holicki, M., Drijkoningen, G., and Wapenaar, K., 2019, Acoustic directional snapshot wavefield decomposition: Geophysical
Prospecting, Vol. 67, 32-51
Download: http://homepage.tudelft.nl/t4n4v/4_Journals/Geophys.Prosp/GP_19a.pdf
-5-
INSTALLATION
-------------
......@@ -99,19 +109,22 @@ The MDD kernels depend on BLAS and LAPACK calls. Free downloads of these librari
https://www.netlib.org/blas/index.html
https://www.netlib.org/lapack/index.html
MKL libraries
-------------
If you are running on x64_86 processors you can download (for free) Intel's highly optimised MKL package:
https://software.intel.com/en-us/mkl/choose-download
Usually the MKL libraries are installed in $MKL_ROOT. If that variable is not set, you can try to find the correct path by searching for one of the libraries:
These Libraries include highly optimised libraries of BLAS, LAPACK, and FFT(W).
Usually the MKL libraries are installed in $MKLROOT. If that variable is not set, you can try to find the correct path by searching for one of the libraries:
find /opt/intel -name libmkl_gf_lp64.so
and adjust MKLROOT in Make_include accordingly.
You can also completely disable the use of MKL by commenting out the MKL parts in Make_include.
SEISMIC UNIX
-----------
If you want to use the .su files with SU from CWP:
......@@ -126,12 +139,28 @@ so the XDRFLAG is empty. If you have made that change you have to remake SU with
make remake
make xtremake
ZFP
---
The fdacrtmc code makes use of ZFP compression to store the snaphots in CPU memory. This package is included in this repository for your convenience. The latest package and detailed explanation can be found on:
The fdacrtmc and the 3D Marchenko code makes use of ZFP compression to store snaphots in CPU memory or do effient IO. This package is included in this repository for your convenience. The latest package and detailed explanation can be found on:
https://github.com/LLNL/zfp
FDACRTMC
--------
fdacrtmc uses FFTW and the wisdom computations are stored on disk for re-usage. This directory is defined in fdacrtmc.h
#ifndef WISDOMDIR
#define WISDOMDIR "/tmp/fftw/"
#endif
For the code to run properly the directory /tmp/fftw/ must exist. It is also possible to compile the fdacrtmc code that defines
WISDOMDIR. In fdacrtmc/Makefile add:
CFLAGS += -DWISDOMDIR="/directory/that/exists"
or you can change the name of WISDOMDIR in fdacrtmc.h
MISC
----
......@@ -147,7 +176,7 @@ The latest version of the source code and manual can be found at:
http://www.xs4all.nl/~janth/Software/Software.html
or at github:
or at github:
git clone https://github.com/JanThorbecke/OpenSource.git
git clone git://github.com/JanThorbecke/OpenSource.git
......
marchenko/Makefile
+ 1
− 0
View file @ 9881d53d
......@@ -42,6 +42,7 @@ SRCP = marchenko_primaries.c \
readData.c \
readShotData.c \
readTinvData.c \
findFirstBreak.c \
writeData.c \
writeDataIter.c \
wallclock_time.c \
......
marchenko/demo/oneD/primaries.scr
+ 3
− 3
View file @ 9881d53d
#!/bin/bash -x
#SBATCH -J marchenko_primaries
#SBATCH --cpus-per-task=40
#SBATCH --cpus-per-task=10
#SBATCH --ntasks=1
#SBATCH --time=2:40:00
......@@ -8,7 +8,7 @@
#cd $SLURM_SUBMIT_DIR
export PATH=$HOME/src/OpenSource/bin:$PATH:
export OMP_NUM_THREADS=40
export OMP_NUM_THREADS=10
#shot record to remove internal multiples
select=451
......@@ -17,7 +17,7 @@ makewave fp=20 dt=0.004 file_out=wave.su nt=2048 t0=0.0 scale=0 scfft=1
../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
nshots=901 verbose=2 istart=40 iend=500 fmax=90 pad=1024 \
niter=31 smooth=10 niterec=2 niterskip=1 shift=20 file_rr=pred_rr.su T=0 file_update=update.su
niter=31 smooth=10 niterec=2 niterskip=50 shift=20 file_rr=pred_rr.su T=0 file_update=update.su
exit;
......
marchenko/demo/oneD/primariesPlane.scr
+ 6
− 6
View file @ 9881d53d
#!/bin/bash -x
#SBATCH -J marchenko_primaries
#SBATCH --cpus-per-task=40
#SBATCH --cpus-per-task=10
#SBATCH --ntasks=1
#SBATCH --time=2:40:00
......@@ -8,20 +8,20 @@
#cd $SLURM_SUBMIT_DIR
export PATH=$HOME/src/OpenSource/bin:$PATH:
export OMP_NUM_THREADS=40
export OMP_NUM_THREADS=10
#shot record to remove internal multiples
select=451
makewave fp=20 dt=0.004 file_out=wave.su nt=1024 t0=0.0 scale=0 scfft=1
#../../marchenko_primaries file_shot=shotsdx5_rp.su plane_wave=1 file_src=wave.su src_angle=10 \
#../../marchenko_primaries file_shot=shotsdx5_rp.su plane_wave=1 file_src=wave.su src_angle=30 \
# nshots=901 verbose=2 istart=276 iend=277 fmax=90 pad=1024 \
# xorig=900 niter=2 smooth=10 niterskip=1 file_iter=iterplane.su shift=20 file_rr=plane2_10_rr.su T=0
# niter=2 smooth=10 niterskip=1 file_iter=iterplane.su shift=20 file_rr=plane2_10_rr.su T=0
#exit
../../marchenko_primaries file_shot=shotsdx5_rp.su plane_wave=1 file_src=wave.su src_angle=-10 \
../../marchenko_primaries file_shot=shotsdx5_rp.su plane_wave=1 file_src=wave.su src_angle=30.0 \
nshots=901 verbose=2 istart=40 iend=1000 fmax=90 pad=1024 \
niter=20 smooth=10 niterskip=40 niterec=0 shift=20 file_rr=plane2_10_rr.su T=0
niter=20 smooth=10 niterskip=60 niterec=2 shift=20 file_rr=plane2_10_rr.su T=0
exit
......
marchenko/findFirstBreak.c 0 → 100644
+ 109
− 0
View file @ 9881d53d
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#ifndef MAX
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#endif
#ifndef MIN
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#endif
#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
int findFirstBreak(float *shot, int nx, int nt, int ishot, float *maxval, int tr, int hw, int verbose)
{
int i, j;
int jmax, tstart, tend;
float xmax, tmax, lmax;
/* find consistent (one event) maximum related to maximum value */
/* find global maximum
xmax=0.0;
for (i = 0; i < nx; i++) {
tmax=0.0;
jmax = 0;
for (j = 0; j < nt; j++) {
lmax = fabs(shot[i*nt+j]);
if (lmax > tmax) {
jmax = j;
tmax = lmax;
if (lmax > xmax) {
ishot = i;
xmax=lmax;
}
}
}
maxval[i] = jmax;
}
*/
/* alternative find maximum at source position ishot */
tmax=0.0;
jmax = 0;
if (tr == 0) {
for (j = 0; j < nt; j++) {
lmax = fabs(shot[ishot*nt+j]);
if (lmax > tmax) {
jmax = j;
tmax = lmax;
if (lmax > xmax) {
xmax=lmax;
}
}
}
}
else {
for (j = nt-1; j >= 0; j--) {
lmax = fabs(shot[ishot*nt+j]);
if (lmax > tmax) {
jmax = j;
tmax = lmax;
if (lmax > xmax) {
xmax=lmax;
}
}
}
}
maxval[ishot] = jmax;
if (verbose >= 3) vmess("Mute max at src-trace %d is sample %d", ishot, maxval[ishot]);
/* search forward in trace direction from maximum in file */
for (i = ishot+1; i < nx; i++) {
tstart = MAX(0, (maxval[i-1]-hw));
tend = MIN(nt-1, (maxval[i-1]+hw));
jmax=tstart;
tmax=0.0;
for(j = tstart; j <= tend; j++) {
lmax = fabs(shot[i*nt+j]);
if (lmax > tmax) {
jmax = j;
tmax = lmax;
}
}
maxval[i] = jmax;
}
/* search backward in trace direction from maximum in file */
for (i = ishot-1; i >=0; i--) {
tstart = MAX(0, (maxval[i+1]-hw));
tend = MIN(nt-1, (maxval[i+1]+hw));
jmax=tstart;
tmax=0.0;
for(j = tstart; j <= tend; j++) {
lmax = fabs(shot[i*nt+j]);
if (lmax > tmax) {
jmax = j;
tmax = lmax;
}
}
maxval[i] = jmax;
}
if (tr != 0) {
for (i = 0; i < nx; i++) maxval[i] = nt-1 - maxval[i];
}
return 0;
}
marchenko/marchenko_primaries.c
+ 45
− 15
View file @ 9881d53d
......@@ -29,6 +29,7 @@ int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx
int readTinvData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nfoc, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose);
int findFirstBreak(float *shot, int nx, int nt, int ishot, float *maxval, int tr, int hw, int verbose);
int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int t0shift, int iter);
int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm);
int readData(FILE *fp, float *data, segy *hdrs, int n1);
......@@ -110,7 +111,7 @@ int main (int argc, char **argv)
FILE *fp_out, *fp_rr, *fp_w, *fp_up;
size_t nread;
int i, j, k, l, ret, nshots, Nfoc, nt, nx, nts, nxs, ngath;
int size, n1, n2, ntap, tap, di, ntraces;
int size, n1, n2, ntap, tap, di, ntraces, tr;
int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
int reci, countmin, mode, ixa, ixb, n2out, verbose, ntfft;
int iter, niter, niterec, recur, niterskip, niterrun, tracf, *muteW;
......@@ -121,9 +122,9 @@ int main (int argc, char **argv)
double t0, t1, t2, t3, t4, tsyn, tread, tfft, tcopy, tii;
double energyMi, *energyM0;
float tt0, d1, d2, f1, f2, fxsb, fxse, ft, fx, *xsyn, dxsrc;
float *M0, *DD, *RR, *SRC, dt, dx, dxs, scl, mem;
float *M0, *DD, *RR, *SRC, dt, dx, dxs, scl, mem, scltap;
float *rtrace, *tmpdata, *f1min, *f1plus, *RMi, *Mi, *trace;
float *Mup, *Msp;
float *Mup, *Msp, *maxval;
float xmin, xmax, scale, tsq;
float Q, f0, *ixmask, *costaper;
float src_velo, src_angle, grad2rad, p, *twplane;
......@@ -228,20 +229,20 @@ int main (int argc, char **argv)
//if(!getparint("xorig", &xorig)) xorig=-(nxs-1)/2;
/* compute time delay for plane-wave responses */
twplane = (float *) calloc(nxs,sizeof(float));
twplane = (float *) calloc(nxs,sizeof(float)); /* initialize with zeros */
if (plane_wave==1) {
grad2rad = 17.453292e-3;
p = sin(src_angle*grad2rad)/src_velo;
if (p < 0.0) {
for (i=0; i<nxs; i++) {
twplane[i] = fabsf((nxs-i-1)*dxs*p)+tt0;
if (verbose >=3) vmess("plane-wave delay-time i=%d x=%f t=%f", i, dxs*(nxs-i-1), twplane[i]);
}
}
else {
for (i=0; i<nxs; i++) {
twplane[i] = (i)*dxs*p+tt0;
//twplane[i] = dxs*(xorig+i)*tan(src_angle*grad2rad)/src_velo;
// fprintf(stderr,"plane-wave i=%d x=%f t=%f %f\n", i, dxs*(xorig+i), twplane[i], (xorig+i)*dxs*p);
if (verbose >=3) vmess("plane-wave delay-time i=%d x=%f t=%f", i, dxs*(i), twplane[i]);
}
}
}
......@@ -367,9 +368,9 @@ int main (int argc, char **argv)
if (tap == 2 || tap == 3) {
tapersh = (float *)malloc(nx*sizeof(float));
for (j = 0; j < ntap; j++)
tapersh[j] = (cos(PI*(j-ntap)/ntap)+1)/2.0;
for (j = ntap; j < nx-ntap; j++)
for (j = 1; j < ntap; j++)
tapersh[j-1] = (cos(PI*(j-ntap)/ntap)+1)/2.0;
for (j = ntap-1; j < nx-ntap; j++)
tapersh[j] = 1.0;
for (j = nx-ntap; j < nx; j++)
tapersh[j] =(cos(PI*(j-(nx-ntap))/ntap)+1)/2.0;
......@@ -382,7 +383,6 @@ int main (int argc, char **argv)
}
}
}
free(tapersh);
}
/* check consistency of header values */
......@@ -413,6 +413,7 @@ int main (int argc, char **argv)
scl = 1.0/((float)2.0*ntfft);
rtrace = (float *)calloc(ntfft,sizeof(float));
ctrace = (complex *)calloc(nfreq+1,sizeof(complex));
for (i = 0; i < xnx[ishot]; i++) {
for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
ctrace[j].r = Refl[ishot*nw*nx+m*nx+i].r*cwave[j].r + Refl[ishot*nw*nx+m*nx+i].i*cwave[j].i;
......@@ -423,9 +424,22 @@ int main (int argc, char **argv)
for (j = 0; j < nts; j++) {
DD[0*nxs*nts+i*nts+j] = -1.0*scl*rtrace[j];
}
/* remove taper at edges for selected shot record */
//if (tap == 2 || tap == 3) scltap= scl/tapersh[i];
//for (j = 0; j < nts; j++) {
// DD[0*nxs*nts+i*nts+j] = -1.0*scltap*rtrace[j];
//}
}
/* construct plane wave from all shot records */
/* set timereversal for searching First Break */
/* for experimenting with non flat truncation windows */
//maxval = (float *)calloc(nxs,sizeof(float));
//tr = 1;
//findFirstBreak(DD, xnx[ishot], nts, ishot, maxval, tr, hw, verbose);
//twplane = (float *) calloc(nxs,sizeof(float));
//for (i=0; i<nxs; i++) twplane[i] = dt*(maxval[i]-maxval[ishot]);
//for (i = 0; i < xnx[ishot]; i++) fprintf(stderr,"maxval[%d] = %f tw=%f\n", i, dt*maxval[i], twplane[i]);
/* construct plane wave (time reversed and multiplid with -1) from all shot records */
if (plane_wave==1) {
for (l=0; l<nshots; l++) {
memset(ctrace, 0, sizeof(complex)*(nfreq+1));
......@@ -451,6 +465,7 @@ int main (int argc, char **argv)
tom = j*deltom*twplane[l];
csum.r = cos(-tom);
csum.i = sin(-tom);
/* Optional add wavelet to SRC-field */
//cwav.r = csum.r*cwave[j].r - csum.i*cwave[j].i;
//cwav.i = csum.i*cwave[j].r + csum.r*cwave[j].i;
ctrace[j] = csum;
......@@ -523,8 +538,21 @@ int main (int argc, char **argv)
vmess("number of time samples fft nt nts = %d %d %d", ntfft, nt, nts);
vmess("time sampling = %e ", dt);
vmess("smoothing taper for time-window = %d ", smooth);
if (plane_wave) {
vmess("*** Plane-wave processing selected *** ");
vmess("plane-wave angle = %f ", src_angle);
vmess("plane-wave velocity = %f ", src_velo);
vmess("plane-wave t0-shift = %f ", tt0);
vmess("plane-wave DD input vector = DDplane%03d.su ", src_angle);
vmess("plane-wave SRC for migration = SRCplane%03d.su ", src_angle);
}
if (file_rr != NULL) vmess("RR output file = %s ", file_rr);
if (file_iter != NULL) vmess("Iterations output file = %s ", file_iter);
if (file_iter != NULL) {
vmess("Iterations output file = %s ", file_iter);
vmess("Initialisation input data = M0_%06d.su ", 1000*istart);
vmess("f1min intermediate array = f1min_%03d'iter'.su ", istart);
vmess("Mi intermediate array = Mi_%03d'iter'.su ", istart);
}
}
/*================ initializations ================*/
......@@ -813,8 +841,10 @@ int main (int argc, char **argv)
for (i = 0; i < npos; i++) {
ix = ixpos[i];
iw = NINT((ii*dt+twplane[ix])/dt);
RR[l*nxs*nts+i*nts+iw] = f1min[l*nxs*nts+i*nts+iw];
if (file_update != NULL) Msp[l*nxs*nts+i*nts+iw] = Mup[l*nxs*nts+i*nts+iw];
if ( iw<nts && iw>=0 ) {
RR[l*nxs*nts+i*nts+iw] = f1min[l*nxs*nts+i*nts+iw];
if (file_update != NULL) Msp[l*nxs*nts+i*nts+iw] = Mup[l*nxs*nts+i*nts+iw];
}
}
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment