Skip to content

GitLab

Commit 37dc3095 authored by Florian Wittkamp's avatar Florian Wittkamp
Browse files

SH: Viscoelastic modelling and free surface 

Updated SH viscoelastic modelling to the new grid.
parent 09f6740a
Show whitespace changes
Inline Side-by-side
src/av_mue.c
View file @ 37dc3095
......@@ -19,9 +19,9 @@
#include "fd.h"
void av_mue(float ** u, float ** uipjp, float ** uip, float ** ujp, float ** rho){
void av_mue(float ** u, float ** uipjp,float ** rho){
extern int NX, NY, INVMAT1, WAVETYPE;
extern int NX, NY, INVMAT1;
int i, j;
float u1, u2, u3, u4;
......@@ -31,17 +31,9 @@ void av_mue(float ** u, float ** uipjp, float ** uip, float ** ujp, float ** rho
for (i=1;i<=NX;i++){
uipjp[j][i]=4.0/((1.0/u[j][i])+(1.0/u[j][i+1])+(1.0/u[j+1][i])+(1.0/u[j+1][i+1]));
if(WAVETYPE==2 || WAVETYPE==3) {
uip[j][i]=0.5*(u[j][i+1]+u[j][i]);
ujp[j][i]=0.5*(u[j+1][i]+u[j][i]);
}
if((u[j][i]==0.0)||(u[j][i+1]==0.0)||(u[j+1][i]==0.0)||(u[j+1][i+1]==0.0)){
uipjp[j][i]=0.0;
if(WAVETYPE==2 || WAVETYPE==3) {
uip[j][i]=0.0;
ujp[j][i]=0.0;
}
}
......@@ -61,18 +53,9 @@ void av_mue(float ** u, float ** uipjp, float ** uip, float ** ujp, float ** rho
u4 = rho[j+1][i+1] * u[j+1][i+1] * u[j+1][i+1];
uipjp[j][i]=4.0/((1.0/u1)+(1.0/u2)+(1.0/u3)+(1.0/u4));
if(WAVETYPE==2 || WAVETYPE==3) {
uip[j][i]=0.5*(u1+u2);
ujp[j][i]=0.5*(u3+u1);
}
if((u1==0.0)||(u2==0.0)||(u3==0.0)||(u4==0.0)){
uipjp[j][i]=0.0;
if(WAVETYPE==2 || WAVETYPE==3) {
uip[j][i]=0.0;
ujp[j][i]=0.0;
}
}
}
}
......
src/av_tau.c
View file @ 37dc3095
......@@ -20,18 +20,14 @@
#include "fd.h"
void av_tau(float **taus, float **tausipjp,float **tausip,float **tausjp){
void av_tau(float **taus, float **tausipjp){
extern int NX, NY, WAVETYPE;
extern int NX, NY;
int i, j;
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
tausipjp[j][i] = 0.25*(taus[j][i]+taus[j][i+1]+taus[j+1][i]+taus[j+1][i+1]);
if(WAVETYPE==2 || WAVETYPE==3) {
tausip[j][i]=0.5*(taus[j][i+1]+taus[j][i]);
tausjp[j][i]=0.5*(taus[j+1][i]+taus[j][i]);
}
}
}
}
src/denise.c
View file @ 37dc3095
......@@ -68,7 +68,6 @@ int main(int argc, char **argv){
float ** gradg, ** gradp,** gradg_rho, ** gradp_rho, ** gradg_u, ** gradp_u, ** gradp_u_z,** gradp_rho_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;
float ** puip,** pujp;
float ***forward_prop_z_xz,***forward_prop_z_yz,***forward_prop_rho_z,**waveconv_mu_z;
float ** uxz, ** uyz;
......@@ -93,7 +92,7 @@ int main(int argc, char **argv){
/* 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;
float **tausjp, **tausip, ***pt=NULL, ***po=NULL; // SH Simulation
float ***pt=NULL, ***po=NULL; // SH Simulation
/* Variables for step length calculation */
int step1, step2, step3=0, itests, iteste, stepmax, countstep;
......@@ -472,10 +471,6 @@ int main(int argc, char **argv){
pu = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
punp1 = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
puipjp = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
if(WAVETYPE==2 || WAVETYPE==3) {
puip = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
pujp = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
}
}
vpmat = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
......@@ -507,8 +502,6 @@ int main(int argc, char **argv){
ptaus = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
ptausipjp = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
if(WAVETYPE==2 || WAVETYPE==3) {
tausip = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
tausjp = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
pt = f3tensor(-nd+1,NY+nd,-nd+1,NX+nd,1,L);
po = f3tensor(-nd+1,NY+nd,-nd+1,NX+nd,1,L);
}
......@@ -1033,9 +1026,9 @@ int main(int argc, char **argv){
/* end of MPI split for processors with ntr>0 */
if(!ACOUSTIC) av_mue(pu,puipjp,puip,pujp,prho);
if(!ACOUSTIC) av_mue(pu,puipjp,prho);
av_rho(prho,prip,prjp);
if (!ACOUSTIC && L) av_tau(ptaus,ptausipjp,tausip,tausjp);
if (!ACOUSTIC && L) av_tau(ptaus,ptausipjp);
/* Preparing memory variables for update_s (viscoelastic) */
......@@ -1295,7 +1288,7 @@ int main(int argc, char **argv){
if (!ACOUSTIC){
if (L){
/* viscoelastic */
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);
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);
......@@ -1753,7 +1746,7 @@ 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, puip, pujp, tausip, tausjp, bip, bjm, cip, cjm, etajm, etaip, 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);
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);
}
} else { /* elastic */
if (WAVETYPE==1 || WAVETYPE==3) {
......@@ -1779,7 +1772,7 @@ int main(int argc, char **argv){
if (!ACOUSTIC){
if (L){
/* viscoelastic */
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);
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);
......@@ -2243,7 +2236,7 @@ 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, puip, pujp, tausip, tausjp, bip, bjm, cip, cjm, etajm, etaip, 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);
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);
}
} else{
/* elastic */
......@@ -2268,7 +2261,7 @@ int main(int argc, char **argv){
if (!ACOUSTIC){
if (L){
/* viscoelastic */
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);
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);
......@@ -2882,7 +2875,7 @@ int main(int argc, char **argv){
waveconv_rho=joint_inversion_grad(waveconv_rho,waveconv_rho_z,JOINT_INVERSION_PSV_SH_ALPHA_RHO,JOINT_INVERSION_PSV_SH_TYPE);
/* Output joint gradient to disk */
sprintf(jac,"%s_joint_u_it%i",JACOBIAN,iter);
sprintf(jac,"%s_joint_vs_it%i",JACOBIAN,iter);
write_matrix_disk(waveconv_u, jac);
/* Output joint gradient to disk */
......@@ -3172,7 +3165,7 @@ int main(int argc, char **argv){
MPI_Barrier(MPI_COMM_WORLD);
if(!ACOUSTIC) av_mue(punp1,puipjp,puip,pujp,prhonp1);
if(!ACOUSTIC) av_mue(punp1,puipjp,prhonp1);
av_rho(prhonp1,prip,prjp);
/* Preparing memory variables for update_s (viscoelastic) */
......@@ -3343,7 +3336,7 @@ 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, puip, pujp, tausip, tausjp, bip, bjm, cip, cjm, etajm, etaip, 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);
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);
}
} else { /* elastic */
if(!ACOUSTIC){
......@@ -3373,7 +3366,7 @@ int main(int argc, char **argv){
if (!ACOUSTIC){
if (L){
/* viscoelastic */
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);
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);
......@@ -4251,15 +4244,9 @@ int main(int argc, char **argv){
free_matrix(gradp_u_z,-nd+1,NY+nd,-nd+1,NX+nd);
free_matrix(gradp_rho_z,-nd+1,NY+nd,-nd+1,NX+nd);
if(L){
free_matrix(tausip,-nd+1,NY+nd,-nd+1,NX+nd);
free_matrix(tausjp,-nd+1,NY+nd,-nd+1,NX+nd);
free_f3tensor(pt,-nd+1,NY+nd,-nd+1,NX+nd,1,L);
free_f3tensor(po,-nd+1,NY+nd,-nd+1,NX+nd,1,L);
}
if(!ACOUSTIC){
free_matrix(puip,-nd+1,NY+nd,-nd+1,NX+nd);
free_matrix(pujp,-nd+1,NY+nd,-nd+1,NX+nd);
}
}
......
src/fd.h
View file @ 37dc3095
......@@ -53,11 +53,11 @@ float norm(float ** waveconv, int iter, int sws);
void av_mat(float ** pi, float ** u,
float ** ppijm, float ** puip, float ** pujm);
void av_mue(float ** u, float ** uipjp, float ** uip, float ** ujp, float ** rho);
void av_mue(float ** u, float ** uipjp,float ** rho);
void av_rho(float **rho, float **rip, float **rjp);
void av_tau(float **taus, float **tausipjp,float **tausip,float **tausjp);
void av_tau(float **taus, float **tausipjp);
float median2d(float **mat, int ny, int nx);
......@@ -283,10 +283,9 @@ void surface(int ndepth, float ** pvx, float ** pvy,
void surface_elastic(int ndepth, float ** vx, float ** vy, float ** sxx, float ** syy,
float ** sxy, float ** pi, float ** u, float ** rho, float * hc);
void surface_elastic_PML(int ndepth, float ** vx, float ** vy, float ** sxx, float ** syy,
float ** sxy, float ** syz, float ** pi, float ** u, float ** rho, float * hc, float * K_x, float * a_x, float * b_x, float ** psi_vxx, float ** ux, float ** uy, float ** uxy, float ** uyz,float ** sxz,float **uxz);
void surface_elastic_PML(int ndepth, float ** vx, float ** vy, float ** sxx, float ** syy, float ** sxy, float ** syz, float ** pi, float ** u, float ** rho, float * hc, float * K_x, float * a_x, float * b_x, float ** psi_vxx, float ** ux, float ** uy, float ** uxy, float ** uyz,float ** sxz,float **uxz);
void surface_PML(int ndepth, float ** vx, float ** vy, float ** sxx, float ** syy, float ** sxy, float ** syz, float ***p, float ***q, float ** ppi, float ** pu, float **prho, float **ptaup, float **ptaus, float *etajm, float *peta, float * hc, float * K_x, float * a_x, float * b_x, float ** psi_vxx,float ** ux, float ** uy, float ** uxy);
void surface_PML(int ndepth, float ** vx, float ** vy, float ** sxx, float ** syy, float ** sxy, float ** syz, float ***p, float ***q, float ** ppi, float ** pu, float **prho, float **ptaup, float **ptaus, float *etajm, float *peta, float * hc, float * K_x, float * a_x, float * b_x, float ** psi_vxx, float ** ux, float ** uy, float ** uxy, float ** uyz,float ** sxz,float **uxz);
void timedomain_filt(float ** data, float fc, int order, int ntr, int ns, int method);
void timedomain_filt_vector(float * data, float fc, int order, int ntr, int ns, int method);
......@@ -313,8 +312,7 @@ void update_s_visc_hc(int nx1, int nx2, int ny1, int ny2,
void update_s_elastic_PML_SH(int nx1, int nx2, int ny1, int ny2, float ** vz, float ** sxz, float ** syz, float ** uxz, float ** uyz, float *hc, int infoout,float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half,
float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half,float ** psi_vzx, float ** psi_vzy,float ** uipjp,float ** u,float ** rho);
void update_s_visc_PML_SH(int nx1, int nx2, int ny1, int ny2, float ** vz, float ** sxz, float ** syz, float ***t, float ***o, float ** uip, float ** ujp, float ** tausip, float ** tausjp, float *bip, float *bjm, float *cip, float *cjm, float *etajm, float *etaip, float *hc, int infoout,float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half,
float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half,float ** psi_vzx, float ** psi_vzy);
void update_s_visc_PML_SH(int nx1, int nx2, int ny1, int ny2, float ** vz, float ** sxz, float ** syz, float ***t, float ***o, float *bip, float *bjm, float *cip, float *cjm, float ***d, float ***dip, float **fipjp, float **f, float *hc, int infoout,float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half,float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half,float ** psi_vzx, float ** psi_vzy);
void update_s_rsg(int nx1, int nx2, int ny1, int ny2,
float ** pvx, float ** pvy, float ** psxx, float ** psyy,
......
src/surface_PML.c
View file @ 37dc3095
......@@ -27,7 +27,7 @@
void surface_PML(int ndepth, float ** vx, float ** vy, float ** sxx, float ** syy,
float ** sxy, float ** syz, float ***p, float ***q, float ** ppi, float ** pu, float **prho, float **ptaup, float **ptaus, float *etajm, float *peta, float * hc, float * K_x, float * a_x, float * b_x, float ** psi_vxx,
float ** ux, float ** uy, float ** uxy){
float ** ux, float ** uy, float ** uxy, float ** uyz,float ** sxz,float **uxz){
int i,j,m,h,h1,l;
......@@ -162,10 +162,15 @@ void surface_PML(int ndepth, float ** vx, float ** vy, float ** sxx, float ** sy
}
if (WAVETYPE==2||WAVETYPE==3){
for (i=1;i<=NX;i++){
for (m=0; m<=fdoh-1; m++) {
/* mirroring of stress components for syz
to make a stress free surface in depth y=jsurf*h */
syz[j-m][i]=-syz[j+m+1][i];
syz[j][i]=0.0;
uyz[j][i]=0.0;
for (m=1; m<=fdoh; m++) {
syz[j-m][i]=-syz[j+m][i];
sxz[j-m][i]=-sxz[j+m-1][i];
uyz[j-m][i]=-uyz[j+m][i];
uxz[j-m][i]=-uxz[j+m-1][i];
}
}
}
......
src/update_s_visc_PML.c
View file @ 37dc3095
......@@ -16,7 +16,6 @@
* along with DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/* $Id: update_s_visc_PML.c,v 1.1.1.1 2011/10/06 22:44:52 groos Exp $*/
/*------------------------------------------------------------------------
* updating stress components at gridpoints [nx1...nx2][ny1...ny2]
* by a staggered grid finite difference scheme of arbitrary (FDORDER) order accuracy in space
......
src/update_s_visc_PML_SH.c
View file @ 37dc3095
......@@ -22,14 +22,13 @@
* and second order accuracy in time
* T. Bohlen
*
* SH-Version
* SH-Version F. Wittkamp
*
* ----------------------------------------------------------------------*/
#include "fd.h"
void update_s_visc_PML_SH(int nx1, int nx2, int ny1, int ny2, float ** vz, float ** sxz, float ** syz, float ***t, float ***o, float ** uip, float ** ujp, float ** tausip, float ** tausjp, float *bip, float *bjm, float *cip, float *cjm, float *etajm, float *etaip, float *hc, int infoout,float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half,
float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half,float ** psi_vzx, float ** psi_vzy){
void update_s_visc_PML_SH(int nx1, int nx2, int ny1, int ny2, float ** vz, float ** sxz, float ** syz, float ***t, float ***o, float *bip, float *bjm, float *cip, float *cjm, float ***d, float ***dip, float **fipjp, float **f, float *hc, int infoout,float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half,float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half,float ** psi_vzx, float ** psi_vzy){
int i,j, m, fdoh, h, h1, l;
......@@ -45,9 +44,6 @@ void update_s_visc_PML_SH(int nx1, int nx2, int ny1, int ny2, float ** vz, floa
float sumo=0.0, sumt=0.0;
/*dhi = DT/DH;*/
dhi=1.0/DH;
fdoh = FDORDER/2;
dthalbe = DT/2.0;
......@@ -60,20 +56,83 @@ void update_s_visc_PML_SH(int nx1, int nx2, int ny1, int ny2, float ** vz, floa
}
fprintf(FP," CHECK MATERIAL PARAMETER CHECK MATERIAL PARAMETER \n");
switch (FDORDER){
case 2:
for (j=ny1;j<=ny2;j++){
for (i=nx1;i<=nx2;i++){
vzx = ( hc[1]*(vz[j][i+1]-vz[j][i]))*dhi;
vzy = ( hc[1]*(vz[j][i]-vz[j-1][i]))*dhi;
/* left boundary */
if((!BOUNDARY) && (POS[1]==0) && (i<=FW)){
psi_vzx[j][i] = b_x_half[i] * psi_vzx[j][i] + a_x_half[i] * vzx;
vzx = vzx / K_x_half[i] + psi_vzx[j][i];
}
/* right boundary */
if((!BOUNDARY) && (POS[1]==NPROCX-1) && (i>=nx2-FW+1)){
h1 = (i-nx2+2*FW);
h = i;
psi_vzx[j][h1] = b_x_half[h1] * psi_vzx[j][h1] + a_x_half[h1] * vzx;
vzx = vzx / K_x_half[h1] + psi_vzx[j][h1];
}
/* top boundary */
if((POS[2]==0) && (!(FREE_SURF)) && (j<=FW)){
psi_vzy[j][i] = b_y[j] * psi_vzy[j][i] + a_y[j] * vzy;
vzy = vzy / K_y[j] + psi_vzy[j][i];
}
/* bottom boundary */
if((POS[2]==NPROCY-1) && (j>=ny2-FW+1)){
h1 = (j-ny2+2*FW);
h = j;
psi_vzy[h1][i] = b_y[h1] * psi_vzy[h1][i] + a_y[h1] * vzy;
vzy = vzy / K_y[h1] + psi_vzy[h1][i];
}
/* computing sums of the old memory variables */
sumt=sumo=0.0;
for (l=1;l<=L;l++){
sumo+=o[j][i][l];
sumt+=t[j][i][l];
}
/* updating components of the stress tensor, partially */
sxz[j][i]+=(fipjp[j][i]*vzx)+(dthalbe*sumo);
syz[j][i]+=(f[j][i]*vzy)+(dthalbe*sumt);
/* now updating the memory-variables and sum them up*/
sumt=sumo=0.0;
for (l=1;l<=L;l++){
o[j][i][l]=bip[l]*(o[j][i][l]*cip[l]-(dip[j][i][l]*vzx));
t[j][i][l]=bjm[l]*(t[j][i][l]*cjm[l]-(d[j][i][l]*vzy));
sumt+=t[j][i][l];
sumo+=o[j][i][l];
}
/* and now the components of the stress tensor are
completely updated */
sxz[j][i]+=(dthalbe*sumo);
syz[j][i]+=(dthalbe*sumt);
}
}
break;
case 4:
for (j=ny1;j<=ny2;j++){
for (i=nx1;i<=nx2;i++){
vzx = ( hc[1]*(vz[j][i+1]-vz[j][i])
+ hc[2]*(vz[j][i+2]-vz[j][i-1]))*dhi;
vzy = ( hc[1]*(vz[j+1][i]-vz[j][i])
+ hc[2]*(vz[j+2][i]-vz[j-1][i]))*dhi;
vzy = ( hc[1]*(vz[j][i]-vz[j-1][i])
+ hc[2]*(vz[j+1][i]-vz[j-2][i]))*dhi;
/* left boundary */
......@@ -104,8 +163,6 @@ void update_s_visc_PML_SH(int nx1, int nx2, int ny1, int ny2, float ** vz, floa
vzy = vzy / K_y[h1] + psi_vzy[h1][i];
}
/* computing sums of the old memory variables */
sumt=sumo=0.0;
for (l=1;l<=L;l++){
......@@ -114,15 +171,15 @@ void update_s_visc_PML_SH(int nx1, int nx2, int ny1, int ny2, float ** vz, floa
}
/* updating components of the stress tensor, partially */
sxz[j][i]+=(uip[j][i]*DT*(1.0+L*tausip[j][i])*vzx)+(dthalbe*sumo);
syz[j][i]+=(ujp[j][i]*DT*(1.0+L*tausjp[j][i])*vzy)+(dthalbe*sumt);
sxz[j][i]+=(fipjp[j][i]*vzx)+(dthalbe*sumo);
syz[j][i]+=(f[j][i]*vzy)+(dthalbe*sumt);
/* now updating the memory-variables and sum them up*/
sumt=sumo=0.0;
for (l=1;l<=L;l++){
o[j][i][l]=bip[l]*(o[j][i][l]*cip[l]-(uip[j][i]*etaip[l]*tausip[j][i]*vzx));
t[j][i][l]=bjm[l]*(t[j][i][l]*cjm[l]-(ujp[j][i]*etajm[l]*tausjp[j][i]*vzy));
o[j][i][l]=bip[l]*(o[j][i][l]*cip[l]-(dip[j][i][l]*vzx));
t[j][i][l]=bjm[l]*(t[j][i][l]*cjm[l]-(d[j][i][l]*vzy));
sumt+=t[j][i][l];
sumo+=o[j][i][l];
}
......@@ -131,18 +188,298 @@ void update_s_visc_PML_SH(int nx1, int nx2, int ny1, int ny2, float ** vz, floa
completely updated */
sxz[j][i]+=(dthalbe*sumo);
syz[j][i]+=(dthalbe*sumt);
//fprintf(FP," Signal: %f\n",sxz[j][i]);
}
}
break;
case 6:
for (j=ny1;j<=ny2;j++){
for (i=nx1;i<=nx2;i++){
vzx = ( hc[1]*(vz[j][i+1]-vz[j][i])
+ hc[2]*(vz[j][i+2]-vz[j][i-1])
+ hc[3]*(vz[j][i+3]-vz[j][i-2]))*dhi;
vzy = ( hc[1]*(vz[j][i]-vz[j-1][i])
+ hc[2]*(vz[j+1][i]-vz[j-2][i])
+ hc[3]*(vz[j+2][i]-vz[j-3][i]))*dhi;
/* left boundary */
if((!BOUNDARY) && (POS[1]==0) && (i<=FW)){
psi_vzx[j][i] = b_x_half[i] * psi_vzx[j][i] + a_x_half[i] * vzx;
vzx = vzx / K_x_half[i] + psi_vzx[j][i];
}
/* right boundary */
if((!BOUNDARY) && (POS[1]==NPROCX-1) && (i>=nx2-FW+1)){
h1 = (i-nx2+2*FW);
h = i;
psi_vzx[j][h1] = b_x_half[h1] * psi_vzx[j][h1] + a_x_half[h1] * vzx;
vzx = vzx / K_x_half[h1] + psi_vzx[j][h1];
}
/* top boundary */
if((POS[2]==0) && (!(FREE_SURF)) && (j<=FW)){
psi_vzy[j][i] = b_y[j] * psi_vzy[j][i] + a_y[j] * vzy;
vzy = vzy / K_y[j] + psi_vzy[j][i];
}
/* bottom boundary */
if((POS[2]==NPROCY-1) && (j>=ny2-FW+1)){
h1 = (j-ny2+2*FW);
h = j;
psi_vzy[h1][i] = b_y[h1] * psi_vzy[h1][i] + a_y[h1] * vzy;
vzy = vzy / K_y[h1] + psi_vzy[h1][i];
}
/* computing sums of the old memory variables */
sumt=sumo=0.0;
for (l=1;l<=L;l++){
sumo+=o[j][i][l];
sumt+=t[j][i][l];
}
/* updating components of the stress tensor, partially */
sxz[j][i]+=(fipjp[j][i]*vzx)+(dthalbe*sumo);
syz[j][i]+=(f[j][i]*vzy)+(dthalbe*sumt);
/* now updating the memory-variables and sum them up*/
sumt=sumo=0.0;
for (l=1;l<=L;l++){
o[j][i][l]=bip[l]*(o[j][i][l]*cip[l]-(dip[j][i][l]*vzx));
t[j][i][l]=bjm[l]*(t[j][i][l]*cjm[l]-(d[j][i][l]*vzy));
sumt+=t[j][i][l];
sumo+=o[j][i][l];
}
/* and now the components of the stress tensor are
completely updated */
sxz[j][i]+=(dthalbe*sumo);
syz[j][i]+=(dthalbe*sumt);
}
}
break;
case 8:
for (j=ny1;j<=ny2;j++){
for (i=nx1;i<=nx2;i++){
vzx = ( hc[1]*(vz[j][i+1]-vz[j][i])
+ hc[2]*(vz[j][i+2]-vz[j][i-1])
+ hc[3]*(vz[j][i+3]-vz[j][i-2])
+ hc[4]*(vz[j][i+4]-vz[j][i-3]))*dhi;
vzy = ( hc[1]*(vz[j][i]-vz[j-1][i])
+ hc[2]*(vz[j+1][i]-vz[j-2][i])
+ hc[3]*(vz[j+2][i]-vz[j-3][i])
+ hc[4]*(vz[j+3][i]-vz[j-4][i]))*dhi;
/* left boundary */
if((!BOUNDARY) && (POS[1]==0) && (i<=FW)){
psi_vzx[j][i] = b_x_half[i] * psi_vzx[j][i] + a_x_half[i] * vzx;
vzx = vzx / K_x_half[i] + psi_vzx[j][i];
}
/* right boundary */
if((!BOUNDARY) && (POS[1]==NPROCX-1) && (i>=nx2-FW+1)){
h1 = (i-nx2+2*FW);
h = i;
psi_vzx[j][h1] = b_x_half[h1] * psi_vzx[j][h1] + a_x_half[h1] * vzx;
vzx = vzx / K_x_half[h1] + psi_vzx[j][h1];
}
/* top boundary */
if((POS[2]==0) && (!(FREE_SURF)) && (j<=FW)){
psi_vzy[j][i] = b_y[j] * psi_vzy[j][i] + a_y[j] * vzy;
vzy = vzy / K_y[j] + psi_vzy[j][i];
}
/* bottom boundary */
if((POS[2]==NPROCY-1) && (j>=ny2-FW+1)){
h1 = (j-ny2+2*FW);
h = j;
psi_vzy[h1][i] = b_y[h1] * psi_vzy[h1][i] + a_y[h1] * vzy;
vzy = vzy / K_y[h1] + psi_vzy[h1][i];
}