22.1. Using sundials4py

At a high level, using SUNDIALS from Python via sundials4py looks a lot like using SUNDIALS from C or C++. Below we overview using sundials4py and discuss the few notable differences.

22.1.1. Installation

You can install sundials4py directly from PyPI using pip:

pip install sundials4py

You can also install sundials4py from git:

pip install git+https://github.com/LLNL/sundials.git

The default build of sundials4py that is distributed as a binary wheel uses double precision real types and 64-bit indices. To install SUNDIALS with different precisions and index sizes, you can build from source wheels instead of using the pre-built binary wheels. When building from source wheels instead of binary wheels, you can customize the SUNDIALS precision (real type) and index type at build time by passing the CMake arguments in an environment variable when running pip. For example:

export CMAKE_ARGS="-DSUNDIALS_PRECISION=SINGLE -DSUNDIALS_INDEX_SIZE=64"
pip install sundials4py --no-binary=sundials4py

Other SUNDIALS options can also be accessed in this way. Review §1.1.3 for more information on the available options.

22.1.2. Modules

After installation, you can import the sundials4py module with

import sundials4py

which includes the following submodules (which may also be individually imported) for accessing specific SUNDIALS features:

  • sundials4py.core contains all the shared SUNDIALS classes and functions as well as many of the native SUNDIALS class implementations:

    • NVector: serial and many-vector

    • SUNMatix: band, dense, and sparse

    • SUNLinearSover: band, dense, PCG, SPBCGS, SPFGMR, SPGMR, and SPTFQMR

    • SUNNonlinearSolver: fixed-point and Newton

    • SUNAdaptController: Soderlind, ImEx-Gus, and MRI H-Tol

    • SUNDomEigEstimator: Power

    • SUNAdjointCheckPointScheme: Fixed

  • sundials4py.arkode contains all of the ARKODE specific classes and functions

  • sundials4py.cvodes contains all of the CVODES specific classes and functions

  • sundials4py.idas contains all of the IDAS specific classes and functions

  • sundials4py.kinsol contains all of the KINSOL specific classes and functions

CVODE and IDA dot not have modules because CVODES and IDAS provide all of the same capabilities plus continuous forward and adjoint sensitivity analysis.

Note

Not all SUNDIALS features are supported by the Python interfaces. In particular, third-party libraries are not yet supported.

22.1.3. Example Usage

We now consider a simple CVODE example to illustrate using sundials4py and highlight some of the differences to using SUNDIALS from C/C++. The items highlighted below similarly apply to using other SUNDIALS packages. For more information on usage differences, continue to the next section. Additional examples can be found in the examples/python directory of the SUNDIALS GitHub repository.

This example demonstrates how to use CVODES to solve the Lotka-Volterra equations, a model of predator-prey dynamics in ecology, given by

\[\begin{split}u' &= p_0 u - p_1 u v \\ v' &= -p_2 v + p_3 u v\end{split}\]

