...
 
......@@ -34,7 +34,7 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float **
float ** psi_vxx, float ** psi_vyy, float ** psi_vxy, float ** psi_vyx){
int i,j, m, fdoh, h, h1, l;
float vxx, vyy, vxy, vyx;
float vxx, vyy;
float dhi, dthalbe;
extern float DT, DH;
extern int MYID, FDORDER, FW, L;
......@@ -62,21 +62,85 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float **
for (j=ny1;j<=ny2;j++){
for (i=nx1;i<=nx2;i++){
vxx = ( hc[1]*(vx[j][i] -vx[j][i-1]))*dhi;
vyy = ( hc[1]*(vy[j][i] -vy[j-1][i]))*dhi;
vyx = ( hc[1]*(vy[j][i+1]-vy[j][i]))*dhi;
/* left boundary */
if((!BOUNDARY) && (POS[1]==0) && (i<=FW)){
vxy = ( hc[1]*(vx[j+1][i]-vx[j][i]))*dhi;
psi_vxx[j][i] = b_x[i] * psi_vxx[j][i] + a_x[i] * vxx;
vxx = vxx / K_x[i] + psi_vxx[j][i];
}
vyy = ( hc[1]*(vy[j][i] -vy[j-1][i]))*dhi;
/* right boundary */
if((!BOUNDARY) && (POS[1]==NPROCX-1) && (i>=nx2-FW+1)){
h1 = (i-nx2+2*FW);
h = i;
psi_vxx[j][h1] = b_x[h1] * psi_vxx[j][h1] + a_x[h1] * vxx;
vxx = vxx / K_x[h1] + psi_vxx[j][h1];
}
/* top boundary */
if((POS[2]==0) && (!(FREE_SURF)) && (j<=FW)){
psi_vyy[j][i] = b_y[j] * psi_vyy[j][i] + a_y[j] * vyy;
vyy = vyy / K_y[j] + psi_vyy[j][i];
}
/* bottom boundary */
if((POS[2]==NPROCY-1) && (j>=ny2-FW+1)){
h1 = (j-ny2+2*FW);
h = j;
psi_vyy[h1][i] = b_y[h1] * psi_vyy[h1][i] + a_y[h1] * vyy;
vyy = vyy / K_y[h1] + psi_vyy[h1][i];
}
/* computing sums of the old memory variables */
sump=0.0;
for (l=1;l<=L;l++){
sump+=p[j][i][l];
}
/* updating components of the stress tensor, partially */
sp[j][i] += (g[j][i]*(vxx+vyy))+(dthalbe*sump);
// u[j][i] = ((g[j][i]/DT)*(vxx+vyy))-((2.0*f[j][i]/DT)*vyy)+(0.5*sump);
/* now updating the memory-variables and sum them up*/
sump=0.0;
for (l=1;l<=L;l++){
p[j][i][l] = bjm[l]*(p[j][i][l]*cjm[l]-(e[j][i][l]*(vxx+vyy)));
sump += p[j][i][l];
}
/* and now the components of the stress tensor are
completely updated */
sp[j][i]+=(dthalbe*sump);
// u[j][i]+=(0.5*sump);
}
}
break;
case 4:
for (j=ny1;j<=ny2;j++){
for (i=nx1;i<=nx2;i++){
vxx = ( hc[1]*(vx[j][i] -vx[j][i-1])
+ hc[2]*(vx[j][i+1]-vx[j][i-2]))*dhi;
vyy = ( hc[1]*(vy[j][i] -vy[j-1][i])
+ hc[2]*(vy[j+1][i]-vy[j-2][i]))*dhi;
/* left boundary */
if((!BOUNDARY) && (POS[1]==0) && (i<=FW)){
psi_vxx[j][i] = b_x[i] * psi_vxx[j][i] + a_x[i] * vxx;
vxx = vxx / K_x[i] + psi_vxx[j][i];
psi_vyx[j][i] = b_x_half[i] * psi_vyx[j][i] + a_x_half[i] * vyx;
vyx = vyx / K_x_half[i] + psi_vyx[j][i];
}
/* right boundary */
......@@ -87,22 +151,90 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float **
psi_vxx[j][h1] = b_x[h1] * psi_vxx[j][h1] + a_x[h1] * vxx;
vxx = vxx / K_x[h1] + psi_vxx[j][h1];
}
/* top boundary */
if((POS[2]==0) && (!(FREE_SURF)) && (j<=FW)){
psi_vyy[j][i] = b_y[j] * psi_vyy[j][i] + a_y[j] * vyy;
vyy = vyy / K_y[j] + psi_vyy[j][i];
}
/* bottom boundary */
if((POS[2]==NPROCY-1) && (j>=ny2-FW+1)){
h1 = (j-ny2+2*FW);
h = j;
psi_vyy[h1][i] = b_y[h1] * psi_vyy[h1][i] + a_y[h1] * vyy;
vyy = vyy / K_y[h1] + psi_vyy[h1][i];
}
/* computing sums of the old memory variables */
sump=0.0;
for (l=1;l<=L;l++){
sump+=p[j][i][l];
}
/* updating components of the stress tensor, partially */
sp[j][i] += (g[j][i]*(vxx+vyy))+(dthalbe*sump);
// u[j][i] = ((g[j][i]/DT)*(vxx+vyy))-((2.0*f[j][i]/DT)*vyy)+(0.5*sump);
/* now updating the memory-variables and sum them up*/
sump=0.0;
for (l=1;l<=L;l++){
p[j][i][l] = bjm[l]*(p[j][i][l]*cjm[l]-(e[j][i][l]*(vxx+vyy)));
sump += p[j][i][l];
}
/* and now the components of the stress tensor are
completely updated */
sp[j][i]+=(dthalbe*sump);
// u[j][i]+=(0.5*sump);
}
}
break;
case 6:
for (j=ny1;j<=ny2;j++){
for (i=nx1;i<=nx2;i++){
vxx = ( hc[1]*(vx[j][i] -vx[j][i-1])
+ hc[2]*(vx[j][i+1]-vx[j][i-2])
+ hc[3]*(vx[j][i+2]-vx[j][i-3]))*dhi;
vyy = ( hc[1]*(vy[j][i] -vy[j-1][i])
+ hc[2]*(vy[j+1][i]-vy[j-2][i])
+ hc[3]*(vy[j+2][i]-vy[j-3][i]))*dhi;
/* left boundary */
if((!BOUNDARY) && (POS[1]==0) && (i<=FW)){
psi_vxx[j][i] = b_x[i] * psi_vxx[j][i] + a_x[i] * vxx;
vxx = vxx / K_x[i] + psi_vxx[j][i];
}
/* right boundary */
if((!BOUNDARY) && (POS[1]==NPROCX-1) && (i>=nx2-FW+1)){
/*psi_vyx[j][h] = b_x_half[h] * psi_vyx[j][h] + a_x_half[h] * vyx;
vyx = vyx / K_x_half[h] + psi_vyx[j][h];*/
h1 = (i-nx2+2*FW);
h = i;
psi_vxx[j][h1] = b_x[h1] * psi_vxx[j][h1] + a_x[h1] * vxx;
vxx = vxx / K_x[h1] + psi_vxx[j][h1];
psi_vyx[j][h1] = b_x_half[h1] * psi_vyx[j][h1] + a_x_half[h1] * vyx;
vyx = vyx / K_x_half[h1] + psi_vyx[j][h1];
}
/* top boundary */
if((POS[2]==0) && (!(FREE_SURF)) && (j<=FW)){
psi_vyy[j][i] = b_y[j] * psi_vyy[j][i] + a_y[j] * vyy;
psi_vxy[j][i] = b_y_half[j] * psi_vxy[j][i] + a_y_half[j] * vxy;
vyy = vyy / K_y[j] + psi_vyy[j][i];
vxy = vxy / K_y_half[j] + psi_vxy[j][i];
}
/* bottom boundary */
......@@ -113,12 +245,82 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float **
psi_vyy[h1][i] = b_y[h1] * psi_vyy[h1][i] + a_y[h1] * vyy;
vyy = vyy / K_y[h1] + psi_vyy[h1][i];
}
/* computing sums of the old memory variables */
sump=0.0;
for (l=1;l<=L;l++){
sump+=p[j][i][l];
}
/* updating components of the stress tensor, partially */
sp[j][i] += (g[j][i]*(vxx+vyy))+(dthalbe*sump);
// u[j][i] = ((g[j][i]/DT)*(vxx+vyy))-((2.0*f[j][i]/DT)*vyy)+(0.5*sump);
/* now updating the memory-variables and sum them up*/
sump=0.0;
for (l=1;l<=L;l++){
p[j][i][l] = bjm[l]*(p[j][i][l]*cjm[l]-(e[j][i][l]*(vxx+vyy)));
sump += p[j][i][l];
}
/* and now the components of the stress tensor are
completely updated */
sp[j][i]+=(dthalbe*sump);
// u[j][i]+=(0.5*sump);
}
}
break;
case 8:
for (j=ny1;j<=ny2;j++){
for (i=nx1;i<=nx2;i++){
vxx = ( hc[1]*(vx[j][i] -vx[j][i-1])
+ hc[2]*(vx[j][i+1]-vx[j][i-2])
+ hc[3]*(vx[j][i+2]-vx[j][i-3])
+ hc[4]*(vx[j][i+3]-vx[j][i-4]))*dhi;
vyy = ( hc[1]*(vy[j][i] -vy[j-1][i])
+ hc[2]*(vy[j+1][i]-vy[j-2][i])
+ hc[3]*(vy[j+2][i]-vy[j-3][i])
+ hc[4]*(vy[j+3][i]-vy[j-4][i]))*dhi;
/* left boundary */
if((!BOUNDARY) && (POS[1]==0) && (i<=FW)){
psi_vxx[j][i] = b_x[i] * psi_vxx[j][i] + a_x[i] * vxx;
vxx = vxx / K_x[i] + psi_vxx[j][i];
}
/* right boundary */
if((!BOUNDARY) && (POS[1]==NPROCX-1) && (i>=nx2-FW+1)){
h1 = (i-nx2+2*FW);
h = i;
/*psi_vxy[j][i] = b_y_half[j] * psi_vxy[j][i] + a_y_half[j] * vxy;
vxy = vxy / K_y_half[j] + psi_vxy[j][i];*/
psi_vxx[j][h1] = b_x[h1] * psi_vxx[j][h1] + a_x[h1] * vxx;
vxx = vxx / K_x[h1] + psi_vxx[j][h1];
}
psi_vxy[h1][i] = b_y_half[h1] * psi_vxy[h1][i] + a_y_half[h1] * vxy;
vxy = vxy / K_y_half[h1] + psi_vxy[h1][i];
/* top boundary */
if((POS[2]==0) && (!(FREE_SURF)) && (j<=FW)){
psi_vyy[j][i] = b_y[j] * psi_vyy[j][i] + a_y[j] * vyy;
vyy = vyy / K_y[j] + psi_vyy[j][i];
}
/* bottom boundary */
if((POS[2]==NPROCY-1) && (j>=ny2-FW+1)){
h1 = (j-ny2+2*FW);
h = j;
psi_vyy[h1][i] = b_y[h1] * psi_vyy[h1][i] + a_y[h1] * vyy;
vyy = vyy / K_y[h1] + psi_vyy[h1][i];
}
/* computing sums of the old memory variables */
......@@ -130,6 +332,9 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float **
/* updating components of the stress tensor, partially */
sp[j][i] += (g[j][i]*(vxx+vyy))+(dthalbe*sump);
// u[j][i] = ((g[j][i]/DT)*(vxx+vyy))-((2.0*f[j][i]/DT)*vyy)+(0.5*sump);
/* now updating the memory-variables and sum them up*/
sump=0.0;
for (l=1;l<=L;l++){
......@@ -137,10 +342,12 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float **
sump += p[j][i][l];
}
/* and now the components of the stress tensor are
completely updated */
sp[j][i]+=(dthalbe*sump);
// u[j][i]+=(0.5*sump);
}
}
break;
......@@ -150,18 +357,12 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float **
for (i=nx1;i<=nx2;i++){
vxx = 0.0;
vyy = 0.0;
vyx = 0.0;
vxy = 0.0;
for (m=1; m<=fdoh; m++) {
vxx += hc[m]*(vx[j][i+m-1] -vx[j][i-m] );
vyy += hc[m]*(vy[j+m-1][i] -vy[j-m][i] );
vyx += hc[m]*(vy[j][i+m] -vy[j][i-m+1]);
vxy += hc[m]*(vx[j+m][i] -vx[j-m+1][i]);
}
vxx *= dhi;
vyy *= dhi;
vyx *= dhi;
vxy *= dhi;
sump=0.0;
for (l=1;l<=L;l++){
......