PaperlayerView exact revision on arXiv

Source-linked AI summary

FARGO3D: A new GPU-oriented MHD code

Pablo Benítez-Llambay, Frédéric Masset

arXiv:1602.02359v1astro-ph.IMastro-ph.EPphysics.comp-ph

TL;DR

FARGO3D addresses the computational cost of three-dimensional protoplanetary-disk and planet-interaction simulations. It combines staggered finite-difference HD/MHD methods, orbital advection and GPU-oriented execution with automatic CPU-to-GPU translation. The code conserves mass and angular momentum to computer accuracy, while single-precision calculations can exhibit resolution-dependent numerical problems.

  • Problem

    Three-dimensional planet-disk calculations are computationally expensive, motivating codes that exploit orbital advection and GPU resources for these studies.

  • Method

    FARGO3D uses staggered, dimensionally split HD/MHD algorithms with orbital advection and automatically translates parallel CPU routines into CUDA for GPU execution.

  • Results

    Mass and angular momentum are conserved to computer accuracy, and the code is presented with tests assessing its HD/MHD implementation and orbital-advection properties.

  • Takeaways & Limitations

    FARGO3D provides a GPU-capable tool for computationally demanding protoplanetary-disk and planet-disk interaction studies while remaining adaptable through CPU-side routine development.

  • Takeaways & Limitations

    Single-precision calculations can develop numerical problems at high resolution, and reducing resolution mitigates but does not solve the issue.

Abstract

from arXiv · show

We present the FARGO3D code, recently publicly released. It is a magnetohydrodynamics code developed with special emphasis on protoplanetary disks physics and planet-disk interactions, and parallelized with MPI. The hydrodynamics algorithms are based on finite difference upwind, dimensionally split methods. The magnetohydrodynamics algorithms consist of the constrained transport method to preserve the divergence-free property of the magnetic field to machine accuracy, coupled to a method of characteristics for the evaluation of electromotive forces and Lorentz forces. Orbital advection is implemented, and an N-body solver is included to simulate planets or stars interacting with the gas. We present our implementation in detail and present a number of widely known tests for comparison purposes. One strength of FARGO3D is that it can run on both "Graphical Processing Units" (GPUs) or "Central Processing unit" (CPUs), achieving large speed up with respect to CPU cores. We describe our implementation choices, which allow a user with no prior knowledge of GPU programming to develop new routines for the CPU, and have them translated automatically for the GPU.

1. INTRODUCTION

FARGO3D is designed for computationally demanding three-dimensional protoplanetary-disk and planet-interaction studies by combining orbital advection with GPU execution. Its staggered, dimensionally split approach emphasizes low memory use, accurate disk equilibria and entropy advection, while CPU routines can be translated automatically to GPU code.

  • Three-dimensional magnetized disk studies are computationally expensive because resolving the vertical dimension is required for magneto-rotational instability.
  • GPU execution can make three-dimensional calculations comparable in cost to two-dimensional calculations on a CPU core.
  • FARGO3D solves hydrodynamics and magnetohydrodynamics on Cartesian, cylindrical or spherical meshes, with special emphasis on protoplanetary disks and forming-planet interactions.
  • CPU-developed routines follow strict syntax rules that expose parallelism and are automatically translated into CUDA, allowing users without GPU expertise to modify the code.
  • The dimensionally split method minimizes temporary arrays and memory use, while staggered meshes support accurate treatment of steady disk and planetary-envelope equilibria.
  • Orbital advection is included in both the hydrodynamic and magnetohydrodynamic implementations, extending the predecessor algorithm to MHD with an upstream-averaged electric field.

2. OVERVIEW

FARGO3D solves hydrodynamic or magnetohydrodynamic equations on several mesh geometries using staggered placement of vector and scalar variables. Its formulation includes rotation, viscosity, magnetic induction, internal-energy evolution and selectable equations of state.

  • FARGO3D solves hydrodynamic or magnetohydrodynamic equations on Eulerian Cartesian, cylindrical or spherical meshes, including lower-dimensional configurations.
  • The governing equations include continuity, Navier–Stokes dynamics, external forces, mesh rotation and viscous stress.
  • The energy equation uses a non-conservative formulation based on volumic internal energy.
  • When MHD is enabled, FARGO3D solves the induction equation and closes the governing equations with an equation of state.
  • The public release provides locally isothermal and adiabatic or ideal equation-of-state forms, and coordinate directions are labeled X, Y and Z across geometries.
  • Vector quantities are stored on cell faces, whereas scalar quantities are cell centered; edge-centered electromotive forces support the staggered magnetic-field representation.

