Commit 8abb7c57 authored by Florian Wittkamp's avatar Florian Wittkamp

Higher temporal order with Adams-Basthforth method

Added a new switch FDORDER_TIME to use a fourth order accurate temporal FD-scheme. This is realized with the Adams-Bashforth method.
Possible values for FDORDER_TIME are 2 and 4, whereas 2 refers to the classical leapfrog scheme.
parent 158387de
......@@ -3,6 +3,7 @@
# Sofi2D specific
sofi2D
sofi2D_bench
snapmerge
guide_sofi2D.pdf
......
......@@ -12,7 +12,8 @@
"FD order" : "comment",
"FDORDER" : "4",
"MAXRELERROR" : "1",
"FDORDER_TIME" : "4",
"MAXRELERROR" : "0",
"2-D Grid" : "comment",
"NX" : "300",
......@@ -48,7 +49,7 @@
"WRITE_MODELFILES" : "2",
"Q-approximation" : "comment",
"L" : "0",
"L" : "1",
"FL1" : "5.0",
"TAU" : "0.00001",
......
......@@ -95,6 +95,7 @@ SOFI2D_SRC= \
operators_v.c \
PML_pro.c \
prepare_update_s.c \
prepare_update_s_4.c \
psource.c \
rd_sour.c \
read_checkpoint.c\
......@@ -115,23 +116,37 @@ SOFI2D_SRC= \
surface.c \
surface_elastic.c \
update_s_elastic_abs.c \
update_s_elastic_abs_4.c \
update_s_elastic_interior.c \
update_s_elastic_interior_4.c \
update_s_elastic_PML.c \
update_s_elastic_PML_4.c \
update_s_visc_abs.c \
update_s_visc_interior.c \
update_s_visc_PML.c \
update_s_visc_abs_4.c \
update_s_visc_interior_4.c \
update_s_visc_PML_4.c \
update_v_abs.c \
update_v_abs_4.c \
update_v_interior.c \
update_v_interior_4.c \
update_v_PML.c \
update_v_PML_4.c \
util.c \
wavefield_update_s_el.c \
wavefield_update_s_el_4.c \
wavefield_update_s_visc.c \
wavefield_update_s_visc_4.c \
wavefield_update_v.c \
wavefield_update_v_4.c \
wavelet.c \
write_par.c \
writedsk.c \
writemod.c \
zero_elastic.c \
zero_elastic_4.c \
zero_visco_4.c \
zero_visc.c \
zero_PML_elastic.c \
zero_PML_visc.c
......@@ -170,6 +185,7 @@ SOFI2D_BENCH_SRC= \
operators_v.c \
PML_pro.c \
prepare_update_s.c \
prepare_update_s_4.c \
psource.c \
rd_sour.c \
read_checkpoint.c\
......@@ -190,23 +206,37 @@ SOFI2D_BENCH_SRC= \
surface.c \
surface_elastic.c \
update_s_elastic_abs.c \
update_s_elastic_abs_4.c \
update_s_elastic_interior.c \
update_s_elastic_interior_4.c \
update_s_elastic_PML.c \
update_s_visc_abs.c \
update_s_visc_interior.c \
update_s_visc_PML.c \
update_s_elastic_PML_4.c \
update_s_visc_abs_4.c \
update_s_visc_interior_4.c \
update_s_visc_PML_4.c \
update_v_abs.c \
update_v_abs_4.c \
update_v_interior.c \
update_v_interior_4.c \
update_v_PML.c \
update_v_PML_4.c \
util.c \
wavefield_update_s_el.c \
wavefield_update_s_el_4.c \
wavefield_update_s_visc.c \
wavefield_update_s_visc_4.c \
wavefield_update_v.c \
wavefield_update_v_4.c \
wavelet.c \
write_par.c \
writedsk.c \
writemod.c \
zero_elastic.c \
zero_elastic_4.c \
zero_visco_4.c \
zero_visc.c \
zero_PML_elastic.c \
zero_PML_visc.c \
......
......@@ -35,7 +35,7 @@ void checkfd ( FILE *fp, float ** prho, float ** ppi, float ** pu,
extern float DH, DT, TS, TIME, TSNAP2;
extern float XREC1, XREC2, YREC1, YREC2;
extern int NX, NY, L, MYID, IDX, IDY, NT, NDT, RSG;
extern int READREC, NPROCX,NPROCY, SRCREC, FREE_SURF, ABS_TYPE, FW, BOUNDARY;
extern int READREC, NPROCX,NPROCY, SRCREC, FREE_SURF, ABS_TYPE, FW, BOUNDARY,FDORDER_TIME;
extern int SNAP, SEISMO, CHECKPTREAD, CHECKPTWRITE, SEIS_FORMAT[6], SNAP_FORMAT, POS[4];
extern char SEIS_FILE[STRING_SIZE], CHECKPTFILE[STRING_SIZE], SNAP_FILE[STRING_SIZE];
extern char SOURCE_FILE[STRING_SIZE], REC_FILE[STRING_SIZE];
......@@ -43,10 +43,11 @@ void checkfd ( FILE *fp, float ** prho, float ** ppi, float ** pu,
/* local variables */
float c, cmax_p=0.0, cmin_p=1e9, cmax_s=0.0, cmin_s=1e9, cwater=1.0, fmax, gamma;
float cmax=0.0, cmin=1e9, sum, dtstab, dhstab, ts, cmax_r, cmin_r;
float cmax=0.0, cmin=1e9, sum, dtstab, dhstab, ts, cmax_r, cmin_r, temporal;
float snapoutx=0.0, snapouty=0.0;
float srec_minx=DH*NX*NPROCX+1, srec_miny=DH*NY*NPROCY+1;
float srec_maxx=-1.0, srec_maxy=-1.0;
float CFL;
const float w=2.0*PI/TS; /*center frequency of source*/
int i, j, k, l, ny1=1, nx, ny, myidcounter, nfw;
......@@ -351,14 +352,14 @@ void checkfd ( FILE *fp, float ** prho, float ** ppi, float ** pu,
MPI_Allreduce ( &cmin,&cmin_r,1,MPI_FLOAT,MPI_MIN,MPI_COMM_WORLD );
cmax=cmax_r;
cmin=cmin_r;
if (FDORDER_TIME==4) {temporal=3.0/2.0;} else {temporal=1.0;}
fmax=2.0/TS;
dhstab = ( cmin/ ( hc[0]*fmax ) );
gamma = fabs ( hc[1] ) + fabs ( hc[2] ) + fabs ( hc[3] ) + fabs ( hc[4] ) + fabs ( hc[5] ) + fabs ( hc[6] );
dtstab = DH/ ( sqrt ( 2 ) *gamma*cmax );
if ( RSG ) dtstab=DH/cmax;
dtstab = DH/ ( sqrt ( 2 ) *gamma*cmax*temporal );
CFL=cmax*DT/DH;
if ( RSG ) dtstab=DH/cmax;
if ( MYID == 0 ) {
fprintf ( fp," Global values for entire model: \n" );
......@@ -388,6 +389,7 @@ void checkfd ( FILE *fp, float ** prho, float ** ppi, float ** pu,
fprintf ( fp," In the current simulation cmax is %8.2f m/s .\n\n",cmax );
fprintf ( fp," DT is the timestep and DH is the grid size.\n\n" );
fprintf ( fp," In this simulation the Courant-Friedrichs-Lewy number is %2.4f.\n",CFL );
fprintf ( fp," In this simulation the stability limit for timestep DT is %e seconds .\n",dtstab );
fprintf ( fp," You have specified DT= %e s.\n", DT );
if ( DT>dtstab )
......
......@@ -26,7 +26,7 @@
void exchange_par(void){
/* declaration of extern variables */
extern int NX, NY, FDORDER, MAXRELERROR, SOURCE_TYPE, SOURCE_SHAPE, SNAP, SNAP_FORMAT, L;
extern int NX, NY, FDORDER, MAXRELERROR, SOURCE_TYPE, SOURCE_SHAPE, SNAP, SNAP_FORMAT, L, FDORDER_TIME;
extern float DH, TIME, DT, TS, *FL, TAU, DAMPING, FPML, NPOWER, K_MAX_CPML, VPPML, PLANE_WAVE_DEPTH, PLANE_WAVE_ANGLE, SRCPOSXYZ[3];
extern float XREC1, XREC2, YREC1, YREC2;
extern float REC_ARRAY_DEPTH, REC_ARRAY_DIST;
......@@ -126,6 +126,8 @@ void exchange_par(void){
idum[34] = WRITE_MODELFILES;
idum[35] = ABS_TYPE;
idum[36] = FDORDER_TIME;
} /** if (MYID == 0) **/
......@@ -219,7 +221,10 @@ void exchange_par(void){
RUN_MULTIPLE_SHOTS = idum[33];
WRITE_MODELFILES = idum[34];
ABS_TYPE = idum[35];
ABS_TYPE = idum[35];
FDORDER_TIME = idum[36];
MPI_Bcast(&FL[1],L,MPI_FLOAT,0,MPI_COMM_WORLD);
......
......@@ -37,7 +37,7 @@
#define fsign(x) ((x<0.0)?(-1):1)
#define PI (3.141592653589793238462643383279502884197169)
#define NPAR 40
#define NPAR 41
//#define STRING_SIZE 74 // previous value, sometimes not enough to handle longer file names
#define STRING_SIZE 256
#define REQUEST_COUNT 4
......@@ -266,7 +266,10 @@ void prepare_update_s ( float *etajm, float *etaip, float *peta, float **fipjp,
float **ptausipjp, float **f, float **g, float *bip, float *bjm,
float *cip, float *cjm, float ***dip, float ***d, float ***e );
void prepare_update_s_4 ( float *etajm, float *etaip, float *peta, float **fipjp, float **pu,
float **puipjp, float **ppi, float **ptaus, float **ptaup,
float **ptausipjp, float **f, float **g, float *bip, float *bjm,
float *cip, float *cjm, float ***dip, float ***d, float ***e );
void update_s_visc_abs ( int nx1, int nx2, int ny1, int ny2, int *gx, int *gy, int nt,
float **vx, float **vy, float **sxx, float **syy,
......@@ -276,6 +279,14 @@ void update_s_visc_abs ( int nx1, int nx2, int ny1, int ny2, int *gx, int *gy, i
float *cjm, float ***d, float ***e, float ***dip,
float ** absorb_coeff,float *hc );
void update_s_visc_abs_4 ( int nx1, int nx2, int ny1, int ny2, int *gx, int *gy, int nt,
float **vx, float **vy, float **sxx, float **syy,
float **sxy, float ***r, float *** p, float ***q,
float **pi,
float ** fipjp, float **f, float **g, float *bip, float *bjm, float *cip,
float *cjm, float ***d, float ***e, float ***dip,
float ** absorb_coeff,float *hc, float ** vxx_1,float ** vxx_2,float ** vxx_3,float ** vxx_4,float ** vyy_1,float ** vyy_2,float ** vyy_3,float ** vyy_4,float ** vxy_1,float ** vxy_2,float ** vxy_3,float ** vxy_4,float ** vyx_1,float ** vyx_2,float ** vyx_3,float ** vyx_4,float ** svx_1,float ** svx_2,float ** svx_3,float ** svx_4,float ** svy_1,float ** svy_2,float ** svy_3,float ** svy_4,float ***pr_2,float ***pr_3,float ***pr_4, float ***pp_2, float ***pp_3, float ***pp_4, float ***pq_2, float ***pq_3, float ***pq_4 );
void update_s_elastic ( int nx1, int nx2, int ny1, int ny2, int nt,
float ** vx, float ** vy, float ** sxx, float ** syy,
float ** sxy, float ** pi, float ** u, float ** uipjp, float ** absorb_coeff,
......@@ -286,11 +297,20 @@ void update_s_elastic_abs ( int nx1, int nx2, int ny1, int ny2, int * gx, int *
float ** sxy, float ** pi, float ** u, float ** uipjp,
float ** absorb_coeff, float *hc);
void update_s_elastic_abs_4 ( int nx1, int nx2, int ny1, int ny2, int * gx, int * gy, int nt,
float ** vx, float ** vy, float ** sxx, float ** syy,
float ** sxy, float ** pi, float ** u, float ** uipjp,
float ** absorb_coeff, float *hc ,float ** vxx_1,float ** vxx_2,float ** vxx_3,float ** vxx_4,float ** vyy_1,float ** vyy_2,float ** vyy_3,float ** vyy_4,float ** vxy_1,float ** vxy_2,float ** vxy_3,float ** vxy_4,float ** vyx_1,float ** vyx_2,float ** vyx_3,float ** vyx_4);
void update_s_elastic_interior ( int nx1, int nx2, int ny1, int ny2, int * gx, int * gy, int nt,
float ** vx, float ** vy, float ** sxx, float ** syy,
float ** sxy, float ** pi, float ** u, float ** uipjp,
float *hc );
void update_s_elastic_interior_4 ( int nx1, int nx2, int ny1, int ny2, int * gx, int * gy, int nt,
float ** vx, float ** vy, float ** sxx, float ** syy,
float ** sxy, float ** pi, float ** u, float ** uipjp, float *hc,float ** vxx_1,float ** vxx_2,float ** vxx_3,float ** vxx_4,float ** vyy_1,float ** vyy_2,float ** vyy_3,float ** vyy_4,float ** vxy_1,float ** vxy_2,float ** vxy_3,float ** vxy_4,float ** vyx_1,float ** vyx_2,float ** vyx_3,float ** vyx_4);
void update_s_elastic_PML ( int nx1, int nx2, int ny1, int ny2, int * gx, int * gy, int nt,
float ** vx, float ** vy, float ** sxx, float ** syy,
float ** sxy, float ** pi, float ** u, float ** uipjp, float *hc,
......@@ -298,6 +318,13 @@ void update_s_elastic_PML ( int nx1, int nx2, int ny1, int ny2, int * gx, int *
float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half,
float ** psi_vxx, float ** psi_vyy, float ** psi_vxy, float ** psi_vyx );
void update_s_elastic_PML_4 ( int nx1, int nx2, int ny1, int ny2, int * gx, int * gy, int nt,
float ** vx, float ** vy, float ** sxx, float ** syy,
float ** sxy, float ** pi, float ** u, float ** uipjp, float *hc,
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 ** psi_vxx, float ** psi_vyy, float ** psi_vxy, float ** psi_vyx,float ** vxx_1,float ** vxx_2,float ** vxx_3,float ** vxx_4,float ** vyy_1,float ** vyy_2,float ** vyy_3,float ** vyy_4,float ** vxy_1,float ** vxy_2,float ** vxy_3,float ** vxy_4,float ** vyx_1,float ** vyx_2,float ** vyx_3,float ** vyx_4 );
void update_s_visc_interior(int nx1, int nx2, int ny1, int ny2, int *gx, int *gy, int nt,
float **vx, float **vy, float **sxx, float **syy,
float **sxy, float ***r, float *** p, float ***q,
......@@ -313,6 +340,20 @@ void update_s_visc_PML ( int nx1, int nx2, int ny1, int ny2, int *gx, int *gy,in
float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half,
float ** psi_vxx, float ** psi_vyy, float ** psi_vxy, float ** psi_vyx );
void update_s_visc_interior_4(int nx1, int nx2, int ny1, int ny2, int *gx, int *gy, int nt,
float **vx, float **vy, float **sxx, float **syy,
float **sxy, float ***r, float *** p, float ***q,
float ** fipjp, float **f, float **g, float *bip, float *bjm, float *cip,
float *cjm, float ***d, float ***e, float ***dip,
float *hc, float ** vxx_1,float ** vxx_2,float ** vxx_3,float ** vxx_4,float ** vyy_1,float ** vyy_2,float ** vyy_3,float ** vyy_4,float ** vxy_1,float ** vxy_2,float ** vxy_3,float ** vxy_4,float ** vyx_1,float ** vyx_2,float ** vyx_3,float ** vyx_4,float ** svx_1,float ** svx_2,float ** svx_3,float ** svx_4,float ** svy_1,float ** svy_2,float ** svy_3,float ** svy_4,float ***pr_2,float ***pr_3,float ***pr_4, float ***pp_2, float ***pp_3, float ***pp_4, float ***pq_2, float ***pq_3, float ***pq_4);
void update_s_visc_PML_4 ( int nx1, int nx2, int ny1, int ny2, int *gx, int *gy,int nt,
float ** vx, float ** vy, float ** sxx, float ** syy,
float ** sxy, float *hc, float ***r, float ***p, float ***q, float **fipjp, float **f, float **g,
float *bip, float *bjm, float *cip, float *cjm, float ***d, float ***e, float ***dip,
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 ** psi_vxx, float ** psi_vyy, float ** psi_vxy, float ** psi_vyx, float ** vxx_1,float ** vxx_2,float ** vxx_3,float ** vxx_4,float ** vyy_1,float ** vyy_2,float ** vyy_3,float ** vyy_4,float ** vxy_1,float ** vxy_2,float ** vxy_3,float ** vxy_4,float ** vyx_1,float ** vyx_2,float ** vyx_3,float ** vyx_4,float ** svx_1,float ** svx_2,float ** svx_3,float ** svx_4,float ** svy_1,float ** svy_2,float ** svy_3,float ** svy_4,float ***pr_2,float ***pr_3,float ***pr_4, float ***pp_2, float ***pp_3, float ***pp_4, float ***pq_2, float ***pq_3, float ***pq_4 );
void PML_pro ( float * d_x, float * K_x, float * alpha_prime_x, float * a_x, float * b_x,
float * d_x_half, float * K_x_half, float * alpha_prime_x_half, float * a_x_half, float * b_x_half,
......@@ -330,18 +371,31 @@ void update_v_abs ( int nx1, int nx2, int ny1, int ny2, int * gx, int * gy, int
float ** sxx, float ** syy, float ** sxy, float **rip, float **rjp,
float ** absorb_coeff,float *hc);
void update_v_abs_4 ( int nx1, int nx2, int ny1, int ny2, int * gx, int * gy, int nt,
float ** vx, float ** vy, float ** sxx, float ** syy, float ** sxy,
float **rip, float **rjp, float ** absorb_coeff,float *hc,float ** svx_1,float ** svx_2,float ** svx_3,float ** svx_4,float ** svy_1,float ** svy_2,float ** svy_3,float ** svy_4);
void update_v_interior(int nx1, int nx2, int ny1, int ny2, int *gx, int *gy, int nt,
float ** vx, float ** vy, float ** sxx, float ** syy,
float ** sxy, float **rho, float **rip, float **rjp,
float ** srcpos_loc, float ** signals, int nsrc,
float *hc);
void update_v_interior_4 ( int nx1, int nx2, int ny1, int ny2, int *gx, int *gy, int nt,
float ** vx, float ** vy, float ** sxx, float ** syy,
float ** sxy, float **rho, float **rip, float **rjp,
float ** srcpos_loc, float ** signals, int nsrc,float *hc,float ** svx_1,float ** svx_2,float ** svx_3,float ** svx_4,float ** svy_1,float ** svy_2,float ** svy_3,float ** svy_4);
void update_v_PML ( int nx1, int nx2, int ny1, int ny2,int * gx, int * gy, int nt, float ** vx, float ** vy,
float ** sxx, float ** syy, float ** sxy, float **rip, float **rjp,
float *hc, 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 ** psi_sxx_x, float ** psi_syy_y,
float ** psi_sxy_y, float ** psi_syx_x );
void update_v_PML_4 ( int nx1, int nx2, int ny1, int ny2, int * gx, int * gy, int nt, float ** vx, float ** vy,
float ** sxx, float ** syy, float ** sxy, float **rip, float **rjp,
float *hc, 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 ** psi_sxx_x, float ** psi_syy_y, float ** psi_sxy_y, float ** psi_sxy_x,float ** svx_1,float ** svx_2,float ** svx_3,float ** svx_4,float ** svy_1,float ** svy_2,float ** svy_3,float ** svy_4);
void wavefield_update_s_el ( int i, int j,float vxx, float vyx,float vxy,float vyy, float **sxy,
float **sxx, float ** syy, float ** pi, float ** u, float ** uipjp );
......@@ -351,9 +405,18 @@ void wavefield_update_s_visc ( int i, int j,float vxx, float vyx,float vxy,fl
float ***q,float **fipjp, float **f, float **g, float *bip,
float *bjm,float *cip, float *cjm, float ***d, float ***e, float ***dip );
void wavefield_update_s_visc_4 ( int i, int j,float vxx, float vyx,float vxy,float vyy, float **sxy,
float **sxx, float ** syy, float ***r, float ***p,
float ***q,float **fipjp, float **f, float **g, float *bip,
float *bjm,float *cip, float *cjm, float ***d, float ***e, float ***dip, float ** vxx_1,float ** vxx_2,float ** vxx_3,float ** vxx_4,float ** vyy_1,float ** vyy_2,float ** vyy_3,float ** vyy_4,float ** vxy_1,float ** vxy_2,float ** vxy_3,float ** vxy_4,float ** vyx_1,float ** vyx_2,float ** vyx_3,float ** vyx_4,float ** svx_1,float ** svx_2,float ** svx_3,float ** svx_4,float ** svy_1,float ** svy_2,float ** svy_3,float ** svy_4,float ***pr_2,float ***pr_3,float ***pr_4, float ***pp_2, float ***pp_3, float ***pp_4, float ***pq_2, float ***pq_3, float ***pq_4 );
void wavefield_update_v ( int i, int j,float sxx_x, float sxy_x,float sxy_y,float syy_y, float **vx,
float **vy, float ** rip, float ** rjp );
void wavefield_update_v_4 ( int i, int j,float sxx_x, float sxy_x,float sxy_y,float syy_y, float **vx,
float **vy, float ** rip, float ** rjp,float ** svx_1,float ** svx_2,float ** svx_3,float ** svx_4,float ** svy_1,float ** svy_2,float ** svy_3,float ** svy_4);
void wavefield_update_s_el_4 ( int i, int j,float vxx, float vyx,float vxy,float vyy, float **sxy,
float **sxx, float ** syy, float ** pi, float ** u, float ** uipjp,float ** vxx_1,float ** vxx_2,float ** vxx_3,float ** vxx_4,float ** vyy_1,float ** vyy_2,float ** vyy_3,float ** vyy_4,float ** vxy_1,float ** vxy_2,float ** vxy_3,float ** vxy_4,float ** vyx_1,float ** vyx_2,float ** vyx_3,float ** vyx_4);
float ** wavelet ( float ** srcpos_loc, int nsrc );
......@@ -373,6 +436,10 @@ void writemod ( char modfile[STRING_SIZE], float ** array, int format );
void zero_elastic ( int nx1, int nx2, int ny1, int ny2, float ** vx, float ** vy, float ** sxx, float ** syy, float ** sxy );
void zero_elastic_4 ( int nx1, int nx2, int ny1, int ny2, float ** vxx_1,float ** vxx_2,float ** vxx_3,float ** vxx_4,float ** vyy_1,float ** vyy_2,float ** vyy_3,float ** vyy_4,float ** vxy_1,float ** vxy_2,float ** vxy_3,float ** vxy_4,float ** vyx_1,float ** vyx_2,float ** vyx_3,float ** vyx_4,float ** svx_1,float ** svx_2,float ** svx_3,float ** svx_4,float ** svy_1,float ** svy_2,float ** svy_3,float ** svy_4);
void zero_visco_4 ( int nx1, int nx2, int ny1, int ny2, float ** vxx_1,float ** vxx_2,float ** vxx_3,float ** vxx_4,float ** vyy_1,float ** vyy_2,float ** vyy_3,float ** vyy_4,float ** vxy_1,float ** vxy_2,float ** vxy_3,float ** vxy_4,float ** vyx_1,float ** vyx_2,float ** vyx_3,float ** vyx_4,float ** svx_1,float ** svx_2,float ** svx_3,float ** svx_4,float ** svy_1,float ** svy_2,float ** svy_3,float ** svy_4,float ***pr_2,float ***pr_3,float ***pr_4, float ***pp_2, float ***pp_3, float ***pp_4, float ***pq_2, float ***pq_3, float ***pq_4);
void zero_visc ( int nx1, int nx2, int ny1, int ny2, float ** vx, float ** vy, float ** sxx, float ** syy, float ** sxy, float *** pr, float *** pp, float *** pq );
void zero_PML_elastic ( int ny1, int ny2, int nx1, int nx2, float ** vx, float ** vy, float ** sxx,
......
......@@ -33,7 +33,7 @@ int RUNMODE, WRITE_MODELFILES=0, ABS_TYPE;
int OUTNTIMESTEPINFO=1; /*every OUTNTIMESTEPINFO th timestep, information on the time step will be given to screen/file */
int SEISMO=0, NDT=1, NSRC=1, SEIS_FORMAT=0, FREE_SURF=0, READMOD=0, READREC=0, SRCREC=0, RSG=0, FW=0;
int NX=1, NY=1, NT=0, SOURCE_TYPE=0, SOURCE_SHAPE=0, SNAP=0, SNAP_FORMAT=0, LOG=0, REC_ARRAY=0;
int L=0, BOUNDARY=0, DC=0, DRX=0, NXG=0, NYG=0, IDX=1, IDY=1, CHECKPTREAD=0, CHECKPTWRITE=0, FDORDER=0, MAXRELERROR=0;
int L=0, BOUNDARY=0, DC=0, DRX=0, NXG=0, NYG=0, IDX=1, IDY=1, CHECKPTREAD=0, CHECKPTWRITE=0, FDORDER=0, FDORDER_TIME=0, MAXRELERROR=0;
int RUN_MULTIPLE_SHOTS=0; /* Added for multiple shots */
char SNAP_FILE[STRING_SIZE]="", SOURCE_FILE[STRING_SIZE]="", SIGNAL_FILE[STRING_SIZE]="";
char MFILE[STRING_SIZE]="", REC_FILE[STRING_SIZE]="", CHECKPTFILE[STRING_SIZE]="";
......
/*------------------------------------------------------------------------
* Copyright (C) 2011 For the list of authors, see file AUTHORS.
*
* This file is part of SOFI2D.
*
* SOFI2D is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, version 2.0 of the License only.
*
* SOFI2D is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with SOFI2D. See file COPYING and/or
* <http://www.gnu.org/licenses/gpl-2.0.html>.
--------------------------------------------------------------------------*/
/* ------------------------------------------------------------------------
* prepare of update of the stress tensor
* ------------------------------------------------------------------------*/
/* ------------------------------------------------------------------------
* ATTENTION: The parameters below will be scaled by factor c* due to
* Adams-Bashforth method, so be aware and only call this function when
* FDORDER_TIME is set to 4
* ------------------------------------------------------------------------*/
#include "fd.h"
void prepare_update_s_4(float *etajm, float *etaip, float *peta, float **fipjp, float **pu,
float **puipjp, float **ppi, float **ptaus, float **ptaup,
float **ptausipjp, float **f, float **g, float *bip, float *bjm,
float *cip, float *cjm, float ***dip, float ***d, float ***e) {
extern int NX, NY, L;
extern float DT;
int i, j, l;
float c1; /* Coefficients for Adam Bashforth */
c1=13.0/12.0;
for (l=1;l<=L;l++){
etajm[l] = peta[l];
etaip[l] = peta[l];
}
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
fipjp[j][i] = puipjp[j][i]*DT*(1.0+L*ptausipjp[j][i]);
f[j][i] = pu[j][i]*DT*(1.0+L*ptaus[j][i]);
g[j][i] = ppi[j][i]*DT*(1.0+L*ptaup[j][i]);
for (l=1;l<=L;l++){
bip[l] = 1.0/(1.0+(c1*etaip[l]*0.5));
bjm[l] = 1.0/(1.0+(c1*etajm[l]*0.5));
cip[l] = 1.0-(c1*etaip[l]*0.5);
cjm[l] = 1.0-(c1*etajm[l]*0.5);
dip[j][i][l] = puipjp[j][i]*etaip[l]*ptausipjp[j][i];
d[j][i][l] = pu[j][i]*etajm[l]*ptaus[j][i];
e[j][i][l] = ppi[j][i]*etajm[l]*ptaup[j][i];
}
}
}
}
\ No newline at end of file
......@@ -30,7 +30,7 @@ char ** varname_list,** value_list;
void read_par_json(FILE *fp, char *fileinp){
/* declaration of extern variables */
extern int NX, NY, FDORDER, MAXRELERROR, SOURCE_TYPE, SOURCE_SHAPE, SNAP, SNAP_FORMAT, L;
extern int NX, NY, FDORDER, FDORDER_TIME, MAXRELERROR, SOURCE_TYPE, SOURCE_SHAPE, SNAP, SNAP_FORMAT, L;
extern int SEISMO, NDT, SEIS_FORMAT, FREE_SURF, READMOD, READREC, SRCREC;
extern int BOUNDARY, REC_ARRAY, DRX, LOG, WRITE_MODELFILES; //RSG
extern int NPROCX, NPROCY, MYID, IDX, IDY, CHECKPTREAD, CHECKPTWRITE, RUN_MULTIPLE_SHOTS, ABS_TYPE, FW;
......@@ -81,6 +81,13 @@ void read_par_json(FILE *fp, char *fileinp){
err("Variable RSG could not be retrieved from the json input file!");*/
if (get_int_from_objectlist("FDORDER",number_readobjects,&FDORDER,varname_list, value_list))
err("Variable FDORDER could not be retrieved from the json input file!");
if (get_int_from_objectlist("FDORDER_TIME",number_readobjects,&FDORDER_TIME,varname_list, value_list)) {
FDORDER_TIME=2;
} else {
if(FDORDER_TIME!=2 && FDORDER_TIME!=4) {
err("Only FDORDER_TIME 2 or 4 are supported!");
}
}
if (get_int_from_objectlist("MAXRELERROR",number_readobjects,&MAXRELERROR,varname_list, value_list))
err("Variable MAXRELERROR could not be retrieved from the json input file!");
if (get_int_from_objectlist("NX",number_readobjects,&NX,varname_list, value_list))
......@@ -261,9 +268,9 @@ void read_par_json(FILE *fp, char *fileinp){
err("Variable TAU could not be retrieved from the json input file!");
if (get_int_from_objectlist("L",number_readobjects,&L,varname_list, value_list))
err("Variable L could not be retrieved from the json input file!");
else {
FL=vector(1,L);
switch(L) {
else {
FL=vector(1,L);
switch(L) {
case 0:
break;
case 5:
......
......@@ -14,8 +14,8 @@
*
* You should have received a copy of the GNU General Public License
* along with SOFI2D. See file COPYING and/or
* <http://www.gnu.org/licenses/gpl-2.0.html>.
--------------------------------------------------------------------------*/
* <http://www.gnu.org/licenses/gpl-2.0.html>.
--------------------------------------------------------------------------*/
/* ----------------------------------------------------------------------
* This is program SOFI2D.
* Parallel 2-D Viscoelastic Finite Difference Seismic Modeling
......@@ -46,925 +46,1111 @@
int main ( int argc, char **argv )
{
/* variables in main */
int ns, nseismograms = 0, nt, nd, fdo3;
int lsnap, nsnap = 0, lsamp = 0, buffsize;
int ntr = 0, ntr_loc = 0, ntr_glob = 0, nsrc = 0, nsrc_loc = 0;
int ishot, nshots; /* Added ishot and nshots for multiple shots */
/*Limits for local grids defined in subgrid_bounds.c */
int * gx=NULL, * gy=NULL;
float memdyn, memmodel, memseismograms, membuffer, memtotal, memcpml=0.0;
float fac1, fac2;
char *buff_addr, ext[10], *fileinp = "", modestr[10], infostr[70];
double time1 = 0.0, time2 = 0.0, time3 = 0.0, time4 = 0.0;
double time5 = 0.0, time6 = 0.0, time7 = 0.0, time8 = 0.0;
double time_av_v_update = 0.0, time_av_s_update = 0.0,
time_av_v_exchange = 0.0, time_av_s_exchange = 0.0, time_av_timestep = 0.0;
float **psxx = NULL, **psxy = NULL, **psyy = NULL;
float **pvx = NULL, **pvy = NULL, ***pr = NULL;
float ***pp = NULL, ***pq = NULL;
float **pu = NULL, **puipjp = NULL, **ptaus = NULL, **ptaup = NULL,
*etaip = NULL, *etajm = NULL, *peta = NULL, **ptausipjp = NULL,
**fipjp = NULL, ***dip = NULL, *bip = NULL, *bjm = NULL, *cip = NULL,
*cjm = NULL, ***d = NULL, ***e = NULL, **f = NULL, **g = NULL;
float **prho = NULL, **prip = NULL, **prjp = NULL, **ppi = NULL;
float **sectionvx = NULL, **sectionvy = NULL, **sectionp = NULL,
**sectioncurl = NULL, **sectiondiv = NULL;
float **absorb_coeff = NULL;
float **srcpos = NULL, **srcpos_loc = NULL, **signals = NULL, *hc = NULL,
**srcpos_current = NULL;
int **recpos = NULL, **recpos_loc = NULL;
float **bufferlef_to_rig = NULL, **bufferrig_to_lef = NULL,
**buffertop_to_bot = NULL, **bufferbot_to_top = NULL;
float ** seismo_fulldata=NULL;
int * recswitch = NULL;
/* PML variables */
float * d_x=NULL, * K_x=NULL, * alpha_prime_x=NULL, * a_x=NULL, * b_x=NULL, * d_x_half=NULL,
* K_x_half=NULL, * alpha_prime_x_half=NULL, * a_x_half=NULL, * b_x_half=NULL,
* d_y=NULL, * K_y=NULL, * alpha_prime_y=NULL, * a_y=NULL, * b_y=NULL, * d_y_half=NULL,
* K_y_half=NULL, * alpha_prime_y_half=NULL, * a_y_half=NULL, * b_y_half=NULL;
float ** psi_sxx_x=NULL, ** psi_syy_y=NULL, ** psi_sxy_y=NULL, ** psi_sxy_x=NULL,
** psi_vxx=NULL, ** psi_vyy=NULL, ** psi_vxy=NULL, ** psi_vyx=NULL, ** psi_vxxs=NULL;
FILE *fpinp;
MPI_Request *req_send, *req_rec;
/* MPI_Status *send_statuses, *rec_statuses; */
/* Initialize MPI environment */
MPI_Init ( &argc, &argv );
MPI_Comm_size ( MPI_COMM_WORLD, &NP );
MPI_Comm_rank ( MPI_COMM_WORLD, &MYID );
if ( MYID == 0 ) {
time1 = MPI_Wtime();
clock();
}
/* print program name, version etc to stdout*/
if ( MYID == 0 )
info ( stdout );
/* =================================================== */
/* check of parameter-file can be opened*/
fileinp = argv[1];
fpinp = fopen ( fileinp, "r" );
if ( fpinp == NULL ) {
if ( MYID == 0 ) {
printf (
"\n==================================================================\n" );
printf ( " Cannot open sofi2D input file %s \n", fileinp );
printf (
"\n==================================================================\n\n" );
err ( " --- " );
return 0;
}
} else {
fscanf ( fpinp, "%s %s = %i", infostr, modestr, &RUNMODE );
fclose ( fpinp );
}
/* Check file system */
check_fs ( stdout, argc, fileinp );
/* =================================================== */
if ( RUNMODE == 0 ) {
/* read standard input file */
if ( strstr ( fileinp, ".json" ) )
//read json formated input file
read_par_json ( stdout, fileinp );
else {
if ( MYID == 0 )
err ( " Old Input files (.inp) are no longer supported. \n Please use .json input files instead. " );
}
}
/*else
auto mode: read input files
read_par_auto ( stdout, fileinp ); */
/* =================================================== */
exchange_par();
/* open log-file (each PE is using different file) */
sprintf ( ext, ".%i", MYID );
strcat ( LOG_FILE, ext );
/* nodes MYIDo writes logging info to LOG_FILE or stdout */
if ( MYID == 0 )
switch ( LOG ) {
case 0:
FP = fopen ( "/dev/null", "w" ); /* no logging information will be output */
break;
case 1:
FP = stdout; /* logging information will be written to standard output */
break;
case 2:
if ( ( FP = fopen ( LOG_FILE, "w" ) ) == NULL )
err ( " Opening log-file failed." );
/* logging information will be written to LOG_FILE */
break;
}
/* all other nodes write logging info to LOG_FILE */
if ( MYID > 0 ) {
if ( ( FP = fopen ( LOG_FILE, "w" ) ) == NULL )
err ( " Opening log-file failed." );
fprintf ( FP, " This is the log-file %s generated by PE %d \n\n",
LOG_FILE, MYID );
}
if ( MYID == 0 )
note ( FP );
/* domain decomposition */
initproc();
NT = iround ( TIME/DT ); /* number of timesteps */
ns = iround ( NT/NDT ); /* number of samples per trace */
lsnap = iround ( TSNAP1/DT ); /* first snapshot at this timestep */
lsamp = NDT;
/* output of parameters to log-file or stdout */
if ( MYID == 0 )
write_par ( FP );
/* For the Rotated Staggered Grid only second order FD operators are implemented
if ( RSG ) {
if ( FDORDER > 2 )
err ( " For the Rotated Staggered Grid only second order FD operators are implemented. Please revise parameter FDORDER in the input file! " );
}*/
/* NXG, NYG denote size of the entire (global) grid */
NXG = NX;
NYG = NY;
/* In the following, NX and NY denote size of the local grid ! */
NX = IENDX;
NY = IENDY;
if ( SEISMO ) {
recpos = receiver ( FP, &ntr );
recswitch = ivector ( 1, ntr );
recpos_loc = splitrec ( recpos, &ntr_loc, ntr, recswitch );
ntr_glob = ntr;
ntr = ntr_loc;
}
/* allocate buffer for seismogram output, merged seismogram section of all PEs */
if ( SEISMO ) seismo_fulldata=matrix ( 1,ntr_glob,1,ns );
/* number of seismogram sections which have to be stored in core memory*/
/* allocation of memory for seismogramm merge */
switch ( SEISMO ) {
case 1: /* particle velocities only */
nseismograms = 3;
break;
case 2: /* pressure only */
nseismograms = 1;
break;
case 3: /* curl and div only */
nseismograms = 2;
break;
case 4: /* everything */
nseismograms = 6;
break;
default:
nseismograms = 1;
break;
}
/*allocate memory for dynamic, static and buffer arrays */
nd = FDORDER / 2 + 1;
fdo3 = 2 * nd;