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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 921 - (hide annotations) (download) (as text)
Mon Nov 6 07:49:06 2006 UTC (13 years, 9 months ago) by johnpye
File MIME type: text/x-csrc
File size: 43319 byte(s)
Fixing up some malloc/free problems with  integrator.c, removing some debug output from slv3.
Working on fixing IDA for systems with inactive variables (ongoing).
1 johnpye 669 /* ASCEND modelling environment
2     Copyright (C) 2006 Carnegie Mellon University
3    
4     This program is free software; you can redistribute it and/or modify
5     it under the terms of the GNU General Public License as published by
6     the Free Software Foundation; either version 2, or (at your option)
7     any later version.
8    
9     This program is distributed in the hope that it will be useful,
10     but WITHOUT ANY WARRANTY; without even the implied warranty of
11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12     GNU General Public License for more details.
13    
14     You should have received a copy of the GNU General Public License
15     along with this program; if not, write to the Free Software
16     Foundation, Inc., 59 Temple Place - Suite 330,
17     Boston, MA 02111-1307, USA.
18     *//**
19     @file
20     Integrator API for ASCEND, for solving systems of ODEs and/or DAEs.
21     *//*
22     by John Pye, May 2006
23     based on parts of Integrators.c in the Tcl/Tk interface directory, heavily
24     modified to provide a non-GUI-specific API and modularised for multiple
25     integration engines.
26     */
27     #include <time.h>
28     #include "integrator.h"
29     #include "lsode.h"
30     #include "ida.h"
31    
32 ben.allan 704 #include "samplelist.h"
33 johnpye 669
34 ben.allan 704 /**
35     Define as TRUE to enable debug message printing.
36     @TODO this needs to go away.
37     */
38     #define INTEG_DEBUG TRUE
39    
40     /**
41     Print a debug message with value if INTEG_DEBUG is TRUE.
42     */
43     #define print_debug(msg, value) \
44     if(INTEG_DEBUG){ CONSOLE_DEBUG(msg, value); }
45    
46     /**
47     Print a debug message string if INTEG_DEBUG is TRUE.
48     */
49     #define print_debugstring(msg) \
50     if(INTEG_DEBUG){ CONSOLE_DEBUG(msg); }
51    
52 johnpye 669 /*------------------------------------------------------------------------------
53 ben.allan 704 The following names are of solver_var children or attributes
54     * we support (at least temporarily) to determine who is a state and
55     * who matching derivative.
56     * These should be supported directly in a future solveratominst.
57     */
58 johnpye 669
59 ben.allan 704 static symchar *g_symbols[4];
60    
61     #define STATEFLAG g_symbols[0]
62 johnpye 741 /*
63     Integer child. 0= algebraic, 1 = state, 2 = derivative, 3 = 2nd deriv etc
64 johnpye 669 independent variable is -1.
65     */
66     #define INTEG_OTHER_VAR -1L
67     #define INTEG_ALGEBRAIC_VAR 0L
68     #define INTEG_STATE_VAR 1L
69 ben.allan 704
70     #define STATEINDEX g_symbols[1]
71     /* Integer child. all variables with the same STATEINDEX value are taken to
72     * be derivatives of the same state variable. We really need a compiler
73     * that maintains this info by backpointers, but oh well until then.
74 johnpye 669 */
75 ben.allan 704 #define OBSINDEX g_symbols[2]
76     /* Integer child. All variables with OBSINDEX !=0 will be recorded in
77     * the blsode output file. Tis someone else's job to grok this output.
78     */
79 johnpye 669
80 ben.allan 704 #define FIXEDSYMBOL g_symbols[3]
81    
82     /** Temporary catcher of dynamic variable and observation variable data */
83     struct Integ_var_t {
84     long index;
85     long type;
86 johnpye 669 struct var_variable *i;
87     struct Integ_var_t *derivative_of;
88     struct var_variable *derivative;
89 johnpye 912 int varindx; /**< index into slv_get_master_vars_list, or -1 if not there */
90 ben.allan 704 int isstate;
91     };
92 johnpye 669
93     /*------------------------------------------------------------------------------
94     forward declarations
95     */
96    
97     /* abstractions of setup/teardown procedures for the specific solvers */
98 johnpye 912 void integrator_create_engine(IntegratorSystem *sys);
99     void integrator_free_engine(IntegratorSystem *sys);
100 johnpye 669
101 johnpye 912 static int integrator_analyse_ode(IntegratorSystem *sys);
102     static int integrator_analyse_dae(IntegratorSystem *sys);
103 johnpye 669
104 johnpye 912 typedef void (IntegratorVarVisitorFn)(IntegratorSystem *sys, struct var_variable *var, const int *varindx);
105     static void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitor);
106 johnpye 669 IntegratorVarVisitorFn integrator_ode_classify_var;
107     IntegratorVarVisitorFn integrator_dae_classify_var;
108     IntegratorVarVisitorFn integrator_classify_indep_var;
109    
110 johnpye 912 static int integrator_sort_obs_vars(IntegratorSystem *sys);
111     static void integrator_print_var_stats(IntegratorSystem *sys);
112     static int integrator_check_indep_var(IntegratorSystem *sys);
113 johnpye 669
114     static int Integ_CmpDynVars(struct Integ_var_t *v1, struct Integ_var_t *v2);
115     static int Integ_CmpObs(struct Integ_var_t *v1, struct Integ_var_t *v2);
116 johnpye 908 static void Integ_SetObsId(struct var_variable *v, long oindex);
117 johnpye 669
118 johnpye 908 static long DynamicVarInfo(struct var_variable *v,long *vindex);
119     static struct var_variable *ObservationVar(struct var_variable *v, long *oindex);
120 johnpye 669 static void IntegInitSymbols(void);
121    
122     /*------------------------------------------------------------------------------
123     INSTANTIATION AND DESTRUCTION
124     */
125    
126     /**
127     Create a new IntegratorSystem and assign a slv_system_t to it.
128     */
129 johnpye 921 IntegratorSystem *integrator_new(slv_system_t slvsys, struct Instance *inst){
130     IntegratorSystem *sys;
131 johnpye 669
132 johnpye 921 if (slvsys == NULL) {
133     ERROR_REPORTER_HERE(ASC_PROG_ERR,"slvsys is NULL!");
134 ben.allan 704 return NULL;
135 johnpye 669 }
136    
137 johnpye 921 sys = ASC_NEW_CLEAR(IntegratorSystem);
138     sys->system = slvsys;
139     sys->instance = inst;
140    
141     sys->states = NULL; sys->derivs = NULL;
142     sys->dynvars = NULL; sys->obslist = NULL; sys->indepvars = NULL;
143    
144     sys->y_id = NULL;
145     sys->obs_id = NULL;
146     sys->y = NULL;
147     sys->ydot = NULL;
148     sys->obs = NULL;
149     return sys;
150 johnpye 669 }
151    
152     /**
153     Carefully trash any data in the IntegratorSystem that we own,
154     then destroy the IntegratorSystem struct.
155    
156     Note that the integrator doesn't own the samplelist.
157    
158     @param sys will be destroyed and set to NULL.
159 ben.allan 704 */
160     void integrator_free(IntegratorSystem *sys){
161 johnpye 669 if(sys==NULL)return;
162    
163     integrator_free_engine(sys);
164 ben.allan 704
165 johnpye 669 if(sys->states != NULL)gl_destroy(sys->states);
166 ben.allan 704 if(sys->derivs != NULL)gl_destroy(sys->derivs);
167 johnpye 669
168     if(sys->dynvars != NULL)gl_free_and_destroy(sys->dynvars); /* we own the objects in dynvars */
169 ben.allan 704 if(sys->obslist != NULL)gl_free_and_destroy(sys->obslist); /* and obslist */
170     if (sys->indepvars != NULL)gl_free_and_destroy(sys->indepvars); /* and indepvars */
171 johnpye 669
172 johnpye 921 /* if(sys->y_id != NULL)ASC_FREE(sys->y_id); */
173     if(sys->obs_id != NULL)ASC_FREE(sys->obs_id);
174 johnpye 669
175 johnpye 921 if(sys->y != NULL && !sys->ycount)ASC_FREE(sys->y);
176     if(sys->ydot != NULL && !sys->ydotcount)ASC_FREE(sys->ydot);
177     if(sys->obs != NULL && !sys->obscount)ASC_FREE(sys->obs);
178 johnpye 669
179     ascfree(sys);
180 ben.allan 704 sys=NULL;
181 johnpye 669 }
182    
183     /**
184     Utility function to retreive pointers to the symbols we'll be looking for
185     in the instance hierarchy.
186     */
187 ben.allan 704 static void IntegInitSymbols(void){
188     STATEFLAG = AddSymbol("ode_type");
189     STATEINDEX = AddSymbol("ode_id");
190 johnpye 669 OBSINDEX = AddSymbol("obs_id");
191     FIXEDSYMBOL = AddSymbol("fixed");
192     }
193    
194     /*------------------------------------------------------------------------------
195     INTEGRATOR ENGINE
196     */
197    
198 johnpye 709 /* return 1 on success */
199 johnpye 912 int integrator_set_engine(IntegratorSystem *sys, IntegratorEngine engine){
200 johnpye 669
201     /* verify integrator type ok. always passes for nonNULL inst. */
202 ben.allan 704 if(engine==INTEG_UNKNOWN){
203 johnpye 669 ERROR_REPORTER_NOLINE(ASC_USER_ERROR
204     ,"Integrator has not been specified (or is unknown)."
205 ben.allan 704 );
206     return 0;
207 johnpye 669 }
208    
209 johnpye 912 if(engine==sys->engine){
210 johnpye 709 return 1;
211 johnpye 669 }
212 johnpye 912 if(sys->engine!=INTEG_UNKNOWN){
213     integrator_free_engine(sys);
214 johnpye 669 }
215 johnpye 912 sys->engine = engine;
216     integrator_create_engine(sys);
217 johnpye 709
218     return 1;
219 johnpye 669 }
220    
221 johnpye 912 IntegratorEngine integrator_get_engine(const IntegratorSystem *sys){
222     return sys->engine;
223 johnpye 669 }
224    
225     /**
226 johnpye 741 Free any engine-specific data that was required for the solution of
227 johnpye 912 this system. Note that this data is pointed to by sys->enginedata.
228 johnpye 669 */
229 johnpye 912 void integrator_free_engine(IntegratorSystem *sys){
230     switch(sys->engine){
231     case INTEG_LSODE: integrator_lsode_free(sys->enginedata); break;
232 johnpye 669 #ifdef ASC_WITH_IDA
233 johnpye 912 case INTEG_IDA: integrator_ida_free(sys->enginedata); break;
234 johnpye 669 #endif
235 johnpye 709 default: break;
236 johnpye 669 }
237 johnpye 912 sys->enginedata=NULL;
238 johnpye 669 }
239    
240     /**
241     Create enginedata memory if required for this solver. This doesn't include
242     allocating computation space, since we assume that this stage all we know
243     is that we want to use a specified integrator engine, not the full details
244     of the problem at hand. Allocating space inside enginedata should be done
245     during the solve stage (and freed inside integrator_free_engine)
246     */
247 johnpye 912 void integrator_create_engine(IntegratorSystem *sys){
248     if(sys->enginedata!=NULL)return;
249     switch(sys->engine){
250     case INTEG_LSODE: integrator_lsode_create(sys); break;
251 johnpye 669 #ifdef ASC_WITH_IDA
252 johnpye 912 case INTEG_IDA: integrator_ida_create(sys); break;
253 johnpye 669 #endif
254 johnpye 709 default: break;
255 johnpye 669 }
256     }
257    
258     /*------------------------------------------------------------------------------
259     ANALYSIS
260    
261     Provide two modes in order to provide analysis suitable for solution of both
262     ODEs (the previous technique) and DAEs (new code). These share a common
263 johnpye 741 "visit" method that needs to eventually be integrated with the code in
264 johnpye 669 <solver/analyze.c>. For the moment, we're just hacking in to the compiler.
265     */
266    
267     /**
268     Locate the independent variable. For the purpose of GUI design, this needs
269     to work independent of the integration engine being used.
270     */
271 johnpye 912 int integrator_find_indep_var(IntegratorSystem *sys){
272 johnpye 669 int result = 0;
273    
274 johnpye 912 if(sys->x != NULL){
275     CONSOLE_DEBUG("sys->x already set");
276 johnpye 669 return 1;
277     }
278 johnpye 912 assert(sys->indepvars==NULL);
279     sys->indepvars = gl_create(10L);
280 johnpye 669
281 ben.allan 704 IntegInitSymbols();
282 johnpye 669
283 johnpye 725 /* CONSOLE_DEBUG("About to visit..."); */
284 johnpye 912 integrator_visit_system_vars(sys,&integrator_classify_indep_var);
285 johnpye 669
286 johnpye 725 /* CONSOLE_DEBUG("..."); */
287 johnpye 669
288 johnpye 912 result = integrator_check_indep_var(sys);
289     gl_free_and_destroy(sys->indepvars);
290     sys->indepvars = NULL;
291 johnpye 669
292 johnpye 854 /* ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Returning result %d",result); */
293 ben.allan 704
294 johnpye 669 return result;
295 johnpye 741 }
296 johnpye 669
297     /**
298     Analyse the system, either as DAE or as an ODE system, depending on the
299     solver engine selected.
300    
301     @return 1 on success
302     */
303 johnpye 912 int integrator_analyse(IntegratorSystem *sys){
304     switch(sys->engine){
305     case INTEG_LSODE: return integrator_analyse_ode(sys);
306 johnpye 669 #ifdef ASC_WITH_IDA
307 johnpye 912 case INTEG_IDA: return integrator_analyse_dae(sys);
308 johnpye 669 #endif
309     case INTEG_UNKNOWN:
310     ERROR_REPORTER_HERE(ASC_PROG_ERR,"No engine selected: can't analyse");
311     default:
312     ERROR_REPORTER_HERE(ASC_PROG_ERR
313     ,"The selected integration engine (%d) is not available"
314 johnpye 912 ,sys->engine
315 johnpye 669 );
316     }
317     return 0;
318     }
319    
320     /**
321     To analyse a DAE we need to identify *ALL* variables in the system
322     Except for the highest-level derivatives of any present?
323     We also need to identify the independent variable (just one).
324    
325     @TODO implement Pantelides algorithm in here?
326     @TODO prevent re-analysis without clearing out the data structures?
327     @return 1 on success
328     */
329 johnpye 912 int integrator_analyse_dae(IntegratorSystem *sys){
330 johnpye 669 struct Integ_var_t *info, *prev;
331     char *varname, *derivname;
332 johnpye 912 struct var_variable **varlist;
333     int nvarlist;
334 johnpye 669 int i, j;
335     int numstates;
336     int numy, nrels;
337     int maxderiv;
338    
339     CONSOLE_DEBUG("Starting DAE analysis");
340 ben.allan 704 IntegInitSymbols();
341 johnpye 669
342 johnpye 912 assert(sys->indepvars==NULL);
343 johnpye 669
344 johnpye 912 sys->indepvars = gl_create(10L); /* t var info */
345     sys->dynvars = gl_create(200L); /* y and ydot var info */
346     sys->obslist = gl_create(100L); /* obs info */
347 johnpye 669
348 johnpye 912 if(sys->indepvars==NULL
349     || sys->dynvars==NULL
350     || sys->obslist==NULL
351 johnpye 669 ){
352     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory");
353     return 0;
354     }
355    
356 johnpye 912 integrator_visit_system_vars(sys,&integrator_dae_classify_var);
357 johnpye 669
358 johnpye 912 CONSOLE_DEBUG("Found %lu observation variables:",gl_length(sys->obslist));
359     for(i=1; i<=gl_length(sys->obslist); ++i){
360     info = (struct Integ_var_t *)gl_fetch(sys->obslist, i);
361     varname = var_make_name(sys->system,info->i);
362 johnpye 669 CONSOLE_DEBUG("observation[%d] = \"%s\"",i,varname);
363     ASC_FREE(varname);
364     }
365    
366 johnpye 912 /* CONSOLE_DEBUG("Checking found vars..."); */
367     if(gl_length(sys->dynvars)==0){
368 johnpye 669 ERROR_REPORTER_HERE(ASC_USER_ERROR,"No solver_var vars found to integrate (check 'ode_type'?).");
369     return 0;
370     }
371    
372 johnpye 912 CONSOLE_DEBUG("Found %lu vars.", gl_length(sys->dynvars));
373 johnpye 741
374 johnpye 669 maxderiv = 0;
375     numstates = 0;
376 johnpye 912 for(i=1; i<=gl_length(sys->dynvars); ++i){
377     info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
378 johnpye 903 if(info->type==1 || info->type==0)numstates++;
379 johnpye 669 if(maxderiv < info->type - 1)maxderiv = info->type - 1;
380 johnpye 912 varname = var_make_name(sys->system,info->i);
381 johnpye 903 /* CONSOLE_DEBUG("var[%d] = \"%s\": ode_index = %ld",i,varname,info->type); */
382 johnpye 669 ASC_FREE(varname);
383     }
384     if(maxderiv == 0){
385 johnpye 903 ERROR_REPORTER_HERE(ASC_USER_ERROR,"No derivatives found (check 'ode_type' values for your vars).");
386 johnpye 669 return 0;
387     }
388     if(numstates == 0){
389 johnpye 903 ERROR_REPORTER_HERE(ASC_USER_ERROR,"No states found (check 'odetype' values for your vars).");
390 johnpye 669 return 0;
391     }
392    
393    
394 johnpye 912 if(!integrator_check_indep_var(sys))return 0;
395 johnpye 669
396 johnpye 912 gl_sort(sys->dynvars,(CmpFunc)Integ_CmpDynVars);
397 johnpye 669
398 johnpye 912 /*
399     for(i=1; i<=gl_length(sys->dynvars); ++i){
400     info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
401     varname = var_make_name(sys->system,info->i);
402     // CONSOLE_DEBUG("var[%d] = \"%s\": ode_type = %ld",i,varname,info->type);
403 johnpye 669 ASC_FREE(varname);
404 johnpye 912 }*/
405 johnpye 669
406     /* link up derivative chains */
407 johnpye 741
408 johnpye 669 prev = NULL;
409 johnpye 912 for(i=1; i<=gl_length(sys->dynvars); ++i){ /* why does gl_list index with base 1??? */
410     info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
411 johnpye 669 info->derivative = NULL;
412    
413 johnpye 912 derivname = var_make_name(sys->system,info->i);
414 johnpye 669 if(prev!=NULL){
415 johnpye 912 varname = var_make_name(sys->system,prev->i);
416 johnpye 669 }else{
417     varname = NULL;
418     }
419 johnpye 903 /* CONSOLE_DEBUG("current = \"%s\", previous = \"%s\"",derivname,varname); */
420 johnpye 669 ASC_FREE(derivname);
421     if(varname)ASC_FREE(varname);
422    
423     if(info->type == INTEG_STATE_VAR || info->type == INTEG_ALGEBRAIC_VAR){
424 johnpye 912 varname = var_make_name(sys->system,info->i);
425 johnpye 903 /* CONSOLE_DEBUG("Var \"%s\" is not a derivative",varname); */
426 johnpye 669 ASC_FREE(varname);
427     info->derivative_of = NULL;
428     info->type = INTEG_STATE_VAR;
429     }else{
430     if(prev==NULL || info->index != prev->index){
431 johnpye 903 /* CONSOLE_DEBUG("current current type = %ld",info->type); */
432     /* if(prev!=NULL){
433 johnpye 709 CONSOLE_DEBUG("current index = %ld, previous = %ld",info->index,prev->index);
434 johnpye 669 }else{
435 johnpye 709 CONSOLE_DEBUG("current index = %ld, current type = %ld",info->index,info->type);
436 johnpye 903 } */
437 johnpye 912 varname = var_make_name(sys->system,info->i);
438 johnpye 669 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Derivative %d of \"%s\" is present without its un-differentiated equivalent"
439     , info->type-1
440     , varname
441     );
442     ASC_FREE(varname);
443     return 0;
444     }else if(info->type != prev->type + 1){
445 johnpye 912 derivname = var_make_name(sys->system,info->i);
446     varname = var_make_name(sys->system,prev->i);
447 johnpye 669 ERROR_REPORTER_HERE(ASC_USER_ERROR
448     ,"Looking at \"%s\", expected a derivative (order %d) of \"%s\"."
449     ,varname
450     ,prev->type+1
451     ,derivname
452     );
453     ASC_FREE(varname);
454     ASC_FREE(derivname);
455     return 0;
456     }else{
457 johnpye 912 varname = var_make_name(sys->system,prev->i);
458     derivname = var_make_name(sys->system,info->i);
459 johnpye 669 CONSOLE_DEBUG("Var \"%s\" is the derivative of \"%s\"",derivname,varname);
460     ASC_FREE(varname);
461     ASC_FREE(derivname);
462     info->derivative_of = prev;
463     numy++;
464     }
465 johnpye 741 }
466 johnpye 669 prev = info;
467     }
468    
469     /* record which vars have derivatives and which don't */
470 johnpye 912 for(i=1; i<=gl_length(sys->dynvars); ++i){
471     info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
472 johnpye 669 if(info->derivative_of){
473     info->derivative_of->derivative = info->i;
474 johnpye 912 /* varname = var_make_name(sys->system,info->derivative_of->i);
475     derivname = var_make_name(sys->system,info->derivative_of->derivative);
476 johnpye 903 CONSOLE_DEBUG("Var \"%s\" is the derivative of \"%s\"",derivname,varname);
477     ASC_FREE(varname);
478     ASC_FREE(derivname); */
479 johnpye 741 }
480 johnpye 669 }
481    
482     CONSOLE_DEBUG("Indentifying states...");
483    
484     /* count numy: either it's a state, or it has a higher-order derivative */
485     numy = 0;
486 johnpye 912 for(i=1; i<=gl_length(sys->dynvars); ++i){
487     info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
488 johnpye 669 if(info->type == INTEG_STATE_VAR || info->type == INTEG_ALGEBRAIC_VAR || info->derivative != NULL){
489 johnpye 912 varname = var_make_name(sys->system,info->i);
490 johnpye 903 if(var_fixed(info->i)){
491     CONSOLE_DEBUG("Var \"%s\" is a FIXED state variable",varname);
492     }else{
493     CONSOLE_DEBUG("Var \"%s\" is a state variable",varname);
494 johnpye 908 }
495 johnpye 669 ASC_FREE(varname);
496     info->isstate = 1;
497     numy++;
498     }else{
499 johnpye 912 varname = var_make_name(sys->system,info->i);
500 johnpye 903 CONSOLE_DEBUG("Var \"%s\" is a NON-STATE derivative",varname);
501     ASC_FREE(varname);
502 johnpye 669 info->isstate = 0;
503     }
504     }
505    
506     /*
507     create lists 'y' and 'ydot'. some elements of ydot don't correspond
508     to variables in our model: these are the algebraic vars.
509     */
510    
511     CONSOLE_DEBUG("Identified %d state variables", numy);
512    
513 johnpye 912 sys->y = ASC_NEW_ARRAY(struct var_variable *,numy);
514     sys->ydot = ASC_NEW_ARRAY(struct var_variable *,numy);
515 johnpye 669
516 johnpye 741 /*
517 johnpye 669 at this point we know there are no missing derivatives etc, so we
518     can use (i-1) as the index into y and ydot. any variable with
519     'derivative_of' set to null is a state variable... but it might already
520     be getting added
521     */
522 johnpye 912 for(j=0, i=1; i<=gl_length(sys->dynvars); ++i){
523     info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
524 johnpye 669 if(!info->isstate)continue;
525 johnpye 917 varname = var_make_name(sys->system,info->i);
526 johnpye 669 if(info->derivative == NULL){
527 johnpye 917 CONSOLE_DEBUG("Pure algebraic: %s",varname);
528 johnpye 912 sys->y[j] = info->i; /* a pure algebraic variable */
529     sys->ydot[j] = NULL;
530 johnpye 669 }else{
531 johnpye 917 CONSOLE_DEBUG("Differential variable: %s",varname);
532 johnpye 912 sys->y[j] = info->i; /* a variable whose derivative is present in the model */
533     sys->ydot[j] = info->derivative;
534 johnpye 669 }
535 johnpye 917 ASC_FREE(varname);
536     ++j;
537 johnpye 669 }
538    
539 johnpye 912 /*
540     set up the y_id table so that given a 'variable number' from relman_diff2,
541     we can work out where that fits in our y and ydot vectors.
542    
543     There's something a bit fishy about fetching the varlist at this late stage...
544     */
545    
546 johnpye 921 varlist = slv_get_solvers_var_list(sys->system);
547     nvarlist = slv_get_num_solvers_vars(sys->system);
548 johnpye 912
549 johnpye 921 CONSOLE_DEBUG("WORKING THROUGH THE SOLVER'S VAR LIST %d",nvarlist);
550 johnpye 912
551     sys->y_id = ASC_NEW_ARRAY(long, nvarlist);
552     for(i=0; i< nvarlist; ++i){
553     sys->y_id[i] = -2;
554     }
555    
556 johnpye 917 CONSOLE_DEBUG("WORKING THROUGH %ld DYNVARS",gl_length(sys->dynvars));
557 johnpye 912
558     for(i=1; i <= gl_length(sys->dynvars); ++i){
559     info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
560     sys->y_id[info->varindx] = i;
561 johnpye 917 varname = var_make_name(sys->system,info->i);
562     if(info->varindx > 0){
563     CONSOLE_DEBUG("Variable dynvars[%d] = y_id[%d] = '%s'",i,info->varindx,varname);
564 johnpye 912 }else{
565 johnpye 917 CONSOLE_DEBUG("VARIABLE dynvars[%d] = y_id[%d] = '%s' NOT FOUND IN VARLIST",i,info->varindx,varname);
566 johnpye 912 }
567 johnpye 917 ASC_FREE(varname);
568 johnpye 912 }
569     for(i=0; i< nvarlist; ++i){
570     if(sys->y_id[i] < 0){
571 johnpye 917 varname = var_make_name(sys->system,varlist[i]);
572 johnpye 921 CONSOLE_DEBUG("UNCONNECTED SOLVER VAR varlist[%d] = '%s' (%s)",i,varname,var_fixed(varlist[i])?"fixed":"not fixed");
573 johnpye 912 ASC_FREE(varname);
574     }
575     }
576    
577     nrels = slv_get_num_solvers_rels(sys->system);
578 johnpye 669 if(numy != nrels){
579     ERROR_REPORTER_HERE(ASC_USER_ERROR
580     ,"System is not square: solver has %d rels, found %d system states"
581     ,nrels, numy
582     );
583     return 0;
584     }
585    
586 johnpye 912 sys->n_y = numy;
587 johnpye 669
588 johnpye 912 if(!integrator_sort_obs_vars(sys))return 0;
589 johnpye 669
590     return 1;
591     }
592    
593 johnpye 912 void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitfn){
594 johnpye 669 struct var_variable **vlist;
595     int i, vlen;
596    
597 ben.allan 704 /* visit all the slv_system_t master var lists to collect vars */
598     /* find the vars mostly in this one */
599 johnpye 912 vlist = slv_get_master_var_list(sys->system);
600     vlen = slv_get_num_master_vars(sys->system);
601 ben.allan 704 for (i=0;i<vlen;i++) {
602 johnpye 912 (*visitfn)(sys, vlist[i], &i);
603 ben.allan 704 }
604    
605 johnpye 669 /*
606     CONSOLE_DEBUG("Checked %d vars",vlen);
607 johnpye 912 integrator_print_var_stats(sys);
608 johnpye 669 */
609    
610 ben.allan 704 /* probably nothing here, but gotta check. */
611 johnpye 912 vlist = slv_get_master_par_list(sys->system);
612     vlen = slv_get_num_master_pars(sys->system);
613 ben.allan 704 for (i=0;i<vlen;i++) {
614 johnpye 912 (*visitfn)(sys, vlist[i], NULL);
615 ben.allan 704 }
616 johnpye 669
617     /*
618 ben.allan 704 CONSOLE_DEBUG("Checked %d pars",vlen);
619 johnpye 912 integrator_print_var_stats(sys);
620 johnpye 669 */
621    
622 ben.allan 704 /* might find t here */
623 johnpye 912 vlist = slv_get_master_unattached_list(sys->system);
624     vlen = slv_get_num_master_unattached(sys->system);
625 ben.allan 704 for (i=0;i<vlen;i++) {
626 johnpye 912 (*visitfn)(sys, vlist[i], NULL);
627 ben.allan 704 }
628 johnpye 669
629     /* CONSOLE_DEBUG("Checked %d unattached",vlen); */
630     }
631     /**
632     @return 1 on success
633     */
634 johnpye 912 int integrator_analyse_ode(IntegratorSystem *sys){
635 ben.allan 704 struct Integ_var_t *v1,*v2;
636 johnpye 709 long half,i,len;
637 ben.allan 704 int happy=1;
638 johnpye 894 char *varname1, *varname2;
639 johnpye 669
640 ben.allan 704 CONSOLE_DEBUG("Starting ODE analysis");
641 johnpye 669 IntegInitSymbols();
642    
643 ben.allan 704 /* collect potential states and derivatives */
644 johnpye 912 sys->indepvars = gl_create(10L); /* t var info */
645     sys->dynvars = gl_create(200L); /* y ydot var info */
646     sys->obslist = gl_create(100L); /* obs info */
647     if (sys->dynvars == NULL
648     || sys->obslist == NULL
649     || sys->indepvars == NULL
650 ben.allan 704 ){
651     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
652     return 0;
653     }
654 johnpye 669
655 johnpye 912 sys->nstates = sys->nderivs = 0;
656 johnpye 669
657 johnpye 912 integrator_visit_system_vars(sys,&integrator_ode_classify_var);
658 johnpye 669
659 johnpye 912 integrator_print_var_stats(sys);
660 johnpye 669
661 johnpye 912 if(!integrator_check_indep_var(sys))return 0;
662 johnpye 669
663 ben.allan 704 /* check sanity of state and var lists */
664    
665 johnpye 912 len = gl_length(sys->dynvars);
666 johnpye 669 half = len/2;
667 johnpye 709 CONSOLE_DEBUG("NUMBER OF DYNAMIC VARIABLES = %ld",half);
668 ben.allan 704
669 johnpye 912 if (len % 2 || len == 0L || sys->nstates != sys->nderivs ) {
670 ben.allan 704 /* list length must be even for vars to pair off */
671     ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"n_y != n_ydot, or no dynamic vars found. Fix your indexing.");
672     return 0;
673     }
674 johnpye 912 gl_sort(sys->dynvars,(CmpFunc)Integ_CmpDynVars);
675     if (gl_fetch(sys->dynvars,len)==NULL) {
676 ben.allan 704 ERROR_REPORTER_NOLINE(ASC_PROG_ERR,"Mysterious NULL found!");
677     return 0;
678     }
679 johnpye 912 sys->states = gl_create(half); /* state vars Integ_var_t references */
680     sys->derivs = gl_create(half); /* derivative var atoms */
681 ben.allan 704 for (i=1;i < len; i+=2) {
682 johnpye 912 v1 = (struct Integ_var_t *)gl_fetch(sys->dynvars,i);
683     v2 = (struct Integ_var_t *)gl_fetch(sys->dynvars,i+1);
684 ben.allan 704 if (v1->type!=1 || v2 ->type !=2 || v1->index != v2->index) {
685 johnpye 912 varname1 = var_make_name(sys->system,v1->i);
686     varname2 = var_make_name(sys->system,v2->i);
687 johnpye 908
688 johnpye 894 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Mistyped or misindexed dynamic variables: %s (%s = %ld,%s = %ld) and %s (%s = %ld,%s = %ld).",
689     varname1, SCP(STATEFLAG),v1->type,SCP(STATEINDEX),v1->index,
690     varname2, SCP(STATEFLAG),v2->type,SCP(STATEINDEX),v2->index
691 ben.allan 704 );
692 johnpye 894 ASC_FREE(varname1);
693     ASC_FREE(varname2);
694 ben.allan 704 happy=0;
695     break;
696     } else {
697 johnpye 912 gl_append_ptr(sys->states,(POINTER)v1);
698     gl_append_ptr(sys->derivs,(POINTER)v2->i);
699 ben.allan 704 }
700     }
701     if (!happy) {
702     return 0;
703     }
704 johnpye 912 sys->n_y = half;
705     sys->y = ASC_NEW_ARRAY(struct var_variable *, half);
706     sys->y_id = ASC_NEW_ARRAY(long, half);
707     sys->ydot = ASC_NEW_ARRAY(struct var_variable *, half);
708     if (sys->y==NULL || sys->ydot==NULL || sys->y_id==NULL) {
709 ben.allan 704 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
710     return 0;
711     }
712     for (i = 1; i <= half; i++) {
713 johnpye 912 v1 = (struct Integ_var_t *)gl_fetch(sys->states,i);
714     sys->y[i-1] = v1->i;
715     sys->y_id[i-1] = v1->index;
716     sys->ydot[i-1] = (struct var_variable *)gl_fetch(sys->derivs,i);
717 ben.allan 704 }
718 johnpye 669
719 johnpye 912 if(!integrator_sort_obs_vars(sys))return 0;
720 ben.allan 704
721     /* don't need the gl_lists now that we have arrays for everyone */
722 johnpye 912 gl_destroy(sys->states);
723     gl_destroy(sys->derivs);
724     gl_free_and_destroy(sys->indepvars); /* we own the objects in indepvars */
725     gl_free_and_destroy(sys->dynvars); /* we own the objects in dynvars */
726     gl_free_and_destroy(sys->obslist); /* and obslist */
727     sys->states = NULL;
728     sys->derivs = NULL;
729     sys->indepvars = NULL;
730     sys->dynvars = NULL;
731     sys->obslist = NULL;
732 johnpye 669
733 ben.allan 704 /* analysis completed OK */
734     return 1;
735 johnpye 669 }
736    
737     /**
738     Reindex observations. Sort if the user mostly numbered. Take natural order
739 ben.allan 704 if user just booleaned.
740 johnpye 741
741 johnpye 669 @return 1 on success
742     */
743 johnpye 912 static int integrator_sort_obs_vars(IntegratorSystem *sys){
744 johnpye 777 int half, i, len = 0;
745 ben.allan 704 struct Integ_var_t *v2;
746 johnpye 669
747 johnpye 912 half = sys->n_y;
748     len = gl_length(sys->obslist);
749 ben.allan 704 /* we shouldn't be seeing NULL here ever except if malloc fail. */
750     if (len > 1L) {
751 johnpye 912 half = ((struct Integ_var_t *)gl_fetch(sys->obslist,1))->index;
752 ben.allan 704 /* half != 0 now because we didn't collect 0 indexed vars */
753     for (i=2; i <= len; i++) {
754 johnpye 912 if (half != ((struct Integ_var_t *)gl_fetch(sys->obslist,i))->index) {
755 ben.allan 704 /* change seen. sort and go on */
756 johnpye 912 gl_sort(sys->obslist,(CmpFunc)Integ_CmpObs);
757 ben.allan 704 break;
758     }
759     }
760 johnpye 669 }
761 ben.allan 704 for (i = half = 1; i <= len; i++) {
762 johnpye 912 v2 = (struct Integ_var_t *)gl_fetch(sys->obslist,i);
763 ben.allan 704 if (v2==NULL) {
764     /* we shouldn't be seeing NULL here ever except if malloc fail. */
765 johnpye 912 gl_delete(sys->obslist,i,0); /* should not be gl_delete(so,i,1) */
766 ben.allan 704 } else {
767     Integ_SetObsId(v2->i,half);
768     v2->index = half++;
769     }
770     }
771 johnpye 669
772 ben.allan 704 /* obslist now uniquely indexed, no nulls */
773     /* make into arrays */
774 johnpye 912 half = gl_length(sys->obslist);
775     sys->obs = ASC_NEW_ARRAY(struct var_variable *,half);
776     sys->obs_id = ASC_NEW_ARRAY(long, half);
777     if ( sys->obs==NULL || sys->obs_id==NULL) {
778 ben.allan 704 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
779     return 0;
780     }
781 johnpye 912 sys->n_obs = half;
782 ben.allan 704 for (i = 1; i <= half; i++) {
783 johnpye 912 v2 = (struct Integ_var_t *)gl_fetch(sys->obslist,i);
784     sys->obs[i-1] = v2->i;
785     sys->obs_id[i-1] = v2->index;
786 ben.allan 704 }
787    
788     return 1;
789 johnpye 669 }
790    
791 johnpye 912 static void integrator_print_var_stats(IntegratorSystem *sys){
792     int v = gl_length(sys->dynvars);
793     int i = gl_length(sys->indepvars);
794 johnpye 669 CONSOLE_DEBUG("Currently %d vars, %d indep",v,i);
795     }
796    
797     /**
798     @return 1 on success
799     */
800 johnpye 912 static int integrator_check_indep_var(IntegratorSystem *sys){
801 johnpye 669 int len, i;
802     struct Integ_var_t *info;
803     char *varname;
804    
805 ben.allan 704 /* check the sanity of the independent variable */
806 johnpye 912 len = gl_length(sys->indepvars);
807 ben.allan 704 if (!len) {
808     ERROR_REPORTER_HERE(ASC_PROG_ERR,"No independent variable found.");
809     return 0;
810     }
811 johnpye 669 if (len > 1) {
812 ben.allan 704 ERROR_REPORTER_START_HERE(ASC_USER_ERROR);
813     FPRINTF(ASCERR,"Excess %ld independent variables found:",
814     len);
815     for(i=1; i <=len;i++) {
816 johnpye 912 info = (struct Integ_var_t *)gl_fetch(sys->indepvars,i);
817 johnpye 669 if(info==NULL)continue;
818    
819 johnpye 912 varname = var_make_name(sys->system,info->i);
820 johnpye 894 FPRINTF(ASCERR," %s",varname);
821 ben.allan 704 ASC_FREE(varname);
822     }
823     FPRINTF(ASCERR , "\nSet the \"%s\" flag on all but one of these to %s >= 0.\n"
824 johnpye 669 , SCP(STATEFLAG),SCP(STATEFLAG)
825     );
826 ben.allan 704 error_reporter_end_flush();
827     return 0;
828     }else{
829 johnpye 912 info = (struct Integ_var_t *)gl_fetch(sys->indepvars,1);
830     sys->x = info->i;
831 johnpye 669 }
832     return 1;
833     }
834    
835     /*------------------------------------------------------------------------------
836     CLASSIFICATION OF VARIABLES (for ANALYSIS step)
837     */
838    
839 johnpye 912 #define INTEG_ADD_TO_LIST(info,TYPE,INDEX,VAR,VARINDX,LIST) \
840 johnpye 669 info = ASC_NEW(struct Integ_var_t); \
841     if(info==NULL){ \
842     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory (INTEG_VAR_NEW)"); \
843     return; \
844     } \
845     info->type=TYPE; \
846     info->index=INDEX; \
847     info->i=VAR; \
848 johnpye 912 if(VARINDX==NULL){ \
849     info->varindx = -1; \
850     }else{ \
851     info->varindx = *VARINDX; \
852     } \
853 johnpye 669 gl_append_ptr(LIST,(void *)info); \
854     info = NULL
855    
856     /**
857     In a DAE, it's either the (single) independent variable, or it's a
858     variable in the model.
859    
860     I'm not sure what we should be doing with variables that are already
861     present as derivatives of other variables, I guess those ones need to be
862     removed from the list in a second pass?
863     */
864 johnpye 912 void integrator_dae_classify_var(IntegratorSystem *sys
865     , struct var_variable *var, const int *varindx
866     ){
867 johnpye 669 struct Integ_var_t *info;
868 ben.allan 704 long type,index;
869 johnpye 669
870     /* filter for recognition of solver_vars */
871 ben.allan 704 var_filter_t vfilt;
872     vfilt.matchbits = VAR_SVAR;
873     vfilt.matchvalue = VAR_SVAR;
874 johnpye 669
875 ben.allan 704 assert(var != NULL && var_instance(var)!=NULL );
876 johnpye 669
877     if( var_apply_filter(var,&vfilt) ) {
878 johnpye 921 if(!var_active(var)){
879     CONSOLE_DEBUG("VARIABLE IS NOT ACTIVE");
880     return;
881     }
882 johnpye 669
883     if(!var_fixed(var)){
884 ben.allan 704 /* get the ode_type and ode_id of this solver_var */
885 johnpye 669 type = DynamicVarInfo(var,&index);
886    
887     if(type==INTEG_OTHER_VAR){
888     /* if the var's type is -1, it's independent */
889 johnpye 912 INTEG_ADD_TO_LIST(info,INTEG_OTHER_VAR,0,var,varindx,sys->indepvars);
890 johnpye 669 }else{
891     if(type < 0)type=0;
892     /* any other type of var is in the DAE system, at least for now */
893 johnpye 912 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
894 johnpye 669 }
895 johnpye 921 }
896     #if 1
897     else{
898 johnpye 903 /* fixed variable, only include it if ode_type == 1 */
899     type = DynamicVarInfo(var,&index);
900     if(type==INTEG_STATE_VAR){
901 johnpye 912 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
902 johnpye 903 }
903 johnpye 669 }
904 johnpye 921 #endif
905 johnpye 669
906 ben.allan 704 /* if the var's obs_id > 0, add it to the observation list */
907 johnpye 669 if(ObservationVar(var,&index) != NULL && index > 0L) {
908 johnpye 912 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->obslist);
909 johnpye 669 }
910 ben.allan 704 }
911 johnpye 669 }
912    
913     /**
914     Inspect a specific variable and work out what type it is (what 'ode_type' it
915     has) and what other variable(s) it corresponds to (ie dydt corresponds to
916     y as a derivative).
917    
918     @TODO add ability to create new variables for 'missing' derivative vars?
919 ben.allan 704 */
920 johnpye 912 void integrator_ode_classify_var(IntegratorSystem *sys, struct var_variable *var
921     , const int *varindx
922     ){
923 ben.allan 704 struct Integ_var_t *info;
924     long type,index;
925 johnpye 669
926 ben.allan 704 var_filter_t vfilt;
927     vfilt.matchbits = VAR_SVAR;
928     vfilt.matchvalue = VAR_SVAR;
929 johnpye 669
930 ben.allan 704 assert(var != NULL && var_instance(var)!=NULL );
931    
932 johnpye 669 if( var_apply_filter(var,&vfilt) ) {
933 ben.allan 704 /* it's a solver var: what type of variable? */
934 johnpye 669 type = DynamicVarInfo(var,&index);
935 ben.allan 704
936 johnpye 669 if(type==INTEG_ALGEBRAIC_VAR){
937     /* no action required */
938     }else if(type==INTEG_OTHER_VAR){
939     /* i.e. independent var */
940 johnpye 912 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->indepvars);
941 johnpye 669 }else if(type>=INTEG_STATE_VAR){
942 johnpye 912 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
943 ben.allan 704 if(type == 1){
944 johnpye 912 sys->nstates++;
945 ben.allan 704 }else if(type == 2){ /* what about higher-order derivatives? -- JP */
946 johnpye 912 sys->nderivs++;
947 johnpye 669 }else{
948     ERROR_REPORTER_HERE(ASC_USER_WARNING,"Higher-order (>=2) derivatives are not supported in ODEs.");
949 johnpye 741 } }
950 johnpye 669
951     if(ObservationVar(var,&index) != NULL && index > 0L) {
952 johnpye 912 INTEG_ADD_TO_LIST(info,0L,index,var,varindx,sys->obslist);
953 ben.allan 704 }
954     }
955 johnpye 669 }
956    
957     /**
958     Look at a variable and determine if it's the independent variable or not.
959 johnpye 741 This is just for the purpose of the integrator_find_indep_var function,
960 johnpye 669 which is a utility function provided for use by the GUI.
961     */
962 johnpye 912 void integrator_classify_indep_var(IntegratorSystem *sys
963     , struct var_variable *var, const int *varindx
964     ){
965 ben.allan 704 struct Integ_var_t *info;
966     long type,index;
967 johnpye 669
968 ben.allan 704 var_filter_t vfilt;
969     vfilt.matchbits = VAR_SVAR;
970 johnpye 669 vfilt.matchvalue = VAR_SVAR;
971    
972 johnpye 725 /* CONSOLE_DEBUG("..."); */
973 johnpye 669
974 ben.allan 704 assert(var != NULL && var_instance(var)!=NULL );
975    
976 johnpye 669 if( var_apply_filter(var,&vfilt) ) {
977     type = DynamicVarInfo(var,&index);
978 ben.allan 704
979 johnpye 669 if(type==INTEG_OTHER_VAR){
980     /* i.e. independent var */
981 johnpye 912 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->indepvars);
982 johnpye 669 }
983     }
984 ben.allan 704 }
985    
986 johnpye 669 /**
987     Look at a variable, and if it is an 'ODE variable' (it has a child instance
988     named 'ode_type') return its type, which will be either:
989     - INTEG_OTHER_VAR (if 'ode_type' is -1)
990     - INTEG_ALGEBRAIC_VAR (if 'ode_type' is zero or any negative value < -1)
991     - INTEG_STATE_VAR (if 'ode_type' is 1)
992     - values 2, 3 or up, indicating derivatives (1st deriv=2, 2nd deriv=3, etc)
993    
994     If the parameter 'index' is not null, the value of 'ode_id' will be stuffed
995     there.
996 ben.allan 704 */
997     static long DynamicVarInfo(struct var_variable *v,long *index){
998     struct Instance *c, *d, *i;
999     i = var_instance(v);
1000 johnpye 669 assert(i!=NULL);
1001     assert(STATEFLAG!=NULL);
1002 ben.allan 704 assert(STATEINDEX!=NULL);
1003     c = ChildByChar(i,STATEFLAG);
1004     d = ChildByChar(i,STATEINDEX);
1005     /* lazy evaluation is important in the following if */
1006     if(c == NULL
1007     || d == NULL
1008     || InstanceKind(c) != INTEGER_INST
1009     || InstanceKind(d) != INTEGER_INST
1010     || !AtomAssigned(c)
1011 johnpye 669 || (!AtomAssigned(d) && GetIntegerAtomValue(c) != INTEG_OTHER_VAR)
1012 ben.allan 704 ){
1013     return INTEG_ALGEBRAIC_VAR;
1014     }
1015     if (index != NULL) {
1016     *index = GetIntegerAtomValue(d);
1017     }
1018     return GetIntegerAtomValue(c);
1019     }
1020 johnpye 669
1021     /**
1022     Looks at the given variable checks if it is an 'observation variable'. This
1023     means that it has its 'obs_id' child instance set to a non-zero value.
1024    
1025     If the variable is an observation variable, its index value ('obs_id') is
1026     stuff into *index (provided index!=NULL), and the pointer to the original
1027     instance is rtruend.
1028    
1029 ben.allan 704 If it's not an observation variable, we return NULL and *index is untouched.
1030     */
1031     static struct var_variable *ObservationVar(struct var_variable *v, long *index){
1032     struct Instance *c,*i;
1033     i = var_instance(v);
1034     assert(i!=NULL);
1035     c = ChildByChar(i,OBSINDEX);
1036     if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
1037     return NULL;
1038     }
1039     if (index != NULL) {
1040     *index = GetIntegerAtomValue(c);
1041     }
1042     return v;
1043     }
1044 johnpye 669
1045     /*------------------------------------------------------------------------------
1046     RUNNING THE SOLVER
1047     */
1048    
1049     /*
1050     Make the call to the actual integrator we've selected, for the range of
1051 johnpye 912 time values specified. The sys contains all the specifics.
1052 johnpye 709
1053     Return 1 on success
1054 ben.allan 704 */
1055 johnpye 912 int integrator_solve(IntegratorSystem *sys, long i0, long i1){
1056 johnpye 669
1057 ben.allan 704 long nstep;
1058 johnpye 741 unsigned long start_index=0, finish_index=0;
1059 johnpye 912 assert(sys!=NULL);
1060 johnpye 669
1061 johnpye 912 nstep = integrator_getnsamples(sys)-1;
1062 ben.allan 704 /* check for at least 2 steps and dimensionality of x vs steps here */
1063 johnpye 669
1064     if (i0<0 || i1 <0) {
1065 ben.allan 704 /* dude, there's no way we're writing interactive stuff here... */
1066 johnpye 669 ERROR_REPORTER_HERE(ASC_PROG_ERROR,"Console input of integration limits has been disabled!");
1067 ben.allan 704 return 0;
1068 johnpye 908 } else {
1069 ben.allan 704 start_index=i0;
1070     finish_index =i1;
1071     if (start_index >= (unsigned long)nstep) {
1072 johnpye 669 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Start point (=%lu) must be an index in the range [0,%li]."
1073     ,start_index,nstep
1074 ben.allan 704 );
1075     return 0;
1076     }
1077     if (finish_index > (unsigned long)nstep) {
1078 johnpye 669 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"End point (=%lu) must be an index in the range [0,%li]."
1079     ,finish_index,nstep
1080 ben.allan 704 );
1081     return 0;
1082     }
1083     }
1084 johnpye 669
1085 ben.allan 704 if(finish_index <= start_index) {
1086 johnpye 669 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"End point comes before start point! (start=%lu, end=%lu)"
1087     ,start_index,finish_index
1088 ben.allan 704 );
1089     return 0;
1090 johnpye 669 }
1091    
1092 johnpye 888 CONSOLE_DEBUG("RUNNING INTEGRATION...");
1093    
1094 ben.allan 704 /* now go and run the integrator */
1095 johnpye 912 switch (sys->engine) {
1096 johnpye 908 case INTEG_LSODE:
1097 johnpye 912 return integrator_lsode_solve(sys, start_index, finish_index);
1098 johnpye 908 break;
1099 johnpye 669 #ifdef ASC_WITH_IDA
1100 johnpye 908 case INTEG_IDA:
1101 johnpye 912 return integrator_ida_solve(sys,start_index, finish_index);
1102 johnpye 908 break;
1103 ben.allan 704 #endif
1104 johnpye 709 default:
1105     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown integrator (invalid, or not implemented yet)");
1106     return 0;
1107 ben.allan 704 }
1108 johnpye 669 }
1109 ben.allan 704
1110 johnpye 669 /*---------------------------------------------------------------
1111     HANDLING THE LIST OF TIMESTEMPS
1112     */
1113    
1114     #define GETTER_AND_SETTER(TYPE,NAME) \
1115 johnpye 912 void integrator_set_##NAME(IntegratorSystem *sys, TYPE val){ \
1116     sys->NAME=val; \
1117 johnpye 669 } \
1118 johnpye 912 TYPE integrator_get_##NAME(IntegratorSystem *sys){ \
1119     return sys->NAME; \
1120 johnpye 669 }
1121    
1122 johnpye 908 GETTER_AND_SETTER(SampleList *,samples) /*;*/
1123     GETTER_AND_SETTER(double,maxstep) /*;*/
1124     GETTER_AND_SETTER(double,minstep) /*;*/
1125     GETTER_AND_SETTER(double,stepzero) /*;*/
1126     GETTER_AND_SETTER(int,maxsubsteps) /*;*/
1127 ben.allan 704 #undef GETTER_AND_SETTER
1128    
1129 johnpye 912 long integrator_getnsamples(IntegratorSystem *sys){
1130     assert(sys!=NULL);
1131     assert(sys->samples!=NULL);
1132     return samplelist_length(sys->samples);
1133 ben.allan 704 }
1134    
1135 johnpye 912 double integrator_getsample(IntegratorSystem *sys, long i){
1136     assert(sys!=NULL);
1137     assert(sys->samples!=NULL);
1138     return samplelist_get(sys->samples,i);
1139 ben.allan 704 }
1140    
1141 johnpye 912 void integrator_setsample(IntegratorSystem *sys, long i,double xi){
1142     assert(sys!=NULL);
1143     assert(sys->samples!=NULL);
1144     samplelist_set(sys->samples,i,xi);
1145 ben.allan 704 }
1146    
1147 johnpye 912 const dim_type *integrator_getsampledim(IntegratorSystem *sys){
1148     assert(sys!=NULL);
1149     assert(sys->samples!=NULL);
1150     return samplelist_dim(sys->samples);
1151 ben.allan 704 }
1152 johnpye 669
1153 johnpye 912 ASC_DLLSPEC(long) integrator_getcurrentstep(IntegratorSystem *sys){
1154     return sys->currentstep;
1155 johnpye 669 }
1156 ben.allan 704
1157 johnpye 669 /*------------------------------------------------------------------------------
1158 ben.allan 704 GET/SET VALUE OF THE INDEP VARIABLE
1159 johnpye 669 */
1160    
1161 johnpye 741 /**
1162 johnpye 669 Retrieve the value of the independent variable (time) from ASCEND
1163     and return it as a double.
1164 johnpye 741 */
1165 johnpye 912 double integrator_get_t(IntegratorSystem *sys){
1166     assert(sys->x!=NULL);
1167     return var_value(sys->x);
1168 ben.allan 704 }
1169 johnpye 669
1170     /**
1171     Set the value of the independent variable (time) in ASCEND.
1172 ben.allan 704 */
1173 johnpye 912 void integrator_set_t(IntegratorSystem *sys, double value){
1174     var_set_value(sys->x, value);
1175 johnpye 815 /* CONSOLE_DEBUG("set_t = %g", value); */
1176 ben.allan 704 }
1177    
1178 johnpye 669 /*------------------------------------------------------------------------------
1179     PASSING DIFFERENTIAL VARIABLES AND THEIR DERIVATIVES TO/FROM THE SOLVER
1180     */
1181     /**
1182     Retrieve the current values of the derivatives of the y-variables
1183     and stick them in the/an array that the integrator will use.
1184    
1185     If the pointer 'y' is NULL, the necessary space is allocated (and
1186     must be freed somewhere else).
1187 ben.allan 704 */
1188 johnpye 912 double *integrator_get_y(IntegratorSystem *sys, double *y) {
1189 ben.allan 704 long i;
1190    
1191     if (y==NULL) {
1192 johnpye 912 y = ASC_NEW_ARRAY_CLEAR(double, sys->n_y+1);
1193 ben.allan 704 /* C y[0] <==> ascend d.y[1] <==> f77 y(1) */
1194     }
1195    
1196 johnpye 912 for (i=0; i< sys->n_y; i++) {
1197     assert(sys->y[i]!=NULL);
1198     y[i] = var_value(sys->y[i]);
1199 johnpye 815 /* CONSOLE_DEBUG("ASCEND --> y[%ld] = %g", i+1, y[i]); */
1200 ben.allan 704 }
1201     return y;
1202     }
1203 johnpye 669
1204     /**
1205     Take the values of the derivatives from the array that the integrator
1206     uses, and use them to update the values of the corresponding variables
1207     in ASCEND.
1208 ben.allan 704 */
1209 johnpye 912 void integrator_set_y(IntegratorSystem *sys, double *y) {
1210 johnpye 669 long i;
1211 ben.allan 704 #ifndef NDEBUG
1212 johnpye 669 char *varname;
1213 ben.allan 704 #endif
1214    
1215 johnpye 912 for (i=0; i < sys->n_y; i++) {
1216     assert(sys->y[i]!=NULL);
1217     var_set_value(sys->y[i],y[i]);
1218 johnpye 669 #ifndef NDEBUG
1219 johnpye 912 varname = var_make_name(sys->system, sys->y[i]);
1220 johnpye 815 /* CONSOLE_DEBUG("y[%ld] = \"%s\" = %g --> ASCEND", i+1, varname, y[i]); */
1221 ben.allan 704 ASC_FREE(varname);
1222 johnpye 669 #endif
1223 ben.allan 704 }
1224     }
1225 johnpye 669
1226     /**
1227     Send the values of the derivatives of the 'y' variables to the solver.
1228     Allocate space for an array if necessary.
1229    
1230 johnpye 912 Any element in sys->ydot that is NULL will be passed over (the value
1231 johnpye 669 won't be modified in dydx).
1232 ben.allan 704 */
1233 johnpye 912 double *integrator_get_ydot(IntegratorSystem *sys, double *dydx) {
1234 ben.allan 704 long i;
1235    
1236     if (dydx==NULL) {
1237 johnpye 912 dydx = ASC_NEW_ARRAY_CLEAR(double, sys->n_y+1);
1238 ben.allan 704 /* C dydx[0] <==> ascend d.dydx[1] <==> f77 ydot(1) */
1239     }
1240    
1241 johnpye 912 for (i=0; i < sys->n_y; i++) {
1242     if(sys->ydot[i]!=NULL){
1243     dydx[i] = var_value(sys->ydot[i]);
1244 ben.allan 704 }
1245 johnpye 815 /* CONSOLE_DEBUG("ASCEND --> ydot[%ld] = %g", i+1, dydx[i]); */
1246 ben.allan 704 }
1247     return dydx;
1248     }
1249    
1250 johnpye 912 void integrator_set_ydot(IntegratorSystem *sys, double *dydx) {
1251 ben.allan 704 long i;
1252     #ifndef NDEBUG
1253 johnpye 908 /* char *varname; */
1254 ben.allan 704 #endif
1255 johnpye 912 for (i=0; i < sys->n_y; i++) {
1256     if(sys->ydot[i]!=NULL){
1257     var_set_value(sys->ydot[i],dydx[i]);
1258 johnpye 669 #ifndef NDEBUG
1259 johnpye 912 /* varname = var_make_name(sys->system, sys->ydot[i]);
1260 johnpye 709 CONSOLE_DEBUG("ydot[%ld] = \"%s\" = %g --> ASCEND", i+1, varname, dydx[i]);
1261 johnpye 903 ASC_FREE(varname); */
1262 johnpye 669 #endif
1263     }
1264     #ifndef NDEBUG
1265 johnpye 903 /*else{
1266 johnpye 709 CONSOLE_DEBUG("ydot[%ld] = %g (internal)", i+1, dydx[i]);
1267 johnpye 903 }*/
1268 johnpye 669 #endif
1269 ben.allan 704 }
1270     }
1271    
1272 johnpye 669 /*-------------------------------------------------------------
1273 ben.allan 704 RETRIEVING OBSERVATION DATA
1274 johnpye 908 */
1275 ben.allan 704
1276 johnpye 908 /**
1277 ben.allan 704 This function takes the inst in the solver and returns the vector of
1278     observation variables that are located in the submodel d.obs array.
1279     */
1280 johnpye 912 double *integrator_get_observations(IntegratorSystem *sys, double *obsi) {
1281 ben.allan 704 long i;
1282    
1283     if (obsi==NULL) {
1284 johnpye 912 obsi = ASC_NEW_ARRAY_CLEAR(double, sys->n_obs+1);
1285 ben.allan 704 }
1286    
1287     /* C obsi[0] <==> ascend d.obs[1] */
1288    
1289 johnpye 912 for (i=0; i < sys->n_obs; i++) {
1290     obsi[i] = var_value(sys->obs[i]);
1291 ben.allan 704 /* CONSOLE_DEBUG("*get_d_obs[%ld] = %g\n", i+1, obsi[i]); */
1292     }
1293     return obsi;
1294 johnpye 669 }
1295    
1296 johnpye 912 struct var_variable *integrator_get_observed_var(IntegratorSystem *sys, const long i){
1297 johnpye 669 assert(i>=0);
1298 johnpye 912 assert(i<sys->n_obs);
1299     return sys->obs[i];
1300 ben.allan 704 }
1301    
1302 johnpye 854 /**
1303     @NOTE Although this shouldn't be required for implementation of solver
1304     engines, this is useful for GUI reporting of integration results.
1305     */
1306 johnpye 912 struct var_variable *integrator_get_independent_var(IntegratorSystem *sys){
1307     return sys->x;
1308 johnpye 854 }
1309    
1310    
1311 johnpye 669 /*----------------------------------------------------
1312 johnpye 741 Build an analytic jacobian for solving the state system
1313 johnpye 669
1314 ben.allan 704 This necessarily ugly piece of code attempts to create a unique
1315     list of relations that explicitly contain the variables in the
1316     given input list. The utility of this information is that we know
1317     exactly which relations must be differentiated, to fill in the
1318     df/dy matrix. If the problem has very few derivative terms, this will
1319     be of great savings. If the problem arose from the discretization of
1320     a pde, then this will be not so useful. The decision wether to use
1321     this function or to simply differentiate the entire relations list
1322 johnpye 741 must be done before calling this function.
1323 johnpye 669
1324 ben.allan 704 Final Note: the callee owns the array, but not the array elements.
1325     */
1326     #define AVG_NUM_INCIDENT 4
1327 johnpye 669
1328 ben.allan 704
1329 johnpye 669 /**
1330     This function helps to arrange the observation variables in a sensible order.
1331     The 'obs_id' child instance of v, if present, is assigned the value of the
1332     given parameter 'index'.
1333 ben.allan 704 */
1334     static void Integ_SetObsId(struct var_variable *v, long index){
1335     struct Instance *c, *i;
1336     i = var_instance(v);
1337     assert(i!=NULL);
1338     c = ChildByChar(i,OBSINDEX);
1339     if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
1340     return;
1341     }
1342     SetIntegerAtomValue(c,index,0);
1343     }
1344    
1345 johnpye 669 /**
1346 johnpye 741 Compares observation structs. NULLs should end up at far end.
1347 ben.allan 704 */
1348     static int Integ_CmpObs(struct Integ_var_t *v1, struct Integ_var_t *v2){
1349     if(v1 == NULL)return 1;
1350     if(v2 == NULL)return -1;
1351     if(v1->index > v2->index)return 1;
1352     if(v1->index == v2->index)return 0;
1353     return -1;
1354     }
1355    
1356     /**
1357     Compares dynamic vars structs. NULLs should end up at far end.
1358     List should be sorted primarily by index and then by type, in order
1359     of increasing value of both.
1360     */
1361     static int Integ_CmpDynVars(struct Integ_var_t *v1, struct Integ_var_t *v2){
1362     if(v1 == NULL)return 1;
1363     if(v2 == NULL)return -1;
1364     if(v1->index > v2->index)return 1;
1365     if(v1->index != v2->index)return -1;
1366     if(v1->type > v2->type)return 1;
1367     return -1;
1368     }
1369 johnpye 669 /*----------------------------
1370     Output handling to the GUI/interface.
1371     */
1372    
1373 johnpye 912 int integrator_set_reporter(IntegratorSystem *sys
1374 johnpye 669 , IntegratorReporter *reporter
1375     ){
1376 johnpye 912 assert(sys!=NULL);
1377     sys->reporter = reporter;
1378 johnpye 854 /* ERROR_REPORTER_HERE(ASC_PROG_NOTE,"INTEGRATOR REPORTER HOOKS HAVE BEEN SET\n"); */
1379 johnpye 709 return 1;
1380 johnpye 669 }
1381    
1382 johnpye 912 int integrator_output_init(IntegratorSystem *sys){
1383     assert(sys!=NULL);
1384     assert(sys->reporter!=NULL);
1385     if(sys->reporter->init!=NULL){
1386 johnpye 669 /* call the specified output function */
1387 johnpye 912 return (*(sys->reporter->init))(sys);
1388 johnpye 669 }
1389     ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter init method");
1390     return 1;
1391     }
1392    
1393 johnpye 912 int integrator_output_write(IntegratorSystem *sys){
1394 johnpye 669 static int reported_already=0;
1395 johnpye 912 assert(sys!=NULL);
1396     if(sys->reporter->write!=NULL){
1397     return (*(sys->reporter->write))(sys);
1398 johnpye 669 }
1399     if(!reported_already){
1400     ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter write method (this message only shown once)");
1401     reported_already=1;
1402     }
1403     return 1;
1404     }
1405    
1406 johnpye 912 int integrator_output_write_obs(IntegratorSystem *sys){
1407 johnpye 669 static int reported_already=0;
1408 johnpye 912 assert(sys!=NULL);
1409     if(sys->reporter->write_obs!=NULL){
1410     return (*(sys->reporter->write_obs))(sys);
1411 johnpye 669 }
1412     if(!reported_already){
1413     ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter write_obs method (this message only shown once)");
1414     reported_already=1;
1415     }
1416     return 1;
1417 johnpye 741 }
1418 johnpye 669
1419 johnpye 912 int integrator_output_close(IntegratorSystem *sys){
1420     assert(sys!=NULL);
1421     if(sys->reporter->close!=NULL){
1422     return (*(sys->reporter->close))(sys);
1423 johnpye 669 }
1424     ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter close method");
1425     return 1;
1426     }
1427    
1428     /**
1429     Decode status codes from the integrator, and output them via FPRINTF.
1430 ben.allan 704 */
1431     int integrator_checkstatus(slv_status_t status) {
1432     if (status.converged) {
1433     return 1;
1434     }
1435     if (status.diverged) {
1436 johnpye 903 FPRINTF(stderr, "The derivative system did not converge. Integration will terminate.");
1437 ben.allan 704 return 0;
1438     }
1439     if (status.inconsistent) {
1440 johnpye 903 FPRINTF(stderr, "A numerically inconsistent state was discovered while "
1441     "calculating derivatives. Integration will terminate.");
1442 ben.allan 704 return 0;
1443     }
1444     if (status.time_limit_exceeded) {
1445 johnpye 854 FPRINTF(stderr, "The time limit was exceeded while calculating "
1446 johnpye 903 "derivatives. Integration will terminate.");
1447 ben.allan 704 return 0;
1448     }
1449     if (status.iteration_limit_exceeded) {
1450 johnpye 854 FPRINTF(stderr, "The iteration limit was exceeded while calculating "
1451 johnpye 903 "derivatives. Integration will terminate.");
1452 ben.allan 704 return 0;
1453     }
1454     if (status.panic) {
1455 johnpye 854 FPRINTF(stderr, "The user patience limit was exceeded while "
1456 johnpye 903 "calculating derivatives. Integration will terminate.");
1457 ben.allan 704 return 0;
1458     }
1459     return 0;
1460 johnpye 669 }

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