3. METHODS

FARGO3D advances hydrodynamic and magnetohydrodynamic fields through operator-split source and transport updates on an Eulerian mesh. Its discretization combines finite differences, finite volumes, upwind methods, MOC, and constrained transport, with geometry-specific transported quantities.

  • Overall algorithm: FARGO3D uses an explicit, operator-split scheme with upwind techniques on an Eulerian mesh.Some equations are discretized as finite differences, while others use finite volume methods.
  • Operator splitting: Operator splitting decomposes the governing operator into partial problems whose successive solutions form the numerical update.The source step maps Q0 to Q1, and the transport step maps Q1 to the final state Q2.
  • Overall algorithm: A complete update applies a source step followed by a transport step, with an additional magnetic-field interaction substep for MHD.The full update is organized as successive substeps in the flow chart.
  • Hydrodynamic formulation: In spherical geometry, the solver represents continuity, momentum, and energy equations using coordinate-specific differential operators and transported variables.Transported variables include density, momentum components, specific angular momentum, and energy when the energy equation is solved.
  • Hydrodynamic formulation: The transport formulation embeds several velocity-product and Coriolis terms in conservative updates, while other terms remain in the source step.This allocation is important for planet-disk applications.

3.3. Source step

The source step applies non-transport terms separately from conservative advection, using staggered fields to evaluate pressure, potential, geometric, viscous, and MHD-related contributions. The transport step uses dimensionally split one-dimensional updates and preserves momentum properties through consistent control volumes.

  • Source step: The source step updates terms excluded from transport, including pressure gradients, gravitational forces, viscosity, pressure work, and MHD Lorentz contributions.FARGO3D separates these effects from the conservative transport update.
  • Source step: Three source substeps apply pressure and gravity, artificial viscosity and heating, then pressure work when the energy equation is solved.These correspond to cells 7a, 7b, and 7c of the update flow chart.
  • Transport step: Dimensionally split transport reduces multidimensional advection to successive one-dimensional flux updates using intermediate-stage quantities.The ordering of these operations is characteristic of dimensionally split schemes.
  • Transport step: The upwind transport method evaluates face-centered quantities at half time using the staggered velocity field rather than a Riemann solver.The interface flux is constructed from the transported quantity and the velocity normal to the face.
  • Momentum conservation: The implementation preserves momentum to machine accuracy by using consistent control volumes and averaging star values instead of mass fluxes.This makes momentum transport exactly consistent with zone-centered quantities.
  • Momentum conservation: Angular momentum remains conserved to machine accuracy in rotating frames, supporting correct isothermal shock jump conditions and vortensity production.The result holds after accounting for pressure-gradient source terms and mesh periodicity.

3.5. Stability

FARGO3D controls explicit integration with CFL-based limits and orbital advection that removes large, nearly uniform azimuthal motion from the stability constraint. The method reduces timestep restrictions but has a documented limitation for steepened interpolation in Keplerian disks.

  • CFL stability: The explicit timestep is limited by a CFL condition so information travels no more than one cell per step.The maximum timestep is chosen heuristically from process-specific limits over the computational domain.
  • CFL stability: MHD timestep limits include fast magnetosonic waves, fluid motion, artificial viscosity, viscosity, and resistivity.The fast-wave speed combines sound and Alfvén velocities.
  • Orbital advection: Orbital advection decomposes azimuthal motion into a large uniform component and a smaller residual velocity, effectively using a nearly corotating frame for each ring.The uniform component is handled by a circular shift, while the residual uses advection machinery.
  • Limitation: The steepened Piecewise Parabolic Advection method is unsuitable for Keplerian disks because it spuriously generates small-scale vortices.The standard fractional-shift procedure requires one additional ghost or buffer-zone layer.
  • Orbital advection: The orbital shift does not contribute to the Courant condition, so the stability criterion uses the residual velocity instead of total azimuthal velocity.The shift is a circular permutation and does not introduce numerical errors.
  • Orbital advection: Orbital advection conserves the transported fields to machine accuracy, including mass and angular momentum.The implementation applies the algorithm to the same variables used in the transport formulation.

