PaperlayerView exact revision on arXiv

Source-linked AI summary

GAMERA: A three-dimensional finite-volume MHD solver for non-orthogonal curvilinear geometries

Binzheng Zhang, Kareem A. Sorathia, John G. Lyon, Viacheslav G. Merkin, Jeffrey S. Garretson, Michael Wiltberger

arXiv:1810.10861v2physics.comp-phphysics.space-ph

TL;DR

Plasma simulations require meshes that conform to system geometry, while the aging LFM code is difficult to adapt to modern architectures. GAMERA rewrites LFM to combine flexible curvilinear-grid MHD with upgraded reconstruction, metric calculations, and implementation, including a design aimed at exascale-era computing and simulations that can run near real time on a laptop.

  • Problem

    The aging LFM code is difficult to adapt to modern and future architectures, motivating a new code for geometrically flexible plasma simulations.

  • Method

    GAMERA rewrites LFM around upgraded high-heritage MHD numerics, flexible multidimensional solvers, high-order reconstruction and grid metrics, and modern computational implementation.

  • Results

    GAMERA provides a portable, general-purpose MHD code designed for arbitrary non-orthogonal curvilinear geometries and exascale-era architectures.

  • Takeaways & Limitations

    GAMERA is intended as a useful tool for applications where geometric flexibility is important, with magnetospheric simulations able to run near real time on a laptop.

  • Takeaways & Limitations

    The predecessor LFM code's age makes it difficult to adapt to modern and future architectures.

Abstract

from arXiv · show

Efficient simulation of plasmas in various contexts often involves the use of meshes that conform to the intrinsic geometry of the system under consideration. We present here a description of a new magnetohydrodynamic code, Gamera (Grid Agnostic MHD for Extended Research Applications), designed to combine geometric flexibility with high-order spatial reconstruction and constrained transport to maintain the divergence-free magnetic field. Gamera carries on the legacy of its predecessor, the LFM (Lyon-Fedder-Mobarry), a research code whose use in space physics has spanned three decades. At the time of its initial development the LFM code had a number of novel features: eighth-order centered spatial differencing, the Partial Donor Cell Method limiter for shock capturing, a non-orthogonal staggered mesh with constrained transport, and conservative averaging-reconstruction for axis singularities. A capability to handle multiple ion species was also added later. Gamera preserves the core numerical philosophy of LFM while also incorporating numerous algorithmic and computational improvements. The upgrades in the numerical schemes include accurate grid metric calculations using high-order Gaussian quadrature techniques, high-order upwind reconstruction, non-clipping options for interface values, and improved treatment of axis singularities. The improvements in the code implementation include the use of data structures and memory access patterns conducive to aligned vector operations and the implementation of hybrid parallelism, using MPI and OMP. Gamera is designed to be a portable and easy-to-use code that implements multi-dimensional MHD simulations in arbitrary non-orthogonal curvilinear geometries on modern supercomputer architectures.

1 Introduction

GAMERA is a rewrite of the widely used LFM code that preserves its high-heritage MHD numerics while preparing simulations for modern architectures and flexible curvilinear grids. It adds higher-order geometric and reconstruction methods, improved singularity treatment, and modern implementation features.

  • 1 Introduction: GAMERA modernizes implementation through software and computational changes intended to support efficient execution on current and future architectures.The paper describes implementation details including vectorization and parallelization, and presents five representative test simulations.
  • 1 Introduction: GAMERA rewrites LFM from scratch to retain and improve its high-heritage MHD numerics while preparing the code for the exascale era.LFM had been widely used in space physics and supported non-orthogonal curvilinear geometry, high-order reconstruction, shock limiting, constrained transport, and multiple ion species.
  • 1 Introduction: Its numerical upgrades include high-order Gaussian-quadrature grid metrics, seventh-order upwind reconstruction, non-clipping limiters, and higher-order axis-singularity treatment.These changes extend the legacy numerical approach while addressing geometric and reconstruction requirements.
  • 1 Introduction: The code supports multi-dimensional MHD simulations on general curvilinear computational grids adapted to specific problems.The design targets flexible two- or three-dimensional simulations in non-orthogonal geometries.
  • 1 Introduction: Earth or planetary magnetospheric simulations can run on a laptop close to real time, illustrating the code’s reported efficiency.The paper presents GAMERA as portable and usable for research applications, with algorithms and implementation details intended to support future extensions.

2 The Basic MHD Equations

