/[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 979 - (show annotations) (download) (as text)
Wed Dec 20 14:34:16 2006 UTC (13 years, 1 month ago) by johnpye
File MIME type: text/x-csrc
File size: 47300 byte(s)
Added simplified ASC_PANIC call that uses var-args, added throughout relation_util.c.
Fixed var_filter_t stuff in djex and fvex.
More assertions in integrator.c
Added output of initial state from lsode.c (hoping that's a good idea?)
Fixed output code from relman_diff2.
Added asc_panic_nofunc for non var-arg CPPs.
Disabled -O3 flag in building C++ API
Added __getitem__ and __getattr__ methods in Simuluation for simplified python syntax (eg M.x instead M.sim.x)
Integrator::analyse throws exceptions on error now.

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\tname\ty_id\ty\tydot\n");
666 fprintf(stderr,"-----\t-----\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 // maybe it's the independent variable
684 if(var==sys->x){
685 if(varindx==NULL){
686 fprintf(stderr,".\t%s\t(indep)\n",varname);
687 }else{
688 fprintf(stderr,"%d\t%s\t(index)\n",*varindx,varname);
689 }
690 ASC_FREE(varname);
691 return;
692 }
693
694 // if it's fixed then it's not really a DAE var
695 if(var_fixed(var)){
696 fprintf(stderr,"%d\t%s\t(fixed)\n",*varindx,varname);
697 ASC_FREE(varname);
698 return;
699 }
700
701 // if there's no varindx, and neither of the above, then there's a problem
702 if(varindx==NULL){
703 fprintf(stderr,".\t%s\t(not visited?)\t",varname);
704 ASC_FREE(varname);
705 return;
706 }
707
708 y_id = sys->y_id[*varindx];
709
710 fprintf(stderr,"%d\t%s\t%ld", *varindx, varname,y_id);
711 ASC_FREE(varname);
712 if(y_id >= 0){
713 fprintf(stderr,"\ty[%ld]\t.\n",y_id);
714 }else{
715 fprintf(stderr,"\t.\tydot[%ld]\n",-y_id-1);
716 }
717
718 ASC_ASSERT_LT(*varindx,1e7L);
719 ASC_ASSERT_LT(y_id, 9999999L);
720 ASC_ASSERT_LT(-9999999L, y_id);
721 }
722
723 void integrator_visit_system_vars(IntegratorSystem *sys,IntegratorVarVisitorFn *visitfn){
724 struct var_variable **vlist;
725 int i, vlen;
726
727 /* visit all the slv_system_t master var lists to collect vars */
728 /* find the vars mostly in this one */
729 vlist = slv_get_solvers_var_list(sys->system);
730 vlen = slv_get_num_solvers_vars(sys->system);
731 for (i=0;i<vlen;i++) {
732 (*visitfn)(sys, vlist[i], &i);
733 }
734
735 /*
736 CONSOLE_DEBUG("Checked %d vars",vlen);
737 integrator_print_var_stats(sys);
738 */
739
740 /* probably nothing here, but gotta check. */
741 vlist = slv_get_master_par_list(sys->system);
742 vlen = slv_get_num_master_pars(sys->system);
743 for (i=0;i<vlen;i++) {
744 (*visitfn)(sys, vlist[i], NULL);
745 }
746
747 /*
748 CONSOLE_DEBUG("Checked %d pars",vlen);
749 integrator_print_var_stats(sys);
750 */
751
752 /* might find t here */
753 vlist = slv_get_master_unattached_list(sys->system);
754 vlen = slv_get_num_master_unattached(sys->system);
755 for (i=0;i<vlen;i++) {
756 (*visitfn)(sys, vlist[i], NULL);
757 }
758
759 /* CONSOLE_DEBUG("Checked %d unattached",vlen); */
760 }
761 /**
762 @return 1 on success
763 */
764 int integrator_analyse_ode(IntegratorSystem *sys){
765 struct Integ_var_t *v1,*v2;
766 long half,i,len;
767 int happy=1;
768 char *varname1, *varname2;
769
770 CONSOLE_DEBUG("Starting ODE analysis");
771 IntegInitSymbols();
772
773 /* collect potential states and derivatives */
774 sys->indepvars = gl_create(10L); /* t var info */
775 sys->dynvars = gl_create(200L); /* y ydot var info */
776 sys->obslist = gl_create(100L); /* obs info */
777 if (sys->dynvars == NULL
778 || sys->obslist == NULL
779 || sys->indepvars == NULL
780 ){
781 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
782 return 0;
783 }
784
785 sys->nstates = sys->nderivs = 0;
786
787 integrator_visit_system_vars(sys,&integrator_ode_classify_var);
788
789 integrator_print_var_stats(sys);
790
791 if(!integrator_check_indep_var(sys))return 0;
792
793 /* check sanity of state and var lists */
794
795 len = gl_length(sys->dynvars);
796 half = len/2;
797 CONSOLE_DEBUG("NUMBER OF DYNAMIC VARIABLES = %ld",half);
798
799 if (len % 2 || len == 0L || sys->nstates != sys->nderivs ) {
800 /* list length must be even for vars to pair off */
801 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"n_y != n_ydot, or no dynamic vars found. Fix your indexing.");
802 return 0;
803 }
804 gl_sort(sys->dynvars,(CmpFunc)Integ_CmpDynVars);
805 if (gl_fetch(sys->dynvars,len)==NULL) {
806 ERROR_REPORTER_NOLINE(ASC_PROG_ERR,"Mysterious NULL found!");
807 return 0;
808 }
809 sys->states = gl_create(half); /* state vars Integ_var_t references */
810 sys->derivs = gl_create(half); /* derivative var atoms */
811 for (i=1;i < len; i+=2) {
812 v1 = (struct Integ_var_t *)gl_fetch(sys->dynvars,i);
813 v2 = (struct Integ_var_t *)gl_fetch(sys->dynvars,i+1);
814 if (v1->type!=1 || v2 ->type !=2 || v1->index != v2->index) {
815 varname1 = var_make_name(sys->system,v1->i);
816 varname2 = var_make_name(sys->system,v2->i);
817
818 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Mistyped or misindexed dynamic variables: %s (%s = %ld,%s = %ld) and %s (%s = %ld,%s = %ld).",
819 varname1, SCP(STATEFLAG),v1->type,SCP(STATEINDEX),v1->index,
820 varname2, SCP(STATEFLAG),v2->type,SCP(STATEINDEX),v2->index
821 );
822 ASC_FREE(varname1);
823 ASC_FREE(varname2);
824 happy=0;
825 break;
826 } else {
827 gl_append_ptr(sys->states,(POINTER)v1);
828 gl_append_ptr(sys->derivs,(POINTER)v2->i);
829 }
830 }
831 if (!happy) {
832 return 0;
833 }
834 sys->n_y = half;
835 sys->y = ASC_NEW_ARRAY(struct var_variable *, half);
836 sys->y_id = ASC_NEW_ARRAY(long, half);
837 sys->ydot = ASC_NEW_ARRAY(struct var_variable *, half);
838 if (sys->y==NULL || sys->ydot==NULL || sys->y_id==NULL) {
839 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
840 return 0;
841 }
842 for (i = 1; i <= half; i++) {
843 v1 = (struct Integ_var_t *)gl_fetch(sys->states,i);
844 sys->y[i-1] = v1->i;
845 sys->y_id[i-1] = v1->index;
846 sys->ydot[i-1] = (struct var_variable *)gl_fetch(sys->derivs,i);
847 }
848
849 if(!integrator_sort_obs_vars(sys))return 0;
850
851 /* FIX all states */
852 for(i=0; i<sys->n_y; ++i){
853 if(!var_fixed(sys->y[i])){
854 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Fixing state %d",i);
855 var_set_fixed(sys->y[i], TRUE);
856 }
857 }
858
859 /* don't need the gl_lists now that we have arrays for everyone */
860 gl_destroy(sys->states);
861 gl_destroy(sys->derivs);
862 gl_free_and_destroy(sys->indepvars); /* we own the objects in indepvars */
863 gl_free_and_destroy(sys->dynvars); /* we own the objects in dynvars */
864 gl_free_and_destroy(sys->obslist); /* and obslist */
865 sys->states = NULL;
866 sys->derivs = NULL;
867 sys->indepvars = NULL;
868 sys->dynvars = NULL;
869 sys->obslist = NULL;
870
871 /* analysis completed OK */
872 return 1;
873 }
874
875 /**
876 Reindex observations. Sort if the user mostly numbered. Take natural order
877 if user just booleaned.
878
879 @return 1 on success
880 */
881 static int integrator_sort_obs_vars(IntegratorSystem *sys){
882 int half, i, len = 0;
883 struct Integ_var_t *v2;
884
885 half = sys->n_y;
886 len = gl_length(sys->obslist);
887 /* we shouldn't be seeing NULL here ever except if malloc fail. */
888 if (len > 1L) {
889 half = ((struct Integ_var_t *)gl_fetch(sys->obslist,1))->index;
890 /* half != 0 now because we didn't collect 0 indexed vars */
891 for (i=2; i <= len; i++) {
892 if (half != ((struct Integ_var_t *)gl_fetch(sys->obslist,i))->index) {
893 /* change seen. sort and go on */
894 gl_sort(sys->obslist,(CmpFunc)Integ_CmpObs);
895 break;
896 }
897 }
898 }
899 for (i = half = 1; i <= len; i++) {
900 v2 = (struct Integ_var_t *)gl_fetch(sys->obslist,i);
901 if (v2==NULL) {
902 /* we shouldn't be seeing NULL here ever except if malloc fail. */
903 gl_delete(sys->obslist,i,0); /* should not be gl_delete(so,i,1) */
904 } else {
905 Integ_SetObsId(v2->i,half);
906 v2->index = half++;
907 }
908 }
909
910 /* obslist now uniquely indexed, no nulls */
911 /* make into arrays */
912 half = gl_length(sys->obslist);
913 sys->obs = ASC_NEW_ARRAY(struct var_variable *,half);
914 sys->obs_id = ASC_NEW_ARRAY(long, half);
915 if ( sys->obs==NULL || sys->obs_id==NULL) {
916 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory.");
917 return 0;
918 }
919 sys->n_obs = half;
920 for (i = 1; i <= half; i++) {
921 v2 = (struct Integ_var_t *)gl_fetch(sys->obslist,i);
922 sys->obs[i-1] = v2->i;
923 sys->obs_id[i-1] = v2->index;
924 }
925
926 return 1;
927 }
928
929 static void integrator_print_var_stats(IntegratorSystem *sys){
930 int v = gl_length(sys->dynvars);
931 int i = gl_length(sys->indepvars);
932 CONSOLE_DEBUG("Currently %d vars, %d indep",v,i);
933 }
934
935 /**
936 Check sanity of the independent variable.
937
938 @return 1 on success
939 */
940 static int integrator_check_indep_var(IntegratorSystem *sys){
941 int len, i;
942 struct Integ_var_t *info;
943 char *varname;
944
945 if(sys->x){
946 CONSOLE_DEBUG("Indep var already assigned");
947 return 1;
948 }
949
950 /* check the sanity of the independent variable */
951 len = gl_length(sys->indepvars);
952 if (!len) {
953 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No independent variable found.");
954 return 0;
955 }
956 if (len > 1) {
957 ERROR_REPORTER_START_HERE(ASC_USER_ERROR);
958 FPRINTF(ASCERR,"Excess %ld independent variables found:",
959 len);
960 for(i=1; i <=len;i++) {
961 info = (struct Integ_var_t *)gl_fetch(sys->indepvars,i);
962 if(info==NULL)continue;
963
964 varname = var_make_name(sys->system,info->i);
965 FPRINTF(ASCERR," %s",varname);
966 ASC_FREE(varname);
967 }
968 FPRINTF(ASCERR , "\nSet the \"%s\" flag on all but one of these to %s >= 0.\n"
969 , SCP(STATEFLAG),SCP(STATEFLAG)
970 );
971 error_reporter_end_flush();
972 return 0;
973 }else{
974 info = (struct Integ_var_t *)gl_fetch(sys->indepvars,1);
975 sys->x = info->i;
976 }
977 asc_assert(gl_length(sys->indepvars)==1);
978 asc_assert(sys->x);
979 return 1;
980 }
981
982 /*------------------------------------------------------------------------------
983 CLASSIFICATION OF VARIABLES (for ANALYSIS step)
984 */
985
986 #define INTEG_ADD_TO_LIST(info,TYPE,INDEX,VAR,VARINDX,LIST) \
987 info = ASC_NEW(struct Integ_var_t); \
988 if(info==NULL){ \
989 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Insufficient memory (INTEG_VAR_NEW)"); \
990 return; \
991 } \
992 info->type=TYPE; \
993 info->index=INDEX; \
994 info->i=VAR; \
995 info->derivative=NULL; \
996 info->derivative_of=NULL; \
997 if(VARINDX==NULL){ \
998 info->varindx = -1; \
999 }else{ \
1000 info->varindx = *VARINDX; \
1001 } \
1002 gl_append_ptr(LIST,(void *)info); \
1003 info = NULL
1004
1005 /**
1006 In a DAE, it's either the (single) independent variable, or it's a
1007 variable in the model.
1008
1009 I'm not sure what we should be doing with variables that are already
1010 present as derivatives of other variables, I guess those ones need to be
1011 removed from the list in a second pass?
1012 */
1013 void integrator_dae_classify_var(IntegratorSystem *sys
1014 , struct var_variable *var, const int *varindx
1015 ){
1016 struct Integ_var_t *info;
1017 long type,index;
1018
1019 /* filter for recognition of solver_vars */
1020 var_filter_t vfilt;
1021 vfilt.matchbits = VAR_SVAR;
1022 vfilt.matchvalue = VAR_SVAR;
1023
1024 assert(var != NULL && var_instance(var)!=NULL );
1025
1026 if( var_apply_filter(var,&vfilt) ) {
1027 if(!var_active(var)){
1028 CONSOLE_DEBUG("VARIABLE IS NOT ACTIVE");
1029 return;
1030 }
1031
1032 /* only non-fixed variables are accepted */
1033 if(!var_fixed(var)){
1034 /* get the ode_type and ode_id of this solver_var */
1035 type = DynamicVarInfo(var,&index);
1036
1037 if(type==INTEG_OTHER_VAR){
1038 /* if the var's type is -1, it's independent */
1039 INTEG_ADD_TO_LIST(info,INTEG_OTHER_VAR,0,var,varindx,sys->indepvars);
1040 }else{
1041 if(type < 0)type=0;
1042 /* any other type of var is in the DAE system, at least for now */
1043 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
1044 }
1045 }
1046 #if 0
1047 else{
1048 /* fixed variable, only include it if ode_type == 1 */
1049 type = DynamicVarInfo(var,&index);
1050 if(type==INTEG_STATE_VAR){
1051 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
1052 }
1053 }
1054 #endif
1055
1056 /* if the var's obs_id > 0, add it to the observation list */
1057 if(ObservationVar(var,&index) != NULL && index > 0L) {
1058 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->obslist);
1059 }
1060 }
1061 }
1062
1063 /**
1064 Inspect a specific variable and work out what type it is (what 'ode_type' it
1065 has) and what other variable(s) it corresponds to (ie dydt corresponds to
1066 y as a derivative).
1067
1068 @TODO add ability to create new variables for 'missing' derivative vars?
1069 */
1070 void integrator_ode_classify_var(IntegratorSystem *sys, struct var_variable *var
1071 , const int *varindx
1072 ){
1073 struct Integ_var_t *info;
1074 long type,index;
1075
1076 var_filter_t vfilt;
1077 vfilt.matchbits = VAR_SVAR;
1078 vfilt.matchvalue = VAR_SVAR;
1079
1080 assert(var != NULL && var_instance(var)!=NULL );
1081
1082 if( var_apply_filter(var,&vfilt) ) {
1083 /* it's a solver var: what type of variable? */
1084 type = DynamicVarInfo(var,&index);
1085
1086 if(type==INTEG_ALGEBRAIC_VAR){
1087 /* no action required */
1088 }else if(type==INTEG_OTHER_VAR){
1089 /* i.e. independent var */
1090 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->indepvars);
1091 }else if(type>=INTEG_STATE_VAR){
1092 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->dynvars);
1093 if(type == 1){
1094 sys->nstates++;
1095 }else if(type == 2){ /* what about higher-order derivatives? -- JP */
1096 sys->nderivs++;
1097 }else{
1098 ERROR_REPORTER_HERE(ASC_USER_WARNING,"Higher-order (>=2) derivatives are not supported in ODEs.");
1099 } }
1100
1101 if(ObservationVar(var,&index) != NULL && index > 0L) {
1102 INTEG_ADD_TO_LIST(info,0L,index,var,varindx,sys->obslist);
1103 }
1104 }
1105 }
1106
1107 /**
1108 Look at a variable and determine if it's the independent variable or not.
1109 This is just for the purpose of the integrator_find_indep_var function,
1110 which is a utility function provided for use by the GUI.
1111 */
1112 void integrator_classify_indep_var(IntegratorSystem *sys
1113 , struct var_variable *var, const int *varindx
1114 ){
1115 struct Integ_var_t *info;
1116 long type,index;
1117 var_filter_t vfilt;
1118 #ifdef CLASSIFY_DEBUG
1119 char *varname;
1120 #endif
1121
1122 assert(var != NULL && var_instance(var)!=NULL );
1123 vfilt.matchbits = VAR_SVAR;
1124 vfilt.matchvalue = VAR_SVAR;
1125
1126 #ifdef CLASSIFY_DEBUG
1127 varname = var_make_name(sys->system,var);
1128 #endif
1129
1130 if( var_apply_filter(var,&vfilt) ) {
1131 type = DynamicVarInfo(var,&index);
1132
1133 if(type==INTEG_OTHER_VAR){
1134 /* i.e. independent var */
1135 #ifdef CLASSIFY_DEBUG
1136 CONSOLE_DEBUG("Var '%s' added to indepvars",varname);
1137 #endif
1138 INTEG_ADD_TO_LIST(info,type,index,var,varindx,sys->indepvars);
1139 #ifdef CLASSIFY_DEBUG
1140 }else{
1141 CONSOLE_DEBUG("Var '%s' not correct ode_type",varname);
1142 #endif
1143 }
1144 #ifdef CLASSIFY_DEBUG
1145 }else{
1146 CONSOLE_DEBUG("Var '%s' failed filter",varname);
1147 #endif
1148 }
1149
1150 #ifdef CLASSIFY_DEBUG
1151 ASC_FREE(varname);
1152 #endif
1153 }
1154
1155 /**
1156 Look at a variable, and if it is an 'ODE variable' (it has a child instance
1157 named 'ode_type') return its type, which will be either:
1158 - INTEG_OTHER_VAR (if 'ode_type' is -1)
1159 - INTEG_ALGEBRAIC_VAR (if 'ode_type' is zero or any negative value < -1)
1160 - INTEG_STATE_VAR (if 'ode_type' is 1)
1161 - values 2, 3 or up, indicating derivatives (1st deriv=2, 2nd deriv=3, etc)
1162
1163 If the parameter 'index' is not null, the value of 'ode_id' will be stuffed
1164 there.
1165 */
1166 static long DynamicVarInfo(struct var_variable *v,long *index){
1167 struct Instance *c, *d, *i;
1168 i = var_instance(v);
1169 assert(i!=NULL);
1170 assert(STATEFLAG!=NULL);
1171 assert(STATEINDEX!=NULL);
1172 c = ChildByChar(i,STATEFLAG);
1173 d = ChildByChar(i,STATEINDEX);
1174 /* lazy evaluation is important in the following if */
1175 if(c == NULL
1176 || d == NULL
1177 || InstanceKind(c) != INTEGER_INST
1178 || InstanceKind(d) != INTEGER_INST
1179 || !AtomAssigned(c)
1180 || (!AtomAssigned(d) && GetIntegerAtomValue(c) != INTEG_OTHER_VAR)
1181 ){
1182 return INTEG_ALGEBRAIC_VAR;
1183 }
1184 if (index != NULL) {
1185 *index = GetIntegerAtomValue(d);
1186 }
1187 return GetIntegerAtomValue(c);
1188 }
1189
1190 /**
1191 Looks at the given variable checks if it is an 'observation variable'. This
1192 means that it has its 'obs_id' child instance set to a non-zero value.
1193
1194 If the variable is an observation variable, its index value ('obs_id') is
1195 stuff into *index (provided index!=NULL), and the pointer to the original
1196 instance is rtruend.
1197
1198 If it's not an observation variable, we return NULL and *index is untouched.
1199 */
1200 static struct var_variable *ObservationVar(struct var_variable *v, long *index){
1201 struct Instance *c,*i;
1202 i = var_instance(v);
1203 assert(i!=NULL);
1204 c = ChildByChar(i,OBSINDEX);
1205 if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
1206 return NULL;
1207 }
1208 if (index != NULL) {
1209 *index = GetIntegerAtomValue(c);
1210 }
1211 return v;
1212 }
1213
1214 /*------------------------------------------------------------------------------
1215 RUNNING THE SOLVER
1216 */
1217
1218 /*
1219 Make the call to the actual integrator we've selected, for the range of
1220 time values specified. The sys contains all the specifics.
1221
1222 Return 1 on success
1223 */
1224 int integrator_solve(IntegratorSystem *sys, long i0, long i1){
1225
1226 long nstep;
1227 unsigned long start_index=0, finish_index=0;
1228 assert(sys!=NULL);
1229
1230 assert(sys->internals);
1231 assert(sys->engine!=INTEG_UNKNOWN);
1232
1233 nstep = integrator_getnsamples(sys)-1;
1234 /* check for at least 2 steps and dimensionality of x vs steps here */
1235
1236 if (i0<0 || i1 <0) {
1237 /* dude, there's no way we're writing interactive stuff here... */
1238 ERROR_REPORTER_HERE(ASC_PROG_ERROR,"Console input of integration limits has been disabled!");
1239 return 0;
1240 } else {
1241 start_index=i0;
1242 finish_index =i1;
1243 if (start_index >= (unsigned long)nstep) {
1244 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Start point (=%lu) must be an index in the range [0,%li]."
1245 ,start_index,nstep
1246 );
1247 return 0;
1248 }
1249 if (finish_index > (unsigned long)nstep) {
1250 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"End point (=%lu) must be an index in the range [0,%li]."
1251 ,finish_index,nstep
1252 );
1253 return 0;
1254 }
1255 }
1256
1257 if(finish_index <= start_index) {
1258 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"End point comes before start point! (start=%lu, end=%lu)"
1259 ,start_index,finish_index
1260 );
1261 return 0;
1262 }
1263
1264 CONSOLE_DEBUG("RUNNING INTEGRATION...");
1265
1266 return (sys->internals->solvefn)(sys,start_index,finish_index);
1267 }
1268
1269 /*---------------------------------------------------------------
1270 HANDLING THE LIST OF TIMESTEMPS
1271 */
1272
1273 #define GETTER_AND_SETTER(TYPE,NAME) \
1274 void integrator_set_##NAME(IntegratorSystem *sys, TYPE val){ \
1275 sys->NAME=val; \
1276 } \
1277 TYPE integrator_get_##NAME(IntegratorSystem *sys){ \
1278 return sys->NAME; \
1279 }
1280
1281 GETTER_AND_SETTER(SampleList *,samples) /*;*/
1282 GETTER_AND_SETTER(double,maxstep) /*;*/
1283 GETTER_AND_SETTER(double,minstep) /*;*/
1284 GETTER_AND_SETTER(double,stepzero) /*;*/
1285 GETTER_AND_SETTER(int,maxsubsteps) /*;*/
1286 #undef GETTER_AND_SETTER
1287
1288 long integrator_getnsamples(IntegratorSystem *sys){
1289 assert(sys!=NULL);
1290 assert(sys->samples!=NULL);
1291 return samplelist_length(sys->samples);
1292 }
1293
1294 double integrator_getsample(IntegratorSystem *sys, long i){
1295 assert(sys!=NULL);
1296 assert(sys->samples!=NULL);
1297 return samplelist_get(sys->samples,i);
1298 }
1299
1300 void integrator_setsample(IntegratorSystem *sys, long i,double xi){
1301 assert(sys!=NULL);
1302 assert(sys->samples!=NULL);
1303 samplelist_set(sys->samples,i,xi);
1304 }
1305
1306 const dim_type *integrator_getsampledim(IntegratorSystem *sys){
1307 assert(sys!=NULL);
1308 assert(sys->samples!=NULL);
1309 return samplelist_dim(sys->samples);
1310 }
1311
1312 ASC_DLLSPEC(long) integrator_getcurrentstep(IntegratorSystem *sys){
1313 return sys->currentstep;
1314 }
1315
1316 /*------------------------------------------------------------------------------
1317 GET/SET VALUE OF THE INDEP VARIABLE
1318 */
1319
1320 /**
1321 Retrieve the value of the independent variable (time) from ASCEND
1322 and return it as a double.
1323 */
1324 double integrator_get_t(IntegratorSystem *sys){
1325 assert(sys->x!=NULL);
1326 return var_value(sys->x);
1327 }
1328
1329 /**
1330 Set the value of the independent variable (time) in ASCEND.
1331 */
1332 void integrator_set_t(IntegratorSystem *sys, double value){
1333 var_set_value(sys->x, value);
1334 /* CONSOLE_DEBUG("set_t = %g", value); */
1335 }
1336
1337 /*------------------------------------------------------------------------------
1338 PASSING DIFFERENTIAL VARIABLES AND THEIR DERIVATIVES TO/FROM THE SOLVER
1339 */
1340 /**
1341 Retrieve the current values of the derivatives of the y-variables
1342 and stick them in the/an array that the integrator will use.
1343
1344 If the pointer 'y' is NULL, the necessary space is allocated (and
1345 must be freed somewhere else).
1346 */
1347 double *integrator_get_y(IntegratorSystem *sys, double *y) {
1348 long i;
1349
1350 if (y==NULL) {
1351 y = ASC_NEW_ARRAY_CLEAR(double, sys->n_y+1);
1352 /* C y[0] <==> ascend d.y[1] <==> f77 y(1) */
1353 }
1354
1355 for (i=0; i< sys->n_y; i++) {
1356 assert(sys->y[i]!=NULL);
1357 y[i] = var_value(sys->y[i]);
1358 /* CONSOLE_DEBUG("ASCEND --> y[%ld] = %g", i+1, y[i]); */
1359 }
1360 return y;
1361 }
1362
1363 /**
1364 Take the values of the derivatives from the array that the integrator
1365 uses, and use them to update the values of the corresponding variables
1366 in ASCEND.
1367 */
1368 void integrator_set_y(IntegratorSystem *sys, double *y) {
1369 long i;
1370 #ifdef SOLVE_DEBUG
1371 char *varname;
1372 #endif
1373
1374 for (i=0; i < sys->n_y; i++) {
1375 assert(sys->y[i]!=NULL);
1376 var_set_value(sys->y[i],y[i]);
1377 #ifdef SOLVE_DEBUG
1378 varname = var_make_name(sys->system, sys->y[i]);
1379 CONSOLE_DEBUG("y[%ld] = %g --> '%s'", i+1, y[i], varname);
1380 ASC_FREE(varname);
1381 #endif
1382 }
1383 }
1384
1385 /**
1386 Send the values of the derivatives of the 'y' variables to the solver.
1387 Allocate space for an array if necessary.
1388
1389 Any element in sys->ydot that is NULL will be passed over (the value
1390 won't be modified in dydx).
1391 */
1392 double *integrator_get_ydot(IntegratorSystem *sys, double *dydx) {
1393 long i;
1394
1395 if (dydx==NULL) {
1396 dydx = ASC_NEW_ARRAY_CLEAR(double, sys->n_y+1);
1397 /* C dydx[0] <==> ascend d.dydx[1] <==> f77 ydot(1) */
1398 }
1399
1400 for (i=0; i < sys->n_y; i++) {
1401 if(sys->ydot[i]!=NULL){
1402 dydx[i] = var_value(sys->ydot[i]);
1403 }
1404 /* CONSOLE_DEBUG("ASCEND --> ydot[%ld] = %g", i+1, dydx[i]); */
1405 }
1406 return dydx;
1407 }
1408
1409 void integrator_set_ydot(IntegratorSystem *sys, double *dydx) {
1410 long i;
1411 #ifdef SOLVE_DEBUG
1412 char *varname;
1413 #endif
1414 for (i=0; i < sys->n_y; i++) {
1415 if(sys->ydot[i]!=NULL){
1416 var_set_value(sys->ydot[i],dydx[i]);
1417 #ifdef SOLVE_DEBUG
1418 varname = var_make_name(sys->system, sys->ydot[i]);
1419 CONSOLE_DEBUG("ydot[%ld] = \"%s\" = %g --> ASCEND", i+1, varname, dydx[i]);
1420 ASC_FREE(varname);
1421 }else{
1422 CONSOLE_DEBUG("ydot[%ld] = %g (internal)", i+1, dydx[i]);
1423 #endif
1424 }
1425 }
1426 }
1427
1428 /**
1429 Retrieve the values of 'ode_atol' properties of each of y-variables,
1430 for use in setting absolute error tolerances for the Integrator.
1431
1432 If the pointer 'atol' is NULL, the necessary space is allocated (and
1433 must be freed somewhere else).
1434 */
1435 double *integrator_get_atol(IntegratorSystem *sys, double *atol){
1436 long i;
1437 char *varname;
1438
1439 if (atol==NULL) {
1440 atol = ASC_NEW_ARRAY_CLEAR(double, sys->n_y);
1441 }
1442
1443 for (i=0; i< sys->n_y; i++) {
1444 assert(sys->y[i]!=NULL);
1445 atol[i] = var_odeatol(sys->y[i]);
1446 assert(atol[i]!=-1);
1447 varname = var_make_name(sys->system,sys->y[i]);
1448 CONSOLE_DEBUG("%s.ode_atol = %8.2e",varname,atol[i]);
1449 ASC_FREE(varname);
1450 }
1451 return atol;
1452 }
1453
1454 /*-------------------------------------------------------------
1455 RETRIEVING OBSERVATION DATA
1456 */
1457
1458 /**
1459 This function takes the inst in the solver and returns the vector of
1460 observation variables that are located in the submodel d.obs array.
1461 */
1462 double *integrator_get_observations(IntegratorSystem *sys, double *obsi) {
1463 long i;
1464
1465 if (obsi==NULL) {
1466 obsi = ASC_NEW_ARRAY_CLEAR(double, sys->n_obs+1);
1467 }
1468
1469 /* C obsi[0] <==> ascend d.obs[1] */
1470
1471 for (i=0; i < sys->n_obs; i++) {
1472 obsi[i] = var_value(sys->obs[i]);
1473 /* CONSOLE_DEBUG("*get_d_obs[%ld] = %g\n", i+1, obsi[i]); */
1474 }
1475 return obsi;
1476 }
1477
1478 struct var_variable *integrator_get_observed_var(IntegratorSystem *sys, const long i){
1479 assert(i>=0);
1480 assert(i<sys->n_obs);
1481 return sys->obs[i];
1482 }
1483
1484 /**
1485 @NOTE Although this shouldn't be required for implementation of solver
1486 engines, this is useful for GUI reporting of integration results.
1487 */
1488 struct var_variable *integrator_get_independent_var(IntegratorSystem *sys){
1489 return sys->x;
1490 }
1491
1492
1493 /*----------------------------------------------------
1494 Build an analytic jacobian for solving the state system
1495
1496 This necessarily ugly piece of code attempts to create a unique
1497 list of relations that explicitly contain the variables in the
1498 given input list. The utility of this information is that we know
1499 exactly which relations must be differentiated, to fill in the
1500 df/dy matrix. If the problem has very few derivative terms, this will
1501 be of great savings. If the problem arose from the discretization of
1502 a pde, then this will be not so useful. The decision wether to use
1503 this function or to simply differentiate the entire relations list
1504 must be done before calling this function.
1505
1506 Final Note: the callee owns the array, but not the array elements.
1507 */
1508 #define AVG_NUM_INCIDENT 4
1509
1510
1511 /**
1512 This function helps to arrange the observation variables in a sensible order.
1513 The 'obs_id' child instance of v, if present, is assigned the value of the
1514 given parameter 'index'.
1515 */
1516 static void Integ_SetObsId(struct var_variable *v, long index){
1517 struct Instance *c, *i;
1518 i = var_instance(v);
1519 assert(i!=NULL);
1520 c = ChildByChar(i,OBSINDEX);
1521 if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {
1522 return;
1523 }
1524 SetIntegerAtomValue(c,index,0);
1525 }
1526
1527 /**
1528 Compares observation structs. NULLs should end up at far end.
1529 */
1530 static int Integ_CmpObs(struct Integ_var_t *v1, struct Integ_var_t *v2){
1531 if(v1 == NULL)return 1;
1532 if(v2 == NULL)return -1;
1533 if(v1->index > v2->index)return 1;
1534 if(v1->index == v2->index)return 0;
1535 return -1;
1536 }
1537
1538 /**
1539 Compares dynamic vars structs. NULLs should end up at far end.
1540 List should be sorted primarily by index and then by type, in order
1541 of increasing value of both.
1542 */
1543 static int Integ_CmpDynVars(struct Integ_var_t *v1, struct Integ_var_t *v2){
1544 if(v1 == NULL)return 1;
1545 if(v2 == NULL)return -1;
1546 if(v1->index > v2->index)return 1;
1547 if(v1->index != v2->index)return -1;
1548 if(v1->type > v2->type)return 1;
1549 return -1;
1550 }
1551 /*----------------------------
1552 Output handling to the GUI/interface.
1553 */
1554
1555 int integrator_set_reporter(IntegratorSystem *sys
1556 , IntegratorReporter *reporter
1557 ){
1558 assert(sys!=NULL);
1559 sys->reporter = reporter;
1560 /* ERROR_REPORTER_HERE(ASC_PROG_NOTE,"INTEGRATOR REPORTER HOOKS HAVE BEEN SET\n"); */
1561 return 1;
1562 }
1563
1564 int integrator_output_init(IntegratorSystem *sys){
1565 assert(sys!=NULL);
1566 assert(sys->reporter!=NULL);
1567 if(sys->reporter->init!=NULL){
1568 /* call the specified output function */
1569 return (*(sys->reporter->init))(sys);
1570 }
1571 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter init method");
1572 return 1;
1573 }
1574
1575 int integrator_output_write(IntegratorSystem *sys){
1576 static int reported_already=0;
1577 assert(sys!=NULL);
1578 if(sys->reporter->write!=NULL){
1579 return (*(sys->reporter->write))(sys);
1580 }
1581 if(!reported_already){
1582 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter write method (this message only shown once)");
1583 reported_already=1;
1584 }
1585 return 1;
1586 }
1587
1588 int integrator_output_write_obs(IntegratorSystem *sys){
1589 static int reported_already=0;
1590 assert(sys!=NULL);
1591 if(sys->reporter->write_obs!=NULL){
1592 return (*(sys->reporter->write_obs))(sys);
1593 }
1594 if(!reported_already){
1595 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter write_obs method (this message only shown once)");
1596 reported_already=1;
1597 }
1598 return 1;
1599 }
1600
1601 int integrator_output_close(IntegratorSystem *sys){
1602 assert(sys!=NULL);
1603 if(sys->reporter->close!=NULL){
1604 return (*(sys->reporter->close))(sys);
1605 }
1606 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No integrator reporter close method");
1607 return 1;
1608 }
1609
1610 /**
1611 Decode status codes from the integrator, and output them via FPRINTF.
1612 */
1613 int integrator_checkstatus(slv_status_t status) {
1614 if (status.converged) {
1615 return 1;
1616 }
1617 if (status.diverged) {
1618 FPRINTF(stderr, "The derivative system did not converge. Integration will terminate.");
1619 return 0;
1620 }
1621 if (status.inconsistent) {
1622 FPRINTF(stderr, "A numerically inconsistent state was discovered while "
1623 "calculating derivatives. Integration will terminate.");
1624 return 0;
1625 }
1626 if (status.time_limit_exceeded) {
1627 FPRINTF(stderr, "The time limit was exceeded while calculating "
1628 "derivatives. Integration will terminate.");
1629 return 0;
1630 }
1631 if (status.iteration_limit_exceeded) {
1632 FPRINTF(stderr, "The iteration limit was exceeded while calculating "
1633 "derivatives. Integration will terminate.");
1634 return 0;
1635 }
1636 if (status.panic) {
1637 FPRINTF(stderr, "The user patience limit was exceeded while "
1638 "calculating derivatives. Integration will terminate.");
1639 return 0;
1640 }
1641 return 0;
1642 }

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