/[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 1268 - (show annotations) (download) (as text)
Sat Feb 3 06:18:11 2007 UTC (17 years, 5 months ago) by johnpye
File MIME type: text/x-csrc
File size: 20023 byte(s)
Stable with 4 nodes now, but only with a small peturbation.
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 for(i=0;i<sys->n_y;++i){
239 asc_assert(sys->ydot[i] == 0);
240 }
241
242 CONSOLE_DEBUG("Passing through chains...");
243
244 /* create the lists y and ydot, ignoring 'bad' vars */
245 for(i=0; i<diffvars->nseqs; ++i){
246 /* CONSOLE_DEBUG("i = %d",i); */
247
248 seq = diffvars->seqs[i];
249 asc_assert(seq.n >= 1);
250 v = seq.vars[0];
251 j = var_sindex(v);
252
253 if(!var_apply_filter(v,&integrator_ida_nonderiv)){
254 continue;
255 }
256
257 sys->y[j] = v;
258 /* VARMSG("'%s' is good non-deriv"); */
259
260 if(seq.n > 1 && var_apply_filter(seq.vars[1],&integrator_ida_deriv)){
261 v = seq.vars[1];
262 asc_assert(var_sindex(v) >= sys->n_y);
263 asc_assert(var_sindex(v)-sys->n_y < sys->n_y);
264
265 sys->ydot[j] = v;
266 sys->n_ydot++;
267 /* VARMSG("'%s' is good deriv"); */
268 }else{
269 asc_assert(sys->ydot[j]==NULL);
270 }
271 }
272
273 CONSOLE_DEBUG("Found %d good non-derivs",j);
274
275 /* create the list y_id by looking at non-NULLs from ydot */
276 sys->y_id = ASC_NEW_ARRAY(int,sys->n_ydot);
277 for(i=0,j=0; i < sys->n_y; ++i){
278 if(sys->ydot[i]==NULL)continue;
279 v = sys->ydot[i]; VARMSG("deriv '%s'...");
280 v = sys->y[i]; VARMSG("diff '%s'...");
281 sys->y_id[var_sindex(sys->ydot[i]) - sys->n_y] = i;
282 j++;
283 }
284
285 asc_assert(j==sys->n_ydot);
286
287 return 0;
288 }
289
290 /**
291 Perform additional problem analysis to prepare problem for integration with
292 IDA.
293
294 We assume that the analyse_generate_diffvars routine has been called, so
295 that we just need to call slv_get_diffvars to access the derivative
296 chains.
297
298 We can also assume that the independent variable has been found.
299
300 See mailing list, ~Jan 2007.
301
302 Note, the stuff for identifying the static and output sub-problems should
303 be part of integrator.c, not this file. We will assume this is handled
304
305 @return 0 on success
306 @see integrator_analyse
307 */
308 int integrator_ida_analyse(IntegratorSystem *sys){
309 int res;
310 const SolverDiffVarCollection *diffvars;
311 char *varname;
312 int i;
313
314 asc_assert(sys->engine==INTEG_IDA);
315
316 #ifdef ANALYSE_DEBUG
317 CONSOLE_DEBUG("Starting IDA analysis");
318 #endif
319
320 res = integrator_ida_check_vars(sys);
321 if(res){
322 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem with vars");
323 return 1;
324 }
325
326 #ifdef ANALYSE_DEBUG
327 CONSOLE_DEBUG("Sorting rels and vars");
328 #endif
329
330 res = integrator_ida_sort_rels_and_vars(sys);
331 if(res){
332 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem sorting rels and vars");
333 return 1;
334 }
335
336 #ifdef ANALYSE_DEBUG
337 CONSOLE_DEBUG("Creating lists");
338 #endif
339
340 CONSOLE_DEBUG("BEFORE MAKING LISTS");
341 integrator_ida_debug(sys,stderr);
342
343 res = integrator_ida_create_lists(sys);
344 if(res){
345 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem creating lists");
346 return 1;
347 }
348
349 #ifdef ANALYSE_DEBUG
350 CONSOLE_DEBUG("Checking lists");
351 #endif
352
353 asc_assert(sys->y);
354 asc_assert(sys->ydot);
355 asc_assert(sys->y_id);
356
357 integrator_ida_debug(sys,stderr);
358
359 if(integrator_ida_check_diffindex(sys)){
360 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error with diffindex");
361 return 360;
362 }
363
364 /* the following causes TestIDA.testincidence3 to fail:
365 if(sys->n_y != slv_get_num_solvers_rels(sys->system)){
366 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Problem is not square");
367 return 2;
368 }*/
369
370 #if 0
371 /* check structural singularity for the two IDACalcIC scenarios */
372
373 /* ...(1) FIX the derivatives */
374 CONSOLE_DEBUG("Checking system with derivatives fixed...");
375 for(i=0;i<n_y;++i){
376 if(sys->ydot[i])var_set_fixed(sys->ydot[i],1);
377 }
378
379 slv_block_partition(sys->system);
380 res = integrator_ida_block_check(sys);
381 if(res)return 100 + res;
382
383 /* ...(2) FREE the derivatives, FIX the diffvars */
384 CONSOLE_DEBUG("Checking system with differential variables fixed...");
385 for(i=0;i<n_y;++i){
386 if(sys->ydot[i]){
387 var_set_fixed(sys->ydot[i],0);
388 var_set_fixed(sys->y[i],1);
389 }
390 }
391 slv_block_partition(sys->system);
392 res = integrator_ida_block_check(sys);
393 if(res)return 200 + res;
394
395 /* ...(3) restore as it was, FREE the diffvars */
396 for(i=0;i<n_y;++i){
397 if(sys->ydot[i]){
398 var_set_fixed(sys->y[i],0);
399 }
400 }
401 #endif
402
403 /* get the the dervative chains from the system */
404 diffvars = system_get_diffvars(sys->system);
405
406 /* check the indep var is same as was located elsewhere */
407 if(diffvars->nindep>1){
408 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Multiple variables specified as independent (ode_type=-1)");
409 return 3;
410 }else if(diffvars->nindep<1){
411 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Independent var not set (ode_type=-1)");
412 return 4;
413 }else if(diffvars->indep[0]!=sys->x){
414 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Indep var doesn't match");
415 return 5;
416 }
417
418 #ifdef ANALYSE_DEBUG
419 CONSOLE_DEBUG("Collecting observed variables");
420 #endif
421
422 /* get the observations */
423 /** @TODO this should use a 'problem provider' API instead of static data
424 from the system object */
425 sys->n_obs = diffvars->nobs;
426 sys->obs = ASC_NEW_ARRAY(struct var_variable *,sys->n_obs);
427 for(i=0;i<sys->n_obs;++i){
428 /* we get them all, regardless of flags etc */
429 sys->obs[i] = diffvars->obs[i];
430 varname = var_make_name(sys->system,sys->obs[i]);
431 CONSOLE_DEBUG("'%s' is observation",varname);
432 ASC_FREE(varname);
433 }
434
435 return 0;
436 }
437
438 /*------------------------------------------------------------------------------
439 ANALYSIS ROUTINE (new implementation)
440 */
441
442 static int integrator_ida_check_partitioning(IntegratorSystem *sys){
443 long i, nv;
444 int err=0;
445 char *varname;
446 struct var_variable **vlist, *v;
447 nv = slv_get_num_solvers_vars(sys->system);
448 vlist = slv_get_solvers_var_list(sys->system);
449 var_filter_t vf = {VAR_SVAR|VAR_ACTIVE|VAR_INCIDENT|VAR_FIXED
450 ,VAR_SVAR|VAR_ACTIVE|VAR_INCIDENT|0 };
451 for(i=0;i<nv;++i){
452 v=vlist[i];
453 if(!var_apply_filter(v,&vf))continue;
454 varname = var_make_name(sys->system,v);
455 if(!var_deriv(v)){
456 fprintf(stderr,"vlist[%ld] = '%s' (nonderiv)\n",i,varname);
457 if(i>=sys->n_y){
458 ERROR_REPORTER_HERE(ASC_PROG_ERR,"non-deriv var '%s' is not at the start",varname);
459 err++;
460 }
461 }else{
462 fprintf(stderr,"vlist[%ld] = '%s' (derivative)\n",i,varname);
463 if(i<sys->n_y){
464 ERROR_REPORTER_HERE(ASC_PROG_ERR,"deriv var '%s' is not at the end (n_y = %d, i = %d)"
465 ,varname, sys->n_y, i
466 );
467 err++;
468 }
469 }
470 ASC_FREE(varname);
471 }
472 if(err){
473 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error in IDA partitioning");
474 return 1;
475 }
476 return 0;
477 }
478
479 int integrator_ida_block_check(IntegratorSystem *sys){
480 int res;
481 int dof;
482 #ifdef ANALYSE_DEBUG
483 int *vlist, *vp, i, nv, nv_ok;
484 char *varname;
485 struct var_variable **solversvars;
486
487 var_filter_t vfilt = {
488 VAR_ACTIVE | VAR_INCIDENT | VAR_FIXED,
489 VAR_ACTIVE | VAR_INCIDENT | 0
490 };
491
492 nv = slv_get_num_solvers_vars(sys->system);
493 solversvars = slv_get_solvers_var_list(sys->system);
494 CONSOLE_DEBUG("-------------- nv = %d -------------",nv);
495 for(nv_ok=0, i=0; i < nv; ++i){
496 if(var_apply_filter(solversvars[i],&vfilt)){
497 varname = var_make_name(sys->system,solversvars[i]);
498 fprintf(stderr,"%s\n",varname);
499 ASC_FREE(varname);
500 nv_ok++;
501 }
502 }
503 CONSOLE_DEBUG("----------- got %d ok -------------",nv_ok);
504 #endif
505
506 if(!slvDOF_status(sys->system, &res, &dof)){
507 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to determine DOF status");
508 return -1;
509 }
510 switch(res){
511 case 1: CONSOLE_DEBUG("System is underspecified (%d degrees of freedom)",dof);break;
512 case 2: CONSOLE_DEBUG("System is square"); return 0; /* all OK */
513 case 3: CONSOLE_DEBUG("System is structurally singular"); break;
514 case 4: CONSOLE_DEBUG("System is overspecified"); break;
515 default:
516 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unrecognised slfDOF_status");
517 return -2;
518 }
519
520 #ifdef ANALYSE_DEBUG
521 /* if it was underspecified, what vars could be fixed? */
522 if(res==1){
523 CONSOLE_DEBUG("Need to FIX %d of the following vars:",dof);
524 solversvars = slv_get_solvers_var_list(sys->system);
525 if(!slvDOF_eligible(sys->system, &vlist)){
526 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to det slvDOF_eligble list");
527 return -3;
528 }
529 for(vp=vlist;*vp!=-1;++vp){
530 varname = var_make_name(sys->system, solversvars[*vp]);
531 CONSOLE_DEBUG("Fixable var: %s",varname);
532 ASC_FREE(varname);
533 }
534 CONSOLE_DEBUG("(Found %d fixable vars)",(int)(vp-vlist));
535 return 1;
536 }
537 #endif
538
539 return res;
540 }
541
542 /* return 0 on succes */
543 static int check_dups(IntegratorSystem *sys, struct var_variable **list,int n,int allownull){
544 int i,j;
545 struct var_variable *v;
546 char *varname;
547 for(i=0; i< n; ++i){
548 v=list[i];
549 if(v==NULL){
550 if(allownull)continue;
551 else return 2;
552 }
553 for(j=0; j<i-1;++j){
554 if(list[j]==NULL)continue;
555 if(v==list[j]){
556 CONSOLE_DEBUG("duplicate found (%d, v = %p)",j,v);
557 varname = var_make_name(sys->system,v);
558 if(varname){
559 CONSOLE_DEBUG("Duplicate of '%s' found",varname);
560 ASC_FREE(varname);
561 }else{
562 CONSOLE_DEBUG("Duplicate found (couldn't retrieve name)");
563 }
564 return 1;
565 }
566 }
567 }
568 return 0;
569 }
570
571 /*
572 We are going to check that
573
574 - n_y in range (0,n_vars)
575 - no duplicates anywhere in the varlist
576 - no duplicates in non-NULL elements of ydot
577
578 - first n_y vars in solver's var list match those in vector y and match the non-deriv filter.
579 - var_sindex for first n_y vars match their position the the solver's var list
580
581 - ydot contains n_ydot non-NULL elements, each match the deriv filter.
582 - vars in solver's var list positions [n_y,n_y+n_ydot) all match deriv filter
583
584 - y_id contains n_ydot elements, with int values in range [0,n_y)
585 - var_sindex(ydot[y_id[i]]) - n_y == i for i in [0,n_ydot)
586
587 - all the vars [n_y+n_ydot,n_vars) fail the non-deriv filter and fail the deriv filter.
588 */
589 static int integrator_ida_check_diffindex(IntegratorSystem *sys){
590 int i, n_vars, n_ydot=0;
591 struct var_variable **list, *v;
592 char *varname;
593 int diffindex;
594 const char *msg;
595
596 CONSOLE_DEBUG("Checking diffindex vector");
597
598 if(sys->y_id == NULL || sys->y == NULL || sys->ydot == NULL){
599 ERROR_REPORTER_HERE(ASC_PROG_ERR,"list(s) NULL");
600 return 1;
601 }
602
603 list = slv_get_solvers_var_list(sys->system);
604 n_vars = slv_get_num_solvers_vars(sys->system);
605
606 if(check_dups(sys, list, n_vars,FALSE)){
607 ERROR_REPORTER_HERE(ASC_PROG_ERR,"duplicates or NULLs in solver's var list");
608 return 1;
609 }
610
611 if(check_dups(sys, sys->ydot, n_vars,TRUE)){
612 ERROR_REPORTER_HERE(ASC_PROG_ERR,"duplicates in ydot vector");
613 return 1;
614 }
615
616 /* check n_y in range */
617 if(sys->n_y <= 0 || sys->n_y >= n_vars){
618 ERROR_REPORTER_HERE(ASC_PROG_ERR,"n_y = %d invalid (n_vars = %d)",sys->n_y, n_vars);
619 return 1;
620 }
621
622 /* check first n_y vars */
623 for(i=0; i < sys->n_y; ++i){
624 v = list[i];
625 if(!var_apply_filter(v, &integrator_ida_nonderiv)){
626 msg = "'%s' (in first n_y vars) fails non-deriv filter"; goto finish;
627 }else if(var_sindex(v) != i){
628 msg = "'%s' has wrong var_sindex"; goto finish;
629 }else if(v != sys->y[i]){
630 msg = "'%s' not matched in y vector"; goto finish;
631 }
632 /* meets filter, matches in y vector, has correct var_sindex. */
633 }
634
635 /* count non-NULLs in ydot */
636 for(i=0; i < sys->n_y; ++i){
637 v = sys->ydot[i];
638 if(v!=NULL){
639 if(var_sindex(v) < sys->n_y){
640 msg = "'%s' has var_sindex < n_y"; goto finish;
641 }else if(!var_apply_filter(v,&integrator_ida_deriv)){
642 msg = "'%s' (in next n_ydot vars) fails deriv filter"; goto finish;
643 }
644 /* lies beyond first n_y vars, match deriv filter */
645 n_ydot++;
646 }
647 }
648
649 /* check that vars [n_y, n_y+n_ydot) in solver's var list match the deriv filter */
650 for(i=sys->n_y; i< sys->n_y + n_ydot; ++i){
651 v = list[i];
652 if(!var_apply_filter(v,&integrator_ida_deriv)){
653 msg = "'%s', in range [n_y,n_y+n_ydot), fails deriv filter"; goto finish;
654 }
655 }
656
657 /* check values in y_id are ints int range [0,n_y), and that they point to correct vars */
658 for(i=0; i<n_ydot; ++i){
659 if(sys->y_id[i] < 0 || sys->y_id[i] >= sys->n_y){
660 ERROR_REPORTER_HERE(ASC_PROG_ERR,"value of y_id[%d] is out of range",i);
661 return 1;
662 }
663 v = sys->ydot[sys->y_id[i]];
664 if(var_sindex(v) - sys->n_y != i){
665 msg = "'%s' not index correctly from y_id"; goto finish;
666 }
667 }
668
669 /* check remaining vars fail both filters */
670 for(i=sys->n_y + n_ydot; i<n_vars; ++i){
671 v = list[i];
672 if(var_apply_filter(v,&integrator_ida_nonderiv)){
673 msg = "Var '%s' at end meets non-deriv filter, but shouldn't"; goto finish;
674 }
675 if(var_apply_filter(v,&integrator_ida_deriv)){
676 CONSOLE_DEBUG("position = %d",i);
677 msg = "Var '%s' at end meets deriv filter, but shouldn't"; goto finish;
678 }
679 }
680
681 return 0;
682 finish:
683 varname = var_make_name(sys->system,v);
684 ERROR_REPORTER_HERE(ASC_PROG_ERR,msg,varname);
685 ASC_FREE(varname);
686 return 1;
687 }
688
689 /**
690 Given a derivative variable, return the index of its corresponding differential
691 variable in the y vector (and equivalently the var_sindex of the diff var)
692 */
693 int integrator_ida_diffindex(const IntegratorSystem *sys, const struct var_variable *deriv){
694 asc_assert(var_sindex(deriv) >= sys->n_y);
695 asc_assert(var_sindex(deriv) < sys->n_y + sys->n_ydot);
696 return sys->y_id[var_sindex(deriv) - sys->n_y];
697 }
698
699
700 /**
701 This function will output the data structures provided to use BY THE
702 SYSTEM -- not the ones we're working with here IN THE SOLVER.
703 */
704 int integrator_ida_analyse_debug(const IntegratorSystem *sys,FILE *fp){
705 return system_diffvars_debug(sys->system,fp);
706 }
707

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