where \(u\) is the prey population, \(v\) is the predator population, \(p_0\) is prey birth rate, \(p_1\) is the predation rate, \(p_2\) is the predator death rate, and \(p_3\) is predator growth rate from predation. We use the parameters \(p = [1.5, 1.0, 3.0, 1.0]\), initial condition \(y(0) = [1.0, 1.0]\), and integration interval \(t \in [0, 10]\).

  1import numpy as np
  2import sys
  3import matplotlib.pyplot as plt
  4from sundials4py.core import *  # Always import the core submodule
  5from sundials4py.cvodes import *  # Import the desired SUNDIALS package
  6
  7
  8class LotkaVolterraODE:
  9    """
 10    Encapsulates the Lotka-Volterra ODE problem.
 11
 12    This class defines the ODE system and provides the functions that CVODE needs to
 13    evolve it in time: the right-hand side (RHS) function and the Jacobian.
 14    """
 15
 16    def __init__(self, p):
 17        """
 18        Initialize with model parameters.
 19
 20        Args:
 21            p: list or array of 4 parameters [p_0, p_1, p_2, p_3]
 22        """
 23        self.p = np.array(p, dtype=sunrealtype)
 24        self.NEQ = 2  # Number of equations in the system
 25
 26    def set_init_cond(self, yvec):
 27        """
 28        Set the initial conditions in the solution vector.
 29
 30        Args:
 31            yvec: N_Vector to store initial values
 32
 33        Returns:
 34            0 on success
 35        """
 36        y = N_VGetArrayPointer(yvec)  # Returns a numpy ndarray view of the data
 37        y[0] = 1.0  # Initial prey population
 38        y[1] = 1.0  # Initial predator population
 39        return 0
 40
 41    def rhs(self, t, yvec, ydotvec, _):
 42        """
 43        Right-hand side function: computes dy/dt = f(t, y).
 44
 45        Args:
 46            t: current time (not used in this autonomous system)
 47            yvec: N_Vector with the current solution values y = [u, v]
 48            ydotvec: N_Vector output vector for time derivatives f = [du/dt, dv/dt]
 49            _: user data pointer MUST NOT be used in Python
 50
 51        Returns:
 52            0 on success
 53            positive value for a recoverable error (integrator will retry the step)
 54            negative value for an unrecoverable error (integration will halt)
 55        """
 56        p = self.p
 57        y = N_VGetArrayPointer(yvec)  # Returns a numpy ndarray view of the data
 58        ydot = N_VGetArrayPointer(ydotvec)
 59
 60        # Compute the derivatives
 61        ydot[0] = p[0] * y[0] - p[1] * y[0] * y[1]  # du/dt: prey dynamics
 62        ydot[1] = -p[2] * y[1] + p[3] * y[0] * y[1]  # dv/dt: predator dynamics
 63        return 0
 64
 65    def jac(self, t, yvec, fyvec, J, _, tmp1, tmp2, tmp3):
 66        """
 67        Jacobian function: computes J = df/dy.
 68
 69        Args:
 70            t: current time
 71            yvec: N_Vector with the current solution values
 72            fyvec: N_Vector with the current RHS values (not used here)
 73            J: SUNmatrix to store the output Jacobian matrix
 74            _: user data pointer MUST NOT be used in Python
 75            tmp1, tmp2, tmp3: temporary workspace N_Vectors (not used here)
 76
 77        Returns:
 78            0 on success
 79            positive value for a recoverable error (integrator will retry the step)
 80            negative value for an unrecoverable error (integration will halt)
 81        """
 82        p = self.p
 83        y = N_VGetArrayPointer(yvec)  # Returns a numpy ndarray view of the data
 84        Jdata = SUNDenseMatrix_Data(J)  # Returns a numpy ndarray view of the data
 85
 86        # Compute partial derivatives
 87        # J[0,0] = d(du/dt)/du, J[0,1] = d(du/dt)/dv
 88        Jdata[0, 0] = p[0] - p[1] * y[1]
 89        Jdata[0, 1] = -p[1] * y[0]
 90
 91        # J[1,0] = d(dv/dt)/du, J[1,1] = d(dv/dt)/dv
 92        Jdata[1, 0] = p[3] * y[1]
 93        Jdata[1, 1] = -p[2] + p[3] * y[0]
 94        return 0
 95
 96
 97def main():
 98    """
 99    Main function demonstrating the SUNDIALS/CVODE workflow.
100
101    A typical workflow for solving an ODE with CVODE is:
102    1. Create the SUNDIALS context
103    2. Create the solution vector and set initial conditions
104    3. Create and initialize CVODE
105    4. Configure CVODE (set tolerances, linear solver, Jacobian, etc.)
106    5. Advance the solution in time
107    6. Retrieve CVODE statistics
108    7. Destroy objects and free memory (handled automatically in Python)
109    """
110
111    # ----------------------------------------------------------------------------------
112    # Step 1: Create the SUNDIALS context
113    # ----------------------------------------------------------------------------------
114
115    # The context provides support for features such as error handling, logging, and
116    # profiling and is required to construct all SUNDIALS objects. SUN_COMM_NULL is used
117    # for serial (non-parallel) problems.
118
119    # In C, sunctx is an output parameter (return-by-pointer):
120    # SUNContext_Create(SUN_COMM_NULL, &sunctx). In Python, it is returned in a tuple
121    # where the first entry is the function return value and the second is the context.
122    status, sunctx = SUNContext_Create(SUN_COMM_NULL)
123    assert status == SUN_SUCCESS
124
125    # ----------------------------------------------------------------------------------
126    # Step 2: Set up the ODE problem
127    # ----------------------------------------------------------------------------------
128
129    # Create the ODE problem with parameters p = [1.5, 1.0, 3.0, 1.0]
130    params = [1.5, 1.0, 3.0, 1.0]
131    ode = LotkaVolterraODE(params)
132
133    # Create a serial N_Vector to hold the solution
134    y = N_VNew_Serial(ode.NEQ, sunctx)
135    assert y is not None
136
137    # Set initial conditions: y(0) = [1.0, 1.0]
138    ode.set_init_cond(y)
139
140    # ----------------------------------------------------------------------------------
141    # Step 3: Create and initialize CVODE
142    # ----------------------------------------------------------------------------------
143
144    # Create a CVODE instance using Backward Differentiation Formulas (BDF)
145    cvode = CVodeCreate(CV_BDF, sunctx)
146    assert cvode is not None
147
148    # Initialize CVODE with the RHS function, initial time (t0=0.0), and initial state
149    status = CVodeInit(cvode.get(), ode.rhs, 0.0, y)  # access CVODE memory with .get()
150    assert status == CV_SUCCESS
151
152    # ----------------------------------------------------------------------------------
153    # Step 4: Set CVODE options (tolerances, linear solver, etc.)
154    # ----------------------------------------------------------------------------------
155
156    # CVODE will adapt the step size and method order so the estimated local error meets
157    # the user defined tolerances.
158    reltol = 1e-6  # Relative tolerance
159    abstol = 1e-10  # Absolute tolerance
160    status = CVodeSStolerances(cvode.get(), reltol, abstol)
161    assert status == CV_SUCCESS
162
163    # CVODE defaults to using a Newton iteration to solve nonlinear systems at each time
164    # step which requires solving a linear system each iteration. For small problems, a
165    # dense direct solver is efficient and easy to use.
166
167    # Create a dense matrix (NEQ x NEQ) to store the Jacobian
168    A = SUNDenseMatrix(ode.NEQ, ode.NEQ, sunctx)
169    assert A is not None
170
171    # Create a dense linear solver that works with the matrix A and vector y
172    LS = SUNLinSol_Dense(y, A, sunctx)
173    assert LS is not None
174
175    # Attach the linear solver to CVODE
176    status = CVodeSetLinearSolver(cvode.get(), LS, A)
177    assert status == CV_SUCCESS
178
179    # Providing an analytical Jacobian can significantly improve performance. If not
180    # provided, CVODE will approximate it using finite differences.
181    status = CVodeSetJacFn(cvode.get(), ode.jac)
182    assert status == CV_SUCCESS
183
184    # ----------------------------------------------------------------------------------
185    # Step 5: Advance the ODE in time
186    # ----------------------------------------------------------------------------------
187
188    # Set the current output time and get access to the underlying array data for output
189    tret = 0.0
190    yarr = N_VGetArrayPointer(y)  # Returns a numpy ndarray view of the data
191
192    # Print header with problem information
193    print("\nLotka-Volterra ODE (CVODE):")
194    print(f"    initial condition, y0 = [1.0, 1.0]")
195    print(f"    parameters = {params}")
196    print(f"    reltol = {reltol}, abstol = {abstol}\n")
197    print("        t           u           v")
198    print("   ---------------------------------")
199    print(f"  {tret:10.6f}  {yarr[0]:10.6f}  {yarr[1]:10.6f}")
200
201    # Lists to store solution for plotting
202    t_vals = [tret]
203    u_vals = [yarr[0]]
204    v_vals = [yarr[1]]
205
206    # Integrate from t=0 to t=10, outputting every 0.05 time units using CV_NORMAL mode.
207    # In normal mode, CVODE will take internal steps until it has reached or just passed
208    # the output time and then return a time interpolated solution at the output time.
209    dtout = 0.05
210    iout = 0
211    while tret < 10.0:
212        # Advance the system and return the solution at requested output time
213        status, tret = CVode(cvode.get(), tret + dtout, y, CV_NORMAL)
214        assert status == CV_SUCCESS
215
216        # Store solution for plotting (every output)
217        t_vals.append(tret)
218        u_vals.append(yarr[0])
219        v_vals.append(yarr[1])
220
221        # Print every 10th output (every 1.0 time unit) to keep the output readable
222        iout += 1
223        if iout % 10 == 0 or tret >= 10.0:
224            print(f"  {tret:10.6f}  {yarr[0]:10.6f}  {yarr[1]:10.6f}")
225    print("   ---------------------------------")
226
227    # ----------------------------------------------------------------------------------
228    # Step 6: Get CVODE statistics
229    # ----------------------------------------------------------------------------------
230
231    # CVODE provides various statistics about the integration that can help assess
232    # performance and diagnose issues. These include values such as the step counts,
233    # function evaluations, and information about the linear solver.
234
235    status, nst = CVodeGetNumSteps(cvode.get())
236    assert status == CV_SUCCESS
237
238    status, netf = CVodeGetNumErrTestFails(cvode.get())
239    assert status == CV_SUCCESS
240
241    status, nfe = CVodeGetNumRhsEvals(cvode.get())
242    assert status == CV_SUCCESS
243
244    status, nsetups = CVodeGetNumLinSolvSetups(cvode.get())
245    assert status == CV_SUCCESS
246
247    status, nje = CVodeGetNumJacEvals(cvode.get())
248    assert status == CV_SUCCESS
249
250    # For C functions with multiple output parameters (return-by-pointer),
251    # CVodeGetNonlinSolvStats(cvode_mem, &nni, &ncfn), the first element is always the
252    # function return value (error code) and the subsequent elements follow the same
253    # ordering as in the C function signature.
254    status, nni, ncfn = CVodeGetNonlinSolvStats(cvode.get())
255    assert status == CV_SUCCESS
256
257    print("\nIntegrator Statistics:")
258    print(f"  Number of steps taken                    = {nst}")
259    print(f"  Number of error test fails               = {netf}")
260    print(f"  Number of RHS evaluations                = {nfe}")
261    print(f"  Number of linear solver setups           = {nsetups}")
262    print(f"  Number of Jacobian evaluations           = {nje}")
263    print(f"  Number of nonlinear solver iterations    = {nni}")
264    print(f"  Number of nonlinear convergence failures = {ncfn}")
265
266    # ----------------------------------------------------------------------------------
267    # Plot the solution
268    # ----------------------------------------------------------------------------------
269
270    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
271
272    # Left plot: Time series of prey and predator populations
273    ax1.plot(t_vals, u_vals, "b-", label="Prey (u)")
274    ax1.plot(t_vals, v_vals, "r--", label="Predator (v)")
275    ax1.set_xlabel("Time")
276    ax1.set_ylabel("Population")
277    ax1.set_title("Lotka-Volterra: Population vs Time")
278    ax1.legend()
279    ax1.grid(alpha=0.3)
280
281    # Right plot: Phase portrait (predator vs prey)
282    ax2.plot(u_vals, v_vals, "g-")
283    ax2.plot(u_vals[0], v_vals[0], "ko", label="Start")
284    ax2.plot(u_vals[-1], v_vals[-1], "ks", label="End")
285    ax2.set_xlabel("Prey (u)")
286    ax2.set_ylabel("Predator (v)")
287    ax2.set_title("Phase Portrait")
288    ax2.legend()
289    ax2.grid(alpha=0.3)
290
291    # Display the plot if not running in test mode
292    if "pytest" not in sys.modules:
293        plt.tight_layout()
294        plt.show()
295    else:
296        plt.close("all")
297
298
299if __name__ == "__main__":
300    main()

