Commit 99111605 authored by laura.gassner's avatar laura.gassner

Add missing higher orders for viscoacoustic modeling with PML

parent b77405bb
...@@ -34,7 +34,7 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** ...@@ -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){ float ** psi_vxx, float ** psi_vyy, float ** psi_vxy, float ** psi_vyx){
int i,j, m, fdoh, h, h1, l; int i,j, m, fdoh, h, h1, l;
float vxx, vyy, vxy, vyx; float vxx, vyy;
float dhi, dthalbe; float dhi, dthalbe;
extern float DT, DH; extern float DT, DH;
extern int MYID, FDORDER, FW, L; extern int MYID, FDORDER, FW, L;
...@@ -62,24 +62,16 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** ...@@ -62,24 +62,16 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float **
for (j=ny1;j<=ny2;j++){ for (j=ny1;j<=ny2;j++){
for (i=nx1;i<=nx2;i++){ for (i=nx1;i<=nx2;i++){
vxx = ( hc[1]*(vx[j][i] -vx[j][i-1]))*dhi; 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 */
vxy = ( hc[1]*(vx[j+1][i]-vx[j][i]))*dhi;
vyy = ( hc[1]*(vy[j][i] -vy[j-1][i]))*dhi;
/* left boundary */
if((!BOUNDARY) && (POS[1]==0) && (i<=FW)){ if((!BOUNDARY) && (POS[1]==0) && (i<=FW)){
psi_vxx[j][i] = b_x[i] * psi_vxx[j][i] + a_x[i] * vxx; psi_vxx[j][i] = b_x[i] * psi_vxx[j][i] + a_x[i] * vxx;
vxx = vxx / K_x[i] + psi_vxx[j][i]; 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 */ /* right boundary */
if((!BOUNDARY) && (POS[1]==NPROCX-1) && (i>=nx2-FW+1)){ if((!BOUNDARY) && (POS[1]==NPROCX-1) && (i>=nx2-FW+1)){
h1 = (i-nx2+2*FW); h1 = (i-nx2+2*FW);
...@@ -87,38 +79,23 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** ...@@ -87,38 +79,23 @@ 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; psi_vxx[j][h1] = b_x[h1] * psi_vxx[j][h1] + a_x[h1] * vxx;
vxx = vxx / K_x[h1] + psi_vxx[j][h1]; vxx = vxx / K_x[h1] + psi_vxx[j][h1];
/*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];*/
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 */ /* top boundary */
if((POS[2]==0) && (!(FREE_SURF)) && (j<=FW)){ if((POS[2]==0) && (!(FREE_SURF)) && (j<=FW)){
psi_vyy[j][i] = b_y[j] * psi_vyy[j][i] + a_y[j] * vyy; 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];
vyy = vyy / K_y[j] + psi_vyy[j][i];
vxy = vxy / K_y_half[j] + psi_vxy[j][i];
} }
/* bottom boundary */ /* bottom boundary */
if((POS[2]==NPROCY-1) && (j>=ny2-FW+1)){ if((POS[2]==NPROCY-1) && (j>=ny2-FW+1)){
h1 = (j-ny2+2*FW); h1 = (j-ny2+2*FW);
h = j; h = j;
psi_vyy[h1][i] = b_y[h1] * psi_vyy[h1][i] + a_y[h1] * vyy; psi_vyy[h1][i] = b_y[h1] * psi_vyy[h1][i] + a_y[h1] * vyy;
vyy = vyy / K_y[h1] + psi_vyy[h1][i]; vyy = vyy / K_y[h1] + psi_vyy[h1][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_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];
} }
/* computing sums of the old memory variables */ /* computing sums of the old memory variables */
...@@ -127,20 +104,250 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** ...@@ -127,20 +104,250 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float **
sump+=p[j][i][l]; sump+=p[j][i][l];
} }
/* updating components of the stress tensor, partially */ /* updating components of the stress tensor, partially */
sp[j][i] += (g[j][i]*(vxx+vyy))+(dthalbe*sump); sp[j][i] += (g[j][i]*(vxx+vyy))-(2.0*f[j][i]*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*/ /* now updating the memory-variables and sum them up*/
sump=0.0; sump=0.0;
for (l=1;l<=L;l++){ for (l=1;l<=L;l++){
p[j][i][l] = bjm[l]*(p[j][i][l]*cjm[l]-(e[j][i][l]*(vxx+vyy))); p[j][i][l] = bjm[l]*(p[j][i][l]*cjm[l]-(e[j][i][l]*(vxx+vyy))+(2.0*d[j][i][l]*vyy));
sump += p[j][i][l]; sump += p[j][i][l];
} }
/* and now the components of the stress tensor are /* and now the components of the stress tensor are
completely updated */ completely updated */
sp[j][i]+=(dthalbe*sump); 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];
}
/* 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))-(2.0*f[j][i]*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))+(2.0*d[j][i][l]*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)){
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))-(2.0*f[j][i]*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))+(2.0*d[j][i][l]*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_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))-(2.0*f[j][i]*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))+(2.0*d[j][i][l]*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; break;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment