Added a version of EasyWave that uses multiple physical GPUs to run the simulation.

cpu/branches/multi-gpu/EasyWave.cu

+#define HEADER "\neasyWave ver.2013-04-11\n"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+#include <utilits.h>
+#include "easywave.h"
+
+#include "ewGpuNode.cuh"
+
+CNode *gNode;
+
+double diff(timespec start, timespec end) {
+
+	timespec temp;
+	if ((end.tv_nsec-start.tv_nsec)<0) {
+		temp.tv_sec = end.tv_sec-start.tv_sec-1;
+		temp.tv_nsec = 1000000000+end.tv_nsec-start.tv_nsec;
+	} else {
+		temp.tv_sec = end.tv_sec-start.tv_sec;
+		temp.tv_nsec = end.tv_nsec-start.tv_nsec;
+	}
+
+	return (double)((double)temp.tv_nsec / 1000000000.0 + (double)temp.tv_sec);
+}
+
+int commandLineHelp( void );
+
+int main( int argc, char **argv )
+{
+  char buf[1024];
+  int ierr,argn,elapsed;
+  int lastProgress,lastPropagation,lastDump;
+  int loop;
+
+  printf(HEADER);
+  Err.setchannel(MSG_OUTFILE);
+
+  // Read parameters from command line and use default
+  ierr = ewParam( argc, argv ); if(ierr) return commandLineHelp();
+
+  // Log command line
+  /* FIXME: buffer overflow */
+  sprintf( buf, "Command line: " );
+  for( argn=1; argn<argc; argn++ ) {
+    strcat( buf, " " );
+    strcat( buf, argv[argn] );
+  }
+  Log.print( "%s", buf );
+
+  if( Par.gpu ) {
+	  gNode = new CGpuNode();
+  } else {
+	  gNode = new CStructNode();
+	  //gNode = new CArrayNode();
+  }
+
+  CNode& Node = *gNode;
+
+  // Read bathymetry
+  ierr = ewLoadBathymetry(); if(ierr) return ierr;
+
+  // 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();
+
+  if( Par.outPropagation ) ewStart2DOutput();
+
+  Node.copyToGPU();
+
+  // Main loop
+  Log.print("Starting main loop...");
+  printf("Starting main loop %ld msec\n", clock()/(CLOCKS_PER_SEC/1000) );
+
+  timespec start, end;
+  clock_gettime(CLOCK_MONOTONIC, &start);
+
+  for( Par.time=0,loop=1,lastProgress=Par.outProgress,lastPropagation=Par.outPropagation,lastDump=0;
+    Par.time<=Par.timeMax; loop++,Par.time+=Par.dt,lastProgress+=Par.dt,lastPropagation+=Par.dt ) {
+
+	/* FIXME: check if Par.poiDt can be used for those purposes */
+    if( Par.filePOIs && Par.poiDt && ((Par.time/Par.poiDt)*Par.poiDt == Par.time) ) {
+    	Node.copyIntermediate();
+    	ewSavePOIs();
+    }
+
+    Node.run();
+
+    elapsed = ((int)clock())/CLOCKS_PER_SEC;
+
+    if( Par.outProgress ) {
+      if( lastProgress >= Par.outProgress ) {
+        printf( "Model time = %s,   elapsed: %ld msec\n", utlTimeSplitString(Par.time), clock()/(CLOCKS_PER_SEC/1000) );
+        Log.print( "Model time = %s,   elapsed: %ld msec", utlTimeSplitString(Par.time), clock()/(CLOCKS_PER_SEC/1000) );
+        lastProgress = 0;
+      }
+    }
+
+    if( Par.outPropagation ) {
+      if( lastPropagation >= Par.outPropagation ) {
+    	Node.copyIntermediate();
+        ewOut2D();
+        lastPropagation = 0;
+      }
+    }
+
+    if( Par.outDump ) {
+      if( (elapsed-lastDump) >= Par.outDump ) {
+    	Node.copyIntermediate();
+        ewDumpPOIs();
+        ewDump2D();
+        lastDump = elapsed;
+      }
+    }
+
+  } // main loop
+  clock_gettime(CLOCK_MONOTONIC, &end);
+  Log.print("Finishing main loop");
+
+  /* TODO: check if theses calls can be combined */
+  Node.copyIntermediate();
+  Node.copyFromGPU();
+
+  // Final output
+  Log.print("Final dump...");
+  ewDumpPOIs();
+  ewDump2D();
+
+  Node.freeMem();
+
+  printf_v("Runtime: %.3lf\n", diff(start, end) * 1000.0);
+
+  delete gNode;
+
+  return 0;
+}
+
+
+//========================================================================
+int commandLineHelp( void )
+{
+  printf( "Usage: easyWave  -grid ...  -source ...  -time ... [optional parameters]\n" );
+  printf( "-grid ...         bathymetry in GoldenSoftware(C) GRD format (text or binary)\n" );
+  printf( "-source ...       input wave either als GRD-file or file with Okada faults\n" );
+  printf( "-time ...         simulation time in [min]\n" );
+  printf( "Optional parameters:\n" );
+  printf( "-step ...         simulation time step, default- estimated from bathymetry\n" );
+  printf( "-coriolis         use Coriolis fource, default- no\n" );
+  printf( "-poi ...          POIs file\n" );
+  printf( "-label ...        model name, default- 'eWave'\n" );
+  printf( "-progress ...     show simulation progress each ... minutes, default- 10\n" );
+  printf( "-propagation ...  write wave propagation grid each ... minutes, default- 5\n" );
+  printf( "-dump ...         make solution dump each ... physical seconds, default- 0\n" );
+  printf( "-nolog            deactivate logging\n" );
+  printf( "-poi_dt_out ...   output time step for mariograms in [sec], default- 30\n" );
+  printf( "-poi_search_dist ...  in [km], default- 10\n" );
+  printf( "-poi_min_depth ...    in [m], default- 1\n" );
+  printf( "-poi_max_depth ...    in [m], default- 10 000\n" );
+  printf( "-poi_report       enable POIs loading report, default- disabled\n" );
+  printf( "-ssh0_rel ...     relative threshold for initial wave, default- 0.01\n" );
+  printf( "-ssh0_abs ...     absolute threshold for initial wave in [m], default- 0\n" );
+  printf( "-ssh_arrival ...  threshold for arrival times in [m], default- 0.001\n" );
+  printf( "                  negative value considered as relative threshold\n" );
+  printf( "-gpu              start GPU version of EasyWave (requires a CUDA capable device)\n" );
+  printf( "-verbose          generate verbose output on stdout\n" );
+  printf( "\nExample:\n" );
+  printf( "\t easyWave -grid gebcoIndonesia.grd  -source fault.inp  -time 120\n\n" );
+
+  return -1;
+}

