Implement flexible output Times

.github/workflows/main.yml

@@ -105,6 +105,13 @@ jobs:
       - name: running test 14, Test on AOI implementation aa wall bnb
         run: |
           sed -i 's/test = 13/test = 14/' BG_param.txt
-          ./BG_Flood  
+          ./BG_Flood   
+      - name: running test 15, Test on flexible time input
+        run: |
+          sed -i 's/test = 14/test = 15/' BG_param.txt
+          ./BG_Flood
+          ./BG_Flood BG_param_test15.txt
+       ##   var=$(ncdump -v h_P0 Test15_zoom2.nc)
+       ## if: echo ${{ ${#var} == 7 }}
 
 

src/Arrays.h

@@ -92,6 +92,8 @@ struct outzoneB
 	std::string outname; // name for the output file (one for each zone)
 	int maxlevel; // maximum level in the zone
 	int minlevel; //minimum level in the zone
+	std::vector<double> OutputT; //Next time for the output of this zone
+	int index_next_OutputT = 0; //Index of next time output
 };
 
 
@@ -192,7 +194,7 @@ struct Model
 	std::map<std::string, std::string> Outvarlongname;
 	std::map<std::string, std::string> Outvarstdname;
 	std::map<std::string, std::string> Outvarunits;
-
+	std::vector<double> OutputT;
 
 	//other output
 	//std::vector< std::vector< Pointout > > TSallout;
@@ -229,6 +231,8 @@ struct Loop
 	int nstep = 0;
 	//useful for calculating avg timestep
 	int nstepout = 0;
+	// Needed to identify next output time
+	int indNextoutputtime = 0;
 
 	// usefull for Time series output
 	int nTSsteps = 0;

src/InitialConditions.cu

@@ -17,6 +17,7 @@
 
 
 #include "InitialConditions.h"
+#include "Input.h"
 
 template <class T> void InitialConditions(Param &XParam, Forcing<float> &XForcing, Model<T> &XModel)
 {
@@ -100,14 +101,14 @@ template <class T> void InitialConditions(Param &XParam, Forcing<float> &XForcin
 
 		// Initialise infiltration to IL where h is already wet
 		initinfiltration(XParam, XModel.blocks, XModel.evolv.h, XModel.il, XModel.hgw);
-
 	}
 
-
 	//=====================================
 	// Initialize output variables
 	initoutput(XParam, XModel);
 
+	// Initialise Output times' vector
+	initOutputTimes(XParam, XModel.OutputT, XModel.blocks);
 }
 template void InitialConditions<float>(Param &XParam, Forcing<float> &XForcing, Model<float> &XModel);
 template void InitialConditions<double>(Param &XParam, Forcing<float> &XForcing, Model<double> &XModel);
@@ -835,6 +836,7 @@ template <class T> void Initoutzone(Param& XParam, BlockP<T>& XBlock)
 		XzoneB.nblk = XParam.nblk;
 		XzoneB.maxlevel = XParam.maxlevel;
 		XzoneB.minlevel = XParam.minlevel;
+		XzoneB.OutputT = { XParam.totaltime, XParam.endtime };
 		AllocateCPU(XParam.nblk, 1, XzoneB.blk);
 		int I = 0;
 		for (int ib = 0; ib < XParam.nblk; ib++)
@@ -1329,8 +1331,7 @@ template <class T> __global__ void calcactiveCellGPU(Param XParam, BlockP<T> XBl
 	}
 }
 
-template <class T>
-void initinfiltration(Param XParam, BlockP<T> XBlock, T* h, T* initLoss ,T* hgw)
+template <class T> void initinfiltration(Param XParam, BlockP<T> XBlock, T* h, T* initLoss ,T* hgw)
 {
 //Initialisation to 0 (cold or hot start)
 
@@ -1354,3 +1355,89 @@ void initinfiltration(Param XParam, BlockP<T> XBlock, T* h, T* initLoss ,T* hgw)
 		}
 	}
 }
