Maestro: A low-Mach number stellar hydrodynamics code

Ann Almgren, John Bell, Andy Nonaka, (CCSE/LBL), Mike Zingale (SUNY SB)

Most of the evolution of leading up to the explosion of a Type Ia supernovae occurs at Mach numbers of 0.01 to 0.1. Simulating these processes with traditional compressible hydro codes is challenging, because of the need to follow the soundwaves. We are developing a new low Mach number algorithm for stellar flows that will allow for the efficient simulation of low speed convective flows in stars.

Introduction to low-Mach number hydrodynamics

Figure 1: Mach number for 3 burning rising bubbles from a compressible code (top) and a low Mach number hydrodynamics method (bottom).

A large number of interesting astrophysical phenomena occur at low Mach numbers. Evolving these flows with a fully compressible simulation code is inefficient, because of the need to follow the sound waves. For an explicit time-discretization (i.e., the new state is expressed solely in terms of the present state), a fundamental limitation exists on the size of the allowable timesteps -- the CFL condition. A timestep is restricted such that information may only propagate across one zone in the computational grid per timestep. In compressible flow, information propagates at the speeds: u, u + c, and u - c, where c is the sound speed. Mathematically the timestep restriction is expressed as

Δ t < min { Δ x / (|u| + c) }
For very low Mach number flows, this is
Δ t ~ Δ x / c
This means that for an interface moving at a Mach number M << 1, it takes 1/M timesteps for that interface to move just one zone!

Our desire is to reformulate the equations of hydrodynamics to filter out sound waves, while retaining the compressibility effects important to the problem at hand. This will result in a timestep constraint of the form

Δ t < min{ Δ x / |u| }
Therefore, our algorithm will require far fewer timesteps (~1/M fewer) to simulate low Mach number flows.

Mathematically, we can think of low Mach number flows as exhibiting instantaneous acoustic equilibriation. Figure 1 shows this graphically. Three temperature perturbations were placed in a stratified atmosphere, seeding nuclear reactions. The heat release makes the bubbles buoyant. The top panel shows the Mach number for the compressible solution---a clear signal resulting from the finite propagation speed of the soundwaves is seen surrounding each bubble. In the low panel, a low Mach number method was used. Owing to the instantaneous acoustic equilibriation, the velocity field is realized throughout the whole domain---even far from the bubbles. However, the flow at the bubbles is identical between the two codes, showing that in this case, the soundwaves are not imporatant to the dynamics.

There are several popular methods for filtering soundwaves, which we summarize below.

Incompressible hydrodynamics

The simplest low Mach number approximation is incompressible hydrodynamics. This approximation is formally the M → 0 limit of the Navier-Stokes equations. In incompressible hydrodynamics, the velocity satisfies a constraint equation

∇ ⋅ U = 0
which acts to instantaneously equilibrate the flow, thereby filtering out soundwaves.

The constraint equation implies that

Dρ/Dt = 0
which says that the density is constant along particle paths. This means that there are no compressibility effects modeled in this approximation.

Anelastic hydrodynamics

In the anelastic approximation small amplitude thermodynamic perturbations are carried with respect to a hydrostatic background. The density perturbation is ignored in the continuity equation, resulting in a constraint equation

∇ ⋅ (ρ0U) = 0
This properly captures the compressibility effects due to the stratification of the background. Because there is no evolution equation for the perturbational density, approximations are made to the buoyancy term in the momentum equation.

The low Mach number combustion model

rt

Figure 2: two- and three-dimensional reactive Rayleigh-Taylor instabilities modeled with the smallscale low Mach number combusion code.

In the low Mach number combustion model, the pressure is decomposed into a dynamic, π, and thermodynamic component, p0, the ratio of which is O(M2). The total pressure is replaced everywhere by the thermodynamic pressure, except in the momentum equation. This decouples the pressure and density and filters out the sound waves. Large amplitude density and temperature fluctuations are allowed. The only requirement is that the total pressure stay close to the background pressure, which is assumed constant. This requirement can be expressed as

