PaperlayerView exact revision on arXiv

Source-linked AI summary

Enabling New Flexibility in the SUNDIALS Suite of Nonlinear and Differential/Algebraic Equation Solvers

David J. Gardner, Daniel R. Reynolds, Carol S. Woodward, Cody J. Balos

arXiv:2011.10073v2cs.MS

TL;DR

SUNDIALS needed more flexible interfaces for application-specific and third-party solvers and data structures as HPC systems became more heterogeneous. The paper redesigns solver and vector abstractions and modernizes Fortran interfaces, retaining existing behavior while enabling customized approaches and reported speedups. The principal scope limitation is that one tolerance adjustment does not balance errors among individual residual or solution components.

  • Problem

    Heterogeneous HPC systems increased demand for varied programming models and algebraic solvers, while existing interfaces limited extensibility and Fortran flexibility.

  • Method

    The paper adds object-oriented matrix, linear-solver, and nonlinear-solver classes, expands vector implementations, and creates modern Fortran 2003 interfaces.

  • Results

    Numerical tests retain existing integrator behavior without increasing simulation cost, while problem-specific nonlinear solvers exploiting system structure and locality produce significant speedups.

  • Takeaways & Limitations

    The infrastructure enables reusable external and application-specific solvers and easier comparison of different integrators and solver approaches within the same problem.

  • Takeaways & Limitations

    The tolerance adjustment accounts only for overall linear-residual magnitude and does not balance errors among specific residual or solution components.

Abstract

from arXiv · show

In recent years, the SUite of Nonlinear and DIfferential/ALgebraic equation Solvers (SUNDIALS) has been redesigned to better enable the use of application-specific and third-party algebraic solvers and data structures. Throughout this work, we have adhered to specific guiding principles that minimized the impact to current users while providing maximum flexibility for later evolution of solvers and data structures. The redesign was done through the addition of new linear and nonlinear solvers classes, enhancements to the vector class, and the creation of modern Fortran interfaces. The vast majority of this work has been performed "behind-the-scenes," with minimal changes to the user interface and no reduction in solver capabilities or performance. These changes allow SUNDIALS users to more easily utilize external solver libraries and create highly customized solvers, enabling greater flexibility on extreme-scale, heterogeneous computational architectures.

1 INTRODUCTION

SUNDIALS provides efficient solvers for time-dependent and nonlinear equations on large-scale HPC systems, but evolving heterogeneous architectures created demand for more adaptable solver interfaces. This paper presents enhancements designed to support new programming environments and architecture-aware, problem-specific solvers while minimizing user disruption.

  • SUNDIALS is a collection of packages for solving time-dependent and nonlinear equations on large-scale HPC systems.
  • Heterogeneous architectures increased demand for interfaces supporting new programming models and algebraic solver packages.Existing interfaces would require extensive redundant code with limited extensibility.
  • The redesign emphasizes adaptability to new machines and programming environments while enabling architecture-aware and problem-specific solver construction.
  • The paper focuses on interfaces and functionality, omitting many implementation details while directing readers to user guides for them.

2 SUNDIALS OVERVIEW

SUNDIALS contains six packages covering ODEs, DAEs, sensitivity analysis, linearly implicit and multirate integration, and nonlinear systems. Their methods share vector, matrix, linear-solver, and nonlinear-solver operations needed by time integration and nonlinear solution.

  • The six packages address ODEs, DAEs, sensitivity analysis, linearly implicit and multirate integration, and nonlinear systems.They are CVODE, CVODES, IDA, IDAS, ARKODE, and KINSOL.
  • CVODE and CVODES: CVODE uses variable-order, variable-step implicit multistep methods, with Adams-Moulton orders 1 to 12 for nonstiff problems and BDF orders 1 to 5 for stiff cases.
  • CVODE and CVODES: CVODES and IDAS extend CVODE and IDA with forward and adjoint sensitivity analysis for ODE and DAE initial value problems.
  • ARKODE: ARKODE provides explicit, implicit, and IMEX additive Runge-Kutta methods, plus multirate methods that advance slow and fast components at different time scales.
  • KINSOL: KINSOL solves root-finding and fixed-point systems using Newton, fixed-point, Picard, and optionally Anderson-accelerated iterations.
  • Implicit integrators solve nonlinear systems, whose generic root-finding or fixed-point forms can use methods shared with KINSOL.

3 SUNDIALS STRUCTURE

