Commit d000f3eb authored by Tilman Steinweg's avatar Tilman Steinweg

started to implement structs in ifos

parent b22757f3
......@@ -24,12 +24,13 @@
#include "fd.h"
double exchange_s(float *** sxx, float *** syy, float *** szz,
double exchange_s(st_stress *stress, float *** syy, float *** szz,
float *** sxy, float *** syz, float *** sxz,
float *** bufferlef_to_rig, float *** bufferrig_to_lef,
float *** buffertop_to_bot, float *** bufferbot_to_top,
float *** bufferfro_to_bac, float *** bufferbac_to_fro, MPI_Request * req_send, MPI_Request * req_rec) {
extern int NX, NY, NZ, POS[4], NPROCX, NPROCY, NPROCZ, BOUNDARY, FDORDER, INDEX[7]; /*MYID,LOG,*/
extern const int TAG1,TAG2,TAG3,TAG4,TAG5,TAG6;
......@@ -157,7 +158,7 @@ float *** bufferfro_to_bac, float *** bufferbac_to_fro, MPI_Request * req_send,
/* storage of left edge of local volume into buffer */
n=1;
for (l=1;l<=FDORDER/2;l++){
bufferlef_to_rig[j][k][n++] = sxx[j][l][k];
bufferlef_to_rig[j][k][n++] = stress->sxx[j][l][k];
}
for (l=1;l<=(FDORDER/2-1);l++){
......@@ -180,7 +181,7 @@ float *** bufferfro_to_bac, float *** bufferbac_to_fro, MPI_Request * req_send,
}
for (l=1;l<=(FDORDER/2-1);l++){
bufferrig_to_lef[j][k][n++] = sxx[j][NX-l+1][k];
bufferrig_to_lef[j][k][n++] = stress->sxx[j][NX-l+1][k];
}
}
}
......@@ -215,7 +216,7 @@ float *** bufferfro_to_bac, float *** bufferbac_to_fro, MPI_Request * req_send,
n=1;
for (l=1;l<=(FDORDER/2);l++){
sxx[j][NX+l][k] = bufferlef_to_rig[j][k][n++];
stress->sxx[j][NX+l][k] = bufferlef_to_rig[j][k][n++];
}
for (l=1;l<=(FDORDER/2-1);l++){
......@@ -236,7 +237,7 @@ float *** bufferfro_to_bac, float *** bufferbac_to_fro, MPI_Request * req_send,
}
for (l=1;l<=(FDORDER/2-1);l++){
sxx[j][1-l][k] = bufferrig_to_lef[j][k][n++];
stress->sxx[j][1-l][k] = bufferrig_to_lef[j][k][n++];
}
}
}
......
......@@ -46,6 +46,12 @@
#define NPROCY_MAX 100
#define NPROCZ_MAX 100
typedef struct Stress{
//stress tensor
// Pointer for dynamic wavefields:
float *** sxx, *** syy, *** szz;
float *** sxy, *** syz, *** sxz;
}st_stress;
/* declaration of functions */
......@@ -127,7 +133,7 @@ int plane_wave(float *** force_points);
void CPML_ini_elastic(int * xb, int * yb, int * zb);
void psource(int nt, float *** sxx, float *** syy, float *** szz, float ** srcpos_loc, float ** signals, int nsrc);
void psource(int nt, st_stress *stress, float *** syy, float *** szz, float **srcpos_loc, float **signals, int nsrc);
void psource_acoustic(int nt, float *** sxx, float ** srcpos_loc, float ** signals, int nsrc, int * stype);
......@@ -148,11 +154,7 @@ float *** bufferlef_to_rig, float *** bufferrig_to_lef,
float *** buffertop_to_bot, float *** bufferbot_to_top,
float *** bufferfro_to_bac, float *** bufferbac_to_fro);
double exchange_s(float *** sxx, float *** syy, float *** szz,
float *** sxy, float *** syz, float *** sxz,
float *** bufferlef_to_rig, float *** bufferrig_to_lef,
float *** buffertop_to_bot, float *** bufferbot_to_top,
float *** bufferfro_to_bac, float *** bufferbac_to_fro, MPI_Request * req_send, MPI_Request * req_rec);
double exchange_s(st_stress *stress, float *** syy, float *** szz, float *** sxy, float *** syz, float *** sxz, float *** bufferlef_to_rig, float *** bufferrig_to_lef, float *** buffertop_to_bot, float *** bufferbot_to_top, float *** bufferfro_to_bac, float *** bufferbac_to_fro, MPI_Request *req_send, MPI_Request *req_rec);
double exchange_s_acoustic(float *** sxx,
float *** bufferlef_to_rig, float *** bufferrig_to_lef,
......@@ -195,10 +197,7 @@ void seismo_acoustic(int lsamp, int ntr, int **recpos, float **sectionvx, float
float **sectionvz, float **sectiondiv, float **sectioncurl, float **sectionp,
float ***vx, float ***vy, float ***vz, float ***sxx, float ***pi);
void seismo(int lsamp, int ntr, int **recpos, float **sectionvx, float **sectionvy,
float **sectionvz, float **sectiondiv, float **sectioncurl, float **sectionp,
float ***vx, float ***vy, float ***vz,
float ***sxx, float ***syy, float ***szz, float ***pi, float ***u);
void seismo(int lsamp, int ntr, int **recpos, float **sectionvx, float **sectionvy, float **sectionvz, float **sectiondiv, float **sectioncurl, float **sectionp, float *** vx, float *** vy, float *** vz, st_stress *stress, float *** syy, float *** szz, float *** pi, float *** u);
void seismo_rsg(int lsamp, int ntr, int **recpos, float **sectionvx, float **sectionvy,
float **sectionvz, float **sectiondiv, float **sectioncurl, float **sectionp,
......@@ -210,10 +209,7 @@ float ***vx, float ***vy, float ***vz, float ***sxx, float ***pi,
int idx, int idy, int idz, int nx1, int ny1, int nz1, int nx2,
int ny2, int nz2);
void snap(FILE *fp, int nt, int nsnap, int format, int type,
float ***vx, float ***vy, float ***vz, float ***sxx, float ***syy, float ***szz, float ***u, float ***pi,
int idx, int idy, int idz, int nx1, int ny1, int nz1, int nx2,
int ny2, int nz2);
void snap(FILE *fp, int nt, int nsnap, int format, int type, float *** vx, float *** vy, float *** vz, st_stress *stress, float *** syy, float *** szz, float *** u, float *** pi, int idx, int idy, int idz, int nx1, int ny1, int nz1, int nx2, int ny2, int nz2);
void snap_rsg(FILE *fp, int nt, int nsnap, int format, int type,
......@@ -234,33 +230,16 @@ 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, 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);
void surface(int ndepth, float *** u, float *** pi, float *** taus, float *** taup, float *eta, st_stress *stress, float *** syy, float *** szz, float *** sxy, float *** syz, float *** rxx, float *** ryy, float *** rzz, float *** vx, float *** vy, float *** vz);
void surface_acoustic(int ndepth, float *** pi, float *** sxx, float *** vx, float *** vy, float *** vz);
/*void timing(double * time_v_update, double * time_s_update, double * time_s_exchange, double * time_v_exchange,
double * time_timestep, int ishot);*/
double update_s(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2,
float *** vx, float *** vy, float *** vz,
float *** sxx, float *** syy, float *** szz, float *** sxy,
float *** syz, float *** sxz, float *** rxx, float *** ryy,
float *** rzz, float *** rxy, float *** ryz, float *** rxz,
float *** pi, float *** u, float *** uipjp, float *** ujpkp, float *** uipkp,
float *** taus, float *** tausipjp, float *** tausjpkp, float *** tausipkp,
float *** taup, float * eta);
double update_s_CPML(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float *** vx, float *** vy, float *** vz,
float *** sxx, float *** syy, float *** szz, float *** sxy, float *** syz, float *** sxz, float *** rxx, float *** ryy,
float *** rzz, float *** rxy, float *** ryz, float *** rxz, float *** pi, float *** u, float *** uipjp, float *** ujpkp, float *** uipkp,
float *** taus, float *** tausipjp, float *** tausjpkp, float *** tausipkp, float *** taup, float * eta,
float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half,
float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half,
float * K_z, float * a_z, float * b_z, float * K_z_half, float * a_z_half, float * b_z_half,
float *** psi_vxx, float *** psi_vyx, float *** psi_vzx, float *** psi_vxy, float *** psi_vyy, float *** psi_vzy, float *** psi_vxz, float *** psi_vyz, float *** psi_vzz);
double update_s(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float *** vx, float *** vy, float *** vz, st_stress *stress, float *** syy, float *** szz, float *** sxy, float *** syz, float *** sxz, float *** rxx, float *** ryy, float *** rzz, float *** rxy, float *** ryz, float *** rxz, float *** pi, float *** u, float *** uipjp, float *** ujpkp, float *** uipkp, float *** taus, float *** tausipjp, float *** tausjpkp, float *** tausipkp, float *** taup, float *eta);
double update_s_CPML(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float *** vx, float *** vy, float *** vz, st_stress *stress, float *** syy, float *** szz, float *** sxy, float *** syz, float *** sxz, float *** rxx, float *** ryy, float *** rzz, float *** rxy, float *** ryz, float *** rxz, float *** pi, float *** u, float *** uipjp, float *** ujpkp, float *** uipkp, float *** taus, float *** tausipjp, float *** tausjpkp, float *** tausipkp, float *** taup, float *eta, float *K_x, float *a_x, float *b_x, float *K_x_half, float *a_x_half, float *b_x_half, float *K_y, float *a_y, float *b_y, float *K_y_half, float *a_y_half, float *b_y_half, float *K_z, float *a_z, float *b_z, float *K_z_half, float *a_z_half, float *b_z_half, float *** psi_vxx, float *** psi_vyx, float *** psi_vzx, float *** psi_vxy, float *** psi_vyy, float *** psi_vzy, float *** psi_vxz, float *** psi_vyz, float *** psi_vzz);
double update_s_acoustic(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2,
float *** vx, float *** vy, float *** vz, float *** sxx, float *** pi);
......@@ -278,22 +257,9 @@ float *** rzz, float *** rxy, float *** ryz, float *** rxz,
float *** pi, float *** u,
float *** taus, float *** taup, float * eta);
double update_v(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2,
int nt, float *** vx, float *** vy, float *** vz,
float *** sxx, float *** syy, float *** szz, float *** sxy,
float *** syz, float *** sxz, float *** rho, float *** rjp, float *** rkp, float *** rip,
float ** srcpos_loc, float ** signals, float ** signaly, float ** signalz, int nsrc, float ***absorb_coeff, int back);
double update_v(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, int nt, float *** vx, float *** vy, float *** vz, st_stress *stress, float *** syy, float *** szz, float *** sxy, float *** syz, float *** sxz, float *** rho, float *** rjp, float *** rkp, float *** rip, float **srcpos_loc, float **signals, float **signaly, float **signalz, int nsrc, float *** absorb_coeff, int back);
double update_v_CPML(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2,
int nt, float *** vx, float *** vy, float *** vz,
float *** sxx, float *** syy, float *** szz, float *** sxy,
float *** syz, float *** sxz, float *** rho, float *** rjp, float *** rkp, float *** rip,
float ** srcpos_loc, float ** signals, float ** signaly, float ** signalz, int nsrc, float *** absorb_coeff, int back,
float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half,
float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half,
float * K_z, float * a_z, float * b_z, float * K_z_half, float * a_z_half, float * b_z_half,
float *** psi_sxx_x, float *** psi_sxy_x, float *** psi_sxz_x, float *** psi_sxy_y, float *** psi_syy_y,
float *** psi_syz_y, float *** psi_sxz_z, float *** psi_syz_z, float *** psi_szz_z);
double update_v_CPML(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, int nt, float *** vx, float *** vy, float *** vz, st_stress *stress, float *** syy, float *** szz, float *** sxy, float *** syz, float *** sxz, float *** rho, float *** rjp, float *** rkp, float *** rip, float **srcpos_loc, float **signals, float **signaly, float **signalz, int nsrc, float *** absorb_coeff, int back, float *K_x, float *a_x, float *b_x, float *K_x_half, float *a_x_half, float *b_x_half, float *K_y, float *a_y, float *b_y, float *K_y_half, float *a_y_half, float *b_y_half, float *K_z, float *a_z, float *b_z, float *K_z_half, float *a_z_half, float *b_z_half, float *** psi_sxx_x, float *** psi_sxy_x, float *** psi_sxz_x, float *** psi_sxy_y, float *** psi_syy_y, float *** psi_syz_y, float *** psi_sxz_z, float *** psi_syz_z, float *** psi_szz_z);
double update_v_acoustic(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2,
int nt, float *** vx, float *** vy, float *** vz,
......@@ -385,7 +351,7 @@ void free_dvector(double *v, int nl, int nh);
void readseis(int ishot, float **section, float **sectionf, int ntr, int ns, int comp);
void readseis_split(int ishot, float **section, int ntr, int *rnum_loc, int ns, int comp);
void residual(float **sectiondata, float **sectiondataf, float **section, float **sectiondiff, int ntr, int ns, float *L2,float *L2f );
void zero_wavefield( int NX, int NY, int NZ, float *** vx, float *** vy, float *** vz,float *** sxx, float *** syy, float *** szz, float *** sxy, float *** syz, float *** sxz, float *** rxx, float *** ryy, float *** rzz, float *** rxy, float *** ryz, float *** rxz, float *** psi_sxx_x, float *** psi_sxy_x, float *** psi_sxz_x, float *** psi_sxy_y, float *** psi_syy_y, float *** psi_syz_y, float *** psi_sxz_z, float *** psi_syz_z, float *** psi_szz_z, float *** psi_vxx, float *** psi_vyx, float *** psi_vzx, float *** psi_vxy, float *** psi_vyy, float *** psi_vzy, float *** psi_vxz, float *** psi_vyz, float *** psi_vzz);
void zero_wavefield( int NX, int NY, int NZ, float *** vx, float *** vy, float *** vz, st_stress *stress, float *** syy, float *** szz, float *** sxy, float *** syz, float *** sxz, float *** rxx, float *** ryy, float *** rzz, float *** rxy, float *** ryz, float *** rxz, float *** psi_sxx_x, float *** psi_sxy_x, float *** psi_sxz_x, float *** psi_sxy_y, float *** psi_syy_y, float *** psi_syz_y, float *** psi_sxz_z, float *** psi_syz_z, float *** psi_szz_z, float *** psi_vxx, float *** psi_vyx, float *** psi_vzx, float *** psi_vxy, float *** psi_vyy, float *** psi_vzy, float *** psi_vxz, float *** psi_vyz, float *** psi_vzz);
void gradient_F(int nx,int ny,int nz,float **** fvx,float **** fvy,float **** fvz,float **** fivx,float **** fivy,float **** fivz,float **** bvx,float **** bvy,float **** bvz,float **** bivx,float **** bivy,float **** bivz,float ***gradlam, float ***gradmu,float ***gradrho, int nt,float *** rho, float *** pi, float *** u, float * finv,int nf, int iteration);
void modelupdate(int nx, int ny, int nz, float ***gradlam, float ***gradmu, float ***gradrho, float *** rho, float *** pi, float *** u, float **bfgsmod1, float step, float *beta, int iteration);
......@@ -399,10 +365,10 @@ void zero_grad(int NX, int NY, int NZ, float *** grad1, float *** grad2, float *
void cpmodel(int nx, int ny, int nz, float ***rho, float ***pi, float ***u,float *** testrho, float *** testpi, float *** testu);
void conjugate(int nx,int ny,int nz, float ***grad1, float ***grad2,float ***grad3, float ***gradprior1, float ***gradprior2, float ***gradprior3, float ***gradprior4, float ***gradprior5, float ***gradprior6,float *beta, int iteration,int cdf);
double update_s_elastic(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float *** vx, float *** vy, float *** vz,float *** sxx, float *** syy, float *** szz, float *** sxy, float *** syz, float *** sxz, float *** pi, float *** u, float *** uipjp, float *** ujpkp, float *** uipkp);
double update_s_elastic(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float *** vx, float *** vy, float *** vz, st_stress *stress, float *** syy, float *** szz, float *** sxy, float *** syz, float *** sxz, float *** pi, float *** u, float *** uipjp, float *** ujpkp, float *** uipkp);
double update_s_CPML_elastic(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float *** vx, float *** vy, float *** vz, float *** sxx, float *** syy, float *** szz, float *** sxy, float *** syz, float *** sxz, float *** pi, float *** u, float *** uipjp, float *** ujpkp, float *** uipkp, float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half, float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half, float * K_z, float * a_z, float * b_z, float * K_z_half, float * a_z_half, float * b_z_half, float *** psi_vxx, float *** psi_vyx, float *** psi_vzx, float *** psi_vxy, float *** psi_vyy, float *** psi_vzy, float *** psi_vxz, float *** psi_vyz, 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);
double update_s_CPML_elastic(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float *** vx, float *** vy, float *** vz, st_stress *stress, float *** syy, float *** szz, float *** sxy, float *** syz, float *** sxz, float *** pi, float *** u, float *** uipjp, float *** ujpkp, float *** uipkp, float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half, float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half, float * K_z, float * a_z, float * b_z, float * K_z_half, float * a_z_half, float * b_z_half, float *** psi_vxx, float *** psi_vyx, float *** psi_vzx, float *** psi_vxy, float *** psi_vyy, float *** psi_vzy, float *** psi_vxz, float *** psi_vyz, float *** psi_vzz);
void surface_elastic(int ndepth, float *** u, float *** pi, st_stress *stress, float *** syy, float *** szz, float *** sxy, float *** syz, float *** vx, float *** vy, float *** vz);
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, float *** vx, float *** vy, float *** vz, float **** Fvx, float **** Fvy,float **** Fvz, float **** Fvxi, float **** Fvyi,float **** Fvzi,float *finv, int nf,int ntast, int pshot, int back);
......
This diff is collapsed.
......@@ -23,10 +23,10 @@
#include "fd.h"
void psource(int nt, float *** sxx, float *** syy, float *** szz,
void psource(int nt, st_stress *stress, float *** syy, float *** szz,
float ** srcpos_loc, float ** signals, int nsrc){
extern float DX, DY, DZ;
extern int NT;
int i, j, k, l;
......@@ -46,7 +46,7 @@ extern int NT;
if(nt==1){amp=signals[l][nt+1]/(2.0*DX*DY*DZ);}
if((nt>1)&&(nt<NT)){amp=(signals[l][nt+1]-signals[l][nt-1])/(2.0*DX*DY*DZ);}
if(nt==NT){amp=-signals[l][nt-1]/(2.0*DX*DY*DZ);}
sxx[j][i][k]+=amp;
stress->sxx[j][i][k]+=amp;
syy[j][i][k]+=amp;
szz[j][i][k]+=amp;
}
......
......@@ -27,8 +27,9 @@
void seismo(int lsamp, int ntr, int **recpos, float **sectionvx, float **sectionvy,
float **sectionvz, float **sectiondiv, float **sectioncurl, float **sectionp,
float ***vx, float ***vy, float ***vz,
float ***sxx, float ***syy, float ***szz, float ***pi, float ***u){
st_stress *stress, float ***syy, float ***szz, float ***pi, float ***u){
extern int SEISMO;
int i, j, k, itr, ins, nxrec, nyrec, nzrec;
float amp, dh24x, dh24y, dh24z, vyx, vxy, vxx, vyy, vzx, vyz, vxz, vzy, vzz;
......@@ -52,7 +53,7 @@ float ***sxx, float ***syy, float ***szz, float ***pi, float ***u){
sectionvy[itr][ins]=vy[nyrec][nxrec][nzrec];
sectionvz[itr][ins]=vz[nyrec][nxrec][nzrec];
break;
case 2 : sectionp[itr][ins]=-sxx[nyrec][nxrec][nzrec]
case 2 : sectionp[itr][ins]=-stress->sxx[nyrec][nxrec][nzrec]
-syy[nyrec][nxrec][nzrec]
-szz[nyrec][nxrec][nzrec];
break;
......@@ -92,7 +93,7 @@ float ***sxx, float ***syy, float ***szz, float ***pi, float ***u){
sectionvx[itr][ins]=vx[nyrec][nxrec][nzrec];
sectionvy[itr][ins]=vy[nyrec][nxrec][nzrec];
sectionvz[itr][ins]=vz[nyrec][nxrec][nzrec];
sectionp[itr][ins]=-sxx[nyrec][nxrec][nzrec]-syy[nyrec][nxrec][nzrec]-szz[nyrec][nxrec][nzrec];
sectionp[itr][ins]=-stress->sxx[nyrec][nxrec][nzrec]-syy[nyrec][nxrec][nzrec]-szz[nyrec][nxrec][nzrec];
i=nxrec; j=nyrec; k=nzrec;
/*vxy=(-vx[j+2][i][k]+27.0*(vx[j+1][i][k]-vx[j][i][k])+vx[j-1][i][k])*(dh24y);
......
......@@ -25,7 +25,7 @@
void snap(FILE *fp, int nt, int nsnap, int format, int type,
float ***vx, float ***vy, float ***vz, float ***sxx, float ***syy, float ***szz,
float ***vx, float ***vy, float ***vz, st_stress *stress, float ***syy, float ***szz,
float ***u, float ***pi,
int idx, int idy, int idz, int nx1, int ny1, int nz1, int nx2,
int ny2, int nz2){
......@@ -121,7 +121,7 @@ int ny2, int nz2){
for (k=ny1;k<=ny2;k+=idy)
for (i=nx1;i<=nx2;i+=idx)
for (j=nz1;j<=nz2;j+=idz){
amp=-sxx[j][i][k]-syy[j][i][k]-szz[j][i][k];
amp=-stress->sxx[j][i][k]-syy[j][i][k]-szz[j][i][k];
writedsk(fpx1,amp,format);
......
......@@ -24,7 +24,7 @@
#include "fd.h"
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 * eta, st_stress *stress, float ***syy, float ***szz, float *** sxy,float *** syz,
float *** rxx, float *** ryy, float ***rzz, float *** vx, float *** vy,
float *** vz){
......@@ -89,7 +89,7 @@ float *** vz){
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]);
h=-(DT*(g-f)*(g-f)*(vxx+vzz)/g)-(DT*(g-f)*vyy);
sxx[j][i][k]+=h-(dthalbe*rxx[j][i][k]);
stress->sxx[j][i][k]+=h-(dthalbe*rxx[j][i][k]);
szz[j][i][k]+=h-(dthalbe*rzz[j][i][k]);
/* updating the memory-variable rxx, rzz at the free surface */
......@@ -101,7 +101,7 @@ float *** vz){
rzz[j][i][k]+=h;
/*completely updating the stresses sxx and szz */
sxx[j][i][k]+=(dthalbe*rxx[j][i][k]);
stress->sxx[j][i][k]+=(dthalbe*rxx[j][i][k]);
szz[j][i][k]+=(dthalbe*rzz[j][i][k]);
......@@ -150,7 +150,7 @@ float *** vz){
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]);
h=-(DT*(g-f)*(g-f)*(vxx+vzz)/g)-(DT*(g-f)*vyy);
sxx[j][i][k]+=h-(dthalbe*rxx[j][i][k]);
stress->sxx[j][i][k]+=h-(dthalbe*rxx[j][i][k]);
szz[j][i][k]+=h-(dthalbe*rzz[j][i][k]);
/* updating the memory-variable rxx, rzz at the free surface */
......@@ -162,7 +162,7 @@ float *** vz){
rzz[j][i][k]+=h;
/*completely updating the stresses sxx and szz */
sxx[j][i][k]+=(dthalbe*rxx[j][i][k]);
stress->sxx[j][i][k]+=(dthalbe*rxx[j][i][k]);
szz[j][i][k]+=(dthalbe*rzz[j][i][k]);
}
}
......@@ -216,7 +216,7 @@ float *** vz){
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]);
h=-(DT*(g-f)*(g-f)*(vxx+vzz)/g)-(DT*(g-f)*vyy);
sxx[j][i][k]+=h-(dthalbe*rxx[j][i][k]);
stress->sxx[j][i][k]+=h-(dthalbe*rxx[j][i][k]);
szz[j][i][k]+=h-(dthalbe*rzz[j][i][k]);
/* updating the memory-variable rxx, rzz at the free surface */
......@@ -228,7 +228,7 @@ float *** vz){
rzz[j][i][k]+=h;
/*completely updating the stresses sxx and szz */
sxx[j][i][k]+=(dthalbe*rxx[j][i][k]);
stress->sxx[j][i][k]+=(dthalbe*rxx[j][i][k]);
szz[j][i][k]+=(dthalbe*rzz[j][i][k]);
}
}
......@@ -286,7 +286,7 @@ float *** vz){
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]);
h=-(DT*(g-f)*(g-f)*(vxx+vzz)/g)-(DT*(g-f)*vyy);
sxx[j][i][k]+=h-(dthalbe*rxx[j][i][k]);
stress->sxx[j][i][k]+=h-(dthalbe*rxx[j][i][k]);
szz[j][i][k]+=h-(dthalbe*rzz[j][i][k]);
/* updating the memory-variable rxx, rzz at the free surface */
......@@ -298,7 +298,7 @@ float *** vz){
rzz[j][i][k]+=h;
/*completely updating the stresses sxx and szz */
sxx[j][i][k]+=(dthalbe*rxx[j][i][k]);
stress->sxx[j][i][k]+=(dthalbe*rxx[j][i][k]);
szz[j][i][k]+=(dthalbe*rzz[j][i][k]);
}
}
......@@ -359,7 +359,7 @@ float *** vz){
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]);
h=-(DT*(g-f)*(g-f)*(vxx+vzz)/g)-(DT*(g-f)*vyy);
sxx[j][i][k]+=h-(dthalbe*rxx[j][i][k]);
stress->sxx[j][i][k]+=h-(dthalbe*rxx[j][i][k]);
szz[j][i][k]+=h-(dthalbe*rzz[j][i][k]);
/* updating the memory-variable rxx, rzz at the free surface */
......@@ -371,7 +371,7 @@ float *** vz){
rzz[j][i][k]+=h;
/*completely updating the stresses sxx and szz */
sxx[j][i][k]+=(dthalbe*rxx[j][i][k]);
stress->sxx[j][i][k]+=(dthalbe*rxx[j][i][k]);
szz[j][i][k]+=(dthalbe*rzz[j][i][k]);
}
}
......@@ -438,7 +438,7 @@ float *** vz){
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]);
h=-(DT*(g-f)*(g-f)*(vxx+vzz)/g)-(DT*(g-f)*vyy);
sxx[j][i][k]+=h-(dthalbe*rxx[j][i][k]);
stress->sxx[j][i][k]+=h-(dthalbe*rxx[j][i][k]);
szz[j][i][k]+=h-(dthalbe*rzz[j][i][k]);
/* updating the memory-variable rxx, rzz at the free surface */
......@@ -450,7 +450,7 @@ float *** vz){
rzz[j][i][k]+=h;
/*completely updating the stresses sxx and szz */
sxx[j][i][k]+=(dthalbe*rxx[j][i][k]);
stress->sxx[j][i][k]+=(dthalbe*rxx[j][i][k]);
szz[j][i][k]+=(dthalbe*rzz[j][i][k]);
}
}
......
......@@ -24,9 +24,10 @@
#include "fd.h"
void surface_elastic(int ndepth, float *** u, float *** pi,
float *** sxx, float ***syy, float ***szz, float *** sxy,float *** syz,
st_stress *stress, float ***syy, float ***szz, float *** sxy,float *** syz,
float *** vx, float *** vy, float *** vz){
int i, k ,j, fdoh,m;
float vxx, vyy, vzz;
float f, g, h;
......@@ -67,7 +68,7 @@ void surface_elastic(int ndepth, float *** u, float *** pi,
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);
sxx[j][i][k]+=h;
stress->sxx[j][i][k]+=h;
szz[j][i][k]+=h;
......@@ -103,7 +104,7 @@ void surface_elastic(int ndepth, float *** u, float *** pi,
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);
sxx[j][i][k]+=h;
stress->sxx[j][i][k]+=h;
szz[j][i][k]+=h;
......@@ -148,7 +149,7 @@ void surface_elastic(int ndepth, float *** u, float *** pi,
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);
sxx[j][i][k]+=h;
stress->sxx[j][i][k]+=h;
szz[j][i][k]+=h;
}
}
......@@ -195,7 +196,7 @@ void surface_elastic(int ndepth, float *** u, float *** pi,
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);
sxx[j][i][k]+=h;
stress->sxx[j][i][k]+=h;
szz[j][i][k]+=h;
......@@ -246,7 +247,7 @@ void surface_elastic(int ndepth, float *** u, float *** pi,
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);
sxx[j][i][k]+=h;
stress->sxx[j][i][k]+=h;
szz[j][i][k]+=h;
}
}
......@@ -301,7 +302,7 @@ void surface_elastic(int ndepth, float *** u, float *** pi,
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);
sxx[j][i][k]+=h;
stress->sxx[j][i][k]+=h;
szz[j][i][k]+=h;
}
}
......
......@@ -26,7 +26,7 @@
#include "fd.h"
double update_s(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float *** vx, float *** vy, float *** vz,
float *** sxx, float *** syy, float *** szz, float *** sxy, float *** syz, float *** sxz, float *** rxx, float *** ryy,
st_stress *stress, float *** syy, float *** szz, float *** sxy, float *** syz, float *** sxz, float *** rxx, float *** ryy,
float *** rzz, float *** rxy, float *** ryz, float *** rxz, float *** pi, float *** u, float *** uipjp, float *** ujpkp, float *** uipkp,
float *** taus, float *** tausipjp, float *** tausjpkp, float *** tausipkp, float *** taup, float * eta){
......@@ -101,7 +101,7 @@ float *** taus, float *** tausipjp, float *** tausjpkp, float *** tausip
sxy[j][i][k]+=(fipjp*vxyyx)+(dthalbe*sumrxy);
syz[j][i][k]+=(fjpkp*vyzzy)+(dthalbe*sumryz);
sxz[j][i][k]+=(fipkp*vxzzx)+(dthalbe*sumrxz);
sxx[j][i][k]+=DT*((g*vxxyyzz)-(f*vyyzz))+(dthalbe*sumrxx);
stress->sxx[j][i][k]+=DT*((g*vxxyyzz)-(f*vyyzz))+(dthalbe*sumrxx);
syy[j][i][k]+=DT*((g*vxxyyzz)-(f*vxxzz))+(dthalbe*sumryy);
szz[j][i][k]+=DT*((g*vxxyyzz)-(f*vxxyy))+(dthalbe*sumrzz);
......@@ -135,7 +135,7 @@ float *** taus, float *** tausipjp, float *** tausjpkp, float *** tausip
sxy[j][i][k]+=(dthalbe*sumrxy);
syz[j][i][k]+=(dthalbe*sumryz);
sxz[j][i][k]+=(dthalbe*sumrxz);
sxx[j][i][k]+=(dthalbe*sumrxx);
stress->sxx[j][i][k]+=(dthalbe*sumrxx);
syy[j][i][k]+=(dthalbe*sumryy);
szz[j][i][k]+=(dthalbe*sumrzz);
......@@ -202,7 +202,7 @@ break;
sxy[j][i][k]+=(fipjp*vxyyx)+(dthalbe*sumrxy);
syz[j][i][k]+=(fjpkp*vyzzy)+(dthalbe*sumryz);
sxz[j][i][k]+=(fipkp*vxzzx)+(dthalbe*sumrxz);
sxx[j][i][k]+=DT*((g*vxxyyzz)-(f*vyyzz))+(dthalbe*sumrxx);
stress->sxx[j][i][k]+=DT*((g*vxxyyzz)-(f*vyyzz))+(dthalbe*sumrxx);
syy[j][i][k]+=DT*((g*vxxyyzz)-(f*vxxzz))+(dthalbe*sumryy);
szz[j][i][k]+=DT*((g*vxxyyzz)-(f*vxxyy))+(dthalbe*sumrzz);
......@@ -236,7 +236,7 @@ break;
sxy[j][i][k]+=(dthalbe*sumrxy);
syz[j][i][k]+=(dthalbe*sumryz);
sxz[j][i][k]+=(dthalbe*sumrxz);
sxx[j][i][k]+=(dthalbe*sumrxx);
stress->sxx[j][i][k]+=(dthalbe*sumrxx);
syy[j][i][k]+=(dthalbe*sumryy);
szz[j][i][k]+=(dthalbe*sumrzz);
......@@ -329,7 +329,7 @@ break;
sxy[j][i][k]+=(fipjp*vxyyx)+(dthalbe*sumrxy);
syz[j][i][k]+=(fjpkp*vyzzy)+(dthalbe*sumryz);
sxz[j][i][k]+=(fipkp*vxzzx)+(dthalbe*sumrxz);
sxx[j][i][k]+=DT*((g*vxxyyzz)-(f*vyyzz))+(dthalbe*sumrxx);
stress->sxx[j][i][k]+=DT*((g*vxxyyzz)-(f*vyyzz))+(dthalbe*sumrxx);
syy[j][i][k]+=DT*((g*vxxyyzz)-(f*vxxzz))+(dthalbe*sumryy);
szz[j][i][k]+=DT*((g*vxxyyzz)-(f*vxxyy))+(dthalbe*sumrzz);
......@@ -363,7 +363,7 @@ break;
sxy[j][i][k]+=(dthalbe*sumrxy);
syz[j][i][k]+=(dthalbe*sumryz);
sxz[j][i][k]+=(dthalbe*sumrxz);
sxx[j][i][k]+=(dthalbe*sumrxx);
stress->sxx[j][i][k]+=(dthalbe*sumrxx);
syy[j][i][k]+=(dthalbe*sumryy);
szz[j][i][k]+=(dthalbe*sumrzz);
......@@ -465,7 +465,7 @@ break;
sxy[j][i][k]+=(fipjp*vxyyx)+(dthalbe*sumrxy);
syz[j][i][k]+=(fjpkp*vyzzy)+(dthalbe*sumryz);
sxz[j][i][k]+=(fipkp*vxzzx)+(dthalbe*sumrxz);
sxx[j][i][k]+=DT*((g*vxxyyzz)-(f*vyyzz))+(dthalbe*sumrxx);
stress->sxx[j][i][k]+=DT*((g*vxxyyzz)-(f*vyyzz))+(dthalbe*sumrxx);
syy[j][i][k]+=DT*((g*vxxyyzz)-(f*vxxzz))+(dthalbe*sumryy);
szz[j][i][k]+=DT*((g*vxxyyzz)-(f*vxxyy))+(dthalbe*sumrzz);
......@@ -499,7 +499,7 @@ break;
sxy[j][i][k]+=(dthalbe*sumrxy);
syz[j][i][k]+=(dthalbe*sumryz);
sxz[j][i][k]+=(dthalbe*sumrxz);
sxx[j][i][k]+=(dthalbe*sumrxx);
stress->sxx[j][i][k]+=(dthalbe*sumrxx);
syy[j][i][k]+=(dthalbe*sumryy);
szz[j][i][k]+=(dthalbe*sumrzz);
......@@ -610,7 +610,7 @@ case 10 :
sxy[j][i][k]+=(fipjp*vxyyx)+(dthalbe*sumrxy);
syz[j][i][k]+=(fjpkp*vyzzy)+(dthalbe*sumryz);
sxz[j][i][k]+=(fipkp*vxzzx)+(dthalbe*sumrxz);
sxx[j][i][k]+=DT*((g*vxxyyzz)-(f*vyyzz))+(dthalbe*sumrxx);
stress->sxx[j][i][k]+=DT*((g*vxxyyzz)-(f*vyyzz))+(dthalbe*sumrxx);
syy[j][i][k]+=DT*((g*vxxyyzz)-(f*vxxzz))+(dthalbe*sumryy);
szz[j][i][k]+=DT*((g*vxxyyzz)-(f*vxxyy))+(dthalbe*sumrzz);
......@@ -644,7 +644,7 @@ case 10 :
sxy[j][i][k]+=(dthalbe*sumrxy);
syz[j][i][k]+=(dthalbe*sumryz);
sxz[j][i][k]+=(dthalbe*sumrxz);
sxx[j][i][k]+=(dthalbe*sumrxx);
stress->sxx[j][i][k]+=(dthalbe*sumrxx);
syy[j][i][k]+=(dthalbe*sumryy);
szz[j][i][k]+=(dthalbe*sumrzz);
...