Fixed bug in initialization. Now handles initial equation better.

git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@2160 f25d12d1-65f4-0310-ae8a-bbce733d8d8e

Compiler/DAELow.rml

@@ -263,6 +263,8 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
   relation abs_incidence_matrix: IncidenceMatrix => IncidenceMatrix
 
+  relation vars_incidence_matrix: IncidenceMatrix => int list 
+
   relation dump_incidence_matrix: IncidenceMatrix => ()
 
   relation dump_incidence_matrix_t: IncidenceMatrixT => ()
@@ -310,6 +312,7 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
   relation make_all_states_algebraic: DAELow => DAELow 
 
   relation get_var: (Exp.ComponentRef, Variables) => (Var list, int list)
+  relation vars_size: Variables => int 
   relation get_var_using_name: (string, Variables) => (Var,int) 
   relation add_var: (Var, Variables) => Variables
   relation get_var_at:(Variables,int) => (Var)
@@ -3919,6 +3922,20 @@ relation abs_incidence_matrix: IncidenceMatrix => IncidenceMatrix =
 	abs_incidence_matrix(m) => res
 end
 
+(** relation: vars_incidence_matrix
+ ** author: PA
+ ** Return all variable indices in the incidence matrix, i.e.
+ ** all elements of the matrix.
+ **)
+
+relation vars_incidence_matrix: IncidenceMatrix => int list  =
+
+  rule	array_list(m) => mlst &
+	Util.list_flatten(mlst) => res
+	--------------------------
+	vars_incidence_matrix(m) => res
+end
+
 
 (* relation: transpose_row
  ** author: PA
@@ -4125,6 +4142,9 @@ relation matching_algorithm: (DAELow,IncidenceMatrix, IncidenceMatrixT,
 	int_string(neqns) => ne &
 	int_gt(nvars,0) => true &
 	int_gt(neqns,0) => true &
+
+	(* Worst case, all eqns are differentiated once.
+	 * Create nvars*2 assignment elements *)       
 	int_add(nvars,nvars) => memsize &
 	assignments_create(nvars,memsize,0) => assign1 &
 	assignments_create(nvars,memsize,0) => assign2 &
@@ -4134,7 +4154,9 @@ relation matching_algorithm: (DAELow,IncidenceMatrix, IncidenceMatrixT,
 	assignments_vector(ass1) => vec1 &
 	assignments_vector(ass2) => vec2 
 	--------------------------
-	matching_algorithm(dae,m,mt,match_opts) => (vec1,vec2,dae,m,mt)
+	matching_algorithm(dae as DAELOW(vars,_,eqns,_,_,_,_,_)
+			   ,m,mt,match_opts) 
+	  => (vec1,vec2,dae,m,mt)
 
 (*  rule	(*Print.print_error_buf "#Error, matching failed\n" &*)
 	list_vector([]) => v1 & 
@@ -5225,25 +5247,27 @@ end
 (** relation: equation_size
  ** author: PA
  **
- ** Returns the number of equations in an EquationArray.
- ** Note: This is not corresponding to the number of equations in a system,
- ** since both algorithms and array equations can count as more than one.
- ** For counting the number of equations in the system, use equation_size.
+ ** Returns the number of equations in an EquationArray, which
+ ** corresponds to the number of equations in a system.
+ ** NOTE: Array equations and algorithms are represented several times
+ ** in the array so the number of elements of the array corresponds to 
+ ** the equation system size.
  **)
 relation equation_size:(EquationArray) => int =
 
   axiom	equation_size(EQUATION_ARRAY(n,_,_)) => n
 end
 
 
-(** relation: variable_size
+
+(** relation: vars_size
  ** author: PA
  **
  ** Returns the number of variables
  **)
-relation variable_size:(Variables) => int =
+relation vars_size:(Variables) => int =
 
-  axiom	variable_size(VARIABLES(_,_,_,_,n)) => n 
+  axiom	vars_size(VARIABLES(_,_,_,_,n)) => n 
 end
 
 

Compiler/Main.rml

@@ -503,7 +503,7 @@ relation simcodegen: (Absyn.Path, (* classname *)
 	Util.string_append_list([cname_str,".makefile"]) => makefilename &
 	SimCodegen.generate_functions(p,dae,indexed_dlow',classname,funcfilename) => libs &
 	SimCodegen.generate_simulation_code(dae,indexed_dlow',ass1,ass2,m,mt,comps,classname,filename,funcfilename) &
-	SimCodegen.generate_init_data(indexed_dlow',classname,init_filename,0.0,1.0,0.01) &
+	SimCodegen.generate_init_data(indexed_dlow',classname,init_filename,0.0,1.0,500.0) &
 	SimCodegen.generate_makefile(makefilename,cname_str,libs)
 	----------------------------------------------------------------------------------
 	simcodegen(classname,p,dae,dlow,ass1,ass2,m,mt,comps)

Compiler/SimCodegen.rml

@@ -188,12 +188,13 @@ relation generate_simulation_code: (DAE.DAElist,
 	  => (c_eventchecking, helpVarInfo) & 
 	  list_length helpVarInfo => n_h &
         generate_global_data(class,dlow,n_o,n_i,n_h) => cglobal &  
-        generate_compute_output(cname,dae,dlow,ass1,ass2,blt_no_states) 
+     	  generate_compute_output(cname,dae,dlow,ass1,ass2,blt_no_states) 
 	  => coutput &
-          generate_compute_residual_state(cname,dae,dlow,ass1,ass2,blt_states) 
+        	 generate_compute_residual_state(cname,dae,dlow,ass1,ass2,blt_states) 
 	  => cstate &
-	  generate_ode_code(dlow,blt_states,ass1,ass2,m,mt,class) => c_ode &
+		 generate_ode_code(dlow,blt_states,ass1,ass2,m,mt,class) => c_ode &
 	generate_initial_value_code(dlow) => s_code &
+
 	generate_when_clauses(cname,dae,dlow,ass1,ass2,comps) => cwhen &
 	generate_zero_crossing(cname,dae,dlow,ass1,ass2,comps,helpVarInfo) => czerocross & 
 	Util.string_append_list(["//Simulation code for ",cname,
@@ -1175,34 +1176,53 @@ relation generate_initial_value_code:(DAELow.DAELow) => string =
 	  => (param_assigns,cg_id') &
 	list_append(initial_eqns1,initial_eqns2) => initial_eqns &
 	DAELow.list_equation(initial_eqns) => initial_eqns' &
+
 	list_array([]) => arr &
 	list_array([]) => arr2 &
 	let empty_eqn = DAELow.EQUATION_ARRAY(0,0,arr) &
 
 	let initial_dae = DAELow.DAELOW(vars,knvars,initial_eqns',empty_eqn,empty_eqn,arr2,al,ev) &
+
+
 	(* States are no longer special, they are just a variable that 
 	 * needs an initial value. *)
-	DAELow.make_all_states_algebraic(initial_dae) => initial_dae' &
+	DAELow.make_all_states_algebraic(initial_dae) 
+	  => (initial_dae' as DAELow.DAELOW(v,kv,e,re,ie,ae,al,ev)) &
+
  	DAELow.incidence_matrix(initial_dae') => m &
+
+	(* Collect only variables part of the initial equations*)
+	DAELow.vars_incidence_matrix(m) => v_indx_lst &
+	Util.list_map_1r(v_indx_lst,DAELow.get_var_at,v) => varlst &
+	DAELow.list_var(varlst) => vars' &
+	DAELow.var_list(v) => allvars &
+	DAELow.var_list(kv) => knvars1 &
+	(* Make all other variables known*)
+	Util.list_setdifference_p(allvars,varlst,DAELow.var_equal) => knvars2 &
+	list_append(knvars1,knvars2) => knvars &
+	DAELow.list_var(knvars) => knvars' &
+
+	let new_initial_dae = DAELow.DAELOW(vars',knvars',e,re,ie,ae,al,ev) &
+ 	DAELow.incidence_matrix(new_initial_dae) => m &
  	DAELow.transpose_matrix(m) => mT &
 	(* once dae has been translated, state or not should not matter *)
 	DAELow.abs_incidence_matrix(m) => m' &
 	DAELow.abs_incidence_matrix(mT) => mT' &
-	(*DAELow.dump initial_dae' &
-	DAELow.dump_incidence_matrix m' &*)
+
 	(* index reduction for intitial equations makes no sense.. *)
-	DAELow.matching_algorithm(initial_dae',m',mT',(DAELow.NO_INDEX_REDUCTION,DAELow.ALLOW_UNDERCONSTRAINED)) 
+	DAELow.matching_algorithm(new_initial_dae,m',mT',(DAELow.NO_INDEX_REDUCTION,DAELow.ALLOW_UNDERCONSTRAINED)) 
 	  => (ass1,ass2,initial_dae',m,mT) &
 	DAELow.strong_components(m',mT',ass1,ass2) => blocks &
-	generate_ode_blocks(initial_dae',ass1,ass2,blocks,cg_id')
+	generate_ode_blocks(new_initial_dae,ass1,ass2,blocks,cg_id')
 	  => (block_code,_,extra_funcs) &
 	Codegen.c_make_function("int","initial_function",[],["double *x",
-							       "double *xd", 
-							       "double *y",
-							       "double *p",
-							       "int nx",
-							       "int ny",
-							       "int np"])
+							     "double *xd", 
+							     "double *y",
+							     "double *p",
+							     "double *t",
+							     "int nx",
+							     "int ny",
+							     "int np"])
 	  => init_func' &
 
 	Codegen.c_add_cleanups(init_func', ["return 0;"]) => init_func &
@@ -1227,6 +1247,7 @@ relation generate_initial_value_code:(DAELow.DAELow) => string =
 							       "double *xd", 
 							       "double *y",
 							       "double *p",
+							       "double *t",
 							       "int nx",
 							       "int ny",
 							       "int np"])
@@ -1242,6 +1263,7 @@ relation generate_initial_value_code:(DAELow.DAELow) => string =
 							     "double *xd", 
 							     "double *y",
 							     "double *p",
+							     "double *t",
 							     "int nx",
 							     "int ny",
 							     "int np"])

c_runtime/simulation_runtime.cpp

@@ -199,7 +199,7 @@ int euler_main(int argc,char** argv) {
   // Calculate initial values from (fixed) start attributes and intial equation
   // sections
   init=1;
-  initial_function(x,xd,y,p,nx,ny,np);
+  initial_function(x,xd,y,p,&t,nx,ny,np);
   init=0;
 
   int npts_per_result=int((stop-start)/(step*(numpoints-2)));
@@ -296,7 +296,7 @@ int dassl_main(int argc, char **argv)
   // Calculate initial values from (fixed) start attributes and intial equation
   // sections
   init=1;
-  initial_function(x,xd,y,p,nx,ny,np);
+  initial_function(x,xd,y,p,&t,nx,ny,np);
   t=start;
   tout = newTime(t, step); // TODO: check time events here. Maybe dassl should not be allowed to simulate past the scheduled time event.
 

c_runtime/simulation_runtime.h

@@ -113,6 +113,7 @@ int functionODE(double *x, double *xd, double *y, double *p,
 // function for calculate initial values from initial equations
 // and fixed start attibutes
 int initial_function(double*x, double *xd, double*y, double*p,
+		     double *t,
 		    int nx, int ny, int np); 
 
 // Adds a result to the simulation result data.