+
+// Create a vector of times steps from the input structure Toutput
+std::vector<double> GetTimeOutput(T_output time_info)
+{
+	//std::vector<double> time_vect;
+	//double time;
+
+	//Add independant values
+	//if (!time_info.val.empty())
+	//{
+	//	time_vect = time_info.val;
+	//}
+
+	//Add timesteps from the vector 
+	//time = time_info.init;
+	//while (time < time_info.end)
+	//{
+	//	time_vect.push_back(time);
+	//	time += time_info.tstep;
+	//}
+
+	//Add last timesteps from the vector definition
+	//time_vect.push_back(time_info.end);
+
+	return(time_info.val);
+}
+
+
+
+// Creation of a vector for times requiering a map output
+// Compilations of vectors and independent times from the general input
+// and the different zones outputs
+template <class T> void initOutputTimes(Param XParam, std::vector<double>& OutputT, BlockP<T>& XBlock)
+{
+	std::vector<double> times;
+	std::vector<double> times_partial;
+
+	times_partial = GetTimeOutput(XParam.Toutput);
+	//printf("Time partial:\n");
+	//for (int k = 0; k < times_partial.size(); k++)
+	//{
+	//	printf("%f, ", times_partial[k]);
+	//}
+	//printf("\n");
+
+	times.insert(times.end(), times_partial.begin(), times_partial.end());
+
+	// if zoneOutputs, add their contribution
+	if (XParam.outzone.size() > 0)
+	{
+		for (int ii = 0; ii < XParam.outzone.size(); ii++)
+		{
+			times_partial = GetTimeOutput(XParam.outzone[ii].Toutput);
+
+			//Add to main vector
+			times.insert(times.end(), times_partial.begin(), times_partial.end());
+			//Sort and remove duplicate before saving in outZone struct
+			std::sort(times_partial.begin(), times_partial.end());
+			times_partial.erase(unique(times_partial.begin(), times_partial.end()), times_partial.end());
+
+			XBlock.outZone[ii].OutputT = times_partial;
+		}
+	}
+	else //If not zoneoutput, output zone saved in zoneoutput structure
+	{
+		std::sort(times_partial.begin(), times_partial.end());
+		times_partial.erase(unique(times_partial.begin(), times_partial.end()), times_partial.end());
+
+		XBlock.outZone[0].OutputT = times_partial;
+	}
+
+	// Sort the times for output
+	std::sort(times.begin(), times.end());
+	times.erase(unique(times.begin(), times.end()), times.end());
+
+	printf("Output Times:\n");
+	for (int k = 0; k < times.size(); k++)
+	{
+		printf("%e, ", times[k]);
+	}
+	printf("\n");
+
+
+
+	OutputT = times;
+}

src/InitialConditions.h

@@ -33,5 +33,6 @@ template <class T> void InitzbgradientGPU(Param XParam, Model<T> XModel);
 
 template <class T> void calcactiveCellCPU(Param XParam, BlockP<T> XBlock, Forcing<float>& XForcing, T* zb);
 
+template <class T> void initOutputTimes(Param XParam, std::vector<double>& OutputT, BlockP<T>& XBlock);
 // End of global definition;
 #endif

src/Input.h

@@ -12,13 +12,24 @@ class TSoutnode {
 	std::string outname;
 };
 
