Skip to content
Snippets Groups Projects
Commit 5ff0629f authored by Jan Thorbecke's avatar Jan Thorbecke
Browse files

stable Vz source implementation in fdelmodc

parent 98a2bf8a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
fdelmodc/applySource.c
+ 8
− 6
View file @ 5ff0629f
......@@ -141,14 +141,16 @@ int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int i
/* Force source */
if (src.type == 6) {
// vx[ix*n1+iz] += src_ampl*(dt/mod.dx)/(l2m[ix*n1+iz]);
vx[ix*n1+iz] += src_ampl*rox[ix*n1+iz]/(l2m[ix*n1+iz]);
//vx[ix*n1+iz] += src_ampl*rox[ix*n1+iz]/(l2m[ix*n1+iz]);
/* stable implementation from "Numerical Techniques for Conservation Laws with Source Terms" by Justin Hudson */
vx[ix*n1+iz] = 0.5*(vx[(ix+1)*n1+iz]+vx[(ix-1)*n1+iz])+src_ampl*rox[ix*n1+iz]/(l2m[ix*n1+iz]);
}
else if (src.type == 7) {
/* old amplitude setting does not obey reciprocity */
// vz[ix*n1+iz] += src_ampl*(dt/mod.dx)/(l2m[ix*n1+iz]);
vz[ix*n1+iz] += src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
// fprintf(stderr,"ix=%d iz=%d l2m=%e rox=%e rho=%e dt=%e dx=%e\n", ix, iz, l2m[ix*n1+iz], roz[ix*n1+iz], dt/(mod.dx*roz[ix*n1+iz]), dt, mod.dx);
/* old amplitude setting does not obey reciprocity */
//vz[ix*n1+iz] += src_ampl*(dt/mod.dx)/(l2m[ix*n1+iz]);
//vz[ix*n1+iz] += src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
/* stable implementation from "Numerical Techniques for Conservation Laws with Source Terms" by Justin Hudson */
vz[ix*n1+iz] = 0.5*(vz[ix*n1+iz-1]+vz[ix*n1+iz+1])+src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
} /* src.type */
......
fdelmodc/defineSource.c
+ 18
− 6
View file @ 5ff0629f
......@@ -113,15 +113,14 @@ int defineSource(wavPar wav, srcPar src, modPar mod, float **src_nwav, int rever
optnscale = optn;
nfreqscale = optnscale/2 + 1;
}
// fprintf(stderr,"define S optn=%d ns=%d %e nt=%d %e\n", optn, wav.ns, wav.ds, optnscale, wav.dt);
ctrace = (complex *)calloc(nfreqscale,sizeof(complex));
trace = (float *)calloc(optnscale,sizeof(float));
fprintf(stderr,"define S optn=%d ns=%d %e nt=%d %e\n", optn, wav.ns, wav.ds, optnscale, wav.dt);
df = 1.0/(optn*wav.ds);
deltom = 2.*M_PI*df;
scl = 1.0/optn;
maxampl=0.0;
iwmax = nfreq;
for (i=0; i<wav.nx; i++) {
......@@ -180,16 +179,29 @@ int defineSource(wavPar wav, srcPar src, modPar mod, float **src_nwav, int rever
/* avoid a (small) spike in the last sample
this is done to avoid diffraction from last wavelet sample
which will act as a pulse */
maxampl=0.0;
if (reverse) {
for (j=0; j<wav.nt; j++) src_nwav[i][j] = scl*(trace[wav.nt-j-1]-trace[0]);
for (j=0; j<wav.nt; j++) {
src_nwav[i][j] = scl*(trace[wav.nt-j-1]-trace[0]);
maxampl = MAX(maxampl,fabs(src_nwav[i][j]));
}
}
else {
for (j=0; j<wav.nt; j++) src_nwav[i][j] = scl*(trace[j]-trace[wav.nt-1]);
for (j=0; j<wav.nt; j++) {
src_nwav[i][j] = scl*(trace[j]-trace[wav.nt-1]);
maxampl = MAX(maxampl,fabs(src_nwav[i][j]));
}
}
vmess("Wavelet sampling FFT-interpolated done for trace %d", i);
if (verbose > 3) vmess("Wavelet sampling (FFT-interpolated) done for trace %d", i);
}
}
/* set values smaller than 1e-5 maxampl to zero */
maxampl *= 1e-5;
for (i=0; i<wav.nx; i++) {
for (j=0; j<wav.nt; j++) {
if (fabs(src_nwav[i][j]) < maxampl) src_nwav[i][j] = 0.0;
}
}
free(ctrace);
free(trace);
......
fdelmodc/demo/fdelmodc_circ_medium.scr
+ 12
− 12
View file @ 5ff0629f
......@@ -7,11 +7,11 @@
# receivers are placed on a circle
export PATH=../../utils:$PATH:
export PATH=../../bin:$PATH:
dt=0.0001
#makewave file_out=wavelet.su dt=0.0020 nt=1024 fp=13 shift=1 w=g2 verbose=1
makewave w=fw fmin=0 flef=6 frig=89 fmax=95 dt=$dt file_out=wavefw.su nt=16384 t0=0.3 scale=0 scfft=1
dt=0.0008
makewave file_out=wave.su dt=$dt nt=1024 fp=25 shift=1 w=g2 verbose=1
#makewave w=fw fmin=0 flef=6 frig=89 fmax=95 dt=$dt file_out=wavefw.su nt=16384 t0=0.3 scale=0 scfft=1
#sufft < wavefw.su | suamp | sugain scale=0.0001 | suxgraph style=normal
makemod file_base=model.su \
......@@ -38,7 +38,7 @@ set -x
dxrcv=10 \
dtrcv=0.004 \
xsrc=0 zsrc=-2000 nshot=1 nsrc=1 \
src_type=1 tmod=4.092 \
src_type=1 tmod=1.020 \
ntaper=100 \
left=2 right=2 bottom=2 top=2
......@@ -46,17 +46,17 @@ export filecp=hom_cp.su
export filero=hom_ro.su
file_rcv=hom.su
../fdelmodc \
fdelmodc \
file_cp=$filecp file_den=$filero \
ischeme=1 \
file_src=wavefw.su verbose=2 \
src_type=7 \
file_src=wave.su verbose=2 \
file_rcv=$file_rcv \
rec_type_vz=0 rec_type_p=1 \
xrcv1=-2000 xrcv2=2000 zrcv1=-2000 zrcv2=-2000 \
dxrcv=10 \
rec_type_vz=1 rec_type_p=1 \
rrcv=500 dphi=1 \
dtrcv=0.004 \
xsrc=0 zsrc=-2000 nshot=1 nsrc=1 \
src_type=1 tmod=4.092 \
xsrc=0 zsrc=0 nshot=1 nsrc=1 \
tmod=1.020 \
ntaper=100 \
left=2 right=2 bottom=2 top=2
......
fdelmodc/demo/modelOilGas.scr
+ 12
− 11
View file @ 5ff0629f
......@@ -7,47 +7,48 @@ rho=2500
dx=2.5
dt=0.0005
makemod sizex=5000 sizez=2500 dx=$dx dz=$dx cp0=$cp ro0=$rho \
orig=-2500,0 file_base=syncl.su \
intt=def x=-2500,0,2500 z=250,250,250 poly=0 cp=2300 ro=2000 \
intt=def x=-2500,-2000,-1000,-800,0,800,2500 z=650,650,700,750,900,750,600 poly=2 cp=2600 ro=2500 \
intt=def x=-2500,0,2500 z=1390,1390,1390 poly=0 cp=2000 ro=2000
makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
makewave w=g1 fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
export OMP_NUM_THREADS=2
zsrc=0
zsrc=1100
zsrc=15
../fdelmodc \
fdelmodc \
file_cp=syncl_cp.su ischeme=1 iorder=4 \
file_den=syncl_ro.su \
file_src=wave.su \
file_rcv=shot_fd.su \
src_type=7 \
src_type=6 \
src_orient=1 \
src_injectionrate=1 \
src_injectionrate=0 \
rec_type_vz=1 \
rec_type_vx=1 \
rec_type_p=1 \
rec_int_vz=2 \
dtrcv=0.004 \
rec_delay=0.1 \
verbose=2 \
tmod=2.55 \
verbose=4 \
tmod=2.0 \
dxrcv=10.0 \
xrcv1=-2250 xrcv2=2250 \
zrcv1=0 zrcv2=0 \
xsrc=0 zsrc=$zsrc \
file_snap=snapF_$zsrc \
tsnap1=0.1 tsnap2=4.0 dtsnap=0.05 dxsnap=$dx dzsnap=$dx \
ntaper=101 \
snapwithbnd=1 \
left=1 right=1 top=1 bottom=1
file_snap=snapF_$zsrc \
tsnap1=0.1 tsnap2=4.0 dtsnap=0.05 dxsnap=$dx dzsnap=$dx \
left=2 right=2 top=2 bottom=2
exit
makemod sizex=5000 sizez=2500 dx=$dx dz=$dx cp0=$cp ro0=$rho \
orig=-2500,0 file_base=hom.su
......
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