Revision of the euler and the event handling for the euler until now …

…only work with "if" equations. git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@4567 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
OpenModelica · Nov 24, 2009 · 3d6fea6 · 3d6fea6
1 parent dc2bb27
commit 3d6fea6
Show file tree

Hide file tree

Showing 4 changed files with 337 additions and 74 deletions.
diff --git a/Compiler/SimCodegen.mo b/Compiler/SimCodegen.mo
@@ -6846,6 +6846,7 @@ algorithm
     local
       Codegen.CFunction func_zc,func_handle_zc,cfunc,cfunc0_1,cfunc0,cfunc_1,cfunc_2,func_zc0,func_handle_zc0,func_handle_zc0_1;
       Codegen.CFunction func_handle_zc0_2,func_handle_zc0_3,func_zc0_1,func_zc_1,func_handle_zc_1,cfuncHelpvars;
+      Codegen.CFunction func_ozc;
       Integer cg_id1,cg_id2,cg_id;
       list<CFunction> extra_funcs1,extra_funcs2,extra_funcs;
       String extra_funcs_str,helpvarUpdateStr,func_str,res,cname;
@@ -6890,7 +6891,14 @@ algorithm
         func_zc0_1 = Codegen.cAddCleanups(func_zc0, {"localData->timeValue = timeBackup;", "return 0;"});
         func_zc_1 = Codegen.cMergeFns({func_zc0_1,func_zc});
         func_handle_zc_1 = Codegen.cMergeFns({func_handle_zc0_3,func_handle_zc});
-        func_str = Codegen.cPrintFunctionsStr({func_zc_1,func_handle_zc_1,cfunc_2});
+
+
+        func_ozc = Codegen.cMakeFunction("int", "function_onlyZeroCrossings", {}, {"double *gout,double *t"});
+        func_ozc = addMemoryManagement(func_ozc);
+        func_ozc = Codegen.cMergeFns({func_ozc,func_zc});
+        func_ozc = Codegen.cAddCleanups(func_ozc, {"return 0;"});
+
+        func_str = Codegen.cPrintFunctionsStr({func_zc_1,func_handle_zc_1,cfunc_2,func_ozc});
         res = Util.stringAppendList({extra_funcs_str,func_str});
       then
         res;

diff --git a/c_runtime/simulation_events.cpp b/c_runtime/simulation_events.cpp
@@ -423,6 +423,190 @@ bool edge(double& var) {
 }
 
 bool change(double& var) {
-  return var && !pre(var) || !var && pre(var);
+  return (var && !pre(var)) || (!var && pre(var));
 }
 
