/[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 894 - (show annotations) (download) (as text)
Mon Oct 23 01:07:58 2006 UTC (15 years, 8 months ago) by johnpye
File MIME type: text/x-csrc
File size: 41075 byte(s)
Improved error message in integrator.c (added names of variables with clashing ode_id).
Removed debug output from observer.py.
Fixed icon paths in modelview (context menu)
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 char *varname1, *varname2;
574
575 CONSOLE_DEBUG("Starting ODE analysis");
576 IntegInitSymbols();
577
578 /* collect potential states and derivatives */
579 blsys->indepvars = gl_create(10L); /* t var info */
580 blsys->dynvars = gl_create(200L); /* y ydot var info */
581 blsys->obslist = gl_create(100L); /* obs info */
582 if (blsys->dynvars == NULL
583 || blsys->obslist == NULL
584 || blsys->indepvars == NULL
585 ){
586 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
587 return 0;
588 }
589
590 blsys->nstates = blsys->nderivs = 0;
591
592 integrator_visit_system_vars(blsys,&integrator_ode_classify_var);
593
594 integrator_print_var_stats(blsys);
595
596 if(!integrator_check_indep_var(blsys))return 0;
597
598 /* check sanity of state and var lists */
599
600 len = gl_length(blsys->dynvars);
601 half = len/2;
602 CONSOLE_DEBUG("NUMBER OF DYNAMIC VARIABLES = %ld",half);
603
604 if (len % 2 || len == 0L || blsys->nstates != blsys->nderivs ) {
605 /* list length must be even for vars to pair off */
606 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"n_y != n_ydot, or no dynamic vars found. Fix your indexing.");
607 return 0;
608 }
609 gl_sort(blsys->dynvars,(CmpFunc)Integ_CmpDynVars);
610 if (gl_fetch(blsys->dynvars,len)==NULL) {
611 ERROR_REPORTER_NOLINE(ASC_PROG_ERR,"Mysterious NULL found!");
612 return 0;
613 }
614 blsys->states = gl_create(half); /* state vars Integ_var_t references */
615 blsys->derivs = gl_create(half); /* derivative var atoms */
616 for (i=1;i < len; i+=2) {
617 v1 = (struct Integ_var_t *)gl_fetch(blsys->dynvars,i);
618 v2 = (struct Integ_var_t *)gl_fetch(blsys->dynvars,i+1);
619 if (v1->type!=1 || v2 ->type !=2 || v1->index != v2->index) {
620 varname1 = var_make_name(blsys->system,v1->i);
621 varname2 = var_make_name(blsys->system,v2->i);
622
623 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Mistyped or misindexed dynamic variables: %s (%s = %ld,%s = %ld) and %s (%s = %ld,%s = %ld).",
624 varname1, SCP(STATEFLAG),v1->type,SCP(STATEINDEX),v1->index,
625 varname2, SCP(STATEFLAG),v2->type,SCP(STATEINDEX),v2->index
626 );
627 ASC_FREE(varname1);
628 ASC_FREE(varname2);
629 happy=0;
630 break;
631 } else {
632 gl_append_ptr(blsys->states,(POINTER)v1);
633 gl_append_ptr(blsys->derivs,(POINTER)v2->i);
634 }
635 }
636 if (!happy) {
637 return 0;
638 }
639 blsys->n_y = half;
640 blsys->y = ASC_NEW_ARRAY(struct var_variable *, half);
641 blsys->y_id = ASC_NEW_ARRAY(long, half);
642 blsys->ydot = ASC_NEW_ARRAY(struct var_variable *, half);
643 if (blsys->y==NULL || blsys->ydot==NULL || blsys->y_id==NULL) {
644 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
645 return 0;
646 }
647 for (i = 1; i <= half; i++) {
648 v1 = (struct Integ_var_t *)gl_fetch(blsys->states,i);
649 blsys->y[i-1] = v1->i;
650 blsys->y_id[i-1] = v1->index;
651 blsys->ydot[i-1] = (struct var_variable *)gl_fetch(blsys->derivs,i);
652 }
653
654 if(!integrator_sort_obs_vars(blsys))return 0;
655
656 /* don't need the gl_lists now that we have arrays for everyone */
657 gl_destroy(blsys->states);
658 gl_destroy(blsys->derivs);
659 gl_free_and_destroy(blsys->indepvars); /* we own the objects in indepvars */
660 gl_free_and_destroy(blsys->dynvars); /* we own the objects in dynvars */
661 gl_free_and_destroy(blsys->obslist); /* and obslist */
662 blsys->states = NULL;
663 blsys->derivs = NULL;
664 blsys->indepvars = NULL;
665 blsys->dynvars = NULL;
666 blsys->obslist = NULL;
667
668 /* analysis completed OK */
669 return 1;
670 }
671
672 /**
673 Reindex observations. Sort if the user mostly numbered. Take natural order
674 if user just booleaned.
675
676 @return 1 on success
677 */
678 static int integrator_sort_obs_vars(IntegratorSystem *blsys){
679 int half, i, len = 0;
680 struct Integ_var_t *v2;
681
682 half = blsys->n_y;
683 len = gl_length(blsys->obslist);
684 /* we shouldn't be seeing NULL here ever except if malloc fail. */
685 if (len > 1L) {
686 half = ((struct Integ_var_t *)gl_fetch(blsys->obslist,1))->index;
687 /* half != 0 now because we didn't collect 0 indexed vars */
688 for (i=2; i <= len; i++) {
689 if (half != ((struct Integ_var_t *)gl_fetch(blsys->obslist,i))->index) {
690 /* change seen. sort and go on */
691 gl_sort(blsys->obslist,(CmpFunc)Integ_CmpObs);
692 break;
693 }
694 }
695 }
696 for (i = half = 1; i <= len; i++) {
697 v2 = (struct Integ_var_t *)gl_fetch(blsys->obslist,i);
698 if (v2==NULL) {
699 /* we shouldn't be seeing NULL here ever except if malloc fail. */
700 gl_delete(blsys->obslist,i,0); /* should not be gl_delete(so,i,1) */
701 } else {
702 Integ_SetObsId(v2->i,half);
703 v2->index = half++;
704 }
705 }
706
707 /* obslist now uniquely indexed, no nulls */
708 /* make into arrays */
709 half = gl_length(blsys->obslist);
710 blsys->obs = ASC_NEW_ARRAY(struct var_variable *,half);
711 blsys->obs_id = ASC_NEW_ARRAY(long, half);
712 if ( blsys->obs==NULL || blsys->obs_id==NULL) {
713 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
714 return 0;
715 }
716 blsys->n_obs = half;
717 for (i = 1; i <= half; i++) {
718 v2 = (struct Integ_var_t *)gl_fetch(blsys->obslist,i);
719 blsys->obs[i-1] = v2->i;
720 blsys->obs_id[i-1] = v2->index;
721 }
722
723 return 1;
724 }
725
726 static void integrator_print_var_stats(IntegratorSystem *blsys){
727 int v = gl_length(blsys->dynvars);
728 int i = gl_length(blsys->indepvars);
729 CONSOLE_DEBUG("Currently %d vars, %d indep",v,i);
730 }
731
732 /**
733 @return 1 on success
734 */
735 static int integrator_check_indep_var(IntegratorSystem *blsys){
736 int len, i;
737 struct Integ_var_t *info;
738 char *varname;
739
740 /* check the sanity of the independent variable */
741 len = gl_length(blsys->indepvars);
742 if (!len) {
743 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No independent variable found.");
744 return 0;
745 }
746 if (len > 1) {
747 ERROR_REPORTER_START_HERE(ASC_USER_ERROR);
748 FPRINTF(ASCERR,"Excess %ld independent variables found:",
749 len);
750 for(i=1; i <=len;i++) {
751 info = (struct Integ_var_t *)gl_fetch(blsys->indepvars,i);
752 if(info==NULL)continue;
753
754 varname = var_make_name(blsys->system,info->i);
755 FPRINTF(ASCERR," %s",varname);
756 ASC_FREE(varname);
757 }
758 FPRINTF(ASCERR , "\nSet the \"%s\" flag on all but one of these to %s >= 0.\n"
759 , SCP(STATEFLAG),SCP(STATEFLAG)
760 );
761 error_reporter_end_flush();
762 return 0;
763 }else{
764 info = (struct Integ_var_t *)gl_fetch(blsys->indepvars,1);
765 blsys->x = info->i;
766 }
767 return 1;
768 }
769
770 /*------------------------------------------------------------------------------
771 CLASSIFICATION OF VARIABLES (for ANALYSIS step)
772 */
773
774 #define INTEG_ADD_TO_LIST(info,TYPE,INDEX,VAR,LIST) \
775 info = ASC_NEW(struct Integ_var_t); \
776 if(info==NULL){ \
777 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory (INTEG_VAR_NEW)"); \
778 return; \
779 } \
780 info->type=TYPE; \
781 info->index=INDEX; \
782 info->i=VAR; \
783 gl_append_ptr(LIST,(void *)info); \
784 info = NULL
785
786 /**
787 In a DAE, it's either the (single) independent variable, or it's a
788 variable in the model.
789
790 I'm not sure what we should be doing with variables that are already
791 present as derivatives of other variables, I guess those ones need to be
792 removed from the list in a second pass?
793 */
794 void integrator_dae_classify_var(IntegratorSystem *blsys, struct var_variable *var){
795 struct Integ_var_t *info;
796 long type,index;
797
798 /* filter for recognition of solver_vars */
799 var_filter_t vfilt;
800 vfilt.matchbits = VAR_SVAR;
801 vfilt.matchvalue = VAR_SVAR;
802
803 assert(var != NULL && var_instance(var)!=NULL );
804
805 if( var_apply_filter(var,&vfilt) ) {
806
807 if(!var_fixed(var)){
808 /* get the ode_type and ode_id of this solver_var */
809 type = DynamicVarInfo(var,&index);
810
811 if(type==INTEG_OTHER_VAR){
812 /* if the var's type is -1, it's independent */
813 INTEG_ADD_TO_LIST(info,INTEG_OTHER_VAR,0,var,blsys->indepvars);
814 }else{
815 if(type < 0)type=0;
816 /* any other type of var is in the DAE system, at least for now */
817 INTEG_ADD_TO_LIST(info,type,index,var,blsys->dynvars);
818 }
819 }
820
821 /* if the var's obs_id > 0, add it to the observation list */
822 if(ObservationVar(var,&index) != NULL && index > 0L) {
823 INTEG_ADD_TO_LIST(info,type,index,var,blsys->obslist);
824 }
825 }
826 }
827
828 /**
829 Inspect a specific variable and work out what type it is (what 'ode_type' it
830 has) and what other variable(s) it corresponds to (ie dydt corresponds to
831 y as a derivative).
832
833 @TODO add ability to create new variables for 'missing' derivative vars?
834 */
835 void integrator_ode_classify_var(IntegratorSystem *blsys, struct var_variable *var){
836 struct Integ_var_t *info;
837 long type,index;
838
839 var_filter_t vfilt;
840 vfilt.matchbits = VAR_SVAR;
841 vfilt.matchvalue = VAR_SVAR;
842
843 assert(var != NULL && var_instance(var)!=NULL );
844
845 if( var_apply_filter(var,&vfilt) ) {
846 /* it's a solver var: what type of variable? */
847 type = DynamicVarInfo(var,&index);
848
849 if(type==INTEG_ALGEBRAIC_VAR){
850 /* no action required */
851 }else if(type==INTEG_OTHER_VAR){
852 /* i.e. independent var */
853 INTEG_ADD_TO_LIST(info,type,index,var,blsys->indepvars);
854 }else if(type>=INTEG_STATE_VAR){
855 INTEG_ADD_TO_LIST(info,type,index,var,blsys->dynvars);
856 if(type == 1){
857 blsys->nstates++;
858 }else if(type == 2){ /* what about higher-order derivatives? -- JP */
859 blsys->nderivs++;
860 }else{
861 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Higher-order (>=2) derivatives are not supported in ODEs.");
862 } }
863
864 if(ObservationVar(var,&index) != NULL && index > 0L) {
865 INTEG_ADD_TO_LIST(info,0L,index,var,blsys->obslist);
866 }
867 }
868 }
869
870 /**
871 Look at a variable and determine if it's the independent variable or not.
872 This is just for the purpose of the integrator_find_indep_var function,
873 which is a utility function provided for use by the GUI.
874 */
875 void integrator_classify_indep_var(IntegratorSystem *blsys, struct var_variable *var){
876 struct Integ_var_t *info;
877 long type,index;
878
879 var_filter_t vfilt;
880 vfilt.matchbits = VAR_SVAR;
881 vfilt.matchvalue = VAR_SVAR;
882
883 /* CONSOLE_DEBUG("..."); */
884
885 assert(var != NULL && var_instance(var)!=NULL );
886
887 if( var_apply_filter(var,&vfilt) ) {
888 type = DynamicVarInfo(var,&index);
889
890 if(type==INTEG_OTHER_VAR){
891 /* i.e. independent var */
892 INTEG_ADD_TO_LIST(info,type,index,var,blsys->indepvars);
893 }
894 }
895 }
896
897 /**
898 Look at a variable, and if it is an 'ODE variable' (it has a child instance
899 named 'ode_type') return its type, which will be either:
900 - INTEG_OTHER_VAR (if 'ode_type' is -1)
901 - INTEG_ALGEBRAIC_VAR (if 'ode_type' is zero or any negative value < -1)
902 - INTEG_STATE_VAR (if 'ode_type' is 1)
903 - values 2, 3 or up, indicating derivatives (1st deriv=2, 2nd deriv=3, etc)
904
905 If the parameter 'index' is not null, the value of 'ode_id' will be stuffed
906 there.
907 */
908 static long DynamicVarInfo(struct var_variable *v,long *index){
909 struct Instance *c, *d, *i;
910 i = var_instance(v);
911 assert(i!=NULL);
912 assert(STATEFLAG!=NULL);
913 assert(STATEINDEX!=NULL);
914 c = ChildByChar(i,STATEFLAG);
915 d = ChildByChar(i,STATEINDEX);
916 /* lazy evaluation is important in the following if */
917 if(c == NULL
918 || d == NULL
919 || InstanceKind(c) != INTEGER_INST
920 || InstanceKind(d) != INTEGER_INST
921 || !AtomAssigned(c)
922 || (!AtomAssigned(d) && GetIntegerAtomValue(c) != INTEG_OTHER_VAR)
923 ){
924 return INTEG_ALGEBRAIC_VAR;
925 }
926 if (index != NULL) {
927 *index = GetIntegerAtomValue(d);
928 }
929 return GetIntegerAtomValue(c);
930 }
931
932 /**
933 Looks at the given variable checks if it is an 'observation variable'. This
934 means that it has its 'obs_id' child instance set to a non-zero value.
935
936 If the variable is an observation variable, its index value ('obs_id') is
937 stuff into *index (provided index!=NULL), and the pointer to the original
938 instance is rtruend.
939
940 If it's not an observation variable, we return NULL and *index is untouched.
941 */
942 static struct var_variable *ObservationVar(struct var_variable *v, long *index){
943 struct Instance *c,*i;
944 i = var_instance(v);
945 assert(i!=NULL);
946 c = ChildByChar(i,OBSINDEX);
947 if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
948 return NULL;
949 }
950 if (index != NULL) {
951 *index = GetIntegerAtomValue(c);
952 }
953 return v;
954 }
955
956 /*------------------------------------------------------------------------------
957 RUNNING THE SOLVER
958 */
959
960 /*
961 Make the call to the actual integrator we've selected, for the range of
962 time values specified. The blsys contains all the specifics.
963
964 Return 1 on success
965 */
966 int integrator_solve(IntegratorSystem *blsys, long i0, long i1){
967
968 long nstep;
969 unsigned long start_index=0, finish_index=0;
970
971 assert(blsys!=NULL);
972
973 nstep = integrator_getnsamples(blsys)-1;
974 /* check for at least 2 steps and dimensionality of x vs steps here */
975
976 if (i0<0 || i1 <0) {
977 /* dude, there's no way we're writing interactive stuff here... */
978 ERROR_REPORTER_HERE(ASC_PROG_ERROR,"Console input of integration limits has been disabled!");
979 return 0;
980 }else{
981 start_index=i0;
982 finish_index =i1;
983 if (start_index >= (unsigned long)nstep) {
984 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Start point (=%lu) must be an index in the range [0,%li]."
985 ,start_index,nstep
986 );
987 return 0;
988 }
989 if (finish_index > (unsigned long)nstep) {
990 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"End point (=%lu) must be an index in the range [0,%li]."
991 ,finish_index,nstep
992 );
993 return 0;
994 }
995 }
996
997 if(finish_index <= start_index) {
998 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"End point comes before start point! (start=%lu, end=%lu)"
999 ,start_index,finish_index
1000 );
1001 return 0;
1002 }
1003
1004 CONSOLE_DEBUG("RUNNING INTEGRATION...");
1005
1006 /* now go and run the integrator */
1007 switch (blsys->engine) {
1008 case INTEG_LSODE: return integrator_lsode_solve(blsys, start_index, finish_index); break;
1009 #ifdef ASC_WITH_IDA
1010 case INTEG_IDA: return integrator_ida_solve(blsys,start_index, finish_index); break;
1011 #endif
1012 default:
1013 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown integrator (invalid, or not implemented yet)");
1014 return 0;
1015 }
1016 }
1017
1018 /*---------------------------------------------------------------
1019 HANDLING THE LIST OF TIMESTEMPS
1020 */
1021
1022 #define GETTER_AND_SETTER(TYPE,NAME) \
1023 void integrator_set_##NAME(IntegratorSystem *blsys, TYPE val){ \
1024 blsys->NAME=val; \
1025 } \
1026 TYPE integrator_get_##NAME(IntegratorSystem *blsys){ \
1027 return blsys->NAME; \
1028 }
1029
1030 GETTER_AND_SETTER(SampleList *,samples);
1031 GETTER_AND_SETTER(double,maxstep);
1032 GETTER_AND_SETTER(double,minstep);
1033 GETTER_AND_SETTER(double,stepzero);
1034 GETTER_AND_SETTER(int,maxsubsteps);
1035 #undef GETTER_AND_SETTER
1036
1037 long integrator_getnsamples(IntegratorSystem *blsys){
1038 assert(blsys!=NULL);
1039 assert(blsys->samples!=NULL);
1040 return samplelist_length(blsys->samples);
1041 }
1042
1043 double integrator_getsample(IntegratorSystem *blsys, long i){
1044 assert(blsys!=NULL);
1045 assert(blsys->samples!=NULL);
1046 return samplelist_get(blsys->samples,i);
1047 }
1048
1049 void integrator_setsample(IntegratorSystem *blsys, long i,double xi){
1050 assert(blsys!=NULL);
1051 assert(blsys->samples!=NULL);
1052 samplelist_set(blsys->samples,i,xi);
1053 }
1054
1055 const dim_type *integrator_getsampledim(IntegratorSystem *blsys){
1056 assert(blsys!=NULL);
1057 assert(blsys->samples!=NULL);
1058 return samplelist_dim(blsys->samples);
1059 }
1060
1061 ASC_DLLSPEC(long) integrator_getcurrentstep(IntegratorSystem *blsys){
1062 return blsys->currentstep;
1063 }
1064
1065 /*------------------------------------------------------------------------------
1066 GET/SET VALUE OF THE INDEP VARIABLE
1067 */
1068
1069 /**
1070 Retrieve the value of the independent variable (time) from ASCEND
1071 and return it as a double.
1072 */
1073 double integrator_get_t(IntegratorSystem *blsys){
1074 assert(blsys->x!=NULL);
1075 return var_value(blsys->x);
1076 }
1077
1078 /**
1079 Set the value of the independent variable (time) in ASCEND.
1080 */
1081 void integrator_set_t(IntegratorSystem *blsys, double value){
1082 var_set_value(blsys->x, value);
1083 /* CONSOLE_DEBUG("set_t = %g", value); */
1084 }
1085
1086 /*------------------------------------------------------------------------------
1087 PASSING DIFFERENTIAL VARIABLES AND THEIR DERIVATIVES TO/FROM THE SOLVER
1088 */
1089 /**
1090 Retrieve the current values of the derivatives of the y-variables
1091 and stick them in the/an array that the integrator will use.
1092
1093 If the pointer 'y' is NULL, the necessary space is allocated (and
1094 must be freed somewhere else).
1095 */
1096 double *integrator_get_y(IntegratorSystem *blsys, double *y) {
1097 long i;
1098
1099 if (y==NULL) {
1100 y = ASC_NEW_ARRAY_CLEAR(double, blsys->n_y+1);
1101 /* C y[0] <==> ascend d.y[1] <==> f77 y(1) */
1102 }
1103
1104 for (i=0; i< blsys->n_y; i++) {
1105 assert(blsys->y[i]!=NULL);
1106 y[i] = var_value(blsys->y[i]);
1107 /* CONSOLE_DEBUG("ASCEND --> y[%ld] = %g", i+1, y[i]); */
1108 }
1109 return y;
1110 }
1111
1112 /**
1113 Take the values of the derivatives from the array that the integrator
1114 uses, and use them to update the values of the corresponding variables
1115 in ASCEND.
1116 */
1117 void integrator_set_y(IntegratorSystem *blsys, double *y) {
1118 long i;
1119 #ifndef NDEBUG
1120 char *varname;
1121 #endif
1122
1123 for (i=0; i < blsys->n_y; i++) {
1124 assert(blsys->y[i]!=NULL);
1125 var_set_value(blsys->y[i],y[i]);
1126 #ifndef NDEBUG
1127 varname = var_make_name(blsys->system, blsys->y[i]);
1128 /* CONSOLE_DEBUG("y[%ld] = \"%s\" = %g --> ASCEND", i+1, varname, y[i]); */
1129 ASC_FREE(varname);
1130 #endif
1131 }
1132 }
1133
1134 /**
1135 Send the values of the derivatives of the 'y' variables to the solver.
1136 Allocate space for an array if necessary.
1137
1138 Any element in blsys->ydot that is NULL will be passed over (the value
1139 won't be modified in dydx).
1140 */
1141 double *integrator_get_ydot(IntegratorSystem *blsys, double *dydx) {
1142 long i;
1143
1144 if (dydx==NULL) {
1145 dydx = ASC_NEW_ARRAY_CLEAR(double, blsys->n_y+1);
1146 /* C dydx[0] <==> ascend d.dydx[1] <==> f77 ydot(1) */
1147 }
1148
1149 for (i=0; i < blsys->n_y; i++) {
1150 if(blsys->ydot[i]!=NULL){
1151 dydx[i] = var_value(blsys->ydot[i]);
1152 }
1153 /* CONSOLE_DEBUG("ASCEND --> ydot[%ld] = %g", i+1, dydx[i]); */
1154 }
1155 return dydx;
1156 }
1157
1158 void integrator_set_ydot(IntegratorSystem *blsys, double *dydx) {
1159 long i;
1160 #ifndef NDEBUG
1161 char *varname;
1162 #endif
1163 for (i=0; i < blsys->n_y; i++) {
1164 if(blsys->ydot[i]!=NULL){
1165 var_set_value(blsys->ydot[i],dydx[i]);
1166 #ifndef NDEBUG
1167 varname = var_make_name(blsys->system, blsys->ydot[i]);
1168 CONSOLE_DEBUG("ydot[%ld] = \"%s\" = %g --> ASCEND", i+1, varname, dydx[i]);
1169 ASC_FREE(varname);
1170 #endif
1171 }
1172 #ifndef NDEBUG
1173 else{
1174 CONSOLE_DEBUG("ydot[%ld] = %g (internal)", i+1, dydx[i]);
1175 }
1176 #endif
1177 }
1178 }
1179
1180 /*-------------------------------------------------------------
1181 RETRIEVING OBSERVATION DATA
1182
1183 This function takes the inst in the solver and returns the vector of
1184 observation variables that are located in the submodel d.obs array.
1185 */
1186 double *integrator_get_observations(IntegratorSystem *blsys, double *obsi) {
1187 long i;
1188
1189 if (obsi==NULL) {
1190 obsi = ASC_NEW_ARRAY_CLEAR(double, blsys->n_obs+1);
1191 }
1192
1193 /* C obsi[0] <==> ascend d.obs[1] */
1194
1195 for (i=0; i < blsys->n_obs; i++) {
1196 obsi[i] = var_value(blsys->obs[i]);
1197 /* CONSOLE_DEBUG("*get_d_obs[%ld] = %g\n", i+1, obsi[i]); */
1198 }
1199 return obsi;
1200 }
1201
1202 struct var_variable *integrator_get_observed_var(IntegratorSystem *blsys, const long i){
1203 assert(i>=0);
1204 assert(i<blsys->n_obs);
1205 return blsys->obs[i];
1206 }
1207
1208 /**
1209 @NOTE Although this shouldn't be required for implementation of solver
1210 engines, this is useful for GUI reporting of integration results.
1211 */
1212 struct var_variable *integrator_get_independent_var(IntegratorSystem *blsys){
1213 return blsys->x;
1214 }
1215
1216
1217 /*----------------------------------------------------
1218 Build an analytic jacobian for solving the state system
1219
1220 This necessarily ugly piece of code attempts to create a unique
1221 list of relations that explicitly contain the variables in the
1222 given input list. The utility of this information is that we know
1223 exactly which relations must be differentiated, to fill in the
1224 df/dy matrix. If the problem has very few derivative terms, this will
1225 be of great savings. If the problem arose from the discretization of
1226 a pde, then this will be not so useful. The decision wether to use
1227 this function or to simply differentiate the entire relations list
1228 must be done before calling this function.
1229
1230 Final Note: the callee owns the array, but not the array elements.
1231 */
1232 #define AVG_NUM_INCIDENT 4
1233
1234
1235 /**
1236 This function helps to arrange the observation variables in a sensible order.
1237 The 'obs_id' child instance of v, if present, is assigned the value of the
1238 given parameter 'index'.
1239 */
1240 static void Integ_SetObsId(struct var_variable *v, long index){
1241 struct Instance *c, *i;
1242 i = var_instance(v);
1243 assert(i!=NULL);
1244 c = ChildByChar(i,OBSINDEX);
1245 if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
1246 return;
1247 }
1248 SetIntegerAtomValue(c,index,0);
1249 }
1250
1251 /**
1252 Compares observation structs. NULLs should end up at far end.
1253 */
1254 static int Integ_CmpObs(struct Integ_var_t *v1, struct Integ_var_t *v2){
1255 if(v1 == NULL)return 1;
1256 if(v2 == NULL)return -1;
1257 if(v1->index > v2->index)return 1;
1258 if(v1->index == v2->index)return 0;
1259 return -1;
1260 }
1261
1262 /**
1263 Compares dynamic vars structs. NULLs should end up at far end.
1264 List should be sorted primarily by index and then by type, in order
1265 of increasing value of both.
1266 */
1267 static int Integ_CmpDynVars(struct Integ_var_t *v1, struct Integ_var_t *v2){
1268 if(v1 == NULL)return 1;
1269 if(v2 == NULL)return -1;
1270 if(v1->index > v2->index)return 1;
1271 if(v1->index != v2->index)return -1;
1272 if(v1->type > v2->type)return 1;
1273 return -1;
1274 }
1275 /*----------------------------
1276 Output handling to the GUI/interface.
1277 */
1278
1279 int integrator_set_reporter(IntegratorSystem *blsys
1280 , IntegratorReporter *reporter
1281 ){
1282 assert(blsys!=NULL);
1283 blsys->reporter = reporter;
1284 /* ERROR_REPORTER_HERE(ASC_PROG_NOTE,"INTEGRATOR REPORTER HOOKS HAVE BEEN SET\n"); */
1285 return 1;
1286 }
1287
1288 int integrator_output_init(IntegratorSystem *blsys){
1289 assert(blsys!=NULL);
1290 assert(blsys->reporter!=NULL);
1291 if(blsys->reporter->init!=NULL){
1292 /* call the specified output function */
1293 return (*(blsys->reporter->init))(blsys);
1294 }
1295 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter init method");
1296 return 1;
1297 }
1298
1299 int integrator_output_write(IntegratorSystem *blsys){
1300 static int reported_already=0;
1301 assert(blsys!=NULL);
1302 if(blsys->reporter->write!=NULL){
1303 return (*(blsys->reporter->write))(blsys);
1304 }
1305 if(!reported_already){
1306 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter write method (this message only shown once)");
1307 reported_already=1;
1308 }
1309 return 1;
1310 }
1311
1312 int integrator_output_write_obs(IntegratorSystem *blsys){
1313 static int reported_already=0;
1314 assert(blsys!=NULL);
1315 if(blsys->reporter->write_obs!=NULL){
1316 return (*(blsys->reporter->write_obs))(blsys);
1317 }
1318 if(!reported_already){
1319 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter write_obs method (this message only shown once)");
1320 reported_already=1;
1321 }
1322 return 1;
1323 }
1324
1325 int integrator_output_close(IntegratorSystem *blsys){
1326 assert(blsys!=NULL);
1327 if(blsys->reporter->close!=NULL){
1328 return (*(blsys->reporter->close))(blsys);
1329 }
1330 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter close method");
1331 return 1;
1332 }
1333
1334 /**
1335 Decode status codes from the integrator, and output them via FPRINTF.
1336 */
1337 int integrator_checkstatus(slv_status_t status) {
1338 if (status.converged) {
1339 return 1;
1340 }
1341 if (status.diverged) {
1342 FPRINTF(stderr, "The derivative system did not converge.\nIntegration "
1343 "will be terminated at the end of the current step.");
1344 return 0;
1345 }
1346 if (status.inconsistent) {
1347 FPRINTF(stderr, "A numerical inconsistency was discovered while "
1348 "calculating derivatives. Integration will be terminated at "
1349 "the end of the current step.");
1350 return 0;
1351 }
1352 if (status.time_limit_exceeded) {
1353 FPRINTF(stderr, "The time limit was exceeded while calculating "
1354 "derivatives.\nIntegration will be terminated at the end of "
1355 "the current step.");
1356 return 0;
1357 }
1358 if (status.iteration_limit_exceeded) {
1359 FPRINTF(stderr, "The iteration limit was exceeded while calculating "
1360 "derivatives.\nIntegration will be terminated at "
1361 "the end of the current step.");
1362 return 0;
1363 }
1364 if (status.panic) {
1365 FPRINTF(stderr, "The user patience limit was exceeded while "
1366 "calculating derivatives.\nIntegration will be terminated "
1367 "at the end of the current step.");
1368 return 0;
1369 }
1370 return 0;
1371 }

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