Time Integration

Time integration is the numerical engine that evolves the physical state of a simulation forward in time. In SOFAx, we focus on stiff, constrained mechanical systems, which necessitate implicit integration schemes to ensure stability and precise constraint satisfaction.

Implicit Euler

sequenceDiagram autonumber actor User as Simulation Loop participant TS as TimeStepper participant Solver as SolverNode participant Newton as NewtonSolver participant State as NodeDOFs participant Constr as ConstraintField note over User,Constr: Goal: Advance simulation from t_n to t_{n+1} using implicit Euler
Kinematics: v_{n+1} = v_n + Δt·a_{n+1}, u_{n+1} = u_n + Δt·v_{n+1} User->>TS: step(scene, dt) TS->>TS: t_{n+1} = t_n + dt loop for each SolverNode rect rgba(227,242,253,0.25) note over TS,Constr: A) Prepare trial state TS->>Solver: get current state (u_n, v_n, a_n) Solver->>State: read DOFs at t_n State-->>TS: (u_n, v_n, a_n) TS->>Solver: get constraint multipliers Solver->>Constr: read λ_n Constr-->>TS: λ_n TS->>TS: init guess: y⁽⁰⁾ = (a_n, λ_n) note over TS: Unknown y = (a_{n+1}, λ_{n+1}) end rect rgba(232,245,233,0.25) note over TS,Newton: B) Solve coupled system F(a_{n+1}, λ_{n+1}) = 0 TS->>Newton: solve(coupled_residual_fn, y⁽⁰⁾, dt, tol) note over Newton: Newton-Krylov on coupled system:
F = [R_u(a,λ), c(u)] Newton->>Newton: iterative solve until convergence Newton-->>TS: (a_{n+1}, λ_{n+1}) = y* end rect rgba(255,243,224,0.25) note over TS,State: C) Update kinematics via implicit Euler TS->>TS: v_{n+1} = v_n + dt · a_{n+1} TS->>TS: u_{n+1} = u_n + dt · v_{n+1} TS->>State: write new DOFs (u_{n+1}, v_{n+1}, a_{n+1}) State-->>TS: DOFs updated TS->>Constr: write new multipliers λ_{n+1} Constr-->>TS: multipliers updated end rect rgba(237,231,246,0.25) note over TS,Solver: D) Update solver state TS->>Solver: increment time: time = t_{n+1} TS->>Solver: increment step: step_id = step_id + 1 Solver-->>TS: state updated end end TS-->>User: scene (updated at t_{n+1}) note over User,Constr: All solver nodes advanced to t_{n+1}
Both DOFs and constraint multipliers updated

Continuous Dynamics

Formally, the dynamics of a constrained mechanical system are described by a Differential-Algebraic Equation (DAE) of index 3. Letting \(q\) denote positions and \(v\) velocities, the continuous evolution is governed by:

\[ \begin{aligned} \dot{q} &= v \\ M(q) \dot{v} &= f(q, v, t) + G(q)^T \lambda \\ \phi(q) &= 0 \end{aligned} \]

Here:

\(M(q)\) is the mass matrix.
\(f(q, v, t)\) collects internal and external forces.
\(\phi(q) = 0\) represents holonomic geometric constraints.
\(\lambda\) are the Lagrange multipliers enforcing these constraints.
\(G(q) = \partial \phi / \partial q\) is the constraint Jacobian.

Discretization Schemes

To approximate the flow map \(\Psi_{t, t+\Delta t}\), we discretize the time derivative. Unlike explicit solvers which predict the future state from the current one, implicit solvers couple the future state unknowns, requiring the solution of a system of equations at each step.

SOFAx formulates the integration step as finding the accelerations \(a_{n+1}\) and multipliers \(\lambda_{n+1}\) that satisfy the discrete equations of motion.

Implicit Euler

The implicit Euler scheme is a first-order, L-stable method. It is highly robust and introduces numerical damping, making it ideal for quasistatic simulations or very stiff systems.

Update Rule:

\[ \begin{aligned} v_{n+1} &= v_n + \Delta t \, a_{n+1} \\ q_{n+1} &= q_n + \Delta t \, v_{n+1} \end{aligned} \]

Newmark-Beta

The Newmark family provides a generalized second-order framework. By tuning parameters \(\beta\) and \(\gamma\), one can select between different integration characteristics.

Update Rule:

\[ \begin{aligned} v_{n+1} &= v_n + \Delta t \left[ (1-\gamma)a_n + \gamma a_{n+1} \right] \\ q_{n+1} &= q_n + \Delta t \, v_n + \Delta t^2 \left[ (\tfrac{1}{2}-\beta)a_n + \beta a_{n+1} \right] \end{aligned} \]

The Coupled Root-Finding Problem

In implicit integration, the unknowns \((a_{n+1}, \lambda_{n+1})\) appear on both sides of the equation (implicitly via forces and constraints). This transforms the integration step into a nonlinear root-finding problem.

We define a coupled residual function \(R(z)\) where \(z = (a_{n+1}, \lambda_{n+1})\):

\[ R(z) = \begin{pmatrix} M a_{n+1} - f(q_{n+1}, v_{n+1}, t_{n+1}) - G(q_{n+1})^T \lambda_{n+1} \\ \phi(q_{n+1}) \end{pmatrix} = 0 \]

Solving \(R(z) = 0\) is the core computational task of the simulator. In SOFAx, this is handled by a Newton-Krylov solver, which linearizes \(R(z)\) and solves the resulting coupled primal-dual linear system:

\[ \begin{bmatrix} A & G^T \\ G & 0 \end{bmatrix} \begin{pmatrix} \delta a \\ \delta \lambda \end{pmatrix} = - R(z_k) \]

This system is solved matrix-free using iterative methods (GMRES), fully preserving the coupled nature of the problem without resorting to decoupled free-motion steps.

Integration with Diffrax

SOFAx assumes a structure compatible with diffrax.diffeqsolve, mapping our mechanical components to standard ODE solver abstractions:

SOFAx Concept	Diffrax Abstraction	Description
System State	`y` (PyTree)	The collection of all DOFs \((q, v)\), typically stored as a PyTree.
Physics / Residual	`AbstractTerm`	The operator defining the continuous dynamics (forces, constraints).
Integrator	`AbstractSolver`	The numerical scheme (Euler, Newmark) that defines the discrete step update.
Newton Solve	`nonlinear_solver`	The root-finding algorithm invoked within the implicit solver's step.

Benefits of this architecture

Composable Solvers: Easily swap between Euler for stability and high-order Runge-Kutta methods for precision.
Adjoint Sensitivity: Leverage diffrax's backpropagation capabilities for gradient-based optimization and learning.
Neural ODEs: Seamlessly mix neural networks (as ODETerms) with physical laws, treating the entire simulation as a differentiable layer.