… get off the grid ...
The simulation code GIZMO is a flexible, multi-method magneto-radiation-hydrodynamics+gravity code. The code lets you solve the hydrodynamic equations using a variety of different methods -- whatever is best for the problem at hand. In particular, it introduces a couple of new, Lagrangian Godunov-type methods, that allow you to solve the fluid equations with a moving particle distribution that is automatically adaptive in resolution and avoids the advection errors, angular momentum conservation errors, and excessive diffusion problems that seriously limit the applicability of “adaptive mesh” (AMR) codes, while simultaneously avoiding the low-order errors inherent to simpler methods like smoothed-particle hydrodynamics (SPH). But, if you want to use SPH -- either in “traditional” form or “modern” (more accurate) forms, or use a mesh, you can do that too with GIZMO! Meanwhile, self-gravity is solved fast, with a BH-Tree (optionally a hybrid PM-Tree for periodic boundaries), and on-the-fly adaptive gravitational softenings.
The code is descended from P-GADGET, itself descended from GADGET-2, and many of the naming conventions remain (for the sake of compatibility with the large library of GADGET work and analysis software). Currently available modules include things like:
- Hydrodynamics: with arbitrarily chosen solvers
- Magnetic fields: ideal or non-ideal MHD, including Ohmic resistivity, ambipolar diffusion, and the Hall effect
- Radiation hydro: full and flexible, with a variety of modular solvers and the ability to trivially add or remove frequencies to the transported set
- Cosmological integrations: large volume and cosmological “zoom in” simulations
- Galaxy/Star/Black hole formation and feedback: explicit modules for star and black-hole formation and feedback in galaxy-scale and cosmological simulations including detailed explicit and common sub-grid models
- Individual Star/Planet sink-particle formation: with explicit proto-stellar evolution and feedback, designed for individual star and planet formation simulations
- Self-interacting and scalar-field dark matter;
- Adaptive gravitational resolution: fully-Lagrangian softening for gravitational forces
- Braginskii conduction & viscosity: fully-anisotropic kinetic MHD, and other non-ideal plasma physics effects
- Subgrid-scale turbulent “eddy diffusion” models and the ability to run arbitrary shearing-box or driven turbulent boxes (large-eddy simulations) with any specified driving spectrum
- Dusty fluids: explicit evolution of aerodynamic particles in the Epstein and Stokes limits, with particle-particle collisions, back-reaction, Lorentz forces on charged particles, and more
- Degenerate/stellar equations-of-state: appropriate Helmholtz equations of state and simple nuclear reaction networks
- Ability to insert abitrary external gravitational fields, integration in non-standard cosmologies, cosmic rays, sink particles, and more
No, the code title is not an acronym, I just liked it. It refers both to the code's multi-purpose applications and to its historical relation ship to GADGET.
GIZMO is hosted on a bitbucket repository here. That’s where the real code action is! The development version is not public, but if you want to use it for a project, please contact me and we can set it up.
GIZMO User’s Guide (in download-able html format). Everything you need to know to run and/or analyze simulations with GIZMO. Includes a suite of examples and test problems as well. If you have questions/bug reports/etc that go beyond the user guide, please join the GIZMO Google Group (which serves as a forum, issue tracker, wiki, and mailing list, all in one).
GIZMO code/methods paper (PDF). This explains the details of the new methods, and explores them in a large suite of test problems, systematically comparing different ways to solve the hydrodynamic equations. Read this before anything else!
GIZMO-MHD code/methods paper (PDF). This explains the magneto-hydrodynamic implementation of the new methods and similarly explores them in a large suite of test problems. Read this if you wish to include magnetic fields in your simulations!
This is a publicly-available version of the code. It has full functionality for the “core” physics of the code: (magneto)hydrodynamics+gravity+cosmology. However, almost all of the additional physics modules (cooling, star formation, black holes, degenerate equations of state, conduction, dust-gas mixtures, radiation) are not available in this version of the code (for various reasons, see the User’s Guide for details). The intention is for this code to be used for purposes of learning and testing, code comparisons/methods studies, and development for new physics (as a basis for many experiments), as well as experiments which rely on ideal gas physics+gravity. If you want or need features in the full code, please contact me.
Some advantages of the new methods in GIZMO:
Angular Momentum Conservation:
Can we follow a gaseous disk?
Let’s consider what should be a straightforward test problem, of enormous astro-physical importance. A cold, inviscid, ideal-gas disk (pressure forces less important than the global gravitational forces on the gas) orbiting in a Keplerian (point-mass) potential on circular orbits, with no gas self-gravity. The gas should be stable and simply orbit forever. But what actually happens when we numerically evolve this?
“Standard” (non-moving) fixed-mesh or adaptive mesh codes:
Left: This movie shows the test problem, in a non-moving “grid” code (here, the code ATHENA). The movie shows the gas density, from 0 (black) to 2 (red). The initial condition has density =1 (green) inside a ring (the hole and outer cutoff are just for numerical convenience). We’ve set things up so it’s the “best case” for this sort of code: the resolution is high (256^2), we use a high-order method, and disk is only 2-dimensional and is perfectly aligned in the x-y Cartesian plane (so we don’t worry about vertical motion). The movie runs from t=0 to t=100 in code units, or about 10 orbits of the middle of the ring.
Almost immediately (within ~1 orbit) the disk starts to break up catastrophically. We can change the details by changing things like the order of the method, but it only makes things worse (the disk can look smoother, but all the mass diffuses out in the same timescale!). The problem is that mass must constantly be advected across the grid cells (at irregular angles, too). AMR-type codes do not conserve angular momentum, and have significant advection errors; so the disk breaks up and can’t be reliably evolved for much more than an orbit. Of course, you can help the problem by going to higher resolution, but the test problem here (because its so simple) is already much higher-resolution than most state-of-the art simulations of real problems! And in 3D, you also have a nasty grid-alignment effect: even at extremely high-res (>512^3) in AMR codes, the disk is torqued until the gas aligns with one of the Cartesian axes.
These are serious problems in almost all cosmological, or binary merger, or proto-planetary disk simulations -- anything where the disks should live a large number of orbits. If you’re lucky and can make a problem where you know the problem geometry exactly ahead of time, you can do better (using e.g. a cylindrical mesh) -- but that severely limits the usefulness of these methods.
The new methods we’ve developed for GIZMO avoid these errors entirely. One of the major motivations for the code was to develop a method which would conserve angular momentum accurately on general problems, and it makes a huge difference!
Right: The same simulation, evolved five times as long (~100 orbits)! The disk evolution is almost perfect: theres some noise, but almost no diffusion of the outer or inner disk, and angular momentum is nearly exactly conserved in both rings and total. Since there’s no grid, there’s also no “grid alignment” in 3D. And the conservation is great even at very low resolution: a 64^2 disk looks almost as good!
What about SPH?
If you want, you can use GIZMO to run your problem with the smoothed-particle hydrodynamics (SPH) method. In principle, SPH conserves angular momentum. But in practice, to make SPH actually work for real problems, you need to add an “artificial viscosity.” And an artificial viscosity prescription which allows you do deal with shocks correctly almost always leads to too much shear viscosity (which should be zero) in a flow like this. That triggers the “viscous instability” where sticky bits of the disk lead to it breaking up again.
Here, we see that: Left panel uses the old-fashioned “standard” artificial viscosity with a Balsara switch (density shown again, at time t~100). The disk breaks up catastrophically in ~2-3 orbits (but for different reasons than in the grid codes). Right panel uses the “state of the art” in SPH, a higher-order artificial viscosity “switch” from Cullen & Dehnen. Its much better than old-fashioned SPH or AMR codes, but the disk still starts to break up in ~5 orbits or so.
The Kelvin-Helmholtz Instability:
Standard Non-Moving Mesh GIZMO Meshless Method
Here’s a side-by-side comparison of a high-resolution (512^2) simulation of the Kelvin-Helmholtz instability, with the same initial conditions, but using two different numerical methods to solve the equations of motion. Two fluids (one light, one heavy) are in relative shear motion. The gas density is plotted, ranging between 1 and 2 (blue and red). On the left, a non-moving grid (from ATHENA). On the right, GIZMO’s new meshless finite volume method. Both develop the instability (correctly) at the same time, but in the code with a non-moving, Cartesian mesh, the contact discontinuities must be transported at irregular angles across the cells, which leads to a lot of (incorrect) excess diffusion even at this very high resolution, until the rolls simply “diffuse into” each other (in other words, the contact discontinuity cannot be maintained on small scales), and the whole instability blurs into a smooth “boundary layer” between the two fluids. With the new method, these discontinuities can be followed in exquisite detail. In fact, we had to massively down-grade the resolution on the right to make this movie (while, because of all the diffusion, there isn’t any detail lost on the left)! Click here for a higher-resolution version of the movie.
Here is a zoom-in of just part of one of the KH “whorls” -- still downgraded in resolution! Follow this link to get the a PDF with the full-resolution image of the final snapshot time.
What about SPH?
It’s well known that old-fashioned “traditional” SPH has serious problems with the KH instability. Here the left and right panels show “traditional” (GASOLINE or GADGET-2 type) and “modern” SPH (labeled TSPH and PSPH, respectively). “Modern” SPH refers to SPH using our more accurate “pressure” formulation of the equations of motion (Hopkins), fully conservative corrections for variable smoothing lengths (Springel), a higher-order smoothing kernel (Dehnen & Aly), artificial conductivity terms (Price), higher-order switches for the diffusion/viscosity/conductivity terms (Cullen & Dehnen), and more accurate integral/matrix-based gradient estimators (Garcia-Senz, Rosswog).
In “traditional” SPH, we see the usual problem. If the seed perturbation is small, the low-order errors of the method suppress its growth (the ‘surface tension’ effect in particular is dramatic) and lead to almost no mode growth. In “modern” SPH, the situation is tremendously improved. However, it still is not perfect: the mode amplitude is not exactly correct, and there is some unphysical “jagged” noise resulting from the artificial conductivity terms. More importantly, even in modern SPH with all the other improvements present, only going to very large SPH neighbor number can “beat down” the zeroth-order errors in quantities like the pressure gradients. Thus, if we re-run with a low initial mode amplitude (say, sub-percent, as opposed to a large ~10% initial perturbation), with the same neighbor number (number of particles the SPH particle must “talk to” each hydro interaction) as our GIZMO methods (32, in this problem), we see dramatic suppression of the KH instability! But going to large neighbor numbers is very expensive. And if we wanted to capture a still-smaller amplitude mode, we’d have to go to higher neighbor number still!
What about AMR?
Finally, just in case you’re wondering, going to Adaptive Mesh Refinement (AMR) schemes doesn’t actually help with the problems we described for grid/non-moving mesh codes. The reason is that, unless you just refine everywhere over and over (i.e. just increase the resolution of the simulation -- that’s cheating!), the contact discontinuities will still end up irregularly distributed over and moving over the cartesian mesh edges. It’s not a matter of missing “more cells” in a dense region. This shows that, by showing a grayscale density projection of one of the same tests run in the AMR code RAMSES. The serious diffusion affecting the rolls is obvious -- in fact, its worse here, because at refinement boundaries AMR codes become effectively lower-order accurate.
Some Other Examples:
The “blob test” - showing a cold cloud being disrupted by a supersonic wind. The two methods (left and right) are two different GIZMO methods -- both capture the correct physical destruction of the cloud in this test.
Fluid Mixing & Galilean Invariance:
Evolution of the Rayleigh-Taylor instability in our new GIZMO method: balancing a heavy fluid against gravity on top of a light fluid. The instability develops and buoyant plumes rise up, triggering secondary Kelvin-Helmholtz instabilities.
This is what the same instability looks like at time =3.6 in other state-of-the-art methods.
Top Left: “Modern” SPH (as described above), but using the same number of neighbor particles as our GIZMO method. Because of low-order errors in SPH, the method simply cannot capture small-amplitude instabilities (or numerically converge) without dramatically increasing the number of neighbors inside the interaction kernel (greatly increasing the CPU cost).
Top Right: Modern SPH, with a higher-order kernel with many more neighbors. Now the instability is captured, but low-order errors still corrupt the short-wavelength modes (giving them their odd boundaries)
Bottom Left: A static (non-moving) mesh result. Here things look quite nice, even though there is more diffusion in the “rolls” of the secondary Kelvin-Helmholtz instabilities that develop
Bottom Right: The same non-moving mesh, if we simply “boost” the system (set all the fluid in uniform motion). The solution should be identical, but since the fluid is now moving across instead of with the grid, the diffusion is dramatically amplified and the instability is suppressed. We have to go to >10,000 times as many cells to recover the same accuracy! GIZMO, being Lagrangian, has no such velocity-dependent errors.
Implosions & Explosions:
Here are two tests involving strong super-sonic shocks (projection of density). Left is the Noh implosion test (fluid rushing to the origin at infinite Mach number). The box is spherical, but we just show a slice through the x-y plane, and show one quadrant. Right is the Sedov-Taylor explosion (the result of a massive point injection of energy; now the full x-y plane is shown).
Top: MFM is one of our new GIZMO methods. Note that both solutions display excellent symmetry, smooth flows and well-resolved shock fronts.
Middle: SPH. The symmetry is statistically good, but there is far more (unphysical) noise, because of SPH’s low-order errors. Also, some “wall-heating” leads to the cavity in the center of the Noh test.
Bottom: Grid methods. Here there is more noise still than MFM, and more obviously, severe “carbuncle” and “grid alignment” instabilities, where the Cartesian grid imprints special directions and locations along the grid axes instead of the correct spherical structure!
Check It Out!