22.1.4. Usage Differences

While sundials4py closely follows the C API, some differences are inevitable due to the differences between Python and C as well as the requirements of the code generation tool used to create the bindings. In this section, we note the most critical differences.

22.1.4.1. View Classes and Memory Management

sundials4py provides natural usage of SUNDIALS objects with object lifetimes managed by the Python garbage collection as with any other Python object. There is only one caveat, the SUNDIALS integrator/solver void* objects (those returned by ARKODE, CVODES, IDAS, and KINSOL Create constructors) are wrapped in “View” classes (behind the scenes) for compatibility with nanobind. These view objects cannot be implicitly converted to the underlying void*. As such, when calling a function which operates on these void* objects, one must extract the void* “capsule” from the view object by calling the view’s get method.

# Create CVODE object (returns void* in C)
cvode = CVodeCreate(CV_BDF, sunctx)

# Notice we need to call cvode.get()
status = CVodeInit(cvode.get(), ode_problem.f, T0, y)

22.1.4.2. Return-by-Pointer Parameters

Functions that return values via pointer arguments in the C API are mapped to Python functions that return a tuple where the first element is the function’s return value (typically an error code) and subsequent elements are the values that would be returned via pointer arguments in C, in the same order as the C function signature.

Example 1: Single Return-by-Pointer Value

