Commit 2cee0b7d authored by laura.gassner's avatar laura.gassner

Merge branch 'develop' into Seitosh/import

parents 546a9051 ec4496c1
......@@ -642,7 +642,7 @@ To remove the contribution of the unknown source time function (STF) from the wa
"TRKILL_STF" : "0",
"TRKILL_FILE_STF" : "./trace_kill/trace_kill",
"STF_FULL" : "0",
"TRKILL_STF_OFFSET" : "0",
"TRKILL_STF_OFFSET_LOWER" : "10",
"TRKILL_STF_OFFSET_UPPER" : "20",
......@@ -654,7 +654,7 @@ Default values are:
INV_STF=0
\end{verbatim}}}
INV\_STF should be switched to 1 if you want to invert for the source time function.
INV\_STF should be switched to 1 if you want to invert for the source time function. If STF\_FULL is set to 1 then the total wavefield is used to invert for the source time function and the time window is ignored.
\newline
An example for the parameter string provided in PARA is:
......@@ -718,7 +718,7 @@ With F\_HIGH\_PASS an additional high pass filter can be applied, where F\_HIGH\
With the parameter PRO (see~\ref{json:abort_criterion}) one has to adjust the criterion that defines at which points the bandwidth of the signals are increased.
With the parameter WRITE\_FILTERED\_DATA it is possible to write the time filtered measured data to disk which are filtered with the same filter as the synthetic data. Therefore this output can be used to visualize the residuals between synthetic and measured data. The filtered data is located in DATA\_DIR and are labeled with "\_measured".
With the parameter WRITE\_FILTERED\_DATA=1 it is possible to write the time filtered measured data to disk which are filtered with the same filter as the synthetic data. Therefore this output can be used to visualize the residuals between synthetic and measured data. The filtered data is located in SEIS\_FILE and the seismograms are labeled with "obs" (synthetic seismograms are labeled with "syn"). By choosing WRITE\_FILTERED\_DATA=2 one can directly output the seismograms that are used for calculating the adjoint sources. Additionally time windowing and integration (see Option VELOCITY) are considered and the label "adj" is added to the seismogram name.
If you are using frequeny filtering (TIME\_FILT==1) during the inversion, you can set a minimum number of iterations per frequency. Within this minimum number of iteration per frequency the abort criterion PRO will receive no consideration.
......@@ -906,4 +906,4 @@ Default values are:
MODEL_FILTER=0
\end{verbatim}}}
If MODEL\_FILTER = 1 vp- and vs-models are smoothed with a 2D median filter after every iterationstep. With FILT\_SIZE you can set the filter length in gridpoints.
\ No newline at end of file
If MODEL\_FILTER = 1 vp- and vs-models are smoothed with a 2D median filter after every iterationstep. With FILT\_SIZE you can set the filter length in gridpoints.
This diff is collapsed.
......@@ -29,7 +29,7 @@ double calc_misfit(float **sectiondata, float **section, int ntr, int ns, int LN
extern int TRKILL, NORMALIZE, F_LOW_PASS, TIMEWIN;
extern char TRKILL_FILE[STRING_SIZE];
extern int VELOCITY;
extern int WRITE_FILTERED_DATA;
int i,j;
float l2;
int h;
......@@ -202,6 +202,15 @@ double calc_misfit(float **sectiondata, float **section, int ntr, int ns, int LN
}
}
if(WRITE_FILTERED_DATA==2){
for(i=1;i<=ntr;i++){
for(j=1;j<=ns;j++){
sectiondata[i][j]=intseis_sectiondata[i][j];
section[i][j]=intseis_section[i][j];
}
}
}
l2=L2;
/* printf("\n MYID = %i IN CALC_MISFIT: L2 = %10.12f \n",MYID,l2); */
......
......@@ -30,7 +30,7 @@ void exchange_par(void){
extern float XREC1, XREC2, YREC1, YREC2, FPML;
extern float REC_ARRAY_DEPTH, REC_ARRAY_DIST, MUN, EPSILON, EPSILON_u, EPSILON_rho;
extern int SEISMO, NDT, NGEOPH, SEIS_FORMAT, FREE_SURF, READMOD, READREC, SRCREC;
extern int BOUNDARY, REC_ARRAY, DRX, FW;
extern int BOUNDARY, REC_ARRAY, DRX, FW, STF_FULL;
extern int SNAPSHOT_START,SNAPSHOT_END,SNAPSHOT_INCR;
extern float TSNAP1, TSNAP2, TSNAPINC, REFREC[4];
extern char MFILE[STRING_SIZE], SIGNAL_FILE[STRING_SIZE],SIGNAL_FILE_SH[STRING_SIZE], LOG_FILE[STRING_SIZE];
......@@ -377,7 +377,7 @@ void exchange_par(void){
idum[116]=TRKILL_STF_OFFSET_INVERT;
idum[117]=JOINT_EQUAL_WEIGHTING;
idum[118]=STF_FULL;
} /** if (MYID == 0) **/
MPI_Barrier(MPI_COMM_WORLD);
......@@ -664,7 +664,7 @@ void exchange_par(void){
TRKILL_STF_OFFSET_INVERT=idum[116];
JOINT_EQUAL_WEIGHTING=idum[117];
STF_FULL=idum[118];
if ( MYID!=0 && L>0 ) {
FL=vector(1,L);
}
......
......@@ -14,7 +14,7 @@ float XREC1, XREC2, YREC1, YREC2;
float REC_ARRAY_DEPTH, REC_ARRAY_DIST;
float REFREC[4]={0.0, 0.0, 0.0, 0.0}, FPML;
int SEISMO, NDT, NGEOPH, NSRC=1, SEIS_FORMAT, FREE_SURF, READMOD, READREC, SRCREC, FW=0;
int NX, NY, NT, SOURCE_SHAPE,SOURCE_SHAPE_SH, SOURCE_TYPE, SNAP, SNAP_FORMAT, REC_ARRAY, RUN_MULTIPLE_SHOTS, NTRG;
int NX, NY, NT, SOURCE_SHAPE,SOURCE_SHAPE_SH, SOURCE_TYPE, SNAP, SNAP_FORMAT, REC_ARRAY, RUN_MULTIPLE_SHOTS, NTRG,STF_FULL;
int L, BOUNDARY, DC, DRX, NXG, NYG, IDX, IDY, FDORDER, MAXRELERROR;
char SNAP_FILE[STRING_SIZE], SOURCE_FILE[STRING_SIZE], SIGNAL_FILE[STRING_SIZE], SIGNAL_FILE_SH[STRING_SIZE];
char MFILE[STRING_SIZE], REC_FILE[STRING_SIZE];
......@@ -147,4 +147,4 @@ int JOINT_EQUAL_WEIGHTING;
float JOINT_INVERSION_PSV_SH_ALPHA_VS;
float JOINT_INVERSION_PSV_SH_ALPHA_RHO;
int SNAPSHOT_START,SNAPSHOT_END,SNAPSHOT_INCR;
\ No newline at end of file
int SNAPSHOT_START,SNAPSHOT_END,SNAPSHOT_INCR;
......@@ -128,6 +128,7 @@ void read_par_json(FILE *fp, char *fileinp){
extern float WOLFE_C1_SL;
extern float WOLFE_C2_SL;
extern int STF_FULL;
/* definition of local variables */
int number_readobjects=0,fserr=0;
......@@ -763,6 +764,9 @@ void read_par_json(FILE *fp, char *fileinp){
TRKILL_STF=0;
fprintf(fp,"Variable TRKILL_STF is set to default value %d.\n",TRKILL_STF);}
else {
if (get_int_from_objectlist("STF_FULL",number_readobjects,&STF_FULL,varname_list, value_list)){
STF_FULL=0;
fprintf(fp,"Variable STF_FULL is set to default value %d.\n",STF_FULL);}
if (TRKILL_STF==1) {
if (get_int_from_objectlist("TRKILL_STF_OFFSET",number_readobjects,&TRKILL_STF_OFFSET,varname_list, value_list)){
TRKILL_STF_OFFSET=0;
......
......@@ -24,7 +24,7 @@
void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectionvz,float **sectionp,float **sectioncurl, float **sectiondiv, int **recpos, int **recpos_loc,int ntr, float ** srcpos, int ishot, int ns, int iter, int type_switch){
/* type_switch:
/* type_switch:
* 1== synthetic data
* 2== measured - synthetic data (residuals)
* 3== filtered measured data
......@@ -32,20 +32,21 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
extern int SEISMO, SEIS_FORMAT, RUN_MULTIPLE_SHOTS, WAVETYPE, VERBOSE,FORWARD_ONLY;
extern char SEIS_FILE[STRING_SIZE];
extern int VELOCITY, WRITE_FILTERED_DATA;
char vxf[STRING_SIZE], vyf[STRING_SIZE],vzf[STRING_SIZE], curlf[STRING_SIZE], divf[STRING_SIZE], pf[STRING_SIZE];
int nsrc=1;
switch (type_switch) {
case 1:
sprintf(vxf,"%s_vx.su.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(vyf,"%s_vy.su.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);
if(WAVETYPE==2 || WAVETYPE==3) {
sprintf(vzf,"%s_vz.su.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(vzf,"%s_vz.su.syn.shot%d.it%d",SEIS_FILE,ishot,iter);
}
sprintf(pf,"%s_p.su.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(divf,"%s_div.su.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(curlf,"%s_curl.su.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(pf,"%s_p.su.syn.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(divf,"%s_div.su.syn.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(curlf,"%s_curl.su.syn.shot%d.it%d",SEIS_FILE,ishot,iter);
break;
case 2:
......@@ -60,14 +61,36 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
break;
case 3:
sprintf(vxf,"%s_vx.su.measured.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(vyf,"%s_vy.su.measured.shot%d.it%d",SEIS_FILE,ishot,iter);
if(WRITE_FILTERED_DATA==1){
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);
if(WAVETYPE==2 || WAVETYPE==3) {
sprintf(vzf,"%s_vz.su.obs.shot%d.it%d",SEIS_FILE,ishot,iter);
}
sprintf(pf,"%s_p.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);
}else if(WRITE_FILTERED_DATA==2){
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);
if(WAVETYPE==2 || WAVETYPE==3) {
sprintf(vzf,"%s_vz.su.obs.adj.shot%d.it%d",SEIS_FILE,ishot,iter);
}
sprintf(pf,"%s_p.su.obs.adj.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(divf,"%s_div.su.obs.adj.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(curlf,"%s_curl.su.obs.adj.shot%d.it%d",SEIS_FILE,ishot,iter);
}
break;
case 4:
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);
if(WAVETYPE==2 || WAVETYPE==3) {
sprintf(vzf,"%s_vz.su.measured.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(vzf,"%s_vz.su.syn.adj.shot%d.it%d",SEIS_FILE,ishot,iter);
}
sprintf(pf,"%s_p.su.measured.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(divf,"%s_div.su.measured.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(curlf,"%s_curl.su.measured.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(pf,"%s_p.su.syn.adj.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(divf,"%s_div.su.syn.adj.shot%d.it%d",SEIS_FILE,ishot,iter);
sprintf(curlf,"%s_curl.su.syn.adj.shot%d.it%d",SEIS_FILE,ishot,iter);
break;
default:
......
......@@ -28,7 +28,7 @@ void stf(FILE *fp, float **sectionvy, float ** sectionvy_obs, float ** sectionvy
/* declaration of global variables */
extern float DT, DH;
extern int SEIS_FORMAT, MYID, NT, SOURCE_SHAPE, TIME_FILT, TIMEWIN, TAPER_STF, ORDER;
extern int SEIS_FORMAT, MYID, NT, SOURCE_SHAPE, TIME_FILT, TIMEWIN, TAPER_STF, ORDER, STF_FULL;
extern char PARA[STRING_SIZE], DATA_DIR[STRING_SIZE];
extern int TRKILL_STF, NORMALIZE, USE_WORKFLOW, WORKFLOW_STAGE;
extern char TRKILL_FILE_STF[STRING_SIZE];
......@@ -160,10 +160,10 @@ void stf(FILE *fp, float **sectionvy, float ** sectionvy_obs, float ** sectionvy
}
/* trace killing ends here */
if(TIMEWIN==1){
if((TIMEWIN==1)&&(STF_FULL==0)){
time_window_glob(sectionvy, iter, ntr_glob, ns, ishot);
time_window_glob(sectionvy_obs, iter, ntr_glob, ns, ishot);
}
}
/* NORMALIZE TRACES */
if(NORMALIZE==1){
......@@ -358,4 +358,4 @@ void stf(FILE *fp, float **sectionvy, float ** sectionvy_obs, float ** sectionvy
free_imatrix(kill_tmp,1,ntr_glob,1,nshots);
free_ivector(kill_vector,1,ntr_glob);
}
}
\ No newline at end of file
}
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