Commit 5acb5fe2 authored by Johannes Spazier's avatar Johannes Spazier

Added original EasyWave version for CPU extended by GPU functionality.

parents
#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;
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...");
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
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();
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( "\nExample:\n" );
printf( "\t easyWave -grid gebcoIndonesia.grd -source fault.inp -time 120\n\n" );
return -1;
}
-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
This diff is collapsed.
#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
// 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;
if( u.xmin == u.xmax || u.ymin == u.ymax ) {
ierr = getDeformArea( 1, u.xmin, u.xmax, u.ymin, u.ymax ); if(ierr) return ierr;
}
if( u.nx*u.ny != 0 && u.dx*u.dy != 0 )
u.resetVal();
else if( u.nx*u.ny != 0 )
u.initialize( u.xmin, u.xmax, u.nx, u.ymin, u.ymax, u.ny );
else if( u.dx*u.dy != 0 )
u.initialize( u.xmin, u.xmax, u.dx, u.ymin, u.ymax, u.dy );
else {
u.dx = u.dy = 0.02; // a bit more than 1 arc minute (0.0166667)
u.initialize( u.xmin, u.xmax, u.dx, u.ymin, u.ymax, u.dy );
}
return 0;
}
//=========================================================================
int cOkadaEarthquake::calculate( cObsArray& arr )
{
int k,n,ierr;
double uzf,ulonf,ulatf;
if( !finalized ) return Err.post("cOkadaEarthquake::calculate: eq not finalized");
for( k=0; k<arr.nPos; k++ ) {
for( arr.obs[k][0]=arr.obs[k][1]=arr.obs[k][2]=0., n=0; n<nfault; n++ ) {
ierr = fault[n].calculate( arr.lon[k], arr.lat[k], uzf, ulonf, ulatf ); if(ierr) return ierr;
arr.obs[k][0] += ulonf;
arr.obs[k][1] += ulatf;
arr.obs[k][2] += uzf;
}
}
return 0;
}
//=========================================================================
// print earthquake parameter data into a new string
char* cOkadaEarthquake::sprint()
{
char *info,buf[256];
int n,bufsize;
bufsize = 64 + nfault*128;
info = new char[bufsize];
sprintf( info, "Mw=%.2f M0=%-g", getMw(), getM0() );
sprintf( buf, "\nNumber of faults: %-3d", nfault ); strcat( info, buf );
for( n=0; n<nfault; n++ ) {
sprintf( buf, "\n Fault %d(%d): length=%-.1f width=%-.1f slip=%-.2f top=%-.2f",
n+1, nfault, fault[n].length/1000, fault[n].width/1000, fault[n].slip, fault[n].getZtop()/1000 );
strcat( info, buf );
}
return info;
}
#ifndef OKADAEARTHQUAKE_H
#define OKADAEARTHQUAKE_H
#include <cOgrd.h>
#include <cSphere.h>
#include "cOkadaFault.h"
class cOkadaEarthquake
{
protected:
int finalized;
int getDeformArea( int round, double& lonmin, double& lonmax, double& latmin, double& latmax );
int setGrid( cOgrd& u );
public:
int nfault; // total nuber of Okada faults
double m0; // total earthquake moment
cOkadaFault *fault; // array of composing faults
cOkadaEarthquake();
~cOkadaEarthquake();
int read( char *fname );
int finalizeInput();
double getM0();
double getMw();
int calculate( double lon, double lat, double& uz );
int calculate( double lon, double lat, double& uz, double& ulon, double &ulat );
int calculate( cObsArray& arr );
int calculate( cOgrd& uZ );
int calculate( cOgrd& uZ, cOgrd& uLon, cOgrd& uLat );
char *sprint();
};
#endif // OKADAEARTHQUAKE_H
#include <stdio.h>
#include <string.h>
#include <utilits.h>
#include "cOkadaFault.h"
static const double Rearth=6384.e+3;
//=========================================================================
// Constructor
cOkadaFault::cOkadaFault()
{
// obligatory user parameters
lat = lon = depth = strike = dip = rake = RealMax;
// optional user parameters
mw = slip = length = width = 0.;
adjust_depth = 0;
// default values
refpos = FLT_POS_C;
mu = 3.5e+10;
// derivative parameters
zbot = sind = cosd = sins = coss = tand = coslat = wp = dslip = sslip = 0;
// flags
checked = 0;
}
//=========================================================================