Commit b115902b authored by Tilman Steinweg's avatar Tilman Steinweg

implemented cpml in surface functions

CPML is now applied to update of the surface function.
This could be necessary because these functions overwrite the horizontal components sxx and szz on the surface.
Tests with sofi2D showd that if CPML is not applied again for the horizontal surface points, it could lead to
instabilities.
parent 6202cc7a
......@@ -62,7 +62,7 @@
"Boundary Conditions" : "comment",
"FREE_SURF" : "1",
"ABS_TYPE" : "2",
"ABS_TYPE" : "1",
"FW" : "20.0",
"DAMPING" : "8.0",
"FPML" : "5.0",
......
......@@ -262,11 +262,14 @@ float **splitsrc(float **srcpos,int *nsrc_loc, int nsrc, int * stype_loc, int *s
void surface(int ndepth, float *** u, float *** pi, float ***taus, float *** taup,
float * eta, float *** sxx, float ***syy, float ***szz, float *** sxy,float *** syz,
float *** rxx, float *** ryy, float ***rzz, float *** vx, float *** vy,
float *** vz);
float *** vz, float * K_x, float * a_x, float * b_x, float * K_z, float * a_z, float * b_z,
float *** psi_vxx, float *** psi_vzz );
void surface_elastic(int ndepth, float *** u, float *** pi,
float *** sxx, float ***syy, float ***szz, float *** sxy,float *** syz,
float *** vx, float *** vy, float *** vz);
float *** vx, float *** vy, float *** vz,
float * K_x, float * a_x, float * b_x, float * K_z, float * a_z, float * b_z,
float *** psi_vxx, float *** psi_vzz );
void surface_acoustic(int ndepth, float *** pi, float *** sxx, float *** vx, float *** vy, float *** vz);
......
......@@ -698,6 +698,8 @@ int main(int argc, char **argv){
model_visco(rho,pi,u,taus,taup,eta); /* viscoelastic modeling, L is specified in input file*/
}
}
if (RUN_MULTIPLE_SHOTS) nshots=nsrc; else nshots=1;
/*printf("\n ------------------ checkfd by MYID %i -------------------- \n", MYID);*/
/* check if the FD run will be stable and free of numerical dispersion */
......@@ -941,8 +943,10 @@ int main(int argc, char **argv){
/* stress free surface ? */
if ((FREE_SURF) && (POS[2]==0)){
if (L) surface(1,u,pi,taus,taup,eta,sxx,syy,szz,sxy,syz,rxx,ryy,rzz,vx,vy,vz);
else surface_elastic(1,u,pi,sxx,syy,szz,sxy,syz,vx,vy,vz);
if (L) surface(1,u,pi,taus,taup,eta,sxx,syy,szz,sxy,syz,rxx,ryy,rzz,vx,vy,vz,K_x,a_x,b_x,
K_z,a_z,b_z,psi_vxx, psi_vzz);
else surface_elastic(1,u,pi,sxx,syy,szz,sxy,syz,vx,vy,vz,K_x,a_x,b_x,
K_z,a_z,b_z,psi_vxx, psi_vzz);
}
/* exchange values of stress at boundaries between PEs */
......
......@@ -25,15 +25,17 @@
void surface(int ndepth, float *** u, float *** pi, float ***taus, float *** taup,
float * eta, float *** sxx, float ***syy, float ***szz, float *** sxy,float *** syz,
float *** rxx, float *** ryy, float ***rzz, float *** vx, float *** vy,
float *** vz){
float *** rxx, float *** ryy, float ***rzz, float *** vx, float *** vy, float *** vz,
float * K_x, float * a_x, float * b_x, float * K_z, float * a_z, float * b_z,
float *** psi_vxx, float *** psi_vzz ){
int i, k ,j, fdoh,m;
int i, k ,j, fdoh,m, h1;
float vxx, vyy, vzz;
float b, d, e, f, g, h, dthalbe; /* variables "dh24x, dh24y, dh24z" removed, not in use */
extern int NX, NZ, L, FDORDER, FDCOEFF;
extern int ABS_TYPE, POS[4],FW, NPROCX, NPROCZ;
register float b1, b2, b3, b4, b5, b6;
extern float DT, DX, DY, DZ;
......@@ -83,6 +85,30 @@ float *** vz){
vyy=(-vy[j+1][i][k]+27.0*(vy[j][i][k]-vy[j-1][i][k])+vy[j-2][i][k])/(dh24x);
vzz=(-vz[j][i][k+1]+27.0*(vz[j][i][k]-vz[j][i][k-1])+vz[j][i][k-2])/(dh24x);*/
if (ABS_TYPE==1){
if((POS[1]==0) && (i<=FW)){
psi_vxx[j][i][k] = b_x[i] * psi_vxx[j][i][k] + a_x[i] * vxx;
vxx = vxx / K_x[i] + psi_vxx[j][i][k];
}
if((POS[1]==NPROCX-1) && (i>=NX-FW+1)){
h1 = i-NX+2*FW;
psi_vxx[j][h1][k] = b_x[h1] * psi_vxx[j][h1][k] + a_x[h1] * vxx;
vxx = vxx /K_x[h1] + psi_vxx[j][h1][k];
}
if((POS[3]==0) && (k<=FW)){
psi_vzz[j][i][k] = b_z[k] * psi_vzz[j][i][k] + a_z[k] * vzz;
vzz = vzz / K_z[k] + psi_vzz[j][i][k];
}
if((POS[3]==NPROCZ-1) && (k>=NZ-FW+1)){
h1 = (k-NZ+2*FW);
psi_vzz[j][i][h1] = b_z[h1] * psi_vzz[j][i][h1] + a_z[h1] * vzz;
vzz = vzz / K_z[h1] + psi_vzz[j][i][h1];
}
}
/* partially updating sxx and szz in the same way*/
f=u[j][i][k]*2.0*(1.0+L*taus[j][i][k]);
......@@ -145,6 +171,30 @@ float *** vz){
vyy = (b1*(vy[j][i][k]-vy[j-1][i][k])+b2*(vy[j+1][i][k]-vy[j-2][i][k]))/DY;
vzz = (b1*(vz[j][i][k]-vz[j][i][k-1])+b2*(vz[j][i][k+1]-vz[j][i][k-2]))/DZ;
if (ABS_TYPE==1){
if((POS[1]==0) && (i<=FW)){
psi_vxx[j][i][k] = b_x[i] * psi_vxx[j][i][k] + a_x[i] * vxx;
vxx = vxx / K_x[i] + psi_vxx[j][i][k];
}
if((POS[1]==NPROCX-1) && (i>=NX-FW+1)){
h1 = i-NX+2*FW;
psi_vxx[j][h1][k] = b_x[h1] * psi_vxx[j][h1][k] + a_x[h1] * vxx;
vxx = vxx /K_x[h1] + psi_vxx[j][h1][k];
}
if((POS[3]==0) && (k<=FW)){
psi_vzz[j][i][k] = b_z[k] * psi_vzz[j][i][k] + a_z[k] * vzz;
vzz = vzz / K_z[k] + psi_vzz[j][i][k];
}
if((POS[3]==NPROCZ-1) && (k>=NZ-FW+1)){
h1 = (k-NZ+2*FW);
psi_vzz[j][i][h1] = b_z[h1] * psi_vzz[j][i][h1] + a_z[h1] * vzz;
vzz = vzz / K_z[h1] + psi_vzz[j][i][h1];
}
}
/* partially updating sxx and szz in the same way*/
f=u[j][i][k]*2.0*(1.0+L*taus[j][i][k]);
g=pi[j][i][k]*(1.0+L*taup[j][i][k]);
......
......@@ -25,17 +25,19 @@
void surface_elastic(int ndepth, float *** u, float *** pi,
float *** sxx, float ***syy, float ***szz, float *** sxy,float *** syz,
float *** vx, float *** vy, float *** vz){
float *** vx, float *** vy, float *** vz,
float * K_x, float * a_x, float * b_x, float * K_z, float * a_z, float * b_z,
float *** psi_vxx, float *** psi_vzz ){
int i, k ,j, fdoh,m;
int i, k ,j, fdoh,m,h1;
float vxx, vyy, vzz;
float f, g, h; /* variables "dthalbe, dh24y, dh24z" removed, not in use */
extern int NX, NZ, FDORDER, FDCOEFF;
extern int NX, NZ, FDORDER, FDCOEFF, FW, NPROCX, NPROCZ, POS[4];
register float b1, b2, b3, b4, b5, b6;
extern float DT, DX, DY, DZ;
extern int ABS_TYPE;
j=ndepth; /* The free surface is located exactly in y=(ndepth-1/2)*dh meter!! */
//dthalbe=DT/2.0;
......@@ -62,7 +64,30 @@ void surface_elastic(int ndepth, float *** u, float *** pi,
vxx = (vx[j][i][k]-vx[j][i-1][k])/DX;
vyy = (vy[j][i][k]-vy[j-1][i][k])/DY;
vzz = (vz[j][i][k]-vz[j][i][k-1])/DZ;
if (ABS_TYPE==1){
if((POS[1]==0) && (i<=FW)){
psi_vxx[j][i][k] = b_x[i] * psi_vxx[j][i][k] + a_x[i] * vxx;
vxx = vxx / K_x[i] + psi_vxx[j][i][k];
}
if((POS[1]==NPROCX-1) && (i>=NX-FW+1)){
h1 = i-NX+2*FW;
psi_vxx[j][h1][k] = b_x[h1] * psi_vxx[j][h1][k] + a_x[h1] * vxx;
vxx = vxx /K_x[h1] + psi_vxx[j][h1][k];
}
if((POS[3]==0) && (k<=FW)){
psi_vzz[j][i][k] = b_z[k] * psi_vzz[j][i][k] + a_z[k] * vzz;
vzz = vzz / K_z[k] + psi_vzz[j][i][k];
}
if((POS[3]==NPROCZ-1) && (k>=NZ-FW+1)){
h1 = (k-NZ+2*FW);
psi_vzz[j][i][h1] = b_z[h1] * psi_vzz[j][i][h1] + a_z[h1] * vzz;
vzz = vzz / K_z[h1] + psi_vzz[j][i][h1];
}
}
f=u[j][i][k]*2.0;
g=pi[j][i][k];
......@@ -99,7 +124,30 @@ void surface_elastic(int ndepth, float *** u, float *** pi,
vyy = (b1*(vy[j][i][k]-vy[j-1][i][k])+b2*(vy[j+1][i][k]-vy[j-2][i][k]))/DY;
vzz = (b1*(vz[j][i][k]-vz[j][i][k-1])+b2*(vz[j][i][k+1]-vz[j][i][k-2]))/DZ;
if (ABS_TYPE==1){
if((POS[1]==0) && (i<=FW)){
psi_vxx[j][i][k] = b_x[i] * psi_vxx[j][i][k] + a_x[i] * vxx;
vxx = vxx / K_x[i] + psi_vxx[j][i][k];
}
if((POS[1]==NPROCX-1) && (i>=NX-FW+1)){
h1 = i-NX+2*FW;
psi_vxx[j][h1][k] = b_x[h1] * psi_vxx[j][h1][k] + a_x[h1] * vxx;
vxx = vxx /K_x[h1] + psi_vxx[j][h1][k];
}
if((POS[3]==0) && (k<=FW)){
psi_vzz[j][i][k] = b_z[k] * psi_vzz[j][i][k] + a_z[k] * vzz;
vzz = vzz / K_z[k] + psi_vzz[j][i][k];
}
if((POS[3]==NPROCZ-1) && (k>=NZ-FW+1)){
h1 = (k-NZ+2*FW);
psi_vzz[j][i][h1] = b_z[h1] * psi_vzz[j][i][h1] + a_z[h1] * vzz;
vzz = vzz / K_z[h1] + psi_vzz[j][i][h1];
}
}
f=u[j][i][k]*2.0;
g=pi[j][i][k];
h=-(DT*(g-f)*(g-f)*(vxx+vzz)/g)-(DT*(g-f)*vyy);
......
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