...
 
Commits (3)
......@@ -295,7 +295,7 @@ void PML_pro(float * d_x, float * K_x, float * alpha_prime_x, float * a_x, float
if(abscissa_in_PML >= 0.0){
abscissa_normalized = abscissa_in_PML / thickness_PML_y;
d_y_half[h] = d0_x * pow(abscissa_normalized,NPOWER);
d_y_half[h] = d0_y * pow(abscissa_normalized,NPOWER);
/* this taken from Gedney page 8.2 */
K_y_half[h] = 1.0 + (K_MAX_CPML - 1.0) * pow(abscissa_normalized,NPOWER);
......
......@@ -28,7 +28,7 @@ void model_visco(float ** rho, float ** pi, float ** u, float ** taus, floa
/*--------------------------------------------------------------------------*/
/* extern variables */
extern float DT, *FL, TAU, DH;
extern float DT, *FL, DH;
extern int NX, NY, NXG, NYG, POS[3], L, MYID;
extern int WRITE_MODELFILES;
extern char MFILE[STRING_SIZE];
......@@ -43,10 +43,10 @@ void model_visco(float ** rho, float ** pi, float ** u, float ** taus, floa
/*-----------------material property definition -------------------------*/
/* parameters for layer 1 */
const float vp1=3500.0, vs1=2000.0, rho1=2000.0, h=10000.0;
const float vp1=3500.0, vs1=2000.0, rho1=2000.0, tp1=0.1, ts1=0.1, h=2000.0;
/* parameters for layer 2 */
const float vp2=5400.0, vs2=3700.0, rho2=2500.0;
const float vp2=4550.0, vs2=2600.0, rho2=2600.0,tp2=0.1, ts2=0.1;
/*-----------------------------------------------------------------------*/
......@@ -63,17 +63,10 @@ void model_visco(float ** rho, float ** pi, float ** u, float ** taus, floa
eta[l]=DT/pts[l];
}
ts=TAU;
tp=TAU;
ws=2.0*PI*FL[1];
sumu=0.0;
sumpi=0.0;
for (l=1;l<=L;l++){
sumu=sumu+((ws*ws*pts[l]*pts[l]*ts)/(1.0+ws*ws*pts[l]*pts[l]));
sumpi=sumpi+((ws*ws*pts[l]*pts[l]*tp)/(1.0+ws*ws*pts[l]*pts[l]));
}
......@@ -86,12 +79,18 @@ void model_visco(float ** rho, float ** pi, float ** u, float ** taus, floa
/* two layer case */
if (y<=h){
Vp=vp1; Vs=vs1; Rhov=rho1; }
Vp=vp1; Vs=vs1; Rhov=rho1; tp=tp1; ts=ts1;}
else{
Vp=vp2; Vs=vs2; Rhov=rho2; }
Vp=vp2; Vs=vs2; Rhov=rho2; tp=tp2;ts=ts2;}
sumu=0.0;
sumpi=0.0;
for (l=1;l<=L;l++){
sumu=sumu+((ws*ws*pts[l]*pts[l]*ts)/(1.0+ws*ws*pts[l]*pts[l]));
sumpi=sumpi+((ws*ws*pts[l]*pts[l]*tp)/(1.0+ws*ws*pts[l]*pts[l]));
}
muv=Vs*Vs*Rhov/(1.0+sumu);
piv=Vp*Vp*Rhov/(1.0+sumpi);
......
......@@ -27,7 +27,7 @@
void matcopy_elastic(float ** rho, float ** pi, float ** u){
extern int MYID, NX, NY, INDEX[5];
extern int MYID, NX, NY, INDEX[5],POS[3], BOUNDARY, NPROCY, NPROCX;
extern const int TAG1,TAG2,TAG5,TAG6;
extern FILE *FP;
......@@ -40,8 +40,8 @@ void matcopy_elastic(float ** rho, float ** pi, float ** u){
bufferlef_to_rig = matrix(0,NY+1,1,3);
bufferrig_to_lef = matrix(0,NY+1,1,3);
buffertop_to_bot = matrix(0,NX+1,1,3);
bufferbot_to_top = matrix(0,NX+1,1,3);
buffertop_to_bot = matrix(1,NX,1,3);
bufferbot_to_top = matrix(1,NX,1,3);
fprintf(FP,"\n\n **Message from matcopy_elastic (written by PE %d):",MYID);
......@@ -51,8 +51,8 @@ void matcopy_elastic(float ** rho, float ** pi, float ** u){
/* if (POS[2]!=0)*/ /* no boundary exchange at top of global grid */
for (i=0;i<=NX+1;i++){
if (POS[2]!=0) /* no boundary exchange at top of global grid */
for (i=1;i<=NX;i++){
/* storage of top of local volume into buffer */
buffertop_to_bot[i][1] = rho[1][i];
buffertop_to_bot[i][2] = pi[1][i];
......@@ -60,8 +60,8 @@ void matcopy_elastic(float ** rho, float ** pi, float ** u){
}
/* if (POS[2]!=NPROCY-1)*/ /* no boundary exchange at bottom of global grid */
for (i=0;i<=NX+1;i++){
if (POS[2]!=NPROCY-1) /* no boundary exchange at bottom of global grid */
for (i=1;i<=NX;i++){
/* storage of bottom of local volume into buffer */
bufferbot_to_top[i][1] = rho[NY][i];
......@@ -80,15 +80,15 @@ void matcopy_elastic(float ** rho, float ** pi, float ** u){
MPI_Recv(&bufferbot_to_top[1][1],NX*3,MPI_FLOAT,INDEX[3],TAG6,MPI_COMM_WORLD,&status);
/* if (POS[2]!=NPROCY-1)*/ /* no boundary exchange at bottom of global grid */
for (i=0;i<=NX+1;i++){
if (POS[2]!=NPROCY-1) /* no boundary exchange at bottom of global grid */
for (i=1;i<=NX;i++){
rho[NY+1][i] = buffertop_to_bot[i][1];
pi[NY+1][i] = buffertop_to_bot[i][2];
u[NY+1][i] = buffertop_to_bot[i][3];
}
/* if (POS[2]!=0)*/ /* no boundary exchange at top of global grid */
for (i=0;i<=NX+1;i++){
if (POS[2]!=0) /* no boundary exchange at top of global grid */
for (i=1;i<=NX;i++){
rho[0][i] = bufferbot_to_top[i][1];
pi[0][i] = bufferbot_to_top[i][2];
u[0][i] = bufferbot_to_top[i][3];
......@@ -97,7 +97,7 @@ void matcopy_elastic(float ** rho, float ** pi, float ** u){
/* if (POS[1]!=0)*/ /* no boundary exchange at left edge of global grid */
if ((POS[1]!=0) || (BOUNDARY!=0)) /* no boundary exchange at left edge of global grid */
for (j=0;j<=NY+1;j++)
{
/* storage of left edge of local volume into buffer */
......@@ -107,7 +107,7 @@ void matcopy_elastic(float ** rho, float ** pi, float ** u){
}
/* if (POS[1]!=NPROCX-1)*/ /* no boundary exchange at right edge of global grid */
if ((POS[1]!=NPROCX-1) || (BOUNDARY!=0)) /* no boundary exchange at right edge of global grid */
for (j=0;j<=NY+1;j++){
/* storage of right edge of local volume into buffer */
bufferrig_to_lef[j][1] = rho[j][NX];
......@@ -125,14 +125,14 @@ void matcopy_elastic(float ** rho, float ** pi, float ** u){
MPI_Recv(&bufferrig_to_lef[0][1],(NY+2)*3,MPI_FLOAT,INDEX[1],TAG2,MPI_COMM_WORLD,&status);
/* if (POS[1]!=NPROCX-1)*/ /* no boundary exchange at right edge of global grid */
if ((POS[1]!=NPROCX-1) || (BOUNDARY!=0)) /* no boundary exchange at right edge of global grid */
for (j=0;j<=NY+1;j++){
rho[j][NX+1] = bufferlef_to_rig[j][1];
pi[j][NX+1] = bufferlef_to_rig[j][2];
u[j][NX+1] = bufferlef_to_rig[j][3];
}
/* if (POS[1]!=0)*/ /* no boundary exchange at left edge of global grid */
if ((POS[1]!=0) || (BOUNDARY!=0)) /* no boundary exchange at left edge of global grid */
for (j=0;j<=NY+1;j++){
rho[j][0] = bufferrig_to_lef[j][1];
pi[j][0] = bufferrig_to_lef[j][2];
......