3.6. Magnetohydrodynamics

FARGO3D evolves magnetic fields with Method of Characteristics and constrained transport, combining characteristic estimates of electromotive forces with edge-centered flux updates. The constrained-transport construction maintains initially divergence-free magnetic fields to machine accuracy, while separate steps handle resistivity and Lorentz forces.

  • Magnetic-field update: The MHD solver uses Method of Characteristics with constrained transport to update magnetic fields and preserve zero divergence.The method evaluates half-timestep edge quantities needed for the induction update.
  • Magnetic-field update: Constrained transport evolves face-centered magnetic fluxes from edge-centered electromotive forces using a discrete circulation.Oppositely oriented edge contributions cancel in the zone-boundary flux sum.
  • Divergence preservation: Initially divergence-free magnetic fields remain divergence-free at all times to machine accuracy under the constrained-transport update.The cancellation follows from sharing electromotive forces consistently on zone edges.
  • MHD source terms: Resistive terms are added before constrained transport, while the Lorentz-force update applies magnetic-tension and magnetic-pressure contributions to velocity.These operations are separated into dedicated MHD substeps.
  • Method of characteristics: The Method of Characteristics reconstructs transverse velocity and magnetic-field values along Alfvén characteristics to evaluate the electromotive force.Characteristic speeds use the normal magnetic field and density, with piecewise-linear reconstruction using van Leer slopes.

3.7. FARGO algorithm in MHD - Orbital Advection

FARGO3D extends orbital advection to MHD using constrained transport, preserving the magnetic field’s divergence-free property. In the no-shear case, the magnetic update reduces to an essentially exact azimuthal shift, while residual arcs use interpolated electromotive forces.

  • Method: FARGO3D implements orbital advection for MHD in a staggered-mesh code using a constrained-transport method.The implementation preserves the divergence-free property of the magnetic field and adapts the MHD orbital-advection method to staggered magnetic-field centering.
  • Method: The orbital-advection equations split the dynamics into residual-velocity evolution and a separate azimuthal, axisymmetric orbital-velocity step.The residual system is solved separately, while the magnetic supplement handles the orbital component.
  • Method: In the shearing-sheet formulation, the transverse magnetic field is advected at uniform speed v0, while shear supplies a source term for the longitudinal component.The formulation uses the divergence-free constraint in deriving the component equations.
  • Method: For uniform flow, upstream averaging of the electric field recovers the exact magnetic-field variation for any finite time step.This holds because the velocity is uniform along streamlines and constant in time, making the upstream average equal to the Eulerian time average.
  • Result: In the no-shear case, constrained transport acts as a circular permutation of cell values with no numerical diffusivity or impact on the Courant condition.When the upstream arc is not an integer number of cells or shear is present, the residual-arc electromotive force is evaluated with piecewise parabolic interpolation.

3.8. Orbital integrator for planets

FARGO3D evolves planetary systems with a fifth-order Cash–Karp integrator whose timestep is governed by the CFL condition, while coupling planetary motion to gas forces evaluated during hydrodynamic updates. The integrator is modular, but its fixed timestep does not guarantee accurate close-encounter solutions.

  • Planetary-system integration: FARGO3D supports an arbitrary number of point-like bodies orbiting a central mass and interacting with the gas.
  • Planetary-system integration: Planet positions and velocities are evolved with a fifth-order Cash–Karp Runge–Kutta method using a fixed timestep governed by the CFL condition.
  • Planetary-system integration: The Cash–Karp integrator is decoupled from the rest of the code and can be replaced by a more sophisticated solver when close encounters or collisions must be detected.
  • Gas–planet coupling: Gas forces on planets are evaluated once per hydrodynamical timestep and used to update planetary velocities, either by direct summation over all cells or by removing the axisymmetric gas-density component before summation.
  • Implementation: The planetary positions and velocities are updated at cell 5a of the algorithm shown in Fig. 2.

4. IMPLEMENTATION CONSIDERATIONS

