A More Rigorous Test Problem For Viscous Hydrodynamics Codes

This paper advocates for a more rigorous test problem for viscous hydrodynamics codes by introducing a nonuniform-density Gaussian velocity shear scenario that exposes errors in the calculation of viscous stress tensors, fluxes, and source terms, while also providing a detailed exposition of the Navier-Stokes equations in various coordinate systems.

The Gist

Imagine you are a chef trying to perfect a recipe for a very delicate soup. To test if your kitchen (your computer code) is working correctly, you usually stir in a single drop of dye and watch it spread out evenly in a bowl of water. If the dye spreads smoothly, you think, "Great! My stirring technique is perfect."

But here's the catch: real soup isn't just water. It has chunks of vegetables, different densities, and ingredients that sink or float. If your kitchen only works for plain water but fails when you add a chunky carrot, your "perfect" recipe is actually flawed.

This paper, written by Alexander Dittmann and Geoffrey Ryan, is essentially a warning to scientists who build computer simulations of fluids (like gas in space or water in a pipe). They are saying: "Stop testing your code with plain water. Start testing it with a chunky soup."

Here is the breakdown of their argument using simple analogies:

1. The Old Test: The "Flat Water" Problem

For years, scientists tested their fluid simulation codes using a problem called the Gaussian Velocity Shear.

• The Analogy: Imagine a calm pond. You create a perfect, smooth ripple in the middle. Because the water is uniform (same density everywhere), the ripple just spreads out slowly and predictably, like butter melting on toast.
• The Flaw: If your computer code has a hidden bug in how it calculates the "thickness" (viscosity) of the fluid, this bug might stay hidden because the water is so uniform. The code can get lucky and still look right.

2. The New Test: The "Density Slope" Problem

The authors propose a much harder test. Instead of a flat pond, imagine a river that is thick and heavy at the bottom and thin and light at the top (like a gradient of honey).

• The Setup: You drop that same smooth ripple into this river.
• The Twist: Because the fluid is heavier in some places than others, the ripple doesn't just spread out; it drifts. It gets pushed sideways by the density difference, like a leaf floating down a stream that is deeper on one side.
• Why it matters: To get this drifting motion right, your computer code has to do complex math. It has to calculate how the "thickness" of the fluid interacts with the "weight" of the fluid. If your code has a bug in these calculations, the ripple will drift in the wrong direction or at the wrong speed.

3. The "Gotcha" Moment

The authors ran this new test on three different popular computer codes (Athena++, Disco, and an older version of Disco).

• The Result: Two of the codes handled the "chunky soup" perfectly. They calculated the drift correctly.
• The Failure: One version of the code (Disco v2016) failed miserably. It looked great when tested on the "flat water" (uniform density), but when they added the density slope, it completely missed the drift.
• The Lesson: If you only test with "flat water," you might think your code is a genius. But the moment you introduce real-world complexity (like the swirling gas around a black hole or a binary star system), your code might crash or give you wrong answers.

4. The Appendix: The "Instruction Manual"

The paper includes a massive appendix that looks like a scary wall of math.

• The Analogy: Think of this as the instruction manual for the kitchen. It lists exactly how to calculate the "stirring forces" (viscous stress) in every possible shape of pot you might use—square pots (Cartesian), round pots (Cylindrical), and ball-shaped pots (Spherical).
• Why it's there: Many scientists use these codes in weird shapes (like rings around stars). The authors wanted to make sure everyone has the exact right formulas so they don't accidentally mess up the math when they switch from a square box to a round ring.

The Bottom Line

The authors are advocating for a "Stress Test" for fluid simulations.

Just as a car manufacturer doesn't just test a car on a flat, empty highway to see if the brakes work, they test it on steep, winding mountain roads. Similarly, scientists shouldn't just test their fluid codes on uniform, boring fluids. They need to test them on fluids with density gradients (where the fluid gets heavier or lighter across space).

If your code can handle the "drifting ripple" in the "chunky soup," then you can trust it to simulate the complex, messy, and beautiful physics of the universe. If it can't, you need to fix your code before you try to simulate a black hole or a galaxy.

Technical Summary

1. Problem Statement

Numerical codes designed to solve the Navier-Stokes equations for viscous hydrodynamics are frequently tested using a standard "Gaussian velocity shear" problem. In this classic test, a Gaussian velocity bump evolves in a fluid with uniform background density.

• The Limitation: The authors argue that this standard test is insufficient because it allows errors in the computation of viscous evolutionary terms to go undetected. Specifically, when density is constant, the coupling between density gradients and the viscous stress tensor is trivialized.
• The Goal: To propose and validate a more stringent test problem that forces codes to correctly calculate the full viscous stress tensor, fluxes, and source terms in the presence of non-uniform density.

2. Methodology

