/[ascend]/trunk/base/generic/integrator/idaanalyse.c
ViewVC logotype

Contents of /trunk/base/generic/integrator/idaanalyse.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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