GAMERA’s default solver uses normalized single-fluid ideal MHD equations in semi-conservative form with a plasma energy equation. The formulation represents mass, momentum, energy, pressure, magnetic field, and ideal electric-field relationships.

  • 2 The Basic MHD Equations: The default GAMERA solver advances single-fluid, normalized ideal MHD equations in semi-conservative form.The default equation set includes the plasma energy equation.
  • 2 The Basic MHD Equations: The formulation defines plasma density, bulk velocity, unit tensor, thermal pressure, magnetic field, and ideal electric field.The electric field is defined as E = −u × B through the ideal Ohm’s law.
  • 2 The Basic MHD Equations: Plasma energy is defined as the sum of kinetic and thermal energy, with the equations normalized according to the paper’s stated convention.The specific-heat ratio is also identified in the equation definitions.

Appendix A.

The appendix formulates ideal MHD using conserved and primitive variable sets, while noting that the plasma-energy formulation trades full conservation for advantages in selected low-β applications.

  • Variable formulation: The MHD system is expressed through conserved-density, fluid, magnetic, primitive, and volume-integrated variable vectors used for evolution and reconstruction.Primitive variables are mainly used in reconstruction, while volume-integrated conserved densities are the quantities evolved by the solver.
  • Energy formulation: The plasma energy equation is advantageous for problems involving planetary magnetospheres with strong background fields and cold ambient plasmas, including β < 10−6.
  • Energy formulation: Using plasma energy instead of total energy makes the equation set non-fully-conservative and can produce non-physical results when non-ideal processes occur.
  • Shock behavior: Despite the reduced conservation accuracy, the plasma-energy formulation preserves Rankine-Hugoniot relations to within numerical truncation error, independently of dissipative electric-field transport through shocks.The resulting shock behavior is correct to a good approximation, although its numerical error exceeds that of the total-energy formulation.

3 Numerical Schemes

GAMERA uses finite-volume discretization and staggered constrained transport to solve MHD equations on arbitrary non-orthogonal curvilinear grids. Its numerical schemes combine high-order reconstruction, geometry-aware metrics, and magnetic-flux evolution while preserving integral magnetic divergence.

  • Grid Discretization: GAMERA represents general hexahedral cells in a logically Cartesian computational space and evolves volume-averaged conserved quantities through surface fluxes.Cell geometry is handled through face areas, face normals, and cell volumes in the finite-volume formulation.
  • Grid Discretization: Cell-centered velocity and magnetic-field vectors use Cartesian components, avoiding geometry-related source terms without requiring coordinate orthogonality.Face-normal vectors and coordinate transforms are used when computing numerical fluxes through cell interfaces.
  • Grid Discretization: Grid metrics, including cell volumes, face-center locations, and face areas, are calculated with high-order Gaussian quadrature for curvilinear cells.The solver uses 12th-order quadrature for the cell volume and related geometric quantities.
  • Magnetic Fluxes and Fields: GAMERA evolves face-centered magnetic fluxes and edge-centered electric fields on a staggered non-orthogonal grid, conserving integral magnetic divergence to round-off error when initially divergence-free.Constrained transport updates fluxes through the integral Faraday law, and the staggered discretization preserves the divergence condition during computation.
  • Spatial Reconstruction: The solver applies high-order interface reconstruction, supports upwind schemes, and uses the Partial Donor Cell Method limiter to preserve monotonicity.The default centered reconstruction is eighth order, while the limiter controls whether high-order interface estimates require limiting.
  • Spatial Reconstruction: The seventh-order reconstruction has a dissipative eighth-order leading truncation-error term rather than a dispersive one, improving on the original LFM scheme.This changes the character of the leading error term from dispersive to dissipative.

3.4 Gas-Kinetic Flux Functions

Gamera computes fluid and magnetic-stress fluxes with Maxwellian-based gas-kinetic functions in local face-normal coordinates, then transforms momentum fluxes back to Cartesian coordinates. The formulation does not require orthogonal physical grids, but very low-beta plasmas require additional Alfvén-speed diffusion.

  • Coordinate treatment: Fluxes are evaluated in a local face-normal coordinate system, with momentum fluxes rotated back to Cartesian coordinates for conservative updates.This rotation works without requiring physical-grid orthogonality.
  • Flux formulation: Maxwellian-based gas-kinetic flux functions separately compute fluid terms ∆UF and magnetic-stress terms ∆UB.The separate treatment simplifies implementation of the semi-relativistic Boris correction.
  • Flux evaluation: The gas-kinetic scheme derives mass, momentum, and energy fluxes from velocity moments of left and right Maxwellian distributions across each interface.Identical left and right states recover the one-dimensional ideal gas-dynamics equations.
  • Limitation: For very low-beta plasmas, using thermal speed in the gas-hydro distribution is inadequate for Alfvén-wave-dominated dynamics, so additional numerical diffusion is required.The stated example regime is β < 10^-4.
  • Alfvén-speed correction: The solver adds Rusanov-type diffusion based on the average interface Alfvén speed when left and right states differ, while smooth regions receive no such added diffusion.This targets discontinuities without introducing the correction when interface states are equal or nearly equal.

