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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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