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

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