Commit d3d7b0ac by Florian Wittkamp

Merge branch 'release/Release_2.0.1'

parents 0e8335b9 2414c395
 ... ... @@ -120,37 +120,38 @@ Default values are: Five built-in wavelets of the seismic source are available. The corresponding source time functions are defined in \texttt{src/wavelet.c}. You may modify the time functions in this file and recompile to include your own analytical wavelet or to modify the shape of the built-in wavelets. \newline Ricker wavelet (SOURCE\_SHAPE=1): SOURCE\_SHAPE=1, Ricker wavelet: r(\tau)=\left(1-2\tau^2\right)\exp(-\tau^2) \quad \mbox{with} \quad \tau=\frac{\pi(t-1.5/f_c-t_d)}{1.0/f_c} \label{eq_ricker} Fuchs-M\"uller wavelet (SOURCE\_SHAPE=2): SOURCE\_SHAPE=2, Fuchs-M\"uller wavelet: f_m(t)=\sin(2\pi(t-t_d)f_c)-0.5\sin(4\pi(t-t_d)f_c) \quad \mbox{if} \quad t\in[t_d,t_d+1/fc] \quad \mbox{else} \quad fm(t)=0 \label{eq_fm} $sin^3$ wavelet (SOURCE\_SHAPE=4): SOURCE\_SHAPE=4, $sin^3$ wavelet: s3(t)=0.75 \pi f_c \sin(\pi(t+t_d)f_c)^3\quad \mbox{if} \quad t \in[t_d,t_d+1/fc] \quad \mbox{else} \quad s3(t)=0 \label{eq_s3} First derivative of a Gaussian function (SOURCE\_SHAPE=5): SOURCE\_SHAPE=5, First derivative of a Gaussian function: f(t)= -2.0 a (t-t_s) \exp(-a (t-t_s)^2)\quad \mbox{with} \quad a=\pi^2 f_c^2 \quad \mbox{and} \quad t_s=1.2/f_c \label{eq_deriv_of_gaussian} Delta pulse (SOURCE\_SHAPE=6): Lowpass filtered delta pulse. Note, that it is not clear if the lowpass filter used in the current version works correctly for a delta pulse.\\ SOURCE\_SHAPE=6, delta pulse: Lowpass filtered delta pulse. Note, that it is not clear if the lowpass filter used in the current version works correctly for a delta pulse.\\ Source time function from SIGNAL\_FILE in su format (SOURCE\_SHAPE=7).\\ % Source time function from SIGNAL\_FILE in su format (SOURCE\_SHAPE=7).\\ In these equations, t denotes time and $f_c=1/TS$ is the center frequency. $t_d$ is a time delay which can be defined for each source position in SOURCE\_FILE. Note that the symmetric (zero phase) Ricker signal is always delayed by $1.0/f_c$, which means that after one period the maximum amplitude is excited at the source location. Three of these 5 source wavelets and the corresponding amplitude spectra for a center frequency of $f_c=50$ Hz and a delay of $t_d=0$ are plotted in Figure \ref{fig_source_wavelets_json}. Note the delay of the Ricker signal described above. The Fuchs-M\"uller wavelet has a slightly higher center frequency and covers a broader frequency range. In these equations, t denotes time and $f_c$ is the center frequency. $t_d$ is a time delay which can be defined for each source position. Note that the symmetric (zero phase) Ricker signal is always delayed by $1.0/f_c$, which means that after one period the maximum amplitude is excited at the source location. Three of these 5 source wavelets and the corresponding amplitude spectra for a center frequency of $f_c=50$ Hz and $t_d=0$ are plotted in Figure \ref{fig_source_wavelets_json}. Note the delay of the Ricker signal described above. The Fuchs-M\"uller wavelet has a slightly higher center frequency and covers a broader frequency range. \newline \begin{figure} \begin{center} ... ... @@ -162,8 +163,9 @@ spectrum, c) phase spectrum. } \label{fig_source_wavelets_json} \end{figure} You may also use your own time function as the source wavelet (for instance the signal of the first arrival recorded by a geophone at near offsets). Specify SOURCE\_SHAPE=3 and save the samples of your source wavelet in ASCII-format in SIGNAL\_FILE. SIGNAL\_FILE should contain one sample per line. It should thus look like: \newpage SOURCE\_SHAPE=3 allows you to use your own time function as the source wavelet stored in ASCII-format in SIGNAL\_FILE. SIGNAL\_FILE should then contain one sample per line. It should thus look like: {\color{blue}{\begin{verbatim} 0.0 ... ... @@ -175,12 +177,17 @@ your source wavelet in ASCII-format in SIGNAL\_FILE. SIGNAL\_FILE should contain The time interval between the samples must equal the time step interval (DT) of the FD simulation (see above)! Therefore it might be necessary to resample/interpolate a given source time function with a smaller sample rate. You may use the matlab script mfiles/resamp.m to resample your external source signal to the required sampling interval. \newline It is also possible to read different external source wavelets for each shot. Specify SOURCE\_SHAPE=7 and save the wavelets in su format in SIGNAL\_FILE.shot. The wavelets in each su file must equal the time step intervel (DT) and the number of time steps of the FD simulation! SOURCE\_SHAPE=7 is used for reading different external source wavelets for each shot. The wavelets in SU-format need to be saved in SIGNAL\_FILE.shot.su. If you want to use the source time function inversion (INV\_STF==1, section \ref{sec:STF}) with an external wavelet, this wavelet needs to be provided at SIGNAL\_FILE.shot\_start.su and the inverted wavelets will be stored in SIGNAL\_FILE.shot.su (or SIGNAL\_FILE...su if you use the WORKFLOW option). The wavelets in each su file must have the same number of samples as specified by TIME/DT! \newline The following source types are availabe: explosive sources that excite compressional waves only (SOURCE\_TYPE=1), and point forces in the x- and y-direction (SOURCE\_TYPE=2,3). The force sources excite both P- and S-waves. The explosive source is located at the same position as the diagonal elements of the stress tensor, i.e. at (i,j) (Figure \ref{fig_cell}). The forces are located at the same position as the corresponding components of particle velocity (Figure \ref{fig_cell}). If (x,y) denotes the position at which the source location is defined in source.dat, then the actual force in x-direction is located at (x+DX/2,y) and the actual force in y-direction is located at (x,y+DY/2). With SOURCE\_TYPE=4 a custom directive force can be defined by a force angle between y and x. The angle of the force must be specified in the SOURCE\_FILE after AMP. This force is not aligned along the main directions. \newline The locations of multiple sources must be defined in an external ASCII file (SOURCE\_FILE) that has the following format: {\color{blue}{\begin{verbatim} ... ... @@ -620,6 +627,7 @@ smaller than PRO the inversion aborts or in case of using frequency filtering (T \newpage \section{Source wavelet inversion} \label{sec:STF} To remove the contribution of the unknown source time function (STF) from the waveform residuals, it is necessary to design a filter which minimizes the misfit to the field recordings and raw synthetics. Therefore, a second forward simulation is applied. The first one is done with the wavelet specified in SOURCE\_SHAPE and the second one with the optimized source wavelet saved in SIGNAL\_FILE (see Section~\ref{sec:sources}). This optimized source wavelet is kept constant within N\_STF or within a frequency range (see below).\\ {\color{blue}{\begin{verbatim} ... ...
 *.dat *.txt *.sh *shot* \ No newline at end of file *.* *~ \ No newline at end of file
 ... ... @@ -17,7 +17,7 @@ -----------------------------------------------------------------------------------------*/ /* ---------------------------------------------------------------------- * This is program IFOS. * This is program IFOS Version 2.0.1 * Inversion of Full Observerd Seismograms * * ----------------------------------------------------------------------*/ ... ... @@ -95,9 +95,9 @@ int main(int argc, char **argv){ int step1, step2, step3=0, itests, iteste, stepmax, countstep; float scalefac; /* Variables for Pseudo-Hessian calculation */ int RECINC, ntr1; int SOURCE_SHAPE_OLD; int SOURCE_SHAPE_OLD=0; int SOURCE_SHAPE_OLD_SH=0; /* Variables for L-BFGS */ int LBFGS_NPAR=3; ... ... @@ -903,11 +903,6 @@ int main(int argc, char **argv){ MPI_Barrier(MPI_COMM_WORLD); /* comunication initialisation for persistent communication */ /*comm_ini(bufferlef_to_rig, bufferrig_to_lef, buffertop_to_bot, bufferbot_to_top, req_send, req_rec);*/ snapseis=1; snapseis1=1; SHOTINC=1; RECINC=1; ... ... @@ -917,7 +912,9 @@ int main(int argc, char **argv){ case 2: FC_EXT=filter_frequencies(&nfrq); FC=FC_EXT[FREQ_NR]; break; } /* Save old SOURCE_SHAPE, which is needed for STF */ SOURCE_SHAPE_OLD = SOURCE_SHAPE; if(WAVETYPE==2 || WAVETYPE==3) SOURCE_SHAPE_OLD_SH=SOURCE_SHAPE_SH; nt_out=10000; if(!VERBOSE) nt_out=1e5; ... ... @@ -1163,7 +1160,8 @@ int main(int argc, char **argv){ for (ishot=1;ishot<=nshots;ishot+=SHOTINC){ SOURCE_SHAPE = SOURCE_SHAPE_OLD; if(WAVETYPE==2 || WAVETYPE==3) SOURCE_SHAPE_SH=SOURCE_SHAPE_OLD_SH; /*------------------------------------------------------------------------------*/ /*----------- Start of inversion of source time function -----------------------*/ /*------------------------------------------------------------------------------*/ ... ... @@ -1189,7 +1187,7 @@ int main(int argc, char **argv){ srcpos_loc = splitsrc(srcpos,&nsrc_loc, nsrc); } if((SOURCE_SHAPE==7)||(SOURCE_SHAPE==3))declare_error("SOURCE_SHAPE==7 or SOURCE_SHAPE==3 isn't possible with INV_STF==1"); if(SOURCE_SHAPE==3) declare_error("SOURCE_SHAPE==3 isn't possible with INV_STF==1"); MPI_Barrier(MPI_COMM_WORLD); ... ... @@ -1199,12 +1197,12 @@ int main(int argc, char **argv){ /* calculate wavelet for each source point P SV */ if(WAVETYPE==1||WAVETYPE==3){ signals=NULL; signals=wavelet(srcpos_loc,nsrc_loc,ishot,0); signals=wavelet(srcpos_loc,nsrc_loc,ishot,0,1); } /* calculate wavelet for each source point SH */ if(WAVETYPE==2||WAVETYPE==3){ signals_SH=NULL; signals_SH=wavelet(srcpos_loc,nsrc_loc,ishot,1); signals_SH=wavelet(srcpos_loc,nsrc_loc,ishot,1,1); } ... ... @@ -1558,10 +1556,19 @@ int main(int argc, char **argv){ srcpos_loc = splitsrc(srcpos,&nsrc_loc, nsrc); } /*-------------------*/ /* Use STF wavelet */ /*-------------------*/ if(INV_STF){ SOURCE_SHAPE=7; if(WAVETYPE==1||WAVETYPE==3) fprintf(FP,"\n Using optimized source time function located in %s.shot%d \n",SIGNAL_FILE,ishot); if(WAVETYPE==2||WAVETYPE==3) fprintf(FP,"\n Using optimized source time function located in %s.shot%d \n",SIGNAL_FILE_SH,ishot); if( WAVETYPE==1 || WAVETYPE==3 ) fprintf(FP,"\n Using optimized source time function located in %s.shot%d \n",SIGNAL_FILE,ishot); if( WAVETYPE==2 || WAVETYPE==3 ) { SOURCE_SHAPE_SH=7; fprintf(FP,"\n Using optimized source time function located in %s.shot%d \n",SIGNAL_FILE_SH,ishot); } } MPI_Barrier(MPI_COMM_WORLD); ... ... @@ -1572,19 +1579,19 @@ int main(int argc, char **argv){ /* calculate wavelet for each source point P SV */ if(WAVETYPE==1||WAVETYPE==3){ signals=NULL; signals=wavelet(srcpos_loc,nsrc_loc,ishot,0); signals=wavelet(srcpos_loc,nsrc_loc,ishot,0,0); } /* calculate wavelet for each source point SH */ if(WAVETYPE==2||WAVETYPE==3){ signals_SH=NULL; signals_SH=wavelet(srcpos_loc,nsrc_loc,ishot,1); signals_SH=wavelet(srcpos_loc,nsrc_loc,ishot,1,0); } /*------------------------------------------------------------------------------*/ /*----------- Start of Time Domain Filtering -----------------------------------*/ /*------------------------------------------------------------------------------*/ if (((TIME_FILT==1) || (TIME_FILT==2))){ if (((TIME_FILT==1) || (TIME_FILT==2)) && (SOURCE_SHAPE!=6) && (INV_STF==0)){ fprintf(FP,"\n Time Domain Filter applied: Lowpass with corner frequency of %.2f Hz, order %d\n",FC,ORDER); /*time domain filtering of the source signal */ ... ... @@ -2043,6 +2050,9 @@ int main(int argc, char **argv){ } /* end ADJOINT_TYPE */ } /* --------------------------------- */ /* read seismic data from SU file vz */ /* --------------------------------- */ if(WAVETYPE==2 || WAVETYPE==3){ inseis(fprec,ishot,sectionread,ntr_glob,ns,10,iter); if ((TIME_FILT==1 )|| (TIME_FILT==2)){ ... ... @@ -3091,6 +3101,9 @@ int main(int argc, char **argv){ calc_mat_change_test(waveconv_up,waveconv_rho_up,waveconv_u_up,prhonp1,prho,ppinp1,ppi,punp1,pu,iter,1,FORWARD_ONLY,alpha_SL,0,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start); alpha_SL_old=1; /* If minimum number of iterations would be enforced, L-BFGS is likely to crash */ min_iter_help=0; } /*-----------------------------------------------------*/ ... ... @@ -3262,12 +3275,12 @@ int main(int argc, char **argv){ /* calculate wavelet for each source point P SV */ if(WAVETYPE==1||WAVETYPE==3){ signals=NULL; signals=wavelet(srcpos_loc,nsrc_loc,ishot,0); signals=wavelet(srcpos_loc,nsrc_loc,ishot,0,0); } /* calculate wavelet for each source point SH */ if(WAVETYPE==2||WAVETYPE==3){ signals_SH=NULL; signals_SH=wavelet(srcpos_loc,nsrc_loc,ishot,1); signals_SH=wavelet(srcpos_loc,nsrc_loc,ishot,1,0); } /*------------------------------------------------------------------------------*/ ... ... @@ -3275,7 +3288,7 @@ int main(int argc, char **argv){ /*------------------------------------------------------------------------------*/ /*time domain filtering of the source signal */ if(WAVETYPE==1||WAVETYPE==3){ if(((TIME_FILT==1) || (TIME_FILT==2))){ if(((TIME_FILT==1) || (TIME_FILT==2)) && (INV_STF==0)){ timedomain_filt(signals,FC,ORDER,nsrc_loc,ns,1); } ... ... @@ -3283,7 +3296,7 @@ int main(int argc, char **argv){ /*time domain filtering of the source signal */ if(WAVETYPE==2||WAVETYPE==3){ if(((TIME_FILT==1) || (TIME_FILT==2))){ if(((TIME_FILT==1) || (TIME_FILT==2)) && (INV_STF==0)){ timedomain_filt(signals_SH,FC,ORDER,nsrc_loc,ns,1); } ... ... @@ -3523,6 +3536,9 @@ int main(int argc, char **argv){ } } /* --------------------------------- */ /* read seismic data from SU file vz */ /* --------------------------------- */ if(WAVETYPE==2 || WAVETYPE==3){ inseis(fprec,ishot,sectionread,ntr_glob,ns,10,iter); if ((TIME_FILT==1 )|| (TIME_FILT==2)){ ... ...
 ... ... @@ -83,6 +83,7 @@ void exchange_par(void){ extern int VERBOSE; extern int TRKILL_STF_OFFSET; extern int TRKILL_STF_OFFSET_INVERT; extern float TRKILL_STF_OFFSET_LOWER; extern float TRKILL_STF_OFFSET_UPPER; ... ... @@ -374,6 +375,7 @@ void exchange_par(void){ idum[114]=TRKILL_OFFSET; idum[115]=TRKILL_STF_OFFSET; idum[116]=TRKILL_STF_OFFSET_INVERT; } /** if (MYID == 0) **/ ... ... @@ -660,6 +662,7 @@ void exchange_par(void){ TRKILL_OFFSET=idum[114]; TRKILL_STF_OFFSET=idum[115]; TRKILL_STF_OFFSET_INVERT=idum[116]; MPI_Bcast(&FL[1],L,MPI_FLOAT,0,MPI_COMM_WORLD); ... ...
 ... ... @@ -180,7 +180,7 @@ void outseis_vector(FILE *fp, FILE *fpdata, int comp, float *section, void inseis(FILE *fp, int comp, float **section, int ntr, int ns, int sws, int iter); void inseis_source_wavelet(float *section, int ns, int ishot, int SH); void inseis_source_wavelet(float *section, int ns, int ishot, int SH, int STF); void taper(float *section, int ns, float fc); ... ... @@ -280,7 +280,7 @@ void surface_elastic_PML(int ndepth, float ** vx, float ** vy, float ** sxx, flo void surface_PML(int ndepth, float ** vx, float ** vy, float ** sxx, float ** syy, float ** sxy, float ** syz, float ***p, float ***q, float ** ppi, float ** pu, float **prho, float **ptaup, float **ptaus, float *etajm, float *peta, float * hc, float * K_x, float * a_x, float * b_x, float ** psi_vxx, float ** ux, float ** uy, float ** uxy, float ** uyz,float ** sxz,float **uxz); void timedomain_filt(float ** data, float fc, int order, int ntr, int ns, int method); void timedomain_filt_vector(float * data, float fc, int order, int ntr, int ns, int method); void timedomain_filt_vector(float * data, float fc, int order, int ns, int method); void time_window(float **sectiondata, int iter, int ntr_glob, int **recpos_loc, int ntr, int ns, int ishot); void time_window_glob(float **sectiondata, int iter, int ntr_glob, int ns, int ishot); ... ... @@ -400,7 +400,7 @@ void update_v_rsg_4th(int nx1, int nx2, int ny1, int ny2, int nt, float ** psxy, float ** prho, float ** srcpos_loc, float ** signals, int nsrc, float ** absorb_coeff); float ** wavelet(float ** srcpos_loc, int nsrc, int ishot,int SH); float ** wavelet(float ** srcpos_loc, int nsrc, int ishot, int SH, int STF); float ** wavelet_stf(int nsrc, int ishot, float ** signals_stf); void writebufs(float ** sxx, float ** syy, ... ...
 ... ... @@ -105,6 +105,7 @@ int TRKILL; char TRKILL_FILE[STRING_SIZE]; int TRKILL_OFFSET; int TRKILL_STF_OFFSET_INVERT; float TRKILL_OFFSET_LOWER; float TRKILL_OFFSET_UPPER; ... ...
 ... ... @@ -26,8 +26,8 @@ void info(FILE *fp){ fprintf(fp," ***********************************************************\n"); fprintf(fp," This is program IFOS2D. Version 2.0 \n"); fprintf(fp," Parallel 2-D elastic Finite Difference FWT code \n"); fprintf(fp," This is program IFOS2D. Version 2.0.1 \n"); fprintf(fp," Parallel 2-D elastic Full Waveform Inversion code. \n"); fprintf(fp," \n"); fprintf(fp," ***********************************************************\n"); fprintf(fp,"\n"); ... ...
 ... ... @@ -22,7 +22,7 @@ #include "fd.h" #include "segy.h" void inseis_source_wavelet(float *section, int ns, int ishot, int SH){ void inseis_source_wavelet(float *section, int ns, int ishot, int SH, int STF){ /* declaration of extern variables */ extern int MYID; ... ... @@ -38,17 +38,25 @@ void inseis_source_wavelet(float *section, int ns, int ishot, int SH){ char data[STRING_SIZE]; FILE *fpdata; if(SH==0) { if(USE_WORKFLOW){ sprintf(data,"%s.stage%d.shot%d.su",SIGNAL_FILE,WORKFLOW_STAGE,ishot); if (STF==0){ /* reading inverted signals */ if(SH==0) { if(USE_WORKFLOW){ sprintf(data,"%s.stage%d.shot%d.su",SIGNAL_FILE,WORKFLOW_STAGE,ishot); } else { sprintf(data,"%s.shot%d.su",SIGNAL_FILE,ishot); } } else { sprintf(data,"%s.shot%d.su",SIGNAL_FILE,ishot); if(USE_WORKFLOW){ sprintf(data,"%s.stage%d.shot%d.su",SIGNAL_FILE_SH,WORKFLOW_STAGE,ishot); } else { sprintf(data,"%s.shot%d.su",SIGNAL_FILE_SH,ishot); } } } else { if(USE_WORKFLOW){ sprintf(data,"%s.stage%d.shot%d.su",SIGNAL_FILE_SH,WORKFLOW_STAGE,ishot); } else { /* reading signals for STF inversion */ if(SH==0) { sprintf(data,"%s.shot%d_start.su",SIGNAL_FILE,ishot); } else { sprintf(data,"%s.shot%d.su",SIGNAL_FILE_SH,ishot); sprintf(data,"%s.shot%d_start.su",SIGNAL_FILE_SH,ishot); } } ... ...