Skip to content

Residual operator

The core numerical object in SOFAx v2 is the residual operator. All nonlinear solves, linearizations, and differentiation are expressed in terms of residual evaluation and its directional derivatives.

Rather than assembling global matrices, the solver evaluates a residual function \(R(x)\) and applies its Jacobian through automatic differentiation.

SolverNode Residual Operator
sequenceDiagram autonumber actor Solver as SolverNode
ResidualEval participant Walk as TreeWalker participant Map as Mapping participant Field as PhysicalField participant Constr as ConstraintField participant Acc as Accumulator note over Solver,Acc: Goal: evaluate solver residual R(x) and constraint violations c(x)
by recursive subtree traversal. All children DOFs are mapped to solver DOFs. Solver->>Walk: residual(x_trial) rect rgba(227,242,253,0.25) note over Walk,Map: A) Downward pass
propagate trial kinematics via P Walk->>Walk: visit(node) loop for each MappedNode Walk->>Map: apply_P(x_parent) Map-->>Walk: x_child (projected DOFs) Walk->>Walk: visit(child) recursively end end rect rgba(232,245,233,0.25) note over Walk,Field: B) Local physics evaluation
compute field contributions Walk->>Walk: get node.fields loop for each PhysicalField Walk->>Field: evaluate(node, x_trial) Field->>Field: compute Ma, Bv, Ku, f Field-->>Acc: FieldContribution Acc->>Acc: sum physics residual:
R_phys += Ma + Bv + Ku - f end end rect rgba(237,231,246,0.25) note over Walk,Constr: B') Local constraint evaluation
accumulate constraint violations loop for each ConstraintField Walk->>Constr: evaluate(node, x_trial) Constr->>Constr: compute constraint c(u) Constr->>Constr: compute constraint gradient ∇c Constr-->>Acc: constraint contribution Acc->>Acc: accumulate c_node(u)
store in constraint buffer end end rect rgba(255,243,224,0.25) note over Walk,Acc: C) Upward pass (Folding)
map child residuals to parent space loop for each MappedNode Walk->>Walk: get R_child, c_child Walk->>Map: apply_PT(R_child) Map-->>Acc: R_child in parent coords Acc->>Acc: R_parent += PT · R_child note over Walk,Acc: Constraints accumulate directly
constraint index, no mapping Acc->>Acc: c_parent += c_child end Acc-->>Walk: (R_subtree, c_subtree) end Walk-->>Solver: (R_phys(x), c(x)) at solver root note over Solver,Acc: Returns: physics residual R_phys(x) and constraint buffer c(x)

Role of the residual

At the solver level, the residual encodes:

  • the discretized equations of motion
  • constraint forces induced by Lagrange multipliers
  • violations of holonomic constraints

A Newton step seeks a root of the coupled system:

\[ F(y) = \begin{bmatrix} R_u(y) \\ \Phi(u) \end{bmatrix} = 0, \qquad y = (a, \lambda). \]

Here, \(R_u\) denotes the physics residual augmented with constraint reaction forces \(G(u)^\top \lambda\), while \(\Phi(u)\) measures constraint violations.

Tree-based evaluation

Residual evaluation is performed by a recursive traversal of the scene graph, starting from a solver node and visiting all its descendants.

The traversal follows a fixed three-phase pattern:

  1. Downward pass
    Trial kinematics are propagated from solver DOFs to all descendant nodes using mapping operators \(P\).

  2. Local evaluation
    Each node evaluates its physical fields and constraint fields in its own local coordinates.

  3. Upward pass
    Residual contributions are mapped back to solver coordinates using transpose operators \(P^\top\) and accumulated.

This folding pattern guarantees that all contributions are ultimately expressed in the solver DOF space.

Downward pass

Given a trial state \(x\) at the solver node:

  • child DOFs are computed as \(x_c = P\,x_p\)
  • the mapping is applied recursively along the scene graph

No allocation or in-place mutation occurs: projected states are passed by value through the traversal, preserving functional purity.

Local physics and constraint evaluation

At each node, physical fields contribute to the residual in vector form:

\[ R_{\mathrm{phys}} = M\,a + B\,v + K\,u - f_{\mathrm{ext}}. \]

Constraint fields locally evaluate:

  • constraint violations \(\Phi(u)\)
  • constraint gradients \(G(u)\)

Each field produces a local contribution object containing only vector terms (residuals, forces, constraint values), never assembled matrices.

Upward pass

For mapped nodes:

  • child residuals are mapped back to parent coordinates as \(P^\top R_c\)
  • constraint contributions are accumulated directly by index (no mapping)

The parent residual is formed by summing:

  • its own local contributions
  • all mapped contributions coming from its children

Differentiability

Because the residual operator is:

  • a pure function
  • expressed using JAX-compatible primitives
  • free of side effects and mutation

it can be differentiated through:

  • unrolled Newton or Krylov iterations
  • implicit differentiation
  • adjoint-based inverse problems

All solver components are built by composing and differentiating this single residual operator.

See also