/[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 1317 - (show annotations) (download) (as text)
Mon Mar 5 14:11:37 2007 UTC (15 years, 3 months ago) by johnpye
File MIME type: text/x-csrc
File size: 24532 byte(s)
Added new 'system_jacobian' function for use by IDA (maybe elsewhere?)
Refactored the matrix output stuff in IDA.
Fixed the index checking in idaanalyse
Still a problem with checking rank of small matrices.
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 linsolqr_system_t L;
324 mtx_region_t R;
325 int res, r;
326 struct SystemJacobianStruct df_dydp, dg_dya;
327
328 CONSOLE_DEBUG("system has total of %d rels and %d vars"
329 ,slv_get_num_solvers_rels(sys->system)
330 ,slv_get_num_solvers_vars(sys->system)
331 );
332
333 CONSOLE_DEBUG("VAR_DERIV = 0x%x = %d",VAR_DERIV, VAR_DERIV);
334 CONSOLE_DEBUG("system_vfilter_deriv.matchbits = 0x%x",system_vfilter_deriv.matchbits);
335 CONSOLE_DEBUG("system_vfilter_deriv.matchvalue= 0x%x",system_vfilter_deriv.matchvalue);
336
337 asc_assert(system_vfilter_deriv.matchbits & VAR_DERIV);
338 asc_assert(system_vfilter_deriv.matchvalue & VAR_DERIV);
339
340 res = system_jacobian(sys->system
341 , &system_rfilter_diff
342 , &system_vfilter_deriv
343 , 1 /* safe evaluation */
344 , &df_dydp
345 );
346
347 if(res){
348 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error calculating df/dyd'");
349 }
350 CONSOLE_DEBUG("df/dyd': nr = %d, nv = %d",df_dydp.n_rels,df_dydp.n_vars);
351
352 res = system_jacobian(sys->system
353 , &system_rfilter_algeb
354 , &system_vfilter_algeb
355 , 1 /* safe evaluation */
356 , &dg_dya
357 );
358
359 if(res){
360 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error calculating dg/dya");
361 }
362 CONSOLE_DEBUG("dg/dya: nr = %d, nv = %d",dg_dya.n_rels,dg_dya.n_vars);
363
364 if((df_dydp.n_rels == 0) ^ (df_dydp.n_vars == 0)){
365 ERROR_REPORTER_HERE(ASC_PROG_ERR,"df/dyd' is a bit ambiguous");
366 }
367
368 if(dg_dya.n_rels <= 0){
369 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"No algebraic equations were found in the DAE system!");
370 }else if(dg_dya.n_rels != dg_dya.n_vars){
371 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"The algebraic part of the DAE jacobian, dg/dya, is not square!");
372 }else{
373 /* check the rank */
374 R.row.low = R.col.low = 0;
375 R.row.high = R.col.high = mtx_order(dg_dya.M) - 1;
376
377 L = linsolqr_create_default();
378 linsolqr_set_matrix(L,dg_dya.M);
379 linsolqr_set_region(L,R);
380 linsolqr_prep(L,linsolqr_fmethod_to_fclass(linsolqr_fmethod(L)));
381 linsolqr_reorder(L, &R, linsolqr_rmethod(L));
382 linsolqr_factor(L,linsolqr_fmethod(L));
383 r = linsolqr_rank(L);
384
385 linsolqr_destroy(L);
386
387 if(r != dg_dya.n_rels){
388 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Your DAE system has an index problem: the matrix dg/dya is not full rank");
389 }
390 }
391
392 if(df_dydp.n_rels <= 0){
393 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"No differential equations were found in the DAE system!");
394 }else if(df_dydp.n_rels != df_dydp.n_vars){
395 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"The differential part of the the jacobian dg/dya is not square!");
396 }else{
397 /* check the rank */
398 R.row.low = R.col.low = 0;
399 R.row.high = R.col.high = mtx_order(df_dydp.M) - 1;
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 != df_dydp.n_rels){
412 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Your DAE system has an index problem: the matrix df/dyd' is not full rank");
413 }
414 }
415
416 if(df_dydp.n_rels + dg_dya.n_rels == 0){
417 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Both df/dyd' and dg/dya were empty!");
418 }
419
420 ASC_FREE(df_dydp.vars);
421 ASC_FREE(df_dydp.rels);
422 ASC_FREE(dg_dya.vars);
423 ASC_FREE(dg_dya.rels);
424 mtx_destroy(dg_dya.M);
425 mtx_destroy(df_dydp.M);
426 return 0;
427 }
428
429 /**
430 Perform additional problem analysis to prepare problem for integration with
431 IDA.
432
433 We assume that the analyse_generate_diffvars routine has been called, so
434 that we just need to call slv_get_diffvars to access the derivative
435 chains.
436
437 We can also assume that the independent variable has been found.
438
439 See mailing list, ~Jan 2007.
440
441 Note, the stuff for identifying the static and output sub-problems should
442 be part of integrator.c, not this file. We will assume this is handled
443
444 @return 0 on success
445 @see integrator_analyse
446 */
447 int integrator_ida_analyse(IntegratorSystem *sys){
448 int res;
449 const SolverDiffVarCollection *diffvars;
450 int i;
451 #ifdef ANALYSE_DEBUG
452 char *varname;
453 #endif
454
455 asc_assert(sys->engine==INTEG_IDA);
456
457 #ifdef ANALYSE_DEBUG
458 CONSOLE_DEBUG("Starting IDA analysis");
459 #endif
460
461 res = integrator_ida_check_vars(sys);
462 if(res){
463 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem with vars");
464 return 1;
465 }
466
467 res = integrator_ida_flag_rels(sys);
468 if(res){
469 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem classifying differential equations");
470 return 1;
471 }
472
473 #ifdef ANALYSE_DEBUG
474 CONSOLE_DEBUG("Sorting rels and vars");
475 #endif
476
477 res = integrator_ida_sort_rels_and_vars(sys);
478 if(res){
479 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem sorting rels and vars");
480 return 1;
481 }
482
483 #ifdef ANALYSE_DEBUG
484 CONSOLE_DEBUG("Creating lists");
485
486 CONSOLE_DEBUG("BEFORE MAKING LISTS");
487 integrator_ida_debug(sys,stderr);
488 #endif
489
490 res = integrator_ida_create_lists(sys);
491 if(res){
492 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem creating lists");
493 return 1;
494 }
495
496 #ifdef ANALYSE_DEBUG
497 CONSOLE_DEBUG("Checking lists");
498
499 asc_assert(sys->y);
500 asc_assert(sys->ydot);
501 asc_assert(sys->y_id);
502
503 integrator_ida_debug(sys,stderr);
504 #endif
505
506 if(integrator_ida_check_diffindex(sys)){
507 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error with diffindex");
508 return 360;
509 }
510
511 /* the following causes TestIDA.testincidence3 to fail:
512 if(sys->n_y != slv_get_num_solvers_rels(sys->system)){
513 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Problem is not square");
514 return 2;
515 }*/
516
517 #if 0
518 /* check structural singularity for the two IDACalcIC scenarios */
519
520 /* ...(1) FIX the derivatives */
521 CONSOLE_DEBUG("Checking system with derivatives fixed...");
522 for(i=0;i<n_y;++i){
523 if(sys->ydot[i])var_set_fixed(sys->ydot[i],1);
524 }
525
526 slv_block_partition(sys->system);
527 res = integrator_ida_block_check(sys);
528 if(res)return 100 + res;
529
530 /* ...(2) FREE the derivatives, FIX the diffvars */
531 CONSOLE_DEBUG("Checking system with differential variables fixed...");
532 for(i=0;i<n_y;++i){
533 if(sys->ydot[i]){
534 var_set_fixed(sys->ydot[i],0);
535 var_set_fixed(sys->y[i],1);
536 }
537 }
538 slv_block_partition(sys->system);
539 res = integrator_ida_block_check(sys);
540 if(res)return 200 + res;
541
542 /* ...(3) restore as it was, FREE the diffvars */
543 for(i=0;i<n_y;++i){
544 if(sys->ydot[i]){
545 var_set_fixed(sys->y[i],0);
546 }
547 }
548 #endif
549
550 res = integrator_ida_check_index(sys);
551 if(res){
552 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Your DAE system has an index problem");
553 return 100 + res;
554 }
555
556 /* get the the dervative chains from the system */
557 diffvars = system_get_diffvars(sys->system);
558
559 /* check the indep var is same as was located earlier */
560 if(diffvars->nindep>1){
561 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Multiple variables specified as independent (ode_type=-1)");
562 return 3;
563 }else if(diffvars->nindep<1){
564 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Independent var not set (ode_type=-1)");
565 return 4;
566 }else if(diffvars->indep[0]!=sys->x){
567 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Indep var doesn't match");
568 return 5;
569 }
570
571 #ifdef ANALYSE_DEBUG
572 CONSOLE_DEBUG("Collecting observed variables");
573 #endif
574
575 /* get the observations */
576 /** @TODO this should use a 'problem provider' API instead of static data
577 from the system object */
578 sys->n_obs = diffvars->nobs;
579 sys->obs = ASC_NEW_ARRAY(struct var_variable *,sys->n_obs);
580 for(i=0;i<sys->n_obs;++i){
581 /* we get them all, regardless of flags etc */
582 sys->obs[i] = diffvars->obs[i];
583 #ifdef ANALYSE_DEBUG
584 varname = var_make_name(sys->system,sys->obs[i]);
585 CONSOLE_DEBUG("'%s' is observation",varname);
586 ASC_FREE(varname);
587 #endif
588 }
589
590 return 0;
591 }
592
593 /*------------------------------------------------------------------------------
594 ANALYSIS ROUTINE (new implementation)
595 */
596
597 static int integrator_ida_check_partitioning(IntegratorSystem *sys){
598 long i, nv;
599 int err=0;
600 char *varname;
601 struct var_variable **vlist, *v;
602 nv = slv_get_num_solvers_vars(sys->system);
603 vlist = slv_get_solvers_var_list(sys->system);
604 var_filter_t vf = {VAR_SVAR|VAR_ACTIVE|VAR_INCIDENT|VAR_FIXED
605 ,VAR_SVAR|VAR_ACTIVE|VAR_INCIDENT|0 };
606 for(i=0;i<nv;++i){
607 v=vlist[i];
608 if(!var_apply_filter(v,&vf))continue;
609 varname = var_make_name(sys->system,v);
610 if(!var_deriv(v)){
611 #ifdef ANALYSE_DEBUG
612 fprintf(stderr,"vlist[%ld] = '%s' (nonderiv)\n",i,varname);
613 #endif
614 if(i>=sys->n_y){
615 ERROR_REPORTER_HERE(ASC_PROG_ERR,"non-deriv var '%s' is not at the start",varname);
616 err++;
617 }
618 }else{
619 #ifdef ANALYSE_DEBUG
620 fprintf(stderr,"vlist[%ld] = '%s' (derivative)\n",i,varname);
621 #endif
622 if(i<sys->n_y){
623 ERROR_REPORTER_HERE(ASC_PROG_ERR,"deriv var '%s' is not at the end (n_y = %d, i = %d)"
624 ,varname, sys->n_y, i
625 );
626 err++;
627 }
628 }
629 ASC_FREE(varname);
630 }
631 if(err){
632 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error in IDA partitioning");
633 return 1;
634 }
635 return 0;
636 }
637
638 int integrator_ida_block_check(IntegratorSystem *sys){
639 int res;
640 int dof;
641 #ifdef ANALYSE_DEBUG
642 int *vlist, *vp, i, nv, nv_ok;
643 char *varname;
644 struct var_variable **solversvars;
645
646 var_filter_t vfilt = {
647 VAR_ACTIVE | VAR_INCIDENT | VAR_FIXED,
648 VAR_ACTIVE | VAR_INCIDENT | 0
649 };
650
651 nv = slv_get_num_solvers_vars(sys->system);
652 solversvars = slv_get_solvers_var_list(sys->system);
653 CONSOLE_DEBUG("-------------- nv = %d -------------",nv);
654 for(nv_ok=0, i=0; i < nv; ++i){
655 if(var_apply_filter(solversvars[i],&vfilt)){
656 varname = var_make_name(sys->system,solversvars[i]);
657 fprintf(stderr,"%s\n",varname);
658 ASC_FREE(varname);
659 nv_ok++;
660 }
661 }
662 CONSOLE_DEBUG("----------- got %d ok -------------",nv_ok);
663 #endif
664
665 if(!slvDOF_status(sys->system, &res, &dof)){
666 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to determine DOF status");
667 return -1;
668 }
669 switch(res){
670 case 1: CONSOLE_DEBUG("System is underspecified (%d degrees of freedom)",dof);break;
671 case 2: CONSOLE_DEBUG("System is square"); return 0; /* all OK */
672 case 3: CONSOLE_DEBUG("System is structurally singular"); break;
673 case 4: CONSOLE_DEBUG("System is overspecified"); break;
674 default:
675 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unrecognised slfDOF_status");
676 return -2;
677 }
678
679 #ifdef ANALYSE_DEBUG
680 /* if it was underspecified, what vars could be fixed? */
681 if(res==1){
682 CONSOLE_DEBUG("Need to FIX %d of the following vars:",dof);
683 solversvars = slv_get_solvers_var_list(sys->system);
684 if(!slvDOF_eligible(sys->system, &vlist)){
685 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to det slvDOF_eligble list");
686 return -3;
687 }
688 for(vp=vlist;*vp!=-1;++vp){
689 varname = var_make_name(sys->system, solversvars[*vp]);
690 CONSOLE_DEBUG("Fixable var: %s",varname);
691 ASC_FREE(varname);
692 }
693 CONSOLE_DEBUG("(Found %d fixable vars)",(int)(vp-vlist));
694 return 1;
695 }
696 #endif
697
698 return res;
699 }
700
701 /* return 0 on succes */
702 static int check_dups(IntegratorSystem *sys, struct var_variable **list,int n,int allownull){
703 int i,j;
704 struct var_variable *v;
705 #ifdef ANALYSE_DEBUG
706 char *varname;
707 #endif
708 for(i=0; i< n; ++i){
709 v=list[i];
710 if(v==NULL){
711 if(allownull)continue;
712 else return 2;
713 }
714 asc_assert(v!=(void *)0x31);
715 for(j=0; j<i-1;++j){
716 if(list[j]==NULL)continue;
717 if(v==list[j]){
718 #ifdef ANALYSE_DEBUG
719 varname = var_make_name(sys->system,v);
720 if(varname){
721 CONSOLE_DEBUG("Duplicate of '%s' found",varname);
722 ASC_FREE(varname);
723 }else{
724 CONSOLE_DEBUG("Duplicate found (couldn't retrieve name)");
725 }
726 ASC_FREE(varname);
727 #endif
728 return 1;
729 }
730 }
731 }
732 return 0;
733 }
734
735 /*
736 We are going to check that
737
738 - n_y in range (0,n_vars)
739 - no duplicates anywhere in the varlist
740 - no duplicates in non-NULL elements of ydot
741
742 - first n_y vars in solver's var list match those in vector y and match the non-deriv filter.
743 - var_sindex for first n_y vars match their position the the solver's var list
744
745 - ydot contains n_ydot non-NULL elements, each match the deriv filter.
746 - vars in solver's var list positions [n_y,n_y+n_ydot) all match deriv filter
747
748 - y_id contains n_ydot elements, with int values in range [0,n_y)
749 - var_sindex(ydot[y_id[i]]) - n_y == i for i in [0,n_ydot)
750
751 - all the vars [n_y+n_ydot,n_vars) fail the non-deriv filter and fail the deriv filter.
752 */
753 static int integrator_ida_check_diffindex(IntegratorSystem *sys){
754 int i, n_vars, n_ydot=0;
755 struct var_variable **list, *v;
756 char *varname;
757 const char *msg;
758
759 #ifdef ANALYSE_DEBUG
760 CONSOLE_DEBUG("Checking diffindex vector");
761 #endif
762
763 if(sys->y_id == NULL || sys->y == NULL || sys->ydot == NULL){
764 ERROR_REPORTER_HERE(ASC_PROG_ERR,"list(s) NULL");
765 return 1;
766 }
767
768 list = slv_get_solvers_var_list(sys->system);
769 n_vars = slv_get_num_solvers_vars(sys->system);
770
771 if(check_dups(sys, list, n_vars,FALSE)){
772 ERROR_REPORTER_HERE(ASC_PROG_ERR,"duplicates or NULLs in solver's var list");
773 return 1;
774 }
775
776 if(check_dups(sys, sys->ydot, sys->n_y,TRUE)){
777 ERROR_REPORTER_HERE(ASC_PROG_ERR,"duplicates in ydot vector");
778 return 1;
779 }
780
781 /* check n_y in range */
782 if(sys->n_y <= 0 || sys->n_y >= n_vars){
783 ERROR_REPORTER_HERE(ASC_PROG_ERR,"n_y = %d invalid (n_vars = %d)",sys->n_y, n_vars);
784 return 1;
785 }
786
787 /* check first n_y vars */
788 for(i=0; i < sys->n_y; ++i){
789 v = list[i];
790 if(!var_apply_filter(v, &integrator_ida_nonderiv)){
791 msg = "'%s' (in first n_y vars) fails non-deriv filter"; goto finish;
792 }else if(var_sindex(v) != i){
793 msg = "'%s' has wrong var_sindex"; goto finish;
794 }else if(v != sys->y[i]){
795 msg = "'%s' not matched in y vector"; goto finish;
796 }
797 /* meets filter, matches in y vector, has correct var_sindex. */
798 }
799
800 /* count non-NULLs in ydot */
801 for(i=0; i < sys->n_y; ++i){
802 v = sys->ydot[i];
803 if(v!=NULL){
804 if(var_sindex(v) < sys->n_y){
805 msg = "'%s' has var_sindex < n_y"; goto finish;
806 }else if(!var_apply_filter(v,&integrator_ida_deriv)){
807 msg = "'%s' (in next n_ydot vars) fails deriv filter"; goto finish;
808 }
809 /* lies beyond first n_y vars, match deriv filter */
810 n_ydot++;
811 }
812 }
813
814 /* check that vars [n_y, n_y+n_ydot) in solver's var list match the deriv filter */
815 for(i=sys->n_y; i< sys->n_y + n_ydot; ++i){
816 v = list[i];
817 if(!var_apply_filter(v,&integrator_ida_deriv)){
818 msg = "'%s', in range [n_y,n_y+n_ydot), fails deriv filter"; goto finish;
819 }
820 }
821
822 /* check values in y_id are ints int range [0,n_y), and that they point to correct vars */
823 for(i=0; i<n_ydot; ++i){
824 if(sys->y_id[i] < 0 || sys->y_id[i] >= sys->n_y){
825 ERROR_REPORTER_HERE(ASC_PROG_ERR,"value of y_id[%d] is out of range",i);
826 return 1;
827 }
828 v = sys->ydot[sys->y_id[i]];
829 if(var_sindex(v) - sys->n_y != i){
830 msg = "'%s' not index correctly from y_id"; goto finish;
831 }
832 }
833
834 /* check remaining vars fail both filters */
835 for(i=sys->n_y + n_ydot; i<n_vars; ++i){
836 v = list[i];
837 if(var_apply_filter(v,&integrator_ida_nonderiv)){
838 msg = "Var '%s' at end meets non-deriv filter, but shouldn't"; goto finish;
839 }
840 if(var_apply_filter(v,&integrator_ida_deriv)){
841 CONSOLE_DEBUG("position = %d",i);
842 msg = "Var '%s' at end meets deriv filter, but shouldn't"; goto finish;
843 }
844 }
845
846 return 0;
847 finish:
848 varname = var_make_name(sys->system,v);
849 ERROR_REPORTER_HERE(ASC_PROG_ERR,msg,varname);
850 ASC_FREE(varname);
851 return 1;
852 }
853
854 /**
855 Given a derivative variable, return the index of its corresponding differential
856 variable in the y vector (and equivalently the var_sindex of the diff var)
857 */
858 int integrator_ida_diffindex(const IntegratorSystem *sys, const struct var_variable *deriv){
859 asc_assert(var_sindex(deriv) >= sys->n_y);
860 asc_assert(var_sindex(deriv) < sys->n_y + sys->n_ydot);
861 return sys->y_id[var_sindex(deriv) - sys->n_y];
862 }
863
864 /**
865 This function will output the data structures provided to use BY THE
866 SYSTEM -- not the ones we're working with here IN THE SOLVER.
867 */
868 int integrator_ida_analyse_debug(const IntegratorSystem *sys,FILE *fp){
869 return system_diffvars_debug(sys->system,fp);
870 }
871

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