Commit 48dc252c authored by nikolaos.athanasopoulos's avatar nikolaos.athanasopoulos
Browse files

Merge branch 'time_window' into 'master'

Time window

See merge request !6
parents 9accb585 108ea90a
......@@ -4,8 +4,9 @@
Procedures in the Fourier domain
--------------------------------
Options and parameters in common for procedures in the Fourier domain:
fpow2 use power of two for number of coefficients
fdiv=d use integer multiple of d for number of coefficients
fpad=f padding factor
tshift=d delay source correction filter wavelet by d (in seconds)
fpow2 use power of two for number of coefficients
fdiv=d use integer multiple of d for number of coefficients
fpad=f padding factor
tshift=d delay source correction filter wavelet by d (in seconds)
irtap=t1,t2,t3,t4 taper impulse response of correction filter
# ----- END OF stfinvfourier_summary_usage.txt -----
......@@ -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.
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.
......@@ -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: \\
\noindent
\begin{tabular}{lllllllllllll}
\# & INV\_VS & INV\_VP & INV\_RHO & PRO & TIME\_FILT & F\_HIGH\_PASS & F\_LOW\_PASS & WAVETYPE & 0 & 0 & EPRECOND & EPSILON\_WE\\
\begin{tabular}{llllllllllIlll}
\# & INV\_VS & INV\_VP & INV\_RHO & PRO & TIME\_FILT & F\_HIGH\_PASS & F\_LOW\_PASS & WAVETYPE & 0 & 0 & EPRECOND & EPSILON\_WE & GAMMA \\
\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.
......@@ -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,13 +654,14 @@ 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:
\begin{itemize}
\item To select frequency domain least squares (fdlsq), apply offset dependent weights and a waterlevel use\\
\textit{fdlsq:exp=1.0:waterlevel=0.01}
\item For tapering the inverted STF add the option \textit{irtap=t1;t2;t3;t4} to the parameter string. The four values define the taper window. %It is important to use ";", because "," is used by the JSON-parser to separate variables.
\end{itemize}
In most cases the frequency domain least squares engine is the best approach to find a suitable wavelet. There are also other possibilities, if you want to use those a detailed look at the libraries and the acompanying documentation provided in the folder /contrib is recommended. Here only the two main parameters for the fdlsq approach are described.
......@@ -718,7 +719,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.
......@@ -730,7 +731,7 @@ If you are using frequeny filtering (TIME\_FILT==1) during the inversion, you ca
"PICKS_FILE" : "./picked_times/picks"
"TWLENGTH_PLUS" : "0.01",
"TWLENGTH_MINUS" : "0.01",
"GAMMA" : "100000",
"GAMMA" : "30",
\end{verbatim}}}
{\color{red}{\begin{verbatim}
......@@ -738,9 +739,9 @@ Default values are:
TIMEWIN=0
\end{verbatim}}}
To apply time windowing in a time series the paramter 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.
......@@ -798,7 +799,7 @@ SWS_TAPER_FILE=0
SWS_TAPER_FILE_PER_SHOT=0
\end{verbatim}}}
Different preconditioning matrices can be created and applied to the gradients (using the function \texttt{taper\_grad.c}). To apply a vertical taper one has to set the switch SWS\_TAPER\_GRAD\_VERT to one and for a horizontaltaper SWS\_TAPER\_GRAD\_HOR has to be 1. The parameters for the vertical and the horizontal window are defined by the input file paramters GRADT1, GRADT2, GRADT3 and GRADT4. Please have a look at the function \texttt{taper\_grad.c} directly to obtain more information about the actual definition of the tapers. It is also possible to apply cylindrical tapers around the source positions. This can be done by either setting the switch SWS\_TAPER\_GRAD\_SOURCES or SWS\_TAPER\_CIRCULAR\_PER\_SHOT to 1. If one uses SWS\_TAPER\_GRAD\_SOURCES=1 only the final gradients (that means the gradients obtained by the summation of the gradients of each shots) are multiplied with a taper that decreases the gradients at all shot positions. Therefore, one looses the update information at the source positions. To avoid this one can use SWS\_TAPER\_CIRCULAR\_PER\_SHOT=1. In this case the gradients of the single shots are preconditioned with a window that only decreases at the current shot position. This is done before the summation of all gradients to keep model update information also at the shot positions. The actual tapers are generated by the function \texttt{taper\_grad.c} and \texttt{taper\_grad\_shot.c}, respectively. The circular taper around the source positions decrease from a value of one at the edge of the taper to a value of zero at the source position. The shape of the decrease can be defined by an error function (SRTSHAPE=1) or a log-function (SRTSHAPE=2). The radius of the taper is defined in meter by SRTRADIUS. Note, that this radius must be at least 5 gridpoints. With the parameter FILTSIZE one can extend the region where the taper is zero around the source. The taper is set to zero in a square region of (2*FILTSIZE+1 times 2*FILTSIZE+1) gridpoints. All preconditioning matrices that are applied are saved in the par directory with the file names taper\_coeff\_vert.bin, taper\_coeff\_horz.bin and taper\_coeff\_sources.bin.\\
Different preconditioning matrices can be created and applied to the gradients (using the function \texttt{taper\_grad.c}). To apply a vertical taper one has to set the switch SWS\_TAPER\_GRAD\_VERT to one and for a horizontaltaper SWS\_TAPER\_GRAD\_HOR has to be 1. The parameters for the vertical and the horizontal window are defined by the input file parameters GRADT1, GRADT2, GRADT3 and GRADT4. Please have a look at the function \texttt{taper\_grad.c} directly to obtain more information about the actual definition of the tapers. It is also possible to apply cylindrical tapers around the source positions. This can be done by either setting the switch SWS\_TAPER\_GRAD\_SOURCES or SWS\_TAPER\_CIRCULAR\_PER\_SHOT to 1. If one uses SWS\_TAPER\_GRAD\_SOURCES=1 only the final gradients (that means the gradients obtained by the summation of the gradients of each shots) are multiplied with a taper that decreases the gradients at all shot positions. Therefore, one looses the update information at the source positions. To avoid this one can use SWS\_TAPER\_CIRCULAR\_PER\_SHOT=1. In this case the gradients of the single shots are preconditioned with a window that only decreases at the current shot position. This is done before the summation of all gradients to keep model update information also at the shot positions. The actual tapers are generated by the function \texttt{taper\_grad.c} and \texttt{taper\_grad\_shot.c}, respectively. The circular taper around the source positions decrease from a value of one at the edge of the taper to a value of zero at the source position. The shape of the decrease can be defined by an error function (SRTSHAPE=1) or a log-function (SRTSHAPE=2). The radius of the taper is defined in meter by SRTRADIUS. Note, that this radius must be at least 5 gridpoints. With the parameter FILTSIZE one can extend the region where the taper is zero around the source. The taper is set to zero in a square region of (2*FILTSIZE+1 times 2*FILTSIZE+1) gridpoints. All preconditioning matrices that are applied are saved in the par directory with the file names taper\_coeff\_vert.bin, taper\_coeff\_horz.bin and taper\_coeff\_sources.bin.\\
To apply an externally defined taper on the gradients in IFOS2D, the parameter SWS\_TAPER\_FILE has to be set to 1. Each model parameter requires a taper file which needs to be located in TAPER\_FILE\_NAME.vp for vp, in TAPER\_FILE\_NAME.vs for vs and in TAPER\_FILE\_NAME.rho for the density.\\
......@@ -844,6 +845,9 @@ For GRAD\_FILT\_WAVELENGTH = 1 (and TIME\_FILT=1) a new wavelength dependent fil
\end{equation}
where F\_LOW\_PASS is the corner frequency of TIME\_FILT and A is a weighting factor.
\ \\
WARNING: The option GRAD\_FILTER uses a median filter that makes use of a quicksort routine. Under some (unknown) conditions this can take a huge amount of time and corrupt the inversion.
\section{Model manipulation}
\subsection{Limits for the model parameters}
......@@ -906,4 +910,7 @@ 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.
\ \\
WARNING: The option MODEL\_FILTER uses a median filter that makes use of a quicksort routine. Under some (unknown) conditions this can take a huge amount of time and corrupt the inversion.
# VS VP RHO PRO F_FILT F_HIGH F_LOW WT J J EPRE EPSI
1 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
3 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
# 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 20
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 5
4 1 0 0 0.01 0 0 0 2 0 0 0 0.005 0
This source diff could not be displayed because it is too large. You can view the blob instead.
......@@ -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 int EPRECOND;
extern float EPSILON_WE;
extern float GAMMA;
extern int GRAD_METHOD;
extern int WORKFLOW_STAGE;
......@@ -126,10 +127,12 @@ void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[ST
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)");
}
EPRECOND=workflow[WORKFLOW_STAGE][12];
EPSILON_WE=workflow[WORKFLOW_STAGE][13];
GAMMA=workflow[WORKFLOW_STAGE][14];
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);
}
}
\ No newline at end of file
}
......@@ -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);
}
......
......@@ -534,6 +534,8 @@ float average_matrix(float ** matrix);
float global_maximum(float ** gradiant_1);
void write_matrix_disk(float ** gradient,char path_name[STRING_SIZE]);
float matrix_product(float ** matrix1, float **matrix2);
void get_local_from_global_matrix(float ** global_matrix,float ** local_matrix);
float ** get_global_from_local_matrix(float ** local_matrix);
/* L-BFGS */
void lbfgs(float **grad1, float **grad2, float **grad3,float Vs_avg,float rho_avg,float Vp_avg, float *bfgsscale, float **bfgsmod, float **bfgsgrad,int bfgsnum,int bfgspar, int iteration, int * LBFGS_iter_start);
......
......@@ -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;
......@@ -21,7 +21,7 @@
#include "fd.h"
void write_matrix_disk(float ** gradient,char path_name[STRING_SIZE]){
void write_matrix_disk(float ** local_matrix,char path_name[STRING_SIZE]){
char joint[225];
FILE *FPjoint;
extern int POS[3],MYID;
......@@ -32,7 +32,7 @@ void write_matrix_disk(float ** gradient,char path_name[STRING_SIZE]){
for (i=1;i<=NX;i=i+IDX){
for (j=1;j<=NY;j=j+IDY){
fwrite(&gradient[j][i],sizeof(float),1,FPjoint);
fwrite(&local_matrix[j][i],sizeof(float),1,FPjoint);
}
}
......@@ -122,6 +122,70 @@ float matrix_product(float ** matrix1, float **matrix2) {
return global_sum;
}
float ** get_global_from_local_matrix(float ** local_matrix) {
extern int NXG, NYG;
extern int NX,NY;
extern int POS[3];
float ** global_matrix=NULL,** global_matrix_temp=NULL;
int i=0,j=0;
int ii=0, jj=0;
/* Allocate global matrix temp */
global_matrix_temp=matrix(1,NYG,1,NXG);
if(global_matrix_temp==NULL) {
declare_error("Allocation of global_matrix_temp in get_global_from_local_matrix failed!");
}
/* Allocate global matrix */
/* You have to deallocate this matrix on our own */
global_matrix=matrix(1,NYG,1,NXG);
if(global_matrix==NULL) {
declare_error("Allocation of global_matrix in get_global_from_local_matrix failed!");
}
/* Store local matrix in global matrix */
for (i=1;i<=NXG;i++){
for (j=1;j<=NYG;j++){
if ( (POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY)) ) {
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
global_matrix_temp[j][i]=local_matrix[jj][ii];
}
}
}
MPI_Allreduce(&global_matrix_temp[1][1],&global_matrix[1][1],NXG*NYG,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
free_matrix(global_matrix_temp,1,NYG,1,NXG);
return global_matrix;
}
void get_local_from_global_matrix(float ** global_matrix,float ** local_matrix) {
extern int NXG, NYG;
extern int NX,NY;
extern int POS[3];
int i=0,j=0;
int ii=0, jj=0;
/* Store local matrix in global matrix */
for (i=1;i<=NXG;i++){
for (j=1;j<=NYG;j++){
if ( (POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY)) ) {
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
local_matrix[jj][ii]=global_matrix[j][i];
}
}
}
}
......@@ -25,6 +25,7 @@
float *rd_sour(int *nts,FILE* fp_source){
extern int VERBOSE;
/* local variables */
float *psource;
int i, c;
......@@ -35,7 +36,7 @@ float *rd_sour(int *nts,FILE* fp_source){
while ((c=fgetc(fp_source)) != EOF)
if (c=='\n') ++(*nts);
rewind(fp_source);
printf(" Number of samples (nts) in source file: %i\n",*nts);
if (VERBOSE==1) printf(" Number of samples (nts) in source file: %i\n",*nts);
psource=vector(1,*nts);
for (i=1;i<=*nts;i++) fscanf(fp_source,"%e\n",&psource[i]);
......
......@@ -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;
......@@ -462,7 +463,9 @@ void read_par_json(FILE *fp, char *fileinp){
declare_error("Variable FORWARD_ONLY could not be retrieved from the json input file!");
else {
if (FORWARD_ONLY==0) { /* FWI is calculated */
/* Overwrite IDX/IDY option from forward modeling (used for snapshots), interpolation for FWI not yet implemented*/
IDX=1;
IDY=1;
/* General inversion parameters */
if (get_int_from_objectlist("ITERMAX",number_readobjects,&ITERMAX,varname_list, value_list))
declare_error("Variable ITERMAX could not be retrieved from the json input file!");
......@@ -763,6 +766,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;
......
......@@ -104,7 +104,7 @@ void readmod_elastic(float ** rho, float ** pi, float ** u){
}
if(feof(fp_vs) && feof(fp_rho)){
if(feof(fp_vs) || feof(fp_rho)){
declare_error("Model file VS or RHO is to small. Check dimensions NX*NY of file.");
}
......@@ -141,7 +141,7 @@ void readmod_elastic(float ** rho, float ** pi, float ** u){
fread(&vs, sizeof(float), 1, fp_vs);
fread(&rho, sizeof(float), 1, fp_rho);
if(!feof(fp_vs) && !feof(fp_rho)){
if(!feof(fp_vs) || !feof(fp_rho)){
declare_error("Model file VS or RHO is to big. Check dimensions NX*NY of file.");
}
fclose(fp_vs);
......
......@@ -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,25 +32,43 @@ 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, ADJOINT_TYPE;;
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);
if(ADJOINT_TYPE==1) {
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(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) {
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:
sprintf(vxf,"%s_vx.su.shot%d_adjoint_src",SEIS_FILE,ishot);
sprintf(vyf,"%s_vy.su.shot%d_adjoint_src",SEIS_FILE,ishot);
if(ADJOINT_TYPE==1) {
sprintf(vxf,"%s_vx.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) {
sprintf(vzf,"%s_vz.su.shot%d_adjoint_src",SEIS_FILE,ishot);
}
......@@ -60,14 +78,63 @@ 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){
if(ADJOINT_TYPE==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(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) {
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){
if(ADJOINT_TYPE==1) {
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(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) {
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:
if(ADJOINT_TYPE==1) {
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(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) {
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:
......@@ -76,8 +143,17 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
}
if(FORWARD_ONLY==1){
sprintf(vxf,"%s_vx.su.shot%d",SEIS_FILE,ishot);
sprintf(vyf,"%s_vy.su.shot%d",SEIS_FILE,ishot);
if(ADJOINT_TYPE==1) {
sprintf(vxf,"%s_vx.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) {
sprintf(vzf,"%s_vz.su.shot%d",SEIS_FILE,ishot);
}
......@@ -89,10 +165,21 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
switch (SEISMO){
case 1 : /* particle velocities only */
if (WAVETYPE==1 || WAVETYPE==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 (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==1) {
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 (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==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(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vz) to\n\t %s \n",0,ntr,vzf);
......@@ -114,10 +201,21 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
case 4 : /* everything */
if (WAVETYPE==1 || WAVETYPE==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 (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==1) {
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 (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==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);
outseis_glob(fp,fopen(pf,"w"), 0, sectionp,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
......@@ -134,10 +232,21 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
case 5 : /* everything except curl and div */
if (WAVETYPE==1 || WAVETYPE==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 (vy) to\n\t %s \n",0,ntr<