Source-linked AI summary

A Two-dimensional HLLC Riemann Solver for Conservation Laws : Application to Euler and MHD Flows

Dinshaw S. Balsara

arXiv:1110.0750v1physics.comp-ph

TL;DR

The paper addresses multidimensional upwinding with a genuinely two-dimensional HLLC Riemann solver for four states meeting at a mesh edge. It introduces a stabilized HLL strongly-interacting state and contact-discontinuity substructure, then uses the solver in a second-order three-dimensional Godunov scheme with routinely used CFL numbers of ~0.7 in two dimensions and about 0.6 in three dimensions.

  • Problem

    Existing multidimensional Riemann-solver constructions can be elaborate or difficult to extend, motivating a general solver for multidimensional fluxes and conservation laws.

  • Method

    The method builds a two-dimensional HLLC solver from an improved HLL solver by modeling the strong-interaction region and inserting a contact discontinuity whose orientation is supplied from the density gradient.

  • Results

    The resulting second-order Godunov scheme is general enough for three-dimensional simulations, with CFL numbers of ~0.7 routinely used in two dimensions and about 0.6 empirically found in three dimensions.

  • Takeaways & Limitations

    The multidimensional HLLC technology provides a cost-competitive basis for higher-order Godunov schemes for Euler and MHD simulations.

  • Takeaways & Limitations

    The solver requires a supplied normal direction for the contact discontinuity, inferred from the density gradient of the four zones meeting at an edge.

Abstract

from arXiv · show

In this paper we present a genuinely two-dimensional HLLC Riemann solver. On logically rectangular meshes, it accepts four input states that come together at an edge and outputs the multi-dimensionally upwinded fluxes in both directions. This work builds on, and improves, our prior work on two-dimensional HLL Riemann solvers. The HLL Riemann solver presented here achieves its stabilization by introducing a constant state in the region of strong interaction, where four one-dimensional Riemann problems interact vigorously with one another. A robust version of the HLL Riemann solver is presented here along with a strategy for introducing sub-structure in the strongly-interacting state. Introducing sub-structure turns the two-dimensional HLL Riemann solver into a two-dimensional HLLC Riemann solver. The sub-structure that we introduce represents a contact discontinuity which can be oriented in any direction relative to the mesh. The Riemann solver presented here is general and can work with any system of conservation laws. We also present a second order accurate Godunov scheme that works in three dimensions and is entirely based on the present multidimensional HLLC Riemann solver technology. The methods presented are cost-competitive with traditional higher order Godunov schemes.

II) Derivation of the Two-Dimensional HLLC Riemann Solver

Section II derives the two-dimensional HLL and HLLC Riemann solvers, presenting the HLL construction before extending it with HLLC structure in the Euler-system setting.

  • II.a) Derivation of an Improved HLL Riemann Solver in Two Dimensions: Sub-section II.a derives a robust two-dimensional HLL Riemann solver.
  • II.b) Derivation of the HLLC Riemann Solver in Two Dimensions: Sub-section II.b extends the two-dimensional HLL solver to include HLLC structure.
  • The derivation is presented within the context of the Euler system.

II.a) Derivation of an Improved HLL Riemann Solver in Two Dimensions

The improved two-dimensional HLL solver models four interacting one-dimensional Riemann problems with a bounded strongly-interacting state and computes multidimensional fluxes from this wave structure.

  • The solver accepts four states meeting at a mesh edge and formulates x- and y-directional fluxes on a logically rectangular mesh.For Euler equations, the conserved variables are density, three momentum components, and energy density.
  • Four one-dimensional Riemann problems interact in a central strongly-interacting state bounded by maximal wave speeds.The construction includes two x-directional and two y-directional Riemann problems around the edge.
  • The HLL model replaces unresolved sub-structure in the strong-interaction region with a single constant state whose extent follows the adopted wave model.
  • Wave-speed bounds are obtained by taking maxima and minima of directional signal speeds, then combining face-specific bounds to enclose the strongly-interacting state.
  • Resolved states and fluxes from directional HLL problems are used to construct transverse numerical fluxes in the other direction.
  • The improved integration includes contributions from resolved one-dimensional states, providing a stabilizing contribution relative to the earlier construction.The paper reports robust performance near ultra-strong shocks when using the lower-dissipation equations and their HLLC extensions.

