From 917d9b0bb8afc13af8e1b2978064a0a4877b3ae8 Mon Sep 17 00:00:00 2001 From: Johannes Spazier Date: Thu, 21 Nov 2013 13:40:26 +0000 Subject: [PATCH] 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. --- cpu/branches/web/EasyWave.cu | 4 -- cpu/branches/web/easywave.h | 5 ++ cpu/branches/web/ewGrid.cpp | 91 ++++++++++++++++++++++++++++++++++- cpu/branches/web/ewSource.cpp | 56 ++++++++++++++++++--- 4 files changed, 142 insertions(+), 14 deletions(-) diff --git a/cpu/branches/web/EasyWave.cu b/cpu/branches/web/EasyWave.cu index 81c6a04..a7c807c 100644 --- a/cpu/branches/web/EasyWave.cu +++ b/cpu/branches/web/EasyWave.cu @@ -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(); diff --git a/cpu/branches/web/easywave.h b/cpu/branches/web/easywave.h index 0354426..d293889 100644 --- a/cpu/branches/web/easywave.h +++ b/cpu/branches/web/easywave.h @@ -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(); diff --git a/cpu/branches/web/ewGrid.cpp b/cpu/branches/web/ewGrid.cpp index fbf705a..8daa857 100644 --- a/cpu/branches/web/ewGrid.cpp +++ b/cpu/branches/web/ewGrid.cpp @@ -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 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; +} -- GitLab