Skip to content

Commit

Permalink
Save wave group to netcdf file and read netcdf file
Browse files Browse the repository at this point in the history 
Need to double check tslen in writebndnc
  • Loading branch information
@CyprienBosserelle
CyprienBosserelle committed Jan 27, 2017
1 parent 3a5cfd7 commit 1882711
Show file tree
Hide file tree
Showing 7 changed files with 446 additions and 58 deletions.
3 changes: 2 additions & 1 deletion Wave_gpu.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1753,7 +1753,8 @@ int main(int argc, char **argv)
if (XParam.wavebndtype == 3)
{
//Reuse XBeach_GPU boundary. Similar to XBeach but in netcdf format and all self explanatory
//Not implemented yet
XParam = read_reuse_bndnc_head(XParam);
wavebnd = read_reuse_bndnc_vec(XParam);
}
if (XParam.wavebndtype == 4)
{
Expand Down
140 changes: 137 additions & 3 deletions Wavestep.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ XBGPUParam waveinitGPU(XBGPUParam Param, std::vector<Wavebndparam> wavebnd)
}
if (Param.wavebndtype >= 2)
{
nwavbnd = ceil(Param.rtlength / Param.dtbc);
nwavbnd = ceil(Param.rtlength / Param.dtbc)+1; // +1 needed here
}

