A Field-length based refinement criterion for adaptive mesh simulations of the interstellar medium
Key Words.:hydrodynamics – instabilities – ISM: general, structure
Context:Adequate modelling of the multiphase interstellar medium requires optically thin radiative cooling, comprising an inherent thermal instability. The size of the occurring condensation and evaporation interfaces is determined by the so-called Field-length, which gives the dimension at which the instability is significantly damped by thermal conduction.
Aims:Our aim is to study the relevance of conduction scale effects in the numerical modelling of a bistable medium and check the applicability of conventional and alternative adaptive mesh techniques.
Methods:The low physical value of the thermal conduction within the ISM defines a multiscale problem, hence promoting the use of adaptive meshes. We here introduce a new refinement strategy that applies the Field condition by Koyama & Inutsuka as a refinement criterion. The described method is very similar to the Jeans criterion for gravitational instability by Truelove and efficiently allows to trace the unstable gas situated at the thermal interfaces.
Results:We present test computations that demonstrate the greater accuracy of the newly proposed refinement criterion in comparison to refinement based on the local density gradient. Apart from its usefulness as a refinement trigger, we do not find evidence in favour of the Field criterion as a prerequisite for numerical stability.
The structure of the turbulent interstellar medium (ISM) of the Milky Way and other disk galaxies is determined by various complex physical processes. In particular, the formation of the neutral atomic phase of the ISM is believed to be regulated by condensation under the action of thermally bistable radiative cooling. The related thermal instability (TI) has been extensively studied as a modifying agent for various driving sources of turbulence in the ISM. Numerical studies include the magneto-rotational instability (MRI) (Piontek & Ostriker 2005, 2007), explosions of SNe (Korpi et al. 1999; de Avillez & Breitschwerdt 2005; Joung & Mac Low 2006; Gressel et al. 2008a, b), and converging flows (Audit & Hennebelle 2005; Vázquez-Semadeni et al. 2006; Heitsch et al. 2008). Adaptive mesh MHD simulations of the last scenario have been performed by Hennebelle et al. (2008), who applied a simple density threshold as a refinement criterion.
The development of TI alone has been explored under various conditions by (Sánchez-Salcedo et al. 2002). Recently, there has been a discussion whether TI in conjunction with thermal conduction can act as an independent source of turbulence (Brandenburg et al. 2007; Koyama & Inutsuka 2006).
In dynamic simulations of the ISM, the radiative cooling is usually treated in the limit of an optically thin plasma, i.e., radiative transport effects are neglected. The cooling rate is then prescribed by in units of . Here we adopt the cooling function from Sánchez-Salcedo et al. (2002) that is a piecewise power-law fit to results by Wolfire et al. (1995).
In this paper, we want to focus on an issue in the numerical treatment of thermally bistable gas that has firstly been brought up by Koyama & Inutsuka (2004): The comprehensive theoretical analysis by Field (1965) reveals that, in the case of vanishing thermal conduction, TI has a finite growth rate in the limit of high wave numbers. This implies that, in the discrete representation of the fluid, numerical fluctuations at grid scale will be prone to artificial amplification through the physical instability unless one does take suitable precautions. The most natural way to guarantee a physically meaningful, converged solution is to introduce a physically motivated dissipation length that then needs to be resolved on the numerical grid. For the case when thermal conduction cannot be neglected, Field (1965) derives a stability criterion for the condensation mode which has the general form with the coefficient of thermal conduction in units of . In the case of a power-law cooling function this reduces to . The wavelength where the criterion is marginally fulfilled marks the scale where thermal conduction efficiently damps out fluctuations generated by the TI. Koyama & Inutsuka (2004) have shown that
has to be resolved with at least three grid cells to attain a converged numerical solution. It occurs naturally to use this criterion as a condition for mesh refinement.
2 Numerical Methods
For our simulations, we make use of the newly developed version 3 of the NIRVANA code which is a general purpose MHD fluid tool designed to solve multiscale, self-gravitating magnetohydrodynamics problems. The NIRVANA code comprises (i) a fully conservative, divergence-free Godunov-type central scheme, Ziegler (2004), (ii) block structured adaptive mesh refinement (AMR), and (iii) a multigrid-like adaptive mesh Poisson solver with elliptic matching, Ziegler (2005).
The central aim of adaptive mesh simulations is to resolve regions of interest at enhanced accuracy. In AMR codes like Enzo (O’Shea et al. 2004), FLASH (Fryxell et al. 2000), NIRVANA, or RAMSES (Teyssier 2002), this can be done by a variety of different criteria: baryon or dark matter overdensity (cosmological structure formation), local Jeans-length (protostellar core collapse), or entropy gradients (accretion shocks). Recently, Iapichino et al. (2008) have introduced a vorticity-based criterion to properly resolve a turbulent wake. Alternatively, independent of the underlying physics, the amount of discrepancy between the solutions on different AMR levels can be used as an indication for refinement (Berger & Colella 1989).
2.1 Classical mesh refinement
As a standard criterion to trace developing structures, one typically monitors local gradients in the fluid variables . To be able to follow moving structures (and to anticipate emerging structures), it is desirable to furthermore consider the temporal change of these gradients. Due to the additional overhead in storing quantities between the timesteps, this is usually discarded. Instead, one can make use of the fact that strong gradients are bordered by sections of enhanced curvature – reflected in the second derivatives . Accordingly, NIRVANA’s usual mesh refinement is determined by combining the criterion of normalised gradients with one of the ratio of second to first derivatives of the conserved variables , i.e.,
where and are the first and second numerical derivatives in the Euclidian norm, and controls the mixture of the two indicators. A new block is created if the criterion exceeds a critical value in a two-cell wide buffer zone around the block to be inserted. Blocks are removed if the criterion falls below a factor times the critical value. By tuning , the user can specify a bias towards the first or second type of criterion; a value of corresponds to the standard estimator used in the FLASH code (Fryxell et al. 2000).
For our runs not based on the Field condition, we adopt refinement for the gas density only and use trigger values of and while we set = . The ratio of grid spacings of refinement level and the base level permits a level-dependent refinement via the choice of the exponent . Here we fix , which makes refinements on higher levels successively easier. Finally, we set the additional filter value to prevent the refinement of small ripples.
While this general form of the refinement criterion is very flexible, the large number of free parameters renders the optimal adjustment to a particular problem a tedious task requiring a certain amount of experience. Moreover, since AMR always implies a trade-off between accuracy and efficiency, the optimal choice is by no means clear-cut. This is illustrated in Fig. 1, where we show the variation of the accuracy and AMR speed-up111For a description of the test case, see Sec. 3.1 below. as a function of the two supplementary parameters and of the classical refinement estimator (2). One can, in principle, use these graphs to find the parameter set which yields the highest execution speed at a given level of fidelity. However, this procedure is limited to simple test cases, and its universality cannot easily be proofed.
2.2 Refinement based on the Field condition
In an alternative approach particularly tailored to trace condensation interfaces, we apply grid refinement wherever the local Field-length is not resolved by at least three grid cells. For fine tuning, this number could, furthermore, be varied as a linear function of the refinement level – although we do not make use of this dispensable feature here. The main advantage of the proposed new method lies in the fact that it is virtually parameter free. This implies that the intricate determination of suitable values for and (as described above) can be avoided.
The evaluation of , given by Eq. (1), is straightforward since the values of and have already been computed for the heat flux term. To minimise overheads, they are stored in auxiliary arrays. It is evident that the definition of the Field-length requires the inclusion of a heat conduction term in the energy equation. Contrary to this, the FLASH code supports refinement based on the local cooling time , which is independent of . If, in our case, is scaled with (see below), the two methods should perform similarly.
3.1 Gaussian density perturbation
As a test problem, we revisit the 1D case of a density perturbation as considered in the model DEN3 from Sánchez-Salcedo et al. (2002), hereafter SSVSG. In this setup, the gas with density initially rests at its equilibrium pressure of and is perturbed by a Gaussian overdensity of amplitude , FWHM of , and centred within a domain of . Unlike in our simulations, SSVSG do not apply thermal conduction. This means that we have to choose high enough resolution, such that the conduction (which mainly becomes important near the grid scale) does not overly affect the features seen in their solution. As a reference, we perform a uni-grid run at a linear resolution of 8192 cells, corresponding to a grid spacing of . The actual level of conduction is set such that is resolved properly in the reference run. This yields a value for which is about a factor of 500 larger than the physical Spitzer value in the limit of vanishing electron density. The thus even smaller structures expected in a physically meaningful scenario illustrate the need for applying multiscale techniques.
Fig. 2 shows a plot of the Field-length for the initial and final density and temperature profiles, where we use three different prescriptions for the conduction: Scaling with is typically used together with constant kinematic viscosity yielding a constant value for the Prandtl-number along with a constant numerical time step, which is computationally favourable. This is also reflected in the fact that the Field-length only shows a weak variation with density. For the case of constant , this variation is stronger and accounts for about one order of magnitude in the problem under consideration. Finally, for a Spitzer-type conductivity with , we obtain a variation of over two orders of magnitude in rendering this case particularly interesting for AMR. In the following, we restrict ourselves to the case of a constant coefficient .
In general, the results are only marginally affected when conduction is included and agree well with the solution of SSVSG. In Fig. 3, the resulting profiles for the run with a 512 cell base-grid plus 4 levels of refinement (based on the Field condition) are plotted. Before we proceed with an analysis of the AMR performance, we want to compare our reference run 8192+0 (which is barely distinguishable from the one shown in Fig. 3) with the profiles from SSVSG: The most noticeable discrepancy is in the overall level of the thermal pressure. Due to the assumed periodicity at the domain boundaries, the mean pressure becomes a highly fluctuating quantity (cf. central panel of Fig. 3). Apart from the different offsets, the pressure profiles reasonably match. Notably, there is a small pressure dip at the interface of the condensed structure in our runs that is not seen in SSVSG. A similar feature is observed in panel (d) of their Fig. 7 (model DEN75), although much spikier. This leaves the possibility that the feature is not resolved by the plot for model DEN3 in Fig. 2 of SSVSG. A comparison run with in fact shows that, in our setup, the dip is broadened by the thermal conduction and remains much narrower otherwise. Later investigations on the topic consistently reveal a similar region of lower pressure (see e.g. Vázquez-Semadeni et al. 2003), which can be attributed to the fact that the pressure decreases with density in the unstable regime. In this respect, the pressure dip reflects the S-shape of the equilibrium cooling curve.
In Fig. 4, we plot the relative errors of three AMR runs with respect to the fully resolved solution. The two panels correspond to the two intervals with the strongest deviations, i.e., an outgoing wave (upper panel) and the condensation interface of the cloud (lower panel). We compare the two runs based on the density gradient criterion (with and ) with the run based on the Field condition. At the time the number of AMR blocks for these runs are , , and , respectively. Although the number of refined cells is comparable, the Field condition yields a relative error that is lower by a factor of two at the condensation front. For () the outgoing wave is resolved at (), while the Field condition, expectedly, is insensitive to this feature and does not refine it. Still the relative error is lowest in this case, indicating that the wave might be overly steepened by the refinement based on the density gradients. This finding is rather surprising and deserves to be looked at in greater detail.
As is illustrated in Fig. 5, thermal fragmentation produces turbulent and extremely filamentary structures. Since turbulence can only be modelled properly in three dimensions, adaptive mesh refinement becomes highly beneficial, if not mandatory. Because classical refinement strategies reach their limits if turbulence is involved (see e.g. Iapichino et al. 2008, for a vorticity-based refinement criterion), one has to seek for new tracers of structural properties. This is particularly true in the case of a thermally unstable medium.
Up to now, the criterion introduced by Koyama & Inutsuka (2004) is widely disregarded by many members of the numerical astrophysics community. Notable exceptions are a TI study by Brandenburg et al. (2007) and the MRI simulations of Piontek & Ostriker (2004). The opposite standpoint, represented by a number of authors (see e.g. Gazol et al. 2005; Joung & Mac Low 2006), is to neglect thermal conduction altogether. The general argumentation is that numerical diffusion defines a “numerical Field-length” that is thought to sufficiently damp small-scale modes of the instability. de Avillez & Breitschwerdt (2004), on the other hand, argue that the microscopic heat conduction is suppressed perpendicular to the magnetic field lines (and thus also isotropically for sufficiently tangled fields) and that turbulent transport takes an important role. Although this is certainly true, it does not aid the discussion of whether explicit conduction should be included.
In our own simulations, we only observe artificial growth of modes near the resolution limit in situations where the Courant number is chosen inappropriately high or the gradient based refinement criterion is improperly adjusted. The absence of unphysical growth at high wavenumbers might thus, in fact, be attributed to the finite level of numerical dissipation, which is never really negligible in three-dimensional simulations.
In this respect, the Field criterion of Koyama & Inutsuka indeed seems of academic interest only. Apart from its actual necessity, we, however, demonstrate that the very condition can successfully be used as a refinement criterion for adaptive mesh simulations of the interstellar medium – an approach that is very similar to the refinement procedure based on the Jeans-length as introduced by Truelove et al. (1997) for self-gravitating clouds. In a simple 1D test case, the new refinement estimator is found to produce more accurate results at comparable numerical cost than conventional criteria.
The overhead of evaluating is low compared to the computation of the numerical fluxes. Within NIRVANA, the determination of the heat flux already requires the evaluation of the gas temperature and the (non-uniform) conductivity coefficient. These fields can therefore be stored temporarily such that the only additional expense comes from evaluating the cooling function.
The definition of , nevertheless, implies the inclusion of thermal conduction, which is not a standard feature in many available codes. An alternative approach (implemented in the FLASH code) is to use the cooling time instead, which is independent of . On the other hand, there is evidence (Piontek, in prep.) that heat conduction is in fact physically relevant for the formation of molecular clouds, providing further motivation for the use of the proposed scheme.
Compared to the classical mesh refinement based on gradients, which has to be fine-tuned via multiple parameters according to different situations, the Field criterion is virtually parameter-free and gives appropriate results irrespective of the particular setup. If, of course, other features like shock fronts are supposed to be adequately resolved, one still has to rely on a combination with conventional refinement criteria; this does, however, not impose any difficulties.
Acknowledgements.This work used the NIRVANA code v3.3 developed by Udo Ziegler at the Astrophysical Institute Potsdam. We thank Enrique Vazquez and Simon Glover for the discussion related to this work and acknowledge the helpful comments by the anonymous referee.
- Audit & Hennebelle (2005) Audit, E. & Hennebelle, P. 2005, A&A, 433, 1
- Berger & Colella (1989) Berger, M. J. & Colella, P. 1989, JCoPh, 82, 64
- Brandenburg et al. (2007) Brandenburg, A., Korpi, M. J., & Mee, A. J. 2007, ApJ, 654, 945
- de Avillez & Breitschwerdt (2004) de Avillez, M. A. & Breitschwerdt, D. 2004, A&A, 425, 899
- de Avillez & Breitschwerdt (2005) de Avillez, M. A. & Breitschwerdt, D. 2005, A&A, 436, 585
- Field (1965) Field, G. B. 1965, ApJ, 142, 531
- Fryxell et al. (2000) Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273
- Gazol et al. (2005) Gazol, A., Vázquez-Semadeni, E., & Kim, J. 2005, ApJ, 630, 911
- Gressel et al. (2008b) Gressel, O., Elstner, D., Ziegler, U., & Rüdiger, G. 2008b, A&A, 486, L35
- Gressel et al. (2008a) Gressel, O., Ziegler, U., Elstner, D., & Rüdiger, G. 2008a, AN, 329, 619
- Heitsch et al. (2008) Heitsch, F., Hartmann, L. W., & Burkert, A. 2008, ApJ, 683, 786
- Hennebelle et al. (2008) Hennebelle, P., Banerjee, R., Vázquez-Semadeni, E., Klessen, R. S., & Audit, E. 2008, A&A, 486, L43
- Iapichino et al. (2008) Iapichino, L., Adamek, J., Schmidt, W., & Niemeyer, J. C. 2008, MNRAS, 388, 1079
- Joung & Mac Low (2006) Joung, M. K. R. & Mac Low, M.-M. 2006, ApJ, 653, 1266
- Korpi et al. (1999) Korpi, M. J., Brandenburg, A., Shukurov, A., Tuominen, I., & Nordlund, Å. 1999, ApJ, 514, L99
- Koyama & Inutsuka (2004) Koyama, H. & Inutsuka, S.-i. 2004, ApJ, 602, L25
- Koyama & Inutsuka (2006) Koyama, H. & Inutsuka, S.-i. 2006, [astro-ph/0605528]
- O’Shea et al. (2004) O’Shea, B. W., Bryan, G., Bordner, J., et al. 2004, [astro-ph/0403044]
- Piontek & Ostriker (2004) Piontek, R. A. & Ostriker, E. C. 2004, ApJ, 601, 905
- Piontek & Ostriker (2005) Piontek, R. A. & Ostriker, E. C. 2005, ApJ, 629, 849
- Piontek & Ostriker (2007) Piontek, R. A. & Ostriker, E. C. 2007, ApJ, 663, 183
- Sánchez-Salcedo et al. (2002) Sánchez-Salcedo, F. J., Vázquez-Semadeni, E., & Gazol, A. 2002, ApJ, 577, 768
- Teyssier (2002) Teyssier, R. 2002, A&A, 385, 337
- Truelove et al. (1997) Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, L179+
- Vázquez-Semadeni et al. (2003) Vázquez-Semadeni, E., Gazol, A., & Passot, T. 2003, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 614, Turbulence and Magnetic Fields in Astrophysics, ed. E. Falgarone & T. Passot, 213–251
- Vázquez-Semadeni et al. (2006) Vázquez-Semadeni, E., Ryu, D., Passot, T., González, R. F., & Gazol, A. 2006, ApJ, 643, 245
- Wolfire et al. (1995) Wolfire, M. G., Hollenbach, D., McKee, C. F., Tielens, A. G. G. M., & Bakes, E. L. O. 1995, ApJ, 443, 152
- Ziegler (2004) Ziegler, U. 2004, JCoPh, 196, 393
- Ziegler (2005) Ziegler, U. 2005, CoPhC, 170, 153