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

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