II.b) Derivation of the HLLC Riemann Solver in Two Dimensions

The two-dimensional HLLC solver extends the HLL construction by subdividing the strongly interacting region with a movable contact discontinuity, producing multidimensionally upwinded fluxes. Its formulation uses directional contact speeds, resolved states, and geometric space-time integration, with additional stabilization for difficult regimes.

  • Contact-discontinuity structure: The HLLC solver subdivides the strongly interacting state into two states separated by a contact discontinuity that can move in any x-y direction.The contact orientation is determined from the larger-scale density gradient.
  • Flux construction: Directional HLLC Riemann solvers are constructed on the four faces, supplying x- and y-fluxes from resolved states and face-dependent maximal wave speeds.The lower and right faces are described explicitly, with analogous constructions for the remaining faces.
  • Geometric integration: The space-time control volume is an inverted rectangular pyramid whose triangular and trapezoidal surface integrals simplify evaluation of the multidimensional fluxes.Sub-portions containing each strongly interacting state provide the control volumes used in the integration.
  • Contact-discontinuity structure: The contact speeds are obtained from velocity components of the two-dimensional HLL strongly interacting state and applied consistently to the resolved states.The x- and y-contact speeds are ratios of the corresponding momentum components to density.
  • Flux construction: The numerical fluxes are assembled at computational-mesh edges from the strongly interacting states and the directional face fluxes.The stated implementation equations yield multidimensionally upwinded edge fluxes in both directions.

IV.b) Construction of a Second Order Accurate Predictor-Corrector Scheme Using the Multidimensional Riemann Solver

The scheme uses the two-dimensional HLLC Riemann solver in both predictor and corrector steps of a second-order predictor-corrector method, with edge-aligned sweeps for three-dimensional meshes. This raises stable CFL limits while retaining competitive computational cost relative to conventional Godunov schemes.

  • Scheme design: The predictor-corrector scheme uses the two-dimensional HLLC Riemann solver exclusively in both time-stepping stages.Using it in the corrector restores isotropy to the overall simulation.
  • Stability: A maximal CFL number of about 0.6 was found empirically for three-dimensional simulations, while two-dimensional problems routinely used approximately 0.7.Some two-dimensional problems ran at larger CFL numbers.
  • Three-dimensional implementation: Three-dimensional implementation processes z-, y-, and x-edge sweeps, obtaining four states at each edge and returning fluxes in two directions.Axis permutations are used for y- and x-sweeps before restoring the flux components.
  • Performance: 90,400 three-dimensional zones per second were updated with the multidimensional solver, versus 212,600 for a conventional Euler scheme.The multidimensional scheme operated at a higher maximal CFL number, about 0.6 versus about 0.3.
  • Performance: 56,700 three-dimensional MHD zones per second were updated with the multidimensional solver, and the method was judged cost-competitive with conventional second-order Godunov schemes.The comparison used a conventional MHD scheme updating 121,500 zones per second.

V.a) Hydrodynamic Tests

The hydrodynamic tests are designed to demonstrate advantages of introducing multidimensional Riemann solvers into flow-solver technology.

  • Test objectives: The paper presents several hydrodynamical test problems to show advantages of multidimensional Riemann solvers in current flow-solver technology.The tests are presented as demonstrations of the multidimensional approach.

V.a.1) Spherical Blast Wave in Three Dimensions

The three-dimensional spherical blast-wave test compares a conventional Godunov scheme with the multidimensional scheme on a resolution-limited mesh. The multidimensional solver produces a more isotropic density distribution and can reveal instability on a smaller mesh.

  • Test setup: A 64^3-zone mesh was used to test a spherical blast wave under numerical-resolution starvation.The domain was [-0.5, 0.5]^3 with continuative boundaries.
  • Results: The multidimensional scheme produced a more isotropic density distribution than the conventional second-order Godunov scheme.The comparison used CFL 0.6 for the multidimensional scheme and CFL 0.3 for the conventional scheme.
  • Results: The multidimensional Riemann solver revealed the instability on a much smaller mesh while using the conventional predictor step.The resulting solution was reported as as good as the higher-CFL result.

