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

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