working on Marchenko planewave demo

extrap3d/lib/Makefile

@@ -37,8 +37,11 @@ CSRC = \
 
 #if there is a Fortran compiler you can also use Remez optimised operators
 # uncomment this line to do so.
+#ifneq ($(strip $(FC)),)
 #FSRC = \
-#	optRemez.f \
+#	optRemez.f 
+#CFLAGS += -DREM 
+#endif
 
 OBJC    = $(CSRC:%.c=%.o)
 OBJF    = $(FSRC:%.f=%.o)

marchenko/demo/planewave/Ini2Focus.scr

+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+dx=2.5
+dt=0.0005
+
+makewave fp=25 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+
+seq -3000 2.5 3000 > srcx.txt
+seq 600 0.40000 1560.1 > srcz.txt
+paste srcx.txt srcz.txt > srcxz.txt
+
+
+cat << EOF > slurm.job 
+#!/bin/bash
+#SBATCH --cpus-per-task=4
+#SBATCH --ntasks=1
+#SBATCH -J mod
+#SBATCH -V
+#SBATCH -p max2h
+
+export PATH=\$HOME/src/OpenSource/bin:\$PATH:
+
+export OMP_NUM_THREADS=4
+varo=SourceDipa.su
+
+dx=2.5
+dt=0.0005
+
+fdelmodc \
+        file_cp=ge_cp.su ischeme=1 iorder=4 \
+        file_den=ge_ro.su \
+        file_src=wave.su \
+        file_rcv=\$varo \
+        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.500 \
+        dxrcv=10.0 \
+		src_txt=srcxz.txt \
+        tsnap1=0.1 tsnap2=0.5 file_snap=snapDipa.su \
+        xrcv1=-3000 xrcv2=3000 \
+        zrcv1=0 zrcv2=0 \
+        left=2 right=2 top=2 bottom=2\
+
+EOF
+
+sbatch slurm.job
+
+

marchenko/demo/planewave/IniFocus.scr

@@ -42,7 +42,7 @@ varq=SourceNew-${angle}a${a}.su
 dx=2.5
 dt=0.0005
 
-/vardim/home/thorbcke/src/OpenSource/fdelmodc/fdelmodc \
+fdelmodc \
         file_cp=ge_cp.su ischeme=1 iorder=4 \
         file_den=ge_ro.su \
         file_src=wave.su \
@@ -67,7 +67,7 @@ dt=0.0005
         left=2 right=2 top=2 bottom=2\
 
 
-/vardim/home/thorbcke/src/OpenSource/fdelmodc/fdelmodc \
+fdelmodc \
         file_cp=ge_cp.su ischeme=1 iorder=4 \
         file_den=ge_ro.su \
         file_src=wave.su \

marchenko/demo/planewave/README