p = p0,
and differentiating this along particle paths leads to a constraint on the velocity field:
∇ ⋅ U = S
This looks like the constraint for incompressible hydrodynamics, but now we have a source term, S, representing the compressibility effects due to the energy generation and thermal diffusion.

Since the background pressure is taken to be constant, we cannot model flows that cover a large fraction of a pressure scale height. However, this method is ideal for exploring the physics of flames. We formulated this algorithm for astrophysical flows and used it to explore the dynamics of Rayleigh-Taylor unstable flame fronts in Type Ia supernovae in two- and three-dimensions:

Our small-scale, low Mach number astrophysical combustion algorithm is described in the following paper:

Adaptive Low Mach Number Simulations of Nuclear Flame Microphysics [Science Direct]
Bell, J. B., Day, M. S., Rendleman, C. A., Woosley, S. E. & Zingale, M. A. 2004, J. Comput. Phys., 195, 2, 677.

Papers describing scientific results with this algorithm include:

Direct Numerical Simulations of Type Ia Supernovae Flames I: The Landau-Darrieus Instability [ads]
Bell, J. B., Day, M. S., Rendleman, C. A., Woosley, S. E. & Zingale, M. A. 2004, ApJ, 606, 1029.
Direct Numerical Simulations of Type Ia Supernovae Flames II: The Rayleigh-Taylor Instability [ads]
Bell, J. B., Day, M. S., Rendleman, C. A., Woosley, S. E. & Zingale, M. A. 2004, ApJ, 608, 883.
Three-Dimensional Numerical Simulations of Rayleigh-Taylor Unstable Flames in Type Ia Supernovae [ads]
Zingale, M., Woosley, S. E., Rendleman, C. A., Day, M. S., & Bell, J. B. 2005, ApJ, 632, 1021.
The physics of flames in Type Ia supernovae [IOP Online]
Zingale, M., Woosley, S. E., Bell, J. B., Day, M. S., & Rendleman, C. A. 2005, J. Phys.: Conf. Ser. 16 405
Turbulence-Flame Interactions in Type Ia Supernovae [ads]
Aspden, A. J., Bell, J. B., Day, M. S., Woosley, S. E., & Zingale, M. 2008, ApJ, 689, 1173.

Low Mach Number Stellar Hydrodynamics Algorithm Description

To model the full star, we begin with the low Mach number combustion model described above, but no longer assume that the background pressure is constant, but instead, is hydrostatically stratified. This results in a constraint on the velocity field of:
∇ ⋅ ( β0U) = &beta0[ S - 1/(Γ1p0) ∂p0/∂t]
Here, β0 is a density-like variable. Again, large amplitude density and temperature variations are allowed. The only restriction is that the pressure remain close to the background pressure. In the presence of heating, the star will expand, and therefore we need to evolve the background state in response to the local heating. As with anelastic hydrodynamics, this constraint incorporates the effects of the background stratification. However, in contrast to anelastic, the right-hand side is non-zero, and represents the compressibility effects due to heating and the change in the background state.

We call the resulting simulation code Maestro. Our algorithm for solving the resulting system of equations is developed in a series of papers:

Low Mach Number Modeling of Type Ia Supernovae. I. Hydrodynamics [ads]
Almgren, A. S., Bell, J. B., Rendleman, C. A., & Zingale, M. 2006, ApJ, 637, 922.
Low Mach Number Modeling of Type Ia Supernovae. II. Energy Evolution [ads]
Almgren, A. S., Bell, J. B., Rendleman, C. A., & Zingale, M. 2006, ApJ, 649, 927.
MAESTRO: A Low Mach Number Stellar Hydrodynamics Code [IoP Online]
Almgren, A. S., Bell, J. B., & Zingale, M. 2007, Proceedings of SciDAC 2007, Journal of Physics: Conference Series, 78, 012085
Low Mach Number Modeling of Type Ia Supernovae. III. Reactions [ads]
Almgren, A. S., Bell, J. B., Nonaka, A., & Zingale, M. 2008, ApJ, 684, 449.
Low Mach Number Modeling of Type Ia Supernovae. IV. White Dwarf Convection [ads]
Zingale, M., Almgren, A. S., Bell, Nonaka, A., & Woosley, S. E. 2009, ApJ, 704, 196.

