3. CAE solver architectures

Large-scale CAE simulation performance is determined less by raw floating-point throughput than by the solver architecture chosen for a given physics, model size, and hardware topology. This chapter examines the four interlocking decisions that define that architecture: how time is advanced (implicit vs. explicit), how the resulting sparse linear systems are solved, how nonlinearity is driven to convergence, and how the model is partitioned across MPI ranks. Each decision constrains the others, so they must be evaluated as a system rather than in isolation.

Implicit and Explicit Time Integration

Explicit schemes such as the central-difference method advance the solution using only diagonal (lumped) mass information. The update for node i is local: accelerations depend on internal and external forces at the current step, so no global matrix solve is required. The price is conditional stability — the critical step is governed by the Courant–Friedrichs–Lewy (CFL) condition, Δt ≤ 2/ω_max, where ω_max is the highest natural frequency of the mesh. For a solid element of characteristic length L_e and dilatational wave speed c, this becomes Δt ≈ L_e/c. A steel model (c ≈ 5100 m/s) with 5 mm elements therefore requires Δt ≈ 1 µs, meaning a 10 ms event needs roughly 10⁴ steps. Explicit methods excel for crash, impact, blast, and short-duration nonlinear events where the physical time scale is itself small.

Implicit schemes (Newmark-β, generalized-α, or backward Euler) evaluate forces at the unknown time level, producing a global system KΔu = r at every step. They are unconditionally stable for linear problems, so Δt is chosen by accuracy rather than stability — often 10²–10⁴ times larger than the explicit limit. The trade-off is that each step requires assembling and factorizing a sparse Jacobian and solving it, which dominates cost for models beyond a few million DOF. The practical rule: use explicit when the event is short and highly nonlinear; use implicit when the event is long, quasi-static, or dominated by low-frequency response.

Sparse Linear Systems

Implicit solvers spend most of their time in the sparse linear solve. Two families dominate. Direct sparse solvers (MUMPS, SuperLU_DIST, PARDISO) perform a fill-reducing ordering (AMD, METIS nested dissection), symbolic factorization, and numeric factorization. They are robust and handle ill-conditioning well, but factorization cost scales roughly as O(n^1.5) to O(n²) in 3D and memory grows with fill-in. Iterative Krylov solvers (CG for symmetric positive-definite systems, GMRES or BiCGStab for nonsymmetric ones) require only matrix–vector products, scaling as O(n) per iteration, but need effective preconditioning. Algebraic multigrid (AMG, e.g., BoomerAMG in hypre) is the standard preconditioner for solid mechanics and Poisson-like problems; incomplete LU (ILU) or block-Jacobi works for more localized systems. A common production strategy is a hybrid: use a direct solver for small or stiff subsystems and an AMG-preconditioned Krylov method for the global system.

Nonlinear Iteration

Geometric, material, and contact nonlinearities require an outer Newton–Raphson loop. At iteration k, the residual is r^k = F_ext − F_int(u^k), and the tangent system K_T^k δu = r^k is solved for the correction δu. Convergence is declared when both the residual norm and the displacement correction norm fall below tolerances (typically 10⁻³–10⁻⁶ relative). Modified Newton, which reuses K_T across several iterations, trades slower convergence for fewer factorizations. Line search or trust-region damping is essential when contact status changes or material softening causes divergence. Arc-length (Riks) methods are used when the load–displacement path includes limit points.

Domain Decomposition Strategies

For models exceeding a single node's memory, the mesh is partitioned into subdomains, one per MPI rank. Two broad approaches exist. Overlapping Schwarz methods (additive Schwarz, RAS) exchange a layer of ghost elements and iterate on local solves; they converge quickly but duplicate interface work. Non-overlapping substructuring methods (FETI, BDDC, and their variants) solve local Neumann problems and enforce interface continuity through a coarse problem and Lagrange multipliers. BDDC/FETI-DP scale to tens of thousands of cores because the coarse space propagates global information in few iterations. The partitioner (METIS, ParMETIS, Zoltan) directly affects communication volume and load balance; a poor partition can double runtime even with perfect solvers.

Putting It Together: A Worked Sizing Example

Consider a 30-million-DOF automotive body-in-white model subjected to a 20 ms crash pulse. The smallest element is 4 mm in aluminum (c ≈ 5100 m/s), giving Δt ≈ 0.78 µs and about 25,600 explicit steps. Each step requires one internal-force evaluation and a diagonal mass update — roughly 30 MFLOP per step, or 0.77 TFLOP total. On 512 cores at 20% efficiency, this completes in minutes. If the same model were run implicitly with 0.2 ms steps (100 steps), each step would require a sparse factorization of a 30-million-DOF Jacobian. With AMG-preconditioned CG and BDDC, each solve might take 30–60 seconds, yielding 50–100 minutes — slower than explicit despite 256× fewer steps. This illustrates why the explicit/implicit choice must be made before tuning solvers.

For a 5-million-DOF static tire-contact problem, the picture reverses. Implicit Newton–Raphson with 8–12 iterations per load step, each solved by a direct sparse factorization (fill-in manageable at this size), completes in a few minutes. Explicit would require millions of steps to reach static equilibrium and is impractical. The domain decomposition here is trivial — a single node may suffice — but if the model grows to 50 million DOF, BDDC with 256 subdomains and a coarse problem of a few thousand DOF becomes necessary to keep iteration counts below 30.

The architecture is not a single choice but a coordinated set. Explicit integration eliminates the linear solve but demands tiny steps; implicit integration enables large steps but requires scalable sparse solvers and robust nonlinear iteration; domain decomposition determines whether those solvers can use thousands of cores effectively. Getting the combination right — and validating it with a scaling study on the actual hardware — is what separates a simulation that finishes overnight from one that never converges.