...
 
Commits (1)
......@@ -375,6 +375,29 @@ Small $Q$ values ($Q<50$) may lead to significant amplitude decay and velocity d
The frequency dependence of attenuation, i.e. $Q$ and phase velocity as a function of frequency, may be calculated using the Matlab functions in the directory mfiles. The Matlab script /mfiles/qplot.m can be used to plot $Q(\omega)$ for different values of L, $f_l$ and $\tau$. The m-file qapprox.m in the same directory finds optimal values for L, $f_l$ and $\tau$ that fit a desired function $Q(\omega)=const$ in a least-squares sense.
\section{Anisotropic modeling and inversion}
\label{anisotropy}
{\color{blue}{\begin{verbatim}
"Anisotropy" : "comment",
"VTI" : "0",
\end{verbatim}}}
With VTI=1, the forward modeling and the inversion is calculated for a vertically transversely isotropic medium. For the forward modeling with WAVETYPE=1, additional models of the Thomsen parameters $\varepsilon$ and $\delta$ are needed with the file name expansions ``.epsilon'' and ``.delta''. For WAVETYPE=2, the input models are density, vertical S-wave velocity and horizontal S-wave velocity (``.rho'', ``.vs'', ``.vshor''). The inversion is only implemented for SH waves (WAVETYPE=2) and should only be used with VELOCITY=1. Additional parameters that are used equivalent to those existing for the isotropic case are
{\color{blue}{\begin{verbatim}
"INV_VSHOR_ITER" : "0",
"VSHORUPPERLIM" : "5000",
"VSHORLOWERLIM" : "0",
\end{verbatim}}}
\noindent
Additionally, with
{\color{blue}{\begin{verbatim}
"GAMMA_VTI" : "0",
\end{verbatim}}}
a maximal absolute value of the Thomsen parameter $\gamma$ can be defined to constrain updates of the velocities (constraint is not active for GAMMA\_VTI=0).
\section{Wavefield snapshots}
{\color{blue}{\begin{verbatim}
"Snapshots" : "comment",
......
/*-----------------------------------------------------------------------------------------
* 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>.
-----------------------------------------------------------------------------------------*/
/*
* Model 1D linear gradient
* Vertically transversely isotropic model
*/
#include "fd.h"
void model_elastic_vti(float ** rho, float ** pi, float ** u, float ** delta, float ** epsilon, float ** vshor){
/*--------------------------------------------------------------------------*/
/* extern variables */
extern int NX, NY, NXG, NYG, POS[3], L, MYID;
extern char MFILE[STRING_SIZE];
extern char INV_MODELFILE[STRING_SIZE];
extern float DH;
/* local variables */
float vp, vs, rhov, grad1, grad2, grad3, grad4, grad5, grad6, y;
float c11v, c13v, c33v, c55v, c66v, deltav, epsilonv, vshorv;
int i, j, ii, jj;
char modfile[STRING_SIZE];
float ** pwavemod=NULL, ** swavemod=NULL;
/* parameter defintion:
* vp = vertical P-wave-velocity, vs = vertical S-wave-velocity, vshor = horizontal S-wave-velocity
* delta, epsilon = Thomsen parameters */
/* parameters for layer 1 */
const float vp1=500.0, vs1=300.0, rho1=1800.0, h=15.0;
const float delta1=0.1, epsilon1=0.2, vshor1=550;
/* parameters for layer 2 due to calculation of grad1, grad2 and grad3*/
const float vp2=1200.0, vs2=700.0, rho2=2000.0;
const float delta2=0.1, epsilon2=0.2, vshor2=1500;
/*-----------------------------------------------------------------------*/
y=h/DH;
if(y==NYG) declare_error(" \n y is equal NYG !! see src/model_grad.c \n ");
grad1=(vp2-vp1)/y;
grad2=(vs2-vs1)/y;
grad3=(rho2-rho1)/y;
grad4=(delta2-delta1)/y;
grad5=(epsilon2-epsilon1)/y;
grad6=(vshor2-vshor1)/y;
/* loop over global grid */
for (i=1;i<=NXG;i++){
for (j=1;j<=NYG;j++){
if(j<=y){
vp=vp1+(j*grad1);
vs=vs1+(j*grad2);
rhov=rho1+(j*grad3);
deltav=delta1+(j*grad4);
epsilonv=epsilon1+(j*grad5);
vshorv=vshor1+(j*grad6);
}
else{
vp=vp2;
vs=vs2;
rhov=rho2;
deltav=delta2;
epsilonv=epsilon2;
vshorv=vshor2;
}
/* 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]=vs;
pi[jj][ii]=vp;
rho[jj][ii]=rhov;
delta[jj][ii]=deltav;
epsilon[jj][ii]=epsilonv;
vshor[jj][ii]=vshorv;
}
}
}
sprintf(modfile,"%s_rho_it_0.bin",INV_MODELFILE);
writemod(modfile,rho,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
sprintf(modfile,"%s_vs_it_0.bin",INV_MODELFILE);
writemod(modfile,u,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
sprintf(modfile,"%s_vp_it_0.bin",INV_MODELFILE);
writemod(modfile,pi,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
sprintf(modfile,"%s_vshor_it_0.bin",INV_MODELFILE);
writemod(modfile,vshor,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
}
......@@ -62,6 +62,9 @@
"NDT" : "1",
"SEIS_FORMAT" : "1",
"SEIS_FILE" : "su/IFOS",
"Anisotropy" : "comment",
"VTI" : "0",
"Q-approximation" : "comment",
"L" : "0",
......
......@@ -62,6 +62,9 @@
"NDT" : "1",
"SEIS_FORMAT" : "1",
"SEIS_FILE" : "su/IFOS",
"Anisotropy" : "comment",
"VTI" : "0",
"Q-approximation" : "comment",
"L" : "0",
......@@ -222,6 +225,8 @@
"VSLOWERLIM" : "0",
"RHOUPPERLIM" : "5000",
"RHOLOWERLIM" : "0",
"VSHORUPPERLIM" : "5000",
"VSHORLOWERLIM" : "0",
"Limited update of model parameters in reference to the starting model" : "comment",
"S" : "0",
......@@ -232,6 +237,9 @@
"Minimum Vp/Vs-ratio" : "comment",
"VP_VS_RATIO" : "0.0",
"Maximal gamma-value (Thomsen parameter for VTI)" : "comment",
"GAMMA_VTI" : "0",
"Definition of smoothing the models vp and vs" : "comment",
"MODEL_FILTER" : "0",
"FILT_SIZE" : "5",
......
This diff is collapsed.
......@@ -25,14 +25,14 @@
void lbfgs_reset(int iter, int N_LBFGS, int NPAR_LBFGS,float ** s_LBFGS1, float ** y_LBFGS1, float * rho_LBFGS1);
void lbfgs_core(int iteration, int N_LBFGS, int NPAR_LBFGS,float ** s_LBFGS, float ** y_LBFGS, float * rho_LBFGS,float *q_LBFGS,float *alpha_LBFGS,float *r_LBFGS);
void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float rho_avg,float Vp_avg, float *rho_LBFGS, float **s_LBFGS, float **y_LBFGS,int N_LBFGS,int NPAR_LBFGS, int iteration, int *LBFGS_iter_start){
void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp, float **grad_vshor, float Vs_avg,float rho_avg,float Vp_avg, float Vshor_avg, float *rho_LBFGS, float **s_LBFGS, float **y_LBFGS,int N_LBFGS,int NPAR_LBFGS, int iteration, int *LBFGS_iter_start){
/* global */
extern int NX,NY,MYID;
extern FILE *FP;
extern int POS[3];
extern char JACOBIAN[STRING_SIZE];
extern int WAVETYPE;
extern int WAVETYPE, VTI;
extern int ACOUSTIC;
extern float LBFGS_SCALE_GRADIENTS;
......@@ -61,6 +61,11 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
write_matrix_disk(grad_vp, jac);
}
if(WAVETYPE==2 && VTI) {
sprintf(jac,"%s_grad1_vshor_it%d",JACOBIAN,iteration);
write_matrix_disk(grad_vshor, jac);
}
/*---------------------*/
/* Experimental */
/*---------------------*/
......@@ -76,7 +81,8 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
for (j=1;j<=NY;j++){
if(!ACOUSTIC) grad_vs[j][i]=grad_vs[j][i]*LBFGS_SCALE_GRADIENTS;
if(NPAR_LBFGS>1) grad_rho[j][i]=grad_rho[j][i]*LBFGS_SCALE_GRADIENTS;
if(NPAR_LBFGS>2) grad_vp[j][i]=grad_vp[j][i]*LBFGS_SCALE_GRADIENTS;
if(NPAR_LBFGS>2 && WAVETYPE!=2) grad_vp[j][i]=grad_vp[j][i]*LBFGS_SCALE_GRADIENTS;
if(NPAR_LBFGS>2 && VTI) grad_vshor[j][i]=grad_vshor[j][i]*LBFGS_SCALE_GRADIENTS;
}
}
}
......@@ -94,7 +100,10 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
l++;
if(!ACOUSTIC) y_LBFGS[w][l]=-grad_vs[j][i]*Vs_avg; /* VS */
if(NPAR_LBFGS>1) y_LBFGS[w][l+NX*NY]=-grad_rho[j][i]*rho_avg; /* RHO */
if(NPAR_LBFGS>2) y_LBFGS[w][l+2*NX*NY]=-grad_vp[j][i]*Vp_avg; /* VP */
if(NPAR_LBFGS>2){
if(WAVETYPE!=2) y_LBFGS[w][l+2*NX*NY]=-grad_vp[j][i]*Vp_avg; /* VP */
else if(VTI) y_LBFGS[w][l+2*NX*NY]=-grad_vshor[j][i]*Vshor_avg; /* VS HORIZONTAL (VTI) */
}
}
}
}
......@@ -141,8 +150,14 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
/* VP */
if(NPAR_LBFGS>2) {
y_LBFGS[w][l+2*NY*NX]+=grad_vp[j][i]*Vp_avg; /* add grad(i) to build grad(i)-grad(i-1) */
q_LBFGS[l+2*NY*NX]=grad_vp[j][i]*Vp_avg; /* Normalisation */
if (WAVETYPE!=2) {
y_LBFGS[w][l+2*NY*NX]+=grad_vp[j][i]*Vp_avg; /* add grad(i) to build grad(i)-grad(i-1) */
q_LBFGS[l+2*NY*NX]=grad_vp[j][i]*Vp_avg; /* Normalisation */
}
else if (VTI) {
y_LBFGS[w][l+2*NY*NX]+=grad_vshor[j][i]*Vshor_avg; /* add grad(i) to build grad(i)-grad(i-1) */
q_LBFGS[l+2*NY*NX]=grad_vshor[j][i]*Vshor_avg; /* Normalisation */
}
}
/* Debugging */
......@@ -194,8 +209,14 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
/* VP */
if(NPAR_LBFGS>2) {
y_LBFGS[w][l+2*NY*NX]=-grad_vp[j][i]*Vp_avg; /* add -grad(i-1) to build grad(i)-grad(i-1) */
grad_vp[j][i]=r_LBFGS[l+2*NY*NX]*Vp_avg; /* Denormalization */
if (WAVETYPE!=2) {
y_LBFGS[w][l+2*NY*NX]=-grad_vp[j][i]*Vp_avg; /* add -grad(i-1) to build grad(i)-grad(i-1) */
grad_vp[j][i]=r_LBFGS[l+2*NY*NX]*Vp_avg; /* Denormalization */
}
else if (VTI) {
y_LBFGS[w][l+2*NY*NX]=-grad_vshor[j][i]*Vshor_avg; /* add -grad(i-1) to build grad(i)-grad(i-1) */
grad_vshor[j][i]=r_LBFGS[l+2*NY*NX]*Vshor_avg; /* Denormalization */
}
}
}
}
......@@ -221,6 +242,11 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
sprintf(jac,"%s_grad2_vp_it%d",JACOBIAN,iteration);
write_matrix_disk(grad_vp, jac);
}
if(WAVETYPE==2 && VTI) {
sprintf(jac,"%s_grad2_vshor_it%d",JACOBIAN,iteration);
write_matrix_disk(grad_vshor, jac);
}
}
void lbfgs_core(int iteration, int N_LBFGS, int NPAR_LBFGS,float ** s_LBFGS, float ** y_LBFGS, float * rho_LBFGS,float *q_LBFGS,float *alpha_LBFGS,float *r_LBFGS) {
......
......@@ -9,6 +9,7 @@
MODEL = ../genmod/1D_linear_gradient_visc.c
MODEL_AC = ../genmod/1D_linear_gradient_ac.c
MODEL_EL = ../genmod/1D_linear_gradient_el.c
MODEL_EL_VTI = ../genmod/1D_linear_gradient_el_vti.c
MODEL_VAC = ../genmod/1D_linear_gradient_viscac.c
EXEC= ../bin
......@@ -72,10 +73,12 @@ IFOS2D= \
window_cos.c \
alloc_sections.c \
calc_mat_change_test.c \
calc_mat_change_test_vti.c \
calc_res.c \
calc_misfit.c \
calc_opt_step.c \
calc_energy.c \
change_parameterization_vti.c \
checkfd.c \
checkfd_ssg_elastic.c \
checkfd_ssg_visc.c \
......@@ -91,6 +94,7 @@ IFOS2D= \
snap_ssg_SH.c \
seismo_ssg.c \
surface_elastic_PML.c \
surface_el_vti_PML.c \
surface_acoustic_PML.c \
surface_PML.c \
update_v_ssg.c \
......@@ -98,15 +102,20 @@ IFOS2D= \
update_v_PML_SH.c \
update_v_acoustic_PML.c \
prepare_update_s.c \
prepare_update_s_vti.c \
update_p_PML.c \
update_s_elastic_ssg.c \
update_s_elastic_PML.c \
update_s_el_vti_PML.c \
update_s_elastic_PML_SH.c \
update_s_el_vti_PML_SH.c \
update_s_visc_PML.c \
update_s_visc_PML_SH.c \
update_s_visc_vti_PML_SH.c \
av_mue.c \
av_rho.c \
av_tau.c \
av_c55c66.c \
median2D.c \
exchange_par.c \
info.c \
......@@ -120,9 +129,12 @@ IFOS2D= \
$(MODEL) \
$(MODEL_AC) \
$(MODEL_EL) \
$(MODEL_EL_VTI) \
$(MODEL_VAC) \
matcopy.c \
matcopy_vti.c \
matcopy_elastic.c \
matcopy_elastic_vti.c \
matcopy_acoustic.c \
mergemod.c \
max_grad.c \
......@@ -134,11 +146,14 @@ IFOS2D= \
output_source_signal.c \
PCG.c \
PCG_SH.c \
PCG_SH_vti.c \
PML_pro.c \
readdsk.c \
read_par_json.c \
readmod.c \
readmod_vti.c \
readmod_elastic.c \
readmod_el_vti.c \
readmod_acoustic.c \
receiver.c \
rd_sour.c \
......
This diff is collapsed.
......@@ -38,7 +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 int GRAD_METHOD;
extern int GRAD_METHOD, WOLFE_CONDITION;
extern int WORKFLOW_STAGE;
/******************/
......@@ -54,6 +54,13 @@ void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[ST
}
}
/* PCG or LBFGS */
if( TIME_FILT > 0 ) {
if( GRAD_METHOD != workflow[WORKFLOW_STAGE][1] ) *LBFGS_iter_start=*iter;
GRAD_METHOD=workflow[WORKFLOW_STAGE][1];
if(workflow[WORKFLOW_STAGE][1]==1) WOLFE_CONDITION=0;
}
/* Inversion of material parameter */
if(workflow[WORKFLOW_STAGE][2]!=-1) {
if(workflow[WORKFLOW_STAGE][2]==1) {
......
/*-----------------------------------------------------------------------------------------
* 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>.
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* averaging of elastic constants c55 and c66
* ----------------------------------------------------------------------*/
#include "fd.h"
void av_c55c66(float **c55, float **c55ipjp, float **c66, float **c66ipjp){
extern int NX, NY;
int i, j;
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
// c55ipjp[j][i] = 0.25*(c55[j][i]+c55[j][i+1]+c55[j+1][i]+c55[j+1][i+1]);
// c66ipjp[j][i] = 0.25*(c66[j][i]+c66[j][i+1]+c66[j+1][i]+c66[j+1][i+1]);
c55ipjp[j][i] = 4.0/(1.0/c55[j][i]+1.0/c55[j][i+1]+1.0/c55[j+1][i]+1.0/c55[j+1][i+1]);
if((c55[j][i]==0.0)||(c55[j][i+1]==0.0)||(c55[j+1][i]==0.0)||(c55[j+1][i+1]==0.0)){
c55ipjp[j][i]=0.0;
}
c66ipjp[j][i] = 4.0/(1.0/c66[j][i]+1.0/c66[j][i+1]+1.0/c66[j+1][i]+1.0/c66[j+1][i+1]);
if((c66[j][i]==0.0)||(c66[j][i+1]==0.0)||(c66[j+1][i]==0.0)||(c66[j+1][i+1]==0.0)){
c66ipjp[j][i]=0.0;
}
}
}
}
This diff is collapsed.
/* 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>.
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* transforms parameters (rho,vp,vs,epsilon,gamma,delta)
* into (rho,c11,c13,c33,c55,c66)
* and averages c55 and c66
* ----------------------------------------------------------------------*/
#include "fd.h"
void change_parameterization_vti(float **rho, float **pi, float **u, float **epsilon, float **delta, float **vshor,
float **c11, float **c13, float **c33, float **c55, float **c66) {
extern int NX,NY;
int i,j;
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
c33[j][i]=pi[j][i]*pi[j][i]*rho[j][i];
c55[j][i]=u[j][i]*u[j][i]*rho[j][i];
c11[j][i]=c33[j][i]*(1+2*epsilon[j][i]);
c13[j][i]=sqrt(2*delta[j][i]*c33[j][i]*(c33[j][i]-c55[j][i])+(c33[j][i]-c55[j][i])*(c33[j][i]-c55[j][i]))-c55[j][i];
c66[j][i]=vshor[j][i]*vshor[j][i]*rho[j][i];
}
}
}
\ No newline at end of file
......@@ -26,11 +26,12 @@
#include <limits.h>
#include "fd.h"
void checkfd(FILE *fp, float ** prho, float ** ppi, float ** pu, float ** ptaus, float ** ptaup, float *peta, float *hc, float **srcpos, int nsrc, int **recpos, int ntr){
void checkfd(FILE *fp, float ** prho, float ** ppi, float ** pu, float ** ptaus, float ** ptaup, float *peta, float *hc, float **srcpos, int nsrc, int **recpos, int ntr,
float ** pepsilon, float ** pdelta, float ** pvshor) {
extern float DH, DT, TS;
extern float XREC1, XREC2, YREC1, YREC2;
extern int NX, NY, MYID, PARAMETERIZATION, FW, L, NT, NDT, ACOUSTIC;
extern int NX, NY, MYID, PARAMETERIZATION, FW, VTI, L, NT, NDT, ACOUSTIC;
extern int READREC, NPROCX,NPROCY, SRCREC, FREE_SURF;
extern int SEISMO, SEIS_FORMAT[6];
extern char SOURCE_FILE[STRING_SIZE], REC_FILE[STRING_SIZE];
......@@ -38,6 +39,7 @@ void checkfd(FILE *fp, float ** prho, float ** ppi, float ** pu, float ** ptaus,
/* local variables */
float c = 0.0, cmax_p=0.0, cmin_p=1e9, cmax_s=0.0, cmin_s=1e9, fmax, gamma;
float chor, cmax_phor=0.0, cmin_phor=1e9, cmax_shor=0, cmin_shor=1e9;
float cmax=0.0, cmin=1e9, dtstab, dhstab, cmax_r, cmin_r;
float sumu, sumpi, ws, ts, ppi_ref = 0.0, pu_ref;
int nfw=iround(FW/DH);
......@@ -89,17 +91,35 @@ void checkfd(FILE *fp, float ** prho, float ** ppi, float ** pu, float ** ptaus,
}
}
} else { /*L=0, elastic*/
for (i=1+nfw;i<=(nx-nfw);i++){
for (j=ny1;j<=(ny-nfw);j++){
if(PARAMETERIZATION==3){
c=sqrt(pu[j][i]/prho[j][i]);
} else {
c=pu[j][i];
if (VTI) {
for (i=1+nfw; i<= (nx-nfw); i++) {
for (j=ny1; j<= (ny-nfw); j++) {
c=pu[j][i]; /* vertical s-velocity */
chor=pvshor[j][i]; /* horizontal s-velocity */
if (cmax_s<c) cmax_s=c;
if (cmin_s>c) cmin_s=c;
if (cmax_shor<chor) cmax_shor=chor;
if (cmin_shor>chor) cmin_shor=chor;
if (cmax_shor>cmax_s) cmax_s=cmax_shor;
if (cmin_shor<cmin_s) cmin_s=cmin_shor;
}
}
}
else {
for (i=1+nfw;i<=(nx-nfw);i++){
for (j=ny1;j<=(ny-nfw);j++){
if (cmax_s<c) cmax_s=c;
if (cmin_s>c) cmin_s=c;
if(PARAMETERIZATION==3){
c=sqrt(pu[j][i]/prho[j][i]);
} else {
c=pu[j][i];
}
if (cmax_s<c) cmax_s=c;
if (cmin_s>c) cmin_s=c;
}
}
}
}
......@@ -137,20 +157,38 @@ void checkfd(FILE *fp, float ** prho, float ** ppi, float ** pu, float ** ptaus,
}
}
else{ /*L=0, elastic*/
for (i=1+nfw;i<=(nx-nfw);i++){
for (j=ny1;j<=(ny-nfw);j++){
if (VTI) {
for (i=1+nfw; i<= (nx-nfw); i++) {
for (j=ny1; j<= (ny-nfw); j++) {
c=ppi[j][i]; /* vertical p-velocity */
chor=c*(1+pepsilon[j][i]); /* horizontal p-velocity */
if (cmax_p<c) cmax_p=c;
if (cmin_p>c) cmin_p=c;
if (cmax_phor<chor) cmax_phor=chor;
if (cmin_phor>chor) cmin_phor=chor;
if (cmax_phor>cmax_p) cmax_p=cmax_phor;
if (cmin_phor<cmin_p) cmin_p=cmin_phor;
}
}
}
else {
for (i=1+nfw;i<=(nx-nfw);i++){
for (j=ny1;j<=(ny-nfw);j++){
if(PARAMETERIZATION==3){
c=sqrt((ppi[j][i]+2.0*pu[j][i])/prho[j][i]);}
if(PARAMETERIZATION==3){
c=sqrt((ppi[j][i]+2.0*pu[j][i])/prho[j][i]);}
if(PARAMETERIZATION==1){
c=ppi[j][i];}
if(PARAMETERIZATION==1){
c=ppi[j][i];}
if (cmax_p<c) cmax_p=c;
if (cmin_p>c) cmin_p=c;
if (cmax_p<c) cmax_p=c;
if (cmin_p>c) cmin_p=c;
}
}
}
}
}
......
......@@ -25,10 +25,10 @@
void exchange_par(void){
/* declaration of extern variables */
extern int NX, NY, FDORDER, MAXRELERROR, SOURCE_SHAPE, SOURCE_TYPE, SNAP, SNAP_FORMAT, L;
extern int NX, NY, FDORDER, MAXRELERROR, SOURCE_SHAPE, SOURCE_TYPE, SNAP, SNAP_FORMAT, L, VTI;
extern float DH, TIME, DT, TS, *FL, TAU, VPPML, PLANE_WAVE_DEPTH, PHI, F_REF;
extern float XREC1, XREC2, YREC1, YREC2, FPML;
extern float REC_ARRAY_DEPTH, REC_ARRAY_DIST, MUN, EPSILON, EPSILON_u, EPSILON_rho;
extern float REC_ARRAY_DEPTH, REC_ARRAY_DIST, MUN, EPSILON, EPSILON_u, EPSILON_rho, EPSILON_vshor;
extern int SEISMO, NDT, NGEOPH, SEIS_FORMAT, FREE_SURF, READMOD, READREC, SRCREC;
extern int BOUNDARY, REC_ARRAY, DRX, FW;
extern int SNAPSHOT_START,SNAPSHOT_END,SNAPSHOT_INCR;
......@@ -51,11 +51,11 @@ void exchange_par(void){
extern float SRTRADIUS;
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 INV_RHO_ITER, INV_VP_ITER, INV_VS_ITER, INV_VSHOR_ITER;
extern int MIN_ITER;
extern int nfstart, nf;
extern int nfstart_jac, nf_jac;
extern float VPUPPERLIM, VPLOWERLIM, VSUPPERLIM, VSLOWERLIM, RHOUPPERLIM, RHOLOWERLIM;
extern float VPUPPERLIM, VPLOWERLIM, VSUPPERLIM, VSLOWERLIM, RHOUPPERLIM, RHOLOWERLIM, VSHORUPPERLIM, VSHORLOWERLIM;
extern float npower, k_max_PML;
extern int INV_STF, N_STF, N_STF_START;
extern char PARA[STRING_SIZE];
......@@ -74,7 +74,7 @@ void exchange_par(void){
extern char MISFIT_LOG_FILE[STRING_SIZE];
extern int VELOCITY;
extern float WATERLEVEL_LNORM8;
extern float VP_VS_RATIO;
extern float VP_VS_RATIO, GAMMA_VTI;
extern int S;
extern float S_VP, S_VS, S_RHO;
extern int GRAD_FILT_WAVELENGTH;
......@@ -162,6 +162,7 @@ void exchange_par(void){
fdum[27] = EPSILON;
fdum[28] = EPSILON_u;
fdum[29] = EPSILON_rho;
fdum[30] = EPSILON_vshor;
fdum[31] = FPML;
fdum[32] = SRTRADIUS;
......@@ -221,6 +222,10 @@ void exchange_par(void){
fdum[69]=TRKILL_OFFSET_UPPER;
fdum[70]=LBFGS_SCALE_GRADIENTS;
fdum[71]=VSHORUPPERLIM;
fdum[72]=VSHORLOWERLIM;
fdum[73]=GAMMA_VTI;
/***********/
/* Integer */
......@@ -253,6 +258,8 @@ void exchange_par(void){
idum[26] = NGEOPH;
idum[27] = NDT;
idum[28] = SEIS_FORMAT;
idum[29] = VTI;
idum[31] = FDORDER;
......@@ -377,6 +384,8 @@ void exchange_par(void){
idum[116]=TRKILL_STF_OFFSET_INVERT;
idum[117]=JOINT_EQUAL_WEIGHTING;
idum[118]=INV_VSHOR_ITER;
} /** if (MYID == 0) **/
......@@ -444,6 +453,7 @@ void exchange_par(void){
EPSILON = fdum[27];
EPSILON_u = fdum[28];
EPSILON_rho = fdum[29];
EPSILON_vshor = fdum[30];
FPML = fdum[31];
SRTRADIUS = fdum[32];
......@@ -506,6 +516,10 @@ void exchange_par(void){
LBFGS_SCALE_GRADIENTS=fdum[70];
VSHORUPPERLIM=fdum[71];
VSHORLOWERLIM=fdum[72];
GAMMA_VTI=fdum[73];
/***********/
/* Integer */
/***********/
......@@ -537,6 +551,8 @@ void exchange_par(void){
NGEOPH = idum[26];
NDT = idum[27];
SEIS_FORMAT = idum[28];
VTI = idum[29];
FDORDER = idum[31];
......@@ -665,6 +681,8 @@ void exchange_par(void){
JOINT_EQUAL_WEIGHTING=idum[117];
INV_VSHOR_ITER=idum[118];
if ( MYID!=0 && L>0 ) {
FL=vector(1,L);
}
......
This diff is collapsed.
......@@ -15,7 +15,7 @@ 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 L, BOUNDARY, DC, DRX, NXG, NYG, IDX, IDY, FDORDER, MAXRELERROR;
int VTI, 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];
char LOG_FILE[STRING_SIZE];
......@@ -45,7 +45,7 @@ int ITERMAX, REC1, REC2, PARAMETERIZATION, FORWARD_ONLY, ADJOINT_TYPE;
int GRAD_METHOD;
float TSHIFT_back;
int MODEL_FILTER, FILT_SIZE;
float EPSILON, MUN, EPSILON_u, EPSILON_rho;
float EPSILON, MUN, EPSILON_u, EPSILON_rho, EPSILON_vshor;
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;
......@@ -55,7 +55,7 @@ char TAPER_FILE_NAME[STRING_SIZE];
int SPATFILTER, SPAT_FILT_SIZE, SPAT_FILT_1, SPAT_FILT_ITER;
int INV_RHO_ITER, INV_VS_ITER, INV_VP_ITER;
int INV_RHO_ITER, INV_VS_ITER, INV_VP_ITER, INV_VSHOR_ITER;
int MIN_ITER;
......@@ -64,7 +64,7 @@ int nfstart, nf;
int nfstart_jac, nf_jac;
float VPUPPERLIM, VPLOWERLIM, VSUPPERLIM, VSLOWERLIM, RHOUPPERLIM, RHOLOWERLIM;
float VPUPPERLIM, VPLOWERLIM, VSUPPERLIM, VSLOWERLIM, RHOUPPERLIM, RHOLOWERLIM, VSHORUPPERLIM, VSHORLOWERLIM;
float npower, k_max_PML;
......@@ -129,6 +129,8 @@ float WATERLEVEL_LNORM8;
float VP_VS_RATIO;
float GAMMA_VTI;
int S;
float S_VP, S_VS, S_RHO;
......
/*-----------------------------------------------------------------------------------------
* 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 pc13blished by
* the Free Software Foc13ndation, 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 Pc13blic 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>.
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* For the averaging of material properties each process requires valc13es
* at the indices 0 and NX+1 etc. These lie on the neighbouring processes.
* Thus, they have to be coc11ed which is done by this function.
*
*
* ----------------------------------------------------------------------*/
#include "fd.h"
void matcopy_elastic_vti(float ** rho, float ** c11, float ** c13, float ** c33, float ** c55, float ** c66){
extern int MYID, NX, NY, INDEX[5],VERBOSE;
extern const int TAG1,TAG2,TAG5,TAG6;
extern FILE *FP;
MPI_Status status;
double time1, time2;
int i, j;
float ** bufferlef_to_rig_1, ** bufferrig_to_lef_1;
float ** buffertop_to_bot_1, ** bufferbot_to_top_1;
bufferlef_to_rig_1 = matrix(0,NY+1,1,6);
bufferrig_to_lef_1 = matrix(0,NY+1,1,6);
buffertop_to_bot_1 = matrix(0,NX+1,1,6);
bufferbot_to_top_1 = matrix(0,NX+1,1,6);
if(VERBOSE){
fprintf(FP,"\n\n **Message from matcopy (written by PE %d):",MYID);
fprintf(FP,"\n Copy material properties at inner boundaries ... \n");
time1=MPI_Wtime();
}
/* if (POS[2]!=0)*/ /* no boundary exchange at top of global grid */
for (i=0;i<=NX+1;i++){
/* storage of top of local volc13me into buffer */
buffertop_to_bot_1[i][1] = rho[1][i];
buffertop_to_bot_1[i][2] = c11[1][i];
buffertop_to_bot_1[i][3] = c13[1][i];
buffertop_to_bot_1[i][4] = c33[1][i];
buffertop_to_bot_1[i][5] = c55[1][i];
buffertop_to_bot_1[i][6] = c66[1][i];
}
/* if (POS[2]!=NPROCY-1)*/ /* no boc13ndary exchange at bottom of global grid */
for (i=0;i<=NX+1;i++){
/* storage of bottom of local volume into buffer */
bufferbot_to_top_1[i][1] = rho[NY][i];
bufferbot_to_top_1[i][2] = c11[NY][i];
bufferbot_to_top_1[i][3] = c13[NY][i];
bufferbot_to_top_1[i][4] = c33[NY][i];
bufferbot_to_top_1[i][5] = c55[NY][i];
bufferbot_to_top_1[i][6] = c66[NY][i];
}
/*=========sending and receiving of the boundaries==========*/
MPI_Bsend(&buffertop_to_bot_1[0][1],(NX+2)*6,MPI_FLOAT,INDEX[3],TAG5,MPI_COMM_WORLD);
MPI_Barrier(MPI_COMM_WORLD);
MPI_Recv(&buffertop_to_bot_1[0][1],(NX+2)*6,MPI_FLOAT,INDEX[4],TAG5,MPI_COMM_WORLD,&status);
MPI_Bsend(&bufferbot_to_top_1[0][1],(NX+2)*6,MPI_FLOAT,INDEX[4],TAG6,MPI_COMM_WORLD);
MPI_Barrier(MPI_COMM_WORLD);
MPI_Recv(&bufferbot_to_top_1[0][1],(NX+2)*6,MPI_FLOAT,INDEX[3],TAG6,MPI_COMM_WORLD,&status);
/* if (POS[2]!=NPROCY-1)*/ /* no boundary exchange at bottom of global grid */
for (i=0;i<=NX+1;i++){
rho[NY+1][i] = buffertop_to_bot_1[i][1];
c11[NY+1][i] = buffertop_to_bot_1[i][2];
c13[NY+1][i] = buffertop_to_bot_1[i][3];
c33[NY+1][i] = buffertop_to_bot_1[i][4];
c55[NY+1][i] = buffertop_to_bot_1[i][5];
c66[NY+1][i] = buffertop_to_bot_1[i][6];
}
/* if (POS[2]!=0)*/ /* no boc13ndary exchange at top of global grid */
for (i=0;i<=NX+1;i++){
rho[0][i] = bufferbot_to_top_1[i][1];
c11[0][i] = bufferbot_to_top_1[i][2];
c13[0][i] = bufferbot_to_top_1[i][3];
c33[0][i] = bufferbot_to_top_1[i][4];
c55[0][i] = bufferbot_to_top_1[i][5];
c66[0][i] = bufferbot_to_top_1[i][6];
}
/* if (POS[1]!=0)*/ /* no boundary exchange at left edge of global grid */
for (j=0;j<=NY+1;j++)
{
/* storage of left edge of local volc13me into buffer */
bufferlef_to_rig_1[j][1] = rho[j][1];
bufferlef_to_rig_1[j][2] = c11[j][1];
bufferlef_to_rig_1[j][3] = c13[j][1];
bufferlef_to_rig_1[j][4] = c33[j][1];
bufferlef_to_rig_1[j][5] = c55[j][1];
bufferlef_to_rig_1[j][6] = c66[j][1];
}
/* if (POS[1]!=NPROCX-1)*/ /* no boundary exchange at right edge of global grid */
for (j=0;j<=NY+1;j++){
/* storage of right edge of local volume into buffer */
bufferrig_to_lef_1[j][1] = rho[j][NX];
bufferrig_to_lef_1[j][2] = c11[j][NX];
bufferrig_to_lef_1[j][3] = c13[j][NX];
bufferrig_to_lef_1[j][4] = c33[j][NX];
bufferrig_to_lef_1[j][5] = c55[j][NX];
bufferrig_to_lef_1[j][6] = c66[j][NX];
}
MPI_Bsend(&bufferlef_to_rig_1[0][1],(NY+2)*6,MPI_FLOAT,INDEX[1],TAG1,MPI_COMM_WORLD);
MPI_Barrier(MPI_COMM_WORLD);
MPI_Recv(&bufferlef_to_rig_1[0][1],(NY+2)*6,MPI_FLOAT,INDEX[2],TAG1,MPI_COMM_WORLD,&status);
MPI_Bsend(&bufferrig_to_lef_1[0][1],(NY+2)*6,MPI_FLOAT,INDEX[2],TAG2,MPI_COMM_WORLD);
MPI_Barrier(MPI_COMM_WORLD);
MPI_Recv(&bufferrig_to_lef_1[0][1],(NY+2)*6,MPI_FLOAT,INDEX[1],TAG2,MPI_COMM_WORLD,&status);
/* if (POS[1]!=NPROCX-1)*/ /* no boundary exchange at right edge of global grid */
for (j=0;j<=NY+1;j++){
rho[j][NX+1] = bufferlef_to_rig_1[j][1];
c11[j][NX+1] = bufferlef_to_rig_1[j][2];
c13[j][NX+1] = bufferlef_to_rig_1[j][3];
c33[j][NX+1] = bufferlef_to_rig_1[j][4];
c55[j][NX+1] = bufferlef_to_rig_1[j][5];
c66[j][NX+1] = bufferlef_to_rig_1[j][6];
}
/* if (POS[1]!=0)*/ /* no boundary exchange at left edge of global grid */
for (j=0;j<=NY+1;j++){
rho[j][0] = bufferrig_to_lef_1[j][1];
c11[j][0] = bufferrig_to_lef_1[j][2];
c13[j][0] = bufferrig_to_lef_1[j][3];
c33[j][0] = bufferrig_to_lef_1[j][4];
c55[j][0] = bufferrig_to_lef_1[j][5];
c66[j][0] = bufferrig_to_lef_1[j][6];
}
if (MYID==0&&VERBOSE){
time2=MPI_Wtime();
fprintf(FP," finished (real time: %4.2f s).\n",time2-time1);
}
free_matrix(bufferlef_to_rig_1,0,NY+1,1,6);
free_matrix(bufferrig_to_lef_1,0,NY+1,1,6);
free_matrix(buffertop_to_bot_1,0,NX+1,1,6);
free_matrix(bufferbot_to_top_1,0,NX+1,1,6);
}
/*-----------------------------------------------------------------------------------------
* 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 pc13blished by
* the Free Software Foc13ndation, 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 Pc13blic 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>.
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* For the averaging of material properties each process requires valc13es
* at the indices 0 and NX+1 etc. These lie on the neighbouring processes.
* Thus, they have to be coc11ed which is done by this function.
*
*
* ----------------------------------------------------------------------*/
#include "fd.h"
void matcopy_vti(float ** rho, float ** c11, float ** c13, float ** c33, float ** c55, float ** c66, float ** taus, float **taup){
extern int MYID, NX, NY, INDEX[5],VERBOSE;
extern const int TAG1,TAG2,TAG5,TAG6;
extern FILE *FP;
MPI_Status status;
double time1, time2;
int i, j;
float ** bufferlef_to_rig_1, ** bufferrig_to_lef_1;
float ** buffertop_to_bot_1, ** bufferbot_to_top_1;
bufferlef_to_rig_1 = matrix(0,NY+1,1,8);
bufferrig_to_lef_1 = matrix(0,NY+1,1,8);
buffertop_to_bot_1 = matrix(0,NX+1,1,8);
bufferbot_to_top_1 = matrix(0,NX+1,1,8);
if(VERBOSE){
fprintf(FP,"\n\n **Message from matcopy (written by PE %d):",MYID);
fprintf(FP,"\n Copy material properties at inner boundaries ... \n");
time1=MPI_Wtime();
}
/* if (POS[2]!=0)*/ /* no boundary exchange at top of global grid */
for (i=0;i<=NX+1;i++){
/* storage of top of local volc13me into buffer */
buffertop_to_bot_1[i][1] = rho[1][i];
buffertop_to_bot_1[i][2] = c11[1][i];
buffertop_to_bot_1[i][3] = c13[1][i];
buffertop_to_bot_1[i][4] = c33[1][i];
buffertop_to_bot_1[i][5] = c55[1][i];
buffertop_to_bot_1[i][6] = c66[1][i];
buffertop_to_bot_1[i][7] = taus[1][i];
buffertop_to_bot_1[i][8] = taup[1][i];
}
/* if (POS[2]!=NPROCY-1)*/ /* no boc13ndary exchange at bottom of global grid */
for (i=0;i<=NX+1;i++){
/* storage of bottom of local volume into buffer */
bufferbot_to_top_1[i][1] = rho[NY][i];
bufferbot_to_top_1[i][2] = c11[NY][i];
bufferbot_to_top_1[i][3] = c13[NY][i];
bufferbot_to_top_1[i][4] = c33[NY][i];
bufferbot_to_top_1[i][5] = c55[NY][i];
bufferbot_to_top_1[i][6] = c66[NY][i];
bufferbot_to_top_1[i][7] = taus[NY][i];
bufferbot_to_top_1[i][8] = taup[NY][i];
}
/*=========sending and receiving of the boundaries==========*/
MPI_Bsend(&buffertop_to_bot_1[0][1],(NX+2)*8,MPI_FLOAT,INDEX[3],TAG5,MPI_COMM_WORLD);
MPI_Barrier(MPI_COMM_WORLD);
MPI_Recv(&buffertop_to_bot_1[0][1],(NX+2)*8,MPI_FLOAT,INDEX[4],TAG5,MPI_COMM_WORLD,&status);
MPI_Bsend(&bufferbot_to_top_1[0][1],(NX+2)*8,MPI_FLOAT,INDEX[4],TAG6,MPI_COMM_WORLD);
MPI_Barrier(MPI_COMM_WORLD);
MPI_Recv(&bufferbot_to_top_1[0][1],(NX+2)*8,MPI_FLOAT,INDEX[3],TAG6,MPI_COMM_WORLD,&status);
/* if (POS[2]!=NPROCY-1)*/ /* no boundary exchange at bottom of global grid */
for (i=0;i<=NX+1;i++){
rho[NY+1][i] = buffertop_to_bot_1[i][1];
c11[NY+1][i] = buffertop_to_bot_1[i][2];
c13[NY+1][i] = buffertop_to_bot_1[i][3];
c33[NY+1][i] = buffertop_to_bot_1[i][4];
c55[NY+1][i] = buffertop_to_bot_1[i][5];
c66[NY+1][i] = buffertop_to_bot_1[i][6];
taus[NY+1][i] = buffertop_to_bot_1[i][7];
taup[NY+1][i] = buffertop_to_bot_1[i][8];
}
/* if (POS[2]!=0)*/ /* no boc13ndary exchange at top of global grid */
for (i=0;i<=NX+1;i++){
rho[0][i] = bufferbot_to_top_1[i][1];
c11[0][i] = bufferbot_to_top_1[i][2];
c13[0][i] = bufferbot_to_top_1[i][3];
c33[0][i] = bufferbot_to_top_1[i][4];
c55[0][i] = bufferbot_to_top_1[i][5];
c66[0][i] = bufferbot_to_top_1[i][6];
taus[0][i] = bufferbot_to_top_1[i][7];
taup[0][i] = bufferbot_to_top_1[i][8];
}
/* if (POS[1]!=0)*/ /* no boundary exchange at left edge of global grid */
for (j=0;j<=NY+1;j++)
{
/* storage of left edge of local volc13me into buffer */
bufferlef_to_rig_1[j][1] = rho[j][1];
bufferlef_to_rig_1[j][2] = c11[j][1];
bufferlef_to_rig_1[j][3] = c13[j][1];
bufferlef_to_rig_1[j][4] = c33[j][1];
bufferlef_to_rig_1[j][5] = c55[j][1];
bufferlef_to_rig_1[j][6] = c66[j][1];
bufferlef_to_rig_1[j][7] = taus[j][1];
bufferlef_to_rig_1[j][8] = taup[j][1];
}
/* if (POS[1]!=NPROCX-1)*/ /* no boundary exchange at right edge of global grid */
for (j=0;j<=NY+1;j++){
/* storage of right edge of local volume into buffer */
bufferrig_to_lef_1[j][1] = rho[j][NX];
bufferrig_to_lef_1[j][2] = c11[j][NX];
bufferrig_to_lef_1[j][3] = c13[j][NX];
bufferrig_to_lef_1[j][4] = c33[j][NX];
bufferrig_to_lef_1[j][5] = c55[j][NX];
bufferrig_to_lef_1[j][6] = c66[j][NX];
bufferrig_to_lef_1[j][7] = taus[j][NX];
bufferrig_to_lef_1[j][8] = taup[j][NX];
}
MPI_Bsend(&bufferlef_to_rig_1[0][1],(NY+2)*8,MPI_FLOAT,INDEX[1],TAG1,MPI_COMM_WORLD);
MPI_Barrier(MPI_COMM_WORLD);
MPI_Recv(&bufferlef_to_rig_1[0][1],(NY+2)*8,MPI_FLOAT,INDEX[2],TAG1,MPI_COMM_WORLD,&status);
MPI_Bsend(&bufferrig_to_lef_1[0][1],(NY+2)*8,MPI_FLOAT,INDEX[2],TAG2,MPI_COMM_WORLD);
MPI_Barrier(MPI_COMM_WORLD);
MPI_Recv(&bufferrig_to_lef_1[0][1],(NY+2)*8,MPI_FLOAT,INDEX[1],TAG2,MPI_COMM_WORLD,&status);
/* if (POS[1]!=NPROCX-1)*/ /* no boundary exchange at right edge of global grid */
for (j=0;j<=NY+1;j++){
rho[j][NX+1] = bufferlef_to_rig_1[j][1];
c11[j][NX+1] = bufferlef_to_rig_1[j][2];
c13[j][NX+1] = bufferlef_to_rig_1[j][3];
c33[j][NX+1] = bufferlef_to_rig_1[j][4];
c55[j][NX+1] = bufferlef_to_rig_1[j][5];
c66[j][NX+1] = bufferlef_to_rig_1[j][6];
taus[j][NX+1] = bufferlef_to_rig_1[j][7];
taup[j][NX+1] = bufferlef_to_rig_1[j][8];
}
/* if (POS[1]!=0)*/ /* no boundary exchange at left edge of global grid */
for (j=0;j<=NY+1;j++){
rho[j][0] = bufferrig_to_lef_1[j][1];
c11[j][0] = bufferrig_to_lef_1[j][2];
c13[j][0] = bufferrig_to_lef_1[j][3];
c33[j][0] = bufferrig_to_lef_1[j][4];
c55[j][0] = bufferrig_to_lef_1[j][5];
c66[j][0] = bufferrig_to_lef_1[j][6];
taus[j][0] = bufferrig_to_lef_1[j][7];
taup[j][0] = bufferrig_to_lef_1[j][8];
}
if (MYID==0&&VERBOSE){
time2=MPI_Wtime();
fprintf(FP," finished (real time: %4.2f s).\n",time2-time1);
}
free_matrix(bufferlef_to_rig_1,0,NY+1,1,8);
free_matrix(bufferrig_to_lef_1,0,NY+1,1,8);
free_matrix(buffertop_to_bot_1,0,NX+1,1,8);
free_matrix(bufferbot_to_top_1,0,NX+1,1,8);
}
/*-----------------------------------------------------------------------------------------
* 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>.
-----------------------------------------------------------------------------------------*/
#include "fd.h"
void prepare_update_s_vti(float *etajm, float *etaip, float *peta, float **fipjp, float **pc55,
float **pc66ipjp, float **ppi, float **prho, float **ptaus, float **ptaup,
float **ptausipjp, float **f, float **g, float *bip, float *bjm,
float *cip, float *cjm, float ***dip, float ***d, float ***e) {
extern int NX, NY, L, PARAMETERIZATION, MYID;
extern float DT, *FL;
int i, j, l;
extern char MFILE[STRING_SIZE];
extern float F_REF;
extern int VERBOSE;
float sumu, sumpi, ws, *pts;
float c55, pi, c66ipjp;
/* vector for maxwellbodies */
pts=vector(1,L);
for (l=1;l<=L;l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
}
/*ws=2.0*PI*FL[1];*/
ws=2.0*PI*F_REF;
if(MYID==0&&VERBOSE)printf("MYID %d: F_REF = %f\n",MYID,F_REF);
sumu=0.0;
sumpi=0.0;
for (l=1;l<=L;l++){
sumu=sumu+((ws*ws*pts[l]*pts[l])/(1.0+ws*ws*pts[l]*pts[l]));
sumpi=sumpi+((ws*ws*pts[l]*pts[l])/(1.0+ws*ws*pts[l]*pts[l]));
}
for (l=1;l<=L;l++){
etajm[l] = peta[l];
etaip[l] = peta[l];
}
if (PARAMETERIZATION==1){
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
c55=pc55[j][i]/(1.0+sumu*ptaus[j][i]);
pi=(ppi[j][i]*ppi[j][i]*prho[j][i])/(1.0+sumpi*ptaup[j][i]);
c66ipjp=pc66ipjp[j][i]/(1.0+sumu*ptausipjp[j][i]);
fipjp[j][i] = c66ipjp*DT*(1.0+L*ptausipjp[j][i]);
f[j][i] = c55*DT*(1.0+L*ptaus[j][i]);
g[j][i] = pi*DT*(1.0+L*ptaup[j][i]);
for (l=1;l<=L;l++){
bip[l] = 1.0/(1.0+(etaip[l]*0.5));
bjm[l] = 1.0/(1.0+(etajm[l]*0.5));
cip[l] = 1.0-(etaip[l]*0.5);
cjm[l] = 1.0-(etajm[l]*0.5);
dip[j][i][l] = c66ipjp*etaip[l]*ptausipjp[j][i];
d[j][i][l] = c55*etajm[l]*ptaus[j][i];
e[j][i][l] = pi*etajm[l]*ptaup[j][i];
}
}
}
}
/*if (PARAMETERIZATION==3){
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
mu=pu[j][i]/(1.0+sumu*ptaus[j][i]);
pi=(ppi[j][i]+2*pu[j][i])/(1.0+sumpi*ptaup[j][i]);
muipjp=puipjp[j][i]/(1.0+sumu*ptausipjp[j][i]);
fipjp[j][i] = muipjp*DT*(1.0+L*ptausipjp[j][i]);
f[j][i] = mu*DT*(1.0+L*ptaus[j][i]);
g[j][i] = pi*DT*(1.0+L*ptaup[j][i]);
for (l=1;l<=L;l++){
bip[l] = 1.0/(1.0+(etaip[l]*0.5));
bjm[l] = 1.0/(1.0+(etajm[l]*0.5));
cip[l] = 1.0-(etaip[l]*0.5);
cjm[l] = 1.0-(etajm[l]*0.5);
dip[j][i][l] = muipjp*etaip[l]*ptausipjp[j][i];
d[j][i][l] = mu*etajm[l]*ptaus[j][i];
e[j][i][l] = pi*etajm[l]*ptaup[j][i];
}
}
}
}*/
free_vector(pts,1,L);
}
......@@ -28,7 +28,7 @@ char ** varname_list,** value_list;
void read_par_json(FILE *fp, char *fileinp){
/* declaration of extern variables */
extern int NX, NY, FDORDER, MAXRELERROR, SOURCE_SHAPE,SOURCE_SHAPE_SH, SOURCE_TYPE, SNAP, SNAP_FORMAT, ACOUSTIC, L, VERBOSE, WAVETYPE,JOINT_INVERSION_PSV_SH_TYPE,JOINT_EQUAL_WEIGHTING;
extern int NX, NY, FDORDER, MAXRELERROR, SOURCE_SHAPE,SOURCE_SHAPE_SH, SOURCE_TYPE, SNAP, SNAP_FORMAT, ACOUSTIC, VTI, L, VERBOSE, WAVETYPE,JOINT_INVERSION_PSV_SH_TYPE,JOINT_EQUAL_WEIGHTING;
extern float DH, TIME, DT, TS, *FL, TAU, VPPML, PLANE_WAVE_DEPTH, PHI, F_REF,JOINT_INVERSION_PSV_SH_ALPHA_VS,JOINT_INVERSION_PSV_SH_ALPHA_RHO;
extern float XREC1, XREC2, YREC1, YREC2, FPML;
extern float REC_ARRAY_DEPTH, REC_ARRAY_DIST;
......@@ -52,11 +52,11 @@ void read_par_json(FILE *fp, char *fileinp){
extern float SRTRADIUS;
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 int INV_RHO_ITER, INV_VS_ITER, INV_VP_ITER, INV_VSHOR_ITER;
extern char INV_MODELFILE[STRING_SIZE];
extern int nfstart, nf;
extern int nfstart_jac, nf_jac;
extern float VPUPPERLIM, VPLOWERLIM, VSUPPERLIM, VSLOWERLIM, RHOUPPERLIM, RHOLOWERLIM;
extern float VPUPPERLIM, VPLOWERLIM, VSUPPERLIM, VSLOWERLIM, RHOUPPERLIM, RHOLOWERLIM, VSHORUPPERLIM, VSHORLOWERLIM;
extern float npower, k_max_PML;
extern int INV_STF, N_STF, N_STF_START;
......@@ -100,7 +100,7 @@ void read_par_json(FILE *fp, char *fileinp){
extern float WATERLEVEL_LNORM8;
extern float VP_VS_RATIO;
extern float VP_VS_RATIO, GAMMA_VTI;
extern int S;
......@@ -444,6 +444,9 @@ void read_par_json(FILE *fp, char *fileinp){
if (get_float_from_objectlist("F_REF",number_readobjects,&F_REF,varname_list, value_list)){
F_REF=-1.0;
fprintf(fp,"Reference frequency for viscoelastic modeling is set to center frequency of the source wavelet.\n");}
if (get_int_from_objectlist("VTI",number_readobjects,&VTI,varname_list, value_list))
declare_error("Variable VTI could not be retrieved from the json input file!");
}
......@@ -658,11 +661,21 @@ void read_par_json(FILE *fp, char *fileinp){
if (get_int_from_objectlist("INV_VS_ITER",number_readobjects,&INV_VS_ITER,varname_list, value_list)){
INV_VS_ITER=0;
fprintf(fp,"Variable INV_VS_ITER is set to default value %d.\n",INV_VS_ITER);}
/* Inversion for gamma */
if (get_int_from_objectlist("INV_VSHOR_ITER",number_readobjects,&INV_VSHOR_ITER,varname_list, value_list)){
INV_VSHOR_ITER=0;
fprintf(fp,"Variable INV_VSHOR_ITER is set to default value %d.\n",INV_VSHOR_ITER);}
/* Vp/Vs-Ratio */
if (get_float_from_objectlist("VP_VS_RATIO",number_readobjects,&VP_VS_RATIO,varname_list, value_list)){
VP_VS_RATIO=0.0;
fprintf(fp,"Variable VP_VS_RATIO is set to default value %4.2f which means that it is disregarded.\n",VP_VS_RATIO);}
/* Limited Thomsen parameter gamma */
if (get_float_from_objectlist("GAMMA_VTI",number_readobjects,&GAMMA_VTI,varname_list, value_list)){
GAMMA_VTI=0.0;
fprintf(fp,"Variable GAMMA_VTI is set to default value %4.2f which means that it is disregarded.\n",GAMMA_VTI);}
/* Limited update of model parameters in reference to the starting model */
if (get_int_from_objectlist("S",number_readobjects,&S,varname_list, value_list)){
......@@ -699,6 +712,12 @@ void read_par_json(FILE *fp, char *fileinp){
if (get_float_from_objectlist("RHOLOWERLIM",number_readobjects,&RHOLOWERLIM,varname_list, value_list)){
RHOLOWERLIM=0.0;
fprintf(fp,"Variable RHOLOWERLIM is set to default value %f.\n",RHOLOWERLIM);}
if (get_float_from_objectlist("VSHORUPPERLIM",number_readobjects,&VSHORUPPERLIM,varname_list, value_list)){
VSHORUPPERLIM=25000.0;
fprintf(fp,"Variable VSHORUPPERLIM is set to default value %f.\n",VSHORUPPERLIM);}
if (get_float_from_objectlist("VSHORLOWERLIM",number_readobjects,&VSHORLOWERLIM,varname_list, value_list)){
VSHORLOWERLIM=0.0;
fprintf(fp,"Variable VSHORLOWERLIM is set to default value %f.\n",VSHORLOWERLIM);}
/* Hessian and Gradient-Method */
......
......@@ -174,7 +174,7 @@ void readmod(float ** rho, float ** pi, float ** u, float ** taus, float **
if (sw_Qs){
if(feof(fp_vs)){
if(feof(fp_qs)){
declare_error("Model file QS is to small. Check dimensions NX*NY of file.");
}
......
/*-----------------------------------------------------------------------------------------
* 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>.
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* Read elastic model properties (vp,vs,density,delta,epsilon,gamma) from files
*
* ----------------------------------------------------------------------*/
/* This file contains function readmod, which has the purpose
to read data from model-files for viscoelastic simulation */
#include "fd.h"
void readmod_el_vti(float ** rho, float ** pi, float ** u, float **delta, float **epsilon, float **vshor){
extern int NX, NY, NXG, NYG, POS[3], MYID, PARAMETERIZATION,WAVETYPE;
extern char MFILE[STRING_SIZE];
extern FILE *FP;
/* local variables */
float rhov, vp = 0.0, vs, deltav, epsilonv, vshorv;
int i, j, ii, jj;
FILE *fp_vs, *fp_vp = NULL, *fp_rho, *fp_epsilon, *fp_delta, *fp_vshor;
char filename[STRING_SIZE];
fprintf(FP,"\n...reading model information from model-files...\n");
/* ----------------------------------- */
/* read density and seismic velocities */
/* ----------------------------------- */
if(PARAMETERIZATION==1){
if(WAVETYPE==1||WAVETYPE==3){
fprintf(FP,"\t Vp:\n\t %s.vp\n\n",MFILE);
sprintf(filename,"%s.vp",MFILE);
fp_vp=fopen(filename,"r");
if (fp_vp==NULL) declare_error(" Could not open model file for Vp ! ");
fprintf(FP,"\t Epsilon:\n\t %s.epsilon\n\n",MFILE);
sprintf(filename,"%s.epsilon",MFILE);
fp_epsilon=fopen(filename,"r");
if (fp_epsilon==NULL) declare_error(" Could not open model file for epsilon ! ");
fprintf(FP,"\t Delta:\n\t %s.delta\n\n",MFILE);
sprintf(filename,"%s.delta",MFILE);
fp_delta=fopen(filename,"r");
if (fp_delta==NULL) declare_error(" Could not open model file for delta ! ");
}
if(WAVETYPE==2||WAVETYPE==3){
fprintf(FP,"\t Vshor:\n\t %s.vshor\n\n",MFILE);
sprintf(filename,"%s.vshor",MFILE);
fp_vshor=fopen(filename,"r");
if (fp_vshor==NULL) declare_error(" Could not open model file for vshor ! ");
}
fprintf(FP,"\t Vs:\n\t %s.vs\n\n",MFILE);
sprintf(filename,"%s.vs",MFILE);
fp_vs=fopen(filename,"r");
if (fp_vs==NULL) declare_error(" Could not open model file for Vs ! ");
fprintf(FP,"\t Density:\n\t %s.rho\n\n",MFILE);
sprintf(filename,"%s.rho",MFILE);
fp_rho=fopen(filename,"r");
if (fp_rho==NULL) declare_error(" Could not open model file for densities ! ");
} else {
/* read density and Lame parameters */
/* ----------------------------------- */
fprintf(FP,"\t Lame parameter lambda:\n\t %s.lam\n\n",MFILE);
sprintf(filename,"%s.lam",MFILE);
fp_vp=fopen(filename,"r");
if (fp_vp==NULL) declare_error(" Could not open model file for Lame parameter lambda ! ");
fprintf(FP,"\t Lame parameter mu:\n\t %s.mu\n\n",MFILE);
sprintf(filename,"%s.mu",MFILE);
fp_vs=fopen(filename,"r");