User-supplied functions

The user-supplied functions for MRIStep consist of:

Additionally, a user may supply a custom set of slow-to-fast coupling coefficients for the MRI method.

ODE right-hand side

The user must supply a function of type ARKRhsFn to specify the “slow” right-hand side of the ODE system:

typedef int (*ARKRhsFn)(realtype t, N_Vector y, N_Vector ydot, void* user_data)

This function computes a portion of the ODE right-hand side for a given value of the independent variable \(t\) and state vector \(y\).

Arguments:
  • t – the current value of the independent variable.
  • y – the current value of the dependent variable vector.
  • ydot – the output vector that forms a portion the ODE RHS \(f(t,y)\).
  • user_data – the user_data pointer that was passed to MRIStepSetUserData().

Return value: An ARKRhsFn should return 0 if successful, a positive value if a recoverable error occurred (in which case MRIStep will attempt to correct), or a negative value if it failed unrecoverably (in which case the integration is halted and ARK_RHSFUNC_FAIL is returned).

Notes: Allocation of memory for ydot is handled within the MRIStep module.

The vector ydot may be uninitialized on input; it is the user’s responsibility to fill this entire vector with meaningful values.

A recoverable failure error return from the ARKRhsFn is typically used to flag a value of the dependent variable \(y\) that is “illegal” in some way (e.g., negative where only a non-negative value is physically meaningful). If such a return is made within an implicit solve, MRIStep may attempt to recover by repeating the nonlinear iteration in order to avoid this recoverable error return. However, since MRIStep currently requires fixed time stepping at the slow time scale, no other recovery mechanisms are available, and MRIStep may halt on a recoverable error flag.

Error message handler function

As an alternative to the default behavior of directing error and warning messages to the file pointed to by errfp (see MRIStepSetErrFile()), the user may provide a function of type ARKErrHandlerFn to process any such messages.

typedef void (*ARKErrHandlerFn)(int error_code, const char* module, const char* function, char* msg, void* user_data)

This function processes error and warning messages from MRIStep and its sub-modules.

Arguments:
  • error_code – the error code.
  • module – the name of the MRIStep module reporting the error.
  • function – the name of the function in which the error occurred.
  • msg – the error message.
  • user_data – a pointer to user data, the same as the eh_data parameter that was passed to MRIStepSetErrHandlerFn().

Return value: An ARKErrHandlerFn function has no return value.

Notes: error_code is negative for errors and positive (ARK_WARNING) for warnings. If a function that returns a pointer to memory encounters an error, it sets error_code to 0.

Error weight function

As an alternative to providing the relative and absolute tolerances, the user may provide a function of type ARKEwtFn to compute a vector ewt containing the weights in the WRMS norm \(\|v\|_{WRMS} = \left(\frac{1}{n} \sum_{i=1}^n \left(ewt_i\; v_i\right)^2 \right)^{1/2}\). These weights will be used in place of those defined in the section Error norms.

typedef int (*ARKEwtFn)(N_Vector y, N_Vector ewt, void* user_data)

This function computes the WRMS error weights for the vector \(y\).

Arguments:
  • y – the dependent variable vector at which the weight vector is to be computed.
  • ewt – the output vector containing the error weights.
  • user_data – a pointer to user data, the same as the user_data parameter that was passed to MRIStepSetUserData().

Return value: An ARKEwtFn function must return 0 if it successfully set the error weights, and -1 otherwise.

Notes: Allocation of memory for ewt is handled within MRIStep.

The error weight vector must have all components positive. It is the user’s responsibility to perform this test and return -1 if it is not satisfied.

Implicit stage prediction function

A user may supply a function to update the prediction for each implicit stage solution. If supplied, this routine will be called after any existing MRIStep predictor algorithm completes, so that the predictor may be modified by the user as desired. In this scenario, a user may provide a function of type ARKStagePredictFn to provide this implicit predictor to MRIStep. This function takes as input the already-predicted implicit stage solution and the corresponding ‘time’ for that prediction; it then updates the prediction vector as desired. If the user-supplied routine will construct a full prediction (and thus the MRIStep prediction is irrelevant), it is recommended that the user not call MRIStepSetPredictorMethod(), thereby leaving the default trivial predictor in place.