C:
int retval;
long int numsteps;
retval = CVodeGetNumSteps(cvode_mem, &numsteps);
printf("Number of steps: %ld\n", numsteps);
Python:
retval, numsteps = CVodeGetNumSteps(cvode_mem.get())
print(f"Number of steps: {numsteps}")

Example 2: Multiple Return-by-Pointer Values

C:
int retval;
long int nni, ncfn;
retval = CVodeGetNonlinSolvStats(cvode_mem, &nni, &ncfn)
printf("Nonlinear iterations: %ld, Nonlinear convergence fails: %ld\n", nni, ncfn);
Python:
retval, nni, ncfn = CVodeGetNonlinSolvStats(cvode_mem.get())
print(f"Nonlinear iterations: {nni}, Nonlinear convergence fails: {ncnf}");

22.1.4.3. Arrays

N_Vector objects in sundials4py are compatible with numpy’s ndarray. Each N_Vector can work on a numpy arrays without copies, and you can access and modify the underlying data directly using N_VGetArrayPointer, which returns a numpy ndarray view of the data.

SUNDIALS matrix types (dense, banded, sparse) are also exposed as Python objects that provide access to their underlying data as numpy arrays (e.g., via SUNDenseMatrix_Data).

Arrays of scalars (e.g., scaling factors passed to N_VLinearCombination) are also represented as numpy arrays.

