Skip to content

GitLab

IFOS2D
IFOS2D
Commit bcca82b5 authored by niklas.thiel's avatar niklas.thiel
Browse files
Options

Merge branch 'develop' into pup 

Merge branch 'develop' of git.scc.kit.edu:GPIAG-Software/IFOS2D into pup

Conflicts:
	src/IFOS2D.c
	src/exchange_par.c
parents d8de885d 675ff853
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 5120 additions and 4821 deletions
README.md
View file @ bcca82b5
# What is IFOS2D?
**IFOS2D** (**I**nversion of **F**ull **O**bserved **S**eismograms) is a 2D elastic full waveform inversion code.
The inversion problem is solved by a conjugate gradient method and the gradients are computed in the time domain with the adjoint method.
The forward modeling is done by a time domain finite difference scheme.
**IFOS2D** (**I**nversion of **F**ull **O**bserved **S**eismograms) is a 2-D elastic full waveform inversion code.
The inversion problem is solved by a conjugate gradient method and the gradients are computed in the time domain by the adjoint state method.
The forward modeling is done by a time domain Finite-Difference scheme.
IFOS2D is the reverse (inverse) of our 2D finite difference forward solver [**SOFI2D**](https://git.scc.kit.edu/GPIAG-Software/SOFI2D).
IFOS2D is the reverse (inverse) of our 2-D Finite-Difference forward solver [**SOFI2D**](https://git.scc.kit.edu/GPIAG-Software/SOFI2D).
The [**manual**](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/wikis/home) is included in the download archive or can be downloaded [here](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/wikis/home).
# Download and Newsletter
You can Download the [**latest Release**](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/tags/Release_2.0) or the current [**Beta-Version**](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/tree/master)
You can download the [**latest Release**](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/tags/Release_2.0.1) or the current [**Beta-Version**](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/tree/develop).
To receive news and updates please [register](http://www.gpi.kit.edu/Software-FWI.php) on the email list [IFOS@lists.kit.edu](http://www.gpi.kit.edu/Software-FWI.php).
Please use this list also to ask questions on using the software or to report problems or bugs.
\ No newline at end of file
Please use this list also to ask questions or to report problems or bugs.
\ No newline at end of file
doc/latex/.gitignore
View file @ bcca82b5
......@@ -2,6 +2,7 @@
# Created by https://www.gitignore.io/api/latex
### LaTeX ###
.texpadtmp/
*.acn
*.acr
*.alg
......
doc/latex/2_Theoretical_Background.tex
View file @ bcca82b5
......@@ -3,7 +3,7 @@
\chapter{Theoretical Background}
%------------------------------------------------------------------------------------------------%
This chapter is mainly based on the dissertation of Daniel Köhn (\cite{koehn:11}), who originally has written this manual.
\section{Equations of motion for an elastic medium}\label{elastic_fd_model}
The propagation of waves in a general elastic medium can be described by a system of coupled linear partial differential equations. They consist of the equations of motion
\EQ{m:1}{\begin{split}
......@@ -27,7 +27,6 @@ This expression is called {\bf{Stress-Displacement}} formulation. Another common
\rm{\frac{\partial \epsilon_{ij}}{\partial t}}&\rm{=\frac{1}{2}\biggl(\frac{\partial v_i}{\partial x_j}+\frac{\partial v_j}{\partial x_i}\biggr)}
\end{split}}
This expression is called {\bf{Stress-Velocity}} formulation. For simple cases \ER{2:20} and \ER{2:20:1} can be solved analytically. More complex problems require numerical solutions. One possible approach for a numerical solution is described in the next section.
\newpage
\section{Solution of the elastic wave equation by finite differences}\label{elastic_FD_Code}
\subsection{Discretization of the wave equation}
For the numerical solution of the elastic equations of motion, Eqs. \ER{2:20} have to be discretized in time and space on a grid. The
......@@ -182,7 +181,7 @@ A comparison between the exponential damping and the PML boundary is shown in Fi
\begin{figure}[ht]
\begin{center}
\includegraphics[width=13cm]{figures/ABS_PML_comp_shots.pdf}
\caption{\label{comp_EXP_PML} Comparison between exponential damping (left column) and PML (right column) absorbing boundary conditions for a homogeneous full space model.}
\caption{\label{comp_EXP_PML} Comparison between exponential damping (left column) and PML (right column) absorbing boundary conditions for a homogeneous full space model. (\cite{koehn:11})}
\end{center}
\end{figure}
\end{enumerate}
......@@ -251,7 +250,7 @@ Holberg) of FD operators. For the Holberg coefficients n is calculated for a min
\epsfig{file=figures/homogenous_grid_n_16_5.pdf, width=6 cm}\epsfig{file=figures/homogenous_grid_n_16_10.pdf, width=6 cm}
\epsfig{file=figures/homogenous_grid_n_4_5.pdf, width=6 cm}\epsfig{file=figures/homogenous_grid_n_4_10.pdf, width=6 cm}
\epsfig{file=figures/homogenous_grid_n_2_5.pdf, width=6 cm}\epsfig{file=figures/homogenous_grid_n_2_10.pdf, width=6 cm}
\caption{\label{grid_disp_pics} The influence of grid dispersion in FD modeling: Spatial sampling of the wavefield using n=16 (top), n=4 (center) and n=2 gridpoints (bottom) per minimum wavelength $\rm{\lambda_{min}}$.}
\caption{\label{grid_disp_pics} The influence of grid dispersion in FD modeling: Spatial sampling of the wavefield using n=16 (top), n=4 (center) and n=2 gridpoints (bottom) per minimum wavelength $\rm{\lambda_{min}}$. (\cite{koehn:11})}
\end{center}
\end{figure}
\clearpage
......@@ -284,7 +283,7 @@ FDORDER & h (Taylor) & h (Holberg) \\ \hline
\begin{center}
\epsfig{file=figures/courandt_1.pdf, width=7 cm}\epsfig{file=figures/courandt_2.pdf, width=7 cm}
\epsfig{file=figures/courandt_3.pdf, width=7 cm}\epsfig{file=figures/courandt_4.pdf, width=7 cm}
\caption{\label{courandt_pics} Temporal evolution of the Courant instability. In the colored areas the wave amplitudes are extremly large.}
\caption{\label{courandt_pics} Temporal evolution of the Courant instability. In the colored areas the wave amplitudes are extremly large. (\cite{koehn:11})}
\end{center}
\end{figure}
\clearpage
doc/latex/3_Adjoint_Problem.tex
View file @ bcca82b5
% ------------------------------------
\chapter{The adjoint problem}
% ------------------------------------
This chapter is mainly based on the dissertation of Daniel Köhn (\cite{koehn:11}), who originally has written this manual.\\
\newline
The aim of full waveform tomography is to find an ''optimum'' model which can explain the data very well. It should not only explain the
first arrivals of specific phases of the seismic wavefield like refractions or reflections, but also the amplitudes which contain information
on the distribution of the elastic material parameters in the underground. To achieve this goal three problems have to be solved:
......@@ -31,7 +32,7 @@ has a special physical meaning. It represents the residual elastic energy contai
\begin{figure}[!bh]
\begin{center}
\includegraphics[width=15 cm]{figures/data_res_sketch_1.pdf}
\caption{Definition of data residuals $\rm{\mathbf{\delta u}}$.}
\caption{Definition of data residuals $\rm{\mathbf{\delta u}}$. (\cite{koehn:11})}
\label{sketch_data_res}
\end{center}
\end{figure}
......@@ -49,7 +50,7 @@ along the search direction $\rm{\mathbf{\delta m}_{1}}$ with the step length $\r
\begin{figure}[!bh]
\begin{center}
\includegraphics[width=12.5cm]{figures/sketch_grad_1.pdf}
\caption{Schematic sketch of the residual energy at one point in space as a function of two model parameters $\rm{m_1}$ and $\rm{m_2}$. The blue dot denotes the starting point in the parameter space, while the red cross marks a minimum of the objective function.}
\caption{Schematic sketch of the residual energy at one point in space as a function of two model parameters $\rm{m_1}$ and $\rm{m_2}$. The blue dot denotes the starting point in the parameter space, while the red cross marks a minimum of the objective function. (\cite{koehn:11})}
\label{sketch_grad}
\end{center}
\end{figure}
......@@ -105,7 +106,7 @@ Eq. \ER{dL2_dm} can be related to the mapping of small changes from the data to
\begin{figure}[ht]
\begin{center}
\includegraphics[width=10cm]{figures/mapping_data_model.pdf}
\caption{Mapping between model and data space and vice versa.}
\caption{Mapping between model and data space and vice versa. (\cite{koehn:11})}
\label{mapping_data_model}
\end{center}
\end{figure}
......@@ -424,7 +425,7 @@ The aim is to find the minimum of this function loacted at the point [1,1] which
\begin{center}
\includegraphics[width=15cm]{figures/Rosenbrock_1.pdf}\\
\includegraphics[width=15cm]{figures/Rosenbrock_2.pdf}\\
\caption{Results of the convergence test for the Rosenbrock function. The minimum is marked with a red cross, the starting point with a blue point. The maximum number of iterations is 16000. The step length $\rm{\mu_n}$ varies between $\rm{2e-3}$ (top) and $\rm{6.1e-3}$ (bottom).}
\caption{Results of the convergence test for the Rosenbrock function. The minimum is marked with a red cross, the starting point with a blue point. The maximum number of iterations is 16000. The step length $\rm{\mu_n}$ varies between $\rm{2e-3}$ (top) and $\rm{6.1e-3}$ (bottom). (\cite{koehn:11})}
\label{Rosenbrock_constant}
\end{center}
\end{figure}
......@@ -432,7 +433,7 @@ gradient of the Rosenbrock function is large. After reaching the narrow valley t
\begin{figure}[ht]
\begin{center}
\includegraphics[width=15cm]{figures/sl_case1_final}
\caption{Line search algorithm to find the optimum step length $\rm{\mu_{opt}}$: The true misfit function (yellow line) is approximated by a parabola fitted by 3 points.}
\caption{Line search algorithm to find the optimum step length $\rm{\mu_{opt}}$: The true misfit function (yellow line) is approximated by a parabola fitted by 3 points. (\cite{koehn:11})}
\label{sl_case1_final}
\end{center}
\end{figure}
......@@ -502,7 +503,7 @@ This approach leads to a smoother decrease of the objective function, but also i
\begin{figure}[ht]
\begin{center}
\includegraphics[width=16cm]{figures/Rosenbrock_3}
\caption{Results of the convergence test for the Rosenbrock function. The minimum is marked by a red cross, the starting point by a blue point. The maximum number of iterations is 4000. The optimum step length is calculated at each iteration by the parabola fitting algorithm. Note the criss-cross pattern of the updates in the narrow valley near the minimum.}
\caption{Results of the convergence test for the Rosenbrock function. The minimum is marked by a red cross, the starting point by a blue point. The maximum number of iterations is 4000. The optimum step length is calculated at each iteration by the parabola fitting algorithm. Note the criss-cross pattern of the updates in the narrow valley near the minimum. (\cite{koehn:11})}
\label{Rosenbrock_variable}
\end{center}
\end{figure}
......@@ -544,7 +545,7 @@ I use the very popular choice $\rm{\beta_n = max[0,\beta_n^{PR}]}$ which provide
\begin{figure}[ht]
\includegraphics[width=17cm]{figures/Rosenbrock_4}\\
\caption{Results of the convergence test for the Rosenbrock function using the conjugate gradient method, where the optimum step length is calculated with the parabolic fitting algorithm. The minimum is marked by a red cross, the starting point by a blue point. The maximum number of iterations is 2000.}
\caption{Results of the convergence test for the Rosenbrock function using the conjugate gradient method, where the optimum step length is calculated with the parabolic fitting algorithm. The minimum is marked by a red cross, the starting point by a blue point. The maximum number of iterations is 2000. (\cite{koehn:11})}
\label{Rosenbrock_cg}
\end{figure}
......
doc/latex/6_Parameter_Definition.tex
View file @ bcca82b5
This diff is collapsed.
doc/latex/manual_IFOS2D.tex
View file @ bcca82b5
......@@ -112,10 +112,12 @@
\section*{Authors}
The IFOS2D code (formerly DENISE) was at first developed by Daniel K\"ohn, Denise De Nil and Andr$\rm{\acute{e}}$ Kurzmann at the Christian-Albrechts-Universit\"at Kiel and TU Bergakademie Freiberg (Germany) from 2005 to 2009.\\
\newline
This documentation was originally written by Daniel K\"ohn.\\ Large parts of the shown theory is extracted from his dissertation \cite{koehn:11}.\\
\newline
The forward code is based on the viscoelastic FD code fdveps (now SOFI2D) by \cite{bohlen:02}.\\
\newline
Different external libraries for timedomain filtering are used.\\
Different external libraries for time domain filtering are used.\\
The copyright of the source codes are held by different persons:\\
\newline
cseife.c, cseife.h, lib\_stfinv, lib\_aff, lib\_fourier:\\
......@@ -160,13 +162,14 @@ The Authors of IFOS2D are listed in file \path{AUTHORS}.
%------------------------------------------------------------------------------------------------%
We thank for constructive discussions and further code improvements:\\
\newline
\newline
Daniel K\"ohn (Christian-Albrechts-Universit\"at Kiel), \\
Anna Przebindowska (Karlsruhe Institute of Technology), \\
Olaf Hellwig (TU Bergakademie Freiberg), \\
Dennis Wilken and Wolfgang Rabbel (Christian-Albrechts-Universit\"at Kiel).
\newline
\noindent The development of the code was suppported by the Christian-Albrechts-Universität Kiel, TU Bergakademie Freiberg, Deutsche Forschungsgemeinschaft (DFG), Bundesministerium für Bildung und Forschung (BMBF), the Wave Inversion Technology (WIT) Consortium and the Verbundnetz-Gas AG (VNG).
\noindent The development of the code was supported by the Christian-Albrechts-Universität Kiel, TU Bergakademie Freiberg, Deutsche Forschungsgemeinschaft (DFG), Bundesministerium für Bildung und Forschung (BMBF), the Wave Inversion Technology (WIT) Consortium and the Verbundnetz-Gas AG (VNG).
\noindent The code was tested and optimized at the computing centres of Kiel University, TU Bergakademie Freiberg, TU Chemnitz, TU Dresden, the Karlsruhe Institute of Technology (KIT) and the Hochleistungsrechenzentrum Nord (HLRN 1+2).
\newline
......
par/.gitignore
View file @ bcca82b5
*.dat
*.txt
*.sh
*shot*
\ No newline at end of file
*.*
*~
\ No newline at end of file
par/in_and_out/IFOS2D_INV_all_parameters.json
View file @ bcca82b5
......@@ -162,17 +162,20 @@
"N_STF" : "10",
"N_STF_START" : "1",
"TAPER_STF" : "0",
"TRKILL_STF" : "0",
"TRKILL_FILE_STF" : "./trace_kill/trace_kill",
"TRKILL_STF_OFFSET" : "0",
"TRKILL_STF_OFFSET_LOWER , TRKILL_STF_OFFSET_UPPER" : "0 , 10",
"TRKILL_FILE_STF" : "./trace_kill/trace_kill",
"TRKILL_STF_OFFSET_INVERT" : "0",
"Frequency filtering during inversion" : "comment",
"TIME_FILT" : "0",
"F_HP" : "1",
"FC_START" : "10.0",
"FC_END" : "75.0",
"FC_INCR" : "10.0",
"F_HIGH_PASS" : "1",
"F_LOW_PASS_START" : "10.0",
"F_LOW_PASS_END" : "75.0",
"F_LOW_PASS_INCR" : "10.0",
"ORDER" : "2",
"ZERO_PHASE" : "0",
"FREQ_FILE" : "frequencies.dat",
......@@ -191,9 +194,10 @@
"Trace killing" : "comment",
"TRKILL" : "0",
"TRKILL_FILE" : "./trace_kill/trace_kill",
"TRKILL_OFFSET" : "0",
"TRKILL_OFFSET_LOWER , TRKILL_OFFSET_UPPER" : "0 , 10",
"TRKILL_FILE" : "./trace_kill/trace_kill",
"Definition of a gradient taper" : "comment",
"SWS_TAPER_GRAD_VERT" : "0",
......
par/in_and_out/toy_example/toy_example_INV.json
View file @ bcca82b5
......@@ -127,8 +127,8 @@
"Frequency filtering during inversion" : "comment",
"TIME_FILT" : "1",
"FC_START" : "10.0",
"FC_END" : "70.0",
"FC_INCR" : "10.0",
"F_LOW_PASS_END" : "70.0",
"F_LOW_PASS_INCR" : "10.0",
"ORDER" : "4",
"Minimum number of iteration per frequency" : "comment",
......
par/in_and_out/toy_example/toy_example_INV_SH.json
View file @ bcca82b5
......@@ -139,9 +139,9 @@
"Frequency filtering during inversion" : "comment",
"TIME_FILT" : "0",
"FC_START" : "5.0",
"FC_END" : "5.0",
"FC_INCR" : "10.0",
"F_LOW_PASS_START" : "5.0",
"F_LOW_PASS_END" : "5.0",
"F_LOW_PASS_INCR" : "10.0",
"ORDER" : "4",
"Minimum number of iteration per frequency" : "comment",
......
par/in_and_out/toy_example/toy_example_ac_INV.json
View file @ bcca82b5
......@@ -136,10 +136,10 @@
"Frequency filtering during inversion" : "comment",
"TIME_FILT" : "0",
"F_HP" : "1",
"FC_START" : "10.0",
"FC_END" : "70.0",
"FC_INCR" : "10.0",
"F_HIGH_PASS" : "1",
"F_LOW_PASS_START" : "10.0",
"F_LOW_PASS_END" : "70.0",
"F_LOW_PASS_INCR" : "10.0",
"ORDER" : "4",
"ZERO_PHASE" : "1",
......
src/IFOS2D.c
View file @ bcca82b5
This diff is collapsed.
src/PCG.c
View file @ bcca82b5
......@@ -24,7 +24,7 @@
#include "fd.h"
void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, float C_vp, float ** gradp, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float FC){
void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, float C_vp, float ** gradp, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS, int PCG_iter_start){
extern int NX, NY, IDX, IDY, SPATFILTER, GRAD_FILTER;
extern int FORWARD_ONLY, SWS_TAPER_GRAD_VERT, SWS_TAPER_GRAD_HOR, SWS_TAPER_GRAD_SOURCES, SWS_TAPER_FILE;
......@@ -41,8 +41,14 @@ void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int
int use_conjugate_2=1;
FILE *FP3, *FP4, *FP6, *FP5 = NULL;
/* Check if conjugate gradient can be used */
if( (iter-PCG_iter_start) < 2 ) {
use_conjugate_2=0;
if( (iter-PCG_iter_start) < 1 ) {
use_conjugate_1=0;
}
}
/* ===================================================================================================================================================== */
/* ===================================================== GRADIENT ZP ================================================================================== */
/* ===================================================================================================================================================== */
......@@ -99,7 +105,7 @@ void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int
spat_filt(waveconv,iter,1);}
/* apply 2D-Gaussian filter*/
if(GRAD_FILTER==1){smooth(waveconv,1,1,Vs_avg,FC);}
if(GRAD_FILTER==1){smooth(waveconv,1,1,Vs_avg,F_LOW_PASS);}
/* output of the preconditioned gradient */
for (i=1;i<=NX;i=i+IDX){
......@@ -331,7 +337,7 @@ void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int
spat_filt(waveconv_u,iter,2);}
/* apply 2D-Gaussian filter*/
if(GRAD_FILTER==1){smooth(waveconv_u,2,1,Vs_avg,FC);}
if(GRAD_FILTER==1){smooth(waveconv_u,2,1,Vs_avg,F_LOW_PASS);}
/* output of the preconditioned gradient */
for (i=1;i<=NX;i=i+IDX){
......@@ -558,7 +564,7 @@ void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int
spat_filt(waveconv_rho,iter,3);}
/* apply 2D-Gaussian filter*/
if(GRAD_FILTER==1){smooth(waveconv_rho,3,1,Vs_avg,FC);}
if(GRAD_FILTER==1){smooth(waveconv_rho,3,1,Vs_avg,F_LOW_PASS);}
/* output of the preconditioned gradient */
for (i=1;i<=NX;i=i+IDX){
......
src/PCG_SH.c
View file @ bcca82b5
......@@ -24,7 +24,7 @@
#include "fd.h"
void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float FC){
void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS, int PCG_iter_start){
extern int NX, NY, IDX, IDY, SPATFILTER, GRAD_FILTER;
extern int FORWARD_ONLY, SWS_TAPER_GRAD_VERT, SWS_TAPER_GRAD_HOR, SWS_TAPER_GRAD_SOURCES, SWS_TAPER_FILE;
......@@ -40,6 +40,14 @@ void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int
int use_conjugate_2=1;
/* Check if conjugate gradient can be used */
if( (iter-PCG_iter_start) < 2 ) {
use_conjugate_2=0;
if( (iter-PCG_iter_start) < 1 ) {
use_conjugate_1=0;
}
}
/* ===================================================================================================================================================== */
/* ===================================================== GRADIENT Zs ================================================================================== */
/* ===================================================================================================================================================== */
......@@ -91,7 +99,7 @@ void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int
spat_filt(waveconv_u,iter,2);}
/* apply 2D-Gaussian filter*/
if(GRAD_FILTER==1){smooth(waveconv_u,2,1,Vs_avg,FC);}
if(GRAD_FILTER==1){smooth(waveconv_u,2,1,Vs_avg,F_LOW_PASS);}
/* output of the preconditioned gradient */
for (i=1;i<=NX;i=i+IDX){
......@@ -320,7 +328,7 @@ void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int
spat_filt(waveconv_rho,iter,3);}
/* apply 2D-Gaussian filter*/
if(GRAD_FILTER==1){smooth(waveconv_rho,3,1,Vs_avg,FC);}
if(GRAD_FILTER==1){smooth(waveconv_rho,3,1,Vs_avg,F_LOW_PASS);}
/* output of the preconditioned gradient */
for (i=1;i<=NX;i=i+IDX){
......
src/apply_workflow.c
View file @ bcca82b5
......@@ -23,7 +23,7 @@
#include "fd.h"
void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[STRING_SIZE],int *iter,float *FC,int wavetype_start, int * change_wavetype_iter, int * LBFGS_iter_start){
void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[STRING_SIZE],int *iter,float *F_LOW_PASS,int wavetype_start, int * change_wavetype_iter, int * LBFGS_iter_start){
/* local variables */
int x;
......@@ -90,12 +90,12 @@ void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[ST
/* Frequency filtering */
if(TIME_FILT==1) {
TIME_FILT=workflow[WORKFLOW_STAGE][6];
if(*FC>workflow[WORKFLOW_STAGE][7]&&(workflow[WORKFLOW_STAGE][6]>0)) {
if(MYID==0)printf("\n Due to the abort criteriom FC is already higher than specified in workflow\n");
if(MYID==0)printf(" therefore instead of %.2f HZ FC=%.2f HZ is used\n",workflow[WORKFLOW_STAGE][7],*FC);
if(*F_LOW_PASS>workflow[WORKFLOW_STAGE][7]&&(workflow[WORKFLOW_STAGE][6]>0)) {
if(MYID==0)printf("\n Due to the abort criteriom F_LOW_PASS is already higher than specified in workflow\n");
if(MYID==0)printf(" therefore instead of %.2f HZ F_LOW_PASS=%.2f HZ is used\n",workflow[WORKFLOW_STAGE][7],*F_LOW_PASS);
} else {
if(*FC!=workflow[WORKFLOW_STAGE][7]) *LBFGS_iter_start=*iter;
*FC=workflow[WORKFLOW_STAGE][7];
if(*F_LOW_PASS!=workflow[WORKFLOW_STAGE][7]) *LBFGS_iter_start=*iter;
*F_LOW_PASS=workflow[WORKFLOW_STAGE][7];
}
} else {
if(MYID==0&&(workflow[WORKFLOW_STAGE][6]>0))printf("\n TIME_FILT cannot be activated due to it is not activated in the JSON File \n");
......
src/calc_energy.c
View file @ bcca82b5
......@@ -26,7 +26,7 @@ double calc_energy(float **sectiondata, int ntr, int ns, float energy, int ntr_g
/* declaration of variables */
extern float DT;
extern int MYID, USE_WORKFLOW, WORKFLOW_STAGE;
extern int TRKILL, NORMALIZE, FC, TIMEWIN;
extern int TRKILL, NORMALIZE, F_LOW_PASS, TIMEWIN;
extern char TRKILL_FILE[STRING_SIZE];
extern int VELOCITY;
int i, j;
......
src/calc_misfit.c
View file @ bcca82b5
......@@ -26,7 +26,7 @@ double calc_misfit(float **sectiondata, float **section, int ntr, int ns, int LN
/* declaration of variables */
extern float DT;
extern int MYID, USE_WORKFLOW, WORKFLOW_STAGE;
extern int TRKILL, NORMALIZE, FC, TIMEWIN;
extern int TRKILL, NORMALIZE, F_LOW_PASS, TIMEWIN;
extern char TRKILL_FILE[STRING_SIZE];
extern int VELOCITY;
......
src/calc_res.c
View file @ bcca82b5
......@@ -26,7 +26,7 @@ double calc_res(float **sectiondata, float **section, float **sectiondiff, float
/* declaration of variables */
extern float DT, WATERLEVEL_LNORM8;
extern int REC1, REC2, MYID, ACOUSTIC;
extern int TRKILL, NORMALIZE, FC, TIMEWIN;
extern int TRKILL, NORMALIZE, F_LOW_PASS, TIMEWIN;
extern char TRKILL_FILE[STRING_SIZE];
extern int VELOCITY, USE_WORKFLOW, WORKFLOW_STAGE;
float RMS, signL1;
......
src/exchange_par.c
View file @ bcca82b5
......@@ -60,7 +60,7 @@ void exchange_par(void){
extern int INV_STF, N_STF, N_STF_START;
extern char PARA[STRING_SIZE];
extern int TIME_FILT, ORDER, ZERO_PHASE,WRITE_FILTERED_DATA;
extern float FC_START, FC_END, FC_INCR, F_HP;
extern float F_LOW_PASS_START, F_LOW_PASS_END, F_LOW_PASS_INCR, F_HIGH_PASS;
extern int LNORM, DTINV;
extern int STEPMAX;
extern float EPS_SCALE, SCALEFAC;
......@@ -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;
......@@ -108,7 +109,6 @@ void exchange_par(void){
extern int EPRECOND_PER_SHOT_SH;
extern float LBFGS_SCALE_GRADIENTS;
extern int LBFGS_SURFACE;
extern int LBFGS_STEP_LENGTH;
extern int N_LBFGS;
......@@ -178,9 +178,9 @@ void exchange_par(void){
fdum[39] = npower;
fdum[40] = k_max_PML;
fdum[42] = FC_START;
fdum[43] = FC_END;
fdum[44] = FC_INCR;
fdum[42] = F_LOW_PASS_START;
fdum[43] = F_LOW_PASS_END;
fdum[44] = F_LOW_PASS_INCR;
fdum[45] = EPS_SCALE;
fdum[46] = SCALEFAC;
......@@ -206,7 +206,7 @@ void exchange_par(void){
fdum[58] = A;
fdum[59] = F_HP;
fdum[59] = F_HIGH_PASS;
fdum[60] = JOINT_INVERSION_PSV_SH_ALPHA_VS;
fdum[61] = JOINT_INVERSION_PSV_SH_ALPHA_RHO;
......@@ -364,7 +364,6 @@ void exchange_par(void){
idum[101]=EPRECOND;
idum[102]=EPRECOND_ITER;
idum[103]=LBFGS_SURFACE;
idum[104]=LBFGS_STEP_LENGTH;
idum[105]=EPRECOND_PER_SHOT;
......@@ -383,14 +382,11 @@ void exchange_par(void){
idum[114]=TRKILL_OFFSET;
idum[115]=TRKILL_STF_OFFSET;
idum[116]=WAVESEP;
idum[116]=TRKILL_STF_OFFSET_INVERT;
idum[117]=WAVESEP;
} /** if (MYID == 0) **/
if (MYID != 0) FL=vector(1,L);
MPI_Barrier(MPI_COMM_WORLD);
MPI_Bcast(&idum,NPAR,MPI_INT,0,MPI_COMM_WORLD);
......@@ -470,9 +466,9 @@ void exchange_par(void){
k_max_PML = fdum[40];
FC_START = fdum[42];
FC_END = fdum[43];
FC_INCR = fdum[44];
F_LOW_PASS_START = fdum[42];
F_LOW_PASS_END = fdum[43];
F_LOW_PASS_INCR = fdum[44];
EPS_SCALE = fdum[45];
SCALEFAC = fdum[46];
......@@ -499,7 +495,7 @@ void exchange_par(void){
A = fdum[58];
F_HP = fdum[59];
F_HIGH_PASS = fdum[59];
JOINT_INVERSION_PSV_SH_ALPHA_VS = fdum[60];
JOINT_INVERSION_PSV_SH_ALPHA_RHO = fdum[61];
......@@ -658,7 +654,6 @@ void exchange_par(void){
EPRECOND=idum[101];
EPRECOND_ITER=idum[102];
LBFGS_SURFACE=idum[103];
LBFGS_STEP_LENGTH=idum[104];
EPRECOND_PER_SHOT= idum[105];
......@@ -679,9 +674,15 @@ void exchange_par(void){
TRKILL_OFFSET=idum[114];
TRKILL_STF_OFFSET=idum[115];
WAVESEP=idum[116];
TRKILL_STF_OFFSET_INVERT=idum[116];
WAVESEP=idum[117];
MPI_Bcast(&FL[1],L,MPI_FLOAT,0,MPI_COMM_WORLD);
if ( MYID!=0 && L>0 ) {
FL=vector(1,L);
}
if ( L>0 ) {
MPI_Bcast(&FL[1],L,MPI_FLOAT,0,MPI_COMM_WORLD);
}
}
src/fd.h
View file @ bcca82b5
......@@ -58,7 +58,7 @@ void window_cos(float **win, int npad, int nsrc, float it1, float it2, float it3
void catseis(float **data, float **fulldata, int *recswitch, int ntr_glob, MPI_Comm newcomm_nodentr);
void stf(FILE *fp, float **sectionvy, float ** sectionvy_obs, float ** sectionvy_conv, float * source_time_function, int **recpos, int **recpos_loc,
int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots, float FC, int SH,int nsrc_glob);
int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots, float F_LOW_PASS, int SH,int nsrc_glob);
int **splitrec(int **recpos,int *ntr_loc, int ntr, int *recswitch);
......@@ -207,15 +207,15 @@ 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);
void output_source_signal(FILE *fp, float **signals, int ns, int seis_form);
void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, float C_vp, float ** gradp, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float FC);
void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, float C_vp, float ** gradp, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS, int PCG_iter_start);
void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float FC);
void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS, int PCG_iter_start);
void PML_pro(float * d_x, float * K_x, float * alpha_prime_x, float * a_x, float * b_x,
float * d_x_half, float * K_x_half, float * alpha_prime_x_half, float * a_x_half, float * b_x_half,
......@@ -308,7 +308,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);
......@@ -428,7 +428,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,
......@@ -454,7 +454,7 @@ void zero_fdveps_visc(int ny1, int ny2, int nx1, int nx2, float ** vx, float **
void FLnode(float ** rho, float ** pi, float ** u, float ** taus, float ** taup, float * eta);
void smooth(float ** mat, int sws, int filter, float Vs_avg, float FC);
void smooth(float ** mat, int sws, int filter, float Vs_avg, float F_LOW_PASS);
/* declaration of functions for parser*/
......@@ -545,7 +545,7 @@ void read_workflow(char file_in[STRING_SIZE],float *** workflow, int *workflow_l
float ** joint_inversion_grad ( float ** gradiant_1,float ** gradiant_2, float alpha, int joint_type);
void snap_SH(FILE *fp,int nt, int nsnap, float ** vz, float **u, float **pi, float *hc,int ishot);
void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[STRING_SIZE],int *iter,float *FC,int wavetype_start, int * change_wavetype_iter, int * LBFGS_iter_start);
void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[STRING_SIZE],int *iter,float *F_LOW_PASS,int wavetype_start, int * change_wavetype_iter, int * LBFGS_iter_start);
void eprecond(float ** W, float ** vx, float ** vy);
void eprecond_SH(float ** W, float ** vz);
......
src/filter_frequencies.c
View file @ bcca82b5
......@@ -31,7 +31,7 @@ extern FILE *FP;
FILE *freqf;
char cline[256];
int IT=1;
float *FC;
float *F_LOW_PASS;
/*----------------------------------- open FREQ_FILE ------------------------------*/
......@@ -48,8 +48,8 @@ fprintf(FP,"\n Reading frequencies from %s \n",FREQ_FILE);
fscanf(freqf,"%i",nfrq);
/*----------------------------alocate frequency array------------------------------*/
FC=vector(1,*nfrq);
//FC = malloc((*nfrq+1) * sizeof(float));
F_LOW_PASS=vector(1,*nfrq);
//F_LOW_PASS = malloc((*nfrq+1) * sizeof(float));
rewind (freqf);
fprintf(FP," Number of freqencies specified: %d \n",*nfrq);
......@@ -57,16 +57,16 @@ fprintf(FP," Number of freqencies specified: %d \n",*nfrq);
/*----------------------Read frequencies from FREQ_FILE----------------------------*/
for (IT=1;IT<=(*nfrq);IT++){
fgets(cline,255,freqf);
fscanf(freqf,"%f",&FC[IT]);
fscanf(freqf,"%f",&F_LOW_PASS[IT]);