10. Then the magnetic field strength in the µ-direction at the cell edge

The magnetic-field update reconstructs edge quantities in non-orthogonal curvilinear geometry, computes edge electric fields in edge-aligned orthogonal frames, and advances magnetic fluxes with Faraday’s law. Constrained transport preserves initially divergence-free fields to round-off error, while the implementation scales to production simulations at O(10k) CPU-cores.

  • Edge stabilization: The resistive electric-field term activates only near detected discontinuities; in smooth regions the edge field is a high-order approximation of −u × B.The additional term damps spurious oscillations associated with Alfvén waves near discontinuities.
  • Edge electric fields: Edge electric fields are computed by transforming velocity and magnetic-field quantities into an orthogonal frame aligned with the cell edge, then applying a Rusanov-like calculation.The edge-centered electric field is multiplied by edge length for the Faraday-law flux update.
  • Magnetic-flux evolution: Initially divergence-free magnetic fields remain divergence-free to round-off error during the Maxwell-solver evolution step.The volume-integrated ∇·B is unchanged by the magnetic-field update.
  • Implementation performance: Gamera supports production simulations at O(10k) CPU-cores and achieves approximately 500k zone-cycles per second per core with 85% of optimal OpenMP scaling.The implementation also shows near-linear scaling until approximately 163 cells per compute core.

4 Test Results

GAMERA’s tests assess high-order reconstruction, accuracy, and robustness across linear advection and multidimensional MHD problems on Cartesian, distorted, and non-orthogonal grids.

  • 4.1 Circular Advection: The 7th-order PDM scheme achieves comparable spatial gradients with 128 × 128 cells, while the 2nd-order TVD scheme requires 1024 × 1024 cells.Despite costing about 24% more per computation than the second-order method, the second-order run requires more than two orders of magnitude higher total computational cost.
  • 4.2 Field-Loop Advection: The field-loop tests retain the loop shape on highly distorted grids, while magnetic-energy retention decreases from approximately 98.5% to 97.7% as distortion increases.Lower reconstruction order broadens the magnetic-field-loop edge through increased numerical diffusion width.
  • 4.3 Non-linear Alfvén Wave: The default 7th-order reconstruction produces an average L1 error approximately ten times lower than 2nd-order reconstruction, and non-clipping reduces it by another factor of about three.The tests also find approximately second-order convergence for the smooth nonlinear Alfvén wave, independent of stencil size.
  • 4.3 Non-linear Alfvén Wave: Wavefronts, pressure distributions, shocks, and blast-wave structures remain similar across Cartesian and distorted non-orthogonal curvilinear grids.Minor differences occur near shock transitions, while smooth regions show very little difference and the blast-wave results show no noticeable grid-alignment effects.

5 Summary

GAMERA reinvents the LFM MHD kernel with upgraded numerical schemes and modern implementation choices for flexible, general-purpose simulations in non-orthogonal curvilinear grids.

  • GAMERA re-engineers the LFM MHD kernel with substantial numerical and software upgrades while retaining its established computational foundation.The paper presents the code as a reinvention of LFM rather than a wholly unrelated solver.
  • 12th-order Gaussian quadrature, high-order upwind reconstruction, an extremum-preserving PDM option, and background-field splitting improve the numerical toolkit.
  • Ring-Average treatment addresses spherical-axis singularities, supporting high-resolution simulations in general geometries with such singularities.
  • Modern Fortran, aligned memory blocks, contiguous access, OpenMP threading, and MPI distribution make the implementation more compatible with HPC architectures.
  • GAMERA targets multidimensional MHD simulations on non-orthogonal curvilinear grids and is intended for applications where geometric flexibility matters.
  • Its applications include global planetary-magnetosphere, inner-heliosphere and solar-wind, and local plasma-physics simulations.

C. The General Form of the PDM Operator

