diff --git a/extrap3d/lib/dataFileIO.c b/extrap3d/lib/dataFileIO.c
index 84b910c2df5d78198459de67cfc679fb188388d3..4612557c0031fadc4aba670f5ae7465f44e33813 100644
--- a/extrap3d/lib/dataFileIO.c
+++ b/extrap3d/lib/dataFileIO.c
@@ -107,7 +107,7 @@ int read_FFT_DataFile(FILE *fp, complex *data, Area vel_area, int nfft, int nw,
}
-int write_FFT_DataFile(FILE *fp, complex *data, Area data_area, int fldr, int nt, int nfft, int nw, int nw_low, float dt, int out_su, int conjg, int verbose)
+int write_FFT_DataFile(FILE *fp, complex *data, Area data_area, int fldr, int nt, int nfft, int nw, int nw_low, float dt, int out_su, int conjg, int freq, int verbose)
{
int ix, iy, iw, i, nx, ny, sxy, nfreq, pos, sign;
int xvmin, yvmin;
@@ -142,37 +142,72 @@ int write_FFT_DataFile(FILE *fp, complex *data, Area data_area, int fldr, int nt
hdr[0].d2 = dx;
hdr[0].scalco = -1000;
- for (iy=0; iy<ny; iy++) {
- for (ix=0; ix<nx; ix++) {
- pos = iy*nx+ix;
- for (i=0; i<nfreq; i++)
- ctrace[i].r = ctrace[i].i = 0.0;
- for (iw=0; iw<nw; iw++) {
- ctrace[iw+nw_low].r = data[iw*sxy+pos].r;
- ctrace[iw+nw_low].i = conjg*data[iw*sxy+pos].i;
- }
- /* transform back to time */
- cr1fft(ctrace, trace, nfft, sign);
-
- /* write to output file */
- if (out_su) {
- hdr[0].gx = (int)(xvmin+ix*dx)*1000;
- hdr[0].gy = (int)(yvmin+iy*dy)*1000;
- hdr[0].f2 = xvmin+ix*dx;
- hdr[0].tracf = ix+1;
- hdr[0].tracl = iy*nx+ix+1;
- nwrite = fwrite(hdr, 1, TRCBYTES, fp);
- assert( nwrite == TRCBYTES );
- nwrite = fwrite(trace, sizeof(float), nt, fp);
- assert( nwrite == nt );
+ if (freq==0) {
+ for (iy=0; iy<ny; iy++) {
+ for (ix=0; ix<nx; ix++) {
+ pos = iy*nx+ix;
+ for (i=0; i<nfreq; i++)
+ ctrace[i].r = ctrace[i].i = 0.0;
+ for (iw=0; iw<nw; iw++) {
+ ctrace[iw+nw_low].r = data[iw*sxy+pos].r;
+ ctrace[iw+nw_low].i = conjg*data[iw*sxy+pos].i;
+ }
+
+ /* transform back to time */
+ cr1fft(ctrace, trace, nfft, sign);
+
+ /* write to output file */
+ if (out_su) {
+ hdr[0].gx = (int)(xvmin+ix*dx)*1000;
+ hdr[0].gy = (int)(yvmin+iy*dy)*1000;
+ hdr[0].f2 = xvmin+ix*dx;
+ hdr[0].tracf = ix+1;
+ hdr[0].tracl = iy*nx+ix+1;
+ nwrite = fwrite(hdr, 1, TRCBYTES, fp);
+ assert( nwrite == TRCBYTES );
+ nwrite = fwrite(trace, sizeof(float), nt, fp);
+ assert( nwrite == nt );
+ }
+ else {
+ nwrite = fwrite(trace, sizeof(float), nt, fp);
+ assert( nwrite == nt );
+ }
}
- else {
- nwrite = fwrite(trace, sizeof(float), nt, fp);
- assert( nwrite == nt );
+ }
+ }
+ else { /* Frequency output slices */
+ hdr[0].trid = 122;
+ hdr[0].fldr = fldr;
+ hdr[0].trwf = ny;
+ hdr[0].ntr = nw*ny;
+ hdr[0].ns = nx;
+ hdr[0].dt = 1000000*dt;
+ hdr[0].f1 = 0.0;
+ hdr[0].f2 = yvmin;
+
+ for (iw=0; iw<nw; iw++) {
+ for (iy=0; iy<ny; iy++) {
+ /* write to output file */
+ if (out_su) {
+ hdr[0].gx = (int)(xvmin)*1000;
+ hdr[0].gy = (int)(yvmin+iy*dy)*1000;
+ hdr[0].f1 = xvmin+ix*dx;
+ hdr[0].f2 = yvmin+iy*dy;
+ hdr[0].tracf = iy+1;
+ hdr[0].tracl = iw*ny+iy+1;
+ nwrite = fwrite(hdr, 1, TRCBYTES, fp);
+ assert( nwrite == TRCBYTES );
+ nwrite = fwrite(&data[iw*sxy+iy*nx], sizeof(complex), nx, fp);
+ assert( nwrite == nx );
+ }
+ else {
+ nwrite = fwrite(&data[iw*sxy+iy*nx], sizeof(complex), nx, fp);
+ assert( nwrite == nx );
+ }
}
}
- }
+ }
fflush(fp);
free(trace);
diff --git a/extrap3d/main/cfpmod3d.c b/extrap3d/main/cfpmod3d.c
index 6882e550f4100f2f1e766f12d684eb44d4f12682..2eb38afd05035008d486a6657c93dbe943fe49f7 100644
--- a/extrap3d/main/cfpmod3d.c
+++ b/extrap3d/main/cfpmod3d.c
@@ -24,7 +24,7 @@ void cr1fft(complex *cdata, float *rdata, int n, int sign);
/****** IO routines *******/
int openVelocityFile(char *file_vel, FILE **fp, Area *vel_area, int verbose);
void readVelocitySlice(FILE *fp, float *velocity, int iz, int nyv, int nxv);
-int write_FFT_DataFile(FILE *fp, complex *data, Area data_area, int fldr, int nt, int nfft, int nw, int nw_low, float dt, int out_su, int conjg, int verbose);
+int write_FFT_DataFile(FILE *fp, complex *data, Area data_area, int fldr, int nt, int nfft, int nw, int nw_low, float dt, int out_su, int conjg, int freq, int verbose);
int write_ImageFile(FILE *fp, float *data, Area data_area, int fldr, int d, int out_su, int verbose);
@@ -129,6 +129,7 @@ char *sdoc[] = {
" filter_inc=1 ............. the increment in dz for the Li filter (0=off)",
" nterms=1 ................. the number of terms in the paraxial expansion",
" OUTPUT DEFINITION ",
+" freq=0 ................... output in frequency slices",
" snap=0 ................... snapshots (not yet implemented)",
" beam=0 ................... beams ",
" verbose=0 ................ >1: shows various parameters and results",
@@ -152,7 +153,7 @@ int main(int argc, char *argv[])
{
FILE *out_file, *vel_file, *src_file, *beam_file;
size_t nread, size, size_out;
- int verbose, method, ntraces, oper_opt, MB, out_su, verb_root;
+ int verbose, method, ntraces, oper_opt, MB, out_su, out_freq, verb_root;
int nxv, nyv, nzv, dstep, fldr;
int nt, err;
int ntap, tap_opt, order, McC, oplx, oply, fine;
@@ -239,6 +240,7 @@ int main(int argc, char *argv[])
if(!getparint("fine", &fine)) fine = 2;
if(!getparint("beam", &beam)) beam = 0;
if(!getparint("add", &add)) add = 0;
+ if(!getparint("freq", &out_freq)) out_freq = 0;
if(!getparint("verbose", &verbose)) verbose = 0;
if(!ISODD(oplx)) oplx += 1;
@@ -651,7 +653,7 @@ int main(int argc, char *argv[])
fldr = is+1;
write_FFT_DataFile(out_file, rec_all, shot_area, fldr,
- nt, nfft, nw, nw_low, dt, out_su, conjg, verbose);
+ nt, nfft, nw, nw_low, dt, out_su, conjg, out_freq, verbose);
fclose(out_file);
if (verbose)
diff --git a/extrap3d/main/extrap3d.c b/extrap3d/main/extrap3d.c
index 17f6dab3b55b671f1d4021c07cd18fbef4184e94..6a4d6d3389122f55fbbc0803794c1ec6a1f15c25 100644
--- a/extrap3d/main/extrap3d.c
+++ b/extrap3d/main/extrap3d.c
@@ -22,7 +22,7 @@ int openVelocityFile(char *file_vel, FILE **fp, Area *vel_area, int verbose);
void readVelocitySlice(FILE *fp, float *velocity, int iz, int nyv, int nxv);
-int write_FFT_DataFile(FILE *fp, complex *data, Area data_area, int fldr, int nt, int nfft, int nw, int nw_low, float dt, int out_su, int conjg, int verbose);
+int write_FFT_DataFile(FILE *fp, complex *data, Area data_area, int fldr, int nt, int nfft, int nw, int nw_low, float dt, int out_su, int conjg, int freq, int verbose);
int read_FFT_DataFile(FILE *fp, complex *data, Area vel_area, int nfft, int nw, int nw_low, int *tr_read_in, int *tr_shot, int *ixmin, int *ixmax, int *iymin, int *iymax, int *sx, int *sy, int conjg, int verbose);
@@ -107,6 +107,7 @@ char *sdoc[] = {
" filter_inc=1 ............. the increment in dz for the Li filter (0=off)",
" nterms=1 ................. the number of terms in the paraxial expansion",
" OUTPUT DEFINITION ",
+" freq=0 ................... output in frequency slices",
" snap=0 ................... snapshots (not yet implemented)",
" beam=0 ................... beams ",
" verbose=0 ................ >1: shows various parameters and results",
@@ -136,7 +137,7 @@ int main(int argc, char *argv[])
size_t nread, bytes, size, trace_sz, size_out;
int verbose, method, ntraces, verb_root;
int nxv, nyv, nzv, binary_file, dstep, id, id1, id2;
- int d, nt, ndepth, i, j, conjg, conjgs, mode, out_su;
+ int d, nt, ndepth, i, j, conjg, conjgs, mode, out_su, out_freq;
int ntap, tap_opt, order, McC, oplx, oply, fine, MB;
int stackmigr, imc, area, ixmin, ixmax, iymin, iymax, ns;
int nfft, nfreq, nw_high, nw_low, nw, sx, sy, ix, iy;
@@ -206,6 +207,7 @@ int main(int argc, char *argv[])
if(!getparfloat("weights", &weights)) weights = 1e-2;
if(!getparint("fine", &fine)) fine = 2;
if(!getparint("beam", &beam)) beam = 0;
+ if(!getparint("freq", &out_freq)) out_freq = 0;
if(!getparint("verbose", &verbose)) verbose = 0;
if(!ISODD(oplx)) oplx += 1;
@@ -602,7 +604,7 @@ int main(int argc, char *argv[])
assert( out_file );
write_FFT_DataFile(out_file, rec_all, shot_area, (is+1),
- nt, nfft, nw, nw_low, dt, out_su, conjg, verb_root);
+ nt, nfft, nw, nw_low, dt, out_su, conjg, out_freq, verb_root);
fclose(out_file);
}
diff --git a/extrap3d/main/zomodel3d.c b/extrap3d/main/zomodel3d.c
index 39f66a2591032f7728c48c8d3d6795d284441fad..1bb9a53509a0b5f72777be8d027e0a30fd952dc1 100644
--- a/extrap3d/main/zomodel3d.c
+++ b/extrap3d/main/zomodel3d.c
@@ -24,7 +24,7 @@ void cr1fft(complex *cdata, float *rdata, int n, int sign);
/****** IO routines *******/
int openVelocityFile(char *file_vel, FILE **fp, Area *vel_area, int verbose);
void readVelocitySlice(FILE *fp, float *velocity, int iz, int nyv, int nxv);
-int write_FFT_DataFile(FILE *fp, complex *data, Area data_area, int fldr, int nt, int nfft, int nw, int nw_low, float dt, int out_su, int conjg, int verbose);
+int write_FFT_DataFile(FILE *fp, complex *data, Area data_area, int fldr, int nt, int nfft, int nw, int nw_low, float dt, int out_su, int conjg, int freq, int verbose);
/***** Operator tables *****/
void tablecalc_2D(int oplx, int oply, int nx, float dx, int ny, float dy,
@@ -64,7 +64,6 @@ char *sdoc[] = {
" Optional parameters:",
" ",
" conjg=0 .................. 1: take complex conjugate of input data",
-" verbose=0 ................ >1: shows various parameters and results",
" VELOCITY MODEL ",
" dxv=dxv .................. stepsize in x-direction of velocity model ",
" dyv=dyv .................. stepsize in y-direction of velocity model ",
@@ -113,6 +112,9 @@ char *sdoc[] = {
" fine=2 ................... fine sampling in operator table",
" filter_inc=1 ............. the increment in dz for the Li filter (0=off)",
" nterms=1 ................. the number of terms in the paraxial expansion",
+" OUTPUT DEFINITION ",
+" freq=0 ................... output in frequency slices",
+" verbose=0 ................ >1: shows various parameters and results",
" ",
" Options for method:",
" - 1 = Direct 2D convolution (Walter's scheme)",
@@ -148,7 +150,7 @@ int main(int argc, char *argv[])
int nfft, nfreq, nw_high, nw_low, nw, ix, iy;
int sxy, iw, sign, ixv, iyv;
int nxy, pos, id, idp, oper_opt;
- int out_su, nterms, filter_inc, interpol;
+ int out_freq, out_su, nterms, filter_inc, interpol;
Area shot_area;
float alpha, weight;
float depth1,depth2;
@@ -217,6 +219,7 @@ int main(int argc, char *argv[])
if(!getparint("fine", &fine)) fine = 2;
if(!getparfloat("dt", &dt)) dt = 0.004;
if(!getparint("nt", &nt)) nt = 512;
+ if(!getparint("freq", &out_freq)) out_freq = 0;
if(!getparint("verbose", &verbose)) verbose = 0;
if(!ISODD(oplx)) oplx += 1;
@@ -699,7 +702,7 @@ int main(int argc, char *argv[])
out_file = fopen( file_out, "w+" ); assert( out_file );
write_FFT_DataFile(out_file, rec_field, shot_area, 1,
- nt, nfft, nw, nw_low, dt, out_su, conjg, verbose);
+ nt, nfft, nw, nw_low, dt, out_su, conjg, out_freq, verbose);
fclose(out_file);
}
diff --git a/fdelmodc/Makefile b/fdelmodc/Makefile
index 543b2122ef369771ec9b5466b9bbe304587f42e6..56b3c4db58b1edee66c9341b593d29f36afe820d 100644
--- a/fdelmodc/Makefile
+++ b/fdelmodc/Makefile
@@ -7,7 +7,7 @@ include ../Make_include
ALLINC = -I.
LIBS += -L$L -lgenfft -lm $(LIBSM)
#LIBS += -L$L -lgenfft -lm -lc
-#OPTC = -g -Wall -fsignaling-nans -O0
+#OPTC = -g -Wall -fsignaling-nans -O0 -fopenmp
#OPTC += -fopenmp -Waddress
#OPTC := $(subst -O3 -ffast-math, -O1 -g ,$(OPTC))
#PGI options for compiler feedback
@@ -27,6 +27,7 @@ SRCC = $(PRG).c \
acoustic6.c \
viscoacoustic4.c \
elastic4.c \
+ elastic4dc.c \
elastic6.c \
viscoelastic4.c \
defineSource.c \
diff --git a/fdelmodc/demo/fdelmodc_doublecouple.scr b/fdelmodc/demo/fdelmodc_doublecouple.scr
index e2eebc45e8bb727b6a94306849af4c3e9e89741b..ffea7d400c08043b744be9e5f73fa048f0f9d309 100755
--- a/fdelmodc/demo/fdelmodc_doublecouple.scr
+++ b/fdelmodc/demo/fdelmodc_doublecouple.scr
@@ -20,12 +20,12 @@ export filero=model_ro.su
export OMP_NUM_THREADS=1
../fdelmodc \
- file_cp=$filecp file_cs=$filecs file_den=$filero \
- ischeme=3 \
+ file_cp=$filecp file_den=$filero \
+ ischeme=5 \
src_type=9 dip=0 tmod=2.0 \
file_src=wavelet.su verbose=2 \
- file_rcv=recS.su \
- file_snap=snapS0.su \
+ file_rcv=rec4S.su \
+ file_snap=snap4S0.su \
xrcv1=0 xrcv2=2100 dxrcv=15 \
zrcv1=400 zrcv2=400 \
rec_type_vx=1 rec_int_vx=1 \
diff --git a/fdelmodc/demo/modelOilGas.scr b/fdelmodc/demo/modelOilGas.scr
index 50055435d7d4cb61e46944c9811a0a899c98f299..5fa26904388a3d67e07cc72a6bc7861d1b5bec16 100755
--- a/fdelmodc/demo/modelOilGas.scr
+++ b/fdelmodc/demo/modelOilGas.scr
@@ -1,6 +1,6 @@
#!/bin/bash
-#export PATH=$HOME/bin:$HOME/src/OpenSource/utils:$PATH:
+export PATH=$HOME/bin:$HOME/src/OpenSource/utils:$PATH:
cp=2000
rho=2500
@@ -16,12 +16,14 @@ makemod sizex=5000 sizez=2500 dx=$dx dz=$dx cp0=$cp ro0=$rho \
makewave w=g1 fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
-export OMP_NUM_THREADS=8
+export OMP_NUM_THREADS=1
zsrc=1100
zsrc=0
-~/bin/fdelmodc \
+which fdelmodc
+
+fdelmodc \
file_cp=syncl_cp.su ischeme=1 iorder=4 \
file_den=syncl_ro.su \
file_src=wave.su \
diff --git a/fdelmodc/demo/staal.scr b/fdelmodc/demo/staal.scr
new file mode 100755
index 0000000000000000000000000000000000000000..a99d0a59fa1444de173a8c3dca12d149049d9712
--- /dev/null
+++ b/fdelmodc/demo/staal.scr
@@ -0,0 +1,20 @@
+#!/bin/bash
+
+makemod file_base=betonplaat3.su cp0=3700 ro0=2400 cs0=3400 verbose=1 sizex=5000 sizez=300 dx=5.0 dz=5.0
+
+makewave file_out=betonplaatwave.su dt=1.0e-4 nt=10000 w=g2 shift=1 fp=6 verbose=1
+
+../fdelmodc file_cp=betonplaat3_cp.su file_cs=betonplaat3_cs.su \
+file_den=betonplaat3_ro.su file_src=betonplaatwave.su \
+file_rcv=betonplaat2recv.su file_snap=betonplaat2snap.su \
+ischeme=3 tmod=5 \
+top=1 left=4 right=1 bottom=1 \
+src_type=3 xsrc=2500 zsrc=0.0 \
+tsnap1=1.0e-2 tsnap2=5 dtsnap=1.0e-2 dxsnap=5.0 dzsnap=5.0 \
+dtrcv=1.0e-2 verbose=1
+
+!sufrac power=1 < betonplaat3recv_rvz.su > sufracbetonplaat3.su
+!sustrip < sufracbetonplaat3.su | b2a n1=2051 > Sufrac_betonplaat3_Tzz.txt
+!suximage<sufracbetonplaat3.su title="FDM shot record Concrete Slab Sourcetype=Tzz [Left=4 Right=1]" label2="Lateral Position [mm]"
+label1="Time [ms]" wbox=2500 hbox=500 cmap=hsv2 legend=1
+
diff --git a/fdelmodc/elastic4dc.c b/fdelmodc/elastic4dc.c
new file mode 100644
index 0000000000000000000000000000000000000000..d119cfb3db178a7f3d5ade881188634191b1c686
--- /dev/null
+++ b/fdelmodc/elastic4dc.c
@@ -0,0 +1,160 @@
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include"fdelmodc.h"
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+
+int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, int izsrc, float *vx, float *vz, float *tzz,
+float *txx, float *txz, float *rox, float *roz, float *l2m, float **src_nwav, int verbose);
+
+int storeSourceOnSurface(modPar mod, srcPar src, bndPar bnd, int ixsrc, int izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, int verbose);
+
+int reStoreSourceOnSurface(modPar mod, srcPar src, bndPar bnd, int ixsrc, int izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, int verbose);
+
+int boundariesP(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, int itime, int verbose);
+
+int boundariesV(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, int itime, int verbose);
+
+int elastic4dc(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, int izsrc, float **src_nwav, float *vx,
+float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float
+*l2m, float *lam, float *mul, int verbose)
+{
+/*********************************************************************
+ COMPUTATIONAL OVERVIEW OF THE 4th ORDER STAGGERED GRID:
+
+ The captial symbols T (=Txx,Tzz) Txz,Vx,Vz represent the actual grid
+ The indices ix,iz are related to the T grid, so the capital T
+ symbols represent the actual modelled grid.
+
+ one cel (iz,ix)
+ |
+ V extra column of vx,txz
+ |
+ ------- V
+ | txz vz| txz vz txz vz txz vz txz vz txz vz txz
+ | |
+ | vx t | vx t vx t vx t vx t vx t vx
+ -------
+ txz vz txz vz txz vz txz vz txz vz txz vz txz
+
+ vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
+ | | | | | | |
+ txz vz txz Vz--Txz-Vz--Txz-Vz Txz-Vz txz vz txz
+ | | | | | | |
+ vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
+ | | | | | | |
+ txz vz txz Vz Txz-Vz Txz-Vz Txz-Vz txz vz txz
+ | | | | | | |
+ vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
+ | | | | | | |
+ txz vz txz Vz Txz-Vz Txz-Vz Txz-Vz txz vz txz
+ | | | | | | |
+ vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
+
+ txz vz txz vz txz vz txz vz txz vz txz vz txz
+
+ vx t vx t vx t vx t vx t vx t vx
+
+ txz vz txz vz txz vz txz vz txz vz txz vz txz <--|
+ |
+ extra row of txz/vz |
+
+ AUTHOR:
+ Jan Thorbecke (janth@xs4all.nl)
+ The Netherlands
+
+***********************************************************************/
+
+ float c1, c2;
+ float dvx, dvz;
+ int ix, iz;
+ int n1;
+
+ c1 = 9.0/8.0;
+ c2 = -1.0/24.0;
+ n1 = mod.naz;
+
+ /* calculate vx for all grid points except on the virtual boundary*/
+#pragma omp for private (ix, iz) nowait schedule(guided,1)
+ for (ix=mod.ioXx; ix<mod.ieXx; ix++) {
+#pragma ivdep
+ for (iz=mod.ioXz; iz<mod.ieXz; iz++) {
+ vx[ix*n1+iz] -= rox[ix*n1+iz]*(
+ c1*(txx[ix*n1+iz] - txx[(ix-1)*n1+iz] +
+ txz[ix*n1+iz+1] - txz[ix*n1+iz]) +
+ c2*(txx[(ix+1)*n1+iz] - txx[(ix-2)*n1+iz] +
+ txz[ix*n1+iz+2] - txz[ix*n1+iz-1]) );
+ }
+ }
+
+ /* calculate vz for all grid points except on the virtual boundary */
+#pragma omp for private (ix, iz) schedule(guided,1)
+ for (ix=mod.ioZx; ix<mod.ieZx; ix++) {
+#pragma ivdep
+ for (iz=mod.ioZz; iz<mod.ieZz; iz++) {
+ vz[ix*n1+iz] -= roz[ix*n1+iz]*(
+ c1*(tzz[ix*n1+iz] - tzz[ix*n1+iz-1] +
+ txz[(ix+1)*n1+iz] - txz[ix*n1+iz]) +
+ c2*(tzz[ix*n1+iz+1] - tzz[ix*n1+iz-2] +
+ txz[(ix+2)*n1+iz] - txz[(ix-1)*n1+iz]) );
+ }
+ }
+
+ /* Add force source */
+ if (src.type > 5) {
+ applySource(mod, src, wav, bnd, itime, ixsrc, izsrc, vx, vz, tzz, txx, txz, rox, roz, l2m, src_nwav, verbose);
+ }
+
+
+ /* boundary condition clears velocities on boundaries */
+ boundariesP(mod, bnd, vx, vz, tzz, txx, txz, rox, roz, l2m, lam, mul, itime, verbose);
+
+ /* calculate Txx/tzz for all grid points except on the virtual boundary */
+#pragma omp for private (ix, iz, dvx, dvz) nowait schedule(guided,1)
+ for (ix=mod.ioPx; ix<mod.iePx; ix++) {
+#pragma ivdep
+ for (iz=mod.ioPz; iz<mod.iePz; iz++) {
+ dvx = c1*(vx[(ix+1)*n1+iz] - vx[ix*n1+iz]) +
+ c2*(vx[(ix+2)*n1+iz] - vx[(ix-1)*n1+iz]);
+ dvz = c1*(vz[ix*n1+iz+1] - vz[ix*n1+iz]) +
+ c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]);
+ txx[ix*n1+iz] -= l2m[ix*n1+iz]*dvx + l2m[ix*n1+iz]*dvz;
+ tzz[ix*n1+iz] -= l2m[ix*n1+iz]*dvz + l2m[ix*n1+iz]*dvx;
+ }
+ }
+
+
+
+ /* calculate Txz for all grid points except on the virtual boundary */
+/*
+#pragma omp for private (ix, iz) schedule(guided,1)
+ for (ix=mod.ioTx; ix<mod.ieTx; ix++) {
+#pragma ivdep
+ for (iz=mod.ioTz; iz<mod.ieTz; iz++) {
+ txz[ix*n1+iz] -= mul[ix*n1+iz]*(
+ c1*(vx[ix*n1+iz] - vx[ix*n1+iz-1] +
+ vz[ix*n1+iz] - vz[(ix-1)*n1+iz]) +
+ c2*(vx[ix*n1+iz+1] - vx[ix*n1+iz-2] +
+ vz[(ix+1)*n1+iz] - vz[(ix-2)*n1+iz]) );
+ }
+ }
+*/
+
+ /* Add stress source */
+ if (src.type < 6) {
+ applySource(mod, src, wav, bnd, itime, ixsrc, izsrc, vx, vz, tzz, txx, txz, rox, roz, l2m, src_nwav, verbose);
+ }
+
+ /* check if there are sources placed on the boundaries */
+ storeSourceOnSurface(mod, src, bnd, ixsrc, izsrc, vx, vz, tzz, txx, txz, verbose);
+
+ /* Free surface: calculate free surface conditions for stresses */
+ boundariesV(mod, bnd, vx, vz, tzz, txx, txz, rox, roz, l2m, lam, mul, itime, verbose);
+
+ /* restore source positions on the edge */
+ reStoreSourceOnSurface(mod, src, bnd, ixsrc, izsrc, vx, vz, tzz, txx, txz, verbose);
+
+ return 0;
+}
diff --git a/fdelmodc/fdelmodc.c b/fdelmodc/fdelmodc.c
index d9a70e2b032510457c9c2fbdd15e39b2350515e1..251937797e7baa88cee21af505fd2f286b31b0d0 100644
--- a/fdelmodc/fdelmodc.c
+++ b/fdelmodc/fdelmodc.c
@@ -46,6 +46,8 @@ int viscoacoustic4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, in
int elastic4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, int izsrc, float **src_nwav, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, int verbose);
+int elastic4dc(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, int izsrc, float **src_nwav, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, int verbose);
+
int viscoelastic4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, int izsrc, float **src_nwav, float *vx, float
*vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, float *ts, float *tep, float
*tes, float *r, float *q, float *p, int verbose);
@@ -96,7 +98,7 @@ char *sdoc[] = {
" dt= ............... read from file_src: if dt is set it will interpolate file_src to dt sampling",
"" ,
" OPTIONAL PARAMETERS:",
-" ischeme=3 ......... 1=acoustic, 2=visco-acoustic 3=elastic, 4=visco-elastic",
+" ischeme=3 ......... 1=acoustic, 2=visco-acoustic 3=elastic, 4=visco-elastic, 5=double-couple",
" tmod=(nt-1)*dt .... total modeling time (nt from file_src)",
" ntaper=0 .......... length of taper in points at edges of model",
" npml=35 ........... length of PML layer in points at edges of model",
@@ -584,6 +586,10 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, izsrc, it, src_nwav, verbose)
vx, vz, tzz, txx, txz, rox, roz, l2m, lam, mul,
tss, tep, tes, r, q, p, verbose);
break;
+ case 5 : /* Elastic FD kernel with S-velocity set to zero*/
+ elastic4dc(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
+ vx, vz, tzz, txx, txz, rox, roz, l2m, lam, mul, verbose);
+ break;
}
/* write samples to file if rec.nt samples are calculated */
diff --git a/fdelmodc/getParameters.c b/fdelmodc/getParameters.c
index 56fb84d5b177a53d60ec8ec468a13ec3de7c1d73..8200cd644d2af60a57415c7f092050e83cf1ebf3 100644
--- a/fdelmodc/getParameters.c
+++ b/fdelmodc/getParameters.c
@@ -78,7 +78,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
if (!getparstring("file_den",&mod->file_ro)) {
verr("parameter file_den required!");
}
- if (mod->ischeme>2) {
+ if (mod->ischeme>2 && mod->ischeme!=5) {
if (!getparstring("file_cs",&mod->file_cs)) {
verr("parameter file_cs required!");
}
@@ -110,7 +110,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
if (nz != n1)
vwarn("nz differs for file_cp and file_den!");
- if (mod->ischeme>2) {
+ if (mod->ischeme>2 && mod->ischeme!=5) {
getModelInfo(mod->file_cs, &n1, &n2, &d1, &d2, &zstart, &xstart, &cs_min, &cs_max, &axis, 1, verbose);
mod->cs_max = cs_max;
mod->cs_min = cs_min;
@@ -123,6 +123,12 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
if (nz != n1)
vwarn("nz differs for file_cp and file_cs!");
}
+ if (mod->ischeme==5) {
+ cs_max=0.0; cs_min=0.0;
+ mod->cs_max = cs_max;
+ mod->cs_min = cs_min;
+ }
+
mod->dz = dz;
mod->dx = dx;
mod->nz = nz;
@@ -223,6 +229,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
if (mod->ischeme == 2) vmess("Visco-Acoustic staggered grid, pressure/velocity");
if (mod->ischeme == 3) vmess("Elastic staggered grid, stress/velocity");
if (mod->ischeme == 4) vmess("Visco-Elastic staggered grid, stress/velocity");
+ if (mod->ischeme == 5) vmess("Acoustic staggered grid, Txx/Tzz/velocity");
if (mod->grid_dir) vmess("Time reversed modelling");
else vmess("Forward modelling");
vmess("*******************************************");
@@ -233,7 +240,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
vmess("zmin = %8.4f zmax = %8.4f", sub_z0, zmax);
vmess("xmin = %8.4f xmax = %8.4f", sub_x0, xmax);
vmess("min(cp) = %9.3f max(cp) = %9.3f", cp_min, cp_max);
- if (mod->ischeme>2) vmess("min(cs) = %9.3f max(cs) = %9.3f", cs_min, cs_max);
+ if (mod->ischeme>2 && mod->ischeme!=5) vmess("min(cs) = %9.3f max(cs) = %9.3f", cs_min, cs_max);
vmess("min(ro) = %9.3f max(ro) = %9.3f", ro_min, ro_max);
if (mod->ischeme==2 || mod->ischeme==4) {
if (mod->file_qp!=NULL) vmess("Qp from file %s ", mod->file_qp);
diff --git a/fdelmodc/readModel.c b/fdelmodc/readModel.c
index 4fc6d24f2f295de162c793930620505dc198dcf9..a3975ef3cb657a813aeaf0b6b91054b321c94c38 100644
--- a/fdelmodc/readModel.c
+++ b/fdelmodc/readModel.c
@@ -85,8 +85,8 @@ int readModel(modPar mod, bndPar bnd, float *rox, float *roz, float *l2m, float
nread = fread(&hdr, 1, TRCBYTES, fpro);
assert(nread == TRCBYTES);
- if (mod.ischeme>2) {
- cs = (float *)malloc(nz*nx*sizeof(float));
+ cs = (float *)calloc(nz*nx,sizeof(float));
+ if (mod.ischeme>2 && mod.ischeme!=5) {
fpcs = fopen( mod.file_cs, "r" );
assert( fpcs != NULL);
nread = fread(&hdr, 1, TRCBYTES, fpcs);
@@ -120,7 +120,7 @@ int readModel(modPar mod, bndPar bnd, float *rox, float *roz, float *l2m, float
assert (nread == hdr.ns);
nread = fread(&ro[i*nz], sizeof(float), hdr.ns, fpro);
assert (nread == hdr.ns);
- if (mod.ischeme>2) {
+ if (mod.ischeme>2 && mod.ischeme!=5) {
nread = fread(&cs[i*nz], sizeof(float), hdr.ns, fpcs);
assert (nread == hdr.ns);
}
@@ -181,7 +181,7 @@ int readModel(modPar mod, bndPar bnd, float *rox, float *roz, float *l2m, float
if (nread==0) break;
nread = fread(&hdr, 1, TRCBYTES, fpro);
if (nread==0) break;
- if (mod.ischeme>2) {
+ if (mod.ischeme>2 && mod.ischeme!=5) {
nread = fread(&hdr, 1, TRCBYTES, fpcs);
if (nread==0) break;
}
@@ -197,7 +197,7 @@ int readModel(modPar mod, bndPar bnd, float *rox, float *roz, float *l2m, float
}
fclose(fpcp);
fclose(fpro);
- if (mod.ischeme>2) fclose(fpcs);
+ if (mod.ischeme>2 && mod.ischeme!=5) fclose(fpcs);
if (fpqp != NULL) fclose(fpqp);
if (fpqs != NULL) fclose(fpqs);
@@ -784,7 +784,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
free(cp);
free(ro);
- if (mod.ischeme>2) free(cs);
+ free(cs);
return 0;
}
diff --git a/marchenko/Makefile b/marchenko/Makefile
index a0017f27943755d76ef6c70647c2e3d282973437..237c134e904d7c4b7311e63065c4c09282ace67f 100644
--- a/marchenko/Makefile
+++ b/marchenko/Makefile
@@ -7,7 +7,7 @@ LIBS += -L$L -lgenfft -lm $(LIBSM)
#ALL: fmute marchenko marchenko2
-ALL: fmute marchenko marchenko_tshift
+ALL: fmute marchenko
SRCJ = fmute.c \
getFileInfo.c \
@@ -37,23 +37,6 @@ SRCH = marchenko.c \
OBJJ = $(SRCJ:%.c=%.o)
-SRCT = marchenko_tshift.c \
- getFileInfo.c \
- readData.c \
- readShotData.c \
- readTinvData.c \
- applyMute_tshift.c \
- writeData.c \
- writeDataIter.c \
- wallclock_time.c \
- name_ext.c \
- verbosepkg.c \
- atopkge.c \
- docpkge.c \
- getpars.c
-
-OBJT = $(SRCT:%.c=%.o)
-
fmute: $(OBJJ)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o fmute $(OBJJ) $(LIBS)
@@ -62,15 +45,12 @@ OBJH = $(SRCH:%.c=%.o)
marchenko: $(OBJH)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko $(OBJH) $(LIBS)
-marchenko_tshift: $(OBJT)
- $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko_tshift $(OBJT) $(LIBS)
-
OBJH2 = $(SRCH2:%.c=%.o)
marchenko2: $(OBJH2)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko2 $(OBJH2) $(LIBS)
-install: fmute marchenko marchenko_tshift
+install: fmute marchenko
cp fmute $B
cp marchenko $B
diff --git a/utils/demo/makemodsmooth.scr b/utils/demo/makemodsmooth.scr
new file mode 100755
index 0000000000000000000000000000000000000000..587e7aa8200cae47b4b735e28b0825fafcc2715d
--- /dev/null
+++ b/utils/demo/makemodsmooth.scr
@@ -0,0 +1,56 @@
+#!/bin/bash
+
+dx=1.0
+dt=0.0001
+dx=2.5
+dt=0.0004
+
+makewave w=fw fmin=0 flef=5 frig=70 fmax=90 dt=$dt file_out=wavefw4.su nt=16384 t0=0.4 scale=0 scfft=1
+makewave w=fw fmin=1.5 flef=4 frig=45 fmax=60 dt=$dt file_out=wavefw4.su nt=16384 t0=0.4 scale=0 scfft=1
+
+../makemod sizex=12000 sizez=4500 dx=$dx dz=$dx cp0=1654 ro0=1000 \
+ orig=-6000,0 file_base=blind.su reflectivity=1 supersmooth=1 sigma=0.8 \
+ intt=def x=-6000,0,6000 z=450,450,450 poly=0 cp=2316 ro=2500 grad=0 gradcp=0 \
+ intt=def x=-6000,0,6000 z=590,860,980 poly=2 cp=1734 ro=1500 grad=0 gradcp=0 gradt=2 \
+ intt=elipse x=4100 z=1000 poly=0 cp=1022 ro=1911 var=300,900 grad=0 gradcp=500 \
+ intt=rough x=-6000,6000 z=1000,1000 poly=0 cp=2222 ro=1111 var=120,4,13 grad=0 gradcp=0 \
+ intt=def x=-6000,2000 z=490,3000 poly=0 cp=2110 ro=2300 grad=0 gradcp=0 \
+ intt=def x=-6000,2000 z=550,3060 poly=0 cp=2438 ro=1300 grad=0 gradcp=0 \
+ intt=def x=-6000,2000 z=650,3160 poly=0 cp=2106 ro=2300 grad=0 gradcp=0 \
+ intt=def x=-6000,2000 z=880,3390 poly=0 cp=2512 ro=1700 grad=0 gradcp=0 \
+ intt=def x=-6000,2000 z=1200,3710 poly=0 cp=2912 ro=1000 grad=0 gradcp=0 \
+ intt=def x=-1000,6000 z=2700,1710 poly=0 cp=2442 ro=2700 grad=0 gradcp=0 \
+ intt=def x=-1000,6000 z=2900,1910 poly=0 cp=2642 ro=2900 grad=0 gradcp=0 \
+ intt=def x=-1000,6000 z=2960,1970 poly=0 cp=2242 ro=2400 grad=0 gradcp=0 \
+ intt=def x=-1000,6000 z=3400,2410 poly=0 cp=3012 ro=2500 grad=0 gradcp=0 \
+ intt=def x=-1145,2860 z=2015,2150 poly=0 cp=2012 ro=1300 grad=0 gradcp=0 \
+ intt=def x=-1145,2860 z=2215,2350 poly=0 cp=2512 ro=1300 grad=0 gradcp=0 \
+ intt=def x=-1145,2860 z=2315,2450 poly=0 cp=2712 ro=2500 grad=0 gradcp=0 \
+ intt=def x=-6000,-3000,0,2000,6000 z=2000,2500,2600,3000,3100 poly=2 cp=2822 ro=1811 grad=0 gradcp=0 \
+ intt=def x=-6000,6000 z=3500,3500 poly=0 cp=3104 ro=1500 grad=0 gradcp=0 \
+ intt=def x=-6000,2000 z=4000,3500 poly=0 cp=3504 ro=2500 grad=0 gradcp=0 \
+ intt=def x=2000,6000 z=3500,3500 poly=0 cp=3504 ro=2500 grad=0 gradcp=0
+
+fdelmodc \
+ file_cp=blind_cp.su ischeme=1 iorder=4 \
+ file_den=blind_ro.su \
+ file_src=wavefw4.su \
+ file_rcv=shot_dx${dx}.su \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=1 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.4 fmax=40 \
+ verbose=2 \
+ tmod=8.588 \
+ dxrcv=12.5 \
+ xrcv1=-6000 xrcv2=6000 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=0 \
+ ntaper=400 \
+ dtsnap=0.1 dxsnap=10 dzsnap=10 tsnap2=8.5 snapwithbnd=1 \
+ tapleft=2 tapright=2 taptop=2 tapbottom=2
+
diff --git a/utils/makemod.c b/utils/makemod.c
index 909a060d4c8e8d47d8ee483f905d996703f35798..d6bfb9782ed83dfbe77ab2b04740e5c539907e85 100644
--- a/utils/makemod.c
+++ b/utils/makemod.c
@@ -966,3 +966,4 @@ int main(int argc, char **argv)
return 0;
}
+