...
 
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",
......
......@@ -40,7 +40,7 @@ int main(int argc, char **argv){
float muss, lamss;
float memdyn, memmodel, memseismograms, membuffer, memtotal, eps_scale;
float fac1, fac2;
float opteps_vp, opteps_vs, opteps_rho, Vp_avg, C_vp, Vs_avg, C_vs, rho_avg, C_rho;
float opteps_vp, opteps_vs, opteps_rho, Vp_avg, C_vp, Vs_avg, C_vs, rho_avg, C_rho, vshor_avg, C_vshor;
float memfwt, memfwt1, memfwtdata;
char *buff_addr, ext[10], *fileinp;
char jac[225];
......@@ -61,8 +61,9 @@ int main(int argc, char **argv){
float ** psxx, ** psxy, ** psyy, ** psxz, ** psyz, **psp, ** ux, ** uy, ** uxy, ** uyx, ** u, ** Vp0, ** uttx, ** utty, ** Vs0, ** Rho0;
float ** pvx, ** pvy, ** pvz, **waveconv, **waveconv_lam, **waveconv_mu, **waveconv_rho, **waveconv_rho_s, **waveconv_u, **waveconvtmp, **wcpart, **wavejac,**waveconv_rho_s_z,**waveconv_u_z,**waveconv_rho_z;
float **waveconv_shot, **waveconv_u_shot, **waveconv_rho_shot, **waveconv_u_shot_z, **waveconv_rho_shot_z;
float **waveconv_c55_shot_z, **waveconv_c66_shot_z, **waveconv_c55_s_z, **waveconv_c66_s_z, **waveconv_vshor_shot_z, **waveconv_vshor_z;
float ** pvxp1, ** pvyp1, ** pvzp1, ** pvxm1, ** pvym1, ** pvzm1;
float ** gradg, ** gradp,** gradg_rho, ** gradp_rho, ** gradg_u, ** gradp_u, ** gradp_u_z,** gradp_rho_z;
float ** gradg, ** gradp,** gradg_rho, ** gradp_rho, ** gradg_u, ** gradp_u, ** gradp_u_z,** gradp_rho_z, **gradp_vshor_z;
float ** prho,** prhonp1, **prip=NULL, **prjp=NULL, **pripnp1=NULL, **prjpnp1=NULL, ** ppi, ** pu, ** punp1, ** puipjp, ** ppinp1;
float ** vpmat, ***forward_prop_x, ***forward_prop_y, ***forward_prop_rho_x, ***forward_prop_u, ***forward_prop_rho_y, ***forward_prop_p;
......@@ -86,6 +87,10 @@ int main(int argc, char **argv){
float ** psi_sxx_x, ** psi_syy_y, ** psi_sxy_y, ** psi_sxy_x, ** psi_vxx, ** psi_vyy, ** psi_vxy, ** psi_vyx, ** psi_vxxs;
float ** psi_sxz_x, ** psi_syz_y, ** psi_vzx, ** psi_vzy;
/* Variables for vertically transversely isotropic modeling */
float **pc11=NULL, **pc13=NULL, **pc33=NULL, **pc55=NULL, **pc55ipjp=NULL, **pc66=NULL, **pc66ipjp=NULL;
float **pdelta=NULL, **pepsilon=NULL, **pgamma=NULL, **pvshornp1=NULL, **pvshor=NULL;
/* Variables for viscoelastic modeling */
float **ptaus=NULL, **ptaup=NULL, *etaip=NULL, *etajm=NULL, *peta=NULL, **ptausipjp=NULL, **fipjp=NULL, ***dip=NULL, *bip=NULL, *bjm=NULL;
float *cip=NULL, *cjm=NULL, ***d=NULL, ***e=NULL, ***pr=NULL, ***pp=NULL, ***pq=NULL, **f=NULL, **g=NULL;
......@@ -115,8 +120,8 @@ int main(int argc, char **argv){
int gradient_optimization=1;
float alpha_SL_min=0, alpha_SL_max=0, alpha_SL=1.0;
float alpha_SL_old;
float ** waveconv_old,** waveconv_u_old,** waveconv_rho_old;
float ** waveconv_up,** waveconv_u_up,** waveconv_rho_up;
float ** waveconv_old,** waveconv_u_old,** waveconv_rho_old,** waveconv_vshor_old;
float ** waveconv_up,** waveconv_u_up,** waveconv_rho_up,** waveconv_vshor_up;
float L2_SL_old=0, L2_SL_new=0;
float c1_SL=1e-4, c2_SL=0.9;
int wolfe_status;
......@@ -388,7 +393,7 @@ int main(int argc, char **argv){
if(GRAD_METHOD==2) {
/* Allocate memory for L-BFGS */
if(WAVETYPE==2) LBFGS_NPAR=2;
if(WAVETYPE==2 && !VTI) LBFGS_NPAR=2;
s_LBFGS=fmatrix(1,N_LBFGS,1,LBFGS_NPAR*NX*NY);
......@@ -488,6 +493,24 @@ int main(int argc, char **argv){
}
}
if (VTI) {
/* dynamic arrays for anisotropic modeling */
pc11 = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
pc13 = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
pc33 = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
pc55 = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
pc55ipjp = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
pc66 = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
pc66ipjp = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
pdelta = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
pepsilon = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
if(WAVETYPE==2 || WAVETYPE==3) {
pgamma = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
pvshornp1 = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
}
pvshor = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
}
if (L) {
/* dynamic (wavefield) arrays for viscoelastic modeling */
pr = f3tensor(-nd+1,NY+nd,-nd+1,NX+nd,1,L);
......@@ -550,7 +573,16 @@ int main(int argc, char **argv){
forward_prop_z_xz = f3tensor(-nd+1,NY+nd,-nd+1,NX+nd,1,NT/DTINV);
forward_prop_z_yz = f3tensor(-nd+1,NY+nd,-nd+1,NX+nd,1,NT/DTINV);
waveconv_rho_shot_z = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
waveconv_u_shot_z = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
if (VTI) {
waveconv_c55_shot_z = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
waveconv_c66_shot_z = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
waveconv_c55_s_z = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
waveconv_c66_s_z = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
waveconv_vshor_shot_z = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
waveconv_vshor_z = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
gradp_vshor_z = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
}
waveconv_u_shot_z = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
waveconv_mu_z = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
waveconv_rho_s_z = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
waveconv_u_z = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
......@@ -571,11 +603,17 @@ int main(int argc, char **argv){
c2_SL=WOLFE_C2_SL;
waveconv_old= matrix(-nd+1,NY+nd,-nd+1,NX+nd);
if(!ACOUSTIC) waveconv_u_old= matrix(-nd+1,NY+nd,-nd+1,NX+nd);
if(!ACOUSTIC){
waveconv_u_old= matrix(-nd+1,NY+nd,-nd+1,NX+nd);
waveconv_vshor_old= matrix(-nd+1,NY+nd,-nd+1,NX+nd);
}
waveconv_rho_old= matrix(-nd+1,NY+nd,-nd+1,NX+nd);
waveconv_up= matrix(-nd+1,NY+nd,-nd+1,NX+nd);
if(!ACOUSTIC) waveconv_u_up= matrix(-nd+1,NY+nd,-nd+1,NX+nd);
if(!ACOUSTIC) {
waveconv_u_up= matrix(-nd+1,NY+nd,-nd+1,NX+nd);
waveconv_vshor_up= matrix(-nd+1,NY+nd,-nd+1,NX+nd);
}
waveconv_rho_up= matrix(-nd+1,NY+nd,-nd+1,NX+nd);
}
......@@ -778,7 +816,12 @@ int main(int argc, char **argv){
if(L){
if(!ACOUSTIC){
if (READMOD){
readmod(prho,ppi,pu,ptaus,ptaup,peta);
if(VTI){
readmod_vti(prho,ppi,pu,pdelta,pepsilon,pvshor,ptaus,ptaup,peta);
}
else {
readmod(prho,ppi,pu,ptaus,ptaup,peta);
}
}else{
model(prho,ppi,pu,ptaus,ptaup,peta);
}
......@@ -792,9 +835,19 @@ int main(int argc, char **argv){
}else{
if(!ACOUSTIC){
if (READMOD){
readmod_elastic(prho,ppi,pu);
if (VTI) {
readmod_el_vti(prho,ppi,pu,pdelta,pepsilon,pvshor);
}
else {
readmod_elastic(prho,ppi,pu);
}
}else{
model_elastic(prho,ppi,pu);
if(VTI){
model_elastic_vti(prho,ppi,pu,pdelta,pepsilon,pvshor);
}
else {
model_elastic(prho,ppi,pu);
}
}
}else{
if (READMOD){
......@@ -806,7 +859,7 @@ int main(int argc, char **argv){
}
/* check if the FD run will be stable and free of numerical dispersion */
checkfd(FP, prho, ppi, pu, ptaus, ptaup, peta, hc, srcpos, nsrc, recpos, ntr_glob);
checkfd(FP, prho, ppi, pu, ptaus, ptaup, peta, hc, srcpos, nsrc, recpos, ntr_glob, pepsilon, pdelta, pvshor);
/* calculate damping coefficients for CPMLs*/
if(FW>0)
......@@ -899,15 +952,28 @@ int main(int argc, char **argv){
they have to be averaged. For this, values lying at 0 and NX+1,
for example, are required on the local grid. These are now copied from the
neighbouring grids */
if (VTI){
change_parameterization_vti(prho, ppi, pu, pepsilon, pdelta, pvshor, pc11, pc13, pc33, pc55, pc66);
}
if (L){
if(!ACOUSTIC){
matcopy(prho,ppi,pu,ptaus,ptaup);
if(VTI){
matcopy_vti(prho, pc11, pc13, pc33, pc55, pc66, ptaus, ptaup);
}
else {
matcopy(prho,ppi,pu,ptaus,ptaup);
}
} else {
matcopy_viscac(prho,ppi,ptaup);
}
}else{
if(!ACOUSTIC){
matcopy_elastic(prho, ppi, pu);
if (VTI) {
matcopy_elastic_vti(prho, pc11, pc13, pc33, pc55, pc66);
}
else {
matcopy_elastic(prho, ppi, pu);
}
}else{
matcopy_acoustic(prho, ppi);
}
......@@ -929,12 +995,18 @@ int main(int argc, char **argv){
if(!ACOUSTIC) av_mue(pu,puipjp,prho);
av_rho(prho,prip,prjp);
if (!ACOUSTIC && L) av_tau(ptaus,ptausipjp);
if (VTI) av_c55c66(pc55,pc55ipjp, pc66,pc66ipjp);
/* Preparing memory variables for update_s (viscoelastic) */
if (L) {
if(!ACOUSTIC){
prepare_update_s(etajm,etaip,peta,fipjp,pu,puipjp,ppi,prho,ptaus,ptaup,ptausipjp,f,g,bip,bjm,cip,cjm,dip,d,e);
if(VTI){
prepare_update_s_vti(etajm,etaip,peta,fipjp,pc55,pc66ipjp,ppi,prho,ptaus,ptaup,ptausipjp,f,g,bip,bjm,cip,cjm,dip,d,e);
}
else {
prepare_update_s(etajm,etaip,peta,fipjp,pu,puipjp,ppi,prho,ptaus,ptaup,ptausipjp,f,g,bip,bjm,cip,cjm,dip,d,e);
}
} else {
prepare_update_p(etajm,peta,ppi,prho,ptaup,g,bjm,cjm,e);
}
......@@ -978,13 +1050,24 @@ int main(int argc, char **argv){
Vp_avg=average_matrix(ppi);
rho_avg=average_matrix(prho);
if(!ACOUSTIC) Vs_avg=average_matrix(pu);
if (VTI) {
vshor_avg=average_matrix(pvshor);
}
if(!ACOUSTIC) if(VERBOSE) printf("MYID = %d \t Vp_avg = %e \t Vs_avg = %e \t rho_avg = %e \n ",MYID,Vp_avg,Vs_avg,rho_avg);
if(!ACOUSTIC) {
if (VTI) {
if (VERBOSE) printf("MYID = %d \t Vp_avg = %e \t Vs_avg = %e \t rho_avg = %e \t vshor_avg = %e \n ",MYID,Vp_avg,Vs_avg,rho_avg,vshor_avg);
}
else if(VERBOSE) printf("MYID = %d \t Vp_avg = %e \t Vs_avg = %e \t rho_avg = %e \n ",MYID,Vp_avg,Vs_avg,rho_avg);
}
else if(VERBOSE) printf("MYID = %d \t Vp_avg = %e \t rho_avg = %e \n ",MYID,Vp_avg,rho_avg);
C_vp = Vp_avg*Vp_avg;
if(!ACOUSTIC) C_vs = Vs_avg*Vs_avg;
C_rho = rho_avg*rho_avg;
if (VTI) {
C_vshor = vshor_avg*vshor_avg;
}
}
......@@ -1057,6 +1140,9 @@ int main(int argc, char **argv){
for (i=1;i<=NX;i=i+IDX){
waveconv_rho_z[j][i]=0.0;
waveconv_u_z[j][i]=0.0;
if (VTI) {
waveconv_vshor_z[j][i]=0.0;
}
}
}
......@@ -1251,18 +1337,33 @@ int main(int argc, char **argv){
}
}
if (WAVETYPE==2 || WAVETYPE==3) {
update_s_visc_PML_SH(1, NX, 1, NY, pvz, psxz, psyz, pt, po, bip, bjm, cip, cjm, d, dip,fipjp, f, hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy);
if(VTI){
update_s_visc_vti_PML_SH(1, NX, 1, NY, pvz, uxz, uyz, psxz, psyz, pt, po, bip, bjm, cip, cjm, d, dip,fipjp, f, hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy, pc55,pc66ipjp,prho);
}
else {
update_s_visc_PML_SH(1, NX, 1, NY, pvz, uxz, uyz, psxz, psyz, pt, po, bip, bjm, cip, cjm, d, dip,fipjp, f, hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy);
}
}
} else { /* elastic */
if (WAVETYPE==1 || WAVETYPE==3) {
if(!ACOUSTIC) {
update_s_elastic_PML(1, NX, 1, NY, pvx, pvy, ux, uy, uxy, uyx, psxx, psyy, psxy, ppi, pu, puipjp, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
if (VTI) {
update_s_el_vti_PML(1, NX, 1, NY, pvx, pvy, ux, uy, uxy, uyx, psxx, psyy, psxy, pc11, pc13, pc33, pc55ipjp, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
}
else {
update_s_elastic_PML(1, NX, 1, NY, pvx, pvy, ux, uy, uxy, uyx, psxx, psyy, psxy, ppi, pu, puipjp, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
}
} else {
update_p_PML(1, NX, 1, NY, pvx, pvy, psp, u, ppi, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
}
}
if (WAVETYPE==2 || WAVETYPE==3) {
update_s_elastic_PML_SH(1, NX, 1, NY, pvz,psxz,psyz,uxz,uyz,hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy,puipjp,pu,prho);
if (VTI) {
update_s_el_vti_PML_SH(1, NX, 1, NY, pvz,psxz,psyz,uxz,uyz,hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy,pc55,pc66ipjp,prho);
}
else {
update_s_elastic_PML_SH(1, NX, 1, NY, pvz,psxz,psyz,uxz,uyz,hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy,puipjp,pu,prho);
}
}
}
......@@ -1279,7 +1380,12 @@ int main(int argc, char **argv){
surface_PML(1, pvx, pvy, psxx, psyy, psxy,psyz, pp, pq, ppi, pu, prho, ptaup, ptaus, etajm, peta, hc, K_x, a_x, b_x, psi_vxxs, ux, uy,uxy,uyz,psxz,uxz);
}else{
/* elastic */
surface_elastic_PML(1, pvx, pvy, psxx, psyy, psxy,psyz, ppi, pu, prho, hc, K_x, a_x, b_x, psi_vxxs, ux, uy, uxy,uyz,psxz,uxz);
if (VTI) {
surface_el_vti_PML(1, pvx, pvy, psxx, psyy, psxy,psyz, pc11, pc13, pc33, prho, hc, K_x, a_x, b_x, psi_vxxs, ux, uy, uxy,uyz,psxz,uxz);
}
else {
surface_elastic_PML(1, pvx, pvy, psxx, psyy, psxy,psyz, ppi, pu, prho, hc, K_x, a_x, b_x, psi_vxxs, ux, uy, uxy,uyz,psxz,uxz);
}
}
} else {
/* viscoelastic and elastic ACOUSTIC */
......@@ -1588,12 +1694,23 @@ int main(int argc, char **argv){
/*initialize gradient matrices for each shot with zeros SH*/
if(WAVETYPE==2 || WAVETYPE==3){
if(FORWARD_ONLY==0){
for(j=1;j<=NY;j=j+IDY){
for(i=1;i<=NX;i=i+IDX){
waveconv_rho_shot_z[j][i]=0.0;
waveconv_u_shot_z[j][i]=0.0;
}
}
if (VTI) {
for(j=1;j<=NY;j=j+IDY){
for(i=1;i<=NX;i=i+IDX){
waveconv_rho_shot_z[j][i]=0.0;
waveconv_c55_shot_z[j][i]=0.0;
waveconv_c66_shot_z[j][i]=0.0;
}
}
}
else {
for(j=1;j<=NY;j=j+IDY){
for(i=1;i<=NX;i=i+IDX){
waveconv_rho_shot_z[j][i]=0.0;
waveconv_u_shot_z[j][i]=0.0;
}
}
}
}
}
......@@ -1694,18 +1811,33 @@ int main(int argc, char **argv){
}
}
if (WAVETYPE==2 || WAVETYPE==3) {
update_s_visc_PML_SH(1, NX, 1, NY, pvz, psxz, psyz, pt, po, bip, bjm, cip, cjm, d, dip,fipjp, f, hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy);
if(VTI){
update_s_visc_vti_PML_SH(1, NX, 1, NY, pvz, uxz, uyz, psxz, psyz, pt, po, bip, bjm, cip, cjm, d, dip,fipjp, f, hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy, pc55,pc66ipjp,prho);
}
else {
update_s_visc_PML_SH(1, NX, 1, NY, pvz, uxz, uyz, psxz, psyz, pt, po, bip, bjm, cip, cjm, d, dip,fipjp, f, hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy);
}
}
} else { /* elastic */
if (WAVETYPE==1 || WAVETYPE==3) {
if(!ACOUSTIC) {
update_s_elastic_PML(1, NX, 1, NY, pvx, pvy, ux, uy, uxy, uyx, psxx, psyy, psxy, ppi, pu, puipjp, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
if (VTI) {
update_s_el_vti_PML(1, NX, 1, NY, pvx, pvy, ux, uy, uxy, uyx, psxx, psyy, psxy, pc11, pc13, pc33, pc55ipjp, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
}
else {
update_s_elastic_PML(1, NX, 1, NY, pvx, pvy, ux, uy, uxy, uyx, psxx, psyy, psxy, ppi, pu, puipjp, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
}
} else {
update_p_PML(1, NX, 1, NY, pvx, pvy, psp, u, ppi, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
}
}
if (WAVETYPE==2 || WAVETYPE==3) {
update_s_elastic_PML_SH(1, NX, 1, NY, pvz,psxz,psyz,uxz,uyz,hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy,puipjp,pu,prho);
if (VTI) {
update_s_el_vti_PML_SH(1, NX, 1, NY, pvz,psxz,psyz,uxz,uyz,hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy,pc55,pc66ipjp,prho);
}
else {
update_s_elastic_PML_SH(1, NX, 1, NY, pvz,psxz,psyz,uxz,uyz,hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy,puipjp,pu,prho);
}
}
}
......@@ -1723,7 +1855,12 @@ int main(int argc, char **argv){
surface_PML(1, pvx, pvy, psxx, psyy, psxy,psyz, pp, pq, ppi, pu, prho, ptaup, ptaus, etajm, peta, hc, K_x, a_x, b_x, psi_vxxs, ux, uy,uxy,uyz,psxz,uxz);
}else{
/* elastic */
surface_elastic_PML(1, pvx, pvy, psxx, psyy, psxy,psyz, ppi, pu, prho, hc, K_x, a_x, b_x, psi_vxxs, ux, uy, uxy,uyz,psxz,uxz);
if (VTI) {
surface_el_vti_PML(1, pvx, pvy, psxx, psyy, psxy,psyz, pc11, pc13, pc33, prho, hc, K_x, a_x, b_x, psi_vxxs, ux, uy, uxy,uyz,psxz,uxz);
}
else {
surface_elastic_PML(1, pvx, pvy, psxx, psyy, psxy,psyz, ppi, pu, prho, hc, K_x, a_x, b_x, psi_vxxs, ux, uy, uxy,uyz,psxz,uxz);
}
}
} else {
/* viscoelastic and elastic ACOUSTIC */
......@@ -2201,16 +2338,31 @@ int main(int argc, char **argv){
}
}
if (WAVETYPE==2 || WAVETYPE==3) {
update_s_visc_PML_SH(1, NX, 1, NY, pvz, psxz, psyz, pt, po, bip, bjm, cip, cjm, d, dip,fipjp, f, hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy);
if(VTI){
update_s_visc_vti_PML_SH(1, NX, 1, NY, pvz, uxz, uyz, psxz, psyz, pt, po, bip, bjm, cip, cjm, d, dip,fipjp, f, hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy, pc55,pc66ipjp,prho);
}
else {
update_s_visc_PML_SH(1, NX, 1, NY, pvz, uxz, uyz, psxz, psyz, pt, po, bip, bjm, cip, cjm, d, dip,fipjp, f, hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy);
}
}
} else{
/* elastic */
if(!ACOUSTIC){
if (WAVETYPE==1 || WAVETYPE==3) {
update_s_elastic_PML(1, NX, 1, NY, pvx, pvy, ux, uy, uxy, uyx, psxx, psyy, psxy, ppi, pu, puipjp, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
if (VTI) {
update_s_el_vti_PML(1, NX, 1, NY, pvx, pvy, ux, uy, uxy, uyx, psxx, psyy, psxy, pc11, pc13, pc33, pc55ipjp, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
}
else {
update_s_elastic_PML(1, NX, 1, NY, pvx, pvy, ux, uy, uxy, uyx, psxx, psyy, psxy, ppi, pu, puipjp, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
}
}
if (WAVETYPE==2 || WAVETYPE==3) {
update_s_elastic_PML_SH(1, NX, 1, NY, pvz,psxz,psyz,uxz,uyz,hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy,puipjp,pu,prho);
if (VTI) {
update_s_el_vti_PML_SH(1, NX, 1, NY, pvz,psxz,psyz,uxz,uyz,hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy,pc55,pc66ipjp,prho);
}
else {
update_s_elastic_PML_SH(1, NX, 1, NY, pvz,psxz,psyz,uxz,uyz,hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy,puipjp,pu,prho);
}
}
} else {
update_p_PML(1, NX, 1, NY, pvx, pvy, psp, u, ppi, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
......@@ -2229,7 +2381,12 @@ int main(int argc, char **argv){
surface_PML(1, pvx, pvy, psxx, psyy, psxy,psyz, pp, pq, ppi, pu, prho, ptaup, ptaus, etajm, peta, hc, K_x, a_x, b_x, psi_vxxs, ux, uy,uxy,uyz,psxz,uxz);
}else{
/* elastic */
surface_elastic_PML(1, pvx, pvy, psxx, psyy, psxy,psyz, ppi, pu, prho, hc, K_x, a_x, b_x, psi_vxxs, ux, uy, uxy,uyz,psxz,uxz);
if (VTI) {
surface_el_vti_PML(1, pvx, pvy, psxx, psyy, psxy,psyz, pc11, pc13, pc33, prho, hc, K_x, a_x, b_x, psi_vxxs, ux, uy, uxy,uyz,psxz,uxz);
}
else {
surface_elastic_PML(1, pvx, pvy, psxx, psyy, psxy,psyz, ppi, pu, prho, hc, K_x, a_x, b_x, psi_vxxs, ux, uy, uxy,uyz,psxz,uxz);
}
}
} else {
/* viscoelastic and elastic ACOUSTIC */
......@@ -2292,9 +2449,15 @@ int main(int argc, char **argv){
}
if (WAVETYPE==2 || WAVETYPE==3) {
waveconv_rho_shot_z[j][i]+=(pvzp1[j][i]*forward_prop_rho_z[j][i][NTDTINV-hin+1]);
muss = prho[j][i] * pu[j][i] * pu[j][i];
waveconv_u_shot_z[j][i]+=(1.0/(muss*muss))*(forward_prop_z_xz[j][i][NTDTINV-hin+1]*psxz[j][i]+forward_prop_z_yz[j][i][NTDTINV-hin+1]*psyz[j][i]);
waveconv_rho_shot_z[j][i]+=(pvzp1[j][i]*forward_prop_rho_z[j][i][NTDTINV-hin+1]);
if (VTI) {
waveconv_c55_shot_z[j][i] += (1.0/(pc55[j][i]*pc55[j][i]))*(forward_prop_z_yz[j][i][NTDTINV-hin+1]*psyz[j][i]);
waveconv_c66_shot_z[j][i] += (1.0/(pc66[j][i]*pc66[j][i]))*(forward_prop_z_xz[j][i][NTDTINV-hin+1]*psxz[j][i]);
}
else {
muss = prho[j][i] * pu[j][i] * pu[j][i];
waveconv_u_shot_z[j][i]+=(1.0/(muss*muss))*(forward_prop_z_xz[j][i][NTDTINV-hin+1]*psxz[j][i]+forward_prop_z_yz[j][i][NTDTINV-hin+1]*psyz[j][i]);
}
}
}
}
......@@ -2439,24 +2602,37 @@ int main(int argc, char **argv){
/* interpolate unknown values */
if((IDXI>1)||(IDYI>1)){
interpol(IDXI,IDYI,waveconv_u_shot_z,1);
if (VTI) {
interpol(IDXI,IDYI,waveconv_c55_shot_z,1);
interpol(IDXI,IDYI,waveconv_c66_shot_z,1);
}
}
/* calculate complete gradient */
for (j=1;j<=NY;j=j+IDY){
for (i=1;i<=NX;i=i+IDX){
/* calculate mu gradient */
waveconv_mu_z[j][i] = - DT * waveconv_u_shot_z[j][i];
if(PARAMETERIZATION==1){
/* calculate Vs gradient */
waveconv_u_shot_z[j][i] = (2.0 * prho[j][i] * pu[j][i] * waveconv_mu_z[j][i]);
}
if(PARAMETERIZATION==3){
/* calculate u gradient */
waveconv_u_shot_z[j][i] = waveconv_mu_z[j][i];
}
if (VTI) {
waveconv_c55_s_z[j][i] = -DT * waveconv_c55_shot_z[j][i];
waveconv_c66_s_z[j][i] = -DT * waveconv_c66_shot_z[j][i];
waveconv_u_shot_z[j][i] = 2.0*prho[j][i]*pu[j][i]*waveconv_c55_s_z[j][i];
}
else {
/* calculate mu gradient */
waveconv_mu_z[j][i] = - DT * waveconv_u_shot_z[j][i];
if(PARAMETERIZATION==1){
/* calculate Vs gradient */
waveconv_u_shot_z[j][i] = (2.0 * prho[j][i] * pu[j][i] * waveconv_mu_z[j][i]);
}
if(PARAMETERIZATION==3){
/* calculate u gradient */
waveconv_u_shot_z[j][i] = waveconv_mu_z[j][i];
}
}
}
}
}
......@@ -2515,22 +2691,45 @@ int main(int argc, char **argv){
/* calculate density gradient rho' */
waveconv_rho_s_z[j][i]= - DT * waveconv_rho_shot_z[j][i];
if(PARAMETERIZATION==1){
/* calculate density gradient */
waveconv_rho_shot_z[j][i] = ( (pu[j][i] * pu[j][i] * waveconv_mu_z[j][i]) + waveconv_rho_s_z[j][i]);
}
if(PARAMETERIZATION==3){
/* calculate density gradient */
waveconv_rho_shot_z[j][i] = waveconv_rho_s_z[j][i];
}
if (VTI) {
waveconv_rho_shot_z[j][i] = pu[j][i]*pu[j][i]*waveconv_c55_s_z[j][i] + pvshor[j][i]*pvshor[j][i]*waveconv_c66_s_z[j][i] + waveconv_rho_s_z[j][i];
}
else {
if(PARAMETERIZATION==1){
/* calculate density gradient */
waveconv_rho_shot_z[j][i] = ( (pu[j][i] * pu[j][i] * waveconv_mu_z[j][i]) + waveconv_rho_s_z[j][i]);
}
if(PARAMETERIZATION==3){
/* calculate density gradient */
waveconv_rho_shot_z[j][i] = waveconv_rho_s_z[j][i];
}
}
}
}
}
/* -------------------------------------------- */
/* calculate gradient direction horizontal Vs (VTI, SH) */
/* -------------------------------------------- */
if (VTI) {
if((WAVETYPE==2 || WAVETYPE==3)&&(FORWARD_ONLY==0)){
for (j=1;j<=NY;j=j+IDY){
for (i=1;i<=NX;i=i+IDX){
waveconv_vshor_shot_z[j][i] = 2.0*prho[j][i]*pvshor[j][i]*waveconv_c66_s_z[j][i];
}
}
}
}
/* -------------------------------------------- */
/* calculate and apply energy preconditioning */
/* -------------------------------------------- */
......@@ -2566,11 +2765,16 @@ int main(int argc, char **argv){
if(EPRECOND_PER_SHOT_SH && (WAVETYPE==2 || WAVETYPE==3)){
fprintf(FP,"\n Applying approx. Hessian for shot %i SH. EPRECOND=%i, EPSILON_WE=%f",ishot,EPRECOND,EPSILON_WE);
for (j=1;j<=NY;j=j+IDY){
for (i=1;i<=NX;i=i+IDX){
We_SH[j][i]=We_SH[j][i]/We_max_SH;
waveconv_u_shot_z[j][i] = waveconv_u_shot_z[j][i]/(We_SH[j][i]);
waveconv_rho_shot_z[j][i] = waveconv_rho_shot_z[j][i]/(We_SH[j][i]);
if (VTI) {
waveconv_vshor_shot_z[j][i] = waveconv_vshor_shot_z[j][i]/(We_SH[j][i]);
}
}
}
}
......@@ -2594,7 +2798,9 @@ int main(int argc, char **argv){
/* applying the preconditioning */
taper_grad_shot(waveconv_u_shot_z,taper_coeff,srcpos,nsrc,recpos,ntr_glob,ishot,1);
taper_grad_shot(waveconv_rho_shot_z,taper_coeff,srcpos,nsrc,recpos,ntr_glob,ishot,1);
if (VTI) {
taper_grad_shot(waveconv_vshor_shot_z,taper_coeff,srcpos,nsrc,recpos,ntr_glob,ishot,1);
}
}
/* end of SWS_TAPER_CIRCULAR_PER_SHOT == 1 */
......@@ -2604,6 +2810,7 @@ int main(int argc, char **argv){
/* applying the preconditioning */
taper_grad_shot(waveconv_u_shot_z,taper_coeff,srcpos,nsrc,recpos,ntr_glob,ishot,3); /*taper vs gradient */
taper_grad_shot(waveconv_rho_shot_z,taper_coeff,srcpos,nsrc,recpos,ntr_glob,ishot,4); /*taper rho gradient */
taper_grad_shot(waveconv_vshor_shot_z,taper_coeff,srcpos,nsrc,recpos,ntr_glob,ishot,5); /*taper vshor gradient */
}
}
......@@ -2657,12 +2864,25 @@ int main(int argc, char **argv){
/* Summing up the gradient for all shots SH */
/* ----------------------------------------- */
if (WAVETYPE==2 || WAVETYPE==3) {
for (j=1;j<=NY;j=j+IDY){
for (i=1;i<=NX;i=i+IDX){
waveconv_u_z[j][i] += waveconv_u_shot_z[j][i];
waveconv_rho_z[j][i] += waveconv_rho_shot_z[j][i];
}
}
if (VTI) {
for (j=1;j<=NY;j=j+IDY){
for (i=1;i<=NX;i=i+IDX){
waveconv_u_z[j][i] += waveconv_u_shot_z[j][i];
waveconv_rho_z[j][i] += waveconv_rho_shot_z[j][i];
waveconv_vshor_z[j][i] += waveconv_vshor_shot_z[j][i];
}
}
}
else{
for (j=1;j<=NY;j=j+IDY){
for (i=1;i<=NX;i=i+IDX){
waveconv_u_z[j][i] += waveconv_u_shot_z[j][i];
waveconv_rho_z[j][i] += waveconv_rho_shot_z[j][i];
}
}
}
}
}
......@@ -2699,6 +2919,7 @@ int main(int argc, char **argv){
}
}
if(!EPRECOND_PER_SHOT_SH && (WAVETYPE==2 || WAVETYPE==3)){
fprintf(FP,"\n Applying approx. Hessian to summed gradient SH. EPRECOND=%i, EPSILON_WE=%f",EPRECOND,EPSILON_WE);
......@@ -2708,6 +2929,9 @@ int main(int argc, char **argv){
We_sum_SH[j][i]=We_sum_SH[j][i]/We_sum_max1;
waveconv_u_z[j][i] = waveconv_u_z[j][i]*We_sum_SH[j][i];
waveconv_rho_z[j][i] = waveconv_rho_z[j][i]*We_sum_SH[j][i];
if (VTI) {
waveconv_vshor_z[j][i] = waveconv_vshor_z[j][i]*We_sum_SH[j][i];
}
}
}
}
......@@ -2740,6 +2964,9 @@ int main(int argc, char **argv){
for (i=1;i<=NX;i=i+IDX){
if(iter<INV_VS_ITER) waveconv_u_z[j][i] = 0.0;
if(iter<INV_RHO_ITER) waveconv_rho_z[j][i] = 0.0;
if (VTI) {
if(iter<INV_VSHOR_ITER) waveconv_vshor_z[j][i] = 0.0;
}
}
}
}
......@@ -2825,20 +3052,33 @@ int main(int argc, char **argv){
printf("sum_killed_traces_testshots=%d\n",sum_killed_traces_testshots);
printf("sum_killed_traces=%d\n",sum_killed_traces);}}
else{
L2t[1]=2.0*(1.0+(L2sum/((float)(NO_OF_TESTSHOTS*ntr_glob-sum_killed_traces_testshots))));
L2t[4]=2.0*(1.0+(L2sum_all_shots/((float)(nsrc_glob*ntr_glob-sum_killed_traces))));
if (MYID==0){
printf("sum_killed_traces_testshots=%d\n",sum_killed_traces_testshots);
printf("sum_killed_traces=%d\n",sum_killed_traces);
printf("ntr_glob=%d\n",ntr_glob);
printf("nsrc_glob=%d\n",nsrc_glob);}}}
L2t[1]=2.0*(1.0+(L2sum/((float)(NO_OF_TESTSHOTS*ntr_glob-sum_killed_traces_testshots))));
L2t[4]=2.0*(1.0+(L2sum_all_shots/((float)(nsrc_glob*ntr_glob-sum_killed_traces))));
if (MYID==0){
printf("sum_killed_traces_testshots=%d\n",sum_killed_traces_testshots);
printf("sum_killed_traces=%d\n",sum_killed_traces);
printf("ntr_glob=%d\n",ntr_glob);
printf("nsrc_glob=%d\n",nsrc_glob);}
}
if(WAVETYPE==2){
L2t[1]=2.0*(1.0+(L2sum_SH/((float)((NO_OF_TESTSHOTS*ntr_glob-sum_killed_traces_testshots)))));
L2t[4]=2.0*(1.0+(L2sum_all_shots_SH/((float)((nsrc_glob*ntr_glob-sum_killed_traces)))));
}
}
else{
if(ADJOINT_TYPE==1){ /* x and y component are used in the inversion */
L2t[1]=2.0*(1.0+(L2sum/((float)NO_OF_TESTSHOTS*(float)ntr_glob*2.0)));
L2t[4]=2.0*(1.0+(L2sum_all_shots/((float)nsrc_glob*(float)ntr_glob*2.0)));}
else{
L2t[1]=2.0*(1.0+(L2sum/((float)NO_OF_TESTSHOTS*(float)ntr_glob)));
L2t[4]=2.0*(1.0+(L2sum_all_shots/((float)nsrc_glob*(float)ntr_glob)));}
if (WAVETYPE==1 || WAVETYPE==3) {
L2t[1]=2.0*(1.0+(L2sum/((float)NO_OF_TESTSHOTS*(float)ntr_glob)));
L2t[4]=2.0*(1.0+(L2sum_all_shots/((float)nsrc_glob*(float)ntr_glob)));
}
}
if (WAVETYPE==2) {
L2t[1]=2.0*(1.0+(L2sum_SH/((float)NO_OF_TESTSHOTS*(float)ntr_glob)));
L2t[4]=2.0*(1.0+(L2sum_all_shots_SH/((float)nsrc_glob*(float)ntr_glob)));
}
}
break;
case 8:
......@@ -2846,8 +3086,14 @@ int main(int argc, char **argv){
L2t[4]=L2sum_all_shots/energy_sum_all_shots;
break;
default:
L2t[1]=L2sum;
L2t[4]=L2sum_all_shots;
if ((WAVETYPE==1) || (WAVETYPE==3)) {
L2t[1]=L2sum;
L2t[4]=L2sum_all_shots;
}
else {
L2t[1]=L2sum_SH;
L2t[4]=L2sum_all_shots_SH;
}
break;
}
......@@ -2888,9 +3134,16 @@ int main(int argc, char **argv){
}
}
if (WAVETYPE==1 || WAVETYPE==3) PCG(waveconv, taper_coeff, nsrc, srcpos, recpos, ntr_glob, iter, C_vp, gradp, nfstart_jac, waveconv_u, C_vs, gradp_u, waveconv_rho, C_rho, gradp_rho,Vs_avg,F_LOW_PASS,PCG_iter_start);
if (WAVETYPE==2 || WAVETYPE==3) PCG_SH(taper_coeff, nsrc, srcpos, recpos, ntr_glob, iter, nfstart_jac, waveconv_u_z, C_vs, gradp_u_z, waveconv_rho_z, C_rho, gradp_rho_z,Vs_avg,F_LOW_PASS,PCG_iter_start);
if (WAVETYPE==1 || WAVETYPE==3) PCG(waveconv, taper_coeff, nsrc, srcpos, recpos, ntr_glob, iter, C_vp, gradp, nfstart_jac, waveconv_u, C_vs, gradp_u, waveconv_rho, C_rho, gradp_rho,Vs_avg,F_LOW_PASS,PCG_iter_start);
if (WAVETYPE==2 || WAVETYPE==3) {
if (VTI) {
PCG_SH_vti(taper_coeff, nsrc, srcpos, recpos, ntr_glob, iter, nfstart_jac, waveconv_u_z, C_vs, gradp_u_z, waveconv_rho_z, C_rho, gradp_rho_z,waveconv_vshor_z,C_vshor,gradp_vshor_z,Vs_avg,F_LOW_PASS,PCG_iter_start);
}
else {
PCG_SH(taper_coeff, nsrc, srcpos, recpos, ntr_glob, iter, nfstart_jac, waveconv_u_z, C_vs, gradp_u_z, waveconv_rho_z, C_rho, gradp_rho_z,Vs_avg,F_LOW_PASS,PCG_iter_start);
}
}
}
/* ---------------------------------------------------------*/
......@@ -2958,6 +3211,9 @@ int main(int argc, char **argv){
if (SWS_TAPER_FILE){ /* read taper from BIN-File*/
taper_grad(waveconv_rho,taper_coeff,srcpos,nsrc,recpos,ntr_glob,6);}
if (SWS_TAPER_FILE && VTI && WAVETYPE==2){
taper_grad(waveconv_vshor_z,taper_coeff,srcpos,nsrc,recpos,ntr_glob,7);}
if(GRAD_FILTER==1 && !ACOUSTIC){smooth(waveconv_u,2,1,Vs_avg,F_LOW_PASS);}
if(GRAD_FILTER==1 && !ACOUSTIC){
......@@ -2965,6 +3221,9 @@ int main(int argc, char **argv){
}else if(GRAD_FILTER==1 && ACOUSTIC){
smooth(waveconv_rho,3,1,Vp_avg,F_LOW_PASS);
}
if(GRAD_FILTER==1 && VTI && WAVETYPE==2){
smooth(waveconv_vshor_z,7,1,vshor_avg,F_LOW_PASS);
}
if(WOLFE_CONDITION) {
for (j=1;j<=NY;j=j+IDY){
......@@ -2980,6 +3239,11 @@ int main(int argc, char **argv){
punp1[j][i] =pu[j][i];
}
if(VTI && WAVETYPE==2){
waveconv_vshor_old[j][i] = waveconv_vshor_z[j][i];
pvshornp1[j][i] = pvshor[j][i];
}
waveconv_rho_old[j][i] = waveconv_rho[j][i];
prhonp1[j][i] = prho[j][i];
}
......@@ -2989,7 +3253,7 @@ int main(int argc, char **argv){
/*---------------------*/
/* L-BFGS */
/*---------------------*/
lbfgs(waveconv_u, waveconv_rho, waveconv,Vs_avg,rho_avg,Vp_avg, rho_LBFGS, s_LBFGS, y_LBFGS, N_LBFGS,LBFGS_NPAR, iter,&LBFGS_iter_start);
lbfgs(waveconv_u, waveconv_rho, waveconv, waveconv_vshor_z, Vs_avg,rho_avg,Vp_avg, vshor_avg, rho_LBFGS, s_LBFGS, y_LBFGS, N_LBFGS,LBFGS_NPAR, iter,&LBFGS_iter_start);
}
......@@ -2997,6 +3261,7 @@ int main(int argc, char **argv){
for (j=1;j<=NY;j=j+IDY){
for (i=1;i<=NX;i=i+IDX){
if(WAVETYPE!=2) waveconv_up[j][i]=waveconv[j][i];
else if (VTI) waveconv_vshor_up[j][i]=waveconv_vshor_z[j][i];
if(!ACOUSTIC) waveconv_u_up[j][i] = waveconv_u[j][i];
waveconv_rho_up[j][i] = waveconv_rho[j][i];
}
......@@ -3039,14 +3304,19 @@ int main(int argc, char **argv){
if(iter>2 && WOLFE_TRY_OLD_STEPLENGTH) alpha_SL=alpha_SL_old;
/* Calculate update */
calc_mat_change_test(waveconv_up,waveconv_rho_up,waveconv_u_up,prhonp1,prho,ppinp1,ppi,punp1,pu,iter,1,FORWARD_ONLY,alpha_SL,1,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
FWI_run=1;
if (VTI) {
calc_mat_change_test_vti(waveconv_up,waveconv_rho_up,waveconv_u_up,waveconv_vshor_up,prhonp1,prho,ppinp1,ppi,punp1,pu,pvshornp1,pvshor,iter,1,FORWARD_ONLY,alpha_SL,1,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,vshor_avg,LBFGS_iter_start);
}
else {
calc_mat_change_test(waveconv_up,waveconv_rho_up,waveconv_u_up,prhonp1,prho,ppinp1,ppi,punp1,pu,iter,1,FORWARD_ONLY,alpha_SL,1,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
}
FWI_run=1;
}
/* Check if current step length satisfy wolfe condition, if not call linesearch for new step length */
if(countstep>0) {
wolfe_status=check_wolfe(alpha_SL, L2_SL_old, L2_SL_new, waveconv_u_old, waveconv_u, waveconv_u_up,waveconv_rho_old, waveconv_rho, waveconv_rho_up,waveconv_old, waveconv, waveconv_up, c1_SL, c2_SL,LBFGS_NPAR);
wolfe_status=check_wolfe(alpha_SL, L2_SL_old, L2_SL_new, waveconv_u_old, waveconv_u, waveconv_u_up,waveconv_rho_old, waveconv_rho, waveconv_rho_up,waveconv_old, waveconv, waveconv_up, waveconv_vshor_old, waveconv_vshor_z, waveconv_vshor_up, c1_SL, c2_SL,LBFGS_NPAR);
if(wolfe_status==0) {
/* Current step length satisfy wolfe condition, abort step length search */
......@@ -3067,8 +3337,13 @@ int main(int argc, char **argv){
/* Current step length do not satisfy wolfe condition, try new step length */
wolfe_linesearch(wolfe_status, &alpha_SL_min, &alpha_SL_max, &alpha_SL);
calc_mat_change_test(waveconv_up,waveconv_rho_up,waveconv_u_up,prhonp1,prho,ppinp1,ppi,punp1,pu,iter,1,FORWARD_ONLY,alpha_SL,1,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
fprintf(FP,"; New steplength=%1.3f",alpha_SL);
if (VTI) {
calc_mat_change_test_vti(waveconv_up,waveconv_rho_up,waveconv_u_up,waveconv_vshor_up,prhonp1,prho,ppinp1,ppi,punp1,pu,pvshornp1,pvshor,iter,1,FORWARD_ONLY,alpha_SL,1,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,vshor_avg,LBFGS_iter_start);
}
else {
calc_mat_change_test(waveconv_up,waveconv_rho_up,waveconv_u_up,prhonp1,prho,ppinp1,ppi,punp1,pu,iter,1,FORWARD_ONLY,alpha_SL,1,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
}
fprintf(FP,"; New steplength=%1.3f",alpha_SL);
FWI_run=1;
}
......@@ -3097,8 +3372,12 @@ int main(int argc, char **argv){
/* No update is done here, however model fils are written to disk for easy post processing */
alpha_SL=0.0;
calc_mat_change_test(waveconv_up,waveconv_rho_up,waveconv_u_up,prhonp1,prho,ppinp1,ppi,punp1,pu,iter,1,FORWARD_ONLY,alpha_SL,0,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
if (VTI) {
calc_mat_change_test_vti(waveconv_up,waveconv_rho_up,waveconv_u_up,waveconv_vshor_up,prhonp1,prho,ppinp1,ppi,punp1,pu,pvshornp1,pvshor,iter,1,FORWARD_ONLY,alpha_SL,0,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,vshor_avg,LBFGS_iter_start);
}
else {
calc_mat_change_test(waveconv_up,waveconv_rho_up,waveconv_u_up,prhonp1,prho,ppinp1,ppi,punp1,pu,iter,1,FORWARD_ONLY,alpha_SL,0,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
}
alpha_SL_old=1;
/* If minimum number of iterations would be enforced, L-BFGS is likely to crash */
......@@ -3118,13 +3397,19 @@ int main(int argc, char **argv){
for (i=1;i<=NX;i=i+IDX){
prho[j][i]=prhonp1[j][i];
if(WAVETYPE!=2) ppi[j][i]=ppinp1[j][i];
else if (VTI) pvshor[j][i]=pvshornp1[j][i];
if(!ACOUSTIC) pu[j][i]=punp1[j][i];
}
}
/* do the final model update */
calc_mat_change_test(waveconv_up,waveconv_rho_up,waveconv_u_up,prhonp1,prho,ppinp1,ppi,punp1,pu,iter,1,FORWARD_ONLY,alpha_SL,0,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
if (VTI) {
calc_mat_change_test_vti(waveconv_up,waveconv_rho_up,waveconv_u_up,waveconv_vshor_up,prhonp1,prho,ppinp1,ppi,punp1,pu,pvshornp1,pvshor,iter,1,FORWARD_ONLY,alpha_SL,0,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,vshor_avg,LBFGS_iter_start);
}
else{
calc_mat_change_test(waveconv_up,waveconv_rho_up,waveconv_u_up,prhonp1,prho,ppinp1,ppi,punp1,pu,iter,1,FORWARD_ONLY,alpha_SL,0,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
}
L2_hist[iter]=L2t[4];
/* write L2 log file */
......@@ -3211,13 +3496,21 @@ int main(int argc, char **argv){
for (itest=itests;itest<=iteste;itest++){
/* calculate change in the material parameters */
calc_mat_change_test(waveconv,waveconv_rho,waveconv_u,prho,prhonp1,ppi,ppinp1,pu,punp1,iter,1,FORWARD_ONLY,eps_scale,1,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
if (VTI) {
calc_mat_change_test_vti(waveconv,waveconv_rho,waveconv_u,waveconv_vshor_z,prho,prhonp1,ppi,ppinp1,pu,punp1,pvshor,pvshornp1,iter,1,FORWARD_ONLY,eps_scale,1,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,vshor_avg,LBFGS_iter_start);
}
else {
calc_mat_change_test(waveconv,waveconv_rho,waveconv_u,prho,prhonp1,ppi,ppinp1,pu,punp1,iter,1,FORWARD_ONLY,eps_scale,1,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
}
/* For the calculation of the material parameters beteween gridpoints
the have to be averaged. For this, values lying at 0 and NX+1,
for example, are required on the local grid. These are now copied from the
neighbouring grids */
if(!ACOUSTIC) {
if (VTI) {
change_parameterization_vti(prhonp1,ppinp1,punp1,pepsilon,pdelta,pvshornp1,pc11,pc13,pc33,pc55,pc66 );
matcopy_elastic_vti(prhonp1,pc11,pc13,pc33,pc55,pc66);
}
matcopy_elastic(prhonp1, ppinp1, punp1); /* no differentiation of elastic and viscoelastic modelling because the viscoelastic parameters did not change during the forward modelling */
}else{
matcopy_acoustic(prhonp1, ppinp1);
......@@ -3227,11 +3520,19 @@ int main(int argc, char **argv){
if(!ACOUSTIC) av_mue(punp1,puipjp,prhonp1);
av_rho(prhonp1,prip,prjp);
if (VTI) {
av_c55c66(pc55,pc55ipjp,pc66,pc66ipjp);
}
/* Preparing memory variables for update_s (viscoelastic) */
if (L) {
if(!ACOUSTIC) {
prepare_update_s(etajm,etaip,peta,fipjp,punp1,puipjp,ppinp1,prhonp1,ptaus,ptaup,ptausipjp,f,g, bip,bjm,cip,cjm,dip,d,e);
if(VTI){
prepare_update_s_vti(etajm,etaip,peta,fipjp,pc55,pc66ipjp,ppinp1,prhonp1,ptaus,ptaup,ptausipjp,f,g,bip,bjm,cip,cjm,dip,d,e);
}
else {
prepare_update_s(etajm,etaip,peta,fipjp,punp1,puipjp,ppinp1,prhonp1,ptaus,ptaup,ptausipjp,f,g, bip,bjm,cip,cjm,dip,d,e);
}
}else{
prepare_update_p(etajm,peta,ppinp1,prhonp1,ptaup,g,bjm,cjm,e);
}
......@@ -3416,7 +3717,12 @@ int main(int argc, char **argv){
}
}
if (WAVETYPE==2 || WAVETYPE==3) {
update_s_visc_PML_SH(1, NX, 1, NY, pvz, psxz, psyz, pt, po, bip, bjm, cip, cjm, d, dip,fipjp, f, hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy);
if(VTI){
update_s_visc_vti_PML_SH(1, NX, 1, NY, pvz, uxz, uyz, psxz, psyz, pt, po, bip, bjm, cip, cjm, d, dip,fipjp, f, hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy, pc55,pc66ipjp,prho);
}
else {
update_s_visc_PML_SH(1, NX, 1, NY, pvz, uxz, uyz, psxz, psyz, pt, po, bip, bjm, cip, cjm, d, dip,fipjp, f, hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy);
}
}
} else { /* elastic */
if(!ACOUSTIC){
......@@ -3424,8 +3730,13 @@ int main(int argc, char **argv){
update_s_elastic_PML(1, NX, 1, NY, pvx, pvy, ux, uy, uxy, uyx, psxx, psyy, psxy, ppinp1, punp1, puipjp, absorb_coeff, prhonp1, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
}
if (WAVETYPE==2 || WAVETYPE==3) {
update_s_elastic_PML_SH(1, NX, 1, NY, pvz,psxz,psyz,uxz,uyz,hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy,puipjp,punp1,prhonp1);
}
if (VTI) {
update_s_el_vti_PML_SH(1, NX, 1, NY, pvz,psxz,psyz,uxz,uyz,hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy,pc55,pc66ipjp,prhonp1);
}
else {
update_s_elastic_PML_SH(1, NX, 1, NY, pvz,psxz,psyz,uxz,uyz,hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy,puipjp,punp1,prhonp1);
}
}
}
else /* acoustic */
update_p_PML(1, NX, 1, NY, pvx, pvy, psp, u, ppinp1, absorb_coeff, prhonp1, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
......@@ -3449,8 +3760,13 @@ int main(int argc, char **argv){
surface_PML(1, pvx, pvy, psxx, psyy, psxy,psyz, pp, pq, ppinp1, punp1, prhonp1, ptaup, ptaus, etajm, peta, hc, K_x, a_x, b_x, psi_vxxs, ux, uy,uxy,uyz,psxz,uxz);
}else{
/* elastic */
surface_elastic_PML(1, pvx, pvy, psxx, psyy, psxy,psyz, ppinp1, punp1, prhonp1, hc, K_x, a_x, b_x, psi_vxxs, ux, uy, uxy,uyz,psxz,uxz);
}
if (VTI) {
surface_el_vti_PML(1, pvx, pvy, psxx, psyy, psxy,psyz, pc11, pc13, pc33, prhonp1, hc, K_x, a_x, b_x, psi_vxxs, ux, uy, uxy,uyz,psxz,uxz);
}
else {
surface_elastic_PML(1, pvx, pvy, psxx, psyy, psxy,psyz, ppinp1, punp1, prhonp1, hc, K_x, a_x, b_x, psi_vxxs, ux, uy, uxy,uyz,psxz,uxz);
}
}
} else {
/* viscoelastic and elastic ACOUSTIC */
surface_acoustic_PML(1, psp);
......@@ -3621,25 +3937,42 @@ int main(int argc, char **argv){
if(ADJOINT_TYPE==1){ /* x and y component are used in the inversion */
L2t[itest]=2.0*(1.0+(L2sum/((float)((NO_OF_TESTSHOTS*ntr_glob-sum_killed_traces_testshots)*2.0))));
if (MYID==0){
printf("sum_killed_traces_testshots=%d\n",sum_killed_traces_testshots);}}
printf("sum_killed_traces_testshots=%d\n",sum_killed_traces_testshots);}
}
else{
L2t[itest]=2.0*(1.0+(L2sum/((float)(NO_OF_TESTSHOTS*ntr_glob-sum_killed_traces_testshots))));
if (MYID==0){
printf("sum_killed_traces_testshots=%d\n",sum_killed_traces_testshots);
printf("ntr_glob=%d\n",ntr_glob);
printf("NO_OF_TESTSHOTS=%d\n",NO_OF_TESTSHOTS);}}}
L2t[itest]=2.0*(1.0+(L2sum/((float)(NO_OF_TESTSHOTS*ntr_glob-sum_killed_traces_testshots))));
if (MYID==0){
printf("sum_killed_traces_testshots=%d\n",sum_killed_traces_testshots);
printf("ntr_glob=%d\n",ntr_glob);
printf("NO_OF_TESTSHOTS=%d\n",NO_OF_TESTSHOTS);}
}
if(WAVETYPE==2){
L2t[itest]=2.0*(1.0+(L2sum_SH/((float)((NO_OF_TESTSHOTS*ntr_glob-sum_killed_traces_testshots)))));
}
}
else{
if(ADJOINT_TYPE==1){ /* x and y component are used in the inversion */
L2t[itest]=2.0*(1.0+(L2sum/((float)NO_OF_TESTSHOTS*(float)ntr_glob*2.0)));}
else{
L2t[itest]=2.0*(1.0+(L2sum/((float)NO_OF_TESTSHOTS*(float)ntr_glob)));}
if (WAVETYPE==1 || WAVETYPE==3) {
L2t[itest]=2.0*(1.0+(L2sum/((float)NO_OF_TESTSHOTS*(float)ntr_glob)));
}
}
if (WAVETYPE==2) {
L2t[itest]=2.0*(1.0+(L2sum_SH/((float)NO_OF_TESTSHOTS*(float)ntr_glob)));
}
}
break;
case 8:
L2t[itest]=L2sum/energy_sum;
break;
default:
L2t[itest]=L2sum;
if ((WAVETYPE==1) || (WAVETYPE==3)) {
L2t[itest]=L2sum;
}
if (WAVETYPE==2) {
L2t[itest]=L2sum_SH;
}
break;
}
......@@ -3819,8 +4152,15 @@ int main(int argc, char **argv){
/* -----------------------------------------------------------------------*/
/* ----------- Do the actual update to the material parameters -----------*/
/* -----------------------------------------------------------------------*/
calc_mat_change_test(waveconv,waveconv_rho,waveconv_u,prho,prhonp1,ppi,ppinp1,pu,punp1,iter,1,FORWARD_ONLY,eps_scale,0,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
fprintf(FP,"=================================================================================================\n");
if (VTI) {
calc_mat_change_test_vti(waveconv,waveconv_rho,waveconv_u,waveconv_vshor_z,prho,prhonp1,ppi,ppinp1,pu,punp1,pvshor,pvshornp1,iter,1,FORWARD_ONLY,eps_scale,0,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,vshor_avg,LBFGS_iter_start);
}
else {
calc_mat_change_test(waveconv,waveconv_rho,waveconv_u,prho,prhonp1,ppi,ppinp1,pu,punp1,iter,1,FORWARD_ONLY,eps_scale,0,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
}
fprintf(FP,"=================================================================================================\n");
} /* end of if(FORWARD_ONLY!=4) */
......@@ -3830,6 +4170,11 @@ int main(int argc, char **argv){
if (FORWARD_ONLY==0 && (opteps_vp>0.0 || WOLFE_CONDITION)){
if(!ACOUSTIC){
if(WAVETYPE==1||WAVETYPE==3) if(MODEL_FILTER)smooth(ppi,4,2,Vs_avg,F_LOW_PASS);
if(WAVETYPE==2||WAVETYPE==3) {
if (VTI) {
if(MODEL_FILTER)smooth(pvshor,8,2,vshor_avg,F_LOW_PASS);
}
}
if(MODEL_FILTER)smooth(pu,5,2,Vs_avg,F_LOW_PASS);
if(MODEL_FILTER)smooth(prho,6,2,Vs_avg,F_LOW_PASS);
}else{
......@@ -4085,6 +4430,10 @@ int main(int argc, char **argv){
free_matrix(uttx,-nd+1,NY+nd,-nd+1,NX+nd);
free_matrix(utty,-nd+1,NY+nd,-nd+1,NX+nd);
}
if (WAVETYPE==2 || WAVETYPE==3) {
free_matrix(uxz,-nd+1,NY+nd,-nd+1,NX+nd);
free_matrix(uyz,-nd+1,NY+nd,-nd+1,NX+nd);
}
}
free_matrix(Vp0,-nd+1,NY+nd,-nd+1,NX+nd);
if(!ACOUSTIC)
......@@ -4128,6 +4477,25 @@ int main(int argc, char **argv){
free_vector(cip,1,L);
free_vector(cjm,1,L);
}
/* free memory for anisotropic modeling variables */
if (VTI) {
free_matrix ( pc11, -nd + 1, NY + nd, -nd + 1, NX + nd );
free_matrix ( pc13, -nd + 1, NY + nd, -nd + 1, NX + nd );
free_matrix ( pc33, -nd + 1, NY + nd, -nd + 1, NX + nd );
free_matrix ( pc55, -nd + 1, NY + nd, -nd + 1, NX + nd );
free_matrix ( pc55ipjp, -nd + 1, NY + nd, -nd + 1, NX + nd );
free_matrix ( pc66, -nd + 1, NY + nd, -nd + 1, NX + nd );
free_matrix ( pc66ipjp, -nd + 1, NY + nd, -nd + 1, NX + nd );
free_matrix ( pdelta, -nd+1,NY+nd,-nd+1,NX+nd);
free_matrix ( pepsilon, -nd+1,NY+nd,-nd+1,NX+nd);
if (WAVETYPE==2 || WAVETYPE==3) {
free_matrix ( pgamma, -nd+1,NY+nd,-nd+1,NX+nd);
free_matrix ( pvshornp1, -nd+1,NY+nd,-nd+1,NX+nd);
}
free_matrix (pvshor, -nd+1,NY+nd,-nd+1,NX+nd);
}
if(FORWARD_ONLY==0){
free_matrix(waveconv,-nd+1,NY+nd,-nd+1,NX+nd);
free_matrix(waveconv_lam,-nd+1,NY+nd,-nd+1,NX+nd);
......@@ -4301,11 +4669,17 @@ int main(int argc, char **argv){