HYD system:

- Bug fix FP angle -90 degree
- added a dynamic stepsize for FP2FP coupling via Courant criteria
- added hit-key figure, how often the limiter in FP2FP coupling is hit; also a percentage evaluation
- added the error negative discharge boundaries, which can not be extracted out of the model to the hydrological bilan

Author: dabachma

source_code/system_hydraulic/source_code/models/coupling/Hyd_Coupling_FP2FP.cpp

@@ -8,6 +8,9 @@ Hyd_Coupling_FP2FP::Hyd_Coupling_FP2FP(void){
 	this->defining_polysegment.set_line_type(_hyd_floodplain_polysegment_type::FP_BOUNDARYLINE);
 	this->area_limiter = 0.0;
 	this->counter_limiter = 0;
+	this->counter_tot_coupling_active=0;
+	this->timestep_opt_hit_limit=0.0;
+	this->timestep_opt_cor=0.0;
 
 	//count the memory
 	Sys_Memory_Count::self()->add_mem(sizeof(Hyd_Coupling_FP2FP)-sizeof(Hyd_Coupling_Point_FP2FP_List), _sys_system_modules::HYD_SYS);
@@ -53,6 +56,9 @@ void Hyd_Coupling_FP2FP::init_coupling(void){
 		this->area_limiter=*this->floodplain_model_1->Param_FP.get_ptr_width_x()* *this->floodplain_model_1->Param_FP.get_ptr_width_y();
 		this->area_limiter = min(this->area_limiter, *this->floodplain_model_2->Param_FP.get_ptr_width_x() * *this->floodplain_model_2->Param_FP.get_ptr_width_y());
 		this->counter_limiter = 0;
+		this->counter_tot_coupling_active = 0;
+		this->timestep_opt_hit_limit = 0.0;
+		this->timestep_opt_cor = 0.0;
 
 		//sort it along the defining line
 		this->list.sort_distance_along_polysegment();
@@ -96,8 +102,12 @@ void Hyd_Coupling_FP2FP::synchronise_models(const double timepoint, const double
 
 
 		//}
+		if (time_check == true) {
+			this->timestep_opt_hit_limit = 0.0;
+			this->timestep_opt_cor = 0.0;
+		}
 
-		this->list.syncronisation_models_bylistpoints(timepoint, delta_t, time_check, internal_counter, this->area_limiter,&this->counter_limiter);
+		this->list.syncronisation_models_bylistpoints(timepoint, delta_t, time_check, internal_counter, this->area_limiter,&this->counter_limiter, &this->counter_tot_coupling_active, &this->timestep_opt_hit_limit, &this->timestep_opt_cor);
 
 	}
 	catch(Error msg){
@@ -145,6 +155,9 @@ void Hyd_Coupling_FP2FP::clone_couplings(Hyd_Coupling_FP2FP *coupling, Hyd_Hydra
 	this->defining_polysegment.clone_polysegment(&coupling->defining_polysegment);
 	this->area_limiter = coupling->area_limiter;
 	this->counter_limiter = 0;
+	this->counter_tot_coupling_active = 0;
+	this->timestep_opt_hit_limit = 0.0;
+	this->timestep_opt_cor = 0.0;
 
 	this->list.set_defining_polysegment(&this->defining_polysegment);
 	this->list.clone_list(&coupling->list, this->floodplain_model_1, this->floodplain_model_2); 
@@ -153,13 +166,25 @@ void Hyd_Coupling_FP2FP::clone_couplings(Hyd_Coupling_FP2FP *coupling, Hyd_Hydra
 //Reset counter limiter
 void Hyd_Coupling_FP2FP::reset_counter_limiter(void) {
 	this->counter_limiter = 0;
+	this->counter_tot_coupling_active = 0;
 }
 //Output number of limiter hits
 void Hyd_Coupling_FP2FP::output_number_limiter_hits(int* total, ostringstream* out) {
-	*out << " Hits FP" <<this->floodplain_model_1->Param_FP.get_floodplain_number()<<" to FP"<< this->floodplain_model_2->Param_FP.get_floodplain_number()<<"   :" << this->counter_limiter << endl;
+	double perc = 0.0;
+	if (this->counter_tot_coupling_active > 0) {
+		perc = ((double)this->counter_limiter / (double)this->counter_tot_coupling_active) * 100.0;
+
+	}
+	*out << " Hits FP" <<this->floodplain_model_1->Param_FP.get_floodplain_number()<<" to FP"<< this->floodplain_model_2->Param_FP.get_floodplain_number()<<"   :" << this->counter_limiter <<" (active: "<< this->counter_tot_coupling_active<< P(1) << " "<< perc<<"%)"<< endl;
 
 	*total = *total + this->counter_limiter;
 
+}
+//Get the minimal timestep for coupling
+double Hyd_Coupling_FP2FP::get_min_timestep(void) {
+
+	return this->list.get_min_timestep();
+
 }
 //______________
 //private

source_code/system_hydraulic/source_code/models/coupling/Hyd_Coupling_FP2FP.h

@@ -61,6 +61,9 @@ class Hyd_Coupling_FP2FP
 	///Output number of limiter hits
 	void output_number_limiter_hits(int* total, ostringstream* out);
 
+	///Get the minimal timestep for coupling
+	double get_min_timestep(void);
+
 
 private:
 	//members
@@ -75,6 +78,13 @@ class Hyd_Coupling_FP2FP
 	///Counter how often the limiter is reached
 	int counter_limiter;
 
+	///Total counter how often a coupling is active
+	int counter_tot_coupling_active;
+	///Optimal timestep to avoid hitting limiter
+	double timestep_opt_hit_limit;
+	///Optimal timestep after Courant-criteria
+	double timestep_opt_cor;
+
 	//methods
 	///Set the defining polysegment with the boundary of the first floodplain model
 	void set_defining_polysegment(void);

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

@@ -1411,6 +1411,24 @@ void Hyd_Coupling_Management::output_number_limiter_hits(void) {
 	}
 
 
+}
+///Get minimal timestep for FP2FP couplinbg
+double Hyd_Coupling_Management::get_min_timestep_fp2fp(void) {
+	double min_t = 1000000000;
+
+	if (this->number_fp2fp > 0) {
+		
+
+
+		for (int i = 0; i < this->number_fp2fp; i++) {
+			min_t = min(this->coupling_fp2fp[i].list.get_min_timestep(), min_t);
+		}
+
+
+	}
+
+	return min_t;
+
 }
 //__________________
 //private

source_code/system_hydraulic/source_code/models/coupling/Hyd_Coupling_Management.h

@@ -246,6 +246,9 @@ class Hyd_Coupling_Management : public _Sys_Common_System
 	///Out number of limiter hits
 	void output_number_limiter_hits(void);
 
+	///Get minimal timestep for FP2FP couplinbg
+	double get_min_timestep_fp2fp(void);
+
 
 private:
 

source_code/system_hydraulic/source_code/models/coupling/Hyd_Coupling_Point_FP2FP.cpp

@@ -194,7 +194,7 @@ void Hyd_Coupling_Point_FP2FP::reset_coupling_discharge(void){
 	}
 }
 //Syncronisation of the coupled models with the couplingspoint
-void Hyd_Coupling_Point_FP2FP::syncronisation_coupled_models(const double timepoint, const double delta_t, const bool time_check, const int internal_counter, const double area_limiter, int *counter_limiter){
+void Hyd_Coupling_Point_FP2FP::syncronisation_coupled_models(const double timepoint, const double delta_t, const bool time_check, const int internal_counter, const double area_limiter, int *counter_limiter, int* counter_tot, double* opt_time_hit, double* opt_time_cor){
 
 	_Hyd_Coupling_Point::syncronisation_coupled_models();
 	this->delta_t=delta_t;
@@ -326,24 +326,18 @@ void Hyd_Coupling_Point_FP2FP::syncronisation_coupled_models(const double timepo
 		//	this->q_list.smooth_discharge.append(this->current_q);
 		//}
 
-		//calculate the flow-velocity
-		if(flow_depth!=0.0){
-			this->coupling_v=-1.0*this->current_q/(this->distance_down*flow_depth);
-		}
-		else{
-			this->coupling_v=0.0;
-		}
 
 
 
 
-		//if (timepoint > 48600 && abs(this->current_q) > 0.0) {
+
+		//if (timepoint > 83000 && abs(this->current_q) > 0.0 && time_check==true) {
 		//	//	//	Daniel test alfnitz
 		//	ostringstream cout;
 		//	cout << "t " << timepoint << " fp_1_elem " << this->fp1_elem_index << " fp_2_elem " << this->fp2_elem_index ;
 		//	cout << " s_1 " << this->fp1_elem->get_z_value() + h_one_buff << " s_2 " << this->fp2_elem->get_z_value() + h_two_buff;
 		//	cout << " wd " << this->distance_down << " flow_h " << flow_depth;
-		//	cout << " cp_Q_buf " << q_buff << " cp_Q_s " << this->current_q<< " Q_max "<< max_Q;
+		//	cout << " cp_Q_buf " << q_buff << " cp_Q_s " << this->current_q;
 		//	if (abs(this->current_q) > 5) {
 		//		cout << "*";
 		//	}
@@ -362,20 +356,64 @@ void Hyd_Coupling_Point_FP2FP::syncronisation_coupled_models(const double timepo
 
 
 		//Test a delimiter
+		//Total overflow to total hits delimiter
+		if (abs(this->current_q) > 0.0) {
+			if (time_check == false) {
+				(*counter_tot)++; //double* opt_time_cor
+			}
+		}
+
 		double max_Q = 0.0;
+		//calculate the delimiter
 		max_Q = (area_limiter* this->distance_down) / (4 * delta_t);
 		if (max_Q < abs(this->current_q)) {
 			if (time_check == false) {
 				(*counter_limiter)++;
 			}
 		}
+		//calculate optimal timestep for not hitting delimiter
+		if (time_check == true) {
+			if (max_Q < abs(this->current_q)) {
+
+				(*opt_time_hit) = (area_limiter * this->distance_down) / (4 * abs(this->current_q));
+			}
+			else {
+
+				(*opt_time_hit) = -1.0;
+			}
+		}
+		//set the delimiter
 		max_Q=min(abs(this->current_q), max_Q);
+	
 		if (this->current_q < 0.0) {
 			this->current_q = -1.0 * max_Q;
 		}
 		else {
 			this->current_q =  max_Q;
 		}
+		
+
+		//calculate the flow-velocity
+		if (flow_depth != 0.0) {
+			this->coupling_v = -1.0 * this->current_q / (this->distance_down * flow_depth);
+		}
+		else {
+			this->coupling_v = 0.0;
+		}
+
+		//calculate optimal timestep with Courant-Criteria
+		if (time_check == true) {
+			double courant = 0.5;
+			if (abs(this->coupling_v) > 0) {
+				(*opt_time_cor) = (courant * this->flow_distance) / abs(this->coupling_v);
+				
+			}
+			else {
+				(*opt_time_cor) = -1.0;
+
+			}
+
+		}
 
 
 		this->calculate_maximum_values(timepoint, this->current_q, &this->max_coupling_v);

source_code/system_hydraulic/source_code/models/coupling/Hyd_Coupling_Point_FP2FP.h

@@ -60,7 +60,7 @@ friend class Hyd_Coupling_Point_FP2FP_List;
 	///Reset the current coupling discharge of the points and the coupled element
 	void reset_coupling_discharge(void);
 	///Syncronisation of the coupled models with the couplingspoint
-	void syncronisation_coupled_models(const double timepoint, const double delta_t, const bool time_check, const int internal_counter, const double area_limiter, int *counter_limiter);
+	void syncronisation_coupled_models(const double timepoint, const double delta_t, const bool time_check, const int internal_counter, const double area_limiter, int *counter_limiter, int *counter_tot, double *opt_time_hit, double *opt_time_cor);
 
 	///Reset the smoothing calculation members
 	void reset_smoothing(void);

source_code/system_hydraulic/source_code/models/coupling/Hyd_Coupling_Point_FP2FP_List.cpp

@@ -5,6 +5,7 @@
 Hyd_Coupling_Point_FP2FP_List::Hyd_Coupling_Point_FP2FP_List(void){
 	this->number=0;
 	this->points=NULL;
+	this->min_timestep = 0.0;
 
 	//count the memory
 	Sys_Memory_Count::self()->add_mem(sizeof(Hyd_Coupling_Point_FP2FP_List), _sys_system_modules::HYD_SYS);
@@ -128,6 +129,7 @@ void Hyd_Coupling_Point_FP2FP_List::set_new_couplingpoint(Hyd_Coupling_Point_FP2
 void Hyd_Coupling_Point_FP2FP_List::delete_list(void){
 	this->delete_couplingpoints();
 	this->number=0;
+	this->min_timestep = 0.0;
 }
 //Sort the points after the distance to the beginning of the defining polysegment
 void Hyd_Coupling_Point_FP2FP_List::sort_distance_along_polysegment(void){
@@ -260,12 +262,54 @@ void Hyd_Coupling_Point_FP2FP_List::reset_points(void){
 	}
 }
 //Syncronisation of the coupled models with the stored couplingspoints in the list
-void Hyd_Coupling_Point_FP2FP_List::syncronisation_models_bylistpoints(const double timepoint, const double delta_t, const bool time_check, const int internal_counter, const double area_limiter, int *counter_limiter){
+void Hyd_Coupling_Point_FP2FP_List::syncronisation_models_bylistpoints(const double timepoint, const double delta_t, const bool time_check, const int internal_counter, const double area_limiter, int *counter_limiter, int* counter_tot, double* opt_time_hit, double* opt_time_cor){
+	Sys_Median_Calculator by_hits;
+	Sys_Median_Calculator by_cor;
+
 	for(int i=0; i<this->number;i++){
 
-		this->points[i].syncronisation_coupled_models(timepoint,delta_t, time_check, internal_counter, area_limiter, counter_limiter);	
+		this->points[i].syncronisation_coupled_models(timepoint,delta_t, time_check, internal_counter, area_limiter, counter_limiter, counter_tot,opt_time_hit, opt_time_cor);
+
+		//claculate the median of the timesteps
+		if (time_check == true) {
+			
+			if (*opt_time_hit > 0.0) {
+				by_hits.add_value(*opt_time_hit);
+			}
+			
+			if (*opt_time_cor > 0.0) {
+				by_cor.add_value(*opt_time_cor);
+				//output
+				/*ostringstream cout;
+				cout << "t " << timepoint << " step " << delta_t << " opt_time " << (*opt_time_cor);
+				cout << endl;
+				Sys_Common_Output::output_hyd->output_txt(&cout, true);*/
+			}
+
+
+
+
+
+		}
+
 
 	}
+
+	if(time_check==true){
+		//output
+		//ostringstream cout;
+		//cout << "t " << timepoint << " step " << delta_t<< " med_hit " << by_hits.get_median()<< " perc25_hit "<< by_hits.get_perc25();
+		//cout << " med_cor " << by_cor.get_median() << " perc_25 " << by_cor.get_perc25() << " s";
+		//cout << endl;
+		//Sys_Common_Output::output_hyd->output_txt(&cout, true);
+		 
+		
+		//Make decision wich one courant crit or hit crit ...median, perc etc.
+		this->min_timestep = by_cor.get_median();
+	
+	}
+
+
 }
 //Clone the coupling point list
 void Hyd_Coupling_Point_FP2FP_List::clone_list(Hyd_Coupling_Point_FP2FP_List *list, Hyd_Model_Floodplain *fp1, Hyd_Model_Floodplain *fp2){
@@ -284,6 +328,10 @@ void Hyd_Coupling_Point_FP2FP_List::clone_list(Hyd_Coupling_Point_FP2FP_List *li
 		this->points[i].set_second_fpelem_pointer(fp2->get_ptr_floodplain_elem(list->points[i].fp2_elem_index));
 	}
 }
+//Get the minimal timestep
+double Hyd_Coupling_Point_FP2FP_List::get_min_timestep(void) {
+	return this->min_timestep;
+}
 //Output the discharge list to file
 //void Hyd_Coupling_Point_FP2FP_List::output_discharge_lists(void){
 //	for(int i=0; i< this->number; i++){

source_code/system_hydraulic/source_code/models/coupling/Hyd_Coupling_Point_FP2FP_List.h

@@ -72,20 +72,25 @@ class Hyd_Coupling_Point_FP2FP_List : public _Hyd_Coupling_Point_List
 	void reset_points(void);
 
 	///Syncronisation of the coupled models with the stored couplingspoints in the list
-	void syncronisation_models_bylistpoints(const double timepoint, const double delta_t, const bool time_check, const int internal_counter, const double area_limiter, int *counter_limiter);
+	void syncronisation_models_bylistpoints(const double timepoint, const double delta_t, const bool time_check, const int internal_counter, const double area_limiter, int *counter_limiter, int* counter_tot, double* opt_time_hit, double* opt_time_cor);
 
 	///Clone the coupling point list
 	void clone_list(Hyd_Coupling_Point_FP2FP_List *list, Hyd_Model_Floodplain *fp1, Hyd_Model_Floodplain *fp2);
 
 	//Output the discharge list to file (just for testing)
 	//void output_discharge_lists(void);
 
+	///Get the minimal timestep
+	double get_min_timestep(void);
+
 
 private:
 
 	//members
 	///Coupling points in list
 	Hyd_Coupling_Point_FP2FP *points;
+	///Minimal timestep
+	double min_timestep;
 
 
 	//methods

source_code/system_hydraulic/source_code/models/floodplain/Hyd_Model_Floodplain.cpp

@@ -1186,6 +1186,16 @@ void Hyd_Model_Floodplain::make_hyd_balance_max(const double time_point){
 	for(int i=0; i< this->number_coup_cond; i++){
 		this->floodplain_elems[this->coup_cond_id[i]].element_type->calculate_hydrolocigal_balance_coupling(delta_t);
 	}
+
+	if (this->Param_FP.get_scheme_info().scheme_type != model::schemeTypes::kDiffusiveCPU) {
+		//check balanc deficit in gpu modells
+		
+
+		this->error_zero_outflow_volume= this->pManager->getDomain()->getScheme()->get_bilan_deficit()*-1.0;
+
+	}
+
+
 }
 //Reset the solver of the model
 void Hyd_Model_Floodplain::reset_solver(void){

source_code/system_hydraulic/source_code/models/floodplain/Hyd_Parse_FP.cpp

@@ -176,6 +176,17 @@ void Hyd_Parse_FP::parse_general(_hyd_keyword_file Key, word Command){
 		}
 		else if (Key == eANGLE){
 			buffer>>this->fp_params.angle;
+
+			if (this->fp_params.angle >= 90.0 - constant::angle_epsilon && this->fp_params.angle <= 90.0 + constant::angle_epsilon) {
+
+				this->fp_params.angle = this->fp_params.angle - constant::angle_epsilon*10;
+			}
+			if (this->fp_params.angle >= -90.0 - constant::angle_epsilon && this->fp_params.angle <= -90.0 + constant::angle_epsilon) {
+
+				this->fp_params.angle = this->fp_params.angle + constant::angle_epsilon*10;
+			}
+
+
 			this->found_angle=true;
 		}
 		else if (Key == eLOWLEFTX){