Merge branch 'master' of https://github.com/JanThorbecke/OpenSource

FFTlib/Makefile

@@ -44,7 +44,7 @@ $(LIB)  : $(ARCH)
 mkdirs:
 	-mkdir -p $(ROOT)/lib
 	-mkdir -p $(ROOT)/include
-	-cp -rp genfft.h ../include/genfft.h
+	-cp -rp genfft.h $(ROOT)/include/genfft.h
 
 bins: 
 	cd test ; $(MAKE)

Make_include_swan

+# Makefile for general rules
+
+# To Change the compile environment to your current system you should set:
+#   -1- ROOT variable to the directory where you found this file
+#   -2- if needed use a different compiler (CC) if gcc is not available
+#   -3- on Solaris system use RANLIB=ranlib which is defined below
+
+# the current directory (in vi ":r!pwd")
+ROOT=/lus/scratch/${USER}/OpenSource
+
+#############################################################################
+# Some convenient abbreviations
+B = $(ROOT)/bin
+I = $(ROOT)/include
+L = $(ROOT)/lib
+
+########################################################################
+# C compiler; change this only if you are using a different C-compiler
+
+#GNU 
+#CC = gcc-mp-5 
+#CC = gcc
+#FC = gfortran
+# Linux gcc version 4.x
+#OPTC = -O3 -ffast-math
+#to include parallelisation with OpenMP
+#OPTC += -fopenmp
+
+# for better performing code you can try:
+#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
+
+#for double precision use FFTlib and emmod
+#OPTC += -DDOUBLE
+#OPTF += -fdefault-double-8 -fdefault-real-8
+
+#Cray
+#CC=cc
+#FC=ftn
+#OPTC = -O2
+#OPTF = -O2
+
+#Intel 
+CC = icc
+FC = ifort
+### Linux
+##OPTC = -O3 -no-prec-div -qopt-report-phase=vec,openmp
+##OPTF = -O3 -no-prec-div -qopt-report-phase=vec,openmp
+OPTC = -O3 -no-prec-div -xSANDYBRIDGE 
+#OPTF = -O3 -no-prec-div -xCORE-AVX2
+##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
+#OPTC = -fast  -Minfo=vect -Mvect=simd:256 -Msafeptr
+#OPTC = -fast  -Minfo=vect -Mvect=simd:256 -Msafeptr -Mprof=lines
+#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
+#CC = opencc
+#FC = openf95
+#OPTC = -Ofast
+#OPTF = -Ofast
+#LDFLAGS = -static -Ofast
+
+#############################################################################
+# BLAS and LAPACK libraries 
+MKLLIB=${MKLROOT}/lib
+#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 on OSX
+#BLAS = -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
+
+#############################################################################
+# FOR FFT LIBRARIES
+#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
+#Intel MKL
+#MKLROOT=/opt/intel/mkl/
+MKLLIB=${MKLROOT}/lib
+OPTC += -DMKL -I$(MKLROOT)/include
+#for GNU compilers
+#FFT = -Wl,-rpath ${MKLLIB} -L${MKLLIB} -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
+#for GNU on OSX
+#FFT = -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
+#FFT = -Wl,-rpath ${MKLLIB} -Wl,--start-group -L${MKLLIB} -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -Wl,--end-group -lpthread -lm -ldl
+#for Intel compilers
+FFT = -Wl,-rpath ${MKLLIB} -L${MKLLIB} -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
+
+#LIBARIES
+LIBS= -L$L -lgenfft $(FFT) $(BLAS)
+
+
+########################################################################
+# standard CFLAGS
+CFLAGS = -I$I -I.
+
+#############################################################################
+# the archiver
+AR = ar
+
+#############################################################################
+# ar FLAGS
+ARFLAGS = rv
+
+#############################################################################
+# ranlib definition
+RANLIB = ar -s
+# on Sun SOLARIS use:
+#RANLIB = ranlib
+
+
+.SUFFIXES : .o .c .cc .f .a .F90
+.c.o	:
+	$(CC) -c $(CFLAGS) $(OPTC) $<
+.c.a	:
+	$(CC) -c $(CFLAGS) $(OPTC) $<
+	$(AR) $(ARFLAGS) $@ $*.o
+	rm -f $*.o
+.o.a	:
+	$(AR) $(ARFLAGS) $@ $*.o
+	rm -f $*.o
+.f.o	:
+	$(FC) -c $(FFLAGS) $(OPTF) $<
+.F90.o	:
+	$(FC) -c $(FFLAGS) $(OPTF) $<
+.f.a	:
+	$(FC) -c $(FFLAGS) -I$I $<
+	$(AR) $(ARFLAGS) $@ $*.o
+	rm -f $*.o
+.cc.a	:
+	$(C++) -c $(C++FLAGS) -I$I $<
+	$(AR) $(ARFLAGS) $@ $*.o
+	rm -f $*.o
+.cc.o	:
+	$(C++) -c $(C++FLAGS) $<
+

