Skip to content
Snippets Groups Projects
Commit 983533c0 authored by Jan Willem Thorbecke's avatar Jan Willem Thorbecke
Browse files

added marchenko primaries that implemets Lele's scheme

parent ff009af1
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1226 additions and 943 deletions
.gitignore
+ 1
− 0
View file @ 983533c0
......@@ -10,6 +10,7 @@ Make_include
fdelmodc_cuda.tgz
marchenko/fmute
marchenko/marchenko
marchenko/marchenko_primaries
raytime/raytime
utils/basop
utils/correigen
......
marchenko/Makefile
+ 4
− 10
View file @ 983533c0
......@@ -5,7 +5,6 @@ include ../Make_include
#LIBS += -L$L -lgenfft -lm $(LIBSM)
#OPTC += -g -O0 -Wall
#ALL: fmute marchenko marchenko2
ALL: fmute marchenko marchenko_primaries
......@@ -21,6 +20,7 @@ SRCJ = fmute.c \
getpars.c
SRCH = marchenko.c \
synthesis.c \
getFileInfo.c \
readData.c \
readShotData.c \
......@@ -36,6 +36,7 @@ SRCH = marchenko.c \
getpars.c
SRCP = marchenko_primaries.c \
synthesis.c \
getFileInfo.c \
readData.c \
readShotData.c \
......@@ -64,23 +65,16 @@ OBJP = $(SRCP:%.c=%.o)
marchenko_primaries: $(OBJP)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko_primaries $(OBJP) $(LIBS)
OBJH2 = $(SRCH2:%.c=%.o)
marchenko2: $(OBJH2)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko2 $(OBJH2) $(LIBS)
install: fmute marchenko marchenko_primaries
cp fmute $B
cp marchenko $B
cp marchenko_primaries $B
# cp marchenko2 $B
clean:
rm -f core fmute $(OBJJ) marchenko $(OBJH) marchenko2 $(OBJH2) marchenko_primaries $(OBJP)
rm -f core fmute $(OBJJ) marchenko $(OBJH) marchenko_primaries $(OBJP)
realclean: clean
rm -f $B/fmute $B/marchenko $B/marchenko2 $B/marchenko_primaries
rm -f $B/fmute $B/marchenko $B/marchenko_primaries
......
marchenko/demo/twoD/marchenko.scr
+ 3
− 3
View file @ 983533c0
......@@ -2,15 +2,15 @@
export PATH=$HOME/src/OpenSource/bin:$PATH:
export OMP_NUM_THREADS=1
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
#apply the Marchenko algorithm
marchenko file_shot=shots/refl_rp.su file_tinv=p0plus.su nshots=601 verbose=1 \
../../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_green=pgreen2.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
......
marchenko/marchenko.c
+ 7
− 343
View file @ 983533c0
......@@ -25,9 +25,6 @@ void omp_set_num_threads(int num_threads);
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#endif
#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
int compareInt(const void *a, const void *b)
{ return (*(int *)a-*(int *)b); }
#ifndef COMPLEX
typedef struct _complexStruct { /* complex number */
......@@ -56,7 +53,7 @@ Nfoc, float *xrcv, float *xsrc, int *xnx, float fxse, float fxsb, float dxs, flo
nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int npos, double *tfft, int *isxcount, int
*reci_xsrc, int *reci_xrcv, float *ixmask, int verbose);
void synthesisPosistions(int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nfoc, float *xrcv, float *xsrc, int *xnx,
void synthesiPosistions(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 nshots, int *ixpos, int *npos, int *isxcount, int countmin, int reci, int verbose);
int linearsearch(int *array, size_t N, int value);
......@@ -227,7 +224,6 @@ int main (int argc, char **argv)
muteW = (int *)calloc(Nfoc*nxs,sizeof(int));
tsynW = (int *)malloc(Nfoc*nxs*sizeof(int)); // time-shift for Giovanni's plane-wave on non-zero times
trace = (float *)malloc(ntfft*sizeof(float));
tapersy = (float *)malloc(nxs*sizeof(float));
xrcvsyn = (float *)calloc(Nfoc*nxs,sizeof(float)); // x-rcv postions of focal points
xsyn = (float *)malloc(Nfoc*sizeof(float)); // x-src position of focal points
zsyn = (float *)malloc(Nfoc*sizeof(float)); // z-src position of focal points
......@@ -235,7 +231,6 @@ int main (int argc, char **argv)
ixpos = (int *)calloc(nxs,sizeof(int)); // x-position of source of shot in G_d domain (nxs with dxs)
Refl = (complex *)malloc(nw*nx*nshots*sizeof(complex));
tapersh = (float *)malloc(nx*sizeof(float));
xrcv = (float *)calloc(nshots*nx,sizeof(float)); // x-rcv postions of shots
xsrc = (float *)calloc(nshots,sizeof(float)); //x-src position of shots
zsrc = (float *)calloc(nshots,sizeof(float)); // z-src position of shots
......@@ -281,17 +276,13 @@ int main (int argc, char **argv)
/* define tapers to taper edges of acquisition */
if (tap == 1 || tap == 3) {
tapersy = (float *)malloc(nxs*sizeof(float));
for (j = 0; j < ntap; j++)
tapersy[j] = (cos(PI*(j-ntap)/ntap)+1)/2.0;
for (j = ntap; j < nxs-ntap; j++)
tapersy[j] = 1.0;
for (j = nxs-ntap; j < nxs; j++)
tapersy[j] =(cos(PI*(j-(nxs-ntap))/ntap)+1)/2.0;
}
else {
for (j = 0; j < nxs; j++) tapersy[j] = 1.0;
}
if (tap == 1 || tap == 3) {
if (verbose) vmess("Taper for operator applied ntap=%d", ntap);
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < nxs; i++) {
......@@ -300,6 +291,7 @@ int main (int argc, char **argv)
}
}
}
free(tapersy);
}
/* check consistency of header values */
......@@ -318,19 +310,14 @@ int main (int argc, char **argv)
readShotData(file_shot, xrcv, xsrc, zsrc, xnx, Refl, nw, nw_low, nshots, nx, nxs, fxsb, dxs, ntfft,
mode, scale, tsq, Q, f0, reci, &nshots_r, isxcount, reci_xsrc, reci_xrcv, ixmask, verbose);
tapersh = (float *)malloc(nx*sizeof(float));
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++)
tapersh[j] = 1.0;
for (j = nx-ntap; j < nx; j++)
tapersh[j] =(cos(PI*(j-(nx-ntap))/ntap)+1)/2.0;
}
else {
for (j = 0; j < nx; j++) tapersh[j] = 1.0;
}
if (tap == 2 || tap == 3) {
if (verbose) vmess("Taper for shots applied ntap=%d", ntap);
for (l = 0; l < nshots; l++) {
for (j = 1; j < nw; j++) {
......@@ -340,8 +327,8 @@ int main (int argc, char **argv)
}
}
}
free(tapersh);
}
free(tapersh);
/* check consistency of header values */
fxf = xsrc[0];
......@@ -408,10 +395,10 @@ int main (int argc, char **argv)
/* dry-run of synthesis to get all x-positions calcalated by the integration */
synthesisPosistions(nx, nt, nxs, nts, dt, xsyn, Nfoc, xrcv, xsrc, xnx, fxse, fxsb,
synthesisPositions(nx, nt, nxs, nts, dt, xsyn, Nfoc, xrcv, xsrc, xnx, fxse, fxsb,
dxs, dxsrc, dx, nshots, ixpos, &npos, isxcount, countmin, reci, verbose);
if (verbose) {
vmess("synthesisPosistions: nshots=%d npos=%d", nshots, npos);
vmess("synthesisPositions: nshots=%d npos=%d", nshots, npos);
}
/*================ set variables for output data ================*/
......@@ -758,331 +745,8 @@ sqrt(energyNi/energyN0));
/*================ free memory ================*/
free(hdrs_out);
free(tapersy);
exit(0);
}
/*================ Convolution and Integration ================*/
void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, 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
nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int npos, double *tfft, int *isxcount, int
*reci_xsrc, int *reci_xrcv, float *ixmask, int verbose)
{
int nfreq, size, inx;
float scl;
int i, j, l, m, iw, ix, k, ixsrc, il, ik;
float *rtrace, idxs;
complex *sum, *ctrace;
int npe;
static int first=1, *ixrcv;
static double t0, t1, t;
size = nxs*nts;
nfreq = ntfft/2+1;
/* scale factor 1/N for backward FFT,
* scale dt for correlation/convolution along time,
* scale dx (or dxsrc) for integration over receiver (or shot) coordinates */
scl = 1.0*dt/((float)ntfft);
#ifdef _OPENMP
npe = omp_get_max_threads();
/* parallelisation is over number of virtual source positions (Nfoc) */
if (npe > Nfoc) {
vmess("Number of OpenMP threads set to %d (was %d)", Nfoc, npe);
omp_set_num_threads(Nfoc);
}
#endif
t0 = wallclock_time();
/* reset output data to zero */
memset(&iRN[0], 0, Nfoc*nxs*nts*sizeof(float));
ctrace = (complex *)calloc(ntfft,sizeof(complex));
if (!first) {
/* transform muted Ni (Top) to frequency domain, input for next iteration */
for (l = 0; l < Nfoc; l++) {
/* set Fop to zero, so new operator can be defined within ixpos points */
memset(&Fop[l*nxs*nw].r, 0, nxs*nw*2*sizeof(float));
for (i = 0; i < npos; i++) {
rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
ix = ixpos[i];
for (iw=0; iw<nw; iw++) {
Fop[l*nxs*nw+iw*nxs+ix].r = ctrace[nw_low+iw].r;
Fop[l*nxs*nw+iw*nxs+ix].i = mode*ctrace[nw_low+iw].i;
}
}
}
}
else { /* only for first call to synthesis using all nxs traces in G_d */
/* transform G_d to frequency domain, over all nxs traces */
first=0;
for (l = 0; l < Nfoc; l++) {
/* set Fop to zero, so new operator can be defined within all ix points */
memset(&Fop[l*nxs*nw].r, 0, nxs*nw*2*sizeof(float));
for (i = 0; i < nxs; i++) {
rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
for (iw=0; iw<nw; iw++) {
Fop[l*nxs*nw+iw*nxs+i].r = ctrace[nw_low+iw].r;
Fop[l*nxs*nw+iw*nxs+i].i = mode*ctrace[nw_low+iw].i;
}
}
}
idxs = 1.0/dxs;
ixrcv = (int *)malloc(nshots*nx*sizeof(int));
for (k=0; k<nshots; k++) {
for (i = 0; i < nx; i++) {
ixrcv[k*nx+i] = NINT((xrcv[k*nx+i]-fxsb)*idxs);
}
}
}
free(ctrace);
t1 = wallclock_time();
*tfft += t1 - t0;
/* Loop over total number of shots */
if (reci == 0 || reci == 1) {
for (k=0; k<nshots; k++) {
if ((xsrc[k] < 0.999*fxsb) || (xsrc[k] > 1.001*fxse)) continue;
ixsrc = NINT((xsrc[k] - fxsb)/dxs);
inx = xnx[k]; /* number of traces per shot */
/*================ SYNTHESIS ================*/
#pragma omp parallel default(none) \
shared(iRN, dx, npe, nw, verbose) \
shared(Refl, Nfoc, reci, xrcv, xsrc, xsyn, fxsb, fxse, nxs, dxs) \
shared(nx, dxsrc, inx, k, nfreq, nw_low, nw_high) \
shared(Fop, size, nts, ntfft, scl, ixrcv, ixsrc) \
private(l, ix, j, m, i, sum, rtrace)
{ /* start of parallel region */
sum = (complex *)malloc(nfreq*sizeof(complex));
rtrace = (float *)calloc(ntfft,sizeof(float));
#pragma omp for schedule(guided,1)
for (l = 0; l < Nfoc; l++) {
/* compute integral over receiver positions */
/* multiply R with Fop and sum over nx */
memset(&sum[0].r,0,nfreq*2*sizeof(float));
for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
for (i = 0; i < inx; i++) {
ix = ixrcv[k*nx+i];
sum[j].r += Refl[k*nw*nx+m*nx+i].r*Fop[l*nw*nxs+m*nxs+ix].r -
Refl[k*nw*nx+m*nx+i].i*Fop[l*nw*nxs+m*nxs+ix].i;
sum[j].i += Refl[k*nw*nx+m*nx+i].i*Fop[l*nw*nxs+m*nxs+ix].r +
Refl[k*nw*nx+m*nx+i].r*Fop[l*nw*nxs+m*nxs+ix].i;
}
}
/* transfrom result back to time domain */
cr1fft(sum, rtrace, ntfft, 1);
/* place result at source position ixsrc; dx = receiver distance */
for (j = 0; j < nts; j++)
iRN[l*size+ixsrc*nts+j] += rtrace[j]*scl*dx;
} /* end of parallel Nfoc loop */
free(sum);
free(rtrace);
#ifdef _OPENMP
#pragma omp single
npe = omp_get_num_threads();
#endif
} /* end of parallel region */
if (verbose>4) vmess("*** Shot gather %d processed ***", k);
} /* end of nshots (k) loop */
} /* end of if reci */
/* if reciprocal traces are enabled start a new loop over reciprocal shot positions */
if (reci != 0) {
for (k=0; k<nxs; k++) {
if (isxcount[k] == 0) continue;
ixsrc = k;
inx = isxcount[ixsrc]; /* number of traces per reciprocal shot */
#pragma omp parallel default(none) \
shared(iRN, dx, nw, verbose) \
shared(Refl, Nfoc, reci, xrcv, xsrc, xsyn, fxsb, fxse, nxs, dxs) \
shared(nx, dxsrc, inx, k, nfreq, nw_low, nw_high) \
shared(reci_xrcv, reci_xsrc, ixmask) \
shared(Fop, size, nts, ntfft, scl, ixrcv, ixsrc) \
private(l, ix, j, m, i, sum, rtrace, ik, il)
{ /* start of parallel region */
sum = (complex *)malloc(nfreq*sizeof(complex));
rtrace = (float *)calloc(ntfft,sizeof(float));
#pragma omp for schedule(guided,1)
for (l = 0; l < Nfoc; l++) {
/* compute integral over (reciprocal) source positions */
/* multiply R with Fop and sum over nx */
memset(&sum[0].r,0,nfreq*2*sizeof(float));
for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
for (i = 0; i < inx; i++) {
il = reci_xrcv[ixsrc*nxs+i];
ik = reci_xsrc[ixsrc*nxs+i];
ix = NINT((xsrc[il] - fxsb)/dxs);
sum[j].r += Refl[il*nw*nx+m*nx+ik].r*Fop[l*nw*nxs+m*nxs+ix].r -
Refl[il*nw*nx+m*nx+ik].i*Fop[l*nw*nxs+m*nxs+ix].i;
sum[j].i += Refl[il*nw*nx+m*nx+ik].i*Fop[l*nw*nxs+m*nxs+ix].r +
Refl[il*nw*nx+m*nx+ik].r*Fop[l*nw*nxs+m*nxs+ix].i;
}
}
/* transfrom result back to time domain */
cr1fft(sum, rtrace, ntfft, 1);
/* place result at source position ixsrc; dxsrc = shot distance */
for (j = 0; j < nts; j++)
iRN[l*size+ixsrc*nts+j] = ixmask[ixsrc]*(iRN[l*size+ixsrc*nts+j]+rtrace[j]*scl*dxsrc);
} /* end of parallel Nfoc loop */
free(sum);
free(rtrace);
} /* end of parallel region */
} /* end of reciprocal shots (k) loop */
} /* end of if reci */
t = wallclock_time() - t0;
if (verbose) {
vmess("OMP: parallel region = %f seconds (%d threads)", t, npe);
}
return;
}
void synthesisPosistions(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 nshots, int *ixpos, int *npos, int *isxcount, int countmin, int reci, int verbose)
{
int i, j, l, ixsrc, ixrcv, dosrc, k, *count;
float x0, x1;
count = (int *)calloc(nxs,sizeof(int)); // number of traces that contribute to the integration over x
/*================ SYNTHESIS ================*/
for (l = 0; l < 1; l++) { /* assuming all focal operators cover the same lateral area */
// for (l = 0; l < Nfoc; l++) {
*npos=0;
if (reci == 0 || reci == 1) {
for (k=0; k<nshots; k++) {
ixsrc = NINT((xsrc[k] - fxsb)/dxs);
if (verbose>=3) {
vmess("source position: %.2f in operator %d", xsrc[k], ixsrc);
vmess("receiver positions: %.2f <--> %.2f", xrcv[k*nx+0], xrcv[k*nx+nx-1]);
vmess("focal point positions: %.2f <--> %.2f", fxsb, fxse);
}
if ((NINT(xsrc[k]-fxse) > 0) || (NINT(xrcv[k*nx+nx-1]-fxse) > 0) ||
(NINT(xrcv[k*nx+nx-1]-fxsb) < 0) || (NINT(xsrc[k]-fxsb) < 0) ||
(NINT(xrcv[k*nx+0]-fxsb) < 0) || (NINT(xrcv[k*nx+0]-fxse) > 0) ) {
vwarn("source/receiver positions are outside synthesis aperture");
vmess("xsrc = %.2f xrcv_1 = %.2f xrvc_N = %.2f", xsrc[k], xrcv[k*nx+0], xrcv[k*nx+nx-1]);
vmess("source position: %.2f in operator %d", xsrc[k], ixsrc);
vmess("receiver positions: %.2f <--> %.2f", xrcv[k*nx+0], xrcv[k*nx+nx-1]);
vmess("focal point positions: %.2f <--> %.2f", fxsb, fxse);
}
if ( (xsrc[k] >= 0.999*fxsb) && (xsrc[k] <= 1.001*fxse) ) {
j = linearsearch(ixpos, *npos, ixsrc);
if (j < *npos) { /* the position (at j) is already included */
count[j] += xnx[k];
}
else { /* add new postion */
ixpos[*npos]=ixsrc;
count[*npos] += xnx[k];
*npos += 1;
}
// vmess("source position %d is inside synthesis model %f *npos=%d count=%d", k, xsrc[k], *npos, count[*npos]);
}
} /* end of nshots (k) loop */
} /* end of reci branch */
/* if reci=1 or reci=2 source-receive reciprocity is used and new (reciprocal-)sources are added */
if (reci != 0) {
for (k=0; k<nxs; k++) { /* check count in total number of shots added by reciprocity */
if (isxcount[k] >= countmin) {
j = linearsearch(ixpos, *npos, k);
if (j < *npos) { /* the position (at j) is already included */
count[j] += isxcount[k];
}
else { /* add new postion */
ixpos[*npos]=k;
count[*npos] += isxcount[k];
*npos += 1;
}
}
else {
isxcount[k] = 0;
}
}
} /* end of reci branch */
} /* end of Nfoc loop */
if (verbose>=4) {
for (j=0; j < *npos; j++) {
vmess("ixpos[%d] = %d count=%d", j, ixpos[j], count[j]);
}
}
free(count);
/* sort ixpos into increasing values */
qsort(ixpos, *npos, sizeof(int), compareInt);
return;
}
int linearsearch(int *array, size_t N, int value)
{
int j;
/* Check is position is already in array */
j = 0;
while (j < N && value != array[j]) {
j++;
}
return j;
}
/*
void update(float *field, float *term, int Nfoc, int nx, int nt, int reverse, int ixpos)
{
int i, j, l, ix;
if (reverse) {
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) {
j = 0;
Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+nts-j];
}
}
}
}
else {
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) {
j = 0;
Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+nts-j];
}
}
}
}
return;
}
*/
marchenko/marchenko_primaries.c
+ 243
− 569
View file @ 983533c0
This diff is collapsed.
marchenko/synthesis.c 0 → 100644
+ 351
− 0
View file @ 983533c0
/*
* Copyright (c) 2017 by the Society of Exploration Geophysicists.
* 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.
*/
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <genfft.h>
int omp_get_max_threads(void);
int omp_get_num_threads(void);
void omp_set_num_threads(int num_threads);
#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 compareInt(const void *a, const void *b)
{ return (*(int *)a-*(int *)b); }
#ifndef COMPLEX
typedef struct _complexStruct { /* complex number */
float r,i;
} complex;
#endif/* complex */
double wallclock_time(void);
void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, 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
nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int npos, double *tfft, int *isxcount, int
*reci_xsrc, int *reci_xrcv, float *ixmask, int verbose);
void synthesisPositions(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 nshots, int *ixpos, int *npos, int *isxcount, int countmin, int reci, int verbose);
int linearsearch(int *array, size_t N, int value);
/*================ Convolution and Integration ================*/
/* 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 iRN 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 *iRN, 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
nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int npos, double *tfft, int *isxcount, int
*reci_xsrc, int *reci_xrcv, float *ixmask, int verbose)
{
int nfreq, size, inx;
float scl;
int i, j, l, m, iw, ix, k, ixsrc, il, ik;
float *rtrace, idxs;
complex *sum, *ctrace;
int npe;
static int first=1, *ixrcv;
static double t0, t1, t;
size = nxs*nts;
nfreq = ntfft/2+1;
/* scale factor 1/N for backward FFT,
* scale dt for correlation/convolution along time,
* scale dx (or dxsrc) for integration over receiver (or shot) coordinates */
scl = 1.0*dt/((float)ntfft);
#ifdef _OPENMP
npe = omp_get_max_threads();
/* parallelisation is over number of shot positions (nshots) */
if (npe > nshots) {
vmess("Number of OpenMP threads set to %d (was %d)", nshots, npe);
omp_set_num_threads(nshots);
}
#endif
t0 = wallclock_time();
/* reset output data to zero */
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
* 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) {
/* transform muted Ni (Top) to frequency domain, input for next iteration */
for (l = 0; l < Nfoc; l++) {
/* set Fop to zero, so new operator can be defined within ixpos points */
memset(&Fop[l*nxs*nw].r, 0, nxs*nw*2*sizeof(float));
for (i = 0; i < npos; i++) {
rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
ix = ixpos[i];
for (iw=0; iw<nw; iw++) {
Fop[l*nxs*nw+iw*nxs+ix].r = ctrace[nw_low+iw].r;
Fop[l*nxs*nw+iw*nxs+ix].i = mode*ctrace[nw_low+iw].i;
}
}
}
}
else { /* only for first call to synthesis using all nxs traces in G_d */
/* transform G_d to frequency domain, over all nxs traces */
first=0;
for (l = 0; l < Nfoc; l++) {
/* set Fop to zero, so new operator can be defined within all ix points */
memset(&Fop[l*nxs*nw].r, 0, nxs*nw*2*sizeof(float));
for (i = 0; i < nxs; i++) {
rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
for (iw=0; iw<nw; iw++) {
Fop[l*nxs*nw+iw*nxs+i].r = ctrace[nw_low+iw].r;
Fop[l*nxs*nw+iw*nxs+i].i = mode*ctrace[nw_low+iw].i;
}
}
}
idxs = 1.0/dxs;
ixrcv = (int *)malloc(nshots*nx*sizeof(int));
for (k=0; k<nshots; k++) {
for (i = 0; i < nx; i++) {
ixrcv[k*nx+i] = NINT((xrcv[k*nx+i]-fxsb)*idxs);
}
}
}
free(ctrace);
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(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)
{ /* start of parallel region */
sum = (complex *)malloc(nfreq*sizeof(complex));
rtrace = (float *)calloc(ntfft,sizeof(float));
#pragma omp for schedule(guided,1)
for (k=0; k<nshots; k++) {
if ((xsrc[k] < 0.999*fxsb) || (xsrc[k] > 1.001*fxse)) continue;
ixsrc = NINT((xsrc[k] - fxsb)/dxs);
inx = xnx[k]; /* number of traces per shot */
for (l = 0; l < Nfoc; l++) {
/* compute integral over receiver positions */
/* multiply R with Fop and sum over nx */
memset(&sum[0].r,0,nfreq*2*sizeof(float));
for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
for (i = 0; i < inx; i++) {
ix = ixrcv[k*nx+i];
sum[j].r += Refl[k*nw*nx+m*nx+i].r*Fop[l*nw*nxs+m*nxs+ix].r -
Refl[k*nw*nx+m*nx+i].i*Fop[l*nw*nxs+m*nxs+ix].i;
sum[j].i += Refl[k*nw*nx+m*nx+i].i*Fop[l*nw*nxs+m*nxs+ix].r +
Refl[k*nw*nx+m*nx+i].r*Fop[l*nw*nxs+m*nxs+ix].i;
}
}
/* transfrom result back to time domain */
cr1fft(sum, rtrace, ntfft, 1);
/* place result at source position ixsrc; dx = receiver distance */
for (j = 0; j < nts; j++)
iRN[l*size+ixsrc*nts+j] += rtrace[j]*scl*dx;
} /* end of parallel Nfoc loop */
if (verbose>4) vmess("*** Shot gather %d processed ***", k);
} /* end of nshots (k) loop */
free(sum);
free(rtrace);
} /* end of parallel region */
} /* end of if reci */
/* if reciprocal traces are enabled start a new loop over reciprocal shot positions */
if (reci != 0) {
#pragma omp parallel default(none) \
shared(iRN, dx, nw, verbose) \
shared(Refl, Nfoc, reci, xrcv, 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) \
private(l, ix, j, m, i, k, sum, rtrace, ik, il, ixsrc, inx)
{ /* start of parallel region */
sum = (complex *)malloc(nfreq*sizeof(complex));
rtrace = (float *)calloc(ntfft,sizeof(float));
#pragma omp for schedule(guided,1)
for (k=0; k<nxs; k++) {
if (isxcount[k] == 0) continue;
ixsrc = k;
inx = isxcount[ixsrc]; /* number of traces per reciprocal shot */
for (l = 0; l < Nfoc; l++) {
/* compute integral over (reciprocal) source positions */
/* multiply R with Fop and sum over nx */
memset(&sum[0].r,0,nfreq*2*sizeof(float));
for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
for (i = 0; i < inx; i++) {
il = reci_xrcv[ixsrc*nxs+i];
ik = reci_xsrc[ixsrc*nxs+i];
ix = NINT((xsrc[il] - fxsb)/dxs);
sum[j].r += Refl[il*nw*nx+m*nx+ik].r*Fop[l*nw*nxs+m*nxs+ix].r -
Refl[il*nw*nx+m*nx+ik].i*Fop[l*nw*nxs+m*nxs+ix].i;
sum[j].i += Refl[il*nw*nx+m*nx+ik].i*Fop[l*nw*nxs+m*nxs+ix].r +
Refl[il*nw*nx+m*nx+ik].r*Fop[l*nw*nxs+m*nxs+ix].i;
}
}
/* transfrom result back to time domain */
cr1fft(sum, rtrace, ntfft, 1);
/* place result at source position ixsrc; dxsrc = shot distance */
for (j = 0; j < nts; j++)
iRN[l*size+ixsrc*nts+j] = ixmask[ixsrc]*(iRN[l*size+ixsrc*nts+j]+rtrace[j]*scl*dxsrc);
} /* end of Nfoc loop */
} /* end of parallel reciprocal shots (k) loop */
free(sum);
free(rtrace);
} /* end of parallel region */
} /* end of if reci */
t = wallclock_time() - t0;
if (verbose>2) {
vmess("OMP: parallel region = %f seconds (%d threads)", t, npe);
}
return;
}
void synthesisPositions(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 nshots, int *ixpos, int *npos, int *isxcount, int countmin, int reci, int verbose)
{
int i, j, l, ixsrc, ixrcv, dosrc, k, *count;
float x0, x1;
count = (int *)calloc(nxs,sizeof(int)); // number of traces that contribute to the integration over x
/*================ SYNTHESIS ================*/
/* assuming all focal operators cover the same lateral area */
*npos=0;
if (reci == 0 || reci == 1) {
for (k=0; k<nshots; k++) {
ixsrc = NINT((xsrc[k] - fxsb)/dxs);
if (verbose>=3) {
vmess("source position: %.2f in operator %d", xsrc[k], ixsrc);
vmess("receiver positions: %.2f <--> %.2f", xrcv[k*nx+0], xrcv[k*nx+nx-1]);
vmess("focal point positions: %.2f <--> %.2f", fxsb, fxse);
}
if ((NINT(xsrc[k]-fxse) > 0) || (NINT(xrcv[k*nx+nx-1]-fxse) > 0) ||
(NINT(xrcv[k*nx+nx-1]-fxsb) < 0) || (NINT(xsrc[k]-fxsb) < 0) ||
(NINT(xrcv[k*nx+0]-fxsb) < 0) || (NINT(xrcv[k*nx+0]-fxse) > 0) ) {
vwarn("source/receiver positions are outside synthesis aperture");
vmess("xsrc = %.2f xrcv_1 = %.2f xrvc_N = %.2f", xsrc[k], xrcv[k*nx+0], xrcv[k*nx+nx-1]);
vmess("source position: %.2f in operator %d", xsrc[k], ixsrc);
vmess("receiver positions: %.2f <--> %.2f", xrcv[k*nx+0], xrcv[k*nx+nx-1]);
vmess("focal point positions: %.2f <--> %.2f", fxsb, fxse);
}
//fprintf(stderr,"k=%d xsrc[k]=%f 0.999*fxsb=%f, 1.001*fxse=%f %f %f\n",k, xsrc[k], 0.999*fxsb, 1.001*fxse, fxsb, fxse);
//if ( (xsrc[k] >= 0.999*fxsb) && (xsrc[k] <= 1.001*fxse) ) {
if ( (ixsrc < nxs) && (ixsrc >= 0) ) {
j = linearsearch(ixpos, *npos, ixsrc);
if (j < *npos) { /* the position (at j) is already included */
count[j] += xnx[k];
}
else { /* add new postion */
ixpos[*npos]=ixsrc;
count[*npos] += xnx[k];
*npos += 1;
}
if (verbose>=3) {
vmess("source position %d is inside synthesis model %f *npos=%d count=%d", k, xsrc[k], *npos, count[*npos]);
vmess("ixpos[%d] = %d ixsrc=%d", *npos-1, ixpos[*npos-1], ixsrc);
}
}
else {
if (verbose>=2) {
vwarn("source position %d is outside synthesis model %f ixsrc=%d", k, xsrc[k], ixsrc);
}
}
} /* end of nshots (k) loop */
} /* end of reci branch */
/* if reci=1 or reci=2 source-receive reciprocity is used and new (reciprocal-)sources are added */
if (reci != 0) {
for (k=0; k<nxs; k++) { /* check count in total number of shots added by reciprocity */
if (isxcount[k] >= countmin) {
j = linearsearch(ixpos, *npos, k);
if (j < *npos) { /* the position (at j) is already included */
count[j] += isxcount[k];
}
else { /* add new postion */
ixpos[*npos]=k;
count[*npos] += isxcount[k];
*npos += 1;
}
}
else {
isxcount[k] = 0;
}
}
} /* end of reci branch */
if (verbose>=4) {
for (j=0; j < *npos; j++) {
vmess("ixpos[%d] = %d count=%d", j, ixpos[j], count[j]);
}
}
free(count);
/* sort ixpos into increasing values */
qsort(ixpos, *npos, sizeof(int), compareInt);
return;
}
int linearsearch(int *array, size_t N, int value)
{
int j;
/* Check is position is already in array */
j = 0;
while (j < N && value != array[j]) {
j++;
}
return j;
}
scripts/SolidEarth2018/ScatterModel/ClassicTimeReverseR1.scr 0 → 100755
+ 207
− 0
View file @ 983533c0
#!/bin/bash
#
#SBATCH -J ClassicTimeReverse
#SBATCH --cpus-per-task=12
#SBATCH --ntasks=1
#SBATCH --time=1:00:00
#SBATCH --hint=nomultithread
export PATH=:$HOME/src/OpenSource/bin:$HOME/bin:$PATH:
which makewave
which makemod
which fdelmodc
cd /vardim/home/thorbcke/data/Kees/Marchenko/SolidEarth/ScatterModel
#makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
#makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=4096 t0=0.3
dx=2.5
dt=0.0005
depth=850
ix1a=1
ix1b=$(echo "scale=0; ${ix1a}+6000/${dx}" | bc -l)
base=$(echo "scale=0; ${depth}/${dx}" | bc -l)
makewave fp=25 dt=$dt file_out=wave.su nt=1024 t0=0.1 scale=1
file_mod=scat
export OMP_NUM_THREADS=12
#forward model of scattered response of source at depth
#fdelmodc \
file_cp=hom_cp.su ischeme=1 iorder=4 \
file_den=scatter_ro.su \
file_src=wave.su \
file_rcv=ctr.su \
src_type=1 \
src_injectionrate=1 \
src_orient=1 \
rec_type_vz=0 \
rec_type_p=1 \
rec_int_vz=2 \
dtrcv=$dt \
rec_delay=0.1 \
verbose=2 \
tmod=4.1955 \
dxrcv=$dx \
xrcv1=-2250 xrcv2=2250 \
zrcv1=0 zrcv2=0 \
xsrc=0 zsrc=$depth \
npml=250 \
left=2 right=2 top=2 bottom=2
#impulse response through scattered medium
#fdelmodc \
file_cp=hom_cp.su ischeme=1 iorder=4 \
file_den=scatter_ro.su \
file_src=wave.su \
file_rcv=ctr_impulse_response.su \
src_type=1 \
src_injectionrate=1 \
src_orient=1 \
rec_type_vz=0 \
rec_type_p=1 \
rec_int_vz=2 \
rec_delay=0.1 \
tmod=4.1940 \
dtrcv=0.004 \
verbose=2 \
dxrcv=10.0 \
xrcv1=-2250 xrcv2=2250 \
zrcv1=$depth zrcv2=$depth \
xsrc=$xsrc zsrc=0 \
ntaper=250 \
left=2 right=2 top=2 bottom=2
#Forward model of homogenoeus response of source at depth
#fdelmodc \
file_cp=hom_cp.su ischeme=1 iorder=4 \
file_den=hom_ro.su \
file_src=wave.su \
file_rcv=shom.su \
src_type=1 \
src_injectionrate=1 \
src_orient=1 \
rec_type_vz=0 \
rec_type_p=1 \
rec_int_vz=2 \
dtrcv=$dt \
rec_delay=0.1 \
verbose=2 \
tmod=4.1955 \
dxrcv=$dx \
xrcv1=-2250 xrcv2=2250 \
zrcv1=0 zrcv2=0 \
xsrc=0 zsrc=$depth \
npml=250 \
left=2 right=2 top=2 bottom=2
app=300
app=2250
filectr=ctr_${app}.su
fileshom=shom_${app}.su
suwind < ctr_rp.su key=gx min=-${app}000 max=${app}000 > $filectr
suwind < shom_rp.su key=gx min=-${app}000 max=${app}000 > $fileshom
suwind key=gx j=100000 itmax=801 min=-${app}000 max=${app}000 < ctr_impulse_response_rp.su | \
supswigp wbox=2 hbox=4 titlesize=-1 labelsize=8 f2=-300 d2=100 f2num=-300 d2num=100 d1num=0.5 axescolor=black > impulse_response_app${app}.eps
suwind key=gx j=100000 itmax=801 min=-${app}000 max=${app}000 < ctr_impulse_response_rp.su | \
supswigp wbox=2 hbox=4 titlesize=-1 labelsize=4 frame=0 axescolor=black > impulse_response_app${app}_noaxis.eps
#Time reverse of scattered field through scattered medium
fdelmodc \
file_cp=hom_cp.su ischeme=1 iorder=4 \
file_den=scatter_ro.su \
file_src=$filectr \
grid_dir=1 \
src_type=1 \
src_injectionrate=1 \
src_orient=1 \
rec_delay=0.0 \
verbose=2 \
tmod=4.51 \
file_snap=${file_mod}_timerev_scat_scat${app}.su \
tsnap1=3.4950 dtsnap=0.05 tsnap2=4.5005 sna_type_vz=0 \
dxsnap=$dx dzsnap=$dx zsnap1=0 zsnap2=1200 xsnap1=-2250 xsnap2=2250 \
npml=250 \
left=2 right=2 top=2 bottom=2
#Time reverse of homogenoeus field through homogenoeus medium
fdelmodc \
file_cp=hom_cp.su ischeme=1 iorder=4 \
file_den=hom_ro.su \
file_src=$fileshom \
grid_dir=1 \
src_type=1 \
src_injectionrate=1 \
src_orient=1 \
rec_delay=0.0 \
verbose=2 \
tmod=4.51 \
file_snap=${file_mod}_timerev_hom_hom${app}.su \
tsnap1=3.4950 dtsnap=0.05 tsnap2=4.5005 sna_type_vz=0 \
dxsnap=$dx dzsnap=$dx zsnap1=0 zsnap2=1200 xsnap1=-2250 xsnap2=2250 \
npml=250 \
left=2 right=2 top=2 bottom=2
# curve=curve1,curve2,curve3 npair=25,25,25 curvecolor=black,black,black curvedash=3,3,3 \
#sumax < ${file_mod}_back_hom_ctr${app}_sp.su mode=abs outpar=nep
#x2end=`cat nep | awk '{print $1}'`
#echo $x2end
for file in ${file_mod}_timerev_hom_hom${app} ${file_mod}_timerev_scat_scat${app}
do
sumax < ${file}_sp.su mode=abs outpar=nep
clip=`cat nep | awk '{print $1/7}'`
echo $file has clip $clip
for fldr in 10 13 16
do
times=$(echo "scale=2; 0.05*(13-${fldr})" | bc -l)
atime=`printf "%4.2f" $times`
suwind key=fldr min=$fldr max=$fldr < ${file}_sp.su | \
supsimage hbox=4 wbox=6.7 labelsize=10 \
x1beg=0 x1end=1200 clip=$clip \
n1tic=4 x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file}_$atime.eps
done
suwind itmin=$base itmax=$base key=fldr min=13 max=13 < ${file}_sp.su | sumax mode=abs outpar=nep
scl=`cat nep | awk '{print 1.0/$1}'`
echo scale for trace = $scl
suwind itmin=$base itmax=$base key=fldr min=13 max=13 < ${file}_sp.su | sugain scale=$scl | sustrip > trace.bin
suaddhead < trace.bin n1=1801 dt=$dx | supsgraph hbox=2 wbox=6 labelsize=10 \
f1=-2250 d1=$dx x1beg=-500 x1end=500 f1num=-500 d1num=500 style=normal > ${file}_z${depth}_t0.eps
suaddhead < trace.bin n1=1801 dt=$dx > ${file}_z${depth}_t0.su
(( imin = base - 50 ))
(( imax = base + 50 ))
echo $base $imin $imax
suwind key=fldr min=13 max=13 < ${file}_sp.su | \
suwind itmin=$imin itmax=$imax key=gx min=-125000 max=125000 | \
sustrip > ${file}_t0.bin
python3 readbin.py ${file}_t0.bin
done
(cat ${file_mod}_timerev_hom_hom${app}_z${depth}_t0.su; cat ${file_mod}_timerev_scat_scat${app}_z${depth}_t0.su ) | \
supsgraph hbox=2 wbox=6 labelsize=10 \
f1=-2250 d1=$dx x1beg=-1000 x1end=1000 f1num=-1000 d1num=500 x2beg=-0.1 \
style=normal linecolor=red,blue,green > ${file_mod}_timerev_z${depth}_t0.eps
rm nep trace.bin
exit;
xgraph < trace.bin n=451 pairs=2 d1=10 title=hom
suwind itmin=75 itmax=75 key=fldr min=13 max=13 < snap_back_ctr_sp.su | sustrip > trace.bin
xgraph < trace.bin n=451 pairs=2 d1=10 title=scatter
scripts/SolidEarth2018/ScatterModel/backpropf1plus.scr 0 → 100755
+ 88
− 0
View file @ 983533c0
#!/bin/bash
#PBS -N fdelmod
#PBS -q verylong
#PBS -l nodes=1
#PBS -k eo
#PBS -j eo
export PATH=:$HOME/src/OpenSource/utils:$HOME/bin64:$PATH:
which makewave
which makemod
tmpdir=/tmp/shotS
cd /home/thorbcke/data/Kees/MultElim/ScatterModel/Redatum
#makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
#makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=4096 t0=0.3
dx=2.5
dt=0.0005
makemod sizex=5000 sizez=1100 dx=$dx dz=$dx cp0=1900 ro0=1200 \
orig=-2500,0 file_base=scatterf1.su verbose=2 \
intt=randdf x=-2500 z=200 cp=2700 ro=4700 var=3000,5 \
intt=def x=-2500,0,2500 z=700,700,700 cp=1900 ro=1200
# intt=diffr x=-90,0,90 z=900,900,900 cp=2300,2300,2300 ro=2800,2800,2800 var=5 \
# verbose=4
makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
export OMP_NUM_THREADS=2
interpolate file_in=f1plus0.su d1out=$dt file_out=f1plusdt.su
$HOME/bin64/fdelmodc \
file_cp=scatterf1_cp.su ischeme=1 iorder=4 \
file_den=scatterf1_ro.su \
file_src=f1plusdt.su \
file_rcv=backprop_f1plusz800.su \
grid_dir=0 \
src_type=1 \
src_injectionrate=1 \
src_orient=1 \
rec_type_vz=0 \
rec_type_p=1 \
rec_int_vz=2 \
dtrcv=0.002 \
rec_delay=0.0 \
verbose=2 \
tmod=4.094 \
dxrcv=10.0 \
xrcv1=-2250 xrcv2=2250 \
zrcv1=800 zrcv2=800 \
xsrc=0 zsrc=0 \
ntaper=250 \
file_snap=back.su tsnap1=1.5005 dtsnap=0.05 tsnap2=2.5005 dxsnap=10 dzsnap=10 zsnap1=0 zsnap2=2000 xsnap1=-2250 xsnap2=2250 \
left=4 right=4 top=4 bottom=4
suwind key=sx min=0 max=0 itmax=1023 < p0plusall.su > p0plus.su
basop file_in=p0plus.su file_out=p0plusr.su choice=5 verbose=1
rotate nrot=512 < p0plusr.su > f1plusi0.su
interpolate file_in=f1plusi0.su d1out=$dt file_out=f1plusi0dt.su
$HOME/bin64/fdelmodc \
file_cp=scatterf1_cp.su ischeme=1 iorder=4 \
file_den=scatterf1_ro.su \
file_src=f1plusi0dt.su \
file_rcv=backprop_f1plusi0z800.su \
grid_dir=0 \
src_type=1 \
src_injectionrate=1 \
src_orient=1 \
rec_type_vz=0 \
rec_type_p=1 \
rec_int_vz=2 \
dtrcv=0.002 \
rec_delay=0.0 \
verbose=2 \
tmod=4.094 \
dxrcv=10.0 \
xrcv1=-2250 xrcv2=2250 \
zrcv1=800 zrcv2=800 \
xsrc=0 zsrc=0 \
ntaper=250 \
file_snap=backi0.su tsnap1=1.5005 dtsnap=0.05 tsnap2=2.5005 dxsnap=10 dzsnap=10 zsnap1=0 zsnap2=2000 xsnap1=-2250 xsnap2=2250 \
left=4 right=4 top=4 bottom=4
scripts/SolidEarth2018/ScatterModel/backpropf2.scr 0 → 100755
+ 81
− 0
View file @ 983533c0
#!/bin/bash
export PATH=$HOME/OpenSource/bin/:$PATH:
cd /vardim/home/thorbcke/data/Kees/Marchenko/SolidEarth/ScatterModel
dx=2.5
dt=0.0005
file_cp=scatter_cp.su
file_ro=scatter_ro.su
export OMP_NUM_THREADS=8
# t=0 focal time is at 2.0445 seconds back=propagating (dtrcv*(ns/2-1)+dt)
# shift f2.su such that t=0 is positioned in the middle of the time axis
# the extra shift of 0.000250 is needed because of the staggered time implementation of the Finite Difference program.
ns=`surange <f2.su | grep ns | awk '{print $2}'`
dtrcv=`surange < f2.su | grep dt | awk '{print $2/1000000.0}'`
suwind key=gx min=-2250000 max=2250000 itmax=1023 < f2.su > nep.su
shift=$(echo "scale=6; ($dtrcv*(($ns)/2.0-1)+$dt)" | bc -l)
echo $shift
basop choice=shift shift=$shift file_in=nep.su verbose=1 > pplus.su
# the f2.su is sampled with 4ms the FD program needs 0.5ms
# time axis is interpolated by making use of FFT's: sinc interpolation
#ftr1d file_in=pplus.su file_out=freq.su
#sushw <freq.su key=nhs,dt a=8192,500 >fr.su
#ftr1d file_in=fr.su n1=8194 file_out=pplusdt.su verbose=1
midsnap=2.04
app=300
suwind < pplus.su key=gx min=-${app}000 max=${app}000 > pplus_${app}.su
#backpropagate f2.su and collect snapshots
fdelmodc \
file_cp=$file_cp ischeme=1 iorder=4 \
file_den=$file_ro \
file_src=pplus_${app}.su \
file_rcv=backprop_f2_z850_${app}.su \
dt=$dt \
grid_dir=0 \
src_type=1 \
src_injectionrate=1 \
src_orient=1 \
rec_type_vz=0 \
rec_type_p=1 \
rec_int_vz=2 \
dtrcv=0.004 \
rec_delay=0.0 \
verbose=2 \
tmod=3.10 \
dxrcv=10.0 \
xrcv1=-3000 xrcv2=3000 \
zrcv1=1100 zrcv2=1100 \
zsrc=0 xsrc=0 \
npml=101 \
file_snap=backpropmar_${app}.su tsnap1=`echo "$midsnap-0.1" | bc -l` dtsnap=0.01 tsnap2=`echo "$midsnap+0.1" | bc -l` dxsnap=5 dzsnap=5 zsnap1=0 zsnap2=2000 xsnap1=-2250 xsnap2=2250 \
sna_type_vz=0 \
sna_type_p=1 \
left=2 right=2 top=2 bottom=2
for file in backpropmar_${app}
do
sumax < ${file}_sp.su mode=abs outpar=nep
clip=`cat nep | awk '{print $1/7}'`
echo $file has clip $clip
for fldr in 10 11 12
do
times=$(echo "scale=2; 0.01*(11-${fldr})" | bc -l)
atime=`printf "%4.2f" $times`
suwind key=fldr min=$fldr max=$fldr < ${file}_sp.su | \
supsimage hbox=4 wbox=6.7 labelsize=10 \
x1beg=0 x1end=1200 clip=$clip \
n1tic=4 x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file}_$atime.eps
done
done
scripts/SolidEarth2018/ScatterModel/corr.scr
+ 9
− 6
View file @ 983533c0
......@@ -29,7 +29,7 @@ sudipfilt < BwithTraceA.su slopes=-0.0007,-0.0001,0,0.0001,0.0007 amps=0,1,1,1,0
#optimal results obtained by placing recievers at level of TraceA (virtual shot position) in time-reversal computation
#fdelmodc \
fdelmodc \
file_cp=hom_cp.su ischeme=1 iorder=4 \
file_den=scatter_ro.su \
file_src=ctr_rp.su \
......@@ -60,10 +60,13 @@ sudipfilt < SI_reference_timerev_rp.su slopes=-0.0007,-0.0001,0,0.0001,0.0007 am
#plot noisetraces
for offset in -1000 -500 0 500 1000
do
suwind key=offset min=$offset max=$offset < noiseTraces.su | \
supswigp wbox=1 hbox=5 titlesize=-1 labelsize=-1 axescolor=white > noiseTraces_off${offset}.eps
done
#for offset in -1000 -500 0 500 1000
#do
#suwind key=offset min=$offset max=$offset < noiseTraces.su | \
# supswigp wbox=1 hbox=5 titlesize=-1 labelsize=-1 axescolor=white > noiseTraces_off${offset}.eps
#done
suwind key=gx j=100000 min=-1000000 max=1000000 < noiseTraces.su | \
supswigp wbox=6.7 hbox=4 titlesize=-1 labelsize=8 f2=-1000 d2=100 f2num=-1000 d2num=500 \
x1beg=0 x1end=1.0 axescolor=black > noiseTraces.eps
scripts/SolidEarth2018/ScatterModel/direct.scr 0 → 100755
+ 36
− 0
View file @ 983533c0
#!/bin/bash
cd /vardim/home/thorbcke/data/Kees/Marchenko/SolidEarth/ScatterModel
dx=2.5
dt=0.0005
makemod sizex=10000 sizez=500 dx=$dx dz=$dx cp0=1900 ro0=1200 \
orig=-5000,0 file_base=homL verbose=2
makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=4096 t0=0.3 scale=0
export OMP_NUM_THREADS=8
fdelmodc \
file_cp=homL_cp.su ischeme=1 iorder=4 \
file_den=homL_ro.su \
file_src=wavefw.su \
file_rcv=direct.su \
src_type=7 \
src_injectionrate=1 \
src_orient=1 \
rec_type_vz=0 \
rec_type_p=1 \
rec_int_vz=2 \
dtrcv=0.004 \
rec_delay=0.3 \
verbose=2 \
tmod=4.394 \
dxrcv=10.0 \
xrcv1=-4500 xrcv2=4500 \
zrcv1=0 zrcv2=0 \
xsrc=0 zsrc=0 \
npml=250 \
left=2 right=2 top=2 bottom=2
scripts/SolidEarth2018/ScatterModel/marchenko.scr 0 → 100755
+ 55
− 0
View file @ 983533c0
#!/bin/bash -x
export PATH=$HOME/src/OpenSource/bin:$PATH:
export OMP_NUM_THREADS=1
depth=850
cd /vardim/home/thorbcke/data/Kees/Marchenko/SolidEarth/ScatterModel
#compute Td
#fdelmodc \
file_cp=hom_cp.su ischeme=1 iorder=4 \
file_den=hom_ro.su \
file_src=wave.su \
file_rcv=Td.su \
src_type=1 \
src_injectionrate=1 \
src_orient=1 \
rec_type_vz=0 \
rec_type_p=1 \
rec_int_vz=2 \
dtrcv=0.004 \
rec_delay=0.1 \
verbose=2 \
tmod=4.1920 \
dxrcv=10 \
xrcv1=-2250 xrcv2=2250 \
zrcv1=0 zrcv2=0 \
xsrc=0 zsrc=$depth \
npml=250 \
left=2 right=2 top=2 bottom=2
#apply the Marchenko algorithm
marchenko file_shot=shotsRefl/refl_rp.su file_tinv=Td_rp.su nshots=451 verbose=2 \
tap=0 niter=15 hw=8 shift=7 smooth=3 \
file_green=pgreen.su file_gplus=Gplus.su file_gmin=Gmin.su \
file_f1plus=f1plus.su file_f1min=f1min.su file_f2=f2.su
exit;
#compare Green's funtions on Marhcenko 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}'`
rf=`cat neprf | awk '{print $1}'`
value=${value/[eE][+][0]/*10^}
mg=${mg/[eE][+][0]/*10^}
rf=${rf/[eE][+][0]/*10^}
rm nep*
scale=$(echo "scale=3; ($rf)/($mg)" | bc -l)
echo $scale
(suwind key=gx min=0 max=0 itmax=511 < pgreen.su | sugain scale=$scale; \
suwind key=gx min=0 max=0 < referenceP_rp.su) | suxgraph
scripts/SolidEarth2018/ScatterModel/model.scr
+ 21
− 12
View file @ 983533c0
......@@ -16,30 +16,39 @@ dt=0.0005
#shots=var=6000,5
#high contrast model
makemod sizex=5000 sizez=1200 dx=$dx dz=$dx cp0=1900 ro0=1200 \
orig=-2500,0 file_base=scatter.su \
intt=randdf x=-2500 z=200 cp=1900 ro=4700 var=4001,11 \
intt=def x=-2500,0,2500 z=700,700,700 cp=1900 ro=1200 \
verbose=4
#low contrast model
#makemod sizex=5000 sizez=1200 dx=$dx dz=$dx cp0=1900 ro0=1200 \
# orig=-2500,0 file_base=scatter.su \
# intt=randdf x=-2500 z=200 cp=1900 ro=2400 var=4001,11 \
# intt=randdf x=-2500 z=200 cp=1900 ro=4700 var=4001,11 \
# intt=def x=-2500,0,2500 z=700,700,700 cp=1900 ro=1200 \
# verbose=4
#low contrast model for Revision 1 in paper
makemod sizex=5000 sizez=1200 dx=$dx dz=$dx cp0=1900 ro0=1200 \
orig=-2500,0 file_base=scatter.su \
intt=randdf x=-2500 z=200 cp=1900 ro=2400 var=4001,11 \
intt=def x=-2500,0,2500 z=700,700,700 cp=1900 ro=1200 \
verbose=4
makemod sizex=5000 sizez=1200 dx=$dx dz=$dx cp0=1900 ro0=1200 \
orig=-2500,0 file_base=hom verbose=2
makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=4096 t0=0.3 scale=0
supsimage hbox=4 wbox=6 labelsize=10 < scatter_ro.su\
x1beg=0 x1end=1000.0 legend=0 threecolor=1 \
#supsimage hbox=4 wbox=6 labelsize=10 < scatter_ro.su\
# x1beg=0 x1end=1200.0 legend=0 threecolor=1 \
# d1s=1.0 d2s=1.0 \
# wrgb=1.0,.5,0 grgb=0,.7,1.0 brgb=0,1.0,0 \
# bps=24 bclip=2400 wclip=1500 blockinterp=1 legend=1 \
# n1tic=5 x2beg=-2500 f2num=-2000 d2num=1000 x2end=2500 > modelscatter_ro.eps
#for R1 on scale with snapshots
supsimage hbox=4 wbox=6.7 labelsize=10 < scatter_ro.su\
legend=0 threecolor=1 \
d1s=1.0 d2s=1.0 \
wrgb=1.0,.5,0 grgb=0,.7,1.0 brgb=0,1.0,0 \
bps=24 bclip=2400 wclip=1500 blockinterp=1 legend=1 \
n1tic=5 x2beg=-2500 f2num=-2000 d2num=1000 x2end=2500 > modelscatter_ro.eps
x1beg=0 x1end=1200 \
bps=24 bclip=2400 wclip=1500 legend=1 \
n1tic=4 x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > modelscatter_ro.eps
exit;
......
scripts/SolidEarth2018/ScatterModel/remove_direct.scr 0 → 100755
+ 51
− 0
View file @ 983533c0
#!/bin/bash
#PBS -N fdelmod
#PBS -q verylong
#PBS -l nodes=1
#PBS -k eo
#PBS -j eo
export PATH=$HOME/bin:$HOME/src/OpenSource/bin:$PATH:
cd /vardim/home/thorbcke/data/Kees/Marchenko/SolidEarth/ScatterModel
dxshot=10
ishot=0
nshots=451
file_R=shotsRefl/refl_rp.su
rm $file_R
while (( ishot < nshots ))
do
(( xsrc = -2250 + ${ishot}*${dxshot} ))
(( sx = ${xsrc}*1000 ))
(( iishot = ${ishot}*${dxshot}/10 ))
(( tr1 = 451 - ${iishot} ))
(( tr2 = ${tr1} + 450 ))
echo xsrc=$xsrc tr1=$tr1 tr2=$tr2
# direct wave
suwind < direct_rp.su key=tracl min=$tr1 max=$tr2 > direct.su
# 2D shot
file_rcv=shotsRefl/shot_${xsrc}_rp.su
#suwind key=tracl min=1 max=601 < $file_rcv > shotz0.su
sudiff $file_rcv direct.su | suwind itmax=1023 > refl.su
# 1D shot
# suwind < shot1d${fast}_rp.su key=tracl min=$tr1 max=$tr2 > shotz0.su
# sudiff shotz0.su direct.su > refl.su
(( ishot = $ishot + 1))
sushw < refl.su key=fldr a=$ishot >> $file_R
done
#rm shotz0.su direct.su refl.su
#make sure nt is number for FFT
#convert file_in=shots/refl_rp.su saminc=4 | suwind itmax=1539 > shots/refl_4ms.su
scripts/SolidEarth2018/ScatterModel/shotsRefl.scr 0 → 100755
+ 69
− 0
View file @ 983533c0
#!/bin/bash
#
export PATH=:$HOME/src/OpenSource/bin:$HOME/bin:$PATH:
cd /vardim/home/thorbcke/data/Kees/Marchenko/SolidEarth/ScatterModel
dt=0.0005
makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=4096 t0=0.3 scale=0
mkdir -p shotsRefl
mkdir -p jobs
zsrc=0
dxshot=10
ishot=0
nshots=451
while (( ishot < nshots ))
do
(( xsrc = -2250 + ${ishot}*${dxshot} ))
echo ishot=$ishot xsrc=$xsrc zsrc=$zsrc
cat << EOF > jobs/job_$ishot.job
#!/bin/bash
#
#SBATCH -J scat_${xsrc}
#SBATCH --cpus-per-task=8
#SBATCH --ntasks=1
#SBATCH --time=0:10:00
cd \$SLURM_SUBMIT_DIR
export PATH=:\$HOME/src/OpenSource/bin:\$HOME/bin:\$PATH:
export OMP_NUM_THREADS=8
file_rcv=shotsRefl/shot_${xsrc}.su
fdelmodc \
file_cp=scatter_cp.su ischeme=1 iorder=4 \
file_den=scatter_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=-2250 xrcv2=2250 \
zrcv1=0 zrcv2=0 \
xsrc=$xsrc zsrc=0 \
ntaper=250 \
left=2 right=2 top=2 bottom=2
EOF
sbatch jobs/job_$ishot.job
(( ishot = $ishot + 1))
done
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