cpu/branches/multi-gpu/Makefile

+-include ../../../Makefile.inc
+
+CC=g++
+NVCC=nvcc
+CFLAGS=-O3 -I.
+
+ARCH?=compute_20
+CODE?=sm_21
+NVFLAGS=$(CFLAGS) -gencode arch=$(ARCH),code=$(CODE)
+
+CPPS=$(wildcard *.cpp)
+CUS=$(wildcard *.cu)
+
+CU_OBJS=$(patsubst %.cu,%.o,$(CUS) )
+
+OBJECTS=$(patsubst %.cpp, %.o, $(CPPS) ) $(patsubst %.cu,%.o,$(CUS) )
+
+all: EasyWave
+
+EasyWave: $(OBJECTS) link.o
+	$(NVCC) -o $@ $^
+
+%.o: %.cpp *.h
+	$(CC) -c $(CFLAGS) -o $@ $<
+	
+%.o: %.cu *.cuh *.h
+	$(NVCC) -dc $(NVFLAGS) -x cu -o $@ $<
+
+link.o: $(CU_OBJS)
+	$(NVCC) -dlink $(NVFLAGS) -o $@ $^
+
+clean:
+	rm -f EasyWave *.o

cpu/branches/multi-gpu/Makefile_debug

+CC=g++
+NVCC=nvcc
+CFLAGS=-O3 -I. -g
+
+CODE=sm_35
+ARCH=compute_35
+NVFLAGS=$(CFLAGS) -gencode arch=$(ARCH),code=$(CODE)
+
+CPPS=$(wildcard *.cpp)
+CUS=$(wildcard *.cu)
+
+CU_OBJS=$(patsubst %.cu,%.o,$(CUS) )
+
+OBJECTS=$(patsubst %.cpp, %.o, $(CPPS) ) $(patsubst %.cu,%.o,$(CUS) )
+
+all: EasyWave
+
+EasyWave: $(OBJECTS) link.o
+	$(NVCC) -g -G -o $@ $^
+
+%.o: %.cpp *.h
+	$(CC) -c $(CFLAGS) -o $@ $<
+	
+%.o: %.cu *.cuh *.h
+	$(NVCC) -G -dc $(NVFLAGS) -x cu -o $@ $<
+
+link.o: $(CU_OBJS)
+	$(NVCC) -G -dlink $(NVFLAGS) -o $@ $^
+
+clean:
+	rm -f EasyWave *.o

