Implement Groundwater module

Makefile

@@ -238,7 +238,7 @@ ifeq ($(SAMPLE_ENABLED),0)
 EXEC ?= @echo "[@]"
 endif
 
-OBJ = BG_Flood.o AdaptCriteria.o Adaptation.o Advection.o  Boundary.o ConserveElevation.o FlowCPU.o FlowGPU.o Friction.o Gradients.o GridManip.o Halo.o InitEvolv.o InitialConditions.o Kurganov.o Mainloop.o Meanmax.o MemManagement.o Mesh.o ReadForcing.o ReadInput.o Read_netcdf.o Reimann.o Setup_GPU.o Testing.o Updateforcing.o Util_CPU.o Write_netcdf.o utctime.o Write_txtlog.o FlowMLGPU.o Multilayer.o
+OBJ = BG_Flood.o AdaptCriteria.o Adaptation.o Advection.o  Boundary.o ConserveElevation.o FlowCPU.o FlowGPU.o Friction.o Gradients.o GridManip.o Halo.o InitEvolv.o InitialConditions.o Kurganov.o Mainloop.o Meanmax.o MemManagement.o Mesh.o ReadForcing.o ReadInput.o Read_netcdf.o Reimann.o Setup_GPU.o Testing.o Updateforcing.o Util_CPU.o Write_netcdf.o utctime.o Write_txtlog.o FlowMLGPU.o Multilayer.o groundwater.o
 
 
 ################################################################################
@@ -361,7 +361,10 @@ FlowMLGPU.o:./src/FlowMLGPU.cu
 Multilayer.o:./src/Multilayer.cu
 	$(EXEC) $(NVCC) $(INCLUDES) $(ALL_CCFLAGS) $(GENCODE_FLAGS) -o $@ -c $<
 
-BG_Flood: BG_Flood.o AdaptCriteria.o Adaptation.o Advection.o Boundary.o ConserveElevation.o FlowCPU.o FlowGPU.o Friction.o Gradients.o GridManip.o Halo.o InitEvolv.o InitialConditions.o Kurganov.o Mainloop.o Meanmax.o MemManagement.o Mesh.o ReadForcing.o ReadInput.o Read_netcdf.o Reimann.o Setup_GPU.o Testing.o Updateforcing.o Util_CPU.o Write_netcdf.o Write_txtlog.o Poly.o Spherical.o utctime.o FlowMLGPU.o Multilayer.o
+groundwater.o:./src/groundwater.cu
+	$(EXEC) $(NVCC) $(INCLUDES) $(ALL_CCFLAGS) $(GENCODE_FLAGS) -o $@ -c $<
+
+BG_Flood: BG_Flood.o AdaptCriteria.o Adaptation.o Advection.o Boundary.o ConserveElevation.o FlowCPU.o FlowGPU.o Friction.o Gradients.o GridManip.o Halo.o InitEvolv.o InitialConditions.o Kurganov.o Mainloop.o Meanmax.o MemManagement.o Mesh.o ReadForcing.o ReadInput.o Read_netcdf.o Reimann.o Setup_GPU.o Testing.o Updateforcing.o Util_CPU.o Write_netcdf.o Write_txtlog.o Poly.o Spherical.o utctime.o FlowMLGPU.o Multilayer.o groundwater.o
 
 	$(EXEC) $(NVCC) $(ALL_LDFLAGS) $(GENCODE_FLAGS) -o $@ $+ $(LIBRARIES)
 	$(EXEC) mkdir -p ../../bin/$(TARGET_ARCH)/$(TARGET_OS)/$(BUILD_TYPE)
@@ -372,7 +375,7 @@ run: build
 
 clean:
 