+/*
+ * All event functions from new, are till now only used in Euler  
+ * 
+*/
+
+//
+// This function checks for Events in Intervall=[oldTime,timeValue]
+// If a zerocrossing Function cause a sign chage, root finding  
+// process will start
+//
+int CheckForNewEvent(int flag) {
+
+	if (flag != INTERVAL){
+	   while(checkForDiscreteVarChanges()) {
+		   saveall();
+		   function_updateDependents();
+	   }
+	}
+
+	function_onlyZeroCrossings(gout,&globalData->timeValue);
+
+	for (long i = 0; i < globalData->nZeroCrossing; i++) {
+
+		if (gout[i] < 0) { // check also zero crossings that are on zero.
+			if (sim_verbose) cout << "gout[" << i << "] = " << gout[i] << endl;
+			if (sim_verbose) {
+				cout << "adding event " << i << " at time: "
+				<< globalData->timeValue << endl;
+			}
+			AddEvent(i);
+		}
+	}
+	if (!EventQueue.empty() && flag == INTERVAL){
+		FindRoot();
+		EventHandle();
+		return 1;
+	}else if(!EventQueue.empty()){
+		EventHandle();
+		return 1;
+	}	
+	return 0;
+}
+
+//
+// This function handle events and change all 
+// needed variables for an event
+//
+void EventHandle(){
+
+	while(!EventQueue.empty()){
+		long event_id;
+
+		event_id = EventQueue.front(); 
+
+		if (sim_verbose) cout << "Handle Event ID: " << event_id << endl;
+		if (zeroCrossingEnabled[event_id] == 1){
+			zeroCrossingEnabled[event_id] = -1;}
+		else if (zeroCrossingEnabled[event_id] == -1){
+			zeroCrossingEnabled[event_id] = 1;}
+		else{
+			zeroCrossingEnabled[event_id] = 1;}
+
+		handleZeroCrossing(event_id);
+	    function_updateDependents();
+	    saveall();
+
+		EventQueue.pop_front();
+	}
+	CheckForNewEvent(NOINTERVAL);
+}
+
+//
+// This function perform a root finding for 
+// Intervall=[oldTime,timeValue]
+//
+void FindRoot(){
+	double EventTime = 0;
+	long int event_id =0;
+
+	double *states_right = new double[globalData->nStates];
+	double *states_left = new double[globalData->nStates];
+
+	//write states to work and backup Array
+	for(int i=0;i<globalData->nStates;i++){
+		states_right[i] = globalData->states[i];
+		states_left[i] = globalData->oldStates[i];
+	}
+
+	// Search for event time with Bisection method
+	EventTime = BiSection(globalData->oldTime, globalData->timeValue, states_right, states_left, &event_id);
+
+	//if (EventTime!=0){ // Found event at EventTime
+
+	globalData->timeValue = EventTime;
+	if (sim_verbose) {
+		cout << "Found event " << event_id << " at time: "<< globalData->timeValue << endl;
+	}
+	AddEvent(event_id);
+
+	for(int i=0;i<globalData->nStates;i++){
+		globalData->states[i] = states_left[i];
+	}
+	functionODE();
+	functionDAE_output();
+	emit();
+
+
+	for(int i=0;i<globalData->nStates;i++){
+		globalData->states[i] = states_right[i];
+	}
+
+	delete[] states_left;
+	delete[] states_right;
+
+	EventHandle();
+}
+
+// 
+// Method to find root in Intervall[oldTime,timeValue]
+//
+double BiSection(double a, double b, double* states_a, double* states_b,long int* event_id )	
+{
+
+	//double TTOL =  DBL_EPSILON;//*fabs(2*b-a)*100;
+	double TTOL = 1e-05;
+	if(TOL!=0) TTOL = TOL;
+	double c;
+
+	if (sim_verbose){
+			cout << "Check Intervall [" << a << "," << b << "]" << endl; 
+			cout << "TTOL is set to: " << TTOL << endl;
+	}
+
+	while ( fabs(b-a) > TTOL){
+
+		c = (a+b)/2.0;
+		globalData->timeValue = c;
+
+		//calculates states at time c 
+		for(int i=0;i<globalData->nStates;i++){
+			globalData->states[i] = (states_a[i] + states_b[i]) / 2.0; 
+		}
+
+		//calculates Values dependents on new states 
+		functionODE();
+		functionDAE_output();
+
+		if ( CheckZeroCrossings(event_id)){ //If Zerocrossing in left Section
+
+			//if (sim_verbose) cout << "\t\tSearch in the left section" << endl;
+			for(int i=0;i<globalData->nStates;i++){
+				states_a[i] = globalData->states[i];
+			}
+			b = c;
+		}else{   //If Zerocrossing in right Section
+
+			//if (sim_verbose) cout << "\t\tSearch in the right section" << endl;
+			for(int i=0;i<globalData->nStates;i++){
+				states_b[i] = globalData->states[i];
+			}
+			a = c;
+		} 
+	}
+
+	c = (a+b)/2.0;
+	return c;
+}
+
+
+//
+// Check if at least one zerocrossing has change sign
+// is used in BiSection
+//
+int CheckZeroCrossings(long int *eventid) {
+  function_onlyZeroCrossings(gout,&globalData->timeValue);
+  for (long i = 0; i < globalData->nZeroCrossing; i++) {
+	  //if (sim_verbose) cout << "check gout[" << i << "] = " << gout[i] << endl;
+	  if (gout[i] < 0) {
+		  *eventid = i;
+		  return 1;
+	  }
+  }
+  return 0;
+}
diff --git a/c_runtime/simulation_events.h b/c_runtime/simulation_events.h
@@ -118,4 +118,22 @@ function_when(int i);
 
 extern long* zeroCrossingEnabled;
 
+int
+function_onlyZeroCrossings(double* gout ,double* t);
+
+int CheckForNewEvent(int flag);
+
+void EventHandle();
+
+void FindRoot();
+
+double BiSection(double, double, double*, double*, long int*);
+
+int CheckZeroCrossings(long int*);
+
+#define INTERVAL 1
+#define NOINTERVAL 0
+
+extern double TOL;
+
 #endif