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

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