efficient implicit timestepping

To address the insufficient time-integration capabilities of the SUNDIALS numerical library in high-performance scientific computing, this project systematically extends its time-stepping solvers. Methodologically, it introduces three novel classes of single-step methods—low-storage Runge–Kutta (LSRK), symplectic structure-preserving block RK, and general-purpose operator-splitting schemes—alongside a new multi-rate adaptive step-size controller and explicit RK-based adjoint sensitivity analysis (filling a longstanding gap). It further enhances nonlinear solvers with Anderson acceleration and improves error handling and logging infrastructure. These contributions significantly improve efficiency, stability, and accuracy for large-scale transient simulations, enabling long-duration, high-fidelity, and multiphysics-coupled modeling. Validation across multiple HPC applications demonstrates speedups of 1.5–3× and markedly improved numerical robustness.

This work addresses the complexity of developing numerical time integrators and the difficulty of algorithm switching on distributed heterogeneous platforms. We propose a deep integration of the Odeint and OpenFPM frameworks, leveraging template metaprogramming and a declarative integration interface to uniformly abstract multistage, multistep, and adaptive ODE solvers as plug-and-play components—enabling algorithm switching with a single line of code while preserving MPI+GPU parallel scalability. To our knowledge, this is the first unified architecture synergistically combining Boost.Odeint, OpenFPM, CUDA, and MPI. Within just 60 lines of core code, we achieve efficient CPU/GPU-accelerated simulation of the 3D Gray–Scott reaction-diffusion model, supporting both exponential and sigmoidal kinetics. Experimental evaluation demonstrates strong scalability and high computational efficiency across heterogeneous clusters.

Explicit Material Point Method (MPM) suffers from numerical instability under large time steps, hindering its integration with partitioned large-step solvers and multi-physics frameworks. Method: This paper proposes a plug-and-play substepping algorithm that encapsulates explicit MPM as a pseudo-implicit scheme without modifying the underlying explicit integrator code. By introducing internal substeps coupled with constraint enforcement and projection operations, the method ensures numerical stability and physical consistency at macroscopic large time steps. Contribution/Results: The approach is inherently compatible with multi-solver coupling, complex constraint handling, and multi-physics integration, significantly enhancing computational robustness and efficiency. Experiments demonstrate that, while preserving accuracy, the method enables time-step enlargement by several-fold—providing a practical pathway for embedding MPM into large-scale, multi-physics simulation frameworks.

For stiff nonlinear dynamical systems, implicit time integration suffers from slow convergence and high computational cost due to repeated nonlinear equation solving at each time step. This paper proposes a deep learning–enhanced hybrid Newton method to address these challenges. The core innovation is an unsupervised, target-oriented learning strategy specifically designed to optimize Newton’s initial guess—requiring no labeled data while enabling neural networks to generate highly accurate starting points. We theoretically derive both a convergence acceleration bound and a generalization error upper bound. By tightly integrating deep neural networks with classical Newton iteration, the method significantly reduces the number of iterations (by 40–60%) on benchmark 1D and 2D stiff problems, while preserving numerical stability and solution accuracy. This work establishes a novel, efficient, and robust solver paradigm for implicit time stepping in stiff nonlinear dynamics.

This work addresses the high computational cost of gradient-based optimization for transient convection–diffusion problems using reduced-order models (ROMs). We propose a “optimize-then-reduce” coupled framework wherein the primal and adjoint systems are simultaneously solved within the reduced space at each time step. To alleviate the dependency on high-dimensional adjoint snapshots, we devise an efficient adjoint snapshot collection strategy. Furthermore, we overcome conventional limitations in adjoint basis construction by introducing an adaptive adjoint basis selection method guided by energy decay, iteration count, and computational time. Numerical experiments demonstrate that the proposed approach significantly reduces the computational overhead of adjoint system solves and overall simulation time, while maintaining controllable error levels. The method enables real-time, high-fidelity optimization in ROM–ROM coupled settings.

This paper addresses the efficiency–accuracy trade-off in quantized state system (QSS) methods for numerically solving ordinary differential equations (ODEs). We propose the Generalized Linear Implicit Quantized State System (GLIQSS) framework—the first extension of Linear Implicit QSS (LIQSS) to non-uniform quantization and higher-order linear implicit integration structures—yielding two novel algorithm families. GLIQSS integrates state quantization, event-driven simulation, and rigorous error and stability analysis, guaranteeing global error bounds and unconditional stability while substantially reducing event-processing overhead. Experimental results on two representative applications demonstrate that GLIQSS achieves significant computational speedups over classical methods such as RK4 and BDF, reduces event-triggering frequency by over 30%, and maintains superior numerical accuracy and robustness.

This work proposes a stabilization framework based on score-based generative models to address non-physical instabilities and structural distortions commonly encountered in the numerical solution of time-dependent partial differential equations (PDEs). For the first time, score-based generative models are introduced into PDE numerical stabilization, where a conditional stabilizing operator with manifold-contracting properties is constructed by learning the physically admissible solution manifold. This operator corrects intermediate solutions during time integration. Numerical experiments demonstrate that the method significantly enhances robustness for convection, Korteweg–de Vries (KdV), nonlinear Schrödinger, and Burgers equations, effectively suppressing spurious oscillations while preserving essential dynamical features.

This work addresses the high memory and computational complexity typically associated with solving three-dimensional partial differential equations on Cartesian grids. By exploiting tensor-product structure, the proposed method decomposes the 3D operator into one-dimensional banded kernels aligned with coordinate axes, thereby avoiding explicit assembly of the global matrix and enabling a matrix-free solution strategy. Within a unified framework that integrates diverse numerical approaches—including Kronecker product algebra, compact finite differences, isogeometric analysis, and direct diagonalization—the study systematically identifies three key techniques: multi-right-hand-side reshaping, sum factorization, and pencil-style MPI decomposition. These innovations collectively enhance hardware affinity and parallel scalability, reducing algorithmic complexity to O(N) and storage requirements to O(Nₓ + Nᵧ + N_z), thus enabling efficient large-scale 3D PDE simulations.

This work addresses the high computational cost of solving high-dimensional partial integro-differential equations (PIDEs), which arises from their nonlocal jump terms. The authors propose a mesh-free iterative neural solver that reformulates PIDE solution as a recursive regression problem over the entire space-time domain. By leveraging a single-jump Monte Carlo sampling strategy, the method implicitly handles the nonlocal integral term, thereby avoiding explicit integration and full residual differentiation. Embedded within the physics-informed neural networks (PINNs) framework and augmented with an iterative learning strategy, the approach efficiently captures the global solution. The method is theoretically guaranteed to converge via contraction mapping for linear PIDEs and demonstrates high accuracy and strong scalability across multiple high-dimensional linear and nonlinear test cases.

This work addresses the high memory overhead and neglect of local structure in gradient computation for implicit nonlinear solvers within differentiable simulation. The authors propose a solver-level differentiation method that constructs an adjoint algorithm symmetric to the forward solve by reverse-scanning a block-structured implicit solver, entirely avoiding the assembly of a global Jacobian matrix. For the first time, adjoint computation is aligned with the block structure of the forward solver, combining vertex-block descent with reverse-colored Gauss–Seidel sweeps to enable efficient backpropagation using only local 3×3 adjoint solves. This approach leverages operator-view approximations of the inverse and its transpose. On a single GPU, it achieves a 33× speedup and 71× reduction in memory compared to unrolled automatic differentiation, enabling, for the first time, differentiable elastic dynamics simulation of million-contact coupled soft bodies with up to 8 million vertices.