Commit ca959455 authored by laura.gassner's avatar laura.gassner
Browse files

Merge branch 'develop' into feature/smoothing

parents 030f098d 20f473aa
......@@ -19,6 +19,10 @@ remainder may consist of several control parameters being separated by colons
or may come along with a parameter value. The value is separated from the
parameter by an equal sign (=).
Where several values in an argument to a parameter must be separated (like in
the 'irtap' option of the Fourier domain procedures) white space ( ), commas
(,), and semicolons (;) are allowed as field delimiters, at your convenience.
Examples:
- To select Fourier domain least squares and shift the returned source
correction filter wavelet by 0.4s and switch on verbose mode, pass the
......
......@@ -4,17 +4,31 @@
# 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)
in order to expose acausal components
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)
in order to expose acausal components
irtap=t1,t2,t3,t4 taper impulse response of correction filter
These options define the number of samples N used for the FFT (Fast Fourier
Transform). This number N should be larger than the number of samples M in the
original input time series to avoid wrap-around. If fpow2 is set, N will be
the next power of 2 larger than M*f. Else if fdiv is set, N will be the next
integer multiple of d larger than M*f.
integer multiple of d larger than M*f. If fdiv is not set explicitely, a
default value for d (commonly 100) is used. If option fpad ist used, N will be
f times larger than without padding. Without explicitely setting fpad, a
default value for f is used (which commonly equals 1.5). When defining the
number of samples N, first padding is considered (fpad), then the either
selection of a power of two (pow2) or the divisor criterion (fdiv) is applied.
The latter is only applied, if pow2 ist not selected.
Input time series with M samples will be padded with (N-M) zeros to create the
time series which actually will be transformed to the Fourier domain. Upon
inverse FFT the additional (N-M) samples of the resulting time series will be
discarded before returning the M remaining samples to the caller. Note, that
this is a form of implicite taper. In particular the caller will not obtain
exactly the filter response, which was used for convolution internally.
The derived correction filter in some cases can contain acausal components.
This means that the impulse response is non-zero for negative time values.
......@@ -22,4 +36,24 @@ Since by definition, the impulse response is output for the time interval of
the input data, these acausal components can remain unnoticed. The option
tshift can be used to shift the impulse response as obtained by inverse FFT in
order to expose acausal components.
A time domain taper can be applied to the impulse response of the correction
filter by using option irtap. Four time values are given in units of seconds:
t1, t2, t3, and t4. They must be in increasing order and (t4-t1) must be
smaller than the total duration of the time series used to represent signals
internally. Times value are allowed to be negative. Time series are understood
to be periodic (due to discrete Fourier transformation). Prior to application
of the correction filter to the time series passed to the algorithm, the
correction filter is transformed to the time domain, tapered, and then
transformed to the Fourier domain again. The values of the taper are:
0 if t < t1
0.5-0.5*cos(pi*(t-t1)/(t2-t1)) if t1 <= t <= t2
1 if t2 < t < t3
0.5+0.5*cos(pi*(t-t3)/(t3-t4)) if t3 <= t <= t4
0 if t > t4
Time values are given in the same unit in which the sampling interval is given
in the input time series. I.e. if sampling interval is specified as a fraction
of seconds (which is standard) then all time values passed as parameters are
also given as fractions or multiples of seconds.
# ----- END OF stfinvfourier_description_usage.txt -----
......@@ -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 -----
......@@ -322,7 +322,9 @@ The locations of the receivers may either be specified in a separate file REC\_F
\end{verbatim}}}
These receiver coordinates in REC\_FILE are shifted by REFREC[1], REFREC[2] into the x- and y-direction, respectively. This allows for completely moving the receiver spread without modifying REC\_FILE. This may be useful for the simulation of moving profiles in reflection seismics.
If READREC=0 the receiver locations must be specified in the parameter file. In this case, it is assumed that the receivers are located along a straight line. The first receiver position is defined by (XREC1, YREC1), and the last receiver position by (XREC1, YREC1). The spacing between receivers is NGEOPH grid points.
If READREC=0 the receiver locations must be specified in the parameter file. In this case, it is assumed that the receivers are located along a straight line. The first receiver position is defined by (XREC1, YREC1), and the last receiver position by (XREC1, YREC1). The spacing between receivers is NGEOPH grid points.
READREC=2 is made for a moving streamer aquisition geometry. The receiver positions for shot 1 must be specified in the parameter file (as described for READREC=0). For all other shots the receiver positions are moved according to the movement of the following source. This implementation allows a memory saving usage of moving geometries. The receivers must be in the same depth and all receivers must be located in the model for all shots.
Receivers are always located on full grid indices, i.e. a receiver that is located between two grid points will be shifted by the FD program to the closest next grid point. It is not yet possible
to output seismograms for arbitrary receiver locations since this would require a certain wavefield interpolation.
......@@ -640,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",
......@@ -652,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.
......@@ -716,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.
......@@ -736,7 +739,7 @@ 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.
......@@ -796,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.\\
......@@ -904,4 +907,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.
/*-----------------------------------------------------------------------------------------
* Copyright (C) 2005 <name of author>
*
* This file is part of DENISE.
*
* DENISE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, version 2.0 of the License only.
*
* DENISE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/* $Id: hh_el.c,v 2.3 2007/08/21 13:16:19 tbohlen Exp $*/
/*
* Model defined by flnode file. In this version a constant Q-model with Q_p=Q_s is used. It
* is defined via the relaxation frequencies and the TAU value given in the input file.
*/
#include "fd.h"
void model(float ** rho, float ** pi, float ** u, float ** taus, float ** taup, float * eta){
/*--------------------------------------------------------------------------*/
/* extern variables */
extern int NX, NY, NXG, NYG, POS[3], MYID, L;
extern float DH, DT, *FL, TAU;
extern char MFILE[STRING_SIZE];
extern char INV_MODELFILE[STRING_SIZE];
/* local variables */
float muv, piv, vp, vs, rhov, ts, tp, *pts;
int i, j, ii, jj, l;
char modfile[STRING_SIZE];
FILE *flfile;
int nodes;
char cline[256];
float *fldepth, *flrho, *flvp, *flvs;
/*---------------------------------------
Trench "Ettlinger Linie" units in m
----------------------------------------*/
/* trench=1 - elliptic shaped trench , trench=2 - V shaped trench */
int trench=2;
float width=10;
float depth=3;
float pos=NXG/2*DH;
// Parameters of trench "Ettlinger Linie" velocities in km/s and density in kg/m^3
float evs=0.0917;
float evp=0.2569;
float erho=1.5;
/* fopen FL-node file*/
flfile=fopen("model_true/flnodes.toy_example","r");
if (flfile==NULL) declare_error(" FL-file could not be opened !");
/* nodes specified in FL-node File*/
nodes=7;
/*-----------No changes below that line needed ---------------------------*/
/* 2 linear lines as border for the v-shape */
float m1=2*depth/width;
float m2=-m1;
float c1=depth*(1-(2*pos/width));
float c2=depth*(1+(2*pos/width));
fldepth=vector(1,nodes);
flrho=vector(1,nodes);
flvp=vector(1,nodes);
flvs=vector(1,nodes);
/* vector for maxwellbodies */
pts=vector(1,L);
for (l=1;l<=L;l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
}
/*read FL nodes from File*/
for (l=1;l<=nodes;l++){
fgets(cline,255,flfile);
if (cline[0]!='#'){
sscanf(cline,"%f%f%f%f",&fldepth[l], &flrho[l], &flvp[l], &flvs[l]);
}
else l=l-1;
}
if(MYID==0){
printf(" ------------------------------------------------------------------ \n\n");
printf(" Information of FL nodes: \n\n");
printf(" \t depth \t rho \t vp \t vs \n\n");
for (l=1;l<=nodes;l++){
printf(" \t %f \t %f \t %f \t %f\n\n",fldepth[l],flrho[l],flvp[l],flvs[l]);
}
printf(" ------------------------------------------------------------------ \n\n");
}
/*-----------------------------------------------------------------------*/
/*-----------------------------------------------------------------------*/
/* loop over global grid */
for (i=1;i<=NXG;i++){
for (l=1;l<nodes;l++){
if (fldepth[l]==fldepth[l+1]){
if ((i==1) && (MYID==0)){
printf("depth: %f m: double node\n",fldepth[l]);}}
else{
for (j=(int)(fldepth[l]/DH)+1;j<=(int)(fldepth[l+1]/DH);j++){
vp=0.0; vs=0.0; rhov=0.0;
vp=(DH*(j-1)-fldepth[l])*(flvp[l+1]-flvp[l])/(fldepth[l+1]-fldepth[l])+flvp[l];
vp=vp*1000.0;
vs=(DH*(j-1)-fldepth[l])*(flvs[l+1]-flvs[l])/(fldepth[l+1]-fldepth[l])+flvs[l];
vs=vs*1000.0;
rhov=(DH*(j-1)-fldepth[l])*(flrho[l+1]-flrho[l])/(fldepth[l+1]-fldepth[l])+flrho[l];
rhov=rhov*1000.0;
muv=vs;
piv=vp;
ts=TAU;
tp=TAU;
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
u[jj][ii]=muv;
rho[jj][ii]=rhov;
pi[jj][ii]=piv;
taus[jj][ii]=ts;
taup[jj][ii]=tp;
}
}
}
}
for (j=(int)(fldepth[nodes]/DH)+1;j<=NYG;j++){
vp=0.0; vs=0.0; rhov=0.0;
vp=flvp[nodes]*1000.0; vs=flvs[nodes]*1000.0; rhov=flrho[nodes]*1000.0;
muv=vs;
piv=vp;
ts=TAU;
tp=TAU;
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
u[jj][ii]=muv;
rho[jj][ii]=rhov;
pi[jj][ii]=piv;
taus[jj][ii]=ts;
taup[jj][ii]=tp;
}
}
/*Trench Ettlinger Linie*/
if (trench==1) {
for (j=1;j<=(int) depth/DH+1;j++){
if ( ( ((pos-i*DH) * (pos-i*DH))/(width/2*width/2) + ((j*DH)*(j*DH))/(depth*depth) ) <=1){
vp=evp*1000.0;
vs=evs*1000.0;
rhov=erho*1000.0;
muv=vs;
piv=vp;
ts=TAU;
tp=TAU;
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
u[jj][ii]=muv;
rho[jj][ii]=rhov;
pi[jj][ii]=piv;
taus[jj][ii]=ts;
taup[jj][ii]=tp;
}
}
}
}
/*Trench V-shape Ettlinger Linie*/
if (trench==2) {
for (j=1;j<=(int) depth/DH+1;j++){
if ( ((i*DH)*m1+c1 >=0) && j*DH <= (i*DH)*m1+c1 && j*DH <= (i*DH)*m2+c2){
vp=evp*1000.0;
vs=evs*1000.0;
rhov=erho*1000.0;
muv=vs;
piv=vp;
ts=TAU;
tp=TAU;
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
u[jj][ii]=muv;
rho[jj][ii]=rhov;
pi[jj][ii]=piv;
taus[jj][ii]=ts;
taup[jj][ii]=tp;
}
}
}
}
}
sprintf(modfile,"%s_rho_it0.bin",INV_MODELFILE);
writemod(modfile,rho,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"%s_rho_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
remove(modfile);
sprintf(modfile,"%s_vs_it0.bin",INV_MODELFILE);
writemod(modfile,u,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"%s_vs_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
remove(modfile);
sprintf(modfile,"%s_vp_it0.bin",INV_MODELFILE);
writemod(modfile,pi,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"%s_vp_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
remove(modfile);
sprintf(modfile,"%s_taus_it0.bin",INV_MODELFILE);
writemod(modfile,taus,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"%s_taus_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
remove(modfile);
sprintf(modfile,"%s_taup_it0.bin",INV_MODELFILE);
writemod(modfile,taup,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"%s_taup_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
remove(modfile);
free_vector(fldepth,1,nodes);
free_vector(flrho,1,nodes);
free_vector(flvp,1,nodes);
free_vector(flvs,1,nodes);
free_vector(pts,1,L);
}
......@@ -7,6 +7,14 @@
# compiling all libraries and IFOS
make install MODEL=../genmod/toy_example_true.c
#-------------------------------------------------------------------------------------
# Insert the trench "Ettlinger Linie" into the toy_example model
# in the file ../genmod/toy_example_true_trench.c the size and position can be adjusted
# Either an elliptic trench or a v-shaped trench can be choosen
#make install MODEL=../genmod/toy_example_true_trench.c
#---------------------------------------------------------------------------------------
# starting IFOS for forward modeling
# (depending on the MPI package you maybe have to adjust the following programme call,
# e.g. when using openMPI you do not have to use the 'lamboot' command)
......
This diff is collapsed.
......@@ -19,7 +19,7 @@ EXEC= ../bin
# LINUX with OpenMPI / IntelMPI and INTEL Compiler
# Use icc whenever possible, this will be much faster than gcc
CC=mpiicc
CC=mpicc
LFLAGS=-lm -lcseife -lstfinv -laff -lfourierxx -lfftw3 -lstdc++
CFLAGS=-O3
SFLAGS=-L./../contrib/libcseife -L./../contrib/bin
......@@ -70,6 +70,7 @@ IFOS2D= \
IFOS2D.c \
stf.c \
window_cos.c \
alloc_sections.c \
calc_mat_change_test.c \
calc_res.c \
calc_misfit.c \
......@@ -205,4 +206,4 @@ clean:
install: clean all
-include $(IFOS_OBJ:.o=.d)
\ No newline at end of file
-include $(IFOS_OBJ:.o=.d)
/*-----------------------------------------------------------------------------------------
* Copyright (C) 2016 For the list of authors, see file AUTHORS.
*
* This file is part of IFOS.
*
* IFOS is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, version 2.0 of the License only.
*
* IFOS is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/* Computation of local source coordinates
*
*
*/
#include "fd.h"
void alloc_sections(int ntr,int ns,float ***sectionvx,float ***sectionvy,float ***sectionvz,float ***sectionp,float ***sectionpnp1,float ***sectionpn,float ***sectioncurl,float ***sectiondiv,
float ***sectionpdata,float ***sectionpdiff,float ***sectionpdiffold,float ***sectionvxdata,float ***sectionvxdiff,float ***sectionvxdiffold,float ***sectionvydata,
float ***sectionvydiff,float ***sectionvydiffold,float ***sectionvzdata,float ***sectionvzdiff,float ***sectionvzdiffold)
{
extern int SEISMO, WAVETYPE;
extern FILE *FP;
switch (SEISMO){
case 1 : /* particle velocities only */
switch (WAVETYPE) {
case 1:
*sectionvx=matrix(1,ntr,1,ns);
*sectionvy=matrix(1,ntr,1,ns);
break;
case 2:
*sectionvz=matrix(1,ntr,1,ns);
break;
case 3:
*sectionvx=matrix(1,ntr,1,ns);
*sectionvy=matrix(1,ntr,1,ns);
*sectionvz=matrix(1,ntr,1,ns);
break;
}
break;
case 2 : /* pressure only */
*sectionp=matrix(1,ntr,1,ns);
*sectionpnp1=matrix(1,ntr,1,ns);
*sectionpn=matrix(1,ntr,1,ns);
break;
case 3 : /* curl and div only */
*sectioncurl=matrix(1,ntr,1,ns);
*sectiondiv=matrix(1,ntr,1,ns);
break;
case 4 : /* everything */
switch (WAVETYPE) {
case 1:
*sectionvx=matrix(1,ntr,1,ns);
*sectionvy=matrix(1,ntr,1,ns);
break;
case 2:
*sectionvz=matrix(1,ntr,1,ns);
break;
case 3:
*sectionvx=matrix(1,ntr,1,ns);
*sectionvy=matrix(1,ntr,1,ns);
*sectionvz=matrix(1,ntr,1,ns);
break;
}
*sectioncurl=matrix(1,ntr,1,ns);
*sectiondiv=matrix(1,ntr,1,ns);
*sectionp=matrix(1,ntr,1,ns);
break;
case 5 : /* everything except curl and div*/
switch (WAVETYPE) {
case 1:
*sectionvx=matrix(1,ntr,1,ns);
*sectionvy=matrix(1,ntr,1,ns);
break;
case 2:
*sectionvz=matrix(1,ntr,1,ns);
break;
case 3:
*sectionvx=matrix(1,ntr,1,ns);
*sectionvy=matrix(1,ntr,1,ns);