/[ascend]/trunk/base/generic/solver/integrator.h
ViewVC logotype

Contents of /trunk/base/generic/solver/integrator.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 945 - (show annotations) (download) (as text)
Sat Nov 25 12:41:03 2006 UTC (14 years ago) by johnpye
File MIME type: text/x-chdr
File size: 16495 byte(s)
Fixed a bug in slv_param_char.
More work on IDA ongoing.
1 /* ASCEND modelling environment
2 Copyright (C) 2006 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330,
17 Boston, MA 02111-1307, USA.
18 *//** @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 slv_parameters_t params; /**< structure containing parameters applicable to this Integrator */
159
160 int nstates; /**< was a local global in integrator.c, moved it here. */
161 int nderivs; /**< ditto, as for nstates */
162
163 /* temporaries for build. these elements will be NULL to clients */
164 struct gl_list_t *indepvars; /**< all apparent independent vars */
165 struct gl_list_t *dynvars; /**< all state and deriv instances plus indices */
166 struct gl_list_t *obslist; /**< observation instance plus indices */
167 struct gl_list_t *states; /**< ordered list of state variables and indices*/
168 struct gl_list_t *derivs; /**< ordered list of derivative (ydot) insts */
169
170 /* persistent, these elements are valid to C clients. */
171 struct var_variable *x; /**< independent variable */
172 struct var_variable **y; /**< array form of states */
173 struct var_variable **ydot; /**< array form of derivatives */
174 struct var_variable **obs; /**< array form of observed variables */
175 long *y_id; /**< array form of y/ydot user indices, for DAEs we use negatives here for derivative vars */
176 long *obs_id; /**< array form of obs user indices */
177 long n_y;
178 long n_obs;
179 long currentstep; /**< current step number (also @see integrator_getnsamples) */
180
181 int ycount; /**< number of external references to y */
182 int ydotcount; /**< number of external references to ydot */
183 int obscount; /**< number of external references to obs */
184
185 /** @TODO move the following to the 'params' structure? Or maybe better not to? */
186 int maxsubsteps; /**< most steps between mesh poins */
187 double stepzero; /**< initial step length, SI units. */
188 double minstep; /**< shortest step length, SI units. */
189 double maxstep; /**< longest step length, SI units. */
190 };
191
192 typedef struct IntegratorSystemStruct IntegratorSystem;
193
194 /*------------------------------------------------------------------------------
195 PUBLIC INTERFACE (for use by the GUI/CLI)]
196 */
197
198 ASC_DLLSPEC(IntegratorSystem *) integrator_new(slv_system_t sys, struct Instance *inst);
199
200 ASC_DLLSPEC(int) integrator_analyse(IntegratorSystem *blsys);
201
202 ASC_DLLSPEC(int) integrator_solve(IntegratorSystem *blsys, long i0, long i1);
203 /**<
204 Takes the type of integrator and sets up the global variables into the
205 current integration instance.
206
207 @return 1 on success, 0 on failure.
208 */
209
210 ASC_DLLSPEC(void) integrator_free(IntegratorSystem *blsys);
211 /**<
212 Deallocates any memory used and sets all integration global points to NULL.
213 */
214
215 ASC_DLLSPEC(int) integrator_set_engine(IntegratorSystem *blsys, IntegratorEngine engine);
216 /**<
217 Sets the engine for this integrator. Checks that the integrator can be used
218 on the given system.
219
220 @return 0 on success,
221 1 if unknown engine,
222 2 if invalid engine (not suitable for this present problem)
223 */
224
225 ASC_DLLSPEC(IntegratorEngine) integrator_get_engine(const IntegratorSystem *blsys);
226 /**<
227 Returns the engine (ID) selected for use in this IntegratorSystem (may return
228 INTEG_UNKNOWN if none or invalid setting).
229 */
230
231 typedef int IntegratorParamsDefaultFn(IntegratorSystem *blsys);
232 /**<
233 Integrators must provide a function like this that can be used to retrieve
234 the default set of parameters.
235 */
236
237 ASC_DLLSPEC(int) integrator_params_get(const IntegratorSystem *blsys, slv_parameters_t *parameters);
238 /**<
239 Copies the current set of Integrator parameters into the indicated parameter set, 'parameters'.
240
241 @TODO test these, not sure if the memory copy stuff is right or not.
242
243 @return 0 on success
244 */
245
246 ASC_DLLSPEC(int) integrator_params_set(IntegratorSystem *blsys, const slv_parameters_t *parameters);
247 /**<
248 Copies the the given parameter set into the Integrator's data structure,
249 which immediately makes the new values available to the integrator engine.
250
251 @TODO test these, not sure if the memory copy stuff is right or not.
252
253 @return 0 on success
254 */
255
256 ASC_DLLSPEC(int) integrator_set_reporter(IntegratorSystem *blsys
257 , IntegratorReporter *r);
258 /**<
259 Use this function from your interface to tell the integrator how it needs
260 to make its output. You need to pass in a struct containing the required
261 function pointers. (We will wrap IntegratorReporter and access it from
262 Python, eventually)
263 */
264
265 ASC_DLLSPEC(int) integrator_find_indep_var(IntegratorSystem *blsys);
266 /**<
267 Attempt to locate the independent variable. It will be stored in the blsys
268 structure if found.
269 @return 1 if found, 0 if not found.
270 */
271
272 extern int integrator_checkstatus(slv_status_t status);
273 /**<
274 * Takes a status and checks for convergence. If converged then return 1
275 * else return 0 and print error message if appropriate.
276 */
277
278 ASC_DLLSPEC(void) integrator_set_samples(IntegratorSystem *blsys, SampleList *samples);
279 /**<
280 Sets values of time samples to the values given (ns of them) and
281 keeps both the dim pointer and vector given. The vector and dimp
282 given may be freed by calling this again, but xsamples owns
283 them until then. If called with ns<1 or values==NULL, will free any
284 previously captured values/dimp. If called with dimp==NULL, will
285 assume WildDimension. Don't call this with a dimp we can't free later.
286 Return is 1 if for some reason we can't set as expected, 0 otherwise.
287 */
288
289 ASC_DLLSPEC(void) integrator_set_stepzero(IntegratorSystem *blsys, double);
290 ASC_DLLSPEC(double) integrator_get_stepzero(IntegratorSystem *blsys);
291 /**<
292 Returns the length of the initial step user specified,
293 or 0.0 if none was set.
294 */
295
296 ASC_DLLSPEC(void) integrator_set_minstep(IntegratorSystem *blsys, double);
297 ASC_DLLSPEC(double) integrator_get_minstep(IntegratorSystem *blsys);
298 /**<
299 Returns the length of the longest allowable step.
300 or 0.0 if none was set by user.
301 */
302
303 ASC_DLLSPEC(void) integrator_set_maxstep(IntegratorSystem *blsys, double);
304 ASC_DLLSPEC(double) integrator_get_maxstep(IntegratorSystem *blsys);
305 /**<
306 Returns the length of the shortest allowable step.
307 or 0.0 if none was set by user.
308 */
309
310 ASC_DLLSPEC(void) integrator_set_maxsubsteps(IntegratorSystem *blsys, int);
311 ASC_DLLSPEC(int) integrator_get_maxsubsteps(IntegratorSystem *blsys);
312 /**<
313 Returns the most internal steps allowed between
314 two time samples, or 0 if none was set by user.
315 */
316
317 ASC_DLLSPEC(long) integrator_getnsamples(IntegratorSystem *blsys);
318 /**<
319 Returns the number of values currently stored in xsamples.
320 */
321
322 ASC_DLLSPEC(long) integrator_getcurrentstep(IntegratorSystem *blsys);
323 /**<
324 Returns the current step number (blsys->currentstep). Should be reset to
325 zero inside your solver's integrator_*_solve method.
326 */
327
328 extern double integrator_getsample(IntegratorSystem *blsys, long i);
329 /**< The following functions fetch and set parts with names specific to
330 the type definitions in ivp.lib, the ASCEND initial value problem
331 model file. They are for use by any ivp solver interface.
332 All names are given at the scope of ivp.
333
334 Returns the value stored in xsamples[i]. Will whine if
335 if xsample[i] does not exist. No, there is no way to get
336 back the pointer to the xsamples vector.
337 */
338
339 extern void integrator_setsample(IntegratorSystem *blsys, long i, double val);
340 /**<
341 Sets the value stored in xsamples[i] to val. Will whine if
342 if xsample[i] does not exist. No, there is no way to get
343 back the pointer to the xsamples vector.
344 */
345
346 extern const dim_type *integrator_getsampledim(IntegratorSystem *blsys);
347 /**<
348 Returns a pointer to the dimensionality of the samples currently
349 stored, or NULL if none are stored. Do not free this pointer: we
350 own it.
351 */
352
353 /*------------------------------------------------------------------------------
354 BACKEND INTERFACE (for use by the integrator engine)
355 */
356
357 /* Problem size parameters. */
358
359 /* Parts of type definition derivatives refinements. */
360
361 ASC_DLLSPEC(double) integrator_get_t(IntegratorSystem *blsys);
362 /**<
363 Gets value of the independent variable value from the ASCEND model
364 (retrieves the value from the ASCEND compiler's instance hierarchy)
365 */
366
367 extern void integrator_set_t(IntegratorSystem *blsys, double value);
368 /**<
369 Sets value of the independent variable in the model (inserts the value in
370 the correct place in the compiler's instance hierarchy via the var_variable
371 API).
372 */
373
374 extern double *integrator_get_y(IntegratorSystem *blsys, double *vector);
375 /**<
376 Gets the value of vector 'y' from the model (what 'y' is depends on your
377 particular integrator).
378
379 If vector given is NULL, allocates vector, which the caller then owns.
380 Vector, if given, should be IntegGet_d_neq()+1 long.
381 */
382
383 extern void integrator_set_y(IntegratorSystem *blsys, double *vector);
384 /**<
385 Sets d.y[] to values in vector.
386 */
387
388 extern double *integrator_get_ydot(IntegratorSystem *blsys, double *vector);
389 /**<
390 Returns the vector ydot (derivatives of the 'states')
391
392 If vector given is NULL, allocates vector, which the caller then owns.
393 Vector, if given, should be IntegGet_d_neq()+1 long.
394 */
395
396 extern void integrator_set_ydot(IntegratorSystem *blsys, double *vector);
397 /**<
398 Sets d.dydx[] to values in vector.
399 */
400
401 ASC_DLLSPEC(double *) integrator_get_observations(IntegratorSystem *blsys, double *vector);
402 /**<
403 Returns the vector d.obs.
404 Vector should be of sufficient length (g_intginst_d_n_obs+1).
405 If NULL is given a vector is allocated and is the caller's
406 responsibility to deallocate.
407 */
408
409 ASC_DLLSPEC(struct var_variable *) integrator_get_observed_var(IntegratorSystem *blsys, const long i);
410 /**<
411 Returns the var_variable contained in the ith position in the observed variable list.
412 */
413
414 ASC_DLLSPEC(struct var_variable *) integrator_get_independent_var(IntegratorSystem *blsys);
415 /**<
416 Return a pointer to the variable identified as the independent variable.
417 */
418
419 extern double *integrator_get_atol(IntegratorSystem *blsys, double *vector);
420 /**<
421 Return abolute tolerance values for the 'ode_atol' values in the MODEL.
422 */
423
424 /*-------------------------------
425 Stuff to facilitate output to the interface
426 */
427
428 /**
429 This call should be used to get the file streams ready, output column
430 headings etc.
431 */
432 extern int integrator_output_init(IntegratorSystem *blsys);
433
434 /**
435 This call should be used to output the immediate integration results from
436 a single timestep. For an ODE solver, this means just the 'y' values but
437 perhaps not the values of the algebraic varaibles, which must be calculated
438 separately. They'll be stored in the Integ_system_t somehow (not yet
439 known how)
440
441 @return 1 on success, return 0 on failure (integration will be cancelled)
442 */
443 extern int integrator_output_write(IntegratorSystem *blsys);
444
445 /**
446 Write out the 'observed values' for the integration. In the case of ODE
447 integration, we assume that the values of all the algabraic variables are
448 also now calculated.
449 */
450 extern int integrator_output_write_obs(IntegratorSystem *blsys);
451
452 /**
453 This call will close file stream and perhaps perform some kind of
454 user notification or screen update, etc.
455 */
456 extern int integrator_output_close(IntegratorSystem *blsys);
457
458 #endif

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