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 |
Struct containin the list of supported integrators |
86 |
*/ |
87 |
typedef enum{ |
88 |
INTEG_UNKNOWN, |
89 |
INTEG_LSODE, |
90 |
INTEG_IDA |
91 |
} IntegratorEngine; |
92 |
|
93 |
/*------------------------ |
94 |
abstraction of the integrator output interface |
95 |
*/ |
96 |
struct IntegratorSystemStruct; |
97 |
|
98 |
/** |
99 |
Initialisation. This hook allows initialisation of the GUI or reporting |
100 |
mechanism to be performed when integration begins |
101 |
*/ |
102 |
typedef int IntegratorOutputInitFn(struct IntegratorSystemStruct *); |
103 |
|
104 |
/** Status report. This hook allows raw data to be output as integration |
105 |
proceeds, or for a GUI to perform status updates. An integrator should check |
106 |
the return status on this one, as this is the suggested way to perform |
107 |
GUI interruption of the integrator. |
108 |
|
109 |
@return 1 on success, 0 on user interrupt |
110 |
*/ |
111 |
typedef int IntegratorOutputWriteFn(struct IntegratorSystemStruct *); |
112 |
|
113 |
/** |
114 |
Observation reporting. This hook should be implemented to record |
115 |
observations in a way that can be presented to the use, recorded in a |
116 |
file, etc. |
117 |
*/ |
118 |
typedef int IntegratorOutputWriteObsFn(struct IntegratorSystemStruct *); |
119 |
|
120 |
/** |
121 |
Finalisation. This hook can be used to terminate recording of observations, |
122 |
close files, terminate GUI status reporting, etc. |
123 |
*/ |
124 |
typedef int IntegratorOutputCloseFn(struct IntegratorSystemStruct *); |
125 |
|
126 |
/** |
127 |
This struct allows arbitrary functions to be used for the reporting |
128 |
of integrator progress. |
129 |
*/ |
130 |
typedef struct{ |
131 |
IntegratorOutputInitFn *init; |
132 |
IntegratorOutputWriteFn *write; |
133 |
IntegratorOutputWriteObsFn *write_obs; |
134 |
IntegratorOutputCloseFn *close; |
135 |
} IntegratorReporter; |
136 |
|
137 |
/*------------------------------------*/ |
138 |
/** |
139 |
Initial value problem description struct. Anyone making a copy of |
140 |
the y, ydot, or obs pointers who plans to free that pointer later |
141 |
should increment the reference count for that pointer so if blsys |
142 |
is destroyed later we don't end up double freeing it. Note this |
143 |
scheme will only work for the first copy and could be eliminated |
144 |
if we decoupled blsode from lsode as we will once things restabilize. |
145 |
|
146 |
JP: adding to this structure the instance being solved, the engine |
147 |
that we're solving it with, and a struct containing the output function ptrs |
148 |
*/ |
149 |
struct IntegratorSystemStruct{ |
150 |
struct Instance *instance; /**< not sure if this one is really necessary... -- JP */ |
151 |
slv_system_t system; /**< the system that we're integrating in ASCEND */ |
152 |
IntegratorEngine engine; /**< enum containing the ID of the integrator engine we're using */ |
153 |
IntegratorReporter *reporter;/**< functions for reporting integration results */ |
154 |
SampleList *samples; /**< pointer to the list of samples. we *don't own* this **/ |
155 |
void *enginedata; /**< space where the integrator engine can store stuff */ |
156 |
void *clientdata; /**< any stuff that the GUI/CLI needs to associate with this */ |
157 |
|
158 |
int nstates; /**< was a local global in integrator.c, moved it here. */ |
159 |
int nderivs; /**< ditto, as for nstates */ |
160 |
|
161 |
/* temporaries for build. these elements will be NULL to clients */ |
162 |
struct gl_list_t *indepvars; /**< all apparent independent vars */ |
163 |
struct gl_list_t *dynvars; /**< all state and deriv instances plus indices */ |
164 |
struct gl_list_t *obslist; /**< observation instance plus indices */ |
165 |
struct gl_list_t *states; /**< ordered list of state variables and indices*/ |
166 |
struct gl_list_t *derivs; /**< ordered list of derivative (ydot) insts */ |
167 |
|
168 |
/* persistent, these elements are valid to C clients. */ |
169 |
struct var_variable *x; /**< independent variable */ |
170 |
struct var_variable **y; /**< array form of states */ |
171 |
struct var_variable **ydot; /**< array form of derivatives */ |
172 |
struct var_variable **obs; /**< array form of observed variables */ |
173 |
long *y_id; /**< array form of y/ydot user indices */ |
174 |
long *obs_id; /**< array form of obs user indices */ |
175 |
long n_y; |
176 |
long n_obs; |
177 |
long currentstep; /**< current step number (also @see integrator_getnsamples) */ |
178 |
int ycount; /**< number of external references to y */ |
179 |
int ydotcount; /**< number of external references to ydot */ |
180 |
int obscount; /**< number of external references to obs */ |
181 |
int maxsubsteps; /**< most steps between mesh poins */ |
182 |
double stepzero; /**< initial step length, SI units. */ |
183 |
double minstep; /**< shortest step length, SI units. */ |
184 |
double maxstep; /**< longest step length, SI units. */ |
185 |
}; |
186 |
|
187 |
typedef struct IntegratorSystemStruct IntegratorSystem; |
188 |
|
189 |
/*------------------------------------------------------------------------------ |
190 |
PUBLIC INTERFACE (for use by the GUI/CLI)] |
191 |
*/ |
192 |
|
193 |
ASC_DLLSPEC(IntegratorSystem *) integrator_new(slv_system_t sys, struct Instance *inst); |
194 |
|
195 |
ASC_DLLSPEC(int) integrator_analyse(IntegratorSystem *blsys); |
196 |
|
197 |
ASC_DLLSPEC(int) integrator_solve(IntegratorSystem *blsys, long i0, long i1); |
198 |
/**< |
199 |
Takes the type of integrator and sets up the global variables into the |
200 |
current integration instance. |
201 |
|
202 |
@return 1 on success, 0 on failure. |
203 |
*/ |
204 |
|
205 |
ASC_DLLSPEC(void) integrator_free(IntegratorSystem *blsys); |
206 |
/**< |
207 |
Deallocates any memory used and sets all integration global points to NULL. |
208 |
*/ |
209 |
|
210 |
ASC_DLLSPEC(int) integrator_set_engine(IntegratorSystem *blsys, IntegratorEngine engine); |
211 |
/** |
212 |
Sets the engine for this integrator. Checks that the integrator can be used |
213 |
on the given system. |
214 |
|
215 |
@return 1 on success, if engine is compatible with the system being integrated. |
216 |
*/ |
217 |
|
218 |
ASC_DLLSPEC(IntegratorEngine) integrator_get_engine(const IntegratorSystem *blsys); |
219 |
/** |
220 |
Returns the engine (ID) selected for use in this IntegratorSystem (may return |
221 |
INTEG_UNKNOWN if none or invalid setting). |
222 |
*/ |
223 |
|
224 |
ASC_DLLSPEC(int) integrator_set_reporter(IntegratorSystem *blsys |
225 |
, IntegratorReporter *r); |
226 |
/**< |
227 |
Use this function from your interface to tell the integrator how it needs |
228 |
to make its output. You need to pass in a struct containing the required |
229 |
function pointers. (We will wrap IntegratorReporter and access it from |
230 |
Python, eventually) |
231 |
*/ |
232 |
|
233 |
ASC_DLLSPEC(int) integrator_find_indep_var(IntegratorSystem *blsys); |
234 |
/**< |
235 |
Attempt to locate the independent variable. It will be stored in the blsys |
236 |
structure if found. |
237 |
@return 1 if found, 0 if not found. |
238 |
*/ |
239 |
|
240 |
extern int integrator_checkstatus(slv_status_t status); |
241 |
/**< |
242 |
* Takes a status and checks for convergence. If converged then return 1 |
243 |
* else return 0 and print error message if appropriate. |
244 |
*/ |
245 |
|
246 |
ASC_DLLSPEC(void) integrator_set_samples(IntegratorSystem *blsys, SampleList *samples); |
247 |
/**< |
248 |
Sets values of time samples to the values given (ns of them) and |
249 |
keeps both the dim pointer and vector given. The vector and dimp |
250 |
given may be freed by calling this again, but xsamples owns |
251 |
them until then. If called with ns<1 or values==NULL, will free any |
252 |
previously captured values/dimp. If called with dimp==NULL, will |
253 |
assume WildDimension. Don't call this with a dimp we can't free later. |
254 |
Return is 1 if for some reason we can't set as expected, 0 otherwise. |
255 |
*/ |
256 |
|
257 |
ASC_DLLSPEC(void) integrator_set_stepzero(IntegratorSystem *blsys, double); |
258 |
ASC_DLLSPEC(double) integrator_get_stepzero(IntegratorSystem *blsys); |
259 |
/**< |
260 |
Returns the length of the initial step user specified, |
261 |
or 0.0 if none was set. |
262 |
*/ |
263 |
|
264 |
ASC_DLLSPEC(void) integrator_set_minstep(IntegratorSystem *blsys, double); |
265 |
ASC_DLLSPEC(double) integrator_get_minstep(IntegratorSystem *blsys); |
266 |
/**< |
267 |
Returns the length of the longest allowable step. |
268 |
or 0.0 if none was set by user. |
269 |
*/ |
270 |
|
271 |
ASC_DLLSPEC(void) integrator_set_maxstep(IntegratorSystem *blsys, double); |
272 |
ASC_DLLSPEC(double) integrator_get_maxstep(IntegratorSystem *blsys); |
273 |
/**< |
274 |
Returns the length of the shortest allowable step. |
275 |
or 0.0 if none was set by user. |
276 |
*/ |
277 |
|
278 |
ASC_DLLSPEC(void) integrator_set_maxsubsteps(IntegratorSystem *blsys, int); |
279 |
ASC_DLLSPEC(int) integrator_get_maxsubsteps(IntegratorSystem *blsys); |
280 |
/**< |
281 |
Returns the most internal steps allowed between |
282 |
two time samples, or 0 if none was set by user. |
283 |
*/ |
284 |
|
285 |
ASC_DLLSPEC(long) integrator_getnsamples(IntegratorSystem *blsys); |
286 |
/**< |
287 |
Returns the number of values currently stored in xsamples. |
288 |
*/ |
289 |
|
290 |
ASC_DLLSPEC(long) integrator_getcurrentstep(IntegratorSystem *blsys); |
291 |
/**< |
292 |
Returns the current step number (blsys->currentstep). Should be reset to |
293 |
zero inside your solver's integrator_*_solve method. |
294 |
*/ |
295 |
|
296 |
extern double integrator_getsample(IntegratorSystem *blsys, long i); |
297 |
/**< The following functions fetch and set parts with names specific to |
298 |
the type definitions in ivp.lib, the ASCEND initial value problem |
299 |
model file. They are for use by any ivp solver interface. |
300 |
All names are given at the scope of ivp. |
301 |
|
302 |
Returns the value stored in xsamples[i]. Will whine if |
303 |
if xsample[i] does not exist. No, there is no way to get |
304 |
back the pointer to the xsamples vector. |
305 |
*/ |
306 |
|
307 |
extern void integrator_setsample(IntegratorSystem *blsys, long i, double val); |
308 |
/**< |
309 |
Sets the value stored in xsamples[i] to val. Will whine if |
310 |
if xsample[i] does not exist. No, there is no way to get |
311 |
back the pointer to the xsamples vector. |
312 |
*/ |
313 |
|
314 |
extern const dim_type *integrator_getsampledim(IntegratorSystem *blsys); |
315 |
/**< |
316 |
Returns a pointer to the dimensionality of the samples currently |
317 |
stored, or NULL if none are stored. Do not free this pointer: we |
318 |
own it. |
319 |
*/ |
320 |
|
321 |
/*------------------------------------------------------------------------------ |
322 |
BACKEND INTERFACE (for use by the integrator engine) |
323 |
*/ |
324 |
|
325 |
/* Problem size parameters. */ |
326 |
|
327 |
/* Parts of type definition derivatives refinements. */ |
328 |
|
329 |
ASC_DLLSPEC(double) integrator_get_t(IntegratorSystem *blsys); |
330 |
/**< |
331 |
Gets value of the independent variable value from the ASCEND model |
332 |
(retrieves the value from the ASCEND compiler's instance hierarchy) |
333 |
*/ |
334 |
|
335 |
extern void integrator_set_t(IntegratorSystem *blsys, double value); |
336 |
/**< |
337 |
Sets value of the independent variable in the model (inserts the value in |
338 |
the correct place in the compiler's instance hierarchy via the var_variable |
339 |
API). |
340 |
*/ |
341 |
|
342 |
extern double *integrator_get_y(IntegratorSystem *blsys, double *vector); |
343 |
/**< |
344 |
Gets the value of vector 'y' from the model (what 'y' is depends on your |
345 |
particular integrator). |
346 |
|
347 |
If vector given is NULL, allocates vector, which the caller then owns. |
348 |
Vector, if given, should be IntegGet_d_neq()+1 long. |
349 |
*/ |
350 |
|
351 |
extern void integrator_set_y(IntegratorSystem *blsys, double *vector); |
352 |
/**< |
353 |
Sets d.y[] to values in vector. |
354 |
*/ |
355 |
|
356 |
extern double *integrator_get_ydot(IntegratorSystem *blsys, double *vector); |
357 |
/**< |
358 |
Returns the vector ydot (derivatives of the 'states') |
359 |
|
360 |
If vector given is NULL, allocates vector, which the caller then owns. |
361 |
Vector, if given, should be IntegGet_d_neq()+1 long. |
362 |
*/ |
363 |
|
364 |
extern void integrator_set_ydot(IntegratorSystem *blsys, double *vector); |
365 |
/**< |
366 |
Sets d.dydx[] to values in vector. |
367 |
*/ |
368 |
|
369 |
ASC_DLLSPEC(double *) integrator_get_observations(IntegratorSystem *blsys, double *vector); |
370 |
/**< |
371 |
Returns the vector d.obs. |
372 |
Vector should be of sufficient length (g_intginst_d_n_obs+1). |
373 |
If NULL is given a vector is allocated and is the caller's |
374 |
responsibility to deallocate. |
375 |
*/ |
376 |
|
377 |
ASC_DLLSPEC(struct var_variable *) integrator_get_observed_var(IntegratorSystem *blsys, const long i); |
378 |
/**< |
379 |
Returns the var_variable contained in the ith position in the observed variable list. |
380 |
*/ |
381 |
|
382 |
ASC_DLLSPEC(struct var_variable *) integrator_get_independent_var(IntegratorSystem *blsys); |
383 |
/**< |
384 |
Return a pointer to the variable identified as the independent variable. |
385 |
*/ |
386 |
|
387 |
/*------------------------------- |
388 |
Stuff to facilitate output to the interface |
389 |
*/ |
390 |
|
391 |
/** |
392 |
This call should be used to get the file streams ready, output column |
393 |
headings etc. |
394 |
*/ |
395 |
extern int integrator_output_init(IntegratorSystem *blsys); |
396 |
|
397 |
/** |
398 |
This call should be used to output the immediate integration results from |
399 |
a single timestep. For an ODE solver, this means just the 'y' values but |
400 |
perhaps not the values of the algebraic varaibles, which must be calculated |
401 |
separately. They'll be stored in the Integ_system_t somehow (not yet |
402 |
known how) |
403 |
*/ |
404 |
extern int integrator_output_write(IntegratorSystem *blsys); |
405 |
|
406 |
/** |
407 |
Write out the 'observed values' for the integration. In the case of ODE |
408 |
integration, we assume that the values of all the algabraic variables are |
409 |
also now calculated. |
410 |
*/ |
411 |
extern int integrator_output_write_obs(IntegratorSystem *blsys); |
412 |
|
413 |
/** |
414 |
This call will close file stream and perhaps perform some kind of |
415 |
user notification or screen update, etc. |
416 |
*/ |
417 |
extern int integrator_output_close(IntegratorSystem *blsys); |
418 |
|
419 |
#endif |