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

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