FARGO3D combines MPI domain decomposition with automatic CPU-to-CUDA translation for GPU execution. Its implementation emphasizes low memory use, tunable kernel performance, and substantial acceleration, while exposing trade-offs in communication and setup-dependent throughput.

  • MPI and domain decomposition: FARGO3D uses C with MPI and slab decomposition in the Y and Z directions, keeping the full X direction on one processing element.The mesh is split radially and vertically in cylindrical coordinates, and in radius and colatitude in spherical coordinates.
  • MPI and domain decomposition: Communications occur twice per time step with three ghost layers, balancing communication frequency against the thickness of exchanged boundary regions.More frequent communications permit thinner ghost layers, while additional physics may require more layers.
  • Automatic GPU translation: FARGO3D automatically converts computationally expensive CPU mesh functions into CUDA kernels and wrapper functions using structured input/output declarations.The parser maps loop indices to CUDA thread and block indices and handles required data movement and kernel setup.
  • Automatic GPU translation: The CPU implementation uses a parallelizable mesh-function structure with explicit INPUT and OUTPUT declarations for fields consumed and produced by each routine.These declarations support the corresponding GPU conversion and data-transfer bookkeeping.
  • Automatic GPU translation: CUDA code is generated during compilation from C comments that follow strict syntax rules separating internal variables, kernel arguments, and pre- and post-operations.The public mesh functions follow these rules and can run on the GPU; the same design could support an OpenCL parser.
  • Performance and memory: Kernel block-size exploration provides a further 10–20 % speed up relative to a unique default size, while one GPU typically achieves 40× the execution speed of a similar-generation CPU.The speedup depends on the problem and mesh size, with larger meshes producing larger speedups.
  • Performance and memory: GPU performance and memory footprint depend strongly on setup geometry and dimensionality, with ghost-zone work causing apparent throughput reductions in some configurations.On a K20x with 6 Gb of RAM, an isothermal three-dimensional Cartesian HD simulation of about 320^3 cells can fit after reserving overhead.

5. TESTS

FARGO3D is evaluated with standard hydrodynamic and MHD tests against analytic solutions and established reference results. The reported results show good agreement across shock, diffusion, rotor, vortex, and magnetic-loop problems, with localized deviations and oscillations noted.

  • Test suite: The test suite covers published hydrodynamic and MHD benchmarks organized broadly by increasing setup dimensionality.The hydrodynamic tests include the Sod shock tube and viscous ring, while the MHD tests include Riemann, resistive, rotor, vortex, and loop problems.
  • Hydrodynamic tests: In the 300-cell Sod shock tube, the shock spans 2–3 zones, the contact discontinuity spans 7–8 zones, and δρ/ρ is offset by about 10^-2 from the exact solution.The authors report good agreement despite a small deviation in specific energy.
  • Hydrodynamic tests: The viscous ring test shows good agreement between numerical density evolution and the analytic expectations for an axisymmetric Keplerian ring.The simulation uses constant viscosity ν0 = 10^-5 on a 512-cell cylindrical mesh.
  • MHD tests: The Brio & Wu MHD Riemann test exhibits overall good behavior, despite small velocity oscillations around 500 ≤ z ≤ 650.The setup uses 800 equally spaced cells and reflecting boundaries.
  • MHD tests: The resistive current-sheet test shows good agreement between the numerical solution and the analytic diffusion solution.With zero velocity, the induction equation reduces to a diffusion equation, allowing direct comparison with the analytic profile.
  • MHD tests: The MHD rotor test validates the solver in Cartesian and cylindrical coordinates, using a dense rapidly spinning cylinder to test strong torsional Alfvén-wave propagation.The reported results are closely comparable with reference results, using a modified setup with a central hole.
  • MHD tests: The magnetic-loop results show magnetic-energy decay comparable to reference solvers for horizontal advection and preserve central velocity-node symmetry in the evolving field lines.At t = 2, the reported magnetic energy is Em ≈ 1.235 · 10^-7; at later times, four islands remain along the discontinuities.

6. DISCUSSION

