12.1. The SUNAdaptController API

Added in version 6.7.0.

The SUNAdaptController base class provides a common API for accuracy-based adaptivity controllers to be used by SUNDIALS integrators. These controllers estimate step sizes (among other things) such that the next step solution satisfies a desired temporal accuracy, while striving to maximize computational efficiency. We note that in the descriptions below, we frequently use dsm to represent temporal error. This is not the raw temporal error estimate; instead, it is a norm of the temporal error estimate after scaling by the user-supplied accuracy tolerances (see (2.19)),

\[\text{dsm} = \left( \frac{1}{N} \sum_{i=1}^N \left(\frac{\text{error}_i}{\text{rtol}\cdot |y_{n-1,i}| + \text{atol}_i}\right)^2\right)^{1/2}.\]

Thus dsm values below one represent errors estimated to be more accurate than needed, whereas errors above one are considered to be larger than allowable.

The base SUNAdaptController class is modeled after SUNDIALS’ other object-oriented classes, in that this class contains a pointer to an implementation-specific content, an ops structure with generic controller operations, and a SUNContext object. Specifically, the type SUNAdaptController is defined as:

typedef struct _generic_SUNAdaptController *SUNAdaptController

and the base class structure is defined as

struct _generic_SUNAdaptController {
     void* content;
     generic_SUNAdaptController_Ops* ops;
     SUNContext sunctx;
 };

Here, _generic_SUNAdaptController_Ops is the pointer to a structure containing function pointers to the various controller operations, and is defined as

struct _generic_SUNAdaptController_Ops {
    SUNAdaptController_Type (*getid)(SUNAdaptController C);
    SUNErrCode (*destroy)(SUNAdaptController C);
    SUNErrCode (*estimatestep)(SUNAdaptController C, sunrealtype h, int p, sunrealtype dsm, sunrealtype* hnew);
    SUNErrCode (*reset)(SUNAdaptController C);
    SUNErrCode (*setdefaults)(SUNAdaptController C);
    SUNErrCode (*write)(SUNAdaptController C, FILE* fptr);
    SUNErrCode (*seterrorbias)(SUNAdaptController C, sunrealtype bias);
    SUNErrCode (*updateh)(SUNAdaptController C, sunrealtype h, sunrealtype dsm);
    SUNErrCode (*space)(SUNAdaptController C, long int *lenrw, long int *leniw);
};

12.1.1. SUNAdaptController Types

The time integrators in SUNDIALS adapt a variety of parameters to achieve accurate and efficient computations. To this end, each SUNAdaptController implementation should note its type, so that integrators will understand the types of adaptivity that the controller is designed to perform. These are encoded in the following set of SUNAdaptController types:

enum SUNAdaptController_Type: The enumerated type SUNAdaptController_Type defines the enumeration constants for SUNDIALS error controller types

enumerator SUN_ADAPTCONTROLLER_NONE: Empty object that performs no control.

enumerator SUN_ADAPTCONTROLLER_H: Controls a single-rate step size.

12.1.2. SUNAdaptController Operations

The base SUNAdaptController class defines and implements all SUNAdaptController functions. Most of these routines are merely wrappers for the operations defined by a particular SUNAdaptController implementation, which are accessed through the ops field of the SUNAdaptController structure. The base SUNAdaptController class provides the constructor

SUNAdaptController SUNAdaptController_NewEmpty(SUNContext sunctx)

This function allocates a new generic SUNAdaptController object and initializes its content pointer and the function pointers in the operations structure to NULL.

Parameters:

sunctx – the SUNContext object (see §1.4)

Returns:

If successful, a generic SUNAdaptController object. If unsuccessful, a NULL pointer will be returned.

Each of the following methods are optional for any specific SUNAdaptController implementation, however some may be required based on the implementation’s SUNAdaptController_Type (see Section §12.1.1). We note these requirements below. Additionally, we note the behavior of the base SUNAdaptController methods when they perform an action other than only a successful return.

SUNAdaptController_Type SUNAdaptController_GetType(SUNAdaptController C)

Returns the type identifier for the controller C. Returned values are given in Section §12.1.1

