Welcome

I am a postdoctoral researcher in the School of Natural Sciences at the Institute for Advanced Study, and a Visiting Research Collaborator in the Department of Astrophysical Sciences at Princeton University.

My research is on numerical methods for astrophysical fluid dynamics. I am especially interested in high-order methods, which reach a given accuracy with far fewer degrees of freedom than the second-order schemes most codes rely on, and in making them robust enough for real astrophysical flows, with their shocks, steep gradients and scales that cannot be fully resolved. I am equally interested in performance: how these methods are implemented and optimized so that they run efficiently on modern supercomputers, GPUs in particular. On the application side, my background is in planet formation and planet–disc interaction.

Previously: postdoctoral researcher at Princeton University (2022–2025) and at the Institute for Computational Science, University of Zurich (2019–2022). PhD in Computational Modelling and Scientific Computing, UAEM, Cuernavaca, Mexico (2019).

Research

High-order methods and their robustness

High-order methods reach a given accuracy with far fewer degrees of freedom than second-order schemes, but at shocks and other discontinuities they produce spurious oscillations, which can drive the density or the pressure negative. I work on the Spectral Difference method, a discontinuous finite-element scheme of arbitrary order, combined with a posteriori limiting: every update is tested against numerical and physical admissibility criteria, and in the cells that fail, the high-order fluxes are replaced by those of a second-order MUSCL–Hancock scheme defined on the sub-cells of each element. The solution keeps its high order where the flow is smooth and remains stable at discontinuities. With Romain Teyssier and Maria Han Veiga, we have developed the method for the Euler equations, for low Mach number flows such as stellar convection, and for the Navier–Stokes equations.

The Spectral Difference method, interactively →
The fallback scheme, interactively →

Multigrid solvers for self-gravity

Self-gravity requires the solution of Poisson's equation for the gravitational potential at every time step. Relaxation methods such as Jacobi or Gauss–Seidel solve it with simple local updates, but they only damp the error efficiently at the scale of the mesh spacing: its smooth, long-wavelength part decays slowly, and the number of iterations needed to converge grows with the resolution. Multigrid applies the same relaxation on a hierarchy of coarser meshes, where those long wavelengths become short and are damped quickly, so that a converged solution is reached in a few cycles, whatever the resolution. I work on the multigrid self-gravity solver of AthenaK, on adaptively refined meshes and on GPUs: the interpolation at the boundaries between refinement levels, which has to be conservative for orbits not to drift, and the performance of the solver on large GPU clusters. It is now used in work on tidal disruption events and on black holes inside massive stars.

Relaxation and multigrid, interactively →

GPU codes at scale

I write and optimize simulation codes for GPUs: spd_K, a Kokkos implementation of the Spectral Difference method, the nested-mesh version of FARGO3D, and contributions to AthenaK and RAMSES. On these machines most of the performance is decided by how a code uses the hardware rather than by the number of operations it performs: the memory access pattern and the data locality of its kernels, the number of kernels launched per step, and the way ghost zones are exchanged between mesh blocks and between ranks.

Planets in gas and dust discs

A planet embedded in a disc heats the gas around it, and that heat changes the force the gas exerts back. With Frédéric Masset I studied these thermal forces, the dynamical friction on luminous bodies, and how accretion of pebbles can drive a low-mass planet onto an eccentric orbit, using nested meshes to resolve scales far below the planet's Hill radius.

Interactive explainers

Three short pages where the algorithm runs live in the browser and you control it.

The Spectral Difference method

Drag the solution inside an element, follow the four steps of an update, then advect a square wave at orders 1 to 8 with and without a posteriori limiting.

The fallback scheme

Follow one stage as every control volume is checked and, where it fails, recomputed at second order; then watch the limiter stay silent until a shock forms.

Publications

First-author papers

With collaborators

Movies

Fluid instabilities