Papers describing scientific results with this algorithm include:

Astrophysical Applications of the Maestro Code [IoP Online]
Zingale, M., Almgren, A. S., Bell, J. B., Malone, C. M., & Nonaka, A. 2008, J. Phys.: Conf. Ser. 125 012013
A New Low Mach Number Approach in Astrophysics [CiSE]
Almgren, A. S., Bell, J. B., Nonaka, A., & Zingale, M. 2009, Computing in Science and Engineering, vol. 11, no. 2, 24

Hydrodynamics Tests

bubble comparison

Figure 3: A comparison of several hydrodynamics algorithms modeling a rising bubble

To test this code, we perform comparisons against fully compressible algorithms in a variety of setups. To test just the hydrodynamics, we put a hot bubble in a white dwarf atmosphere (taken to be plane-parallel), and watch the evolution as it buoyantly rises. Figure 3 shows a comparison between two fully compressible codes (PPM and the unsplit method of Colella 1990), our low Mach number method, traditional anelastic, and incompressible:

This bubble reached speeds of up to Mach 0.2. We see that the low Mach number method shows good agreement with the two compressible codes. The anelastic method does not get the rise velocity correct, because of approximations to the buoyancy term. Finally, the incompressible method completely fails to capture the right behavior (as expected).

Energy Evolution in Plane-Parallel Geometries

1-d hydrostatic readjustment test

Figure 4: Density vs. height for a 1-d heating problem using a fully compressible code (solid line), a low Mach number method that allows the base state to expand (dotted), and a low Mach number method without expansion (dashed).

When there is heating, a star or atmosphere will expand. If the expansion is significant, for the low Mach number to remain valid, we much adjust the hydrostatic base state in response to the large scale heating.

Figure 4 shows the result of heating a white dwarf atmosphere (again in a plane-parallel geometry) for 5 s. The change in density is shown. Here, the solid black line is the fully compressible (PPM) solution. The dotted line that lies directly on top of it is the low Mach number model where we've allowed the base state to adjust -- as we see, we have almost perfect agreement. For contrast, the dashed line is the low Mach number model where we did not allow the base state to evolve. We see that this sharply disagrees with the correct solution.

An analogous procedure can be derived for spherical problems.

Reactions

The algorithm has been extended to include heat release from reactions. Strang spliting is used to incorporate the reactions in a second-order accurate manner. Comparisons between fully compressible and our low Mach number method demonstrate that we accurately capture the dynamics of reacting bubbles in conditions appropriate to a white dwarf (see Figure 5).

reaction comparison

Figure 5: A comparison of compressible and low Mach number methods of a reacting bubble test, showing the temperature (top) and ash mass fraction (bottom).

Applications

Type Ia supernovae

Maestro has been extended to handle spherical geometries (represented in Cartesian coordinates). The resulting simulation is being used to model convection in white dwarfs leading up to Type Ia supernovae.

Figure 6: the radial velocity field (red: outward flowing; blue: inward flowing) in a convecting white dwarf shown at different times.

The image above shows the radial velocity in a convecting white dwarf leading up to the point of ignition of a Type Ia supernovae. The red contours are radially-outward moving fluid and the blue contours are inward moving fluid. At early times, we see the characteristic dipole pattern seen in previous studies. However, we note that these are the first such calculations to model the entire star.

Acknowledgements

This work was supported by the Applied Mathematics Program of the DOE Office of Mathematics, Information, and Computational Sciences under the U.S. Department of Energy under contract No. DE-AC02-05CH11231 to LBNL and by a DOE Office of Nuclear Physics Outstanding Junior Investigator award (No. DE-FG02-06ER41448) to SUNY Stony Brook.