/[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 1049 - (show annotations) (download) (as text)
Fri Jan 5 13:45:13 2007 UTC (16 years, 8 months ago) by johnpye
File MIME type: text/x-chdr
File size: 19053 byte(s)
Changed integrator_solve to return 0 on success
Changed integrator_checkstatus to return 0 on success
Some flow-through changes as a result of these.
Fixed problem with reporting of errors in do_method (gtkbrowser.py)
Expanded error msg in integrator.cpp
Fixed catching of Integrator C++ exceptions in integratorreporter.py.
Removed M.checkStructuralSingularity() in TestSteam.testdsgsat -- causes crash (why?)
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 0 on success, else error
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, and prints status info if the
354 system is not converged.
355 @return 0 if status says system has converged, non-zero otherwise.
356 */
357
358 ASC_DLLSPEC(void) integrator_set_samples(IntegratorSystem *blsys, SampleList *samples);
359 /**<
360 Sets values of time samples to the values given (ns of them) and
361 keeps both the dim pointer and vector given. The vector and dimp
362 given may be freed by calling this again, but xsamples owns
363 them until then. If called with ns<1 or values==NULL, will free any
364 previously captured values/dimp. If called with dimp==NULL, will
365 assume WildDimension. Don't call this with a dimp we can't free later.
366 Return is 1 if for some reason we can't set as expected, 0 otherwise.
367 */
368
369 ASC_DLLSPEC(void) integrator_set_stepzero(IntegratorSystem *blsys, double);
370 ASC_DLLSPEC(double) integrator_get_stepzero(IntegratorSystem *blsys);
371 /**<
372 Returns the length of the initial step user specified,
373 or 0.0 if none was set.
374 */
375
376 ASC_DLLSPEC(void) integrator_set_minstep(IntegratorSystem *blsys, double);
377 ASC_DLLSPEC(double) integrator_get_minstep(IntegratorSystem *blsys);
378 /**<
379 Returns the length of the longest allowable step.
380 or 0.0 if none was set by user.
381 */
382
383 ASC_DLLSPEC(void) integrator_set_maxstep(IntegratorSystem *blsys, double);
384 ASC_DLLSPEC(double) integrator_get_maxstep(IntegratorSystem *blsys);
385 /**<
386 Returns the length of the shortest allowable step.
387 or 0.0 if none was set by user.
388 */
389
390 ASC_DLLSPEC(void) integrator_set_maxsubsteps(IntegratorSystem *blsys, int);
391 ASC_DLLSPEC(int) integrator_get_maxsubsteps(IntegratorSystem *blsys);
392 /**<
393 Returns the most internal steps allowed between
394 two time samples, or 0 if none was set by user.
395 */
396
397 ASC_DLLSPEC(long) integrator_getnsamples(IntegratorSystem *blsys);
398 /**<
399 Returns the number of values currently stored in xsamples.
400 */
401
402 ASC_DLLSPEC(long) integrator_getcurrentstep(IntegratorSystem *blsys);
403 /**<
404 Returns the current step number (blsys->currentstep). Should be reset to
405 zero inside your solver's integrator_*_solve method.
406 */
407
408 extern double integrator_getsample(IntegratorSystem *blsys, long i);
409 /**< The following functions fetch and set parts with names specific to
410 the type definitions in ivp.lib, the ASCEND initial value problem
411 model file. They are for use by any ivp solver interface.
412 All names are given at the scope of ivp.
413
414 Returns the value stored in xsamples[i]. Will whine if
415 if xsample[i] does not exist. No, there is no way to get
416 back the pointer to the xsamples vector.
417 */
418
419 extern void integrator_setsample(IntegratorSystem *blsys, long i, double val);
420 /**<
421 Sets the value stored in xsamples[i] to val. Will whine if
422 if xsample[i] does not exist. No, there is no way to get
423 back the pointer to the xsamples vector.
424 */
425
426 extern const dim_type *integrator_getsampledim(IntegratorSystem *blsys);
427 /**<
428 Returns a pointer to the dimensionality of the samples currently
429 stored, or NULL if none are stored. Do not free this pointer: we
430 own it.
431 */
432
433 /*------------------------------------------------------------------------------
434 BACKEND INTERFACE (for use by the integrator engine)
435 */
436
437 /* Problem size parameters. */
438
439 /* Parts of type definition derivatives refinements. */
440
441 ASC_DLLSPEC(double) integrator_get_t(IntegratorSystem *blsys);
442 /**<
443 Gets value of the independent variable value from the ASCEND model
444 (retrieves the value from the ASCEND compiler's instance hierarchy)
445 */
446
447 extern void integrator_set_t(IntegratorSystem *blsys, double value);
448 /**<
449 Sets value of the independent variable in the model (inserts the value in
450 the correct place in the compiler's instance hierarchy via the var_variable
451 API).
452 */
453
454 extern double *integrator_get_y(IntegratorSystem *blsys, double *vector);
455 /**<
456 Gets the value of vector 'y' from the model (what 'y' is depends on your
457 particular integrator).
458
459 If vector given is NULL, allocates vector, which the caller then owns.
460 Vector, if given, should be IntegGet_d_neq()+1 long.
461 */
462
463 extern void integrator_set_y(IntegratorSystem *blsys, double *vector);
464 /**<
465 Sets d.y[] to values in vector.
466 */
467
468 extern double *integrator_get_ydot(IntegratorSystem *blsys, double *vector);
469 /**<
470 Returns the vector ydot (derivatives of the 'states')
471
472 If vector given is NULL, allocates vector, which the caller then owns.
473 Vector, if given, should be IntegGet_d_neq()+1 long.
474 */
475
476 extern void integrator_set_ydot(IntegratorSystem *blsys, double *vector);
477 /**<
478 Sets d.dydx[] to values in vector.
479 */
480
481 ASC_DLLSPEC(double *) integrator_get_observations(IntegratorSystem *blsys, double *vector);
482 /**<
483 Returns the vector d.obs.
484 Vector should be of sufficient length (g_intginst_d_n_obs+1).
485 If NULL is given a vector is allocated and is the caller's
486 responsibility to deallocate.
487 */
488
489 ASC_DLLSPEC(struct var_variable *) integrator_get_observed_var(IntegratorSystem *blsys, const long i);
490 /**<
491 Returns the var_variable contained in the ith position in the observed variable list.
492 */
493
494 ASC_DLLSPEC(struct var_variable *) integrator_get_independent_var(IntegratorSystem *blsys);
495 /**<
496 Return a pointer to the variable identified as the independent variable.
497 */
498
499 extern double *integrator_get_atol(IntegratorSystem *blsys, double *vector);
500 /**<
501 Return abolute tolerance values for the 'ode_atol' values in the MODEL.
502 */
503
504 /*-------------------------------
505 Stuff to facilitate output to the interface
506 */
507
508 /**
509 This call should be used to get the file streams ready, output column
510 headings etc.
511 */
512 extern int integrator_output_init(IntegratorSystem *blsys);
513
514 /**
515 This call should be used to output the immediate integration results from
516 a single timestep. For an ODE solver, this means just the 'y' values but
517 perhaps not the values of the algebraic varaibles, which must be calculated
518 separately. They'll be stored in the Integ_system_t somehow (not yet
519 known how)
520
521 @return 1 on success, return 0 on failure (integration will be cancelled)
522 */
523 extern int integrator_output_write(IntegratorSystem *blsys);
524
525 /**
526 Write out the 'observed values' for the integration. In the case of ODE
527 integration, we assume that the values of all the algabraic variables are
528 also now calculated.
529 */
530 extern int integrator_output_write_obs(IntegratorSystem *blsys);
531
532 /**
533 This call will close file stream and perhaps perform some kind of
534 user notification or screen update, etc.
535 */
536 extern int integrator_output_close(IntegratorSystem *blsys);
537
538 #endif

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