Kelvin–Helmholtz instability. Density, fourth-order Spectral Difference with a posteriori limiting, 800×800 elements (3200² degrees of freedom). The method of Paper I; see where the limiter acts.
Rayleigh–Taylor instability. Density with a second-order Godunov method at 512, 1024, 2048 and 4096 cells.
Viscous Kelvin–Helmholtz instability. The dye in the test of Lecoanet et al. (2016) at Reynolds number 10⁵, with both streams at the same density: explicit viscosity and diffusion give the problem a well-defined solution to converge to. Fourth-order Spectral Difference with the Navier–Stokes terms of Paper III, run with spd_K at 2048×4096 degrees of freedom; one of the two shear layers is shown.
The same test with a density jump. One stream is now twice as dense as the other, and nothing else changes. The billow rolls up as before and then breaks into small-scale structure. This case is far harder to converge, and at this resolution the details of the break-up are not.
Taylor–Green vortex in three dimensions. A volume rendering of the vorticity magnitude at Reynolds number 1600 and Mach 0.1: the initial vortex is stretched into sheets, the sheets roll up into tubes, and the tubes break down into turbulence. Fourth-order Spectral Difference with the Navier–Stokes terms of Paper III, 128³ degrees of freedom; the colour scale follows the growth of the vorticity.
A hot bubble rising at low Mach number, and the turbulence it leaves. A convective layer between two stably stratified ones, with gravity pointing down. The colour is the temperature relative to the hydrostatic equilibrium, within ±0.4 per cent: the bubble starts as a perturbation of one per cent, rises at a few per cent of the sound speed, rolls up into a vortex pair and breaks apart under the stable layer. Its debris stirs the whole convective layer into turbulence, while gravity waves grow in the stable layers until they saturate this colour scale. Well-balanced fourth-order Spectral Difference scheme, 768×384 degrees of freedom, from the tests of Paper II.

Shocks at high order

Double Mach reflection, and where the limiter acts. A Mach 10 shock runs along a wall. Top: density with the fourth-order Spectral Difference scheme, 800×200 elements. Bottom: the control volumes that fell back to second order (red, drawn enlarged). They trace the shocks and stay under half a per cent of the cells; everywhere else the scheme remains fourth order. From the tests of Paper I.
The same problem at orders 2, 4 and 8. A close-up of the Mach stem and the wall jet, with the same number of degrees of freedom in the three runs (3200×800). The shear layer stays smooth at second order and rolls up at fourth and eighth.
Mach 3 wind tunnel with a step. Density with the same scheme, 600×200 elements: the bow shock, its reflections off the walls, and the shear layer that rolls up along the top.

Planets, discs and dynamical friction

A planet crossing a dust ring. Dust surface density around a planet of four Earth masses on an eccentric orbit, in the frame that follows the planet. At each periastron it dives into the ring, accretes, and pulls out a trail of dust. FARGO3D with three nested mesh levels, from Velasco-Romero et al. (2024).
The same run in three dimensions. A volume rendering of the dust density, with the vertical direction stretched five times: the dust layer is only about a hundredth of the orbital radius thick. Each passage of the planet lifts a thin shell of dust several scale heights out of the midplane, and the shell then drifts along the ring.
The same ring, a planet twice as massive. With eight Earth masses the planet disturbs the ring much more: it leaves dust-poor lanes behind it, and the mark of each passage is still visible when it returns.
Gas flowing past a massive, luminous body. A passive tracer at three luminosities, increasing from left to right. FARGO3D with nested meshes, from the work on dynamical friction with radiative feedback.

Self-gravity with multigrid, in AthenaK

Refinement must not perturb the orbit. Gravitational potential of a binary, with mesh blocks that refine around each mass and follow it. Left: mass-conserving prolongation at coarse–fine boundaries; the masses stay on the analytic orbit (dashed). Right: plain interpolation; the orbit drifts. How the solver works →
Collapse and fragmentation of a rotating cloud. The Boss–Bodenheimer test: a rotating isothermal cloud with a small two-armed perturbation collapses into a bar, and the bar fragments. Face-on column density over 2.5 free-fall times; the white dots are sink particles, which form in the densest gas, accrete and merge. AthenaK with multigrid self-gravity.
Adaptive refinement in action. The same run with its mesh drawn. Every square is a block of 16³ cells; a block splits in two along each direction wherever the gas would otherwise not resolve its Jeans length, down to five levels. The fine blocks follow the bar, then each fragment along its orbit, and merge back into coarse ones behind them.
Turbulent fragmentation under self-gravity. Column density of a turbulent, self-gravitating cloud at 2048³ cells, seen along two axes. Filaments that exceed their local Jeans mass collapse; the multigrid solver supplies the potential at every step.
Around a fragmenting cloud. Volume rendering of the gas density in a turbulent, magnetised, self-gravitating box at 1024³ cells, when the first cores have collapsed: filaments in purple and orange, the densest cores in green and white.

