Implementation of VTI forward solver and VTI inversion for SH component

doc/latex/6_Parameter_Definition.tex

@@ -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",

genmod/1D_linear_gradient_el_vti.c

+/*-----------------------------------------------------------------------------------------
+ * 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);
+}
+
+
+

par/in_and_out/IFOS2D_FW_all_parameters.json

@@ -62,6 +62,9 @@
 			"NDT" : "1",
 			"SEIS_FORMAT" : "1",
 			"SEIS_FILE" : "su/IFOS",
+			
+"Anisotropy" : "comment",
+            "VTI" : "0",				
 
 "Q-approximation" : "comment",
 			"L" : "0",

par/in_and_out/IFOS2D_INV_all_parameters.json

@@ -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",

src/LBFGS.c

@@ -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) {

src/Makefile

@@ -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 \

src/apply_workflow.c

@@ -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) {

src/av_c55c66.c

+/*-----------------------------------------------------------------------------------------
+ * 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;
+	    }
+        }
+    }
+}

src/change_parameterization_vti.c

+/* 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

src/checkfd.c

@@ -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;
+				}
 			}
-		}
-
+		}	
 	}
 
 

src/exchange_par.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;