typedef int (*ARKStagePredictFn)(realtype t, N_Vector zpred, void* user_data)

This function updates the prediction for the implicit stage solution.

Arguments:
  • t – the current value of the independent variable.
  • zpred – the MRIStep-predicted stage solution on input, and the user-modified predicted stage solution on output.
  • user_data – a pointer to user data, the same as the user_data parameter that was passed to MRIStepSetUserData().

Return value: An ARKStagePredictFn function should return 0 if it successfully set the upcoming stable step size, and a non-zero value otherwise.

Notes: This may be useful if there are bound constraints on the solution, and these should be enforced prior to beginning the nonlinear or linear implicit solver algorithm.

Rootfinding function

If a rootfinding problem is to be solved during the integration of the ODE system, the user must supply a function of type ARKRootFn.

typedef int (*ARKRootFn)(realtype t, N_Vector y, realtype* gout, void* user_data)

This function implements a vector-valued function \(g(t,y)\) such that the roots of the nrtfn components \(g_i(t,y)\) are sought.

Arguments:
  • t – the current value of the independent variable.
  • y – the current value of the dependent variable vector.
  • gout – the output array, of length nrtfn, with components \(g_i(t,y)\).
  • user_data – a pointer to user data, the same as the user_data parameter that was passed to MRIStepSetUserData().

Return value: An ARKRootFn function should return 0 if successful or a non-zero value if an error occurred (in which case the integration is halted and MRIStep returns ARK_RTFUNC_FAIL).

Notes: Allocation of memory for gout is handled within MRIStep.

Jacobian construction (matrix-based linear solvers)

If a matrix-based linear solver module is used (i.e., a non-NULL SUNMatrix object was supplied to MRIStepSetLinearSolver() in section A skeleton of the user’s main program), the user may provide a function of type ARKLsJacFn to provide the Jacobian approximation or ARKLsLinSysFn to provide an approximation of the linear system \(A = I - \gamma J\).

typedef int (*ARKLsJacFn)(realtype t, N_Vector y, N_Vector fy, SUNMatrix Jac, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)

This function computes the Jacobian matrix \(J = \frac{\partial f^S}{\partial y}\) (or an approximation to it).

Arguments:
  • t – the current value of the independent variable.
  • y – the current value of the dependent variable vector, namely the predicted value of \(y(t)\).
  • fy – the current value of the vector \(f^S(t,y)\).
  • Jac – the output Jacobian matrix.
  • user_data – a pointer to user data, the same as the user_data parameter that was passed to MRIStepSetUserData().
  • tmp1, tmp2, tmp3 – pointers to memory allocated to variables of type N_Vector which can be used by an ARKLsJacFn as temporary storage or work space.

Return value: An ARKLsJacFn function should return 0 if successful, a positive value if a recoverable error occurred (in which case MRIStep will attempt to correct, while ARKLS sets last_flag to ARKLS_JACFUNC_RECVR), or a negative value if it failed unrecoverably (in which case the integration is halted, MRIStepEvolve() returns ARK_LSETUP_FAIL and ARKLS sets last_flag to ARKLS_JACFUNC_UNRECVR).

Notes: Information regarding the structure of the specific SUNMatrix structure (e.g.~number of rows, upper/lower bandwidth, sparsity type) may be obtained through using the implementation-specific SUNMatrix interface functions (see the section Matrix Data Structures for details).

When using a linear solver of type SUNLINEARSOLVER_DIRECT, prior to calling the user-supplied Jacobian function, the Jacobian matrix \(J(t,y)\) is zeroed out, so only nonzero elements need to be loaded into Jac.

If the user’s ARKLsJacFn function uses difference quotient approximations, then it may need to access quantities not in the argument list. These include the current step size, the error weights, etc. To obtain these, the user will need to add a pointer to the ark_mem structure to their user_data, and then use the MRIStepGet* functions listed in Optional output functions. The unit roundoff can be accessed as UNIT_ROUNDOFF, which is defined in the header file sundials_types.h.