theta = (DECNUM *)malloc(ntheta*sizeof(DECNUM));
Expand Down Expand Up @@ -126,7 +126,7 @@ XBGPUParam waveinitGPU(XBGPUParam Param, std::vector<Wavebndparam> wavebnd)
if (Param.wavebndtype == 3)
{
// Reuse XBeach_GPU style wave boundary. same as normal XBeach but as a self documented netcdf file
//Not implemented yet
read_reuse_bndnc(Param, 0, Trep, qfile, Stfile);
}
if (Param.wavebndtype == 4)
{
Expand All @@ -144,7 +144,61 @@ XBGPUParam waveinitGPU(XBGPUParam Param, std::vector<Wavebndparam> wavebnd)
//void GenWGnLBW(XBGPUParam Param, int nf, int ndir, double * HRfreq, double * HRdir, double * HRSpec, float Trep, double * qfile, double * Stfile)
GenWGnLBW(Param, nfHR, ndHR, HRfreq, HRdir, HRSpec, Trep, qfile, Stfile);
//create2dnc(nfHR, ndHR, HRfreq[1] - HRfreq[0], HRdir[1] - HRdir[0], 0.0, HRfreq, HRdir, HRSpec);

//////////////////////////////////////
//Save to Netcdf file
//////////////////////////////////////
double * yyfx, *ttfx, *thetafx;
double * qxtemp, *qytemp, *eetemp;
int tslenbc = nwavbnd;

qxtemp = (double *)malloc(ny*tslenbc*sizeof(double));
qytemp = (double *)malloc(ny*tslenbc*sizeof(double));
eetemp = (double *)malloc(ny*Param.ntheta*tslenbc*sizeof(double));

yyfx = (double *)malloc(ny*sizeof(double));
ttfx = (double *)malloc(tslenbc*sizeof(double));
thetafx = (double *)malloc(Param.ntheta*sizeof(double));

for (int j = 0; j < Param.ny; j++)
{
yyfx[j] = j*Param.dx;
}

for (int m = 0; m < tslenbc; m++)
{
ttfx[m] = m*Param.dtbc;
}

for (int itheta = 0; itheta < Param.ntheta; itheta++)
{
thetafx[itheta] = itheta*(Param.dtheta) + Param.thetamin + 0.5f*Param.dtheta;
}

//for Debugging
//create2dnc(ny, tslen, Param.dx, dtin, 0.0, yyfx, tin, qx);
//create3dnc(ny, Param.ntheta, tslen, Param.dx, Param.dtheta, dtin, 0.0, yyfx, thetafx, tin, zeta);

for (int j = 0; j < Param.ny; j++)
{
for (int m = 0; m < tslenbc; m++)
{
qxtemp[m + j*tslenbc] = qfile[j + 0 * ny + m*ny * 4];
qytemp[m + j*tslenbc] = qfile[j + 1 * ny + m*ny * 4];

for (int itheta = 0; itheta < Param.ntheta; itheta++)
{
eetemp[m + j*tslenbc + itheta*ny*tslenbc] = Stfile[j + itheta*ny + m*ny*Param.ntheta];
}
}
}


createbndnc(tslenbc, ny, Param.ntheta, Param.dx, Param.dtheta, 0.0, ttfx, yyfx, thetafx, eetemp, qxtemp, qytemp);
//Trep = 15.0;//
free(HRSpec);
free(HRfreq);
free(HRdir);



Expand Down Expand Up @@ -382,10 +436,90 @@ void wavebnd(XBGPUParam Param, std::vector<Wavebndparam> wavebndvec)

if (Param.wavebndtype == 2)
{
//Read new STfile and qfile
//Read new STfile and qfile XBeach style
readXbbndstep(Param, wavebndvec, WAVstepinbnd - 1, Trep, qfile, Stfile);

}

if (Param.wavebndtype == 3)
{
// Reuse XBeach_GPU style wave boundary. same as normal XBeach but as a self documented netcdf file
read_reuse_bndnc(Param, WAVstepinbnd - 1, Trep, qfile, Stfile);
}

if (Param.wavebndtype == 4)
{
//JONSWAP
//First generate a Highres 2D spec
double * HRfreq;
double * HRdir;
double * HRSpec;

int nfHR, ndHR;

makjonswap(Param, wavebndvec, WAVstepinbnd - 1, nfHR, ndHR, HRfreq, HRdir, HRSpec);

//Then generate wave group timeseries based on that spectra
//void GenWGnLBW(XBGPUParam Param, int nf, int ndir, double * HRfreq, double * HRdir, double * HRSpec, float Trep, double * qfile, double * Stfile)
GenWGnLBW(Param, nfHR, ndHR, HRfreq, HRdir, HRSpec, Trep, qfile, Stfile);
//////////////////////////////////////
//Save to Netcdf file
//////////////////////////////////////

// Stfile is not ordered teh way we want to save it to file so we need a temporary storage to rearange
double * yyfx, *ttfx, *thetafx;
double * qxtemp, *qytemp, *eetemp;
int tslenbc = (int)ceil(Param.rtlength / Param.dtbc)+1;

qxtemp = (double *)malloc(ny*tslenbc*sizeof(double));
qytemp = (double *)malloc(ny*tslenbc*sizeof(double));
eetemp = (double *)malloc(ny*Param.ntheta*tslenbc*sizeof(double));

yyfx = (double *)malloc(ny*sizeof(double));
ttfx = (double *)malloc(tslenbc*sizeof(double));
thetafx = (double *)malloc(Param.ntheta*sizeof(double));

for (int j = 0; j < Param.ny; j++)
{
yyfx[j] = j*Param.dx;
}

for (int m = 0; m < tslenbc; m++)
{
ttfx[m] = m*Param.dtbc;
}

for (int itheta = 0; itheta < Param.ntheta; itheta++)
{
thetafx[itheta] = itheta*(Param.dtheta) + Param.thetamin + 0.5f*Param.dtheta;
}

//for Debugging
//create2dnc(ny, tslen, Param.dx, dtin, 0.0, yyfx, tin, qx);
//create3dnc(ny, Param.ntheta, tslen, Param.dx, Param.dtheta, dtin, 0.0, yyfx, thetafx, tin, zeta);

for (int j = 0; j < Param.ny; j++)
{
for (int m = 0; m < tslenbc; m++)
{
qxtemp[m + j*tslenbc] = qfile[j + 0 * ny + m*ny * 4];
qytemp[m + j*tslenbc] = qfile[j + 1 * ny + m*ny * 4];

for (int itheta = 0; itheta < Param.ntheta; itheta++)
{
eetemp[m + j*tslenbc + itheta*ny*tslenbc] = Stfile[j + itheta*ny + m*ny*Param.ntheta];
}
}
}


writebndnc(tslenbc, ny, Param.ntheta, Param.dx, Param.dtheta, wavebndvec[WAVstepinbnd-1].time, ttfx, yyfx, thetafx, eetemp, qxtemp, qytemp);
//Trep = 15.0;//
free(HRSpec);
free(HRfreq);
free(HRdir);

}
}

