HYD_system: bug fixes in

-bilan gpu;
- mixed gpu-cpu FP

Author: dabachma

source_code/system_hydraulic/source_code/models/_Hyd_Model.cpp

@@ -265,8 +265,8 @@ void _Hyd_Model::init_solver_gpu(Hyd_Param_Global* global_params) {
 	pManager->setSelectedDevice(scheme_info.selected_device);							// Set GPU device to Use. Important: Has to be called after setExecutor. Default is the faster one.
 	pManager->setSimulationLength(simulationLength);									// Set Simulation Length
 	pManager->setOutputFrequency(outputFrequency);										// Set Output Frequency
-	pManager->setFloatPrecision(model::floatPrecision::kDouble);						// Set Precision
-	
+	pManager->setFloatPrecision(model::floatPrecision::kDouble);						// Set Precision 
+	//pManager->setFloatPrecision(model::floatPrecision::kSingle);
 	//Create the domain
 	CDomainCartesian* ourCartesianDomain = pManager->getDomain();
 	ourCartesianDomain->setCellResolution(*myFloodplain->Param_FP.get_ptr_width_x(), *myFloodplain->Param_FP.get_ptr_width_y());

source_code/system_hydraulic/source_code/models/coupling/Hyd_Coupling_Management.cpp

@@ -1421,7 +1421,9 @@ double Hyd_Coupling_Management::get_min_timestep_fp2fp(void) {
 
 
 		for (int i = 0; i < this->number_fp2fp; i++) {
-			min_t = min(this->coupling_fp2fp[i].list.get_min_timestep(), min_t);
+			if (this->coupling_fp2fp[i].list.get_min_timestep()>0) {
+				min_t = min(this->coupling_fp2fp[i].list.get_min_timestep(), min_t);
+			}
 		}
 
 

source_code/system_hydraulic/source_code/models/floodplain/elements/Hyd_Element_Floodplain_Type_Standard.cpp

@@ -254,15 +254,15 @@ void Hyd_Element_Floodplain_Type_Standard::calculate_v_out(const double sin_valu
 		this->v_x_out = (this->v_x + this->element_x_minus->element_type->get_flowvelocity_vx())*0.5;
 	}
 	else {
-		this->v_x_out = (this->v_x + 0.0)*0.5;
+		this->v_x_out = (this->v_x)*0.5;
 	}
 
 	//y-direction
 	if (this->element_y_minus != NULL) {
 		this->v_y_out = (this->v_y + this->element_y_minus->element_type->get_flowvelocity_vy())*0.5;
 	}
 	else {
-		this->v_y_out = (this->v_y + 0.0)*0.5;
+		this->v_y_out = (this->v_y)*0.5;
 	}
 
 	//turning
@@ -1358,6 +1358,7 @@ void Hyd_Element_Floodplain_Type_Standard::calculate_wet_dry_duration(const doub
 	if(this->was_wet_flag==false && this->h_value>wet_boundary){
 		this->was_wet_flag=true;
 	}
+
 	if(this->wet_flag==false && this->h_value>wet_boundary){
 		this->wet_flag=true;
 		this->time_wet_start=time_point;
@@ -1369,6 +1370,7 @@ void Hyd_Element_Floodplain_Type_Standard::calculate_wet_dry_duration(const doub
 		this->wet_duration=this->wet_duration+(time_point-this->time_wet_start);
 		this->time_wet_start=time_point;
 	}
+
 	//set the duration
 	if(this->wet_flag==true && this->h_value<wet_boundary){
 		this->wet_flag=false;

source_code/system_hydraulic/source_code/solvergpu/source_code/CSchemeGodunov.cpp

@@ -59,6 +59,7 @@ CSchemeGodunov::CSchemeGodunov(void)
 	this->ulReductionGlobalSize = 0;
 	this->ulReductionWorkgroupSize = 0;
 	this->dAvgTimestep = 0.0;
+	this->sum_deficit_bilan = 0.0;
 
 	this->dCurrentTime = 0.0;
 	this->dCurrentTimestepMovAvg = 0.0;
@@ -572,10 +573,14 @@ void CSchemeGodunov::prepare1OMemory(){
 	this->oclBufferCellManning = new COCLBuffer("Manning coefficients", oclModel, true, true);
 	if (this->bUseOptimizedBoundary == false) {
 		this->oclBufferCellBoundary = new COCLBuffer("Boundary Values", oclModel, false, true);
+		this->oclBufferBilanDef = new COCLBuffer("BilanDef", oclModel, false, true);
+		this->sum_deficit_bilan = 0.0;
 	}
 	else {
 		this->oclBufferCouplingIDs = new COCLBuffer("Coupling IDs", oclModel, true, true);
 		this->oclBufferCouplingValues = new COCLBuffer("Coupling Values", oclModel, false, true);
+		this->oclBufferBilanDef = new COCLBuffer("BilanDef", oclModel, false, true);
+		this->sum_deficit_bilan = 0.0;
 	}
 	this->oclBufferUsePoleni = new COCLBuffer("Poleni Booleans", oclModel, true, true);
 	this->oclBuffer_opt_zxmax = new COCLBuffer("opt_zxmax Values", oclModel, true, true);
@@ -584,15 +589,19 @@ void CSchemeGodunov::prepare1OMemory(){
 	this->oclBuffer_opt_cy = new COCLBuffer("opt_cy Values", oclModel, true, true);
 	this->oclBufferCellBed = new COCLBuffer("Bed elevations", oclModel, true, true);
 
+	
+
 	this->oclBufferCellStates->setPointer(pCellStates, ucFloatSize * 4 * pDomain->getCellCount());
 	this->oclBufferCellStatesAlt->setPointer(pCellStates, ucFloatSize * 4 * pDomain->getCellCount());
 	this->oclBufferCellManning->setPointer(pManningValues, ucFloatSize * pDomain->getCellCount());
 	if (this->bUseOptimizedBoundary == false) {
 		this->oclBufferCellBoundary->setPointer(pBoundaryValues, ucFloatSize * pDomain->getCellCount());
+		this->oclBufferBilanDef->setPointer(pBoundaryValues, ucFloatSize * pDomain->getCellCount());
 	}
 	else {
 		this->oclBufferCouplingIDs->setPointer(pCouplingIDs, sizeof(cl_ulong) * this->ulCouplingArraySize);
 		this->oclBufferCouplingValues->setPointer(pCouplingValues, ucFloatSize * this->ulCouplingArraySize);
+		this->oclBufferBilanDef->setPointer(pCouplingValues, ucFloatSize * this->ulCouplingArraySize);
 	}
 	this->oclBufferUsePoleni->setPointer(pPoleniValues, sizeof(sUsePoleni) * pDomain->getCellCount());
 	this->oclBuffer_opt_zxmax->setPointer(pOpt_zxmax, ucFloatSize * pDomain->getCellCount());
@@ -606,10 +615,12 @@ void CSchemeGodunov::prepare1OMemory(){
 	this->oclBufferCellManning->createBuffer();
 	if (this->bUseOptimizedBoundary == false) {
 		this->oclBufferCellBoundary->createBuffer();
+		this->oclBufferBilanDef->createBuffer();
 	}
 	else {
 		this->oclBufferCouplingIDs->createBuffer();
 		this->oclBufferCouplingValues->createBuffer();
+		this->oclBufferBilanDef->createBuffer();
 	}
 	this->oclBufferUsePoleni->createBuffer();
 	this->oclBuffer_opt_zxmax->createBuffer();
@@ -628,7 +639,7 @@ void CSchemeGodunov::prepare1OMemory(){
 	this->oclBufferTimeTarget = new COCLBuffer("Target time (sync)", oclModel, false, true, ucFloatSize, true);
 	this->oclBufferTimestepMovAvg = new COCLBuffer("Timestep Moving Average", oclModel, false, true, ucFloatSize, true);
 	this->oclBufferTimeHydrological = new COCLBuffer("Time (hydrological)", oclModel, false, true, ucFloatSize, true);
-	this->oclBufferBilanDef= new COCLBuffer("BilanDef", oclModel, false, true, ucFloatSize, true);
+	
 
 	// We duplicate the time and timestep variables if we're using single-precision so we have copies in both formats
 	if (cModel->getFloatPrecision() == model::floatPrecision::kSingle)
@@ -638,23 +649,23 @@ void CSchemeGodunov::prepare1OMemory(){
 		*(this->oclBufferTimestepMovAvg->getHostBlock<float*>()) = 0.0f;
 		*(this->oclBufferTimeHydrological->getHostBlock<float*>()) = 0.0f;
 		*(this->oclBufferTimeTarget->getHostBlock<float*>()) = 0.0f;
-		*(this->oclBufferBilanDef->getHostBlock<float*>()) = 0.0f;
+		
 	}
 	else {
 		*(this->oclBufferTime->getHostBlock<double*>()) = this->dCurrentTime;
 		*(this->oclBufferTimestep->getHostBlock<double*>()) = this->dCurrentTimestep;
 		*(this->oclBufferTimestepMovAvg->getHostBlock<double*>()) = 0.0;
 		*(this->oclBufferTimeHydrological->getHostBlock<double*>()) = 0.0;
 		*(this->oclBufferTimeTarget->getHostBlock<double*>()) = 0.0;
-		*(this->oclBufferBilanDef->getHostBlock<double*>()) = 0.0;
+		
 	}
 	//create the buffers
 	this->oclBufferTimestep->createBuffer();
 	this->oclBufferTime->createBuffer();
 	this->oclBufferTimestepMovAvg->createBuffer();
 	this->oclBufferTimeHydrological->createBuffer();
 	this->oclBufferTimeTarget->createBuffer();
-	this->oclBufferBilanDef->createBuffer();
+	
 
 	// --
 	// Timestep reduction global array
@@ -791,6 +802,7 @@ void CSchemeGodunov::release1OResources()
 	if (this->oclBufferCellBoundary != NULL)				delete oclBufferCellBoundary;
 	if (this->oclBufferCouplingIDs != NULL)				delete oclBufferCouplingIDs;
 	if (this->oclBufferCouplingValues != NULL)			delete oclBufferCouplingValues;
+	if (this->oclBufferBilanDef != NULL)   delete this->oclBufferBilanDef;
 	if (this->oclBufferUsePoleni != NULL)					delete oclBufferUsePoleni;
 	if (this->oclBuffer_opt_zxmax != NULL)				delete oclBuffer_opt_zxmax;
 	if (this->oclBuffer_opt_cx != NULL)					delete oclBuffer_opt_cx;
@@ -820,6 +832,7 @@ void CSchemeGodunov::release1OResources()
 	oclBufferCellBoundary = NULL;
 	oclBufferCouplingIDs = NULL;
 	oclBufferCouplingValues = NULL;
+	this->oclBufferBilanDef = NULL;
 	oclBufferUsePoleni = NULL;
 	oclBuffer_opt_zxmax = NULL;
 	oclBuffer_opt_cx = NULL;
@@ -853,10 +866,12 @@ void	CSchemeGodunov::prepareSimulation()
 	oclBufferCellManning->queueWriteAll();
 	if (this->bUseOptimizedBoundary == false) {
 		oclBufferCellBoundary->queueWriteAll();
+		this->oclBufferBilanDef->queueWriteAll();
 	}
 	else {
 		oclBufferCouplingIDs->queueWriteAll();
 		oclBufferCouplingValues->queueWriteAll();
+		this->oclBufferBilanDef->queueWriteAll();
 	}
 	oclBufferUsePoleni->queueWriteAll();
 	oclBuffer_opt_zxmax->queueWriteAll();
@@ -867,7 +882,7 @@ void	CSchemeGodunov::prepareSimulation()
 	oclBufferTimestep->queueWriteAll();
 	oclBufferTimestepMovAvg->queueWriteAll();
 	oclBufferTimeHydrological->queueWriteAll();
-	this->oclBufferBilanDef->queueWriteAll();
+	
 
 	this->pDomain->getDevice()->blockUntilFinished();
 
@@ -1036,7 +1051,7 @@ void CSchemeGodunov::Threaded_runBatch()
 			if (this->dCurrentTime < dTargetTime - 1e-8) {
 				oclKernelResetCounters->scheduleExecution();
 				for (unsigned int i = 0; i < uiQueueAmount; i++) {
-
+					///here the clacluation is started
 					this->scheduleIteration();
 					uiIterationsSinceTargetChanged++;
 					ulCurrentCellsCalculated += this->pDomain->getCellCount();
@@ -1052,7 +1067,56 @@ void CSchemeGodunov::Threaded_runBatch()
 			oclBufferBatchSkipped->queueReadAll();
 			oclBufferBatchSuccessful->queueReadAll();
 			oclBufferBatchTimesteps->queueReadAll();
-			this->oclBufferBilanDef->queueReadAll();
+			this->sum_deficit_bilan = 0.0;
+			if (this->bUseOptimizedBoundary == false) {
+				this->oclBufferBilanDef->queueReadAll();
+				for (int i = 0; i < pDomain->getCellCount();i++) {
+					if (cModel->getFloatPrecision() == model::floatPrecision::kSingle)
+					{
+						float* data = (this->oclBufferBilanDef->getHostBlock<float*>());
+						
+						this->sum_deficit_bilan = this->sum_deficit_bilan+data[i];
+						
+					}
+					else {
+						
+						double* data = (this->oclBufferBilanDef->getHostBlock<double*>());
+
+						this->sum_deficit_bilan = this->sum_deficit_bilan+data[i];
+						
+					}
+					
+
+
+				}
+
+
+
+			}
+			else {
+				this->oclBufferBilanDef->queueReadAll();
+				for (int i = 0; i < this->ulCouplingArraySize; i++) {
+					if (cModel->getFloatPrecision() == model::floatPrecision::kSingle)
+					{
+						float* data = (this->oclBufferBilanDef->getHostBlock<float*>());
+
+						this->sum_deficit_bilan = this->sum_deficit_bilan+ data[i];
+
+					}
+					else {
+
+						double* data = (this->oclBufferBilanDef->getHostBlock<double*>());
+
+						this->sum_deficit_bilan = this->sum_deficit_bilan+ data[i];
+
+					}
+
+
+				}
+			}
+
+
+			
 
 			this->pDomain->getDevice()->blockUntilFinished();
 
@@ -1082,7 +1146,7 @@ void CSchemeGodunov::Threaded_runBatch()
 	this->bThreadTerminated = true;
 }
 
-//Runs the actual simulation until completion or error
+//Runs the actual simulation until completion or error here the simulation is set together via the kernels this connects the clc-methods as kernels 
 void	CSchemeGodunov::scheduleIteration() {
 	//std::cout << "start iteration" << std::endl;
 	COCLBuffer* bufferSrc = NULL;
@@ -1106,6 +1170,7 @@ void	CSchemeGodunov::scheduleIteration() {
 	oclKernelTimestepReduction->assignArgument(0, bufferDst);
 
 	profiless
+	//here the timesteps is claculated depending to the scheme, e.g. function pro_cacheDisabled in CLSchemePromaides.clc/.h
 	oclKernelFullTimestep->scheduleExecution();
 	profilese
 
@@ -1115,11 +1180,14 @@ void	CSchemeGodunov::scheduleIteration() {
 	}
 
 	profilebs
+	//here the´bounadries are set e.g. function pro_cacheDisabled in CLSchemePromaides.clc/.h
 	oclKernelBoundary->scheduleExecution();
 	profilebe
 
 	profilers
-	if (this->bDynamicTimestep) { oclKernelTimestepReduction->scheduleExecution(); }
+	if (this->bDynamicTimestep) { 
+		oclKernelTimestepReduction->scheduleExecution(); 
+	}
 	profilere
 
 	profilets
@@ -1220,15 +1288,7 @@ void CSchemeGodunov::setDebugger(unsigned int debugX, unsigned int debugY) {
 //Function to get the bilan deficit
 double CSchemeGodunov::get_bilan_deficit(void) {
 	if (this->oclBufferBilanDef != NULL) {
-		
-
-		if (cModel->getFloatPrecision() == model::floatPrecision::kSingle)
-		{
-			return *(this->oclBufferBilanDef->getHostBlock<float*>());
-		}
-		else {
-			return *(this->oclBufferBilanDef->getHostBlock<double*>());
-		}
+		return this->sum_deficit_bilan;
 	}
 	else {
 		return 0.0;

source_code/system_hydraulic/source_code/solvergpu/source_code/CSchemeGodunov.h

@@ -101,6 +101,7 @@ class CSchemeGodunov : public CScheme
 		unsigned int		uiDebugCellX;											// Debug info cell X
 		unsigned int		uiDebugCellY;											// Debug info cell Y
 		unsigned int		uiTimestepReductionWavefronts;							// Number of wavefronts used in reduction
+		double sum_deficit_bilan; //Sum of Bilan deficitin
 
 		// Private functions
 		virtual void		prepareCode();											// Prepare the code required

source_code/system_hydraulic/source_code/solvergpu/source_code/opencl/CLBoundaries.clc

@@ -41,6 +41,9 @@ __kernel void bdy_Promaides (
 	if (dLclTimestep <= 0.0 || pCellData.y <= -9999.0 || pCellData.x == -9999.0){
 		return;
 	}
+	if(dRate==0.0){
+		return;
+	}
 
 	//if (ulIdx == 7116 ){
 	//	//printf("Postion: X: %f Y:%f \n",lIdxX*1.0 , lIdxY*1.0 );
@@ -56,11 +59,11 @@ __kernel void bdy_Promaides (
 	//}
 
 	//calculate the deficit
-	if(pCellData.x + dRate * dLclTimestep<dCellBedElev){
-		
-		(*pBilanDef)=(*pBilanDef)+((dRate-((dCellBedElev-pCellData.x)/dLclTimestep))*dLclTimestep)*-1.0;
-		//printf("deficit  %f %f %f\n",(*pBilanDef), dLclTimestep, dRate);
-		
+	if(dRate<0.0){
+		if(pCellData.x + dRate * dLclTimestep<dCellBedElev){
+			pBilanDef[ulIdx]=pBilanDef[ulIdx]+((dRate-((dCellBedElev-pCellData.x)/dLclTimestep))*dLclTimestep)*-1.0;
+			//printf("deficit  %f %f %f %f\n",pBilanDef[ulIdx], dLclTimestep, dRate, pCellData.x );
+		}
 	}
 
 	// Apply the value...
@@ -107,12 +110,18 @@ __kernel void bdy_Promaides_by_id (
 	if (dLclTimestep <= 0.0){
 		return;
 	}
+	if(dRate==0.0){
+		return;
+	}
 
 	//calculate the deficit
-	if(pCellData.x + dRate * dLclTimestep<dCellBedElev){
+	if(dRate<0.0){
+		if(pCellData.x + dRate * dLclTimestep<dCellBedElev){
+		
+			pBilanDef[lId]=pBilanDef[ulIdx]+((dRate-((dCellBedElev-pCellData.x)/dLclTimestep))*dLclTimestep)*-1.0;
+			//printf("deficit  %f %f %f %f\n",pBilanDef[lId], dLclTimestep, dRate, pCellData.x );
 
-		(*pBilanDef)=(*pBilanDef)+((dRate-((dCellBedElev-pCellData.x)/dLclTimestep))*dLclTimestep)*-1.0;
-		//printf("deficit  %f %f %f\n",(*pBilanDef), dLclTimestep, dRate);
+		}
 
 	}
 

source_code/system_hydraulic/source_code/solvergpu/source_code/opencl/CLSchemePromaides.clc

@@ -15,6 +15,8 @@
 /*
  *  Calculate everything without using LDS caching
  */
+ //______________________
+ //calculation of the overfow for the four neihgboring elements; it is i decide her if Manning (normal flow) or poleni (weir-flow over 1d-dikeline) is used
 __kernel REQD_WG_SIZE_FULL_TS
 void pro_cacheDisabled ( 
 			__constant	cl_double *  				dTimestep,					// Timestep
@@ -315,7 +317,8 @@ void pro_cacheDisabled (
 	pCellStateDst[ ulIdx ] = pCellData;
 }
 
-
+//_________________________________________________
+//calculation the overflow via Manning-equation (normal flow)
 cl_double2 manning_Solver(
 	cl_double timestep,		// Current Time-step
 	cl_double manning,     	// Manning of Main Cell
@@ -416,6 +419,8 @@ cl_double2 manning_Solver(
 
 	}
 
+//_________________________________________________
+//calculation the overflow between elements via Poleni weir equation (1d-dike line element)
 cl_double2 poleni_Solver(
 	cl_double timestep,		// Current Time-step
 	cl_double manning,     	// Manning of Main Cell