Parameters:

C – the SUNAdaptController object.

Returns:

SUNAdaptController_Type type identifier.

Usage:

SUNAdaptController_Type id = SUNAdaptController_GetType(C);

SUNErrCode SUNAdaptController_Destroy(SUNAdaptController C)

Deallocates the controller C. If this method is not provided by the implementation, the base class method will free both the content and ops objects – this should be sufficient unless a controller implementation performs dynamic memory allocation of its own (note that the SUNDIALS-provided SUNAdaptController implementations do not need to supply this routine).

Parameters:

C – the SUNAdaptController object.

Returns:

SUNErrCode indicating success or failure.

Usage:

retval = SUNAdaptController_Destroy(C);

SUNErrCode SUNAdaptController_EstimateStep(SUNAdaptController C, sunrealtype h, int p, sunrealtype dsm, sunrealtype *hnew)

Estimates a single-rate step size. This routine is required for controllers of type SUN_ADAPTCONTROLLER_H. If this is not provided by the implementation, the base class method will set *hnew = h and return.

Parameters:

C – the SUNAdaptController object.
h – the step size from the previous step attempt.
p – the current order of accuracy for the time integration method.
dsm – the local temporal estimate from the previous step attempt.
hnew – (output) the estimated step size.

Returns:

SUNErrCode indicating success or failure.

Usage:

retval = SUNAdaptController_EstimateStep(C, hcur, p, dsm, &hnew);

SUNErrCode SUNAdaptController_Reset(SUNAdaptController C)

Resets the controller to its initial state, e.g., if it stores a small number of previous dsm or h values.

Parameters:

C – the SUNAdaptController object.

Returns:

SUNErrCode indicating success or failure.

Usage:

retval = SUNAdaptController_Reset(C);

SUNErrCode SUNAdaptController_SetDefaults(SUNAdaptController C)

Sets the controller parameters to their default values.

Parameters:

C – the SUNAdaptController object.

Returns:

SUNErrCode indicating success or failure.

Usage:

retval = SUNAdaptController_SetDefaults(C);

SUNErrCode SUNAdaptController_Write(SUNAdaptController C, FILE *fptr)

Writes all controller parameters to the indicated file pointer.

Parameters:

C – the SUNAdaptController object.
fptr – the output stream to write the parameters to.

Returns:

SUNErrCode indicating success or failure.

Usage:

retval = SUNAdaptController_Write(C, stdout);

SUNErrCode SUNAdaptController_SetErrorBias(SUNAdaptController C, sunrealtype bias)

Sets an error bias factor for scaling the local error factors. This is typically used to slightly exaggerate the temporal error during the estimation process, leading to a more conservative estimated step size.

Parameters:

C – the SUNAdaptController object.
bias – the error bias factor – an input \(\leq 0\) indicates to use the default value for the controller.

Returns:

SUNErrCode indicating success or failure.

Usage:

retval = SUNAdaptController_SetErrorBias(C, 1.2);

SUNErrCode SUNAdaptController_UpdateH(SUNAdaptController C, sunrealtype h, sunrealtype dsm)

Notifies a controller of type SUN_ADAPTCONTROLLER_H that a successful time step was taken with stepsize h and local error factor dsm, indicating that these can be saved for subsequent controller functions. This is typically relevant for controllers that store a history of either step sizes or error estimates for performing the estimation process.

Parameters:

C – the SUNAdaptController object.
h – the successful step size.
dsm – the successful temporal error estimate.

Returns:

SUNErrCode indicating success or failure.

Usage:

retval = SUNAdaptController_UpdateH(C, h, dsm);

SUNErrCode SUNAdaptController_Space(SUNAdaptController C, long int *lenrw, long int *leniw)

Informative routine that returns the memory requirements of the SUNAdaptController object.

Parameters:

C – the SUNAdaptController object..
lenrw – (output) number of sunsunrealtype words stored in the controller.
leniw – (output) number of sunindextype words stored in the controller. This may also include pointers, int and long int words.

Returns:

SUNErrCode indicating success or failure.

Usage:

retval = SUNAdaptController_Space(C, &lenrw, &leniw);

