High-order 2D Hydrodynamics

This example is a Tungsten implementation of the Laghos algorithm for high-order Lagrangian compressible hydrodynamics in two dimensions, targeting Cerebras wafer-scale hardware. The computational mesh moves with the material, coupling its deformation to momentum and internal energy.

Governing equations

The inviscid Euler equations in the Lagrangian framework are

$$\begin{aligned} \frac{1}{\rho}\frac{d\rho}{dt} &= -\nabla\cdot v, \\ \rho\frac{dv}{dt} &= \nabla\cdot\sigma, \\ \rho\frac{de}{dt} &= \sigma:\nabla v, \\ \frac{dx}{dt} &= v, \\ p &= \operatorname{EOS}(\rho,e), \\ \sigma &= -pI. \end{aligned}$$

Here, $t$ is time, $x$ is the current position of a material point, $v$ is its velocity, $\rho$ is density, $e$ is internal energy per unit mass, $p$ is pressure, and $\sigma$ is the stress tensor. In 2D, $x$ and $v$ are two-component vectors and $I$ is the $2\times2$ identity matrix. The colon denotes the matrix inner product.

A material point identifies the same infinitesimal fluid parcel throughout its motion. Its material coordinates $X$ are fixed labels in the undeformed reference configuration, whereas $x=x(X,t)$ is the point's current position in the deformed configuration. Spatial derivatives $\nabla=\nabla_x$ and $\nabla\cdot$ act with respect to the current (deformed) spatial coordinates $x$, not $X$. The material derivative $d/dt=\partial/\partial t+v\cdot\nabla_x$ follows that material point, holding $X$ fixed; here $\partial/\partial t$ is taken at fixed $x$.

The example uses the ideal-gas equation of state $p=(\gamma-1)\rho e$, where $\gamma$ is the adiabatic index.

Spatial FEM semi-discretization

The moving quadrilateral mesh supports two tensor-product finite element spaces: a continuous kinematic space of degree $P+1$ in each coordinate for position and velocity, and a discontinuous thermodynamic space of degree $P$ for internal energy. The spatial discretization produces

$$\begin{aligned} M_V\frac{d\mathbf{v}}{dt} &= -F\mathbf{1}, \\ M_E\frac{d\mathbf{e}}{dt} &= F^T\mathbf{v}, \\ \frac{d\mathbf{x}}{dt} &= \mathbf{v}. \end{aligned}$$

The bold variables are finite element coefficient vectors, and $\mathbf{1}$ represents the constant 1 thermodynamic function. $M_V$ is the kinematic mass operator, $M_E$ is the element-block-diagonal thermodynamic mass operator, and the force operator $F$ couples the spaces through pressure forces and their work. The mass operators are constant in material coordinates; $F$ changes with geometry and the thermodynamic state. Density is recovered from mass conservation, $\rho\det J=\rho_0\det J_0$, where $J$ is the reference-cell-to-current-cell Jacobian and subscript $0$ denotes the initial state.

Time integration

The state advances with a prescribed timestep using the second-order, two-stage Runge-Kutta method with averaging, RK2+Averaging. Midpoint and full-step updates use averaged velocities to advance position and internal energy, with force data refreshed for the corresponding stage states.

Kernel composition

The implementation provides the following composable kernels and supporting operations:

• Kinematic mass action: evaluates $M_V\mathbf{u}$ for a kinematic vector $\mathbf{u}$.
• Force action: evaluates $F\mathbf{1}$ for the momentum equation.
• Force-transpose action: evaluates $F^T\mathbf{u}$ for the internal-energy equation.
• Thermodynamic inverse mass action: applies $M_E^{-1}$ using precomputed cell-local dense inverse blocks.
• Force-data update: recomputes geometric Jacobians, density, pressure, and the quadrature-point data defining $F$ for the current stage.
• Global DOT product: sums local contributions, counting each shared degree of freedom once. An all-reduce collects row contributions toward the east, reduces them southward to the south-east corner of the processor layout, and broadcasts the result back to all workers.
• Supporting operations: diagonal preconditioning, vector updates, shared-DOF gather-scatter through neighbor exchanges, and essential boundary-condition enforcement.

The mass and force actions use partial-assembly matrix-free methods: geometry and physics are evaluated at quadrature points, while sum factorization expresses 2D basis evaluations as sequences of 1D contractions. Gather-scatter assembles contributions at shared kinematic degrees of freedom without storing global sparse matrices.

The kinematic mass action, global DOT product, diagonal preconditioning, vector updates, and boundary enforcement are combined into a preconditioned conjugate gradient (CG) solver. It approximates the inverse action $M_V^{-1}\mathbf{b}$ needed for velocity updates, without explicitly assembling or inverting the global mass matrix.

Finally, the explicit RK2+Averaging loop combines CG with the force-data update, force action, force-transpose action, thermodynamic inverse mass action, boundary enforcement, and state-vector updates. Together these components advance position, velocity, and internal energy through the midpoint and full-step stages.