V.a.2) Multidimensional Riemann Problems in Two Dimensions

Two-dimensional multidimensional Riemann problems test shock, slip-line, and contact interactions at high resolution and CFL 0.9. The solver captures strong-interaction structures and reveals instability on relatively small meshes.

  • Problem setup: The tests use multidimensional Riemann problems whose initial conditions were catalogued from earlier studies.The problems are used to calibrate numerical schemes and assess multidimensional behavior.
  • First Riemann problem: The first problem starts from two weak shocks and two slip lines and evolves into Mach reflections and contact surfaces.A pronounced density valley moves toward the intersection of the four shocks.
  • First Riemann problem: The first problem was run at CFL 0.9 on a 1000×1000 mesh and produced a solution dominated by strong interaction.The simulation used an r=3 WENO reconstruction with the multidimensional Riemann solver.
  • First Riemann problem: The multidimensional solver showed Kelvin-Helmholtz instability emerging on a much smaller mesh when true multidimensionality and sub-structure were included.The instability is contrasted with its emergence on a larger mesh using a conventional scheme.
  • Second Riemann problem: A second problem produces a double Mach reflection and a shock propagating at 45° to the mesh.It was run at CFL 0.9 on a 1000×1000 mesh and the mushroom cap was captured crisply.

V.a.3) Double Mach Reflection Problem

The double Mach reflection test shows that the multidimensional scheme resolves the flow structures and Mach-stem instabilities at second-order resolution.

  • V.a.3) Double Mach Reflection Problem: The second-order 1920×480 simulation captures Kelvin–Helmholtz instability around the Mach stem.A 2400×600 simulation further improves the result, while a cited comparison required fourth order at 1920×480.

V.b) MHD Tests

The MHD tests use a divergence-free scheme combining multidimensionally upwinded electric fields with the presented multidimensional Riemann solver.

  • V.b) MHD Tests: The MHD tests use a divergence-free scheme with multidimensionally upwinded electric fields obtained through the multidimensional Riemann solver.The electric-field construction follows Balsara and Spicer's suggestion.

V.b.1) Long Term Decay of Alfven Waves in Two Dimensions

The torsional Alfvén-wave test evaluates long-term numerical dissipation through amplitude decay, while noting that nonlinear stabilization prevents quantitative dissipation measurement.

  • V.b.1) Long Term Decay of Alfven Waves in Two Dimensions: The test measures numerical dissipation by tracking long-term decay of torsional Alfvén-wave amplitudes in a periodic two-dimensional domain.The z-velocity and z-magnetic-field amplitudes decay numerically, whereas the analytical solution preserves the wave.
  • V.b.1) Long Term Decay of Alfven Waves in Two Dimensions: Nonlinear stabilization prevents this test from quantitatively measuring numerical dissipation.The amplitude decay remains a qualitative reflection of dissipation in the scheme.
  • V.b.1) Long Term Decay of Alfven Waves in Two Dimensions: Figure 9 compares z-velocity and z-magnetic-field decay for multidimensional and one-dimensional Riemann-solver schemes with different limiters.The plotted curves use the multidimensional solver with an r=3 WENO or MC limiter and compare against one-dimensional HLLC results.

V.b.2) Field Loop Advection in Two Dimensions

The field-loop advection test evaluates magnetic-loop preservation after one orbit and finds a more isotropic, less diffusive result with the multidimensional HLLC scheme.

  • V.b.2) Field Loop Advection in Two Dimensions: The test advects a low-magnetic-pressure field loop diagonally across a periodic 128x64 domain.The loop traverses a domain spanning [-1,1]×[-0.5,0.5].
  • V.b.2) Field Loop Advection in Two Dimensions: After one orbit, the multidimensional HLLC scheme produces an extremely isotropic loop with minimal boundary diffusion and no magnetic-pressure oscillations.The conventional second-order Godunov result is described as clearly inferior despite using a lower CFL number and a linearized Riemann solver.
  • V.b.2) Field Loop Advection in Two Dimensions: Extension of the multidimensional Riemann-solver technology beyond second order is left for future work.The paper states that this task will be undertaken later.
  • V.b.2) Field Loop Advection in Two Dimensions: Using the multidimensional solver only in the corrector step retains the test's solution quality with a two-dimensional CFL limit of ~0.45.The predictor step can remain conventional while producing results as good as the full multidimensional-solver result.