12.1.3. C/C++ API Usage

Specific SUNDIALS adaptivity controller modules can be used in C and C++ programs by including the corresponding header file for that module, e.g. sunadaptcontroller/sunadaptcontroller_XYZ.h.

Example usage (here SUNAdaptController_XYZ is a placeholder for an actual SUNAdaptController constructor):

#include <stdio.h>
#include <stdlib.h>
#include <sundials/sundials_context.h>
#include <sundials/sundials_types.h>
#include <sunadaptcontroller/sunadaptcontroller_XYZ.h>

int main()
{
    /* Create a SUNContext object */
    SUNContext sunctx = ...;

    /* Create a SUNAdaptController object */
    SUNAdaptController C = SUNAdaptController_XYZ(sunctx);

    /* Use the control object */

    /* Destroy the control object */
    retval = SUNAdaptController_Destroy(C);

    return 0;
}

12.2. The SUNAdaptController_Soderlind Module

The Soderlind implementation of the SUNAdaptController class, SUNAdaptController_Soderlind, implements a general structure for temporal control proposed by G. Soderlind in [113], [114], and [115]. This controller has the form

\[h' = h_n \varepsilon_n^{-k_1/(p+1)} \varepsilon_{n-1}^{-k_2/(p+1)} \varepsilon_{n-2}^{-k_3/(p+1)} \left(\dfrac{h_n}{h_{n-1}}\right)^{k_4} \left(\dfrac{h_{n-1}}{h_{n-2}}\right)^{k_5}\]

with default parameter values \(k_1 = 1.25\), \(k_2 = 0.5\), \(k_3 = -0.75\), \(k_4 = 0.25\), and \(k_5 = 0.75\), where \(p\) is the global order of the time integration method. In this estimate, a floor of \(\varepsilon_* > 10^{-10}\) is enforced to avoid division-by-zero errors. During the first two steps (when \(\varepsilon_{n-2}\), \(\varepsilon_{n-1}\), \(h_{n-2}\), and \(h_{n-2}\) may be unavailable), the corresponding terms are merely omitted during estimation of \(h'\).

The SUNAdaptController_Soderlind controller is implemented as a derived SUNAdaptController class, and defines its content field as:

struct _SUNAdaptControllerContent_Soderlind {
  sunrealtype k1;
  sunrealtype k2;
  sunrealtype k3;
  sunrealtype k4;
  sunrealtype k5;
  sunrealtype bias;
  sunrealtype ep;
  sunrealtype epp;
  sunrealtype hp;
  sunrealtype hpp;
  int firststeps;
};

These entries of the content field contain the following information:

k1, k2, k3, k4, k5 - controller parameters above.
bias - error bias factor, that converts from an input temporal error estimate via \(\varepsilon = \text{bias}*\text{dsm}\).
ep, epp - storage for the two previous error estimates, \(\varepsilon_{n-1}\) and \(\varepsilon_{n-2}\).
hp, hpp - storage for the previous two step sizes, \(h_{n-1}\) and \(h_{n-2}\).
firststeps - counter to handle first two steps (where previous step sizes and errors are unavailable).

The header file to be included when using this module is sunadaptcontroller/sunadaptcontroller_soderlind.h.

We note that through appropriate selection of the parameters \(k_*\), this controller may create a wide range of proposed temporal adaptivity controllers, including the PID, PI, I, as well as Gustafsson’s explicit and implicit controllers, [58] and [59]. As a convenience, utility routines to create these controllers and set their parameters (as special cases of the SUNAdaptController_Soderlind) are provided.

The SUNAdaptController_Soderlind class provides implementations of all operations relevant to a SUN_ADAPTCONTROLLER_H controller listed in §12.1.2. This class also provides the following additional user-callable routines:

SUNAdaptController SUNAdaptController_Soderlind(SUNContext sunctx)

This constructor creates and allocates memory for a SUNAdaptController_Soderlind object, and inserts its default parameters.

Parameters:

sunctx – the current SUNContext object.

Returns:

if successful, a usable SUNAdaptController object; otherwise it will return NULL.

Usage:

SUNAdaptController C = SUNAdaptController_Soderlind(sunctx);

SUNErrCode SUNAdaptController_SetParams_Soderlind(SUNAdaptController C, sunrealtype k1, sunrealtype k2, sunrealtype k3, sunrealtype k4, sunrealtype k5)

This user-callable function provides control over the relevant parameters above. This should be called before the time integrator is called to evolve the problem.

Parameters:

C – the SUNAdaptController_Soderlind object.
k1 – parameter used within the controller time step estimate.
k2 – parameter used within the controller time step estimate.
k3 – parameter used within the controller time step estimate.
k4 – parameter used within the controller time step estimate.
k5 – parameter used within the controller time step estimate.

Returns:

SUNErrCode indicating success or failure.

Usage:

/* Specify parameters for Soderlind's H_{0}312 controller */
retval = SUNAdaptController_SetParams_Soderlind(C, 0.25, 0.5, 0.25, -0.75, -0.25);

SUNAdaptController SUNAdaptController_PID(SUNContext sunctx)

This constructor creates and allocates memory for a SUNAdaptController_Soderlind object, set up to replicate a PID controller, and inserts its default parameters \(k_1=0.58\), \(k_2=-0.21\), \(k_3=0.1\), and \(k_4=k_5=0\).

Parameters:

sunctx – the current SUNContext object.

Returns:

if successful, a usable SUNAdaptController object; otherwise it will return NULL.

Usage:

SUNAdaptController C = SUNAdaptController_PID(sunctx);

SUNErrCode SUNAdaptController_SetParams_PID(SUNAdaptController C, sunrealtype k1, sunrealtype k2, sunrealtype k3)

This user-callable function provides control over the relevant parameters above for a PID controller, setting \(k_4 = k_5 = 0\). This should be called before the time integrator is called to evolve the problem.

Parameters:

C – the SUNAdaptController_Soderlind object.
k1 – parameter used within the controller time step estimate.
k2 – parameter used within the controller time step estimate.
k3 – parameter used within the controller time step estimate.

Returns:

SUNErrCode indicating success or failure.

Usage:

retval = SUNAdaptController_SetParams_PID(C, 0.58, -0.21, 0.1);

SUNAdaptController SUNAdaptController_PI(SUNContext sunctx)

This constructor creates and allocates memory for a SUNAdaptController_Soderlind object, set up to replicate a PI controller, and inserts its default parameters \(k_1=0.8\), \(k_2=-0.31\), and \(k_3=k_4=k_5=0\).

Parameters:

sunctx – the current SUNContext object.

Returns:

if successful, a usable SUNAdaptController object; otherwise it will return NULL.

Usage:

SUNAdaptController C = SUNAdaptController_PI(sunctx);

SUNErrCode SUNAdaptController_SetParams_PI(SUNAdaptController C, sunrealtype k1, sunrealtype k2)

This user-callable function provides control over the relevant parameters above for a PI controller, setting \(k_3 = k_4 = k_5 = 0\). This should be called before the time integrator is called to evolve the problem.

Parameters:

C – the SUNAdaptController_Soderlind object.
k1 – parameter used within the controller time step estimate.
k2 – parameter used within the controller time step estimate.

Returns:

SUNErrCode indicating success or failure.

Usage:

retval = SUNAdaptController_SetParams_PI(C, 0.8, -0.31);

SUNAdaptController SUNAdaptController_I(SUNContext sunctx)

This constructor creates and allocates memory for a SUNAdaptController_Soderlind object, set up to replicate an I controller, and inserts its default parameters \(k_1=1.0\) and \(k_2=k_3=k_4=k_5=0\).

Parameters:

sunctx – the current SUNContext object.

Returns:

if successful, a usable SUNAdaptController object; otherwise it will return NULL.

Usage:

SUNAdaptController C = SUNAdaptController_I(sunctx);

SUNErrCode SUNAdaptController_SetParams_I(SUNAdaptController C, sunrealtype k1)

This user-callable function provides control over the relevant parameters above for an I controller, setting \(k_2 = k_3 = k_4 = k_5 = 0\). This should be called before the time integrator is called to evolve the problem.

Parameters:

C – the SUNAdaptController_Soderlind object.
k1 – parameter used within the controller time step estimate.

Returns:

SUNErrCode indicating success or failure.

Usage:

retval = SUNAdaptController_SetParams_I(C, 1.0);

SUNAdaptController SUNAdaptController_ExpGus(SUNContext sunctx)

This constructor creates and allocates memory for a SUNAdaptController_Soderlind object, set up to replicate Gustafsson’s explicit controller [58], and inserts its default parameters \(k_1=0.635\), \(k_2=-0.268\), and \(k_3=k_4=k_5=0\).

Parameters:

sunctx – the current SUNContext object.

Returns:

if successful, a usable SUNAdaptController object; otherwise it will return NULL.

Usage:

SUNAdaptController C = SUNAdaptController_ExpGus(sunctx);

SUNErrCode SUNAdaptController_SetParams_ExpGus(SUNAdaptController C, sunrealtype k1_hat, sunrealtype k2_hat)

This user-callable function provides control over the relevant parameters above for the explicit Gustafsson controller, setting \(k_3 = k_4 = k_5 = 0\). This should be called before the time integrator is called to evolve the problem.

Note

Gustafsson’s explicit controller has the form

\[h' = h_n \varepsilon_n^{-\hat{k}_1/(p+1)} \left(\frac{\varepsilon_n}{\varepsilon_{n-1}}\right)^{-\hat{k}_2/(p+1)}.\]

The inputs to this function correspond to the values of \(\hat{k}_1\) and \(\hat{k}_2\), which are internally transformed into the Soderlind coeficients \(k_1 = \hat{k}_1+\hat{k}_2\) and \(k_2 = -\hat{k}_2\).

Parameters:

C – the SUNAdaptController_Soderlind object.
k1_hat – parameter used within the explicit Gustafsson controller time step estimate.
k2_hat – parameter used within the explicit Gustafsson controller time step estimate.

Returns:

SUNErrCode indicating success or failure.

Usage:

retval = SUNAdaptController_SetParams_ExpGus(C, 0.367, 0.268);

SUNAdaptController SUNAdaptController_ImpGus(SUNContext sunctx)

This constructor creates and allocates memory for a SUNAdaptController_Soderlind object, set up to replicate Gustafsson’s implicit controller [59], and inserts its default parameters \(k_1=1.93\), \(k_2=-0.95\), \(k_4=1\), and \(k_3=k_5=0\).

Parameters:

sunctx – the current SUNContext object.

Returns:

if successful, a usable SUNAdaptController object; otherwise it will return NULL.

Usage:

SUNAdaptController C = SUNAdaptController_ImpGus(sunctx);

SUNErrCode SUNAdaptController_SetParams_ImpGus(SUNAdaptController C, sunrealtype k1_hat, sunrealtype k2_hat)

This user-callable function provides control over the relevant parameters above for the implicit Gustafsson controller, setting \(k_3 = k_4 = k_5 = 0\). This should be called before the time integrator is called to evolve the problem.

Note

Gustafsson’s implicit controller has the form

\[h' = h_n \varepsilon_n^{-\hat{k}_1/(p+1)} \left(\frac{\varepsilon_n}{\varepsilon_{n-1}}\right)^{-\hat{k}_2/(p+1)} \left(\frac{h_n}{h_{n-1}}\right).\]

The inputs to this function correspond to the values of \(\hat{k}_1\) and \(\hat{k}_2\), which are internally transformed into the Soderlind coeficients \(k_1 = \hat{k}_1+\hat{k}_2\), \(k_2 = -\hat{k}_2\), and \(k_4=1\).

Parameters:

C – the SUNAdaptController_Soderlind object.
k1_hat – parameter used within the implicit Gustafsson controller time step estimate.
k2_hat – parameter used within the implicit Gustafsson controller time step estimate.

Returns:

SUNErrCode indicating success or failure.

Usage:

retval = SUNAdaptController_SetParams_ImpGus(C, 0.98, 0.95);

12.3. The SUNAdaptController_ImExGus Module

The ImEx Gustafsson implementation of the SUNAdaptController class, SUNAdaptController_ImExGus, implements a combination of two adaptivity controllers proposed by K. Gustafsson. His “explicit” controller was proposed in [58], is primarily useful with explicit Runge–Kutta methods, and has the form

(12.1)\[\begin{split}h' \;=\; \begin{cases} h_1\; \varepsilon_1^{-1/(p+1)}, &\quad\text{on the first step}, \\ h_n\; \varepsilon_n^{-k_1^E/(p+1)}\; \left(\dfrac{\varepsilon_n}{\varepsilon_{n-1}}\right)^{k_2^E/(p+1)}, & \quad\text{on subsequent steps}. \end{cases}\end{split}\]

Similarly, Gustafsson’s “implicit” controller was proposed in [59] with the form

(12.2)\[\begin{split}h' = \begin{cases} h_1 \varepsilon_1^{-1/(p+1)}, &\quad\text{on the first step}, \\ h_n \left(\dfrac{h_n}{h_{n-1}}\right) \varepsilon_n^{-k_1^I/(p+1)} \left(\dfrac{\varepsilon_n}{\varepsilon_{n-1}}\right)^{-k_2^I/(p+1)}, & \quad\text{on subsequent steps}. \end{cases}\end{split}\]

In the above formulas, the default values of \(k_1^E\), \(k_2^E\), \(k_1^I\), and \(k_2^I\) are 0.367, 0.268, 0.98, and 0.95, respectively, and \(p\) is the global order of the time integration method. In these estimates, a floor of \(\varepsilon_* > 10^{-10}\) is enforced to avoid division-by-zero errors.

The SUNAdaptController_ImExGus controller implements both formulas (12.1) and (12.2), and sets its recommended step size as the minimum of these two. It is implemented as a derived SUNAdaptController class, and defines its content field as:

struct _SUNAdaptControllerContent_ImExGus {
  sunrealtype k1e;
  sunrealtype k2e;
  sunrealtype k1i;
  sunrealtype k2i;
  sunrealtype bias;
  sunrealtype ep;
  sunrealtype hp;
  sunbooleantype firststep;
};

These entries of the content field contain the following information:

k1e, k2e - explicit controller parameters used in (12.1).
k1i, k2i - implicit controller parameters used in (12.2).
bias - error bias factor, that converts from an input temporal error estimate via \(\varepsilon = \text{bias}*\text{dsm}\).
ep - storage for the previous error estimate, \(\varepsilon_{n-1}\).
hp - storage for the previous step size, \(h_{n-1}\).
firststep - flag indicating whether a step has completed successfully, allowing the formulas above to transition between \(h_1\) and \(h_n\).

The header file to be included when using this module is sunadaptcontroller/sunadaptcontroller_imexgus.h.

The SUNAdaptController_ImExGus class provides implementations of all operations relevant to a SUN_ADAPTCONTROLLER_H controller listed in §12.1.2. The SUNAdaptController_ImExGus class also provides the following additional user-callable routines:

SUNAdaptController SUNAdaptController_ImExGus(SUNContext sunctx)

This constructor creates and allocates memory for a SUNAdaptController_ImExGus object, and inserts its default parameters.

Parameters:

sunctx – the current SUNContext object.

Returns:

if successful, a usable SUNAdaptController object; otherwise it will return NULL.

Usage:

SUNAdaptController C = SUNAdaptController_ImExGus(sunctx);

SUNErrCode SUNAdaptController_SetParams_ImExGus(SUNAdaptController C, sunrealtype k1e, sunrealtype k2e, sunrealtype k1i, sunrealtype k2i)

This user-callable function provides control over the relevant parameters above. This should be called before the time integrator is called to evolve the problem.

Parameters:

C – the SUNAdaptController_ImExGus object.
k1e – parameter used within the controller time step estimate.
k2e – parameter used within the controller time step estimate.
k1i – parameter used within the controller time step estimate.
k2i – parameter used within the controller time step estimate.

Returns:

SUNErrCode indicating success or failure.

Usage:

retval = SUNAdaptController_SetParams_ImExGus(C, 0.4, 0.3, -1.0, 1.0);