Commit 573b32d7 authored by tilman.metz's avatar tilman.metz

fixed instability issues with the free surface

the free surface functions recalculate the updates on the first gridpoint.
PML needs to be applied here ad the edges.
parent a63f40a9
......@@ -264,7 +264,7 @@ int **splitrec(int **recpos,int *ntr_loc, int ntr, int *rnum_loc);
float **splitsrc(float **srcpos,int *nsrc_loc, int nsrc, int *snum_loc);
void surface(int ndepth, st_model *mod, st_stress *stress, st_visc_mem *mem, st_velocity *vel);
void surface(int ndepth, st_model *mod, st_pml_coeff *pml_coeff, st_pml_wfd *pml_wfd, st_visc_mem *mem, st_stress *stress, st_velocity *vel);
void surface_acoustic(int ndepth, st_model *mod, st_stress *stress, st_velocity *vel);
......@@ -383,7 +383,7 @@ void conjugate(int nx, int ny, int nz, st_gradient *grad, st_gradient *grad_prio
double update_s_elastic(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, st_velocity *vel, st_stress *stress, st_model *mod, st_model_av *mod_av);
double update_s_CPML_elastic(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, st_velocity *vel, st_stress *stress, st_model *mod, st_model_av *mod_av, st_pml_coeff *pml_coeff, st_pml_wfd *pml_wfd);
void surface_elastic(int ndepth, st_model *mod, st_stress *stress, st_velocity *vel);
void surface_elastic(int ndepth, st_model *mod, st_pml_coeff *pml_coeff, st_pml_wfd *pml_wfd, st_stress *stress, st_velocity *vel);
void readinv(float *finv, int *nf, int *groupnum,int *itpergroup,int nfmax);
void discfourier(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, int nt, st_velocity *vel, st_freq_velocity *fourier_vel, float *finv, int nf, int ntast, int pshot1, int back);
......
......@@ -810,8 +810,8 @@ MPI_Barrier(MPI_COMM_WORLD);
/* stress free surface ? */
if ((FREE_SURF) && (POS[2]==0))
surface_elastic(1,mod,stress,vel);
/*surface(1,mod,stress,visco_mem,vx,vy,vz);*/
surface_elastic(1,mod,pml_coeff,pml_wfd,stress,vel);
/*surface(1,mod,pml_coeff,pml_wfd,visco_mem,stress,vel);*/
/* store amplitudes at receivers in sectionvx-sectionvz */
if ((SEISMO) && (ntr>0) && (nt==lsamp)){
......@@ -956,8 +956,8 @@ MPI_Barrier(MPI_COMM_WORLD);
/* stress free surface ? */
if ((FREE_SURF) && (POS[2]==0))
/*surface(1,mod,stress,visco_mem,vel);*/
surface_elastic(1,mod,stress,vel);
/*surface(1,mod,pml_coeff,pml_wfd,stress,visco_mem,vel);*/
surface_elastic(1,mod,pml_coeff,pml_wfd,stress,vel);
/* store amplitudes at receivers in sectionvx-sectionvz */
if ((SEISMO) && (ntr>0) && (nt==lsamp)){
seismo(nlsamp,ntr,acq->recpos_loc,section,vel,stress,mod);
......@@ -1106,8 +1106,8 @@ MPI_Barrier(MPI_COMM_WORLD);
time_s_exchange[nt]+=exchange_s(stress,stressbuff, sreq_send, sreq_rec);
if ((FREE_SURF) && (POS[2]==0))
/*surface(1,testmod,stress,visco_mem,vel);*/
surface_elastic(1,testmod,stress,vel);
/*surface(1,testmod,pml_coeff,pml_wfd,stress,visco_mem,vel);*/
surface_elastic(1,testmod,pml_coeff,pml_wfd,stress,vel);
if ((SEISMO) && (ntr>0) && (nt==lsamp)){
seismo(nlsamp,ntr,acq->recpos_loc,section,vel,stress,testmod);
nlsamp++;
......
......@@ -23,15 +23,15 @@
#include "fd.h"
void surface(int ndepth, st_model *mod, st_stress *stress,
st_visc_mem *mem, st_velocity *vel){
void surface(int ndepth, st_model *mod, st_pml_coeff *pml_coeff,st_pml_wfd *pml_wfd, st_visc_mem *mem, st_stress *stress,
st_velocity *vel){
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; /*dh24x*/
extern int NX, NZ, L, FDORDER, FDCOEFF;
extern int NX, NZ, L, FDORDER, FW, ABS_TYPE, FDCOEFF, POS[4],NPROCX,NPROCZ;
register float b1, b2, b3, b4, b5, b6;
extern float DT, DX, DY, DZ;
......@@ -82,6 +82,29 @@ st_visc_mem *mem, st_velocity *vel){
vyy=(-vel->vy[j+1][i][k]+27.0*(vel->vy[j][i][k]-vel->vy[j-1][i][k])+vel->vy[j-2][i][k])/(dh24x);
vzz=(-vel->vz[j][i][k+1]+27.0*(vel->vz[j][i][k]-vel->vz[j][i][k-1])+vel->vz[j][i][k-2])/(dh24x);*/
if (ABS_TYPE==1){
if((POS[1]==0) && (i<=FW)){
pml_wfd->psi_vxx[j][i][k] = pml_coeff->b_x[i] * pml_wfd->psi_vxx[j][i][k] + pml_coeff->a_x[i] * vxx;
vxx = vxx / pml_coeff->K_x[i] + pml_wfd->psi_vxx[j][i][k];
}
if((POS[1]==NPROCX-1) && (i>=NX-FW+1)){
h1 = i-NX+2*FW;
pml_wfd->psi_vxx[j][h1][k] = pml_coeff->b_x[h1] * pml_wfd->psi_vxx[j][h1][k] + pml_coeff->a_x[h1] * vxx;
vxx = vxx /pml_coeff->K_x[h1] + pml_wfd->psi_vxx[j][h1][k];
}
if((POS[3]==0) && (k<=FW)){
pml_wfd->psi_vzz[j][i][k] = pml_coeff->b_z[k] * pml_wfd->psi_vzz[j][i][k] + pml_coeff->a_z[k] * vzz;
vzz = vzz / pml_coeff->K_z[k] + pml_wfd->psi_vzz[j][i][k];
}
if((POS[3]==NPROCZ-1) && (k>=NZ-FW+1)){
h1 = (k-NZ+2*FW);
pml_wfd->psi_vzz[j][i][h1] = pml_coeff->b_z[h1] * pml_wfd->psi_vzz[j][i][h1] + pml_coeff->a_z[h1] * vzz;
vzz = vzz / pml_coeff->K_z[h1] + pml_wfd->psi_vzz[j][i][h1];
}
}
/* partially updating sxx and szz in the same way*/
f=mod->u[j][i][k]*2.0*(1.0+L*mod->taus[j][i][k]);
......@@ -144,6 +167,30 @@ st_visc_mem *mem, st_velocity *vel){
vyy = (b1*(vel->vy[j][i][k]-vel->vy[j-1][i][k])+b2*(vel->vy[j+1][i][k]-vel->vy[j-2][i][k]))/DY;
vzz = (b1*(vel->vz[j][i][k]-vel->vz[j][i][k-1])+b2*(vel->vz[j][i][k+1]-vel->vz[j][i][k-2]))/DZ;
if (ABS_TYPE==1){
if((POS[1]==0) && (i<=FW)){
pml_wfd->psi_vxx[j][i][k] = pml_coeff->b_x[i] * pml_wfd->psi_vxx[j][i][k] + pml_coeff->a_x[i] * vxx;
vxx = vxx / pml_coeff->K_x[i] + pml_wfd->psi_vxx[j][i][k];
}
if((POS[1]==NPROCX-1) && (i>=NX-FW+1)){
h1 = i-NX+2*FW;
pml_wfd->psi_vxx[j][h1][k] = pml_coeff->b_x[h1] * pml_wfd->psi_vxx[j][h1][k] + pml_coeff->a_x[h1] * vxx;
vxx = vxx /pml_coeff->K_x[h1] + pml_wfd->psi_vxx[j][h1][k];
}
if((POS[3]==0) && (k<=FW)){
pml_wfd->psi_vzz[j][i][k] = pml_coeff->b_z[k] * pml_wfd->psi_vzz[j][i][k] + pml_coeff->a_z[k] * vzz;
vzz = vzz / pml_coeff->K_z[k] + pml_wfd->psi_vzz[j][i][k];
}
if((POS[3]==NPROCZ-1) && (k>=NZ-FW+1)){
h1 = (k-NZ+2*FW);
pml_wfd->psi_vzz[j][i][h1] = pml_coeff->b_z[h1] * pml_wfd->psi_vzz[j][i][h1] + pml_coeff->a_z[h1] * vzz;
vzz = vzz / pml_coeff->K_z[h1] + pml_wfd->psi_vzz[j][i][h1];
}
}
/* partially updating sxx and szz in the same way*/
f=mod->u[j][i][k]*2.0*(1.0+L*mod->taus[j][i][k]);
g=mod->pi[j][i][k]*(1.0+L*mod->taup[j][i][k]);
......
......@@ -23,16 +23,16 @@
#include "fd.h"
void surface_elastic(int ndepth, st_model *mod,
void surface_elastic(int ndepth, st_model *mod, st_pml_coeff *pml_coeff,st_pml_wfd *pml_wfd,
st_stress *stress, st_velocity *vel){
int i, k ,j, fdoh,m;
int i, k ,j, fdoh,m,h1;
float vxx, vyy, vzz;
float f, g, h;
extern int NX, NZ, FDORDER, FDCOEFF;
extern int NX, NZ, FDORDER, FDCOEFF,FW,NPROCX,NPROCZ, ABS_TYPE, POS[4];
register float b1, b2, b3, b4, b5, b6;
extern float DT, DX, DY, DZ;
......@@ -62,7 +62,30 @@ void surface_elastic(int ndepth, st_model *mod,
vxx = (vel->vx[j][i][k]-vel->vx[j][i-1][k])/DX;
vyy = (vel->vy[j][i][k]-vel->vy[j-1][i][k])/DY;
vzz = (vel->vz[j][i][k]-vel->vz[j][i][k-1])/DZ;
if (ABS_TYPE==1){
if((POS[1]==0) && (i<=FW)){
pml_wfd->psi_vxx[j][i][k] = pml_coeff->b_x[i] * pml_wfd->psi_vxx[j][i][k] + pml_coeff->a_x[i] * vxx;
vxx = vxx / pml_coeff->K_x[i] + pml_wfd->psi_vxx[j][i][k];
}
if((POS[1]==NPROCX-1) && (i>=NX-FW+1)){
h1 = i-NX+2*FW;
pml_wfd->psi_vxx[j][h1][k] = pml_coeff->b_x[h1] * pml_wfd->psi_vxx[j][h1][k] + pml_coeff->a_x[h1] * vxx;
vxx = vxx /pml_coeff->K_x[h1] + pml_wfd->psi_vxx[j][h1][k];
}
if((POS[3]==0) && (k<=FW)){
pml_wfd->psi_vzz[j][i][k] = pml_coeff->b_z[k] * pml_wfd->psi_vzz[j][i][k] + pml_coeff->a_z[k] * vzz;
vzz = vzz / pml_coeff->K_z[k] + pml_wfd->psi_vzz[j][i][k];
}
if((POS[3]==NPROCZ-1) && (k>=NZ-FW+1)){
h1 = (k-NZ+2*FW);
pml_wfd->psi_vzz[j][i][h1] = pml_coeff->b_z[h1] * pml_wfd->psi_vzz[j][i][h1] + pml_coeff->a_z[h1] * vzz;
vzz = vzz / pml_coeff->K_z[h1] + pml_wfd->psi_vzz[j][i][h1];
}
}
f=mod->u[j][i][k]*2.0;
g=mod->pi[j][i][k];
......@@ -99,7 +122,31 @@ void surface_elastic(int ndepth, st_model *mod,
vyy = (b1*(vel->vy[j][i][k]-vel->vy[j-1][i][k])+b2*(vel->vy[j+1][i][k]-vel->vy[j-2][i][k]))/DY;
vzz = (b1*(vel->vz[j][i][k]-vel->vz[j][i][k-1])+b2*(vel->vz[j][i][k+1]-vel->vz[j][i][k-2]))/DZ;
if (ABS_TYPE==1){
if((POS[1]==0) && (i<=FW)){
pml_wfd->psi_vxx[j][i][k] = pml_coeff->b_x[i] * pml_wfd->psi_vxx[j][i][k] + pml_coeff->a_x[i] * vxx;
vxx = vxx / pml_coeff->K_x[i] + pml_wfd->psi_vxx[j][i][k];
}
if((POS[1]==NPROCX-1) && (i>=NX-FW+1)){
h1 = i-NX+2*FW;
pml_wfd->psi_vxx[j][h1][k] = pml_coeff->b_x[h1] * pml_wfd->psi_vxx[j][h1][k] + pml_coeff->a_x[h1] * vxx;
vxx = vxx /pml_coeff->K_x[h1] + pml_wfd->psi_vxx[j][h1][k];
}
if((POS[3]==0) && (k<=FW)){
pml_wfd->psi_vzz[j][i][k] = pml_coeff->b_z[k] * pml_wfd->psi_vzz[j][i][k] + pml_coeff->a_z[k] * vzz;
vzz = vzz / pml_coeff->K_z[k] + pml_wfd->psi_vzz[j][i][k];
}
if((POS[3]==NPROCZ-1) && (k>=NZ-FW+1)){
h1 = (k-NZ+2*FW);
pml_wfd->psi_vzz[j][i][h1] = pml_coeff->b_z[h1] * pml_wfd->psi_vzz[j][i][h1] + pml_coeff->a_z[h1] * vzz;
vzz = vzz / pml_coeff->K_z[h1] + pml_wfd->psi_vzz[j][i][h1];
}
}
f=mod->u[j][i][k]*2.0;
g=mod->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