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

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