  Many complex systems can be accurately modeled as a set of coupled
time-dependent partial differential equations (PDEs). However, solving such
equations can be prohibitively expensive, easily taxing the world's largest
supercomputers. One pragmatic strategy for attacking such problems is to split
the PDEs into components that can more easily be solved in isolation. This
operator splitting approach is used ubiquitously across scientific domains, and
in many cases leads to a set of ordinary differential equations (ODEs) that
need to be solved as part of a larger "outer-loop" time-stepping approach. The
SUNDIALS library provides a plethora of robust time integration algorithms for
solving ODEs, and the U.S. Department of Energy Exascale Computing Project
(ECP) has supported its extension to applications on exascale-capable computing
hardware. In this paper, we highlight some SUNDIALS capabilities and its
deployment in combustion and cosmology application codes (Pele and Nyx,
respectively) where operator splitting gives rise to numerous, small ODE
systems that must be solved concurrently.
