---

A&A 506, 1065-1070 (2009)

Two-dimensional adaptive mesh refinement simulations of colliding flows

M. Niklaus^1,3 - W. Schmidt^2,3 - J. C. Niemeyer^2,3

1 - Deutsches Fernerkundungsdatenzentrum, Deutsches Zentrum für Luft- und Raumfahrt, Oberpfaffenhofen, Germany
2 - Institut für Astrophysik, Universität Göttingen, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany
3 - Lehrstuhl für Astronomie, Institut für Theoretische Physik und Astrophysik, Universität Würzburg, Am Hubland, 97074 Würzburg, Germany

Received 13 May 2009 / Accepted 5 August 2009

Abstract
Context. Colliding flows are a commonly used scenario for the formation of molecular clouds in numerical simulations. Turbulence is produced by cooling, because of the thermal instability of the warm neutral medium.
Aims. We carried out a two-dimensional numerical study of colliding flows to test whether statistical properties inferred from adaptive mesh refinement (AMR) simulations are robust with respect to the applied refinement criteria.
Methods. We compare probability density functions of various quantities, as well as the clump statistics and fractal dimension of the density fields in AMR simulations to a static-grid simulation. The static grid with 2048² cells matches the resolution of the most refined subgrids in the AMR simulations.
Results. The density statistics are reproduced fairly well by AMR. Refinement criteria based on the cooling time or the turbulence intensity appear to be superior to the standard technique of refinement by overdensity. Nevertheless, substantial differences in the flow structure become apparent.
Conclusions. In general, it is difficult to separate numerical effects from genuine physical processes in AMR simulations.

Key words: hydrodynamics - turbulence - instabilities - ISM: kinematics and dynamics - methods: numerical - ISM: clouds

1 Introduction

Computational fluid dynamics in astrophysics rely on numerical methods that are capable of covering a huge range of scales. Apart from smoothed particle hydrodynamics (Monaghan 1992), adaptive mesh refinement (AMR) has been applied to variety of problems. This method was developed by Berger & Oliger (1984) and Berger & Colella (1989). Among the widely used, publicly available AMR codes for astrophysical fluid dynamics are FLASH (Fryxell et al. 2000), Enzo (O'Shea et al. 2004) and Ramses (Teyssier 2002). Although there are comparative studies of AMR vs. SPH (for example, O'Shea et al. 2005; Agertz et al. 2007; Commerçon et al. 2008), the degree of reliance of AMR in comparison to non-adaptive methods has met only little attention so far.

Especially for turbulent flows, it is a non-trivial question whether the solutions obtained from AMR simulations agree with the correct solutions of the fluid dynamical equations at a given resolution level. For this reason, we systematically compare AMR and static-grid simulations for a particular test problem in this article. We chose a scenario that has been investigated in the context of molecular cloud formation, namely, the frontal collision of opposing flows of warm atomic hydrogen at supersonic speed (Hennebelle et al. 2008; Walder & Folini 2000; Hennebelle & Audit 2007a; Heitsch et al. 2006; Vázquez-Semadeni et al. 2007). Because of the cooling instability at densities and temperatures of a few thousand Kelvin, the gas becomes highly turbulent at the collision interface. Since the instabilities develop on length scales much smaller than the integral scale, this problem is computationally extremely demanding. The two-dimensional resolution study of Hennebelle & Audit (2007a) showed that the properties of the turbulent multi-phase medium evolving in these simulations is highly resolution-dependent, and numerical convergence is seen only at resolutions well above 1000². In three-dimensional simulations, such high resolutions are infeasible if static grids are used. Consequently, Hennebelle et al. (2008) and Banerjee et al. (2008) applied refinement by fixed density thresholds and refinement by Jeans mass, respectively, in their three-dimensional high-resolution AMR simulations.

In this article, we consider two-dimensional colliding flows without self-gravity and magnetic fields for a systematic comparison of AMR simulations to a reference simulation on a static grid. We analyze both statistical properties and the morphology of the gas fragmentation due to the cooling instability. This work is organized as follows: in Sect. 2 the numerical methods are described and the setup of the simulations will be presented in detail. In Sect. 3, we compare the results from the different simulations. Section 4 concludes this paper with a summary of the main results and general remarks on AMR.

2 Numerical methods and simulation setup

