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

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