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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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