/[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 940 - (show annotations) (download) (as text)
Fri Nov 24 08:07:33 2006 UTC (13 years, 10 months ago) by johnpye
File MIME type: text/x-chdr
File size: 15284 byte(s)
Fixed 'testLSODE' in test.py, added 'IntegratorReporterNull' that makes no output (for testing purposes)
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, for DAEs we use negatives here for derivative vars */
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 @return 1 on success, return 0 on failure (integration will be cancelled)
405 */
406 extern int integrator_output_write(IntegratorSystem *blsys);
407
408 /**
409 Write out the 'observed values' for the integration. In the case of ODE
410 integration, we assume that the values of all the algabraic variables are
411 also now calculated.
412 */
413 extern int integrator_output_write_obs(IntegratorSystem *blsys);
414
415 /**
416 This call will close file stream and perhaps perform some kind of
417 user notification or screen update, etc.
418 */
419 extern int integrator_output_close(IntegratorSystem *blsys);
420
421 #endif

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