/[ascend]/trunk/base/generic/integrator/idaanalyse.c
ViewVC logotype

Annotation of /trunk/base/generic/integrator/idaanalyse.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1337 - (hide annotations) (download) (as text)
Sat Mar 10 12:33:15 2007 UTC (13 years, 6 months ago) by jpye
File MIME type: text/x-csrc
File size: 26045 byte(s)
Added support for detection of boundary crossing in the IDA integrator.
Cleared up some code comments and formatting.
1 johnpye 1239 #include "idaanalyse.h"
2    
3     #include <utilities/ascPanic.h>
4 jpye 1337 #include <utilities/error.h>
5 johnpye 1239
6     #ifdef ASC_IDA_NEW_ANALYSE
7 johnpye 1316 # include <linear/linsolqr.h>
8    
9     # include <system/diffvars.h>
10     # include <system/slv_stdcalls.h>
11     # include <system/block.h>
12     # include <system/diffvars_impl.h>
13 johnpye 1317 # include <system/jacobian.h>
14 jpye 1337 # include <system/cond_config.h>
15 johnpye 1309 # include <solver/slvDOF.h>
16 johnpye 1239 #endif
17    
18 johnpye 1309 #include "ida_impl.h"
19    
20 johnpye 1271 /* #define ANALYSE_DEBUG */
21 johnpye 1240
22 johnpye 1249 #define VARMSG(MSG) \
23     varname = var_make_name(sys->system,v); \
24     CONSOLE_DEBUG(MSG,varname); \
25     ASC_FREE(varname)
26    
27 johnpye 1239 static int integrator_ida_check_partitioning(IntegratorSystem *sys);
28     static int integrator_ida_check_diffindex(IntegratorSystem *sys);
29 johnpye 1249 /* static int integrator_ida_rebuild_diffindex(IntegratorSystem *sys); */
30 johnpye 1239
31 johnpye 1240 const var_filter_t integrator_ida_nonderiv = {
32 johnpye 1246 VAR_SVAR | VAR_ACTIVE | VAR_FIXED | VAR_DERIV,
33     VAR_SVAR | VAR_ACTIVE | 0 | 0
34 johnpye 1240 };
35    
36     const var_filter_t integrator_ida_deriv = {
37     VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | VAR_FIXED | VAR_DERIV,
38     VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | 0 | VAR_DERIV
39     };
40    
41     const rel_filter_t integrator_ida_rel = {
42     REL_INCLUDED | REL_EQUALITY | REL_ACTIVE,
43     REL_INCLUDED | REL_EQUALITY | REL_ACTIVE
44     };
45    
46     /**
47 johnpye 1252 This is the first step in the DAE analysis process. We inspect the
48     derivative chains from the solver's analyse() routine.
49 johnpye 1240 Check derivative chains: if a var or derivative is inelligible,
50     follow down the chain fixing & setting zero any derivatives.
51     */
52     static int integrator_ida_check_vars(IntegratorSystem *sys){
53     const SolverDiffVarCollection *diffvars;
54 johnpye 1246 char *varname;
55 johnpye 1240 int n_y = 0;
56     int i, j;
57     struct var_variable *v;
58     SolverDiffVarSequence seq;
59 johnpye 1246 int vok;
60 johnpye 1240
61 johnpye 1271 #ifdef ANALYSE_DEBUG
62 johnpye 1241 CONSOLE_DEBUG("BEFORE CHECKING VARS");
63     integrator_ida_analyse_debug(sys,stderr);
64 johnpye 1271 #endif
65 johnpye 1241
66 johnpye 1240 /* we shouldn't have allocated these yet: just be sure */
67     asc_assert(sys->y==NULL);
68     asc_assert(sys->ydot==NULL);
69     asc_assert(sys->n_y==0);
70    
71     /* get the the dervative chains from the system */
72 johnpye 1249 diffvars = system_get_diffvars(sys->system);
73 johnpye 1240 if(diffvars==NULL){
74     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");
75     return 1;
76     }
77 johnpye 1246
78 johnpye 1271 #ifdef ANALYSE_DEBUG
79 johnpye 1246 system_var_list_debug(sys->system);
80 johnpye 1271 #endif
81 johnpye 1240
82     /* add the variables from the derivative chains */
83     for(i=0; i<diffvars->nseqs; ++i){
84 johnpye 1246 seq = diffvars->seqs[i];
85     asc_assert(seq.n >= 1);
86     v = seq.vars[0];
87 johnpye 1241
88 johnpye 1265 /* update the var_fixed flag */
89     (void)var_fixed(v);
90 johnpye 1246 vok = var_apply_filter(v,&integrator_ida_nonderiv);
91 johnpye 1265
92 johnpye 1246 if(vok && !var_incident(v)){
93 johnpye 1270 /* VARMSG("good var '%s' is not incident"); */
94 johnpye 1246 /* var meets our filter, but perhaps it's not incident? */
95     if(seq.n == 1 || var_apply_filter(seq.vars[1],&integrator_ida_nonderiv)){
96 johnpye 1270 /* VARMSG("DEACTIVATING NON-INCIDENT VAR '%s' (NO DERIVATIVE)"); */
97 johnpye 1246 var_set_active(v,0);
98     vok = 0;
99     }else{
100 johnpye 1270 /* VARMSG("'%s' has a derivative that's OK"); */
101 johnpye 1246 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.");
102     return 1;
103     }
104     }
105 johnpye 1241
106 johnpye 1246 if(!vok){
107 johnpye 1270 /* VARMSG("'%s' fails non-deriv filter");
108 johnpye 1265 if(var_fixed(v)){
109     CONSOLE_DEBUG("(var is fixed");
110     }
111     CONSOLE_DEBUG("passes nonderiv? %s (flags = 0x%x)"
112     , (var_apply_filter(v,&integrator_ida_nonderiv) ? "TRUE" : "false")
113 johnpye 1270 , var_flags(v)
114     );
115     */
116 johnpye 1240 for(j=1;j<seq.n;++j){
117     v = seq.vars[j];
118     var_set_active(v,FALSE);
119     var_set_value(v,0);
120 johnpye 1270 /* VARMSG("Derivative '%s' SET ZERO AND INACTIVE"); */
121 johnpye 1240 }
122     continue;
123 johnpye 1246 }
124 johnpye 1240
125 johnpye 1254 /* VARMSG("We will use var '%s'"); */
126 johnpye 1240 if(seq.n > 2){
127     varname = var_make_name(sys->system,v);
128 johnpye 1242 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Too-long derivative chain with var '%s'",varname);
129 johnpye 1240 ASC_FREE(varname);
130     return 2;
131     }else if(seq.n==2){
132     /* diff var */
133     if(var_apply_filter(seq.vars[1],&integrator_ida_deriv)){
134     /* add the diff & deriv vars to the lists */
135     n_y++;
136 johnpye 1270 /* VARMSG("Added diff var '%s'");
137     v = seq.vars[1]; VARMSG("and its derivative '%s'"); */
138 johnpye 1240 continue;
139     }
140 johnpye 1270 /* VARMSG("Diff var '%s' being converted to alg var...");
141     v = seq.vars[1]; VARMSG("...because deriv var '%s' fails filter"); */
142 johnpye 1240 /* fall through */
143     }
144    
145 johnpye 1270 /* VARMSG("Adding '%s' to algebraic"); */
146 johnpye 1240 n_y++;
147     }
148    
149     /* we assert that all vars in y meet the integrator_ida_nonderiv filter */
150     /* we assert that all vars in ydot meet the integrator_ida_deriv filter */
151    
152 johnpye 1271 #ifdef ANALYSE_DEBUG
153 johnpye 1240 CONSOLE_DEBUG("Found %d good non-derivative vars", n_y);
154 johnpye 1271 #endif
155 johnpye 1240 sys->n_y = n_y;
156 johnpye 1241
157 johnpye 1240 return 0;
158     }
159    
160 johnpye 1309 /**
161     Flag relations that contain derivatives as 'REL_DIFFERENTIAL'
162     */
163     static int integrator_ida_flag_rels(IntegratorSystem *sys){
164     int i, n, c, nd=0;
165     struct rel_relation **rels;
166     n = slv_get_num_solvers_rels(sys->system);
167     rels = slv_get_solvers_rel_list(sys->system);
168     for(i=0;i<n;++i){
169     c = rel_classify_differential(rels[i]);
170     if(c){
171     nd++;
172     /* CONSOLE_DEBUG("Rel %d is DIFFERENTIAL", i); */
173     }
174     }
175     CONSOLE_DEBUG("Found %d differential equations (so %d algebraic)",nd, n - nd);
176     sys->n_diffeqs = nd;
177     return 0;
178     }
179 johnpye 1240
180     /**
181     Sort the lists. First we will put the non-derivative vars, then we will
182     put the derivative vars. Then we will put all the others.
183     */
184     static int integrator_ida_sort_rels_and_vars(IntegratorSystem *sys){
185     int ny1, nydot, nr;
186    
187 johnpye 1241
188 johnpye 1271 #ifdef ANALYSE_DEBUG
189 johnpye 1241 CONSOLE_DEBUG("BEFORE SORTING RELS AND VARS");
190     integrator_ida_analyse_debug(sys,stderr);
191 johnpye 1271 #endif
192 johnpye 1241
193 johnpye 1240 /* we should not have allocated y or ydot yet */
194     asc_assert(sys->y==NULL && sys->ydot==NULL);
195    
196     /* but we should have found some variables (and know how many) */
197     asc_assert(sys->n_y);
198    
199     if(system_cut_vars(sys->system, 0, &integrator_ida_nonderiv, &ny1)){
200     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting non-derivs");
201     return 1;
202     }
203    
204 johnpye 1271 #ifdef ANALYSE_DEBUG
205 johnpye 1240 CONSOLE_DEBUG("cut_vars: ny1 = %d, sys->n_y = %d",ny1,sys->n_y);
206 johnpye 1271 #endif
207 johnpye 1240 asc_assert(ny1 == sys->n_y);
208    
209     if(system_cut_vars(sys->system, ny1, &integrator_ida_deriv, &nydot)){
210     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting derivs");
211     return 1;
212     }
213    
214     if(system_cut_rels(sys->system, 0, &integrator_ida_rel, &nr)){
215     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting derivs");
216     return 1;
217     }
218    
219     if(ny1 != nr){
220     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem is not square (ny = %d, nr = %d)",ny1,nr);
221     return 2;
222     }
223    
224     return 0;
225     }
226    
227     /**
228     Build up the lists 'ydot' and 'y_id.
229     'ydot' allows us to find the derivative of a variable given its sindex.
230     'y_id' allows us to locate the variable of a derivative given its sindex.
231    
232     Hence
233     y_id[var_sindex(derivvar)-sys->n_y] = var_sindex(diffvar)
234     ydot[var_sindex(diffvar)] = derivvar (or NULL if diffvar is algebraic)
235    
236     Note that there is 'shifting' required when looking up y_id.
237    
238     Note that we want to get rid of 'y' but for the moment we are also assigning
239     that. We want to be rid of it because if the vars are ordered correctly
240     it shouldn't be needed.
241     */
242     static int integrator_ida_create_lists(IntegratorSystem *sys){
243     const SolverDiffVarCollection *diffvars;
244     int i, j;
245     struct var_variable *v;
246    
247     SolverDiffVarSequence seq;
248    
249     asc_assert(sys->y==NULL);
250     asc_assert(sys->ydot==NULL);
251     asc_assert(sys->y_id== NULL);
252    
253     /* get the the dervative chains from the system */
254 johnpye 1249 diffvars = system_get_diffvars(sys->system);
255 johnpye 1240 if(diffvars==NULL){
256     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");
257     return 1;
258     }
259    
260     /* allocate space (more than enough) */
261     sys->y = ASC_NEW_ARRAY(struct var_variable *,sys->n_y);
262     sys->ydot = ASC_NEW_ARRAY_CLEAR(struct var_variable *,sys->n_y);
263 johnpye 1250 sys->n_ydot = 0;
264 johnpye 1268 for(i=0;i<sys->n_y;++i){
265     asc_assert(sys->ydot[i] == 0);
266     }
267 johnpye 1240
268 johnpye 1271 #ifdef ANALYSE_DEBUG
269 johnpye 1241 CONSOLE_DEBUG("Passing through chains...");
270 johnpye 1271 #endif
271 johnpye 1250 /* create the lists y and ydot, ignoring 'bad' vars */
272 johnpye 1240 for(i=0; i<diffvars->nseqs; ++i){
273 johnpye 1268 /* CONSOLE_DEBUG("i = %d",i); */
274 johnpye 1241
275 johnpye 1240 seq = diffvars->seqs[i];
276     asc_assert(seq.n >= 1);
277     v = seq.vars[0];
278 johnpye 1250 j = var_sindex(v);
279 johnpye 1241
280 johnpye 1240 if(!var_apply_filter(v,&integrator_ida_nonderiv)){
281     continue;
282     }
283    
284 johnpye 1250 sys->y[j] = v;
285 johnpye 1268 /* VARMSG("'%s' is good non-deriv"); */
286 johnpye 1241
287 johnpye 1240 if(seq.n > 1 && var_apply_filter(seq.vars[1],&integrator_ida_deriv)){
288 johnpye 1250 v = seq.vars[1];
289     asc_assert(var_sindex(v) >= sys->n_y);
290     asc_assert(var_sindex(v)-sys->n_y < sys->n_y);
291 johnpye 1241
292 johnpye 1250 sys->ydot[j] = v;
293     sys->n_ydot++;
294 johnpye 1268 /* VARMSG("'%s' is good deriv"); */
295 johnpye 1240 }else{
296     asc_assert(sys->ydot[j]==NULL);
297     }
298 johnpye 1250 }
299 johnpye 1240
300 johnpye 1271 #ifdef ANALYSE_DEBUG
301 johnpye 1250 CONSOLE_DEBUG("Found %d good non-derivs",j);
302 johnpye 1271 #endif
303 johnpye 1250 /* create the list y_id by looking at non-NULLs from ydot */
304     sys->y_id = ASC_NEW_ARRAY(int,sys->n_ydot);
305     for(i=0,j=0; i < sys->n_y; ++i){
306     if(sys->ydot[i]==NULL)continue;
307 johnpye 1270 /* v = sys->ydot[i]; VARMSG("deriv '%s'..."); */
308     /* v = sys->y[i]; VARMSG("diff '%s'..."); */
309 johnpye 1250 sys->y_id[var_sindex(sys->ydot[i]) - sys->n_y] = i;
310 johnpye 1240 j++;
311     }
312 johnpye 1250
313     asc_assert(j==sys->n_ydot);
314 johnpye 1240
315     return 0;
316     }
317    
318 johnpye 1309 /**
319 johnpye 1317 Check for index-1 DAE system. We plan to do this by checking that df/dyd'
320     and dg/dya are both non-singular (and of course square). Valid systems
321     can (we think?) be written that don't meet these requirements but they
322     will have index problems and may not solve with IDA.
323 johnpye 1309 */
324     int integrator_ida_check_index(IntegratorSystem *sys){
325 jpye 1324 #if 1
326 johnpye 1309 linsolqr_system_t L;
327 jpye 1324 mtx_range_t range;
328 johnpye 1309 mtx_region_t R;
329 johnpye 1317 int res, r;
330     struct SystemJacobianStruct df_dydp, dg_dya;
331 johnpye 1309
332 johnpye 1317 CONSOLE_DEBUG("system has total of %d rels and %d vars"
333     ,slv_get_num_solvers_rels(sys->system)
334     ,slv_get_num_solvers_vars(sys->system)
335     );
336    
337     CONSOLE_DEBUG("VAR_DERIV = 0x%x = %d",VAR_DERIV, VAR_DERIV);
338     CONSOLE_DEBUG("system_vfilter_deriv.matchbits = 0x%x",system_vfilter_deriv.matchbits);
339     CONSOLE_DEBUG("system_vfilter_deriv.matchvalue= 0x%x",system_vfilter_deriv.matchvalue);
340    
341     asc_assert(system_vfilter_deriv.matchbits & VAR_DERIV);
342     asc_assert(system_vfilter_deriv.matchvalue & VAR_DERIV);
343    
344     res = system_jacobian(sys->system
345     , &system_rfilter_diff
346     , &system_vfilter_deriv
347     , 1 /* safe evaluation */
348     , &df_dydp
349     );
350    
351     if(res){
352     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error calculating df/dyd'");
353 johnpye 1309 }
354 johnpye 1317 CONSOLE_DEBUG("df/dyd': nr = %d, nv = %d",df_dydp.n_rels,df_dydp.n_vars);
355 johnpye 1309
356 johnpye 1317 res = system_jacobian(sys->system
357     , &system_rfilter_algeb
358     , &system_vfilter_algeb
359     , 1 /* safe evaluation */
360     , &dg_dya
361     );
362    
363     if(res){
364     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error calculating dg/dya");
365 johnpye 1309 }
366 johnpye 1317 CONSOLE_DEBUG("dg/dya: nr = %d, nv = %d",dg_dya.n_rels,dg_dya.n_vars);
367    
368     if((df_dydp.n_rels == 0) ^ (df_dydp.n_vars == 0)){
369     ERROR_REPORTER_HERE(ASC_PROG_ERR,"df/dyd' is a bit ambiguous");
370 johnpye 1309 }
371 johnpye 1317
372     if(dg_dya.n_rels <= 0){
373     ERROR_REPORTER_HERE(ASC_PROG_WARNING,"No algebraic equations were found in the DAE system!");
374     }else if(dg_dya.n_rels != dg_dya.n_vars){
375     ERROR_REPORTER_HERE(ASC_PROG_WARNING,"The algebraic part of the DAE jacobian, dg/dya, is not square!");
376     }else{
377     /* check the rank */
378 jpye 1324 range.low = 0; range.high = mtx_order(dg_dya.M) - 1;
379     R.row = range; R.col = range;
380 johnpye 1309
381 johnpye 1317 L = linsolqr_create_default();
382     linsolqr_set_matrix(L,dg_dya.M);
383     linsolqr_set_region(L,R);
384     linsolqr_prep(L,linsolqr_fmethod_to_fclass(linsolqr_fmethod(L)));
385     linsolqr_reorder(L, &R, linsolqr_rmethod(L));
386     linsolqr_factor(L,linsolqr_fmethod(L));
387     r = linsolqr_rank(L);
388 johnpye 1309
389 johnpye 1317 linsolqr_destroy(L);
390    
391     if(r != dg_dya.n_rels){
392     ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Your DAE system has an index problem: the matrix dg/dya is not full rank");
393     }
394     }
395 johnpye 1309
396 jpye 1324 ASC_FREE(dg_dya.vars);
397     ASC_FREE(dg_dya.rels);
398     mtx_destroy(dg_dya.M);
399    
400 johnpye 1317 if(df_dydp.n_rels <= 0){
401     ERROR_REPORTER_HERE(ASC_PROG_WARNING,"No differential equations were found in the DAE system!");
402     }else if(df_dydp.n_rels != df_dydp.n_vars){
403     ERROR_REPORTER_HERE(ASC_PROG_WARNING,"The differential part of the the jacobian dg/dya is not square!");
404 jpye 1324 ASC_FREE(df_dydp.vars);
405     ASC_FREE(df_dydp.rels);
406     mtx_destroy(df_dydp.M);
407     return 1;
408 johnpye 1317 }else{
409     /* check the rank */
410 jpye 1324 range.low = 0; range.high = mtx_order(df_dydp.M) - 1;
411     R.row = range; R.col = range;
412 johnpye 1309
413 johnpye 1317 L = linsolqr_create_default();
414 jpye 1324 linsolqr_set_matrix(L,df_dydp.M);
415 johnpye 1317 linsolqr_set_region(L,R);
416     linsolqr_prep(L,linsolqr_fmethod_to_fclass(linsolqr_fmethod(L)));
417     linsolqr_reorder(L, &R, linsolqr_rmethod(L));
418     linsolqr_factor(L,linsolqr_fmethod(L));
419     r = linsolqr_rank(L);
420 johnpye 1309
421 johnpye 1317 linsolqr_destroy(L);
422    
423     if(r != df_dydp.n_rels){
424     ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Your DAE system has an index problem: the matrix df/dyd' is not full rank");
425     }
426 johnpye 1309 }
427 johnpye 1317
428     if(df_dydp.n_rels + dg_dya.n_rels == 0){
429     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Both df/dyd' and dg/dya were empty!");
430     }
431    
432     ASC_FREE(df_dydp.vars);
433     ASC_FREE(df_dydp.rels);
434     mtx_destroy(df_dydp.M);
435     return 0;
436 jpye 1320 #else
437     ERROR_REPORTER_HERE(ASC_PROG_ERR,"check_index disabled");
438     return 0;
439     #endif
440 johnpye 1309 }
441    
442 johnpye 1239 /**
443     Perform additional problem analysis to prepare problem for integration with
444     IDA.
445    
446     We assume that the analyse_generate_diffvars routine has been called, so
447     that we just need to call slv_get_diffvars to access the derivative
448     chains.
449    
450     We can also assume that the independent variable has been found.
451    
452     See mailing list, ~Jan 2007.
453    
454     Note, the stuff for identifying the static and output sub-problems should
455     be part of integrator.c, not this file. We will assume this is handled
456    
457     @return 0 on success
458     @see integrator_analyse
459     */
460     int integrator_ida_analyse(IntegratorSystem *sys){
461 johnpye 1240 int res;
462 johnpye 1239 const SolverDiffVarCollection *diffvars;
463 johnpye 1271 int i;
464     #ifdef ANALYSE_DEBUG
465 johnpye 1239 char *varname;
466 johnpye 1271 #endif
467 johnpye 1239
468     asc_assert(sys->engine==INTEG_IDA);
469    
470 jpye 1337 CONSOLE_DEBUG("System contains a total of %d bnds and %d reals"
471     ,slv_get_num_solvers_bnds(sys->system)
472     ,slv_get_num_solvers_rels(sys->system)
473     );
474    
475     /* set the active flags on variables depending on the state of WHENs */
476     CONSOLE_DEBUG("Currently %d rels active",slv_count_solvers_rels(sys->system, &integrator_ida_rel));
477    
478     reanalyze_solver_lists(sys->system);
479    
480     CONSOLE_DEBUG("After analysing WHENs, there are %d rels active"
481     ,slv_count_solvers_rels(sys->system, &integrator_ida_rel)
482     );
483    
484    
485 johnpye 1239 #ifdef ANALYSE_DEBUG
486     CONSOLE_DEBUG("Starting IDA analysis");
487     #endif
488    
489 jpye 1337 /* set the flags on differential and derivative and algebraic vars */
490 johnpye 1240 res = integrator_ida_check_vars(sys);
491     if(res){
492     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem with vars");
493 johnpye 1239 return 1;
494     }
495 johnpye 1240
496 johnpye 1309 res = integrator_ida_flag_rels(sys);
497     if(res){
498     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem classifying differential equations");
499     return 1;
500     }
501    
502 johnpye 1240 #ifdef ANALYSE_DEBUG
503     CONSOLE_DEBUG("Sorting rels and vars");
504 johnpye 1239 #endif
505    
506 johnpye 1240 res = integrator_ida_sort_rels_and_vars(sys);
507     if(res){
508     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem sorting rels and vars");
509     return 1;
510 johnpye 1239 }
511    
512 johnpye 1240 #ifdef ANALYSE_DEBUG
513     CONSOLE_DEBUG("Creating lists");
514 johnpye 1239
515 johnpye 1249 CONSOLE_DEBUG("BEFORE MAKING LISTS");
516     integrator_ida_debug(sys,stderr);
517 johnpye 1271 #endif
518 johnpye 1249
519 johnpye 1240 res = integrator_ida_create_lists(sys);
520     if(res){
521     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem creating lists");
522 johnpye 1239 return 1;
523     }
524    
525 jpye 1337 CONSOLE_DEBUG("After ida_create_lists, there are %d rels active"
526     ,slv_count_solvers_rels(sys->system, &integrator_ida_rel)
527     );
528    
529 johnpye 1240 #ifdef ANALYSE_DEBUG
530 johnpye 1249 CONSOLE_DEBUG("Checking lists");
531 johnpye 1239
532 johnpye 1240 asc_assert(sys->y);
533     asc_assert(sys->ydot);
534     asc_assert(sys->y_id);
535 johnpye 1239
536 johnpye 1249 integrator_ida_debug(sys,stderr);
537 johnpye 1271 #endif
538 johnpye 1239
539     if(integrator_ida_check_diffindex(sys)){
540     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error with diffindex");
541     return 360;
542     }
543    
544 jpye 1337
545     CONSOLE_DEBUG("After ida_check_diffindex, there are %d rels active"
546     ,slv_count_solvers_rels(sys->system, &integrator_ida_rel)
547     );
548    
549 johnpye 1239 /* the following causes TestIDA.testincidence3 to fail:
550     if(sys->n_y != slv_get_num_solvers_rels(sys->system)){
551     ERROR_REPORTER_HERE(ASC_USER_ERROR,"Problem is not square");
552     return 2;
553     }*/
554    
555     #if 0
556     /* check structural singularity for the two IDACalcIC scenarios */
557    
558     /* ...(1) FIX the derivatives */
559     CONSOLE_DEBUG("Checking system with derivatives fixed...");
560     for(i=0;i<n_y;++i){
561     if(sys->ydot[i])var_set_fixed(sys->ydot[i],1);
562     }
563    
564     slv_block_partition(sys->system);
565     res = integrator_ida_block_check(sys);
566     if(res)return 100 + res;
567    
568     /* ...(2) FREE the derivatives, FIX the diffvars */
569     CONSOLE_DEBUG("Checking system with differential variables fixed...");
570     for(i=0;i<n_y;++i){
571     if(sys->ydot[i]){
572     var_set_fixed(sys->ydot[i],0);
573     var_set_fixed(sys->y[i],1);
574     }
575     }
576     slv_block_partition(sys->system);
577     res = integrator_ida_block_check(sys);
578     if(res)return 200 + res;
579    
580     /* ...(3) restore as it was, FREE the diffvars */
581     for(i=0;i<n_y;++i){
582     if(sys->ydot[i]){
583     var_set_fixed(sys->y[i],0);
584     }
585     }
586     #endif
587    
588 johnpye 1309 res = integrator_ida_check_index(sys);
589     if(res){
590 jpye 1324 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Your DAE system has an index problem");
591 johnpye 1309 return 100 + res;
592     }
593    
594 jpye 1337 CONSOLE_DEBUG("After ida_check_index, there are %d rels active"
595     ,slv_count_solvers_rels(sys->system, &integrator_ida_rel)
596     );
597    
598    
599 johnpye 1240 /* get the the dervative chains from the system */
600 johnpye 1249 diffvars = system_get_diffvars(sys->system);
601 johnpye 1240
602 johnpye 1309 /* check the indep var is same as was located earlier */
603 johnpye 1239 if(diffvars->nindep>1){
604     ERROR_REPORTER_HERE(ASC_USER_ERROR,"Multiple variables specified as independent (ode_type=-1)");
605     return 3;
606     }else if(diffvars->nindep<1){
607     ERROR_REPORTER_HERE(ASC_USER_ERROR,"Independent var not set (ode_type=-1)");
608     return 4;
609     }else if(diffvars->indep[0]!=sys->x){
610     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Indep var doesn't match");
611     return 5;
612     }
613    
614 johnpye 1249 #ifdef ANALYSE_DEBUG
615     CONSOLE_DEBUG("Collecting observed variables");
616     #endif
617    
618 johnpye 1239 /* get the observations */
619     /** @TODO this should use a 'problem provider' API instead of static data
620     from the system object */
621     sys->n_obs = diffvars->nobs;
622     sys->obs = ASC_NEW_ARRAY(struct var_variable *,sys->n_obs);
623     for(i=0;i<sys->n_obs;++i){
624 johnpye 1267 /* we get them all, regardless of flags etc */
625 johnpye 1239 sys->obs[i] = diffvars->obs[i];
626 johnpye 1271 #ifdef ANALYSE_DEBUG
627 johnpye 1239 varname = var_make_name(sys->system,sys->obs[i]);
628     CONSOLE_DEBUG("'%s' is observation",varname);
629     ASC_FREE(varname);
630 johnpye 1271 #endif
631 johnpye 1239 }
632    
633 jpye 1337 CONSOLE_DEBUG("rels matchbits: 0x%x",integrator_ida_rel.matchbits);
634     CONSOLE_DEBUG("rels matchvalue: 0x%x",integrator_ida_rel.matchvalue);
635    
636     CONSOLE_DEBUG("At the end of ida_analyse, there are %d rels active"
637     ,slv_count_solvers_rels(sys->system, &integrator_ida_rel)
638     );
639    
640 johnpye 1239 return 0;
641     }
642    
643     /*------------------------------------------------------------------------------
644     ANALYSIS ROUTINE (new implementation)
645     */
646    
647     static int integrator_ida_check_partitioning(IntegratorSystem *sys){
648     long i, nv;
649     int err=0;
650     char *varname;
651     struct var_variable **vlist, *v;
652     nv = slv_get_num_solvers_vars(sys->system);
653     vlist = slv_get_solvers_var_list(sys->system);
654     var_filter_t vf = {VAR_SVAR|VAR_ACTIVE|VAR_INCIDENT|VAR_FIXED
655     ,VAR_SVAR|VAR_ACTIVE|VAR_INCIDENT|0 };
656     for(i=0;i<nv;++i){
657     v=vlist[i];
658     if(!var_apply_filter(v,&vf))continue;
659     varname = var_make_name(sys->system,v);
660     if(!var_deriv(v)){
661 johnpye 1271 #ifdef ANALYSE_DEBUG
662 johnpye 1239 fprintf(stderr,"vlist[%ld] = '%s' (nonderiv)\n",i,varname);
663 johnpye 1271 #endif
664 johnpye 1239 if(i>=sys->n_y){
665     ERROR_REPORTER_HERE(ASC_PROG_ERR,"non-deriv var '%s' is not at the start",varname);
666     err++;
667     }
668     }else{
669 johnpye 1271 #ifdef ANALYSE_DEBUG
670 johnpye 1239 fprintf(stderr,"vlist[%ld] = '%s' (derivative)\n",i,varname);
671 johnpye 1271 #endif
672 johnpye 1239 if(i<sys->n_y){
673     ERROR_REPORTER_HERE(ASC_PROG_ERR,"deriv var '%s' is not at the end (n_y = %d, i = %d)"
674     ,varname, sys->n_y, i
675     );
676     err++;
677     }
678     }
679     ASC_FREE(varname);
680     }
681     if(err){
682     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error in IDA partitioning");
683     return 1;
684     }
685     return 0;
686     }
687    
688     int integrator_ida_block_check(IntegratorSystem *sys){
689     int res;
690     int dof;
691     #ifdef ANALYSE_DEBUG
692 johnpye 1246 int *vlist, *vp, i, nv, nv_ok;
693 johnpye 1239 char *varname;
694     struct var_variable **solversvars;
695    
696     var_filter_t vfilt = {
697     VAR_ACTIVE | VAR_INCIDENT | VAR_FIXED,
698     VAR_ACTIVE | VAR_INCIDENT | 0
699     };
700    
701     nv = slv_get_num_solvers_vars(sys->system);
702     solversvars = slv_get_solvers_var_list(sys->system);
703 johnpye 1246 CONSOLE_DEBUG("-------------- nv = %d -------------",nv);
704 johnpye 1239 for(nv_ok=0, i=0; i < nv; ++i){
705     if(var_apply_filter(solversvars[i],&vfilt)){
706     varname = var_make_name(sys->system,solversvars[i]);
707     fprintf(stderr,"%s\n",varname);
708     ASC_FREE(varname);
709     nv_ok++;
710     }
711     }
712 johnpye 1246 CONSOLE_DEBUG("----------- got %d ok -------------",nv_ok);
713 johnpye 1239 #endif
714    
715     if(!slvDOF_status(sys->system, &res, &dof)){
716     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to determine DOF status");
717     return -1;
718     }
719     switch(res){
720     case 1: CONSOLE_DEBUG("System is underspecified (%d degrees of freedom)",dof);break;
721     case 2: CONSOLE_DEBUG("System is square"); return 0; /* all OK */
722     case 3: CONSOLE_DEBUG("System is structurally singular"); break;
723     case 4: CONSOLE_DEBUG("System is overspecified"); break;
724     default:
725     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unrecognised slfDOF_status");
726     return -2;
727     }
728    
729     #ifdef ANALYSE_DEBUG
730     /* if it was underspecified, what vars could be fixed? */
731     if(res==1){
732     CONSOLE_DEBUG("Need to FIX %d of the following vars:",dof);
733     solversvars = slv_get_solvers_var_list(sys->system);
734     if(!slvDOF_eligible(sys->system, &vlist)){
735     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to det slvDOF_eligble list");
736     return -3;
737     }
738     for(vp=vlist;*vp!=-1;++vp){
739     varname = var_make_name(sys->system, solversvars[*vp]);
740     CONSOLE_DEBUG("Fixable var: %s",varname);
741     ASC_FREE(varname);
742     }
743 johnpye 1240 CONSOLE_DEBUG("(Found %d fixable vars)",(int)(vp-vlist));
744 johnpye 1239 return 1;
745     }
746     #endif
747    
748     return res;
749     }
750    
751 johnpye 1249 /* return 0 on succes */
752 johnpye 1268 static int check_dups(IntegratorSystem *sys, struct var_variable **list,int n,int allownull){
753 johnpye 1249 int i,j;
754     struct var_variable *v;
755 johnpye 1271 #ifdef ANALYSE_DEBUG
756 johnpye 1268 char *varname;
757 johnpye 1271 #endif
758 johnpye 1249 for(i=0; i< n; ++i){
759     v=list[i];
760     if(v==NULL){
761     if(allownull)continue;
762     else return 2;
763     }
764 johnpye 1271 asc_assert(v!=(void *)0x31);
765 johnpye 1249 for(j=0; j<i-1;++j){
766     if(list[j]==NULL)continue;
767 johnpye 1268 if(v==list[j]){
768 johnpye 1271 #ifdef ANALYSE_DEBUG
769 johnpye 1268 varname = var_make_name(sys->system,v);
770     if(varname){
771     CONSOLE_DEBUG("Duplicate of '%s' found",varname);
772     ASC_FREE(varname);
773     }else{
774     CONSOLE_DEBUG("Duplicate found (couldn't retrieve name)");
775     }
776 johnpye 1270 ASC_FREE(varname);
777 johnpye 1271 #endif
778 johnpye 1268 return 1;
779     }
780 johnpye 1249 }
781     }
782     return 0;
783     }
784    
785     /*
786     We are going to check that
787    
788     - n_y in range (0,n_vars)
789     - no duplicates anywhere in the varlist
790     - no duplicates in non-NULL elements of ydot
791    
792     - first n_y vars in solver's var list match those in vector y and match the non-deriv filter.
793     - var_sindex for first n_y vars match their position the the solver's var list
794    
795     - ydot contains n_ydot non-NULL elements, each match the deriv filter.
796     - vars in solver's var list positions [n_y,n_y+n_ydot) all match deriv filter
797    
798     - y_id contains n_ydot elements, with int values in range [0,n_y)
799     - var_sindex(ydot[y_id[i]]) - n_y == i for i in [0,n_ydot)
800    
801     - all the vars [n_y+n_ydot,n_vars) fail the non-deriv filter and fail the deriv filter.
802     */
803 johnpye 1239 static int integrator_ida_check_diffindex(IntegratorSystem *sys){
804 johnpye 1249 int i, n_vars, n_ydot=0;
805     struct var_variable **list, *v;
806     char *varname;
807     const char *msg;
808 johnpye 1239
809 johnpye 1271 #ifdef ANALYSE_DEBUG
810 johnpye 1239 CONSOLE_DEBUG("Checking diffindex vector");
811 johnpye 1271 #endif
812 johnpye 1239
813 johnpye 1249 if(sys->y_id == NULL || sys->y == NULL || sys->ydot == NULL){
814     ERROR_REPORTER_HERE(ASC_PROG_ERR,"list(s) NULL");
815 johnpye 1239 return 1;
816     }
817 johnpye 1249
818     list = slv_get_solvers_var_list(sys->system);
819     n_vars = slv_get_num_solvers_vars(sys->system);
820    
821 johnpye 1268 if(check_dups(sys, list, n_vars,FALSE)){
822 johnpye 1249 ERROR_REPORTER_HERE(ASC_PROG_ERR,"duplicates or NULLs in solver's var list");
823     return 1;
824     }
825 johnpye 1239
826 johnpye 1270 if(check_dups(sys, sys->ydot, sys->n_y,TRUE)){
827 johnpye 1249 ERROR_REPORTER_HERE(ASC_PROG_ERR,"duplicates in ydot vector");
828     return 1;
829     }
830 johnpye 1239
831 johnpye 1249 /* check n_y in range */
832     if(sys->n_y <= 0 || sys->n_y >= n_vars){
833     ERROR_REPORTER_HERE(ASC_PROG_ERR,"n_y = %d invalid (n_vars = %d)",sys->n_y, n_vars);
834     return 1;
835     }
836    
837     /* check first n_y vars */
838     for(i=0; i < sys->n_y; ++i){
839     v = list[i];
840     if(!var_apply_filter(v, &integrator_ida_nonderiv)){
841     msg = "'%s' (in first n_y vars) fails non-deriv filter"; goto finish;
842 johnpye 1250 }else if(var_sindex(v) != i){
843     msg = "'%s' has wrong var_sindex"; goto finish;
844 johnpye 1249 }else if(v != sys->y[i]){
845     msg = "'%s' not matched in y vector"; goto finish;
846     }
847     /* meets filter, matches in y vector, has correct var_sindex. */
848     }
849    
850     /* count non-NULLs in ydot */
851     for(i=0; i < sys->n_y; ++i){
852     v = sys->ydot[i];
853     if(v!=NULL){
854     if(var_sindex(v) < sys->n_y){
855     msg = "'%s' has var_sindex < n_y"; goto finish;
856     }else if(!var_apply_filter(v,&integrator_ida_deriv)){
857     msg = "'%s' (in next n_ydot vars) fails deriv filter"; goto finish;
858 johnpye 1239 }
859 johnpye 1249 /* lies beyond first n_y vars, match deriv filter */
860     n_ydot++;
861 johnpye 1239 }
862     }
863 johnpye 1249
864     /* check that vars [n_y, n_y+n_ydot) in solver's var list match the deriv filter */
865     for(i=sys->n_y; i< sys->n_y + n_ydot; ++i){
866     v = list[i];
867     if(!var_apply_filter(v,&integrator_ida_deriv)){
868     msg = "'%s', in range [n_y,n_y+n_ydot), fails deriv filter"; goto finish;
869     }
870 johnpye 1239 }
871 johnpye 1249
872     /* check values in y_id are ints int range [0,n_y), and that they point to correct vars */
873     for(i=0; i<n_ydot; ++i){
874     if(sys->y_id[i] < 0 || sys->y_id[i] >= sys->n_y){
875     ERROR_REPORTER_HERE(ASC_PROG_ERR,"value of y_id[%d] is out of range",i);
876     return 1;
877     }
878     v = sys->ydot[sys->y_id[i]];
879     if(var_sindex(v) - sys->n_y != i){
880     msg = "'%s' not index correctly from y_id"; goto finish;
881     }
882     }
883    
884     /* check remaining vars fail both filters */
885     for(i=sys->n_y + n_ydot; i<n_vars; ++i){
886     v = list[i];
887     if(var_apply_filter(v,&integrator_ida_nonderiv)){
888     msg = "Var '%s' at end meets non-deriv filter, but shouldn't"; goto finish;
889     }
890     if(var_apply_filter(v,&integrator_ida_deriv)){
891     CONSOLE_DEBUG("position = %d",i);
892     msg = "Var '%s' at end meets deriv filter, but shouldn't"; goto finish;
893     }
894     }
895    
896     return 0;
897     finish:
898     varname = var_make_name(sys->system,v);
899     ERROR_REPORTER_HERE(ASC_PROG_ERR,msg,varname);
900     ASC_FREE(varname);
901     return 1;
902 johnpye 1239 }
903    
904     /**
905 johnpye 1250 Given a derivative variable, return the index of its corresponding differential
906     variable in the y vector (and equivalently the var_sindex of the diff var)
907 johnpye 1239 */
908     int integrator_ida_diffindex(const IntegratorSystem *sys, const struct var_variable *deriv){
909     asc_assert(var_sindex(deriv) >= sys->n_y);
910 johnpye 1250 asc_assert(var_sindex(deriv) < sys->n_y + sys->n_ydot);
911     return sys->y_id[var_sindex(deriv) - sys->n_y];
912 johnpye 1239 }
913    
914 johnpye 1241 /**
915     This function will output the data structures provided to use BY THE
916     SYSTEM -- not the ones we're working with here IN THE SOLVER.
917     */
918 johnpye 1239 int integrator_ida_analyse_debug(const IntegratorSystem *sys,FILE *fp){
919 johnpye 1247 return system_diffvars_debug(sys->system,fp);
920 johnpye 1239 }
921    

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