cpu/branches/multi-gpu/cOgrd.cpp

+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <utilits.h>
+#include <cOgrd.h>
+
+#define DEFAULT_NOVAL 0.0
+#define ROWWISE 1
+#define COLWISE 2
+#define TOP2BOT -1
+#define BOT2TOP 1
+
+
+int iabs( int i )
+{
+  return( (i<0) ? (-i) : i );
+}
+
+
+//=========================================================================
+// Zero Constructor
+cOgrd::cOgrd( )
+{
+  xmin = xmax = dx = 0.;
+  nx = 0;
+  ymin = ymax = dy = 0.;
+  ny = 0;
+  nnod = 0;
+  noval = DEFAULT_NOVAL;
+  val = NULL;
+}
+
+
+//=========================================================================
+// Copy constructor similar to operator =
+cOgrd::cOgrd( const cOgrd& grd )
+{
+  xmin=grd.xmin; xmax=grd.xmax; dx=grd.dx; nx=grd.nx;
+  ymin=grd.ymin; ymax=grd.ymax; dy=grd.dy; ny=grd.ny;
+  nnod=grd.nnod;
+  noval=grd.noval;
+
+  val = new double[nnod];
+  memcpy( val, grd.val, nnod*sizeof(double) );
+}
+
+
+//=========================================================================
+cOgrd::cOgrd( double xmin0, double xmax0, int nx0, double ymin0, double ymax0, int ny0 )
+{
+  noval = DEFAULT_NOVAL;
+  val = NULL;
+  initialize( xmin0, xmax0, nx0, ymin0, ymax0, ny0 );
+}
+
+
+//=========================================================================
+cOgrd::cOgrd( double xmin0, double xmax0, double dx0, double ymin0, double ymax0, double dy0 )
+{
+  noval = DEFAULT_NOVAL;
+  val = NULL;
+  initialize( xmin0, xmax0, dx0, ymin0, ymax0, dy0 );
+}
+
+
+//=========================================================================
+// Destructor
+cOgrd::~cOgrd()
+{
+  if(val) delete [] val;
+}
+
+
+//=========================================================================
+// Initialize with number of nodes
+int cOgrd::initialize( double xmin0, double xmax0, int nx0, double ymin0, double ymax0, int ny0 )
+{
+
+  xmin = xmin0;
+  xmax = xmax0;
+  nx = nx0;
+  dx = (xmax-xmin)/(nx-1);
+
+  ymin = ymin0;
+  ymax = ymax0;
+  ny = ny0;
+  dy = (ymax-ymin)/(ny-1);
+
+  nnod = nx*ny;
+
+  if( val ) delete [] val;
+
+  val = new double[nnod];
+  for( int l=0; l<nnod; l++ ) val[l]=noval;
+
+  return 0;
+}
+
+
+//=========================================================================
+// Initialize with grid step
+int cOgrd::initialize( double xmin0, double xmax0, double dx0, double ymin0, double ymax0, double dy0 )
+{
+
+  xmin = xmin0;
+  xmax = xmax0;
+  dx = dx0;
+  nx = (int)((xmax-xmin)/dx + 1.e-10) + 1;
+
+  ymin = ymin0;
+  ymax = ymax0;
+  dy = dy0;
+  ny = (int)((ymax-ymin)/dy + 1.e-10) + 1;
+
+  nnod = nx*ny;
+
+  if( val ) delete [] val;
+
+  val = new double[nnod];
+  for( int l=0; l<nnod; l++ ) val[l]=noval;
+
+  return 0;
+}
+
+
+//=========================================================================
+void cOgrd::setNoval( double val )
+{
+  noval = val;
+}
+
+
+//=========================================================================
+double cOgrd::getNoval()
+{
+  return noval;
+}
+
+
+//=========================================================================
+// only read header from GRD-file
+int cOgrd::readHeader( const char *grdfile )
+{
+  FILE *fp;
+  char dsaa_label[8];
+  int ierr,isBin;
+  short shval;
+  float fval;
+  double dval;
+
+  if( (fp = fopen( grdfile, "rb" )) == NULL ) return Err.post(Err.msgOpenFile(grdfile));
+  memset( dsaa_label, 0, 5 );
+  ierr = fread( dsaa_label, 4, 1, fp );
+  if( !strcmp( dsaa_label,"DSAA" ) )
+    isBin = 0;
+  else if( !strcmp( dsaa_label,"DSBB" ) )
+    isBin = 1;
+  else {
+    fclose(fp);
+    return Err.post(Err.msgReadFile( grdfile, 1, "not a GRD-file" ));
+  }
+  fclose(fp);
+
+  if( isBin ) {
+    fp = fopen( grdfile, "rb" );
+    ierr = fread( dsaa_label, 4, 1, fp );
+    ierr = fread( &shval, sizeof(short), 1, fp ); nx = shval;
+    ierr = fread( &shval, sizeof(short), 1, fp ); ny = shval;
+    ierr = fread( &xmin, sizeof(double), 1, fp ); ierr = fread( &xmax, sizeof(double), 1, fp );
+    ierr = fread( &ymin, sizeof(double), 1, fp ); ierr = fread( &ymax, sizeof(double), 1, fp );
+    ierr = fread( &dval, sizeof(double), 1, fp ); ierr = fread( &dval, sizeof(double), 1, fp ); // zmin zmax
+  }
+  else {
+    fp = fopen( grdfile, "rt" );
+    ierr = fscanf( fp, "%s", dsaa_label );
+    ierr = fscanf( fp, " %d %d ", &nx, &ny );
+    ierr = fscanf( fp, " %lf %lf ", &xmin, &xmax );
+    ierr = fscanf( fp, " %lf %lf ", &ymin, &ymax );
+    ierr = fscanf( fp, " %*s %*s " );   // zmin, zmax
+  }
+
+  fclose( fp );
+
+  nnod = nx*ny;
+  dx = (xmax-xmin)/(nx-1);
+  dy = (ymax-ymin)/(ny-1);
+
+  return 0;
+}
+
+
+//=========================================================================
+// read shape from GRD-file, reset values to zero
+int cOgrd::readShape( const char *grdfile )
+{
+
+  int ierr = readGRD( grdfile ); if(ierr) return ierr;
+
+  for( int l=0; l<nnod; l++ ) val[l]=noval;
+
+  return 0;
+}
+
+
+//=========================================================================
+// Grid initialization from Golden Software GRD-file
+int cOgrd::readGRD( const char *grdfile )
+{
+  FILE *fp;
+  char dsaa_label[8];
+  int i,j,ierr,isBin;
+  short shval;
+  float fval;
+  double dval;
+
+
+  // check if bathymetry file is in ascii or binary format
+  if( (fp = fopen( grdfile, "rb" )) == NULL ) return Err.post(Err.msgOpenFile(grdfile));
+  memset( dsaa_label, 0, 5 );
+  ierr = fread( dsaa_label, 4, 1, fp );
+  if( !strcmp( dsaa_label,"DSAA" ) )
+    isBin = 0;
+  else if( !strcmp( dsaa_label,"DSBB" ) )
+    isBin = 1;
+  else {
+    fclose(fp);
+    return Err.post(Err.msgReadFile( grdfile, 1, "not GRD-file" ));
+  }
+  fclose(fp);
+
+
+  // Read Surfer GRD-file
+  if( isBin ) {
+    fp = fopen( grdfile, "rb" );
+    ierr = fread( dsaa_label, 4, 1, fp );
+    ierr = fread( &shval, sizeof(short), 1, fp ); nx = shval;
+    ierr = fread( &shval, sizeof(short), 1, fp ); ny = shval;
+  }
+  else {
+    fp = fopen( grdfile, "rt" );
+    ierr = fscanf( fp, "%s", dsaa_label );
+    ierr = fscanf( fp, " %d %d ", &nx, &ny );
+  }
+
+  nnod = nx*ny;
+
+  if( val ) delete [] val;
+  val = new double[nnod];
+  for( int l=0; l<nnod; l++ ) val[l]=noval;
+
+  if( isBin ) {
+    ierr = fread( &xmin, sizeof(double), 1, fp ); ierr = fread( &xmax, sizeof(double), 1, fp );
+    ierr = fread( &ymin, sizeof(double), 1, fp ); ierr = fread( &ymax, sizeof(double), 1, fp );
+    ierr = fread( &dval, sizeof(double), 1, fp ); ierr = fread( &dval, sizeof(double), 1, fp ); // zmin zmax
+  }
+  else {
+    ierr = fscanf( fp, " %lf %lf ", &xmin, &xmax );
+    ierr = fscanf( fp, " %lf %lf ", &ymin, &ymax );
+    ierr = fscanf( fp, " %*s %*s " );   // zmin, zmax
+  }
+
+  dx = (xmax-xmin)/(nx-1);
+  dy = (ymax-ymin)/(ny-1);
+
+  for( j=0; j<ny; j++ ) {
+    for( i=0; i<nx; i++ ) {
+
+      if( isBin )
+        ierr = fread( &fval, sizeof(float), 1, fp );
+      else
+        ierr = fscanf( fp, " %f ", &fval );
+
+      val[idx(i,j)] = (double)fval;
+    }
+  }
+
+  fclose( fp );
+
+  return 0;
+}
+
+
+
+//=========================================================================
+// Grid initialization from XYZ-file
+int cOgrd::readXYZ( const char *fname )
+{
+  #define MAXRECLEN 254
+  FILE *fp;
+  char record[254];
+  int i,j,line;
+  int ordering,ydirection;
+  int nxny;
+  double x0,x,y0,y,z;
+
+
+  // Open plain XYZ-file
+  if( (fp = fopen( fname, "rt" )) == NULL )
+    return utlPostError( utlErrMsgOpenFile(fname) );
+
+  // Check xyz-file format, define number of nodes and min-max
+  rewind( fp );
+  for( line=0,nxny=0,xmin=ymin=RealMax,xmax=ymax=-RealMax; utlReadNextDataRecord( fp, record, &line ) != EOF; ) {
+    if( sscanf( record, "%lf %lf %lf", &x, &y, &z ) != 3 ) return Err.post(Err.msgReadFile( fname, line, "X Y Z" ));
+    nxny++;
+    if( x < xmin ) xmin = x;
+    if( x > xmax ) xmax = x;
+    if( y < ymin ) ymin = y;
+    if( y > ymax ) ymax = y;
+  }
+
+  // Read first two lines and define if file is row- or column-wise
+  rewind( fp );
+  line = 0;
+  utlReadNextDataRecord( fp, record, &line );
+  sscanf( record, "%lf %lf %lf", &x0, &y0, &z );
+  utlReadNextDataRecord( fp, record, &line );
+  sscanf( record, "%lf %lf %lf", &x, &y, &z );
+  if( x0==x && y0!=y )
+    ordering = COLWISE;
+  else if( x0!=x && y0==y )
+    ordering = ROWWISE;
+  else
+    return Err.post(Err.msgReadFile( fname, line, "Cannot recognise data ordering" ));
+
+  // Define nx and ny
+  rewind( fp );
+  line = 0;
+  utlReadNextDataRecord( fp, record, &line );
+  sscanf( record, "%lf %lf %lf", &x0, &y0, &z );
+  if( ordering == ROWWISE ) {
+    nx = 1;
+    while( utlReadNextDataRecord( fp, record, &line ) != EOF ) {
+      sscanf( record, "%lf %lf %lf", &x, &y, &z );
+      if( y != y0 ) break;
+      nx++;
+    }
+    ny = nxny / nx;
+    if( y > y0 )
+      ydirection = BOT2TOP;
+    else
+      ydirection = TOP2BOT;
+
+  }
+  else if( ordering == COLWISE ) {
+    ny = 1;
+    while( utlReadNextDataRecord( fp, record, &line ) != EOF ) {
+      sscanf( record, "%lf %lf %lf", &x, &y, &z );
+      if( x != x0 ) break;
+      if( ny == 1 ) {
+        if( y > y0 )
+          ydirection = BOT2TOP;
+        else
+          ydirection = TOP2BOT;