Arbitrary-order DG solver for the 3D Euler equations

Finite volume methods are robust and well understood, but raising their order of accuracy means widening the stencil, which becomes awkward on unstructured meshes and worse in parallel. Discontinuous Galerkin keeps the stencil local: accuracy comes from the polynomial order inside each element rather than from reaching further into its neighbours.

Locality is the practical advantage. Each element communicates only across the faces it shares, whatever the order, which is the property that lets the solver scale across processors and makes a deforming mesh tractable rather than catastrophic.

The elements are triangles, mapped isoparametrically: the same shape functions that interpolate the solution also describe the element geometry, so curved and distorted elements are handled without a separate geometric treatment.

Each element computes the Jacobian of that mapping, its determinant, and the gradients of the shape functions in physical coordinates. The determinant also provides the area, and a near-zero value is trapped as a degenerate element rather than quietly producing nonsense downstream. Small guards like this are most of what separates a solver that runs from a solver that can be trusted.

Polynomial order is a parameter rather than a decision baked into the discretisation. The same code runs at low order where the flow is smooth and robustness matters, and at high order where resolving fine structure cheaply is worth more than robustness.

This matters specifically for the eigenmode work. The quantities of interest there are derived from the flow rather than read off it directly, and at low order the truncation error can be the same size as the signal being measured.

A swimming body deforms, so the mesh must deform with it. That breaks the assumption most solvers rest on, that geometry is fixed, and it introduces a conservation requirement of its own. In an arbitrary Lagrangian-Eulerian formulation the scheme has to account for the volume each face sweeps as it moves. Get it wrong and a uniform flow on a moving mesh stops being uniform: the solver manufactures disturbances out of its own grid motion.

Satisfying the geometric conservation law is unglamorous and load-bearing. It is the difference between a solver that produces plausible pictures and one that produces numbers worth publishing.

Fortran remains the right tool for dense numerical kernels with predictable memory access, and it sits directly on top of BLAS and LAPACK, which do the heavy lifting in the element-local operations. The choice is about the maturity of the numerical libraries and the compiler's willingness to vectorise, not about nostalgia.