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

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