/[ascend]/trunk/solvers/ida/idaanalyse.c
ViewVC logotype

Contents of /trunk/solvers/ida/idaanalyse.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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