/[ascend]/branches/ksenija2/solvers/ida/idaanalyse.c
ViewVC logotype

Contents of /branches/ksenija2/solvers/ida/idaanalyse.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2855 - (show annotations) (download) (as text)
Sat Mar 21 02:53:28 2015 UTC (3 years, 11 months ago) by jpye
File MIME type: text/x-csrc
File size: 30116 byte(s)
fixing error introduced earlier!

1 /* ASCEND modelling environment
2 Copyright (C) 2006-2011 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>.
16 *//** @file
17 Analysis routines for the ASCEND wrapper of the IDA integrator.
18 These functions perform sorting of variables and relations and create
19 additional lists of variables as required for use by our ida.c code.
20 */
21 #include "idaanalyse.h"
22 #include "idatypes.h"
23 #include "idaio.h"
24
25 #include <ascend/general/panic.h>
26 #include <ascend/utilities/error.h>
27
28 #include <ascend/linear/linsolqr.h>
29
30 #include <ascend/system/diffvars.h>
31 #include <ascend/system/slv_stdcalls.h>
32 #include <ascend/system/block.h>
33 #include <ascend/system/diffvars.h>
34 #include <ascend/system/diffvars_impl.h>
35 #include <ascend/system/jacobian.h>
36 #include <ascend/system/cond_config.h>
37 #include <ascend/solver/slvDOF.h>
38
39 #define ANALYSE_DEBUG
40
41 /*
42 define DERIV_WITHOUT_DIFF to enable experimental handling of derivatives
43 for which corresponding differential vars were not found to be incident.
44
45 (this is important! it means that if we have eg "Qdot = der(Q)" in our model,
46 but Q itself isn't referred to anywhere. We want IDA to include Q in our
47 system in that case, because we want to calculate the integral of Qdot, Q,
48 as time marches forward. -- JP)
49 */
50 #define DERIV_WITHOUT_DIFF
51
52 #define VARMSG(MSG) \
53 varname = var_make_name(integ->system,v); \
54 CONSOLE_DEBUG(MSG,varname); \
55 ASC_FREE(varname)
56
57 //static int integrator_ida_check_partitioning(IntegratorSystem *integ);
58 static int integrator_ida_check_diffindex(IntegratorSystem *integ);
59 /* static int integrator_ida_rebuild_diffindex(IntegratorSystem *integ); */
60
61 const var_filter_t integrator_ida_nonderiv = {
62 VAR_SVAR | VAR_ACTIVE | VAR_FIXED | VAR_DERIV,
63 VAR_SVAR | VAR_ACTIVE | 0 | 0
64 };
65
66 const var_filter_t integrator_ida_deriv = {
67 VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | VAR_FIXED | VAR_DERIV,
68 VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | 0 | VAR_DERIV
69 };
70
71 const rel_filter_t integrator_ida_rel = {
72 REL_INCLUDED | REL_EQUALITY | REL_ACTIVE,
73 REL_INCLUDED | REL_EQUALITY | REL_ACTIVE
74 };
75
76 /**
77 This is the first step in the DAE analysis process. We inspect the
78 derivative chains from the solver's analyse() routine.
79 Check derivative chains: if a var or derivative is inelligible,
80 follow down the chain fixing & setting zero any derivatives.
81 */
82 static int integrator_ida_check_vars(IntegratorSystem *integ){
83 const SolverDiffVarCollection *diffvars;
84 char *varname;
85 int n_y = 0;
86 int i, j, k, n;
87 struct var_variable *v;
88 SolverDiffVarSequence seq;
89 int vok;
90 struct var_variable **vlist, **oldvl;
91 int var_found;
92
93 #ifdef ANALYSE_DEBUG
94 CONSOLE_DEBUG("BEFORE CHECKING VARS");
95 system_diffvars_debug(integ->system,stderr);
96 #endif
97
98 /* we shouldn't have allocated these yet: just be sure */
99 asc_assert(integ->y==NULL);
100 asc_assert(integ->ydot==NULL);
101 asc_assert(integ->n_y==0);
102
103 /* get the the dervative chains from the system */
104 diffvars = system_get_diffvars(integ->system);
105 if(diffvars==NULL){
106 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");
107 return 1;
108 }
109
110 #ifdef ANALYSE_DEBUG
111 CONSOLE_DEBUG("AFTER system_get_diffvars...");
112 system_var_list_debug(integ->system);
113 #endif
114
115 /* add the variables from the derivative chains */
116 for(i=0; i<diffvars->nseqs; ++i){
117 seq = diffvars->seqs[i];
118 asc_assert(seq.n >= 1);
119 v = seq.vars[0];
120
121 /* update the var_fixed flag */
122 (void)var_fixed(v);
123 vok = var_apply_filter(v,&integrator_ida_nonderiv);
124 if(!var_incident(v) || !var_active(v)){
125 /* VARMSG("good var '%s' is not incident"); */
126 /* var meets our filter, but perhaps it's not incident? */
127 if(seq.n == 1 || var_apply_filter(seq.vars[1],&integrator_ida_nonderiv)){
128 /* VARMSG("DEACTIVATING NON-INCIDENT VAR '%s' (NO DERIVATIVE)"); */
129 var_set_active(v,0);
130 vok = 0;
131 }else{
132 var_found = 0;
133 #ifdef DERIV_WITHOUT_DIFF
134 VARMSG("'%s' has a derivative present, so needs to be included in the system");
135 #ifdef ANALYSE_DEBUG
136 CONSOLE_DEBUG("That var %s active",(var_active(v) ? "is" : "is NOT"));
137 #endif
138 var_set_active(v,1);
139 vok = var_apply_filter(v,&integrator_ida_nonderiv);
140 if (vok) {
141 if (!var_incident(v)) var_set_incident(v,1);
142 n = slv_get_num_solvers_vars(integ->system);
143 vlist = ASC_NEW_ARRAY(struct var_variable*,n + 2);
144 oldvl = slv_get_solvers_var_list(integ->system);
145 for(k = 0; k < n; k++) {
146 vlist[k] = oldvl[k];
147 if (vlist[k]==v) var_found = 1;
148 }
149 vlist[n] = v;
150 vlist[n+1] = NULL;
151 if (var_found) {
152 slv_set_solvers_var_list(integ->system,oldvl,n);
153 ascfree(vlist);
154 }
155 else {
156 slv_set_solvers_var_list(integ->system,vlist,n+1);
157 ascfree(oldvl);
158 }
159 }
160 #else
161 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Non-incident var with an incident derivative. ASCEND can't handle this case at the moment, but we hope to fix it.");
162 return 1;
163 #endif
164 }
165 }
166
167 if(!vok){
168 #ifdef ANALYSE_DEBUG
169 if(var_fixed(v)){
170 VARMSG("Fixed variable '%s' fails non-deriv filter");
171 }else{
172 VARMSG("'%s' fails non-deriv filter");
173 }
174 CONSOLE_DEBUG("passes nonderiv? %s (flags = 0x%x)"
175 , (var_apply_filter(v,&integrator_ida_nonderiv) ? "TRUE" : "false")
176 , var_flags(v)
177 );
178 #endif
179 for(j=1;j<seq.n;++j){
180 v = seq.vars[j];
181 var_set_active(v,FALSE);
182 var_set_value(v,0);
183 #ifdef ANALYSE_DEBUG
184 VARMSG("Derivative '%s' SET ZERO AND INACTIVE");
185 #endif
186 }
187 continue;
188 }
189
190 /* VARMSG("We will use var '%s'"); */
191 if(seq.n > 2){
192 varname = var_make_name(integ->system,v);
193 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Too-long derivative chain with var '%s'",varname);
194 ASC_FREE(varname);
195 return 2;
196 }else if(seq.n==2){
197 /* diff var */
198 if(var_apply_filter(seq.vars[1],&integrator_ida_deriv)){
199 /* add the diff & deriv vars to the lists */
200 n_y++;
201 /* VARMSG("Added diff var '%s'");
202 v = seq.vars[1]; VARMSG("and its derivative '%s'"); */
203 continue;
204 }
205 /* VARMSG("Diff var '%s' being converted to alg var...");
206 v = seq.vars[1]; VARMSG("...because deriv var '%s' fails filter"); */
207 /* fall through */
208 }
209
210 /* VARMSG("Adding '%s' to algebraic"); */
211 n_y++;
212 }
213
214 #ifdef ANALYSE_DEBUG
215 system_var_list_debug(integ->system);
216 #endif
217
218 /* we assert that all vars in y meet the integrator_ida_nonde6: v_wind
219 riv filter */
220 /* we assert that all vars in ydot meet the integrator_ida_deriv filter */
221
222 #ifdef ANALYSE_DEBUG
223 CONSOLE_DEBUG("Found %d good non-derivative vars", n_y);
224 #endif
225 integ->n_y = n_y;
226
227 return 0;
228 }
229
230 /**
231 Flag relations that contain derivatives as 'REL_DIFFERENTIAL'
232
233 @TODO what to do about relations that make reference to 't'?
234 */
235 static int integrator_ida_flag_rels(IntegratorSystem *integ){
236 int i, n, c, nd=0;
237 struct rel_relation **rels;
238 n = slv_get_num_solvers_rels(integ->system);
239 rels = slv_get_solvers_rel_list(integ->system);
240 for(i=0;i<n;++i){
241 c = rel_classify_differential(rels[i]);
242 if(c){
243 nd++;
244 /* CONSOLE_DEBUG("Rel %d is DIFFERENTIAL", i); */
245 }
246 }
247 #ifdef ANALYSE_DEBUG
248 CONSOLE_DEBUG("Found %d differential equations (so %d algebraic)",nd, n - nd);
249 #endif
250 integ->n_diffeqs = nd;
251 return 0;
252 }
253
254 /**
255 Sort the lists. First we will put the non-derivative vars, then we will
256 put the derivative vars. Then we will put all the others.
257
258 See http://ascend4.org/File:Diffvars.png
259 */
260 static int integrator_ida_sort_rels_and_vars(IntegratorSystem *integ){
261 int ny1, nydot, nr;
262
263
264 #ifdef ANALYSE_DEBUG
265 CONSOLE_DEBUG("BEFORE SORTING RELS AND VARS");
266 system_diffvars_debug(integ->system,stderr);
267 #endif
268
269 /* we should not have allocated y or ydot yet */
270 asc_assert(integ->y==NULL && integ->ydot==NULL);
271
272 /* but we should have found some variables (and know how many) */
273 asc_assert(integ->n_y);
274
275 if(system_cut_vars(integ->system, 0, &integrator_ida_nonderiv, &ny1)){
276 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting non-derivs");
277 return 1;
278 }
279
280 #ifdef ANALYSE_DEBUG
281 CONSOLE_DEBUG("Cut %d non-derivative vars to start of list. cf integ->n_y = %d",ny1,integ->n_y);
282 #endif
283 asc_assert(ny1 == integ->n_y);
284
285 ERROR_REPORTER_HERE(ASC_USER_NOTE,"moving derivs to start of remainder\n");
286 if(system_cut_vars(integ->system, ny1, &integrator_ida_deriv, &nydot)){
287 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting derivs");
288 return 1;
289 }
290
291 if(system_cut_rels(integ->system, 0, &integrator_ida_rel, &nr)){
292 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting derivs");
293 return 1;
294 }
295
296 if(ny1 != nr){
297 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem is not square (ny = %d, nr = %d)",ny1,nr);
298 return 2;
299 }
300
301 return 0;
302 }
303
304 /**
305 Build up the lists 'ydot' and 'y_id.
306 'ydot' allows us to find the derivative of a variable given its sindex.
307 'y_id' allows us to locate the variable of a derivative given its sindex.
308
309 Hence
310 y_id[var_sindex(derivvar)-integ->n_y] = var_sindex(diffvar)
311 ydot[var_sindex(diffvar)] = derivvar (or NULL if diffvar is algebraic)
312
313 Note that there is 'shifting' required when looking up y_id.
314
315 Note that we want to get rid of 'y' but for the moment we are also assigning
316 that. We want to be rid of it because if the vars are ordered correctly
317 it shouldn't be needed.
318 */
319 static int integrator_ida_create_lists(IntegratorSystem *integ){
320 const SolverDiffVarCollection *diffvars;
321 int i, j;
322 struct var_variable *v;
323
324 SolverDiffVarSequence seq;
325
326 asc_assert(integ->y==NULL);
327 asc_assert(integ->ydot==NULL);
328 asc_assert(integ->y_id== NULL);
329
330 /* get the the dervative chains from the system */
331 diffvars = system_get_diffvars(integ->system);
332 if(diffvars==NULL){
333 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");
334 return 1;
335 }
336
337 /* allocate space (more than enough) */
338 integ->y = ASC_NEW_ARRAY(struct var_variable *,integ->n_y);
339 integ->ydot = ASC_NEW_ARRAY_CLEAR(struct var_variable *,integ->n_y);
340 integ->n_ydot = 0;
341 for(i=0;i<integ->n_y;++i){
342 asc_assert(integ->ydot[i] == 0);
343 }
344
345 #ifdef ANALYSE_DEBUG
346 CONSOLE_DEBUG("Passing through chains...");
347 #endif
348 /* create the lists y and ydot, ignoring 'bad' vars */
349 for(i=0; i<diffvars->nseqs; ++i){
350 /* CONSOLE_DEBUG("i = %d",i); */
351
352 seq = diffvars->seqs[i];
353 asc_assert(seq.n >= 1);
354 v = seq.vars[0];
355 j = var_sindex(v);
356
357 if(!var_apply_filter(v,&integrator_ida_nonderiv)){
358 continue;
359 }
360
361 integ->y[j] = v;
362 /* VARMSG("'%s' is good non-deriv"); */
363
364 if(seq.n > 1 && var_apply_filter(seq.vars[1],&integrator_ida_deriv)){
365 v = seq.vars[1];
366 asc_assert(var_sindex(v) >= integ->n_y);
367 asc_assert(var_sindex(v)-integ->n_y < integ->n_y);
368
369 integ->ydot[j] = v;
370 integ->n_ydot++;
371 /* VARMSG("'%s' is good deriv"); */
372 }else{
373 asc_assert(integ->ydot[j]==NULL);
374 }
375 }
376
377 #ifdef ANALYSE_DEBUG
378 CONSOLE_DEBUG("Found %d good non-derivs",j);
379 #endif
380 /* create the list y_id by looking at non-NULLs from ydot */
381 integ->y_id = ASC_NEW_ARRAY(int,integ->n_ydot);
382 for(i=0,j=0; i < integ->n_y; ++i){
383 if(integ->ydot[i]==NULL)continue;
384 /* v = integ->ydot[i]; VARMSG("deriv '%s'..."); */
385 /* v = integ->y[i]; VARMSG("diff '%s'..."); */
386 integ->y_id[var_sindex(integ->ydot[i]) - integ->n_y] = i;
387 j++;
388 }
389
390 asc_assert(j==integ->n_ydot);
391
392 return 0;
393 }
394
395 /**
396 Check for index-1 DAE system. We plan to do this by checking that df/dyd'
397 and dg/dya are both non-singular (and of course square). Valid systems
398 can (we think?) be written that don't meet these requirements but they
399 will have index problems and may not solve with IDA.
400 */
401 int integrator_ida_check_index(IntegratorSystem *integ){
402 #if 1
403 linsolqr_system_t L;
404 mtx_range_t range;
405 mtx_region_t R;
406 int res, r;
407 struct SystemJacobianStruct df_dydp, dg_dya;
408
409 #ifdef ANALYSE_DEBUG
410 CONSOLE_DEBUG("system has total of %d rels and %d vars"
411 ,slv_get_num_solvers_rels(integ->system)
412 ,slv_get_num_solvers_vars(integ->system)
413 );
414
415 CONSOLE_DEBUG("VAR_DERIV = 0x%x = %d",VAR_DERIV, VAR_DERIV);
416 CONSOLE_DEBUG("system_vfilter_deriv.matchbits = 0x%x",system_vfilter_deriv.matchbits);
417 CONSOLE_DEBUG("system_vfilter_deriv.matchvalue= 0x%x",system_vfilter_deriv.matchvalue);
418 #endif
419
420 asc_assert(system_vfilter_deriv.matchbits & VAR_DERIV);
421 asc_assert(system_vfilter_deriv.matchvalue & VAR_DERIV);
422
423 #ifdef ANALYSE_DEBUG
424 CONSOLE_DEBUG("system has %d vars matching deriv filter",slv_count_solvers_vars(integ->system, &system_vfilter_deriv));
425 #endif
426
427 res = system_jacobian(integ->system
428 , &system_rfilter_diff
429 , &system_vfilter_deriv
430 , 1 /* safe evaluation */
431 , &df_dydp
432 );
433
434 if(res){
435 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error calculating df/dyd'");
436 }
437
438 #ifdef ANALYSE_DEBUG
439 CONSOLE_DEBUG("df/dyd': nr = %d, nv = %d",df_dydp.n_rels,df_dydp.n_vars);
440 #endif
441
442 res = system_jacobian(integ->system
443 , &system_rfilter_algeb
444 , &system_vfilter_algeb
445 , 1 /* safe evaluation */
446 , &dg_dya
447 );
448
449 if(res){
450 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error calculating dg/dya");
451 }
452
453 #ifdef ANALYSE_DEBUG
454 CONSOLE_DEBUG("dg/dya: nr = %d, nv = %d",dg_dya.n_rels,dg_dya.n_vars);
455 #endif
456
457 if((df_dydp.n_rels == 0) ^ (df_dydp.n_vars == 0)){
458 ERROR_REPORTER_HERE(ASC_PROG_ERR,"df/dyd' is a bit ambiguous");
459 }
460
461 if(dg_dya.n_rels <= 0){
462 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"No algebraic equations were found in the DAE system!");
463 }else if(dg_dya.n_rels != dg_dya.n_vars){
464 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"The algebraic part of the DAE jacobian, dg/dya, is not square!");
465 }else{
466 /* check the rank */
467 range.low = 0; range.high = mtx_order(dg_dya.M) - 1;
468 R.row = range; R.col = range;
469
470 L = linsolqr_create_default();
471 linsolqr_set_matrix(L,dg_dya.M);
472 linsolqr_set_region(L,R);
473 linsolqr_prep(L,linsolqr_fmethod_to_fclass(linsolqr_fmethod(L)));
474 linsolqr_reorder(L, &R, linsolqr_rmethod(L));
475 linsolqr_factor(L,linsolqr_fmethod(L));
476 r = linsolqr_rank(L);
477
478 linsolqr_destroy(L);
479
480 if(r != dg_dya.n_rels){
481 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Your DAE system has an index problem: the matrix dg/dya is not full rank");
482 }
483 }
484
485 ASC_FREE(dg_dya.vars);
486 ASC_FREE(dg_dya.rels);
487 mtx_destroy(dg_dya.M);
488
489 if(df_dydp.n_rels <= 0){
490 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"No differential equations were found in the DAE system!");
491 }else if(df_dydp.n_rels != df_dydp.n_vars){
492 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"The differential part of the the jacobian dg/dya is not square!");
493 ASC_FREE(df_dydp.vars);
494 ASC_FREE(df_dydp.rels);
495 mtx_destroy(df_dydp.M);
496 return 1;
497 }else{
498 /* check the rank */
499 range.low = 0; range.high = mtx_order(df_dydp.M) - 1;
500 R.row = range; R.col = range;
501
502 L = linsolqr_create_default();
503 linsolqr_set_matrix(L,df_dydp.M);
504 linsolqr_set_region(L,R);
505 linsolqr_prep(L,linsolqr_fmethod_to_fclass(linsolqr_fmethod(L)));
506 linsolqr_reorder(L, &R, linsolqr_rmethod(L));
507 linsolqr_factor(L,linsolqr_fmethod(L));
508 r = linsolqr_rank(L);
509
510 linsolqr_destroy(L);
511
512 if(r != df_dydp.n_rels){
513 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Your DAE system has an index problem: the matrix df/dyd' is not full rank");
514 }
515 }
516
517 if(df_dydp.n_rels + dg_dya.n_rels == 0){
518 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Both df/dyd' and dg/dya were empty!");
519 }
520
521 ASC_FREE(df_dydp.vars);
522 ASC_FREE(df_dydp.rels);
523 mtx_destroy(df_dydp.M);
524 return 0;
525 #else
526 ERROR_REPORTER_HERE(ASC_PROG_ERR,"check_index disabled");
527 return 0;
528 #endif
529 }
530
531 /*------------------------------------------------------------------------------
532 ANALYSIS ROUTINE (new implementation)
533 */
534
535 /**
536 Perform additional problem analysis to prepare problem for integration with
537 IDA.
538
539 We assume that the analyse_generate_diffvars routine has been called, so
540 that we just need to call slv_get_diffvars to access the derivative
541 chains.
542
543 We can also assume that the independent variable has been found.
544
545 See mailing list, ~Jan 2007.
546
547 Note, the stuff for identifying the static and output sub-problems should
548 be part of integrator.c, not this file. We will assume this is handled
549
550 @return 0 on success
551 @see integrator_analyse
552 */
553 int integrator_ida_analyse(IntegratorSystem *integ){
554 int res;
555 const SolverDiffVarCollection *diffvars;
556 int i;
557 #ifdef ANALYSE_DEBUG
558 char *varname;
559 #endif
560
561 asc_assert(integ->engine==INTEG_IDA);
562
563 #ifdef ANALYSE_DEBUG
564 CONSOLE_DEBUG("System contains a total of %d bnds and %d rels"
565 ,slv_get_num_solvers_bnds(integ->system)
566 ,slv_get_num_solvers_rels(integ->system)
567 );
568
569 /* set the active flags on variables depending on the state of WHENs */
570 CONSOLE_DEBUG("Currently %d rels active",slv_count_solvers_rels(integ->system, &integrator_ida_rel));
571 #endif
572
573 reanalyze_solver_lists_cont(integ->system);
574
575 #ifdef ANALYSE_DEBUG
576 CONSOLE_DEBUG("After analysing WHENs, there are %d rels active"
577 ,slv_count_solvers_rels(integ->system, &integrator_ida_rel)
578 );
579 #endif
580
581
582 #ifdef ANALYSE_DEBUG
583 CONSOLE_DEBUG("Starting IDA analysis");
584 #endif
585
586 /* set the flags on differential and derivative and algebraic vars */
587 res = integrator_ida_check_vars(integ);
588 if(res){
589 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem with vars");
590 return 1;
591 }
592
593 res = integrator_ida_flag_rels(integ);
594 if(res){
595 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem classifying differential equations");
596 return 1;
597 }
598
599 #ifdef ANALYSE_DEBUG
600 CONSOLE_DEBUG("Sorting rels and vars");
601 #endif
602
603 res = integrator_ida_sort_rels_and_vars(integ);
604 if(res){
605 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem sorting rels and vars");
606 return 1;
607 }
608
609 #ifdef ANALYSE_DEBUG
610 CONSOLE_DEBUG("Creating lists");
611
612 CONSOLE_DEBUG("BEFORE MAKING LISTS");
613 integrator_ida_debug(integ,stderr);
614 #endif
615
616 res = integrator_ida_create_lists(integ);
617 if(res){
618 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem creating lists");
619 return 1;
620 }
621
622 #ifdef ANALYSE_DEBUG
623 CONSOLE_DEBUG("After ida_create_lists, there are %d rels active"
624 ,slv_count_solvers_rels(integ->system, &integrator_ida_rel)
625 );
626 CONSOLE_DEBUG("Checking lists");
627 #endif
628
629 #ifdef ANALYSE_DEBUG
630 asc_assert(integ->y);
631 asc_assert(integ->ydot);
632 asc_assert(integ->y_id);
633
634 integrator_ida_debug(integ,stderr);
635 #endif
636
637 if(integrator_ida_check_diffindex(integ)){
638 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error with diffindex");
639 return 360;
640 }
641
642
643 #ifdef ANALYSE_DEBUG
644 CONSOLE_DEBUG("After ida_check_diffindex, there are %d rels active"
645 ,slv_count_solvers_rels(integ->system, &integrator_ida_rel)
646 );
647 #endif
648
649 /* the following causes TestIDA.testincidence3 to fail:
650 if(integ->n_y != slv_get_num_solvers_rels(integ->system)){
651 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Problem is not square");
652 return 2;
653 }*/
654
655 #if 0
656 /* WARNING: uncommenting this section seems to break the solver!!!! -- JP */
657 /* check structural singularity for the two IDACalcIC scenarios */
658
659 /* ...(1) FIX the derivatives */
660 CONSOLE_DEBUG("Checking system with derivatives fixed...");
661 for(i=0;i<integ->n_y;++i){
662 if(integ->ydot[i])var_set_fixed(integ->ydot[i],1);
663 }
664
665 slv_block_partition(integ->system);
666 res = integrator_ida_block_check(integ);
667 if(res)return 100 + res;
668
669 /* ...(2) FREE the derivatives, FIX the diffvars */
670 CONSOLE_DEBUG("Checking system with differential variables fixed...");
671 for(i=0;i<integ->n_y;++i){
672 if(integ->ydot[i]){
673 var_set_fixed(integ->ydot[i],0);
674 var_set_fixed(integ->y[i],1);
675 }
676 }
677 slv_block_partition(integ->system);
678 res = integrator_ida_block_check(integ);
679 if(res)return 200 + res;
680
681 /* ...(3) restore as it was, FREE the diffvars */
682 for(i=0;i<integ->n_y;++i){
683 if(integ->ydot[i]){
684 var_set_fixed(integ->y[i],0);
685 }
686 }
687 #endif
688
689 res = integrator_ida_check_index(integ);
690 if(res){
691 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Your DAE system has an index problem");
692 return 100 + res;
693 }
694
695 #ifdef ANALYSE_DEBUG
696 CONSOLE_DEBUG("After ida_check_index, there are %d rels active"
697 ,slv_count_solvers_rels(integ->system, &integrator_ida_rel)
698 );
699 #endif
700
701
702 /* get the the dervative chains from the system */
703 diffvars = system_get_diffvars(integ->system);
704
705 /* check the indep var is same as was located earlier */
706 if(diffvars->nindep>1){
707 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Multiple variables specified as independent (ode_type=-1)");
708 return 3;
709 }else if(diffvars->nindep<1){
710 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Independent var not set (ode_type=-1)");
711 return 4;
712 }else if(diffvars->indep[0]!=integ->x){
713 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Indep var doesn't match");
714 return 5;
715 }
716
717 #ifdef ANALYSE_DEBUG
718 CONSOLE_DEBUG("Collecting observed variables");
719 #endif
720
721 /* get the observations */
722 /** @TODO this should use a 'problem provider' API instead of static data
723 from the system object */
724 integ->n_obs = diffvars->nobs;
725 integ->obs = ASC_NEW_ARRAY(struct var_variable *,integ->n_obs);
726 for(i=0;i<integ->n_obs;++i){
727 /* we get them all, regardless of flags etc */
728 integ->obs[i] = diffvars->obs[i];
729 #ifdef ANALYSE_DEBUG
730 varname = var_make_name(integ->system,integ->obs[i]);
731 CONSOLE_DEBUG("'%s' is observation",varname);
732 ASC_FREE(varname);
733 #endif
734 }
735
736 integ->n_dobs = diffvars->ndobs;
737 integ->dobs = ASC_NEW_ARRAY(struct dis_discrete *,integ->n_dobs);
738 for(i=0;i<integ->n_dobs;++i){
739 /* we get them all, regardless of flags etc */
740 integ->dobs[i] = diffvars->dobs[i];
741 #ifdef ANALYSE_DEBUG
742 varname = dis_make_name(integ->system,integ->dobs[i]);
743 CONSOLE_DEBUG("'%s' is observation",varname);
744 ASC_FREE(varname);
745 #endif
746 }
747
748 #ifdef ANALYSE_DEBUG
749 CONSOLE_DEBUG("rels matchbits: 0x%x",integrator_ida_rel.matchbits);
750 CONSOLE_DEBUG("rels matchvalue: 0x%x",integrator_ida_rel.matchvalue);
751
752 CONSOLE_DEBUG("At the end of ida_analyse, there are %d rels active"
753 ,slv_count_solvers_rels(integ->system, &integrator_ida_rel)
754 );
755 #endif
756
757 return 0;
758 }
759
760 #if 0 /* disused function, apparently -- JP */
761 static int integrator_ida_check_partitioning(IntegratorSystem *integ){
762 long i, nv;
763 int err=0;
764 char *varname;
765 struct var_variable **vlist, *v;
766 nv = slv_get_num_solvers_vars(integ->system);
767 vlist = slv_get_solvers_var_list(integ->system);
768 var_filter_t vf = {VAR_SVAR|VAR_ACTIVE|VAR_INCIDENT|VAR_FIXED
769 ,VAR_SVAR|VAR_ACTIVE|VAR_INCIDENT|0 };
770 for(i=0;i<nv;++i){
771 v=vlist[i];
772 if(!var_apply_filter(v,&vf))continue;
773 varname = var_make_name(integ->system,v);
774 if(!var_deriv(v)){
775 #ifdef ANALYSE_DEBUG
776 fprintf(stderr,"vlist[%ld] = '%s' (nonderiv)\n",i,varname);
777 #endif
778 if(i>=integ->n_y){
779 ERROR_REPORTER_HERE(ASC_PROG_ERR,"non-deriv var '%s' is not at the start",varname);
780 err++;
781 }
782 }else{
783 #ifdef ANALYSE_DEBUG
784 fprintf(stderr,"vlist[%ld] = '%s' (derivative)\n",i,varname);
785 #endif
786 if(i<integ->n_y){
787 ERROR_REPORTER_HERE(ASC_PROG_ERR,"deriv var '%s' is not at the end (n_y = %d, i = %d)"
788 ,varname, integ->n_y, i
789 );
790 err++;
791 }
792 }
793 ASC_FREE(varname);
794 }
795 if(err){
796 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error in IDA partitioning");
797 return 1;
798 }
799 return 0;
800 }
801 #endif
802
803 int integrator_ida_block_check(IntegratorSystem *integ){
804 int res;
805 int dof;
806 #ifdef ANALYSE_DEBUG
807 int *vlist, *vp, i, nv, nv_ok;
808 char *varname;
809 struct var_variable **solversvars;
810
811 var_filter_t vfilt = {
812 VAR_ACTIVE | VAR_INCIDENT | VAR_FIXED,
813 VAR_ACTIVE | VAR_INCIDENT | 0
814 };
815
816 nv = slv_get_num_solvers_vars(integ->system);
817 solversvars = slv_get_solvers_var_list(integ->system);
818 CONSOLE_DEBUG("-------------- nv = %d -------------",nv);
819 for(nv_ok=0, i=0; i < nv; ++i){
820 if(var_apply_filter(solversvars[i],&vfilt)){
821 varname = var_make_name(integ->system,solversvars[i]);
822 fprintf(stderr,"%s\n",varname);
823 ASC_FREE(varname);
824 nv_ok++;
825 }
826 }
827 CONSOLE_DEBUG("----------- got %d ok -------------",nv_ok);
828 #endif
829
830 if(!slvDOF_status(integ->system, &res, &dof)){
831 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to determine DOF status");
832 return -1;
833 }
834 switch(res){
835 case 1: CONSOLE_DEBUG("System is underspecified (%d degrees of freedom)",dof);break;
836 case 2: CONSOLE_DEBUG("System is square"); return 0; /* all OK */
837 case 3: CONSOLE_DEBUG("System is structurally singular"); break;
838 case 4: CONSOLE_DEBUG("System is overspecified"); break;
839 default:
840 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unrecognised slfDOF_status");
841 return -2;
842 }
843
844 #ifdef ANALYSE_DEBUG
845 /* if it was underspecified, what vars could be fixed? */
846 if(res==1){
847 CONSOLE_DEBUG("Need to FIX %d of the following vars:",dof);
848 solversvars = slv_get_solvers_var_list(integ->system);
849 if(!slvDOF_eligible(integ->system, &vlist)){
850 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to det slvDOF_eligble list");
851 return -3;
852 }
853 for(vp=vlist;*vp!=-1;++vp){
854 varname = var_make_name(integ->system, solversvars[*vp]);
855 CONSOLE_DEBUG("Fixable var: %s",varname);
856 ASC_FREE(varname);
857 }
858 CONSOLE_DEBUG("(Found %d fixable vars)",(int)(vp-vlist));
859 return 1;
860 }
861 #endif
862
863 return res;
864 }
865
866 /* return 0 on succes */
867 static int check_dups(IntegratorSystem *integ, struct var_variable **list,int n,int allownull){
868 int i,j;
869 struct var_variable *v;
870 #ifdef ANALYSE_DEBUG
871 char *varname;
872 #endif
873 for(i=0; i< n; ++i){
874 v=list[i];
875 if(v==NULL){
876 if(allownull)continue;
877 else return 2;
878 }
879 asc_assert(v!=(void *)0x31);
880 for(j=0; j<i-1;++j){
881 if(list[j]==NULL)continue;
882 if(v==list[j]){
883 #ifdef ANALYSE_DEBUG
884 varname = var_make_name(integ->system,v);
885 if(varname){
886 CONSOLE_DEBUG("Duplicate of '%s' found",varname);
887 ASC_FREE(varname);
888 }else{
889 CONSOLE_DEBUG("Duplicate found (couldn't retrieve name)");
890 }
891 ASC_FREE(varname);
892 #endif
893 return 1;
894 }
895 }
896 }
897 return 0;
898 }
899
900 /*
901 We are going to check that
902
903 - n_y in range (0,n_vars)
904 - no duplicates anywhere in the varlist
905 - no duplicates in non-NULL elements of ydot
906
907 - first n_y vars in solver's var list match those in vector y and match the non-deriv filter.
908 - var_sindex for first n_y vars match their position the the solver's var list
909
910 - ydot contains n_ydot non-NULL elements, each match the deriv filter.
911 - vars in solver's var list positions [n_y,n_y+n_ydot) all match deriv filter
912
913 - y_id contains n_ydot elements, with int values in range [0,n_y)
914 - var_sindex(ydot[y_id[i]]) - n_y == i for i in [0,n_ydot)
915
916 - all the vars [n_y+n_ydot,n_vars) fail the non-deriv filter and fail the deriv filter.
917 */
918 static int integrator_ida_check_diffindex(IntegratorSystem *integ){
919 int i, n_vars, n_ydot=0;
920 struct var_variable **list, *v;
921 char *varname;
922 const char *msg;
923
924 #ifdef ANALYSE_DEBUG
925 CONSOLE_DEBUG("Checking diffindex vector");
926 #endif
927
928 if(integ->y_id == NULL || integ->y == NULL || integ->ydot == NULL){
929 ERROR_REPORTER_HERE(ASC_PROG_ERR,"list(s) NULL");
930 return 1;
931 }
932
933 list = slv_get_solvers_var_list(integ->system);
934 n_vars = slv_get_num_solvers_vars(integ->system);
935
936 if(check_dups(integ, list, n_vars,FALSE)){
937 ERROR_REPORTER_HERE(ASC_PROG_ERR,"duplicates or NULLs in solver's var list");
938 return 1;
939 }
940
941 if(check_dups(integ, integ->ydot, integ->n_y,TRUE)){
942 ERROR_REPORTER_HERE(ASC_PROG_ERR,"duplicates in ydot vector");
943 return 1;
944 }
945
946 /* check n_y in range */
947 if(integ->n_y <= 0 || integ->n_y >= n_vars){
948 ERROR_REPORTER_HERE(ASC_PROG_ERR,"n_y = %d invalid (n_vars = %d)",integ->n_y, n_vars);
949 return 1;
950 }
951
952 /* check first n_y vars */
953 for(i=0; i < integ->n_y; ++i){
954 v = list[i];
955 if(!var_apply_filter(v, &integrator_ida_nonderiv)){
956 msg = "'%s' (in first n_y vars) fails non-deriv filter"; goto finish;
957 }else if(var_sindex(v) != i){
958 msg = "'%s' has wrong var_sindex"; goto finish;
959 }else if(v != integ->y[i]){
960 msg = "'%s' not matched in y vector"; goto finish;
961 }
962 /* meets filter, matches in y vector, has correct var_sindex. */
963 }
964
965 /* count non-NULLs in ydot */
966 for(i=0; i < integ->n_y; ++i){
967 v = integ->ydot[i];
968 if(v!=NULL){
969 if(var_sindex(v) < integ->n_y){
970 msg = "'%s' has var_sindex < n_y"; goto finish;
971 }else if(!var_apply_filter(v,&integrator_ida_deriv)){
972 msg = "'%s' (in next n_ydot vars) fails deriv filter"; goto finish;
973 }
974 /* lies beyond first n_y vars, match deriv filter */
975 n_ydot++;
976 }
977 }
978
979 /* check that vars [n_y, n_y+n_ydot) in solver's var list match the deriv filter */
980 for(i=integ->n_y; i< integ->n_y + n_ydot; ++i){
981 v = list[i];
982 if(!var_apply_filter(v,&integrator_ida_deriv)){
983 msg = "'%s', in range [n_y,n_y+n_ydot), fails deriv filter"; goto finish;
984 }
985 }
986
987 /* check values in y_id are ints int range [0,n_y), and that they point to correct vars */
988 for(i=0; i<n_ydot; ++i){
989 if(integ->y_id[i] < 0 || integ->y_id[i] >= integ->n_y){
990 ERROR_REPORTER_HERE(ASC_PROG_ERR,"value of y_id[%d] is out of range",i);
991 return 1;
992 }
993 v = integ->ydot[integ->y_id[i]];
994 if(var_sindex(v) - integ->n_y != i){
995 msg = "'%s' not index correctly from y_id"; goto finish;
996 }
997 }
998
999 /* check remaining vars fail both filters */
1000 for(i=integ->n_y + n_ydot; i<n_vars; ++i){
1001 v = list[i];
1002 if(var_apply_filter(v,&integrator_ida_nonderiv)){
1003 msg = "Var '%s' at end meets non-deriv filter, but shouldn't"; goto finish;
1004 }
1005 if(var_apply_filter(v,&integrator_ida_deriv)){
1006 CONSOLE_DEBUG("position = %d",i);
1007 msg = "Var '%s' at end meets deriv filter, but shouldn't"; goto finish;
1008 }
1009 }
1010
1011 return 0;
1012 finish:
1013 varname = var_make_name(integ->system,v);
1014 ERROR_REPORTER_HERE(ASC_PROG_ERR,msg,varname);
1015 ASC_FREE(varname);
1016 return 1;
1017 }
1018
1019 /**
1020 Given a derivative variable, return the index of its corresponding differential
1021 variable in the y vector (and equivalently the var_sindex of the diff var)
1022 */
1023 int integrator_ida_diffindex(const IntegratorSystem *integ, const struct var_variable *deriv){
1024 asc_assert(var_sindex(deriv) >= integ->n_y);
1025 asc_assert(var_sindex(deriv) < integ->n_y + integ->n_ydot);
1026 return integ->y_id[var_sindex(deriv) - integ->n_y];
1027 }
1028
1029 /**
1030 A less assertive version of diffindex...
1031 */
1032 int integrator_ida_diffindex1(const IntegratorSystem *integ, const struct var_variable *deriv){
1033 if(var_sindex(deriv) >= integ->n_y)return -1;
1034 if(var_sindex(deriv) < integ->n_y + integ->n_ydot)return -2;
1035 return integ->y_id[var_sindex(deriv) - integ->n_y];
1036 }
1037

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22