Commit 2607899f authored by Tilman Steinweg's avatar Tilman Steinweg

imolemented structs

parent d000f3eb
......@@ -84,12 +84,7 @@
#include "fd.h"
void CPML_coeff(float * K_x, float * alpha_prime_x, float * a_x, float * b_x,
float * K_x_half, float * alpha_prime_x_half, float * a_x_half, float * b_x_half,
float * K_y, float * alpha_prime_y, float * a_y, float * b_y,
float * K_y_half, float * alpha_prime_y_half, float * a_y_half, float * b_y_half,
float * K_z, float * alpha_prime_z, float * a_z, float * b_z,
float * K_z_half, float * alpha_prime_z_half, float * a_z_half, float * b_z_half)
void CPML_coeff(st_pml_coeff *pml_coeff)
{
extern float DX, DY, DZ , VPPML, DT, FPML;
......@@ -126,7 +121,7 @@ void CPML_coeff(float * K_x, float * alpha_prime_x, float * a_x, float * b_x,
for (i=1;i<=FW+1;i++){
K_x[i] = 1.0;
pml_coeff->K_x[i] = 1.0;
/* define damping profile at the grid points */
position_in_PML = (FW-i+1)*DX; /*distance to boundary in meter */
......@@ -135,16 +130,16 @@ void CPML_coeff(float * K_x, float * alpha_prime_x, float * a_x, float * b_x,
d_x[i] = d0_x *(a*position_norm+b*pow(position_norm,npower)+c*pow(position_norm,(4)));
/* this taken from Gedney page 8.2 */
K_x[i] = 1.0 + (k_max_PML - 1.0) * pow(position_norm,npower);
alpha_prime_x[i] = alpha_max_PML * (1.0 - position_norm);
pml_coeff->K_x[i] = 1.0 + (k_max_PML - 1.0) * pow(position_norm,npower);
pml_coeff->alpha_prime_x[i] = alpha_max_PML * (1.0 - position_norm);
if(alpha_prime_x[i] < 0.0){ fprintf(FP,"ERROR:alpha_prime_x[i] < 0.0, i %d", i);}
if(pml_coeff->alpha_prime_x[i] < 0.0){ fprintf(FP,"ERROR:alpha_prime_x[i] < 0.0, i %d", i);}
b_x[i] = exp(- (d_x[i] / K_x[i] + alpha_prime_x[i]) * DT);
pml_coeff->b_x[i] = exp(- (d_x[i] / pml_coeff->K_x[i] + pml_coeff->alpha_prime_x[i]) * DT);
/* avoid division by zero outside the PML */
if(abs(d_x[i]) > 1.0e-6){ a_x[i] = d_x[i] * (b_x[i] - 1.0) / (K_x[i] * (d_x[i] + K_x[i] * alpha_prime_x[i]));}
else a_x[i]=0.0;
if(abs(d_x[i]) > 1.0e-6){ pml_coeff->a_x[i] = d_x[i] * (pml_coeff->b_x[i] - 1.0) / (pml_coeff->K_x[i] * (d_x[i] + pml_coeff->K_x[i] * pml_coeff->alpha_prime_x[i]));}
else pml_coeff->a_x[i]=0.0;
if(i<=FW){
......@@ -154,35 +149,35 @@ void CPML_coeff(float * K_x, float * alpha_prime_x, float * a_x, float * b_x,
i1=i;
K_x_half[i1] = 1.0;
pml_coeff->K_x_half[i1] = 1.0;
d_x_half[i1] = d0_x * (a*position_norm+b*pow(position_norm,npower)+c*pow(position_norm,(4)));
if(position_in_PML < 0.0) {fprintf(FP,"ERROR: Position in PML (x-boundary) smaller 0");}
/* this taken from Gedney page 8.2 */
K_x_half[i1] = 1.0 + (k_max_PML - 1.0) * pow(position_norm,npower);
alpha_prime_x_half[i1] = alpha_max_PML * (1.0 - position_norm);
pml_coeff->K_x_half[i1] = 1.0 + (k_max_PML - 1.0) * pow(position_norm,npower);
pml_coeff->alpha_prime_x_half[i1] = alpha_max_PML * (1.0 - position_norm);
/* just in case, for -5 at the end */
if(alpha_prime_x_half[i1] < 0.0) {fprintf(FP,"ERROR:alpha_prime_x_half[i] < 0.0, i %d", i);}
if(pml_coeff->alpha_prime_x_half[i1] < 0.0) {fprintf(FP,"ERROR:alpha_prime_x_half[i] < 0.0, i %d", i);}
b_x_half[i1] = exp(- (d_x_half[i1] / K_x_half[i1] + alpha_prime_x_half[i1]) * DT);
pml_coeff->b_x_half[i1] = exp(- (d_x_half[i1] / pml_coeff->K_x_half[i1] + pml_coeff->alpha_prime_x_half[i1]) * DT);
if(abs(d_x_half[i1]) > 1.0e-6){ a_x_half[i1] = d_x_half[i1] * (b_x_half[i1] - 1.0) / (K_x_half[i1] * (d_x_half[i1] + K_x_half[i1] * alpha_prime_x_half[i1]));}
if(abs(d_x_half[i1]) > 1.0e-6){ pml_coeff->a_x_half[i1] = d_x_half[i1] * (pml_coeff->b_x_half[i1] - 1.0) / (pml_coeff->K_x_half[i1] * (d_x_half[i1] + pml_coeff->K_x_half[i1] * pml_coeff->alpha_prime_x_half[i1]));}
/* right boundary --> mirroring left boundary*/
l = 2* FW -i+1;
if(i>1){
K_x[l+1]=K_x[i];
b_x[l+1] = b_x[i];
if(abs(d_x[i]) > 1.0e-6){ a_x[l+1] = a_x[i]; }
pml_coeff->K_x[l+1]=pml_coeff->K_x[i];
pml_coeff->b_x[l+1] = pml_coeff->b_x[i];
if(abs(d_x[i]) > 1.0e-6){ pml_coeff->a_x[l+1] = pml_coeff->a_x[i]; }
}
K_x_half[l]=K_x_half[i];
b_x_half[l] = b_x_half[i]; /*half a grid point in +x)*/
if(abs(d_x[i]) > 1.0e-6){ a_x_half[l] = a_x_half[i]; }
pml_coeff->K_x_half[l]=pml_coeff->K_x_half[i];
pml_coeff->b_x_half[l] = pml_coeff->b_x_half[i]; /*half a grid point in +x)*/
if(abs(d_x[i]) > 1.0e-6){ pml_coeff->a_x_half[l] = pml_coeff->a_x_half[i]; }
}
}
......@@ -194,7 +189,7 @@ void CPML_coeff(float * K_x, float * alpha_prime_x, float * a_x, float * b_x,
for (i=1;i<=FW+1;i++){
K_y[i] = 1.0;
pml_coeff->K_y[i] = 1.0;
/* define damping profile at the grid points */
position_in_PML = (FW-i+1)*DY; /*distance to boundary in meter */
......@@ -203,17 +198,17 @@ void CPML_coeff(float * K_x, float * alpha_prime_x, float * a_x, float * b_x,
d_y[i] = d0_y * (a*position_norm+b*pow(position_norm,npower)+c*pow(position_norm,(4)));
/* this taken from Gedney page 8.2 */
K_y[i] = 1.0 + (k_max_PML - 1.0) * pow(position_norm,npower);
alpha_prime_y[i] = alpha_max_PML * (1.0 - position_norm);
pml_coeff->K_y[i] = 1.0 + (k_max_PML - 1.0) * pow(position_norm,npower);
pml_coeff->alpha_prime_y[i] = alpha_max_PML * (1.0 - position_norm);
/* just in case, for -5 at the end */
if(alpha_prime_y[i] < 0.0){ fprintf(FP,"ERROR:alpha_prime_y[i] < 0.0, i %d", i);}
if(pml_coeff->alpha_prime_y[i] < 0.0){ fprintf(FP,"ERROR:alpha_prime_y[i] < 0.0, i %d", i);}
b_y[i] = exp(- (d_y[i] / K_y[i] + alpha_prime_y[i]) * DT);
pml_coeff->b_y[i] = exp(- (d_y[i] / pml_coeff->K_y[i] + pml_coeff->alpha_prime_y[i]) * DT);
/* avoid division by zero outside the PML */
if(abs(d_y[i]) > 1.0e-6){ a_y[i] = d_y[i] * (b_y[i] - 1.0) / (K_y[i] * (d_y[i] + K_y[i] * alpha_prime_y[i]));}
else a_x[i]=0.0;
if(abs(d_y[i]) > 1.0e-6){ pml_coeff->a_y[i] = d_y[i] * (pml_coeff->b_y[i] - 1.0) / (pml_coeff->K_y[i] * (d_y[i] + pml_coeff->K_y[i] *pml_coeff->alpha_prime_y[i]));}
else pml_coeff->a_x[i]=0.0;
if(i<=FW){
......@@ -222,33 +217,33 @@ void CPML_coeff(float * K_x, float * alpha_prime_x, float * a_x, float * b_x,
position_norm = position_in_PML / (FW*DY);
i1=i;
K_y_half[i1] = 1.0;
pml_coeff->K_y_half[i1] = 1.0;
d_y_half[i1] = d0_y * (a*position_norm+b*pow(position_norm,npower)+c*pow(position_norm,(4)));
if(position_in_PML < 0.0) {fprintf(FP,"ERROR: Position in PML (y-boundary) <0");}
/* this taken from Gedney page 8.2 */
K_y_half[i1] = 1.0 + (k_max_PML - 1.0) * pow(position_norm,npower);
alpha_prime_y_half[i1] = alpha_max_PML * (1.0 - position_norm);
pml_coeff->K_y_half[i1] = 1.0 + (k_max_PML - 1.0) * pow(position_norm,npower);
pml_coeff->alpha_prime_y_half[i1] = alpha_max_PML * (1.0 - position_norm);
if(alpha_prime_y_half[i1] < 0.0) {fprintf(FP,"ERROR:alpha_prime_y_half[i] < 0.0, i %d", i);}
b_y_half[i1] = exp(- (d_y_half[i1] / K_y_half[i1] + alpha_prime_y_half[i1]) * DT);
if(pml_coeff->alpha_prime_y_half[i1] < 0.0) {fprintf(FP,"ERROR:alpha_prime_y_half[i] < 0.0, i %d", i);}
pml_coeff->b_y_half[i1] = exp(- (d_y_half[i1] / pml_coeff->K_y_half[i1] + pml_coeff->alpha_prime_y_half[i1]) * DT);
if(abs(d_y_half[i1]) > 1.0e-6){ a_y_half[i1] = d_y_half[i1] * (b_y_half[i1] - 1.0) / (K_y_half[i1] * (d_y_half[i1] + K_y_half[i1] * alpha_prime_y_half[i1]));}
if(abs(d_y_half[i1]) > 1.0e-6){ pml_coeff->a_y_half[i1] = d_y_half[i1] * (pml_coeff->b_y_half[i1] - 1.0) / (pml_coeff->K_y_half[i1] * (d_y_half[i1] + pml_coeff->K_y_half[i1] * pml_coeff->alpha_prime_y_half[i1]));}
/* right boundary --> mirroring left boundary*/
l = 2* FW -i+1;
if(i>1){
K_y[l+1] = K_y[i];
b_y[l+1] = b_y[i];
if(abs(d_y[i]) > 1.0e-6){ a_y[l+1] = a_y[i]; }
pml_coeff->K_y[l+1] = pml_coeff->K_y[i];
pml_coeff->b_y[l+1] = pml_coeff->b_y[i];
if(abs(d_y[i]) > 1.0e-6){ pml_coeff->a_y[l+1] = pml_coeff->a_y[i]; }
}
K_y_half[l]=K_y_half[i];
b_y_half[l] = b_y_half[i]; /*half a grid point in +x)*/
if(abs(d_y[i]) > 1.0e-6){ a_y_half[l] = a_y_half[i]; }
pml_coeff->K_y_half[l]=pml_coeff->K_y_half[i];
pml_coeff->b_y_half[l] = pml_coeff->b_y_half[i]; /*half a grid point in +x)*/
if(abs(d_y[i]) > 1.0e-6){ pml_coeff->a_y_half[l] = pml_coeff->a_y_half[i]; }
}
}
......@@ -258,7 +253,7 @@ void CPML_coeff(float * K_x, float * alpha_prime_x, float * a_x, float * b_x,
for (i=1;i<=FW+1;i++){
K_z[i] = 1.0;
pml_coeff->K_z[i] = 1.0;
/* define damping profile at the grid points */
position_in_PML = (FW-i+1)*DZ; /*distance to boundary in meter */
......@@ -267,16 +262,16 @@ void CPML_coeff(float * K_x, float * alpha_prime_x, float * a_x, float * b_x,
d_z[i] = d0_z * (a*position_norm+b*pow(position_norm,npower)+c*pow(position_norm,(4)));
/* this taken from Gedney page 8.2 */
K_z[i] = 1.0 + (k_max_PML - 1.0) * pow(position_norm,npower);
alpha_prime_z[i] = alpha_max_PML * (1.0 - position_norm);
pml_coeff->K_z[i] = 1.0 + (k_max_PML - 1.0) * pow(position_norm,npower);
pml_coeff->alpha_prime_z[i] = alpha_max_PML * (1.0 - position_norm);
/* just in case, for -5 at the end */
if(alpha_prime_z[i] < 0.0){ fprintf(FP,"ERROR:alpha_prime_z[i] < 0.0, i %d", i);}
if(pml_coeff->alpha_prime_z[i] < 0.0){ fprintf(FP,"ERROR:alpha_prime_z[i] < 0.0, i %d", i);}
b_z[i] = exp(- (d_z[i] / K_z[i] + alpha_prime_z[i]) * DT);
pml_coeff->b_z[i] = exp(- (d_z[i] / pml_coeff->K_z[i] + pml_coeff->alpha_prime_z[i]) * DT);
/* avoid division by zero outside the PML */
if(abs(d_z[i]) > 1.0e-6){ a_z[i] = d_z[i] * (b_z[i] - 1.0) / (K_z[i] * (d_z[i] + K_z[i] * alpha_prime_z[i]));}
if(abs(d_z[i]) > 1.0e-6){ pml_coeff->a_z[i] = d_z[i] * (pml_coeff->b_z[i] - 1.0) / (pml_coeff->K_z[i] * (d_z[i] + pml_coeff->K_z[i] *pml_coeff->alpha_prime_z[i]));}
if(i<=FW){
/* define damping profile at half the grid points (half a grid point in -x)*/
......@@ -285,34 +280,34 @@ void CPML_coeff(float * K_x, float * alpha_prime_x, float * a_x, float * b_x,
i1=i;
K_z_half[i1] = 1.0;
pml_coeff->K_z_half[i1] = 1.0;
d_z_half[i1] = d0_z *(a*position_norm+b* pow(position_norm,npower)+c*pow(position_norm,(4)));
if(position_in_PML < 0.0) {fprintf(FP,"ERROR: Position in PML (y-boundary) <0");}
/* this taken from Gedney page 8.2 */
K_z_half[i1] = 1.0 + (k_max_PML - 1.0) * pow(position_norm,npower);
alpha_prime_z_half[i1] = alpha_max_PML * (1.0 - position_norm);
pml_coeff->K_z_half[i1] = 1.0 + (k_max_PML - 1.0) * pow(position_norm,npower);
pml_coeff->alpha_prime_z_half[i1] = alpha_max_PML * (1.0 - position_norm);
if(alpha_prime_z_half[i1] < 0.0) {fprintf(FP,"ERROR:alpha_prime_z_half[i] < 0.0, i %d", i);}
if(pml_coeff->alpha_prime_z_half[i1] < 0.0) {fprintf(FP,"ERROR:alpha_prime_z_half[i] < 0.0, i %d", i);}
b_z_half[i1] = exp(- (d_z_half[i1] / K_z_half[i1] + alpha_prime_z_half[i1]) * DT);
pml_coeff->b_z_half[i1] = exp(- (d_z_half[i1] / pml_coeff->K_z_half[i1] + pml_coeff->alpha_prime_z_half[i1]) * DT);
if(abs(d_z_half[i1]) > 1.0e-6){ a_z_half[i1] = d_z_half[i1] * (b_z_half[i1] - 1.0) / (K_z_half[i1] * (d_z_half[i1] + K_z_half[i1] * alpha_prime_z_half[i1]));}
if(abs(d_z_half[i1]) > 1.0e-6){ pml_coeff->a_z_half[i1] = d_z_half[i1] * (pml_coeff->b_z_half[i1] - 1.0) / (pml_coeff->K_z_half[i1] * (d_z_half[i1] + pml_coeff->K_z_half[i1] * pml_coeff->alpha_prime_z_half[i1]));}
/* right boundary --> mirroring left boundary*/
l = 2* FW -i+1;
if(i>1){
K_z[l+1] = K_z[i];
b_z[l+1] = b_z[i];
if(abs(d_z[i]) > 1.0e-6){ a_z[l+1] = a_z[i]; }
pml_coeff->K_z[l+1] = pml_coeff->K_z[i];
pml_coeff->b_z[l+1] = pml_coeff->b_z[i];
if(abs(d_z[i]) > 1.0e-6){ pml_coeff->a_z[l+1] = pml_coeff->a_z[i]; }
}
K_z_half[l]=K_z_half[i];
b_z_half[l] = b_z_half[i]; /*half a grid point in +x)*/
if(abs(d_z[i]) > 1.0e-6){ a_z_half[l] = a_z_half[i]; }
pml_coeff->K_z_half[l]=pml_coeff->K_z_half[i];
pml_coeff->b_z_half[l] = pml_coeff->b_z_half[i]; /*half a grid point in +x)*/
if(abs(d_z[i]) > 1.0e-6){ pml_coeff->a_z_half[l] = pml_coeff->a_z_half[i]; }
}
}
free_vector(d_x,1,2*FW);
......
......@@ -25,10 +25,7 @@
#include "fd.h"
void av_mat(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup,
float *** uipjp, float *** ujpkp, float *** uipkp, float *** tausipjp,
float *** tausjpkp, float *** tausipkp, float *** rjp, float *** rkp, float *** rip ){
void av_mat(st_model *mod,st_model_av *mod_av){
extern int NX, NY, NZ, MYID,L,VERBOSE;
......@@ -48,20 +45,20 @@ float *** tausjpkp, float *** tausipkp, float *** rjp, float *** rkp, float
/* harmonic averaging of shear modulus */
uipjp[j][i][k]=4.0/((1.0/u[j][i][k])+(1.0/u[j][i+1][k])+(1.0/u[j+1][i+1][k])+(1.0/u[j+1][i][k]));
ujpkp[j][i][k]=4.0/((1.0/u[j][i][k])+(1.0/u[j][i][k+1])+(1.0/u[j+1][i][k+1])+(1.0/u[j+1][i][k]));
uipkp[j][i][k]=4.0/((1.0/u[j][i][k])+(1.0/u[j][i][k+1])+(1.0/u[j][i+1][k+1])+(1.0/u[j][i+1][k]));
mod_av->uipjp[j][i][k]=4.0/((1.0/mod->u[j][i][k])+(1.0/mod->u[j][i+1][k])+(1.0/mod->u[j+1][i+1][k])+(1.0/mod->u[j+1][i][k]));
mod_av->ujpkp[j][i][k]=4.0/((1.0/mod->u[j][i][k])+(1.0/mod->u[j][i][k+1])+(1.0/mod->u[j+1][i][k+1])+(1.0/mod->u[j+1][i][k]));
mod_av->uipkp[j][i][k]=4.0/((1.0/mod->u[j][i][k])+(1.0/mod->u[j][i][k+1])+(1.0/mod->u[j][i+1][k+1])+(1.0/mod->u[j][i+1][k]));
/* arithmetic averaging of TAU for S-waves and density */
if (L){
tausipjp[j][i][k]=0.25*(taus[j][i][k]+taus[j][i+1][k]+taus[j+1][i+1][k]+taus[j+1][i][k]);
tausjpkp[j][i][k]=0.25*(taus[j][i][k]+taus[j+1][i][k]+taus[j+1][i][k+1]+taus[j][i][k+1]);
tausipkp[j][i][k]=0.25*(taus[j][i][k]+taus[j][i+1][k]+taus[j][i+1][k+1]+taus[j][i][k+1]);
mod_av->tausipjp[j][i][k]=0.25*(mod->taus[j][i][k]+mod->taus[j][i+1][k]+mod->taus[j+1][i+1][k]+mod->taus[j+1][i][k]);
mod_av->tausjpkp[j][i][k]=0.25*(mod->taus[j][i][k]+mod->taus[j+1][i][k]+mod->taus[j+1][i][k+1]+mod->taus[j][i][k+1]);
mod_av->tausipkp[j][i][k]=0.25*(mod->taus[j][i][k]+mod->taus[j][i+1][k]+mod->taus[j][i+1][k+1]+mod->taus[j][i][k+1]);
}
rjp[j][i][k]=0.5*(rho[j][i][k]+rho[j+1][i][k]);
rkp[j][i][k]=0.5*(rho[j][i][k]+rho[j][i][k+1]);
rip[j][i][k]=0.5*(rho[j][i][k]+rho[j][i+1][k]);
mod_av->rjp[j][i][k]=0.5*(mod->rho[j][i][k]+mod->rho[j+1][i][k]);
mod_av->rkp[j][i][k]=0.5*(mod->rho[j][i][k]+mod->rho[j][i][k+1]);
mod_av->rip[j][i][k]=0.5*(mod->rho[j][i][k]+mod->rho[j][i+1][k]);
}
......
......@@ -37,8 +37,7 @@ void err2(char errformat[],char errfilename[]){
err(outtxt);
}
void checkfd(FILE *fp, float *** prho, float *** ppi, float *** pu,
float *** ptaus, float *** ptaup, float *peta, float **srcpos, int nsrc, int **recpos, int ntr){
void checkfd(FILE *fp, st_model *mod, float **srcpos, int nsrc, int **recpos, int ntr){
extern float DX, DY, DZ, DT, TS, TIME, TSNAP2;
extern float XREC1, XREC2, YREC1, YREC2, ZREC1, ZREC2;
......@@ -80,12 +79,12 @@ float *** ptaus, float *** ptaup, float *peta, float **srcpos, int nsrc, int **r
for (k=1+nfw;k<=(nz-nfw);k++){
for (i=1+nfw;i<=(nx-nfw);i++){
for (j=ny1;j<=(ny-nfw);j++){
if (prho[j][i][k]==0.0)
printf("prho=0 detected from MYID=%d at i=%d,j=%d,k=%d\n",MYID,i,j,k);
if (mod->rho[j][i][k]==0.0)
printf("rho=0 detected from MYID=%d at i=%d,j=%d,k=%d\n",MYID,i,j,k);
else{
if(L){
c=sqrt(pu[j][i][k]*(1.0+L*ptaus[j][i][k])/prho[j][i][k]);}
else c=sqrt(pu[j][i][k]/prho[j][i][k]);
c=sqrt(mod->u[j][i][k]*(1.0+L*mod->taus[j][i][k])/mod->rho[j][i][k]);}
else c=sqrt(mod->u[j][i][k]/mod->rho[j][i][k]);
}
if (cmax_s<c) cmax_s=c;
......@@ -93,13 +92,13 @@ float *** ptaus, float *** ptaup, float *peta, float **srcpos, int nsrc, int **r
frequency of the source */
sum=0.0;
for (l=1;l<=L;l++){
ts=DT/peta[l];
sum=sum+((w*w*ptaus[j][i][k]*ts*ts)/(1.0+w*w*ts*ts));
ts=DT/mod->eta[l];
sum=sum+((w*w*mod->taus[j][i][k]*ts*ts)/(1.0+w*w*ts*ts));
}
if (prho[j][i][k]==0.0)
printf("prho=0 detected from MYID=%d at i=%d,j=%d,k=%d\n",MYID,i,j,k);
if (mod->rho[j][i][k]==0.0)
printf("rho=0 detected from MYID=%d at i=%d,j=%d,k=%d\n",MYID,i,j,k);
else
c=sqrt(pu[j][i][k]*(1.0+sum)/prho[j][i][k]);
c=sqrt(mod->u[j][i][k]*(1.0+sum)/mod->rho[j][i][k]);
if ((c>cwater)&&(cmin_s>c)) cmin_s=c;
}
}
......@@ -112,24 +111,24 @@ float *** ptaus, float *** ptaup, float *peta, float **srcpos, int nsrc, int **r
for (k=1+nfw;k<=(nz-nfw);k++){
for (i=1+nfw;i<=(nx-nfw);i++){
for (j=ny1;j<=(ny-nfw);j++){
if (prho[j][i][k]==0.0)
printf("prho=0 detected from MYID=%d at i=%d,j=%d,k=%d\n",MYID,i,j,k);
if (mod->rho[j][i][k]==0.0)
printf("rho=0 detected from MYID=%d at i=%d,j=%d,k=%d\n",MYID,i,j,k);
else{
if(L) c=sqrt(ppi[j][i][k]*(1.0+L*ptaup[j][i][k])/prho[j][i][k]);
else c=sqrt(ppi[j][i][k]/prho[j][i][k]);
if(L) c=sqrt(mod->pi[j][i][k]*(1.0+L*mod->taup[j][i][k])/mod->rho[j][i][k]);
else c=sqrt(mod->pi[j][i][k]/mod->rho[j][i][k]);
}
if (cmax_p<c) cmax_p=c;
/* find minimum model phase velocity of P-waves at center
frequency of the source */
sum=0.0;
for (l=1;l<=L;l++){
ts=DT/peta[l];
sum=sum+((w*w*ptaup[j][i][k]*ts*ts)/(1.0+w*w*ts*ts));
ts=DT/mod->eta[l];
sum=sum+((w*w*mod->taup[j][i][k]*ts*ts)/(1.0+w*w*ts*ts));
}
if (prho[j][i][k]==0.0)
printf("prho=0 detected from MYID=%d at i=%d,j=%d,k=%d\n",MYID,i,j,k);
if (mod->rho[j][i][k]==0.0)
printf("rho=0 detected from MYID=%d at i=%d,j=%d,k=%d\n",MYID,i,j,k);
else
c=sqrt(ppi[j][i][k]*(1.0+sum)/prho[j][i][k]);
c=sqrt(mod->pi[j][i][k]*(1.0+sum)/mod->rho[j][i][k]);
if (cmin_p>c) cmin_p=c;
}
}
......
......@@ -26,10 +26,7 @@
#include "fd.h"
void comm_ini(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){
void comm_ini(st_buffer *buffer, MPI_Request *req_send, MPI_Request *req_rec){
extern int NX, NY, NZ, INDEX[7], FDORDER;
......@@ -50,22 +47,24 @@ MPI_Request *req_send, MPI_Request *req_rec){
nf2=nf1-1;
MPI_Bsend_init(&bufferlef_to_rig[1][1][1],NY*NZ*nf1,MPI_FLOAT,INDEX[1],TAG1,MPI_COMM_WORLD,&req_send[0]);
MPI_Bsend_init(&bufferrig_to_lef[1][1][1],NY*NZ*nf2,MPI_FLOAT,INDEX[2],TAG2,MPI_COMM_WORLD,&req_send[1]);
MPI_Bsend_init(&bufferfro_to_bac[1][1][1],NX*NY*nf1,MPI_FLOAT,INDEX[5],TAG3,MPI_COMM_WORLD,&req_send[2]);
MPI_Bsend_init(&bufferbac_to_fro[1][1][1],NX*NY*nf2,MPI_FLOAT,INDEX[6],TAG4,MPI_COMM_WORLD,&req_send[3]);
MPI_Bsend_init(&buffertop_to_bot[1][1][1],NX*NZ*nf1,MPI_FLOAT,INDEX[3],TAG5,MPI_COMM_WORLD,&req_send[4]);
MPI_Bsend_init(&bufferbot_to_top[1][1][1],NX*NZ*nf2,MPI_FLOAT,INDEX[4],TAG6,MPI_COMM_WORLD,&req_send[5]);
MPI_Bsend_init(&buffer->lef_to_rig[1][1][1],NY*NZ*nf1,MPI_FLOAT,INDEX[1],TAG1,MPI_COMM_WORLD,&req_send[0]);
MPI_Bsend_init(&buffer->rig_to_lef[1][1][1],NY*NZ*nf2,MPI_FLOAT,INDEX[2],TAG2,MPI_COMM_WORLD,&req_send[1]);
MPI_Bsend_init(&buffer->fro_to_bac[1][1][1],NX*NY*nf1,MPI_FLOAT,INDEX[5],TAG3,MPI_COMM_WORLD,&req_send[2]);
MPI_Bsend_init(&buffer->bac_to_fro[1][1][1],NX*NY*nf2,MPI_FLOAT,INDEX[6],TAG4,MPI_COMM_WORLD,&req_send[3]);
MPI_Bsend_init(&buffer->top_to_bot[1][1][1],NX*NZ*nf1,MPI_FLOAT,INDEX[3],TAG5,MPI_COMM_WORLD,&req_send[4]);
MPI_Bsend_init(&buffer->bot_to_top[1][1][1],NX*NZ*nf2,MPI_FLOAT,INDEX[4],TAG6,MPI_COMM_WORLD,&req_send[5]);
/* initialising of receive of buffer arrays. Same arrays for send and receive
are used. Thus, before starting receive, it is necessary to check if
Bsend has copied data into local buffers, i.e. has completed.
*/
MPI_Recv_init(&bufferlef_to_rig[1][1][1],NY*NZ*nf1,MPI_FLOAT,INDEX[2],TAG1,MPI_COMM_WORLD,&req_rec[0]);
MPI_Recv_init(&bufferrig_to_lef[1][1][1],NY*NZ*nf2,MPI_FLOAT,INDEX[1],TAG2,MPI_COMM_WORLD,&req_rec[1]);
MPI_Recv_init(&bufferfro_to_bac[1][1][1],NX*NY*nf1,MPI_FLOAT,INDEX[6],TAG3,MPI_COMM_WORLD,&req_rec[2]);
MPI_Recv_init(&bufferbac_to_fro[1][1][1],NX*NY*nf2,MPI_FLOAT,INDEX[5],TAG4,MPI_COMM_WORLD,&req_rec[3]);
MPI_Recv_init(&buffertop_to_bot[1][1][1],NX*NZ*nf1,MPI_FLOAT,INDEX[4],TAG5,MPI_COMM_WORLD,&req_rec[4]);
MPI_Recv_init(&bufferbot_to_top[1][1][1],NX*NZ*nf2,MPI_FLOAT,INDEX[3],TAG6,MPI_COMM_WORLD,&req_rec[5]);
MPI_Recv_init(&buffer->lef_to_rig[1][1][1],NY*NZ*nf1,MPI_FLOAT,INDEX[2],TAG1,MPI_COMM_WORLD,&req_rec[0]);
MPI_Recv_init(&buffer->rig_to_lef[1][1][1],NY*NZ*nf2,MPI_FLOAT,INDEX[1],TAG2,MPI_COMM_WORLD,&req_rec[1]);
MPI_Recv_init(&buffer->fro_to_bac[1][1][1],NX*NY*nf1,MPI_FLOAT,INDEX[6],TAG3,MPI_COMM_WORLD,&req_rec[2]);
MPI_Recv_init(&buffer->bac_to_fro[1][1][1],NX*NY*nf2,MPI_FLOAT,INDEX[5],TAG4,MPI_COMM_WORLD,&req_rec[3]);
MPI_Recv_init(&buffer->top_to_bot[1][1][1],NX*NZ*nf1,MPI_FLOAT,INDEX[4],TAG5,MPI_COMM_WORLD,&req_rec[4]);
MPI_Recv_init(&buffer->bot_to_top[1][1][1],NX*NZ*nf2,MPI_FLOAT,INDEX[3],TAG6,MPI_COMM_WORLD,&req_rec[5]);
}
......@@ -7,10 +7,7 @@
#include "fd.h"
void comm_ini_s(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){
void comm_ini_s(st_buffer *buffer, MPI_Request *req_send, MPI_Request *req_rec){
extern int NX, NY, NZ, INDEX[7], FDORDER;
......@@ -31,22 +28,22 @@ MPI_Request *req_send, MPI_Request *req_rec){
nf2=nf1-1;
MPI_Bsend_init(&bufferlef_to_rig[1][1][1],NY*NZ*nf2,MPI_FLOAT,INDEX[1],TAG1,MPI_COMM_WORLD,&req_send[0]);
MPI_Bsend_init(&bufferrig_to_lef[1][1][1],NY*NZ*nf1,MPI_FLOAT,INDEX[2],TAG2,MPI_COMM_WORLD,&req_send[1]);
MPI_Bsend_init(&bufferfro_to_bac[1][1][1],NX*NY*nf2,MPI_FLOAT,INDEX[5],TAG3,MPI_COMM_WORLD,&req_send[2]);
MPI_Bsend_init(&bufferbac_to_fro[1][1][1],NX*NY*nf1,MPI_FLOAT,INDEX[6],TAG4,MPI_COMM_WORLD,&req_send[3]);
MPI_Bsend_init(&buffertop_to_bot[1][1][1],NX*NZ*nf2,MPI_FLOAT,INDEX[3],TAG5,MPI_COMM_WORLD,&req_send[4]);
MPI_Bsend_init(&bufferbot_to_top[1][1][1],NX*NZ*nf1,MPI_FLOAT,INDEX[4],TAG6,MPI_COMM_WORLD,&req_send[5]);
MPI_Bsend_init(&buffer->lef_to_rig[1][1][1],NY*NZ*nf2,MPI_FLOAT,INDEX[1],TAG1,MPI_COMM_WORLD,&req_send[0]);
MPI_Bsend_init(&buffer->rig_to_lef[1][1][1],NY*NZ*nf1,MPI_FLOAT,INDEX[2],TAG2,MPI_COMM_WORLD,&req_send[1]);
MPI_Bsend_init(&buffer->fro_to_bac[1][1][1],NX*NY*nf2,MPI_FLOAT,INDEX[5],TAG3,MPI_COMM_WORLD,&req_send[2]);
MPI_Bsend_init(&buffer->bac_to_fro[1][1][1],NX*NY*nf1,MPI_FLOAT,INDEX[6],TAG4,MPI_COMM_WORLD,&req_send[3]);
MPI_Bsend_init(&buffer->top_to_bot[1][1][1],NX*NZ*nf2,MPI_FLOAT,INDEX[3],TAG5,MPI_COMM_WORLD,&req_send[4]);
MPI_Bsend_init(&buffer->bot_to_top[1][1][1],NX*NZ*nf1,MPI_FLOAT,INDEX[4],TAG6,MPI_COMM_WORLD,&req_send[5]);
/* initialising of receive of buffer arrays. Same arrays for send and receive
are used. Thus, before starting receive, it is necessary to check if
Bsend has copied data into local buffers, i.e. has completed.
*/
MPI_Recv_init(&bufferlef_to_rig[1][1][1],NY*NZ*nf2,MPI_FLOAT,INDEX[2],TAG1,MPI_COMM_WORLD,&req_rec[0]);
MPI_Recv_init(&bufferrig_to_lef[1][1][1],NY*NZ*nf1,MPI_FLOAT,INDEX[1],TAG2,MPI_COMM_WORLD,&req_rec[1]);
MPI_Recv_init(&bufferfro_to_bac[1][1][1],NX*NY*nf2,MPI_FLOAT,INDEX[6],TAG3,MPI_COMM_WORLD,&req_rec[2]);
MPI_Recv_init(&bufferbac_to_fro[1][1][1],NX*NY*nf1,MPI_FLOAT,INDEX[5],TAG4,MPI_COMM_WORLD,&req_rec[3]);
MPI_Recv_init(&buffertop_to_bot[1][1][1],NX*NZ*nf2,MPI_FLOAT,INDEX[4],TAG5,MPI_COMM_WORLD,&req_rec[4]);
MPI_Recv_init(&bufferbot_to_top[1][1][1],NX*NZ*nf1,MPI_FLOAT,INDEX[3],TAG6,MPI_COMM_WORLD,&req_rec[5]);
MPI_Recv_init(&buffer->lef_to_rig[1][1][1],NY*NZ*nf2,MPI_FLOAT,INDEX[2],TAG1,MPI_COMM_WORLD,&req_rec[0]);
MPI_Recv_init(&buffer->rig_to_lef[1][1][1],NY*NZ*nf1,MPI_FLOAT,INDEX[1],TAG2,MPI_COMM_WORLD,&req_rec[1]);
MPI_Recv_init(&buffer->fro_to_bac[1][1][1],NX*NY*nf2,MPI_FLOAT,INDEX[6],TAG3,MPI_COMM_WORLD,&req_rec[2]);
MPI_Recv_init(&buffer->bac_to_fro[1][1][1],NX*NY*nf1,MPI_FLOAT,INDEX[5],TAG4,MPI_COMM_WORLD,&req_rec[3]);
MPI_Recv_init(&buffer->top_to_bot[1][1][1],NX*NZ*nf2,MPI_FLOAT,INDEX[4],TAG5,MPI_COMM_WORLD,&req_rec[4]);
MPI_Recv_init(&buffer->bot_to_top[1][1][1],NX*NZ*nf1,MPI_FLOAT,INDEX[3],TAG6,MPI_COMM_WORLD,&req_rec[5]);
}
......@@ -24,7 +24,7 @@
#include "fd.h"
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){
void conjugate(int nx,int ny,int nz, st_gradient *grad, st_gradient *grad_prior1,st_gradient *grad_prior2, float *beta, int iteration,int cdf){
extern FILE *FP;
......@@ -43,12 +43,12 @@ void conjugate(int nx,int ny,int nz, float ***grad1, float ***grad2,float ***gra
for (j=1;j<=ny;j++){
for (i=1;i<=nx;i++){
for (k=1;k<=nz;k++){
gradprior1[j][i][k]=grad1[j][i][k];
gradprior2[j][i][k]=grad2[j][i][k];
gradprior3[j][i][k]=grad3[j][i][k];
gradprior4[j][i][k]=grad1[j][i][k];
gradprior5[j][i][k]=grad2[j][i][k];
gradprior6[j][i][k]=grad3[j][i][k];
grad_prior1->vp[j][i][k]=grad->vp[j][i][k];
grad_prior1->vs[j][i][k]=grad->vs[j][i][k];
grad_prior1->rho[j][i][k]=grad->rho[j][i][k];
grad_prior2->vp[j][i][k]=grad->vp[j][i][k];
grad_prior2->vs[j][i][k]=grad->vs[j][i][k];
grad_prior2->rho[j][i][k]=grad->rho[j][i][k];
}
}
}
......@@ -59,17 +59,17 @@ void conjugate(int nx,int ny,int nz, float ***grad1, float ***grad2,float ***gra
for (j=1;j<=ny;j++){
for (i=1;i<=nx;i++){
for (k=1;k<=nz;k++){
betadum[0]+=grad1[j][i][k]*(grad1[j][i][k]-gradprior1[j][i][k]);
betadum[1]+=gradprior1[j][i][k]*gradprior1[j][i][k];
betadum[2]+=grad2[j][i][k]*(grad2[j][i][k]-gradprior2[j][i][k]);
betadum[3]+=gradprior2[j][i][k]*gradprior2[j][i][k];
betadum[4]+=grad3[j][i][k]*(grad3[j][i][k]-gradprior3[j][i][k]);
betadum[5]+=gradprior3[j][i][k]*gradprior3[j][i][k];
betadum[0]+=grad->vp[j][i][k]*(grad->vp[j][i][k]-grad_prior1->vp[j][i][k]);
betadum[1]+=grad_prior1->vp[j][i][k]*grad_prior1->vp[j][i][k];
betadum[2]+=grad->vs[j][i][k]*(grad->vs[j][i][k]-grad_prior1->vs[j][i][k]);
betadum[3]+=grad_prior1->vs[j][i][k]*grad_prior1->vs[j][i][k];
betadum[4]+=grad->rho[j][i][k]*(grad->rho[j][i][k]-grad_prior1->rho[j][i][k]);
betadum[5]+=grad_prior1->rho[j][i][k]*grad_prior1->rho[j][i][k];
/*save preconditioned gradient for next iteration*/
gradprior1[j][i][k]=grad1[j][i][k];