+Code to reproduce the example in section `Plane-Wave Marchenko algorithm` of the manual.
+
+Description of files:
+
+
+1) ./model.scr computes the gridded velocity/density model 
+
+2) ./shots.scr computes 601 shots using slurm arrays: ~100 s runtime for each shot
+
+3) ./direct.scr compute 1 shot that contains direct wave only
+   
+4) ./remove_direct.scr removes the direct wave from the data created in 2 and places all shots in one file shots/refl_rp.su 
+
+   Figures of the model and the middle shot are generated by
+   ./epsModel.scr model
+
+5) IniFocus.scr compute plane wave responses at 600 m depth with angles -5 0 and 5 degrees.
+   To compute the plane wave respone with the sources all startin at t=0 on a slanted line in the model run
+   ./Ini2Focus.scr
+   To generate the 12 snapshot figures
+   ./epsModel.scr snapshots
+
+6) To generate the Marchenko results for plane wave with angle=0 (horizontal)
+
+./marchenkoPlaneIter.scr plane0
+./epsModel.scr plane0
+
+[ 7) Marchenko with and without time-shifts in a one-D model: requires shift to build in again in marchenko code: does not generate
+correct results at the moment
+   cd shifts
+   compute 1D model, f_d and marchenko with and without shift
+   ./model.scr
+   ./initialFocus.scr
+   ./marchenko.scr
+   ./eps.scr
+]
+
+
+2) itertions.scr computes the intermediate results of the multiple attenutation scheme and produces all output files that are used in the manuscript.
+   - runtime on 4 cores is 
+3) epsPrimaries.scr selected output from step 2) are converted to .eps pictures that are used in the Figures to explain the method.
+   To reproduce the postscript files of the manuscript SU postscript plotting programs are required.
+-3) epsModel.scr to generate the postscript files for the numerical model
+
+optional scripts not needed to reproduce the figures:
++) primaries.scr computes the internal multiple attenuated (middle) shot record.
+   - runtime on 4 cores is ~500 s.
++) primariesPlane.scr: computes the internal moval scheme for plane-waves (see Meles 2020)
++) clean: remove all produced files and start with a clean directory
+
+
+--------------------------
+To reproduce the Figures in the Manuscript:
+--------------------------
+* Figure 2: Model + Initial wavefield
+
+==> run model.scr to generate the data .su files: this will take 4-5 minutes. The files generate are:
+	- hom_cp.su, hom_ro.su
+	- model10_cp.su, model10_ro.su
+	- shot5_fd_rp.su
+	- shot5_hom_fd_rp.su
+	- shot5_rp.su
+	- wavefw.su
+
+==> run './epsPrimaries.scr Figure2' to generate the postscript files of Figure 2 
+
+model_cp_line.eps 	=> Figure 2a 
+model_ro_line.eps 	=> Figure 2b
+shotx0_rp.eps 		=> Figure 2c
+
+It also produces two extra pictures of the wavelet used in the FD modelling:
+wavefw_freq.eps
+wavefw.eps 
+
+--------------------------
+* Figure 3: First Iteration
+
+==> run './iterations.scr Figure34910' to compute the intermediate results of the first iterations of the Marchenko Primaries algorithm.
+This will take 15 seconds. The generated files are:
+
+	- M0_276000.su
+	- Mi_2760##.su
+	- k1min_2760##.su
+	- v1plus_2760##.su
+	- iter_2760##.su (not used)
+	- pred_rr_276.su (not used)
+	- DDshot_450.su (not used): selected shot record convolved with file_src=wave.su
+
+	where ## ranges from 01 to 34
+
+To generate the postscript files for Figure 3:
+
+==> run './epsPrimaries.scr Figure3' 
+
+This will produce the following files:
+
+shotx0_rp.eps 		=> Figure 2c == Figure 3a
+M0_276000_flip.eps	=> Figure 3b
+fconvN0fulltime.eps => Figure 3c
+fconvN0flip.eps 	=> Figure 3d
+Mi_276001.eps		=> Figure 3e
+
+
+--------------------------
+* Figure 4 second iteration
+To generate the postscript files for Figure 4:
+==> run './epsPrimaries.scr Figure4' 
+
+This will produce the following files:
+
+fconvN1fulltime.eps	=> Figure 4c
+fconvN1flip.eps		=> Figure 4d
+Mi_276002.eps		=> Figure 4e
+
+The window time function in Figure 5 is not reproduced. 
+
+--------------------------
+* Figure 6 v1plus and convergence
+==> run './iterations.scr Figure6' to compute the marchenko results with ii=200
+
+To generate the postscript files for Figure 6:
+==> run './epsPrimaries.scr Figure6' 
+
+This will produce the following files:
+v1plus_200001.eps		=> Figure 6a
+v1plus_max.eps			=> Figure 6b
+k1min_200030.eps		=> Figure 6b
+
+--------------------------
+* Figure 8 To compute the convergence for a strong contrast medium:
+cd strongContrast
+==> run ./model.scr
+==> run ./iterations.scr Figure8
+
+To generate the postscript files for Figure 8:
+==> run './epsPrimaries.scr Figure8' 
+
+This will produce the following files:
+v1plusStrong_max.eps		=> Figure 8
+
+Don't forget to go back to the main directory with the regular contrast results
+cd ../
+
+
+--------------------------
+* Figure 9 iterations M_i
+
+To generate the postscript files for Figure 9:
+==> run './epsPrimaries.scr Figure9' 
+
+This will produce the following files:
+Mi_276002.eps	=> Figure 9b
+Mi_276004.eps	=> Figure 9d
+Mi_276012.eps	=> Figure 9f
+Mi_276020.eps	=> Figure 9h
+Mi_276001.eps	=> Figure 9a
+Mi_276003.eps	=> Figure 9c
+Mi_276011.eps	=> Figure 9e
+Mi_276019.eps	=> Figure 9g
+
+
+--------------------------
+* Figure 10 iterations M_i and k_1^-
+
+To generate the postscript files for Figure 10:
+==> run './epsPrimaries.scr Figure10' 
+
+This will produce the following files:
+Mi_276002flip.eps	=> Figure 10a
+k1min_276002.eps	=> Figure 10b
+Mi_276004flip.eps	=> Figure 10c
+k1min_276004.eps	=> Figure 10d
+Mi_276010flip.eps	=> Figure 10e
+k1min_276010.eps	=> Figure 10f
+Mi_276020flip.eps	=> Figure 10g
+k1min_276020.eps	=> Figure 10h
+
+
+--------------------------
+* Figure 11 iterations k_1^- for different ii 246:316:10
+
+To generate the data 
+==> run ./iterations.scr Figure11
+this will take ~2 minutes and generate a lot of files
+
+To generate the postscript files for Figure 11:
+==> run './epsPrimaries.scr Figure11' 
+
+This will produce the following files:
+k1min_246032.eps	=> Figure 11a
+k1min_256032.eps	=> Figure 11b
+k1min_266032.eps	=> Figure 11c
+k1min_276032.eps	=> Figure 11d
+k1min_286032.eps	=> Figure 11e
+k1min_296032.eps	=> Figure 11f
+k1min_306032.eps	=> Figure 11g
+k1min_316032.eps	=> Figure 11h
+
+
+--------------------------
+* Figure 13 iterations M_i and k_1^- for ii-276 T-MME scheme
+
+To generate the data 
+==> run ./iterations.scr Figure13
+this will take ~15 seconds
+****NOTE this will overwrite the results of the MME-scheme !
+
+To generate the postscript files for Figure 13:
+==> run './epsPrimaries.scr Figure13' 
+
+This will produce the following files:
+Mi_276002T.eps	=> Figure 13a
+Mi_276004T.eps	=> Figure 13b
+Mi_276001T.eps	=> Figure 13c
+Mi_276003T.eps	=> Figure 13d
+k1min_276002T.eps	=> Figure 13e
+k1min_276004T.eps	=> Figure 13f
+k1min_276010T.eps	=> Figure 13g
+k1min_276020T.eps	=> Figure 13h
+
+
+
+

