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

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