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

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.

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

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:

rt

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

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.
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
Zingale, M., Almgren, A. S., Bell, Nonaka, A., & Woosley, S. E. 2009, ApJ, submitted.

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

Verification and Validation

Hydrodynamics

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. The figure below 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:

bubble comparison

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

When there is heating, the base state will adjust. The figure below shows the result of heating a white dwarf atmosphere (again in a plane-parallel geometry) for 5 s. The change in density is shown:

1-d hydrostatic readjustment test

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.

We can look at the evolution in two-dimensions where there is localized heating, as shown below. Here, three localized heating sources and one broad heating source at constant height were used.

2-d heating test

The red contours are the fully compressible solution and the green are the low Mach number stellar hydrodynamics solution. Again we see excellent agreement in the development of the buoyant plumes.

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:

reaction comparison

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.

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.