SUNDIALS organizes shared numerical operations behind abstract vector, matrix, linear-solver, and nonlinear-solver classes. The redesigned structure removes package-specific duplication and supports reusable, external, and application-specific solver implementations through object-oriented interfaces.

  • SUNDIALS packages require vector, nonlinear-system, linear-system, and matrix operations, motivating four corresponding abstract base classes.
  • The four base classes are N_Vector, SUNMatrix, SUNLinearSolver, and SUNNonlinearSolver.
  • Abstract interfaces hide implementation-specific data structures, layouts, parallelization methods, and algebraic solvers from SUNDIALS packages.
  • Previous SUNDIALS structure: Earlier linear-solver interfaces duplicated package-specific code, while nonlinear solvers were embedded in integrators and lacked external or user-defined solver interfaces.
  • Redesign goals: The redesign collapses linear-solver interfaces, unifies nonlinear solvers, preserves integrator-specific controls, and streamlines use of external or application-specific solvers.
  • SUNDIALS Class Structure: Each base class combines derived-class member data with a virtual method table of function pointers.
  • SUNDIALS Class Structure: Pure virtual functions define required operations, while optional functions can be overridden or used to improve performance.
  • SUNDIALS Class Structure: Constructors, cloning utilities, and method tables let derived classes provide custom implementations with minimal impact on current users.

4 ENHANCED VECTOR STRUCTURES

SUNDIALS enhances its vector structures with fused operations and many-vector implementations to reduce computational overhead and support heterogeneous, multiphysics simulations. These extensions preserve customization while allowing data and computation to be partitioned across distinct resources and communicators.

  • New Vector Operations: Fused vector operations combine repeated computations to reduce memory accesses, kernel launches, and MPI reductions without changing the total floating-point work.They increase arithmetic intensity and reduce auxiliary overhead.
  • New Vector Operations: Optional fused operations let custom N_Vector implementations provide only the routines most impactful for their use case while retaining equivalent base-class behavior.The operations can be enabled or disabled at runtime through VMT function pointers.
  • New Vector Operations: Standalone CPU tests found fused streaming operations beneficial for many-vector workloads above approximately 1,000 unknowns per MPI task, while fused reductions helped with fewer than approximately 130K unknowns per MPI task.These thresholds are machine dependent.
  • Many-vector: Many-vector implementations combine distinct vector instances, potentially using different parallelization approaches, into one cohesive vector without directly manipulating their data.Underlying vector instances perform the computations, preserving control over data placement and resources.
  • Many-vector: Many-vector designs support partitioning data across heterogeneous resources, separate MPI communicators for coupled simulations, user-supplied subvector routines, and new partitioned or block solver implementations.They are intended to function seamlessly with existing SUNDIALS packages and classes.
  • Many-vector: Reactive-flow data can partition fluid variables into parallel vectors and chemical species into task-specific serial, OpenMP, or GPU vectors.Many-vector use cases also include combining multiphysics simulations with distinct MPI intracommunicators and constructing hierarchical communication patterns.
  • Many-vector: The MPIManyVector WRMS algorithm uses 1 global reduction instead of s+1 reductions when each WRMS_Norm call requires an MPI_Allreduce.This reduces global communication for the collection of subvectors.
  • Many-vector: MPIPlusX provides a thin wrapper around MPIManyVector to simplify construction of hybrid MPI+X vectors using an MPI communicator and an on-node N_Vector.The on-node vector may use SUNDIALS or user-defined parallelism such as OpenMP, CUDA, RAJA, OpenCL, or OpenACC.

13 MPI_Allreduce(&local_sum, &global_sum, 1, MPI_DOUBLE, MPI_SUM, inter_comm)

The MPIPlusX implementation wraps MPIManyVector to simplify hybrid MPI+X vector construction and provide utility functions for the wrapped on-node vector.

  • MPIPlusX: MPIPlusX is a thin wrapper of MPIManyVector that simplifies construction and provides utility functions for the wrapped on-node vector.Construction supplies an MPI communicator and an on-node N_Vector.
  • MPIPlusX: The wrapped on-node vector can support GPU-based implementations whose performance is compared beneath the MPIPlusX vector.The passage refers readers to a separate performance comparison of various GPU on-node vectors.
  • MPIPlusX: The construction pattern combines an MPI communicator with an MPI-unaware on-node N_Vector to create a hybrid MPI+X vector.MPIManyVector supplies the underlying mechanism.

5 LINEAR SOLVERS

SUNDIALS introduces a generic SUNLinearSolver interface that supports application-specific and third-party linear solvers while preserving integrator-specific optimizations. The class accommodates direct, iterative, matrix-based, and matrix-free approaches, including scaling, preconditioning, and a reusable hypre integration demonstrated on a two-dimensional heat equation.

  • SUNLinearSolver class: SUNLinearSolver requires implementations of SUNLinSolGetType and SUNLinSolSolve, while optional routines provide integrators with efficiency-enhancing controls and solver artifacts.SUNLinSolGetType identifies the solver category, and SUNLinSolSolve computes a solution to A x = b.
  • SUNLinearSolver class: The solver type distinguishes matrix-based direct, matrix-based iterative, and matrix-free iterative methods, enabling SUNDIALS to select solver-specific optimizations and tolerances.For iterative methods, SUNDIALS selects tolerances intended to support outer Newton iterations without unnecessary linear iterations.
  • Integrator-specific optimizations: Optional SUNLinSolSetup amortizes costly matrix factorizations or preconditioner setups by allowing them to be performed less frequently than linear solves.Without this routine, setup costs are incurred at every SUNLinSolSolve call, which occurs more frequently.
  • Matrix interfaces: A custom matrix-free solver receives a function for computing A x → z, whereas matrix-based solvers use compatible SUNMatrix implementations with required and optional operations.The custom hypre example defines SUNMatrix and SUNLinearSolver data structures and implements only the operations needed for its iterative matrix-based solver.
  • Iterative solver support: Scaling and preconditioning routines support transformed iterative systems that balance component errors and accelerate convergence.SUNDIALS supplies scaling matrices and preconditioning operators; unsupported scaling can be approximated through an adjusted tolerance, but this does not balance individual components.
  • Demonstration problem: In the heat-equation demonstration, CVODE with a matrix-based hypre solver is most efficient at the desired accuracy, while the new infrastructure preserves reusable interfaces across integrators.The comparison includes CVODE and ARKODE with SUNDIALS and hypre solvers; the redesign removes duplicate interfaces and leaves performance within observed variability.