-	rm -f BG_Flood AdaptCriteria.o Adaptation.o Advection.o BG_Flood.o Boundary.o ConserveElevation.o FlowCPU.o FlowGPU.o Friction.o Gradients.o GridManip.o Halo.o InitEvolv.o InitialConditions.o Kurganov.o Mainloop.o Meanmax.o MemManagement.o Mesh.o ReadForcing.o ReadInput.o Read_netcdf.o Reimann.o Setup_GPU.o Testing.o Updateforcing.o Util_CPU.o Write_netcdf.o Write_txtlog.o Poly.o Spherical.o utctime.o FlowMLGPU.o Multilayer.o
+	rm -f BG_Flood AdaptCriteria.o Adaptation.o Advection.o BG_Flood.o Boundary.o ConserveElevation.o FlowCPU.o FlowGPU.o Friction.o Gradients.o GridManip.o Halo.o InitEvolv.o InitialConditions.o Kurganov.o Mainloop.o Meanmax.o MemManagement.o Mesh.o ReadForcing.o ReadInput.o Read_netcdf.o Reimann.o Setup_GPU.o Testing.o Updateforcing.o Util_CPU.o Write_netcdf.o Write_txtlog.o Poly.o Spherical.o utctime.o FlowMLGPU.o Multilayer.o groundwater.o
 
 	rm -rf ../../bin/$(TARGET_ARCH)/$(TARGET_OS)/$(BUILD_TYPE)/template
 

src/Arrays.h

@@ -174,7 +174,21 @@ struct AdvanceP
 	T* dhv;
 };
 
-
+template <class T>
+struct GoundWaterP
+{
+	// Groundwater Physical Parameters
+	T* zs;// groundwater celing elevation
+	T* h;// groundwater depth (z_gw - zaqb)
+	T* K;
+	T* fs;
+	T* Sy;
+	T* zb;// Elevation of aquifer bottom
+
+	// Groundwater Fluxes
+	T* Qx;
+	T* Qy;
+};
 struct outP
 {
 	float* z;
@@ -315,6 +329,8 @@ struct Model
 
 	//GroundWater elevation (due to the accumulation of water by infiltration during the simulation)
 	T* hgw;
+
+	GoundWaterP<T> gw;
 
 	// Used for external forcing too
 	// May need a better placeholder

src/FlowGPU.cu

@@ -1,4 +1,5 @@
 #include "FlowGPU.h"
+#include "groundwater.h"
 
 /**
  * @brief Main GPU flow solver for the flood model.
@@ -79,6 +80,8 @@ template <class T> void FlowGPU(Param XParam, Loop<T>& XLoop, Forcing<float> XFo
 	reset_var <<< gridDim, blockDim, 0 >>> (XParam.halowidth,XModel.blocks.active,XLoop.hugeposval,XModel.time.dtmax);
 	CUDA_CHECK(cudaDeviceSynchronize());
 
+
+
 	//============================================
 	// Calculate gradient for evolving parameters for predictor step
 	gradientGPUnew(XParam, XModel.blocks, XModel.evolv, XModel.grad, XModel.zb);
@@ -288,7 +291,11 @@ template <class T> void FlowGPU(Param XParam, Loop<T>& XLoop, Forcing<float> XFo
 		CUDA_CHECK(cudaDeviceSynchronize());
 	}
 
-	if (XParam.infiltration)
+	if (XParam.groundwater)
+	{
+		GroundwaterStepGPU(XParam, XLoop, XModel);
+	}
+	else if (XParam.infiltration)
 	{
 		AddinfiltrationImplicitGPU <<< gridDim, blockDim, 0 >>> (XParam, XLoop, XModel.blocks, XModel.il, XModel.cl, XModel.evolv, XModel.hgw);
 		CUDA_CHECK(cudaDeviceSynchronize());

src/Forcing.h

@@ -217,6 +217,36 @@ struct Forcing
 	Default: (see constant in parameters)
 	*/
 
+	StaticForcingP<T> K_gw;
+	/*Hydraulic conductivity map (m/s)
+	Ex: K_gw=K_gw.nc?K;
+	Default: (see constant in parameters)
+	*/
+
+	StaticForcingP<T> fs_gw;
+	/*Saturated infiltration rate map (m/s)
+	Ex: fs_gw=fs_gw.nc?fs;
+	Default: (see constant in parameters)
+	*/
+
+	StaticForcingP<T> Sy_gw;
+	/*Specific yield map
+	Ex: Sy_gw=Sy_gw.nc?Sy;
+	Default: (see constant in parameters)
+	*/
+
+	StaticForcingP<T> zb_gw;
+	/*Aquifer depth map (m)
+	Ex: Aquifer_Depth=Aquifer_Depth.nc?depth;
+	Default: (see constant in parameters)
+	*/
+
+	StaticForcingP<T> zs_gw_init;
+	/*Initial groundwater elevation map (m)
+	Ex: zgw_init=zgw_init.nc?zgw;
+	Default: (see constant in parameters)
+	*/
+
 	std::vector<StaticForcingP<int>> targetadapt;
 
 	std::vector<deformmap<T>> deform;

src/InitEvolv.cu

@@ -548,6 +548,34 @@ int AddZSoffset(Param XParam, BlockP<T> XBlock, EvolvingP<T> &XEv, T*zb)
 	return success;
 }
 
