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 |
*//** @file |
19 |
Implementation of non-GUI-specific parts of the Integration Interface. |
20 |
|
21 |
Here we are integrating a model with independent variable t for a specified |
22 |
set of timesteps, and and observing the values of certain of the |
23 |
independent variables. |
24 |
|
25 |
You can specify the values of the time steps that you want via the |
26 |
'samplelist' interface. |
27 |
|
28 |
You can have your results reported however you like using the |
29 |
'IntegratorReporter' interface. |
30 |
|
31 |
There's nothing here yet to explicitly support DAEs -- that's the next task. |
32 |
|
33 |
(old annotation:) |
34 |
The following functions fetch and set parts with names specific to |
35 |
the type definitions in ivp.lib, the ASCEND initial value problem |
36 |
model file. They are for use by any ivp solver interface. |
37 |
All names are given at the scope of ivp. |
38 |
|
39 |
*//* |
40 |
by John Pye, May 2006. |
41 |
|
42 |
Derived from tcltk/.../Integrators.h (with heavy editing). |
43 |
Removed all global variables and created a 'IntegratorReporter' interface |
44 |
for reporting results of integration to whatever interface you're using. |
45 |
*/ |
46 |
|
47 |
#ifndef ASC_INTEGRATOR_H |
48 |
#define ASC_INTEGRATOR_H |
49 |
|
50 |
#include <utilities/ascConfig.h> |
51 |
#include <utilities/ascMalloc.h> |
52 |
#include <general/list.h> |
53 |
#include <general/dstring.h> |
54 |
#include <compiler/compiler.h> |
55 |
#include <compiler/instance_enum.h> |
56 |
|
57 |
#include <compiler/fractions.h> |
58 |
#include <compiler/dimen.h> |
59 |
#include <compiler/units.h> |
60 |
#include <compiler/module.h> |
61 |
#include <compiler/library.h> |
62 |
#include <compiler/expr_types.h> |
63 |
#include <compiler/child.h> |
64 |
#include <compiler/type_desc.h> |
65 |
#include <compiler/atomvalue.h> |
66 |
#include <compiler/instance_name.h> |
67 |
#include <compiler/instquery.h> |
68 |
#include <compiler/parentchild.h> |
69 |
#include <compiler/symtab.h> |
70 |
#include <compiler/instance_io.h> |
71 |
|
72 |
#include "slv_types.h" |
73 |
#include "mtx.h" |
74 |
#include "var.h" |
75 |
#include "rel.h" |
76 |
#include "discrete.h" |
77 |
#include "conditional.h" |
78 |
#include "logrel.h" |
79 |
#include "bnd.h" |
80 |
#include "slv_common.h" |
81 |
#include "samplelist.h" |
82 |
|
83 |
/*---------------------------*/ |
84 |
|
85 |
#ifdef ASC_WITH_IDA |
86 |
# define IDA_OPTIONAL S I(IDA,integrator_ida_internals) |
87 |
#else |
88 |
# define IDA_OPTIONAL |
89 |
#endif |
90 |
|
91 |
/* we add IDA to the list of integrators at build time, if it is selected */ |
92 |
#define INTEG_LIST \ |
93 |
I(LSODE ,integrator_lsode_internals) \ |
94 |
IDA_OPTIONAL \ |
95 |
S I(AWW ,integrator_aww_internals) |
96 |
|
97 |
/** |
98 |
Struct containin the list of supported integrators |
99 |
*/ |
100 |
typedef enum{ |
101 |
#define I(N,P) INTEG_##N |
102 |
#define S , |
103 |
INTEG_LIST |
104 |
#undef I |
105 |
#undef S |
106 |
,INTEG_UNKNOWN /* must be last in the list*/ |
107 |
} IntegratorEngine; |
108 |
|
109 |
|
110 |
typedef struct{ |
111 |
IntegratorEngine id; |
112 |
const char *name; |
113 |
} IntegratorLookup; |
114 |
|
115 |
/*------------------------ |
116 |
abstraction of the integrator output interface |
117 |
*/ |
118 |
struct IntegratorSystemStruct; |
119 |
|
120 |
/** |
121 |
Initialisation. This hook allows initialisation of the GUI or reporting |
122 |
mechanism to be performed when integration begins |
123 |
*/ |
124 |
typedef int IntegratorOutputInitFn(struct IntegratorSystemStruct *); |
125 |
|
126 |
/** Status report. This hook allows raw data to be output as integration |
127 |
proceeds, or for a GUI to perform status updates. An integrator should check |
128 |
the return status on this one, as this is the suggested way to perform |
129 |
GUI interruption of the integrator. |
130 |
|
131 |
@return 1 on success, 0 on user interrupt |
132 |
*/ |
133 |
typedef int IntegratorOutputWriteFn(struct IntegratorSystemStruct *); |
134 |
|
135 |
/** |
136 |
Observation reporting. This hook should be implemented to record |
137 |
observations in a way that can be presented to the use, recorded in a |
138 |
file, etc. |
139 |
*/ |
140 |
typedef int IntegratorOutputWriteObsFn(struct IntegratorSystemStruct *); |
141 |
|
142 |
/** |
143 |
Finalisation. This hook can be used to terminate recording of observations, |
144 |
close files, terminate GUI status reporting, etc. |
145 |
*/ |
146 |
typedef int IntegratorOutputCloseFn(struct IntegratorSystemStruct *); |
147 |
|
148 |
/** |
149 |
This struct allows arbitrary functions to be used for the reporting |
150 |
of integrator progress. |
151 |
*/ |
152 |
typedef struct{ |
153 |
IntegratorOutputInitFn *init; |
154 |
IntegratorOutputWriteFn *write; |
155 |
IntegratorOutputWriteObsFn *write_obs; |
156 |
IntegratorOutputCloseFn *close; |
157 |
} IntegratorReporter; |
158 |
|
159 |
/*------------------------------------*/ |
160 |
/* INTEGRATOR INTERNALS - define the routines that consitute a specific integrator */ |
161 |
|
162 |
/* forward decl */ |
163 |
struct IntegratorSystemStruct; |
164 |
|
165 |
typedef void IntegratorCreateFn(struct IntegratorSystemStruct *blsys); |
166 |
/**< Integrators must provide a routine here that allocates the necessary |
167 |
data structures. |
168 |
*/ |
169 |
|
170 |
typedef int IntegratorParamsDefaultFn(struct IntegratorSystemStruct *blsys); |
171 |
/**< |
172 |
Integrators must provide a function like this that can be used to retrieve |
173 |
the default set of parameters. |
174 |
*/ |
175 |
|
176 |
typedef int IntegratorAnalyseFn(struct IntegratorSystemStruct *blsys); |
177 |
|
178 |
|
179 |
typedef int IntegratorSolveFn(struct IntegratorSystemStruct *blsys |
180 |
, unsigned long start_index, unsigned long finish_index); |
181 |
/**< |
182 |
Integrators must provide a function like this that actually runs the |
183 |
integration. |
184 |
*/ |
185 |
|
186 |
typedef void IntegratorFreeFn(void *enginedata); |
187 |
/**< |
188 |
Integrators must provide a function like this that frees internal |
189 |
data that they have allocated in their 'enginedata' structure. |
190 |
*/ |
191 |
|
192 |
typedef struct{ |
193 |
IntegratorCreateFn *createfn; |
194 |
IntegratorParamsDefaultFn *paramsdefaultfn; |
195 |
IntegratorAnalyseFn *analysefn; |
196 |
IntegratorSolveFn *solvefn; |
197 |
IntegratorFreeFn *freefn; |
198 |
IntegratorEngine engine; |
199 |
const char name[]; |
200 |
} IntegratorInternals; |
201 |
|
202 |
/*------------------------------------*/ |
203 |
/** |
204 |
Initial Value Problem description struct. Anyone making a copy of |
205 |
the y, ydot, or obs pointers who plans to free that pointer later |
206 |
should increment the reference count for that pointer so if blsys |
207 |
is destroyed later we don't end up double freeing it. Note this |
208 |
scheme will only work for the first copy and could be eliminated |
209 |
if we decoupled blsode from lsode as we will once things restabilize. |
210 |
|
211 |
JP: adding to this structure the instance being solved, the engine |
212 |
that we're solving it with, and a struct containing the output function ptrs |
213 |
*/ |
214 |
struct IntegratorSystemStruct{ |
215 |
struct Instance *instance; /**< not sure if this one is really necessary... -- JP */ |
216 |
slv_system_t system; /**< the system that we're integrating in ASCEND */ |
217 |
IntegratorEngine engine; /**< enum containing the ID of the integrator engine we're using */ |
218 |
const IntegratorInternals *internals;/**< pointers to the various functions belonging to this integrator */ |
219 |
IntegratorReporter *reporter;/**< functions for reporting integration results */ |
220 |
SampleList *samples; /**< pointer to the list of samples. we *don't own* this **/ |
221 |
void *enginedata; /**< space where the integrator engine can store stuff */ |
222 |
void *clientdata; /**< any stuff that the GUI/CLI needs to associate with this */ |
223 |
|
224 |
slv_parameters_t params; /**< structure containing parameters applicable to this Integrator */ |
225 |
|
226 |
int nstates; /**< was a local global in integrator.c, moved it here. */ |
227 |
int nderivs; /**< ditto, as for nstates */ |
228 |
|
229 |
/* temporaries for build. these elements will be NULL to clients */ |
230 |
struct gl_list_t *indepvars; /**< all apparent independent vars */ |
231 |
struct gl_list_t *dynvars; /**< all state and deriv instances plus indices */ |
232 |
struct gl_list_t *obslist; /**< observation instance plus indices */ |
233 |
struct gl_list_t *states; /**< ordered list of state variables and indices*/ |
234 |
struct gl_list_t *derivs; /**< ordered list of derivative (ydot) insts */ |
235 |
|
236 |
/* persistent, these elements are valid to C clients. */ |
237 |
struct var_variable *x; /**< independent variable */ |
238 |
struct var_variable **y; /**< array form of states */ |
239 |
struct var_variable **ydot; /**< array form of derivatives */ |
240 |
struct var_variable **obs; /**< array form of observed variables */ |
241 |
long *y_id; /**< array form of y/ydot user indices, for DAEs we use negatives here for derivative vars */ |
242 |
long *obs_id; /**< array form of obs user indices */ |
243 |
long n_y; |
244 |
long n_obs; |
245 |
long currentstep; /**< current step number (also @see integrator_getnsamples) */ |
246 |
|
247 |
int ycount; /**< number of external references to y */ |
248 |
int ydotcount; /**< number of external references to ydot */ |
249 |
int obscount; /**< number of external references to obs */ |
250 |
|
251 |
/** @TODO move the following to the 'params' structure? Or maybe better not to? */ |
252 |
int maxsubsteps; /**< most steps between mesh poins */ |
253 |
double stepzero; /**< initial step length, SI units. */ |
254 |
double minstep; /**< shortest step length, SI units. */ |
255 |
double maxstep; /**< longest step length, SI units. */ |
256 |
}; |
257 |
|
258 |
typedef struct IntegratorSystemStruct IntegratorSystem; |
259 |
|
260 |
/*------------------------------------------------------------------------------ |
261 |
PUBLIC INTERFACE (for use by the GUI/CLI)] |
262 |
*/ |
263 |
|
264 |
ASC_DLLSPEC(IntegratorSystem *) integrator_new(slv_system_t sys, struct Instance *inst); |
265 |
|
266 |
ASC_DLLSPEC(int) integrator_analyse(IntegratorSystem *blsys); |
267 |
|
268 |
/** |
269 |
These routines will hopefully be sharable by different Integrator engines. |
270 |
The idea is that if we change from the 'ode_id' and 'ode_type' syntax, |
271 |
the only place in the Integrator code that would be affected is here. |
272 |
However for the moment, we allow Integrators to specify how they would |
273 |
like to analyse the system, and for that, we export these routines so that |
274 |
they can be referred to in lsode.h, ida.h, etc. |
275 |
*/ |
276 |
ASC_DLLSPEC(int) integrator_analyse_ode(IntegratorSystem *blsys); |
277 |
ASC_DLLSPEC(int) integrator_analyse_dae(IntegratorSystem *blsys); |
278 |
/**< see integrator_analyse_ode */ |
279 |
|
280 |
ASC_DLLSPEC(int) integrator_solve(IntegratorSystem *blsys, long i0, long i1); |
281 |
/**< |
282 |
Takes the type of integrator and sets up the global variables into the |
283 |
current integration instance. |
284 |
|
285 |
@return 1 on success, 0 on failure. |
286 |
*/ |
287 |
|
288 |
ASC_DLLSPEC(void) integrator_free(IntegratorSystem *blsys); |
289 |
/**< |
290 |
Deallocates any memory used and sets all integration global points to NULL. |
291 |
*/ |
292 |
|
293 |
ASC_DLLSPEC(const IntegratorLookup *) integrator_get_engines(); |
294 |
/**< |
295 |
Return a {INTEG_UNKNOWN,NULL} terminated list of integrator currently |
296 |
available. At present this is determined at compile time but will be |
297 |
changed to being determined at load time if necessary. |
298 |
*/ |
299 |
|
300 |
ASC_DLLSPEC(int) integrator_set_engine(IntegratorSystem *blsys, IntegratorEngine engine); |
301 |
/**< |
302 |
Sets the engine for this integrator. Checks that the integrator can be used |
303 |
on the given system. |
304 |
|
305 |
@return 0 on success, |
306 |
1 if unknown engine, |
307 |
2 if invalid engine (not suitable for this present problem) |
308 |
*/ |
309 |
|
310 |
ASC_DLLSPEC(IntegratorEngine) integrator_get_engine(const IntegratorSystem *blsys); |
311 |
/**< |
312 |
Returns the engine (ID) selected for use in this IntegratorSystem (may return |
313 |
INTEG_UNKNOWN if none or invalid setting). |
314 |
*/ |
315 |
|
316 |
ASC_DLLSPEC(int) integrator_params_get(const IntegratorSystem *blsys, slv_parameters_t *parameters); |
317 |
/**< |
318 |
Copies the current set of Integrator parameters into the indicated parameter set, 'parameters'. |
319 |
|
320 |
@TODO test these, not sure if the memory copy stuff is right or not. |
321 |
|
322 |
@return 0 on success |
323 |
*/ |
324 |
|
325 |
ASC_DLLSPEC(int) integrator_params_set(IntegratorSystem *blsys, const slv_parameters_t *parameters); |
326 |
/**< |
327 |
Copies the the given parameter set into the Integrator's data structure, |
328 |
which immediately makes the new values available to the integrator engine. |
329 |
|
330 |
@TODO test these, not sure if the memory copy stuff is right or not. |
331 |
|
332 |
@return 0 on success |
333 |
*/ |
334 |
|
335 |
ASC_DLLSPEC(int) integrator_set_reporter(IntegratorSystem *blsys |
336 |
, IntegratorReporter *r); |
337 |
/**< |
338 |
Use this function from your interface to tell the integrator how it needs |
339 |
to make its output. You need to pass in a struct containing the required |
340 |
function pointers. (We will wrap IntegratorReporter and access it from |
341 |
Python, eventually) |
342 |
*/ |
343 |
|
344 |
ASC_DLLSPEC(int) integrator_find_indep_var(IntegratorSystem *blsys); |
345 |
/**< |
346 |
Attempt to locate the independent variable. It will be stored in the blsys |
347 |
structure if found. |
348 |
@return 1 if found, 0 if not found. |
349 |
*/ |
350 |
|
351 |
extern int integrator_checkstatus(slv_status_t status); |
352 |
/**< |
353 |
* Takes a status and checks for convergence. If converged then return 1 |
354 |
* else return 0 and print error message if appropriate. |
355 |
*/ |
356 |
|
357 |
ASC_DLLSPEC(void) integrator_set_samples(IntegratorSystem *blsys, SampleList *samples); |
358 |
/**< |
359 |
Sets values of time samples to the values given (ns of them) and |
360 |
keeps both the dim pointer and vector given. The vector and dimp |
361 |
given may be freed by calling this again, but xsamples owns |
362 |
them until then. If called with ns<1 or values==NULL, will free any |
363 |
previously captured values/dimp. If called with dimp==NULL, will |
364 |
assume WildDimension. Don't call this with a dimp we can't free later. |
365 |
Return is 1 if for some reason we can't set as expected, 0 otherwise. |
366 |
*/ |
367 |
|
368 |
ASC_DLLSPEC(void) integrator_set_stepzero(IntegratorSystem *blsys, double); |
369 |
ASC_DLLSPEC(double) integrator_get_stepzero(IntegratorSystem *blsys); |
370 |
/**< |
371 |
Returns the length of the initial step user specified, |
372 |
or 0.0 if none was set. |
373 |
*/ |
374 |
|
375 |
ASC_DLLSPEC(void) integrator_set_minstep(IntegratorSystem *blsys, double); |
376 |
ASC_DLLSPEC(double) integrator_get_minstep(IntegratorSystem *blsys); |
377 |
/**< |
378 |
Returns the length of the longest allowable step. |
379 |
or 0.0 if none was set by user. |
380 |
*/ |
381 |
|
382 |
ASC_DLLSPEC(void) integrator_set_maxstep(IntegratorSystem *blsys, double); |
383 |
ASC_DLLSPEC(double) integrator_get_maxstep(IntegratorSystem *blsys); |
384 |
/**< |
385 |
Returns the length of the shortest allowable step. |
386 |
or 0.0 if none was set by user. |
387 |
*/ |
388 |
|
389 |
ASC_DLLSPEC(void) integrator_set_maxsubsteps(IntegratorSystem *blsys, int); |
390 |
ASC_DLLSPEC(int) integrator_get_maxsubsteps(IntegratorSystem *blsys); |
391 |
/**< |
392 |
Returns the most internal steps allowed between |
393 |
two time samples, or 0 if none was set by user. |
394 |
*/ |
395 |
|
396 |
ASC_DLLSPEC(long) integrator_getnsamples(IntegratorSystem *blsys); |
397 |
/**< |
398 |
Returns the number of values currently stored in xsamples. |
399 |
*/ |
400 |
|
401 |
ASC_DLLSPEC(long) integrator_getcurrentstep(IntegratorSystem *blsys); |
402 |
/**< |
403 |
Returns the current step number (blsys->currentstep). Should be reset to |
404 |
zero inside your solver's integrator_*_solve method. |
405 |
*/ |
406 |
|
407 |
extern double integrator_getsample(IntegratorSystem *blsys, long i); |
408 |
/**< The following functions fetch and set parts with names specific to |
409 |
the type definitions in ivp.lib, the ASCEND initial value problem |
410 |
model file. They are for use by any ivp solver interface. |
411 |
All names are given at the scope of ivp. |
412 |
|
413 |
Returns the value stored in xsamples[i]. Will whine if |
414 |
if xsample[i] does not exist. No, there is no way to get |
415 |
back the pointer to the xsamples vector. |
416 |
*/ |
417 |
|
418 |
extern void integrator_setsample(IntegratorSystem *blsys, long i, double val); |
419 |
/**< |
420 |
Sets the value stored in xsamples[i] to val. Will whine if |
421 |
if xsample[i] does not exist. No, there is no way to get |
422 |
back the pointer to the xsamples vector. |
423 |
*/ |
424 |
|
425 |
extern const dim_type *integrator_getsampledim(IntegratorSystem *blsys); |
426 |
/**< |
427 |
Returns a pointer to the dimensionality of the samples currently |
428 |
stored, or NULL if none are stored. Do not free this pointer: we |
429 |
own it. |
430 |
*/ |
431 |
|
432 |
/*------------------------------------------------------------------------------ |
433 |
BACKEND INTERFACE (for use by the integrator engine) |
434 |
*/ |
435 |
|
436 |
/* Problem size parameters. */ |
437 |
|
438 |
/* Parts of type definition derivatives refinements. */ |
439 |
|
440 |
ASC_DLLSPEC(double) integrator_get_t(IntegratorSystem *blsys); |
441 |
/**< |
442 |
Gets value of the independent variable value from the ASCEND model |
443 |
(retrieves the value from the ASCEND compiler's instance hierarchy) |
444 |
*/ |
445 |
|
446 |
extern void integrator_set_t(IntegratorSystem *blsys, double value); |
447 |
/**< |
448 |
Sets value of the independent variable in the model (inserts the value in |
449 |
the correct place in the compiler's instance hierarchy via the var_variable |
450 |
API). |
451 |
*/ |
452 |
|
453 |
extern double *integrator_get_y(IntegratorSystem *blsys, double *vector); |
454 |
/**< |
455 |
Gets the value of vector 'y' from the model (what 'y' is depends on your |
456 |
particular integrator). |
457 |
|
458 |
If vector given is NULL, allocates vector, which the caller then owns. |
459 |
Vector, if given, should be IntegGet_d_neq()+1 long. |
460 |
*/ |
461 |
|
462 |
extern void integrator_set_y(IntegratorSystem *blsys, double *vector); |
463 |
/**< |
464 |
Sets d.y[] to values in vector. |
465 |
*/ |
466 |
|
467 |
extern double *integrator_get_ydot(IntegratorSystem *blsys, double *vector); |
468 |
/**< |
469 |
Returns the vector ydot (derivatives of the 'states') |
470 |
|
471 |
If vector given is NULL, allocates vector, which the caller then owns. |
472 |
Vector, if given, should be IntegGet_d_neq()+1 long. |
473 |
*/ |
474 |
|
475 |
extern void integrator_set_ydot(IntegratorSystem *blsys, double *vector); |
476 |
/**< |
477 |
Sets d.dydx[] to values in vector. |
478 |
*/ |
479 |
|
480 |
ASC_DLLSPEC(double *) integrator_get_observations(IntegratorSystem *blsys, double *vector); |
481 |
/**< |
482 |
Returns the vector d.obs. |
483 |
Vector should be of sufficient length (g_intginst_d_n_obs+1). |
484 |
If NULL is given a vector is allocated and is the caller's |
485 |
responsibility to deallocate. |
486 |
*/ |
487 |
|
488 |
ASC_DLLSPEC(struct var_variable *) integrator_get_observed_var(IntegratorSystem *blsys, const long i); |
489 |
/**< |
490 |
Returns the var_variable contained in the ith position in the observed variable list. |
491 |
*/ |
492 |
|
493 |
ASC_DLLSPEC(struct var_variable *) integrator_get_independent_var(IntegratorSystem *blsys); |
494 |
/**< |
495 |
Return a pointer to the variable identified as the independent variable. |
496 |
*/ |
497 |
|
498 |
extern double *integrator_get_atol(IntegratorSystem *blsys, double *vector); |
499 |
/**< |
500 |
Return abolute tolerance values for the 'ode_atol' values in the MODEL. |
501 |
*/ |
502 |
|
503 |
/*------------------------------- |
504 |
Stuff to facilitate output to the interface |
505 |
*/ |
506 |
|
507 |
/** |
508 |
This call should be used to get the file streams ready, output column |
509 |
headings etc. |
510 |
*/ |
511 |
extern int integrator_output_init(IntegratorSystem *blsys); |
512 |
|
513 |
/** |
514 |
This call should be used to output the immediate integration results from |
515 |
a single timestep. For an ODE solver, this means just the 'y' values but |
516 |
perhaps not the values of the algebraic varaibles, which must be calculated |
517 |
separately. They'll be stored in the Integ_system_t somehow (not yet |
518 |
known how) |
519 |
|
520 |
@return 1 on success, return 0 on failure (integration will be cancelled) |
521 |
*/ |
522 |
extern int integrator_output_write(IntegratorSystem *blsys); |
523 |
|
524 |
/** |
525 |
Write out the 'observed values' for the integration. In the case of ODE |
526 |
integration, we assume that the values of all the algabraic variables are |
527 |
also now calculated. |
528 |
*/ |
529 |
extern int integrator_output_write_obs(IntegratorSystem *blsys); |
530 |
|
531 |
/** |
532 |
This call will close file stream and perhaps perform some kind of |
533 |
user notification or screen update, etc. |
534 |
*/ |
535 |
extern int integrator_output_close(IntegratorSystem *blsys); |
536 |
|
537 |
#endif |