/[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 928 - (hide annotations) (download) (as text)
Wed Nov 22 10:32:18 2006 UTC (15 years, 8 months ago) by johnpye
File MIME type: text/x-csrc
File size: 40845 byte(s)
Commented out some stream redirection stuff for simplicity.
The CUnit test suite now works as expected (but without output suppression, for the moment).
Some more effort on IDA (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 928 struct Integ_var_t *derivative;
87     struct Integ_var_t *derivative_of;
88 johnpye 669 struct var_variable *i;
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 johnpye 928 int yindex;
338 johnpye 669 int maxderiv;
339    
340     CONSOLE_DEBUG("Starting DAE analysis");
341 ben.allan 704 IntegInitSymbols();
342 johnpye 669
343 johnpye 912 assert(sys->indepvars==NULL);
344 johnpye 669
345 johnpye 912 sys->indepvars = gl_create(10L); /* t var info */
346     sys->dynvars = gl_create(200L); /* y and ydot var info */
347     sys->obslist = gl_create(100L); /* obs info */
348 johnpye 669
349 johnpye 912 if(sys->indepvars==NULL
350     || sys->dynvars==NULL
351     || sys->obslist==NULL
352 johnpye 669 ){
353     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory");
354     return 0;
355     }
356    
357 johnpye 912 integrator_visit_system_vars(sys,&integrator_dae_classify_var);
358 johnpye 669
359 johnpye 912 CONSOLE_DEBUG("Found %lu observation variables:",gl_length(sys->obslist));
360     for(i=1; i<=gl_length(sys->obslist); ++i){
361     info = (struct Integ_var_t *)gl_fetch(sys->obslist, i);
362     varname = var_make_name(sys->system,info->i);
363 johnpye 669 CONSOLE_DEBUG("observation[%d] = \"%s\"",i,varname);
364     ASC_FREE(varname);
365     }
366    
367 johnpye 912 /* CONSOLE_DEBUG("Checking found vars..."); */
368     if(gl_length(sys->dynvars)==0){
369 johnpye 669 ERROR_REPORTER_HERE(ASC_USER_ERROR,"No solver_var vars found to integrate (check 'ode_type'?).");
370     return 0;
371     }
372    
373 johnpye 912 CONSOLE_DEBUG("Found %lu vars.", gl_length(sys->dynvars));
374 johnpye 741
375 johnpye 669 maxderiv = 0;
376     numstates = 0;
377 johnpye 912 for(i=1; i<=gl_length(sys->dynvars); ++i){
378     info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
379 johnpye 903 if(info->type==1 || info->type==0)numstates++;
380 johnpye 669 if(maxderiv < info->type - 1)maxderiv = info->type - 1;
381 johnpye 928 /* varname = var_make_name(sys->system,info->i);
382     CONSOLE_DEBUG("var[%d] = \"%s\": ode_index = %ld",i,varname,info->type);
383     ASC_FREE(varname); */
384 johnpye 669 }
385     if(maxderiv == 0){
386 johnpye 903 ERROR_REPORTER_HERE(ASC_USER_ERROR,"No derivatives found (check 'ode_type' values for your vars).");
387 johnpye 669 return 0;
388     }
389 johnpye 928 if(maxderiv > 1){
390     ERROR_REPORTER_HERE(ASC_USER_ERROR,"Higher-order derivatives found. You must provide a reduced order formulation for your model.");
391 johnpye 669 return 0;
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 928 fprintf(stderr,"\n\n\nSORTED VARS\n");
399 johnpye 912 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 johnpye 928 CONSOLE_DEBUG("var[%d] = \"%s\": ode_type = %ld",i,varname,info->type);
403 johnpye 669 ASC_FREE(varname);
404 johnpye 928 }
405 johnpye 669
406 johnpye 928 /* link up variables with their derivatives */
407 johnpye 669 prev = NULL;
408 johnpye 912 for(i=1; i<=gl_length(sys->dynvars); ++i){ /* why does gl_list index with base 1??? */
409     info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
410 johnpye 928
411 johnpye 669 if(info->type == INTEG_STATE_VAR || info->type == INTEG_ALGEBRAIC_VAR){
412 johnpye 912 varname = var_make_name(sys->system,info->i);
413 johnpye 928 CONSOLE_DEBUG("Var \"%s\" is an algebraic variable",varname);
414 johnpye 669 ASC_FREE(varname);
415 johnpye 928 info->type = INTEG_STATE_VAR;
416 johnpye 669 info->derivative_of = NULL;
417     }else{
418     if(prev==NULL || info->index != prev->index){
419 johnpye 928 /* derivative, but without undifferentiated var present in model */
420 johnpye 912 varname = var_make_name(sys->system,info->i);
421 johnpye 669 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Derivative %d of \"%s\" is present without its un-differentiated equivalent"
422     , info->type-1
423     , varname
424     );
425     ASC_FREE(varname);
426     return 0;
427     }else if(info->type != prev->type + 1){
428 johnpye 928 /* derivative, but missing the next-lower-order derivative */
429 johnpye 912 derivname = var_make_name(sys->system,info->i);
430     varname = var_make_name(sys->system,prev->i);
431 johnpye 669 ERROR_REPORTER_HERE(ASC_USER_ERROR
432     ,"Looking at \"%s\", expected a derivative (order %d) of \"%s\"."
433     ,varname
434     ,prev->type+1
435     ,derivname
436     );
437     ASC_FREE(varname);
438     ASC_FREE(derivname);
439     return 0;
440     }else{
441 johnpye 928 /* variable with derivative */
442 johnpye 912 varname = var_make_name(sys->system,prev->i);
443     derivname = var_make_name(sys->system,info->i);
444 johnpye 669 CONSOLE_DEBUG("Var \"%s\" is the derivative of \"%s\"",derivname,varname);
445     ASC_FREE(varname);
446     ASC_FREE(derivname);
447     info->derivative_of = prev;
448     }
449 johnpye 741 }
450 johnpye 669 prev = info;
451     }
452    
453 johnpye 928 /* record which vars have derivatives and which don't, and count 'states' */
454     numy = 0;
455 johnpye 912 for(i=1; i<=gl_length(sys->dynvars); ++i){
456     info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
457 johnpye 669 if(info->derivative_of){
458 johnpye 928 info->derivative_of->derivative = info;
459     }else{
460 johnpye 669 numy++;
461     }
462     }
463    
464 johnpye 928 /* allocate storage for the 'y' and 'ydot' arrays */
465 johnpye 912 sys->y = ASC_NEW_ARRAY(struct var_variable *,numy);
466     sys->ydot = ASC_NEW_ARRAY(struct var_variable *,numy);
467 johnpye 928 sys->y_id = ASC_NEW_ARRAY(int, slv_get_num_master_vars(sys->system));
468 johnpye 669
469 johnpye 928 /* now add variables and their derivatives to 'ydot' and 'y' */
470     yindex = 0;
471    
472     for(i=1; i<=gl_length(sys->dynvars); ++i){
473 johnpye 912 info = (struct Integ_var_t *)gl_fetch(sys->dynvars, i);
474 johnpye 928 if(info->derivative_of)continue;
475     if(info->derivative){
476     sys->y[yindex] = info->i;
477     sys->ydot[yindex] = info->derivative->i;
478     if(info->varindx >= 0){
479     sys->y_id[info->varindx] = yindex;
480     CONSOLE_DEBUG("y_id[%d] = %d",info->varindx,yindex);
481     }
482     if(info->derivative->varindx >= 0){
483     sys->y_id[info->derivative->varindx] = -1-yindex;
484     CONSOLE_DEBUG("y_id[%d] = %d",info->derivative->varindx,-1-yindex);
485     }
486 johnpye 669 }else{
487 johnpye 928 sys->y[yindex] = info ->i;
488     sys->ydot[yindex] = NULL;
489     if(info->varindx >= 0){
490     sys->y_id[info->varindx] = yindex;
491     CONSOLE_DEBUG("y_id[%d] = %d",info->varindx,yindex);
492     }
493 johnpye 669 }
494 johnpye 928 yindex++;
495 johnpye 669 }
496    
497 johnpye 912 nrels = slv_get_num_solvers_rels(sys->system);
498 johnpye 669 if(numy != nrels){
499     ERROR_REPORTER_HERE(ASC_USER_ERROR
500     ,"System is not square: solver has %d rels, found %d system states"
501     ,nrels, numy
502     );
503     return 0;
504     }
505    
506 johnpye 928 CONSOLE_DEBUG("THERE ARE %d VARIABLES IN THE INTEGRATION SYSTEM",numy);
507    
508 johnpye 912 sys->n_y = numy;
509 johnpye 669
510 johnpye 912 if(!integrator_sort_obs_vars(sys))return 0;
511 johnpye 669
512     return 1;
513     }
514    
515 johnpye 912 void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitfn){
516 johnpye 669 struct var_variable **vlist;
517     int i, vlen;
518    
519 ben.allan 704 /* visit all the slv_system_t master var lists to collect vars */
520     /* find the vars mostly in this one */
521 johnpye 912 vlist = slv_get_master_var_list(sys->system);
522     vlen = slv_get_num_master_vars(sys->system);
523 ben.allan 704 for (i=0;i<vlen;i++) {
524 johnpye 912 (*visitfn)(sys, vlist[i], &i);
525 ben.allan 704 }
526    
527 johnpye 669 /*
528     CONSOLE_DEBUG("Checked %d vars",vlen);
529 johnpye 912 integrator_print_var_stats(sys);
530 johnpye 669 */
531    
532 ben.allan 704 /* probably nothing here, but gotta check. */
533 johnpye 912 vlist = slv_get_master_par_list(sys->system);
534     vlen = slv_get_num_master_pars(sys->system);
535 ben.allan 704 for (i=0;i<vlen;i++) {
536 johnpye 912 (*visitfn)(sys, vlist[i], NULL);
537 ben.allan 704 }
538 johnpye 669
539     /*
540 ben.allan 704 CONSOLE_DEBUG("Checked %d pars",vlen);
541 johnpye 912 integrator_print_var_stats(sys);
542 johnpye 669 */
543    
544 ben.allan 704 /* might find t here */
545 johnpye 912 vlist = slv_get_master_unattached_list(sys->system);
546     vlen = slv_get_num_master_unattached(sys->system);
547 ben.allan 704 for (i=0;i<vlen;i++) {
548 johnpye 912 (*visitfn)(sys, vlist[i], NULL);
549 ben.allan 704 }
550 johnpye 669
551     /* CONSOLE_DEBUG("Checked %d unattached",vlen); */
552     }
553     /**
554     @return 1 on success
555     */
556 johnpye 912 int integrator_analyse_ode(IntegratorSystem *sys){
557 ben.allan 704 struct Integ_var_t *v1,*v2;
558 johnpye 709 long half,i,len;
559 ben.allan 704 int happy=1;
560 johnpye 894 char *varname1, *varname2;
561 johnpye 669
562 ben.allan 704 CONSOLE_DEBUG("Starting ODE analysis");
563 johnpye 669 IntegInitSymbols();
564    
565 ben.allan 704 /* collect potential states and derivatives */
566 johnpye 912 sys->indepvars = gl_create(10L); /* t var info */
567     sys->dynvars = gl_create(200L); /* y ydot var info */
568     sys->obslist = gl_create(100L); /* obs info */
569     if (sys->dynvars == NULL
570     || sys->obslist == NULL
571     || sys->indepvars == NULL
572 ben.allan 704 ){
573     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
574     return 0;
575     }
576 johnpye 669
577 johnpye 912 sys->nstates = sys->nderivs = 0;
578 johnpye 669
579 johnpye 912 integrator_visit_system_vars(sys,&integrator_ode_classify_var);
580 johnpye 669
581 johnpye 912 integrator_print_var_stats(sys);
582 johnpye 669
583 johnpye 912 if(!integrator_check_indep_var(sys))return 0;
584 johnpye 669
585 ben.allan 704 /* check sanity of state and var lists */
586    
587 johnpye 912 len = gl_length(sys->dynvars);
588 johnpye 669 half = len/2;
589 johnpye 709 CONSOLE_DEBUG("NUMBER OF DYNAMIC VARIABLES = %ld",half);
590 ben.allan 704
591 johnpye 912 if (len % 2 || len == 0L || sys->nstates != sys->nderivs ) {
592 ben.allan 704 /* list length must be even for vars to pair off */
593     ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"n_y != n_ydot, or no dynamic vars found. Fix your indexing.");
594     return 0;
595     }
596 johnpye 912 gl_sort(sys->dynvars,(CmpFunc)Integ_CmpDynVars);
597     if (gl_fetch(sys->dynvars,len)==NULL) {
598 ben.allan 704 ERROR_REPORTER_NOLINE(ASC_PROG_ERR,"Mysterious NULL found!");
599     return 0;
600     }
601 johnpye 912 sys->states = gl_create(half); /* state vars Integ_var_t references */
602     sys->derivs = gl_create(half); /* derivative var atoms */
603 ben.allan 704 for (i=1;i < len; i+=2) {
604 johnpye 912 v1 = (struct Integ_var_t *)gl_fetch(sys->dynvars,i);
605     v2 = (struct Integ_var_t *)gl_fetch(sys->dynvars,i+1);
606 ben.allan 704 if (v1->type!=1 || v2 ->type !=2 || v1->index != v2->index) {
607 johnpye 912 varname1 = var_make_name(sys->system,v1->i);
608     varname2 = var_make_name(sys->system,v2->i);
609 johnpye 908
610 johnpye 894 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Mistyped or misindexed dynamic variables: %s (%s = %ld,%s = %ld) and %s (%s = %ld,%s = %ld).",
611     varname1, SCP(STATEFLAG),v1->type,SCP(STATEINDEX),v1->index,
612     varname2, SCP(STATEFLAG),v2->type,SCP(STATEINDEX),v2->index
613 ben.allan 704 );
614 johnpye 894 ASC_FREE(varname1);
615     ASC_FREE(varname2);
616 ben.allan 704 happy=0;
617     break;
618     } else {
619 johnpye 912 gl_append_ptr(sys->states,(POINTER)v1);
620     gl_append_ptr(sys->derivs,(POINTER)v2->i);
621 ben.allan 704 }
622     }
623     if (!happy) {
624     return 0;
625     }
626 johnpye 912 sys->n_y = half;
627     sys->y = ASC_NEW_ARRAY(struct var_variable *, half);
628     sys->y_id = ASC_NEW_ARRAY(long, half);
629     sys->ydot = ASC_NEW_ARRAY(struct var_variable *, half);
630     if (sys->y==NULL || sys->ydot==NULL || sys->y_id==NULL) {
631 ben.allan 704 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
632     return 0;
633     }
634     for (i = 1; i <= half; i++) {
635 johnpye 912 v1 = (struct Integ_var_t *)gl_fetch(sys->states,i);
636     sys->y[i-1] = v1->i;
637     sys->y_id[i-1] = v1->index;
638     sys->ydot[i-1] = (struct var_variable *)gl_fetch(sys->derivs,i);
639 ben.allan 704 }
640 johnpye 669
641 johnpye 912 if(!integrator_sort_obs_vars(sys))return 0;
642 ben.allan 704
643     /* don't need the gl_lists now that we have arrays for everyone */
644 johnpye 912 gl_destroy(sys->states);
645     gl_destroy(sys->derivs);
646     gl_free_and_destroy(sys->indepvars); /* we own the objects in indepvars */
647     gl_free_and_destroy(sys->dynvars); /* we own the objects in dynvars */
648     gl_free_and_destroy(sys->obslist); /* and obslist */
649     sys->states = NULL;
650     sys->derivs = NULL;
651     sys->indepvars = NULL;
652     sys->dynvars = NULL;
653     sys->obslist = NULL;
654 johnpye 669
655 ben.allan 704 /* analysis completed OK */
656     return 1;
657 johnpye 669 }
658    
659     /**
660     Reindex observations. Sort if the user mostly numbered. Take natural order
661 ben.allan 704 if user just booleaned.
662 johnpye 741
663 johnpye 669 @return 1 on success
664     */
665 johnpye 912 static int integrator_sort_obs_vars(IntegratorSystem *sys){
666 johnpye 777 int half, i, len = 0;
667 ben.allan 704 struct Integ_var_t *v2;
668 johnpye 669
669 johnpye 912 half = sys->n_y;
670     len = gl_length(sys->obslist);
671 ben.allan 704 /* we shouldn't be seeing NULL here ever except if malloc fail. */
672     if (len > 1L) {
673 johnpye 912 half = ((struct Integ_var_t *)gl_fetch(sys->obslist,1))->index;
674 ben.allan 704 /* half != 0 now because we didn't collect 0 indexed vars */
675     for (i=2; i <= len; i++) {
676 johnpye 912 if (half != ((struct Integ_var_t *)gl_fetch(sys->obslist,i))->index) {
677 ben.allan 704 /* change seen. sort and go on */
678 johnpye 912 gl_sort(sys->obslist,(CmpFunc)Integ_CmpObs);
679 ben.allan 704 break;
680     }
681     }
682 johnpye 669 }
683 ben.allan 704 for (i = half = 1; i <= len; i++) {
684 johnpye 912 v2 = (struct Integ_var_t *)gl_fetch(sys->obslist,i);
685 ben.allan 704 if (v2==NULL) {
686     /* we shouldn't be seeing NULL here ever except if malloc fail. */
687 johnpye 912 gl_delete(sys->obslist,i,0); /* should not be gl_delete(so,i,1) */
688 ben.allan 704 } else {
689     Integ_SetObsId(v2->i,half);
690     v2->index = half++;
691     }
692     }
693 johnpye 669
694 ben.allan 704 /* obslist now uniquely indexed, no nulls */
695     /* make into arrays */
696 johnpye 912 half = gl_length(sys->obslist);
697     sys->obs = ASC_NEW_ARRAY(struct var_variable *,half);
698     sys->obs_id = ASC_NEW_ARRAY(long, half);
699     if ( sys->obs==NULL || sys->obs_id==NULL) {
700 ben.allan 704 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
701     return 0;
702     }
703 johnpye 912 sys->n_obs = half;
704 ben.allan 704 for (i = 1; i <= half; i++) {
705 johnpye 912 v2 = (struct Integ_var_t *)gl_fetch(sys->obslist,i);
706     sys->obs[i-1] = v2->i;
707     sys->obs_id[i-1] = v2->index;
708 ben.allan 704 }
709    
710     return 1;
711 johnpye 669 }
712    
713 johnpye 912 static void integrator_print_var_stats(IntegratorSystem *sys){
714     int v = gl_length(sys->dynvars);
715     int i = gl_length(sys->indepvars);
716 johnpye 669 CONSOLE_DEBUG("Currently %d vars, %d indep",v,i);
717     }
718    
719     /**
720 johnpye 928 Check sanity of the independent variable.
721    
722 johnpye 669 @return 1 on success
723     */
724 johnpye 912 static int integrator_check_indep_var(IntegratorSystem *sys){
725 johnpye 669 int len, i;
726     struct Integ_var_t *info;
727     char *varname;
728    
729 ben.allan 704 /* check the sanity of the independent variable */
730 johnpye 912 len = gl_length(sys->indepvars);
731 ben.allan 704 if (!len) {
732     ERROR_REPORTER_HERE(ASC_PROG_ERR,"No independent variable found.");
733     return 0;
734     }
735 johnpye 669 if (len > 1) {
736 ben.allan 704 ERROR_REPORTER_START_HERE(ASC_USER_ERROR);
737     FPRINTF(ASCERR,"Excess %ld independent variables found:",
738     len);
739     for(i=1; i <=len;i++) {
740 johnpye 912 info = (struct Integ_var_t *)gl_fetch(sys->indepvars,i);
741 johnpye 669 if(info==NULL)continue;
742    
743 johnpye 912 varname = var_make_name(sys->system,info->i);
744 johnpye 894 FPRINTF(ASCERR," %s",varname);
745 ben.allan 704 ASC_FREE(varname);
746     }
747     FPRINTF(ASCERR , "\nSet the \"%s\" flag on all but one of these to %s >= 0.\n"
748 johnpye 669 , SCP(STATEFLAG),SCP(STATEFLAG)
749     );
750 ben.allan 704 error_reporter_end_flush();
751     return 0;
752     }else{
753 johnpye 912 info = (struct Integ_var_t *)gl_fetch(sys->indepvars,1);
754     sys->x = info->i;
755 johnpye 669 }
756     return 1;
757     }
758    
759     /*------------------------------------------------------------------------------
760     CLASSIFICATION OF VARIABLES (for ANALYSIS step)
761     */
762    
763 johnpye 912 #define INTEG_ADD_TO_LIST(info,TYPE,INDEX,VAR,VARINDX,LIST) \
764 johnpye 669 info = ASC_NEW(struct Integ_var_t); \
765     if(info==NULL){ \
766     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory (INTEG_VAR_NEW)"); \
767     return; \
768     } \
769     info->type=TYPE; \
770     info->index=INDEX; \
771     info->i=VAR; \
772 johnpye 912 if(VARINDX==NULL){ \
773     info->varindx = -1; \
774     }else{ \
775     info->varindx = *VARINDX; \
776     } \
777 johnpye 669 gl_append_ptr(LIST,(void *)info); \
778     info = NULL
779    
780     /**
781     In a DAE, it's either the (single) independent variable, or it's a
782     variable in the model.
783    
784     I'm not sure what we should be doing with variables that are already
785     present as derivatives of other variables, I guess those ones need to be
786     removed from the list in a second pass?
787     */
788 johnpye 912 void integrator_dae_classify_var(IntegratorSystem *sys
789     , struct var_variable *var, const int *varindx
790     ){
791 johnpye 669 struct Integ_var_t *info;
792 ben.allan 704 long type,index;
793 johnpye 669
794     /* filter for recognition of solver_vars */
795 ben.allan 704 var_filter_t vfilt;
796     vfilt.matchbits = VAR_SVAR;
797     vfilt.matchvalue = VAR_SVAR;
798 johnpye 669
799 ben.allan 704 assert(var != NULL && var_instance(var)!=NULL );
800 johnpye 669
801     if( var_apply_filter(var,&vfilt) ) {
802 johnpye 921 if(!var_active(var)){
803     CONSOLE_DEBUG("VARIABLE IS NOT ACTIVE");
804     return;
805     }
806 johnpye 928
807     /* only non-fixed variables are accepted */
808 johnpye 669 if(!var_fixed(var)){
809 ben.allan 704 /* get the ode_type and ode_id of this solver_var */
810 johnpye 669 type = DynamicVarInfo(var,&index);
811    
812     if(type==INTEG_OTHER_VAR){
813     /* if the var's type is -1, it's independent */
814 johnpye 912 INTEG_ADD_TO_LIST(info,INTEG_OTHER_VAR,0,var,varindx,sys->indepvars);
815 johnpye 669 }else{
816     if(type < 0)type=0;
817     /* any other type of var is in the DAE system, at least for now */
818 johnpye 912 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
819 johnpye 669 }
820 johnpye 921 }
821 johnpye 928 #if 0
822     else{
823 johnpye 903 /* fixed variable, only include it if ode_type == 1 */
824     type = DynamicVarInfo(var,&index);
825     if(type==INTEG_STATE_VAR){
826 johnpye 912 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
827 johnpye 903 }
828 johnpye 669 }
829 johnpye 921 #endif
830 johnpye 669
831 ben.allan 704 /* if the var's obs_id > 0, add it to the observation list */
832 johnpye 669 if(ObservationVar(var,&index) != NULL && index > 0L) {
833 johnpye 912 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->obslist);
834 johnpye 669 }
835 ben.allan 704 }
836 johnpye 669 }
837    
838     /**
839     Inspect a specific variable and work out what type it is (what 'ode_type' it
840     has) and what other variable(s) it corresponds to (ie dydt corresponds to
841     y as a derivative).
842    
843     @TODO add ability to create new variables for 'missing' derivative vars?
844 ben.allan 704 */
845 johnpye 912 void integrator_ode_classify_var(IntegratorSystem *sys, struct var_variable *var
846     , const int *varindx
847     ){
848 ben.allan 704 struct Integ_var_t *info;
849     long type,index;
850 johnpye 669
851 ben.allan 704 var_filter_t vfilt;
852     vfilt.matchbits = VAR_SVAR;
853     vfilt.matchvalue = VAR_SVAR;
854 johnpye 669
855 ben.allan 704 assert(var != NULL && var_instance(var)!=NULL );
856    
857 johnpye 669 if( var_apply_filter(var,&vfilt) ) {
858 ben.allan 704 /* it's a solver var: what type of variable? */
859 johnpye 669 type = DynamicVarInfo(var,&index);
860 ben.allan 704
861 johnpye 669 if(type==INTEG_ALGEBRAIC_VAR){
862     /* no action required */
863     }else if(type==INTEG_OTHER_VAR){
864     /* i.e. independent var */
865 johnpye 912 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->indepvars);
866 johnpye 669 }else if(type>=INTEG_STATE_VAR){
867 johnpye 912 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
868 ben.allan 704 if(type == 1){
869 johnpye 912 sys->nstates++;
870 ben.allan 704 }else if(type == 2){ /* what about higher-order derivatives? -- JP */
871 johnpye 912 sys->nderivs++;
872 johnpye 669 }else{
873     ERROR_REPORTER_HERE(ASC_USER_WARNING,"Higher-order (>=2) derivatives are not supported in ODEs.");
874 johnpye 741 } }
875 johnpye 669
876     if(ObservationVar(var,&index) != NULL && index > 0L) {
877 johnpye 912 INTEG_ADD_TO_LIST(info,0L,index,var,varindx,sys->obslist);
878 ben.allan 704 }
879     }
880 johnpye 669 }
881    
882     /**
883     Look at a variable and determine if it's the independent variable or not.
884 johnpye 741 This is just for the purpose of the integrator_find_indep_var function,
885 johnpye 669 which is a utility function provided for use by the GUI.
886     */
887 johnpye 912 void integrator_classify_indep_var(IntegratorSystem *sys
888     , struct var_variable *var, const int *varindx
889     ){
890 ben.allan 704 struct Integ_var_t *info;
891     long type,index;
892 johnpye 669
893 ben.allan 704 var_filter_t vfilt;
894     vfilt.matchbits = VAR_SVAR;
895 johnpye 669 vfilt.matchvalue = VAR_SVAR;
896    
897 johnpye 725 /* CONSOLE_DEBUG("..."); */
898 johnpye 669
899 ben.allan 704 assert(var != NULL && var_instance(var)!=NULL );
900    
901 johnpye 669 if( var_apply_filter(var,&vfilt) ) {
902     type = DynamicVarInfo(var,&index);
903 ben.allan 704
904 johnpye 669 if(type==INTEG_OTHER_VAR){
905     /* i.e. independent var */
906 johnpye 912 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->indepvars);
907 johnpye 669 }
908     }
909 ben.allan 704 }
910    
911 johnpye 669 /**
912     Look at a variable, and if it is an 'ODE variable' (it has a child instance
913     named 'ode_type') return its type, which will be either:
914     - INTEG_OTHER_VAR (if 'ode_type' is -1)
915     - INTEG_ALGEBRAIC_VAR (if 'ode_type' is zero or any negative value < -1)
916     - INTEG_STATE_VAR (if 'ode_type' is 1)
917     - values 2, 3 or up, indicating derivatives (1st deriv=2, 2nd deriv=3, etc)
918    
919     If the parameter 'index' is not null, the value of 'ode_id' will be stuffed
920     there.
921 ben.allan 704 */
922     static long DynamicVarInfo(struct var_variable *v,long *index){
923     struct Instance *c, *d, *i;
924     i = var_instance(v);
925 johnpye 669 assert(i!=NULL);
926     assert(STATEFLAG!=NULL);
927 ben.allan 704 assert(STATEINDEX!=NULL);
928     c = ChildByChar(i,STATEFLAG);
929     d = ChildByChar(i,STATEINDEX);
930     /* lazy evaluation is important in the following if */
931     if(c == NULL
932     || d == NULL
933     || InstanceKind(c) != INTEGER_INST
934     || InstanceKind(d) != INTEGER_INST
935     || !AtomAssigned(c)
936 johnpye 669 || (!AtomAssigned(d) && GetIntegerAtomValue(c) != INTEG_OTHER_VAR)
937 ben.allan 704 ){
938     return INTEG_ALGEBRAIC_VAR;
939     }
940     if (index != NULL) {
941     *index = GetIntegerAtomValue(d);
942     }
943     return GetIntegerAtomValue(c);
944     }
945 johnpye 669
946     /**
947     Looks at the given variable checks if it is an 'observation variable'. This
948     means that it has its 'obs_id' child instance set to a non-zero value.
949    
950     If the variable is an observation variable, its index value ('obs_id') is
951     stuff into *index (provided index!=NULL), and the pointer to the original
952     instance is rtruend.
953    
954 ben.allan 704 If it's not an observation variable, we return NULL and *index is untouched.
955     */
956     static struct var_variable *ObservationVar(struct var_variable *v, long *index){
957     struct Instance *c,*i;
958     i = var_instance(v);
959     assert(i!=NULL);
960     c = ChildByChar(i,OBSINDEX);
961     if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
962     return NULL;
963     }
964     if (index != NULL) {
965     *index = GetIntegerAtomValue(c);
966     }
967     return v;
968     }
969 johnpye 669
970     /*------------------------------------------------------------------------------
971     RUNNING THE SOLVER
972     */
973    
974     /*
975     Make the call to the actual integrator we've selected, for the range of
976 johnpye 912 time values specified. The sys contains all the specifics.
977 johnpye 709
978     Return 1 on success
979 ben.allan 704 */
980 johnpye 912 int integrator_solve(IntegratorSystem *sys, long i0, long i1){
981 johnpye 669
982 ben.allan 704 long nstep;
983 johnpye 741 unsigned long start_index=0, finish_index=0;
984 johnpye 912 assert(sys!=NULL);
985 johnpye 669
986 johnpye 912 nstep = integrator_getnsamples(sys)-1;
987 ben.allan 704 /* check for at least 2 steps and dimensionality of x vs steps here */
988 johnpye 669
989     if (i0<0 || i1 <0) {
990 ben.allan 704 /* dude, there's no way we're writing interactive stuff here... */
991 johnpye 669 ERROR_REPORTER_HERE(ASC_PROG_ERROR,"Console input of integration limits has been disabled!");
992 ben.allan 704 return 0;
993 johnpye 908 } else {
994 ben.allan 704 start_index=i0;
995     finish_index =i1;
996     if (start_index >= (unsigned long)nstep) {
997 johnpye 669 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Start point (=%lu) must be an index in the range [0,%li]."
998     ,start_index,nstep
999 ben.allan 704 );
1000     return 0;
1001     }
1002     if (finish_index > (unsigned long)nstep) {
1003 johnpye 669 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"End point (=%lu) must be an index in the range [0,%li]."
1004     ,finish_index,nstep
1005 ben.allan 704 );
1006     return 0;
1007     }
1008     }
1009 johnpye 669
1010 ben.allan 704 if(finish_index <= start_index) {
1011 johnpye 669 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"End point comes before start point! (start=%lu, end=%lu)"
1012     ,start_index,finish_index
1013 ben.allan 704 );
1014     return 0;
1015 johnpye 669 }
1016    
1017 johnpye 888 CONSOLE_DEBUG("RUNNING INTEGRATION...");
1018    
1019 ben.allan 704 /* now go and run the integrator */
1020 johnpye 912 switch (sys->engine) {
1021 johnpye 908 case INTEG_LSODE:
1022 johnpye 912 return integrator_lsode_solve(sys, start_index, finish_index);
1023 johnpye 908 break;
1024 johnpye 669 #ifdef ASC_WITH_IDA
1025 johnpye 908 case INTEG_IDA:
1026 johnpye 912 return integrator_ida_solve(sys,start_index, finish_index);
1027 johnpye 908 break;
1028 ben.allan 704 #endif
1029 johnpye 709 default:
1030     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown integrator (invalid, or not implemented yet)");
1031     return 0;
1032 ben.allan 704 }
1033 johnpye 669 }
1034 ben.allan 704
1035 johnpye 669 /*---------------------------------------------------------------
1036     HANDLING THE LIST OF TIMESTEMPS
1037     */
1038    
1039     #define GETTER_AND_SETTER(TYPE,NAME) \
1040 johnpye 912 void integrator_set_##NAME(IntegratorSystem *sys, TYPE val){ \
1041     sys->NAME=val; \
1042 johnpye 669 } \
1043 johnpye 912 TYPE integrator_get_##NAME(IntegratorSystem *sys){ \
1044     return sys->NAME; \
1045 johnpye 669 }
1046    
1047 johnpye 908 GETTER_AND_SETTER(SampleList *,samples) /*;*/
1048     GETTER_AND_SETTER(double,maxstep) /*;*/
1049     GETTER_AND_SETTER(double,minstep) /*;*/
1050     GETTER_AND_SETTER(double,stepzero) /*;*/
1051     GETTER_AND_SETTER(int,maxsubsteps) /*;*/
1052 ben.allan 704 #undef GETTER_AND_SETTER
1053    
1054 johnpye 912 long integrator_getnsamples(IntegratorSystem *sys){
1055     assert(sys!=NULL);
1056     assert(sys->samples!=NULL);
1057     return samplelist_length(sys->samples);
1058 ben.allan 704 }
1059    
1060 johnpye 912 double integrator_getsample(IntegratorSystem *sys, long i){
1061     assert(sys!=NULL);
1062     assert(sys->samples!=NULL);
1063     return samplelist_get(sys->samples,i);
1064 ben.allan 704 }
1065    
1066 johnpye 912 void integrator_setsample(IntegratorSystem *sys, long i,double xi){
1067     assert(sys!=NULL);
1068     assert(sys->samples!=NULL);
1069     samplelist_set(sys->samples,i,xi);
1070 ben.allan 704 }
1071    
1072 johnpye 912 const dim_type *integrator_getsampledim(IntegratorSystem *sys){
1073     assert(sys!=NULL);
1074     assert(sys->samples!=NULL);
1075     return samplelist_dim(sys->samples);
1076 ben.allan 704 }
1077 johnpye 669
1078 johnpye 912 ASC_DLLSPEC(long) integrator_getcurrentstep(IntegratorSystem *sys){
1079     return sys->currentstep;
1080 johnpye 669 }
1081 ben.allan 704
1082 johnpye 669 /*------------------------------------------------------------------------------
1083 ben.allan 704 GET/SET VALUE OF THE INDEP VARIABLE
1084 johnpye 669 */
1085    
1086 johnpye 741 /**
1087 johnpye 669 Retrieve the value of the independent variable (time) from ASCEND
1088     and return it as a double.
1089 johnpye 741 */
1090 johnpye 912 double integrator_get_t(IntegratorSystem *sys){
1091     assert(sys->x!=NULL);
1092     return var_value(sys->x);
1093 ben.allan 704 }
1094 johnpye 669
1095     /**
1096     Set the value of the independent variable (time) in ASCEND.
1097 ben.allan 704 */
1098 johnpye 912 void integrator_set_t(IntegratorSystem *sys, double value){
1099     var_set_value(sys->x, value);
1100 johnpye 815 /* CONSOLE_DEBUG("set_t = %g", value); */
1101 ben.allan 704 }
1102    
1103 johnpye 669 /*------------------------------------------------------------------------------
1104     PASSING DIFFERENTIAL VARIABLES AND THEIR DERIVATIVES TO/FROM THE SOLVER
1105     */
1106     /**
1107     Retrieve the current values of the derivatives of the y-variables
1108     and stick them in the/an array that the integrator will use.
1109    
1110     If the pointer 'y' is NULL, the necessary space is allocated (and
1111     must be freed somewhere else).
1112 ben.allan 704 */
1113 johnpye 912 double *integrator_get_y(IntegratorSystem *sys, double *y) {
1114 ben.allan 704 long i;
1115    
1116     if (y==NULL) {
1117 johnpye 912 y = ASC_NEW_ARRAY_CLEAR(double, sys->n_y+1);
1118 ben.allan 704 /* C y[0] <==> ascend d.y[1] <==> f77 y(1) */
1119     }
1120    
1121 johnpye 912 for (i=0; i< sys->n_y; i++) {
1122     assert(sys->y[i]!=NULL);
1123     y[i] = var_value(sys->y[i]);
1124 johnpye 815 /* CONSOLE_DEBUG("ASCEND --> y[%ld] = %g", i+1, y[i]); */
1125 ben.allan 704 }
1126     return y;
1127     }
1128 johnpye 669
1129     /**
1130     Take the values of the derivatives from the array that the integrator
1131     uses, and use them to update the values of the corresponding variables
1132     in ASCEND.
1133 ben.allan 704 */
1134 johnpye 912 void integrator_set_y(IntegratorSystem *sys, double *y) {
1135 johnpye 669 long i;
1136 ben.allan 704 #ifndef NDEBUG
1137 johnpye 669 char *varname;
1138 ben.allan 704 #endif
1139    
1140 johnpye 912 for (i=0; i < sys->n_y; i++) {
1141     assert(sys->y[i]!=NULL);
1142     var_set_value(sys->y[i],y[i]);
1143 johnpye 669 #ifndef NDEBUG
1144 johnpye 912 varname = var_make_name(sys->system, sys->y[i]);
1145 johnpye 815 /* CONSOLE_DEBUG("y[%ld] = \"%s\" = %g --> ASCEND", i+1, varname, y[i]); */
1146 ben.allan 704 ASC_FREE(varname);
1147 johnpye 669 #endif
1148 ben.allan 704 }
1149     }
1150 johnpye 669
1151     /**
1152     Send the values of the derivatives of the 'y' variables to the solver.
1153     Allocate space for an array if necessary.
1154    
1155 johnpye 912 Any element in sys->ydot that is NULL will be passed over (the value
1156 johnpye 669 won't be modified in dydx).
1157 ben.allan 704 */
1158 johnpye 912 double *integrator_get_ydot(IntegratorSystem *sys, double *dydx) {
1159 ben.allan 704 long i;
1160    
1161     if (dydx==NULL) {
1162 johnpye 912 dydx = ASC_NEW_ARRAY_CLEAR(double, sys->n_y+1);
1163 ben.allan 704 /* C dydx[0] <==> ascend d.dydx[1] <==> f77 ydot(1) */
1164     }
1165    
1166 johnpye 912 for (i=0; i < sys->n_y; i++) {
1167     if(sys->ydot[i]!=NULL){
1168     dydx[i] = var_value(sys->ydot[i]);
1169 ben.allan 704 }
1170 johnpye 815 /* CONSOLE_DEBUG("ASCEND --> ydot[%ld] = %g", i+1, dydx[i]); */
1171 ben.allan 704 }
1172     return dydx;
1173     }
1174    
1175 johnpye 912 void integrator_set_ydot(IntegratorSystem *sys, double *dydx) {
1176 ben.allan 704 long i;
1177     #ifndef NDEBUG
1178 johnpye 908 /* char *varname; */
1179 ben.allan 704 #endif
1180 johnpye 912 for (i=0; i < sys->n_y; i++) {
1181     if(sys->ydot[i]!=NULL){
1182     var_set_value(sys->ydot[i],dydx[i]);
1183 johnpye 669 #ifndef NDEBUG
1184 johnpye 912 /* varname = var_make_name(sys->system, sys->ydot[i]);
1185 johnpye 709 CONSOLE_DEBUG("ydot[%ld] = \"%s\" = %g --> ASCEND", i+1, varname, dydx[i]);
1186 johnpye 903 ASC_FREE(varname); */
1187 johnpye 669 #endif
1188     }
1189     #ifndef NDEBUG
1190 johnpye 903 /*else{
1191 johnpye 709 CONSOLE_DEBUG("ydot[%ld] = %g (internal)", i+1, dydx[i]);
1192 johnpye 903 }*/
1193 johnpye 669 #endif
1194 ben.allan 704 }
1195     }
1196    
1197 johnpye 669 /*-------------------------------------------------------------
1198 ben.allan 704 RETRIEVING OBSERVATION DATA
1199 johnpye 908 */
1200 ben.allan 704
1201 johnpye 908 /**
1202 ben.allan 704 This function takes the inst in the solver and returns the vector of
1203     observation variables that are located in the submodel d.obs array.
1204     */
1205 johnpye 912 double *integrator_get_observations(IntegratorSystem *sys, double *obsi) {
1206 ben.allan 704 long i;
1207    
1208     if (obsi==NULL) {
1209 johnpye 912 obsi = ASC_NEW_ARRAY_CLEAR(double, sys->n_obs+1);
1210 ben.allan 704 }
1211    
1212     /* C obsi[0] <==> ascend d.obs[1] */
1213    
1214 johnpye 912 for (i=0; i < sys->n_obs; i++) {
1215     obsi[i] = var_value(sys->obs[i]);
1216 ben.allan 704 /* CONSOLE_DEBUG("*get_d_obs[%ld] = %g\n", i+1, obsi[i]); */
1217     }
1218     return obsi;
1219 johnpye 669 }
1220    
1221 johnpye 912 struct var_variable *integrator_get_observed_var(IntegratorSystem *sys, const long i){
1222 johnpye 669 assert(i>=0);
1223 johnpye 912 assert(i<sys->n_obs);
1224     return sys->obs[i];
1225 ben.allan 704 }
1226    
1227 johnpye 854 /**
1228     @NOTE Although this shouldn't be required for implementation of solver
1229     engines, this is useful for GUI reporting of integration results.
1230     */
1231 johnpye 912 struct var_variable *integrator_get_independent_var(IntegratorSystem *sys){
1232     return sys->x;
1233 johnpye 854 }
1234    
1235    
1236 johnpye 669 /*----------------------------------------------------
1237 johnpye 741 Build an analytic jacobian for solving the state system
1238 johnpye 669
1239 ben.allan 704 This necessarily ugly piece of code attempts to create a unique
1240     list of relations that explicitly contain the variables in the
1241     given input list. The utility of this information is that we know
1242     exactly which relations must be differentiated, to fill in the
1243     df/dy matrix. If the problem has very few derivative terms, this will
1244     be of great savings. If the problem arose from the discretization of
1245     a pde, then this will be not so useful. The decision wether to use
1246     this function or to simply differentiate the entire relations list
1247 johnpye 741 must be done before calling this function.
1248 johnpye 669
1249 ben.allan 704 Final Note: the callee owns the array, but not the array elements.
1250     */
1251     #define AVG_NUM_INCIDENT 4
1252 johnpye 669
1253 ben.allan 704
1254 johnpye 669 /**
1255     This function helps to arrange the observation variables in a sensible order.
1256     The 'obs_id' child instance of v, if present, is assigned the value of the
1257     given parameter 'index'.
1258 ben.allan 704 */
1259     static void Integ_SetObsId(struct var_variable *v, long index){
1260     struct Instance *c, *i;
1261     i = var_instance(v);
1262     assert(i!=NULL);
1263     c = ChildByChar(i,OBSINDEX);
1264     if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
1265     return;
1266     }
1267     SetIntegerAtomValue(c,index,0);
1268     }
1269    
1270 johnpye 669 /**
1271 johnpye 741 Compares observation structs. NULLs should end up at far end.
1272 ben.allan 704 */
1273     static int Integ_CmpObs(struct Integ_var_t *v1, struct Integ_var_t *v2){
1274     if(v1 == NULL)return 1;
1275     if(v2 == NULL)return -1;
1276     if(v1->index > v2->index)return 1;
1277     if(v1->index == v2->index)return 0;
1278     return -1;
1279     }
1280    
1281     /**
1282     Compares dynamic vars structs. NULLs should end up at far end.
1283     List should be sorted primarily by index and then by type, in order
1284     of increasing value of both.
1285     */
1286     static int Integ_CmpDynVars(struct Integ_var_t *v1, struct Integ_var_t *v2){
1287     if(v1 == NULL)return 1;
1288     if(v2 == NULL)return -1;
1289     if(v1->index > v2->index)return 1;
1290     if(v1->index != v2->index)return -1;
1291     if(v1->type > v2->type)return 1;
1292     return -1;
1293     }
1294 johnpye 669 /*----------------------------
1295     Output handling to the GUI/interface.
1296     */
1297    
1298 johnpye 912 int integrator_set_reporter(IntegratorSystem *sys
1299 johnpye 669 , IntegratorReporter *reporter
1300     ){
1301 johnpye 912 assert(sys!=NULL);
1302     sys->reporter = reporter;
1303 johnpye 854 /* ERROR_REPORTER_HERE(ASC_PROG_NOTE,"INTEGRATOR REPORTER HOOKS HAVE BEEN SET\n"); */
1304 johnpye 709 return 1;
1305 johnpye 669 }
1306    
1307 johnpye 912 int integrator_output_init(IntegratorSystem *sys){
1308     assert(sys!=NULL);
1309     assert(sys->reporter!=NULL);
1310     if(sys->reporter->init!=NULL){
1311 johnpye 669 /* call the specified output function */
1312 johnpye 912 return (*(sys->reporter->init))(sys);
1313 johnpye 669 }
1314     ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter init method");
1315     return 1;
1316     }
1317    
1318 johnpye 912 int integrator_output_write(IntegratorSystem *sys){
1319 johnpye 669 static int reported_already=0;
1320 johnpye 912 assert(sys!=NULL);
1321     if(sys->reporter->write!=NULL){
1322     return (*(sys->reporter->write))(sys);
1323 johnpye 669 }
1324     if(!reported_already){
1325     ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter write method (this message only shown once)");
1326     reported_already=1;
1327     }
1328     return 1;
1329     }
1330    
1331 johnpye 912 int integrator_output_write_obs(IntegratorSystem *sys){
1332 johnpye 669 static int reported_already=0;
1333 johnpye 912 assert(sys!=NULL);
1334     if(sys->reporter->write_obs!=NULL){
1335     return (*(sys->reporter->write_obs))(sys);
1336 johnpye 669 }
1337     if(!reported_already){
1338     ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter write_obs method (this message only shown once)");
1339     reported_already=1;
1340     }
1341     return 1;
1342 johnpye 741 }
1343 johnpye 669
1344 johnpye 912 int integrator_output_close(IntegratorSystem *sys){
1345     assert(sys!=NULL);
1346     if(sys->reporter->close!=NULL){
1347     return (*(sys->reporter->close))(sys);
1348 johnpye 669 }
1349     ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter close method");
1350     return 1;
1351     }
1352    
1353     /**
1354     Decode status codes from the integrator, and output them via FPRINTF.
1355 ben.allan 704 */
1356     int integrator_checkstatus(slv_status_t status) {
1357     if (status.converged) {
1358     return 1;
1359     }
1360     if (status.diverged) {
1361 johnpye 903 FPRINTF(stderr, "The derivative system did not converge. Integration will terminate.");
1362 ben.allan 704 return 0;
1363     }
1364     if (status.inconsistent) {
1365 johnpye 903 FPRINTF(stderr, "A numerically inconsistent state was discovered while "
1366     "calculating derivatives. Integration will terminate.");
1367 ben.allan 704 return 0;
1368     }
1369     if (status.time_limit_exceeded) {
1370 johnpye 854 FPRINTF(stderr, "The time limit was exceeded while calculating "
1371 johnpye 903 "derivatives. Integration will terminate.");
1372 ben.allan 704 return 0;
1373     }
1374     if (status.iteration_limit_exceeded) {
1375 johnpye 854 FPRINTF(stderr, "The iteration limit was exceeded while calculating "
1376 johnpye 903 "derivatives. Integration will terminate.");
1377 ben.allan 704 return 0;
1378     }
1379     if (status.panic) {
1380 johnpye 854 FPRINTF(stderr, "The user patience limit was exceeded while "
1381 johnpye 903 "calculating derivatives. Integration will terminate.");
1382 ben.allan 704 return 0;
1383     }
1384     return 0;
1385 johnpye 669 }

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