dense:

A user-supplied dense Jacobian function must load the N by N dense matrix Jac with an approximation to the Jacobian matrix \(J(t,y)\) at the point \((t,y)\). The accessor macros SM_ELEMENT_D and SM_COLUMN_D allow the user to read and write dense matrix elements without making explicit references to the underlying representation of the SUNMATRIX_DENSE type. SM_ELEMENT_D(J, i, j) references the (i,j)-th element of the dense matrix J (for i, j between 0 and N-1). This macro is meant for small problems for which efficiency of access is not a major concern. Thus, in terms of the indices \(m\) and \(n\) ranging from 1 to N, the Jacobian element \(J_{m,n}\) can be set using the statement SM_ELEMENT_D(J, m-1, n-1) = \(J_{m,n}\). Alternatively, SM_COLUMN_D(J, j) returns a pointer to the first element of the j-th column of J (for j ranging from 0 to N-1), and the elements of the j-th column can then be accessed using ordinary array indexing. Consequently, \(J_{m,n}\) can be loaded using the statements col_n = SM_COLUMN_D(J, n-1); col_n[m-1] = \(J_{m,n}\). For large problems, it is more efficient to use SM_COLUMN_D than to use SM_ELEMENT_D. Note that both of these macros number rows and columns starting from 0. The SUNMATRIX_DENSE type and accessor macros are documented in section The SUNMATRIX_DENSE Module.

band:

A user-supplied banded Jacobian function must load the band matrix Jac with the elements of the Jacobian \(J(t,y)\) at the point \((t,y)\). The accessor macros SM_ELEMENT_B, SM_COLUMN_B, and SM_COLUMN_ELEMENT_B allow the user to read and write band matrix elements without making specific references to the underlying representation of the SUNMATRIX_BAND type. SM_ELEMENT_B(J, i, j) references the (i,j)-th element of the band matrix J, counting from 0. This macro is meant for use in small problems for which efficiency of access is not a major concern. Thus, in terms of the indices \(m\) and \(n\) ranging from 1 to N with \((m, n)\) within the band defined by mupper and mlower, the Jacobian element \(J_{m,n}\) can be loaded using the statement SM_ELEMENT_B(J, m-1, n-1) \(= J_{m,n}\). The elements within the band are those with -mupper \(\le m-n \le\) mlower. Alternatively, SM_COLUMN_B(J, j) returns a pointer to the diagonal element of the j-th column of J, and if we assign this address to realtype *col_j, then the i-th element of the j-th column is given by SM_COLUMN_ELEMENT_B(col_j, i, j), counting from 0. Thus, for \((m,n)\) within the band, \(J_{m,n}\) can be loaded by setting col_n = SM_COLUMN_B(J, n-1); SM_COLUMN_ELEMENT_B(col_n, m-1, n-1) \(= J_{m,n}\) . The elements of the j-th column can also be accessed via ordinary array indexing, but this approach requires knowledge of the underlying storage for a band matrix of type SUNMATRIX_BAND. The array col_n can be indexed from -mupper to mlower. For large problems, it is more efficient to use SM_COLUMN_B and SM_COLUMN_ELEMENT_B than to use the SM_ELEMENT_B macro. As in the dense case, these macros all number rows and columns starting from 0. The SUNMATRIX_BAND type and accessor macros are documented in section The SUNMATRIX_BAND Module.

sparse:

A user-supplied sparse Jacobian function must load the compressed-sparse-column (CSC) or compressed-sparse-row (CSR) matrix Jac with an approximation to the Jacobian matrix \(J(t,y)\) at the point \((t,y)\). Storage for Jac already exists on entry to this function, although the user should ensure that sufficient space is allocated in Jac to hold the nonzero values to be set; if the existing space is insufficient the user may reallocate the data and index arrays as needed. The amount of allocated space in a SUNMATRIX_SPARSE object may be accessed using the macro SM_NNZ_S or the routine SUNSparseMatrix_NNZ(). The SUNMATRIX_SPARSE type is further documented in the section The SUNMATRIX_SPARSE Module.

