/[ascend]/trunk/ascend/integrator/integrator.h
ViewVC logotype

Contents of /trunk/ascend/integrator/integrator.h

Parent Directory Parent Directory | Revision Log Revision Log


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

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