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

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