typedef int (*ARKLsLinSysFn)(realtype t, N_Vector y, N_Vector fy, SUNMatrix A, SUNMatrix M, booleantype jok, booleantype *jcur, realtype gamma, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)

This function computes the linear system matrix \(A = I - \gamma J\) (or an approximation to it).

Arguments:
  • t – the current value of the independent variable.
  • y – the current value of the dependent variable vector, namely the predicted value of \(y(t)\).
  • fy – the current value of the vector \(f^S(t,y)\).
  • A – the output linear system matrix.
  • M – the argument will be NULL since MRIStep does not support non-identity mass matrices.
  • jok – is an input flag indicating whether the Jacobian-related data needs to be updated. The jok argument provides for the reuse of Jacobian data. When jok = SUNFALSE, the Jacobian-related data should be recomputed from scratch. When jok = SUNTRUE the Jacobian data, if saved from the previous call to this function, can be reused (with the current value of gamma). A call with jok = SUNTRUE can only occur after a call with jok = SUNFALSE.
  • jcur – is a pointer to a flag which should be set to SUNTRUE if Jacobian data was recomputed, or set to SUNFALSE if Jacobian data was not recomputed, but saved data was still reused.
  • gamma – the scalar \(\gamma\) appearing in the Newton matrix given by \(A=I-\gamma J\).
  • user_data – a pointer to user data, the same as the user_data parameter that was passed to MRIStepSetUserData().
  • tmp1, tmp2, tmp3 – pointers to memory allocated to variables of type N_Vector which can be used by an ARKLsLinSysFn as temporary storage or work space.

Return value: An ARKLsLinSysFn function should return 0 if successful, a positive value if a recoverable error occurred (in which case MRIStep will attempt to correct, while ARKLS sets last_flag to ARKLS_JACFUNC_RECVR), or a negative value if it failed unrecoverably (in which case the integration is halted, MRIStepEvolve() returns ARK_LSETUP_FAIL and ARKLS sets last_flag to ARKLS_JACFUNC_UNRECVR).

Jacobian-vector product (matrix-free linear solvers)

When using a matrix-free linear solver module for the implicit stage solves (i.e., a NULL-valued SUNMATRIX argument was supplied to MRIStepSetLinearSolver() in the section A skeleton of the user’s main program), the user may provide a function of type ARKLsJacTimesVecFn in the following form, to compute matrix-vector products \(Jv\). If such a function is not supplied, the default is a difference quotient approximation to these products.

typedef int (*ARKLsJacTimesVecFn)(N_Vector v, N_Vector Jv, realtype t, N_Vector y, N_Vector fy, void* user_data, N_Vector tmp)

This function computes the product \(Jv = \left(\frac{\partial f^S}{\partial y}\right)v\) (or an approximation to it).

Arguments:
  • v – the vector to multiply.
  • Jv – the output vector computed.
  • t – the current value of the independent variable.
  • y – the current value of the dependent variable vector.
  • fy – the current value of the vector \(f^S(t,y)\).
  • user_data – a pointer to user data, the same as the user_data parameter that was passed to MRIStepSetUserData().
  • tmp – pointer to memory allocated to a variable of type N_Vector which can be used as temporary storage or work space.

Return value: The value to be returned by the Jacobian-vector product function should be 0 if successful. Any other return value will result in an unrecoverable error of the generic Krylov solver, in which case the integration is halted.

Notes: If the user’s ARKLsJacTimesVecFn function uses difference quotient approximations, it may need to access quantities not in the argument list. These include the current step size, the error weights, etc. To obtain these, the user will need to add a pointer to the ark_mem structure to their user_data, and then use the MRIStepGet* functions listed in Optional output functions. The unit roundoff can be accessed as UNIT_ROUNDOFF, which is defined in the header file sundials_types.h.

Jacobian-vector product setup (matrix-free linear solvers)

If the user’s Jacobian-times-vector routine requires that any Jacobian-related data be preprocessed or evaluated, then this needs to be done in a user-supplied function of type ARKLsJacTimesSetupFn, defined as follows:

typedef int (*ARKLsJacTimesSetupFn)(realtype t, N_Vector y, N_Vector fy, void* user_data)

This function preprocesses and/or evaluates any Jacobian-related data needed by the Jacobian-times-vector routine.

Arguments:
  • t – the current value of the independent variable.
  • y – the current value of the dependent variable vector.
  • fy – the current value of the vector \(f^S(t,y)\).
  • user_data – a pointer to user data, the same as the user_data parameter that was passed to MRIStepSetUserData().

Return value: The value to be returned by the Jacobian-vector setup function should be 0 if successful, positive for a recoverable error (in which case the step will be retried), or negative for an unrecoverable error (in which case the integration is halted).

Notes: Each call to the Jacobian-vector setup function is preceded by a call to the implicit ARKRhsFn user function with the same \((t,y)\) arguments. Thus, the setup function can use any auxiliary data that is computed and saved during the evaluation of the implicit ODE right-hand side.

If the user’s ARKLsJacTimesSetupFn function uses difference quotient approximations, it may need to access quantities not in the argument list. These include the current step size, the error weights, etc. To obtain these, the user will need to add a pointer to the ark_mem structure to their user_data, and then use the MRIStepGet* functions listed in Optional output functions. The unit roundoff can be accessed as UNIT_ROUNDOFF, which is defined in the header file sundials_types.h.

Preconditioner solve (iterative linear solvers)

If a user-supplied preconditioner is to be used with a SUNLinSol solver module, then the user must provide a function of type ARKLsPrecSolveFn to solve the linear system \(Pz=r\), where \(P\) corresponds to either a left or right preconditioning matrix. Here \(P\) should approximate (at least crudely) the Newton matrix \(A=I-\gamma J\), where \(J = \frac{\partial f^S}{\partial y}\) If preconditioning is done on both sides, the product of the two preconditioner matrices should approximate \(A\).

typedef int (*ARKLsPrecSolveFn)(realtype t, N_Vector y, N_Vector fy, N_Vector r, N_Vector z, realtype gamma, realtype delta, int lr, void* user_data)

This function solves the preconditioner system \(Pz=r\).

Arguments:
  • t – the current value of the independent variable.
  • y – the current value of the dependent variable vector.
  • fy – the current value of the vector \(f^S(t,y)\).
  • r – the right-hand side vector of the linear system.
  • z – the computed output solution vector.
  • gamma – the scalar \(\gamma\) appearing in the Newton matrix given by \(A=I-\gamma J\).
  • delta – an input tolerance to be used if an iterative method is employed in the solution. In that case, the residual vector \(Res = r-Pz\) of the system should be made to be less than delta in the weighted \(l_2\) norm, i.e. \(\left(\sum_{i=1}^n \left(Res_i * ewt_i\right)^2 \right)^{1/2} < \delta\), where \(\delta =\) delta. To obtain the N_Vector ewt, call MRIStepGetErrWeights().
  • lr – an input flag indicating whether the preconditioner solve is to use the left preconditioner (lr = 1) or the right preconditioner (lr = 2).
  • user_data – a pointer to user data, the same as the user_data parameter that was passed to MRIStepSetUserData().

Return value: The value to be returned by the preconditioner solve function is a flag indicating whether it was successful. This value should be 0 if successful, positive for a recoverable error (in which case the step will be retried), or negative for an unrecoverable error (in which case the integration is halted).

Preconditioner setup (iterative linear solvers)

If the user’s preconditioner routine requires that any data be preprocessed or evaluated, then these actions need to occur within a user-supplied function of type ARKLsPrecSetupFn.

typedef int (*ARKLsPrecSetupFn)(realtype t, N_Vector y, N_Vector fy, booleantype jok, booleantype* jcurPtr, realtype gamma, void* user_data)

This function preprocesses and/or evaluates Jacobian-related data needed by the preconditioner.

Arguments:
  • t – the current value of the independent variable.
  • y – the current value of the dependent variable vector.
  • fy – the current value of the vector \(f^S(t,y)\).
  • jok – is an input flag indicating whether the Jacobian-related data needs to be updated. The jok argument provides for the reuse of Jacobian data in the preconditioner solve function. When jok = SUNFALSE, the Jacobian-related data should be recomputed from scratch. When jok = SUNTRUE the Jacobian data, if saved from the previous call to this function, can be reused (with the current value of gamma). A call with jok = SUNTRUE can only occur after a call with jok = SUNFALSE.
  • jcurPtr – is a pointer to a flag which should be set to SUNTRUE if Jacobian data was recomputed, or set to SUNFALSE if Jacobian data was not recomputed, but saved data was still reused.
  • gamma – the scalar \(\gamma\) appearing in the Newton matrix given by \(A=I-\gamma J\).
  • user_data – a pointer to user data, the same as the user_data parameter that was passed to MRIStepSetUserData().

Return value: The value to be returned by the preconditioner setup function is a flag indicating whether it was successful. This value should be 0 if successful, positive for a recoverable error (in which case the step will be retried), or negative for an unrecoverable error (in which case the integration is halted).

Notes: The operations performed by this function might include forming a crude approximate Jacobian, and performing an LU factorization of the resulting approximation to \(A = I - \gamma J\).

Each call to the preconditioner setup function is preceded by a call to the implicit ARKRhsFn user function with the same \((t,y)\) arguments. Thus, the preconditioner setup function can use any auxiliary data that is computed and saved during the evaluation of the ODE right-hand side.

This function is not called in advance of every call to the preconditioner solve function, but rather is called only as often as needed to achieve convergence in the Newton iteration.

If the user’s ARKLsPrecSetupFn function uses difference quotient approximations, it may need to access quantities not in the call list. These include the current step size, the error weights, etc. To obtain these, the user will need to add a pointer to the ark_mem structure to their user_data, and then use the MRIStepGet* functions listed in Optional output functions. The unit roundoff can be accessed as UNIT_ROUNDOFF, which is defined in the header file sundials_types.h.

Vector resize function

For simulations involving changes to the number of equations and unknowns in the ODE system (e.g. when using spatial adaptivity in a PDE simulation), the MRIStep integrator may be “resized” between integration steps, through calls to the MRIStepResize() function. Typically, when performing adaptive simulations the solution is stored in a customized user-supplied data structure, to enable adaptivity without repeated allocation/deallocation of memory. In these scenarios, it is recommended that the user supply a customized vector kernel to interface between SUNDIALS and their problem-specific data structure. If this vector kernel includes a function of type ARKVecResizeFn to resize a given vector implementation, then this function may be supplied to MRIStepResize() so that all internal MRIStep vectors may be resized, instead of deleting and re-creating them at each call. This resize function should have the following form:

typedef int (*ARKVecResizeFn)(N_Vector y, N_Vector ytemplate, void* user_data)

This function resizes the vector y to match the dimensions of the supplied vector, ytemplate.

Arguments:
  • y – the vector to resize.
  • ytemplate – a vector of the desired size.
  • user_data – a pointer to user data, the same as the resize_data parameter that was passed to MRIStepResize().

Return value: An ARKVecResizeFn function should return 0 if it successfully resizes the vector y, and a non-zero value otherwise.

Notes: If this function is not supplied, then MRIStep will instead destroy the vector y and clone a new vector y off of ytemplate.

Pre inner integrator communication function

The user may supply a function of type MRIStepPreInnerFn that will be called before each inner integration to perform any communication or memory transfers of forcing data supplied by the outer integrator to the inner integrator for the inner integration.

typedef int (*MRIStepPreInnerFn)(realtype t, N_Vector* f, int num_vecs, void* user_data)
Arguments:
  • t – the current value of the independent variable.
  • f – an N_Vector array of outer forcing vectors.
  • num_vecs – the number of vectors in the N_Vector array.
  • user_data – the user_data pointer that was passed to MRIStepSetUserData().

Return value: An MRIStepPreInnerFn function should return 0 if successful, a positive value if a recoverable error occurred, or a negative value if an unrecoverable error occurred. As the MRIStep module only supports fixed step sizes at this time any non-zero return value will halt the integration.

