Integrator API

From ASCEND
Jump to: navigation, search

The Integrator API is the interface through which DAE and ODE solvers for dynamic modelling must communicate with ASCEND.

Process flow

This section describes the process by which an Integrator is instantiated and used in ASCEND, including some detail about the associate GUI actions, outputs, and so on.

Registration

An integrator is built as an external library that can be either loaded automatically by ASCEND, or manually using an IMPORT statement in the user's model code.

These external libraries are loaded using 'dlopen' and then a register function is called. This function should return a data structure containing a set of function pointers, an integrator ID number, and a string integrator name, as shown:

static const IntegratorInternals integrator_lsode_internals = {
	integrator_lsode_create
	,integrator_lsode_params_default
	,integrator_analyse_ode
	,integrator_lsode_solve
	,integrator_lsode_write_matrix
	,NULL /* debugfn */

	,integrator_lsode_free
	,INTEG_LSODE
	,"LSODE";
};

extern ASC_EXPORT int lsode_register(void){

	CONSOLE_DEBUG("Registering LSODE...");
	return integrator_register(&integrator_lsode_internals);
}

The integrator ID number (INTEG_LSODE in this case, defined as a number) is bit of a fossil; we aim to remove this ID number eventually. The integrator_register function checks that no integrator with the same name or ID has already been registered. If that's OK, the integrator is added to the list of available integrators.

Build

When the PyGTK GUI user clicks the integrate button, the system object (slv_system_t) is built if it hasn't already been built. This means that the system_build function (common with all problems, eg NLAs) is called, which in turn results in calls to analyze_make_problem. This last analyze_make_problem includes calls to analyze_make_master_lists and system_generate_diffvars.

In analyze_make_master_lists there is a 'visit' process that uses the function classify_instance to identify the following types of variables, and add them to some intermediate lists in the problem_t structure:

  • variables that need to be 'observed' when integrator output is reported
  • independent variables (hopefully only one, ode_type == -1)
  • differential/derivative variables (ode_type >= 1)
  • algebraic variables (ode_type == 0)

The function system_generate_diffvars is also called during analyze_make_problem uses a tricky sorting algorithm to arrange the diffvars into a suitable order, so that 'chains' of derivatives can be identified. These chains are split into separate ones for each variable set (eg {y, dydt, d2ydt2}), and stored in a big SolverDiffVarCollection object, which is ultimately stored in the system object (as the diffvars element).

Selection

The user must select which integrator engine they want to use to solve a problem (although a GUI can of course offer a default choice).

A list of all available (loaded) integrator enginers can be obtained with the function integrator_get_engines, which returns a list of IntegratorInternals from which all the details of the loaded integrators can be determined, including the names.

An integrator engine is selected by the user with the function integrator_set_engine, which uses the short string name for the engine as a parameter. If the integrator name exists in the engine list, integrator_free_engine is called on the previously-selected engine if necessary, then the internals field value is updated with the new IntegratorInternals pointer, and then integrator_create_engine function is called to allocate space as required.

Parameter and timestep input

Currently, integrators need some additional parameters that are not stored within the model itself. These include number of time-steps, duration of integration, etc.

During this stage, the PyGTK GUI makes a call to the function integrator_find_indep_var, to determine the current value of the independent variable, so that its value can be inserted into the parameter dialog box.

The user inputs number of timesteps and the spacing of the timesteps (linear, logarithmic,...), which is used to create a suitable SampleList (see ascend/integrator/samplelist.h) containing the values of the independent variable for which solution values are to be output/recorded. Actually, the stuff that assigns the SampleList contents is C++ code from pygtk/integrator.h.

Associated with the integrator is an IntegratorParamsDefaultFn, which will create the list of parameters that configures the integration engine. The integration engine provides the initial list; the GUI is then used by the user to make any adjustments to the predefined defaults.

We would like to migrate towards a solution that allows default values for these parameters to be specified within the model file, for easier model re-usability.

Analysis

When all the timestep (SampleList) information is specified and the integrator parameters are set, the integrator_analyse function is called. This function allows the integrator engine to do any necessary preparatory analysis of the model. This may include re-ordering the equations and variables for optimal solving, or creating internal data structures relating derivatives and their base variables.

