About Static Pressure Condensation

Static pressure condensation is the single most worthwhile structural optimization in Sentinel's forward solve, and it is unusual in that it makes the answer more accurate while making it cheaper. This page is about why both of those are true, and why they are true of this discretization in particular rather than of mixed formulations in general.

For the derivation see Static Pressure Condensation in the Mathematical Reference; for defaults, benchmarks, and how to turn it off see How to Tune Performance.

Where the pressure unknowns come from

Soft tissue is nearly incompressible: its bulk modulus exceeds its shear modulus by roughly six orders of magnitude. Discretizing that directly causes volumetric locking — the elements become artificially stiff because they cannot satisfy the near-constant-volume constraint with displacement degrees of freedom alone.

The standard remedy is a mixed formulation, which introduces pressure as an independent unknown. Sentinel does this, adding four element-local pressure degrees of freedom per element, and the element stiffness acquires a saddle-point block structure with a displacement block, a pressure block, and the coupling between them.

That solves the locking problem and creates a different one. Saddle-point systems are harder to solve than positive-definite ones, and this particular saddle-point system is badly scaled — which is the subject of the next two sections.

Why elimination is exact here

Eliminating unknowns from a linear system by forming a Schur complement is ordinary algebra, available in any system. What varies is the price. Usually the Schur complement is dense, or much denser than what you started with, and the elimination costs more than it saves.

Two features of this discretization make it nearly free.

The pressure basis is discontinuous and element-local. A pressure degree of freedom belongs to exactly one element and is shared with no neighbour. That makes the global pressure–pressure block block-diagonal, with one independent 4×4 block per element, so the elimination is a 4×4 solve inside each element rather than a global operation. It also means the element-level Schur complement summed over elements reproduces the global Schur complement exactly — there is no approximation being made anywhere.

No new sparsity appears. The coupling that gets eliminated is confined to the displacement degrees of freedom of a single element, and those already couple to each other. So the condensed matrix has the same sparsity pattern as the pure-displacement stiffness. The system shrinks from $3n_{\text{nodes}}+4n_e$ unknowns to $3n_{\text{nodes}}$ — about 12% fewer on the MGH brain mesh — with no added fill in the factorization.

A smaller matrix with the same sparsity and no approximation is not a trade-off. It is why condensation is on by default rather than offered as an option.

The conditioning argument, which matters more than the size

The 12% size reduction is the less interesting half.

The penalty block scales as $1/{\kappa }^{\ast }$ with a bulk modulus around 2 GPa, while the displacement block scales as ${\mu }^{\ast }$, on the order of kilopascals. The two blocks therefore sit about twelve orders of magnitude apart on the diagonal, and applying boundary conditions widens the spread further. A direct factorization of that system in double precision spends several of its sixteen significant digits coping with the scale disparity.

Condensation does not rescale the problem. It removes the offending block entirely, leaving a displacement stiffness that is well conditioned on its own terms. The source of the ill-conditioning is gone rather than compensated for.

This is why the accuracy result comes out the way it does. Against the raw double-precision saddle-point solve, the condensed reconstruction differs by about ${10}^{-2}$ in relative L2 — and the natural reading, that condensation introduced the error, is backwards. An independent well-conditioned path (single-precision factorization with Jacobi equilibration and double-precision iterative refinement) lands on the same answer as condensation, agreeing with it to round-off, and deviates from the raw saddle solve by the same amount. Two independent methods agreeing with each other and disagreeing with a third is the usual signature of the third being wrong. The raw saddle solve is the inaccurate one; it lost the digits its conditioning predicted. The measured comparison is tabulated in How to Tune Performance.

It is worth being explicit about the general lesson, because it recurs: a formulation that is mathematically equivalent in exact arithmetic can be substantially better or worse in floating point, and "we solved the equations as written" is not the same as "we solved the problem accurately".

Where it does not apply

Condensation is not universal, and the boundaries are structural rather than cautious.

It needs an invertible pressure block, which the penalty-regularized incompressible models (1, 2, 5, and 6) provide. Model 7 uses a pure Lagrange multiplier formulation whose pressure block is zero — singular, and not eliminable this way.

It also assumes you do not need the pressure solution. The forward displacement, the objective, and the gradients with respect to shear modulus and density all depend only on the displacement. The gradients that genuinely read the pressure are the bulk-modulus gradient and the anisotropic models' pressure-coupling terms. So condensation applies cleanly to the standard viscoelastic MRE reconstruction — complex shear modulus, optionally density, bulk modulus held fixed — and that is exactly the case it is enabled for automatically.

If you reconstruct the bulk modulus, the pressure unknowns are load-bearing and condensation is not available. That is a real restriction, but it coincides with a reconstruction that is rarely wanted.

Related

Condensation is one of two responses to the saddle-point structure. The other is to solve the saddle system iteratively with a preconditioner strong enough to cope with it, which is the subject of Iterative Solvers on the GPU — worth reading alongside this page, because it explains why the preconditioner has to be as strong as it does, and condensation is precisely the route that avoids needing one.

• Static Pressure Condensation — the derivation.

• How to Tune Performance — defaults, validation table, benchmarks.

• About the MRE Inverse Problem — why the mixed formulation is there at all.