Example: Accessing and modifying an N_Vector

y_nvec = N_VNew_Serial(10, sunctx)
y = N_VGetArrayPointer(y_nvec)
y[:] = np.linspace(0, 1, 10)  # Set values using numpy

Example: Using a matrix as a numpy array

mat = SUNDenseMatrix(3, 3, sunctx)
arr = SUNDenseMatrix_Data(mat)
arr[:] = np.eye(3)  # Set to identity matrix

This allows you to use numpy operations for vector and matrix data, and to pass numpy arrays to and from SUNDIALS routines efficiently and without unnecessary copies.

22.1.4.4. User-Supplied Callback Functions

SUNDIALS packages and several modules/classes require user-supplied callback functions to define problem-specific behavior, such as the right-hand side of an ODE or a nonlinear system function. In sundials4py, you can provide these as standard Python functions or lambdas.

The callback signatures follow the C API. As such, N_Vector arguments are passed as N_Vector objects and the underlying ndarray must be extracted in the user code. The only caveat is that return-by-pointer parameters are removed from the signature, and instead become return values (mirroring how return-by-pointer parameters for other functions are handled)

Warning

The C function signatures for most callbacks include a void* user_data argument. In Python, this argument must be present in the signature, but it should be ignored to avoid catastrophic errors. We recommend using _ as the parameter name in the callback signature to indicate this argument is unused.

Example: ODE right-hand side for ARKStep

# The C signature is:
# int(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data)
def rhs(t, y_nvector, ydot_nvector, _): # note _ in place of user_data
   # Compute ydot = f(t, y)
   y = N_VGetArrayPointer(y_nvector)
   ydot = N_VGetArrayPointer(ydot_nvector)
   ydot[:] = -y
   return 0

ark = ARKStepCreate(rhs, None, t0, y, sunctx)

Example: Nonlinear system for KINSOL

# The C signature is:
# int(N_Vector u, N_Vector g, void* user_data)
def fp_function(u_nvector, g_nvector, _): # note _ in place of user_data
   # Compute g = F(u)
   u = N_VGetArrayPointer(u_nvector)
   g = N_VGetArrayPointer(g_nvector)
   g[:] = u**2 - 1
   return 0

kin = KINCreate(sunctx)
KINInit(kin.get(), fp_function, u)

Example: ARKODE LSRKStep dominant eigenvalue estimation function with return-by-pointer parameters

# The C signature is:
# int(sunrealtype t, N_Vector y, N_Vector fn,
#     sunrealtype* lambdaR, sunrealtype* lambdaI,
#     void* user_data,
#     N_Vector temp1, N_Vector temp2, N_Vector temp3)
def dom_eig(t, yvec, fnvec, _, temp1, temp2, temp3): # note the _ in place of user_data
     lamdbaR = L
     lamdbaI = 0.0
     # lambdaR and lambdaI should be returned in the order that they appear
     # as parameters in the C API and follow the error code to return
     return 0, lamdbaR, lamdbaI

22.1.4.5. Error Codes

The named SUN_ERR_* code constants are not available in Python. However, all negative values of SUNErrCode are still errors, zero is success, and positive values are warnings. As such, users Users can call SUNGetErrMsg from Python with the returned SUNErrCode to get further information about an error.