From 0583597e7c5dfe76cf7f5ebf5aaa3305d8726eea Mon Sep 17 00:00:00 2001 From: JanThorbecke <janth@xs4all.nl> Date: Mon, 17 Aug 2020 14:33:57 +0200 Subject: [PATCH] preparing for primaries manuscript R1 --- INSTALL | 11 +++++ MDD/deconvolve.c | 2 +- Make_include_template | 75 ++++++++++++--------------------- README.md | 70 +++++++++--------------------- REPRODUCE | 52 +++++++++++++++++++++++ marchenko/marchenko_primaries.c | 14 ++++-- marchenko/synthesis.c | 10 +++-- 7 files changed, 127 insertions(+), 107 deletions(-) create mode 100644 REPRODUCE diff --git a/INSTALL b/INSTALL index 5383758..bc66f15 100644 --- a/INSTALL +++ b/INSTALL @@ -11,6 +11,9 @@ cp Make_include_template Make_include 2) Check the compiler and CFLAGS options in the file Make_include and adapt to the system you are using. The default options are set for a the GNU C-compiler on a Linux system. A Fortran compiler is only needed for the MDD package. A C++ compiler is not needed to compile the code. The Makefile has been tested with GNU make. If the fortran compiler (FC=) is not set the MDD package will not be build. +The Make_include file contains example settings for GNU, Intel, AMD, PGI and Cray compilers. If you want to use one of these compilers don't +forget to comment the compilers you do not want to use. + 3) If the compiler options are set in the Make_include file you can type ``` make clean @@ -59,8 +62,13 @@ Usually the MKL libraries are installed in $MKLROOT. If that variable is not set find /opt/intel -name libmkl_gf_lp64.so and adjust MKLROOT in Make_include accordingly. + If MKLROOT is not set in your environment and also not set in Make_include then MKL will not be used and default options for BLAS, LAPACK, and FFT are chosen. +The Finite Different based RTM code (FDACRTMC) will not be installed if MKL is not available. This code depends of FFTW and we do want to +include the full FFTW source code package to this. Our aim is to keep the compilation of the code in the GitHub repository +as simple as possible with the smallest number of external libraries. + SEISMIC UNIX ----------- @@ -94,6 +102,7 @@ https://www.netlib.org/lapack/index.html FDACRTMC -------- fdacrtmc uses FFTW and the wisdom computations are stored on disk for re-usage. This directory is defined in fdacrtmc.h +MKL includes an FFTW interface and in presence of MKL (and MKLROOT defined) the FDACRTMC package will also get compiled. #ifndef WISDOMDIR #define WISDOMDIR "/tmp/fftw/" @@ -114,6 +123,7 @@ Other make commands which can be useful: make clean : removes all object files, but leaves libraries and executables make realclean: removes also object files, libraries and executables. + COMPILATION PROBLEMS -------------------- If you encounter missing trigometric / mathematical functions during compilation, for example; @@ -124,6 +134,7 @@ add '-lm -lc' around line 109 in Make_include: LIBS= -L$L -lgenfft $(BLAS) -lm -lc + UPDATES AND LATEST VERSION -------------------------- The latest version of the source code and manual can be found at: diff --git a/MDD/deconvolve.c b/MDD/deconvolve.c index 3e677c6..9d7d81e 100644 --- a/MDD/deconvolve.c +++ b/MDD/deconvolve.c @@ -97,7 +97,7 @@ void Aprod2_(int *m, int *n, float *x, float *y, size_t *indw); void computeMatrixInverse(complex *matrix, int nxm, int rthm, float eps_a, float eps_r, float numacc, int eigenvalues, float *eigen, int iw, int verbose); -extern complex *cB; +complex *cB; int deconvolve(complex *cA, complex *cB, complex *cC, complex *oBB, int nfreq, int nblock, size_t nstationA, size_t nstationB, float eps_a, float eps_r, float numacc, int eigenvalues, float *eigen, int rthm, int mdd, int conjgA, int conjgB, int lsqr_iter, float lsqr_damp, int k_iter, float TCscl, int verbose) { diff --git a/Make_include_template b/Make_include_template index ab90618..9528e90 100644 --- a/Make_include_template +++ b/Make_include_template @@ -6,7 +6,7 @@ # -3- on Solaris system use RANLIB=ranlib which is defined below # the current directory (in vi ":r!pwd") -ROOT=/Users/jan/src/OpenSource +ROOT=/vardim/home/thorbcke/src/OpenSource ############################################################################# # Some convenient abbreviations @@ -23,15 +23,14 @@ CC = gcc FC = gfortran # Linux gcc version 4.x OPTC = -O3 -ffast-math -#to include parallelisation with OpenMP +##to include parallelisation with OpenMP OPTC += -fopenmp - # for better performing code you can try: -#OPTC = -O3 -fno-trapping-math -ffast-math -funroll-all-loops -fomit-frame-pointer -ftree-vectorizer-verbose=1 -# Linux gcc very old version 3.x -#OPTC = -O3 -ffast-math -funroll-all-loops +#OPTC = -O3 -fno-trapping-math -ffast-math -funroll-all-loops -mavx -fomit-frame-pointer -mfpmath=sse -ftree-vectorizer-verbose=1 +# Linux gcc version 3.x +#OPTC = -O3 -ffast-math -funroll-all-loops -mfpmath=sse -#enable double precision only for use in FFTlib and emmod +#for double precision use FFTlib and emmod #OPTC += -DDOUBLE #OPTF += -fdefault-double-8 -fdefault-real-8 @@ -49,12 +48,9 @@ OPTC += -fopenmp ##OPTF = -O3 -no-prec-div -qopt-report-phase=vec,openmp #OPTC = -O3 -no-prec-div -xCORE-AVX2 #OPTF = -O3 -no-prec-div -xCORE-AVX2 -##to include parallelisation with OpenMP +###to include parallelisation with OpenMP #OPTC += -qopenmp -# Apple OSX intel 11.1.076 snow leopard 10.6.2 -#OPTC = -O3 -msse3 -no-prec-div -vec-report2 -fno-builtin-__sprintf_chk - #PGI #CC = pgcc #FC = pgf90 @@ -63,18 +59,12 @@ OPTC += -fopenmp #OPTF = -fast #LDFLAGS = -fast -Minfo=vect -Mvect=simd:256 -Msafeptr -#Pathscale -#CC = cc -#FC = ftn -#OPTC = -Ofast -OPT:Ofast -fno-math-errno -#OPTF = -Ofast -OPT:Ofast -fno-math-errno - #Apple OSX clang/gcc (10.9) llvm #CC = clang #OPTC = -Ofast #LDFLAGS = -Ofast -#AMD Open64 +#AMD Open64 this will be replaced with the latest AOCC compilers #CC = opencc #FC = openf95 #OPTC = -Ofast @@ -82,40 +72,27 @@ OPTC += -fopenmp #LDFLAGS = -static -Ofast ############################################################################# -# BLAS and LAPACK libraries and FFT LIBRARIES +# BLAS and LAPACK libraries #MKLROOT=/opt/intel/mkl/ -MKLLIB=${MKLROOT}/lib -OPTC += -DMKL -I$(MKLROOT)/include -I$(MKLROOT)/include/fftw -#for GNU compilers -#you might need to add intel64 to : ${MKLROOT}/lib/intel64 -BLAS = -L${MKLROOT}/lib/ -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl -#for GNU/clang on OSX -#BLAS = -Wl,-rpath ${MKLLIB} -L${MKLROOT}/lib/ -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl -#on linux you want to use groups and MKL is in lib/intel64 -#MKLLIB=${MKLROOT}/lib/intel64 -#BLAS = -Wl,-rpath ${MKLLIB} -Wl,--start-group -L${MKLLIB} -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -Wl,--end-group -lpthread -lm -ldl -#for intel compilers -#BLAS = -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 - -############################################################################# -# AMD ACML 4.4.0 -#AMDROOT = /home/thorbcke/amdsdk/v1.0/acml/open64_64 -#OPTC += -DACML440 -I$(AMDROOT)/include -#BLAS = -L$(AMDROOT)/lib -lacml -lfortran -lffio -lrt -lm - -# FOR ZFP LIBRARIES -ZFPROOT = ${ROOT}/zfp -ZFP = -L${ZFPROOT}/lib -lzfp -ZFPINCL = ${ZFPROOT}/include -OPTC += -I${ZFPINCL} +MKLROOT= +ifneq ($(strip $(MKLROOT)),) + ifeq ($(CC),icc) + #for intel compilers + OPTC += -DMKL -I$(MKLROOT)/include -I$(MKLROOT)/include/fftw + BLAS = -mkl + else + #for GNU and other compilers + #on linux you want to use groups and MKL is in lib/intel64 + MKLLIB=${MKLROOT}/lib/intel64 + OPTC += -DMKL -I$(MKLROOT)/include -I$(MKLROOT)/include/fftw + BLAS = -Wl,-rpath ${MKLLIB} -Wl,--start-group -L${MKLLIB} -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -Wl,--end-group -lpthread -lm -ldl +endif +else + BLAS = -llapack -lblas +endif #LIBARIES -LIBS= -L$L -lgenfft $(FFT) $(BLAS) $(ZFP) - +LIBS= -L$L -lgenfft $(BLAS) ######################################################################## diff --git a/README.md b/README.md index 28e547d..e436901 100644 --- a/README.md +++ b/README.md @@ -39,56 +39,38 @@ Jan Thorbecke, Evert Slob, Joeri Brackenhoff, Joost van der Neut, and Kees Wapen 2017, Geophysics, Vol. 82, no. 6 (November-December); p. WB29--WB45, doi: 10.1190/GEO2017-0108.1 Download: https://janth.home.xs4all.nl/Publications/Articles/ThorbeckeGPY2017.pdf --3- When you are using the marchenko_primaries algorithm developed by Lele Zhang please refer to the following paper +-3- If you used the code to construct homogenoeus Green's functions, please refeer to this paper in your related publications: + +Virtual acoustics in inhomogeneous media with single-sided access: + Wapenaar, K., Brackenhoff, J., Thorbecke, J., van der Neut, J., Slob, E., and Verschuur, E., +2018, Scientific Reports, Vol. 8, 2497. +Download: http://homepage.tudelft.nl/t4n4v/4_Journals/Nature/SR_18.pdf + +-4- When you are using the marchenko_primaries algorithm developed by Lele Zhang please refer to the following paper Free-surface and internal multiple elimination in one step without adaptive subtraction Lele Zhang and Evert Slob 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: +-5- 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 ------------- +See the seperate INSTALL file. -1) To compile and link the code you first have to set the ROOT variable in the Make_include file which can be found in the directory where you have found this README. -You can use Make_include_template as a first start: - -``` -cp Make_include_template Make_include -``` - -2) Check the compiler and CFLAGS options in the file Make_include and adapt to the system you are using. The default options are set for a the GNU C-compiler on a Linux system. A Fortran or C++ compiler are not needed to compile the code. The Makefile has been tested with GNU make. -3) If the compiler options are set in the Make_include file you can type -``` -make clean -make -``` -and the Makefile will make: +REPRODUCING +---------- +Almost all Figures in the papers mentioned above can be reproduced with the sofwtare in this repository. Please see the file REPRODUCE for further instructions -- FFT library -- fdelmodc / 3D -- marchenko / 3D -- utils -The libraries will be placed in the lib/ directory and the executables in the bin/ directory. - -To use the executables don't forget to include the pathname in your PATH: - -``` -bash/sh: -export PATH='path_to_this_directory'/bin:$PATH: -csh: -setenv PATH 'path_to_this_directory'/bin:$PATH: -``` Finite Difference Modeling: FDELMODC ------------------------------------ @@ -111,7 +93,7 @@ To reproduce the Figures shown in the GEOPHYICS paper "Implementation of the Mar To reproduce the Figures shown in the Scientific Reports paper "Virtual acoustics in inhomogeneous media with single-sided access" the scripts in marchenko/demo/ScientificReports directory can be used. The README in this directory gives more instructions and guidelines. -To reproduce the Figures shown in the GEOPHYICS paper "Implementation of the Marchenko Multiple Elimination algorithm" the scripts in marchenko/demo/oneD directory can be used. The README_PRIMARIES in this directory gives more instructions and guidelines. +To reproduce the Figures shown in the GEOPHYICS paper "Implementation of the Marchenko Multiple Elimination algorithm" the scripts in marchenko/demo/mme directory can be used. The README_PRIMARIES in this directory gives more instructions and guidelines. MDD --- @@ -135,6 +117,8 @@ 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. +In case MKL is installed on your system there is no need to install the netlib packages mentioned for MDD. + SEISMIC UNIX ----------- @@ -172,25 +156,9 @@ CFLAGS += -DWISDOMDIR="/directory/that/exists" or you can change the name of WISDOMDIR in fdacrtmc.h - -MISC ----- -Other make commands which can be useful: - -make clean : removes all object files, but leaves libraries and executables -make realclean: removes also object files, libraries and executables. - -COMPILATION PROBLEMS --------------------- -If you encounter missing trigometric / mathematical functions during compilation, for example; - -defineSource.c:(.text+0x2356): undefined reference to sin getParameters.o: In function getParameters: - -add '-lm -lc' around line 109 in Make_include: - -LIBS= -L$L -lgenfft $(BLAS) -lm -lc - - +The Finite Different based RTM code (FDACRTMC) will not be installed if MKL is not available. This code depends on FFTW and we do want to +include the full FFTW source code package to this. Our aim is to keep the compilation of the code in the GitHub repository +as simple as possible with the smallest number of external libraries. UPDATES AND LATEST VERSION -------------------------- diff --git a/REPRODUCE b/REPRODUCE new file mode 100644 index 0000000..524471a --- /dev/null +++ b/REPRODUCE @@ -0,0 +1,52 @@ + +REFERENCES +--------- +-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 + +*** To reproduce the Figures shown in the GEOPHYICS manuscript "Finite-difference modeling experiments for seismic interferometry" the scripts in fdelmodc/FiguresPaper directory can be used. Please read the README in the FiguresPaper directory for more instructions and guidelines. + + +-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 +2017, Geophysics, Vol. 82, no. 6 (November-December); p. WB29--WB45, doi: 10.1190/GEO2017-0108.1 +Download: https://janth.home.xs4all.nl/Publications/Articles/ThorbeckeGPY2017.pdf + +*** To reproduce the Figures shown in the GEOPHYICS paper "Implementation of the Marchenko method" the scripts in marchenko/demo/oneD directory can be used. The README in this directory gives more instructions and guidelines. + + +-3- If you used the code to construct homogenoeus Green's functions, please refeer to this paper in your related publications: + +Virtual acoustics in inhomogeneous media with single-sided access: + Wapenaar, K., Brackenhoff, J., Thorbecke, J., van der Neut, J., Slob, E., and Verschuur, E., +2018, Scientific Reports, Vol. 8, 2497. +Download: http://homepage.tudelft.nl/t4n4v/4_Journals/Nature/SR_18.pdf + +*** To reproduce the Figures shown in the Scientific Reports paper "Virtual acoustics in inhomogeneous media with single-sided access" the scripts in marchenko/demo/ScientificReports directory can be used. The README in this directory gives more instructions and guidelines. + +-4- When you are using the marchenko_primaries algorithm developed by Lele Zhang please refer to the following paper + +Free-surface and internal multiple elimination in one step without adaptive subtraction +Lele Zhang and Evert Slob +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 + +*** To reproduce the Figures shown in the GEOPHYICS paper "Implementation of the Marchenko Multiple Elimination algorithm" the scripts in marchenko/demo/mme directory can be used. The README_PRIMARIES in this directory gives more instructions and guidelines. + + +-5- 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 + diff --git a/marchenko/marchenko_primaries.c b/marchenko/marchenko_primaries.c index 1eddba6..3f6f9dc 100644 --- a/marchenko/marchenko_primaries.c +++ b/marchenko/marchenko_primaries.c @@ -60,7 +60,7 @@ char *sdoc[] = { " ", " INTEGRATION ", " ishot=nshots/2 ........... shot number(s) to remove internal multiples ", -" file_tinv= ............... shot-record (from R) to remove internal multiples", +" file_tinv= ............... shot-record to remove internal multiples", " file_src= ................ optional source wavelet to convolve selected ishot(s)", " COMPUTATION", " tap=0 .................... lateral taper R_ishot(1), file_shot(2), or both(3)", @@ -187,6 +187,8 @@ int main (int argc, char **argv) if (reci && ntap) vwarn("tapering influences the reciprocal result"); +/* defines the smooth transition zone of the time-window */ + smooth = MIN(smooth, shift); if (smooth) { costaper = (float *)malloc(smooth*sizeof(float)); @@ -297,6 +299,7 @@ int main (int argc, char **argv) } /*================ Read focusing operator(s) ================*/ +/* this is an optional functionality, typically one uses one of the shots from R */ if (file_tinv != NULL) { /* M0 is named DD */ muteW = (int *)calloc(Nfoc*nxs,sizeof(int)); @@ -504,7 +507,7 @@ int main (int argc, char **argv) dxs = dx; } - /* define tapers to taper edges of acquisition */ + /* define optional tapers to taper edges of acquisition */ if (tap == 1 || tap == 3) { tapersy = (float *)malloc(nxs*sizeof(float)); for (j = 0; j < ntap; j++) @@ -619,6 +622,7 @@ int main (int argc, char **argv) } perc=(iend-istart)/100;if(!perc)perc=1; + /* for the plane-wave option write the constructed plane wave to disk */ if (plane_wave) { writeDataIter("SRCplane.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, NINT(src_angle)); } @@ -656,6 +660,7 @@ int main (int argc, char **argv) */ /*================ start loop over number of time-samples ================*/ +/* for each time sample in this loop the Marchenko equations are solved */ for (ii=istart; ii<iend; ii++) { @@ -669,6 +674,8 @@ int main (int argc, char **argv) t5 = wallclock_time(); memset(M0, 0, Nfoc*nxs*ntfft*sizeof(float)); memset(v1plus, 0, Nfoc*nxs*ntfft*sizeof(float)); + + /* M0 equation (3) in the Geophysics implementation paper */ /* once every 'niterskip' time-steps start from fresh M0 and do niter (~20) iterations */ if ( ((ii-istart)%niterskip==0) || (ii==istart) ) { niterrun=niter; @@ -733,12 +740,13 @@ int main (int argc, char **argv) /*================ number of Marchenko iterations ================*/ + /* niterrun is m in equation (1) in the Geophysics implementation paper */ for (iter=0; iter<niterrun; iter++) { t2 = wallclock_time(); /*================ construction of Mi(-t) = - \int R(x,t) Mi(t) ================*/ - + /* synthesis process is the compute kernel in equations (4) and (5) in the Geophysics implementation paper */ synthesisp(Refl, Fop, Mi, RMi, nx, nt, nacq, nts, dt, xsyn, Nfoc, xrcv, xsrc, xnx, fxse, fxsb, dxs, dxsrc, dx, ntfft, nw, nw_low, nw_high, mode, reci, nshots, ixpos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv, ixmask, verbose); diff --git a/marchenko/synthesis.c b/marchenko/synthesis.c index 7531c25..39a3150 100644 --- a/marchenko/synthesis.c +++ b/marchenko/synthesis.c @@ -41,10 +41,14 @@ float fxse, float fxsb, float dxs, float dxsrc, float dx, int nshots, int *ixpos int linearsearch(int *array, size_t N, int value); /*================ Convolution and Integration ================*/ -/* Refl has the full acquisition grid R(x_r, x_s) +/* The synthesis process is the compute kernel in equations (4) and (5) in the Geophysics implementation paper + * It computes M_{2m}(x'_0,x''_0,t,t_2) = \int_{t'=0}^{+\infty} \int_{\partial \setD_0} R(x'''_0, x'_0,t') + * & M_{2m-1}(x'''_0,x''_0,t-t',t_2) d x'''_0 dt' + * + * Refl has the full acquisition grid R(x_r, x_s) * Fop has the acquisition grid of the operator, ideally this should be equal to the acquisition grid of Refl, - * so all traces can be used to compute R*Fop. - * The output RNi has the traces in the grid of Fop, these are the x_s positions of R(x_r,x_s) */ + * so all traces can be used to compute R*Fop. + * The output RNi has the traces in the grid of Fop, these are the x_s positions of R(x_r,x_s) */ void synthesis(complex *Refl, complex *Fop, float *Top, float *RNi, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nfoc, float *xrcv, float *xsrc, int *xnx, float fxse, float fxsb, float dxs, float dxsrc, float dx, int ntfft, int -- GitLab