Skip to content
Snippets Groups Projects
Commit 67742d01 authored by JanThorbecke's avatar JanThorbecke
Browse files

plane waves FD for negative times

parent 44281295
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
MDD/Makefile
+ 1
− 1
View file @ 67742d01
......@@ -16,7 +16,7 @@ ALLINC = -I.
#FFLAGS += -g -mcmodel=medium
#LIBS += -lblas -llapack -L$L -lgenfft $(LIBSM) -lc -lm
all: mdd
all: mdd
PRG = mdd
......
MDD/deconvolve.c
+ 5
− 3
View file @ 67742d01
......@@ -131,7 +131,7 @@ int deconvolve(complex *cA, complex *cB, complex *cC, complex *oBB, int nfreq, i
if (conjgA) NC = nshots;
else NC = nstationB;
// if (verbose) fprintf(stderr,"transa=%s transb=%s %d %d %d\n", transa, transb, NA, NB, nshots);
if (verbose>1) vmess("MMD tranpose options transa=%s transb=%s N(rec)A=%d N(rec)B=%d LDA(nshots)=%d\n", transa, transb, NA, NB, nshots);
#pragma omp for schedule(static) \
private(iw, iwnA, iwnB, iwAB, iwBB)
......@@ -141,10 +141,12 @@ private(iw, iwnA, iwnB, iwAB, iwBB)
iwnB = iw*nstationB*nshots;
iwAB = iw*NC*NC;
if (mdd==0) { /* Correlation */
/* cblas_cgemm(CblasRowMajor,CblasNoTrans, CblasConjTrans, NA, NB, nshots, &alpha.r,
/*
cblas_cgemm(CblasColMajor,CblasNoTrans, CblasNoTrans, NA, NB, nshots, &alpha.r,
&cA[iwnA].r, NA,
&cB[iwnB].r, NB, &beta.r,
&cC[iwAB].r, NC); */
&cC[iwAB].r, NC);
*/
cgemm_(transa, transb, &NA, &NB, &nshots, &alpha.r,
&cA[iwnA].r, &NA,
&cB[iwnB].r, &NB, &beta.r,
......
fdelmodc/demo/compare_green_inmanual.scr
+ 0
− 1
View file @ 67742d01
......@@ -127,7 +127,6 @@ sudiff trace_green.su trace_fd.su | basop choice=shift shift=-0.1 | sugain scale
# rvz
(suwind key=tracl min=101 max=101 < shot_fd_dip_rvz.su | basop choice=shift shift=-$shift; suwind key=tracl min=101 max=101 < shot_green_dip_rvz.su ) | basop choice=shift shift=-0.1 | supsgraph style=normal labelsize=10 wbox=4 hbox=2 label1='time in seconds' label2="Amplitude" linewidth=0.1 linecolor=red,green d2num=6e-10 f2num=-1.2e-9 x2end=2e-9 x2beg=-1.3e-9 > dip_rvz.eps
(suwind key=tracl min=101 max=101 < shot_fd_dip_rvz.su | basop choice=shift shift=-$shift; suwind key=tracl min=101 max=101 < shot_green_dip_rvz.su ) | basop choice=shift shift=-0.1 | sugain scale=-1 | supsgraph style=normal labelsize=10 wbox=2 hbox=2 label1='time in seconds' label2="Amplitude" linewidth=0.1 linecolor=red,green d2num=5e-9 f2num=9.5e-7 x2end=9.78e-7 x2beg=9.5e-7 x1beg=0.2430 x1end=0.24601 f1=0 f1num=0.243 d1num=0.002 > dip_zoom_rvz.eps
......
fdelmodc/fdelmodc.c
+ 4
− 8
View file @ 67742d01
......@@ -521,12 +521,12 @@ int main(int argc, char **argv)
its=-1;
it0=0;
it1=mod.nt;
it1=mod.nt+NINT(-mod.t0/mod.dt);
its=1;
}
else {
it0=0;
it1=mod.nt;
it1=mod.nt+NINT(-mod.t0/mod.dt);
its=1;
}
perc=it1/100;if(!perc)perc=1;
......@@ -607,20 +607,18 @@ 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)/rec.skipdt)+1)%rec.nt );
writeToFile = ! ( (((it-rec.delay+NINT(mod.t0/mod.dt))/rec.skipdt)+1)%rec.nt );
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 */
//isam = (it-rec.delay-itwritten)/rec.skipdt+1;
/* negative time correction with mod.t0 for dipping plane waves modeling */
isam = (it-rec.delay-itwritten+NINT(mod.t0/mod.dt))/rec.skipdt+1;
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,
rec_vx, rec_vz, rec_txx, rec_tzz, rec_txz,
rec_p, rec_pp, rec_ss, rec_udp, rec_udvz, verbose);
/* at the end of modeling a shot, write receiver array to output file(s) */
if (writeToFile && (it+rec.skipdt <= it1-1) ) {
fileno = ( ((it-rec.delay)/rec.skipdt)+1)/rec.nt;
......@@ -686,8 +684,6 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, izsrc, it, src_nwav, verbose)
}
}
}
/* negative time correction with mod.t0 for dipping plane waves modeling */
isam += NINT(-mod.t0/mod.dt)/rec.skipdt;
writeRec(rec, mod, bnd, wav, ixsrc, izsrc, isam+1, ishot, fileno,
rec_vx, rec_vz, rec_txx, rec_tzz, rec_txz,
rec_p, rec_pp, rec_ss, rec_udp, rec_udvz, verbose);
......
fdelmodc/getParameters.c
+ 2
− 1
View file @ 67742d01
......@@ -884,7 +884,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
src_iz0=MAX(0,NINT((zsrc-sub_z0)/dz));
src_iz0=MIN(src_iz0,nz);
/* for a source defined on mutliple gridpoint calculate p delay factor */
/* plane-wave defined on multiple gridpoints calculate p delay factor */
src->x = (int *)malloc(nsrc*sizeof(int));
src->z = (int *)malloc(nsrc*sizeof(int));
......@@ -892,6 +892,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
src->tend = (float *)malloc(nsrc*sizeof(float));
grad2rad = 17.453292e-3;
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;
}
......
fdelmodc/getRecTimes.c
+ 29
− 29
View file @ 67742d01
......@@ -86,7 +86,7 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
C10 = tzz[(ix+1)*n1+iz];
C01 = tzz[ix*n1+iz+1];
C11 = tzz[(ix+1)*n1+iz+1];
rec_p[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
rec_p[irec*rec.nt+isam] += C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
C01*(1.0-rdx)*rdz + C11*rdx*rdz;
}
if (rec.type.txx) {
......@@ -94,7 +94,7 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
C10 = txx[(ix+1)*n1+iz];
C01 = txx[ix*n1+iz+1];
C11 = txx[(ix+1)*n1+iz+1];
rec_txx[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
rec_txx[irec*rec.nt+isam] += C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
C01*(1.0-rdx)*rdz + C11*rdx*rdz;
}
if (rec.type.tzz) {
......@@ -102,7 +102,7 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
C10 = tzz[(ix+1)*n1+iz];
C01 = tzz[ix*n1+iz+1];
C11 = tzz[(ix+1)*n1+iz+1];
rec_tzz[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
rec_tzz[irec*rec.nt+isam] += C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
C01*(1.0-rdx)*rdz + C11*rdx*rdz;
}
if (rec.type.txz) {
......@@ -110,7 +110,7 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
C10 = txz[(ix2+1)*n1+iz2];
C01 = txz[ix2*n1+iz2+1];
C11 = txz[(ix2+1)*n1+iz2+1];
rec_txz[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
rec_txz[irec*rec.nt+isam] += C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
C01*(1.0-rdx)*rdz + C11*rdx*rdz;
}
if (rec.type.pp) {
......@@ -122,7 +122,7 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
vz[ix*n1+iz2+1]-vz[ix*n1+iz+1])/mod.dx;
C11 = (vx[(ix2+1)*n1+iz+1]-vx[(ix+1)*n1+iz+1] +
vz[(ix+1)*n1+iz2+1]-vz[(ix+1)*n1+iz+1])/mod.dx;
rec_pp[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
rec_pp[irec*rec.nt+isam] += C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
C01*(1.0-rdx)*rdz + C11*rdx*rdz;
}
if (rec.type.ss) {
......@@ -134,7 +134,7 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
(vz[ix2*n1+iz2+1]-vz[ix*n1+iz2+1]))/mod.dx;;
C11 = (vx[(ix2+1)*n1+iz2+1]-vx[(ix2+1)*n1+iz+1] -
(vz[(ix2+1)*n1+iz2+1]-vz[(ix+1)*n1+iz2+1]))/mod.dx;
rec_ss[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
rec_ss[irec*rec.nt+isam] += C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
C01*(1.0-rdx)*rdz + C11*rdx*rdz;
}
if (rec.type.vz) {
......@@ -142,7 +142,7 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
C10 = vz[(ix+1)*n1+iz2];
C01 = vz[ix*n1+iz2+1];
C11 = vz[(ix+1)*n1+iz2+1];
rec_vz[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
rec_vz[irec*rec.nt+isam] += C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
C01*(1.0-rdx)*rdz + C11*rdx*rdz;
}
if (rec.type.vx) {
......@@ -150,7 +150,7 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
C10 = vx[(ix2+1)*n1+iz];
C01 = vx[ix2*n1+iz+1];
C11 = vx[(ix2+1)*n1+iz+1];
rec_vx[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
rec_vx[irec*rec.nt+isam] += C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
C01*(1.0-rdx)*rdz + C11*rdx*rdz;
}
......@@ -173,10 +173,10 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
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]*(dvx+dvz);
rec_p[irec*rec.nt+isam] = 0.5*field;
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]);
rec_p[irec*rec.nt+isam] += 0.5*(tzz[ix*n1+iz1]+tzz[ix*n1+iz]);
}
}
else if (rec.int_p == 2) {
......@@ -191,40 +191,40 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
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]*(dvx+dvz);
rec_p[irec*rec.nt+isam] = 0.5*field;
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]);
rec_p[irec*rec.nt+isam] += 0.5*(tzz[ix1*n1+iz]+tzz[ix*n1+iz]);
}
}
else {
rec_p[irec*rec.nt+isam] = tzz[ix*n1+iz];
rec_p[irec*rec.nt+isam] += tzz[ix*n1+iz];
}
}
if (rec.type.txx) rec_txx[irec*rec.nt+isam] = txx[ix*n1+iz];
if (rec.type.tzz) rec_tzz[irec*rec.nt+isam] = tzz[ix*n1+iz];
if (rec.type.txx) rec_txx[irec*rec.nt+isam] += txx[ix*n1+iz];
if (rec.type.tzz) rec_tzz[irec*rec.nt+isam] += tzz[ix*n1+iz];
if (rec.type.txz) { /* time interpolation to be done */
if (rec.int_vz == 2 || rec.int_vx == 2) {
rec_txz[irec*rec.nt+isam] = 0.25*(
rec_txz[irec*rec.nt+isam] += 0.25*(
txz[ix*n1+iz2]+txz[ix2*n1+iz2]+
txz[ix*n1+iz]+txz[ix2*n1+iz]);
}
else {
rec_txz[irec*rec.nt+isam] = txz[ix2*n1+iz2];
rec_txz[irec*rec.nt+isam] += txz[ix2*n1+iz2];
}
}
if (rec.type.pp) {
rec_pp[irec*rec.nt+isam] = (vx[ix2*n1+iz]-vx[ix*n1+iz] +
rec_pp[irec*rec.nt+isam] += (vx[ix2*n1+iz]-vx[ix*n1+iz] +
vz[ix*n1+iz2]-vz[ix*n1+iz])/mod.dx;
}
if (rec.type.ss) {
rec_ss[irec*rec.nt+isam] = (vx[ix2*n1+iz2]-vx[ix2*n1+iz] -
rec_ss[irec*rec.nt+isam] += (vx[ix2*n1+iz2]-vx[ix2*n1+iz] -
(vz[ix2*n1+iz2]-vz[ix*n1+iz2]))/mod.dx;
}
if (rec.type.vz) {
/* interpolate vz to vx position to the right and above of vz */
if (rec.int_vz == 1) {
rec_vz[irec*rec.nt+isam] = 0.25*(
rec_vz[irec*rec.nt+isam] += 0.25*(
vz[ix*n1+iz2]+vz[ix1*n1+iz2]+
vz[ix*n1+iz] +vz[ix1*n1+iz]);
}
......@@ -237,14 +237,14 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
field += vz[ix*n1+iz2] - 0.5*roz[ix*n1+iz2]*(
c1*(tzz[ix*n1+iz2] - tzz[ix*n1+iz2-1]) +
c2*(tzz[ix*n1+iz2+1] - tzz[ix*n1+iz2-2]));
rec_vz[irec*rec.nt+isam] = 0.5*field;
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]);
rec_vz[irec*rec.nt+isam] += 0.5*(vz[ix*n1+iz2]+vz[ix*n1+iz]);
}
}
else {
rec_vz[irec*rec.nt+isam] = vz[ix*n1+iz2];
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]);
}
......@@ -252,7 +252,7 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
if (rec.type.vx) {
/* interpolate vx to vz position to the left and below of vx */
if (rec.int_vx == 1) {
rec_vx[irec*rec.nt+isam] = 0.25*(
rec_vx[irec*rec.nt+isam] += 0.25*(
vx[ix2*n1+iz]+vx[ix2*n1+iz1]+
vx[ix*n1+iz]+vx[ix*n1+iz1]);
}
......@@ -265,14 +265,14 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
field += vx[ix2*n1+iz] - 0.5*rox[ix2*n1+iz]*(
c1*(tzz[ix2*n1+iz] - tzz[(ix2-1)*n1+iz]) +
c2*(tzz[(ix2+1)*n1+iz] - tzz[(ix2-2)*n1+iz]));
rec_vx[irec*rec.nt+isam] = 0.5*field;
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]);
rec_vx[irec*rec.nt+isam] += 0.5*(vx[ix2*n1+iz]+vx[ix*n1+iz]);
}
}
else {
rec_vx[irec*rec.nt+isam] = vx[ix2*n1+iz];
rec_vx[irec*rec.nt+isam] += vx[ix2*n1+iz];
}
}
}
......@@ -297,8 +297,8 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
}
for (ix=0; ix<mod.nax; ix++) {
/* -2- compute average in time and depth to get Vz at same depth and time as P */
rec_udvz[ix*rec.nt+isam] = 0.25*(vz[ix*n1+iz2]+vz[ix*n1+iz]+vz_t[mod.nax+ix]+vz_t[ix]);
rec_udp[ix*rec.nt+isam] = tzz[ix*n1+iz];
rec_udvz[ix*rec.nt+isam] += 0.25*(vz[ix*n1+iz2]+vz[ix*n1+iz]+vz_t[mod.nax+ix]+vz_t[ix]);
rec_udp[ix*rec.nt+isam] += tzz[ix*n1+iz];
}
free(vz_t);
}
......
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