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

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