/[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 854 - (show annotations) (download) (as text)
Wed Sep 20 13:36:40 2006 UTC (13 years, 9 months ago) by johnpye
File MIME type: text/x-csrc
File size: 40830 byte(s)
First tentative version in 'integration reporting':
Values of observed variables from the simulation are added to an Observer table after simulation completes.
This is not very efficiently coded at this stage but is a start.
Also some minor changes to text and comments in some base/generic files.
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, i, len = 0;
674 struct Integ_var_t *v2;
675
676 half = blsys->n_y;
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 unsigned long start_index=0, finish_index=0;
964
965 assert(blsys!=NULL);
966
967 nstep = integrator_getnsamples(blsys)-1;
968 /* check for at least 2 steps and dimensionality of x vs steps here */
969
970 if (i0<0 || i1 <0) {
971 /* dude, there's no way we're writing interactive stuff here... */
972 ERROR_REPORTER_HERE(ASC_PROG_ERROR,"Console input of integration limits has been disabled!");
973 return 0;
974 }else{
975 start_index=i0;
976 finish_index =i1;
977 if (start_index >= (unsigned long)nstep) {
978 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Start point (=%lu) must be an index in the range [0,%li]."
979 ,start_index,nstep
980 );
981 return 0;
982 }
983 if (finish_index > (unsigned long)nstep) {
984 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"End point (=%lu) must be an index in the range [0,%li]."
985 ,finish_index,nstep
986 );
987 return 0;
988 }
989 }
990
991 if(finish_index <= start_index) {
992 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"End point comes before start point! (start=%lu, end=%lu)"
993 ,start_index,finish_index
994 );
995 return 0;
996 }
997
998 /* now go and run the integrator */
999 switch (blsys->engine) {
1000 case INTEG_LSODE: return integrator_lsode_solve(blsys, start_index, finish_index); break;
1001 #ifdef ASC_WITH_IDA
1002 case INTEG_IDA: return integrator_ida_solve(blsys,start_index, finish_index); break;
1003 #endif
1004 default:
1005 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown integrator (invalid, or not implemented yet)");
1006 return 0;
1007 }
1008 }
1009
1010 /*---------------------------------------------------------------
1011 HANDLING THE LIST OF TIMESTEMPS
1012 */
1013
1014 #define GETTER_AND_SETTER(TYPE,NAME) \
1015 void integrator_set_##NAME(IntegratorSystem *blsys, TYPE val){ \
1016 blsys->NAME=val; \
1017 } \
1018 TYPE integrator_get_##NAME(IntegratorSystem *blsys){ \
1019 return blsys->NAME; \
1020 }
1021
1022 GETTER_AND_SETTER(SampleList *,samples);
1023 GETTER_AND_SETTER(double,maxstep);
1024 GETTER_AND_SETTER(double,minstep);
1025 GETTER_AND_SETTER(double,stepzero);
1026 GETTER_AND_SETTER(int,maxsubsteps);
1027 #undef GETTER_AND_SETTER
1028
1029 long integrator_getnsamples(IntegratorSystem *blsys){
1030 assert(blsys!=NULL);
1031 assert(blsys->samples!=NULL);
1032 return samplelist_length(blsys->samples);
1033 }
1034
1035 double integrator_getsample(IntegratorSystem *blsys, long i){
1036 assert(blsys!=NULL);
1037 assert(blsys->samples!=NULL);
1038 return samplelist_get(blsys->samples,i);
1039 }
1040
1041 void integrator_setsample(IntegratorSystem *blsys, long i,double xi){
1042 assert(blsys!=NULL);
1043 assert(blsys->samples!=NULL);
1044 samplelist_set(blsys->samples,i,xi);
1045 }
1046
1047 const dim_type *integrator_getsampledim(IntegratorSystem *blsys){
1048 assert(blsys!=NULL);
1049 assert(blsys->samples!=NULL);
1050 return samplelist_dim(blsys->samples);
1051 }
1052
1053 ASC_DLLSPEC(long) integrator_getcurrentstep(IntegratorSystem *blsys){
1054 return blsys->currentstep;
1055 }
1056
1057 /*------------------------------------------------------------------------------
1058 GET/SET VALUE OF THE INDEP VARIABLE
1059 */
1060
1061 /**
1062 Retrieve the value of the independent variable (time) from ASCEND
1063 and return it as a double.
1064 */
1065 double integrator_get_t(IntegratorSystem *blsys){
1066 assert(blsys->x!=NULL);
1067 return var_value(blsys->x);
1068 }
1069
1070 /**
1071 Set the value of the independent variable (time) in ASCEND.
1072 */
1073 void integrator_set_t(IntegratorSystem *blsys, double value){
1074 var_set_value(blsys->x, value);
1075 /* CONSOLE_DEBUG("set_t = %g", value); */
1076 }
1077
1078 /*------------------------------------------------------------------------------
1079 PASSING DIFFERENTIAL VARIABLES AND THEIR DERIVATIVES TO/FROM THE SOLVER
1080 */
1081 /**
1082 Retrieve the current values of the derivatives of the y-variables
1083 and stick them in the/an array that the integrator will use.
1084
1085 If the pointer 'y' is NULL, the necessary space is allocated (and
1086 must be freed somewhere else).
1087 */
1088 double *integrator_get_y(IntegratorSystem *blsys, double *y) {
1089 long i;
1090
1091 if (y==NULL) {
1092 y = ASC_NEW_ARRAY_CLEAR(double, blsys->n_y+1);
1093 /* C y[0] <==> ascend d.y[1] <==> f77 y(1) */
1094 }
1095
1096 for (i=0; i< blsys->n_y; i++) {
1097 assert(blsys->y[i]!=NULL);
1098 y[i] = var_value(blsys->y[i]);
1099 /* CONSOLE_DEBUG("ASCEND --> y[%ld] = %g", i+1, y[i]); */
1100 }
1101 return y;
1102 }
1103
1104 /**
1105 Take the values of the derivatives from the array that the integrator
1106 uses, and use them to update the values of the corresponding variables
1107 in ASCEND.
1108 */
1109 void integrator_set_y(IntegratorSystem *blsys, double *y) {
1110 long i;
1111 #ifndef NDEBUG
1112 char *varname;
1113 #endif
1114
1115 for (i=0; i < blsys->n_y; i++) {
1116 assert(blsys->y[i]!=NULL);
1117 var_set_value(blsys->y[i],y[i]);
1118 #ifndef NDEBUG
1119 varname = var_make_name(blsys->system, blsys->y[i]);
1120 /* CONSOLE_DEBUG("y[%ld] = \"%s\" = %g --> ASCEND", i+1, varname, y[i]); */
1121 ASC_FREE(varname);
1122 #endif
1123 }
1124 }
1125
1126 /**
1127 Send the values of the derivatives of the 'y' variables to the solver.
1128 Allocate space for an array if necessary.
1129
1130 Any element in blsys->ydot that is NULL will be passed over (the value
1131 won't be modified in dydx).
1132 */
1133 double *integrator_get_ydot(IntegratorSystem *blsys, double *dydx) {
1134 long i;
1135
1136 if (dydx==NULL) {
1137 dydx = ASC_NEW_ARRAY_CLEAR(double, blsys->n_y+1);
1138 /* C dydx[0] <==> ascend d.dydx[1] <==> f77 ydot(1) */
1139 }
1140
1141 for (i=0; i < blsys->n_y; i++) {
1142 if(blsys->ydot[i]!=NULL){
1143 dydx[i] = var_value(blsys->ydot[i]);
1144 }
1145 /* CONSOLE_DEBUG("ASCEND --> ydot[%ld] = %g", i+1, dydx[i]); */
1146 }
1147 return dydx;
1148 }
1149
1150 void integrator_set_ydot(IntegratorSystem *blsys, double *dydx) {
1151 long i;
1152 #ifndef NDEBUG
1153 char *varname;
1154 #endif
1155 for (i=0; i < blsys->n_y; i++) {
1156 if(blsys->ydot[i]!=NULL){
1157 var_set_value(blsys->ydot[i],dydx[i]);
1158 #ifndef NDEBUG
1159 varname = var_make_name(blsys->system, blsys->ydot[i]);
1160 CONSOLE_DEBUG("ydot[%ld] = \"%s\" = %g --> ASCEND", i+1, varname, dydx[i]);
1161 ASC_FREE(varname);
1162 #endif
1163 }
1164 #ifndef NDEBUG
1165 else{
1166 CONSOLE_DEBUG("ydot[%ld] = %g (internal)", i+1, dydx[i]);
1167 }
1168 #endif
1169 }
1170 }
1171
1172 /*-------------------------------------------------------------
1173 RETRIEVING OBSERVATION DATA
1174
1175 This function takes the inst in the solver and returns the vector of
1176 observation variables that are located in the submodel d.obs array.
1177 */
1178 double *integrator_get_observations(IntegratorSystem *blsys, double *obsi) {
1179 long i;
1180
1181 if (obsi==NULL) {
1182 obsi = ASC_NEW_ARRAY_CLEAR(double, blsys->n_obs+1);
1183 }
1184
1185 /* C obsi[0] <==> ascend d.obs[1] */
1186
1187 for (i=0; i < blsys->n_obs; i++) {
1188 obsi[i] = var_value(blsys->obs[i]);
1189 /* CONSOLE_DEBUG("*get_d_obs[%ld] = %g\n", i+1, obsi[i]); */
1190 }
1191 return obsi;
1192 }
1193
1194 struct var_variable *integrator_get_observed_var(IntegratorSystem *blsys, const long i){
1195 assert(i>=0);
1196 assert(i<blsys->n_obs);
1197 return blsys->obs[i];
1198 }
1199
1200 /**
1201 @NOTE Although this shouldn't be required for implementation of solver
1202 engines, this is useful for GUI reporting of integration results.
1203 */
1204 struct var_variable *integrator_get_independent_var(IntegratorSystem *blsys){
1205 return blsys->x;
1206 }
1207
1208
1209 /*----------------------------------------------------
1210 Build an analytic jacobian for solving the state system
1211
1212 This necessarily ugly piece of code attempts to create a unique
1213 list of relations that explicitly contain the variables in the
1214 given input list. The utility of this information is that we know
1215 exactly which relations must be differentiated, to fill in the
1216 df/dy matrix. If the problem has very few derivative terms, this will
1217 be of great savings. If the problem arose from the discretization of
1218 a pde, then this will be not so useful. The decision wether to use
1219 this function or to simply differentiate the entire relations list
1220 must be done before calling this function.
1221
1222 Final Note: the callee owns the array, but not the array elements.
1223 */
1224 #define AVG_NUM_INCIDENT 4
1225
1226
1227 /**
1228 This function helps to arrange the observation variables in a sensible order.
1229 The 'obs_id' child instance of v, if present, is assigned the value of the
1230 given parameter 'index'.
1231 */
1232 static void Integ_SetObsId(struct var_variable *v, long index){
1233 struct Instance *c, *i;
1234 i = var_instance(v);
1235 assert(i!=NULL);
1236 c = ChildByChar(i,OBSINDEX);
1237 if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
1238 return;
1239 }
1240 SetIntegerAtomValue(c,index,0);
1241 }
1242
1243 /**
1244 Compares observation structs. NULLs should end up at far end.
1245 */
1246 static int Integ_CmpObs(struct Integ_var_t *v1, struct Integ_var_t *v2){
1247 if(v1 == NULL)return 1;
1248 if(v2 == NULL)return -1;
1249 if(v1->index > v2->index)return 1;
1250 if(v1->index == v2->index)return 0;
1251 return -1;
1252 }
1253
1254 /**
1255 Compares dynamic vars structs. NULLs should end up at far end.
1256 List should be sorted primarily by index and then by type, in order
1257 of increasing value of both.
1258 */
1259 static int Integ_CmpDynVars(struct Integ_var_t *v1, struct Integ_var_t *v2){
1260 if(v1 == NULL)return 1;
1261 if(v2 == NULL)return -1;
1262 if(v1->index > v2->index)return 1;
1263 if(v1->index != v2->index)return -1;
1264 if(v1->type > v2->type)return 1;
1265 return -1;
1266 }
1267 /*----------------------------
1268 Output handling to the GUI/interface.
1269 */
1270
1271 int integrator_set_reporter(IntegratorSystem *blsys
1272 , IntegratorReporter *reporter
1273 ){
1274 assert(blsys!=NULL);
1275 blsys->reporter = reporter;
1276 /* ERROR_REPORTER_HERE(ASC_PROG_NOTE,"INTEGRATOR REPORTER HOOKS HAVE BEEN SET\n"); */
1277 return 1;
1278 }
1279
1280 int integrator_output_init(IntegratorSystem *blsys){
1281 CONSOLE_DEBUG("...");
1282 assert(blsys!=NULL);
1283 assert(blsys->reporter!=NULL);
1284 if(blsys->reporter->init!=NULL){
1285 /* call the specified output function */
1286 return (*(blsys->reporter->init))(blsys);
1287 }
1288 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter init method");
1289 return 1;
1290 }
1291
1292 int integrator_output_write(IntegratorSystem *blsys){
1293 static int reported_already=0;
1294 assert(blsys!=NULL);
1295 if(blsys->reporter->write!=NULL){
1296 return (*(blsys->reporter->write))(blsys);
1297 }
1298 if(!reported_already){
1299 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter write method (this message only shown once)");
1300 reported_already=1;
1301 }
1302 return 1;
1303 }
1304
1305 int integrator_output_write_obs(IntegratorSystem *blsys){
1306 static int reported_already=0;
1307 assert(blsys!=NULL);
1308 if(blsys->reporter->write_obs!=NULL){
1309 return (*(blsys->reporter->write_obs))(blsys);
1310 }
1311 if(!reported_already){
1312 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter write_obs method (this message only shown once)");
1313 reported_already=1;
1314 }
1315 return 1;
1316 }
1317
1318 int integrator_output_close(IntegratorSystem *blsys){
1319 assert(blsys!=NULL);
1320 if(blsys->reporter->close!=NULL){
1321 return (*(blsys->reporter->close))(blsys);
1322 }
1323 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter close method");
1324 return 1;
1325 }
1326
1327 /**
1328 Decode status codes from the integrator, and output them via FPRINTF.
1329 */
1330 int integrator_checkstatus(slv_status_t status) {
1331 if (status.converged) {
1332 return 1;
1333 }
1334 if (status.diverged) {
1335 FPRINTF(stderr, "The derivative system did not converge.\nIntegration "
1336 "will be terminated at the end of the current step.");
1337 return 0;
1338 }
1339 if (status.inconsistent) {
1340 FPRINTF(stderr, "A numerical inconsistency was discovered while "
1341 "calculating derivatives. Integration will be terminated at "
1342 "the end of the current step.");
1343 return 0;
1344 }
1345 if (status.time_limit_exceeded) {
1346 FPRINTF(stderr, "The time limit was exceeded while calculating "
1347 "derivatives.\nIntegration will be terminated at the end of "
1348 "the current step.");
1349 return 0;
1350 }
1351 if (status.iteration_limit_exceeded) {
1352 FPRINTF(stderr, "The iteration limit was exceeded while calculating "
1353 "derivatives.\nIntegration will be terminated at "
1354 "the end of the current step.");
1355 return 0;
1356 }
1357 if (status.panic) {
1358 FPRINTF(stderr, "The user patience limit was exceeded while "
1359 "calculating derivatives.\nIntegration will be terminated "
1360 "at the end of the current step.");
1361 return 0;
1362 }
1363 return 0;
1364 }

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