marchenko/demo/planewave/clean

+#!/bin/bash
+
+rm *.su *.bin *.eps nep line* *.asci src*.txt
+

marchenko/demo/planewave/epsModel.scr

@@ -3,6 +3,13 @@
 srcname=Source
 srcname=SourceNew
 
+ntfft=1024
+ns=602
+nsd=351 #number of samples to display
+(( nstart = ns - nsd ))
+
+if [[ "$1" == "model" ]];
+then
 supsimage hbox=4 wbox=6 labelsize=12 < ge_cp.su \
         x1beg=0 x1end=1400.0 d1num=200 lstyle=vertright legend=1 threecolor=0 \
         label1="depth (m)" label2="lateral distance (m)" \
@@ -15,14 +22,6 @@ supsimage hbox=4 wbox=6 labelsize=12 < ge_ro.su \
 
 suwind < shots/refl_rp.su key=sx min=0 max=0 > shot0_rp.su
 
-ntfft=1024
-ns=602
-nsd=351 #number of samples to display
-(( nstart = ns - nsd ))
-
-sumax < ${srcname}0a120_rp.su mode=abs outpar=nep
-clipS=`cat nep | awk '{print $1/7}'`
-
 file=shot0_rp.su
 sumax < $file mode=abs outpar=nep
 clipf1=`cat nep | awk '{print $1/7}'`
@@ -31,6 +30,33 @@ supsimage < $file hbox=6 wbox=4 labelsize=8 linewidth=0.0 \
       n1tic=2 x1end=1.4 f1=0 \
       label1="time [s]" label2="lateral distance [m]" \
       f2=-3000 f2num=-3000 d2num=1500 clip=$clipf1 > $file_base.eps
