Skip to content

Commit

Permalink
COmpleted Sl Bnd update
Browse files Browse the repository at this point in the history 
This was compiled and tested
  • Loading branch information
@CyprienBosserelle
CyprienBosserelle committed Nov 16, 2016
1 parent e8d13ed commit 3601e0c
Show file tree
Hide file tree
Showing 5 changed files with 88 additions and 38 deletions.
6 changes: 3 additions & 3 deletions Flow_CPU.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ extern "C" int pplus2C(int ix, int nx)



extern "C" void ubndCPU(int nx, int ny, DECNUM dx, DECNUM dt,DECNUM g, DECNUM rho,DECNUM totaltime,DECNUM wavbndtime,DECNUM rt,DECNUM slbndtime, DECNUM rtsl,DECNUM zsbndold,DECNUM zsbndnew,DECNUM Trep,DECNUM * qbndold, DECNUM * qbndnew,DECNUM *&zs, DECNUM * &uu,DECNUM * &vv, DECNUM *vu, DECNUM * umean, DECNUM * vmean,DECNUM * zb,DECNUM * cg,DECNUM * hum, DECNUM * zo, DECNUM *Fx,DECNUM *&hh)
extern "C" void ubndCPU(int nx, int ny, DECNUM dx, DECNUM dt,DECNUM g, DECNUM rho,DECNUM totaltime,DECNUM wavbndtime,DECNUM rt,DECNUM zsbnd, DECNUM Trep,DECNUM * qbndold, DECNUM * qbndnew,DECNUM *&zs, DECNUM * &uu,DECNUM * &vv, DECNUM *vu, DECNUM * umean, DECNUM * vmean,DECNUM * zb,DECNUM * cg,DECNUM * hum, DECNUM * zo, DECNUM *Fx,DECNUM *&hh)
{
int ix=0;

Expand All @@ -92,7 +92,7 @@ extern "C" void ubndCPU(int nx, int ny, DECNUM dx, DECNUM dt,DECNUM g, DECNUM rh
DECNUM epsi=0.005; //Not used!
DECNUM ur,uumean,vvmean,urr,alphanew;
DECNUM dbetadx,dbetady,dvudy,dhdx;
DECNUM qx,qy,zsbnd;
DECNUM qx,qy;
DECNUM order=2.0f;
DECNUM ccg=cg[i];
DECNUM cats=4; // number of wave period to average the current from
Expand All @@ -103,7 +103,7 @@ extern "C" void ubndCPU(int nx, int ny, DECNUM dx, DECNUM dt,DECNUM g, DECNUM rh

qx=(qbndold[iy]+(totaltime-wavbndtime+rt)*(qbndnew[iy]-qbndold[iy])/rt)*taper;
qy=(qbndold[iy+ny]+(totaltime-wavbndtime+rt)*(qbndnew[iy+ny]-qbndold[iy+ny])/rt)*taper;
zsbnd=zsbndold+(totaltime-rtsl)*(zsbndnew-zsbndold)/(slbndtime-rtsl);
//zsbnd=zsbndold+(totaltime-rtsl)*(zsbndnew-zsbndold)/(slbndtime-rtsl);

ht=zsbnd+zb[i];
htr=zsbnd+zb[xplus+iy*nx];
Expand Down
6 changes: 3 additions & 3 deletions Flow_kernel.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -144,7 +144,7 @@ __global__ void ubnd(int nx, int ny, DECNUM dx, DECNUM dt, DECNUM g, DECNUM rho,


}
__global__ void ubnd1D(int nx, int ny, DECNUM dx, DECNUM dt, DECNUM g, DECNUM rho, DECNUM totaltime, DECNUM wavbndtime, DECNUM rt, DECNUM slbndtime, DECNUM rtsl, DECNUM zsbndold, DECNUM zsbndnew, DECNUM Trep, DECNUM * qbndold, DECNUM * qbndnew, DECNUM *zs, DECNUM * uu, DECNUM * vv, DECNUM *vu, DECNUM * umean, DECNUM * vmean, DECNUM * zb, DECNUM * cg, DECNUM * hum, DECNUM * zo, DECNUM *Fx, DECNUM *hh)
__global__ void ubnd1D(int nx, int ny, DECNUM dx, DECNUM dt, DECNUM g, DECNUM rho, DECNUM totaltime, DECNUM wavbndtime, DECNUM rt, DECNUM zsbnd, DECNUM Trep, DECNUM * qbndold, DECNUM * qbndnew, DECNUM *zs, DECNUM * uu, DECNUM * vv, DECNUM *vu, DECNUM * umean, DECNUM * vmean, DECNUM * zb, DECNUM * cg, DECNUM * hum, DECNUM * zo, DECNUM *Fx, DECNUM *hh)
{
unsigned int ix = blockIdx.x*blockDim.x + threadIdx.x;
unsigned int iy = blockIdx.y*blockDim.y + threadIdx.y;
Expand All @@ -164,7 +164,7 @@ __global__ void ubnd1D(int nx, int ny, DECNUM dx, DECNUM dt, DECNUM g, DECNUM rh
DECNUM epsi = 0.005; //Not used!
DECNUM ur, uumean, vvmean, urr, alphanew;

DECNUM qx, qy, zsbnd;
DECNUM qx, qy;

DECNUM cats = 4.0f; // number of wave period to average the current from
DECNUM factime = 1.0f/cats/Trep*dt;
Expand All @@ -175,7 +175,7 @@ __global__ void ubnd1D(int nx, int ny, DECNUM dx, DECNUM dt, DECNUM g, DECNUM rh
qx = (qbndold[iy] + (totaltime - wavbndtime + rt)*(qbndnew[iy] - qbndold[iy]) / rt)*taper;
qy = (qbndold[iy + ny] + (totaltime - wavbndtime + rt)*(qbndnew[iy + ny] - qbndold[iy + ny]) / rt)*taper;

zsbnd = zsbndold + (totaltime - rtsl)*(zsbndnew - zsbndold) / (slbndtime - rtsl);
//zsbnd = zsbndold + (totaltime - rtsl)*(zsbndnew - zsbndold) / (slbndtime - rtsl);

ht = zsbnd + zb[i];
htr = zsbnd + zb[xplus + iy*nx];
Expand Down
83 changes: 54 additions & 29 deletions Wave_gpu.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -267,7 +267,7 @@ exit(EXIT_FAILURE);


// Main loop that actually runs the model
void mainloopGPU(XBGPUParam Param)
void mainloopGPU(XBGPUParam Param, std::vector<SLBnd> slbnd)
{
double dt = Param.dt;
int nx, ny;
Expand All @@ -286,7 +286,7 @@ void mainloopGPU(XBGPUParam Param)
wdt = dt; // previous timestep

//Calculate timestep
FLOWDT << <gridDim, blockDim, 0 >> >(nx, ny, Param.dx, 0.8f, dtflow_g, hh_g);
FLOWDT << <gridDim, blockDim, 0 >> >(nx, ny, Param.dx, Param.CFL, dtflow_g, hh_g);
CUDA_CHECK(cudaThreadSynchronize());


Expand All @@ -300,15 +300,15 @@ void mainloopGPU(XBGPUParam Param)
//CUDA_CHECK(cudaMemcpy(arrmax, arrmax_g, nx*ny*sizeof(DECNUM), cudaMemcpyDeviceToHost));
CUDA_CHECK(cudaMemcpy(arrmin, arrmin_g, nx*ny*sizeof(DECNUM), cudaMemcpyDeviceToHost));

dt = arrmin[0] * 0.5f;
dt = arrmin[0] * 0.5; // Not sure why this is but it is in the original XBeach!!



if ((Param.swave == 1 ) && totaltime>0.0)
{
float dtwave;
double dtwave;
// Make sure the CFL condition for flow do not violate CFL condition for Waves
WAVEDT << <gridDim, blockDim, 0 >> >(nx, ny, ntheta, 0.8f, dtheta, dtflow_g, ctheta_g);
WAVEDT << <gridDim, blockDim, 0 >> >(nx, ny, ntheta, Param.CFL, dtheta, dtflow_g, ctheta_g);
CUDA_CHECK(cudaThreadSynchronize());


Expand All @@ -322,7 +322,7 @@ void mainloopGPU(XBGPUParam Param)
//CUDA_CHECK(cudaMemcpy(arrmax, arrmax_g, nx*ny*sizeof(DECNUM), cudaMemcpyDeviceToHost));
CUDA_CHECK(cudaMemcpy(arrmin, arrmin_g, nx*ny*sizeof(DECNUM), cudaMemcpyDeviceToHost));

dtwave = arrmin[0]*0.5f;
dtwave = arrmin[0]*0.5;

dt = min(dt, dtwave);
}
Expand Down Expand Up @@ -357,7 +357,7 @@ void mainloopGPU(XBGPUParam Param)

if (Param.flow == 1)
{
flowbnd(Param);// Calculate the flow boundary for this step
flowbnd(Param,slbnd);// Calculate the flow boundary for this step
}

if (Param.swave == 1 )
Expand Down Expand Up @@ -509,7 +509,7 @@ void mainloopGPU(XBGPUParam Param)
}


void mainloopCPU(XBGPUParam Param)
void mainloopCPU(XBGPUParam Param,std::vector<SLBnd> slbnd)
{
printf("Computing CPU mode\n");

Expand All @@ -533,7 +533,7 @@ void mainloopCPU(XBGPUParam Param)

if (Param.flow == 1)
{
flowbnd(Param);// Calculate the flow boundary for this step
flowbnd(Param,slbnd);// Calculate the flow boundary for this step
}
if (Param.swave == 1)
{
Expand Down Expand Up @@ -590,25 +590,42 @@ void mainloopCPU(XBGPUParam Param)



void flowbnd(XBGPUParam Param)
void flowbnd(XBGPUParam Param,std::vector<SLBnd> slbnd)
{

double zsbndi;
int stepinbnd;
int nx, ny;

nx = Param.nx;
ny = Param.ny;
//update sl bnd


if (totaltime >= slbndtime)
for (int i = 1; i++; i < slbnd.size())
{
// find next timestep
stepinbnd = i;
double difft = slbnd[i].time - totaltime;
if (difft >= 0.0)
{
zsbndi = interptime(slbnd[i].wlev, slbnd[i - 1].wlev, slbnd[i].time - slbnd[i - 1].time, totaltime - slbnd[i - 1].time);
break;
}
}

//std::cout << "i= " << stepinbnd << "; " << zsbndi << "\n" << std::endl;



zsbndold = zsbndnew;
rtsl = slbndtime;
fscanf(fsl, "%f\t%f", &slbndtime, &zsbndnew);
//if (totaltime >= slbndtime)
//{

// zsbndold = zsbndnew;
// rtsl = slbndtime;
// fscanf(fsl, "%f\t%f", &slbndtime, &zsbndnew);
//slbndtime=+rtsl;
//zsbnd=zsbndold+(t-slbndtime+rtsl)*(zsbndnew-zsbndold)/rtsl;
}
//}



Expand All @@ -618,7 +635,7 @@ void flowbnd(XBGPUParam Param)
{
for (int ni = 0; ni < ny; ni++)
{
zsbnd[ni] = zsbndold + ((float) totaltime - rtsl)*(zsbndnew - zsbndold) / (slbndtime - rtsl);
zsbnd[ni] = zsbndi;//zsbndold + ((float) totaltime - rtsl)*(zsbndnew - zsbndold) / (slbndtime - rtsl);
}
}

Expand All @@ -630,13 +647,13 @@ void flowbnd(XBGPUParam Param)
dim3 gridDim(ceil((nx*1.0f) / blockDim.x), ceil((ny*1.0f) / blockDim.y), 1);
// FLow abs_2d should be here not at the flow step
// Set weakly reflective offshore boundary
ubnd1D << <gridDim, blockDim, 0 >> >(nx, ny, Param.dx, Param.dt, Param.g, Param.rho, (float)totaltime, wavbndtime, dtwavbnd, slbndtime, rtsl, zsbndold, zsbndnew, Trep, qbndold_g, qbndnew_g, zs_g, uu_g, vv_g, vu_g, umeanbnd_g, vmeanbnd_g, zb_g, cg_g, hum_g, cfm_g, Fx_g, hh_g);
ubnd1D << <gridDim, blockDim, 0 >> >(nx, ny, Param.dx, Param.dt, Param.g, Param.rho, (float)totaltime, wavbndtime, dtwavbnd, zsbndi, Trep, qbndold_g, qbndnew_g, zs_g, uu_g, vv_g, vu_g, umeanbnd_g, vmeanbnd_g, zb_g, cg_g, hum_g, cfm_g, Fx_g, hh_g);
//CUT_CHECK_ERROR("ubnd execution failed\n");
CUDA_CHECK(cudaThreadSynchronize());
}
else
{
ubndCPU(nx, ny, Param.dx, Param.dt, Param.g, Param.rho, (float)totaltime, wavbndtime, dtwavbnd, slbndtime, rtsl, zsbndold, zsbndnew, Trep, qbndold_g, qbndnew_g, zs_g, uu_g, vv_g, vu_g, umeanbnd_g, vmeanbnd_g, zb_g, cg_g, hum_g, cfm_g, Fx_g, hh_g);
ubndCPU(nx, ny, Param.dx, Param.dt, Param.g, Param.rho, (float)totaltime, wavbndtime, dtwavbnd, zsbndi, Trep, qbndold_g, qbndnew_g, zs_g, uu_g, vv_g, vu_g, umeanbnd_g, vmeanbnd_g, zb_g, cg_g, hum_g, cfm_g, Fx_g, hh_g);

}
}
Expand Down Expand Up @@ -1205,6 +1222,8 @@ int main(int argc, char **argv)

XBGPUParam XParam;

std::vector<SLBnd> slbnd;

std::ifstream fs("XBG_param.txt");

if (fs.fail()){
Expand Down Expand Up @@ -1327,14 +1346,20 @@ int main(int argc, char **argv)

printf("Opening sea level bnd...");

//std::vector<SLBnd> readWLfile(std::string WLfilename);
slbnd = readWLfile(XParam.slbnd);

fsl = fopen(XParam.slbnd.c_str(), "r");

fscanf(fsl, "%f\t%f", &rtsl, &zsbndold);
//fsl = fopen(XParam.slbnd.c_str(), "r");

//fscanf(fsl, "%f\t%f", &rtsl, &zsbndold);
rtsl = slbnd[0].time;
zsbndold = slbnd[0].wlev;

//Note: the first rtsl should be 0
fscanf(fsl, "%f\t%f", &slbndtime, &zsbndnew);
//fscanf(fsl, "%f\t%f", &slbndtime, &zsbndnew);
slbndtime = slbnd[1].time;
zsbndnew = slbnd[1].wlev;

printf("done\n");

//zsbnd sea level in bnd file
Expand Down Expand Up @@ -2060,7 +2085,7 @@ int main(int argc, char **argv)
// Calculate initial maximum timestep


FLOWDT << <gridDim, blockDim, 0 >> >(nx, ny, dx, 0.25f, dtflow_g, hh_g);
FLOWDT << <gridDim, blockDim, 0 >> >(nx, ny, dx, 0.5*XParam.CFL, dtflow_g, hh_g);
CUDA_CHECK(cudaThreadSynchronize());


Expand Down Expand Up @@ -2088,9 +2113,9 @@ int main(int argc, char **argv)



dt = arrmin[0];
dt = arrmin[0]*0.5;

printf("Initial timestep: dt=%f\n", arrmin[0]);
printf("Initial timestep: dt=%f\n", dt);


double tiny = 0.00000001;
Expand Down Expand Up @@ -2170,18 +2195,18 @@ int main(int argc, char **argv)
// Run the model
if (XParam.GPUDEVICE >= 0)
{
mainloopGPU(XParam);
mainloopGPU(XParam, slbnd);
}
else
{
mainloopCPU(XParam);
mainloopCPU(XParam, slbnd);
}




//close the bnd files and clean up a bit
fclose(fsl);
//fclose(fsl);
fclose(fwind);

endcputime = clock();
Expand Down
6 changes: 4 additions & 2 deletions XBeachGPU.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -134,7 +134,7 @@ extern "C" int signC(DECNUM x);
extern "C" void set_bndCPU(int nx, int ny,DECNUM Trep,int ntheta,DECNUM * theta,DECNUM *&sigm);
extern "C" void dispersion_initCPU(int nx,int ny,DECNUM twopi,DECNUM g,DECNUM aphi,DECNUM bphi,DECNUM * sigm,DECNUM * hh,DECNUM * &cg);
extern "C" void updatezomCPU(int nx, int ny,DECNUM cf,DECNUM cf2,DECNUM fw,DECNUM fw2,DECNUM * structdepth, DECNUM * &cfm,DECNUM * &fwm);
extern "C" void ubndCPU(int nx, int ny, DECNUM dx, DECNUM dt,DECNUM g, DECNUM rho,DECNUM totaltime,DECNUM wavbndtime,DECNUM rt,DECNUM slbndtime, DECNUM rtsl,DECNUM zsbndold,DECNUM zsbndnew,DECNUM Trep,DECNUM * qbndold, DECNUM * qbndnew,DECNUM *&zs, DECNUM * &uu,DECNUM * &vv, DECNUM *vu, DECNUM * umean, DECNUM * vmean,DECNUM * zb,DECNUM * cg,DECNUM * hum, DECNUM * zo, DECNUM *Fx,DECNUM *&hh);
extern "C" void ubndCPU(int nx, int ny, DECNUM dx, DECNUM dt,DECNUM g, DECNUM rho,DECNUM totaltime,DECNUM wavbndtime,DECNUM rt,DECNUM zsbnd,DECNUM Trep,DECNUM * qbndold, DECNUM * qbndnew,DECNUM *&zs, DECNUM * &uu,DECNUM * &vv, DECNUM *vu, DECNUM * umean, DECNUM * vmean,DECNUM * zb,DECNUM * cg,DECNUM * hum, DECNUM * zo, DECNUM *Fx,DECNUM *&hh);
extern "C" void offshorebndWavCPU(int nx, int ny,int ntheta,DECNUM totaltime,DECNUM Trep,DECNUM *St,DECNUM *&sigm, DECNUM *&ee);
extern "C" void sanityCPU(int nx, int ny,DECNUM eps,DECNUM * hh,DECNUM * sigm, int ntheta,DECNUM * &ee);
extern "C" void dispersionCPU(int nx,int ny,DECNUM twopi,DECNUM g,DECNUM aphi,DECNUM bphi,DECNUM * sigm,DECNUM * hh,DECNUM * &k,DECNUM * &c,DECNUM * &kh,DECNUM * &sinh2kh,DECNUM * &cg);
Expand Down Expand Up @@ -176,6 +176,8 @@ void split(const std::string &s, char delim, std::vector<std::string> &elems);
std::vector<std::string> split(const std::string &s, char delim);
std::string trim(const std::string& str, const std::string& whitespace);
XBGPUParam checkparamsanity(XBGPUParam param);
std::vector<SLBnd> readWLfile(std::string WLfilename);
double interptime(double next, double prev, double timenext, double time);

// General functions
//void CUDA_CHECK(cudaError CUDerr);
Expand All @@ -188,7 +190,7 @@ template <class T> const T& max (const T& a, const T& b);
extern "C"
void waveinitGPU(XBGPUParam Param);
void wavebnd(XBGPUParam Param);
void flowbnd(XBGPUParam Param);
void flowbnd(XBGPUParam Param, std::vector<SLBnd> slbnd);
void wavestep(XBGPUParam Param);
void flowstep(XBGPUParam Param);
void sedimentstep(XBGPUParam Param);
Expand Down
25 changes: 24 additions & 1 deletion read_input.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -141,6 +141,7 @@ std::vector<SLBnd> readWLfile(std::string WLfilename)
}

std::string line;
std::vector<std::string> lineelements;
SLBnd slbndline;
while (std::getline(fs, line))
{
Expand All @@ -152,14 +153,22 @@ std::vector<SLBnd> readWLfile(std::string WLfilename)
//Data should be in teh format :
//BASIN,CY,YYYYMMDDHH,TECHNUM/MIN,TECH,TAU,LatN/S,LonE/W,VMAX,MSLP,TY,RAD,WINDCODE,RAD1,RAD2,RAD3,RAD4,RADP,RRP,MRD,GUSTS,EYE,SUBREGION,MAXSEAS,INITIALS,DIR,SPEED,STORMNAME,DEPTH,SEAS,SEASCODE,SEAS1,SEAS2,SEAS3,SEAS4,USERDEFINED,userdata

slbndline = readBSHline(line);
//by default we expect tab delimitation
lineelements = split(line, '\t');
slbndline.time = std::stod(lineelements[0]);
slbndline.wlev = std::stod(lineelements[1]);

//slbndline = readBSHline(line);
slbnd.push_back(slbndline);
//std::cout << line << std::endl;
}

}
fs.close();

//std::cout << slbnd[0].wlev << std::endl;


return slbnd;
}

Expand Down Expand Up @@ -539,6 +548,15 @@ XBGPUParam readparamstr(std::string line, XBGPUParam param)
///////////////////////////////////////////////////////
// Input and output files
//

parameterstr = "outfile =";
parametervalue = findparameter(parameterstr, line);
if (!parametervalue.empty())
{
param.outfile = parametervalue;
//std::cerr << "Bathymetry file found!" << std::endl;
}

parameterstr = "SedThkfile =";
parametervalue = findparameter(parameterstr, line);
if (!parametervalue.empty())
Expand Down Expand Up @@ -701,6 +719,11 @@ std::string trim(const std::string& str, const std::string& whitespace)

return str.substr(strBegin, strRange);
}

double interptime(double next, double prev, double timenext, double time)
{
return prev + (time) / (timenext)*(next - prev);
}
/*
extern "C" void readbathy(int &nx, int &ny, float &dx, float &grdalpha, float *&zb)
{
Expand Down

0 comments on commit 3601e0c

Please sign in to comment.