Makefile

@@ -4,7 +4,7 @@ all: mkdirs
 	cd FFTlib		; $(MAKE)
 	cd fdelmodc		; $(MAKE) install
 	cd utils		; $(MAKE) install
-	cd marchenko	; $(MAKE) install
+	cd marchenko		; $(MAKE) install
 	cd corrvir		; $(MAKE) install
 	cd raytime		; $(MAKE) install
 	cd MDD			; $(MAKE) install
@@ -18,7 +18,7 @@ clean:
 	cd FFTlib 		; $(MAKE) $@
 	cd fdelmodc		; $(MAKE) $@
 	cd utils		; $(MAKE) $@
-	cd marchenko	; $(MAKE) $@
+	cd marchenko		; $(MAKE) $@
 	cd corrvir		; $(MAKE) $@
 	cd raytime		; $(MAKE) $@
 	cd MDD			; $(MAKE) $@
@@ -27,7 +27,7 @@ realclean:
 	cd FFTlib 		; $(MAKE) $@
 	cd fdelmodc		; $(MAKE) $@
 	cd utils		; $(MAKE) $@
-	cd marchenko	; $(MAKE) $@
+	cd marchenko		; $(MAKE) $@
 	cd corrvir		; $(MAKE) $@
 	cd raytime		; $(MAKE) $@
 	cd MDD			; $(MAKE) $@

README

@@ -65,10 +65,6 @@ export PATH='path_to_this_directory'/bin:$PATH:
 csh:
 setenv PATH 'path_to_this_directory'/bin:$PATH:
 
-*** FOR SU USERS:
-If you want to use the .su files with SU from CWP:
-==> Please make sure that SU is compiled without XDR (in $CWPROOT/Makefile.config make sure that XDRFLAG is NOT set). The SU output files of fdelmodc are all based on local IEEE data.
-*** 
 
 > make clean
 
@@ -100,10 +96,25 @@ 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
 
-If you are running on Intel processors you can download (for free) Intel's highly optimised MKL package:
+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:
+
+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 these parts in Make_include. 
+
+SEISMIC UNIX
+-----------
+If you want to use the .su files with SU from CWP:
+git clone https://github.com/JohnWStockwellJr/SeisUnix
+
+==> Please make sure that SU is compiled without XDR (in $CWPROOT/Makefile.config make sure that XDRFLAG is NOT set). The SU output files of fdelmodc are all based on local IEEE data.
 
 MISC
 ----

marchenko/demo/WS15/README.1

