Commit 7e8cc012 authored by Tilman Steinweg's avatar Tilman Steinweg

renamed err function

(err is part of an standard library)
parent 74df1237
......@@ -219,7 +219,7 @@ void check_fs(FILE *fp, int argc, char *fileinp)
{
fprintf(fp, "\n");
sprintf(errmsg, "\n in: <check_fs.c> \n");
err(errmsg);
declare_error(errmsg);
}
}
......@@ -27,9 +27,8 @@
#include "fd.h"
void checkfd ( FILE *fp, float ** prho, float ** ppi, float ** pu,
float ** ptaus, float ** ptaup, float *peta, float *hc, float **srcpos, int nsrc, int **recpos, int ntr )
{
void checkfd(FILE *fp, float **prho, float **ppi, float **pu,
float **ptaus, float **ptaup, float *peta, float *hc, float **srcpos, int nsrc, int **recpos, int ntr) {
extern float DH, DT, TS, TIME, TSNAP2;
......@@ -47,7 +46,7 @@ void checkfd ( FILE *fp, float ** prho, float ** ppi, float ** pu,
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;
float CFL;
const float w=2.0*PI/TS; /*center frequency of source*/
int i, j, k, l, ny1=1, nx, ny, myidcounter, nfw;
......@@ -59,11 +58,11 @@ void checkfd ( FILE *fp, float ** prho, float ** ppi, float ** pu,
nx=NX;
ny=NY;
fprintf ( fp,"\n **********************************************************" );
fprintf ( fp,"\n ************ CHECKS OF INPUT FILE PARAMETERS ************" );
fprintf ( fp,"\n **********************************************************\n\n" );
fprintf ( fp,"\n **Message from checkfd_ssg (printed by PE %d):\n",MYID );
fprintf ( fp,"\n\n ------------------ CHECK OUTPUT FILES --------------------------\n" );
fprintf(fp,"\n **********************************************************");
fprintf(fp,"\n ************ CHECKS OF INPUT FILE PARAMETERS ************");
fprintf(fp,"\n **********************************************************\n\n");
fprintf(fp,"\n **Message from checkfd_ssg (printed by PE %d):\n",MYID);
fprintf(fp,"\n\n ------------------ CHECK OUTPUT FILES --------------------------\n");
/* The original checks might delete files accidentally that would not be overwritten anyway.
and did not test accessibility from all CPUs which may be vary, especially in distributed clusters */
......@@ -73,69 +72,98 @@ void checkfd ( FILE *fp, float ** prho, float ** ppi, float ** pu,
/* only MYID is performing the file checks, if to many PEs try to do it simultaneously,
* the code and/or FILESYSTEM and/or MPI implementation will cause segmentation faults
*/
if ( ( SNAP>0 ) && ( MYID==0 ) ) {
switch ( SNAP_FORMAT ) {
case 1:
sprintf ( file_ext, ".su" );
strcpy ( xmod, "ab" );
break;
case 2:
sprintf ( file_ext, ".asc" );
strcpy ( xmod, "a" );
break;
case 3:
sprintf ( file_ext, ".bin" );
strcpy ( xmod, "ab" );
break;
default:
err ( " Sorry. Snapshot format (SNAP_FORMAT) unknown. \n" );
break;
if ((SNAP>0) && (MYID==0)) {
switch (SNAP_FORMAT) {
case 1:
sprintf(file_ext, ".su");
strcpy(xmod, "ab");
break;
case 2:
sprintf(file_ext, ".asc");
strcpy(xmod, "a");
break;
case 3:
sprintf(file_ext, ".bin");
strcpy(xmod, "ab");
break;
default:
declare_error(" Sorry. Snapshot format (SNAP_FORMAT) unknown. \n");
break;
}
fprintf ( fp," Check accessibility for snapshot files ... \n" );
switch ( SNAP ) {
case 1 :
sprintf ( xfile,"%s%s.vx.%i.%i",SNAP_FILE,file_ext,POS[1],POS[2] );
fprintf ( fp," Check accessibility for snapshot file %s... \n",xfile );
if ( ( fpcheck=fopen ( xfile,xmod ) ) ==NULL ) err2 ( " PE0 cannot write snapshots to %s!",xfile );
else fclose ( fpcheck );
sprintf ( xfile,"%s%s.vy.%i.%i",SNAP_FILE,file_ext,POS[1],POS[2] );
fprintf ( fp," Check accessibility for snapshot file %s... \n",xfile );
if ( ( fpcheck=fopen ( xfile,xmod ) ) ==NULL ) err2 ( " PE0 cannot write snapshots to %s!",xfile );
else fclose ( fpcheck );
break;
case 2 :
sprintf ( xfile,"%s%s.p.%i.%i",SNAP_FILE,file_ext,POS[1],POS[2] );
fprintf ( fp," Check accessibility for snapshot file %s... \n",xfile );
if ( ( fpcheck=fopen ( xfile,xmod ) ) ==NULL ) err2 ( " PE0 cannot write snapshots to %s!",xfile );
else fclose ( fpcheck );
break;
case 4 :
sprintf ( xfile,"%s%s.vx.%i.%i",SNAP_FILE,file_ext,POS[1],POS[2] );
fprintf ( fp," Check accessibility for snapshot file %s... \n",xfile );
if ( ( fpcheck=fopen ( xfile,xmod ) ) ==NULL ) err2 ( " PE0 cannot write snapshots to %s!",xfile );
else fclose ( fpcheck );
sprintf ( xfile,"%s%s.vy.%i.%i",SNAP_FILE,file_ext,POS[1],POS[2] );
fprintf ( fp," Check accessibility for snapshot file %s... \n",xfile );
if ( ( fpcheck=fopen ( xfile,xmod ) ) ==NULL ) err2 ( " PE0 cannot write snapshots to %s!",xfile );
else fclose ( fpcheck );
sprintf ( xfile,"%s%s.p.%i.%i",SNAP_FILE,file_ext,POS[1],POS[2] );
fprintf ( fp," Check accessibility for snapshot file %s... \n",xfile );
if ( ( fpcheck=fopen ( xfile,xmod ) ) ==NULL ) err2 ( " PE0 cannot write snapshots to %s!",xfile );
else fclose ( fpcheck );
case 3 :
sprintf ( xfile,"%s%s.div.%i.%i",SNAP_FILE,file_ext,POS[1],POS[2] );
fprintf ( fp," Check accessibility for snapshot file %s... \n",xfile );
if ( ( fpcheck=fopen ( xfile,xmod ) ) ==NULL ) err2 ( " PE0 cannot write snapshots to %s!",xfile );
else fclose ( fpcheck );
sprintf ( xfile,"%s%s.curl.%i.%i",SNAP_FILE,file_ext,POS[1],POS[2] );
fprintf ( fp," Check accessibility for snapshot file %s... \n",xfile );
if ( ( fpcheck=fopen ( xfile,xmod ) ) ==NULL ) err2 ( " PE0 cannot write snapshots to %s!",xfile );
else fclose ( fpcheck );
break;
fprintf(fp," Check accessibility for snapshot files ... \n");
switch (SNAP) {
case 1 :
sprintf(xfile,"%s%s.vx.%i.%i",SNAP_FILE,file_ext,POS[1],POS[2]);
fprintf(fp," Check accessibility for snapshot file %s... \n",xfile);
if ((fpcheck=fopen(xfile,xmod)) ==NULL) err2(" PE0 cannot write snapshots to %s!",xfile);
else fclose(fpcheck);
sprintf(xfile,"%s%s.vy.%i.%i",SNAP_FILE,file_ext,POS[1],POS[2]);
fprintf(fp," Check accessibility for snapshot file %s... \n",xfile);
if ((fpcheck=fopen(xfile,xmod)) ==NULL) err2(" PE0 cannot write snapshots to %s!",xfile);
else fclose(fpcheck);
break;
case 2 :
sprintf(xfile,"%s%s.p.%i.%i",SNAP_FILE,file_ext,POS[1],POS[2]);
fprintf(fp," Check accessibility for snapshot file %s... \n",xfile);
if ((fpcheck=fopen(xfile,xmod)) ==NULL) err2(" PE0 cannot write snapshots to %s!",xfile);
else fclose(fpcheck);
break;
case 4 :
sprintf(xfile,"%s%s.vx.%i.%i",SNAP_FILE,file_ext,POS[1],POS[2]);
fprintf(fp," Check accessibility for snapshot file %s... \n",xfile);
if ((fpcheck=fopen(xfile,xmod)) ==NULL) err2(" PE0 cannot write snapshots to %s!",xfile);
else fclose(fpcheck);
sprintf(xfile,"%s%s.vy.%i.%i",SNAP_FILE,file_ext,POS[1],POS[2]);
fprintf(fp," Check accessibility for snapshot file %s... \n",xfile);
if ((fpcheck=fopen(xfile,xmod)) ==NULL) err2(" PE0 cannot write snapshots to %s!",xfile);
else fclose(fpcheck);
sprintf(xfile,"%s%s.p.%i.%i",SNAP_FILE,file_ext,POS[1],POS[2]);
fprintf(fp," Check accessibility for snapshot file %s... \n",xfile);
if ((fpcheck=fopen(xfile,xmod)) ==NULL) err2(" PE0 cannot write snapshots to %s!",xfile);
else fclose(fpcheck);
case 3 :
sprintf(xfile,"%s%s.div.%i.%i",SNAP_FILE,file_ext,POS[1],POS[2]);
fprintf(fp," Check accessibility for snapshot file %s... \n",xfile);
if ((fpcheck=fopen(xfile,xmod)) ==NULL) err2(" PE0 cannot write snapshots to %s!",xfile);
else fclose(fpcheck);
sprintf(xfile,"%s%s.curl.%i.%i",SNAP_FILE,file_ext,POS[1],POS[2]);
fprintf(fp," Check accessibility for snapshot file %s... \n",xfile);
if ((fpcheck=fopen(xfile,xmod)) ==NULL) err2(" PE0 cannot write snapshots to %s!",xfile);
else fclose(fpcheck);
break;
}
}
......@@ -145,100 +173,145 @@ void checkfd ( FILE *fp, float ** prho, float ** ppi, float ** pu,
/* only MYID is performing the file checks, if to many PEs try to do it simultaneously,
* the code and/or FILESYSTEM and/or MPI implementation will cause segmentation faults
*/
if ( ( SEISMO>0 ) && ( MYID==0 ) ) {
fprintf ( fp," Check accessibility for seismogram files of each PE ... \n" );
fprintf ( fp," However, the list below refers only to PE0 ... \n" );
switch ( SEIS_FORMAT[0] ) {
case 1:
sprintf ( file_ext,"su" );
break;
case 2:
sprintf ( file_ext,"txt" );
break;
case 3:
sprintf ( file_ext,"bin" );
break;
if ((SEISMO>0) && (MYID==0)) {
fprintf(fp," Check accessibility for seismogram files of each PE ... \n");
fprintf(fp," However, the list below refers only to PE0 ... \n");
switch (SEIS_FORMAT[0]) {
case 1:
sprintf(file_ext,"su");
break;
case 2:
sprintf(file_ext,"txt");
break;
case 3:
sprintf(file_ext,"bin");
break;
}
if ( SEIS_FORMAT[0]==2 ) strcpy ( xmod,"a" );
else strcpy ( xmod,"w+b" );
if (SEIS_FORMAT[0]==2) strcpy(xmod,"a");
else strcpy(xmod,"w+b");
/*MYID=0 is checking if all seismogram file written by other PEs can be written
* assuming that all PEs can address the files ystem equally
*/
for ( myidcounter=0; myidcounter< ( NPROCX*NPROCY ); myidcounter++ ) {
switch ( SEISMO ) {
case 1: /* particle velocities only */
sprintf ( xfile,"%s_vx.%s.%d",SEIS_FILE,file_ext,myidcounter );
/*in case of number of PE's=500, there will be 500 messages, too many to be displayed! */
if ( MYID==myidcounter ) fprintf ( fp," Check accessibility for seismogram file %s... \n",xfile );
sprintf ( errormessage,"PE %i cannot write seismogram file %s!",MYID,xfile );
//if (access(xfile,W_OK|X_OK)==-1) err(errormessage);
if ( ( fpcheck=fopen ( xfile,xmod ) ) ==NULL ) err ( errormessage );
else fclose ( fpcheck );
remove ( xfile );
sprintf ( xfile,"%s_vy.%s.%d",SEIS_FILE,file_ext,myidcounter );
if ( MYID==myidcounter ) fprintf ( fp," Check accessibility for seismogram file %s... \n",xfile );
sprintf ( errormessage,"PE %i cannot write seismogram file %s!",MYID,xfile );
if ( ( fpcheck=fopen ( xfile,xmod ) ) ==NULL ) err ( errormessage );
else fclose ( fpcheck );
remove ( xfile );
for (myidcounter=0; myidcounter< (NPROCX*NPROCY); myidcounter++) {
break;
case 2 : /* pressure only */
sprintf ( xfile,"%s_p.%s.%d",SEIS_FILE,file_ext,myidcounter );
if ( MYID==myidcounter ) fprintf ( fp," Check accessibility for seismogram file %s... \n",xfile );
sprintf ( errormessage,"PE %i cannot write seismogram file %s!",MYID,xfile );
if ( ( fpcheck=fopen ( xfile,xmod ) ) ==NULL ) err ( errormessage );
else fclose ( fpcheck );
remove ( xfile );
switch (SEISMO) {
case 1: /* particle velocities only */
sprintf(xfile,"%s_vx.%s.%d",SEIS_FILE,file_ext,myidcounter);
break;
/*in case of number of PE's=500, there will be 500 messages, too many to be displayed! */
if (MYID==myidcounter) fprintf(fp," Check accessibility for seismogram file %s... \n",xfile);
case 4 : /* everything */
sprintf ( xfile,"%s_vx.%s.%d",SEIS_FILE,file_ext,myidcounter );
/*in case of number of PE's=500, there will be 500 messages, too many to be displayed! */
if ( MYID==myidcounter ) fprintf ( fp," Check accessibility for seismogram file %s... \n",xfile );
sprintf ( errormessage,"PE %i cannot write seismogram file %s!",MYID,xfile );
//if (access(xfile,W_OK|X_OK)==-1) err(errormessage);
if ( ( fpcheck=fopen ( xfile,xmod ) ) ==NULL ) err ( errormessage );
else fclose ( fpcheck );
remove ( xfile );
sprintf ( xfile,"%s_vy.%s.%d",SEIS_FILE,file_ext,myidcounter );
if ( MYID==myidcounter ) fprintf ( fp," Check accessibility for seismogram file %s... \n",xfile );
//if (access(xfile,W_OK|X_OK)==-1) err(errormessage);
sprintf ( errormessage,"PE %i cannot write seismogram file %s!",MYID,xfile );
if ( ( fpcheck=fopen ( xfile,xmod ) ) ==NULL ) err ( errormessage );
else fclose ( fpcheck );
remove ( xfile );
sprintf ( xfile,"%s_p.%s.%d",SEIS_FILE,file_ext,myidcounter );
if ( MYID==myidcounter ) fprintf ( fp," Check accessibility for seismogram file %s... \n",xfile );
sprintf ( errormessage,"PE %i cannot write seismogram file %s!",MYID,xfile );
if ( ( fpcheck=fopen ( xfile,xmod ) ) ==NULL ) err ( errormessage );
else fclose ( fpcheck );
remove ( xfile );
case 3 : /* curl and div only */
sprintf ( xfile,"%s_div.%s.%d",SEIS_FILE,file_ext,myidcounter );
if ( MYID==myidcounter ) fprintf ( fp," Check accessibility for seismogram file %s... \n",xfile );
sprintf ( errormessage,"PE %i cannot write seismogram file %s!",MYID,xfile );
if ( ( fpcheck=fopen ( xfile,xmod ) ) ==NULL ) err ( errormessage );
else fclose ( fpcheck );
remove ( xfile );
sprintf ( xfile,"%s_curl.%s.%d",SEIS_FILE,file_ext,myidcounter );
if ( MYID==myidcounter ) fprintf ( fp," Check accessibility for seismogram file %s... \n",xfile );
sprintf ( errormessage,"PE %i cannot write seismogram file %s!",MYID,xfile );
if ( ( fpcheck=fopen ( xfile,xmod ) ) ==NULL ) err ( errormessage );
else fclose ( fpcheck );
remove ( xfile );
sprintf(errormessage,"PE %i cannot write seismogram file %s!",MYID,xfile);
break;
//if (access(xfile,W_OK|X_OK)==-1) declare_error(errormessage);
if ((fpcheck=fopen(xfile,xmod)) ==NULL) declare_error(errormessage);
else fclose(fpcheck);
remove(xfile);
sprintf(xfile,"%s_vy.%s.%d",SEIS_FILE,file_ext,myidcounter);
if (MYID==myidcounter) fprintf(fp," Check accessibility for seismogram file %s... \n",xfile);
sprintf(errormessage,"PE %i cannot write seismogram file %s!",MYID,xfile);
if ((fpcheck=fopen(xfile,xmod)) ==NULL) declare_error(errormessage);
else fclose(fpcheck);
remove(xfile);
break;
case 2 : /* pressure only */
sprintf(xfile,"%s_p.%s.%d",SEIS_FILE,file_ext,myidcounter);
if (MYID==myidcounter) fprintf(fp," Check accessibility for seismogram file %s... \n",xfile);
sprintf(errormessage,"PE %i cannot write seismogram file %s!",MYID,xfile);
if ((fpcheck=fopen(xfile,xmod)) ==NULL) declare_error(errormessage);
else fclose(fpcheck);
remove(xfile);
break;
case 4 : /* everything */
sprintf(xfile,"%s_vx.%s.%d",SEIS_FILE,file_ext,myidcounter);
/*in case of number of PE's=500, there will be 500 messages, too many to be displayed! */
if (MYID==myidcounter) fprintf(fp," Check accessibility for seismogram file %s... \n",xfile);
sprintf(errormessage,"PE %i cannot write seismogram file %s!",MYID,xfile);
//if (access(xfile,W_OK|X_OK)==-1) declare_error(errormessage);
if ((fpcheck=fopen(xfile,xmod)) ==NULL) declare_error(errormessage);
else fclose(fpcheck);
remove(xfile);
sprintf(xfile,"%s_vy.%s.%d",SEIS_FILE,file_ext,myidcounter);
if (MYID==myidcounter) fprintf(fp," Check accessibility for seismogram file %s... \n",xfile);
//if (access(xfile,W_OK|X_OK)==-1) declare_error(errormessage);
sprintf(errormessage,"PE %i cannot write seismogram file %s!",MYID,xfile);
if ((fpcheck=fopen(xfile,xmod)) ==NULL) declare_error(errormessage);
else fclose(fpcheck);
remove(xfile);
sprintf(xfile,"%s_p.%s.%d",SEIS_FILE,file_ext,myidcounter);
if (MYID==myidcounter) fprintf(fp," Check accessibility for seismogram file %s... \n",xfile);
sprintf(errormessage,"PE %i cannot write seismogram file %s!",MYID,xfile);
if ((fpcheck=fopen(xfile,xmod)) ==NULL) declare_error(errormessage);
else fclose(fpcheck);
remove(xfile);
case 3 : /* curl and div only */
sprintf(xfile,"%s_div.%s.%d",SEIS_FILE,file_ext,myidcounter);
if (MYID==myidcounter) fprintf(fp," Check accessibility for seismogram file %s... \n",xfile);
sprintf(errormessage,"PE %i cannot write seismogram file %s!",MYID,xfile);
if ((fpcheck=fopen(xfile,xmod)) ==NULL) declare_error(errormessage);
else fclose(fpcheck);
remove(xfile);
sprintf(xfile,"%s_curl.%s.%d",SEIS_FILE,file_ext,myidcounter);
if (MYID==myidcounter) fprintf(fp," Check accessibility for seismogram file %s... \n",xfile);
sprintf(errormessage,"PE %i cannot write seismogram file %s!",MYID,xfile);
if ((fpcheck=fopen(xfile,xmod)) ==NULL) declare_error(errormessage);
else fclose(fpcheck);
remove(xfile);
break;
}
}
......@@ -246,24 +319,29 @@ void checkfd ( FILE *fp, float ** prho, float ** ppi, float ** pu,
/*Checking CHECKPOINT Output */
/*-------------------------- */
if ( CHECKPTREAD>0 ) {
strcpy ( xmod,"rb" );
sprintf ( xfile,"%s.%d",CHECKPTFILE,MYID );
fprintf ( fp," Check readability for checkpoint files %s... \n",xfile );
if ( ( ( fpcheck=fopen ( xfile,xmod ) ) ==NULL ) && ( MYID==0 ) ) err ( " PE 0 cannot read checkpoints!" );
else fclose ( fpcheck );
if (CHECKPTREAD>0) {
strcpy(xmod,"rb");
sprintf(xfile,"%s.%d",CHECKPTFILE,MYID);
fprintf(fp," Check readability for checkpoint files %s... \n",xfile);
if (((fpcheck=fopen(xfile,xmod)) ==NULL) && (MYID==0)) declare_error(" PE 0 cannot read checkpoints!");
else fclose(fpcheck);
}
if ( ( CHECKPTWRITE>0 ) ) {
strcpy ( xmod,"ab" );
sprintf ( xfile,"%s.%d",CHECKPTFILE,MYID );
fprintf ( fp," Check writability for checkpoint files %s... \n",xfile );
if ( ( ( fpcheck=fopen ( xfile,xmod ) ) ==NULL ) && ( MYID==0 ) ) err ( " PE 0 cannot write checkpoints!" );
else fclose ( fpcheck ); /* Is there any reason to remove it? */
if ((CHECKPTWRITE>0)) {
strcpy(xmod,"ab");
sprintf(xfile,"%s.%d",CHECKPTFILE,MYID);
fprintf(fp," Check writability for checkpoint files %s... \n",xfile);
if (((fpcheck=fopen(xfile,xmod)) ==NULL) && (MYID==0)) declare_error(" PE 0 cannot write checkpoints!");
else fclose(fpcheck); /* Is there any reason to remove it? */
}
fprintf ( fp," Accessibility of output files from PE %d has been checked successfully.\n", MYID );
fprintf(fp," Accessibility of output files from PE %d has been checked successfully.\n", MYID);
fprintf ( fp,"\n\n --------- DETERMININATION OF MIN AND MAX VELOCITIES -----------\n" );
fprintf(fp,"\n\n --------- DETERMININATION OF MIN AND MAX VELOCITIES -----------\n");
/* low Q frame not yet applied as a absorbing boundary */
/* if (!FREE_SURF) ny1=1+FW;*/
......@@ -277,175 +355,214 @@ void checkfd ( FILE *fp, float ** prho, float ** ppi, float ** pu,
/* find maximum model phase velocity of shear waves at infinite
frequency within the whole model */
if ( L>0 ) { /*viscoelastic*/
for ( i=1+nfw; i<= ( nx-nfw ); i++ ) {
for ( j=ny1; j<= ( ny-nfw ); j++ ) {
c=sqrt ( pu[j][i]* ( 1.0/prho[j][i] ) * ( 1.0+L*ptaus[j][i] ) );
if (L>0) { /*viscoelastic*/
for (i=1+nfw; i<= (nx-nfw); i++) {
for (j=ny1; j<= (ny-nfw); j++) {
c=sqrt(pu[j][i]* (1.0/prho[j][i]) * (1.0+L*ptaus[j][i]));
/*if c is close to zero (water, air), c will be ignored for finding
cmax,cmin*/
if ( ( cmax_s<c ) && ( c>cwater ) ) cmax_s=c;
if ((cmax_s<c) && (c>cwater)) cmax_s=c;
/* find minimum model phase velocity of shear waves at center
frequency of the source */
sum=0.0;
for ( l=1; l<=L; l++ ) {
for (l=1; l<=L; l++) {
ts=DT/peta[l];
sum=sum+ ( ( w*w*ptaus[j][i]*ts*ts ) / ( 1.0+w*w*ts*ts ) );
sum=sum+ ((w*w*ptaus[j][i]*ts*ts) / (1.0+w*w*ts*ts));
}
c=sqrt ( ( pu[j][i]/prho[j][i] ) );
if ( ( cmin_s>c ) && ( c>cwater ) ) cmin_s=c;
c=sqrt((pu[j][i]/prho[j][i]));
if ((cmin_s>c) && (c>cwater)) cmin_s=c;
}
}
} else { /* L=0, elastic */
for ( i=1+nfw; i<= ( nx-nfw ); i++ ) {
for ( j=ny1; j<= ( ny-nfw ); j++ ) {
c=sqrt ( pu[j][i]/prho[j][i] );
for (i=1+nfw; i<= (nx-nfw); i++) {
for (j=ny1; j<= (ny-nfw); j++) {
c=sqrt(pu[j][i]/prho[j][i]);
/*if c is close to zero (water, air), c will be ignored for finding cmax,cmin*/
if ( ( c>cwater ) && ( cmax_s<c ) ) cmax_s=c;
if ( ( c>cwater ) && ( cmin_s>c ) ) cmin_s=c;
if ((c>cwater) && (cmax_s<c)) cmax_s=c;
if ((c>cwater) && (cmin_s>c)) cmin_s=c;
}
}
}
/* find maximum model phase velocity of P-waves at infinite
frequency within the whole model */
if ( L>0 ) { /*viscoelastic*/
for ( i=1+nfw; i<= ( nx-nfw ); i++ ) {
for ( j=ny1; j<= ( ny-nfw ); j++ ) {
c=sqrt ( ppi[j][i]* ( 1.0/prho[j][i] ) * ( 1.0+L*ptaup[j][i] ) );
if ( ( c>cwater ) && ( cmax_p<c ) ) cmax_p=c;
if (L>0) { /*viscoelastic*/
for (i=1+nfw; i<= (nx-nfw); i++) {
for (j=ny1; j<= (ny-nfw); j++) {
c=sqrt(ppi[j][i]* (1.0/prho[j][i]) * (1.0+L*ptaup[j][i]));
if ((c>cwater) && (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++ ) {
for (l=1; l<=L; l++) {
ts=DT/peta[l];
sum=sum+ ( ( w*w*ptaup[j][i]*ts*ts ) / ( 1.0+w*w*ts*ts ) );
sum=sum+ ((w*w*ptaup[j][i]*ts*ts) / (1.0+w*w*ts*ts));
}
c=sqrt ( ( ppi[j][i]/prho[j][i] ) );
if ( ( c>cwater ) && ( cmin_p>c ) ) cmin_p=c;
c=sqrt((ppi[j][i]/prho[j][i]));
if ((c>cwater) && (cmin_p>c)) cmin_p=c;
}
}
} else { /* L=0, elastic */
for ( i=1+nfw; i<= ( nx-nfw ); i++ ) {
for ( j=ny1; j<= ( ny-nfw ); j++ ) {
c=sqrt ( ppi[j][i]/prho[j][i] );
for (i=1+nfw; i<= (nx-nfw); i++) {
for (j=ny1; j<= (ny-nfw); j++) {
c=sqrt(ppi[j][i]/prho[j][i]);
/*if c is close to zero (water, air), c will be ignored for finding cmax,cmin*/
if ( ( c>cwater ) && ( cmax_p<c ) ) cmax_p=c;
if ( ( c>cwater ) && ( cmin_p>c ) ) cmin_p=c;
if ((c>cwater) && (cmax_p<c)) cmax_p=c;
if ((c>cwater) && (cmin_p>c)) cmin_p=c;
}
}
}
fprintf ( fp," Minimum and maximum P-wave and S-wave velocities within subvolumes: \n " );
fprintf ( fp," MYID\t Vp_min(f=fc) \t Vp_max(f=inf) \t Vs_min(f=fc) \t Vs_max(f=inf) \n" );
fprintf ( fp," %d \t %5.2f \t %5.2f \t %5.2f \t %5.2f \n", MYID, cmin_p, cmax_p, cmin_s, cmax_s );
fprintf(fp," Minimum and maximum P-wave and S-wave velocities within subvolumes: \n ");
fprintf(fp," MYID\t Vp_min(f=fc) \t Vp_max(f=inf) \t Vs_min(f=fc) \t Vs_max(f=inf) \n");
fprintf(fp," %d \t %5.2f \t %5.2f \t %5.2f \t %5.2f \n", MYID, cmin_p, cmax_p, cmin_s, cmax_s);
fprintf(fp," Note : if any P- or S-wave velocity is set below 1.0 m/s to simulate water or air,\n");
fprintf(fp," this minimum velocity will be ignored for determining stable DH and DT.\n\n");
fprintf ( fp," Note : if any P- or S-wave velocity is set below 1.0 m/s to simulate water or air,\n" );
fprintf ( fp," this minimum velocity will be ignored for determining stable DH and DT.\n\n" );
if (cmax_s>cmax_p) cmax=cmax_s;
if ( cmax_s>cmax_p ) cmax=cmax_s;
else cmax=cmax_p;
if ( cmin_s<cmin_p ) cmin=cmin_s;
if (cmin_s<cmin_p) cmin=cmin_s;
else cmin=cmin_p;
/* find global maximum for Vp and global minimum for Vs*/
MPI_Allreduce ( &cmax,&cmax_r,1,MPI_FLOAT,MPI_MAX,MPI_COMM_WORLD );
MPI_Allreduce ( &cmin,&cmin_r,1,MPI_FLOAT,MPI_MIN,MPI_COMM_WORLD );
MPI_Allreduce(&cmax,&cmax_r,1,MPI_FLOAT,MPI_MAX,MPI_COMM_WORLD);
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;}
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*temporal );
CFL=cmax*DT/DH;
if ( RSG ) dtstab=DH/cmax;
if ( MYID == 0 ) {
fprintf ( fp," Global values for entire model: \n" );
fprintf ( fp," V_max= %5.2f m/s \t V_min= %5.2f m/s \n\n", cmax,cmin );
fprintf ( fp,"\n\n ------------------ CHECK FOR GRID DISPERSION --------------------\n" );
fprintf ( fp," To satisfactorily limit grid dispersion the number of gridpoints \n" );
fprintf ( fp," per minimum wavelength (of S-waves) should be 6 (better more).\n" );
fprintf ( fp," Here the minimum wavelength is assumed to be minimum model phase velocity \n" );
fprintf ( fp," (of S-waves) at maximum frequency of the source\n" );
fprintf ( fp," devided by maximum frequency of the source.\n" );
fprintf ( fp," Maximum frequency of the source is approximately %6.3f Hz\n",2.0/TS );
fprintf ( fp," The minimum wavelength (P- or S-wave) in the following simulation will\n" );
fprintf ( fp," be %5.3f meter.\n", cmin/fmax );
fprintf ( fp," Thus, the recommended value for DH is %5.3f meter.\n", dhstab );
fprintf ( fp," You have specified DH= %5.3f meter.\n\n", DH );
if ( DH>dhstab )
warning ( " Grid dispersion will influence wave propagation, choose smaller grid spacing (DH)." );
fprintf ( fp," \n\n ----------------------- CHECK FOR STABILITY ---------------------\n" );
fprintf ( fp," The following simulation is stable provided that\n\n" );
if ( RSG ) fprintf ( fp," \t p=cmax*DT/DH < 1,\n\n" );
else fprintf ( fp," \t p=cmax*DT/DH < 1/(sqrt(2)*gamma),\n\n" );
fprintf ( fp," where cmax is the maximum phase velocity at infinite frequency\n" );
if ( RSG==0 ) fprintf ( fp," and gamma = sum(|FD coeff.|)\n" );
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 )
err ( " The simulation will get unstable, choose smaller DT. " );
else fprintf ( fp," The simulation will be stable.\n" );
fprintf ( fp," \n\n --------------------- CHECK FOR INPUT ERRORS ---------------------\n" );
fprintf ( fp," Checking the time of wave propagation. \n" );
if ( SNAP ) {
fprintf ( fp," Checking the snapshot parameters. \n" );
if (TSNAP2>TIME) {
sprintf(errormessage,"\nTSNAP2 = %e (last snapshot) > Time of wave propagation %e. TSNAP2 was changed to be equal to TIME.\n",TSNAP2, TIME);
TSNAP2=TIME;
if (MYID==0)
warning(errormessage); /* if TSNAP2>simulation TIME, snapmerge will generate "additional" snapshots out of nowhere, thus, snapshot files size blow up */
}
snapoutx=NX/ ( float ) IDX;
snapouty=NY/ ( float ) IDY;
fprintf ( fp," Output of snapshot gridpoints per node (NX/NPROCX/IDX) %8.2f .\n", snapoutx );
fprintf ( fp," Output of snapshot gridpoints per node (NY/NPROCY/IDY) %8.2f .\n", snapouty );
if ( snapoutx- ( int ) snapoutx>0 )
err ( "\n\n Ratio NX-NPROCX-IDX must be whole-numbered \n\n" );
if ( snapouty- ( int ) snapouty>0 )
err ( "\n\n Ratio NY-NPROCY-IDY must be whole-numbered \n\n" );
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*temporal);
CFL=cmax*DT/DH;