The PDM operator combines reconstruction order with a parameter A that controls numerical diffusion and limiter aggressiveness. Its implementation supports high-order interface values while recovering the diffusive donor-cell method when A = 0.

  • C. The General Form of the PDM Operator: For GAMERA, the default PDM settings are A = 4.0 and Courant number NCFL = 0.3 because larger A values provide little additional square-wave resolution.The supplied text reports little improvement for A > 4 and identifies these as the solver defaults.
  • C. The General Form of the PDM Operator: The PDM limiter can combine with arbitrarily high-order reconstruction to calculate left- and right-state interface values according to the velocity direction.The high-order interface value is used in the state calculation, with the sweep direction given by sign(v).
  • C. The General Form of the PDM Operator: When A = 0, the PDM operator becomes the first-order donor-cell method, identified as highly diffusive.This limiting case is stated explicitly for the general operator.
  • C. The General Form of the PDM Operator: A controls numerical diffusion: increasing A reduces diffusion but makes limiting more aggressive and lowers the allowable Courant number.The operator permits A > 0, while the monotonicity constraint links A to the Courant number.

D. The Choice of Reconstruction Order

GAMERA chooses seventh-order upwind reconstruction as its default because it reduces numerical diffusion with reasonable computational cost. The choice targets three-dimensional global-scale plasma simulations where both accuracy and resources matter.

  • D. The Choice of Reconstruction Order: GAMERA uses seventh-order upwind reconstruction by default, while retaining the original LFM eighth-order centered reconstruction for reference.The default reflects a deliberate update from the predecessor code rather than removal of the earlier scheme.
  • D. The Choice of Reconstruction Order: High-order reconstruction is selected to achieve low numerical diffusion with a reasonable amount of computing resource.This trade-off motivates the default scheme for large simulations.
  • D. The Choice of Reconstruction Order: The reconstruction choice is motivated by three-dimensional global-scale simulations of large space-plasma systems, including planetary magnetospheres and the heliosphere.The section frames computational efficiency as important for these applications.

1-D Linear Advection of Four Shapes in Non-uniform Grid

A non-uniform-grid advection test evaluates reconstruction orders using Gaussian, square-wave, triangular, and half-circle profiles. Higher order substantially improves smooth-profile resolution, while gains for square-wave transitions level off beyond fifth to seventh order.

  • 1-D Linear Advection of Four Shapes in Non-uniform Grid: The test uses a non-uniform grid, periodic boundaries, and four shapes chosen to probe smooth features, value discontinuities, gradient changes, and curvature discontinuities.The half-circle is especially challenging because of abrupt and gradual gradient changes on limited cells.
  • 1-D Linear Advection of Four Shapes in Non-uniform Grid: The first-order scheme severely smears all four profiles, whereas seventh order resolves the narrow Gaussian peak within nine cells.The test evolves four shapes to t = 2.0 on a non-uniform grid.
  • 1-D Linear Advection of Four Shapes in Non-uniform Grid: Gaussian-peak resolution improves through seventh order, but improvements above seventh order are small.The reported saturation concerns the narrow Gaussian profile.
  • 1-D Linear Advection of Four Shapes in Non-uniform Grid: Square-wave resolution improves with order, but its sharp right-side transition remains four cells wide above fifth order.The text characterizes further improvement in square-wave resolution as small.
  • 1-D Linear Advection of Four Shapes in Non-uniform Grid: The seventh-order scheme costs about 10% more total computing time than the fifth-order scheme, leading GAMERA to select it as the default.The timing includes reconstruction and limiting.

The Effective Diffusion Coefficient

The effective diffusion analysis quantifies numerical diffusion by fitting numerical step-advection solutions to an analytical advection-diffusion solution. Diffusion decreases strongly with reconstruction order, while improvements become small above seventh order.

  • The Effective Diffusion Coefficient: The analysis estimates numerical diffusion by fitting numerical solutions of a step-function advection problem to an analytical advection-diffusion solution.The setup uses a uniform domain with 256 cells and estimates D from the fitted solution.
  • The Effective Diffusion Coefficient: The effective diffusion coefficient decreases from 0.0604 to 0.0051 as reconstruction order increases from first to twelfth order.The values summarize the reconstruction schemes tested against the analytical advection-diffusion solution.
  • The Effective Diffusion Coefficient: At orders higher than seventh, the numerical diffusion width remains four cells, indicating small improvement in square-wave slope resolution.This quantitative result agrees with the four-shape advection simulations.
  • The Effective Diffusion Coefficient: Seventh-order reconstruction increases the effective grid Reynolds number by approximately two orders of magnitude relative to first order, with only small further gains above seventh order.The relative Reynolds number is normalized by the first-order result.
  • The Effective Diffusion Coefficient: The efficiency measure R/Order peaks at seventh order, supporting the choice of seventh-order reconstruction.The table reports diffusion coefficient, diffusion width, normalized Reynolds number, and R/Order across first- through twelfth-order schemes.
Loading 1810.10861v2…