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

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