Skip to content

GitLab

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

Changed naming convention of individual spatial tapers for gradients, updated... 

Changed naming convention of individual spatial tapers for gradients, updated to usage with WORKFLOW
parent 4a288b60
Hide whitespace changes
Inline Side-by-side
Showing with 142 additions and 348 deletions
doc/latex/6_Parameter_Definition.tex
View file @ c56120d8
......@@ -439,7 +439,7 @@ INV_VP_ITER=0
INV_VS_ITER=0
\end{verbatim}}}
This section covers some general inversion parameters. The maximum number of iterations is defined by ITERMAX. The switch FORWARD\_ONLY controls if only the forward modeling code should be used (FORWARD\_ONLY=1), e.\,g. to calculate synthetic seismograms or a complete FWT run (FORWARD\_ONLY=0). The seismic sections of the real data need to be located in DATA\_DIR and should have the ending \_x.su.shot<shotnumber> for the x-component and so on. As noted in section \ref{model parametrizations} the gradients can be expressed for different model parameterizations. The switch PARAMETERIZATION defines which parameterization should be used, seismic velocities and density (Vp,Vs,rho, PARAMETERIZATION=1), seismic impedances (Zp,Zs,rho, PARAMETERIZATION=2) or Lam$\rm{\acute{e}}$ parameters ($\rm{\lambda,\mu,\rho}$, PARAMETERIZATION=3). Please use PARAMETERIZATION>1 with care, as current developers are only working with PARAMETERIZATION=1.
This section covers some general inversion parameters. The maximum number of iterations is defined by ITERMAX. The switch FORWARD\_ONLY controls if only the forward modeling code should be used (FORWARD\_ONLY=1), e.\,g. to calculate synthetic seismograms or a complete FWT run (FORWARD\_ONLY=0). The seismic sections of the real data need to be located in DATA\_DIR and should have the ending \_vx.su.shot<shotnumber> for the x-component and so on. As noted in section \ref{model parametrizations} the gradients can be expressed for different model parameterizations. The switch PARAMETERIZATION defines which parameterization should be used, seismic velocities and density (Vp,Vs,rho, PARAMETERIZATION=1), seismic impedances (Zp,Zs,rho, PARAMETERIZATION=2) or Lam$\rm{\acute{e}}$ parameters ($\rm{\lambda,\mu,\rho}$, PARAMETERIZATION=3). Please use PARAMETERIZATION>1 with care, as current developers are only working with PARAMETERIZATION=1.
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.
......@@ -768,9 +768,7 @@ Moreover, you can directly define a offset range that should be killed. Simply s
"FILTSIZE" : "1",
"SWS_TAPER_FILE" : "0",
"SWS_TAPER_FILE_PER_SHOT" : "0",
"TAPER_FILE_NAME" : "taper.bin",
"TAPER_FILE_NAME_U" : "taper_u.bin",
"TAPER_FILE_NAME_RHO" : "taper_rho.bin",
"TAPER_FILE_NAME" : "taper",
\end{verbatim}}}
{\color{red}{\begin{verbatim}
......@@ -785,9 +783,11 @@ SWS_TAPER_FILE_PER_SHOT=0
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.\\
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 for vp, in TAPER\_FILE\_NAME\_U for vs and in TAPER\_FILE\_NAME\_RHO for the density.\\
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.\\
It is also possible to apply externally defined taper files to the gradients of the single shots before summation of these gradients. This can be done by setting SWS\_TAPER\_FILE\_PER\_SHOT to one. IFOS2D expects the tapers in TAPER\_FILE\_NAME.shot<shotnumber> for the vp gradients, in TAPER\_FILE\_NAME\_U.shot<shotnumber> for the vs gradients and in TAPER\_FILE\_NAME\_RHO.shot<shotnumber> for the density gradients.
It is also possible to apply externally defined taper files to the gradients of the single shots before summation of these gradients. This can be done by setting SWS\_TAPER\_FILE\_PER\_SHOT to one. IFOS2D expects the tapers in TAPER\_FILE\_NAME.shot<shotnumber>.vp for the vp gradients, in TAPER\_FILE\_NAME.shot<shotnumber>.vs for the vs gradients and in TAPER\_FILE\_NAME.shot<shotnumber>.rho for the density gradients.\\
If you use the option "Workflow" (section \ref{sec:workflow}) it is possible to define separate files for each workflow stage. They have to be named TAPER\_FILE\_NAME(.shot<shotnumber>)\_<workflowstage>.vp etc. If you don't provide separate files for some or all stages the standard file (TAPER\_FILE\_NAME(.shot<shotnumber>).vp etc.) is used. This means a standard file should exist in any case.
\subsection{Spatial filtering of the gradients}
{\color{blue}{\begin{verbatim}
......
genmod/toy_example_ac_start.c
View file @ c56120d8
/*-----------------------------------------------------------------------------------------
* Copyright (C) 2005 <name of author>
* Copyright (C) 2016 For the list of authors, see file AUTHORS.
*
* This file is part of DENISE.
* This file is part of IFOS.
*
* DENISE is free software: you can redistribute it and/or modify
* 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.
*
* DENISE is distributed in the hope that it will be useful,
* 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 DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
* along with IFOS. 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.
*/
------------------------------------------------------------------------*/
#include "fd.h"
#include "../src/fd.h"
......@@ -31,7 +30,7 @@ void model_acoustic(float **rho, float **pi){
extern float DH, DT;
extern char MFILE[STRING_SIZE];
extern char INV_MODELFILE[STRING_SIZE];
extern char TAPER_FILE_NAME[STRING_SIZE], TAPER_FILE_NAME_RHO[STRING_SIZE];
extern char TAPER_FILE_NAME[STRING_SIZE];
extern int SWS_TAPER_FILE;
/* local variables */
......@@ -261,7 +260,7 @@ void model_acoustic(float **rho, float **pi){
free_vector(flvp,1,nodes);
sprintf(modfile,"%s_rho_it0.bin",INV_MODELFILE);
sprintf(modfile,"%s_rho_it0.bin",INV_MODELFILE);
writemod(modfile,rho,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
......@@ -278,20 +277,20 @@ void model_acoustic(float **rho, float **pi){
remove(modfile);
if(SWS_TAPER_FILE){
sprintf(modfile,"%s",TAPER_FILE_NAME);
sprintf(modfile,"%s.vp",TAPER_FILE_NAME);
writemod(modfile,taper,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"%s.%i%i",TAPER_FILE_NAME,POS[1],POS[2]);
sprintf(modfile,"%s.vp.%i%i",TAPER_FILE_NAME,POS[1],POS[2]);
remove(modfile);
sprintf(modfile,"%s",TAPER_FILE_NAME_RHO);
sprintf(modfile,"%s.rho",TAPER_FILE_NAME);
writemod(modfile,taper,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"%s.%i%i",TAPER_FILE_NAME_RHO,POS[1],POS[2]);
sprintf(modfile,"%s.rho.%i%i",TAPER_FILE_NAME,POS[1],POS[2]);
remove(modfile);
}
......
par/in_and_out/IFOS2D_INV_all_parameters.json
View file @ c56120d8
......@@ -201,9 +201,7 @@
"FILTSIZE" : "1",
"SWS_TAPER_FILE" : "0",
"SWS_TAPER_FILE_PER_SHOT" : "0",
"TAPER_FILE_NAME" : "taper.bin",
"TAPER_FILE_NAME_U" : "taper_u.bin",
"TAPER_FILE_NAME_RHO" : "taper_rho.bin",
"TAPER_FILE_NAME" : "taper",
"Definition of spatial filtering of the gradients" : "comment",
"SPATFILTER" : "0",
......
par/in_and_out/toy_example/toy_example_ac_INV.json
View file @ c56120d8
......@@ -171,9 +171,7 @@
"SRTRADIUS" : "3.0",
"FILTSIZE" : "1",
"SWS_TAPER_FILE" : "1",
"TAPER_FILE_NAME" : "taper_te_ac.bin",
"TAPER_FILE_NAME_U" : "taper_te_ac_u.bin",
"TAPER_FILE_NAME_RHO" : "taper_te_ac_rho.bin",
"TAPER_FILE_NAME" : "taper_te_ac",
"Upper and lower limits for model parameters" : "comment",
"VPUPPERLIM" : "5000",
......
src/exchange_par.c
View file @ c56120d8
......@@ -16,11 +16,8 @@
* along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/**/
/*------------------------------------------------------------------------
* Exchange FD-Parameters between PEs
* last update 29/06/2002
*
* ----------------------------------------------------------------------*/
#include "fd.h"
......@@ -52,7 +49,7 @@ void exchange_par(void){
extern int SWS_TAPER_GRAD_VERT, SWS_TAPER_GRAD_HOR, SWS_TAPER_GRAD_SOURCES, SWS_TAPER_CIRCULAR_PER_SHOT, SRTSHAPE, FILTSIZE;
extern int SWS_TAPER_FILE, SWS_TAPER_FILE_PER_SHOT;
extern float SRTRADIUS;
extern char TAPER_FILE_NAME[STRING_SIZE], TAPER_FILE_NAME_U[STRING_SIZE], TAPER_FILE_NAME_RHO[STRING_SIZE];
extern char TAPER_FILE_NAME[STRING_SIZE];
extern int SPATFILTER, SPAT_FILT_SIZE, SPAT_FILT_1, SPAT_FILT_ITER;
extern int INV_RHO_ITER, INV_VP_ITER, INV_VS_ITER;
extern int MIN_ITER;
......@@ -409,8 +406,6 @@ void exchange_par(void){
MPI_Bcast(&PICKS_FILE,STRING_SIZE,MPI_CHAR,0,MPI_COMM_WORLD);
MPI_Bcast(&TAPER_FILE_NAME,STRING_SIZE,MPI_CHAR,0,MPI_COMM_WORLD);
MPI_Bcast(&TAPER_FILE_NAME_U,STRING_SIZE,MPI_CHAR,0,MPI_COMM_WORLD);
MPI_Bcast(&TAPER_FILE_NAME_RHO,STRING_SIZE,MPI_CHAR,0,MPI_COMM_WORLD);
MPI_Barrier(MPI_COMM_WORLD);
......
src/globvar.h
View file @ c56120d8
......@@ -51,7 +51,7 @@ int TESTSHOT_START, TESTSHOT_END, TESTSHOT_INCR, NO_OF_TESTSHOTS;
int SWS_TAPER_GRAD_VERT, SWS_TAPER_GRAD_HOR, SWS_TAPER_GRAD_SOURCES, SWS_TAPER_CIRCULAR_PER_SHOT, SRTSHAPE, FILTSIZE;
int SWS_TAPER_FILE, SWS_TAPER_FILE_PER_SHOT;
float SRTRADIUS;
char TAPER_FILE_NAME[STRING_SIZE], TAPER_FILE_NAME_U[STRING_SIZE], TAPER_FILE_NAME_RHO[STRING_SIZE];
char TAPER_FILE_NAME[STRING_SIZE];
int SPATFILTER, SPAT_FILT_SIZE, SPAT_FILT_1, SPAT_FILT_ITER;
......
src/read_par_json.c
View file @ c56120d8
......@@ -51,7 +51,7 @@ void read_par_json(FILE *fp, char *fileinp){
extern int SWS_TAPER_GRAD_VERT, SWS_TAPER_GRAD_HOR, SWS_TAPER_GRAD_SOURCES, SWS_TAPER_CIRCULAR_PER_SHOT, SRTSHAPE, FILTSIZE;
extern int SWS_TAPER_FILE, SWS_TAPER_FILE_PER_SHOT;
extern float SRTRADIUS;
extern char TAPER_FILE_NAME[STRING_SIZE], TAPER_FILE_NAME_U[STRING_SIZE], TAPER_FILE_NAME_RHO[STRING_SIZE];
extern char TAPER_FILE_NAME[STRING_SIZE];
extern int SPATFILTER, SPAT_FILT_SIZE, SPAT_FILT_1, SPAT_FILT_ITER;
extern int INV_RHO_ITER, INV_VS_ITER, INV_VP_ITER;
extern char INV_MODELFILE[STRING_SIZE];
......@@ -563,26 +563,20 @@ void read_par_json(FILE *fp, char *fileinp){
}
if (get_int_from_objectlist("SWS_TAPER_FILE",number_readobjects,&SWS_TAPER_FILE,varname_list, value_list)){
SWS_TAPER_FILE=0;
fprintf(fp,"Variable SWS_TAPER_FILE is set to default value %d.\n",SWS_TAPER_FILE);}
fprintf(fp,"Variable SWS_TAPER_FILE is set to default value %d.\n",SWS_TAPER_FILE);
}
if (SWS_TAPER_FILE==1){
if (get_string_from_objectlist("TAPER_FILE_NAME",number_readobjects,TAPER_FILE_NAME,varname_list, value_list))
err("Variable TAPER_FILE_NAME could not be retrieved from the json input file!");
if(!ACOUSTIC){
if (get_string_from_objectlist("TAPER_FILE_NAME_U",number_readobjects,TAPER_FILE_NAME_U,varname_list, value_list))
err("Variable TAPER_FILE_NAME_U could not be retrieved from the json input file!");
}
if (get_string_from_objectlist("TAPER_FILE_NAME_RHO",number_readobjects,TAPER_FILE_NAME_RHO,varname_list, value_list))
err("Variable TAPER_FILE_NAME_RHO could not be retrieved from the json input file!");}
}
if (get_int_from_objectlist("SWS_TAPER_FILE_PER_SHOT",number_readobjects,&SWS_TAPER_FILE_PER_SHOT,varname_list, value_list)){
SWS_TAPER_FILE_PER_SHOT=0;
fprintf(fp,"Variable SWS_TAPER_FILE_PER_SHOT is set to default value %d.\n",SWS_TAPER_FILE_PER_SHOT);}
fprintf(fp,"Variable SWS_TAPER_FILE_PER_SHOT is set to default value %d.\n",SWS_TAPER_FILE_PER_SHOT);
}
if (SWS_TAPER_FILE_PER_SHOT==1){
if (get_string_from_objectlist("TAPER_FILE_NAME",number_readobjects,TAPER_FILE_NAME,varname_list, value_list))
err("Variable TAPER_FILE_NAME could not be retrieved from the json input file!");
if (get_string_from_objectlist("TAPER_FILE_NAME_U",number_readobjects,TAPER_FILE_NAME_U,varname_list, value_list))
err("Variable TAPER_FILE_NAME_U could not be retrieved from the json input file!");
if (get_string_from_objectlist("TAPER_FILE_NAME_RHO",number_readobjects,TAPER_FILE_NAME_RHO,varname_list, value_list))
err("Variable TAPER_FILE_NAME_RHO could not be retrieved from the json input file!");}
}
......
src/taper_grad.c
View file @ c56120d8
......@@ -16,13 +16,13 @@
* along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/*
/*------------------------------------------------------------------------
* taper gradient with a gaussian frame to damp inversion artefacts near the sources and receivers
sws == 1 vertical taper (for tomography geometry)
sws == 2 horizontal taper (for reflection geometry)
sws == 3 local radial taper at the source and receiver positions
sws == 4 custom taper
*/
------------------------------------------------------------------------*/
#include "fd.h"
......@@ -31,11 +31,12 @@ void taper_grad(float ** waveconv,float ** taper_coeff, float **srcpos, int nsho
/* extern variables */
extern float DH;
extern float DH;
extern int FREE_SURF, NX, NY, NXG, NYG;
extern int NPROCX, NPROCY, MYID, POS[3];
extern FILE *FP;
extern char TAPER_FILE_NAME[STRING_SIZE], TAPER_FILE_NAME_U[STRING_SIZE], TAPER_FILE_NAME_RHO[STRING_SIZE];
extern char TAPER_FILE_NAME[STRING_SIZE];
extern int USE_WORKFLOW, WORKFLOW_STAGE;
/* local variables */
int i, j, h, ifw, ii, jj, n, xb, yb, xe, ye, taperlength,taperlength2, VTON, SRTON;
......@@ -53,10 +54,6 @@ void taper_grad(float ** waveconv,float ** taper_coeff, float **srcpos, int nsho
extern int VERBOSE;
FILE *fp_taper;
/*SRTSHAPE=2;
SRTRADIUS=25.0;
filtsize=2;*/
/* =============== */
/* Vertical taper */
/* =============== */
......@@ -68,9 +65,6 @@ void taper_grad(float ** waveconv,float ** taper_coeff, float **srcpos, int nsho
taperlength2=GRADT4-GRADT3;
ifw = GRADT4-GRADT1;
/*printf("%d \t %d \t %d \t %d \n",GRADT1, GRADT2, GRADT3, GRADT4);
printf("%d \t %d \t %d \n",taperlength, taperlength2,ifw);*/
if (MYID==0&&VERBOSE)
{
fprintf(FP,"\n **Message from taper_grid (printed by PE %d):\n",MYID);
......@@ -136,11 +130,7 @@ void taper_grad(float ** waveconv,float ** taper_coeff, float **srcpos, int nsho
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii = i-POS[1]*NX;
jj = j-POS[2]*NY;
/*waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
+ 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
+ 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;*/
waveconvtmp[jj][ii] = 0.0;
}
}
......@@ -158,11 +148,7 @@ void taper_grad(float ** waveconv,float ** taper_coeff, float **srcpos, int nsho
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii = i-POS[1]*NX;
jj = j-POS[2]*NY;
/*waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
+ 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
+ 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;*/
waveconvtmp[jj][ii] = 0.0;
}
......@@ -232,10 +218,6 @@ void taper_grad(float ** waveconv,float ** taper_coeff, float **srcpos, int nsho
if(j>GRADT2){
/*grad_tap=((float)(j*DH))*((float)(j*DH))*((float)(j*DH));*/
/* grad_tap=((float)(j*DH))*((float)(j*DH));*/
/*grad_tap=((float)(j*DH));*/ /*Version of Daniel */
grad_tap=1.0;
}
......@@ -264,48 +246,6 @@ void taper_grad(float ** waveconv,float ** taper_coeff, float **srcpos, int nsho
}
/* apply filter at shot and receiver points */
/*for (n=1;n<=nshots;n++)
{
i = iround(srcpos[1][n]/DH);
j = iround(srcpos[2][n]/DH);
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii = i-POS[1]*NX;
jj = j-POS[2]*NY;
waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
+ 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
+ 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;
}
}
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
waveconv[j][i] = waveconvtmp[j][i];
}
}
for (n=1;n<=ntr;n++)
{
i = recpos[1][n];
j = recpos[2][n];
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii = i-POS[1]*NX;
jj = j-POS[2]*NY;
waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
+ 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
+ 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;
}
}
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
waveconv[j][i] = waveconvtmp[j][i];
}
}*/
sprintf(modfile,"taper_coeff_hor.bin");
writemod(modfile,taper_coeff,3);
MPI_Barrier(MPI_COMM_WORLD);
......@@ -399,9 +339,6 @@ void taper_grad(float ** waveconv,float ** taper_coeff, float **srcpos, int nsho
}
}
/* for (iy=1;iy<=NYG;iy++)
for (ix=1;ix<=NXG;ix++) msum[iy][ix] += mm[iy][ix];
*/
for (iy=1;iy<=NYG;iy++)
for (ix=1;ix<=NXG;ix++)
if (msum[iy][ix] > mm[iy][ix])
......@@ -409,33 +346,6 @@ void taper_grad(float ** waveconv,float ** taper_coeff, float **srcpos, int nsho
}
/***********************/
/* loop over receivers */
/*for (n=1;n<=ntr;n++) {
for (iy=1;iy<=NYG;iy++)
for (ix=1;ix<=NXG;ix++) mm[iy][ix] = 1.0;
i = recpos[1][n];
j = recpos[2][n];
for (iy=1;iy<=srctaper_gridpt;iy++) {
for (ix=1;ix<=srctaper_gridpt;ix++) {
xx = i + ix - ijc;
yy = j + iy - ijc;
if ((xx<1) || (xx>NXG) || (yy<1) || (yy>NYG)) continue;
mm[yy][xx] = m[iy][ix];
}
}*/
/* for (iy=1;iy<=NYG;iy++) Die kommenden zwei Zeilen wurden von Daniel auskommentiert.
for (ix=1;ix<=NXG;ix++) msum[iy][ix] += mm[iy][ix];
*/
/* for (iy=1;iy<=NYG;iy++)
for (ix=1;ix<=NXG;ix++)
if (msum[iy][ix] > mm[iy][ix])
msum[iy][ix] = mm[iy][ix];
}*/
minm = minimum_m(msum,NXG,NYG);
for (iy=1;iy<=NYG;iy++)
for (ix=1;ix<=NXG;ix++) msum[iy][ix] -= minm;
......@@ -448,8 +358,6 @@ void taper_grad(float ** waveconv,float ** taper_coeff, float **srcpos, int nsho
if ((POS[1]==((ix-1)/NX)) && (POS[2]==((iy-1)/NY))){
ii = ix-POS[1]*NX;
jj = iy-POS[2]*NY;
/*Diese Zeile wurde von Daniel auskommentiert: taper_coeff[jj][ii] = taper_coeff_glob[iy][ix] * ((float)(iy*DH));*/
/*taper_coeff[jj][ii] = ((float)(iy*DH)) * ((float)(iy*DH)) * ((float)(iy*DH));*/
taper_coeff[jj][ii]=msum[iy][ix]/maxm;
}
}
......@@ -489,29 +397,6 @@ void taper_grad(float ** waveconv,float ** taper_coeff, float **srcpos, int nsho
}
/*for (n=1;n<=ntr;n++)
{
i1 = recpos[1][n];
j1 = recpos[2][n];
for (i=i1-FILTSIZE;i<=i1+FILTSIZE;i++){
for (j=j1-FILTSIZE;j<=j1+FILTSIZE;j++){
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii = i-POS[1]*NX;
jj = j-POS[2]*NY; */
/* Die kommenden 4 Zeilen wurden von Daniel auskommentiert. waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
+ 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
+ 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;*/
/* waveconvtmp[jj][ii] = 0.0;
}
}
}
}*/
/* apply taper on local gradient */
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
......@@ -527,10 +412,6 @@ void taper_grad(float ** waveconv,float ** taper_coeff, float **srcpos, int nsho
free_matrix(taper_coeff_glob,1,NYG,1,NXG);
free_matrix(waveconvtmp,0,NX+1,0,NY+1);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"taper_coeff_sources.bin");
......@@ -546,29 +427,55 @@ void taper_grad(float ** waveconv,float ** taper_coeff, float **srcpos, int nsho
if((sws>=4)&&(sws<=6)){
if (MYID==0&&VERBOSE)
{
fprintf(FP,"\n **Message from taper_grid (printed by PE %d):\n",MYID);
fprintf(FP," Coefficients for gradient taper are now calculated.\n");
}
if (MYID==0&&VERBOSE)
{
fprintf(FP,"\n **Message from taper_grad (printed by PE %d):\n",MYID);
fprintf(FP," Coefficients for gradient tapers are now read in.\n");
}
if(sws==4){
sprintf(taper_file,"%s",TAPER_FILE_NAME);
fp_taper=fopen(taper_file,"r");
}
if(USE_WORKFLOW){
sprintf(taper_file,"%s_%i.vp",TAPER_FILE_NAME,WORKFLOW_STAGE);
fp_taper=fopen(taper_file,"r");
if(fp_taper==NULL){
sprintf(taper_file,"%s.vp",TAPER_FILE_NAME);
fp_taper=fopen(taper_file,"r");
}
}else{
sprintf(taper_file,"%s.vp",TAPER_FILE_NAME);
fp_taper=fopen(taper_file,"r");
}
}
if(sws==5){
sprintf(taper_file,"%s",TAPER_FILE_NAME_U);
fp_taper=fopen(taper_file,"r");
}
if(USE_WORKFLOW){
sprintf(taper_file,"%s_%i.vs",TAPER_FILE_NAME,WORKFLOW_STAGE);
fp_taper=fopen(taper_file,"r");
if(fp_taper==NULL){
sprintf(taper_file,"%s.vs",TAPER_FILE_NAME);
fp_taper=fopen(taper_file,"r");
}
}else{
sprintf(taper_file,"%s.vs",TAPER_FILE_NAME);
fp_taper=fopen(taper_file,"r");
}
}
if(sws==6){
sprintf(taper_file,"%s",TAPER_FILE_NAME_RHO);
fp_taper=fopen(taper_file,"r");
if(USE_WORKFLOW){
sprintf(taper_file,"%s_%i.rho",TAPER_FILE_NAME,WORKFLOW_STAGE);
fp_taper=fopen(taper_file,"r");
if(fp_taper==NULL){
sprintf(taper_file,"%s.rho",TAPER_FILE_NAME);
fp_taper=fopen(taper_file,"r");
}
}else{
sprintf(taper_file,"%s.rho",TAPER_FILE_NAME);
fp_taper=fopen(taper_file,"r");
}
}
if(fp_taper==NULL) {
err("Taper file could not opened");
}
err("Taper file could not be opened");
}
/* loop over global grid */
for (i=1;i<=NXG;i++){
......@@ -604,7 +511,4 @@ void taper_grad(float ** waveconv,float ** taper_coeff, float **srcpos, int nsho
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
}
}
}
\ No newline at end of file
src/taper_grad_shot.c
View file @ c56120d8
<
......@@ -16,13 +16,13 @@
* along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/*
/*------------------------------------------------------------------------
* taper gradient with a gaussian frame to damp inversion artefacts near the sources and receivers
sws == 1 vertical taper (for tomography geometry)
sws == 2 horizontal taper (for reflection geometry)
sws == 3 local radial taper at the source and receiver positions
sws == 4
*/
------------------------------------------------------------------------*/
#include "fd.h"
......@@ -35,8 +35,10 @@ void taper_grad_shot(float ** waveconv,float ** taper_coeff, float **srcpos, int
extern int FREE_SURF, NX, NY, NXG, NYG;
extern int NPROCX, NPROCY, MYID, POS[3];
extern FILE *FP;
extern char TAPER_FILE_NAME[STRING_SIZE], TAPER_FILE_NAME_U[STRING_SIZE], TAPER_FILE_NAME_RHO[STRING_SIZE];
extern char TAPER_FILE_NAME[STRING_SIZE];
extern int VERBOSE;
extern int USE_WORKFLOW, WORKFLOW_STAGE;
/* local variables */
int i, j, h, ifw, ii, jj, n, xb, yb, xe, ye, taperlength,taperlength2, VTON, SRTON;
int ijc, iy, ix, iii, jjj, xx, yy, srctaper_gridpt, i1, j1;
......@@ -52,8 +54,6 @@ void taper_grad_shot(float ** waveconv,float ** taper_coeff, float **srcpos, int
FILE *fp_taper;
if(sws==1){
/* =================================== */
......@@ -138,42 +138,11 @@ void taper_grad_shot(float ** waveconv,float ** taper_coeff, float **srcpos, int
}
}
/* for (iy=1;iy<=NYG;iy++)
for (ix=1;ix<=NXG;ix++) msum[iy][ix] += mm[iy][ix];
*/
for (iy=1;iy<=NYG;iy++)
for (ix=1;ix<=NXG;ix++)
if (msum[iy][ix] > mm[iy][ix])
msum[iy][ix] = mm[iy][ix];
/* }*/
/***********************/
/* loop over receivers */
/*for (n=1;n<=ntr;n++) {
for (iy=1;iy<=NYG;iy++)
for (ix=1;ix<=NXG;ix++) mm[iy][ix] = 1.0;
i = recpos[1][n];
j = recpos[2][n];
for (iy=1;iy<=srctaper_gridpt;iy++) {
for (ix=1;ix<=srctaper_gridpt;ix++) {
xx = i + ix - ijc;
yy = j + iy - ijc;
if ((xx<1) || (xx>NXG) || (yy<1) || (yy>NYG)) continue;
mm[yy][xx] = m[iy][ix];
}
}*/
/* for (iy=1;iy<=NYG;iy++) Die kommenden zwei Zeilen wurden von Daniel auskommentiert.
for (ix=1;ix<=NXG;ix++) msum[iy][ix] += mm[iy][ix];
*/
/* for (iy=1;iy<=NYG;iy++)
for (ix=1;ix<=NXG;ix++)
if (msum[iy][ix] > mm[iy][ix])
msum[iy][ix] = mm[iy][ix];
}*/
minm = minimum_m(msum,NXG,NYG);
for (iy=1;iy<=NYG;iy++)
......@@ -203,26 +172,15 @@ void taper_grad_shot(float ** waveconv,float ** taper_coeff, float **srcpos, int
/* apply filter at shot and receiver points */
/*for (n=1;n<=nshots;n++)
{*/
n=ishot;
i1 = iround(srcpos[1][n]/DH);
j1 = iround(srcpos[2][n]/DH);
/*fprintf(FP,"\n Shot %d (printed by PE %d):\n",n,MYID);
fprintf(FP,"\n i1: %d (printed by PE %d):\n",i1,MYID);
fprintf(FP,"\n j1: %d (printed by PE %d):\n",j1,MYID);*/
for (i=i1-FILTSIZE;i<=i1+FILTSIZE;i++){
for (j=j1-FILTSIZE;j<=j1+FILTSIZE;j++){
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii = i-POS[1]*NX;
jj = j-POS[2]*NY;
/*printf("\n ii: %d (printed by PE %d):\n",ii,MYID);
printf("\n jj: %d (printed by PE %d):\n",jj,MYID);*/
/*waveconvtmp[jj][ii] = 1*waveconv[jj][ii]