Creating a new branch to implement some CPU related optimizations.

code/branches/cpu/EasyWave.cu

+/*
+ * EasyWave - A realtime tsunami simulation program with GPU support.
+ * Copyright (C) 2014  Andrey Babeyko, Johannes Spazier
+ * GFZ German Research Centre for Geosciences (http://www.gfz-potsdam.de)
+ *
+ * Parts of this program (especially the GPU extension) were developed
+ * within the context of the following publicly funded project:
+ * - TRIDEC, EU 7th Framework Programme, Grant Agreement 258723
+ *   (http://www.tridec-online.eu)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Affero General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * This program 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 Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#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.copyPOIs();
+    	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;
+      }
+    }
+
+    fflush(stdout); 
+
+    if( Par.outDump ) {
+      if( (elapsed-lastDump) >= Par.outDump ) {
+    	Node.copyIntermediate();
+    	/* FIXME: needs tArr as well */
+        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;
+}

code/branches/cpu/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

code/branches/cpu/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

code/branches/cpu/cOgrd.h

+/*
+ * EasyWave - A realtime tsunami simulation program with GPU support.
+ * Copyright (C) 2014  Andrey Babeyko, Johannes Spazier
+ * GFZ German Research Centre for Geosciences (http://www.gfz-potsdam.de)
+ *
+ * Parts of this program (especially the GPU extension) were developed
+ * within the context of the following publicly funded project:
+ * - TRIDEC, EU 7th Framework Programme, Grant Agreement 258723
+ *   (http://www.tridec-online.eu)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Affero General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * This program 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 Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef OGRD_H
+#define OGRD_H
+
+class cOgrd
+{
+
+protected:
+
+
+public:
+  int nx,ny;
+  int nnod;
+  double noval;
+  double xmin,xmax;
+  double ymin,ymax;
+  double dx,dy;
+  double *val;
+
+  cOgrd();
+  cOgrd( const cOgrd& );
+  cOgrd( double xmin0, double xmax0, int nx0, double ymin0, double ymax0, int ny0 );
+  cOgrd( double xmin0, double xmax0, double dx0, double ymin0, double ymax0, double dy0 );
+  ~cOgrd();
+
+  void setNoval( double val );
+  double getNoval();
+  int initialize( double xmin0, double xmax0, int nx0, double ymin0, double ymax0, int ny0 );
+  int initialize( double xmin0, double xmax0, double dx0, double ymin0, double ymax0, double dy0 );
+  int readShape( const char *grdfile );
+  int readHeader( const char *grdfile );
+  int readGRD( const char *grdfile );
+  int readXYZ( const char *xyzfile );
+  int readRasterStream( FILE *fp, int ordering, int ydirection );
+  cOgrd* extract( int i1, int i2, int j1, int j2 );
+  cOgrd* extract( double x1, double x2, double y1, double y2 );
+  int idx( int i, int j );
+  double& operator() ( int i, int j );
+  double& operator() ( int idx );
+  cOgrd& operator= ( const cOgrd& grd );
+  cOgrd& operator*= ( double coeff );
+  cOgrd& operator+= ( cOgrd& grd );
+  void getIJ( int idx, int& i, int& j );
+  int getIJ( double x, double y, int& i, int& j );
+  double getX( int i, int j );
+  double getX( int idx );
+  double getY( int i, int j );
+  double getY( int idx );
+  double getVal( int idx );
+  double getVal( int i, int j );
+  double getVal( double x, double y );
+  double getMaxVal();
+  double getMaxVal( int& i, int& j );
+  double getMinVal();
+  double getMinVal( int& i, int& j );
+  double getMaxAbsVal();
+  double getMaxAbsVal( int& i, int& j );
+  double getMaxAbsValBnd();
+  void setVal( double value, int i, int j );
+  void setVal( double value, int idx );
+  void reset();
+  void resetVal();
+  int getIntersectionRegion( const cOgrd &grd, int& imin, int& imax, int& jmin, int& jmax );
+  int interpolateFrom( cOgrd &grd, int resetValues );
+  int getNearestIdx( double x, double y );
+  int getNearestIdx( double x, double y, double rangemin, double rangemax );
+  int equalTo( cOgrd& grd );
+  int isSameShapeTo( cOgrd& grd );
+  void smooth( int radius );
+  int writeGRD( const char *fname );
+  int writeGRDbin( const char *fname );
+  int writeXYZ( const char *fname );
+};
+
+
+#endif  // OGRD_H

code/branches/cpu/cOkadaEarthquake.cpp

+/*
+ * EasyWave - A realtime tsunami simulation program with GPU support.
+ * Copyright (C) 2014  Andrey Babeyko, Johannes Spazier
+ * GFZ German Research Centre for Geosciences (http://www.gfz-potsdam.de)
+ *
+ * Parts of this program (especially the GPU extension) were developed
+ * within the context of the following publicly funded project:
+ * - TRIDEC, EU 7th Framework Programme, Grant Agreement 258723
+ *   (http://www.tridec-online.eu)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Affero General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * This program 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 Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+// cOkadaEarthquake.cpp: implementation of the cOkadaEarthquake class
+//
+//=========================================================================
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <utilits.h>
+#include "cOkadaEarthquake.h"
+
+
+
+//=========================================================================
+cOkadaEarthquake::cOkadaEarthquake()
+{
+  finalized = 0;
+  m0 = 0;
+  nfault = 1;
+  fault = new cOkadaFault [1];
+}
+
+
+//=========================================================================
+cOkadaEarthquake::~cOkadaEarthquake()
+{
+  if( fault != NULL ) delete [] fault;
+}
+
+
+//=========================================================================
+// Read fault(s) from a file
+int cOkadaEarthquake::read( char *fname )
+{
+  FILE *fp;
+  char record[256];
+  int ierr;
+
+  nfault = utlGetNumberOfRecords( fname );
+  if( fault != NULL ) delete [] fault;
+  fault = new cOkadaFault [nfault];
+
+  if( (fp = fopen(fname,"rt")) == NULL ) return Err.post(Err.msgOpenFile(fname));
+
+  for( int n=0, line=0; n<nfault; n++ ) {
+    ierr = utlReadNextRecord(fp, record, &line); if( ierr == EOF ) return Err.post(Err.msgReadFile(fname, line, "Unexpected EOF"));
+    ierr = fault[n].read( record ); if(ierr) return ierr;
+  }
+
+  fclose( fp );
+
+  return 0;
+}
+
+
+//=========================================================================
+// check all ruptures for parameter integrity
+int cOkadaEarthquake::finalizeInput()
+{
+  char msg[32];
+  int ierr,n;
+
+
+  for( m0=0., n=0; n<nfault; n++ ) {
+
+    ierr = fault[n].check();
+
+    if( ierr ) {
+      Err.post( "Fault number: %d", n+1 );
+      return( 10*n + ierr );
+    }
+
+    m0 += fault[n].getM0();
+
+  }
+
+  finalized = 1;
+
+  return 0;
+}
+
+
+//=========================================================================
+double cOkadaEarthquake::getM0()
+{
+  if( !finalized ) { Err.post( "cOkadaEarthquake::getM0: eq not finalized" ); return -RealMax; }
+
+  return m0;
+}
+
+//=========================================================================
+double cOkadaEarthquake::getMw()
+{
+  if( !finalized ) { Err.post( "cOkadaEarthquake::getMw: eq not finalized" ); return -RealMax; }
+
+  return( 2./3.*(log10(m0)-9.1) );
+}
+
+
+//=========================================================================
+// get ruptured area to calculate surface displacement
+int cOkadaEarthquake::getDeformArea( int round, double& lonmin, double& lonmax, double& latmin, double& latmax )
+{
+  int ierr;
+  double xmin,xmax,ymin,ymax;
+
+
+  if( !finalized ) return Err.post("cOkadaEarthquake::calculate: eq not finalized");
+
+  lonmin = latmin = RealMax;
+  lonmax = latmax = -RealMax;
+
+  for( int n=0; n<nfault; n++ ) {
+
+    ierr = fault[n].getDeformArea( xmin, xmax, ymin, ymax ); if(ierr) return ierr;
+
+    if( xmin < lonmin ) lonmin = xmin; if( xmax > lonmax ) lonmax = xmax;
+    if( ymin < latmin ) latmin = ymin; if( ymax > latmax ) latmax = ymax;
+  }
+
+  if( round ) {
+    lonmin = floor( lonmin ); lonmax = ceil( lonmax );
+    latmin = floor( latmin ); latmax = ceil( latmax );
+  }
+
+  return 0;
+}
+
+
+//=========================================================================
+// Calculate surface displacements by summing effect from all ruptures
+int cOkadaEarthquake::calculate( double lon, double lat, double& uz, double& ulon, double &ulat )
+{
+  int n,ierr;
+  double uzf,ulonf,ulatf;
+
+
+  if( !finalized ) return Err.post("cOkadaEarthquake::calculate: eq not finalized");
+
+  for( ulon=ulat=uz=0, n=0; n<nfault; n++ ) {
+    ierr = fault[n].calculate( lon, lat, uzf, ulonf, ulatf ); if(ierr) return ierr;
+    uz+= uzf; ulon += ulonf; ulat += ulatf;
+  }
+
+  return 0;
+}
+
+
+//=========================================================================
+// Calculate surface displacements by summing effect from all ruptures
+int cOkadaEarthquake::calculate( double lon, double lat, double& uz )
+{
+  int n,ierr;
+  double uzf;
+
+
+  if( !finalized ) return Err.post("cOkadaEarthquake::calculate: eq not finalized");
+
+  for( uz=0, n=0; n<nfault; n++ ) {
+    ierr = fault[n].calculate( lon, lat, uzf ); if(ierr) return ierr;
+    uz += uzf;
+  }
+
+  return 0;
+}
+
+
+//=========================================================================
+// Calculate surface displacements by summing effect from all ruptures
+int cOkadaEarthquake::calculate( cOgrd& uZ )
+{
+  int i,j,n,ierr;
+  double uzf;
+
+
+  if( !finalized ) return Err.post("cOkadaEarthquake::calculate: eq not finalized");
+  ierr = setGrid( uZ ); if( ierr ) return ierr;
+
+  // calculate displacenents on a grid
+  for( i=0; i<uZ.nx; i++ ) {
+    for( j=0; j<uZ.ny; j++ ) {
+
+      for( n=0; n<nfault; n++ ) {
+        ierr = fault[n].calculate( uZ.getX(i,j), uZ.getY(i,j), uzf ); if(ierr) return ierr;
+        uZ(i,j) = uZ(i,j) + uzf;
+      }
+
+    }
+  }
+
+  return 0;
+}
+
+
+//=========================================================================
+// Calculate surface displacements by summing effect from all ruptures
+int cOkadaEarthquake::calculate( cOgrd& uZ, cOgrd& uLon, cOgrd& uLat )
+{
+  int i,j,n,ierr;
+  double uzf,ulonf,ulatf;
+
+  if( !finalized ) return Err.post("cOkadaEarthquake::calculate: eq not finalized");
+  ierr = setGrid( uZ ); if(ierr) return ierr;
+
+  uLon = uZ;
+  uLat = uZ;
+
+  // calculate displacenents on a grid
+  for( i=0; i<uZ.nx; i++ ) {
+    for( j=0; j<uZ.ny; j++ ) {
+
+      for( n=0; n<nfault; n++ ) {
+        ierr = fault[n].calculate( uZ.getX(i,j), uZ.getY(i,j), uzf, ulonf, ulatf ); if(ierr) return ierr;
+        uZ(i,j) = uZ(i,j) + uzf;
+        uLon(i,j) = uLon(i,j) + ulonf;
+        uLat(i,j) = uLat(i,j) + ulatf;
+      }
+
+    }
+  }
+
+  return 0;
+}
+
+
+//=========================================================================
+// Check and prepare grid for calculations
+int cOkadaEarthquake::setGrid( cOgrd& u )
+{
+  int ierr;
+