MMS: Verification of the Incompressible Navier-Stokes Equations

Good question. MMS (Method of Manufactured Solutions) is the ultimate tool for code verification. The key point is not "fake solution" but "a solution of your own choosing."

Complex equations like the Navier-Stokes equations generally do not have analytical solutions. But with MMS, by arbitrarily declaring any smooth function to be "the solution," you can rigorously check if the solver's implementation is mathematically correct.

That's the trick of MMS. For example, if we denote the differential operator of the Navier-Stokes equations as $\mathcal{L}$, then originally we look for $u$ that satisfies

$$ \mathcal{L}(u) = 0 $$

. But with MMS, you choose any $u_\text{mfg}$ you like and substitute it into the equation. Naturally, $\mathcal{L}(u_\text{mfg}) \neq 0$, so we add that residual as a source term $f$ on the right-hand side:

$$ \mathcal{L}(u_\text{mfg}) = f $$

In other words, we artificially create a problem where "in a world with this source term $f$, $u_\text{mfg}$ is the exact solution." We make the solver solve this problem and check if it matches $u_\text{mfg}$.

It does. Because we check the order of accuracy when systematically refining the mesh. If discretized with a 2nd-order scheme, the error should become $1/4$ when the mesh width $h$ is halved. If this breaks, there's a bug somewhere in the code.

For example, there was a case with an automotive exterior CFD code where the discretization of the convection term had a sign error. Standard benchmarks (like Lid-driven cavity) gave results that "looked about right," but MMS revealed the order of convergence dropped to 1st order, leading to the bug being found.

The incompressible NS equations consist of two parts: the momentum conservation equation and the continuity equation (mass conservation). In 2D, they are written as follows:

Theoretically, anything, but there are practical rules to follow:

1. It must satisfy the continuity equation: If $\nabla \cdot \mathbf{u}_\text{mfg} \neq 0$, you'd be giving contradictory input to an incompressible solver.
2. It must be sufficiently smooth: If not differentiable, you cannot verify high-order accuracy schemes.
3. It must "activate" all terms: For example, a constant manufactured solution would make both convection and viscous terms zero, rendering the test meaningless.

A typical example is the Taylor-Green vortex. This is actually an exact solution to the NS equations (holds without a source term), but it's perfect for explaining MMS.

Just actually compute the partial derivatives:

$$ \frac{\partial u_\text{mfg}}{\partial x} = \pi\sin(\pi x)\sin(\pi y)\,e^{-2\pi^2\nu t} $$ $$ \frac{\partial v_\text{mfg}}{\partial y} = -\pi\sin(\pi x)\sin(\pi y)\,e^{-2\pi^2\nu t} $$

They sum to zero. It perfectly satisfies the continuity equation. This manufactured solution is designed so that the structures of $u$ and $v$ are "sign inversion + sin/cos swap," ensuring the incompressibility condition is always met.

Exactly. Calculate each term in order. Let's show it for the $x$-direction momentum equation.

First, the unsteady term:

$$ \frac{\partial u_\text{mfg}}{\partial t} = 2\pi^2\nu\cos(\pi x)\sin(\pi y)\,e^{-2\pi^2\nu t} $$

Convection term ($u \cdot \partial u/\partial x + v \cdot \partial u/\partial y$):

$$ u\frac{\partial u}{\partial x} = -\pi\cos(\pi x)\sin(\pi y)\cdot\pi\sin(\pi x)\sin(\pi y)\cdot e^{-4\pi^2\nu t} $$ $$ v\frac{\partial u}{\partial y} = \pi\sin(\pi x)\cos(\pi y)\cdot(-\pi)\cos(\pi x)\cos(\pi y)\cdot e^{-4\pi^2\nu t} $$

Viscous term:

$$ \nu\nabla^2 u = \nu\cdot 2\pi^2\cos(\pi x)\sin(\pi y)\,e^{-2\pi^2\nu t} $$

The unsteady term and the viscous term cancel, and for the Taylor-Green vortex, the source term becomes $f_x = 0$. That is, this manufactured solution is also an exact solution of the NS equations.

It can. Even with $f=0$, the fact that "the exact solution is known" remains unchanged, so convergence order verification can be done without issue. However, in practice, more general manufactured solutions are often used. For example:

$$ u_\text{mfg} = \sin(ax)\sin(by) + c $$ $$ v_\text{mfg} = \cos(ax)\cos(by) + d $$

You can choose something like this that does not satisfy the continuity equation and adjust the coefficients $a, b$. In this case, the source term becomes non-zero, so you can also test the source term implementation.

Here we use the Taylor-Green vortex for educational purposes, but in practice, the golden rule is to use more complex manufactured solutions (combinations of polynomials + trigonometric functions) to activate all discretized terms.

Systematically refine the mesh and calculate the error $e$ on each mesh. The observed convergence order between two meshes is:

Let's look at an actual numerical example. When refining the mesh as $h$, $h/2$, $h/4$, $h/8$:

Let's also write the formula for the $L_2$ norm:

Since FVM deals with cell-averaged values, the source term is also integrated over the cell volume. The discretized momentum equation for cell $i$ is:

$$ \sum_f \phi_f u_f - \sum_f \nu (\nabla u)_f \cdot \mathbf{S}_f = -\frac{1}{\rho}\sum_f p_f \mathbf{S}_f + \int_{V_i} f \, dV $$

Here, $\phi_f = \mathbf{u}_f \cdot \mathbf{S}_f$ is the face flux, $\mathbf{S}_f$ is the face normal vector. The volume integral of the source term $\int_{V_i} f \, dV$ is approximated by the cell center value $f_i \cdot V_i$ (single-point approximation at cell center).

Good point. The cell center single-point approximation is 2nd-order accurate, so it's fine for verifying 2nd-order schemes. However, if you want to verify 3rd-order or higher accuracy, you need to integrate the source term with higher precision using Gauss quadrature. This is an easily overlooked point, so be careful.

There are two main methods:

• Method 1: Explicitly add the source term -- Add a user-defined source term function to the right-hand side of the momentum equation. In OpenFOAM, use fvOptions; in Fluent, use the UDF macro DEFINE_SOURCE.
• Method 2: Add as a body force -- Give the source term in the same form as the gravity term. Implementation is simple, but sometimes the source term for the pressure term needs to be handled separately.

In practice, Method 1 is common. Implement the source term function either hard-coded or as a user function taking coordinates and time as arguments.

It is desirable to calculate three types of error norms: $L_1$, $L_2$, $L_\infty$:

I wouldn't say it's not enough, but there are cases where the $L_\infty$ (maximum error) norm breaks down. For example, a bug where accuracy drops only near the boundary might be obscured by averaging in $L_2$, but $L_\infty$ would reveal it instantly. If all three norms show the correct convergence order, the code's reliability is quite high.

In practice, the golden rule is to use SymPy (Python's symbolic computation library) for automatic derivation. Keep manual calculations only for verification. Especially, sign errors frequently occur when expanding the nonlinear convection term $(u \cdot \nabla)u$, so I strongly recommend cross-checking with symbolic computation. Reference code: