Well-balanced treatment of gravity in astrophysical fluid dynamics simulations at low Mach numbers
P. V. F. Edelmann, L. Horst, J. P. Berberich, R. Andrassy, J. Higl, C. Klingenberg, F. K. Roepke
AAstronomy & Astrophysics manuscript no. wellbalancing © ESO 2021March 1, 2021
Well-balanced treatment of gravity in astrophysical fluid dynamicssimulations at low Mach numbers
P. V. F. Edelmann , L. Horst , J. P. Berberich , R. Andrassy , J. Higl , C. Klingenberg , and F. K. Röpke , X Computational Physics (XCP) Division and Center for Theoretical Astrophysics (CTA), Los Alamos National Laboratory, LosAlamos, NM 87545, USAe-mail: [email protected] Heidelberger Institut für Theoretische Studien, Schloss-Wolfsbrunnenweg 35, D-69118 Heidelberg, Germany Department of Mathematics, Würzburg University, Emil-Fischer-Str. 40, D-97074 Würzburg, Germany Zentrum für Astronomie der Universität Heidelberg, Institut für Theoretische Astrophysik, Philosophenweg 12, D-69120 Heidel-berg, GermanyReceived xxxx xx, xxxx / accepted xxxx xx, xxxx ABSTRACT
Context.
Accurate simulations of flows in stellar interiors are crucial to improving our understanding of stellar structure and evolution.Because the typically slow flows are but tiny perturbations on top of a close balance between gravity and pressure gradient, suchsimulations place heavy demands on numerical hydrodynamics schemes.
Aims.
We demonstrate how discretization errors on grids of reasonable size can lead to spurious flows orders of magnitude faster thanthe physical flow. Well-balanced numerical schemes can deal with this problem.
Methods.
Three such schemes are applied in the implicit, finite-volume code SLH in combination with a low-Mach-number numericalflux function. We compare how the schemes perform in four numerical experiments addressing some of the challenges imposed bytypical problems in stellar hydrodynamics.
Results.
We find that the α - β and Deviation well-balancing methods can accurately maintain hydrostatic solutions provided that grav-itational potential energy is included in the total energy balance. They accurately conserve minuscule entropy fluctuations advectedin an isentropic stratification, which enables the methods to reproduce the expected scaling of convective flow speed with the heatingrate. The Deviation method also substantially increases accuracy of maintaining stationary orbital motions in a Keplerian disk on longtime scales. The Cargo–LeRoux method fares substantially worse in our tests, although its simplicity may still o ff er some merits incertain situations. Conclusions.
Overall, we find the well-balanced treatment of gravity in combination with low Mach number flux functions essentialto reproducing correct physical solutions to challenging stellar slow-flow problems on a ff ordable collocated grids. Key words. hydrodynamics – methods: numerical – convection
1. Introduction
Astrophysical modeling often involves self-gravitating fluids.They are commonly described by the equations of fluid dynam-ics – viscous e ff ects are negligible in most astrophysical sys-tems and therefore the non-viscous Euler equations are used –with a gravitational source term. Such systems can attain station-ary equilibrium configurations in which a pressure gradient bal-ances gravity, i.e. hydrostatic equilibrium. A prominent exampleare stars, modeled in classical approaches as spherically sym-metric gaseous objects. Apart from this dimensional reduction,the assumption of hydrostatic equilibrium considerably simpli-fies the modeling of the – in reality rather complex – structureof stars. The resulting equations of stellar structure (e.g., Kip-penhahn et al. 2012) enable successful qualitative modeling ofthe evolution of stars through di ff erent stages. The price for thissuccess is a parametrization of multidimensional and dynami-cal processes that limits the predictive power of such theoreticalmodels and requires their calibration with observations. Recentattempts to simulate inherently multi-dimensional and dynami-cal processes, such as convection in stellar interiors (e.g., Brown-ing et al. 2004; Meakin & Arnett 2006, 2007; Woodward et al.2013; Rogers et al. 2013; Viallet et al. 2013; Pratt et al. 2016; Müller et al. 2016; Cristini et al. 2017; Edelmann et al. 2019;Horst et al. 2020), try to overcome this shortcoming.Such simulations pose a number of challenges to the under-lying numerical techniques. Not only is the range of relevant spa-tial and temporal scales excessive, but the flows of interest arisein a configuration that is often close to hydrostatic equilibrium.This has two implications: (i) The schemes must be able to pre-serve hydrostatic equilibrium in stable setups over a long periodof time compared to the typical time scales of the flows of inter-est. (ii) The flow speed v expected to arise from a small perturba-tion of the equilibrium configuration should be slow comparedto the speed of sound c , thus the corresponding Mach number, M ≡ v / c , is expected to be low.If we decide to avoid approximating the equations and in-clude all e ff ects of compressibility, aspect (ii) above calls forspecial low-Mach-number solvers in numerical fluid dynamicscombined with time-implicit discretization to enable time stepsdetermined from the actual fluid velocity instead of the – forthe problem at hand irrelevant – speed of sound as requiredby the CFL stability criterion (Courant et al. 1928) of time-explicit schemes. Several suitable methods are implemented inthe S even -L eague H ydro (SLH) code. Numerical and theoreti- Article number, page 1 of 24 a r X i v : . [ a s t r o - ph . S R ] F e b & A proofs: manuscript no. wellbalancing cal details are discussed in Barsukow et al. (2017); Barsukowet al. (2017); Edelmann & Röpke (2016); Miczek et al. (2015);Edelmann (2014); Miczek (2013) and examples of applicationto astrophysical problems can be found in Horst et al. (2020);Röpke et al. (2018); Edelmann et al. (2017); Michel (2019); Bo-laños Rosales (2016).Aspect (i), however, also requires attention. The conditionfor hydrostatic equilibrium is part of the equations of stellarstructure, that are discretized and numerically solved in classi-cal stellar evolution modeling approaches. In contrast, hydro-static equilibrium is only a special solution to the full gravo-hydrodynamic system at the level of the partial di ff erential equa-tions, but it is not guaranteed that discretizations of these equa-tions can reproduce the physically correct equilibrium state. Thisis in particular the case because gravitational source terms areusually treated in an operator-splitting approach, resulting in dif-ferent discretizations of the pressure and gravity terms. Astro-physical fluid dynamics simulations often employ finite-volumeschemes, in which hydrodynamical flows are modeled with aGodunov-type flux across cell interfaces. Hydrodynamical quan-tities are therefore determined at these locations. The gravita-tional source term, in contrast, is discretized in a completely dif-ferent and independent way. In a second-order code, for exam-ple, it is often calculated using cell-averaged densities assignedto cell centers. In general, this procedure does not lead to an ex-act cancellation of gravity and pressure gradient in hydrostaticconfigurations. Spurious motions are introduced that mask thedelicate low-Mach-number flows arising from perturbations ofthis equilibrium, such as, for instance, convection driven by nu-clear energy release.To overcome the problem of aspect (i), so-called well-balancing methods have been introduced, which are numeri-cal methods that ensure exact preservation of certain stationarystates. Methods of this type have predominantly been developedfor the simulation of shallow-water-type models in order to re-solve stationary solutions such as the lake-at-rest solution with-out numerical artifacts (e.g., Brufau et al. 2002; Audusse et al.2004; Bermudez & Vázquez 1994; LeVeque 1998; Desveauxet al. 2016b; Touma & Klingenberg 2015; Castro & Semplice2018; Barsukow & Berberich 2020). These stationary states canbe described using an algebraic relation, which favors the devel-opment of well-balanced methods. In the simulation of hydrody-namics under the influence of a gravitational field, the situationis di ff erent, since hydrostatic solutions are described by a dif-ferential equation that admits a large variety of solutions thatdepend on temperature and chemical composition profiles, aswell as the equation of state (EoS). In practice, the concrete hy-drostatic profile is determined by equations describing physicalprocesses other than hydrodynamics and gravity, such as ther-mal and chemical transport and the change in energy and speciesabundance due to reactions.Di ff erent approaches can be used to deal with this: The ma-jority of well-balanced methods for the Euler equations withgravity, for example Chandrashekar & Klingenberg (2015);Desveaux et al. (2016a); Touma et al. (2016) and referencestherein, are designed to only balance certain classes of hydro-static states, often isothermal, polytropic, or isentropic stratifica-tions, under the assumption of an ideal gas EoS. However, formany astrophysical applications, in particular, cases involvinglate stellar evolutionary stages and massive stars, nonideal ef-fects of the gas may be important. In stellar interiors, the mostimportant additions to the ideal gas EoS are radiation pressureand electron degeneracy e ff ects. This requires a more complex– often in parts tabulated – EoS to properly describe the ther- modynamical properties of the gas. We discuss an example ofsuch an EoS in Section 2.2.2. Well-balanced methods whichare capable of balancing hydrostatic states for general EoS havebeen introduced by Cargo & Le Roux (1994); Käppeli & Mishra(2014, 2016); Grosheintz-Laval & Käppeli (2019); Berberichet al. (2018, 2019); Berberich et al. (2020); Berberich et al.(2021); Berberich (2021).Most methods that have been discussed in the astrophysi-cal context and literature (e.g., Zingale et al. 2002; Perego et al.2016; Käppeli et al. 2011; Käppeli & Mishra 2016; Popov et al.2019) balance a second-order approximation of the hydrostaticstate rather than the hydrostatic state itself. Another recent ap-proach is the well-balanced, all-Mach-number scheme by Pa-dioleau et al. (2019). None of these publications tested a low-Mach-number, well-balanced method in more than one spatialdimension in a stable stratification over long time scales. As weshow in this paper, long-term stability cannot be automaticallyinferred from one-dimensional (1D) tests, yet it is of fundamen-tal importance for applications in stellar astrophysics.Using a staggered grid, which in this context means storingpressure on the cell interfaces instead of the cell centers, can al-leviate some of the problems of well-balancing the atmosphere,as shown, for example, in the MUSIC code (Go ff rey et al. 2017,Sec. 6). For this approach, it still has to be shown that convectivevelocities scale correctly with the strength of the driving forceat low Mach numbers, which we found very challenging in ourapproach, see Section 5.3.The methods introduced in Berberich et al. (2018, 2019,2021) can balance any hydrostatic stratification exactly. The onlyassumption is that the hydrostatic solution to be balanced isknown a priori. This poses no severe restriction for many as-trophysical applications where the initial condition is often con-structed under the assumption of hydrostatic equilibrium. Anexample are simulations of stellar convection, where the initialmodel is commonly derived from classical stellar evolution cal-culations that by construction impose hydrostatic equilibrium.Here, we discuss three possible well-balancing methods thatfollow rather di ff erent approaches. The first method extends thework of Cargo & Le Roux (1994) which only applied to 1Dsetups into the three-dimensional (3D) case and achieves well-balancing by modifying the pressure part of a general EoS. Werefer to this as the Cargo–LeRoux (CL) well-balancing method.The other two methods modify how variables are extrapolated tothe cell interfaces. We refer to them as the α - β well-balancing(Berberich et al. 2018, 2019) and the Deviation well-balancingmethod (Berberich et al. 2021). For these three schemes, we de-scribe their theoretical background and study their impact on theaccuracy of solutions to a set of simplified test problems, whichare designed to resemble typical situations in astrophysics.The structure of the paper is as follows: Section 2 reviews thebasic set of equations of fluid dynamics and their implications.It also introduces the notation that is used in the subsequent sec-tions. In Section 3 we discuss the discretization of these equa-tions and describe the AUSM + -up flux used in the later tests. Thewell-balancing schemes are introduced in Sections 4.1 to 4.3. InSection 5 we test the applicability of the well-balancing methodsand their performance in an extensive suite of simple applicationexamples. Conclusions are drawn in Section 6.
2. Equations of compressible, ideal hydrodynamics
This section introduces the general set of equations that aresolved with the SLH code in their formulation in general coordi-
Article number, page 2 of 24. V. F. Edelmann et al.: Well-balancing in astrophysical hydrodynamics nates. The following subsections closely follow the presentationof Miczek (2013).
We employ curvilinear coordinates x = ( x , y , z ) = ( x , x , x )with a smooth mapping x : R → R , x (cid:55)→ x ( ξ ) . (1)to Cartesian coordinates ξ = ( ξ, η, ζ ) = ( ξ , ξ , ξ ). The rea-soning here is that the coordinates ξ simplify the computations,while the coordinates x are adapted to the physical object, suchas a spherical star.The compressible Euler equations on curvilinear coordinatesthen read J ∂ u ∂ t + A ξ ∂ f ξ ∂ξ + A η ∂ f η ∂η + A ζ ∂ f ζ ∂ζ = J s , (2)with the vector u of conserved variables and the fluxes f ξ l givenby u = (cid:37)(cid:37) u (cid:37) v (cid:37) wE , f ξ l = (cid:37) n T ξ l v (cid:37) u n T ξ l v + ( n ξ l ) x p (cid:37) v n T ξ l v + ( n ξ l ) y p (cid:37) w n T ξ l v + ( n ξ l ) z p n T ξ i v ( E + p ) , (3)for l = , ,
3. Here, density and pressure are denoted by (cid:37) and p ,respectively. The velocity vector, expressed through its curvilin-ear components, reads v = ( u , v , w ) and enters the equation forthe total energy density E = (cid:37)(cid:15) + (cid:37) | v | + (cid:37)φ with the specificenergy (cid:15) and the gravitational potential φ .The Euler equations (2) in their curvilinear form in dependon the derivatives of the coordinate transformation. Its Jacobideterminant is J = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∂ x ∂ξ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:88) l , m , n = (cid:15) lmn ∂ x l ∂ξ ∂ x m ∂η ∂ x n ∂ζ , (4)where (cid:15) lmn is the three-dimensional Levi-Civita symbol. The nor-mal vector n ξ l and interface area A ξ l in ξ l -direction are n ξ l = JA ξ l ∂ξ l ∂ x ∂ξ l ∂ y ∂ξ l ∂ z , A ξ l = (cid:115)(cid:32) J ∂ξ l ∂ x (cid:33) + (cid:32) J ∂ξ l ∂ y (cid:33) + (cid:32) J ∂ξ l ∂ z (cid:33) . (5)External forces that enter Eq. (2) are – with an exception dis-cussed below – collectively denoted by the source term s . Whilethere are di ff erent possible contributions, for example energygeneration due to nuclear burning, gravity inevitably appears inany astrophysical setup. At the same time it might pose di ffi cul-ties for hydrodynamical codes to maintain hydrostatic solutionsto Eq. (2) (see Section 2.3) if it includes a strong gravitationalsource term as is common in the interior of stars. When gravityis the only source term, the expression for s reads s = − (cid:37) ∂φ∂ x − (cid:37) ∂φ∂ y − (cid:37) ∂φ∂ z . (6) This adds gravitational forces to the momentum part of Eq. (2).The presence of a gravitational field also a ff ects the evolutionof total energy. We include this e ff ect in the definition of thetotal energy rather than the source term, since this treatment sig-nificantly improved our accuracy in numerical experiments. Wefound this to be crucial in simulations of low Mach number con-vection.The source term adds gravitational force to the momentumequations. For our current treatment and test setups, the gravi-tational potential φ is fixed in time and changes the mass distri-bution during the simulation is excluded, meaning self-gravityis neglected. It is a reasonable simplification if the setup is veryclose to hydrostatic equilibrium and the change in the mass dis-tribution is negligible. Such an approximation, however, fails forstellar core simulations at later evolutionary stages, where asym-metries in the mass distribution may arise due to violent convec-tive motions. However, these setups are not the typical use casesof the well-balancing techniques presented in this paper as devi-ations from hydrostatic equilibrium are non-negligible.In our notation lower indices do not indicate partial deriva-tives to avoid confusion. The common choice to close the Euler system Eq. (2) is using anEoS. There are few physically relevant EoS which can be givenin a short, explicit analytical form. Two of these are discussedin the following. We assume that all components of the gas arein local thermodynamic equilibrium, that is they can all be de-scribed with a common temperature.
The ideal gas is one of the simplest EoS, yet with a wide rangeof applications. It describes an ensemble of randomly moving,noninteracting particles in thermodynamic equilibrium. It is anacceptable model for terrestrial gases, such as air, for which theinteractions between the particles are small. It serves well in thecase of a fully ionized plasma, such as in the interior of stars,as long as the e ff ects of degeneracy and radiation pressure aresmall.The ideal gas pressure is given by p ( (cid:37), (cid:15) ) = p ( (cid:37), T ( (cid:37), (cid:15) )) = R µ (cid:37) T ( (cid:37), (cid:15) ) , (7)with the temperature T ( (cid:37), (cid:15) ) = ( γ − µ R · (cid:15)(cid:37) . (8)The gas constant R for the ideal gas is8 .
314 462 618 153 24 × erg K − mol − . The specific heatratio γ depends on the internal degrees of freedom in the under-lying gas mixture, typical values are 5 / / temperature T instead of the explicit dependence of (cid:15) . Weuse this form of the EoS to formulate hydrostatic equilibria withcertain temperature profiles. While the ideal gas is a useful approximation of the EoS of stellarinteriors, it does not capture the e ff ect of partially or fully degen-erate electrons or of radiation pressure. A commonly used EoS Article number, page 3 of 24 & A proofs: manuscript no. wellbalancing that includes these e ff ects to great precision is the HelmholtzEoS (Timmes & Swesty 2000). It relies on a interpolation ofthe Helmholtz free energy from tabulated values using biquinticHermite polynomials. All other quantities are then derived fromexpressions involving derivatives of the Helmholtz free energy.This approach ensures that all thermodynamic consistency rela-tions are fulfilled automatically. This EoS has a wide range ofapplicability and serves as a typical example of a general tab-ulated EoS, contrasting our approaches to some well-balancedmethods relying on using a specific EoS, such as the ideal gas. Except for the very late stages of stellar evolution, stars can beconsidered as gaseous spheres, which change only over verylong time scales, much longer that those of the fluid motions.Dynamical processes acting in the interiors, as for example con-vective motions, however, evolve on much shorter time scales.Thus, in hydrodynamical simulations that aim to follow such fastprocesses, a star can, to first order, be described as a static strat-ification with pressure and density profiles constant in time, thatis v ≡ , (cid:37) ( t , x ) = (cid:37) ( x ) , and p ( t , x ) = p ( x ) . (9)These conditions reduce the first and the last part of Eq. (2) tothe trivial relations ∂ t (cid:37) = ∂ t ( (cid:37) E ) = . (10)The momentum equations lead to the hydrostatic equation ∇ p ( (cid:37), T ) = − (cid:37) ∇ φ. (11)This equation is invariant under transformations between di ff er-ent sets of curvilinear coordinates. A pair of constant-in-timefunctions (cid:37) and p , which satisfy Eq. (11) together with the cho-sen EoS is called hydrostatic solution or hydrostatic equilibrium .Since the EoS usually depends on temperature, there is in manycases a whole continuum of hydrostatic solutions rather thanuniqueness. Depending on the stratification, perturbations to the hydrostaticsolution may lead to dynamical phenomena. One important ex-ample is convection, where hydrostatic equilibrium is not per-fectly fulfilled anymore but deviations are small.The criterion for stability against convection is typically de-rived by considering the behavior of a small fluid element beingperturbed from the surrounding stratification. The frequency atwhich the element oscillates around its equilibrium position χ is called the Brunt–Väisälä frequency N . Its square is given by N = ∂φ∂χ (cid:37) ext (cid:32) ∂(cid:37) int ∂χ − ∂(cid:37) ext ∂χ (cid:33)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) χ , (12)where χ denotes the vertical coordinate , (cid:37) int is the density ofthe small fluid element, and (cid:37) ext is the density of the backgroundstratification. It is assumed that the fluid element changes its stateadiabatically, that is without exchanging heat with its surround-ing, and the derivative ∂(cid:37) int /∂χ is interpreted as the adiabaticchange of density while maintaining pressure equilibrium withthe background stratification at height χ . For the full derivation i.e. the direction opposing the vector of gravitational acceleration of Eq. (12) we refer the reader to any textbook on stellar astro-physics (e.g., Maeder 2009; Kippenhahn et al. 2012).It is common to express the gradients in Eq. (12) in termsof di ff erent variables. In the case of homogeneous composition( µ ( x ) = const.) Eq. (12) is equivalent to N = − ∂φ∂χ ∂∂χ (cid:32) sc p (cid:33) = − ∂φ∂χ δ T (cid:34) ∂ T ∂χ − (cid:32) ∂ T ∂χ (cid:33) ad (cid:35) , (13)with specific entropy s and specific heat at constant pressure c p and the equation of state derivative, δ = − ∂ ln (cid:37)∂ ln T , (14)which is 1 in the ideal gas case. The subscript “ad” denotes theadiabatic derivative as mentioned above.Another common form of this equation is using a variant ofthe temperature gradients expressed using pressure as a coordi-nate, ∇ = ∂ ln T ∂ ln p . (15)Using this definition Eq. (12) is equivalent to N = − ∂φ∂χ δ H p ( ∇ − ∇ ad ) , (16)with the pressure scale height, H p = − ∂χ∂ ln p = − p ∂χ∂ p . (17)We call a hydrostatic equilibrium stable with respect to convec-tion or convectively stable , if N ≥
0. Otherwise we call it un-stable with respect to convection or convectively unstable . Thisis a local definition, which means that a hydrostatic solution canbe convectively stable in one region and convectively unstable inanother. A suitable reference time for convectively stable setupsis the minimal Brunt–Väisälä time t BV = min x ∈ Ω t locBV ( x ) = π max x ∈ Ω N ( x ) . (18)It seems to be a natural time scale for the evolution of small per-turbations as explained for example in Berberich et al. (2019).For any hydrostatic solution, the value of N can be either cal-culated analytically (for simple EoS, like the ideal gas EoS) ornumerically for more complex EoS.Another useful timescale is the sound crossing time t SC through the domain. Similar to Käppeli & Mishra (2016), wedefine t SC in the direction of coordinate ξ l as t SC = ξ i , i (cid:44) l (cid:90) ξ lU ξ lL d ξ l c (cid:0) ξ , ξ , ξ (cid:1) , (19)with ξ lL and ξ lR being the lower and upper boundaries of the do-main in that direction. The speed of sound c is calculated usingthe equation of state. An expression for the sound speed using ageneral equation of state is given by c ( (cid:37), (cid:15) ) = (cid:115) ∂ p ( (cid:37), (cid:15) ) ∂(cid:37) + ∂ p ( (cid:37), (cid:15) ) ∂(cid:15) · (cid:15) + p ( (cid:37), (cid:15) ) (cid:37) , (20)where the function for the pressure p comes from an equation ofstate (see Section 2.2). Article number, page 4 of 24. V. F. Edelmann et al.: Well-balancing in astrophysical hydrodynamics
3. Discretization
Analytic solutions to Eq. (2) as for example given by the hy-drostatic solution of Section 2.3 are exceptions and require spe-cial initial conditions. To obtain more general solutions, whichalso allow for more complex dynamics such as turbulent convec-tion developed from perturbations in a hydrostatic stratification,Eq. (2) needs to be solved numerically. While there are severaldi ff erent numerical approaches, this section focuses on the meth-ods that are employed by the SLH code. For a more general in-troduction on this topic, see Toro (2009). For the numerical solution, the underlying equations have to bediscretized on a, possibly curvilinear, mesh that resembles thephysical spatial domain. This grid is then mapped to Cartesiancoordinates on which the computations are conducted. A set ofintegers ( i , j , k ) denotes the center of the ( i , j , k )-th cell while,for example, ( i + / , j , k ) denotes the interface between cell( i , j , k ) and ( i + , j , k ). The semi-discrete finite-volume schemeis obtained by integrating Eq. (2) over the cell volume in com-putational space, leading to V i jk ∂ U i jk ∂ t = − A i + , j , k ( ˆ F ξ ) i + , j , k + A i − , j , k ( ˆ F ξ ) i − , j , k − A i , j + , k ( ˆ F η ) i , j + , k + A i , j − , k ( ˆ F η ) i , j − , k − A i , j , k + ( ˆ F ζ ) i , j , k + + A i , j , k − ( ˆ F ζ ) i , j , k − + V i jk ˆ S i jk , (21)where V i jk is the volume of the corresponding cell in physicalspace and A i + , j , k , A i , j + , k , and A i , j , k + are the interface areas ofthe interfaces in ξ , η , and ζ -direction respectively. Details on thecomputation of cell volumes and interface areas are computedto second order following Kifonidis & Müller (2012). The cell-averaged source term is approximated to second order byˆ S i jk = (cid:37) i jk g x ) i jk ( g y ) i jk ( g z ) i jk , (22)where (cid:37) i jk is the cell-averaged value of density and the cell-centered gravitational acceleration ( g χ ) i jk = − ∂φ∂χ (cid:12)(cid:12)(cid:12)(cid:12) x ijk is com-puted analytically from the given gravitational potential φ .In Eq. (21), ˆ F ξ l is an approximation of the interface flux for l = , ,
3. There is some freedom in constructing the approxi-mate flux function that calculates ˆ F ξ l and many approaches canbe found in the literature. However, the specific choice is cru-cial for the accuracy of the numerical solution. This is furtherdiscussed in Section 3.2 in the context of flows at low Machnumbers M = | v | c , (23)where c is the speed of sound given by Eq. (20).The values that enter the approximate flux function need tobe reconstructed from the center of the cells to the correspond-ing interfaces. The reconstruction and the evaluation of the fluxis done for each coordinate direction separately, before the re-sulting fluxes over the surfaces are added for each cell. Table 1.
Parameter a l used in the computation of the signal velocity inEq. (27). condition a l n ξ l · v > n ξ l · g > + n ξ l · v ≤ n ξ l · g ≤ − n ξ l · v > n ξ l · g ≤ (cid:16) n ξ l · v (cid:17) n ξ l · g + c CFL ( ∆ ξ l ) ≤ + n ξ l · v > n ξ l · g ≤ (cid:16) n ξ l · v (cid:17) n ξ l · g + c CFL ( ∆ ξ l ) > − n ξ l · v ≤ n ξ l · g > (cid:16) n ξ l · v (cid:17) n ξ l · g + c CFL ( ∆ ξ l ) > − n ξ l · v ≤ n ξ l · g > (cid:16) n ξ l · v (cid:17) n ξ l · g + c CFL ( ∆ ξ l ) ≤ + ffi ciently accurate ODE solver thisdiscretization yields a second-order accurate scheme as has beennumerically verified by Berberich et al. (2019). For the testsin this article we mainly use the implicit second-order accu-rate three step Runge–Kutta method ESDIRK23 of Hosea &Shampine (1996).We chose an advective CFL time step (CFL u ), not strictly forreasons of stability, but as a good compromise between accuracyand e ffi ciency. In curvilinear coordinates it takes the form ∆ t CFL u = c CFL min i jkl ( ∆ ξ l ) i jk | n ξ l · v | i jk , (24)with a constant c CFL of order unity and an estimate of the celllength in direction ξ l given by( ∆ ξ ) i jk = V i jk (cid:16) A i − , j , k + A i + , j , k (cid:17) , (25)and accordingly for ( ∆ ξ ) i jk and ( ∆ ξ ) i jk .This time step criterion generally works well when the flowis fully developed, but it has problems when the Mach numberson the grid are very small (e.g. in the beginning of a simula-tion with zero initial velocities), because this yields very largeor infinite time steps. As a way to prevent this, Miczek (2013)suggests to include the free-fall signal velocity in the time stepcalculation. The so-called CFL ug time step is then given by ∆ t CFL ug = c CFL min i jkl ( ∆ ξ l ) i jk s i jkl , (26)with the signal velocity s i jkl = (cid:18) a l | n ξ l · v | + (cid:113) | n ξ l · v | + a l c CFL ( ∆ ξ l ) n ξ l · g (cid:19) i jk . (27)The parameter a l selects the right branch of the quadratic solu-tion and is given in Table 1.For Mach numbers close to M =
1, however, it is usuallymore e ffi cient to use explicit time stepping. For this we use thethird order accurate RK3 scheme of Shu & Osher (1988) with aCFL uc time step controlled by the fluid velocity and sound speed.It is given by ∆ t CFL uc = c CFL N dim min i jkl ( ∆ ξ l ) i jk | n ξ l · v | i jk + c i jk , (28)where N dim denotes the spatial dimensionality of the equations.In contrast to the previous criteria, this is a strict stability crite-rion for the explicit time stepping. Article number, page 5 of 24 & A proofs: manuscript no. wellbalancing
A fundamental part of the discretization is the choice of a nu-merical two-state flux. These fluxes give approximate solutionsof the two-state Riemann problem at the cell interfaces. Choos-ing di ff erent numerical fluxes yields di ff erent properties for thescheme. One of the most commonly used numerical flux func-tions for compressible Euler equations is the Roe flux (Roe1981), which has been developed to accurately capture shocks.However, at low-Mach numbers, the accuracy of the Roe fluxand similar methods su ff er from excessive Mach number de-pendent di ff usion, the origin of which is an upwind term in theschemes that is needed for numerical stability (e.g., Turkel 1987;Miczek et al. 2015). To correct this behavior, a number of lowMach number fixes have been proposed that aim on reducing theexcessive di ff usion and make it independent of the Mach number(e.g., Turkel 1987; Li & Gu 2008; Rieper 2011; Oßwald et al.2015; Miczek et al. 2015; Barsukow et al. 2017; Berberich &Klingenberg 2020).One peculiarity of astrophysical setups compared to, for ex-ample, setups in the engineering community is the presence ofstrong stratifications where pressure and density may change byorders of magnitudes within the computational domain. In suchsetups, the reduction of di ff usion comes with the risk of reducingstability and many of the schemes found in the literature developinstabilities. The SLH code is designed in a modular fashion thatfacilitates the implementation and testing of di ff erent types offlux functions. In numerical tests we find the so-called AUSM + -up method to yield appropriate results in the low-Mach regimein combination with the well-balancing method discussed here.The basic construction and a modification for improved low-Mach behavior is is discussed in the next section. + -up While the Roe solver and related methods are based on the ap-proximate solution of the Riemann problem at the cell inter-faces, a di ff erent approach is followed in the class of AdvectiveUpstream Splitting Methods (AUSM), which have been first in-troduced by Liou & Ste ff en (1993). In Liou (1996) the AUSMscheme was extended to AUSM + , the idea of which we willbriefly describe in the following. To be consistent with the orig-inal publication, we use dimensionalized quantities.The central idea is to split the analytical flux function f χ ofEq. (2) into a pressure and a mass flux via f x i = p e i + + ˙ m i ψ, (29)with˙ m i = (cid:37) v i , ψ = uvwE + p (cid:37) , i ∈ [1 , ,
3] (30)and the i -th canonical basis vector in the five-dimensional fluxvector space e i . This formulation is given for Cartesian coordi-nates, a transformation to curvilinear coordinates is possible.The pressure and mass flux of Eq. (29) are discretized sepa-rately which results in the numerical flux functionˆ F x i ( U L , U R ) = p / ( U L , U R ) e i + + ˙ m / ( U L , U R ) ψ up ( U L , U R ) , (31) where the upwind term ψ up is given by ψ up ( U L , U R ) = (cid:40) ψ ( U L ) if ˙ m / ( U L , U R ) ≥ ,ψ ( U R ) otherwise . (32)The core properties of this numerical flux function are deter-mined by the definition of the interface values p / and ˙ m / ofthe pressure p and the mass flux ˙ m i . With the initially proposeddefinitions of Liou (1996), this flux function is not capable ofresolving low Mach number flows. However, Liou (2006) ex-tended the AUSM + scheme to AUSM + -up with enhanced lowMach number capability.For AUSM + -up, the interface pressure is defined as p / = P + (5) ( M L ) p L + P − (5) ( M R ) p R − K u P + (5) ( M L ) P + (5) ( M R )( (cid:37) L + (cid:37) R )( f a c / )( u R − u L ) , (33)where P ± (5) are fifth degree polynomial functions, c / is an ap-proximation for the interface speed of sound, and K u is a constantthat can be set to a value between zero and unity. We refer thereader to Liou (2006) for the detailed definitions of the terms.The third term on the right hand side of Eq. (33) that includesvelocity-di ff usion is called u -term and is designed to reduce thenumerical dissipation at low Mach numbers. It involves a scalingfactor f a defined as f a = M o (2 − M o ) , (34)with M o = min [1 , max (Ma , M cut )] , (35)where M cut is a cut-o ff Mach number that ensures that f a doesnot approach 0 in the limit of very small Mach numbers. Thisis necessary to prevent singularities as the inverse of f a entersinto the mass-flux part [see Eq. (36)] in the original definitionof AUSM + -up. However, as described below, the SLH code setsthe scaling in these two parts independently such that M cut canbe theoretically set to zero. In SLH, for implementation reasonsthe value is set to a small value, typically to 10 − , to avoid diver-gence at smaller Mach numbers. This could easily be changed,but does not have a practical influence on our calculations. Themass flux in AUSM + -up is given by˙ m / = c / M / (cid:40) (cid:37) L if ˙ M / > ,(cid:37) R otherwise , (36)with the interface Mach number M / = M + (4) ( M L ) + M − (4) ( M R ) − K p f p a max (cid:16) − σ ¯ M , (cid:17) p R − p L (cid:37) / c / . (37)Here, M ± (4) are fourth degree polynomial functions and ¯ M = ( u L + u R ) / (2 c / ). K p and σ are constants between zero and unity.The last term is called the p -term and is a pressure di ff usionterm that was introduced to ensure pressure–velocity coupling atlow speeds (Edwards & Liou 1998). In the original AUSM + -upscheme the Mach number dependent scaling of the p -term, f p a ,was chosen identical to that of the u -term, f a . Similar to Miczek(2013) we chose independent cut-o ff values for f a and f p a = M po (2 − M po ) , M po = min (cid:104) , max (cid:16) M , M p cut (cid:17)(cid:105) . (38) Article number, page 6 of 24. V. F. Edelmann et al.: Well-balancing in astrophysical hydrodynamics
This way, f p a can be defined with a significantly higher cut-o ff Mach number (typically around 10 − ) compared to f a . This pre-vents stability issues, which can occur at locally very low Machnumbers when 1 / f p a becomes exceedingly large. We use thismodified scheme in the presented tests and refer to it as AUSM + -up. Still even this scheme is struggling with maintaining a simplehydrostatic stratification as shown in Section 5.1.Note that the role of the e ff ect of the pressure di ff usion is al-tered by combining AUSM + -up with any of the well-balancingmethods described in the following: Well-balancing techniqueslead to an exact reconstruction of the hydrostatic pressure. Onlythe non-hydrostatic pressure deviations are captured by the re-construction and given to the numerical flux. Hence, when well-balancing is applied, the pressure di ff usion acts on the non-hydrostatic pressure only.
4. Well-balancing methods
To illustrate the issue with configurations close to hydrostaticequilibrium in finite-volume codes we consider the e ff ect of onetime step on an initial configuration in perfect hydrostatic equi-librium. For simplicity we discretize the time derivative usingthe forward Euler method and only consider the one-dimensionalEuler equations. From the analytic result we expect that this stepwill not alter the states and any subsequent steps will also keepthe hydrostatic equilibrium intact. To this end we associate a dis-crete stationary solution that provides a good approximation tothe hydrostatic equilibrium. If starting from a discrete hydro-static equilibrium, the solution of the time evolutionary prob-lem does not change for a numerical scheme, we call it well-balanced .The density, momentum, and total energy after one time stepof length ∆ t (denoted by the superscript “1”) are calculated fromthe previous values (denoted by the superscript “0”), the inter-face fluxes ˆ F , and the cell-centered source term ˆ S . The result forcell i is (cid:37) i = (cid:37) i − ∆ t ∆ x (cid:20)(cid:18) ˆ F i + (cid:19) − (cid:18) ˆ F i − (cid:19) (cid:21) , (39)( (cid:37) u ) i = ( (cid:37) u ) i − ∆ t ∆ x (cid:20)(cid:18) ˆ F i + (cid:19) − (cid:18) ˆ F i − (cid:19) (cid:21) + ∆ t (cid:16) ˆ S i (cid:17) , (40) E i = E i − ∆ t ∆ x (cid:20)(cid:18) ˆ F i + (cid:19) − (cid:18) ˆ F i − (cid:19) (cid:21) , (41)where the indexes outside the parentheses stand for the vec-tor component of the flux or source term. To be well-balanced,a scheme needs to guarantee that the step leaves the state un-changed, which leads to the conditions0 = (cid:18) ˆ F i + (cid:19) − (cid:18) ˆ F i − (cid:19) , (42)0 = (cid:18) ˆ F i + (cid:19) − (cid:18) ˆ F i − (cid:19) + ∆ x (cid:16) ˆ S i (cid:17) , (43)0 = (cid:18) ˆ F i + (cid:19) − (cid:18) ˆ F i − (cid:19) . (44)As the fluxes are evaluated at di ff erent states, none of these con-ditions is fulfilled automatically. A consistent numerical flux willautomatically satisfy Eqs. (42) and (44) because the 1- and 3-components of the fluxes are zero for u =
0. The case of Eq. (43)is less straightforward. Here, the discretization of the source termmust be constructed to match the flux di ff erence in the hydro-static case. The methods described in the following subsectionsachieve this in di ff erent ways. One method to turn almost any hydrodynamics scheme intoa well-balanced scheme was suggested by Cargo & Le Roux(1994) (see also Le Roux 1999). The only prerequisites it re-quires is support for a general equation of state and flux func-tions that resolve contact discontinuities. For completeness wedescribe the original method, which is only applicable in theone-dimensional case, and turn to the multidimensional exten-sion in Section 4.1.2. The one-dimensional Euler equations withgravity read ∂∂ t (cid:37)(cid:37) uE (cid:48) + ∂∂ x (cid:37) u (cid:37) u + pu ( E (cid:48) + p ) = (cid:37) g (cid:37) gu , (45)with g being the, possibly negative, gravitational acceleration inthe x -direction. E (cid:48) = (cid:37)(cid:15) + (cid:37) | v | is the total energy excluding po-tential energy. It is only used in the one-dimensional method inthis subsection. The multidimensional extension in Section 4.1.2follows a slightly di ff erent principle using the total energy E in-cluding the potential as defined in Section 2.Cargo & Le Roux (1994) suggest the introduction of a po-tential q defined by its spatial and temporal derivatives ∂ q ∂ x = (cid:37) g , ∂ q ∂ t = − (cid:37) ug . (46)Numerically, this potential is treated like a composition vari-able, meaning that its time evolution is determined by the advec-tion equation, ∂ ( (cid:37) q ) ∂ t + ∂ ( (cid:37) qu ) ∂ x = . (47)This ensures that the conditions of Eq. (46) are fulfilled at alltimes if they are satisfied initially.Expressing the right side of Eq. (45) using q yields ∂∂ t (cid:37)(cid:37) uE (cid:48) + ∂∂ x (cid:37) u (cid:37) u + pu ( E (cid:48) + p ) = ∂ q ∂ x − ∂ q ∂ t . (48)Collecting derivatives with respect to the same variable and in-serting 0 = q − q results in ∂∂ t (cid:37)(cid:37) uE (cid:48) + q + ∂∂ x (cid:37) u (cid:37) u + p − qu ( E (cid:48) + q + p − q ) = . (49)At this point we introduce a modified EoS based on an arbitraryoriginal EoS by defining a new pressure, Π , and total energy pervolume, F , as Π = p − q , F (cid:48) = E (cid:48) + q . (50)With this definition Eq. (49) takes the form of the homogeneousEuler equations for a modified EoS, ∂∂ t (cid:37)(cid:37) uF (cid:48) + ∂∂ x (cid:37) u (cid:37) u + Π u ( F (cid:48) + Π ) = . (51)The physical meaning of q is that of hydrostatic pressure, ex-cept for an arbitrary constant o ff set. This means that the modifiedpressure Π of an atmosphere in perfect hydrostatic equilibrium Article number, page 7 of 24 & A proofs: manuscript no. wellbalancing is spatially constant, making the solution of the Euler equationstrivial.This new EoS does not change the speed of sound. This caneasily be seen from rewriting the expression for the speed ofsound for a general EoS (20) in terms of Π and F (cid:48) . Because q does not depend on any other thermodynamic variables, therewritten expression takes the same form as the original. An obvious multidimensional extension of the aforementionedmethod would be introduce a potential q with the defining prop-erties ∇ q = (cid:37) g , ∂∂ t q = − (cid:37) g · u . (52)While the algebra shown in Eqs. (48) to (51) is still valid forthis new potential, the problem with this definition is the mereexistence of this potential q . From Eq. (52) follows that ∇ × ( ∇ q ) = ∇ × (cid:37) g = , (53)and using the fact that ∇ × g =
0, due to g being derived from agravitational potential, we can simplify this equation to ∇ × (cid:37) g = (cid:37) ∇ × g − g × ∇ (cid:37) = − g × ∇ (cid:37) = . (54)This means that a potential only exists if the cross product g ×∇ (cid:37) vanishes. This only happens for any of the following three con-ditions: the trivial case of gravity being globally zero, the case ofconstant density, and the case where the gradient of density hasthe same equilibrium with its surroundings. Thus, the approachof Eq. (52) is not suitable.In order to construct a general, but approximate, multidi-mensional extension of the well-balanced method from Sec-tion 4.1.1, we restrict ourselves to problems in which (cid:37) does notvary strongly on equipotential surfaces of the gravitational po-tential φ . We denote an average on these equipotential surfaces,which we call horizontal average , with the operator (cid:104)·(cid:105) , so wecan define an averaged density (cid:37) = (cid:104) (cid:37) (cid:105) . (55)By definition of the horizontal average, the gradient of (cid:37) is par-allel to g . This allows us to define the potential by ∇ q = (cid:37) g , ∂∂ t q = , (56)for which Eq. (54) is fulfilled automatically.The fluxes (Eq. (3)) and the source term (Eq. (6)) in the com-pressible Euler equations (Eq. (2)) can then be rewritten as f ξ l = (cid:37) n T ξ l v (cid:37) u n T ξ l v + ( n ξ l ) x Π (cid:37) v n T ξ l v + ( n ξ l ) y Π (cid:37) w n T ξ l v + ( n ξ l ) z Π n T ξ i v ( F + Π ) , s = − ( (cid:37) − (cid:37) ) ∂φ∂ x − ( (cid:37) − (cid:37) ) ∂φ∂ y − ( (cid:37) − (cid:37) ) ∂φ∂ z . (57)There are two fundamental di ff erences between this form andEq. (51). First, F = E + q is defined including the potential en-ergy. This is di ff erent from the one-dimensional case in Eq. (50). q is now temporally constant, but because of the di ff erent defi-nition of the total energy, the source term in the energy equationstill vanishes. Second, the source term in the momentum equa-tion does not completely vanish anymore. It is now proportional to the local deviation of density from its horizontal average. Eventhough this scheme loses some of the advantageous properties ofthe one-dimensional version, it is still a significant improvementfor multidimensional simulations of stratified atmospheres.Because it only requires to a slight modification of the EoS,it is expected that the Cargo–LeRoux well-balancing method canbe implemented into an existing hydrodynamic code rather eas-ily, provided it supports a general EoS. α - β method Another approach to balance any hydrostatic solution is the α - β well-balancing method presented by Berberich et al. (2018).Similar to the Cargo–LeRoux method, the hydrostatic solutionneeds to be known a priori. An advantage of the α - β well-balancing however is, that it permits an arbitrary, even multidi-mensional, structure of the hydrostatic target solution. Hence,the α - β method is able to deal with the example of a low-density bubble in pressure-equilibrium with its surrounding ofSection 4.1 that could not be balanced with the Cargo–LeRouxmethod.The key idea of this method is to replace the physical valuesof (cid:37) and p by their respective relative deviation prior to the re-construction at the cell interfaces. The reconstructed values arethen multiplied by the known hydrostatic solution at the inter-faces. The second-order approximation of gravity is constructedsuch that the source term exactly cancels the flux over the inter-face and hence that the initial hydrostatic stratification is main-tained to machine precision.For convenience, the working principle of the α - β well-balancing method is demonstrated for the one-dimensional Eulerequations in Cartesian coordinates. The extension to higher di-mensional curvilinear grids follows the same principle and canbe found in Berberich et al. (2019) together with a mathemati-cally more rigorous formulation of the method.For a given gravitational potential φ ( x ), we denote ˜ (cid:37) and ˜ p asthe solution to the hydrostatic equation [Eq. (11)] in one dimen-sion, that is ∂ ˜ p ∂ x = − ˜ (cid:37) ∂φ ( x ) ∂ x . (58)The solutions are written as˜ (cid:37) = (cid:37) α ( x ) , ˜ p = p β ( x ) , (59)where α ( x ) , β ( x ) are dimensionless profiles and (cid:37) , p carry thephysical dimension of density and pressure, respectively. It isassumed that the profiles in Eq. (59) are known at least at coor-dinates that coincide with cell centers and interfaces.The numerical solution of the one-dimensional Euler equa-tion in the finite volume approach requires the reconstruction ofthe quantities from each cell center i of the computational do-main to the respective cell interfaces i + . In general, there issome freedom in choosing the set of quantities that is recon-structed in addition to the velocities. This is exploited by the α - β well-balancing method, which considers the relative deviationfrom the hydrostatic solution. The set of quantities at cell i thatare reconstructed is hence chosen to be W i = (cid:16) W (cid:37) i , W ui , W pi (cid:17) = (cid:32) (cid:37) i α ( x i ) , u i , p i β ( x i ) (cid:33) . (60)There is no restriction on the specific choice of the reconstruc-tion scheme that calculates the value of W i at the interface, thatis W i ± . After reconstruction, the variables are transformed back Article number, page 8 of 24. V. F. Edelmann et al.: Well-balancing in astrophysical hydrodynamics to their physical counterpart by multiplication with the knownhydrostatic solution at the interfaces, (cid:37) L / Ri + = α (cid:16) x i + (cid:17) (cid:18) W (cid:37) i + (cid:19) L / R , u L / Ri + = (cid:18) W ui + (cid:19) L / R , p L / Ri + = β (cid:16) x i + (cid:17) (cid:18) W pi + (cid:19) L / R . (61)Here, L / R denotes the value at the interface when reconstruct-ing from the left or right side, respectively. The values given byEq. (61) enter into the numerical flux function (see Section 3.2).If density (cid:37) and pressure p on the computational grid corre-spond to the hydrostatic solution Eq. (59) and u ≡
0, it followsthat the quantities reconstructed from the left and right side areequal for all cell interfaces and hence U Li + = U Ri + = U i + . (62)If the left and right interface state are the same, any numericalflux function has to equal the analytical flux function to ensurethe consistency of the method. Thus,ˆ F x (cid:18) U Li + , U Ri + (cid:19) = f x (cid:16) U i + (cid:17) = p ( x i + )0 , (63)which immediately follows from Eq. (3) for vanishing velocities.In order to maintain hydrostatic equilibrium, the residual fluxof Eq. (63) has to be balanced exactly by the source term s inEq. (2). To achieve this, the α - β method expresses the gravita-tional potential in Eq. (6) with the aid of the hydrostatic equationEq. (58) as − ∂φ∂ x = p (cid:37) α ( x ) ∂β ( x ) ∂ x . (64)The one-dimensional source term for gravity [see Eq. (22)] isthen given byˆ S i = s i , s i = p (cid:37) β i + − β i − ∆ x (cid:37) i α i , (65)which is a second-order accurate discretization. If the states onthe computational grid correspond to the hydrostatic solution,then Eq. (65) reduces to s i = p (cid:16) x i + (cid:17) − ˜ p (cid:16) x i − (cid:17) = ˆ F x (cid:18) U Li + , U Ri + (cid:19) − ˆ F x (cid:18) U Li − , U Ri − (cid:19) (66)and the discretized source term exactly cancels the interfacefluxes. This leads to zero residual and thus to a well-balancedscheme. The well-balanced property for the one-dimensional α - β method is formally shown in Berberich et al. (2018). In the following we give a short description of the simple andgeneral well-balanced method introduced in Berberich et al.(2021). For more details we refer the reader to this reference.The core of the method is the target solution ˜ u , which must be known a priori. It has to be a stationary solution to the Eulerequations (2), that is it has to satisfy the relation A ξ ∂ f ξ ( ˜ u ) ∂ξ + A η ∂ f η ( ˜ u ) ∂η + A ζ ∂ f ζ ( ˜ u ) ∂ζ = J s ( ˜ u ) . (67)It is noteworthy that, in contrast to other well-balancing meth-ods, it can include a nonzero velocity. In the numerical applica-tions in Section 5 we are going to store the hydrostatic or sta-tionary solution which shall be well-balanced in ˜ u . SubtractingEq. (67) from Eq. (2) yields the evolution equation J ∂ ( ∆ u ) ∂ t + A ξ (cid:32) ∂ f ξ ( ˜ u + ∆ u ) ∂ξ − ∂ f ξ ( ˜ u ) ∂ξ (cid:33) + A η (cid:32) ∂ f η ( ˜ u + ∆ u ) ∂η − ∂ f η ( ˜ u ) ∂η (cid:33) + A ζ (cid:32) ∂ f ζ ( ˜ u + ∆ u ) ∂ζ − ∂ f ζ ( ˜ u ) ∂ζ (cid:33) = J ( s ( ˜ u + ∆ u ) − s ( ˜ u )) (68)for the deviations ∆ u = u − ˜ u from the target solution ˜ u . Inorder to obtain a well-balanced scheme which exactly maintainsthe stationary target solution ˜ u , we discretize Eq. (68) instead ofEq. (2). This yields V i jk ∂ ( ∆ U ) i jk ∂ t = − A i + , j , k ( ˆ F dev ξ ) i + , j , k + A i − , j , k ( ˆ F dev ξ ) i − , j , k − A i , j + , k ( ˆ F dev η ) i , j + , k + A i , j − , k ( ˆ F dev η ) i , j − , k − A i , j , k + ( ˆ F dev ζ ) i , j , k + + A i , j , k − ( ˆ F dev ζ ) i , j , k − + V i jk ˆ S dev i jk , (69)where the numerical flux di ff erences (cid:16) ˆ F dev ξ (cid:17) i + , j , k = (cid:16) ˆ F ξ (cid:17) i + , j , k − f ξ (cid:104) ˜ u (cid:16) x i + , j , k (cid:17)(cid:105) , (cid:16) ˆ F dev η (cid:17) i , j + , k = (cid:16) ˆ F η (cid:17) i , j + , k − f η (cid:104) ˜ u (cid:16) x i , j + , k (cid:17)(cid:105) , (cid:16) ˆ F dev ζ (cid:17) i , j , k + = (cid:16) ˆ F ζ (cid:17) i , j , k + − f ζ (cid:104) ˜ u (cid:16) x i , j , k + (cid:17)(cid:105) are the di ff erences between the numerical fluxes evaluated at thestates ˜ u + ∆ U L / R and the exact fluxes evaluated at the values ofthe target solution at the interface centers. The interface values ∆ U L / R are obtained via reconstruction. To reconstruct in the setof variables u other = T ( u ) that is di ff erent from conserved vari-ables u , the reconstruction is applied to the transformed devia-tions ∆ U other i jk = T (cid:16) ˜ u (cid:16) x i jk (cid:17) + ∆ U i jk (cid:17) − T (cid:16) ˜ u (cid:16) x i jk (cid:17)(cid:17) . (70)After reconstruction, the interface values are transformed back.For the left interface values, this reads U Li + , j , k = T − (cid:18) T (cid:16) ˜ u (cid:16) x i + , j , k (cid:17)(cid:17) + ∆ U other i + , j , k (cid:19) , (71)and right interface values are calculated likewise. The sourceterm di ff erence discretization in Eq. (69) isˆ S dev i jk = ˆ S (cid:16) ˜ u (cid:16) x i jk (cid:17) + ∆ U i jk (cid:17) − ˆ S (cid:16) ˜ u (cid:16) x i jk (cid:17)(cid:17) , (72)where ˆ S ( U ) means that the source term discretization Eq. (22)is evaluated at the state U . It has been shown in Berberich et al.(2021) that this modification of the scheme (21) renders it well-balanced in the sense that the residual vanishes if u = ˜ u and thatthe method can be applied in arbitrarily high order finite-volumecodes. Article number, page 9 of 24 & A proofs: manuscript no. wellbalancing
5. Numerical tests and application examples
This section presents simple test simulations that verify thatthe presented well-balancing methods in combination with alow-Mach flux function are stable and do enable correct rep-resentation of slow flows in stellar-type stratifications. This isdone for steady-state and dynamical test problems. Unless statedotherwise, time integration is performed with the implicit ES-DIRK23 scheme in combination with the CFL ug time step cri-terion [Eq. (26)]. The corresponding value of c CFL is stated foreach of the simulations individually. For the numerical flux, theAUSM + -up scheme (Section 3.2.1) is used with f p a = . f a = − . Interface values are calculated from cell averagesby applying linear reconstruction without any slope or flux lim-iters. The set of reconstructed quantities is either (cid:37) - p or (cid:37) - T ,depending of the well-balancing method, and will be stated ex-plicitly for the respective simulations. The choice of appropriateboundary conditions depends on the specific setup, hence theyare given for each setup individually. All of our test cases aretwo-dimensional for computational reasons, although the meth-ods are equally valid in three spatial dimensions. Qualitativelydi ff erent results for the 3D case are not expected. A useful test problem for the quality of a well-balanced schemeis a stably stratified atmosphere. A zero-velocity initial conditionshould be maintained perfectly, at least down to rounding errors.However, in a not well-balanced scheme the pressure-gradientforce and gravity do not cancel exactly and the systems ends upin a state with small, but nonzero, acceleration. Depending on thedetails of the numerical flux, this acceleration may prevent thesimulation of flows at very low Mach numbers or it may causeunphysical convection in formally stable stratifications. In anycase, keeping a hydrostatic stratification stable is a necessary,but not a su ffi cient condition for a well-balanced scheme to beuseful for the simulation of stratified atmospheres.We start with the one-dimensional (1D), isothermal test caseof Käppeli & Mishra (2016) on the domain [0 , x -direction and the gravitational potential isgiven by φ = s x , (73)with a steepness s of 1 cm s − . The constant temperature profileis convectively stable for any equation of state with a positivevalue of δ [Eq. (14)]. The density and pressure profiles fulfillingEq. (11) are given by, (cid:37) = (cid:37) exp (cid:32) − (cid:37) p φ (cid:33) , p = p exp (cid:32) − (cid:37) p φ (cid:33) . (74)This expression holds for any form of the gravitational potential,not just Eq. (73). The reference density and pressure are set to (cid:37) = − and p = γ = / ff erent flux functions. The simulations use linear reconstruc-tion in (cid:37) - p variables, c CFL = . t = t BV . The figure shows that runs without t / t BV − − − − − − m a x ( M a ) isothermal T profile, % - p linear reconstruction Roe CL WBRoe α - β WBRoe Deviation WBRoe no WB AUSM + -up CL WBAUSM + -up α - β WBAUSM + -up Deviation WBAUSM + -up no WB0 2000 4000 6000 8000 10000 12000 14000 16000 t / t SC Fig. 1.
Time evolution of the maximum Mach number in a 1D atmo-sphere with an isothermal temperature profile [Eq. (74)] and a lineargravitational potential [Eq. (73)]. The colors indicate di ff erent well-balancing (WB) methods, the markers di ff erent flux functions. CLstands for Cargo–LeRoux. Time is given in units of the Brunt–Väisälätime t BV = .
93 s [Eq. (18)] and sound-crossing time t SC = .
10 s[Eq. (19)]. The curves have been slightly smoothed for better visibil-ity. t / t BV − − − − − − m a x ( M a ) a nd m a x ( ∆ % ) isothermal T profile, % - p linear reconstruction Roe CL WBRoe α - β WBRoe Deviation WBRoe no WB AUSM + -up CL WBAUSM + -up α - β WBAUSM + -up Deviation WBAUSM + -up no WB0 2000 4000 6000 8000 10000 12000 14000 16000 t / t SC Fig. 2.
Same as Fig. 1, but for a 2D atmosphere. The solid lines repre-sent the maximum Mach numbers on the grid and the dashed lines thehorizontal density fluctuations according to Eq. (75). The curves havebeen slightly smoothed for better visibility. well-balancing immediately reach Mach numbers between 10 − and 10 − , while all of the tested well-balancing methods manageto keep the Mach number below 10 − . The choice of flux func-tion does not qualitatively a ff ect the results here, in particular itdoes not play a role if we use a low Mach number flux, such asAUSM + -up, or a standard flux, such as Roe’s approximate Rie-mann solver (e.g., Toro 2009). Article number, page 10 of 24. V. F. Edelmann et al.: Well-balancing in astrophysical hydrodynamics
Certain phenomena, such as convection, are inherently mul-tidimensional and cannot develop in 1D geometry. Therefore werepeat this test using a two-dimensional (2D) atmosphere. Forsimplicity we keep gravity pointing in the negative x -direction.The horizontal boundaries are periodic. Figure 2 shows the evo-lution of the maximum Mach number and horizontal densityfluctuation ∆ (cid:37) over time. The latter is defined as ∆ (cid:37) = (cid:37) − (cid:104) (cid:37) (cid:105) y , (75)with (cid:104)·(cid:105) y denoting an average over the y direction.Similar to the 1D case, we see that the simulations withoutwell-balancing immediately reach Mach numbers around 10 − .All simulations using any well-balanced method start from verylow Mach numbers ( ≈ − ), but the ones using the low Machnumber AUSM + -up flux show an exponential growth over thenext few thousand t BV . The growth of the Mach number is linkedto that of ∆ (cid:37) . This growth even a ff ects the Mach numbers inthe simulation without well-balancing, where the Mach numberincreases further after about 4000 t BV . The well-balanced sim-ulations using the standard Roe flux do not show this behaviorand retain the very low Mach numbers. We found that other lowMach number fluxes, such as the one by Miczek et al. (2015) andLi & Gu (2008), show a similar growth although the rate variesbetween the di ff erent schemes. While the standard flux seems toperform better in this static test case, it is ultimately not suitedfor the simulation of dynamic phenomena at low Mach numbers,such as convection or waves, due to its high numerical dissipa-tion.Figure 3 shows the typical pattern of the exponentially grow-ing perturbation in both Mach number and ∆ (cid:37) . In both quantitieswe see a resolved pattern in the horizontal direction, but a grid-level oscillation in the vertical direction. One hypothesis for thisbehavior is that it is due to unresolved internal gravity waves, seeSection 5.3.We run the same tests for polytropic stratifications. Here theprofiles of density, pressure, and temperature are (cid:37) = (cid:37) θ ν − , p = p θ νν − , T = p µ(cid:37) R θ = T θ, (76)with θ = − ν − ν (cid:37) p φ. (77)We set µ = − . The polytropic index ν determines theslopes of the profiles. If ν is less than the adiabatic exponent γ ,the atmosphere is stable. If ν > γ , it is unstable. If both are equalthe atmosphere is isentropic and therefore marginally stable.Figure 4 shows the maximum Mach number and ∆ (cid:37) for anisentropic atmosphere. Here the Brunt–Väisälä time t BV is notwell defined because N =
0. Instead we use the sound-crossingtime t SC = .
28 s for reference. In contrast to the isothermal testwe see that Cargo–LeRoux well-balancing behaves quite similarto the not well-balanced case. Both quickly reach Mach numbersof about 10 − with flow patterns that resemble 2D convection. α - β and Deviation well-balancing stay at very low Mach num-bers below 10 − , independently of the choice of flux function.This result suggests that the issue with the exponential growth ofperturbations in combination with low Mach number flux func-tions is more severe in more stable stratifications.As a last test we run a polytropic stratification with ν = . < γ = /
3, which is stably stratified, but less so than theisothermal case. The relevant time scales here are t BV = . t SC = .
13 s. The results are shown in Fig. 5. Here we see . . . . . . . . . x / c m .
00 0 .
25 0 .
50 0 .
75 1 .
00 1 .
25 1 .
50 1 .
75 2 . y / cm0 . . . . . . . . . x / c m − − M a × − − . − . . . . ∆ % × − Fig. 3.
Mach number ( top panel ) and horizontal density fluctuation ∆ (cid:37) ( bottom panel ) of the 2D isothermal simulation using the AUSM + -upflux and Deviation well-balancing at time t = t BV . Gravity is di-rected in negative x -direction. that the cases without well-balancing and with Cargo–LeRouxwell-balancing fail to preserve the hydrostatic equilibrium incombination with the low Mach number flux. For the other twowell-balancing methods there is a much clearer di ff erence in thegrowth rate of the perturbation with Deviation well-balancinggrowing significantly more slowly, reaching only Mach numbersof 10 − after 2500 t BV .The spurious growth we found in combination with low-Mach-number fluxes is not necessarily a problem in real-worldapplications. We considered very long time scales of thousandsof sound-crossing and Brunt–Väisälä times. This is much longerthan the time scales typically used to test well-balancing meth-ods. Käppeli & Mishra (2016), for example, ran their hydrostaticsetup only for 2 t SC . While this would definitely show any ma-jor issues with the basic well-balanced property of the scheme,it does not reveal the long-term growth of instabilities in the sta-ble region. This is something to bear in mind when applyingsuch methods to partly convectively stable configurations, suchas stars. Article number, page 11 of 24 & A proofs: manuscript no. wellbalancing t / t SC − − − − − − − m a x ( M a ) a nd m a x ( ∆ % ) isentropic profile, % - p linear reconstruction Roe CL WBRoe α - β WBRoe Deviation WBRoe no WB AUSM + -up CL WBAUSM + -up α - β WBAUSM + -up Deviation WBAUSM + -up no WB Fig. 4.
Same as Fig. 2, but for an isentropic stratification. The solidlines represent the maximum Mach numbers on the grid and dashedlines the horizontal density fluctuations according to Eq. (75). Time isgiven in units of the sound-crossing time t SC = .
28 s. The curves havebeen slightly smoothed for better visibility. The Roe solver runs withoutwell-balancing and Cargo–LeRoux well-balancing experienced convec-tion at high Mach numbers close to the upper boundary causing rapiddivergence leading to the stop of the simulation. t / t BV − − − − − − m a x ( M a ) a nd m a x ( ∆ % ) polytropic ( ν = .
6) profile, % - p linear reconstruction Roe CL WBRoe α - β WBRoe Deviation WBRoe no WB AUSM + -up CL WBAUSM + -up α - β WBAUSM + -up Deviation WBAUSM + -up no WB0 5000 10000 15000 20000 t / t SC Fig. 5.
Same as Fig. 2, but for a polytropic stratification with ν = . γ = /
3. The solid lines represent the max-imum Mach numbers on the grid and dashed lines the horizontal den-sity fluctuations according to Eq. (75). Time is given in units of Brunt–Väisälä time t BV = . t SC = .
13 s. Thecurves have been slightly smoothed for better visibility.
Convection in the stellar interior is usually slow and almost adi-abatic. A typical convection zone is nearly isentropic, althoughthe stratification in the star’s gravitational field can span orders of magnitude in pressure and density. The buoyant accelerationof a fluid parcel is given by its entropy fluctuation with respectto the mean entropy at any given radius. Entropy fluctuations areconstantly created by sources of heating or cooling. However,once the fluid parcel has left the heating or cooling layer andtravels through the rest of the convection zone, the parcel’s en-tropy must be preserved except for those parts that have mixedwith their surroundings. If the flow is slow, all entropy fluctu-ations inside the convection zone are also small and it becomesnumerically challenging to preserve their exact values when den-sity and pressure change by large factors as the fluctuations areadvected along the gravity vector.We test the numerical schemes’ entropy-preservation prop-erties under the conditions described above by simulating thebuoyant rise of a “bubble” with an adjustable initial entropy fluc-tuation embedded in a layer of constant entropy. The layer isstrongly stratified in pressure and density due to the presence ofgravity. We call this setup “hot bubble”, because we use positiveinitial entropy fluctuations. Negative initial entropy fluctuationswould make the bubble fall, but everything else would work inthe same way. Our setup is similar to that used by Almgren et al.(2006) to test their MAESTRO code, although their equation ofstate and stratification di ff er from our setup. We also test ourmethods down to much lower Mach numbers.We construct the stratification in a two-dimensional box10 cm wide and spanning the vertical range from y = y = y max = . × cm. To avoid any influence of boundaryconditions, we use periodic boundaries and a periodic profile ofgravitational acceleration g y = g sin( k y y ) , (78)where g is a constant to be specified later and k y = π/ y max . Thebox is filled with an ideal gas with the ratio of specific heats γ = / µ = − . At y =
0, we setthe pressure to p = Ba and the temperature to T =
300 K.The stratification is isentropic, that is (cid:37) ( y ) = (cid:32) p ( y ) A (cid:33) /γ , (79)where A = A = const . everywhere outside of the bubble. Insidethe bubble, we perturb the entropy via A = A + (cid:32) ∆ AA (cid:33) t = cos (cid:32) π rr (cid:33) , (80)where ( ∆ A / A ) t = is the bubble’s amplitude and r = (cid:104) ( x − x ) + ( y − y ) (cid:105) / is the distance from the bub-ble’s center with x = × cm and y = . × cm. Thebubble has a radius of r = . × cm. We do not perturbthe hydrostatic pressure stratification, which is given by theexpression p ( y ) = p − γ + (cid:32) − γ (cid:33) g A γ k y (cid:104) − cos( k y y ) (cid:105) γγ − . (81)The amplitude g of the gravity profile sets the ratio ofthe maximum to the minimum pressure in the periodic pres-sure profile. To make the problem numerically challenging,we use a pressure ratio of 100, which is achieved with g = − . × cm s − . This stratification is strongerthan that in the convective core of a typical massive main-sequence star, in which the pressure changes by a factor of a few. Article number, page 12 of 24. V. F. Edelmann et al.: Well-balancing in astrophysical hydrodynamics
On the other hand, the relative pressure drop from the bottom ofthe solar convective envelope to the photosphere is ≈ .We start with a moderate initial entropy perturbation of( ∆ A / A ) t = = − , which makes the bubble rise at moderatelylow Mach numbers of a few times 10 − . This allows us to per-form simulations with all three well-balancing methods as wellas simulations without any well-balancing at modest grid reso-lution of 128 × . α - β methods, which require (cid:37) - p reconstruction, wetest both (cid:37) - p and (cid:37) - T reconstruction. The central, most buoyant,part of the bubble accelerates fastest. The bubble gets deformedinto a mushroom-like shape with two trailing vortices and it ex-pands as it rises into layers of lower pressure. Ideally, the initiallypositive entropy fluctuations ∆ A / A should mix with the isen-tropic ( ∆ A / A =
0) background stratification, creating smallerbut still positive entropy fluctuations. The entropy fluctuationsmay locally increase a bit as kinetic energy is slowly dissipatedinto heat, but there is no physical way for them to become nega-tive. Any negative entropy fluctuations in the numerical solutionresult from numerical errors.Figure 6 shows that both the absence of well-balancing andthe Cargo–LeRoux method generate large areas of negative en-tropy fluctuations comparable to or even larger in absolute valuethan the bubble’s initial amplitude. Large-scale positive entropyfluctuations also occur far from the bubble and they clearly donot result from hydrodynamic mixing. They rather seem to becaused by entropy nonconservation as the bubble pushes the sur-rounding stratification upwards and downwards. In the no-well-balancing case, errors in ∆ A / A are smaller when (cid:37) - T reconstruc-tion is employed as compared with (cid:37) - p reconstruction. This maybe due to the fact that the pressure changes by a factor of 100 inthe computational domain whereas the temperature only changesby a factor of 6 .
3. In any case, α - β and Deviation well-balancingclearly provide far superior results with only a mild and highlylocalized undershoot in ∆ A / A above the bubble. With the Devi-ation method, this success is independent of the choice of recon-struction.All of the methods tested converge upon grid refinement,although not in the same way, see the series of runs shown inFig. 7. In this series, we keep the 0 . × . α - β andDeviation methods. The errors are still substantial even on thefinest (256 × α - β and Deviation methods does not decreasein amplitude but it does decrease in spatial extent as the gridis refined. Surprisingly, we do not observe any significant en-tropy fluctuations far from the bubble even on the coarsest grid(64 ×
96) with these two methods.As the bubble’s initial amplitude ( ∆ A / A ) t = is decreased, thetypical Mach number in the flow field decreases asMa ∝ (cid:32) ∆ AA (cid:33) / t = . (82)This scaling results from the fact that the bubble’s acceleration isproportional to ∆ (cid:37)/(cid:37) ∝ ∆ A / A and the velocity an object in uni-formly accelerated motion reaches over a fixed distance (i.e. untilthe bubble has reached the same evolutionary stage) is propor-tional to the square root of the acceleration. Although the bub-ble’s acceleration is not constant, the bubble always evolves in the same way, just slower, when its initial amplitude is decreasedand the scaling still holds.The scaling is also demonstrated in Fig. 8, which shows aseries of runs performed on a 128 ×
192 grid with ( ∆ A / A ) t = ranging from 10 − down to 10 − . As ( ∆ A / A ) t = decreases, weincrease time steps in this order: 0 . α - β and Deviation well-balancing methods in this experiment, be-cause entropy-conservation errors quickly dominate the solutionwhen the initial amplitude is decreased and the Cargo–LeRouxor no well-balancing method is used. Because the Mach numberis expected to scale according to Eq. 82, we scale the time whenthe simulation is stopped with ( ∆ A / A ) − / t = . This way, we com-pare the results when the bubble has reached approximately thesame height and evolutionary stage as the previously discussedcase with ( ∆ A / A ) t = = − at t =
300 s (Fig. 6). Both the α - β and Deviation well-balancing methods provide essentially thesame solutions, reproducing the expected Mach-number scalingdown to Ma ∼ − . The amplitude of the entropy undershootabove the bubble is 24% lower when the α - β method is usedas compared with the Deviation method. The most extreme runwith ( ∆ A / A ) t = = − reaches the limits of our current im-plementation and we can see numerical noise developing in thestratification, see Figs. A.1 and A.2. Figure 8 also shows that theminimum and maximum entropy fluctuations in the evolved flowscale in proportion to the initial amplitude of the bubble. The previous section demonstrated the capabilities of the α - β and Deviation well-balancing schemes to evolve the rise of abubble with an entropy excess in an isentropic stratification atlow Mach numbers. This can be interpreted as a test for the fun-damental mechanism of convection. The purpose of this sectionis to also assess the benefit from well-balancing techniques forfully developed convection in a realistic stellar scenario.The prototype of a convective region in stellar interiors in-cludes a steady heating source in the form of nuclear burningthat injects energy over a long period in time. If radiative trans-port of energy is not e ffi cient enough, the temperature gradientwill steepen until it reaches the adiabatic temperature gradient.According to Eq. (16), this region will become unstable and con-vection will set in. As convection is very e ffi cient in transportingenergy, a common assumption is that the stratification settles toan nearly adiabatic temperature gradient in the steady state.In stars, a convective region is typically adjacent to con-vectively stable regions. The mixing processes across the in-terfaces of convectively stable and unstable regions have a pro-found impact on stellar evolution, yet it is particularly di ffi cult toparametrize these processes in one-dimensional stellar evolutioncodes. A common strategy for improving the current 1D descrip-tion is to investigate the dynamics at the interfaces of convectionzones by means of multidimensional hydrodynamical simula-tions (Jones et al. 2017; Cristini et al. 2019; Pratt et al. 2020;Higl, J. et al. 2021; Horst et al. 2021).The Mach number of convection deep in the stellar interior isestimated to be of the order of 10 − in early stages of stellar evo-lution. Accurate modeling of the early phases is, however, cru-cial as it determines the whole subsequent evolution of the starand inaccuracies will propagate to later stages. Thus, numericalexperiments that address convection in stellar interiors rely onschemes that accurately maintain hydrostatic stratifications and Article number, page 13 of 24 & A proofs: manuscript no. wellbalancing t = s min = -4.99e-04max = Cargo-LeRoux WB % - p rec. min = -4.41e-05max = α - β WB % - p rec. min = -5.14e-05max = Deviation WB % - p rec. min = -4.95e-05max = Deviation WB % - T rec. min = -4.78e-04max = no WB % - p rec. min = -6.37e-05max = no WB % - T rec. t = s min = -2.01e-03max = = -8.35e-05max = = -1.10e-04max = = -1.03e-04max = = -1.64e-03max = = -3.92e-04max = Fig. 6.
Time evolution (top to bottom) of the hot-bubble problem when solved using di ff erent well-balancing and reconstruction schemes (leftto right) on a 128 ×
192 grid. In all of the cases, entropy fluctuations ∆ A / A are shown on the same color scale ranging from − − (dark blue)through 0 (white) to 10 − (dark red). The minimum and maximum values of ∆ A / A in the whole simulation box are given in each panel’s inset. Theamplitude of the initial entropy perturbation is ( ∆ A / A ) t = = − . that are able to follow convection at low Mach numbers for asu ffi cient amount of time.The initial stratification for the test series presented in thissection consists of a convective region with a stable layer on top.One possibility to set up this configuration would be to use re-alistic initial conditions from a 1D stellar evolution code. How-ever, such 1D profiles usually require some extra treatment be-fore they can be used for hydrodynamical simulations, for exam-ple a smoothing of sharp composition gradients or to properlyimpose a flat entropy profile in the convection zone. Instead, weuse analytical initial conditions to be able to test the numericalmethods under well-defined but realistic conditions. This way,any numerical artifacts that may arise can solely be attributed tothe methods applied rather than to inadequate initial conditions.To construct the initial hydrostatic stratification, we followthe procedure described by Edelmann et al. (2017). It imposesthe profile of the superadiabaticity ∆ ∇ = ∇ − ∇ ad while integrat-ing the equation of hydrostatic equilibrium [Eq. (11)]. Accord-ing to Eq. (16), ∆ ∇ determines the sign of the Brunt–Väisälä fre-quency. It is therefore possible to precisely control which regionsare convectively stable ( ∆ ∇ <
0) and unstable ( ∆ ∇ >
0) as wellas the respective transitions between these regions. A marginallystable stratification is imposed inside the convection zone by set-ting ∆ ∇ CZ =
0. In the stable region, we impose ∆ ∇ SZ = −∇ ad ,which corresponds to an isothermal stratification.To connect the two regions, we use a sinusoidal transition,which ensures that the transition between the two ∆ ∇ -values iswell-defined and can be resolved numerically. The profile witha transition between the convection and stable zones then takesthe form ∆ ∇ ( y ) = ∆ ∇ CZ + (cid:20) + sin (cid:18) π η ( y , K , y CZ,SZ ) (cid:19)(cid:21) ( ∆ ∇ SZ − ∆ ∇ CZ ) , (83) with η ( y , K , y i ) = − K ( y − y i ) < − , K ( y − y i ) > , K ( y − y i ) otherwise, (84)where the constant K determines the steepness of the transition.It is set such that the transition is resolved by at least 20 gridcells. The coordinate y CZ , SZ denotes the middle of the transi-tion starting at y CZ , SZ − / K and ending at y CZ , SZ + / K . Thevalue of y CZ , SZ is given in Table 2. The computational domainspans 2 y CZ , SZ in the horizontal direction. The profiles of temper-ature T , pressure p and density (cid:37) then follow from integratingthe equation of hydrostatic equilibrium Eq. (11) as described byEdelmann et al. (2017) with a spatially constant gravitational ac-celeration of | g | = . × cm s − . The initial values requiredfor the integration are listed in Table 2. They are representative ofthe conditions expected in the convective core of a 25 M (cid:12) mainsequence star with a mean atomic weight of A = . Z = . ∆ ∇ , T , and (cid:37) as wellas the fact that the transition between the convective and stablezone is well resolved even on the coarsest SLH grid. Becausecores of massive stars are hot, a considerable fraction of the totalpressure is contributed by the radiation field. The relative impor-tance of radiation pressure is given in the lower panel of Fig. 9. Itranges from roughly 20 % in the bottom region to about 60 % inthe top region where the density is low. This test setup is there-fore also an example where the ideal gas EoS is not su ffi cientto describe the thermodynamic behavior of the gas and whichrequires well-balancing methods that can handle general EoS.To trigger convection in the initially marginally stable con-vection zone, a heat source is placed at the bottom of the convec-tion zone that continuously injects energy into the system with Article number, page 14 of 24. V. F. Edelmann et al.: Well-balancing in astrophysical hydrodynamics × min = -1.45e-02max = Cargo-LeRoux WB % - p rec. min = -6.48e-05max = α - β WB % - p rec. min = -1.03e-04max = Deviation WB % - p rec. min = -9.16e-05max = Deviation WB % - T rec. min = -6.56e-03max = no WB % - p rec. min = -1.39e-03max = no WB % - T rec. × min = -2.01e-03max = = -8.35e-05max = = -1.10e-04max = = -1.03e-04max = = -1.64e-03max = = -3.92e-04max = × min = -8.11e-04max = = -9.19e-05max = = -1.09e-04max = = -1.05e-04max = = -6.55e-04max = = -1.74e-04max = Fig. 7.
Same as Fig. 6, but showing the solutions’ resolution dependence (top to bottom) at t =
300 s. a sinusoidal profile peaking at the bottom of the domain. Theheating profile is given by˙ e h ( y ) = ˙ e a sin (cid:18) π (cid:2) + η ( y , K , (cid:3)(cid:19) erg s − cm − , (85)with a = π K ∆ y / π K ∆ y , K = . y CZ , SZ , (86)where ∆ y denotes the grid spacing in the vertical direction. Asintroduced here, ˙ e is dimensionless. The factor a ensures thatthe total heating rate is independent of grid resolution. In a sta-tionary state, the convective velocity is expected to scale withthe heating rate,Ma rms ∝ ˙ e / . (87)This scaling can be motivated using dimensional arguments(e.g., see Jones et al. 2017) or the mixing-length theory (Kip-penhahn et al. 2012) and has been confirmed in numerical exper-iments (e.g., Cristini et al. 2019; Edelmann et al. 2019; Andrassyet al. 2020).To test the low-Mach capabilities of the methods that are pre-sented in this paper, the amplitude of the heating factor ˙ e is de-creased successively while the rms Mach number measured from Table 2.
Parameters of the convection test setup.
Quantity Value T . × K p . × g cm − s − | g | . × cm s − y bot , y top
0, 2 . × cm y CZ , SZ .
66 y top x right − x left y CZ , SZ ∆ ∇ CZ , ∆ ∇ SZ , −∇ ad A , Z .
3, 1 . e is chosen such that the resulting Mach numberis in a regime where also simulations without well-balancing fol-low the scaling. A deviation from the expected scaling at lowervalues of ˙ e and correspondingly lower Mach number then in-dicates that numerical errors have become significant and themethod has reached its limit of applicability. The scaling testis performed for all available well-balancing methods and alsoin the absence of well-balancing. The domain is discretized by72 ×
144 cells. This rather coarse resolution is chosen to assessthe ability of the well-balancing method to balance hydrostatic
Article number, page 15 of 24 & A proofs: manuscript no. wellbalancing − − − − − − ( ∆ A / A ) t = − − − − − − − − min( ∆ A / A )max( ∆ A / A ) ∝ ( ∆ A / A ) t = max(Ma) ∝ ( ∆ A / A ) / t = Fig. 8.
Dependence of the maximum Mach number Ma and minimumand maximum entropy fluctuations ∆ A / A on the bubble’s initial am-plitude ( ∆ A / A ) t = . All measurements are taken when the bubble hasevolved to a stage comparable to that shown in the bottom row ofFig. 6. Open and filled symbols correspond to runs with α - β and Devia-tion well-balancing, respectively. The lines show that the Mach numberscales according to Eq. (82) and the minimum and maximum entropyfluctuations scale in proportion to the bubble’s initial amplitude. − . − . − . . ∆ ∇ ˙ e s , arb. units0 . . . . . y / cm × % / g c m − . . . p r a d / p t o t Fig. 9.
Initial stratification of the convection setup. The red curve illus-trates the position and shape of the energy injection which has non-zerovalues only at the bottom of the convection zone. Its actual amplitude isset for the di ff erent simulations individually. Dots denote the positionsof cell centers on a grid with 144 vertical cells, the lowest resolutionused in the SLH simulations presented in this section. stratifications even with a moderate number of cells. While anyconsistent method will ultimately be able to follow low-Mach-number flows given a su ffi ciently fine grid, this is computation-ally not feasible, especially in 3D. The simulations use periodicboundary conditions in the horizontal direction. At the top andbottom of the domain solid-wall boundaries are employed thatdo not allow mass to enter or leave the domain. For the timestepping, we set c CFL = .
5. The initial stratification is perturbedby random density fluctuations at the O (cid:16) − (cid:17) level to facilitatethe growth of the convective instability. Additionally, the heatingprofile Eq. (85) is modulated by a sinusoidal function along thehorizontal direction to break the initial horizontal symmetry. Thewavelength is one fourth of the horizontal extent and the ampli-tude is set to 0 .
01 ˙ e h ( y ). The modulation is switched o ff after theflow has reached a substantial fraction of its final speed.For all simulations, a grid file is saved every 200 time steps.For each saved grid, the mass-weighted rms Mach number Ma rms heating rate ˙ e − − M a r m s Ma rms ∝ ˙ e / Deviation WB % - T rec. α - β WB % - p rec.Cargo-LeRoux WB % - p rec.no WB % - p rec.no WB % - T rec. Fig. 10.
Root mean square Mach number Ma rms as a function of theheating rate ˙ e . The dashed line represents the expected scaling ac-cording to Eq. (87). All values correspond to time averages over t ∈ [ t ( N τ = , t ( N τ = rms . is calculated. To this end, only the fixed region from y bot to y CZ , SZ is considered, although the position of the boundary betweenconvective and stable zone is dynamical and may change overtime. However, for the simulation presented in this section, thise ff ect is negligible and a fixed region is chosen for convenience.For all simulations, Ma rms is then averaged over the same timespan in terms of the convective turnover time τ conv , which wedefine as τ conv = y CZ , SZ v rms , (88)where v rms is the rms velocity. This ensures that the stochasticfluctuations, which have di ff erent typical timescales for di ff erentflow speeds, are accounted for in a similar way for all simu-lations. Due to the steady heat injection and the small amountof numerical dissipation in 2D simulations, the value of τ conv slightly decreases over time as velocities slightly increase. Thus,the average of τ conv depends on the time interval considered andis not clear how to chose the proper time interval for the di ff er-ent simulations. Instead, the number of turnover times N τ as afunction of physical time t , N τ ( t ) = (cid:90) t τ conv ( t (cid:48) ) d t (cid:48) , (89)is used to determine the respective time intervals and the aver-aging is done in the time interval t ∈ [ t ( N τ = , t ( N τ = rms as functions of the heating rate are depictedin Fig. 10. For the two highest values of ˙ e , that is at Ma rms around 9 × − and 4 × − , all methods agree and Ma rms fol-lows the scaling given by Eq. (87). This is not the case at lowerheating rates. For ˙ e (cid:46) , the simulation using the Cargo–LeRoux well-balancing method considerably deviates from theexpected scaling by giving a value of Ma rms that does not cor-relate with the energy input anymore but stays rather constantlyslightly below 4 × − . Almost identical results are found whenno well-balancing is applied in combination with (cid:37) - p as recon-struction variables. At this point, it seems that Cargo–LeRouxwell-balancing is not able to improve the behavior at lower Machnumbers. The results slightly improve if (cid:37) - T are used for recon-struction and no well-balancing is used. In this case, the Machnumber settles slightly below Ma rms = × − for an energyinput rate lower than ˙ e = . The reason for this could be the Article number, page 16 of 24. V. F. Edelmann et al.: Well-balancing in astrophysical hydrodynamics non-negligible contribution of radiation-pressure to total pres-sure. As p rad ∝ T , a more accurate reconstruction of T couldlead to an interface pressure that is closer to the hydrostatic so-lution and hence artifacts from imperfect balancing are reduced.By definition, α - β well-balancing requires (cid:37) - p to be recon-structed while the variables can be chosen freely for Devia-tion well-balancing. Hence we have chosen (cid:37) - T for the Devi-ation runs for comparison. However, no major di ff erences canbe seen between α - β and Deviation in Fig. 10. Both Ma rms pro-files closely follow the scaling down to the smallest value of ˙ e .Due to the hydrostatic solution’s being stored at cell interfacesin these two methods, the particular choice of variables for re-construction seem to be less important. In contrast, the Cargo–LeRoux method reconstructs the potential q at the interfaces us-ing values at the cell centers, which can introduce an error in thetotal energy over many time steps. A possible fix for this wouldbe to store q at the cell interfaces, however this is currently notimplemented in our code.To also add a qualitative visual verification of the arising con-vection, the flow patterns are shown in Fig. 11 in the middle ofthe time frame considered. For the highest heating rate, all sim-ulations show the typical flow morphology of two-dimensionalconvection: A pair of large eddies form with a size that is deter-mined by the vertical extent of the convection zone. At lowerheating rates, the flow pattern remains basically the same forDeviation and α - β well-balancing, which strengthens our con-fidence in these solutions. In contrast, the flow patterns ob-tained with Cargo–LeRoux well-balancing or without any well-balancing are dominated by incoherent, small-scale structures.While for the “no WB” run with (cid:37) - T reconstruction, the Machnumbers achieved are somewhat smaller, the general flow pat-tern is similar. We assume that these motions are caused by theimperfect hydrostatic equilibrium.Convective regions are known to excite internal gravitywaves (IGW) in adjacent stable layers (see e.g., Rogers et al.2013; Edelmann et al. 2019; Horst et al. 2020). Assuming thatIGW are predominantly generated at periods close to the convec-tive turnover time τ conv (c.f. Edelmann et al. 2019), we can es-timate their vertical wavelength λ v using the dispersion relationof IGW in the Boussinesq approximation (see, e.g., Sutherland2010), λ v = λ h (cid:112) [ N τ conv / (2 π )] − , (90)where λ h is the horizontal wavelength. It follows that the verti-cal wavelength decreases with increasing τ conv for λ h = const . and the waves will become unresolved when the heating rate isdecreased on the same computational grid.We estimate the vertical wavelengths according to Eq. (90)exemplarily for the simulations with Deviation well-balancingshown in Figs. 10 and 11. To calculate λ v , the value of theBrunt–Väisälä frequency is taken at the top of the box domainand λ h = y CZ , SZ , which corresponds to the maximal horizon-tal wavelength that fits into the computational domain. The ratioof λ v to the vertical grid spacing is shown in Fig. 12 as blackcrosses. At all but the two highest heating rates, the verticalwavelength is less than two cells and it follows from the Nyquistsampling criterion that such waves cannot be represented on ourcoarse grid. Indeed, for these runs strange patterns appear in thestable zone as can be seen in Fig. 11. A peculiar pattern is visiblein stable zone for Deviation well-balancing at the lowest heatingrate. While its origin is not completely understood it looks simi-lar to the other runs for (cid:37) - p reconstruction or a higher grid reso-lution. Our estimate gives typical wavelengths of two to six cells in the two simulations with the highest heating rates. Hence, theflow patterns in Fig. 11 for these two simulation probably in-clude some real internal gravity waves, but they are still domi-nated by artifacts.To confirm that our interpretation is correct, we run an-other series of simulations with increasing vertical resolution at˙ e = . Deviation well-balancing is used in this experiment.Because of the higher computational costs, there is not enoughdata to perform meaningful averages and Ma rms is extracted for asingle snapshot as soon as convection has developed in the wholeconvective zone. The results are shown as blue crosses in Fig. 12.Because Ma rms is measured at earlier times compared with thecorresponding simulation in the heating series, the Ma rms of thelowest-resolution run does not coincide with the (red circled)black cross corresponding to the same simulation. We see thatthe expected typical vertical wavelength is resolved by more thaneight cells on the finest grid. Fig. 13 shows the flow pattern forthe di ff erent grid resolutions. In the center and right panel of thelower row, the resolution in the convective zone corresponds tothe resolution of the upper right panel (288 × /
32 of the starting resolution. This saves computingresources and also illustrates that the patterns do not depend onthe resolution in the convective zone. The transition is shown inFig. B.1.Our resolution study indicates that with increasing verticalresolution in the stable zone artifacts are diminished. For thehighest resolution fine wave patterns are observed. As grid res-olution increases, nearly horizontal wave patterns first appearclose to the convective boundary, where N gradually increasesfrom zero to relatively large values higher up in the stable zone.While our findings can be explained by unresolved IGW, wecannot exclude that at least to some extent also numerical arti-facts contribute to the flows in the stable zone. Nevertheless, ourtests with the simple convective box show that, as soon as we re-solve the IGW su ffi ciently, artifacts tend to disappear and any in-stabilities possibly still present do not visibly dominate the flow.At the same time, this illustrates that, depending on the actualstellar profile, rather high grid resolution is needed to properlyresolve the waves. Some astrophysical problems involve stationary solutions thatare not at rest, for example a rotating star that is partially stabi-lized by the centrifugal force (e.g., Tassoul 2000; Maeder 2009).Another case is the Keplerian motion around a central gravita-tional mass m in its gravitational potential φ ( r ) = − Gmr . Gaburroet al. (2018) describe a non-dimensional test setup of a circulardisk with (cid:37) = p = u ( x , y ) = − sin (cid:2) atan2 ( y , x ) (cid:3)(cid:115) Gmr ( x , y ) (91) v ( x , y ) = cos (cid:2) atan2 ( y , x ) (cid:3)(cid:115) Gmr ( x , y ) , (92)where r = (cid:112) x + y and atan2 ( y , x ) is the typical shortcut forchoosing the quadrant of arctan ( y / x ) correctly. For conveniencewe set G = m =
1. We simulate the Keplerian disk from ra-dius r = Article number, page 17 of 24 & A proofs: manuscript no. wellbalancing ˙ e = D e v i a ti on W B % - T r ec . Ma = . × − ˙ e = = . × − ˙ e = = . × − ˙ e = = . × − ˙ e = = . × − ˙ e = = . × − α - β W B % - p r ec . Ma = . × − = . × − = . × − = . × − = . × − = . × − C a r go - L e R oux W B % - p r ec . Ma = . × − = . × − = . × − = . × − = . × − = . × − no W B % - p r ec . Ma = . × − = . × − = . × − = . × − = . × − = . × − no W B % - T r ec . Ma = . × − = . × − = . × − = . × − = . × − = . × − − − − − − M a Fig. 11.
Mach number of the flow for di ff erent values of the heating rate ˙ e (left to right) and di ff erent well-balancing methods (top to bottom) .The dashed black lines denote the boundaries of the convection zone at y = y top , see Table 2 for an overview of simulation parameters. The insetsdisplays the rms Mach number Ma rms within the convection zone for the snapshot shown. All snapshots are taken at t ( N τ = . − − − Ma rms N u m b e r o f v e r ti ca l ce ll s λ h / ∆ y Fig. 12.
Expected typical vertical wavelength [Eq. (90)] in grid cells ofinternal gravity waves as a function of Ma rms in two series of simulationswith the Deviation well-balancing method. Black crosses correspond tosimulations at fixed resolution but increasing heating rate. Blue crossescorrespond to simulations at a fixed heating rate but di ff erent resolution.The two encircled data points result from the same simulation but Ma rms has been determined at a di ff erent time. cells as well as on a finer grid with 100 radial and 350 angularcells. Polar coordinates are the appropriate choice for the prob-lem’s geometry and should therefore lead to the least amountof numerical errors. We use periodic and solid-wall boundariesin the angular and radial directions, respectively. The flow hasa maximum Mach number of 1 at the inner domain boundarydropping to ≈ . ffi cient in this Mach number range. The timestep of explicit runs is set by the CFL uc criterion [Eq. (28)] with c CFL = .
4. For the lower-resolution setup we also perform sim- ulations with implicit time stepping and the CFL ug time step cri-terion [Eq. (26)] with c CFL = .
4. We find that the results ofthe implicit runs are identical to the explicit time stepping (seedots in upper panel of Fig. 14). Since the disk is isobaric, weuse (cid:37) - p reconstruction for this test case. Here we only comparesimulations without well-balancing to runs with Deviation well-balancing, since the other methods presented in this work are notcapable of stabilizing a target solution with a non-zero velocityfield (see Section 4.3).In order to asses the stability of the setup we add adensity perturbation with (cid:37) = x + . + y < . and follow its evolution up to 10 000 or-bital periods, where the orbital period is taken at the central pointof the perturbation. A perfect solution will maintain the initialradial density distribution of the perturbation at all times. How-ever, due to its radial extent there will be a phase shift betweenthe innermost and outermost regions of the perturbation. Further-more, numerical di ff usion will spread the perturbation predomi-nantly in the direction of its movement. Therefore the perturba-tion evolves into a homogeneous ring orbiting the central objectwith a density that corresponds to the initial angle-averaged den-sity (cid:37) .In Fig. 14 we show how the profile of (cid:37) − ×
70 cells the run without well-balancingonly shows small deviations from the target solution after 10 or-bital periods. However, we already see that mass starts to ac-cumulate in the center. The mass can not leave the domain andaccrete onto the central object due to our wall boundaries. Overtime, more and more mass flows to the inner boundary. After 100
Article number, page 18 of 24. V. F. Edelmann et al.: Well-balancing in astrophysical hydrodynamics ∆ y SZ / × ∆ y SZ / × ∆ y SZ / × ∆ y SZ / × ∆ y SZ /
16 288 × ∆ y SZ /
32 288 × − − − − − − − M a Fig. 13.
Mach number for a heating rate of ˙ e = and Deviation well-balancing at di ff erent resolutions. The upper left box in each panelindicates the change in the vertical resolution relative to the resolution used for the Mach number scaling. The lower right box gives the totalnumber of cells of the particular simulation. All snapshots are taken as soon as convection has fully developed. The dashed white line illustratesthe profile of the Brunt–Väisälä frequency as a function of height with arbitrary units on the x -axis. − − − − − − − − − % − No WBDeviation WB1 . . . . . . − − − − − − − % −
10 orbits100 orbits1000 orbits10000 orbits
Fig. 14.
Time evolution of the angle-averaged density profiles (cid:37) − ×
70 cells, where the dots on top of the lines representcorresponding results computed with implicit time stepping. The bot-tom panel shows the same quantity computed on a much finer grid of100 ×
350 cells. The black dashed lines in each panel show the respectiveinitial density profiles. orbits the center has approximately the same density as the ini-tial perturbation, and after 1000 orbits the profile has completelyshifted towards the inner boundary. This is also clearly visible in the two-dimensional density distribution shown in Fig. C.1. WithDeviation well-balancing the profile remains mostly stable up to1000 orbits. The perturbation is slowly di ff using symmetricallyaround the initial peak, maintaining the general shape of the ini-tial density profile. After 10 000 orbits we see that also the runwith Deviation well-balancing has become noticeably asymmet-ric towards the center. The regions where the density falls bellowthe initial density profile are most likely related to the fact thatwe do not use flux limiters for these runs. Undershoots at steepgradients are a common consequence of this omission.Increasing the resolution improves the stability of the runswithout well-balancing (see bottom panel in Fig. 14). At a reso-lution of 100 ×
350 cells the distribution is still almost identicalto the initial perturbation even after 100 orbits. Only after 1000orbits we start to notice the accumulation of mass at the innerboundary similar to the low-resolution case. After 10 000 orbits,however, the distribution has again completely shifted towardsthe center and the initial perturbation is not recognizable anylonger. The high-resolution run with Deviation well-balancing,on the other hand, is almost identical to the target solution evenafter 1000 orbits. The increased radial resolution has reduced theradial di ff usion observed in the low-resolution runs. The densityundershoots become more noticeable at the 10 000 orbit mark.There we again identify a tendency for a slow drift towards thecenter. However, thanks to the use of well-balancing the shapeof the initial perturbation is still retained approximately.This test shows that the flexibility of the Deviation well-balancing method also allows to maintain stable configurationsother than hydrostatic equilibrium. This is particularly impor-tant for the long-term evolution of such systems. Without well-balancing stability can only be achieved by increasing grid reso-lution, which leads to substantially higher computational cost. Article number, page 19 of 24 & A proofs: manuscript no. wellbalancing
6. Summary and conclusions
We have presented the Deviation, the α - β , and the Cargo–LeRoux well-balancing methods that aim to improve the abil-ity of finite-volume codes to maintain hydrostatic stratificationseven at moderate grid resolution. The performance of thesemethods were assessed in a set of test simulations of static anddynamical setups. Special emphasis was given to flows at lowMach numbers. They are particularly challenging to evolve be-cause they require special low-Mach hydrodynamic flux solvers,which in turn come with reduced dissipation and hence are proneto numerical instabilities. Also, it seems natural that slight devi-ations from hydrostatic equilibrium lead to low-Mach flows, asis the case, for example, in stellar convection. All simulationswere performed with the time-implicit SLH code that solves thefully compressible Euler equations using a modified version ofthe low-Mach AUSM + -up hydrodynamic flux solver. Our ex-perience shows that the inclusion of gravitational potential en-ergy in the total energy is essential to correctly representing slowflows in stratified atmospheres in the cell-centered discretizationof gravity implemented in SLH; see Mullen et al. (2020) for analternative approach. To the best of our knowledge, the presentstudy is the first to reach Mach numbers as low as 2 × − instratified convection using the fully compressible Euler equa-tions.The first test of the well-balancing schemes is to evolve a 1Dhydrostatic atmosphere in time at low resolution for an isother-mal, an isentropic, and a polytropic stratification. In all cases,the absence of a well-balancing scheme quickly led to spuri-ous velocities at significant amplitudes. The application of anyof the considered well-balancing methods removed this problemand managed to keep the flow below Mach numbers of 10 − forvery long times. Repeating the same test in 2D revealed that lowMach number flux functions, such as AUSM + -up, are subjectto an exponential growth of the Mach number and horizontaldensity fluctuations, which is not physically expected in a sta-bly stratified atmosphere. This e ff ect became less pronouncedthe closer the stratification was to the marginally stable, isen-tropic profile. In the isentropic case α - β and Deviation well-balancing in combination with the low-Mach flux remained sta-ble (Ma (cid:46) − ) for long times, while with Cargo–LeRoux well-balancing the setup developed a flow at Mach numbers of about10 − . These examples showed that it is important to test well-balanced schemes in more than one spatial dimension and formore than just a few sound-crossing times, as only then slowlygrowing instabilities become noticeable, especially in very sta-ble stratifications.While it is a necessary condition to maintain an initiallystatic setup, only dynamical setups are of actual interest inmultidimensional simulations. In a second test, we consideredthe rise of a hot bubble in a periodic background stratificationspanning two orders of magnitude in pressure at constant en-tropy. We tuned the bubble’s entropy excess to reach di ff erentrising speeds. With Cargo–LeRoux well-balancing and withoutany well-balancing, unphysical entropy fluctuations appeared inlarge parts of the atmosphere and a relatively fine grid (256 × ∼ × − seemed to be close to a limit of practical appli-cability of these two methods in such a strong stratification. The α - β and Deviation methods fared much better with no entropychanges far from the bubble and only a slight entropy undershootright above the bubble even on a coarse (64 ×
96) grid. Equallygood results were obtained with the initial amplitude decreased by a factor of 10 , leading to Ma ∼ × − . Numerical e ff ectsstarted to dominate only at Ma ∼ × − after the amplitudewas decreased by another factor of 10 .We proceeded with a setup involving a convection zone witha stable zone on top. The stratification was chosen to be repre-sentative of core convection in a 25 M (cid:12) main-sequence star. Ra-diation pressure was a substantial fraction of the total pressure,testing the methods’ capability to deal with a general EoS. Weused volume heating of adjustable amplitude to drive the con-vection. With Cargo–LeRoux well-balancing and without well-balancing, the rms Mach number of the convective flow ceasedto correlate with the heating rate at Ma rms ∼ × − and theflow became dominated by small-scale structures of numericalorigin. The lowest reachable value of Ma rms dropped by abouta factor of two when switching from (cid:37) - p to (cid:37) - T reconstructionin the absence of well-balancing. Only the α - β and Deviationmethods were able to reproduce the expected scaling of Ma rms with the heating rate [Eq. (87)] down to the slowest flows tested(Ma rms ∼ × − ). We observed spurious patterns in the sta-ble layer and demonstrated that they were caused by unresolvedinternal gravity waves, whose vertical wavelength becomes ex-tremely short with increasing period. The patterns disappearedwhen the typical vertical wavelengths of waves with periodsclose to the convective turnover frequency became resolved by 8to 10 cells (grid of 288 × orbits with a moderate resolution of 100 ×
350 cells.In summary, our results show that well-balancing can sub-stantially reduce the grid resolution needed to correctly followtiny perturbations in situations in which the stationary back-ground state involves the balance of two large opposing forces.We obtained comparable results with the α - β and Deviationmethods, both far surpassing the Cargo–LeRoux method in ac-curacy in the low-Mach-number regime. The α - β and Deviationmethods are also expected to be more accurate than the majorityof other well-balanced methods present in astrophysical litera-ture (e.g., Zingale et al. 2002; Perego et al. 2016; Käppeli et al.2011; Käppeli & Mishra 2016; Padioleau et al. 2019), becausethey exactly balance the stationary solution rather than an ap-proximation to it. Although an analytical prescription is often notavailable, the stationary solution can be computed to an arbitrarydegree of accuracy in many astrophysical applications. We pre-fer to use the Deviation method, because it is more general anddoes not impose any restrictions on reconstruction variables. Themethod can be applied both to nearly hydrostatic cases and tocases in which rotation becomes important. The latter can havestrong impact on stellar evolution (Maeder & Meynet 2000) aswell as on the propagation of IGW in stars (Rogers et al. 2013).The Cargo–LeRoux method still has its place, though. Horstet al. (2021) show for the case of convective helium-shell burn-ing that this method considerably reduces numerical di ff usion ofthe hydrostatic background stratification as compared to the un-balanced case if the gravitational energy is not included in thetotal energy. Therefore we consider the Cargo–LeRoux method Article number, page 20 of 24. V. F. Edelmann et al.: Well-balancing in astrophysical hydrodynamics as a valid method that leads to improvements already with littledevelopment e ff ort.This paper focuses on setups with simple grid geometries,but the SLH code can use general curvilinear grids. This allowsus to adapt grid geometry to the problem at hand. Standard spher-ical grids adapt to the geometry of slowly or non-rotating stars,but they su ff er from singularities at the center and at the poles.A star can be inscribed in a simple Cartesian grid (e.g. Wood-ward et al. 2015), but this requires to also impose a sphericalboundary condition on the grid. When done on the cell level, thesphere is rough on small scales and can generate spurious vor-ticity. We have implemented the cubed-sphere grid proposed byCalhoun et al. (2008) in our SLH code. The grid is logically rect-angular but geometrically spherical with a smooth outer bound-ary, see Fig. 15 for a 2D version. Discretization errors along thegrid’s strongly deformed diagonals make it almost impossibleto follow nearly hydrostatic flows without any well-balancingmethod. Figure 15 shows an SLH simulation with a convectiveshell embedded between two stably stratified zones. With De-viation well-balancing, we obtain the expected convective shellwith convection-generated IGW propagating in the stable zones.If we turn the well-balancing o ff , a convection-like flow of nu-merical origin appears and mixes the convection zone with thewhole inner stable zone. This result is promising, but limits ofapplicability of the Deviation method on the cubed-sphere gridare still to be investigated. Acknowledgements.
PVFE was supported by the US Department of Energythrough the Los Alamos National Laboratory (LANL). LANL is operated byTriad National Security, LLC, for the National Nuclear Security Administrationof U.S. Department of Energy (Contract No. 89233218CNA000001). LH, JPB,JH, RA, and FKR acknowledge support by the Klaus Tschira Foundation. Thework of CK and FKR is supported by the German Research Foundation (DFG)through grants KL 566 / / References
Almgren, A. S., Bell, J. B., Rendleman, C. A., & Zingale, M. 2006, ApJ, 637,922Andrassy, R., Herwig, F., Woodward, P., & Ritter, C. 2020, MNRAS, 491, 972Audusse, E., Bouchut, F., Bristeau, M.-O., Klein, R., & Perthame, B. 2004,SIAM Journal on Scientific Computing, 25, 2050Barsukow, W. & Berberich, J. P. 2020, Submitted to Journal on Scientific Com-putingBarsukow, W., Edelmann, P. V. F., Klingenberg, C., Miczek, F., & Röpke, F. K.2017, Journal of Scientific Computing, 72, 623Barsukow, W., Edelmann, P. V. F., Klingenberg, C., & Röpke, F. K. 2017, inESAIM: Proceedings and Surveys, Vol. 58, Workshop on low velocity flows,Paris, 5–6 November 2015, ed. S. Dellacherie & et al., 27–39Berberich, J. P. 2021, doctoralthesis, Universität WürzburgBerberich, J. P., Chandrashekar, P., & Klingenberg, C. 2018, in Theory, Numericsand Applications of Hyperbolic Problems I, Springer Proceedings in Mathe-matics & Statistics 236, ed. C. Klingenberg & M. WestdickenbergBerberich, J. P., Chandrashekar, P., & Klingenberg, C. 2021, Computers & Flu-ids, 104858Berberich, J. P., Chandrashekar, P., Klingenberg, C., & Röpke, F. K. 2019, Com-munications in Computational PhysicsBerberich, J. P., Käppeli, R., Chandrashekar, P., & Klingenberg, C. 2020, Ac-cepted in Communications in Computational Physics, arXiv:2005.01811Berberich, J. P. & Klingenberg, C. 2020, Accepted for publication in: SEMASIMAI Series: Numerical methods for hyperbolic problems Numhyp 2019Bermudez, A. & Vázquez, M. E. 1994, Computers & Fluids, 23, 1049Bolaños Rosales, A. 2016, Dissertation, Julius-Maximilians-UniversitätWürzburgBrowning, M. K., Brun, A. S., & Toomre, J. 2004, ApJ, 601, 512Brufau, P., Vázquez-Cendón, M., & García-Navarro, P. 2002, International Jour-nal for Numerical Methods in Fluids, 39, 247Calhoun, D. A., Helzel, C., & Leveque, R. J. 2008, SIAM Review, 50, 723Cargo, P. & Le Roux, A. 1994, Comptes rendus de l’Académie des sciences.Série 1, Mathématique, 318, 73 Castro, M. J. & Semplice, M. 2018, International Journal for Numerical Methodsin FluidsChandrashekar, P. & Klingenberg, C. 2015, SIAM Journal on Scientific Comput-ing, 37, B382Courant, R., Friedrichs, K. O., & Lewy, H. 1928, Math. Ann., 100, 32Cristini, A., Hirschi, R., Meakin, C., et al. 2019, MNRAS, 484, 4645Cristini, A., Meakin, C., Hirschi, R., et al. 2017, MNRAS, 471, 279Desveaux, V., Zenk, M., Berthon, C., & Klingenberg, C. 2016a, InternationalJournal for Numerical Methods in Fluids, 81, 104, fld.4177Desveaux, V., Zenk, M., Berthon, C., & Klingenberg, C. 2016b, Mathematics ofComputation, 85, 1571Edelmann, P. V. F. 2014, Dissertation, Technische Universität MünchenEdelmann, P. V. F., Ratnasingam, R. P., Pedersen, M. G., et al. 2019, ApJ, 876, 4Edelmann, P. V. F. & Röpke, F. K. 2016, in JUQUEEN Extreme Scaling Work-shop 2016, ed. D. Brömmel, W. Frings, & B. J. N. Wylie, JSC Internal ReportNo. FZJ-JSC-IB-2016-01, 63–67Edelmann, P. V. F., Röpke, F. K., Hirschi, R., Georgy, C., & Jones, S. 2017,A&A, 604, A25Edwards, J. R. & Liou, M.-S. 1998, AIAA journal, 36, 1610Gaburro, E., Castro, M. J., & Dumbser, M. 2018, MNRAS, 477, 2251Go ff rey, T., Pratt, J., Viallet, M., et al. 2017, A&A, 600, A7Grosheintz-Laval, L. & Käppeli, R. 2019, Journal of Computational Physics,378, 324Higl, J., Müller, E., & Weiss, A. 2021, A&A, 646, A133Horst, L., Edelmann, P. V. F., Andrássy, R., et al. 2020, A&A, 641, A18Horst, L., Hirschi, R., Edelmann, P. V. F., Andrássy, R., & Röpke, F. K. 2021,submitted to A&AHosea, M. & Shampine, L. 1996, Applied Numerical Mathematics, 20, 21 ,method of Lines for Time-Dependent ProblemsJones, S., Andrassy, R., Sandalski, S., et al. 2017, MNRAS, 465, 2991Käppeli, R. & Mishra, S. 2014, Journal of Computational Physics, 259, 199Käppeli, R. & Mishra, S. 2016, Astronomy & Astrophysics, 587, A94Käppeli, R., Whitehouse, S. C., Scheidegger, S., Pen, U. L., & Liebendörfer, M.2011, ApJS, 195, 20Kifonidis, K. & Müller, E. 2012, A&A, 544, A47Kippenhahn, R., Weigert, A., & Weiss, A. 2012, Stellar Structure and Evolution(Berlin Heidelberg: Springer-Verlag)Le Roux, A. Y. 1999, ESAIM: Proc., 6, 75LeVeque, R. J. 1998, Journal of computational physics, 146, 346Li, X.-s. & Gu, C.-w. 2008, Journal of Computational Physics, 227, 5144Liou, M.-S. 1996, Journal of Computational Physics, 129, 364Liou, M.-S. 2006, Journal of Computational Physics, 214, 137Liou, M.-S. & Ste ff en, C. J. J. 1993, Journal of Computational Physics, 107, 23Maeder, A. 2009, Physics, Formation and Evolution of Rotating Stars, Astron-omy and Astrophysics Library (Springer Berlin Heidelberg)Maeder, A. & Meynet, G. 2000, ARA&A, 38, 143Meakin, C. A. & Arnett, D. 2006, ApJ, 637, L53Meakin, C. A. & Arnett, D. 2007, ApJ, 667, 448Michel, A. 2019, Dissertation, Ruprecht-Karls-Universität HeidelbergMiczek, F. 2013, Dissertation, Technische Universität MünchenMiczek, F., Röpke, F. K., & Edelmann, P. V. F. 2015, A&A, 576, A50Mullen, P. D., Hanawa, T., & Gammie, C. F. 2020, arXiv e-prints,arXiv:2012.01340Müller, B., Viallet, M., Heger, A., & Janka, H.-T. 2016, ApJ, 833, 124Oßwald, K., Siegmund, A., Birken, P., Hannemann, V., & Meister, A. 2015, In-ternational Journal for Numerical Methods in Fluids, fld.4175Padioleau, T., Tremblin, P., Audit, E., Kestener, P., & Kokh, S. 2019, The Astro-physical Journal, 875, 128Perego, A., Cabezón, R. M., & Käppeli, R. 2016, Astrophysical Journal, Supple-ment, 223, 22Popov, M. V., Walder, R., Folini, D., et al. 2019, A&A, 630, A129Pratt, J., Bara ff e, I., Go ff rey, T., et al. 2020, A&A, 638, A15Pratt, J., Bara ff e, I., Go ff rey, T., et al. 2016, A&A, 593, A121Rieper, F. 2011, Journal of Computational Physics, 230, 5263Roe, P. L. 1981, Journal of Computational Physics, 43, 357Rogers, T. M., Lin, D. N. C., McElwaine, J. N., & Lau, H. H. B. 2013, ApJ, 772,21Röpke, F. K., Berberich, J., Edelmann, P. F. V., et al. 2018, in NIC Series, Vol. 49,NIC Symposium 2018, NIC Symposium 2018, Jülich (Germany), 22 Feb2018 - 23 Feb 2018 (Jülich: Forschungszentrum Jülich GmbH, Zentralbib-liothek), 115 – 122Shu, C.-W. & Osher, S. 1988, Journal of Computational Physics, 77, 439Sutherland, B. 2010, Internal Gravity Waves (Cambridge University Press)Tassoul, J.-L. 2000, Stellar RotationTimmes, F. X. & Swesty, F. D. 2000, ApJS, 126, 501Toro, E. F. 2009, Riemann Solvers and Numerical Methods for Fluid Dynamics:A Practical Introduction (Berlin Heidelberg: Springer)Touma, R. & Klingenberg, C. 2015, Applied Numerical Mathematics, 97, 42Touma, R., Koley, U., & Klingenberg, C. 2016, SIAM Journal on Scientific Com-puting, 38, B773Turkel, E. 1987, Journal of computational physics, 72, 277Viallet, M., Meakin, C., Arnett, D., & Mocák, M. 2013, ApJ, 769, 1Woodward, P. R., Herwig, F., & Lin, P.-H. 2013, ArXiv e-printsWoodward, P. R., Herwig, F., & Lin, P.-H. 2015, ApJ, 798, 49Zingale, M., Dursi, L. J., ZuHone, J., et al. 2002, ApJS, 143, 539 Article number, page 21 of 24 & A proofs: manuscript no. wellbalancing − − x / cm − − y / c m cubed-sphere grid − − x / cmno WB − − x / cmDeviation WB Fig. 15. Left panel:
Geometry of the cubed-sphere grid. The number of cells has been reduced compared to grids used in the other two panels toease the identification of individual cells and their shapes.
Middle and right panel:
Shell convection in a setup comparable to that in Section 5.3but with a shallower stratification and ideal gas EoS. The middle panel shows the color coded Mach number in a run without well-balancing.Numerical discretization errors quickly lead to spurious flows in the central region and the maximum Mach number reaches 3 × − . In the runwith Deviation well-balancing shown in the right panel the convective shell is clearly maintained. The maximum Mach number is 2 × − . α - β W B , % - p r ec . min = -8.35e-05max = ( ∆ A / A ) t = = t = +
02 s min = -8.33e-07max = ( ∆ A / A ) t = = t = +
03 s min = -8.33e-09max = ( ∆ A / A ) t = = t = +
04 s min = -8.35e-11max = ( ∆ A / A ) t = = t = +
05 s min = -2.83e-13max = ( ∆ A / A ) t = = t = +
06 s D e v i a ti on W B , % - p r ec . min = -1.10e-04max = = -1.10e-06max = = -1.10e-08max = = -1.10e-10max = = -7.29e-13max = Fig. A.1.
Same as Fig. 6, but showing the dependence of the solution on the initial amplitude ( ∆ A / A ) t = (left to right) on the 128 ×
192 grid forthe two best well-balancing methods (rows). The final time of the simulations is scaled with ( ∆ A / A ) / t = allowing the bubbles to reach the sameevolutionary stage, see Section 5.2. Appendix A: Hot bubbleAppendix B: Simple convective box setupAppendix C: Keplerian disk
Article number, page 22 of 24. V. F. Edelmann et al.: Well-balancing in astrophysical hydrodynamics α - β W B , % - p r ec . max = ( ∆ A / A ) t = = t = +
02 s max = ( ∆ A / A ) t = = t = +
03 s max = ( ∆ A / A ) t = = t = +
04 s max = ( ∆ A / A ) t = = t = +
05 s max = ( ∆ A / A ) t = = t = +
06 s D e v i a ti on W B , % - p r ec . max = = = = = Fig. A.2.
Mach-number distributions in the solutions shown in Fig. A.1. The color scheme ranges from 0 (black) to the maximum value (white)reached in the simulation box, which is also indicated at the top of the panels. . . . . . . y / cm × − . − . − . . ∆ ∇ ∆ y / c m × ∆ y SZ / ∆ y SZ / Fig. B.1.
Varying vertical grid spacing as a function of the vertical coor-dinate y for the simulations shown in the center and right panels of thelower row in Fig. 13. The superadiabaticity is shown as a blue line, anegative value indicates a convectively stable stratification. The spacingchanges smoothly to a finer resolution slightly before the transition tothe stable zone starts. Article number, page 23 of 24 & A proofs: manuscript no. wellbalancing no W B , e xp l ( % − = . % − = .
10 orbits ( % − = . % − = − .
100 orbits ( % − = . % − = − . ( % − = . % − = − . ( % − = . % − = − . D e v i a ti on W B , e xp l ( % − = . % − = . % − = . % − = − . % − = . % − = − . % − = . % − = − . % − = . % − = − . no W B , ( % − = . % − = . % − = . % − = − . % − = . % − = − . % − = . % − = − . % − = . % − = − . D e v i a ti on W B , ( % − = . % − = . % − = . % − = − . % − = . % − = − . % − = . % − = − . % − = . % − = − . no W B , i m p l ( % − = . % − = . % − = . % − = − . % − = . % − = − . % − = . % − = − . % − = . % − = − . D e v i a ti on W B , i m p l ( % − = . % − = . % − = . % − = − . % − = . % − = − . % − = . % − = − . % − = . % − = − . Fig. C.1.
Snapshots of the density distribution during the evolution of the Keplerian disk. Shown are runs with Deviation well-balancing andwithout well-balancing using explicit time stepping and resolutions of 20 ×
70 cells as well as 100 ×
350 grid cells. The bottom two rows show thesame setup evolved with implicit time stepping at a resolution of 20 ×
70 cells. To emphasize the deviation from the initial background density of (cid:37) = (cid:37) −−