21.2. Using sundials4py

At a high level, using SUNDIALS from Python via sundials4py looks a lot like using SUNDIALS from C or C++.

The few notable differences are:

21.2.1. View Classes and Memory Management

sundials4py provides natural usage of SUNDIALS objects with natural Python object lifetimes managed by the Python garbage collection as with any other Python object. There is only one caveat that arises with void* objects/variables due to restrictions in nanobind: the SUNDIALS integrator/solver memory void* objects are (behind the scenes) wrapped in “View” classes. 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:

from sundials4py.core import *
from sundials4py.cvode import *

sunctx = SUNContext_Create(SUN_COMM_NULL)
cvode = CVodeCreate(CV_BDF, sunctx)
# notice we need to call cvode.get()
status = CVodeInit(cvode.get(), lambda t, y, ydot, _: ode_problem.f(t, y, ydot), T0, y)

21.2.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:

First element: The function’s return value (typically an error code).
Subsequent elements: 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 CVodeGetNumSteps(void *cvode_mem, long int *numsteps);

Python:

retval, numsteps = CVodeGetNumSteps(cvode_mem.get())
print(f"Number of steps: {numsteps}")

Example 2: Multiple Return-by-Pointer Values

C:

int CVodeGetIntegratorStats(void *cvode_mem,
                           long int *nsteps,
                           long int *nfevals,
                           long int *nlinsetups,
                           long int *netfails);

Python:

retval, nsteps, nfevals, nlinsetups, netfails = CVodeGetIntegratorStats(cvode_mem.get())
print(f"Steps: {nsteps}, Function evals: {nfevals}, Linear setups: {nlinsetups}, Error test fails: {netfails}")

21.2.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 = NVectorView.Create(N_VNew_Serial(10, sunctx.get()))
y = N_VGetArrayPointer(y_nvec.get())
y[:] = np.linspace(0, 1, 10)  # Set values using numpy

Example: Using a matrix as a numpy array

mat = SUNMatrixView.Create(SUNDenseMatrix(3, 3, sunctx.get()))
arr = SUNDenseMatrix_Data(mat.get())
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.

21.2.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. Some things to note:

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)
Most callback signatures include a void* user_data argument. In Python, this argument must be present in the signature, but it should be ignored.

Example: ODE right-hand side for ARKStep

def rhs(t, y_nvector, ydot_nvector, _):
   # Compute ydot = f(t, y)
   y = N_VGetArrayPointer(y_nvector)
   ydot = N_VGetArrayPointer(ydot_nvector)
   ydot[:] = -y
   return 0

ark = ARKodeView.Create(ARKStepCreate(rhs, None, t0, y.get(), sunctx.get()))

Example: Nonlinear system for KINSOL

def fp_function(u_nvector, g_nvector, _):
   # Compute g = F(u)
   u = N_VGetArrayPointer(u_nvector)
   g = N_VGetArrayPointer(g_nvector)
   g[:] = u**2 - 1
   return 0

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

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, _):
     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

Warning

The user_data argument should always be None or _ on the Python side. If it is listed otherwise then it should be ignored to avoid causing catastrophic errors.

21.2.5. Examples

Examples can be found in bindings/sundials4py/examples.