Dimensional analysis of the Navier-Stokes equations

This demo demonstrates the automated dimensional analysis in dolfiny.units on a transient incompressible Navier-Stokes lid-driven cavity problem. The demo does not solve the nonlinear cavity problem. Instead, it assembles the Jacobian of the normalized residual at the zero state and studies how its matrix condition number varies with the Euler number.

We start from dimensional reference quantities, check the weak form, and normalize the residual by the convective scale. The resulting Reynolds, Euler, Froude, and Strouhal dimensionless numbers emerge directly from the UFL form.

In particular this demo emphasises

1. automated Buckingham Pi analysis and dimensional checks for a UFL model,

2. extraction of Reynolds, Euler, Froude, and Strouhal numbers by normalising the weak form, and

3. how the conditioning of the normalised Navier-Stokes Jacobian changes under pressure scaling.

Model¶

On the unit square $\Omega = (0, 1)^2$ we consider velocity $\mathbf{v}$ and pressure $p$ such that the incompressible Navier-Stokes equations in the strong form read

Here $\Omega$ is the spatial domain, $\partial_t$ denotes the time derivative, $\nabla$ denotes the spatial gradient or divergence as appropriate, and $\mathbf{I}$ is the identity tensor. Furthermore, $D(\mathbf{v}) = \text{sym}(\nabla \mathbf{v})$ is the strain-rate tensor, $\rho$ is the density, $\nu$ is the kinematic viscosity, $g_\text{ref}$ is a gravity scale, and $\mathbf{b} = (0, -1)$ is a dimensionless body-force direction.

Reference quantities and nondimensional groups¶

The dimensional model is described by the reference scales. The subscript ref marks prescribed dimensional scales used to normalize the equations.

Each scale is created as a Quantity(mesh, scale, unit, symbol). Quantity is a dolfinx.fem.Constant with extra metadata: besides the scalar value used in assembly, it stores the unit, SymPy-based symbol, and derived dimension used for Buckingham Pi analysis and dimensional checks.

Buckingham Pi analysis below reports the basis (so-called Pi groups) of the nullspace of the dimension matrix. For this problem, there are four Pi groups:

A call to dolfiny.units.buckingham_pi_analysis with a list of Quantity objects reports the overview of all dimensional quantities in the model, the dimension matrix, and the Pi groups. Please note, that Pi groups are found using a Matrix.nullspace() call in SymPy. For matrices over the field of rational numbers, the nullspace is usually represented with small integer coefficients, which is favorable for interpretability.

Weak form and dimensional checks¶

We write the residual in term-by-term form so that each contribution can be transformed and factorized independently. We want to find $(\mathbf{v}, p)$ such that the residual $F = 0$ for all test functions $(\delta \mathbf{v}, \delta p)$. Here $\mathbf{v}_0$ denotes the velocity at the previous time level, and $N_\Delta$ is the number of time steps per reference time so that $\Delta t = t_\text{ref} / N_\Delta$:

with

Here $A : B$ denotes the Frobenius product of tensors and $\text{d}x$ denotes integration over $\Omega$. The last term is identically zero. It is kept only so that DOLFINx allocates a pressure-pressure sparsity block, allowing the pressure pinning condition to modify a diagonal entry there. The factor $1 / (p_\text{ref} t_\text{ref})$ is only a dimensional placeholder inside this zero term.

The dictionary mapping below connects the dimensionless coordinates, unknowns, and test functions to their dimensional counterparts: coordinates scale with $l_\text{ref}$, velocity with $v_\text{ref}$, and pressure with $p_\text{ref}$.

The mapping applies the following transformations:

After the transformation, we proceed with the steps: factorization and normalization. Factorization performs a homogeneous factorization of the weak form and extracts the homogeneous factors. Normalization then divides each term by a chosen reference term, which is the convection term in this case. The resulting dimensionless form contains only the Pi groups as coefficients.

Lid-driven cavity boundary conditions¶

The normalized residual uses standard lid-driven cavity boundary data: no-slip on the left, right, and bottom boundaries, unit tangential velocity on the lid, and one pressure degree of freedom fixed at the corner point $(0, 0)$ to remove the constant nullspace. This keeps the conditioning study focused on the reference scaling rather than on an undetermined pressure level. The gravity term is retained only so that the dimensional analysis also exposes the corresponding Froude dimensionless number.

Matrix conditioning under pressure scaling¶

We assemble the Jacobian matrix $A$ of the mixed residual at the zero state $(\mathbf{v}, p) = (\mathbf{0}, 0)$ for a sequence of pressure scales. The Jacobian is assembled into a PETSc MATAIJ matrix, converted to SciPy CSR format with dolfiny.la.petsc_to_scipy, and then densified only for the call to numpy.linalg.cond, which estimates the spectral condition number $\kappa(A) = \|A\|_2 \|A^{-1}\|_2$.

In the figure below, we observe that the condition number grows as the Euler number $\mathrm{Eu}$ decreases. This is expected since $\mathrm{Eu}$ is the coefficient of the pressure term in the normalized form, and the pressure term is responsible for the incompressibility constraint, and in turn the invertibility of the entire saddle-point Jacobian. In addition, for $\mathrm{Eu} \gtrsim 1.7$ the condition number is at its lowest value and does not vary. This is because the pressure coefficient has a sufficiently large norm to ensure good conditioning, and what determines the condition number in this regime is the lowest singular value of the velocity-velocity block, see Section 4.1 of Habera & Zilian (2026) for more details.