V.b.3) Three Dimensional MHD Blast Wave with very low Plasma Beta

The multidimensional HLLC scheme is applied to a stringent three-dimensional, very-low-plasma-beta MHD blast problem. It maintains positive pressure while supporting large time steps and accurate multidirectional magnetic-field propagation.

  • Three-dimensional MHD blast wave: CFL number 0.6 was used to evolve the blast to time 0.014 with the multidimensional Riemann-solver-based scheme.
  • Three-dimensional MHD blast wave: The scheme remains robust in a three-dimensional MHD blast with very low plasma-beta, where conventional second-order Godunov schemes can produce negative pressure.The pressure is computed by subtracting kinetic and magnetic energies from total energy density, making positivity challenging when velocities and magnetic fields are large.
  • Three-dimensional MHD blast wave: The method improves pressure positivity by accurately propagating magnetic fields in any required direction on the computational mesh.
  • Three-dimensional MHD blast wave: The solver accepts four states at an edge and returns two multidimensionally upwinded fluxes, incorporating a contact discontinuity with arbitrary orientation.
  • Three-dimensional MHD blast wave: The associated second-order method supports larger CFL numbers than conventional schemes and propagates strong discontinuities more isotropically.

Appendix:

The appendix specifies how space-time integration is selected and performed over triangular faces of the inverted-pyramid interaction region. A decision tree automates the choice among integration types.

  • Integration cases: Type I, II, and III cases partition triangular faces according to wave-speed positions and intersection conditions in space-time.
  • Integration cases: Type II integration covers the triangular region bounded by the origin and vertices involving the right and lower wave speeds.
  • Integration cases: Type 0 integration skips the lower triangular face, whereas Type III integrates over the entire lower face when the relevant wave-speed ordering permits it.
  • Automated selection: The decision tree assigns integration types to the left, right, down, and upper faces for each wave-speed configuration.
  • Automated selection: Figure 5 provides the wave-speed labels, state ordering, and supersonic-limit states needed to identify triangular-face integration limits.

Appendix B

This appendix develops the space-time geometry and consistency relations underlying the two-dimensional HLLC solver. The construction represents interacting states and an oriented contact discontinuity within an inverted rectangular pyramid.

  • State consistency: The rectangular area containing the HLLC state is related to the HLL state through a consistency condition involving the sub-areas and state vectors.
  • State consistency: Once the first HLLC interacting state is obtained, the consistency relation determines the second state when needed.
  • Contact geometry: The contact discontinuity is characterized by its motion in the x and y directions through the corresponding contact speeds.
  • Space-time construction: Four neighboring zones meet at an edge, and waves propagate outward along the x and y axes to bound the affected rectangular domain.
  • Space-time construction: The control-volume construction integrates the two-dimensional conservation law over rectangular space-time prisms associated with the edge interaction.
  • Contact geometry: The HLLC strongly-interacting states form an inverted rectangular pyramid in space-time, with the contact discontinuity dividing its base projection.

RD U

The appendix documents figures and diagnostics for hydrodynamic and MHD test problems. The plots compare density, wave decay, and magnetic-field structure across multidimensional and conventional schemes.

  • Figure guides: Figure 5 organizes wave-speed labels and states between waves, supporting identification of integration limits and supersonic states.
  • Hydrodynamic tests: Figures 6 and 7 show density contours for hydrodynamic blast and two-dimensional Riemann problems, including comparisons along grid-aligned and 45° directions.
  • Hydrodynamic tests: Figures 8c and 8d magnify the region around the Mach stem from the two preceding computations.
  • MHD tests: Figure 9 compares Alfven-wave decay using different limiters and multidimensional versus one-dimensional Riemann-solver schemes.
  • MHD tests: Figure 10 compares magnetic-field magnitude after one complete field-loop orbit for conventional second-order Godunov methods and the multidimensional HLLC scheme.
  • MHD tests: Figures 11a–11d display midplane contours of density, pressure, velocity magnitude, and magnetic-field magnitude for the MHD blast.
Loading 1110.0750v1…