Commit fbe07179 authored by nikolaos.athanasopoulos's avatar nikolaos.athanasopoulos

Merge branch 'time_window' into 'develop'

Time window

See merge request !7
parents 7db85202 108ea90a
...@@ -452,7 +452,7 @@ This section covers some general inversion parameters. The maximum number of ite ...@@ -452,7 +452,7 @@ This section covers some general inversion parameters. The maximum number of ite
If models are read from binary files appropriate file extensions are required for the different models (see section \ref{gen_of_mod}). Depending on the data different components of the seismic sections can be back propagated. For two component data (x- and y-component) set ADJOINT\_TYPE=1, only the y-component (ADJOINT\_TYPE=2) and only the x-component (ADJOINT\_TYPE=3). For the inversion of pressure seismograms ADJOINT\_TYPE=4 has to be used. If models are read from binary files appropriate file extensions are required for the different models (see section \ref{gen_of_mod}). Depending on the data different components of the seismic sections can be back propagated. For two component data (x- and y-component) set ADJOINT\_TYPE=1, only the y-component (ADJOINT\_TYPE=2) and only the x-component (ADJOINT\_TYPE=3). For the inversion of pressure seismograms ADJOINT\_TYPE=4 has to be used.
During the inversion the misfit values are saved in a log file specified in MISFIT\_LOG\_FILE. The log file consists of eight or nine columns and each line corresponds to one iteration step. The used step length is written in the first column. In the second to fourth column the three test step lengths used for the step length estimation are saved. The corresponding misfit values for these test step lengthes and the test shots are written to column five to seven. Column eight corresponds to the total misfit for all shots and if you use frequency filtering then the ninth column corresponds to the corner frequency of the lowpass filter used in the inversion step. During the inversion the misfit values are saved in a log file specified in MISFIT\_LOG\_FILE. The log file consists of nine or ten columns and each line corresponds to one iteration step. The used step length is written in the first column. In the second to fourth column the three test step lengths used for the step length estimation are saved. The corresponding misfit values for these test step lengthes and the test shots are written to column five to seven. Column eight corresponds to the total misfit for all shots and if you use frequency filtering then the ninth column corresponds to the corner frequency of the lowpass filter used in the inversion step. If TIMEWIN=1 then the tenth column shows the GAMMA value that was used at the current iteration.
In general IFOS2D tries to minimize the misfit in the particle displacement between the observed data and the synthetic data. If you set the switch VELOCITY to 1 the misfit in the particle velocity seismograms is minimized. In general IFOS2D tries to minimize the misfit in the particle displacement between the observed data and the synthetic data. If you set the switch VELOCITY to 1 the misfit in the particle velocity seismograms is minimized.
...@@ -489,12 +489,12 @@ Default values are: ...@@ -489,12 +489,12 @@ Default values are:
With the use of a workflow file, you can define different FWI stages. For instance, one FWI stage could refer to one corner frequency of a low-pass filter or to a specific time window. Every line in the workflow file corresponds to one FWI stage. To use a workflow file the switch USE\_WORKFLOW have to be set to 1. The algorithm will automatically change to the next line of the workflow file if the abort criterium of the current line is reached or if no step length could be found which reduces the misfit. The structure of the variables inside the workflow file is as follow: \\ With the use of a workflow file, you can define different FWI stages. For instance, one FWI stage could refer to one corner frequency of a low-pass filter or to a specific time window. Every line in the workflow file corresponds to one FWI stage. To use a workflow file the switch USE\_WORKFLOW have to be set to 1. The algorithm will automatically change to the next line of the workflow file if the abort criterium of the current line is reached or if no step length could be found which reduces the misfit. The structure of the variables inside the workflow file is as follow: \\
\noindent \noindent
\begin{tabular}{lllllllllllll} \begin{tabular}{llllllllllIlll}
\# & INV\_VS & INV\_VP & INV\_RHO & PRO & TIME\_FILT & F\_HIGH\_PASS & F\_LOW\_PASS & WAVETYPE & 0 & 0 & EPRECOND & EPSILON\_WE\\ \# & INV\_VS & INV\_VP & INV\_RHO & PRO & TIME\_FILT & F\_HIGH\_PASS & F\_LOW\_PASS & WAVETYPE & 0 & 0 & EPRECOND & EPSILON\_WE & GAMMA \\
\end{tabular} \end{tabular}
\ \\ \ \\
The first column is the number of the line. With INV\_* etc. you can activate the inversion for VS, VP or RHO, respectively. The abort criterium in percent for this FWI stage will be the declared in the variable PRO. With TIME\_FILT you can activate the frequency filtering with the corner frequency F\_HIGH\_PASS of the high-pass and F\_LOW\_PASS of the low-pass filter. WAVETYPE can be used to switch between PSV and SH modeling, however it is recommended to use PSV modeling (WAVETYPE==1) due to the current development on SH modeling. The following to zeros are placeholders for an upcoming update. With EPRECOND and EPSILON\_WE you can control the approx. Hessian. Please note, that all features which are used eg. TIME\_FILT (see section \ref{sec:filtering}) within the workflow have to be activated in the .JSON file. The first column is the number of the line. With INV\_* etc. you can activate the inversion for VS, VP or RHO, respectively. The abort criterium in percent for this FWI stage will be the declared in the variable PRO. With TIME\_FILT you can activate the frequency filtering with the corner frequency F\_HIGH\_PASS of the high-pass and F\_LOW\_PASS of the low-pass filter. WAVETYPE can be used to switch between PSV and SH modeling, however it is recommended to use PSV modeling (WAVETYPE==1) due to the current development on SH modeling. The following two zeros are placeholders for an upcoming update. With EPRECOND and EPSILON\_WE you can control the approx. Hessian. Please note, that all features which are used eg. With GAMMA you can define the exponential taper applied to both synthetic and observed data. TIME\_FILT (see section \ref{sec:filtering}) within the workflow have to be activated in the .JSON file.
For an example of a workflow file, have a look in the par/ folder. For an example of a workflow file, have a look in the par/ folder.
...@@ -731,7 +731,7 @@ If you are using frequeny filtering (TIME\_FILT==1) during the inversion, you ca ...@@ -731,7 +731,7 @@ If you are using frequeny filtering (TIME\_FILT==1) during the inversion, you ca
"PICKS_FILE" : "./picked_times/picks" "PICKS_FILE" : "./picked_times/picks"
"TWLENGTH_PLUS" : "0.01", "TWLENGTH_PLUS" : "0.01",
"TWLENGTH_MINUS" : "0.01", "TWLENGTH_MINUS" : "0.01",
"GAMMA" : "100000", "GAMMA" : "30",
\end{verbatim}}} \end{verbatim}}}
{\color{red}{\begin{verbatim} {\color{red}{\begin{verbatim}
...@@ -741,7 +741,7 @@ Default values are: ...@@ -741,7 +741,7 @@ Default values are:
To apply time windowing in a time series the parameter TIMEWIN must set to 1. A automatic picker routine is not integrated. The point in time (picked time) for each source must be specified in separate files. The folder and file name can be set with the parameter PICKS\_FILE. The files must be named like this PICKS\_FILE\_<sourcenumber>.dat. So the number of sources must be equal to the number of files. Each file must contain the picked times for every receiver. To apply time windowing in a time series the parameter TIMEWIN must set to 1. A automatic picker routine is not integrated. The point in time (picked time) for each source must be specified in separate files. The folder and file name can be set with the parameter PICKS\_FILE. The files must be named like this PICKS\_FILE\_<sourcenumber>.dat. So the number of sources must be equal to the number of files. Each file must contain the picked times for every receiver.
The parameters TWLENGTH\_PLUS and TWLENGTH\_MINUS specify the length of the time window after (PLUS) and before (MINUS) the picked time. The unit is seconds. The damping factor GAMMA must be set individually. The parameters TWLENGTH\_PLUS and TWLENGTH\_MINUS specify the length of the time window after (PLUS) and before (MINUS) the picked time. The unit is seconds. The damping factor GAMMA must be set individually. In case where the workflow is used you can specify different values for GAMMA per stage.
When using TW\_IND = 1 three columns are expected in PICK\_FILE for each trace with the picked time (PT), time window length plus (TWLP) and time window length minus (TWLM). In this case TWLENGTH\_PLUS and TWLENGTH\_MINUS are ignored. It is also possible to use two time windows for each trace which can be utilized with TW\_IND = 2. PICK\_FILE must then contain six columns with values for PT\_1, TWLP\_1, TWLM\_1, PT\_2, TWLP\_2 and TWLM\_2 for each trace. When using TW\_IND = 1 three columns are expected in PICK\_FILE for each trace with the picked time (PT), time window length plus (TWLP) and time window length minus (TWLM). In this case TWLENGTH\_PLUS and TWLENGTH\_MINUS are ignored. It is also possible to use two time windows for each trace which can be utilized with TW\_IND = 2. PICK\_FILE must then contain six columns with values for PT\_1, TWLP\_1, TWLM\_1, PT\_2, TWLP\_2 and TWLM\_2 for each trace.
......
# VS VP RHO PRO F_FILT F_HIGH F_LOW WT J J EPRE EPSI # VS VP RHO PRO F_FILT F_HIGH F_LOW WT J J EPRE EPSI GAMMA
1 1 0 0 0.01 0 0 0 2 0 0 0 0.005 1 1 0 0 0.01 0 0 0 2 0 0 0 0.005 20
2 1 0 0 0.01 0 0 0 2 0 0 0 0.005 2 1 0 0 0.01 0 0 0 2 0 0 0 0.005 10
3 1 0 0 0.01 0 0 0 2 0 0 0 0.005 3 1 0 0 0.01 0 0 0 2 0 0 0 0.005 5
4 1 0 0 0.01 0 0 0 2 0 0 0 0.005 4 1 0 0 0.01 0 0 0 2 0 0 0 0.005 0
...@@ -957,8 +957,8 @@ int main(int argc, char **argv){ ...@@ -957,8 +957,8 @@ int main(int argc, char **argv){
if(PARAMETERIZATION==2){ if(PARAMETERIZATION==2){
Vp0[j][i] = sqrt((ppi[j][i]+2.0*pu[j][i])*prho[j][i]); Vp0[j][i] = sqrt((ppi[j][i]+2.0*pu[j][i])/prho[j][i]);
Vs0[j][i] = sqrt((pu[j][i])*prho[j][i]); Vs0[j][i] = sqrt((pu[j][i])/prho[j][i]);
Rho0[j][i] = prho[j][i]; Rho0[j][i] = prho[j][i];
} }
...@@ -1003,9 +1003,9 @@ int main(int argc, char **argv){ ...@@ -1003,9 +1003,9 @@ int main(int argc, char **argv){
/* Write header for misfit log file */ /* Write header for misfit log file */
if(GRAD_METHOD==1&&VERBOSE) { if(GRAD_METHOD==1&&VERBOSE) {
if (TIME_FILT==0){ if (TIME_FILT==0){
fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \n");} fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \t GAMMA \n");}
else{ else{
fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \t F_LOW_PASS \n"); fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \t F_LOW_PASS \t GAMMA \n");
} }
} }
...@@ -3108,9 +3108,9 @@ int main(int argc, char **argv){ ...@@ -3108,9 +3108,9 @@ int main(int argc, char **argv){
if(wolfe_SLS_failed) { if(wolfe_SLS_failed) {
if (TIME_FILT==0){ if (TIME_FILT==0){
if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \n",0.0,iter,wolfe_sum_FWI,0.0,countstep-1,L2_SL_old,L2_SL_old);} if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \t %f\n",0.0,iter,wolfe_sum_FWI,0.0,countstep-1,L2_SL_old,L2_SL_old, GAMMA);}
else{ else{
if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \t %f\n",0.0,iter,wolfe_sum_FWI,0.0,countstep-1,L2_SL_old,L2_SL_old,F_LOW_PASS); if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \t %f \t %f \n",0.0,iter,wolfe_sum_FWI,0.0,countstep-1,L2_SL_old,L2_SL_old,F_LOW_PASS, GAMMA);
} }
if(WAVETYPE==3 && MYID==0){ if(WAVETYPE==3 && MYID==0){
...@@ -3153,9 +3153,9 @@ int main(int argc, char **argv){ ...@@ -3153,9 +3153,9 @@ int main(int argc, char **argv){
float diff=0.0; float diff=0.0;
diff=fabs((L2_hist[iter-2]-L2_hist[iter])/L2_hist[iter-2]); diff=fabs((L2_hist[iter-2]-L2_hist[iter])/L2_hist[iter-2]);
if (TIME_FILT==0){ if (TIME_FILT==0){
if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \n",alpha_SL,iter,wolfe_sum_FWI,diff,countstep-1,L2_SL_old,L2_SL_new);} if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \t %f \n",alpha_SL,iter,wolfe_sum_FWI,diff,countstep-1,L2_SL_old,L2_SL_new, GAMMA);}
else{ else{
if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \t %f\n",alpha_SL,iter,wolfe_sum_FWI,diff,countstep-1,L2_SL_old,L2_SL_new,F_LOW_PASS); if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \t %f \t %f\n",alpha_SL,iter,wolfe_sum_FWI,diff,countstep-1,L2_SL_old,L2_SL_new,F_LOW_PASS, GAMMA);
} }
if(WAVETYPE==3 && MYID==0){ if(WAVETYPE==3 && MYID==0){
...@@ -3825,9 +3825,9 @@ int main(int argc, char **argv){ ...@@ -3825,9 +3825,9 @@ int main(int argc, char **argv){
if(MYID==0&&(!VERBOSE)){printf("L2-Norm[4] = %e (final L2 of all shots)\n",L2t[4]);} if(MYID==0&&(!VERBOSE)){printf("L2-Norm[4] = %e (final L2 of all shots)\n",L2t[4]);}
printf("MYID = %d \t epst1[1] = %e \t epst1[2] = %e \t epst1[3] = %e \n",MYID,epst1[1],epst1[2],epst1[3]); printf("MYID = %d \t epst1[1] = %e \t epst1[2] = %e \t epst1[3] = %e \n",MYID,epst1[1],epst1[2],epst1[3]);
if (TIME_FILT==0){ if (TIME_FILT==0){
fprintf(FPL2,"%e \t %e \t %e \t %e \t %e \t %e \t %e \t %e \n",opteps_vp,epst1[1],epst1[2],epst1[3],L2t[1],L2t[2],L2t[3],L2t[4]);} fprintf(FPL2,"%e \t %e \t %e \t %e \t %e \t %e \t %e \t %e \t %f \n",opteps_vp,epst1[1],epst1[2],epst1[3],L2t[1],L2t[2],L2t[3],L2t[4],GAMMA);}
else{ else{
fprintf(FPL2,"%e \t %e \t %e \t %e \t %e \t %e \t %e \t %e \t %f\n",opteps_vp,epst1[1],epst1[2],epst1[3],L2t[1],L2t[2],L2t[3],L2t[4],F_LOW_PASS);} fprintf(FPL2,"%e \t %e \t %e \t %e \t %e \t %e \t %e \t %e \t %f \t %f \n",opteps_vp,epst1[1],epst1[2],epst1[3],L2t[1],L2t[2],L2t[3],L2t[4],F_LOW_PASS,GAMMA);}
if(WAVETYPE==3 && MYID==0){ if(WAVETYPE==3 && MYID==0){
fprintf(FPL2_JOINT,"%d \t %f \t %f\n",iter,L2sum_all_shots/energy_sum_all_shots,L2sum_all_shots_SH/energy_sum_all_shots_SH); fprintf(FPL2_JOINT,"%d \t %f \t %f\n",iter,L2sum_all_shots/energy_sum_all_shots,L2sum_all_shots_SH/energy_sum_all_shots_SH);
} }
......
...@@ -38,6 +38,7 @@ void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[ST ...@@ -38,6 +38,7 @@ void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[ST
extern float JOINT_INVERSION_PSV_SH_ALPHA_RHO; extern float JOINT_INVERSION_PSV_SH_ALPHA_RHO;
extern int EPRECOND; extern int EPRECOND;
extern float EPSILON_WE; extern float EPSILON_WE;
extern float GAMMA;
extern int GRAD_METHOD; extern int GRAD_METHOD;
extern int WORKFLOW_STAGE; extern int WORKFLOW_STAGE;
...@@ -126,8 +127,10 @@ void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[ST ...@@ -126,8 +127,10 @@ void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[ST
if(EPRECOND==0 && workflow[WORKFLOW_STAGE][12]!=0){ if(EPRECOND==0 && workflow[WORKFLOW_STAGE][12]!=0){
if(MYID==0) printf(" WARNING: EPRECOND have to be set >0 in JSON (if so, ignore this message)"); if(MYID==0) printf(" WARNING: EPRECOND have to be set >0 in JSON (if so, ignore this message)");
} }
EPRECOND=workflow[WORKFLOW_STAGE][12]; EPRECOND=workflow[WORKFLOW_STAGE][12];
EPSILON_WE=workflow[WORKFLOW_STAGE][13]; EPSILON_WE=workflow[WORKFLOW_STAGE][13];
GAMMA=workflow[WORKFLOW_STAGE][14];
if(*LBFGS_iter_start==*iter && GRAD_METHOD==2){ if(*LBFGS_iter_start==*iter && GRAD_METHOD==2){
if(MYID==0)printf("\n L-BFGS will be used from iteration %d on.",*LBFGS_iter_start+1); if(MYID==0)printf("\n L-BFGS will be used from iteration %d on.",*LBFGS_iter_start+1);
......
...@@ -32,15 +32,24 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio ...@@ -32,15 +32,24 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
extern int SEISMO, SEIS_FORMAT, RUN_MULTIPLE_SHOTS, WAVETYPE, VERBOSE,FORWARD_ONLY; extern int SEISMO, SEIS_FORMAT, RUN_MULTIPLE_SHOTS, WAVETYPE, VERBOSE,FORWARD_ONLY;
extern char SEIS_FILE[STRING_SIZE]; extern char SEIS_FILE[STRING_SIZE];
extern int VELOCITY, WRITE_FILTERED_DATA; extern int VELOCITY, WRITE_FILTERED_DATA, ADJOINT_TYPE;;
char vxf[STRING_SIZE], vyf[STRING_SIZE],vzf[STRING_SIZE], curlf[STRING_SIZE], divf[STRING_SIZE], pf[STRING_SIZE]; char vxf[STRING_SIZE], vyf[STRING_SIZE],vzf[STRING_SIZE], curlf[STRING_SIZE], divf[STRING_SIZE], pf[STRING_SIZE];
int nsrc=1; int nsrc=1;
switch (type_switch) { switch (type_switch) {
case 1: case 1:
if(ADJOINT_TYPE==1) {
sprintf(vxf,"%s_vx.su.syn.shot%d.it%d",SEIS_FILE,ishot,iter); sprintf(vxf,"%s_vx.su.syn.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(vyf,"%s_vy.su.syn.shot%d.it%d",SEIS_FILE,ishot,iter); sprintf(vyf,"%s_vy.su.syn.shot%d.it%d",SEIS_FILE,ishot,iter);
}
if(ADJOINT_TYPE==2){
sprintf(vyf,"%s_vy.su.syn.shot%d.it%d",SEIS_FILE,ishot,iter);
}
if(ADJOINT_TYPE==3){
sprintf(vxf,"%s_vx.su.syn.shot%d.it%d",SEIS_FILE,ishot,iter);
}
if(WAVETYPE==2 || WAVETYPE==3) { if(WAVETYPE==2 || WAVETYPE==3) {
sprintf(vzf,"%s_vz.su.syn.shot%d.it%d",SEIS_FILE,ishot,iter); sprintf(vzf,"%s_vz.su.syn.shot%d.it%d",SEIS_FILE,ishot,iter);
} }
...@@ -50,8 +59,16 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio ...@@ -50,8 +59,16 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
break; break;
case 2: case 2:
if(ADJOINT_TYPE==1) {
sprintf(vxf,"%s_vx.su.shot%d_adjoint_src",SEIS_FILE,ishot); sprintf(vxf,"%s_vx.su.shot%d_adjoint_src",SEIS_FILE,ishot);
sprintf(vyf,"%s_vy.su.shot%d_adjoint_src",SEIS_FILE,ishot); sprintf(vyf,"%s_vy.su.shot%d_adjoint_src",SEIS_FILE,ishot);
}
if(ADJOINT_TYPE==2){
sprintf(vyf,"%s_vy.su.shot%d_adjoint_src",SEIS_FILE,ishot);
}
if(ADJOINT_TYPE==3){
sprintf(vxf,"%s_vx.su.shot%d_adjoint_src",SEIS_FILE,ishot);
}
if(WAVETYPE==2 || WAVETYPE==3) { if(WAVETYPE==2 || WAVETYPE==3) {
sprintf(vzf,"%s_vz.su.shot%d_adjoint_src",SEIS_FILE,ishot); sprintf(vzf,"%s_vz.su.shot%d_adjoint_src",SEIS_FILE,ishot);
} }
...@@ -62,8 +79,17 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio ...@@ -62,8 +79,17 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
case 3: case 3:
if(WRITE_FILTERED_DATA==1){ if(WRITE_FILTERED_DATA==1){
if(ADJOINT_TYPE==1) {
sprintf(vxf,"%s_vx.su.obs.shot%d.it%d",SEIS_FILE,ishot,iter); sprintf(vxf,"%s_vx.su.obs.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(vyf,"%s_vy.su.obs.shot%d.it%d",SEIS_FILE,ishot,iter); sprintf(vyf,"%s_vy.su.obs.shot%d.it%d",SEIS_FILE,ishot,iter);
}
if(ADJOINT_TYPE==2){
sprintf(vyf,"%s_vy.su.obs.shot%d.it%d",SEIS_FILE,ishot,iter);
}
if(ADJOINT_TYPE==3){
sprintf(vxf,"%s_vx.su.obs.shot%d.it%d",SEIS_FILE,ishot,iter);
}
if(WAVETYPE==2 || WAVETYPE==3) { if(WAVETYPE==2 || WAVETYPE==3) {
sprintf(vzf,"%s_vz.su.obs.shot%d.it%d",SEIS_FILE,ishot,iter); sprintf(vzf,"%s_vz.su.obs.shot%d.it%d",SEIS_FILE,ishot,iter);
} }
...@@ -71,8 +97,17 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio ...@@ -71,8 +97,17 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
sprintf(divf,"%s_div.su.obs.shot%d.it%d",SEIS_FILE,ishot,iter); sprintf(divf,"%s_div.su.obs.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(curlf,"%s_curl.su.obs.shot%d.it%d",SEIS_FILE,ishot,iter); sprintf(curlf,"%s_curl.su.obs.shot%d.it%d",SEIS_FILE,ishot,iter);
}else if(WRITE_FILTERED_DATA==2){ }else if(WRITE_FILTERED_DATA==2){
if(ADJOINT_TYPE==1) {
sprintf(vxf,"%s_vx.su.obs.adj.shot%d.it%d",SEIS_FILE,ishot,iter); sprintf(vxf,"%s_vx.su.obs.adj.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(vyf,"%s_vy.su.obs.adj.shot%d.it%d",SEIS_FILE,ishot,iter); sprintf(vyf,"%s_vy.su.obs.adj.shot%d.it%d",SEIS_FILE,ishot,iter);
}
if(ADJOINT_TYPE==2){
sprintf(vyf,"%s_vy.su.obs.adj.shot%d.it%d",SEIS_FILE,ishot,iter);
}
if(ADJOINT_TYPE==3){
sprintf(vxf,"%s_vx.su.obs.adj.shot%d.it%d",SEIS_FILE,ishot,iter);
}
if(WAVETYPE==2 || WAVETYPE==3) { if(WAVETYPE==2 || WAVETYPE==3) {
sprintf(vzf,"%s_vz.su.obs.adj.shot%d.it%d",SEIS_FILE,ishot,iter); sprintf(vzf,"%s_vz.su.obs.adj.shot%d.it%d",SEIS_FILE,ishot,iter);
} }
...@@ -83,8 +118,17 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio ...@@ -83,8 +118,17 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
break; break;
case 4: case 4:
if(ADJOINT_TYPE==1) {
sprintf(vxf,"%s_vx.su.syn.adj.shot%d.it%d",SEIS_FILE,ishot,iter); sprintf(vxf,"%s_vx.su.syn.adj.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(vyf,"%s_vy.su.syn.adj.shot%d.it%d",SEIS_FILE,ishot,iter); sprintf(vyf,"%s_vy.su.syn.adj.shot%d.it%d",SEIS_FILE,ishot,iter);
}
if(ADJOINT_TYPE==2){
sprintf(vyf,"%s_vy.su.syn.adj.shot%d.it%d",SEIS_FILE,ishot,iter);
}
if(ADJOINT_TYPE==3){
sprintf(vxf,"%s_vx.su.syn.adj.shot%d.it%d",SEIS_FILE,ishot,iter);
}
if(WAVETYPE==2 || WAVETYPE==3) { if(WAVETYPE==2 || WAVETYPE==3) {
sprintf(vzf,"%s_vz.su.syn.adj.shot%d.it%d",SEIS_FILE,ishot,iter); sprintf(vzf,"%s_vz.su.syn.adj.shot%d.it%d",SEIS_FILE,ishot,iter);
} }
...@@ -99,8 +143,17 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio ...@@ -99,8 +143,17 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
} }
if(FORWARD_ONLY==1){ if(FORWARD_ONLY==1){
if(ADJOINT_TYPE==1) {
sprintf(vxf,"%s_vx.su.shot%d",SEIS_FILE,ishot); sprintf(vxf,"%s_vx.su.shot%d",SEIS_FILE,ishot);
sprintf(vyf,"%s_vy.su.shot%d",SEIS_FILE,ishot); sprintf(vyf,"%s_vy.su.shot%d",SEIS_FILE,ishot);
}
if(ADJOINT_TYPE==2){
sprintf(vyf,"%s_vy.su.shot%d",SEIS_FILE,ishot);
}
if(ADJOINT_TYPE==3){
sprintf(vxf,"%s_vx.su.shot%d",SEIS_FILE,ishot);
}
if(WAVETYPE==2 || WAVETYPE==3) { if(WAVETYPE==2 || WAVETYPE==3) {
sprintf(vzf,"%s_vz.su.shot%d",SEIS_FILE,ishot); sprintf(vzf,"%s_vz.su.shot%d",SEIS_FILE,ishot);
} }
...@@ -112,11 +165,22 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio ...@@ -112,11 +165,22 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
switch (SEISMO){ switch (SEISMO){
case 1 : /* particle velocities only */ case 1 : /* particle velocities only */
if (WAVETYPE==1 || WAVETYPE==3) { if (WAVETYPE==1 || WAVETYPE==3) {
if(ADJOINT_TYPE==1) {
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vx) to\n\t %s \n",0,ntr,vxf); if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vx) to\n\t %s \n",0,ntr,vxf);
outseis_glob(fp,fopen(vxf,"w"),1,sectionvx,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1); outseis_glob(fp,fopen(vxf,"w"),1,sectionvx,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vy) to\n\t %s \n",0,ntr,vyf); if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vy) to\n\t %s \n",0,ntr,vyf);
outseis_glob(fp,fopen(vyf,"w"),2,sectionvy,recpos,recpos_loc, ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1); outseis_glob(fp,fopen(vyf,"w"),2,sectionvy,recpos,recpos_loc, ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
} }
if(ADJOINT_TYPE==2){
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vy) to\n\t %s \n",0,ntr,vyf);
outseis_glob(fp,fopen(vyf,"w"),2,sectionvy,recpos,recpos_loc, ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
}
if(ADJOINT_TYPE==3){
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vx) to\n\t %s \n",0,ntr,vxf);
outseis_glob(fp,fopen(vxf,"w"),1,sectionvx,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
}
}
if(WAVETYPE==2 || WAVETYPE==3) { if(WAVETYPE==2 || WAVETYPE==3) {
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vz) to\n\t %s \n",0,ntr,vzf); if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vz) to\n\t %s \n",0,ntr,vzf);
outseis_glob(fp,fopen(vzf,"w"),2,sectionvz,recpos,recpos_loc, ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1); outseis_glob(fp,fopen(vzf,"w"),2,sectionvz,recpos,recpos_loc, ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
...@@ -137,10 +201,21 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio ...@@ -137,10 +201,21 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
case 4 : /* everything */ case 4 : /* everything */
if (WAVETYPE==1 || WAVETYPE==3) { if (WAVETYPE==1 || WAVETYPE==3) {
if(ADJOINT_TYPE==1) {
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vx) to\n\t %s \n",0,ntr,vxf); if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vx) to\n\t %s \n",0,ntr,vxf);
outseis_glob(fp,fopen(vxf,"w"),1,sectionvx,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1); outseis_glob(fp,fopen(vxf,"w"),1,sectionvx,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vy) to\n\t %s \n",0,ntr,vyf); if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vy) to\n\t %s \n",0,ntr,vyf);
outseis_glob(fp,fopen(vyf,"w"),2,sectionvy,recpos,recpos_loc, ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1); outseis_glob(fp,fopen(vyf,"w"),2,sectionvy,recpos,recpos_loc, ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
}
if(ADJOINT_TYPE==2){
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vy) to\n\t %s \n",0,ntr,vyf);
outseis_glob(fp,fopen(vyf,"w"),2,sectionvy,recpos,recpos_loc, ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
}
if(ADJOINT_TYPE==3){
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vx) to\n\t %s \n",0,ntr,vxf);
outseis_glob(fp,fopen(vxf,"w"),1,sectionvx,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
}
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms of pressure to\n\t %s \n",0,ntr,pf); if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms of pressure to\n\t %s \n",0,ntr,pf);
outseis_glob(fp,fopen(pf,"w"), 0, sectionp,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1); outseis_glob(fp,fopen(pf,"w"), 0, sectionp,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
...@@ -157,10 +232,21 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio ...@@ -157,10 +232,21 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
case 5 : /* everything except curl and div */ case 5 : /* everything except curl and div */
if (WAVETYPE==1 || WAVETYPE==3) { if (WAVETYPE==1 || WAVETYPE==3) {
if(ADJOINT_TYPE==1) {
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vx) to\n\t %s \n",0,ntr,vxf); if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vx) to\n\t %s \n",0,ntr,vxf);
outseis_glob(fp,fopen(vxf,"w"),1,sectionvx,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1); outseis_glob(fp,fopen(vxf,"w"),1,sectionvx,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vy) to\n\t %s \n",0,ntr,vyf); if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vy) to\n\t %s \n",0,ntr,vyf);
outseis_glob(fp,fopen(vyf,"w"),2,sectionvy,recpos,recpos_loc, ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1); outseis_glob(fp,fopen(vyf,"w"),2,sectionvy,recpos,recpos_loc, ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
}
if(ADJOINT_TYPE==2){
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vy) to\n\t %s \n",0,ntr,vyf);
outseis_glob(fp,fopen(vyf,"w"),2,sectionvy,recpos,recpos_loc, ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
}
if(ADJOINT_TYPE==3){
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vx) to\n\t %s \n",0,ntr,vxf);
outseis_glob(fp,fopen(vxf,"w"),1,sectionvx,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
}
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms of pressure to\n\t %s \n",0,ntr,pf); if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms of pressure to\n\t %s \n",0,ntr,pf);
outseis_glob(fp,fopen(pf,"w"), 0, sectionp,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1); outseis_glob(fp,fopen(pf,"w"), 0, sectionp,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
} }
......
...@@ -139,10 +139,10 @@ void time_window(float **sectiondata, int iter, int ntr_glob, int **recpos_loc, ...@@ -139,10 +139,10 @@ void time_window(float **sectiondata, int iter, int ntr_glob, int **recpos_loc,
time = (float)(j * DT); time = (float)(j * DT);
dumpa = (time-picked_times_m[1][i]-picked_times_m[2][i]); dumpa = (time-picked_times_m[1][i]-picked_times_m[2][i]);
taper = exp(-GAMMA*dumpa*dumpa); taper = exp(-GAMMA*dumpa);
dumpb = (time-picked_times_m[1][i]+picked_times_m[3][i]); dumpb = (time-picked_times_m[1][i]+picked_times_m[3][i]);
taper1 = exp(-GAMMA*dumpb*dumpb); taper1 = exp(-GAMMA*dumpb);
if(time>=picked_times_m[1][i]+picked_times_m[2][i]){ if(time>=picked_times_m[1][i]+picked_times_m[2][i]){
sectiondata[i][j] = sectiondata[i][j] * taper;} sectiondata[i][j] = sectiondata[i][j] * taper;}
...@@ -160,10 +160,10 @@ void time_window(float **sectiondata, int iter, int ntr_glob, int **recpos_loc, ...@@ -160,10 +160,10 @@ void time_window(float **sectiondata, int iter, int ntr_glob, int **recpos_loc,
time = (float)(j * DT); time = (float)(j * DT);
dumpa = (time-picked_times_m[1][i]-picked_times_m[2][i]); dumpa = (time-picked_times_m[1][i]-picked_times_m[2][i]);
taper = exp(-GAMMA*dumpa*dumpa); taper = exp(-GAMMA*dumpa);
dumpb = (time-picked_times_m[1][i]+picked_times_m[3][i]); dumpb = (time-picked_times_m[1][i]+picked_times_m[3][i]);
taper1 = exp(-GAMMA*dumpb*dumpb); taper1 = exp(-GAMMA*dumpb);
dummysection[i][j] = sectiondata[i][j]; dummysection[i][j] = sectiondata[i][j];
...@@ -174,10 +174,10 @@ void time_window(float **sectiondata, int iter, int ntr_glob, int **recpos_loc, ...@@ -174,10 +174,10 @@ void time_window(float **sectiondata, int iter, int ntr_glob, int **recpos_loc,
dummysection[i][j] = sectiondata[i][j] * taper1;} dummysection[i][j] = sectiondata[i][j] * taper1;}
dumpc = (time-picked_times_m[4][i]-picked_times_m[5][i]); dumpc = (time-picked_times_m[4][i]-picked_times_m[5][i]);
taper2 = exp(-GAMMA*dumpc*dumpc); taper2 = exp(-GAMMA*dumpc);
dumpd = (time-picked_times_m[4][i]+picked_times_m[6][i]); dumpd = (time-picked_times_m[4][i]+picked_times_m[6][i]);
taper3 = exp(-GAMMA*dumpd*dumpd); taper3 = exp(-GAMMA*dumpd);
if(time>=picked_times_m[4][i]+picked_times_m[5][i]){ if(time>=picked_times_m[4][i]+picked_times_m[5][i]){
sectiondata[i][j] = sectiondata[i][j] * taper2;} sectiondata[i][j] = sectiondata[i][j] * taper2;}
...@@ -196,10 +196,10 @@ void time_window(float **sectiondata, int iter, int ntr_glob, int **recpos_loc, ...@@ -196,10 +196,10 @@ void time_window(float **sectiondata, int iter, int ntr_glob, int **recpos_loc,
time = (float)(j * DT); time = (float)(j * DT);
dumpa = (time-picked_times[i]-TWLENGTH_PLUS); dumpa = (time-picked_times[i]-TWLENGTH_PLUS);
taper = exp(-GAMMA*dumpa*dumpa); taper = exp(-GAMMA*dumpa);
dumpb = (time-picked_times[i]+TWLENGTH_MINUS); dumpb = (time-picked_times[i]+TWLENGTH_MINUS);
taper1 = exp(-GAMMA*dumpb*dumpb); taper1 = exp(-GAMMA*dumpb);
if(time>=picked_times[i]+TWLENGTH_PLUS){ if(time>=picked_times[i]+TWLENGTH_PLUS){
sectiondata[i][j] = sectiondata[i][j] * taper;} sectiondata[i][j] = sectiondata[i][j] * taper;}
......
...@@ -101,10 +101,10 @@ void time_window_glob(float **sectiondata, int iter, int ntr_glob, int ns, int i ...@@ -101,10 +101,10 @@ void time_window_glob(float **sectiondata, int iter, int ntr_glob, int ns, int i
time = (float)(j * DT); time = (float)(j * DT);
dumpa = (time-picked_times_m[1][i]-picked_times_m[2][i]); dumpa = (time-picked_times_m[1][i]-picked_times_m[2][i]);
taper = exp(-GAMMA*dumpa*dumpa); taper = exp(-GAMMA*dumpa);
dumpb = (time-picked_times_m[1][i]+picked_times_m[3][i]); dumpb = (time-picked_times_m[1][i]+picked_times_m[3][i]);
taper1 = exp(-GAMMA*dumpb*dumpb); taper1 = exp(-GAMMA*dumpb);
if(time>=picked_times_m[1][i]+picked_times_m[2][i]){ if(time>=picked_times_m[1][i]+picked_times_m[2][i]){
sectiondata[i][j] = sectiondata[i][j] * taper;} sectiondata[i][j] = sectiondata[i][j] * taper;}
...@@ -123,10 +123,10 @@ void time_window_glob(float **sectiondata, int iter, int ntr_glob, int ns, int i ...@@ -123,10 +123,10 @@ void time_window_glob(float **sectiondata, int iter, int ntr_glob, int ns, int i
time = (float)(j * DT); time = (float)(j * DT);
dumpa = (time-picked_times_m[1][i]-picked_times_m[2][i]); dumpa = (time-picked_times_m[1][i]-picked_times_m[2][i]);
taper = exp(-GAMMA*dumpa*dumpa); taper = exp(-GAMMA*dumpa);
dumpb = (time-picked_times_m[1][i]+picked_times_m[3][i]); dumpb = (time-picked_times_m[1][i]+picked_times_m[3][i]);
taper1 = exp(-GAMMA*dumpb*dumpb); taper1 = exp(-GAMMA*dumpb);
dummysection[i][j] = sectiondata[i][j]; dummysection[i][j] = sectiondata[i][j];
...@@ -137,10 +137,10 @@ void time_window_glob(float **sectiondata, int iter, int ntr_glob, int ns, int i ...@@ -137,10 +137,10 @@ void time_window_glob(float **sectiondata, int iter, int ntr_glob, int ns, int i
dummysection[i][j] = sectiondata[i][j] * taper1;} dummysection[i][j] = sectiondata[i][j] * taper1;}
dumpc = (time-picked_times_m[4][i]-picked_times_m[5][i]); dumpc = (time-picked_times_m[4][i]-picked_times_m[5][i]);
taper2 = exp(-GAMMA*dumpc*dumpc); taper2 = exp(-GAMMA*dumpc);
dumpd = (time-picked_times_m[4][i]+picked_times_m[6][i]); dumpd = (time-picked_times_m[4][i]+picked_times_m[6][i]);
taper3 = exp(-GAMMA*dumpd*dumpd); taper3 = exp(-GAMMA*dumpd);
if(time>=picked_times_m[4][i]+picked_times_m[5][i]){ if(time>=picked_times_m[4][i]+picked_times_m[5][i]){
sectiondata[i][j] = sectiondata[i][j] * taper2;} sectiondata[i][j] = sectiondata[i][j] * taper2;}
...@@ -159,10 +159,10 @@ void time_window_glob(float **sectiondata, int iter, int ntr_glob, int ns, int i ...@@ -159,10 +159,10 @@ void time_window_glob(float **sectiondata, int iter, int ntr_glob, int ns, int i
time = (float)(j * DT); time = (float)(j * DT);
dumpa = (time-picked_times[i]-TWLENGTH_PLUS); dumpa = (time-picked_times[i]-TWLENGTH_PLUS);
taper = exp(-GAMMA*dumpa*dumpa); taper = exp(-GAMMA*dumpa);
dumpb = (time-picked_times[i]+TWLENGTH_MINUS); dumpb = (time-picked_times[i]+TWLENGTH_MINUS);
taper1 = exp(-GAMMA*dumpb*dumpb); taper1 = exp(-GAMMA*dumpb);
if(time>=picked_times[i]+TWLENGTH_PLUS){ if(time>=picked_times[i]+TWLENGTH_PLUS){
sectiondata[i][j] = sectiondata[i][j] * taper;} sectiondata[i][j] = sectiondata[i][j] * taper;}
......
Markdown is supported
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