Commit 917d9b0b authored by Johannes Spazier's avatar Johannes Spazier

This is an alpha version that dynamically creates a subset of the input grid...

This is an alpha version that dynamically creates a subset of the input grid according to the source location.
The number of steps that are necessary to run a complete simulation is hardcoded for now. Maybe there is a better
approach.

This feature is needed to use EasyWave from the web where a user can specify any source location and we have to
choose the right grid to answer his request.
parent 08d748cd
......@@ -64,10 +64,6 @@ int main( int argc, char **argv )
// Read points of interest
ierr = ewLoadPOIs(); if(ierr) return ierr;
// Init tsunami with faults or uplift-grid
ierr = ewSource(); if(ierr) return ierr;
Log.print( "Read source from %s", Par.fileSource );
// Write model parameters into the log
ewLogParams();
......
......@@ -70,6 +70,9 @@ extern int Jmax;
#define getLon(i) (LonMin+(i-1)*DLon)
#define getLat(j) (LatMin+(j-1)*DLat)
#define getI(lon) round( ( (lon-LonMin) / DLon) )
#define getJ(lat) round( ( (lat-LatMin) / DLat) )
int ewLoadBathymetry();
int ewParam( int argc, char **argv );
void ewLogParams(void);
......@@ -78,6 +81,8 @@ int ewSource();
int ewStep();
int ewStepCor();
int setWaveHeights();
int ewStart2DOutput();
int ewOut2D();
int ewDump2D();
......
......@@ -56,7 +56,6 @@ int ewLoadBathymetry()
// try to allocate memory for GRIDNODE structure and for caching arrays
printf_v("Size: %d %d %luMB\n", NLon, NLat, sizeof(float)*MAX_VARS_PER_NODE*NLon*NLat/1024/1024);
if( Node.mallocMem() ) return Err.post( Err.msgAllocateMem() );
if( isBin ) {
ierr = fread( &LonMin, sizeof(double), 1, fp ); ierr = fread( &LonMax, sizeof(double), 1, fp );
......@@ -74,12 +73,72 @@ int ewLoadBathymetry()
Dx = Re * g2r( DLon ); // in m along the equator
Dy = Re * g2r( DLat );
// Init tsunami with faults or uplift-grid
ierr = ewSource(); if(ierr) return ierr;
Log.print( "Read source from %s", Par.fileSource );
int rgL, rgR, rgT, rgB;
/* TODO: How do we know how many steps are necessary? Assume one grid extension each 5 steps for now. */
int numIter = ceil( Par.timeMax / Par.dt / 5.0 );
rgL = My_max( 0, Imin - numIter );
rgR = My_min( Imax + numIter, NLon - 1 );
rgT = My_max( 0, Jmin - numIter );
rgB = My_min( Jmax + numIter, NLat - 1 );
double LonMinSec = getLon( rgL + 1 );
double LonMaxSec = getLon( rgR + 1 );
double LatMinSec = getLat( rgT + 1 );
double LatMaxSec = getLat( rgB + 1 );
/* Use the following lines for a static range that equals "gebco_08_atlantic_2.grd". */
/*double LonMinSec = -100.004167;
double LonMaxSec = 19.995833;
double LatMinSec = -89.995833;
double LatMaxSec = 89.987500;
rgL = getI( LonMinSec );
rgR = getI( LonMaxSec );
rgT = getJ( LatMinSec );
rgB = getJ( LatMaxSec );*/
LonMin = LonMinSec;
LonMax = LonMaxSec;
LatMin = LatMinSec;
LatMax = LatMaxSec;
printf_v("rgL, rgR, rgT, rgB: %u %u %u %u\n", rgL, rgR, rgT, rgB );
/* upon here NLot and NLat are set according to the desired section */
int NLonBase = NLon;
int NLatBase = NLat;
NLon = rgR - rgL + 1;
NLat = rgB - rgT + 1;
printf_v("NLon, NLat: %u %u\n", NLon, NLat );
if( Node.mallocMem() ) return Err.post( Err.msgAllocateMem() );
/* NOTE: optimal would be reading everything in one step, but that does not work because rows and columns are transposed
* (only possible with binary data at all) - use temporary buffer for now (consumes additional memory!) */
float *buf = new float[ NLat*NLon ];
ierr = fread( buf, sizeof(float), NLat*NLon, fp );
//ierr = fread( buf, sizeof(float), NLat*NLon, fp );
fseek( fp, rgT * NLonBase * sizeof(float), SEEK_CUR );
for( j = 1; j <= NLat; j++ ) {
fseek ( fp, rgL * sizeof(float), SEEK_CUR );
ierr = fread( buf + (j-1) * NLon, sizeof(float), NLon, fp );
fseek ( fp, ( (NLonBase - 1) - rgR) * sizeof(float), SEEK_CUR );
}
double zmin = 0.0;
double zmax = 0.0;
for( i=1; i<=NLon; i++ ) {
for( j=1; j<=NLat; j++ ) {
......@@ -91,6 +150,9 @@ int ewLoadBathymetry()
else
ierr = fscanf( fp, " %f ", &fval );
if( fval > zmax ) zmax = fval;
if( fval < zmin ) zmin = fval;
Node(m, iTopo) = fval;
Node(m, iTime) = -1;
Node(m, iD) = -fval;
......@@ -104,6 +166,29 @@ int ewLoadBathymetry()
}
}
FILE *fout = fopen( "range.grd", "wb");
fwrite( "DSBB", 1, 4, fout );
shval = NLon;
fwrite( &shval, sizeof(unsigned short), 1, fout );
shval = NLat;
fwrite( &shval, sizeof(unsigned short), 1, fout );
fwrite( &LonMinSec, sizeof(double), 1, fout );
fwrite( &LonMaxSec, sizeof(double), 1, fout );
fwrite( &LatMinSec, sizeof(double), 1, fout );
fwrite( &LatMaxSec, sizeof(double), 1, fout );
fwrite( &zmin, sizeof(double), 1, fout );
fwrite( &zmax, sizeof(double), 1, fout );
for( j=0; j<NLat; j++ ) {
fwrite( buf + j * NLon, sizeof(float), NLon, fout );
}
fclose( fout );
delete[] buf;
for( k=1; k<MAX_VARS_PER_NODE-2; k++ ) {
......@@ -199,5 +284,7 @@ int ewLoadBathymetry()
if( Node(idx(j,NLon), iD) != 0 ) C4[j] = 1./sqrt(Gravity*Node(idx(j,NLon), iD));
}
setWaveHeights();
return 0;
}
......@@ -10,6 +10,8 @@ int Imax;
int Jmin;
int Jmax;
cOgrd uZ;
#define SRC_GRD 1
#define SRC_FLT 2
......@@ -20,10 +22,8 @@ int ewSource()
int i,j,ierr,srcType;
double lon,lat,dz,absuzmax,absuzmin;
FILE *fp;
cOkadaEarthquake eq;
cOgrd uZ;
CNode& Node = *gNode;
cOkadaEarthquake eq;
// check input file type: GRD or fault
if( (fp = fopen( Par.fileSource, "rb" )) == NULL ) return Err.post( Err.msgOpenFile(Par.fileSource) );
......@@ -160,11 +160,7 @@ int ewSource()
lon = getLon(i);
lat = getLat(j);
if( Node(idx(j,i), iD) != 0. )
dz = Node(idx(j,i), iH) = uZ.getVal( lon,lat );
else
dz = Node(idx(j,i), iH) = 0.;
dz = uZ.getVal( lon,lat );
if( fabs(dz) > Par.sshClipThreshold ) {
Imin = My_min( Imin, i );
......@@ -183,5 +179,49 @@ int ewSource()
Jmin = My_max( Jmin - 2, 2 );
Jmax = My_min( Jmax + 2, NLat-1 );
printf_v("Imin, Imax, Jmin, Jmax: %u %u %u %u\n", Imin, Imax, Jmin, Jmax);
return 0;
}
int setWaveHeights() {
int i, j;
double lon, lat, dz;
CNode& Node = *gNode;
Imin = NLon; Imax = 1; Jmin = NLat; Jmax = 1;
for( i=1; i<=NLon; i++ ) {
for( j=1; j<=NLat; j++ ) {
lon = getLon(i);
lat = getLat(j);
//lon = (LonMinSec+(i-1)*DLon);
//lat = (LatMinSec+(j-1)*DLat);
if( Node(idx(j,i), iD) != 0. )
dz = Node(idx(j,i), iH) = uZ.getVal( lon,lat );
else
dz = Node(idx(j,i), iH) = 0.;
if( fabs(dz) > Par.sshClipThreshold ) {
Imin = My_min( Imin, i );
Imax = My_max( Imax, i );
Jmin = My_min( Jmin, j );
Jmax = My_max( Jmax, j );
}
}
}
Imin = My_max( Imin - 2, 2 );
Imax = My_min( Imax + 2, NLon-1 );
Jmin = My_max( Jmin - 2, 2 );
Jmax = My_min( Jmax + 2, NLat-1 );
printf_v("Imin, Imax, Jmin, Jmax: %u %u %u %u\n", Imin, Imax, Jmin, Jmax);
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment