diff --git a/FFTlib/Makefile b/FFTlib/Makefile
index b6e4faf67ca1a6e79d25089228b2937f4e131c2b..bf53a5af2097ef8891af3fdb526acfea4c47070b 100644
--- a/FFTlib/Makefile
+++ b/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)
diff --git a/Make_include_swan b/Make_include_swan
new file mode 100644
index 0000000000000000000000000000000000000000..b37cffbbd5a282e872b94c0ba45d5400cddb6cb7
--- /dev/null
+++ b/Make_include_swan
@@ -0,0 +1,165 @@
+# 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) $<
+
diff --git a/Makefile b/Makefile
index 74e45defbad119dcd86d34c28ca92712ca38f366..ebce8175a55a82205d6f42cbfefad78c3c684530 100644
--- a/Makefile
+++ b/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) $@
diff --git a/README b/README
index fb6b56317cde12ca8e35b6eadf52418f0c5fad2e..f7643f03892752fe21ca916c1498ea2cbcf9af2e 100644
--- a/README
+++ b/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
----
diff --git a/marchenko/demo/WS15/README.1 b/marchenko/demo/WS15/README.1
new file mode 100644
index 0000000000000000000000000000000000000000..dcc4e80c872e44bf38dde9b205c336a03a2475fe
--- /dev/null
+++ b/marchenko/demo/WS15/README.1
@@ -0,0 +1,66 @@
+
+*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
+
diff --git a/marchenko/demo/WS15/README.2 b/marchenko/demo/WS15/README.2
new file mode 100644
index 0000000000000000000000000000000000000000..aa9603a8f39ff8773bbfd22d07e2a8b03efda6f3
--- /dev/null
+++ b/marchenko/demo/WS15/README.2
@@ -0,0 +1,36 @@
+
+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).
+
+
+
diff --git a/marchenko/demo/WS15/README.3 b/marchenko/demo/WS15/README.3
new file mode 100644
index 0000000000000000000000000000000000000000..37f0f302e0cc317af54a956bd8f75e61d3262c83
--- /dev/null
+++ b/marchenko/demo/WS15/README.3
@@ -0,0 +1 @@
+plane waves
diff --git a/marchenko/demo/WS15/README.4 b/marchenko/demo/WS15/README.4
new file mode 100644
index 0000000000000000000000000000000000000000..b6da83df289fa2f49cabe03af2abe301a4403847
--- /dev/null
+++ b/marchenko/demo/WS15/README.4
@@ -0,0 +1,3 @@
+
+
+Primaries only
diff --git a/marchenko/demo/WS15/job.pbs b/marchenko/demo/WS15/job.pbs
new file mode 100755
index 0000000000000000000000000000000000000000..13eef276b548ae634b570172e99460fa976f56cc
--- /dev/null
+++ b/marchenko/demo/WS15/job.pbs
@@ -0,0 +1,25 @@
+#!/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"
+
diff --git a/marchenko/demo/WS15/setup.sh b/marchenko/demo/WS15/setup.sh
new file mode 100755
index 0000000000000000000000000000000000000000..40c56979d607430257b6d88c8e98b520551e6c94
--- /dev/null
+++ b/marchenko/demo/WS15/setup.sh
@@ -0,0 +1,12 @@
+
+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
+
diff --git a/marchenko/demo/invisible/model.scr b/marchenko/demo/invisible/model.scr
index cbd96ebdb49565f710ebc13bbfdf7ac1a0befbf2..22f93aa9ca612ca90d5bd22537b707bd41abf582 100755
--- a/marchenko/demo/invisible/model.scr
+++ b/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 \
diff --git a/marchenko/demo/oneD/marchenko.scr b/marchenko/demo/oneD/marchenko.scr
index 04e495d53dcb305750132b3fa55dc4606d444b46..69de0f8f9a28b37cd89c756eb2d9ca7822585830 100755
--- a/marchenko/demo/oneD/marchenko.scr
+++ b/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 \
diff --git a/marchenko/demo/twoD/eps.scr b/marchenko/demo/twoD/eps.scr
new file mode 100755
index 0000000000000000000000000000000000000000..f66a8bbdb3a6325a78e581233dc0a573628543e8
--- /dev/null
+++ b/marchenko/demo/twoD/eps.scr
@@ -0,0 +1,53 @@
+#!/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
+
+
diff --git a/marchenko/demo/twoD/initialFocus.scr b/marchenko/demo/twoD/initialFocus.scr
new file mode 100755
index 0000000000000000000000000000000000000000..fcc72f51099edf986a557647b32e95d0dce6fcb4
--- /dev/null
+++ b/marchenko/demo/twoD/initialFocus.scr
@@ -0,0 +1,60 @@
+#!/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
+
+
+
diff --git a/marchenko/demo/twoD/initialFocus_pbs.scr b/marchenko/demo/twoD/initialFocus_pbs.scr
index eb4e0c1d48c1a8533905a0ba3a0f6092dc48897d..22df20cfb31e96084159304d169828fba99cfc3d 100755
--- a/marchenko/demo/twoD/initialFocus_pbs.scr
+++ b/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 \
diff --git a/marchenko/demo/twoD/marchenko.scr b/marchenko/demo/twoD/marchenko.scr
index 24131122208386b9ce51a263e9902c84ae5da060..88827395a712e353357876737b1c24bb2646b9f7 100755
--- a/marchenko/demo/twoD/marchenko.scr
+++ b/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}'`
diff --git a/marchenko/demo/twoD/primaries.scr b/marchenko/demo/twoD/primaries.scr
index 7ec27e0c712832d720e734493f8cbc397f73a645..d001b743b833a14ad141e024d9be64a56c504b29 100755
--- a/marchenko/demo/twoD/primaries.scr
+++ b/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 \
diff --git a/marchenko/demo/twoD/referenceShot.scr b/marchenko/demo/twoD/referenceShot.scr
index 4c015f6baa98c092e3210801e3a6caff68fd34ad..abd00a179eb84fca06d2d7aefed28e2847dcb1de 100755
--- a/marchenko/demo/twoD/referenceShot.scr
+++ b/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
diff --git a/marchenko/fmute.c b/marchenko/fmute.c
index f5c7668a16dc382b586c2981427888c6fa256f66..4882a0222d3bb436e280fcb042ce78fc51e1ebca 100644
--- a/marchenko/fmute.c
+++ b/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;
diff --git a/marchenko/synthesis.c b/marchenko/synthesis.c
index 4a8831315428229705ea54f843d5ee0458da2edd..2118510d70aac94a0bd077cbab3b499152ce1ddb 100644
--- a/marchenko/synthesis.c
+++ b/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) \