+
+*1* Source the setup.sh script that you can find on 
+
+source /home/users/jan/WS15setup.sh
+
+The setup.sh script completes the following tasks:
+
+-a-Set PATH to find Seismic Unix and Marchenko code 
+export CWPROOT=/home/users/jan/SeisUnix/
+export PATH=.:$CWPROOT/bin:/lus/scratch/$USER/OpenSource/bin:$PATH:
+
+-b- Create working directory 
+/lus/scratch/$USER
+
+-c- set up compilation environment with the Intel compiler
+
+-d- Get the Marchenko software (and other utilities modeling code))
+git clone https://github.com/JanThorbecke/OpenSource.git
+
+if the git clone takes too much time you can copy the source code:
+
+cp -rp /lus/snx11029/jan/OpenSource /lus/scratch/$USER/
+ 
+
+*2* after the set-up the code can be compiled 
+
+cd OpenSource
+cp Make_include_swan Make_include
+make clean
+make 
+
+After succesfull complication the code can be used for running the exercises and demo's
+
+
+*3* Running jobs 
+
+To run example jobs must be submitted to the Workload Manager PBSpro.
+An example script can be found on:
+
+/lus/scratch/$USER/OpenSource/WS15/job.pbs
+
+in job.pbs leave the number of nodes set to 1; the code is not MPI parallel and cannot use more than 1 node. The number of OMP_NUM_THREADS can be changed and is currently set to 40.
+Add the ...scr command (that can be found in the demo directories) as the last line in the aprun command line and submit the job to the queue with 'qsub job.pbs'
+you can check the status of your job in the queue with 'qstat'
+
+
+*4* Display of results:
+
+suximage < result.su perc=99 &
+
+if X11 display is too slow:
+
+supsimage < result.su perc=99 > results.eps
+convert results.eps ~/results.png
+
+copy png file to your laptop you are working on
+
+scp trxx@swan.cray.com:~/*png . 
+
+
+*5* Usefull job monitor commands specific for Cray
+
+apstat : status of applications running
+xtnodestat : placement of application on system
+xtprocadmin -A : details of installed nodes
+

marchenko/demo/WS15/README.2

+
+Marchenko demo/twoD: reproduce the figures from the paper, but with a 2D model.
+
+cd /lus/scratch/$USER/OpenSource/marchenko/demo/twoD
+
+# Copy pre-computed shots
+cp -rp /lus/scratch/jan/OpenSource/marchenko/demo/twoD/shots .
+cp /lus/scratch/$USER/OpenSource/marchenko/demo/WS15/job.pbs .
+
+Adapt job.pbs to 
+
+Compute the intial focusing operator placed at a depth of 1100 m by running 
+
+initialFocus.scr
+
+To compare the Green's function computed by marchenko generate the reference output
+
+adapt job.pbs to run referenceShot.scr
+
+run the marchenko program interactively (no need to submit job to queue);
+ 
+./marchenko.scr
+
+Compare the Green's function computed by marchenko with the reference (adapt job.pbs to run referenceShot.scr)
+
+The script eps.scr generate .eps files and .png files that can be copied back to your local laptop for display.
+
+You can experiment with 
+-amplitude of R: scale with a constant ampltiude 
+-investigate the wavelet that is used to model R
+-number of iterations, 
+-intial focusing position, 
+-investigate at intermediate results (see ../oneD/marchenkoIter.scr as an example).
+
+
+

marchenko/demo/WS15/README.3

+plane waves

marchenko/demo/WS15/README.4

+
+ 
+Primaries only

marchenko/demo/WS15/job.pbs

+#!/bin/bash 
+#PBS -N Marchenko
+#PBS -j oe
+#PBS -l place=scatter,select=1
+#PBS -S /bin/bash
+#PBS -V
+#PBS -q bw44-sm
+# PBS -q sk40-sm
+#PBS -l walltime=00:10:00
+
+set -x
+cd $PBS_O_WORKDIR
+ulimit -s unlimited
+
+export KMP_AFFINITY=disabled
+export OMP_NUM_THREADS=40
+
+starttime=`date +%s%N`
+
+aprun -n1 -d $OMP_NUM_THREADS 
+
+endtime=`date +%s%N`
+runtime=$(echo "scale=9; 1.0*10^(-9)*(${endtime}-${starttime})" | bc -l)
+echo "Runtime = $runtime seconds" 
+

marchenko/demo/WS15/setup.sh

+
+export CWPROOT=/home/users/jan/SeisUnix/
+export PATH=.:$CWPROOT/bin:/lus/scratch/$USER/OpenSource/bin:$PATH:
+
+module swap PrgEnv-cray PrgEnv-intel
+module list
+
+mkdir -p /lus/scratch/$USER
+
+cd /lus/scratch/$USER
+#git clone https://github.com/JanThorbecke/OpenSource.git
+

marchenko/demo/invisible/model.scr

@@ -25,7 +25,7 @@ makemod sizex=6000 sizez=1500 dx=$dx dz=$dx cp0=1000  ro0=1000 \
 #define wavelet for modeling R
 makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=4096 t0=0.3 scale=0 scfft=1
 
-export OMP_NUM_THREADS=8
+export OMP_NUM_THREADS=12
 
 #Model shot record in middle of model
 fdelmodc \

marchenko/demo/oneD/marchenko.scr

@@ -2,10 +2,13 @@
 
 export PATH=$HOME/src/OpenSource/bin:$PATH:
 export OMP_NUM_THREADS=1
+export MKL_NUM_THREADS=1
 
 #mute all events below the first arrival to get the intial focusing field
 fmute file_shot=iniFocus_rp.su file_out=p0plus.su above=-1 shift=-8 verbose=1 check=0 hw=8
 
+which marchenko
+ldd marchenko
 #apply the Marchenko algorithm
 marchenko file_shot=shotsdx5_rp.su file_tinv=p0plus.su nshots=901 verbose=2 \
 	tap=0 niter=8 hw=8 shift=12 smooth=3 \

marchenko/demo/twoD/eps.scr

+#!/bin/bash
+
+#model
+supsimage hbox=4 wbox=6 labelsize=12 < syncl_cp.su \
+        x1beg=0 x1end=1800.0 d1num=200 lstyle=vertright legend=1 threecolor=1 \
+        label1="depth (m)" label2="lateral distance (m)" wrgb=1.0,.5,0 grgb=0,.7,1.0 brgb=0,1.0,0 \
+        n1tic=5 x2beg=-3000 f2num=-3000 d2num=1000 x2end=3000 > model_cp.eps
+
+#intial focus operator
+file=iniFocus_z1100_x0_rp.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1/2}'`
+
+supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file \
+    label1="time (s)" label2="lateral distance (m)" \
+    n1tic=2 x1beg=0 x1end=2.0 d1num=0.4 \
+    f2=-3000 f2num=-3000 d2num=1000 clip=$clipref > ${file_base}.eps
+convert -quality 90 -antialias ${file_base}.eps ${file_base}.jpg
+
+supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < p0plus.su \
+    label1="time (s)" label2="lateral distance (m)" \
+    n1tic=2 x1beg=0 x1end=2.0 d1num=0.4 \
+    f2=-3000 f2num=-3000 d2num=1000 clip=$clipref > ${file_base}_muted.eps
+convert ${file_base}_muted.eps ${file_base}_muted.png
+
+
+#reference result
+file=referenceP_rp.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1/4}'`
+
+supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file \
+    label1="time (s)" label2="lateral distance (m)" \
+    n1tic=2 x1beg=0 x1end=2.0 d1num=0.4 \
+    f2=-3000 f2num=-3000 d2num=1000 clip=$clipref > ${file_base}.eps
+convert ${file_base}.eps ${file_base}.png
+
+
+#Marchencko computed Greens function
+file=pgreen.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1/4}'`
+
+supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file \
+    label1="time (s)" label2="lateral distance (m)" \
+    n1tic=2 x1beg=0 x1end=2.0 d1num=0.4 \
+    f2=-3000 f2num=-3000 d2num=1000 clip=$clipref > ${file_base}.eps
+convert ${file_base}.eps ${file_base}.png
+
+

marchenko/demo/twoD/initialFocus.scr

+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+dx=2.5
+dt=0.0005
+
+makemod sizex=6000 sizez=2000 dx=$dx dz=$dx cp0=1900  ro0=1200 \
+        orig=-3000,0 file_base=synclDown.su verbose=2 \
+        intt=def x=-3000,500,3000 z=195,195,195 poly=1 cp=1950 ro=3700 \
+        intt=def x=-3000,3000 z=600,600 poly=0 cp=2050 ro=1750 \
+        intt=def x=-3000,3000 z=680,680 poly=0 cp=2150 ro=2220 \
+        intt=def x=-3000,3000 z=780,780 poly=0 cp=2230 ro=1700 \
+        intt=def x=-3000,-2200,-1500,0,1300,2100,3000 z=520,580,680,840,680,600,500 poly=2 cp=2400 ro=2800 \
+
+makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+
+dxshot=10
+ishot=300
+nshots=301
+
+export KMP_AFFINITY=disabled
+export OMP_NUM_THREADS=16
+mkdir -p shots
+mkdir -p jobs
+
+while (( ishot < nshots ))
+do
+
+	(( xsrc = -3000 + ${ishot}*${dxshot} ))
+	echo xsrc=$xsrc
+	file_rcv=iniFocus_z1100_x${xsrc}.su
+	
+fdelmodc \
+	file_cp=synclDown_cp.su ischeme=1 iorder=4 \
+   	file_den=synclDown_ro.su \
+   	file_src=wave.su \
+   	file_rcv=$file_rcv \
+   	src_type=1 \
+	src_orient=1 \
+	src_injectionrate=1 \
+   	rec_type_vz=0 \
+   	rec_type_p=1 \
+   	rec_int_vz=2 \
+	rec_delay=0.1 \
+   	dtrcv=0.004 \
+   	verbose=2 \
+   	tmod=2.100 \
+   	dxrcv=10.0 \
+   	xrcv1=-3000 xrcv2=3000 \
+   	zrcv1=0 zrcv2=0 \
+   	xsrc=$xsrc zsrc=1100 \
+   	ntaper=200 \
+   	left=2 right=2 top=2 bottom=2
+
+	(( ishot = $ishot + 1))
+done
+
+
+

marchenko/demo/twoD/initialFocus_pbs.scr

@@ -33,17 +33,18 @@ do
 cat << EOF > jobs/pbs_$ishot.job 
 #!/bin/bash
 #
-#PBS -q medium
+#PBS -q bw44-sm
 #PBS -N mod_${xsrc}
 #PBS -j eo 
-#PBS -m n 
-#PBS -l nodes=1
+#PBS -l place=scatter,select=1
+#PBS -l walltime=00:10:00
 #PBS -V
 
 export PATH=\$HOME/src/OpenSource/bin:\$PATH:
 cd \$PBS_O_WORKDIR
 
-export OMP_NUM_THREADS=4
+export KMP_AFFINITY=disabled
+export OMP_NUM_THREADS=16
 	
 fdelmodc \
 	file_cp=synclDown_cp.su ischeme=1 iorder=4 \

marchenko/demo/twoD/marchenko.scr

@@ -5,15 +5,15 @@ export PATH=$HOME/src/OpenSource/bin:$PATH:
 export OMP_NUM_THREADS=40
 
 #mute all events below the first arrival to get the intial focusing field
-fmute file_shot=shots/iniFocus_z1100_x0_rp.su file_out=p0plus.su above=-1 shift=-10 verbose=1 check=1 hw=4
+fmute file_shot=iniFocus_z1100_x0_rp.su file_out=p0plus.su above=-1 shift=-10 verbose=1 check=1 hw=4
 
 #apply the Marchenko algorithm
-../../marchenko file_shot=shots/refl_rp.su file_tinv=p0plus.su nshots=601 verbose=1 \
-	tap=0 niter=15 hw=8 shift=7 smooth=3 \
-    file_green=pgreen2.su file_gplus=Gplus0.su file_gmin=Gmin0.su  \
+marchenko file_shot=shots/refl_rp.su file_tinv=p0plus.su nshots=601 verbose=1 \
+    tap=0 niter=15 hw=8 shift=7 smooth=3 \
+    file_green=pgreen.su file_gplus=Gplus0.su file_gmin=Gmin0.su  \
     file_f1plus=f1plus0.su file_f1min=f1min0.su file_f2=f2.su
 
-#compare Green's funtions on Marhcenko and reference result
+#compare Green's funtions between Marchenko and reference result
 suwind key=gx min=0 max=0 itmax=511 < pgreen.su | sumax mode=abs outpar=nepmg
 suwind key=gx min=0 max=0 itmax=511 < referenceP_rp.su | sumax mode=abs outpar=neprf
 mg=`cat nepmg | awk '{print $1}'`

marchenko/demo/twoD/primaries.scr

@@ -15,14 +15,14 @@ export OMP_NUM_THREADS=40
 select=301
 
 makewave fp=20 dt=0.004 file_out=wave.su nt=1024 t0=0.0 scale=0 scfft=1
-suwind key=fldr min=$select max=$select < shots/refl_rp.su > shotsx0.su
-fconv file_in1=shotsx0.su file_in2=wave.su file_out=shotw.su
 
 marchenko_primaries file_shot=shots/refl_rp.su ishot=$select file_src=wave.su \
 	nshots=601 verbose=2 istart=40 iend=650 fmax=90 \
 	niter=22 niterskip=50 shift=20 file_rr=pred_rr.su T=0
 
 #alternative use shotw as input, must first be multiplied by -1 (see theory)
+#suwind key=fldr min=$select max=$select < shots/refl_rp.su > shotsx0.su
+#fconv file_in1=shotsx0.su file_in2=wave.su file_out=shotw.su
 #sugain scale=-1 < shotw.su > pplus.su
 #
 #marchenko_primaries file_shot=shots/refl_rp.su file_tinv=pplus.su \

marchenko/demo/twoD/referenceShot.scr

@@ -9,7 +9,7 @@ dt=0.0005
 
 makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
 
-export OMP_NUM_THREADS=8
+export OMP_NUM_THREADS=12
 
 fdelmodc \
     file_cp=syncl_cp.su ischeme=1 iorder=4 \
@@ -17,21 +17,21 @@ fdelmodc \
     file_src=wave.su \
     file_rcv=referenceP.su \
     src_type=1 \
-	src_orient=1 \
-	src_injectionrate=1 \
+    src_orient=1 \
+    src_injectionrate=1 \
     rec_type_ud=1 \
     rec_type_p=1 \
     rec_int_vz=2 \
     dtrcv=0.004 \
-	rec_delay=0.1 \
+    rec_delay=0.1 \
     verbose=2 \
     tmod=2.144 \
     dxrcv=10.0 \
     xrcv1=-3000 xrcv2=3000 \
     zrcv1=0 zrcv2=0 \
     xsrc=0 zsrc=1100 \
-	file_snap=backpropref.su tsnap1=0.1 dtsnap=0.010 tsnap2=2.100 \
-	dxsnap=10 dzsnap=10 zsnap1=0 zsnap2=2000 xsnap1=-2250 xsnap2=2250 sna_type_vz=0 \
+    file_snap=backpropref.su tsnap1=0.1 dtsnap=0.010 tsnap2=2.100 \
+    dxsnap=10 dzsnap=10 zsnap1=0 zsnap2=2000 xsnap1=-2250 xsnap2=2250 sna_type_vz=0 \
     ntaper=200 \
     left=2 right=2 top=2 bottom=2
 

marchenko/fmute.c

@@ -79,6 +79,7 @@ int main (int argc, char **argv)
 
     if(!getparstring("file_mute", &file_mute)) file_mute=NULL;
     if(!getparstring("file_shot", &file_shot)) file_shot=NULL;
+	assert(file_shot!=NULL);
     if(!getparstring("file_out", &file_out)) file_out=NULL;
     if(!getparint("ntmax", &ntmax)) ntmax = 1024;
     if(!getparint("nxmax", &nxmax)) nxmax = 512;

marchenko/synthesis.c

@@ -87,7 +87,7 @@ nw, int nw_low, int nw_high,  int mode, int reci, int nshots, int *ixpos, int np
     memset(&iRN[0], 0, Nfoc*nxs*nts*sizeof(float));
     ctrace = (complex *)calloc(ntfft,sizeof(complex));
 
-	/* this first check is done to support an acquisition geometry that has more receiver than source
+/* this first check is done to support an acquisition geometry that has more receiver than source
  * postions. In the first iteration the int R(x_r,x_s) Fop(x_r) d x_r results in a grid on x_s. 
  * so for the next interations onlt x_s traces have to be computed on Fop */
     if (!first) {
@@ -131,14 +131,13 @@ nw, int nw_low, int nw_high,  int mode, int reci, int nshots, int *ixpos, int np
     t1 = wallclock_time();
     *tfft += t1 - t0;
 
-/* Loop over total number of shots */
     if (reci == 0 || reci == 1) {
 
 /*================ SYNTHESIS ================*/
 
 #pragma omp parallel default(none) \
  shared(iRN, dx, npe, nw, verbose, nshots, xnx) \
- shared(Refl, Nfoc, reci, xrcv, xsrc, xsyn, fxsb, fxse, nxs, dxs) \
+ shared(Refl, Nfoc, reci, xsrc, xsyn, fxsb, fxse, nxs, dxs) \
  shared(nx, dxsrc, nfreq, nw_low, nw_high) \
  shared(Fop, size, nts, ntfft, scl, ixrcv) \
  private(l, ix, j, m, i, sum, rtrace, k, ixsrc, inx)
@@ -146,6 +145,7 @@ nw, int nw_low, int nw_high,  int mode, int reci, int nshots, int *ixpos, int np
         sum   = (complex *)malloc(nfreq*sizeof(complex));
         rtrace = (float *)calloc(ntfft,sizeof(float));
 
+/* Loop over total number of shots */
 #pragma omp for schedule(guided,1)
         for (k=0; k<nshots; k++) {
             if ((xsrc[k] < 0.999*fxsb) || (xsrc[k] > 1.001*fxse)) continue;
@@ -191,7 +191,7 @@ nw, int nw_low, int nw_high,  int mode, int reci, int nshots, int *ixpos, int np
 
 #pragma omp parallel default(none) \
  shared(iRN, dx, nw, verbose) \
- shared(Refl, Nfoc, reci, xrcv, xsrc, xsyn, fxsb, fxse, nxs, dxs) \
+ shared(Refl, Nfoc, reci, xsrc, xsyn, fxsb, fxse, nxs, dxs) \
  shared(nx, dxsrc, nfreq, nw_low, nw_high) \
  shared(reci_xrcv, reci_xsrc, ixmask, isxcount) \
  shared(Fop, size, nts, ntfft, scl, ixrcv) \