Merge branch 'Test_new_BND_A100' of https://github.com/CyprienBossere…

…lle/BG_Flood into Test_new_BND_A100

.github/workflows/main.yml

@@ -91,7 +91,15 @@ jobs:
         run: |
           sed -i 's/test = 10/test = 11/' BG_param.txt
           ./BG_Flood
-          
 
-     
+      - name: running test 12, Test on calendar time
+        run: |
+          sed -i 's/test = 11/test = 12/' BG_param.txt
+          ./BG_Flood
 
+      - name: running test 13, Test on multi-Bathy / Multi-Roughness input
+        run: |
+          sed -i 's/test = 12/test = 13/' BG_param.txt
+          ./BG_Flood          
+
+     

README.md

@@ -11,5 +11,9 @@ The adaptive version of the model is currently under heavy development. Makefile
 
 [![Documentation (Work in progress)](https://img.shields.io/badge/doc-Doxygen-blue)](https://aliceharang.github.io/)
 
+## Reference
+
+Bosserelle C., Lane E., Harang A., (2021) BG-Flood: A GPU adaptive, open-source, general inundation hazard model. Proceedings of the Australasian Coasts & Ports 2021 Conference. 
+
 ## Testing
 CI test (development branch):![Test](https://github.com/CyprienBosserelle/BG_Flood/actions/workflows//main.yml/badge.svg?branch=development)

src/Adaptation.cu

@@ -72,7 +72,7 @@ template <class T> void Adaptation(Param& XParam, Forcing<float> XForcing, Model
 		//=====================================
 		// Initialise Friction map
 
-		if (!XForcing.cf.inputfile.empty())
+		if (!XForcing.cf.empty())
 		{
 			interp2BUQ(XParam, XModel.blocks, XForcing.cf, XModel.cf);
 		}
@@ -87,8 +87,22 @@ template <class T> void Adaptation(Param& XParam, Forcing<float> XForcing, Model
 		// Initialise the continuous losses map
 		if (XParam.infiltration)
 		{
-			interp2BUQ(XParam, XModel.blocks, XForcing.il, XModel.il);
-			interp2BUQ(XParam, XModel.blocks, XForcing.cl, XModel.cl);
+			if (!XForcing.il.inputfile.empty())
+			{
+				interp2BUQ(XParam, XModel.blocks, XForcing.il, XModel.il);
+			}
+			else
+			{
+				InitArrayBUQ(XParam, XModel.blocks, (T)XParam.il, XModel.il);
+			}
+			if (!XForcing.cl.inputfile.empty())
+			{
+				interp2BUQ(XParam, XModel.blocks, XForcing.cl, XModel.cl);
+			}
+			else
+			{
+				InitArrayBUQ(XParam, XModel.blocks, (T)XParam.cl, XModel.cl);
+			}
 			// Set edges of friction map
 			setedges(XParam, XModel.blocks, XModel.il);
 			setedges(XParam, XModel.blocks, XModel.cl);

src/Forcing.h

@@ -79,7 +79,6 @@ template <class T>
 struct StaticForcingP : public inputmap
 {
 	T *val;
-
 
 };
 
@@ -196,8 +195,9 @@ struct Forcing
 	Default: None but input NECESSARY
 	*/
 
-	StaticForcingP<T> cf;
+	std::vector<StaticForcingP<T>> cf;
 	/*Bottom friction coefficient map (associated to the chosen bottom friction model)
+	A list of roughness map can be provide. At any grid point, the last one defined will be used.
 	Ex: cf=0.001;
 	Ex: cf=bottom_friction.nc?bfc;
 	Default: (see constant in parameters)

src/InitialConditions.cu

@@ -32,7 +32,7 @@ template <class T> void InitialConditions(Param &XParam, Forcing<float> &XForcin
 	//=====================================
 	// Initialise Friction map
 
-	if (!XForcing.cf.inputfile.empty())
+	if (!XForcing.cf.empty())
 	{
 		interp2BUQ(XParam, XModel.blocks, XForcing.cf, XModel.cf);
 	}
@@ -77,8 +77,22 @@ template <class T> void InitialConditions(Param &XParam, Forcing<float> &XForcin
 
 	if (XParam.infiltration)
 	{
-		interp2BUQ(XParam, XModel.blocks, XForcing.il, XModel.il);
-		interp2BUQ(XParam, XModel.blocks, XForcing.cl, XModel.cl);
+		if (!XForcing.il.inputfile.empty())
+		{
+			interp2BUQ(XParam, XModel.blocks, XForcing.il, XModel.il);
+		}
+		else
+		{
+			InitArrayBUQ(XParam, XModel.blocks, (T)XParam.il, XModel.il);
+		}
+		if (!XForcing.cl.inputfile.empty())
+		{
+			interp2BUQ(XParam, XModel.blocks, XForcing.cl, XModel.cl);
+		}
+		else
+		{
+			InitArrayBUQ(XParam, XModel.blocks, (T)XParam.cl, XModel.cl);
+		}
 		// Set edges of friction map
 		setedges(XParam, XModel.blocks, XModel.il);
 		setedges(XParam, XModel.blocks, XModel.cl);

src/Mainloop.cu

@@ -486,7 +486,7 @@ template <class T> void CrashDetection(Param& XParam, Loop<T> XLoop, Model<T> XM
 		log(" \n ");
 		log("\t\tModel CRASHED: time steps (" + std::to_string(XLoop.dt) + ") inferior to " + std::to_string(XParam.dtmin) + "\n");
 
-		std::vector<std::string> outvi = { "zb","h","zs","u","v","ho", "vo", "uo", "zso" };
+		std::vector<std::string> outvi = { "zb","h","zs","u","v","ho","vo","uo","zso" };
 		std::vector<std::string> outvold = XParam.outvars;
 		if (XParam.GPUDEVICE >= 0)
 		{

src/Meanmax.cu

@@ -1,7 +1,7 @@
 
 #include "Meanmax.h"
 
-
+// Calculate the mean and or max of requested variable on the GPU/CPU
 template <class T> void Calcmeanmax(Param XParam, Loop<T>& XLoop, Model<T> XModel, Model<T> XModel_g)
 {
 	dim3 blockDim(XParam.blkwidth, XParam.blkwidth, 1);
@@ -520,4 +520,4 @@ template <class T> __host__ void addwettime_CPU(Param XParam, BlockP<T> XBlock,
 		}
 	}
 
-}
+}

src/Param.h

@@ -22,6 +22,8 @@ class Param {
 	int frictionmodel = 0; // Bottom friction model (-1: Manning model, 0: quadratic, 1: Smart model)
 	double cf = 0.0001; // Bottom friction coefficient for flow model (if constant)
 	double Cd = 0.002; // Wind drag coefficient
+	double il = 0.0; //Initial Loss (if constant)
+	double cl = 0.0; //Continuous Loss (if constant)
 	bool windforcing = false; //not working yet
 	bool atmpforcing = false;
 	bool rainforcing = false;

src/ReadForcing.cu

@@ -149,16 +149,30 @@ void readforcing(Param & XParam, Forcing<T> & XForcing)
 	//==================
 	// Friction maps 
 
-	if (!XForcing.cf.inputfile.empty())
+	if (!XForcing.cf.empty())
 	{
-		XForcing.cf.denanval = 0.0000001;
+		//denanval = 0.0000001;
 		log("\nRead Roughness map (cf) data...");
 		// roughness map was specified!
-		readstaticforcing(XForcing.cf);
-
+		//readstaticforcing(XForcing.cf);
 		//log("...done");
+		// Here we are not using the automated denaning because we want to preserve the Nan in all but the "main/first" listed roughness map. 
+		// This mean that subsequently listed roughness map can have large NAN holes in them.
+		for (int ib = 0; ib < XForcing.cf.size(); ib++)
+		{
+
+			readstaticforcing(XForcing.cf[ib]);
+			if (ib == 0) // Fill Nan for only the first map listed, the others will use values from original bathy topo.
+			{
+				denan(XForcing.cf[ib].nx, XForcing.cf[ib].ny, T(0.0000001), XForcing.cf[ib].val);
+			}
+		}
 	}
 
+
+
+
+
 	//==================
 	// Rain losses maps
 
@@ -546,6 +560,7 @@ std::string readCRSfrombathy(std::string crs_ref, StaticForcingP<float>& Sforcin
 	char* crs_wkt;
 	std::string crs_ref2;
 
+	crs_wkt = "";
 
 	if (!Sforcing.inputfile.empty())
 	{
@@ -633,7 +648,7 @@ std::string readCRSfrombathy(std::string crs_ref, StaticForcingP<float>& Sforcin
 			{
 				printf("CRS_info detected but not understood reverting to default CRS\n Rename attribute in grid-mapping variable\n");
 
-				crs_wkt = "";
+				//crs_wkt = ""; //Move to the top of the file for initialisation
 			}
 
 		}

src/ReadInput.cu

@@ -247,12 +247,34 @@ Param readparamstr(std::string line, Param param)
 	parametervalue = findparameter(paramvec, line);
 	if (!parametervalue.empty())
 	{
-		if (std::isdigit(parametervalue[0]) == true)
+		if (std::any_of(parametervalue.begin(), parametervalue.end(), ::isalpha) == false) //(std::isdigit(parametervalue[0]) == true)
 		{
 			param.cf = std::stod(parametervalue);
 		}
 	}
 
+	paramvec = { "il","Rain_il","initialloss" };
+	parametervalue = findparameter(paramvec, line);
+	if (!parametervalue.empty())
+	{
+		if (std::any_of(parametervalue.begin(), parametervalue.end(), ::isalpha) == false) //(std::isdigit(parametervalue[0]) == true)
+		{
+			param.il = std::stod(parametervalue);
+			param.infiltration = true;
+		}
+	}
+
+	paramvec = { "cl","Rain_cl","continuousloss" };
+	parametervalue = findparameter(paramvec, line);
+	if (!parametervalue.empty())
+	{
+		if (std::any_of(parametervalue.begin(), parametervalue.end(), ::isalpha) == false) //(std::isdigit(parametervalue[0]) == true)
+		{
+			param.cl = std::stod(parametervalue);
+			param.infiltration = true;
+		}
+	}
+
 	paramvec = { "VelThreshold","vthresh","vmax","velmax" };
 	parametervalue = findparameter(paramvec, line);
 	if (!parametervalue.empty())
@@ -424,7 +446,7 @@ Param readparamstr(std::string line, Param param)
 			//Verify that the variable name makes sense?
 			//Need to add more here
 
-			std::vector<std::string> SupportedVarNames = { "zb", "zs", "u", "v", "h", "hmean", "zsmean", "umean", "vmean", "hUmean", "Umean", "hmax", "zsmax", "umax", "vmax", "hUmax", "Umax", "twet", "dhdx","dhdy","dzsdx","dzsdy","dzbdx","dzbdy","dudx","dudy","dvdx","dvdy","Fhu","Fhv","Fqux","Fqvy","Fquy","Fqvx","Su","Sv","dh","dhu","dhv","cf", "Patm", "datmpdx", "datmpdy", "il", "cl", "hgw"};
+			std::vector<std::string> SupportedVarNames = { "zb","zs","u","v","h","hmean","zsmean","umean","vmean","hUmean","Umean","hmax","zsmax","umax","vmax","hUmax","Umax","twet","dhdx","dhdy","dzsdx","dzsdy","dzbdx","dzbdy","dudx","dudy","dvdx","dvdy","Fhu","Fhv","Fqux","Fqvy","Fquy","Fqvx","Su","Sv","dh","dhu","dhv","cf","Patm","datmpdx","datmpdy","il","cl","hgw"};
 
 			std::string vvar = trim(vars[nv], " ");
 			for (int isup = 0; isup < SupportedVarNames.size(); isup++)
@@ -778,8 +800,6 @@ Forcing<T> readparamstr(std::string line, Forcing<T> forcing)
 
 	paramvec = { "AOI","aoipoly" };
 	parametervalue = findparameter(paramvec, line);
-	//parameterstr = "bathy";
-	//parametervalue = findparameter(parameterstr, line);
 	if (!parametervalue.empty())
 	{
 		forcing.AOI.file= parametervalue;
@@ -951,18 +971,27 @@ Forcing<T> readparamstr(std::string line, Forcing<T> forcing)
 	parametervalue = findparameter(paramvec, line);
 	if (!parametervalue.empty())
 	{
-		if (std::isdigit(parametervalue[0]) == false)
+		if (std::any_of(parametervalue.begin(), parametervalue.end(), ::isalpha)) //(std::isdigit(parametervalue[0]) == false)
 		{
-			forcing.cf = readfileinfo(parametervalue, forcing.cf);
+			//forcing.cf = readfileinfo(parametervalue, forcing.cf);
+			StaticForcingP<float> infoRoughness;
+			forcing.cf.push_back(readfileinfo(parametervalue, infoRoughness));
 		}
 	}
 
+
+	//if (!parametervalue.empty())
+	//{
+	//
+		//std::cerr << "Bathymetry file found!" << std::endl;
+	//}
+
 	// Rain losses, initial and continuous loss
 	paramvec = { "il","Rain_il","initialloss" };
 	parametervalue = findparameter(paramvec, line);
 	if (!parametervalue.empty())
 	{
-		if (std::isdigit(parametervalue[0]) == false)
+		if (std::any_of(parametervalue.begin(), parametervalue.end(), ::isalpha)) //(std::isdigit(parametervalue[0]) == false)
 		{
 			forcing.il = readfileinfo(parametervalue, forcing.il);
 		}
@@ -971,7 +1000,7 @@ Forcing<T> readparamstr(std::string line, Forcing<T> forcing)
 	parametervalue = findparameter(paramvec, line);
 	if (!parametervalue.empty())
 	{
-		if (std::isdigit(parametervalue[0]) == false)
+		if (std::any_of(parametervalue.begin(), parametervalue.end(), ::isalpha)) //(std::isdigit(parametervalue[0]) == false)
 		{
 			forcing.cl = readfileinfo(parametervalue, forcing.cl);
 		}
@@ -1397,7 +1426,7 @@ void checkparamsanity(Param& XParam, Forcing<float>& XForcing)
 	}
 
 	//Check that the Initial Loss/ Continuing Loss model is used if il, cl or hgw output are asked by user.
-	if (XForcing.il.inputfile.empty() || XForcing.cl.inputfile.empty())
+	if (!XParam.infiltration) // (XForcing.il.inputfile.empty() && XForcing.cl.inputfile.empty() && (XParam.il == 0.0) && (XParam.cl == 0.0))
 	{
 		std::vector<std::string> namestr = { "il","cl","hgw"};
 		for (int ii = 0; ii < namestr.size(); ii++)

src/Read_netcdf.cu

@@ -279,7 +279,8 @@ void readgridncsize(const std::string ncfilestr, const std::string varstr, std::
 		//size_t start[] = { 0 };
 		//size_t count[] = { (size_t)nt };
 		//status = nc_get_vara_double(ncid, varid, start, count, ttempvar);
-		status = readnctime2(ncid, reftime, nt, ttempvar);
+
+		status = readnctime2(ncid, coordname, reftime, nt, ttempvar);
 
 		to = ttempvar[0];
 		tmax= ttempvar[nt-1];
@@ -426,7 +427,7 @@ int readnctime(std::string filename, double * &time)
 	return status;
 }
 
-int readnctime2(int ncid,std::string refdate,size_t nt, double*& time)
+int readnctime2(int ncid,char * timecoordname,std::string refdate,size_t nt, double*& time)
 {
 
 	int status, varid;
@@ -453,13 +454,13 @@ int readnctime2(int ncid,std::string refdate,size_t nt, double*& time)
 	///	//varstr = "time";
 	///}
 
-	// Warning this could be more robust by taking the unlimited dimention if time does not exist!
-	std::string Varname = "time";
+	// Warning this could be more robust by taking the unlimited dimension if time does not exist!
+	//std::string Varname = "time";
 
 	///status = nc_open(ncfilestr.c_str(), 0, &ncid);
 	///if (status != NC_NOERR) handle_ncerror(status);
 
-	status = nc_inq_varid(ncid, Varname.c_str(), &varid);
+	status = nc_inq_varid(ncid, timecoordname, &varid);
 	if (status != NC_NOERR) handle_ncerror(status);
 
 	// inquire unit of time

src/Read_netcdf.h

@@ -37,7 +37,7 @@ template<class T> void readgridncsize(T& Imap);
 
 int readvarinfo(std::string filename, std::string Varname, size_t *&ddimU);
 int readnctime(std::string filename, double * &time);
-int readnctime2(int ncid, std::string refdate, size_t nt, double*& time);
+int readnctime2(int ncid,char* timecoordname, std::string refdate, size_t nt, double*& time);
 template <class T> int readncslev1(std::string filename, std::string varstr, size_t indx, size_t indy, size_t indt, bool checkhh, double eps, T * &zsa);
 template <class T> int readvardata(std::string filename, std::string Varname, int step, T*& vardata, bool flipx, bool flipy);
 //template <class T> int readhotstartfile(Param XParam, int * leftblk, int *rightblk, int * topblk, int* botblk, double * blockxo, double * blockyo, T * &zs, T * &zb, T * &hh, T *&uu, T * &vv);