Skip to content

GitLab

IFOS2D
IFOS2D
Commit f3c72393 authored by laura.gassner's avatar laura.gassner
Browse files
Options

Add missing higher orders for viscoacoustic modeling with PML

parent b77405bb
Hide whitespace changes
Inline Side-by-side
Showing with 242 additions and 41 deletions
src/update_p_visc_PML.c
View file @ f3c72393
......@@ -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,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 (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;
vxy = ( hc[1]*(vx[j+1][i]-vx[j][i]))*dhi;
vyy = ( hc[1]*(vy[j][i] -vy[j-1][i]))*dhi;
/* left boundary */
/* 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 */
/* right boundary */
if((!BOUNDARY) && (POS[1]==NPROCX-1) && (i>=nx2-FW+1)){
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 **
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][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)){
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];
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 */
/* bottom boundary */
if((POS[2]==NPROCY-1) && (j>=ny2-FW+1)){
h1 = (j-ny2+2*FW);
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];
/*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];
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 */
......@@ -127,9 +104,13 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float **
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++){
......@@ -139,8 +120,234 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float **
/* and now the components of the stress tensor are
completely updated */
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];
}
/* 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 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))+(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_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;
......@@ -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++){
......
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