The simulations presented in this article are accomplished using the open source code Enzo (O'Shea et al. 2004; Bryan & Norman 1997). The compressible Euler equations are solved by means of the staggered grid, finite difference method Zeus (Stone & Norman 1992a,b; Stone et al. 1992). We included the cooling function defined by Audit & Hennebelle (2005) in these equations:

The primitive variables are the mass density , the velocity  and the specific total energy e of the fluid. The total energy per unit mass is given by

where is the adiabatic exponent and the pressure P is related to the mass density  and the temperature T via the ideal gas law:

The constants , and denote the Boltzmann constant, the mean molecular weight and the mass of the hydrogen atom, respectively. The gas is assumed to be a perfect gas with and .

The cooling function of Audit & Hennebelle (2005) includes the cooling by fine-structure lines of CII and OI, further the cooling by H (Ly-line) and the electron recombination onto positively charged grains. The heating is due to the photo-electric effect on small grains and polycyclic aromatic hydrocarbons (PAH) caused by the far-ultraviolet galactic background radiation. For more information about this cooling function see Bakes & Tielens (1994); Wolfire et al. (1995,2003); Spitzer (1978) and Habing (1968). The pressure-equillibrium curve resulting from the cooling function is plotted as black curve in Fig. 1. For the numerical solution of the fluid dynamical equations, we used the radiative cooling routine implemented into Enzo. For each hydrodynamical time step, the state variables are iterated over several subcycles, and the resulting total energy increment for the whole time step is added.

For our numerical study, the two-dimensional setup of Audit & Hennebelle (2005) and Hennebelle & Audit (2007b) was adopted with small modifications. The initial gas content corresponds to the warm neutral material (WNM) in the ISM, i.e., the temperature is T=7100 K, the pressure is P=7  10^-13 erg cm^-3 and the number density of neutral hydrogen is n=0.71 cm^-3. From the left and the right boundaries, warm gas with identical thermodynamic properties flows into the computational domain, where the cosine-shaped inflow velocity profile is modulated with small perturbations, realised by randomly shifted phases. These phase shifts are kept constant for the different simulations, so that initial conditions are exactly the same for all runs to ensure comparability. Following Hennebelle et al. (2008), the top and bottom boundary conditions are periodic. The physical dimensions of the computational domain are 20  20 pc. The two inflows of gas collide in the middle of the domain. The supersonic collision causes a steep rise in the gas density that triggers the thermal instability, and gas undergoes a transition into the phase of the cold neutral material (CNM) in the ISM. In this phase, the gas has temperatures in the range 30-100 K and number densities within 20-50 cm^-3 (Ferrière 2001). The thermal instability produces highly turbulent structures (see Fig. 2) with Mach numbers up to 20. The challenge for AMR is to track these turbulent structures as accurately as possible.

A reference simulation was run with a static grid of 2048² cells. Then the same setup was evolved in AMR simulations with a root-grid resolution of 128² cells and 4 levels of refinement. The resolution between adjacent refinement levels increases by a factor of 2. Hence, the effective resolution at the highest level of refinement is 2048². In these simulations, we employed three different types of refinement criteria:

1.

refinement by overdensity (OD);

2.

refinement by cooling time (CT);

3.

refinement by rate of compression and enstrophy (RCEN).

The first two criteria are widely used in astrophysical AMR simulations. For refinement by overdensity, the mass density must exceed the initial density on the root grid by a certain factor. This overdensity, in turn, defines the initial density for refinement at the first level of refinement and so on. We chose three different values for the overdensity factor, namley, twice the initial density (default OD), as well as three times (OD-3) and fourtimes (OD-4) the initial density. For criterion CT, on the other hand, refinement is triggered for a grid cell if the cooling time becomes less than the sound crossing time over the cell width. Refinement by the rate of compression and the enstrophy uses yet another technique. It was introduced by Schmidt et al. (2009) for the simulation of supersonic isothermal turbulence. The control variables for refinement are the enstrophy and the rate of compression. The enstrophy is given by one-half of the square of the vorticity, while the rate of compression is defined by the substantial time-derivative of the negative divergence of the velocity. The expression used by Schmidt et al. (2009) to evaluate the rate of compression (see Eq. (12) in this paper) is easily generalized to the non-isothermal case, where the speed of sound is not a constant. To trigger refinement by RCEN, dynamic thresholds are calculated from statistical moments of the control variables: a grid cell is flagged for refinement if the local fluctuation of a control variable becomes greater than the maximum of the average and the standard deviation of the variable. On the root grid, averages and variances are computed globally, whereas averaging is constrained to individual grid patches at higher levels of refinement.

For comparison of the simulation results, we calculated probability density functions (pdf) of several quantities. To analyze the gas fragmentation in each simulation, we adapted the clumpfind algorithm implemented by Padoan et al. (2007) to non-isothermal problems. The algorithm identifies the smallest dense regions that fulfill the Jeans criterion for gravitationally unstable gas. Since the clump samples found on the two-dimensional grids used in our simulations are insufficient for the calculation of clump mass spectra, only the total number and the mean size of the clumps are used for quantitative comparisons. In addition, we computed the fractal dimension of gas at densities higher than (corresponding to the minimum density of gas in the cold phase) by means of the box-counting method (Federrath et al. 2009).

3 Results

Due to the gradual accumulation of gas in the simulation domain, no strict statistical equilibrium is approached. For this reason, we evolved the flow until noticeable small-scale structure has developed and the separation of the gas into two phases has emerged. As shown in Fig. 1, two distinct phases are found at time . At this time, the central flow region is in a turbulent state (a contour plot of the mass density of the gas is shown in Fig. 2). Thus, we carry out our analysis for . While the main fraction of the gas is situated in the warm phase with temperatures between 5000 and 10 000 K and low densities  cm^-3, the cold gas with temperatures between 30 K and 100 K and densities in the range 30-350 cm^-3 can be found close to the equilibrium curve.

The pdfs of the mass density and the temperature obtained from different AMR simulations are plotted in Fig. 3. In principle, all refinement criteria reproduce the distributions found in the static-grid simulation quite well, although there is a trend of slightly more cold gas at the cost of warm gas. The discrepancy is more pronounced for refinement by over-density (OD) compared to the other criteria, and it becomes worse for the higher density thresholds (OD-3 and OD-4; not shown in the figure). Nevertheless, it appears that the thermodynamic properties of the gas are quite robust in AMR simulations.

The gravitationally unstable clumps of gas identified by the clumpfind algorithm in the static-grid simulation at time are depicted in Fig. 4a. The corresponding results for the AMR runs are shown in Figs. 4a-d. Table 1 lists the total number and the mean size of the clumps for each simulation. Also listed are the fractal dimensions of the gas regions with number density cm^-3, which are plotted in Figs. 4e-h.

Table 1:   Number N and average size of the clumps in the CNM; fractal dimension D of gas with number density greater than 20 cm^-3.

For refinement by OD, the fragmentation of the CNM is severely underestimated. The number of clumps is roughly half of the number in the static-grid simulation, and the clumps are typically larger. The lower degree of cold gas fragmentation results in a smaller fractal dimension (also see Fig. 4f). If the criteria OD-3 and OD-4 are applied, the number of clumps decreases further, while their average size increases. In the case of criterion OD-4, a slightly higher fractal dimension is obtained, because the cold phase tends to fill broad, area-filling regions. The cooling time criterion CT yields an amount of dense clumps with an average size that compares well to the reference simulation (see Fig. 4c), although the degree of fragmentation appears to be overestimated slightly. However, we found that this overestimation decreases with the further evolution of the colliding flows and, thus, appears to be transient. Refinement by RCEN also reproduces the number of clumps and the fractal dimension of dense gas very well. However, there are some anomalously big clumps, which contribute to an average clump size that is systematically too large. In the plot showing gas at density  cm^-3 (see Fig. 4h), on the other hand, such anomalous structures are not visible. Although refinement by RCEN does not overproduce gas in the cold phase (as one can see from the excellent agreement of the density and temperature pdfs in Figs. 3c and f), there appears to be a bias toward bigger clumps with this refinement method.

In contrast to the phase separation and gas fragmentation, striking deviations of the turbulent flow properties in the AMR vs. static grid simulations become apparent. Generally, a lot of turbulent small-scale structure is missing in the AMR simulations. Even for the criterion RCEN, which is based on control variables related to turbulence, this is apparent from the contour plots of the squared vorticity modulus shown in Fig. 5. Basically, the perturbations of the velocity field imposed at the inflow boundaries are quickly smoothed out in AMR simulations, so that turbulence is only produced by secondary (e.g., Kelvin-Helmholtz) instabilities at the collision interface in the central region of the computational domain. The reason is that all AMR cirteria, including RCEN, select relatively large fluctuations, whereas smaller perturbations are suppressed. On a static grid, on the other hand, the perturbations are transported from the boundaries to the centre and actively contribute to the production of turbulence. Consequently, small eddies are present in almost the whole domain in this case. Accordingly, the probability distribution of vorticity is markedly different (see Fig. 6). In contrast, Schmidt et al. (2009) found very close agreement of the vorticity pdfs in a static-grid and an AMR simulation with criterion RCEN for turbulence in a periodic box with large-scale forcing. Our results thus indicate that the merits of different refinement schemes are non-universal but rather depend on the properties of individual flow structures.

4 Conclusions

We performed two-dimensional simulations of colliding flows of warm atomic hydrogen with a radiative cooling function as source term in the energy equation. The goal of our study was the systematic comparison of AMR simulations, where different criteria for refinement were applied, to a reference simulation on a static grid.

While the probability distributions of mass density and temperature are well reproduced in AMR simulations, regardless of the refinement technique, differences become apparent in the fragmentation properties of the cold gas phase. As indicators, we used the total number of clumps and their average size. The clumps were identified by a clumpfind algorithm. In addition, we calculated the fractal dimension of dense gas, assuming a number density threshold of 20 cm^-3. Remarkably, the largest deviations from the clump statistics and fractal dimension extracted from the static-grid simulation, were encountered for refinement by overdensity, which is a commonly used refinement criterion in astrophysical AMR simulations. The deviations increase with the chosen density threshold. In this regard, it is important to note that Hennebelle et al. (2008) applied a density-based refinement criterion, where the thresholds were chosen even higher than those considered in our study. Good agreement, on the other hand, was obtained if the cooling time or enstrophy in combination with the rate of compression (the negative rate of change of the velocity divergence) were applied.

Table 2:   CPU time for the simulation runs.

Substantial problems with AMR became apparent with regard to turbulent flow properties. Basically, none of our AMR runs were able to reproduce even remotely the small-scale structure of turbulence and the probability distributions of turbulent flow variables such as the vorticity modulus. This defficiency can be attributed to the selection effects introduced by adaptive techniques. The definition of thresholds for triggering refinement either selects strong local fluctuations (for example, large shear that gives rise to Kelvin-Helmholtz instabilities) or large-scale perturbations such as accumulations of mass that become Jeans-unstable in self-gravitating gas. In this respect, the test problem we investigated in this work is particularly tough, because turbulence stems from small-scale instabilities that are seeded by weak initial perturbations. The varying grid resolution in AMR simulations inevitably modulate the growth of these instabilities and, as a consequence, the production of turbulence is suppressed. This defficiency might be overcome by the application of a subgrid scale model, which transports turbulent energy contained in small eddies that are resolved on finer grids across coarser grid regions (see Maier et al. 2009).

The key point of using AMR is the computational cost for a given effective resolution. Indeed, Table 2 demonstrates that a substantial reduction of computation time is achieved with AMR, especially if refinement by overdensity is applied. So, AMR is essentially a trade-off, where fast computation has to be weighted carefully against the question whether the essential physics of the specific problem is captured. Regarding three-dimensional simulations with high spatial resolution the specific impacts of the numerics are hard to predict, since the direct comparison to a high resolution computation on a static grid is not feasible and thus not available so far. Nevertheless, the physical results could be affected as much as in our two-dimensional simulations. So, besides the enormeous reduction in computational time especially in 3D, caution is to be recommended when discussing the results gained from AMR simultions and the conclusions drawn from them.

Acknowledgements

We thank Patrick Hennebelle and Edouard Audit for providing the cooling function that was used for this numerical study. We also thank Paolo Padoan for sharing his clumpfind tool. Computations described in this work were performed using the Enzo code developed by the Laboratory for Computational Astrophysics at the University of California in San Diego (http://lca.ucsd.edu).

References

• Agertz, O., Moore, B., Stadel, J., et al. 2007, MNRAS, 380, 963 [NASA ADS] [CrossRef]
• Audit, E., & Hennebelle, P. 2005, A&A, 433, 1 [NASA ADS] [CrossRef] [EDP Sciences]
• Bakes, E. L. O., & Tielens, A. G. G. M. 1994, ApJ, 427, 822 [NASA ADS] [CrossRef]
• Banerjee, R., Vazquez-Semadeni, E., Hennebelle, P., & Klessen, R. 2008 [arXiv:0808.0986]
• Berger, M. J., & Colella, P. 1989, J. Comp. Phys., 82, 64 [NASA ADS] [CrossRef]
• Berger, M. J., & Oliger, J. 1984, J. Comp. Phys., 53, 484 [NASA ADS] [CrossRef]
• Bryan, G. L., & Norman, M. L. 1997, in Computational Astrophysics; 12th Kingston Meeting on Theoretical Astrophysics, ed. D. A. Clarke, & M. J. West, ASP Conf. Ser., 123, 363
• Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2008, A&A, 482, 371 [NASA ADS] [CrossRef] [EDP Sciences]
• Federrath, C., Klessen, R. S., & Schmidt, W. 2009, ApJ, 692, 364 [NASA ADS] [CrossRef]
• Ferrière, K. M. 2001, Rev. Mod. Phys., 73, 1031 [NASA ADS] [CrossRef]
• Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273 [NASA ADS] [CrossRef]
• Habing, H. J. 1968, Bull. Astron. Inst. Netherlands, 19, 421 [NASA ADS]
• Heitsch, F., Slyz, A. D., Devriendt, J. E. G., Hartmann, L. W., & Burkert, A. 2006, ApJ, 648, 1052 [NASA ADS] [CrossRef]
• Hennebelle, P., & Audit, E. 2007a, A&A, 465, 431 [NASA ADS] [CrossRef] [EDP Sciences]
• Hennebelle, P., & Audit, E. 2007b, A&A, 465, 431 [NASA ADS] [CrossRef] [EDP Sciences]
• Hennebelle, P., Banerjee, R., Vázquez-Semadeni, E., Klessen, R. S., & Audit, E. 2008, A&A, 486, L43 [NASA ADS] [CrossRef] [EDP Sciences]
• Maier, A., Iapichino, L., Schmidt, W., & Niemeyer, J. 2009, ApJ, submitted
• Monaghan, J. J. 1992, ARA&A, 30, 543 [NASA ADS] [CrossRef]
• O'Shea, B. W., Bryan, G., Bordner, J., et al. 2004, arXiv Astrophysics e-prints
• O'Shea, B. W., Nagamine, K., Springel, V., Hernquist, L., & Norman, M. L. 2005, ApJS, 160, 1 [NASA ADS] [CrossRef]
• Padoan, P., Nordlund, Å., Kritsuk, A. G., Norman, M. L., & Li, P. S. 2007, ApJ, 661, 972 [NASA ADS] [CrossRef]
• Schmidt, W., Federrath, C., Hupp, M., Kern, S., & Niemeyer, J. C. 2009, A&A, 494, 127 [NASA ADS] [CrossRef] [EDP Sciences]
• Spitzer, L. 1978, Physical Processes in the Interstellar Medium, 1st edn. (Wiley-Interscience)
• Stone, J. M., & Norman, M. L. 1992a, ApJS, 80, 753 [NASA ADS] [CrossRef]
• Stone, J. M., & Norman, M. L. 1992b, ApJS, 80, 791 [NASA ADS] [CrossRef]
• Stone, J. M., Mihalas, D., & Norman, M. L. 1992, ApJS, 80, 819 [NASA ADS] [CrossRef]
• Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences]
• Vázquez-Semadeni, E., Gómez, G. C., Jappsen, A. K., et al. 2007, ApJ, 657, 870 [NASA ADS] [CrossRef]
• Walder, R., & Folini, D. 2000, Ap&SS, 274, 343 [NASA ADS] [CrossRef]
• Wolfire, M. G., Hollenbach, D., McKee, C. F., Tielens, A. G. G. M., & Bakes, E. L. O. 1995, ApJ, 443, 152 [NASA ADS] [CrossRef]
• Wolfire, M. G., McKee, C. F., Hollenbach, D., & Tielens, A. G. G. M. 2003, ApJ, 587, 278 [NASA ADS] [CrossRef]

 

All Tables

Table 1:   Number N and average size of the clumps in the CNM; fractal dimension D of gas with number density greater than 20 cm^-3.

Table 2:   CPU time for the simulation runs.

All Figures

In the text

In the text