A default analysis function, integrator_analyse_ode is provided (ascend/integrator/integrator.c) that is used by the LSODE and DOPRI5 integrators, but other integrators such as IDA can provide their own function (see integrator_ida_analyse in solvers/ida/idaanalyse.c.

In integrator_analyse_ode, the call integrator_visit_system_vars(sys,&integrator_ode_classify_var) is used to visit all of the variables in the system and classify any that are derivatives or variables for which derivatives are present. The function DynamicVarInfo works out the type of each variable (state/derivative/algebraic/independent) and its 'ode_id', then integrator_ode_classify_var does the job of building up suitable lists of each kind of variable for later use in the analysis process.

If analysis returns an error of any kind, that typically means that the user has configured their problem incorrectly. In this case, the PyGTK GUI opens a dialog and fills it with the output of the function integrator_debug. This function can be used to give suggestions to the user about what may need to be fixed (although the usual error reported mechanism can also be used).

See also IDA/Analysis: the IDA solver does things a bit differently.

Solution

Once the analysis has been completed, the PyGTK GUI creates an IntegratorReporterPython dialog that graphically shows the progress of timesteps as they proceed. The IntegratorReporterPython dialog makes a call to integrator_solve, after first setting up hook functions using integrator_set_reporter that allow intermediate results to be recorded.

integrator_solve makes a call to the specific integration routine provided by the selected integration engine. It is the job of this function to integrate up to each timestep listed in the SampleList and to make calls to the IntegrationReporter at each such point. Some integration libraries will have difficulty providing results for all requested timesteps; in this case, so don't assume that the observer output matches the SampleList in all cases.

In the PyGTK GUI, IntegrationReporterPython looks at the model at records any 'observed' variables (those with varname.obs_id set) into a data array in memory in Python. In the Tcl/Tk GUI, the 'observed' values are recorded to a text file that is later loaded into the GUI.

Results output

When integration has completed, the GUI receives an error message if any problem has occurred.

If all went well, the value of the data array is used (in the PyGTK GUI) to fill the contents of an Observer data table, which is then displayed to the user. This array will contain the timesteps that were requested, together with the values of any variables that had their 'obs_id' value set.

Data structures

This section outlines the key data structures used in communication between ASCEND and the integrator engines.

IntegratorInternals

This structure contains the bare-bones description of an integrator engine, including the hook functions that are used to 'drive' it, as well as a string and integer identifier for the engine.

IntegratorSystem

This structure contains the data common to all integrator engines. These include

  • pointer to the slv_system_t that is to be integrated (contains variable lists, etc).
  • pointer to the IntegratorInternals
  • pointer to an IntegratorReporter
  • pointer to a SampleList
  • pointer to internal data owned by the selected integrator engine
  • pointer to interface data owned by the GUI or other controlling process.
  • solver parameters data structure
  • independent variable
  • list of 'state' variables y.
  • list of 'derivative' variables ydot.
  • list of observed variables obs
  • other stuff

IntegratorReporter

This structure provides functions that will communicate integration results back to the GUI or other controlling process.

  • init, a function that will initialise the reporter, perhaps opening up a window or allocating some memory.
  • write, a function that should provide a quick status report, perhaps just telling the user the current value of the independent variable.
  • write_obs, a function that should examine the system and extract the value of any variables that need to be kept. This function is called each time the integrator reaches one of its designated timesteps as specified in the SampleList.
  • close, a function that can be used to finalise the reporting of results back to the user, deallocating memory, etc.

SampleList

This structure contains a list of timesteps (values of the independent variable, to be precise) for which output data is to be returned by the integrator. This list must be prepared ahead of time, before the integrator_solve function is called. In contains

  • ns, the number of sample values in the list
  • x, an array of double-precision values
  • d, a dim_type containing the DIMENSION of the provided sample values (for cross-checking with the dimensionality of the INDEPENDENT variable. The idea is that the GUI user can provide values of incorrect units (eg time = 10 cm) and we want to be able to catch that before it causes a problem).

See also