Skip to content
Snippets Groups Projects
Commit 0b299db6 authored by JanThorbecke's avatar JanThorbecke
Browse files

marchenko_primaries added, memory leak in fdacrtm

parent e8ae2f41
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
fdacrtmc/extractMigrationSnapshots.c
+ 3
− 5
View file @ 0b299db6
......@@ -60,7 +60,6 @@ int storeVerticalParticleVelocitySnapshot(modPar *mod, wavPar *wav, migPar *mig)
else farr=(float*)malloc(mig->sizem*sizeof(float));
// Store Vertical Particle Velocity Field - Interpolate to Tzz(P)-Grid
mig->wav[mig->it].vz=(float*)malloc(mig->sizem*sizeof(float));
for(ix=0,ix1=mod->ioZx+mig->x1;ix<mig->nx;ix++,ix1+=mig->skipdx){
for(iz=0,iz1=mod->ioZz+mig->z1;iz<mig->nz;iz++,iz1+=mig->skipdz){
farr[ix*mig->nz+iz]=(wav->vz[ix1*mod->naz+iz1]+wav->vz[ix1*mod->naz+iz1+1])/2.0;
......@@ -77,9 +76,6 @@ int storeVerticalParticleVelocitySnapshot(modPar *mod, wavPar *wav, migPar *mig)
int extractMigrationSnapshots(modPar *mod, wavPar *wav, migPar *mig, decompPar *decomp){
size_t ix, ix1, iz, iz1;
//if(mig->it<100){mig->it++;return(0);}
switch(mig->mode){
case 1: /* Conventional Migration - Store Pressure Field (Tzz) */
// P (Tzz)
......@@ -90,8 +86,10 @@ int extractMigrationSnapshots(modPar *mod, wavPar *wav, migPar *mig, decompPar *
storePressureSnapshot(mod,wav,mig);
switch(mig->orient){ //Migration Orientation
case 0: //Image Wavefields not travelling in the same direction
// Vertical Particle Velocity - Interpolate to Tzz(P)-Grid
storeVerticalParticleVelocitySnapshot(mod,wav,mig);
// Horizontal Particle Velocity - Interpolate to Tzz(P)-Grid
storeVerticalParticleVelocitySnapshot(mod,wav,mig); //TODO: What does this case actually do?
storeHorizontalParticleVelocitySnapshot(mod,wav,mig); //TODO: What does this case actually do?
break;
case 1: //Up-Down Imaging
// Vertical Particle Velocity - Interpolate to Tzz(P)-Grid
......
fdacrtmc/rtmImagingCondition.c
+ 4
− 2
View file @ 0b299db6
......@@ -66,6 +66,8 @@ int rtmImagingCondition(modPar *mod, wavPar *wav, migPar *mig, decompPar *decomp
}
}
}
free(mig->wav[mig->it].vx);
free(mig->wav[mig->it].vz);
break;
case 1: // Up-Down Imaging
for(ix=0,ix1=mod->ioPx+mig->x1;ix<mig->nx;ix++,ix1+=mig->skipdx){
......@@ -78,6 +80,7 @@ int rtmImagingCondition(modPar *mod, wavPar *wav, migPar *mig, decompPar *decomp
if(pzf*pzb>0) mig->image[ix*mig->nz+iz]+=mig->dt*(wav->tzz[ix1*mod->naz+iz1]*mig->wav[mig->it].tzz[ix*mig->nz+iz]);
}
}
free(mig->wav[mig->it].vz);
break;
case 2: // Left-Right Imaging
for(ix=0,ix1=mod->ioPx+mig->x1;ix<mig->nx;ix++,ix1+=mig->skipdx){
......@@ -90,6 +93,7 @@ int rtmImagingCondition(modPar *mod, wavPar *wav, migPar *mig, decompPar *decomp
if(pxf*pxb<0) mig->image[ix*mig->nz+iz]+=mig->dt*(wav->tzz[ix1*mod->naz+iz1]*mig->wav[mig->it].tzz[ix*mig->nz+iz]);
}
}
free(mig->wav[mig->it].vx);
break;
case 3: //Normal Imaging
//Not yet implemented
......@@ -103,8 +107,6 @@ int rtmImagingCondition(modPar *mod, wavPar *wav, migPar *mig, decompPar *decomp
break;
}
free(mig->wav[mig->it].tzz);
free(mig->wav[mig->it].vx);
free(mig->wav[mig->it].vz);
break;
case 3: /* Wavefield decomposition zero-lag pressure cross-correlation */
mig->it--;break;
......
marchenko/Makefile
+ 1
− 0
View file @ 0b299db6
......@@ -4,6 +4,7 @@ include ../Make_include
#LIBS += -L$L -lgenfft -lm $(LIBSM)
#OPTC += -g -O0 -Wall
#OPTC += -g -traceback #-check-pointers=rw -check-pointers-undimensioned
ALL: fmute marchenko marchenko_primaries
......
marchenko/demo/oneD/primaries.scr
+ 6
− 4
View file @ 0b299db6
......@@ -8,16 +8,18 @@
#cd $SLURM_SUBMIT_DIR
export PATH=$HOME/src/OpenSource/bin:$PATH:
export OMP_NUM_THREADS=20
export OMP_NUM_THREADS=40
#shot record to remove internal multiples
select=451
select=51
makewave fp=20 dt=0.004 file_out=wave.su nt=2048 t0=0.0 scale=0 scfft=1
../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
nshots=601 verbose=2 istart=40 iend=500 fmax=90 pad=1024 \
niter=25 smooth=10 niterskip=600 shift=20 file_rr=pred_rr.su T=0
nshots=901 verbose=2 istart=40 iend=500 fmax=90 pad=1024 \
niter=20 smooth=10 niterskip=1 shift=20 file_rr=pred_rr.su T=0
exit;
#for reference original shot record from Reflection matrix
suwind key=fldr min=$select max=$select < shotsdx5_rp.su itmax=2047 > shotsx0.su
......
marchenko/marchenko_primaries.c
+ 18
− 12
View file @ 0b299db6
......@@ -187,6 +187,8 @@ int main (int argc, char **argv)
if (!getparfloat("dt", &dt)) dt = d1;
if(!getparint("istart", &istart)) istart=20;
if(!getparint("iend", &iend)) iend=nt;
iend = MIN(iend, nt-shift-1);
istart = MIN(MAX(1,istart),iend);
if (file_tinv == NULL) {/* 'G_d' is one of the shot records */
if(!getparint("ishot", &ishot)) ishot=1+(nshots-1)/2;
......@@ -513,10 +515,10 @@ int main (int argc, char **argv)
G_d[l*nxs*nts+i*nts+j] = DD[l*nxs*nts+i*nts+j];
}
/* apply mute window for samples above nts-ii */
for (j = 0; j < nts-ii+T*shift-smooth; j++) {
for (j = 0; j < MIN(nts, nts-ii+T*shift-smooth); j++) {
G_d[l*nxs*nts+i*nts+j] = 0.0;
}
for (j = nts-ii+T*shift-smooth, k=1; j < nts-ii+T*shift; j++, k++) {
for (j = nts-ii+T*shift-smooth, k=1; j < MIN(nts, nts-ii+T*shift); j++, k++) {
G_d[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
}
}
......@@ -533,7 +535,9 @@ int main (int argc, char **argv)
//}
}
}
writeDataIter("G_d.su", G_d, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii);
if (file_iter != NULL) {
writeDataIter("G_d.su", G_d, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii);
}
}
else { /* use f1min from previous iteration as starting point and do niterec iterations */
niterrun=niterec;
......@@ -548,10 +552,10 @@ int main (int argc, char **argv)
G_d[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j] - f1min[l*nxs*nts+i*nts+nts-j];
}
/* apply mute window for samples above nts-ii */
for (j = 0; j < nts-ii+T*shift-smooth; j++) {
for (j = 0; j < MIN(nts,nts-ii+T*shift-smooth); j++) {
G_d[l*nxs*nts+i*nts+j] = 0.0;
}
for (j = nts-ii+T*shift-smooth, k=1; j < nts-ii+T*shift; j++, k++) {
for (j = nts-ii+T*shift-smooth, k=1; j < MIN(nts, nts-ii+T*shift); j++, k++) {
G_d[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
}
}
......@@ -597,17 +601,17 @@ int main (int argc, char **argv)
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) {
/* apply mute window for samples after ii */
for (j = ii-T*shift+smooth; j < nts; j++) {
for (j = MAX(0,ii-T*shift+smooth); j < nts; j++) {
Ni[l*nxs*nts+i*nts+j] = 0.0;
}
for (j = ii-T*shift+smooth, k=0; j < ii; j++, k++) {
for (j = MAX(0,ii-T*shift+smooth), k=0; j < ii; j++, k++) {
Ni[l*nxs*nts+i*nts+j] *= costaper[k];
}
/* apply mute window for delta function at t=0*/
for (j = 0; j < shift-smooth; j++) {
Ni[l*nxs*nts+i*nts+j] = 0.0;
}
for (j = shift-smooth, k=1; j < shift; j++, k++) {
for (j = MAX(0,shift-smooth), k=1; j < shift; j++, k++) {
Ni[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
}
f1plus[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
......@@ -616,9 +620,11 @@ int main (int argc, char **argv)
}
}
}
/*
if (file_iter != NULL) {
writeDataIter("f1plus.su", f1plus, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter);
}
*/
}
else {/* odd iterations: => f_1^-(t) */
for (l = 0; l < Nfoc; l++) {
......@@ -645,7 +651,7 @@ int main (int argc, char **argv)
for (j = nts-shift+smooth; j < nts; j++) {
Ni[l*nxs*nts+i*nts+j] = 0.0;
}
for (j = nts-shift, k=0; j < nts-shift+smooth; j++, k++) {
for (j = nts-shift, k=0; j < MIN(nts, nts-shift+smooth); j++, k++) {
Ni[l*nxs*nts+i*nts+j] *= costaper[k];
}
/* apply mute window for samples above nts-ii */
......@@ -654,10 +660,10 @@ int main (int argc, char **argv)
// Ni[l*nxs*nts+i*nts+j] = 0.0;
//}
/* apply mute window for samples above nts-ii */
for (j = 0; j < nts-ii+T*shift-smooth; j++) {
for (j = 0; j < MIN(nts,nts-ii+T*shift-smooth); j++) {
Ni[l*nxs*nts+i*nts+j] = 0.0;
}
for (j = nts-ii+T*shift-smooth, k=1; j < nts-ii+T*shift; j++, k++) {
for (j = nts-ii+T*shift-smooth, k=1; j < MIN(nts,nts-ii+T*shift); j++, k++) {
Ni[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
}
}
......@@ -689,7 +695,7 @@ int main (int argc, char **argv)
if((ii-istart)==10)t4=wallclock_time();
if((ii-istart)==20){
t4=(wallclock_time()-t4)*((iend-istart)/10.0);
fprintf(stderr,"\r %s: Estimated total compute time = %.2fs.\n %s: Progress: %3d%%",xargv[0],t4,xargv[0],(ii-istart)/((iend-istart)/100));
fprintf(stderr,"\r %s: Estimated total compute time = %.2fs.\n %s: Progress: %3d%%",xargv[0],t4,xargv[0],(ii-istart)/((iend-istart)/100.0));
}
//t4=wallclock_time();
tii=(t4-t1)*((float)(iend-istart)/(ii-istart+1.0))-(t4-t1);
......
marchenko/synthesis.c
+ 2
− 2
View file @ 0b299db6
......@@ -95,7 +95,7 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
/* 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));
memset(&Fop[l*nxs*nw], 0, nxs*nw*sizeof(complex));
for (i = 0; i < npos; i++) {
rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
ix = ixpos[i];
......@@ -111,7 +111,7 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
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));
memset(&Fop[l*nxs*nw], 0, nxs*nw*sizeof(complex));
for (i = 0; i < nxs; i++) {
rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
for (iw=0; iw<nw; iw++) {
......
utils/syn2d.c
+ 2
− 2
View file @ 0b299db6
......@@ -242,7 +242,6 @@ int main (int argc, char **argv)
if (verbose>=2) vmess("dimensions used: %d x %d",ntmax,nxmax);
}
if (!getparfloat("dt", &dt)) dt = d1;
fprintf(stderr,"dt=%e\n", dt);
size = ntmax * nxmax;
xrcv = (float *)malloc(nxmax*sizeof(float));
......@@ -501,6 +500,7 @@ int main (int argc, char **argv)
for (i = 0; i < n2; i++) {
hdrs_out[i].ns = n1;
hdrs_out[i].dt = (unsigned short)(dt*1000000);
hdrs_out[i].trid = hdrs_in[i].trid;
hdrs_out[i].f1 = f1;
hdrs_out[i].f2 = f2;
......@@ -510,7 +510,7 @@ int main (int argc, char **argv)
hdrs_out[i].trwf = n2out;
hdrs_out[i].scalco = -1000;
hdrs_out[i].sx = NINT(xsyn[l]*1000);
hdrs_out[i].gx = NINT(1000*(f2+i*d2));
hdrs_out[i].gx = NINT(1000.0*(f2+i*d2));
hdrs_out[i].offset = (long)NINT((f2+i*d2) - xsyn[l]);
hdrs_out[i].scalel = -1000;
hdrs_out[i].selev = NINT(zsyn[l]*1000);
......
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