The tests validate FARGO3D across hydrodynamic and magnetohydrodynamic disk problems, while discussion focuses on orbital advection, turbulence, spectral resolution, and precision effects.

  • Vortex in a two-dimensional Keplerian disk: Orbital advection dramatically improves vortensity conservation at all tested resolutions in the two-dimensional Keplerian vortex test.Compared with PLUTO, FARGO3D performs slightly better without orbital advection and slightly worse with it.
  • MRI in unstratified disks: The standard advection scheme is more diffusive away from corotation, yielding lower α than orbital advection except at r = 3.The two values coincide at corotation.
  • MRI in unstratified disks: The MRI calculations reach saturated turbulence later without orbital advection, while the orbital-advection case produces a higher α plateau.The standard case levels off at a slightly smaller α and saturates later.
  • MRI in unstratified disks: The outer-ring magnetic spectrum in the FARGO-MHD calculation resolves smaller scales by roughly a factor of 2, while the inner-ring spectra are similar.The orbital-advection case also shows systematically more high-azimuthal-wavenumber signal in the first ring.
  • Precision effects: Single precision introduces vorticity errors comparable to variations between neighboring radial bins, making vorticity and vortensity non-monotonic.The resulting vortensity extrema can trigger Rossby wave instability at many radii.

7. PERSPECTIVES

The perspectives section describes ongoing extensions to FARGO3D, including radiative transfer, nested meshes, and additional multifluid and non-ideal-MHD physics.

  • Ongoing physics extensions: FARGO3D includes a gray, flux-limited-diffusion radiative-transfer module with separate gas and photon temperatures.The module targets radiative processes relevant to protoplanetary disks and forming planets.
  • Resolution and mesh structure: Nested meshes are used to increase resolution near embedded planets beyond what uniform global meshes can provide.The approach is motivated by the high resolution required for embedded-planet calculations.
  • Ongoing physics extensions: Multifluid capability, Hall effects, and ambipolar diffusion have been implemented and are slated for future publications.

APPENDIX GEOMETRIC COEFFICIENTS

The geometric-coefficient implementation reduces GPU memory use by deriving zone volumes and face surfaces from one-dimensional arrays.

  • Geometric coefficients: One-dimensional arrays store geometric factors whose products provide face surfaces and zone volumes during mesh updates.The design aims to let these arrays fit in fast GPU constant memory for relatively modest setups.

Surfaces

The surface relationships define geometric coefficients for cell faces, including the lower face perpendicular to each coordinate direction.

  • Surfaces: The coefficients s_j/k and the indexed surface S^A(j,k) describe the lower face area of a cell perpendicular to direction A.

Volumes

FARGO3D uses direct expressions for inverse cell volumes because divisions are computationally costly. The coefficients depend on the geometry and are tabulated for the implemented cases.

  • The implementation computes the inverse of a cell volume directly rather than obtaining it through division.This choice reflects the stated high computational cost of divisions.
  • The resulting expressions reproduce the exact surfaces and inverse volumes in all cases.

VISCOUS STRESS TENSOR

The viscous stress tensor is implemented across Cartesian, cylindrical, and spherical geometries, with components placed on a staggered grid and used directly to update velocity. The implementation also decomposes operations to expose parallelism and avoid race conditions.

  • VISCOUS STRESS TENSOR: The viscous stress tensor is expressed separately for Cartesian, cylindrical, and spherical geometries, matching the code implementation.The listed expressions are written exactly as implemented.
  • VISCOUS STRESS TENSOR: Table 5 supplies one-dimensional coefficients for the three geometries, with coordinate meanings defined by geometry-specific conventions.The coefficient depends on either j or k, and the min/max subscript identifies the lower or upper cell-side location.
  • VISCOUS STRESS TENSOR: Diagonal stress components are zone-centered, while cross components lie at edge midpoints according to their coordinate orientation.For example, τxy is located at the midpoint of the lowest-x, lowest-y edge along z within a cell.
  • VISCOUS STRESS TENSOR: After calculation, stress-tensor components directly update the corresponding velocity components.The update is identified as cell 7d in the paper’s algorithmic flow and is written through partial-derivative equations for the three geometries.
  • VISCOUS STRESS TENSOR: Artificial viscous pressure and its gradient are split across two files to avoid race conditions and expose parallelism in each substep.The decomposition is described as common in the implementation.
  • VISCOUS STRESS TENSOR: The implementation uses separate source files for viscous-tensor calculations and velocity updates in Cartesian, cylindrical, and spherical coordinates.The files are listed for the tensor and corresponding velocity-update operations.
Loading 1602.02359v1…