Notes: In a heterogeneous computing environment if any data copies between the host and device vector data are necessary, this is where that should occur.

Post inner integrator communication function

The user may supply a function of type MRIStepPostInnerFn that will be called after each inner integration to perform any communication or memory transfers of state data supplied by the inner integrator to the outer integrator for the outer integration.

typedef int (*MRIStepPostInnerFn)(realtype t, N_Vector y, void* user_data)
Arguments:
  • t – the current value of the independent variable.
  • y – the current value of the dependent variable vector.
  • user_data – the user_data pointer that was passed to MRIStepSetUserData().

Return value: An MRIStepPostInnerFn function should return 0 if successful, a positive value if a recoverable error occurred, or a negative value if an unrecoverable error occurred. As the MRIStep module only supports fixed step sizes at this time any non-zero return value will halt the integration.

Notes: In a heterogeneous computing environment if any data copies between the host and device vector data are necessary, this is where that should occur.

tocdepth:3

MRI Coupling Coefficients Data Structure

As mentioned in the section MRIStep User-callable functions, a user may supply a custom set of coupling coefficients between fast and slow time scales (through MRIStepSetCoupling()). MRIStep uses a custom data type, MRIStepCoupling, to store these coefficients, and provides several related utility routines. As described in the Section MRIStep – Multirate infinitesimal step methods, the coupling from slow to fast time scales in MRI methods may be encoded by a vector of slow ‘stage time’ abcissae, \(c^S \in \mathbb{R}^{s+1}\) and a set of coupling matrices \(\Gamma^{\{k\}}\in\mathbb{R}^{(s+1)\times(s+1)}\). We therefore define the MRIStepCoupling structure to be

typedef MRIStepCouplingMem *MRIStepCoupling

where MRIStepCouplingMem is the structure

struct MRIStepCouplingMem {

  int nmat;
  int stages;
  int q;
  int p;
  realtype ***G;
  realtype *c;

};

Here,

  • nmat corresponds to the number of \(\Gamma^{\{k\}}\) matrices used for coupling,
  • stages is the number of entries in c, i.e., \(s+1\) above,
  • q and p indicate the orders of accuracy for both the MRI method and the embedding, respectively,
  • G is a three-dimensional array with dimensions [nmat][stages][stages] containing the set of \(\Gamma^{\{k\}}\) matrices, and
  • c is an array of length stages containing slow abcissae \(c^S\) for the MRI method.

MRIStepCoupling functions

Function name Description
MRIStepCoupling_LoadTable() Loads a pre-defined MRIStepCoupling table
MRIStepCoupling_Alloc() Allocate an empty MRIStepCoupling table
MRIStepCoupling_Create() Create a new MRIStepCoupling table from coefficients
MRIStepCoupling_MIStoMRI() Create a new MRIStepCoupling table from a slow Butcher table
MRIStepCoupling_Copy() Create a copy of a MRIStepCoupling table
MRIStepCoupling_Space() Get the MRIStepCoupling table real and integer workspace sizes
MRIStepCoupling_Free() Deallocate a MRIStepCoupling table
MRIStepCoupling_Write() Write the MRIStepCoupling table to an output file
MRIStepCoupling MRIStepCoupling_LoadTable(int imethod)

Retrieves a specified MRIStepCoupling table. The prototype for this function, as well as the integer names for each provided method, are defined in top of the header file arkode/arkode_mristep.h. For further information on the current set of coupling tables and their corresponding identifiers, see MRIStepCoupling tables.

Arguments:
  • itable – MRIStepCoupling table identifier to load.
Return value:
  • MRIStepCoupling structure if successful.
  • NULL pointer if itable was invalid or an allocation error occured.
MRIStepCoupling MRIStepCoupling_Alloc(int nmat, int stages)

Allocates an empty MRIStepCoupling table.

Arguments:
  • nmat – number of \(\Gamma^{\{k\}}\) matrices in the coupling table.
  • stages – number of stages in the coupling table.
