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/00046361/201322858  
Published online  26 February 2014 
A fast, robust, and simple implicit method for adaptive timestepping on adaptive meshrefinement grids
^{1}
École Normale Supérieure de Lyon, CRAL, UMR CNRS 5574, Université Lyon
I,
46 Allée d’Italie,
69364
Lyon Cedex 07,
France
email:
benoit.commercon@enslyon.fr
^{2}
Laboratoire de radioastronomie (UMR 8112 CNRS), École Normale
Supérieure et Observatoire de Paris, 24 rue Lhomond, 75231
Paris Cedex 05,
France
^{3}
LESIA−Observatoire de Paris, CNRS, UPMC, Université ParisDiderot,
5 place Jules
Janssen, 92195
Meudon,
France
^{4}
Universität Zürich, Institut für Theoretische Physik,
Winterthurerstrasse
190, 8057
Zürich,
Switzerland
Received:
16
October
2013
Accepted:
17
December
2013
Context. Implicit solvers present strong limitations when used on supercomputing facilities and in particular for adaptive meshrefinement codes.
Aims. We present a new method for implicit adaptive timestepping on adaptive meshrefinement grids. We implement it in the radiationhydrodynamics solver we designed for the RAMSES code for astrophysical purposes and, more particularly, for protostellar collapse.
Methods. We briefly recall the radiationhydrodynamics equations and the adaptive timestepping methodology used for hydrodynamical solvers. We then introduce the different types of boundary conditions (Dirichlet, Neumann, and Robin) that are used at the interface between levels and present our implementation of the new method in the RAMSES code. The method is tested against classical diffusion and radiationhydrodynamics tests, after which we present an application for protostellar collapse.
Results. We show that using Dirichlet boundary conditions at level interfaces is a good compromise between robustness and accuracy and that it can be used in structure formation calculations. The gain in computational time over our former unique time step method ranges from factors of 5 to 50 depending on the level of adaptive timestepping and on the problem. We successfully compare the old and new methods for protostellar collapse calculations that involve highly non linear physics.
Conclusions. We have developed a simple but robust method for adaptive timestepping of implicit scheme on adaptive meshrefinement grids. It can be applied to a wide variety of physical problems that involve diffusion processes.
Key words: methods: numerical / stars: formation / radiation: dynamics / hydrodynamics
© 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 diffusionlike problems. This coupling of explicit and implicit solvers is relatively straightforward and wellstudied on uniform grids (e.g., Turner & Stone 2001), but becomes far more difficult on complex grids, such as those generated by adaptive meshrefinement (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 secondorder parabolic partial differential equations can be generalized as a diffusion problem following (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 (2)with Δx the size of the discretized mesh. The explicit CFL condition for diffusion equation scales as Δx^{2} and is thus far more stringent than the classical CFL derived for the stability of the hyperbolic system formed by the Euler equations (Δt_{hyp} < Δx/c_{s}, with c_{s} 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 timestepping (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 levelbylevel 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 diffusionlike equation in the RAMSES code for radiation hydrodynamics (RHD) using a twotemperature 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 fluxlimited 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 wellknown 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
Fig. 1 Twodimensional 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 fluxlimited 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 (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, E_{T} is the total energy E_{T} = ρϵ + 1/2ρu^{2} + E_{r} (ϵ is the internal specific energy), κ_{P} is the Planck opacity, E_{r} is the radiative energy, and P_{r} 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 timestepping 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, (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., . 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}Δt^{0}, where Δt^{0} is the time step of the coarser level. The scheme is only firstorder 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 that crosses the interface with a cell i at a coarser level ℓ − 1 during the two fine time steps is (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 E_{r}(6)The finite volume discretization in the x direction^{1} of Eq. (6) using an implicit scheme gives (7)where V_{i} is the volume of cell i, S_{i ± 1/2} the surface of exchange between cell i and cell i ± 1, K_{i ± 1/2} is the mean diffusion coefficient computed at the cell interface (e.g., K_{i ± 1/2} = (K_{i ± 1} + K_{i})/2). Equation (7) can be written in the form (8)where C_{i ± 1/2} = K_{i ± 1/2}S_{i ± 1/2}Δt/(ΔxV_{i}). Equation (8) forms a matrix system, Ax = b, where matrix A has to be inverted to get the new value of the radiative energy (the solution vector x). The C_{i ± 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 (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 (10)We define as the mean diffusion coefficient at cells interface, and . 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 In the previous implementation of Paper I, the cells at level ℓ + 1 are not interpolated at level ℓ and the C^{ℓ → ℓ + 1} coefficient equals 2^{3−ndim}A^{ℓ}/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 C_{i−1/2} = C^{ℓ → ℓ+1} or C_{i−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., . Equation (8) then reads where we moved the terms corresponding to E_{r, i−1} on the randhandside (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 F_{i−1/2}, i.e., . Equation (8) reads 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 . Equation (8) reads where α is an adhoc 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 coarsetofine and the finetocoarse. 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. Finetocoarse 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 (14)The contribution of is then moved to the RHS of the matrix system. At the end of the diffusion update, the flux is stored at the interface boundary following Eq. (5), ensuring energy conservation.
3.3.2. Coarsetofine 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., . 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 coarsetofine 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 .
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) (15)where we assume that and that the diffusion coefficient K is constant. Similarly, the flux that crosses the same surface during level ℓ − 1 update equals (16)Energy conservation requires that . In our implementation, the energy mismatch is 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 () and also directly to the first intermediate flux (). Second, it mainly depends on the energy change on the coarse cell itself during the coarse update (). 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 twotemperature 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 levelbylevel 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 L_{2} norm of the residual r^{(j)}/r_{0} < ϵ_{conv} (where r_{0} 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. More details on the iterative solver and the two equations implicit system solver can be found in Paper I.
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 firstorder 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 CrankNicolson integrator scheme.
Secondly, Dirichlet boundary condition at the finetocoarse 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, , has been computed using a desynchronized value of the coarse level energy, so that it actually does not correspond to the correct flux . 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 coarsetofine 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.
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 (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 E_{r,L/2}Δx = E_{0} = 10^{5}. 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 2E_{r,i} − E_{r,i ± 1}/(E_{r,i} + E_{r,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 Δt^{0} = 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 (firstorder). 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.
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
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 secondorder accuracy in space of the scheme. Commerçon et al. (2011c) showed that the combination of AMR and a unique time step is globally secondorder 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., E_{r}(x = 0) = 4 and E_{r}(x = 1) = 0.5. The initial radiative energy is initialized as a step function, i.e., E_{r,x < 0.5} = 4 and E_{r,x > 0.5} = 0.5. The Rosseland opacity is a nonlinear function of E_{r}, 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 E_{ana}(x) is given by (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 L_{2} error, calculated as (21)where N_{eff} is the number of cells given by the effective resolution at the maximum level of refinement and N_{cell} is the number of cells in the calculations. The L_{2} 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 secondorder accuracy in space even with AMR and ATS. The use of the AMR weakens the slope to firstorder compared to the uniform grid one when more than two levels of refinement are used, since our numerical scheme is only firstorder 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 secondorder breaks. Secondorder 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).
Fig. 6 L_{2} 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 (firstorder accuracy in space) and the dashed line the (Δx)^{2} slope (expected secondorder 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
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 semianalytic solutions of the gas (square) and radiative (diamond) temperatures from Lowrie & Edwards (2008) are overplotted. 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). 
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 timestepping (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}. 
Fig. 9 Boss & Bodenheimer test. Density (top row) and gas temperature (bottom row) maps in the equatorial plane at time t_{0} + 1.73 kyr for the DTA (left) and t_{0} + 1.69 kyr for the DTU (right), with t_{0} the time at which the density first exceeds 10^{13} g cm^{3}. The black contours in the density maps shows the contour at ρ = 10^{12} cm^{3} (~3.87 × 10^{12} g cm^{3}) to identify the fragments. 
Fig. 10 Boss & Bodenheimer test. Distribution of the gas in the densitytemperature 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 densitytemperature plane. The oblique lines show the isoJeans 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 semianalytic 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 semianalytic solution can be directly compared with numerical results. We follow the initial setup of Zhang et al. (2011) for the sub and supercritical 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 = T_{r}. 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}, u_{L} = 2.3545 × 10^{5} cm s^{1} and T_{L} = 100 K for the left state, and ρ_{R} = 1.2479 × 10^{12} g cm^{3}, u_{R} = 1.03 × 10^{5} cm s^{1} and T_{R} = 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 semianalytic 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}, u_{R} = 1.63 × 10^{5} cm s^{1} and T_{R} = 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 semianalytic 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 meshrefinement is particularly wellsuited 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 timeconsuming 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 (T_{0} = 10 K) dense core. The ratio between the thermal and gravitational energies, α_{therm} = 5R_{0}k_{B}T_{0}/2GM_{0}μm_{H}, is 0.5 (free fall time t_{ff} ~ 62.3 kyr), which gives an initial radius R_{0} ~ 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 32^{3} 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 wellknown 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 . 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 (2^{17}, ~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 t_{0} + 1.73 kyr for the DTA and t_{0} + 1.69 kyr for the DTU (with t_{0} 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 10^{12} 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 densitytemperature 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 meshrefinement grids in the context of diffusion problems. The method can deal with an adaptive timestepping strategy such as those used by hydrodynamical solvers. Our method has been successfully implemented in the RAMSES code for radiation hydrodynamics using the fluxlimited 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 secondorder 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 supertimestepping (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 timestepping can be applied to the study of astrophysical structures in which other diffusionlike problems such as the propagation of cosmic rays and the anisotropic electronic conduction are involved.
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
 Alexiades, V., Amiez, G., & Gremaud, P.A. 1996, Com. Num. Meth. Eng., 12, 12 [Google Scholar]
 Almgren, A. S., Beckner, V. E., Bell, J. B., et al. 2010, ApJ, 715, 1221 [NASA ADS] [CrossRef] [Google Scholar]
 Bate, M. R. 2009, MNRAS, 392, 1363 [NASA ADS] [CrossRef] [Google Scholar]
 Boss, A. P., & Bodenheimer, P. 1979, ApJ, 234, 289 [NASA ADS] [CrossRef] [Google Scholar]
 Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2008, A&A, 482, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2010, A&A, 510, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Commerçon, B., Audit, E., Chabrier, G., & Chièze, J.P. 2011a, A&A, 530, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Commerçon, B., Hennebelle, P., & Henning, T. 2011b, ApJ, 742, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011c, A&A, 529, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Commerçon, B., Launhardt, R., Dullemond, C., & Henning, T. 2012, A&A, 545, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Guillet, T., & Teyssier, R. 2011, J. Comput. Phys., 230, 4756 [Google Scholar]
 Hincelin, U., Wakelam, V., Commerçon, B., Hersant, F., & Guilloteau, S. 2013, ApJ, 775, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Howell, L. H., & Greenough, J. A. 2003, J. Comput. Phys., 184, 53 [Google Scholar]
 Larson, R. B. 1969, MNRAS, 145, 271 [NASA ADS] [CrossRef] [Google Scholar]
 Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Lowrie, R. B., & Edwards, J. D. 2008, Shock Waves, 18, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Mihalas, D., & Mihalas, B. W. 1984, Foundations of radiation hydrodynamics (New York: Oxford University Press) [Google Scholar]
 Minerbo, G. N. 1978, J. Quant. Spectr. Rad. Transf., 20, 541 [Google Scholar]
 Offner, S. S. R., Klein, R. I., McKee, C. F., & Krumholz, M. R. 2009, ApJ, 703, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tang, W. P. 1992, SIAM J. Sci. Stat. Comput., 13, 573 [CrossRef] [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 The Enzo Collaboration, Bryan, G. L., Norman, M. L., et al. 2013, ApJS, submitted [arXiv:1307.2265] [Google Scholar]
 Turner, N. J., & Stone, J. M. 2001, ApJS, 135, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Vaytet, N., Audit, E., Chabrier, G., Commerçon, B., & Masson, J. 2012, A&A, 543, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zhang, W., Howell, L., Almgren, A., Burrows, A., & Bell, J. 2011, ApJS, 196, 20 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1 Twodimensional 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 
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 
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 
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 
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 
Fig. 6 L_{2} 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 (firstorder accuracy in space) and the dashed line the (Δx)^{2} slope (expected secondorder accuracy). 

In the text 
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 semianalytic solutions of the gas (square) and radiative (diamond) temperatures from Lowrie & Edwards (2008) are overplotted. 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 
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 timestepping (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 
Fig. 9 Boss & Bodenheimer test. Density (top row) and gas temperature (bottom row) maps in the equatorial plane at time t_{0} + 1.73 kyr for the DTA (left) and t_{0} + 1.69 kyr for the DTU (right), with t_{0} the time at which the density first exceeds 10^{13} g cm^{3}. The black contours in the density maps shows the contour at ρ = 10^{12} cm^{3} (~3.87 × 10^{12} g cm^{3}) to identify the fragments. 

In the text 
Fig. 10 Boss & Bodenheimer test. Distribution of the gas in the densitytemperature 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 densitytemperature plane. The oblique lines show the isoJeans mass, ranging from 10^{5} M_{⊙} to 10 M_{⊙}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.