EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 563 (March 2014) A&A, 563 (2014) A11 Full HTML
Free Access
Issue
A&A
Volume 563, March 2014
Article Number A11
Number of page(s) 11
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201322858
Published online 26 February 2014

© ESO, 2014

1. Introduction

The study of structure formation in the Universe involves multiscale highly nonlinear physics, such as hydrodynamics, radiative transfer, gravity, and magnetic fields. Numerical experiments are the best laboratory for studying these structures, but they remain challenging. Thanks to the formidable development of supercomputing facilities, these numerical experiments can integrate many different physical processes and use thousands of processors to achieve unprecedented numerical resolution. Nevertheless, efficient scaling often becomes problematic because of the variety of dynamically important physical processes involved. For instance, some physical processes, such as radiative transfer, involve dynamical timescales that are much shorter than in hydrodynamics. If hydrodynamics and radiative transfer are coupled in a unique nonrelativistic system of equations, the time step at which this system can be integrated is limited by the one derived from radiation transport at the speed of light. Implicit methods have thus been developed and coupled to hydrodynamical solvers to handle the short characteristic timescales of physical processes such as the diffusion. In general, hydrodynamical codes use an operator splitting approach with explicit solvers to integrate the Euler equations and implicit solvers to deal with diffusion-like problems. This coupling of explicit and implicit solvers is relatively straightforward and well-studied on uniform grids (e.g., Turner & Stone 2001), but becomes far more difficult on complex grids, such as those generated by adaptive mesh-refinement (AMR, see the RAMSEScode Teyssier 2002; Commerçon et al. 2011c).

To illustrate the main difficulties of designing an implicit method for AMR grids, let us first consider the simple heat equation in one dimension. The second-order parabolic partial differential equations can be generalized as a diffusion problem following ∂tU−K∂xx2U=0,\begin{equation} \partial_t U-K\partial^2_{xx}U=0, \label{heat} \end{equation}(1)where U(x,t) is function of position x and time t, and K is the diffusion coefficient. Numerically, equation (1) can be integrated with explicit or implicit methods to advance from time level n to time level n + 1. Explicit discretization of (1) leads to the Courant Friedrich Levy (CFL) stability condition Δtexp<Δx22K,\begin{equation} \Delta t_\mathrm{exp}<\frac{\Delta x^2}{2K}, \label{CFL} \end{equation}(2)with Δx the size of the discretized mesh. The explicit CFL condition for diffusion equation scales as Δx2 and is thus far more stringent than the classical CFL derived for the stability of the hyperbolic system formed by the Euler equations (Δthyp < Δx/cs, with cs the gas sound speed). This can lead to extensive computing time when both hyperbolic and parabolic equations are treated simultaneously. For that reason, Eq. (1) is often integrated using an implicit scheme, which is unconditionally stable. The implicit scheme requires solving a matrix system of equation Ax = b, where matrix A has to be inverted to get the vector solution x. While matrix inversion does not present strong conceptual difficulties for uniform grids using preconditioned iterative methods, the problem becomes challenging when the grid is complex like the one generated by RAMSES .

Another feature of AMR codes is the use of adaptive time-stepping (ATS) in their hydrodynamical solvers (Teyssier 2002; Almgren et al. 2010; Bryan et al. 2013) to speed calculations up. For implicit solvers, ATS is not as straightforward as for an explicit scheme, since an operator integrated with an implicit scheme can affect all the grids of the computational domain. Some authors have designed the ATS method for implicit schemes derived from diffusion equations (Howell & Greenough 2003; Zhang et al. 2011). In these methods, the diffusion equation is updated on a level-by-level basis, and the total radiative energy is conserved by storing flux at the level interfaces.

In a previous paper, Commerçon et al. (2011c, hereafter Paper I) proposed a method that integrates a diffusion-like equation in the RAMSES code for radiation hydrodynamics (RHD) using a two-temperature approach. The method in Paper I uses a unique time step for all the levels and does not take advantage of the ATS method developed for the hydrodynamical solver in RAMSES . The purpose of this paper is to present a new ATS method for implicit solvers on AMR grids in order to speed up the solver of Paper I. We seek to keep the method as simple as possible to allow quick implementation in other codes using ATS for their hydrodynamical solver.

The paper is organized as follows. In Sect. 2, we recall the RHD equations we use and briefly present the flux-limited diffusion solver we designed in Paper I. The new implicit solver for the RAMSES code is presented in Sect. 3. In Sect. 4, the method is then tested against well-known test cases for diffusion and RHD. As a final test, RHD dense core collapse calculations with very high resolution are performed in Sect. 5. Section 6 summarizes our work and the main results, and presents our perspectives.

2. Basic concepts

thumbnail Fig. 1

Two-dimensional sketch of possible AMR grid configurations at the interface between two levels. Left: configuration 1 where cell i is at level ℓ and cell i − 1 at level ℓ + 1. Right: configuration 2 where cell i is at level ℓ and cell i − 1 at level ℓ − 1. In both sketches, the vertical red line represents the surface over which energy is exchanged when level ℓ is updated. Similarly, the horizontal dashed red line represents the gradient over which the flux is computed using the energy marked by a red star (the value is interpolated using the value represented by the attached blue circle).

2.1. Radiation hydrodynamics in RAMSES

