/[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 977 - (show annotations) (download) (as text)
Wed Dec 20 00:39:52 2006 UTC (12 years, 11 months ago) by johnpye
File MIME type: text/x-chdr
File size: 19023 byte(s)
Abstracted the internal integrator calls into a struct IntegratorInternals.
Fixed up compile-time list of integrators.
If IDA is not available, then 'INTEG_IDA' will not be defined.
Added ASC_ASSERT_RANGE for assertions x in [low,high).
Changed calling convention for integrator_get_engines().
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

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