if (Param.wavebndtype == 1)
Expand Down
6 changes: 6 additions & 0 deletions XBeachGPU.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -217,6 +217,8 @@ extern "C" void read3Dnc(int nx, int ny,int ntheta,char ncfile[],DECNUM * &ee);
extern "C" void read2Dnc(int nx, int ny,char ncfile[],DECNUM * &hh);

void createbndnc(int tslen, int ny, int ntheta, double dy, double dtheta, double totaltime, double * timevec, double *yy, double *theta, double * ee, double * qx, double * qy);
void writebndnc(int tslen, int ny, int ntheta, double dy, double dtheta, double totaltime, double * timevec, double *yy, double *theta, double * ee, double * qx, double * qy);
void read_reuse_bndnc(XBGPUParam Param, int rec, float &Trep, double * &qfile, double * &Stfile);

void GenCstWave(XBGPUParam Param, std::vector<Wavebndparam> wavebnd, float * theta, double * &Stfile, double * &qfile, double * &Tpfile);

Expand All @@ -229,6 +231,9 @@ extern "C" void readbndhead(const char * wavebndfile,DECNUM &thetamin,DECNUM &th

std::vector<Wavebndparam> ReadJSWPBnd(XBGPUParam XParam);

XBGPUParam read_reuse_bndnc_head(XBGPUParam Param);
std::vector<Wavebndparam> read_reuse_bndnc_vec(XBGPUParam Param);

//Below is for the new CPU routine
extern "C" int mminusC(int ix,int nx);
extern "C" int pplusC(int ix, int nx);
Expand Down Expand Up @@ -285,6 +290,7 @@ std::vector<WindBnd> readWNDfile(std::string WNDfilename, double grdalpha);
std::vector<Wavebndparam> ReadCstBnd(XBGPUParam XParam);
double interptime(double next, double prev, double timenext, double time);
double interp1D(int nx, double *x, double *y, double xx);
double interp1DMono(int nx, double *x, double *y, double xx);
double Interp2(int nx, int ny, double *x, double *y, double *z, double xx, double yy);

extern "C" void readbathyHead(std::string filename, int &nx, int &ny, double &dx, double &grdalpha);
Expand Down
68 changes: 14 additions & 54 deletions makjonswap.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,8 @@ void makjonswap(XBGPUParam Param, std::vector<Wavebndparam> wavebnd, int step, i
double scoeff = wavebnd[step].s;


//printf("fp=%f\n",fp);
printf("Generating JONSWAP spectrum: Hs=%f, Tp=%f, Dp=%f, gam=%f, s=%f\n",Hs,Tp,Dp,gam,scoeff);
write_text_to_log_file("Generating JONSWAP spectrum: Hs=" + std::to_string(Hs) + " Tp=" + std::to_string(Tp) + " Dp=" + std::to_string(Dp) + " gam=" + std::to_string(gam) + " s=" + std::to_string(scoeff) + ";");
////

double fnyq = 3.0f*fp;
Expand Down Expand Up @@ -222,6 +223,11 @@ void makjonswap(XBGPUParam Param, std::vector<Wavebndparam> wavebnd, int step, i

void GenWGnLBW(XBGPUParam Param, int nf, int ndir,double * HRfreq,double * HRdir, double * HRSpec, float &Trep, double * &qfile, double * &Stfile)
{
// Generating Boundary condition: Energy from wave group and Long bound waves

printf("Generating Boundary condition: Energy from wave group and Long bound waves.\n");
write_text_to_log_file("Generating Boundary condition: Energy from wave group and Long bound waves");

int ny = Param.ny;
int K; //Number of wave components

Expand Down Expand Up @@ -414,6 +420,8 @@ void GenWGnLBW(XBGPUParam Param, int nf, int ndir,double * HRfreq,double * HRdir
//}
//interp1D(double *x, double *y, double xx)
//thetagen[i] = interptime(HRdir[dprev + 1], HRdir[dprev], dtheta, dtheta*((number - cdf[dprev]) / (cdf[dprev + 1] - cdf[dprev])));

//Cannot use interp1DMono here because cdf is not monotonic
thetagen[i] = interp1D(ndir, cdf, HRdir, number);

//printf("thetagen[i]=%f\n", thetagen[i]);
Expand Down Expand Up @@ -470,7 +478,7 @@ void GenWGnLBW(XBGPUParam Param, int nf, int ndir,double * HRfreq,double * HRdir

//First assume that internal and bc - writing time step is the same
double dtin = Param.dtbc;
int tslenbc = ceil(Param.rtlength / Param.dtbc);// (int)(Param.rtlength / Param.dtbc) + 1;
int tslenbc = ceil(Param.rtlength / Param.dtbc) + 1;// (int)(Param.rtlength / Param.dtbc) + 1;

//Check whether the internal frequency is high enough to describe the highest frequency
//wave train returned from frange(which can be used in the boundary conditions)
Expand Down Expand Up @@ -565,7 +573,7 @@ void GenWGnLBW(XBGPUParam Param, int nf, int ndir,double * HRfreq,double * HRdir
double dummy;
for (int i = 0; i < K; i++)
{
dummy = interp1D(nf, HRfreq, Sf, fgen[i]);
dummy = interp1DMono(nf, HRfreq, Sf, fgen[i]);
vargenq[i] = min(vargen[i] ,dummy);
}

Expand Down Expand Up @@ -944,7 +952,7 @@ void GenWGnLBW(XBGPUParam Param, int nf, int ndir,double * HRfreq,double * HRdir
//interpolate to boundary timeseries
for (int m = 0; m < tslenbc; m++)
{
Stfile[j + itheta*ny + m*ny*Param.ntheta] = interp1D(tslen, tin, E_tdir, m*Param.dtbc);
Stfile[j + itheta*ny + m*ny*Param.ntheta] = interp1DMono(tslen, tin, E_tdir, m*Param.dtbc);
}

}
Expand Down Expand Up @@ -1323,63 +1331,15 @@ void GenWGnLBW(XBGPUParam Param, int nf, int ndir,double * HRfreq,double * HRdir

for (int m = 0; m < tslenbc; m++)
{
qfile[j + 0 * ny + m*ny * 4] = interp1D(tslen, tin, qtempx, m*Param.dtbc);
qfile[j + 1 * ny + m*ny * 4] = interp1D(tslen, tin, qtempy, m*Param.dtbc);
qfile[j + 0 * ny + m*ny * 4] = interp1DMono(tslen, tin, qtempx, m*Param.dtbc);
qfile[j + 1 * ny + m*ny * 4] = interp1DMono(tslen, tin, qtempy, m*Param.dtbc);
qfile[j + 2 * ny + m*ny * 4] = 0.0;
qfile[j + 3 * ny + m*ny * 4] = 0.0;// qtot[j + m*ny];
}

}

//////////////////////////////////////
//Save to Netcdf file
//////////////////////////////////////
double * yyfx, *ttfx, *thetafx;
double * qxtemp, *qytemp, *eetemp;

qxtemp = (double *)malloc(ny*tslenbc*sizeof(double));
qytemp = (double *)malloc(ny*tslenbc*sizeof(double));
eetemp = (double *)malloc(ny*Param.ntheta*tslenbc*sizeof(double));

yyfx = (double *)malloc(ny*sizeof(double));
ttfx = (double *)malloc(tslenbc*sizeof(double));
thetafx = (double *)malloc(Param.ntheta*sizeof(double));

for (int j = 0; j < Param.ny; j++)
{
yyfx[j] = j*Param.dx;
}

for (int m = 0; m < tslenbc; m++)
{
ttfx[m] = m*Param.dtbc;
}

for (int itheta = 0; itheta < Param.ntheta; itheta++)
{
thetafx[itheta] = itheta*(Param.dtheta) + Param.thetamin + 0.5f*Param.dtheta;
}

//for Debugging
create2dnc(ny, tslen, Param.dx, dtin, 0.0, yyfx, tin, qx);
create3dnc(ny, Param.ntheta, tslen, Param.dx, Param.dtheta, dtin, 0.0, yyfx, thetafx, tin, zeta);

for (int j = 0; j < Param.ny; j++)
{
for (int m = 0; m < tslenbc; m++)
{
qxtemp[m + j*tslenbc] = qfile[j + 0 * ny + m*ny * 4];
qytemp[m + j*tslenbc] = qfile[j + 1 * ny + m*ny * 4];

for (int itheta = 0; itheta < Param.ntheta; itheta++)
{
eetemp[m + j*tslenbc + itheta*ny*tslenbc] = Stfile[j + itheta*ny + m*ny*Param.ntheta];
}
}
}


createbndnc(tslenbc, ny, Param.ntheta, Param.dx, Param.dtheta, 0.0, ttfx, yyfx,thetafx, eetemp, qxtemp, qytemp);

//////////////////////////////////////
//Clean up and desallocate all that memory
Expand Down
3 changes: 3 additions & 0 deletions read_input.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@
using DECNUM = float;



XBGPUParam readXbbndhead(XBGPUParam Param)
{
std::ifstream fs(Param.wavebndfile);
Expand Down Expand Up @@ -107,6 +108,8 @@ std::vector<Wavebndparam> readXbbndfile(XBGPUParam Param)
return wavebndvec;
}



extern "C" void readbndhead(const char * wavebnd, DECNUM &thetamin, DECNUM &thetamax, DECNUM &dtheta, DECNUM &dtwavbnd, int &nwavbnd)
{
FILE * fwav;
Expand Down
32 changes: 32 additions & 0 deletions tools.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,8 @@ double interptime(double next, double prev, double timenext, double time)

double interp1D(int nx, double *x, double *y, double xx)
{
// non monotonic 1D interpolation
// This can be pretty slow so use interp1DMono for faster interpolation if x is monotonic (i.e. the spacing is constant)
double yy;
double prevx = x[0];
double nextx = x[1];
Expand Down Expand Up @@ -58,6 +60,36 @@ double interp1D(int nx, double *x, double *y, double xx)
yy = prevy + (xx - prevx) / (nextx - prevx)*(nexty - prevy);
return yy;
}

double interp1DMono(int nx, double *x, double *y, double xx)
{
// interpolation if x is monotonically increasing (i.e. the spacing is constant)
double yy;
double prevx = x[0];
double nextx = x[1];
double prevy = y[0];
double nexty = y[1];
double dx = nextx - prevx;
int indx = 0;
int indxp;

double diffx = max(xx-x[0],0.0);//This has to be positive!

indx = (int)floor(diffx / dx);



indxp = (int)min(indx*1.0 + 1, nx*1.0 - 1);
prevx = x[indx];
nextx = x[indxp];
prevy = y[indx];
nexty = y[indxp];


yy = prevy + (xx - prevx) / (nextx - prevx)*(nexty - prevy);
return yy;
}

double Interp2(int nx, int ny, double *x, double *y, double *z, double xx, double yy)
//double BilinearInterpolation(double q11, double q12, double q21, double q22, double x1, double x2, double y1, double y2, double x, double y)
{
Expand Down
Loading

0 comments on commit 1882711

Please sign in to comment.