In Paper I, we presented an implementation of a RHD solver into RAMSES using the flux-limited diffusion approximation (FLD, e.g. Minerbo 1978; Levermore & Pomraning 1981). The RHD equations with all the radiative quantities estimated in the comoving frame then read {∂tρ+∇[ρu]=0∂tρu+∇[ρu⊗u+PI]=−λ∇Er∂tET+∇[u(ET+P)]=−Pr∇:u−λu∇Er+∇·(cλρκR∇Er)∂tEr+∇[uEr]=−Pr∇:u+∇·(cλρκR∇Er)+κPρc(aRT4−Er),\begin{equation} \left\{ \begin{array}{llll} \partial_t \rho + \nabla \left[\rho\vec{u} \right] & = & 0 \\ \partial_t \rho \vec{u} + \nabla \left[\rho \vec{u}\otimes \vec{u} + P \mathbb{I} \right]& =& - \lambda\nabla E_\mathrm{r} \\ \partial_t E_\mathrm{T} + \nabla \left[\vec{u}\left( E_\mathrm{T} + P_\mathrm{} \right)\right] &= & - \mathbb{P}_\mathrm{r}\nabla:\vec{u} - \lambda \vec{u} \nabla E_\mathrm{r} \\ & & + \nabla \cdot\left(\frac{c\lambda}{\rho \kappa_\mathrm{R}} \nabla E_\mathrm{r}\right) \\ \partial_t E_\mathrm{r} + \nabla \left[\vec{u}E_\mathrm{r}\right] &=& - \mathbb{P}_\mathrm{r}\nabla:\vec{u} + \nabla \cdot\left(\frac{c\lambda}{\rho \kappa_\mathrm{R}} \nabla E_\mathrm{r}\right) \\ & & + \kappa_\mathrm{P}\rho c(a_\mathrm{R}T^4 - E_\mathrm{r}), \end{array} \right. \end{equation}(3)where ρ is the material density,u is the velocity, P is the thermal pressure, κR is the Rosseland mean opacity, λ is the radiative flux limiter, ET is the total energy ET = ρϵ + 1/2ρu2 + Er (ϵ is the internal specific energy), κP is the Planck opacity, Er is the radiative energy, and Pr is the radiation pressure.

The method presented in Paper I is based on an operator splitting scheme, where the hydrodynamical part is integrated using the hyperbolic explicit solver of RAMSES and the radiative energy diffusion and coupling between matter and radiation terms are integrated using an implicit scheme. The method uses a conjugate gradient algorithm in which all the levels of the AMR grid are coupled so that calculations advance in time following the CFL conditions of the finest level of refinement. The main limitation of this method is that it uses a unique time step, and calculations can become very expensive in numerical experiments involving a large hierarchy of AMR levels. In the following, we present an implementation of ATS for the implicit method presented in Paper I.

2.2. Adaptive time-stepping for hydrodynamics

Most astrophysical problems deal with a large range of physical scales, as for instance in star formation, where scales for the size of the cloud (unit of a parsec) and of the protostar (unit of solar radius, R⊙ ~ 10-8 pc) have to be considered simultaneously. This produces a high level of hierarchy in AMR grids, each level ℓ having a size of Δxℓ. For the classical Euler system of equations (conservation of mass, momentum and total energy), a stability condition can be calculated for each level following the CFL condition, namely, Δt=CΔx|u|+cs,\begin{equation} \Delta t=C\frac{\Delta x}{|\vec{u}|+c_\mathrm{s}}, \end{equation}(4)where C < 1 is the CFL number and | u | is the fluid velocity norm.

To integrate the Euler equations on AMR grids, a unique time step can be used, meaning that all the levels evolve with the same time step, given by the CFL condition on the finest level ℓmax. This method is very powerful, but can be expensive when a high number of AMR levels is used. In the ATS scheme, each level evolves with its own time step which considerably reduces the CPU time. RAMSES uses ATS for the hydrodynamics solver (Teyssier 2002), following the rule that the time step of a level ℓ equals always the sum of two time steps on the finer level, i.e., Δtℓ=Δt1ℓ+1+Δt2ℓ+1\hbox{$\Delta t^\ell=\Delta t^{\ell+1}_1+\Delta t^{\ell+1}_2$}. Ideally, if all the levels of the grid hierarchy are effective and if the problem is isothermal and has a uniform velocity, we have Δtℓmax = 2ℓmaxΔt0, where Δt0 is the time step of the coarser level. The scheme is only first-order accurate at level interfaces, but the errors are localized only at level interfaces and generally the ratio between the interface surface and the volume of the computational domain is relatively small (otherwise a uniform grid would be used, Teyssier 2002). In addition, truncation errors can propagate only if waves move from coarse to fine grids, which is not the case in accretion flows (the accretion shock moves from fine to coarse). In the RAMSES implementation, the fine levels are updated first. When a cell i + 1 at level ℓ is updated, the flux Fi+1/2n+Δtℓ−1\hbox{$F^{n+\Delta t^{\ell-1}}_{i+1/2}$} that crosses the interface with a cell i at a coarser level ℓ − 1 during the two fine time steps Δt1ℓ+Δt2ℓ\hbox{$\Delta t_1^\ell+\Delta t_2^\ell$} is Fi+1/2n+Δtℓ−1=1Δt1ℓ+Δt2ℓ(Δt1ℓFi+1/2n+Δt1ℓ+Δt2ℓFi+1/2n+Δt1ℓ+Δt2ℓ),\begin{equation} F^{n+\Delta t^{\ell-1}}_{i+1/2}=\frac{1}{\Delta t_1^\ell+\Delta t_2^\ell}\left(\Delta t_1^\ell F^{n+\Delta t_1^\ell}_{i+1/2} + \Delta t_2^\ell F^{n+\Delta t_1^\ell+\Delta t_2^\ell}_{i+1/2}\right),\label{mean_flux} \end{equation}(5)such that the total conservative variables (mass, momentum and total energy) are conserved. The aim of this paper is to couple an implicit solver to the hydrodynamical ATS scheme used in RAMSES .

3. Numerical method

3.1. Definitions

Let us consider the following diffusion equation on the radiative energy Er∂Er∂t=∇.K∇Er,\begin{equation} \frac{\partial E_\mathrm{r}}{\partial t}= \nabla. K\nabla E_\mathrm{r}, \label{Er_diff} \end{equation}(6)The finite volume discretization in the x direction1 of Eq. (6) using an implicit scheme gives Er,in+1−Er,inΔtVi=Ki+1/2Er,i+1n+1−Er,in+1ΔxSi+1/2\begin{eqnarray} \frac{E_{\mathrm{r},i}^{n+1}-E_{\mathrm{r},i}^{n}}{\Delta t}V_i & = & K_{i+1/2}\frac{E_{\mathrm{r},i+1}^{n+1}-E_{\mathrm{r},i}^{n+1}}{\Delta x}S_{i+1/2} \nonumber\\ & & - K_{i-1/2}\frac{E_{\mathrm{r},i}^{n+1}-E_{\mathrm{r},i-1}^{n+1}}{\Delta x}S_{i-1/2},\label{Er_imp} \end{eqnarray}(7)where Vi is the volume of cell i, Si ± 1/2 the surface of exchange between cell i and cell i ± 1, Ki ± 1/2 is the mean diffusion coefficient computed at the cell interface (e.g., Ki ± 1/2 = (Ki ± 1 + Ki)/2). Equation (7) can be written in the form −Ci−1/2Er,i−1n+1+(1+Ci−1/2+Ci+1/2)Er,in+1−Ci+1/2Er,i+1n+1=Er,in,\begin{equation} -C_{i-1/2}E_{\mathrm{r},i-1}^{n+1} +(1 + C_{i-1/2} + C_{i+1/2})E_{\mathrm{r},i}^{n+1} -C_{i+1/2} E_{\mathrm{r},i+1}^{n+1} = E_{\mathrm{r},i}^{n}, \label{Er_imp1} \end{equation}(8)where Ci ± 1/2 = Ki ± 1/2Si ± 1/2Δt/(ΔxVi). Equation (8) forms a matrix system, Ax = b, where matrix A has to be inverted to get the new value of the radiative energy Ern+1\hbox{$E_{\mathrm{r}}^{n+1}$} (the solution vector x). The Ci ± 1/2 coefficient depends on the grid configuration, namely if the neighboring cells i −1 and i +1 are or not at the same level of refinement as cell i. In the first configuration, Fig. 1a, cell i is at level ℓ and cell i − 1 at level ℓ + 1. In the simplest case, the neighboring cells at level ℓ + 1 are interpolated on a coarser cell at level ℓ, so that the Cℓ → ℓ + 1 coefficient calculations is reduced to the one on a uniform mesh. We thus have Ci−1/2ℓ→ℓ+1=Ki−1/2(Δxℓ)ndim−1ΔtℓΔxℓ(Δxℓ)ndim,=Ki−1/2Δtℓ(Δxℓ)2,\begin{eqnarray} C^{\ell\rightarrow\ell+1}_{i-1/2} & = & \frac{K_{i-1/2}(\Delta x^{\ell})^{ndim-1}\Delta t^\ell}{\Delta x^{\ell}(\Delta x^{\ell})^{ndim}},\nonumber\\ & = & \frac{K_{i-1/2}\Delta t^\ell}{(\Delta x^{\ell})^2}, \end{eqnarray}(9)where ndim is the number of dimensions of the problem. In the opposite case, Fig. 1b, where cell i is at level ℓ and cell i − 1 at level ℓ − 1, we have Ci−1/2ℓ→ℓ−1=Ki−1/2(Δxℓ)ndim−1Δtℓ3Δxℓ2(Δxℓ)ndim,=23Ki−1/2Δtℓ(Δxℓ)2·\begin{eqnarray} C^{\ell\rightarrow\ell-1}_{i-1/2} & = & \frac{K_{i-1/2}(\Delta x^{\ell})^{ndim-1}\Delta t^\ell}{\frac{3\Delta x^{\ell}}{2}(\Delta x^{\ell})^{ndim}}, \nonumber\\ & = & \frac{2}{3}\frac{K_{i-1/2}\Delta t^\ell}{(\Delta x^{\ell})^2}\cdot \end{eqnarray}(10)We define ˜K\hbox{$\tilde{K}$} as the mean diffusion coefficient at cells interface, and Aℓ=˜KΔtℓ/(Δxℓ)2\hbox{$A^\ell=\tilde{K}\Delta t^\ell/(\Delta x^\ell)^2$}. For a given level ℓ, we thus have three types of coefficient C, namely Cℓ → ℓ+1, Cℓ → ℓ, and, Cℓ → ℓ−1, depending on the grid configuration at cells interface Cℓ→ℓ−1=23Aℓ,Cℓ→ℓ=Aℓ,Cℓ→ℓ+1=Aℓ.\begin{eqnarray} C^{\ell\rightarrow\ell-1} & = &\frac{2}{3}A^\ell,\\ C^{\ell\rightarrow\ell\hspace{8pt}} & = &A^\ell,\\ C^{\ell\rightarrow\ell+1} & = & A^\ell. \end{eqnarray}In the previous implementation of Paper I, the cells at level ℓ + 1 are not interpolated at level ℓ and the Cℓ → ℓ + 1 coefficient equals 23−ndimAℓ/3.

In the following, we solve the diffusion equation for each level ℓ independently from the other levels. In the case where a cell of level ℓ is at an interface with a coarser or a finer level, we need to specify a boundary condition at the level interfaces to solve the matrix system given by Eq. (8). We now study different types of boundary condition at level interfaces that can be used to compute the corresponding flux between cells i and i − 1, so that Ci−1/2 = Cℓ → ℓ+1 or Ci−1/2 = Cℓ → ℓ−1. If cell i is at level ℓ, cell i − 1 can be either at level ℓ + 1 or at level ℓ − 1. We assume that cells i and i + 1 are at the same level ℓ.

3.2. Different types of boundary conditions

3.2.1. Dirichlet boundary condition

The Dirichlet boundary is an imposed boundary condition, i.e., Er,i−1n+1=Er,i−1n\hbox{$E^{n+1}_{\mathrm{r},i-1} = E^{n}_{\mathrm{r},i-1}$}. Equation (8) then reads (1+Ci−1/2+Ci+1/2)Er,in+1−Ci+1/2Er,i+1n+1=Er,in+Ci−1/2Er,i−1n,\begin{eqnarray*} (1 + C_{i-1/2} + C_{i+1/2})E_{\mathrm{r},i}^{n+1} -C_{i+1/2} E_{\mathrm{r},i+1}^{n+1} = E_{\mathrm{r},i}^{n} + C_{i-1/2}E_{\mathrm{r},i-1}^{n},\nonumber \end{eqnarray*}where we moved the terms corresponding to Er, i−1 on the rand-hand-side (RHS) of the matrix system, i.e., in the b vector.

3.2.2. Neumann boundary condition

The Neumann boundary corresponds to an imposed flux Fi−1/2, i.e., Ci−1/2(Er,in+1−Er,i−1n+1)=Fi−1/2\hbox{$C_{i-1/2}(E^{n+1}_{\mathrm{r},i}-E^{n+1}_{\mathrm{r},i-1}) = F_{i-1/2}$}. Equation (8) reads (1+Ci+1/2)Er,in+1−Ci+1/2Er,i+1n+1=Er,in+Fi−1/2,\begin{eqnarray*} (1 + C_{i+1/2})E_{\mathrm{r},i}^{n+1} - C_{i+1/2} E_{\mathrm{r},i+1}^{n+1} = E_{\mathrm{r},i}^{n} + F_{i-1/2},\nonumber \end{eqnarray*}where the imposed flux is also computed in the RHS of the matrix system.

3.2.3. Robin boundary condition

The Robin boundary corresponds to a mix between the Dirichlet and the Neumann boundary conditions, i.e., the energy exchange at level interfaces corresponds to (1−α)Fi−1/2+α(Er,in+1−Er,i−1n)/Δx\hbox{$(1-\alpha)F_{i-1/2} + \alpha (E^{n+1}_{\mathrm{r},i}-E^{n}_{\mathrm{r},i-1})/\Delta x$}. Equation (8) reads (1+αCi−1/2+Ci+1/2)Er,in+1−Ci+1/2Er,i+1n+1=Er,in+αCi−1/2Er,i−1n+(1−α)Fi−1/2,\begin{eqnarray*} (1 + \alpha C_{i-1/2} + C_{i+1/2})E_{\mathrm{r},i}^{n+1} - C_{i+1/2} E_{\mathrm{r},i+1}^{n+1} & = & E_{\mathrm{r},i}^{n} + \alpha C_{i-1/2}E_{\mathrm{r},i-1}^{n}\nonumber\\ & & + (1-\alpha)F_{i-1/2},\nonumber \end{eqnarray*}where α is an ad-hoc parameter, that gives the weight of each type of boundary conditions (0 ≤ α ≤ 1). Robin boundary conditions can be not only used as physical boundary conditions, but also as virtual boundary conditions when solving large matrix system with a parallel algorithm on subdomains (e.g., Tang 1992).

3.3. Implementation in RAMSES

We distinguish two types of interface: the coarse-to-fine and the fine-to-coarse. When the neighboring cells are at the same level, no special trick for the flux calculation is required. As we already mentioned, the fine levels are updated first in the RAMSES ATS algorithm. We define Δx as the size of the the grid mesh at level ℓ (the mesh size is uniform in all direction). In the following, we assume that ATS is used for all levels by default.

3.3.1. Fine-to-coarse interface

We consider the case where cell i − 1 is at level ℓ − 1. In this case, we use a simple Dirichlet boundary condition, where the radiative energy at the coarser level is set to be constant during the two time steps of level ℓ. To compute the radiative energy gradient at the interface, we also assume that the radiative energy is uniform within the neighboring coarse cell (no interpolation, see Fig. 1b), so that ˜Fi−1/2=−Ki−1/2Er,in+1−Er,i−1n32Δx·\begin{equation} \tilde{F}_{i-1/2} = -K_{i-1/2}\frac{E_{\mathrm{r},i}^{n+1}-E_{\mathrm{r},i-1}^{n}}{\frac{3}{2}\Delta x}\cdot \end{equation}(14)The contribution of Er,i−1n\hbox{$E_{\mathrm{r},i-1}^{n}$} is then moved to the RHS of the matrix system. At the end of the diffusion update, the flux ˜Fi−1/2\hbox{$\tilde{F}_{i-1/2}$} is stored at the interface boundary following Eq. (5), ensuring energy conservation.

3.3.2. Coarse-to-fine interface

In this case, we allow the use of the three different types of boundary conditions presented previously. As will be explained below, each has its pros and cons. In the case of the Neumann boundary condition, the flux that is imposed at the level interfaces is given by the flux that has been stored during the update of level ℓ + 1. This ensures energy conservation but can lead to negative energy problems (see Sect. 3.3.5). For the Dirichlet boundary type, a restriction operation is performed on the neighbor parent cell (oct) using the updated value of the leaf cells at level ℓ + 1, i.e., Er,i−1n+Δtℓ+1=Er,i−1n+Δt1ℓ+Δt2ℓ\hbox{$E_{\mathrm{r},i-1}^{n+\Delta t^{\ell+1}}=E_{\mathrm{r},i-1}^{n+\Delta t_1^{\ell}+\Delta t_2^{\ell}}$}. The Robin condition uses a mix between the stored flux and the updated neighbor energy, with parameter α being a user defined parameter. (1 − α) gives the relative amount of the energy that is conserved at the interface.

3.3.3. Energy loss with Dirichlet BC at coarse-to-fine interface

In the case of a Dirichlet boundary condition, the error made on energy conservation at the interface can be computed analytically. We assume that the radiative energy is uniform within the oct at level ℓ + 1, so that the restriction operation on the oct gives (Er,in+Δtℓ−1˜)ℓ−1=(Er,in+Δtℓ−1)ℓ\hbox{$(\tilde{E}_{\mathrm{r},i}^{n+\Delta t^{\ell-1}})^{\ell-1}=(E_{\mathrm{r},i}^{n+\Delta t^{\ell-1}})^\ell$}.

The flux that crosses the surface Sℓ = (Δxℓ)ndim−1 between cell i and cell i − 1 during the update of level ℓ is given by Eq. (5) Fi−1/2ℓ→ℓ−1=12K32Δxℓ(Er,in+Δtℓ+Er,in+2Δtℓ−Er,i−1n),\begin{equation} F^{\ell\rightarrow\ell-1}_{i-1/2}=\frac{1}{2}\frac{K}{\frac{3}{2}\Delta x^\ell}\left( E_{\mathrm{r},i}^{n+\Delta t^{\ell}}+ E_{\mathrm{r},i}^{n+2\Delta t^{\ell}}-E_{\mathrm{r},i-1}^{n}\right), \end{equation}(15)where we assume that Δt1ℓ=Δt2ℓ=1/2Δt1ℓ−1\hbox{$\Delta t_1^{\ell}=\Delta t_2^{\ell}=1/2\Delta t_1^{\ell-1}$} and that the diffusion coefficient K is constant. Similarly, the flux that crosses the same surface during level ℓ − 1 update equals Fi−1/2ℓ−1→ℓ=K2Δxℓ(Er,i−1n+2Δtℓ−Er,in+2Δtℓ).\begin{equation} F^{\ell-1\rightarrow\ell}_{i-1/2}=\frac{K}{2\Delta x^\ell}\left( E_{\mathrm{r},i-1}^{n+2\Delta t^{\ell}}-E_{\mathrm{r},i}^{n+2\Delta t^{\ell}}\right). \end{equation}(16)Energy conservation requires that Fi−1/2ℓ−1→ℓ=−Fi−1/2ℓ→ℓ−1\hbox{$F^{\ell-1\rightarrow\ell}_{i-1/2}=-F^{\ell\rightarrow\ell-1}_{i-1/2}$}. In our implementation, the energy mismatch is ΔEr=K6Δxℓ(Er,in+Δtℓ−Er,in+2Δtℓ+Er,in+Δtℓ−Er,i−1n)+K2Δxℓ(Er,i−1n+2Δtℓ−Er,i−1n).\begin{eqnarray} \Delta E_\mathrm{r}&=&\frac{K}{6\Delta x^\ell}\left( E_{\mathrm{r},i}^{n+\Delta t^{\ell}} - E_{\mathrm{r},i}^{n+2\Delta t^{\ell}} +E_{\mathrm{r},i}^{n+\Delta t^{\ell}}-E_{\mathrm{r},i-1}^{n}\right) \\ & & +\frac{K}{2\Delta x^\ell} \left( E_{\mathrm{r},i-1}^{n+2\Delta t^{\ell}} - E_{\mathrm{r},i-1}^{n} \right). \end{eqnarray}We see that three quantities contribute to the energy loss. First, it is proportional to the rate of change of the energy during the fine level updates (Er,in+Δtℓ−Er,in+2Δtℓ\hbox{$E_{\mathrm{r},i}^{n+\Delta t^{\ell}} - E_{\mathrm{r},i}^{n+2\Delta t^{\ell}} $}) and also directly to the first intermediate flux (Er,in+Δtℓ−Er,i−1n\hbox{$E_{\mathrm{r},i}^{n+\Delta t^{\ell}}-E_{\mathrm{r},i-1}^{n}$}). Second, it mainly depends on the energy change on the coarse cell itself during the coarse update (Er,i−1n+2Δtℓ−Er,i−1n\hbox{$E_{\mathrm{r},i-1}^{n+2\Delta t^{\ell}} - E_{\mathrm{r},i-1}^{n}$}). Energy conservation can then be strongly violated in case of large energy gradients and energy change (early time of an energy pulse propagation).

3.3.4. Implicit update

We follow the same iterative method as in Paper I to solve the implicit system of the coupled equations governing the radiative and gas internal energies using the two-temperature approach. As in Paper I the coupled system of equations is reduced to a single equation on the radiative energy thanks to the linearization of the emission term. The only difference is that the iterative method is called on a level-by-level basis after the hyperbolic update of each AMR level. We use a conjugate gradient algorithm with a diagonal preconditioner. The stopping criterion is based on the L2 norm of the residual r(j)/r0 < ϵconv (where r0 is the initial residual) and on the L∞ norm of the relative change of the radiative energy between two iterations (j) and (j − 1), i.e., max|Er(j)−Er(j−1)|/Er(j−1)<ϵconv\hbox{$|E_\mathrm{r}^{(j)}-E_\mathrm{r}^{(j-1)}|/E_\mathrm{r}^{(j-1)}<\epsilon_\mathrm{conv}$}. More details on the iterative solver and the two equations implicit system solver can be found in Paper I.

thumbnail Fig. 2

Illustration of the negative energy problem. The grid meshes are represented by the dashed lines and the energy gradient by the red curve. The stored flux at the level interfaces (dashed blue) are represented by the blue arrow. In this case, more energy goes out than comes in when updating the coarse level.

3.3.5. Limitations and comparison to other methods

The first limitation of our scheme is that it uses a fully implicit method, which is first-order in time so that it is generally dominated by the truncation error due to the time discretization. This could be improved in the future by using a Crank-Nicolson integrator scheme.

Secondly, Dirichlet boundary condition at the fine-to-coarse interface is a common approximation made by various authors (Howell & Greenough 2003; Zhang et al. 2011). It is a relatively crude approximation since it can lead to flux over-(under-) estimate at level interfaces. This flux is stored at the end of the fine level update in order to allow for energy conservation when updating the coarse level. Nevertheless, this flux, ˜Fi−1/2\hbox{$\tilde{F}_{i-1/2}$}, has been computed using a desynchronized value of the coarse level energy, so that it actually does not correspond to the correct flux Fi−1/2n+1\hbox{$F^{n+1}_{i-1/2}$}. Energy conservation is then ensured using artificial flux corrections. This inaccuracy can be improved using a multilevel solver as proposed by Howell & Greenough (2003). However, for the sake of simplicity, we have decided not to use this in the present work, a choice justified by the strong performance of our method in the tests below. Another drawback of flux storage at level interfaces is that the location where energy is stored is certainly not the correct one. Contrary to the case of the ATS method for the explicit hyperbolic solver, information can propagate across many cells during a single time step with the implicit scheme. For instance, energy could have been transported much further than the cell boundary to even coarser levels. On the other hand, we showed that using a Dirichlet boundary condition at the coarse-to-fine interface leads to unavoidable loss (or gain) of energy. This can be improved using the Robin boundary condition.

Sometimes we experienced severe problems using the flux conservation method in the case where a coarse level ℓ − 1 is surrounded by finer level ℓ with a gradient of energy following the sketch presented in Fig. 2. The flux stored at the left interface is lower that the one on the right. Using an implicit scheme can lead to large fluxes, and it is straightforward to see that energy conservation can lead to negative energy on level ℓ − 1 when it will be updated. (More energy goes out than comes in.) This is problem dependent, but unavoidable. Using Dirichlet boundary conditions avoids negative energies and is in that sense more robust.

thumbnail Fig. 3

Linear diffusion test. Top: radiative energy profile at time t = 2 × 10-13 for three calculations using α = 1 (black), α = 0.5 (red), and α = 0 (blue), i.e., using, respectively, Dirichlet, Robin, and Neumann boundary conditions at level interfaces. The analytical solution is represented by the dashed line. The right axis shows the AMR levels, i.e., the effective resolution profile (dotted line). Bottom: corresponding relative error as a function of the distance.

4. Tests

In this section, we perform a suite of numerical experiments to test and validate our method. We first test the diffusion operator using AMR, and then we perform a full RHD test.

4.1. Linear diffusion test

This first test is the same as in Paper I, Sect. 4.1. It consists in letting an energy pulse diffuse in a uniform medium. We only consider the radiative energy diffusion operator in the optically thick limit ∂Er∂t−∇·(c3ρκR∇Er)=0.\begin{equation} \frac{\partial E_\mathrm{r}}{\partial t} - \nabla \cdot\left(\frac{c}{3 \rho \kappa_\mathrm{R}} \nabla E_\mathrm{r}\right) = 0. \end{equation}(19)We consider a box of length L = 1. The initial radiative energy corresponds to a delta function, namely it is equal to 1 everywhere in the box, except at the center where it equals Er,L/2Δx = E0 = 105. We choose ρκR = 1. We apply periodic boundary conditions to ensure energy conservation. We use three levels of refinement with a coarse grid of 32 cells (effective resolution of 256 cells at the maximum level of refinement ℓ = 3). The mesh is refined when the radiative energy relative gradient exceeds 25% in a cell (e.g., if in cell i we have 2|Er,i − Er,i ± 1|/(Er,i + Er,i ± 1) > 0.25). Each additional level uses a time step twice as smaller as the coarser one. The coarsest time step is kept constant to Δt0 = 3.125 × 10-15 which gives a diffusion CFL of ~4 on the maximum level of refinement. We present three calculations using different values of α corresponding to different boundary condition at level interfaces, namely α = 1 for Dirichlet, α = 0.5 for Robin, and α = 0 for Neumann. The convergence criterion ϵconv is set to 10-8.

Figure 3, top panel, shows the radiative energy profiles (left axis) for the three calculations using α = 1 (black), α = 0.5 (red), and α = 0 (blue) at time t = 2 × 10-13, as well as the AMR levels (right axis, dotted line). The three calculations match the analytic solution (dashed line) remarkably well, even when energy is not conserved at the interfaces. Differences only appear at the center of the domain as illustrated in the bottom panel showing the corresponding relative error profiles. The errors remain of the order of a few per cent and below except in the tail of the diffusion patch. At this location, the increase in the relative error is the same in the three models as a consequence of the truncation error due to the time discretization (first-order). Figure 4 shows the total energy conservation as a function of time for the three calculations. As expected, energy is perfectly conserved in the case where α = 0, but increases by a few per cent (up to ~8%) in the cases where α ≠ 0. Using α = 0.5 results in better energy conservation compared the Dirichlet boundary condition. Finally, it is interesting to note that once the radiative energy profile becomes smoother (t > 10-12), the energy gain (or loss) is stabilized, indicating that energy conservation is much improved thereafter. This first test indicates that using Dirichlet boundary condition and ATS is reasonable even for extreme initial conditions such as a Dirac pulse. In the following, all the calculations we present have been run using the Dirichlet approximation when ATS is used.

thumbnail Fig. 4

Linear diffusion test: energy conservation as a function of time for the three calculations using α = 1 (black), α = 0.5 (red), and α = 0 (blue).

4.2. Equilibrium test with nonlinear diffusion coefficient

thumbnail Fig. 5

Equilibrium test with nonlinear diffusion coefficient. Top: radiative energy profile for the DTA (adaptive time step, black) and DTU (unique time step, red) models, and for the stationary analytic solution (dashed line). The right axis indicates the AMR levels of the DTA and DTU models (dotted line). Bottom: relative error profiles for the DTA (black), DTU (red) and uniform grid (256 cells, blue) models.

In this second test, we check that the introduction of subcycling does not change the global second-order accuracy in space of the scheme. Commerçon et al. (2011c) showed that the combination of AMR and a unique time step is globally second-order accurate in space.

We consider a uniform density (ρ = 1) within a box of L = 1. We impose two different radiation energies at the domain boundaries, i.e., Er(x = 0) = 4 and Er(x = 1) = 0.5. The initial radiative energy is initialized as a step function, i.e., Er,x < 0.5 = 4 and Er,x    >    0.5 = 0.5. The Rosseland opacity is a nonlinear function of Er, κR=1010Era\hbox{$\kappa_\mathrm{R}=10^{10}E_\mathrm{r}^{a}$} with a = 3/2. The system relaxes towards a steady state solution so that we can test the accuracy of the scheme without any limitation due to truncation errors in time. The analytic stationary solution Eana(x) is given by Eana(x)=[(Er(x=1)a+1−Er(x=0)a+1)x+Er(x=0)a+1]1/(n+1).\begin{equation} E_\mathrm{ana}(x)\!=\! \left[\left(E_\mathrm{r}(x=1)^{a+1}-E_\mathrm{r}(x=0)^{a+1} \right)x +E_\mathrm{r}(x=0)^{a+1} \right]^{1/(n+1)}. \end{equation}(20)The coarser grid comprises 8 cells, and we allow for up to four levels of refinement (ℓmax = 4, up to an effective resolution of 128 cells). The convergence criterion ϵconv is set to 10-8 and the mesh is refined where the radiative energy relative gradient exceeds 10%. We perform calculations using α = 1 and ATS (DTA model), the unique time step method of Commerçon et al. (2011c), DTU model, and for uniform grids ranging from 8 to 128 cells. In the DTA and DTU runs, we let vary the maximum level of refinement to achieve effective resolutions comparable to the uniform grid calculations.

Figure 5 (top) shows the radiative energy profile for the DTA (black solid line) and DTU (red solid line) calculations when four levels of refinement are allowed and the analytic equilibrium solution (dashed line). The DTA and DTU calculations give nearly identical results, showing that the stationary state has been reached. Figure 5 (bottom) shows the corresponding relative errors for the DTA and DTU models, and for a simulation run with a uniform grid of 256 cells. The relative errors is of the order of a few per cent, except close to the right boundary, where the energy gradients are the strongest.

Figure 6 shows the norm of the L2 error, calculated as L2=1Neff∑1Ncell[(Er,i−Eana(xi))Δxi]2,\begin{equation} {L}_2 = \sqrt{\frac{1}{N_\mathrm{eff}}\sum_1^{N_\mathrm{cell}}{\left[\left(E_{\mathrm{r},i} - E_{\mathrm{ana}}(x_i) \right)\Delta x_i\right]^2}}, \end{equation}(21)where Neff is the number of cells given by the effective resolution at the maximum level of refinement and Ncell is the number of cells in the calculations. The L2 norm is plotted as a function of the minimum grid size reached using AMR for the DTA (black square) and DTU (red triangle) models, and against the mesh size of the uniform grid models (blue cross). The first important result is that the scheme remains globally close to second-order accuracy in space even with AMR and ATS. The use of the AMR weakens the slope to first-order compared to the uniform grid one when more than two levels of refinement are used, since our numerical scheme is only first-order accurate in space at level interfaces. The ratio between the number of level interfaces and the total number of cells in the computational domain is high (1/7) and the errors are dominated by the one of the coarser level which explains that the second-order breaks. Second-order accuracy can be recovered when the coarse grid resolution is doubled (Guillet & Teyssier 2011). It is also worth mentioning that the error in the DTA model remains very close to the DTU one, which indicates that the error introduced by the lack of energy conservation at level boundaries is limited. There are only 28 cells in the AMR DTA and DTU models for the most resolved calculations (128 cells effective resolution).

thumbnail Fig. 6

L2 norm of the error as a function of the minimum grid size for the three models: DTA (black square), DTU (red triangle), and uniform resolution (blue cross). The dotted line gives the slope that is proportional to Δx (first-order accuracy in space) and the dashed line the (Δx)2 slope (expected second-order accuracy).

The calculations run with α = 0 crashed when we used more than two levels of refinement. This is due to negative energies that appear at the beginning of the calculations, when the mesh is refined only at the center and is thus in a situation similar to that of Fig. 2 close to the left box boundary. This shows the limitations of using Neumann boundary conditions at level interfaces. Last but not least, we see that the error only depends weakly on the method used to compute the flux at level boundaries, and the Dirichlet method, which does not conserve energy, gives remarkably good results.

4.3. Radiative shocks

thumbnail Fig. 7

Radiative shock test. Left: gas (blue line) and radiative (red line) temperature as a function of the distance to the shock in the subcritical radiative shock (M = 2) at time t = 0.1 (top panel). The semi-analytic solutions of the gas (square) and radiative (diamond) temperatures from Lowrie & Edwards (2008) are over-plotted. The AMR level is shown (right axis, dotted line). The bottom panel shows the corresponding relative error for the gas (blue) and radiative (red) temperatures. Right: same as left panels for the supercritical radiative shock (M = 5).

thumbnail Fig. 8

Spherical collapse test. Density, temperature, and velocity profiles as a function of the radius and distribution of temperature as a function of the density for two calculations using adaptive time-stepping (DTA, solid black line) and with a unique time step (DTU, dotted red line) when the central density is ρc ~ 1 × 10-10 g cm-3.

thumbnail Fig. 9

Boss & Bodenheimer test. Density (top row) and gas temperature (bottom row) maps in the equatorial plane at time t0 + 1.73 kyr for the DTA (left) and t0 + 1.69 kyr for the DTU (right), with t0 the time at which the density first exceeds 10-13 g cm-3. The black contours in the density maps shows the contour at ρ = 1012 cm-3 (~3.87 × 10-12 g cm-3) to identify the fragments.

thumbnail Fig. 10

Boss & Bodenheimer test. Distribution of the gas in the density-temperature plane for the DTA (left) and DTU (right) calculations and at the same time as in Fig. 9. The color coding indicates the mass contained within each bin in the density-temperature plane. The oblique lines show the iso-Jeans mass, ranging from 10-5  M⊙ to 10 M⊙.

Radiative shocks are good laboratories to test our RHD method. Classical analysis of radiative shocks can be found in Mihalas & Mihalas (1984). We choose initial conditions following Lowrie & Edwards (2008) who describe a semi-analytic method for the exact solution of radiative shock profiles with gray nonequilibrium diffusion in an optically thick medium. This setup has the advantage of resulting in a stationary shock and the semi-analytic solution can be directly compared with numerical results. We follow the initial setup of Zhang et al. (2011) for the sub- and super-critical radiative shock corresponding to Mach numbers M of 2 and 5, respectively.

The initial setup consists of a one dimensional region made of two uniform states which satisfy the jump relation for a radiating fluid in an optically thick medium (e.g., Mihalas & Mihalas 1984). The boundary conditions are imposed at the initial state values throughout the calculation time. We use an ideal gas equation of state, an adiabatic index γ = 5/3, a mean molecular weight μ = 1, and an optically thick medium, i.e., λ = 1/3. Matter and radiation are assumed to be initially in equilibrium, i.e., T = Tr. The Planck and Rosseland opacities are set to κP = 3.93 × 10-5 cm-1 and κR = 0.848902 cm-1. The initial grid is made of 32 cells and allows for maximum six additional levels of refinement. The refinement criteria are based on the gradients of the density and the radiative energy. We use the hll Riemann solver for the hydrodynamics with a CFL factor of 0.5. We use Dirichlet boundary conditions at level interfaces (α = 0). The convergence criterion ϵconv is set to 10-8.

For the subcritical shock (M = 2), the one dimensional region ranges from −1000 cm to 1000 cm and the discontinuity between the two initial states is located at x = 0 cm. The two states characteristics are: ρL = 5.45887 × 10-13 g cm-3, uL = 2.3545 × 105 cm s-1 and TL = 100 K for the left state, and ρR = 1.2479 × 10-12 g cm-3, uR = 1.03 × 105 cm s-1 and TR = 207.757 K for the right state. A cell is refined when the relative gradient of density or radiative energy exceeds 5%, which provides good resolution at the shock front and in the radiative precursor. Figure 7 (left) shows the gas temperature (squares) and radiative temperature (diamonds) profiles (top panel), and the corresponding relative errors (bottom panel). The semi-analytic solutions are also plotted (blue line for the gas temperature and red line for the radiative temperature). We note that we needed to shift slightly the position x = 0 since the shock moved a little from the initial x = 0 discontinuity during a short period of adjustment to reach the steady state as the shock structure develops from the initial step profile. The AMR level (right axis) is plotted in dotted line. Once the steady state is reached, only four levels of refinement are used (130 cells in the AMR hierarchy), for an effective resolution of 512 cells (3.9 cm). The agreement between the numerical and the analytic solutions is very good as indicated in the relative error plots which shows errors below 1% except at the location of the gas temperature spike. The shock is captured within only two cells thanks to the AMR capabilities. Compared to our previous unique time step method, the gain in CPU time is about a factor 50.

For the supercritical shock (M = 5), the one dimensional domain ranges from −4000 cm to 4000 cm and the discontinuity between the two initial states is placed at x = 0 cm. The left state values are identical to those of the subcritical shock and the right state values read: ρR = 1.964050 × 10-12 g cm-3, uR = 1.63    ×    105 cm s-1 and TR = 855.72 K. A cell is refined when the gradients of density and radiative energy exceed 1% and 10% respectively. Figure 7 (right) shows the profiles of radiative and gas temperatures (and the corresponding relative errors) following the same nomenclature as Fig. 7 (left). In this case, the effective resolution reached at the shock location is ~0.98 cm (2048 cells) and the total number of cell in the computational domain is 225. The numerical and semi-analytic solutions agree again very well with relative errors of a few per cent at most. Compared to the unique time step method, the gain in CPU time is about a factor ~30 in this case.

5. Application to star formation

Adaptive mesh-refinement is particularly well-suited for complex problems using deep hierarchies of levels such as those generated in collapse calculations. Moreover, in the star formation framework, radiative transfer plays a key role in the thermal behavior of the collapsing gas and alters dramatically the fragmentation of prestellar cores (e.g., Bate 2009; Offner et al. 2009; Commerçon et al. 2011b). The method we presented in Paper I has been used with success in star formation studies (Commerçon et al. 2010, 2012; Hincelin et al. 2013) and successfully compared to 1D spherical calculations of Commerçon et al. (2011a). Nevertheless, this method was time-consuming because of the unique time stepping. The improvement we present in this paper with ATS is thus of prime importance for star formation purposes.

5.1. Spherical collapse

We present in this section a spherical collapse test of an isolated and gravitationally unstable 1 M⊙ sphere of molecular gas, similar to the one presented in Paper I. We wish to compare the results obtained with ATS and with a unique time step. We consider a uniform density (ρ0 = 1.14 × 10-18 g cm-3) and temperature (T0 = 10 K) dense core. The ratio between the thermal and gravitational energies, αtherm = 5R0kBT0/2GM0μmH, is 0.5 (free fall time tff ~ 62.3 kyr), which gives an initial radius R0 ~ 0.024 pc. The ratio of specific heats is γ = 5/3 and the mean molecular weight is μ = 2.33. For the opacity, we use the gray opacities of Semenov et al. (2003) as tabulated in Vaytet et al. (2012) for homogeneous spherical dust grains and normal iron content in the silicates (Fe/(Fe+Mg) = 0.3). The coarsest grid is made of 323 cells and we use a refinement criterion based on the local Jeans length, which ensures that the latter is always resolved by at least eight cells. We use the Minerbo (1978) flux limiter, the hll Riemann solver, and a hydrodynamic CFL factor of 0.8. For the simulations with ATS, we use Dirichlet boundary conditions at level interfaces (α = 0). The convergence criterion ϵconv is set to 10-4 for the iterative solver. We run the calculations until the late evolution of the first Larson core (Larson 1969).

Figure 8 shows the density, temperature, and velocity profiles as a function of the radius, and the corresponding distribution of temperature as a function of density for the calculations with ATS (DTA, solid black line) and with a unique time step (DTU, dotted red line) when the central density is ρc ~ 1 × 10-10 g cm-3. The two calculations agree remarkably well for all the different quantities given that the hydrodynamic and the radiation solvers are subcycled in the DTA calculations. We only see a few differences at the tail of the radiative precursor in the temperature profiles (more extended in the DTU model). The first core mass and radius are respectively 6.24 × 10-2  M⊙ and 14.2 AU for the DTA and 6.37 × 10-2  M⊙ and 14.2 AU for the DTU. The acceleration in term of CPU time thanks to ATS is about a factor 25 and the calculations have been run on eight processors.

5.2. Boss & Bodenheimer test

This last test revisits a well-known numerical exercise on protostellar collapse. It is based on the early work by Boss & Bodenheimer (1979). It consists of the collapse of a uniform density and temperature sphere in solid body rotation with an azimuthal density perturbation of amplitude A. This type of test invokes many physical processes: hydrodynamics, gravity, and radiative transfer. In addition, the high nonlinearity of the problem tends to shorten the horizon of predictability so that the comparison between two methods is challenging (e.g., Commerçon et al. 2008). We choose the same ratio of thermal to gravitational energy as in the last section, αtherm = 0.5, a perturbation amplitude A = 0.1, and the angular velocity Ω0 is set by the ratio of rotational to gravitational energy βrot=R03Ω02/3GM0=0.4\hbox{$\beta_\mathrm{rot}=R_0^3\Omega_0^2/3GM_0 =0.4$}. The model has a high initial rotation, which favors the formation of a large disk that is prone to fragmentation. We use the same initial and numerical parameters as in the previous test, except for the refinement criterion which was increased to ten points per Jeans length. The maximum effective resolution that is reached is 131 072 cells (217, ~0.15 AU), corresponding to twelve additional levels of refinement. We ran two calculations, one with a unique time step (DTU) and another one with ATS (DTA). In the DTA calculations, the five first levels share the same time steps and the seven other use ATS.

Figure 9 shows the density and temperature maps for the two models, at time t0 + 1.73 kyr for the DTA and t0 + 1.69 kyr for the DTU (with t0 the time at which the density first exceeds 10-13 g cm-3). The qualitative agreement between the two calculations is good, given that in this comparison of methods, not only is the radiative solver subcycled, but also the hydrodynamics and the gravity solvers. The collapsing cores yield the same number of fragments (ten at this time). The mass within the fragments, i.e., where the density exceeds 1012 cm-3, is 0.098 M⊙ for DTA and 0.1 M⊙ for DTU. The biggest fragments have a mass of 3.45 × 10-2  M⊙ for DTA and 3.56 × 10-2  M⊙ for DTU. The temperature maps show also similar features, such as temperature spikes in the shocked region and heated regions around the fragments. Figure 10 shows the density-temperature distribution in the two calculations at the same time as in Fig. 9. The color coding indicates the mass contained within each bin in the density–temperature distribution. Again, the agreement between the two methods is good, in particular for the green area that represents most of the mass contained within or around the fragments. The typical Jeans mass of the fragments are also similar. In this last test, the gain in term of CPU time is about a factor 5.

6. Conclusion and future work

We have presented in this paper a new method for implicit solvers on adaptive mesh-refinement grids in the context of diffusion problems. The method can deal with an adaptive time-stepping strategy such as those used by hydrodynamical solvers. Our method has been successfully implemented in the RAMSES code for radiation hydrodynamics using the flux-limited diffusion approximation. The principle of this new solver is to consider each level of the AMR hierarchy independently from the others and to use simple recipes for the boundary conditions at level interfaces (Dirichlet, Neumann, and Robin conditions). We have demonstrated that each of the different boundary conditions has its pros and cons. In particular, the Neumann conditions, which ensures energy conservation, can lead to negative energies and errors in the deposit of energy that is stored at level interfaces. On the opposite, we showed that the simple Dirichlet condition is much more robust even if it does not strictly conserve energy. We tested our method against classical numerical exercises (diffusion test, radiative shocks) and compared our results with analytic solutions. The new method is close to second-order accuracy in space and the error only depends weakly on the type of boundary condition used at level interfaces. We applied the method to a star formation test case and successfully compared the new results to the ones obtained using the unique time step method presented in Paper I. The gain in CPU time can vary from a factor 5 to a factor 50, depending on the problem.

This new method makes use of a simple conjugate gradient algorithm as an iterative solver to integrate the diffusion operator. Since all the levels evolve independently from the others, we plan to allow the use of super-time-stepping (Alexiades et al. 1996; Commerçon et al. 2011c) and of explicit time integration depending on the ratio between the Courant condition for the diffusion and the one for the hydrodynamics. Concerning radiative transfer, the method presented in this paper is limited to gray radiation. Extension towards multigroup radiative transfer is in progress (e.g., Vaytet et al. 2012).

Last but not least, the implicit adaptive time-stepping can be applied to the study of astrophysical structures in which other diffusion-like problems such as the propagation of cosmic rays and the anisotropic electronic conduction are involved.

1

We assume that the radiative energy is uniform in the y and z directions for a given position on the x axis, i.e., plane-parallel approximation.

Acknowledgments

We thank the anonymous referee for his/her comments that have improved the quality of the paper. B.C. thanks Neil Vaytet and Patrick Hennebelle for useful comments and discussions. The work presented has been supported by the French ANR Retour Postdoc program.

References

  1. Alexiades, V., Amiez, G., & Gremaud, P.-A. 1996, Com. Num. Meth. Eng., 12, 12 [Google Scholar]
  2. Almgren, A. S., Beckner, V. E., Bell, J. B., et al. 2010, ApJ, 715, 1221 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bate, M. R. 2009, MNRAS, 392, 1363 [NASA ADS] [CrossRef] [Google Scholar]
  4. Boss, A. P., & Bodenheimer, P. 1979, ApJ, 234, 289 [NASA ADS] [CrossRef] [Google Scholar]
  5. Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2008, A&A, 482, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2010, A&A, 510, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Commerçon, B., Audit, E., Chabrier, G., & Chièze, J.-P. 2011a, A&A, 530, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Commerçon, B., Hennebelle, P., & Henning, T. 2011b, ApJ, 742, L9 [NASA ADS] [CrossRef] [Google Scholar]
  9. Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011c, A&A, 529, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Commerçon, B., Launhardt, R., Dullemond, C., & Henning, T. 2012, A&A, 545, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Guillet, T., & Teyssier, R. 2011, J. Comput. Phys., 230, 4756 [Google Scholar]
  12. Hincelin, U., Wakelam, V., Commerçon, B., Hersant, F., & Guilloteau, S. 2013, ApJ, 775, 44 [NASA ADS] [CrossRef] [Google Scholar]
  13. Howell, L. H., & Greenough, J. A. 2003, J. Comput. Phys., 184, 53 [Google Scholar]
  14. Larson, R. B. 1969, MNRAS, 145, 271 [NASA ADS] [CrossRef] [Google Scholar]
  15. Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [NASA ADS] [CrossRef] [Google Scholar]
  16. Lowrie, R. B., & Edwards, J. D. 2008, Shock Waves, 18, 129 [NASA ADS] [CrossRef] [Google Scholar]
  17. Mihalas, D., & Mihalas, B. W. 1984, Foundations of radiation hydrodynamics (New York: Oxford University Press) [Google Scholar]
  18. Minerbo, G. N. 1978, J. Quant. Spectr. Rad. Transf., 20, 541 [Google Scholar]
  19. Offner, S. S. R., Klein, R. I., McKee, C. F., & Krumholz, M. R. 2009, ApJ, 703, 131 [NASA ADS] [CrossRef] [Google Scholar]
  20. Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Tang, W. P. 1992, SIAM J. Sci. Stat. Comput., 13, 573 [CrossRef] [Google Scholar]
  22. Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. The Enzo Collaboration, Bryan, G. L., Norman, M. L., et al. 2013, ApJS, submitted [arXiv:1307.2265] [Google Scholar]
  24. Turner, N. J., & Stone, J. M. 2001, ApJS, 135, 95 [NASA ADS] [CrossRef] [Google Scholar]
  25. Vaytet, N., Audit, E., Chabrier, G., Commerçon, B., & Masson, J. 2012, A&A, 543, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Zhang, W., Howell, L., Almgren, A., Burrows, A., & Bell, J. 2011, ApJS, 196, 20 [NASA ADS] [CrossRef] [Google Scholar]

All Figures

thumbnail Fig. 1

Two-dimensional sketch of possible AMR grid configurations at the interface between two levels. Left: configuration 1 where cell i is at level ℓ and cell i − 1 at level ℓ + 1. Right: configuration 2 where cell i is at level ℓ and cell i − 1 at level ℓ − 1. In both sketches, the vertical red line represents the surface over which energy is exchanged when level ℓ is updated. Similarly, the horizontal dashed red line represents the gradient over which the flux is computed using the energy marked by a red star (the value is interpolated using the value represented by the attached blue circle).
In the text
thumbnail Fig. 2

Illustration of the negative energy problem. The grid meshes are represented by the dashed lines and the energy gradient by the red curve. The stored flux at the level interfaces (dashed blue) are represented by the blue arrow. In this case, more energy goes out than comes in when updating the coarse level.
In the text
thumbnail Fig. 3

Linear diffusion test. Top: radiative energy profile at time t = 2 × 10-13 for three calculations using α = 1 (black), α = 0.5 (red), and α = 0 (blue), i.e., using, respectively, Dirichlet, Robin, and Neumann boundary conditions at level interfaces. The analytical solution is represented by the dashed line. The right axis shows the AMR levels, i.e., the effective resolution profile (dotted line). Bottom: corresponding relative error as a function of the distance.
In the text
thumbnail Fig. 4

Linear diffusion test: energy conservation as a function of time for the three calculations using α = 1 (black), α = 0.5 (red), and α = 0 (blue).
In the text
thumbnail Fig. 5

Equilibrium test with nonlinear diffusion coefficient. Top: radiative energy profile for the DTA (adaptive time step, black) and DTU (unique time step, red) models, and for the stationary analytic solution (dashed line). The right axis indicates the AMR levels of the DTA and DTU models (dotted line). Bottom: relative error profiles for the DTA (black), DTU (red) and uniform grid (256 cells, blue) models.
In the text
thumbnail Fig. 6

L2 norm of the error as a function of the minimum grid size for the three models: DTA (black square), DTU (red triangle), and uniform resolution (blue cross). The dotted line gives the slope that is proportional to Δx (first-order accuracy in space) and the dashed line the (Δx)2 slope (expected second-order accuracy).
In the text
thumbnail Fig. 7

Radiative shock test. Left: gas (blue line) and radiative (red line) temperature as a function of the distance to the shock in the subcritical radiative shock (M = 2) at time t = 0.1 (top panel). The semi-analytic solutions of the gas (square) and radiative (diamond) temperatures from Lowrie & Edwards (2008) are over-plotted. The AMR level is shown (right axis, dotted line). The bottom panel shows the corresponding relative error for the gas (blue) and radiative (red) temperatures. Right: same as left panels for the supercritical radiative shock (M = 5).
In the text
thumbnail Fig. 8

Spherical collapse test. Density, temperature, and velocity profiles as a function of the radius and distribution of temperature as a function of the density for two calculations using adaptive time-stepping (DTA, solid black line) and with a unique time step (DTU, dotted red line) when the central density is ρc ~ 1 × 10-10 g cm-3.
In the text
thumbnail Fig. 9

Boss & Bodenheimer test. Density (top row) and gas temperature (bottom row) maps in the equatorial plane at time t0 + 1.73 kyr for the DTA (left) and t0 + 1.69 kyr for the DTU (right), with t0 the time at which the density first exceeds 10-13 g cm-3. The black contours in the density maps shows the contour at ρ = 1012 cm-3 (~3.87 × 10-12 g cm-3) to identify the fragments.
In the text
thumbnail Fig. 10

Boss & Bodenheimer test. Distribution of the gas in the density-temperature plane for the DTA (left) and DTU (right) calculations and at the same time as in Fig. 9. The color coding indicates the mass contained within each bin in the density-temperature plane. The oblique lines show the iso-Jeans mass, ranging from 10-5  M⊙ to 10 M⊙.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.