Phurbas: An Adaptive, Lagrangian, Meshless, Magnetohydrodynamics Code. I. Algorithm

We present an algorithm for simulating the equations of ideal magnetohydrodynamics and other systems of differential equations on an unstructured set of points represented by sample particles. Local, third-order, least-squares, polynomial interpolations (Moving Least Squares interpolations) are calculated from the field values of neighboring particles to obtain field values and spatial derivatives at the particle position. Field values and particle positions are advanced in time with a second order predictor-corrector scheme. The particles move with the fluid, so the time step is not limited by the Eulerian Courant-Friedrichs-Lewy condition. Full spatial adaptivity is implemented to ensure the particles fill the computational volume, which gives the algorithm substantial flexibility and power. A target resolution is specified for each point in space, with particles being added and deleted as needed to meet this target. Particle addition and deletion is based on a local void and clump detection algorithm. Dynamic artificial viscosity fields provide stability to the integration. The resulting algorithm provides a robust solution for modeling flows that require Lagrangian or adaptive discretizations to resolve. This paper derives and documents the Phurbas algorithm as implemented in Phurbas version 1.1. A following paper presents the implementation and test problem results.

💡 Research Summary

**
The paper introduces Phurbas, a novel adaptive, Lagrangian, mesh‑free algorithm for solving the equations of ideal magnetohydrodynamics (MHD) and related systems of partial differential equations. Traditional Eulerian grid codes suffer from two major drawbacks: the Courant‑Friedrichs‑Lewy (CFL) time‑step restriction depends on both the signal speed and the bulk flow speed relative to the grid, and numerical diffusion is not Galilean invariant, leading to excessive diffusion in high‑velocity or highly sheared flows. Smoothed Particle Hydrodynamics (SPH) overcomes the first issue by moving with the fluid, but it lacks zeroth‑order consistency, suffers from noise in constant fields, and is limited by fixed particle masses and isotropic resolution.

Phurbas addresses these limitations by representing the fluid with a set of moving sample points (particles) that carry the primitive variables (density ρ, velocity V, magnetic field B, internal energy σ). The core numerical technique is a third‑order Moving Least Squares (MLS) interpolation: for each particle, a local polynomial fit is constructed from neighboring particles using a weighted least‑squares minimization. This fit yields the field values and their first and second spatial derivatives at the particle location with guaranteed zeroth‑order (and higher) consistency, eliminating the systematic errors that plague SPH.

Time integration is performed with a second‑order predictor‑corrector scheme applied directly to the Lagrangian form of the MHD equations. Because particles move with the local fluid velocity, the CFL condition depends only on the local signal speed (e.g., the fast magnetosonic speed), allowing substantially larger time steps in problems with strong bulk motions such as Keplerian disks or supersonic jets.

Adaptivity is achieved through a two‑stage process. The user specifies a target spatial resolution (particle number density) as a function of position. During each global step the algorithm evaluates a “void‑clump” criterion: if the distance to the nearest neighbor exceeds a prescribed threshold, a new particle is inserted; if particles are too close, the most redundant particle is removed. New particles inherit their physical quantities from the MLS interpolation of surrounding particles, while deletions are performed so as to conserve mass, momentum, and energy as closely as possible. This dynamic refinement/coarsening maintains the desired resolution without a fixed mesh.

To ensure numerical stability, Phurbas incorporates dynamic artificial viscosity, bulk viscosity, and thermal diffusion terms. These coefficients are locally adaptive, scaling with the divergence and shear of the flow, thereby damping shocks and steep gradients without overwhelming the physical dynamics. The magnetic field divergence constraint (∇·B = 0) is handled by a hyperbolic‑parabolic cleaning method (Dedner et al. 2002), which propagates divergence errors away as damped waves and can be supplemented by additional cleaning steps if needed.

The algorithmic workflow for each global time step is:

1. Initialize or update particle positions and primitive variables.
2. Perform MLS interpolation for each particle to obtain spatial derivatives of all fields.
3. Compute artificial dissipation terms and evaluate the right‑hand sides of the Lagrangian MHD equations, including second‑order time derivatives (provided in the appendix).
4. Execute the predictor step to obtain provisional values.
5. Apply the corrector step using the provisional values to achieve second‑order accuracy.
6. Conduct the void‑clump detection, inserting or deleting particles as required.
7. Re‑evaluate the CFL condition and adjust the global time step.
8. Apply divergence cleaning to the magnetic field.
9. Advance to the next step.

By combining high‑order consistent spatial reconstruction, Lagrangian time integration, and fully automatic adaptive particle management, Phurbas delivers a mesh‑free framework that is both accurate and efficient for a wide range of astrophysical fluid problems. It retains Galilean invariance, reduces numerical diffusion, and focuses resolution where the fluid concentrates or where the user desires, making it especially suitable for simulations of collapsing clouds, accretion disks, turbulent magnetized plasmas, and other multi‑scale, high‑Mach‑number phenomena. The present paper details the mathematical derivation and algorithmic design; a companion paper (McNally et al. 2012) presents implementation specifics and a comprehensive suite of verification tests.