+exit
+fi
+
+if [[ "$1" == "snapshots" ]];
+then
+#plane-wave snapshots
+for file in snap5a120_sp.su snapDipa_sp.su snap-5a120_sp.su 
+do 
+  sumax < $file mode=abs outpar=nep
+  clipf1=`cat nep | awk '{print $1/7}'`
+
+  file_base=${file%.su}
+  for fldr in 1 2 3 4
+  do
+    suwind < $file key=fldr min=$fldr max=$fldr > shottmp.su
+    supsimage < shottmp.su hbox=4 wbox=4 labelsize=8 linewidth=0.0 \
+      n1tic=2 f1=0 \
+      label1="depth [m]" label2="lateral distance [m]" \
+      f2=-3000 f2num=-3000 d2num=1500 clip=$clipf1 > ${file_base}_$fldr.eps
+
+   done
+done
+exit
+fi
+
+sumax < ${srcname}0a120_rp.su mode=abs outpar=nep
+clipS=`cat nep | awk '{print $1/7}'`
 
 #mute windows for angle=5 and angle=-5
 for file in mute05_001.su mute05_002.su mute0-5_001.su mute0-5_002.su #mute0-5s_001.su mute0-5s_002.su
@@ -77,24 +103,6 @@ do
 
 done
 
-#plane-wave snapshots
-for file in snap5a120_sp.su snapDipa_sp.su snap-5a120_sp.su snap-5Aa120_sp.su
-do 
-  sumax < $file mode=abs outpar=nep
-  clipf1=`cat nep | awk '{print $1/7}'`
-
-  file_base=${file%.su}
-  for fldr in 1 2 3 4
-  do
-    suwind < $file key=fldr min=$fldr max=$fldr > shottmp.su
-    supsimage < shottmp.su hbox=4 wbox=4 labelsize=8 linewidth=0.0 \
-      n1tic=2 f1=0 \
-      label1="depth [m]" label2="lateral distance [m]" \
-      f2=-3000 f2num=-3000 d2num=1500 clip=$clipf1 > ${file_base}_$fldr.eps
-
-   done
-done
-
 #iterrations and f1plus and f1min
 
 for file in f1plus${srcname}5a120.su f1min${srcname}5a120.su 

marchenko/demo/planewave/marchenkoPlaneIter.scr

@@ -4,12 +4,25 @@ ROOT=/vardim/home/thorbcke/R1/
 ROOT=/vardim/home/thorbcke/src/
 export PATH=$HOME/src/OpenSource/bin:$PATH:
 export OMP_NUM_THREADS=4
-exe=$ROOT/OpenSource/marchenko/marchenkowd
-exe=$ROOT/OpenSource/marchenko/marchenko
+#exe=$ROOT/OpenSource/marchenko/marchenkowd
 #exe=$ROOT/OpenSource/bin/marchenko
+exe=$ROOT/OpenSource/bin/marchenko
 
-#5 degrees dip without (set src_angle=0) tshift for Gmin in marchenko program (recompiled)
-base=Source5a120_shiftedIter
+#zero angle
+if [[ "$1" == "plane0" ]];
+then
+src_angle=0
+inifile=SourceNew${src_angle}a120_rp.su
+base=${inifile%_rp.su}
+fmute file_shot=$inifile file_out=fdplus${base}.su above=-1 shift=-6 verbose=1 check=0 hw=2
+$exe file_shot=shots/refl_rp.su file_tinv=fdplus${base}.su nshots=601 verbose=2 \
+        tap=3 niter=16 hw=2 shift=7 smooth=3 plane_wave=0 src_angle=0 rotate=0 \
+        file_green=green${base}.su file_gplus=Gplus${base}.su file_gmin=Gmin${base}.su  \
+        file_f1plus=f1plus${base}.su file_f1min=f1min${base}.su file_f2=f2${base}.su 
+exit
+fi
+
+#5 degrees dip without (set src_angle=0) 
 base=Source5a120old
 
 src_angle=5
@@ -25,15 +38,6 @@ $exe file_shot=shots/refl_rp.su file_tinv=fdplus${base}.su nshots=601 verbose=2
         file_f1plus=f1plus${base}.su file_f1min=f1min${base}.su file_f2=f2${base}.su file_iter=iter${base}.su 
 
 #exit