+template <class T>
+int Inith(Param XParam, BlockP<T> XBlock, T* h_gw,T* zs_gw,T* zb_gw)
+{
+	int success = 1;
+	int ib;
+	for (int ibl = 0; ibl < XParam.nblk; ibl++)
+	{
+		ib = XBlock.active[ibl];
+		for (int j = 0; j < XParam.blkwidth; j++)
+		{
+			for (int i = 0; i < XParam.blkwidth; i++)
+			{
+				int n = memloc(XParam, i, j, ib);
+
+				zs_gw[n] = max(zs_gw[n], zb_gw[n]);
+
+				h_gw[n] = utils::max(zs_gw[n] - zb_gw[n], T(0.0));
+
+			}
+
+		}
+	}
+
+	return success;
+}
+template int Inith<float>(Param XParam, BlockP<float> XBlock, float* h_gw, float* zs_gw, float* zb_gw);
+template int Inith<double>(Param XParam, BlockP<double> XBlock, double* h_gw, double* zs_gw, double* zb_gw);
+
 
 /**
  * @brief Read BG_Flood hotstart file and extract block attributes.

src/InitEvolv.h

@@ -22,5 +22,7 @@ template <class T> int AddZSoffset(Param XParam, BlockP<T> XBlock, EvolvingP<T>&
 
 template <class T> int readhotstartfile(Param XParam, BlockP<T> XBlock, EvolvingP<T>& XEv, T*& zb);
 
+template <class T> int Inith(Param XParam, BlockP<T> XBlock, T* h_gw, T* zs_gw, T* zb_gw);
+
 // End of global definition;
 #endif

src/InitialConditions.cu

@@ -113,6 +113,55 @@ template <class T> void InitialConditions(Param &XParam, Forcing<float> &XForcin
 		initinfiltration(XParam, XModel.blocks, XModel.evolv.h, XModel.il, XModel.hgw);
 	}
 
+	//=====================================
+	// Initialise groundwater parameters map
+	if (XParam.groundwater)
+	{
+		log("\tInitializing groundwater parameters");
+		if (!XForcing.K_gw.inputfile.empty())
+			interp2BUQ(XParam, XModel.blocks, XForcing.K_gw, XModel.gw.K);
+		else
+			InitArrayBUQ(XParam, XModel.blocks, (T)XParam.K_gw, XModel.gw.K);
+
+		if (!XForcing.fs_gw.inputfile.empty())
+			interp2BUQ(XParam, XModel.blocks, XForcing.fs_gw, XModel.gw.fs);
+		else
+			InitArrayBUQ(XParam, XModel.blocks, (T)XParam.fs_gw, XModel.gw.fs);
+
+		if (!XForcing.Sy_gw.inputfile.empty())
+			interp2BUQ(XParam, XModel.blocks, XForcing.Sy_gw, XModel.gw.Sy);
+		else
+			InitArrayBUQ(XParam, XModel.blocks, (T)XParam.Sy_gw, XModel.gw.Sy);
+
+		if (!XForcing.zb_gw.inputfile.empty())
+			interp2BUQ(XParam, XModel.blocks, XForcing.zb_gw, XModel.gw.zb);
+		else
+			InitArrayBUQ(XParam, XModel.blocks, (T)XParam.Aquifer_Depth, XModel.gw.zb);
+
+		if (!XForcing.zs_gw_init.inputfile.empty())
+			interp2BUQ(XParam, XModel.blocks, XForcing.zs_gw_init, XModel.gw.zs);
+		else if (!std::isnan(XParam.hgw_init))
+			InitArrayBUQ(XParam, XModel.blocks, (T)XParam.hgw_init, XModel.gw.zs);
+		else
+		{
+			// Default: hgw = zb (saturated to the surface)
+			Copy2CartCPU(XParam.nblk, XParam.blksize, XModel.gw.zs, XModel.zb);
+		}
+
+		Inith(XParam, XModel.blocks, XModel.gw.h, XModel.gw.zs, XModel.gw.zb);
+
+		// Initialise fluxes to 0
+		InitArrayBUQ(XParam, XModel.blocks, T(0.0), XModel.gw.Qx);
+		InitArrayBUQ(XParam, XModel.blocks, T(0.0), XModel.gw.Qy);
+
+		// Set edges for all groundwater parameters
+		setedges(XParam, XModel.blocks, XModel.gw.K);
+		setedges(XParam, XModel.blocks, XModel.gw.fs);
+		setedges(XParam, XModel.blocks, XModel.gw.Sy);
+		setedges(XParam, XModel.blocks, XModel.gw.zb);
+		setedges(XParam, XModel.blocks, XModel.gw.h);
+	}
+
 	//=====================================
 	// Initialize output variables
 	initoutput(XParam, XModel);
@@ -891,6 +940,31 @@ template<class T> void Initmaparray(Model<T>& XModel)
 	XModel.OutputVarMap["twet"] = XModel.wettime;
 	XModel.Outvarlongname["twet"] = "time since the cell has been wet";
 	XModel.Outvarunits["twet"] = "s";
+
+	XModel.OutputVarMap["K_gw"] = XModel.gw.K;
+	XModel.Outvarlongname["K_gw"] = "Groundwater hydraulic conductivity";
+	XModel.Outvarunits["K_gw"] = "m s-1";
+
+	XModel.OutputVarMap["fs_gw"] = XModel.gw.fs;
+	XModel.Outvarlongname["fs_gw"] = "Saturated infiltration rate";
+	XModel.Outvarunits["fs_gw"] = "m s-1";
+
+	XModel.OutputVarMap["Sy_gw"] = XModel.gw.Sy;
+	XModel.Outvarlongname["Sy_gw"] = "Specific yield";
+	XModel.Outvarunits["Sy_gw"] = "dimensionless";
+
+	XModel.OutputVarMap["zb_gw"] = XModel.gw.zb;
+	XModel.Outvarlongname["zb_gw"] = "Aquifer bed elevation above datum ";
+	XModel.Outvarunits["zb_gw"] = "m";
+
+	XModel.OutputVarMap["zs_gw"] = XModel.gw.zs;
+	XModel.Outvarlongname["zs_gw"] = "Aquifer water surface elevation above datum ";
+	XModel.Outvarunits["zs_gw"] = "m";
+
+	XModel.OutputVarMap["h_gw"] = XModel.gw.h;
+	XModel.Outvarlongname["h_gw"] = "Aquifer thickness ";
+	XModel.Outvarunits["h_gw"] = "m";
+
 	//XModel.OutputVarMap["vort"] = XModel.vort;
 }
 

src/MemManagement.cu

@@ -294,6 +294,18 @@ void AllocateCPU(int nblk, int blksize, Param XParam, Model<T>& XModel)
 		AllocateCPU(nblk, blksize, XModel.hgw);
 	}
 
+	if (XParam.groundwater)
+	{
+		AllocateCPU(nblk, blksize, XModel.gw.h);
+		AllocateCPU(nblk, blksize, XModel.gw.zs);
+		AllocateCPU(nblk, blksize, XModel.gw.K);
+		AllocateCPU(nblk, blksize, XModel.gw.fs);
+		AllocateCPU(nblk, blksize, XModel.gw.Sy);
+		AllocateCPU(nblk, blksize, XModel.gw.zb);
+		AllocateCPU(nblk, blksize, XModel.gw.Qx);
+		AllocateCPU(nblk, blksize, XModel.gw.Qy);
+	}
+
 	if (XParam.outmax)
 	{
 		AllocateCPU(nblk, blksize, XModel.evmax);
@@ -540,7 +552,18 @@ void ReallocArray(int nblk, int blksize, Param XParam, Model<T>& XModel)
 		ReallocArray(nblk, blksize, XModel.il);
 		ReallocArray(nblk, blksize, XModel.cl);
 		ReallocArray(nblk, blksize, XModel.hgw);
+	}
 
+	if (XParam.groundwater)
+	{
+		ReallocArray(nblk, blksize, XModel.gw.K);
+		ReallocArray(nblk, blksize, XModel.gw.fs);
+		ReallocArray(nblk, blksize, XModel.gw.Sy);
+		ReallocArray(nblk, blksize, XModel.gw.zb);
+		ReallocArray(nblk, blksize, XModel.gw.zs);
+		ReallocArray(nblk, blksize, XModel.gw.h);
+		ReallocArray(nblk, blksize, XModel.gw.Qx);
+		ReallocArray(nblk, blksize, XModel.gw.Qy);
 	}
 
 	if (XParam.outmax)
@@ -860,6 +883,18 @@ void AllocateGPU(int nblk, int blksize, Param XParam, Model<T>& XModel)
 		AllocateGPU(nblk, blksize, XModel.hgw);
 	}
 
+	if (XParam.groundwater)
+	{
+		AllocateGPU(nblk, blksize, XModel.gw.h);
+		AllocateGPU(nblk, blksize, XModel.gw.zs);
+		AllocateGPU(nblk, blksize, XModel.gw.K);
+		AllocateGPU(nblk, blksize, XModel.gw.fs);
+		AllocateGPU(nblk, blksize, XModel.gw.Sy);
+		AllocateGPU(nblk, blksize, XModel.gw.zb);
+		AllocateGPU(nblk, blksize, XModel.gw.Qx);
+		AllocateGPU(nblk, blksize, XModel.gw.Qy);
+	}
+
 	if (XParam.outmax)
 	{
 		AllocateGPU(nblk, blksize, XModel.evmax);

src/Param.h

@@ -24,10 +24,17 @@ class Param {
 	double Cd = 0.002; // Wind drag coefficient
 	double il = 0.0; //Initial Loss value (if constant)
 	double cl = 0.0; //Continuous Loss value (if constant)
+
+	double K_gw = 0.0; // Hydraulic conductivity (m/s)
+	double fs_gw = 0.0; // Saturated infiltration rate (m/s)
+	double Sy_gw = 0.05; // Specific yield
+	double Aquifer_Depth = 10.0; // Aquifer depth (m)
+	double hgw_init = nan(""); // Initial groundwater elevation (m)
 	bool windforcing = false; //not working yet
 	bool atmpforcing = false;
 	bool rainforcing = false;
 	bool infiltration = false;
+	bool groundwater = false;
 
 	bool conserveElevation = false; //Switch to force the conservation of zs instead of h at the interface between coarse and fine blocks
 	bool wetdryfix = true; // Switch to remove wet/dry instability (i.e. true reoves instability and false leaves the model as is)

src/ReadForcing.cu

@@ -381,6 +381,32 @@ void readforcing(Param & XParam, Forcing<T> & XForcing)
 		}
 	}
 
+	//======================
+	// groundwater file
+
+	if (XParam.groundwater)
+	{
+		//XForcing.K_gw = readforcinghead(XForcing.K_gw);
+		XForcing.K_gw.denanval = 0.0;
+		readstaticforcing(XForcing.K_gw);
+
+		XForcing.fs_gw.denanval = 0.0;
+		readstaticforcing(XForcing.fs_gw);
+
+		XForcing.Sy_gw.denanval = 0.0;
+		readstaticforcing(XForcing.Sy_gw);
+
+		XForcing.zb_gw.denanval = 0.0;
+		readstaticforcing(XForcing.zb_gw);
+
+		XForcing.zs_gw_init.denanval = 0.0;
+		readstaticforcing(XForcing.zs_gw_init);
+		//XForcing.fs_gw = readforcinghead(XForcing.fs_gw);
+		//XForcing.Sy_gw = readforcinghead(XForcing.Sy_gw);
+		//XForcing.Aquifer_Depth = readforcinghead(XForcing.Aquifer_Depth);
+		//XForcing.hgw_init = readforcinghead(XForcing.hgw_init);
+	}
+
 	//======================
 	// Polygon data
 	if (!XForcing.AOI.file.empty())
@@ -452,6 +478,7 @@ template <class T> void readstaticforcing(int step,T& Sforcing)
 		denan(Sforcing.nx, Sforcing.ny, float(Sforcing.denanval), Sforcing.val);
 
 	}
+
 	else
 	{
 		//Error message