diff --git a/fdelmodc/Makefile b/fdelmodc/Makefile
index ed9746c4d455bdabf0df42ac83adbe7a9b9361f0..4153b8cfc59c08f824661f5b94cb461f77118712 100644
--- a/fdelmodc/Makefile
+++ b/fdelmodc/Makefile
@@ -9,11 +9,13 @@ ALLINC = -I.
#OPTC = -g -Wall -fsignaling-nans -O0 -fopenmp
#OPTC += -fopenmp -Waddress
#OPTC += -qopt-report=5 -g
+#OPTC = -fp-trap=invalid,overflow -traceback -g -O1 -qopenmp
#OPTC := $(subst -O3 -ffast-math, -O1 -g ,$(OPTC))
#PGI options for compiler feedback
#OPTC += -Mprof=lines
#OPTC += -qopt-report
#LDFLAGS += -Mprof=lines
+#LDFLAGS += -Wl,-z,muldefs
all: fdelmodc
diff --git a/fdelmodc/applySource.c b/fdelmodc/applySource.c
index ea1a6135a2799208f2dac312f8efe8a960de8b7b..77dd4529b83010ce4adf257ebcd8cc3e9b28fe99 100644
--- a/fdelmodc/applySource.c
+++ b/fdelmodc/applySource.c
@@ -103,7 +103,8 @@ int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int i
id2 = id1+1;
/* delay not reached or no samples left in source wavelet? */
- if ( (time < 0.0) || ( (itime*dt) >= src.tend[isrc]) ) continue;
+ if ( (time < 0.0) || ((itime*dt) >= src.tend[isrc]) || id2 > wav.nt ) continue;
+ //if (isrc==0) fprintf(stderr,"isrc=%d id1=%d time=%f wav.nt=%d tend=%f\n", isrc, id1, time, wav.nt, src.tend[isrc]);
// fprintf(stderr,"isrc=%d ix=%d iz=%d src.x=%d src.z=%d\n", isrc, ix, iz, src.x[isrc], src.z[isrc]);
diff --git a/fdelmodc/boundaries.c b/fdelmodc/boundaries.c
index ffda5cfa0eacda0d59f35c396390ce7de54d5286..55edb021074d95afc6247c2ae7a5492d81acc206 100644
--- a/fdelmodc/boundaries.c
+++ b/fdelmodc/boundaries.c
@@ -134,7 +134,7 @@ int boundariesP(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float
/************************************************************/
npml=bnd.npml; /* lenght of pml in grid-points */
- if ( (npml != 0) && (itime==0) && pml) {
+ if ( (npml != 0) && (allocated==0) && pml) {
#pragma omp master
{
if (allocated) {
@@ -1170,7 +1170,7 @@ int boundariesV(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float
/************************************************************/
npml=bnd.npml; /* lenght of pml in grid-points */
- if ( (npml != 0) && (itime==0) && pml) {
+ if ( (npml != 0) && (allocated==0) && pml) {
#pragma omp master
{
if (allocated) {
diff --git a/fdelmodc/defineSource.c b/fdelmodc/defineSource.c
index 5eaabdad90396856780778981ec0b33b6e2ea783..5d481b0a4de00aaff046702bef8abb0970b797b1 100644
--- a/fdelmodc/defineSource.c
+++ b/fdelmodc/defineSource.c
@@ -128,7 +128,7 @@ int defineSource(wavPar wav, srcPar src, modPar mod, recPar rec, float **src_nwa
randomWavelet(wav, src, &src_nwav[i][0], src.tbeg[i], src.tend[i], verbose);
}
else {
- memset(&ctrace[0].r,0,nfreqscale*sizeof(complex));
+ memset(&ctrace[0],0,nfreqscale*sizeof(complex));
memset(&trace[0],0,optnscale*sizeof(float));
memcpy(&trace[0],&src_nwav[i][0],n1*sizeof(float));
rc1fft(trace,ctrace,optn,-1);
diff --git a/fdelmodc/fdelmodc.c b/fdelmodc/fdelmodc.c
index 4bfb972c3abb84b4392a0f2654e4ce19131e254f..524a21609384d5d1b8025a2c16028c0981c08971 100644
--- a/fdelmodc/fdelmodc.c
+++ b/fdelmodc/fdelmodc.c
@@ -55,7 +55,7 @@ int elastic6(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsr
float *roz, float *l2m, float *lam, float *mul, int verbose);
int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *vx, float *vz, float *tzz, float *txx,
- float *txz, float *l2m, float *rox, float *roz,
+ float *txz, float *l2m, float *lam, float *rox, float *roz,
float *rec_vx, float *rec_vz, float *rec_txx, float *rec_tzz, float *rec_txz,
float *rec_p, float *rec_pp, float *rec_ss, float *rec_udp, float *rec_udvz, int verbose);
@@ -364,7 +364,7 @@ int main(int argc, char **argv)
src_nwav[0] = (float *)calloc(wav.nt*wav.nx,sizeof(float));
assert(src_nwav[0] != NULL);
for (i=0; i<wav.nx; i++) {
- src_nwav[i] = (float *)(src_nwav[0] + wav.nt*i);
+ src_nwav[i] = (float *)(src_nwav[0] + (long)(wav.nt*i));
}
}
@@ -540,7 +540,7 @@ shared (tss, tep, tes, r, q, p) \
shared (tinit, it0, it1, its) \
shared(beam_vx, beam_vz, beam_txx, beam_tzz, beam_txz, beam_p, beam_pp, beam_ss) \
shared(rec_vx, rec_vz, rec_txx, rec_tzz, rec_txz, rec_p, rec_pp, rec_ss) \
-shared (tt, t2, t3) \
+shared (tt, t2, t3, isam) \
shared (shot, bnd, mod, src, wav, rec, ixsrc, izsrc, it, src_nwav, verbose)
{
if (it==it0 && verbose>2) {
@@ -607,7 +607,12 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, izsrc, it, src_nwav, verbose)
if ( (((it-rec.delay) % rec.skipdt)==0) && (it >= rec.delay) ) {
int writeToFile, itwritten;
- writeToFile = ! ( (((it-rec.delay+NINT(mod.t0/mod.dt))/rec.skipdt)+1)%rec.nt );
+ if ((((it-rec.delay+NINT(mod.t0/mod.dt))/rec.skipdt)+1)!=0) {
+ writeToFile = ! ( (((it-rec.delay+NINT(mod.t0/mod.dt))/rec.skipdt)+1)%rec.nt );
+ }
+ else { /* when negative times passes zero-time */
+ writeToFile = 0;
+ }
itwritten = fileno*(rec.nt)*rec.skipdt;
/* Note that time step it=0 (t=0 for t**-fields t=-1/2 dt for v*-field) is not recorded */
/* negative time correction with mod.t0 for dipping plane waves modeling */
@@ -615,7 +620,7 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, izsrc, it, src_nwav, verbose)
if (isam < 0) isam = rec.nt+isam;
/* store time at receiver positions */
getRecTimes(mod, rec, bnd, it, isam, vx, vz, tzz, txx, txz,
- l2m, rox, roz,
+ l2m, lam, rox, roz,
rec_vx, rec_vz, rec_txx, rec_tzz, rec_txz,
rec_p, rec_pp, rec_ss, rec_udp, rec_udvz, verbose);
@@ -640,15 +645,16 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, izsrc, it, src_nwav, verbose)
beam_p, beam_pp, beam_ss, verbose);
}
}
+#pragma omp barrier
#pragma omp master
{
if (verbose) {
- if(!((it1-it)%perc)) fprintf(stderr,"\b\b\b\b%3d%%",it*100/mod.nt);
+ if(!((it1-it)%perc)) fprintf(stderr,"\b\b\b\b%3d%%",it*100/(it1-it0));
if(it==100)t3=wallclock_time();
if(it==500){
- t3=(wallclock_time()-t3)*(mod.nt/400.0);
- fprintf(stderr,"\r %s: Estimated total compute time for this shot = %.2fs.\n %s: Progress: %3d%%",xargv[0],t3,xargv[0],it/(mod.nt/100));
+ t3=(wallclock_time()-t3)*((it1-it0)/400.0);
+ fprintf(stderr,"\r %s: Estimated total compute time for this shot = %.2fs.\n %s: Progress: %3d%%",xargv[0],t3,xargv[0],it/((it1-it0)/100));
}
/*
diff --git a/fdelmodc/getParameters.c b/fdelmodc/getParameters.c
index c779f9942c3257d220ca760e8e1195fad83bd7e5..2297ea7c6f475cde147726e92b512eb24f8439e4 100644
--- a/fdelmodc/getParameters.c
+++ b/fdelmodc/getParameters.c
@@ -58,7 +58,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
int i, j;
int cfree;
int tapleft,tapright,taptop,tapbottom;
- int nxsrc, nzsrc;
+ int nxsrc, nzsrc, nhx;
int largeSUfile;
int is,ntraces,length_random;
float rand;
@@ -881,6 +881,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
}
src_ix0=MAX(0,NINT((xsrc-sub_x0)/dx));
src_ix0=MIN(src_ix0,nx);
+ is0 = -1*floor((nsrc-1)/2);
src_iz0=MAX(0,NINT((zsrc-sub_z0)/dz));
src_iz0=MIN(src_iz0,nz);
@@ -894,9 +895,9 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
p = sin(src_angle*grad2rad)/src_velo;
/* t=0 is defined at position src_ix0 */
for (is=0; is<nsrc; is++) {
- src->tbeg[is] = (is-src_ix0)*dx*p;
+ src->tbeg[is] = (is+is0)*dx*p;
}
- mod->t0 = MIN(-src_ix0*dx*p,(nsrc-src_ix0)*dx*p);
+ mod->t0 = MIN(is0*dx*p,(nsrc+is0)*dx*p);
/*
if (p < 0.0) {
for (is=0; is<nsrc; is++) {
@@ -916,7 +917,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
}
is0 = -1*floor((nsrc-1)/2);
- is0 = -src_ix0;
+ //is0 = -src_ix0;
for (is=0; is<nsrc; is++) {
src->x[is] = is0 + is;
src->z[is] = 0;
@@ -1000,8 +1001,8 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
/* vmess("Memory requirement for sources = %.2f MB.",sizeof(float)*(wav->nx*(wav->nt/(1024.0*1024.0))));*/
vmess("Memory requirement for sources = %.2f MB.",sizeof(float)*(nsamp/(1024.0*1024.0)));
if (src->plane) {
- vmess("Computed p-value = %f.",p);
- vmess("Maximum negative time delay is %f\n", mod->t0);
+ vmess("plane-wave: Computed p-value = %f.",p);
+ vmess("plane-wave: Maximum negative time delay is %f", mod->t0);
}
}
if (src->random) {
diff --git a/fdelmodc/getRecTimes.c b/fdelmodc/getRecTimes.c
index 545e59bfb45bf5635dee4abf37b1c931d3dffd6b..a5d3bf5faaade14f886653ce5ef363fda6f83627 100644
--- a/fdelmodc/getRecTimes.c
+++ b/fdelmodc/getRecTimes.c
@@ -17,7 +17,7 @@
* The Netherlands
**/
-int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *vx, float *vz, float *tzz, float *txx, float *txz, float *l2m, float *rox, float *roz, float *rec_vx, float *rec_vz, float *rec_txx, float *rec_tzz, float *rec_txz, float *rec_p, float *rec_pp, float *rec_ss, float *rec_udp, float *rec_udvz, int verbose)
+int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *vx, float *vz, float *tzz, float *txx, float *txz, float *l2m, float *lam, float *rox, float *roz, float *rec_vx, float *rec_vz, float *rec_txx, float *rec_tzz, float *rec_txz, float *rec_p, float *rec_pp, float *rec_ss, float *rec_udp, float *rec_udvz, int verbose)
{
int n1, ibndx, ibndz;
int irec, ix, iz, ix2, iz2, ix1, iz1;
@@ -176,7 +176,17 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
rec_p[irec*rec.nt+isam] += 0.5*field;
}
else {
- rec_p[irec*rec.nt+isam] += 0.5*(tzz[ix*n1+iz1]+tzz[ix*n1+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]);
+ field = tzz[ix*n1+iz] + 0.5*(l2m[ix*n1+iz]*dvz + lam[ix*n1+iz]*dvx);
+ dvx = c1*(vx[(ix+1)*n1+iz1] - vx[ix*n1+iz1]) +
+ c2*(vx[(ix+2)*n1+iz1] - vx[(ix-1)*n1+iz1]);
+ dvz = c1*(vz[ix*n1+iz1+1] - vz[ix*n1+iz1]) +
+ c2*(vz[ix*n1+iz1+2] - vz[ix*n1+iz1-1]);
+ field += tzz[ix*n1+iz1] + 0.5*(l2m[ix*n1+iz1]*dvz + lam[ix*n1+iz1]*dvx);
+ rec_p[irec*rec.nt+isam] += 0.5*field;
}
}
else if (rec.int_p == 2) {
@@ -194,7 +204,17 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
rec_p[irec*rec.nt+isam] += 0.5*field;
}
else {
- rec_p[irec*rec.nt+isam] += 0.5*(tzz[ix1*n1+iz]+tzz[ix*n1+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]);
+ field = tzz[ix*n1+iz] + 0.5*(l2m[ix*n1+iz]*dvz + lam[ix*n1+iz]*dvx);
+ dvx = c1*(vx[(ix1+1)*n1+iz] - vx[ix1*n1+iz]) +
+ c2*(vx[(ix+2)*n1+iz] - vx[(ix1-1)*n1+iz]);
+ dvz = c1*(vz[ix1*n1+iz+1] - vz[ix1*n1+iz]) +
+ c2*(vz[ix1*n1+iz+2] - vz[ix1*n1+iz-1]);
+ field += tzz[ix1*n1+iz] + 0.5*(l2m[ix1*n1+iz]*dvz + lam[ix1*n1+iz]*dvx);
+ rec_p[irec*rec.nt+isam] += 0.5*field;
}
}
else {
@@ -240,13 +260,21 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
rec_vz[irec*rec.nt+isam] += 0.5*field;
}
else {
- rec_vz[irec*rec.nt+isam] += 0.5*(vz[ix*n1+iz2]+vz[ix*n1+iz]);
+ field = vz[ix*n1+iz] - 0.5*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]));
+ field += vz[ix*n1+iz2] - 0.5*roz[ix*n1+iz2]*(
+ c1*(tzz[ix*n1+iz2] - tzz[ix*n1+iz2-1] +
+ txz[(ix+1)*n1+iz2] - txz[ix*n1+iz2]) +
+ c2*(tzz[ix*n1+iz2+1] - tzz[ix*n1+iz2-2] +
+ txz[(ix+2)*n1+iz2] - txz[(ix-1)*n1+iz2]));
+ rec_vz[irec*rec.nt+isam] += 0.5*field;
}
}
else {
rec_vz[irec*rec.nt+isam] += vz[ix*n1+iz2];
- //rec_vz[irec*rec.nt+isam] = vz[ix*n1+iz];
- //fprintf(stderr,"isam=%d vz[%d]=%e vz[%d]=%e vz[%d]=%e \n",isam, iz-1,vz[ix*n1+iz-1],iz,vz[ix*n1+iz], iz+1, vz[ix*n1+iz+1]);
}
}
if (rec.type.vx) {
@@ -268,7 +296,17 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
rec_vx[irec*rec.nt+isam] += 0.5*field;
}
else {
- rec_vx[irec*rec.nt+isam] += 0.5*(vx[ix2*n1+iz]+vx[ix*n1+iz]);
+ field = vx[ix*n1+iz] - 0.5*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]) );
+ field = vx[ix2*n1+iz] - 0.5*rox[ix2*n1+iz]*(
+ c1*(txx[ix2*n1+iz] - txx[(ix2-1)*n1+iz] +
+ txz[ix2*n1+iz+1] - txz[ix2*n1+iz]) +
+ c2*(txx[(ix2+1)*n1+iz] - txx[(ix2-2)*n1+iz] +
+ txz[ix2*n1+iz+2] - txz[ix2*n1+iz-1]) );
+ rec_vx[irec*rec.nt+isam] += 0.5*field;
}
}
else {
diff --git a/marchenko/applyMute_tshift.c b/marchenko/applyMute_tshift.c
index 26cd141e72f2d170a7bfa809dbaab1a641c2fd4b..8230fb831d767944e4f4d835e597119cfb97696c 100644
--- a/marchenko/applyMute_tshift.c
+++ b/marchenko/applyMute_tshift.c
@@ -219,3 +219,10 @@ void timeShift(float *data, int nsam, int nrec, float dt, float shift, float fmi
return;
}
+
+int map(int j, int nt)
+{
+ if (j<0) return nt-j;
+ else if (j>=nt) return j-nt;
+ else return j;
+}
diff --git a/marchenko/marchenko.c b/marchenko/marchenko.c
index 1021ecf7769a0caf2cb9bf8d79aed8e9fadd0d3f..508302fcb793b05b6cfff3a5d0ec857f3f936414 100644
--- a/marchenko/marchenko.c
+++ b/marchenko/marchenko.c
@@ -261,30 +261,10 @@ int main (int argc, char **argv)
tshift = fabs((nxs-1)*dxs*p);
/* compute mute window for plane waves */
- //for (i=0; i<nxs; i++) fprintf(stderr,"i=%d window=%d\n", i, muteW[i]);
- //itmin = nt;
- //for (i=0; i<nxs; i++) itmin = MIN (itmin, muteW[i]);
-
- /* for negative angles tshift is negative and */
- if (src_angle < 0.0) {
- itmin = NINT(tshift/dt);
- //for (i=0; i<nxs; i++) muteW[i] = MAX(0, muteW[i]-itmin);
- for (i=0; i<nxs; i++) muteW[i] = muteW[i]-itmin;
- timeShift(G_d, nts, nxs, dt, tshift, fmin, fmax);
- }
-
- for (i=0; i<nxs; i++) tsynW[i] = NINT(i*dxs*p/dt);
- //for (i=0; i<nxs; i++) tsynW[i] = 0.0;
+ for (i=0; i<nxs; i++) tsynW[i] = NINT((i-(nxs-1)/2)*dxs*p/dt);
if (Nfoc!=1) verr("For plane-wave focusing only one function can be computed at the same time");
- //fprintf(stderr,"itmin=%d tshift=%f =%d \n", itmin, tshift, NINT(tshift/dt));
- //for (i=0; i<nxs; i++) fprintf(stderr,"i=%d window=%f\n", i, tsynW[i]*dt);
-/* // TESTING SHIFT
- tshift=0.3;
- for (i=0; i<nxs; i++) tsynW[i] = NINT(0.3/dt);
-*/
}
else { /* just fill with zero's */
- //itmin=0;
for (i=0; i<nxs*Nfoc; i++) {
tsynW[i] = 0;
}
@@ -619,15 +599,7 @@ int main (int argc, char **argv)
}
/* Apply mute with window for Gmin */
if ( plane_wave==1 ) {
- /* for plane wave with angle tshift downward */
applyMute_tshift(Gmin, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, 1, tsynW);
- if (src_angle > 0.0) {
- if (verbose>1) vmess("Gmin planewave tshift=%f", tshift);
- timeShift(Gmin, nts, npos, dt, tshift, fmin, fmax);
- }
- else {
- if (verbose>1) vmess("Gmin NO planewave tshift");
- }
}
else {
applyMute(Gmin, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
diff --git a/marchenko/marchenko_primaries.c b/marchenko/marchenko_primaries.c
index 7a4d68aff0ec4047dd205245b9cad6aaefbe116e..90f3c815ffc19d004369d9d6714390ca0a91efe7 100644
--- a/marchenko/marchenko_primaries.c
+++ b/marchenko/marchenko_primaries.c
@@ -190,6 +190,7 @@ int main (int argc, char **argv)
if(!getparfloat("src_velo",&src_velo)) src_velo=1500.;
if(!getparfloat("t0",&tt0)) tt0=0.1;
if( (niterskip>1) && ISODD(niter) ) niter++;
+ if (niterskip==0) niterskip=1;
if (T>0) {
T=-1;