| 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 |
slv_parameters_t params; /**< structure containing parameters applicable to this Integrator */ |
| 159 |
|
| 160 |
int nstates; /**< was a local global in integrator.c, moved it here. */ |
| 161 |
int nderivs; /**< ditto, as for nstates */ |
| 162 |
|
| 163 |
/* temporaries for build. these elements will be NULL to clients */ |
| 164 |
struct gl_list_t *indepvars; /**< all apparent independent vars */ |
| 165 |
struct gl_list_t *dynvars; /**< all state and deriv instances plus indices */ |
| 166 |
struct gl_list_t *obslist; /**< observation instance plus indices */ |
| 167 |
struct gl_list_t *states; /**< ordered list of state variables and indices*/ |
| 168 |
struct gl_list_t *derivs; /**< ordered list of derivative (ydot) insts */ |
| 169 |
|
| 170 |
/* persistent, these elements are valid to C clients. */ |
| 171 |
struct var_variable *x; /**< independent variable */ |
| 172 |
struct var_variable **y; /**< array form of states */ |
| 173 |
struct var_variable **ydot; /**< array form of derivatives */ |
| 174 |
struct var_variable **obs; /**< array form of observed variables */ |
| 175 |
long *y_id; /**< array form of y/ydot user indices, for DAEs we use negatives here for derivative vars */ |
| 176 |
long *obs_id; /**< array form of obs user indices */ |
| 177 |
long n_y; |
| 178 |
long n_obs; |
| 179 |
long currentstep; /**< current step number (also @see integrator_getnsamples) */ |
| 180 |
|
| 181 |
int ycount; /**< number of external references to y */ |
| 182 |
int ydotcount; /**< number of external references to ydot */ |
| 183 |
int obscount; /**< number of external references to obs */ |
| 184 |
|
| 185 |
/** @TODO move the following to the 'params' structure? Or maybe better not to? */ |
| 186 |
int maxsubsteps; /**< most steps between mesh poins */ |
| 187 |
double stepzero; /**< initial step length, SI units. */ |
| 188 |
double minstep; /**< shortest step length, SI units. */ |
| 189 |
double maxstep; /**< longest step length, SI units. */ |
| 190 |
}; |
| 191 |
|
| 192 |
typedef struct IntegratorSystemStruct IntegratorSystem; |
| 193 |
|
| 194 |
/*------------------------------------------------------------------------------ |
| 195 |
PUBLIC INTERFACE (for use by the GUI/CLI)] |
| 196 |
*/ |
| 197 |
|
| 198 |
ASC_DLLSPEC(IntegratorSystem *) integrator_new(slv_system_t sys, struct Instance *inst); |
| 199 |
|
| 200 |
ASC_DLLSPEC(int) integrator_analyse(IntegratorSystem *blsys); |
| 201 |
|
| 202 |
ASC_DLLSPEC(int) integrator_solve(IntegratorSystem *blsys, long i0, long i1); |
| 203 |
/**< |
| 204 |
Takes the type of integrator and sets up the global variables into the |
| 205 |
current integration instance. |
| 206 |
|
| 207 |
@return 1 on success, 0 on failure. |
| 208 |
*/ |
| 209 |
|
| 210 |
ASC_DLLSPEC(void) integrator_free(IntegratorSystem *blsys); |
| 211 |
/**< |
| 212 |
Deallocates any memory used and sets all integration global points to NULL. |
| 213 |
*/ |
| 214 |
|
| 215 |
ASC_DLLSPEC(int) integrator_set_engine(IntegratorSystem *blsys, IntegratorEngine engine); |
| 216 |
/**< |
| 217 |
Sets the engine for this integrator. Checks that the integrator can be used |
| 218 |
on the given system. |
| 219 |
|
| 220 |
@return 0 on success, |
| 221 |
1 if unknown engine, |
| 222 |
2 if invalid engine (not suitable for this present problem) |
| 223 |
*/ |
| 224 |
|
| 225 |
ASC_DLLSPEC(IntegratorEngine) integrator_get_engine(const IntegratorSystem *blsys); |
| 226 |
/**< |
| 227 |
Returns the engine (ID) selected for use in this IntegratorSystem (may return |
| 228 |
INTEG_UNKNOWN if none or invalid setting). |
| 229 |
*/ |
| 230 |
|
| 231 |
typedef int IntegratorParamsDefaultFn(IntegratorSystem *blsys); |
| 232 |
/**< |
| 233 |
Integrators must provide a function like this that can be used to retrieve |
| 234 |
the default set of parameters. |
| 235 |
*/ |
| 236 |
|
| 237 |
ASC_DLLSPEC(int) integrator_params_get(const IntegratorSystem *blsys, slv_parameters_t *parameters); |
| 238 |
/**< |
| 239 |
Copies the current set of Integrator parameters into the indicated parameter set, 'parameters'. |
| 240 |
|
| 241 |
@TODO test these, not sure if the memory copy stuff is right or not. |
| 242 |
|
| 243 |
@return 0 on success |
| 244 |
*/ |
| 245 |
|
| 246 |
ASC_DLLSPEC(int) integrator_params_set(IntegratorSystem *blsys, const slv_parameters_t *parameters); |
| 247 |
/**< |
| 248 |
Copies the the given parameter set into the Integrator's data structure, |
| 249 |
which immediately makes the new values available to the integrator engine. |
| 250 |
|
| 251 |
@TODO test these, not sure if the memory copy stuff is right or not. |
| 252 |
|
| 253 |
@return 0 on success |
| 254 |
*/ |
| 255 |
|
| 256 |
ASC_DLLSPEC(int) integrator_set_reporter(IntegratorSystem *blsys |
| 257 |
, IntegratorReporter *r); |
| 258 |
/**< |
| 259 |
Use this function from your interface to tell the integrator how it needs |
| 260 |
to make its output. You need to pass in a struct containing the required |
| 261 |
function pointers. (We will wrap IntegratorReporter and access it from |
| 262 |
Python, eventually) |
| 263 |
*/ |
| 264 |
|
| 265 |
ASC_DLLSPEC(int) integrator_find_indep_var(IntegratorSystem *blsys); |
| 266 |
/**< |
| 267 |
Attempt to locate the independent variable. It will be stored in the blsys |
| 268 |
structure if found. |
| 269 |
@return 1 if found, 0 if not found. |
| 270 |
*/ |
| 271 |
|
| 272 |
extern int integrator_checkstatus(slv_status_t status); |
| 273 |
/**< |
| 274 |
* Takes a status and checks for convergence. If converged then return 1 |
| 275 |
* else return 0 and print error message if appropriate. |
| 276 |
*/ |
| 277 |
|
| 278 |
ASC_DLLSPEC(void) integrator_set_samples(IntegratorSystem *blsys, SampleList *samples); |
| 279 |
/**< |
| 280 |
Sets values of time samples to the values given (ns of them) and |
| 281 |
keeps both the dim pointer and vector given. The vector and dimp |
| 282 |
given may be freed by calling this again, but xsamples owns |
| 283 |
them until then. If called with ns<1 or values==NULL, will free any |
| 284 |
previously captured values/dimp. If called with dimp==NULL, will |
| 285 |
assume WildDimension. Don't call this with a dimp we can't free later. |
| 286 |
Return is 1 if for some reason we can't set as expected, 0 otherwise. |
| 287 |
*/ |
| 288 |
|
| 289 |
ASC_DLLSPEC(void) integrator_set_stepzero(IntegratorSystem *blsys, double); |
| 290 |
ASC_DLLSPEC(double) integrator_get_stepzero(IntegratorSystem *blsys); |
| 291 |
/**< |
| 292 |
Returns the length of the initial step user specified, |
| 293 |
or 0.0 if none was set. |
| 294 |
*/ |
| 295 |
|
| 296 |
ASC_DLLSPEC(void) integrator_set_minstep(IntegratorSystem *blsys, double); |
| 297 |
ASC_DLLSPEC(double) integrator_get_minstep(IntegratorSystem *blsys); |
| 298 |
/**< |
| 299 |
Returns the length of the longest allowable step. |
| 300 |
or 0.0 if none was set by user. |
| 301 |
*/ |
| 302 |
|
| 303 |
ASC_DLLSPEC(void) integrator_set_maxstep(IntegratorSystem *blsys, double); |
| 304 |
ASC_DLLSPEC(double) integrator_get_maxstep(IntegratorSystem *blsys); |
| 305 |
/**< |
| 306 |
Returns the length of the shortest allowable step. |
| 307 |
or 0.0 if none was set by user. |
| 308 |
*/ |
| 309 |
|
| 310 |
ASC_DLLSPEC(void) integrator_set_maxsubsteps(IntegratorSystem *blsys, int); |
| 311 |
ASC_DLLSPEC(int) integrator_get_maxsubsteps(IntegratorSystem *blsys); |
| 312 |
/**< |
| 313 |
Returns the most internal steps allowed between |
| 314 |
two time samples, or 0 if none was set by user. |
| 315 |
*/ |
| 316 |
|
| 317 |
ASC_DLLSPEC(long) integrator_getnsamples(IntegratorSystem *blsys); |
| 318 |
/**< |
| 319 |
Returns the number of values currently stored in xsamples. |
| 320 |
*/ |
| 321 |
|
| 322 |
ASC_DLLSPEC(long) integrator_getcurrentstep(IntegratorSystem *blsys); |
| 323 |
/**< |
| 324 |
Returns the current step number (blsys->currentstep). Should be reset to |
| 325 |
zero inside your solver's integrator_*_solve method. |
| 326 |
*/ |
| 327 |
|
| 328 |
extern double integrator_getsample(IntegratorSystem *blsys, long i); |
| 329 |
/**< The following functions fetch and set parts with names specific to |
| 330 |
the type definitions in ivp.lib, the ASCEND initial value problem |
| 331 |
model file. They are for use by any ivp solver interface. |
| 332 |
All names are given at the scope of ivp. |
| 333 |
|
| 334 |
Returns the value stored in xsamples[i]. Will whine if |
| 335 |
if xsample[i] does not exist. No, there is no way to get |
| 336 |
back the pointer to the xsamples vector. |
| 337 |
*/ |
| 338 |
|
| 339 |
extern void integrator_setsample(IntegratorSystem *blsys, long i, double val); |
| 340 |
/**< |
| 341 |
Sets the value stored in xsamples[i] to val. Will whine if |
| 342 |
if xsample[i] does not exist. No, there is no way to get |
| 343 |
back the pointer to the xsamples vector. |
| 344 |
*/ |
| 345 |
|
| 346 |
extern const dim_type *integrator_getsampledim(IntegratorSystem *blsys); |
| 347 |
/**< |
| 348 |
Returns a pointer to the dimensionality of the samples currently |
| 349 |
stored, or NULL if none are stored. Do not free this pointer: we |
| 350 |
own it. |
| 351 |
*/ |
| 352 |
|
| 353 |
/*------------------------------------------------------------------------------ |
| 354 |
BACKEND INTERFACE (for use by the integrator engine) |
| 355 |
*/ |
| 356 |
|
| 357 |
/* Problem size parameters. */ |
| 358 |
|
| 359 |
/* Parts of type definition derivatives refinements. */ |
| 360 |
|
| 361 |
ASC_DLLSPEC(double) integrator_get_t(IntegratorSystem *blsys); |
| 362 |
/**< |
| 363 |
Gets value of the independent variable value from the ASCEND model |
| 364 |
(retrieves the value from the ASCEND compiler's instance hierarchy) |
| 365 |
*/ |
| 366 |
|
| 367 |
extern void integrator_set_t(IntegratorSystem *blsys, double value); |
| 368 |
/**< |
| 369 |
Sets value of the independent variable in the model (inserts the value in |
| 370 |
the correct place in the compiler's instance hierarchy via the var_variable |
| 371 |
API). |
| 372 |
*/ |
| 373 |
|
| 374 |
extern double *integrator_get_y(IntegratorSystem *blsys, double *vector); |
| 375 |
/**< |
| 376 |
Gets the value of vector 'y' from the model (what 'y' is depends on your |
| 377 |
particular integrator). |
| 378 |
|
| 379 |
If vector given is NULL, allocates vector, which the caller then owns. |
| 380 |
Vector, if given, should be IntegGet_d_neq()+1 long. |
| 381 |
*/ |
| 382 |
|
| 383 |
extern void integrator_set_y(IntegratorSystem *blsys, double *vector); |
| 384 |
/**< |
| 385 |
Sets d.y[] to values in vector. |
| 386 |
*/ |
| 387 |
|
| 388 |
extern double *integrator_get_ydot(IntegratorSystem *blsys, double *vector); |
| 389 |
/**< |
| 390 |
Returns the vector ydot (derivatives of the 'states') |
| 391 |
|
| 392 |
If vector given is NULL, allocates vector, which the caller then owns. |
| 393 |
Vector, if given, should be IntegGet_d_neq()+1 long. |
| 394 |
*/ |
| 395 |
|
| 396 |
extern void integrator_set_ydot(IntegratorSystem *blsys, double *vector); |
| 397 |
/**< |
| 398 |
Sets d.dydx[] to values in vector. |
| 399 |
*/ |
| 400 |
|
| 401 |
ASC_DLLSPEC(double *) integrator_get_observations(IntegratorSystem *blsys, double *vector); |
| 402 |
/**< |
| 403 |
Returns the vector d.obs. |
| 404 |
Vector should be of sufficient length (g_intginst_d_n_obs+1). |
| 405 |
If NULL is given a vector is allocated and is the caller's |
| 406 |
responsibility to deallocate. |
| 407 |
*/ |
| 408 |
|
| 409 |
ASC_DLLSPEC(struct var_variable *) integrator_get_observed_var(IntegratorSystem *blsys, const long i); |
| 410 |
/**< |
| 411 |
Returns the var_variable contained in the ith position in the observed variable list. |
| 412 |
*/ |
| 413 |
|
| 414 |
ASC_DLLSPEC(struct var_variable *) integrator_get_independent_var(IntegratorSystem *blsys); |
| 415 |
/**< |
| 416 |
Return a pointer to the variable identified as the independent variable. |
| 417 |
*/ |
| 418 |
|
| 419 |
extern double *integrator_get_abstol(IntegratorSystem *blsys, double *vector); |
| 420 |
/**< |
| 421 |
Return abolute tolerance values for the 'ode_atol' values in the MODEL. |
| 422 |
*/ |
| 423 |
|
| 424 |
/*------------------------------- |
| 425 |
Stuff to facilitate output to the interface |
| 426 |
*/ |
| 427 |
|
| 428 |
/** |
| 429 |
This call should be used to get the file streams ready, output column |
| 430 |
headings etc. |
| 431 |
*/ |
| 432 |
extern int integrator_output_init(IntegratorSystem *blsys); |
| 433 |
|
| 434 |
/** |
| 435 |
This call should be used to output the immediate integration results from |
| 436 |
a single timestep. For an ODE solver, this means just the 'y' values but |
| 437 |
perhaps not the values of the algebraic varaibles, which must be calculated |
| 438 |
separately. They'll be stored in the Integ_system_t somehow (not yet |
| 439 |
known how) |
| 440 |
|
| 441 |
@return 1 on success, return 0 on failure (integration will be cancelled) |
| 442 |
*/ |
| 443 |
extern int integrator_output_write(IntegratorSystem *blsys); |
| 444 |
|
| 445 |
/** |
| 446 |
Write out the 'observed values' for the integration. In the case of ODE |
| 447 |
integration, we assume that the values of all the algabraic variables are |
| 448 |
also now calculated. |
| 449 |
*/ |
| 450 |
extern int integrator_output_write_obs(IntegratorSystem *blsys); |
| 451 |
|
| 452 |
/** |
| 453 |
This call will close file stream and perhaps perform some kind of |
| 454 |
user notification or screen update, etc. |
| 455 |
*/ |
| 456 |
extern int integrator_output_close(IntegratorSystem *blsys); |
| 457 |
|
| 458 |
#endif |
Parent Directory
| 
Revision Log