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/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;
}