Black-hole accretion, in AthenaK

An accretion flow and its jets. Volume rendering of a general-relativistic MHD simulation, within 64 gravitational radii of the black hole at t = 8000 M: gas density in orange, and in blue the funnel where the magnetic field dominates over the gas. Eight nested mesh levels of 512³ cells each.

Gallery

FARGO3D with nested meshes

A planet of eight Earth masses crossing a dust ring, with the nested meshes drawn

A planet of eight Earth masses crossing a dust ring: FARGO3D with three nested mesh levels (lines every eight cells)

Flow around a luminous perturber

Dynamical friction on a luminous body: FARGO3D with nested meshes

Nested meshes around a planet

Nested meshes around a pebble-accreting planet

Spectral difference: spd and spd_K

Kelvin–Helmholtz instability

Kelvin–Helmholtz instability: Spectral Difference, 4th order, 3200×3200 degrees of freedom

Zoom into a Kelvin–Helmholtz billow down to the sub-element control volumes

Zooming into a billow, down to the control volumes inside each element

Dye concentration in a viscous Kelvin–Helmholtz billow with a density jump, after it has broken up

Viscous Kelvin–Helmholtz instability with a density jump, Re = 105: Spectral Difference, 4th order, 2048×4096 degrees of freedom, with spd_K

Double Mach reflection: the Mach stem and the wall jet

Double Mach reflection of a Mach 10 shock: Spectral Difference, 4th order, 6400×1600 degrees of freedom

Turbulent convection

Turbulent convection at low Mach number: Spectral Difference, 4th order, 512×1024 degrees of freedom

Rayleigh–Taylor instability at four resolutions

Rayleigh–Taylor instability: Godunov 2nd order at 512, 1024, 2048 and 4096 cells

AthenaK

Dye concentration in a Kelvin–Helmholtz billow with a density jump

Viscous Kelvin–Helmholtz instability with a density jump, Re = 105: AthenaK, 8192×16384 cells

Density of a GRMHD accretion flow in the meridional plane: the disc and the evacuated funnel

GRMHD accretion onto a black hole: density in the meridional plane, ±64 M, at t = 8000 M. AthenaK, eight nested levels of 5123 cells

Density of a GRMHD accretion flow in the equatorial plane

The same flow in the equatorial plane, ±64 M: spiral streams and low-density cavities

Inner region of a GRMHD accretion flow in the equatorial plane

The inner ±16 M of the same flow, down to the black hole

Column density of a self-gravitating magnetised turbulent cloud

Gravo-magneto-turbulent fragmentation with sink particles: column density. AthenaK with multigrid self-gravity, 10243 cells

Density slice of a self-gravitating magnetised turbulent cloud

A density slice through the same run

Density slice of a turbulent self-gravitating box with its adaptive mesh, and a zoom onto a collapsing core

Adaptive refinement onto a collapsing core: a density slice of a self-gravitating turbulent box (left) and the dashed region enlarged (right). AthenaK, 10243 base mesh, blocks of 643 cells, three refinement levels in this plane

Density and magnetic field lines around an accreting magnetised star

Accretion onto a magnetised star with an aligned dipole: density and field lines in the meridional plane. AthenaK, eight refinement levels

Closed dipole loops and open field lines around an accreting magnetised star

Close-up of the magnetosphere: closed dipole loops, the disc and the open polar field

Equatorial density around an accreting magnetised star

The same run in the equatorial plane: the cavity carved by the magnetosphere

Software