+// Flexible definition of time outputs
+class T_output {
+public: 
+	double init = NAN;
+	double tstep = NAN;
+	double end = NAN;
+	std::vector<std::string> inputstr;
+	std::vector<double> val;
+};
+
 // Special output zones for nc files, informatin given by the user
 class outzoneP {
 public:
 	//std::vector<int> blocks; // one zone will spread across multiple blocks (entire blocks containing a part of the area will be output)
 	double xstart, xend, ystart, yend; // definition of the zone needed for special nc output (rectangular zone) by the user
 	//double xo, xmax, yo, ymax; // Real zone for output (because we output full blocks)
 	std::string outname; // name for the output file (one for each zone)
+	T_output Toutput; // time for outputs for the zone
 };
 
 class Flowin {

src/Mainloop.cu

@@ -12,11 +12,12 @@ template <class T> void MainLoop(Param &XParam, Forcing<float> XForcing, Model<T
 	//Define some useful variables 
 	Initmeanmax(XParam, XLoop, XModel, XModel_g);
 
-
-
+	// Check for map output (output initialisation if needed)
+	mapoutput(XParam, XLoop, XModel, XModel_g);
 
 	log("\t\tCompleted");
 	log("Model Running...");
+
 	while (XLoop.totaltime < XParam.endtime)
 	{
 		// Bnd stuff here
@@ -170,7 +171,7 @@ template <class T> Loop<T> InitLoop(Param &XParam, Model<T> &XModel)
 	Loop<T> XLoop;
 	XLoop.atmpuni = T(XParam.Paref);
 	XLoop.totaltime = XParam.totaltime;
-	XLoop.nextoutputtime = XParam.totaltime + XParam.outputtimestep;
+	XLoop.nextoutputtime = XModel.OutputT[0];
 
 	// Prepare output files
 	InitSave2Netcdf(XParam, XModel);
@@ -230,32 +231,57 @@ template <class T> void updateBnd(Param XParam, Loop<T> XLoop, Forcing<float> XF
 
 
 
-template <class T> void mapoutput(Param XParam, Loop<T> &XLoop,Model<T> XModel, Model<T> XModel_g)
+template <class T> void mapoutput(Param XParam, Loop<T> &XLoop, Model<T>& XModel, Model<T> XModel_g)
 {
 	XLoop.nstepout++;
-
-	if (XLoop.nextoutputtime - XLoop.totaltime <= XLoop.dt * T(0.5) && XParam.outputtimestep > 0.0)
+	double tiny = 0.0000001;
+	/*
+	if  (abs(XLoop.nextoutputtime - XParam.totaltime) < tiny)
 	{
-		char buffer[256];
-		sprintf(buffer, "%e", XParam.outputtimestep / XLoop.nstepout);
-		std::string str(buffer);
+		if (XParam.GPUDEVICE >= 0)
+		{
+			for (int ivar = 0; ivar < XParam.outvars.size(); ivar++)
+			{
+				CUDA_CHECK(cudaMemcpy(XModel.OutputVarMap[XParam.outvars[ivar]], XModel_g.OutputVarMap[XParam.outvars[ivar]], XParam.nblkmem * XParam.blksize * sizeof(T), cudaMemcpyDeviceToHost));
+			}
+		}
+
+		SaveInitialisation2Netcdf(XParam, XModel);
+
+		XLoop.indNextoutputtime++;
+		if (XLoop.indNextoutputtime < XModel.OutputT.size())
+		{
+			XLoop.nextoutputtime = min(XModel.OutputT[XLoop.indNextoutputtime], XParam.endtime);
+
+		}
 
-		log("Output to map. Totaltime = "+ std::to_string(XLoop.totaltime) +" s; Mean dt = " + str + " s");
+		XLoop.nstepout = 0;
+
+	}
+	*/
+	if (XLoop.nextoutputtime - XLoop.totaltime <= XLoop.dt * tiny)
+	{
 		if (XParam.GPUDEVICE >= 0)
 		{
 			for (int ivar = 0; ivar < XParam.outvars.size(); ivar++)
 			{
 				CUDA_CHECK(cudaMemcpy(XModel.OutputVarMap[XParam.outvars[ivar]], XModel_g.OutputVarMap[XParam.outvars[ivar]], XParam.nblkmem * XParam.blksize * sizeof(T), cudaMemcpyDeviceToHost));
 			}
 		}
-		
+
 		Save2Netcdf(XParam, XLoop, XModel);
 
+		XLoop.indNextoutputtime++;
+		if (XLoop.indNextoutputtime < XModel.OutputT.size())
+		{
+			XLoop.nextoutputtime = min(XModel.OutputT[XLoop.indNextoutputtime], XParam.endtime);
 
-		XLoop.nextoutputtime = min(XLoop.nextoutputtime + XParam.outputtimestep, XParam.endtime);
+		}
 
 		XLoop.nstepout = 0;
 	}
+
+
 }
 
 template <class T> void pointoutputstep(Param XParam, Loop<T> &XLoop, Model<T> XModel, Model<T> XModel_g)