The authors derive an analytical solution for a non-uniform density Gaussian velocity shear and use it to benchmark three different numerical implementations (Athena++ and two versions of the Disco code).

A. Analytical Derivation

The test setup involves a one-dimensional shear flow in the $y$-direction, uniform in $y$ and $z$, with a specific density profile:

• Density Profile: $\rho(x) = \rho_0 \exp(-\kappa x)$, where $\kappa$ controls the density gradient.
• Assumptions:
  • Pressure is related to density via a spatially dependent sound speed ($P = c_s^2(x)\rho$).
  • The fluid is static in the $z$-direction.
  • Bulk viscosity is zero.
  • Isothermal equation of state is assumed to prevent viscous heating from complicating the evolution (fixing $c_s \ll v$).
• The Evolution Equation: Under these conditions, the evolution of the transverse velocity $v_y$ reduces to a modified diffusion-advection equation:
  $$\partial_t v_y + \nu\kappa \partial_x v_y - \nu \partial^2_x v_y = 0$$
  where $\nu$ is the kinematic viscosity.
• The Analytical Solution: For $t > 0$, the solution is a Gaussian profile that spreads and drifts:
  $$v_y(x, t) = \frac{A}{\sqrt{4\pi\nu t}} \exp\left( -\frac{(x - \kappa\nu t)^2}{4\nu t} \right)$$
  • Key Feature: Unlike the uniform density case ($\kappa=0$), the center of the velocity profile shifts over time with velocity $v_{drift} = \nu\kappa$ along the density gradient. This drift is a direct consequence of the correct calculation of viscous fluxes in a non-uniform medium.

B. Numerical Implementation

The authors tested this problem using:

1. Athena++: Implemented in both Cartesian (uniform mesh) and Cylindrical Polar (log-uniform radial mesh) coordinates.
2. Disco: Implemented in Cylindrical Polar coordinates (linear radial spacing). They tested two viscosity schemes:
   • Disco (v2016): An approximate scheme assuming globally uniform dynamic viscosity.
   • Disco (2021): A complete shear viscosity implementation.

• Parameters: Simulations used $\kappa=0$ (standard test) and $\kappa=15$ (rigorous test), with initial Gaussian width $\sigma_0=0.1$ and amplitude $v_0=10^{-3}$.

3. Key Contributions

1. Proposal of a New Benchmark: The paper introduces the "Non-uniform Density Gaussian Shear" test as a superior standard for validating viscous hydrodynamics codes.
2. Analytical Solution for Drifting Shear: They provide the closed-form analytical solution (Eq. 11) which demonstrates that density gradients induce a drift in the velocity profile, a phenomenon absent in constant-density tests.
3. Appendix of Curvilinear Forms: The authors provide a comprehensive derivation of the compressible Navier-Stokes equations in curvilinear coordinates (Cartesian, Cylindrical, and Spherical). This includes explicit expressions for the fluxes and source terms required for finite-volume numerical schemes, addressing a common pain point for code developers working in non-Cartesian geometries.

4. Results

The authors compared the numerical results against the analytical solution at $t=0.02$:

• Uniform Density ($\kappa=0$): All codes (Athena++, Disco v2016, Disco 2021) converged at second order and agreed well with the analytical solution. This confirms that standard tests are insufficient for detecting specific viscosity implementation errors.
• Non-Uniform Density ($\kappa=15$):
  • Athena++ and Disco (2021): Both codes accurately reproduced the analytical solution, correctly capturing the drift of the velocity profile from $x=0.5$ to $x=0.8$.
  • Disco (v2016): This code failed significantly. It did not converge to the analytical solution and produced large errors.
• Diagnosis: The failure of Disco (v2016) was traced to its assumption that dynamic viscosity is globally uniform. In a non-uniform density field, this assumption leads to incorrect calculation of the viscous stress tensor and fluxes. The new test successfully exposed this flaw, which the $\kappa=0$ test had missed.

5. Significance

• Code Validation: The paper demonstrates that relying solely on uniform-density tests can lead to false confidence in a code's accuracy. The proposed test is "more discriminating" and essential for verifying that a code correctly handles the coupling between density gradients and viscous stresses.
• Astrophysical Relevance: The test is particularly relevant for simulating astrophysical systems with strong density gradients, such as circumbinary disks and gap-opening planets, where the flawed viscosity implementation in Disco (v2016) had previously led to erroneous results.
• Resource for Developers: The detailed Appendix serves as a critical reference for developers implementing viscous hydrodynamics in complex coordinate systems (polar/spherical), ensuring that the geometric source terms and fluxes are derived correctly.

In conclusion, Dittmann and Ryan advocate for the immediate adoption of the non-uniform density shear test to ensure the robustness of viscous hydrodynamics codes, particularly for applications involving strong density gradients.