/[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 1241 - (show annotations) (download) (as text)
Fri Jan 26 11:59:45 2007 UTC (15 years, 5 months ago) by johnpye
File MIME type: text/x-csrc
File size: 16371 byte(s)
Fixed a bug with idaanalyse, hunting another one.
1 #include "idaanalyse.h"
2
3 #include <utilities/ascPanic.h>
4
5 #ifdef ASC_IDA_NEW_ANALYSE
6 # include <solver/diffvars.h>
7 # include <solver/slv_stdcalls.h>
8 # include <solver/block.h>
9 #include <solver/slvDOF.h>
10 #endif
11
12 #define ANALYSE_DEBUG
13
14 static int integrator_ida_check_partitioning(IntegratorSystem *sys);
15 static int integrator_ida_check_diffindex(IntegratorSystem *sys);
16 static int integrator_ida_rebuild_diffindex(IntegratorSystem *sys);
17
18 const var_filter_t integrator_ida_nonderiv = {
19 VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | VAR_FIXED | VAR_DERIV,
20 VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | 0 | 0
21 };
22
23 const var_filter_t integrator_ida_deriv = {
24 VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | VAR_FIXED | VAR_DERIV,
25 VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | 0 | VAR_DERIV
26 };
27
28 const rel_filter_t integrator_ida_rel = {
29 REL_INCLUDED | REL_EQUALITY | REL_ACTIVE,
30 REL_INCLUDED | REL_EQUALITY | REL_ACTIVE
31 };
32
33 /**
34 Check derivative chains: if a var or derivative is inelligible,
35 follow down the chain fixing & setting zero any derivatives.
36 */
37 static int integrator_ida_check_vars(IntegratorSystem *sys){
38 const SolverDiffVarCollection *diffvars;
39 char *varname, *derivname;
40 int n_y = 0;
41 int i, j;
42 struct var_variable *v;
43 SolverDiffVarSequence seq;
44
45 CONSOLE_DEBUG("BEFORE CHECKING VARS");
46 integrator_ida_analyse_debug(sys,stderr);
47
48 /* we shouldn't have allocated these yet: just be sure */
49 asc_assert(sys->y==NULL);
50 asc_assert(sys->ydot==NULL);
51 asc_assert(sys->n_y==0);
52
53 /* get the the dervative chains from the system */
54 diffvars = slv_get_diffvars(sys->system);
55 if(diffvars==NULL){
56 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");
57 return 1;
58 }
59
60 /* add the variables from the derivative chains */
61 for(i=0; i<diffvars->nseqs; ++i){
62
63
64 CONSOLE_DEBUG("BEFORE CHAIN %d",i);
65 integrator_ida_analyse_debug(sys,stderr);
66
67 seq = diffvars->seqs[i];
68 asc_assert(seq.n >= 1);
69 v = seq.vars[0];
70 if(!var_apply_filter(v,&integrator_ida_nonderiv)){
71 varname = var_make_name(sys->system,v);
72 CONSOLE_DEBUG("'%s' fails non-deriv filter",varname);
73 ASC_FREE(varname);
74 for(j=1;j<seq.n;++j){
75 v = seq.vars[j];
76 var_set_active(v,FALSE);
77 var_set_value(v,0);
78 varname = var_make_name(sys->system,v);
79 CONSOLE_DEBUG("Derivative '%s' SET ZERO AND INACTIVE",varname);
80 ASC_FREE(varname);
81 }
82 continue;
83 }
84
85 varname = var_make_name(sys->system,v);
86 CONSOLE_DEBUG("We will use var '%s'",varname);
87 ASC_FREE(varname);
88 if(seq.n > 2){
89 varname = var_make_name(sys->system,v);
90 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Too-long derivative chain with var '%s'");
91 ASC_FREE(varname);
92 return 2;
93 }else if(seq.n==2){
94 /* diff var */
95 if(var_apply_filter(seq.vars[1],&integrator_ida_deriv)){
96 /* add the diff & deriv vars to the lists */
97 n_y++;
98 varname = var_make_name(sys->system,v);
99 derivname = var_make_name(sys->system,seq.vars[1]);
100 CONSOLE_DEBUG("Added '%s' and '%s'",varname,derivname);
101 ASC_FREE(varname);
102 ASC_FREE(derivname);
103 continue;
104 }
105 varname = var_make_name(sys->system,v);
106 derivname = var_make_name(sys->system,seq.vars[1]);
107 CONSOLE_DEBUG("Derivative '%s' of '%s' fails filter; convert '%s' to algebraic",derivname,varname,varname);
108 ASC_FREE(varname);
109 ASC_FREE(derivname);
110 /* fall through */
111 }
112
113 varname = var_make_name(sys->system,v);
114 CONSOLE_DEBUG("Adding '%s' to algebraic",varname);
115 ASC_FREE(varname);
116 n_y++;
117 }
118
119 /* we assert that all vars in y meet the integrator_ida_nonderiv filter */
120 /* we assert that all vars in ydot meet the integrator_ida_deriv filter */
121
122 CONSOLE_DEBUG("Found %d good non-derivative vars", n_y);
123 sys->n_y = n_y;
124
125 CONSOLE_DEBUG("AFTER CHECKING VARS");
126 integrator_ida_analyse_debug(sys,stderr);
127
128 return 0;
129 }
130
131
132 /**
133 Sort the lists. First we will put the non-derivative vars, then we will
134 put the derivative vars. Then we will put all the others.
135 */
136 static int integrator_ida_sort_rels_and_vars(IntegratorSystem *sys){
137 int ny1, nydot, nr;
138
139
140 CONSOLE_DEBUG("BEFORE SORTING RELS AND VARS");
141 integrator_ida_analyse_debug(sys,stderr);
142
143 /* we should not have allocated y or ydot yet */
144 asc_assert(sys->y==NULL && sys->ydot==NULL);
145
146 /* but we should have found some variables (and know how many) */
147 asc_assert(sys->n_y);
148
149 if(system_cut_vars(sys->system, 0, &integrator_ida_nonderiv, &ny1)){
150 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting non-derivs");
151 return 1;
152 }
153
154 CONSOLE_DEBUG("cut_vars: ny1 = %d, sys->n_y = %d",ny1,sys->n_y);
155 asc_assert(ny1 == sys->n_y);
156
157 if(system_cut_vars(sys->system, ny1, &integrator_ida_deriv, &nydot)){
158 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting derivs");
159 return 1;
160 }
161
162 if(system_cut_rels(sys->system, 0, &integrator_ida_rel, &nr)){
163 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting derivs");
164 return 1;
165 }
166
167 if(ny1 != nr){
168 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem is not square (ny = %d, nr = %d)",ny1,nr);
169 return 2;
170 }
171
172 return 0;
173 }
174
175 /**
176 Build up the lists 'ydot' and 'y_id.
177 'ydot' allows us to find the derivative of a variable given its sindex.
178 'y_id' allows us to locate the variable of a derivative given its sindex.
179
180 Hence
181 y_id[var_sindex(derivvar)-sys->n_y] = var_sindex(diffvar)
182 ydot[var_sindex(diffvar)] = derivvar (or NULL if diffvar is algebraic)
183
184 Note that there is 'shifting' required when looking up y_id.
185
186 Note that we want to get rid of 'y' but for the moment we are also assigning
187 that. We want to be rid of it because if the vars are ordered correctly
188 it shouldn't be needed.
189 */
190 static int integrator_ida_create_lists(IntegratorSystem *sys){
191 const SolverDiffVarCollection *diffvars;
192 int i, j;
193 struct var_variable *v;
194 char *varname;
195
196 SolverDiffVarSequence seq;
197
198 asc_assert(sys->y==NULL);
199 asc_assert(sys->ydot==NULL);
200 asc_assert(sys->y_id== NULL);
201
202 CONSOLE_DEBUG("BEFORE MAKING LISTS");
203 integrator_ida_analyse_debug(sys,stderr);
204
205 /* get the the dervative chains from the system */
206 diffvars = slv_get_diffvars(sys->system);
207 if(diffvars==NULL){
208 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");
209 return 1;
210 }
211
212 /* allocate space (more than enough) */
213 sys->y = ASC_NEW_ARRAY(struct var_variable *,sys->n_y);
214 sys->y_id = ASC_NEW_ARRAY_CLEAR(int,sys->n_y);
215 sys->ydot = ASC_NEW_ARRAY_CLEAR(struct var_variable *,sys->n_y);
216 j = 0;
217
218 CONSOLE_DEBUG("Passing through chains...");
219
220 /* add the variables from the derivative chains */
221 for(i=0; i<diffvars->nseqs; ++i){
222
223 CONSOLE_DEBUG("i = %d",i);
224
225 seq = diffvars->seqs[i];
226 asc_assert(seq.n >= 1);
227 v = seq.vars[0];
228
229 varname = var_make_name(sys->system, v);
230 CONSOLE_DEBUG("alg '%s'",varname);
231 ASC_FREE(varname);
232
233
234 if(!var_apply_filter(v,&integrator_ida_nonderiv)){
235 continue;
236 }
237
238 varname = var_make_name(sys->system, v);
239 CONSOLE_DEBUG("alg '%s' is GOOD",varname);
240 ASC_FREE(varname);
241
242 if(seq.n > 1 && var_apply_filter(seq.vars[1],&integrator_ida_deriv)){
243 asc_assert(var_sindex(seq.vars[1]) >= sys->n_y);
244 asc_assert(var_sindex(seq.vars[1])-sys->n_y < sys->n_y);
245
246 varname = var_make_name(sys->system, seq.vars[1]);
247 CONSOLE_DEBUG("diff '%s' IS GOOD",varname);
248 ASC_FREE(varname);
249
250 sys->y_id[var_sindex(seq.vars[1]) - sys->n_y] = j;
251 sys->ydot[j] = seq.vars[1];
252 }else{
253 asc_assert(sys->ydot[j]==NULL);
254 }
255
256 sys->y[j] = v;
257 j++;
258 }
259
260 CONSOLE_DEBUG("AFTER MAKING LISTS");
261 integrator_ida_analyse_debug(sys,stderr);
262
263 return 0;
264 }
265
266 /**
267 Perform additional problem analysis to prepare problem for integration with
268 IDA.
269
270 We assume that the analyse_generate_diffvars routine has been called, so
271 that we just need to call slv_get_diffvars to access the derivative
272 chains.
273
274 We can also assume that the independent variable has been found.
275
276 See mailing list, ~Jan 2007.
277
278 Note, the stuff for identifying the static and output sub-problems should
279 be part of integrator.c, not this file. We will assume this is handled
280
281 @return 0 on success
282 @see integrator_analyse
283 */
284 int integrator_ida_analyse(IntegratorSystem *sys){
285 int res;
286 const SolverDiffVarCollection *diffvars;
287 char *varname;
288 int i;
289
290 asc_assert(sys->engine==INTEG_IDA);
291
292 #ifdef ANALYSE_DEBUG
293 CONSOLE_DEBUG("Starting IDA analysis");
294 #endif
295
296 res = integrator_ida_check_vars(sys);
297 if(res){
298 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem with vars");
299 return 1;
300 }
301
302 #ifdef ANALYSE_DEBUG
303 CONSOLE_DEBUG("Sorting rels and vars");
304 #endif
305
306 res = integrator_ida_sort_rels_and_vars(sys);
307 if(res){
308 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem sorting rels and vars");
309 return 1;
310 }
311
312 #ifdef ANALYSE_DEBUG
313 CONSOLE_DEBUG("Creating lists");
314 #endif
315
316 res = integrator_ida_create_lists(sys);
317 if(res){
318 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem creating lists");
319 return 1;
320 }
321
322 #ifdef ANALYSE_DEBUG
323 CONSOLE_DEBUG("Checking");
324 #endif
325
326 asc_assert(sys->y);
327 asc_assert(sys->ydot);
328 asc_assert(sys->y_id);
329
330 integrator_ida_analyse_debug(sys,stderr);
331
332 if(integrator_ida_check_diffindex(sys)){
333 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error with diffindex");
334 return 360;
335 }
336
337 /* the following causes TestIDA.testincidence3 to fail:
338 if(sys->n_y != slv_get_num_solvers_rels(sys->system)){
339 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Problem is not square");
340 return 2;
341 }*/
342
343 #if 0
344 /* check structural singularity for the two IDACalcIC scenarios */
345
346 /* ...(1) FIX the derivatives */
347 CONSOLE_DEBUG("Checking system with derivatives fixed...");
348 for(i=0;i<n_y;++i){
349 if(sys->ydot[i])var_set_fixed(sys->ydot[i],1);
350 }
351
352 slv_block_partition(sys->system);
353 res = integrator_ida_block_check(sys);
354 if(res)return 100 + res;
355
356 /* ...(2) FREE the derivatives, FIX the diffvars */
357 CONSOLE_DEBUG("Checking system with differential variables fixed...");
358 for(i=0;i<n_y;++i){
359 if(sys->ydot[i]){
360 var_set_fixed(sys->ydot[i],0);
361 var_set_fixed(sys->y[i],1);
362 }
363 }
364 slv_block_partition(sys->system);
365 res = integrator_ida_block_check(sys);
366 if(res)return 200 + res;
367
368 /* ...(3) restore as it was, FREE the diffvars */
369 for(i=0;i<n_y;++i){
370 if(sys->ydot[i]){
371 var_set_fixed(sys->y[i],0);
372 }
373 }
374 #endif
375
376 /* get the the dervative chains from the system */
377 diffvars = slv_get_diffvars(sys->system);
378
379 /* check the indep var is same as was located elsewhere */
380 if(diffvars->nindep>1){
381 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Multiple variables specified as independent (ode_type=-1)");
382 return 3;
383 }else if(diffvars->nindep<1){
384 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Independent var not set (ode_type=-1)");
385 return 4;
386 }else if(diffvars->indep[0]!=sys->x){
387 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Indep var doesn't match");
388 return 5;
389 }
390
391 /* get the observations */
392 /** @TODO this should use a 'problem provider' API instead of static data
393 from the system object */
394 sys->n_obs = diffvars->nobs;
395 sys->obs = ASC_NEW_ARRAY(struct var_variable *,sys->n_obs);
396 for(i=0;i<sys->n_obs;++i){
397 sys->obs[i] = diffvars->obs[i];
398 varname = var_make_name(sys->system,sys->obs[i]);
399 CONSOLE_DEBUG("'%s' is observation",varname);
400 ASC_FREE(varname);
401 }
402
403 /* - 'y' list as [ya|yd] */
404 /* - sparsity pattern for dF/dy and dF/dy' */
405 /* - sparsity pattern for union of above */
406 /* - block decomposition based on above */
407 /* - block decomposition results in reordering of y and y' */
408 /* - boundaries (optional) */
409 /* ERROR_REPORTER_HERE(ASC_PROG_ERR,"Implementation incomplete");
410 return -1; */
411
412 return 0;
413 }
414
415 /*------------------------------------------------------------------------------
416 ANALYSIS ROUTINE (new implementation)
417 */
418
419 static int integrator_ida_check_partitioning(IntegratorSystem *sys){
420 long i, nv;
421 int err=0;
422 char *varname;
423 struct var_variable **vlist, *v;
424 nv = slv_get_num_solvers_vars(sys->system);
425 vlist = slv_get_solvers_var_list(sys->system);
426 var_filter_t vf = {VAR_SVAR|VAR_ACTIVE|VAR_INCIDENT|VAR_FIXED
427 ,VAR_SVAR|VAR_ACTIVE|VAR_INCIDENT|0 };
428 for(i=0;i<nv;++i){
429 v=vlist[i];
430 if(!var_apply_filter(v,&vf))continue;
431 varname = var_make_name(sys->system,v);
432 if(!var_deriv(v)){
433 fprintf(stderr,"vlist[%ld] = '%s' (nonderiv)\n",i,varname);
434 if(i>=sys->n_y){
435 ERROR_REPORTER_HERE(ASC_PROG_ERR,"non-deriv var '%s' is not at the start",varname);
436 err++;
437 }
438 }else{
439 fprintf(stderr,"vlist[%ld] = '%s' (derivative)\n",i,varname);
440 if(i<sys->n_y){
441 ERROR_REPORTER_HERE(ASC_PROG_ERR,"deriv var '%s' is not at the end (n_y = %d, i = %d)"
442 ,varname, sys->n_y, i
443 );
444 err++;
445 }
446 }
447 ASC_FREE(varname);
448 }
449 if(err){
450 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error in IDA partitioning");
451 return 1;
452 }
453 return 0;
454 }
455
456 int integrator_ida_block_check(IntegratorSystem *sys){
457 int res;
458 int dof;
459 #ifdef ANALYSE_DEBUG
460 long *vlist, *vp, i, nv, nv_ok;
461 char *varname;
462 struct var_variable **solversvars;
463
464 var_filter_t vfilt = {
465 VAR_ACTIVE | VAR_INCIDENT | VAR_FIXED,
466 VAR_ACTIVE | VAR_INCIDENT | 0
467 };
468
469 nv = slv_get_num_solvers_vars(sys->system);
470 solversvars = slv_get_solvers_var_list(sys->system);
471 CONSOLE_DEBUG("-------------- nv = %ld -------------",nv);
472 for(nv_ok=0, i=0; i < nv; ++i){
473 if(var_apply_filter(solversvars[i],&vfilt)){
474 varname = var_make_name(sys->system,solversvars[i]);
475 fprintf(stderr,"%s\n",varname);
476 ASC_FREE(varname);
477 nv_ok++;
478 }
479 }
480 CONSOLE_DEBUG("----------- got %ld ok -------------",nv_ok);
481 #endif
482
483 if(!slvDOF_status(sys->system, &res, &dof)){
484 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to determine DOF status");
485 return -1;
486 }
487 switch(res){
488 case 1: CONSOLE_DEBUG("System is underspecified (%d degrees of freedom)",dof);break;
489 case 2: CONSOLE_DEBUG("System is square"); return 0; /* all OK */
490 case 3: CONSOLE_DEBUG("System is structurally singular"); break;
491 case 4: CONSOLE_DEBUG("System is overspecified"); break;
492 default:
493 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unrecognised slfDOF_status");
494 return -2;
495 }
496
497 #ifdef ANALYSE_DEBUG
498 /* if it was underspecified, what vars could be fixed? */
499 if(res==1){
500 CONSOLE_DEBUG("Need to FIX %d of the following vars:",dof);
501 solversvars = slv_get_solvers_var_list(sys->system);
502 if(!slvDOF_eligible(sys->system, &vlist)){
503 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to det slvDOF_eligble list");
504 return -3;
505 }
506 for(vp=vlist;*vp!=-1;++vp){
507 varname = var_make_name(sys->system, solversvars[*vp]);
508 CONSOLE_DEBUG("Fixable var: %s",varname);
509 ASC_FREE(varname);
510 }
511 CONSOLE_DEBUG("(Found %d fixable vars)",(int)(vp-vlist));
512 return 1;
513 }
514 #endif
515
516 return res;
517 }
518
519 static int integrator_ida_check_diffindex(IntegratorSystem *sys){
520 int i, nv, err = 0;
521 struct var_variable **vlist;
522 int diffindex;
523
524 CONSOLE_DEBUG("Checking diffindex vector");
525
526 if(sys->y_id == NULL){
527 CONSOLE_DEBUG("y_id not allocated");
528 return 1;
529 }
530
531 vlist = slv_get_solvers_var_list(sys->system);
532 nv = slv_get_num_solvers_vars(sys->system);
533
534 for(i=0; i<nv; ++i){
535 if(var_deriv(vlist[i])){
536 if(i < sys->n_y){
537 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Deriv in wrong position");
538 err++;
539 }else{
540 if(var_sindex(vlist[i])!=i){
541 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Deriv var with incorrect sindex");
542 err++;
543 }
544 diffindex = sys->y_id[i - sys->n_y];
545 if(diffindex >= sys->n_y){
546 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Diff y_id[%d] value too big",i-sys->n_y);
547 err++;
548 }
549 if(var_sindex(vlist[diffindex])!=diffindex){
550 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Diff var with incorrect sindex");
551 err++;
552 }
553 }
554 }else{
555 if(i!=var_sindex(vlist[i])){
556 ERROR_REPORTER_HERE(ASC_PROG_ERR,"var_sindex incorrect");
557 err++;
558 }
559 }
560 }
561 if(err){
562 CONSOLE_DEBUG("Errors found in diffindex");
563 }
564 return err;
565 }
566
567 /**
568 @TODO provide a macro implementation of this
569 */
570 int integrator_ida_diffindex(const IntegratorSystem *sys, const struct var_variable *deriv){
571 int diffindex;
572 #ifdef DIFFINDEX_DEBUG
573 asc_assert( var_deriv (deriv));
574 asc_assert( var_active (deriv));
575 asc_assert( var_incident (deriv));
576 asc_assert( var_svar (deriv));
577
578 asc_assert(!var_fixed (deriv));
579
580 asc_assert(var_sindex(deriv) >= sys->n_y);
581 asc_assert(diffindex == var_sindex(sys->y[diffindex]));
582 #endif
583 asc_assert(var_sindex(deriv) >= sys->n_y);
584 diffindex = sys->y_id[var_sindex(deriv) - sys->n_y];
585
586 return diffindex;
587 }
588
589
590 /**
591 This function will output the data structures provided to use BY THE
592 SYSTEM -- not the ones we're working with here IN THE SOLVER.
593 */
594 int integrator_ida_analyse_debug(const IntegratorSystem *sys,FILE *fp){
595 return analyse_diffvars_debug(sys->system,fp);
596 }
597

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