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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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