6 NONLINEAR SOLVERS

SUNDIALS introduces a generic nonlinear-solver interface that supports required and optional operations, enabling custom and third-party solvers within its integrators. A task-local solver exploiting spatial locality substantially outperforms the global Newton approach in the demonstration.

  • SUNNonlinearSolver class: The SUNNonlinearSolver class supplies SUNDIALS integrators with problem-specific and third-party nonlinear solvers through required and optional methods.Required methods describe the system type, provide its function, and perform the solve; optional methods support efficiency, configuration, and solver information.
  • SUNNonlinearSolver class: Integrators can pass setup and solve functions for linear systems, while nonlinear solvers that handle these operations internally may omit them.The optional SetLSetupFn and SetLSolveFn methods support integrator-specific actions and reuse SUNLinearSolver objects.
  • Integration and flexibility: The redesigned integrators support generic nonlinear solvers, including PETSc SNES interfaces and custom implementations.Existing solver implementations were consolidated into shared derived classes, expanding the available solver choices.
  • Demonstration problem: The demonstration applies a custom solver to independent task-local nonlinear systems in a stiff one-dimensional advection-reaction problem.The problem uses ARKODE’s third-order IMEX method, with explicit advection and implicit reaction terms; local solves reduce nearly all communication during implicit solves.
  • Demonstration results: 4x to almost 8x speedups are achieved by the custom task-local solver compared with the global Newton solve.The comparison uses average execution times from 20 runs across varying problem sizes on Lassen CPUs.

7 FORTRAN 2003 INTERFACES

SUNDIALS replaces limited Fortran 77 bindings with Fortran 2003 interfaces that expose all six packages and most core classes. The interfaces are generated around ISO C interoperability, closely mirror the C API, and support custom solver implementations.

  • New interfaces: The new Fortran 2003 interfaces provide access to all features in all six SUNDIALS packages and most core classes.They are created using SWIG-Fortran with minimal additional interface code for ISO C interoperability.
  • Maintenance: Generated interface glue code keeps the Fortran 2003 bindings low-maintenance and readily extensible when C APIs change.The glue code is regenerated before package releases and could be integrated into continuous delivery workflows.
  • Usage: The Fortran 2003 API closely resembles the C API while allowing users to write idiomatic Fortran with few caveats.The C API documentation therefore serves as documentation for the Fortran 2003 interface, with notes on data types and subtleties.
  • Demonstration: A Fortran implementation of the task-local nonlinear-solver demonstration follows the same derived-class steps as the C version.The example stores custom solver data in a user-defined module, while a new solver-data type offers a thread-safe alternative.

8 CONCLUSIONS

SUNDIALS’s redesigned infrastructure targets flexible, sustainable integration of application-specific and external solvers across increasingly complex HPC environments. The additions preserve existing integrator behavior while enabling performance gains from problem-specific solver approaches.

  • Design goals: New matrix, linear-solver, and nonlinear-solver classes address growing demands for varied programming models and algebraic solvers on HPC systems.The design emphasizes greater flexibility, sustainability, and minimal changes to user codes.
  • Infrastructure: The new classes improve interoperability with application-specific and external data structures and solvers across SUNDIALS packages.The many-vector and MPI+X vector support accommodates multiphysics simulations, heterogeneous architectures, and data in different memory spaces.
  • Infrastructure: New fused vector operations support hardware utilization and reduced-communication algorithms, while modern Fortran interfaces expose the class flexibility.The many-vector construct also supports vectors operating on potentially different memory spaces.
  • Results: Numerical tests retain existing time-integrator behavior without increasing simulation cost and show significant speedups for structure-exploiting nonlinear solvers.The infrastructure also enables direct performance comparisons between fundamentally different integrators and solver approaches within one problem.
  • Future work: Future work will add vectors and matrices for more on-node parallelism, investigate fused-operation algorithms, expand native solver interfaces, and exploit many-vector structure.These plans target broader support for application-specific parallelism and structural features.
Loading 2011.10073v2…