Return value:
  • MRIStepCoupling structure if successful.
  • NULL pointer if stages was invalid or an allocation error occured.
MRIStepCoupling MRIStepCoupling_Create(int nmat, int stages, int q, int p, realtype *G, realtype *c)

Allocates an MRIStepCoupling table and fills it with the given values.

Arguments:
  • nmat – number of \(\Gamma^{\{k\}}\) matrices in the coupling table.
  • stages – number of stages in the MRI method.
  • q – global order of accuracy for the MRI method.
  • p – global order of accuracy for the embedded MRI method.
  • G – array of coefficients defining the \(\Gamma^{\{k\}}\) matrices. This should be stored as a 1D array of size nmat*stages*stages, in row-major order.
  • c – array (of length stages) of slow abcissae for the MRI method.
Return value:
  • MRIStepCoupling structure if successful.
  • NULL pointer if stages was invalid or an allocation error occured.

Notes: As embeddings are not currently supported in MRIStep, then p should be equal to zero.

MRIStepCoupling MRIStepCoupling_MIStoMRI(ARKodeButcherTable B, int q, int p)

Creates an MRI coupling table for a traditional MIS method based on the slow Butcher table B, following the formula shown in (11).

Arguments:
  • B – the ARKodeButcherTable for the ‘slow’ MIS method.
  • q – the overall order of the MIS/MRI method.
  • p – the overall order of the MIS/MRI embedding.
Return value:
  • MRIStepCoupling structure if successful.
  • NULL pointer an allocation error occured.

Notes: The \(s\)-stage slow Butcher table must have an explicit first stage (i.e., \(c_1=0\) and \(A_{1,j}=0\) for \(1\le j\le s\)) and sorted abcissae (i.e., \(c_{i} \ge c_{i-1}\) for \(2\le i\le s\)).

Since an MIS method is at most third order accurate, and even then only if it meets certain compatibility criteria (see (12)), the values of q and p may differ from the method and embedding orders of accuracy for the Runge–Kutta encoded in B, which is why these arguments should be supplied separately.

As embeddings are not currently supported in MRIStep, then p should be equal to zero.

MRIStepCoupling MRIStepCoupling_Copy(MRIStepCoupling C)

Creates copy of the given coupling table.

Arguments:
  • C – the coupling table to copy.
Return value:
  • MRIStepCoupling structure if successful.
  • NULL pointer an allocation error occured.
void MRIStepCoupling_Space(MRIStepCoupling C, sunindextype *liw, sunindextype *lrw)

Get the real and integer workspace size for an MRIStepCoupling table.

Arguments:
  • C – the coupling table.
  • lenrw – the number of realtype values in the coupling table workspace.
  • leniw – the number of integer values in the coupling table workspace.
Return value:
  • ARK_SUCCESS if successful.
  • ARK_MEM_NULL if the Butcher table memory was NULL.
void MRIStepCoupling_Free(MRIStepCoupling C)

Deallocate the coupling table memory.

Arguments:
  • C – the MRIStepCoupling table.
void MRIStepCoupling_Write(MRIStepCoupling C, FILE *outfile)

Write the coupling table to the provided file pointer.

Arguments:
  • C – the MRIStepCoupling table.
  • outfile – pointer to use for printing the table.

Notes: The outfile argument can be stdout or stderr, or it may point to a specific file created using fopen.

MRIStepCoupling tables

MRIStep currently includes two classes of MRI coupling tables: those that encode methods that are explicit at the slow time scale, and those that are diagonally-implicit and solve-decoupled at the slow time scale. We list the current identifiers, multirate order of accuracy, and relevant references for each in the tables below. For the implicit methods, we also list the number of implicit solves per step that are required at the slow time scale.

Explicit MRI coupling tables:

Table name Order Reference
MIS_KW3 3 [SKAW2009]
MRI_GARK_ERK45a 4 [S2019]

Diagonally-implicit, solve-decoupled MRI coupling tables:

Table name Order Implicit Solves Reference
MRI_GARK_IRK21a 2 1 [S2019]
MRI_GARK_ESDIRK34a 4 3 [S2019]