-#zero angle
-src_angle=0
-inifile=SourceNew${src_angle}a120_rp.su
-base=${inifile%_rp.su}
-fmute file_shot=$inifile file_out=fdplus${base}.su above=-1 shift=-6 verbose=1 check=0 hw=2
-$exe file_shot=shots/refl_rp.su file_tinv=fdplus${base}.su nshots=601 verbose=2 \
-        tap=3 niter=16 hw=2 shift=7 smooth=3 plane_wave=0 src_angle=0 rotate=0 \
-        file_green=green${base}.su file_gplus=Gplus${base}.su file_gmin=Gmin${base}.su  \
-        file_f1plus=f1plus${base}.su file_f1min=f1min${base}.su file_f2=f2${base}.su 
 
 for angle in 5 #10
 do

marchenko/demo/planewave/remove_direct.scr

@@ -23,14 +23,12 @@ do
 	suwind < direct_rp.su key=tracl min=$tr1 max=$tr2 > direct.su
 
 	file_rcv=shots/shots_${xsrc}_rp.su
-	suwind key=tracl min=1 max=601 < $file_rcv > shotz0.su
-
-	sudiff shotz0.su direct.su > refl.su
 
 	(( ishot = $ishot + 1))
 
-	sushw < refl.su key=fldr a=$ishot | \
-	suwind itmax=1023 >> shots/refl_rp.su
+	sudiff $file_rcv direct.su | \
+    	sushw key=fldr a=$ishot | \
+	    suwind itmax=1023 >> shots/refl_rp.su
 
 done
 

marchenko/demo/planewave/shots.scr

@@ -20,51 +20,45 @@ ishot=0
 nshots=601
 zsrc=0
 
-while (( ishot < nshots ))
-do
-
-		(( xsrc = -3000 + ${ishot}*${dxshot} ))
-
-		echo ishot=$ishot xsrc=$xsrc zsrc=$zsrc
-
-  cat << EOF > jobs/pbs_$ishot.job 
+cat << EOF > jobs/slurm.job 
 #!/bin/bash
 #SBATCH --cpus-per-task=4
 #SBATCH --ntasks=1
-#SBATCH -J mod_${xsrc}
+#SBATCH --array=0-600
+#SBATCH -J mod2D
 #SBATCH -V
+#SBATCH -p max30m
 
+(( xsrc = -3000 + \${SLURM_ARRAY_TASK_ID}*${dxshot} ))
 export PATH=\$HOME/src/OpenSource/bin:\$PATH:
 
-	export OMP_NUM_THREADS=4
-	file_rcv=shots/shots_${xsrc}.su
-
-	fdelmodc \
-   		file_cp=ge_cp.su ischeme=1 iorder=4 \
-   		file_den=ge_ro.su \
-   		file_src=wavefw.su \
-   		file_rcv=\$file_rcv \
-		src_type=7 \
-		src_orient=1 \
-		src_injectionrate=1 \
-   		rec_type_vz=1 \
-   		rec_type_p=1 \
-   		rec_int_vz=2 \
-		rec_delay=0.3 \
-   		dtrcv=0.004 \
-   		verbose=2 \
-   		tmod=4.394 \
-   		dxrcv=10.0 \
-   		xrcv1=-3000 xrcv2=3000 \
-   		zrcv1=0 zrcv2=0 \
-   		xsrc=$xsrc zsrc=$zsrc \
-   		ntaper=200 \
-   		left=2 right=2 top=2 bottom=2
+echo ishot=\$SLURM_ARRAY_TASK_ID xsrc=\$xsrc zsrc=$zsrc
+export OMP_NUM_THREADS=4
+file_rcv=shots/shots_\${xsrc}.su
+
+fdelmodc \
+  	file_cp=ge_cp.su ischeme=1 iorder=4 \
+   	file_den=ge_ro.su \
+   	file_src=wavefw.su \
+   	file_rcv=\$file_rcv \
+	src_type=7 \
+	src_orient=1 \
+	src_injectionrate=1 \
+   	rec_type_vz=0 \
+   	rec_type_p=1 \
+   	rec_int_vz=2 \
+	rec_delay=0.3 \
+   	dtrcv=0.004 \
+   	verbose=2 \
+   	tmod=4.394 \
+   	dxrcv=10.0 \
+   	xrcv1=-3000 xrcv2=3000 \
+   	zrcv1=0 zrcv2=0 \
+   	xsrc=\$xsrc zsrc=$zsrc \
+   	ntaper=200 \
+   	left=2 right=2 top=2 bottom=2
 EOF
 
-sbatch jobs/pbs_$ishot.job
-
-   		(( ishot = $ishot + 1))
+#sbatch jobs/slurm.job
 
-done