Max Duarte
email: max.duarte@cd-adapco.com

Richard R Dobbins
email: richard.dobbins@yale.edu

Mitchell D Smooke
email: mitchell.smooke@yale.edu

High order implicit time integration schemes on multiresolution adaptive grids for stiff PDE's

Keywords: High order time discretization, multiresolution finite volume scheme, stiff PDEs, implicit Runge-Kutta schemes, error control AMS subject classifications: 65M50, 65M08, 65G20, 65L04, 65L06, 65M20

We consider high order, implicit Runge-Kutta schemes to solve time-dependent stiff PDEs on dynamically adapted grids generated by multiresolution analysis for unsteady problems disclosing localized fronts. The multiresolution finite volume scheme yields highly compressed representations within a user-defined accuracy tolerance, hence strong reductions of computational requirements to solve large, coupled nonlinear systems of equations. SDIRK and RadauIIA Runge-Kutta schemes are implemented with particular interest in those with L-stability properties and accuracy-based time-stepping capabilities. Numerical evidence is provided of the computational efficiency of the numerical strategy to cope with highly unsteady problems modeling various physical scenarios with a broad spectrum of time and space scales.

Introduction

In many scientific applications, from biomedical models to combustion or air pollution modeling, stiff differential equations must be solved to carry out numerical simulations. A straightforward notion of stiffness was given by Hairer & Wanner by stating that stiff equations are problems for which explicit methods don't work [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF]. The latter is mainly due to the broad spectrum of physical or numerical time scales that a numerical solver must deal with whenever a stiff timedependent differential equation is being solved. Robust and stable methods are thus required for stiff problems in order to properly handle and damp out fast transients. In the past decades high order implicit Runge-Kutta schemes with excellent stability properties were developed and widely investigated to solve stiff problems modeled by ODEs (see [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF] § IV and references therein). The same methods can be naturally considered to solve stiff problems originating from time-dependent PDEs discretized in space. However, the high performance of implicit Runge-Kutta methods for stiff ODEs is often adversely affected by the size of the systems of nonlinear equations arising in the case of semi-discrete PDEs. In particular phenomena involving localized fronts, as considered in this work, commonly require fine spatial representations, hence potentially large systems of equations. Significant effort is required to achieve numerical implementations that solve the corresponding algebraic problems at reasonable computational expenses in terms of both CPU time and memory.

Low order implicit schemes were already successfully used to simulate very complex problems modeled by stiff PDEs. This is the case, for instance, for the numerical simulation of combustion flames accounting for detailed chemical kinetics and multi-species transport (see, e.g., [START_REF] Bennett | Computational and experimental study of the effects of adding dimethyl ether and ethanol to nonpremixed ethylene/air flames[END_REF][START_REF] Smooke | The computation of laminar flames[END_REF] and references therein). The computational performance of high order Runge-Kutta methods, implemented in well-established production aerodynamics codes, was also assessed in the context of laminar and turbulent compressible flows [START_REF] Bijl | Implicit time integration schemes for the unsteady compressible Navier-Stokes equations: Laminar flow[END_REF][START_REF] Carpenter | Fourth-order Runge-Kutta schemes for fluid mechanics applications[END_REF]. In particular Jacobian-free Newton-Krylov methods were investigated in conjunction with high order implicit schemes [START_REF] Isono | A Runge-Kutta-Newton-Krylov algorithm for fourth-order implicit time marching applied to unsteady flows[END_REF][START_REF] Bijl | Iterative solution techniques for unsteady flow computations using higher order time integration schemes[END_REF] to further reduce the computational requirements (see, e.g., [START_REF] Knoll | Jacobian-free Newton-Krylov methods: A survey of approaches and applications[END_REF] for a review on this subject). Similarly, high order space discretization schemes have also been implemented in this context, reducing in practice the computational stencils, hence the size of the nonlinear systems [START_REF] Noskov | An implicit compact scheme solver with application to chemically reacting flows[END_REF][START_REF] Noskov | An implicit compact scheme solver for twodimensional multicomponent flows[END_REF]. Easing the computational load is also achievable by designing efficient parallelization techniques as developed, for instance, in [START_REF] Tosatto | Parallelization strategies for an implicit Newton-based reactive flow solver[END_REF] for reactive flow solvers. Taking into account that grid adaptation techniques for unsteady problems disclosing localized fronts are specifically designed to yield high data compression, we exploit this capability here to efficiently implement implicit integration schemes for stiff PDEs. This strategy was already adopted, for instance, in [START_REF] Bennett | Local rectangular refinement with application to axisymmetric laminar flames[END_REF][START_REF] Bennett | Local rectangular refinement with application to nonreacting and reacting fluid flow problems[END_REF], together with low order implicit solvers.

Among the many adaptive meshing approaches developed in the literature, we consider in this work adaptive multiresolution schemes based on [START_REF] Harten | Adaptive multiresolution schemes for shock computations[END_REF][START_REF] Harten | Multiresolution algorithms for the numerical solution of hyperbolic conservation laws[END_REF], namely the multiresolution finite volume scheme introduced in [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF] for conservation laws. Besides the inherent advantages of grid adaptation, multiresolution techniques rely on biorthogonal wavelet decomposition and thus offer a rigorous mathematical framework for adaptive meshing schemes [START_REF] Cohen | Wavelet Methods in Numerical Analysis[END_REF][START_REF] Müller | Adaptive Multiscale Schemes for Conservation Laws[END_REF]. Consequently, not only approximation errors coming from grid adaptation can be tracked, but general and robust solvers can be implemented since the wavelet decomposition is independent of any physical particularity of the problem and accounts only for the spatial regularity of the discrete variables at a given simulation time. Adaptive multiresolution schemes have been successfully implemented for the simulation of compressible fluids (see, e.g., [START_REF] Brix | Adaptive multiresolution methods: Practical issues on data structures, implementation and parallelization[END_REF][START_REF] Domingues | Adaptive multiresolution methods[END_REF] and references therein), as well as for the numerical solution of time-dependent, parabolic [START_REF] Roussel | A conservative fully adaptive multiresolution algorithm for parabolic PDEs[END_REF][START_REF] Bürger | Fully adaptive multiresolution schemes for strongly degenerate parabolic equations in one space dimension[END_REF] and stiff parabolic [START_REF] Duarte | New resolution strategy for multi-scale reaction waves using time operator splitting, space adaptive multiresolution and dedicated high order implicit/explicit time integrators[END_REF][START_REF] Dumont | Simulation of human ischemic stroke in realistic 3D geometry[END_REF] PDEs. Nevertheless, to the best of our knowledge this is the fist attempt to implement high order implicit time integration schemes in the context of the adaptive multiresolution finite volume method to solve stiff PDEs.

The paper is organized as follows. We give in Section 2 a short introduction on multiresolution finite volume schemes and implicit Runge-Kutta schemes, in particular of SDIRK-(Singly Diagonally Implicit Runge-Kutta) and RadauIIA-type. Some key aspects of the numerical implementation of time implicit schemes on multiresolution grids are detailed in Section 3. Finally, the numerical solution of several stiff time-dependent PDEs is investigated in Section 4.

Numerical methodology

Let us consider a parabolic, time-dependent PDE,

∂ t u = F (∂ 2 x u, ∂ x u, u), t > t 0 , x ∈ R d , u(t 0 , x) = u 0 (x), t = t 0 , x ∈ R d , (1)
where u : R × R d → R m and F : R m → R m , for a model with m variables. For many physically inspired systems, the right hand side F can in general be written as

F (∂ 2 x u, ∂ x u, u) = F (u) = F 1 (u) + F 2 (u) + . . . , (2)
where the F i (u), i = 1, . . ., stand for different physical processes. For instance, a scalar nonlinear reaction-diffusion equation with u : R × R d → R would be given by

F 1 (u) = -∂ x • (D(u)∂ x u) and F 2 (u) = f (u)
for some diffusion coefficient, D : R → R, and a nonlinear function, f : R → R.

Multiresolution analysis

Without loss of generality we perform a finite volume discretization of problem (1) with m = 1 for the sake of simplicity. According to the multiresolution finite volume scheme [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF], we consider a set of nested dyadic grids over a computational domain Ω ⊂ R d as follows. Each cell Ω λ , λ ∈ S j , is the union of 2 d finer cells of equal size Ω µ , µ ∈ S j+1 . The sets S j and S j+1 are thus consecutive embedded grids over Ω, where j = 0, 1, . . . , J, corresponds to the grid-level, from the coarsest to the finest grid, i.e., j equal to 0 and J, respectively. We denote U j := (u λ ) λ∈Sj as the spatial representation of u on the grid S j , where u λ represents the cell-average of u : R × R d → R in Ω λ :

u λ := |Ω λ | -1 Ω λ u(t, x) dx, x ∈ R d . (3)
Data at different levels of discretization are related by two inter-level transformations: the projection and prediction operators, briefly defined in Appendix A. Based on these operations, a multiresolution analysis allows one to define a one-to-one correspondence between two consecutive grid-levels:

U j ←→ (U j-1 , D j ), (4)
where the D j array gathers the so-called details. The latter can be seen as estimators of the local spatial regularity of a given discretized function, in this case U j , and represent the information lost when coarsening the spatial grid, in this case from U j to U j-1 . By iteration of this decomposition, we get a multi-scale representation of U J in terms of

M J := (U 0 , D 1 , D 2 , • • • , D J ): M : U J -→ M J , (5)
and similarly, its inverse M -1 . This multi-scale transform amounts to a representation of U J in a wavelet space spanned by a biorthogonal wavelet basis. Further details can be found in [START_REF] Cohen | Wavelet Methods in Numerical Analysis[END_REF][START_REF] Müller | Adaptive Multiscale Schemes for Conservation Laws[END_REF]. While the transformation ( 5) is exact and the multi-scale representation can be performed back and forth, real computational benefit is achieved by introducing a thresholding operator, as shown in Appendix A. This operator basically discards cells of smooth regularity whose values can be recomputed within an accuracy tolerance whenever needed. As a result a multiresolution approximation U J is obtained. Defining the following normalized 2 -norm:

U J 2 2 := 2 -dJ λ∈S J (u λ ) 2 ,
which corresponds to the L 2 -norm of a piecewise constant function, it can be shown that [START_REF] Duarte | A numerical strategy to discretize and solve the Poisson equation on dynamically adapted multiresolution grids for time-dependent streamer discharge simulations[END_REF]

U J -U J 2 ≤ Cη MR , (6)
where η MR corresponds to an accuracy tolerance 1 . So much is true for steady problems. When solving time-dependent problems, the same behavior is expected in terms of numerical errors introduced by the multiresolution approximation. (The spatially adapted grid is fixed during a given time integration step.) The latter was mathematically proved for hyperbolic problems in an L 1 -norm for both classical and inhomogeneous conservation laws in [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF] and [START_REF] Hovhannisyan | On the stability of fully adaptive multiscale schemes for conservation laws using approximate flux and source reconstruction strategies[END_REF], respectively. Moreover, numerical evidence proves similar behaviors for timedependent, parabolic [START_REF] Roussel | A conservative fully adaptive multiresolution algorithm for parabolic PDEs[END_REF][START_REF] Bürger | Fully adaptive multiresolution schemes for strongly degenerate parabolic equations in one space dimension[END_REF][START_REF] Bendahmane | Adaptive multiresolution schemes with local time stepping for two-dimensional degenerate reaction-diffusion systems[END_REF] and stiff parabolic [START_REF] Duarte | New resolution strategy for multi-scale reaction waves using time operator splitting, space adaptive multiresolution and dedicated high order implicit/explicit time integrators[END_REF][START_REF] Duarte | A new numerical strategy with space-time adaptivity and error control for multi-scale streamer discharge simulations[END_REF][START_REF] Duarte | Time-space adaptive numerical methods for the simulation of combustion fronts[END_REF] PDEs.

Implicit Runge-Kutta schemes

Let us now consider problem (1) discretized on an adapted grid obtained by multiresolution analysis:

d t U = F(U), t > t 0 , U(t 0 ) = U 0 , t = t 0 . (7)
For the ease of reading we denote U J simply as U of size m × N , where N corresponds to the number of cells in the adapted grid and thus F : R m×N → R m×N .

Given a time step ∆t we consider an implicit Runge-Kutta (IRK) scheme of order p for the numerical integration of the semi-discretized problem [START_REF] Bihari | Multiresolution schemes for the numerical solution of 2-D conservation laws I[END_REF]. An s-stage Runge-Kutta scheme is in general defined through a set of arrays b, c ∈ R s , such that b = (b 1 , . . . , b s ) T and c = (c 1 , . . . , c s ) T , and a matrix A ∈ M s (R) such that A = (a ij ) 1≤i,j≤s . These coefficients define the stability properties and the order conditions of the method, and are usually arranged in a Butcher tableau according to c A b T .

In practice, given a set of arrays z 1 , . . . , z s ∈ R m×N , we have to solve the nonlinear system

   z 1 . . . z s    = A    ∆tF(t 0 + c 1 ∆t, U 0 + z 1 ) . . . ∆tF(t 0 + c s ∆t, U 0 + z s )    , ( 8
)
where A is a square block-matrix of size s × m × N built with the coefficients (a ij ) 1≤i,j≤s (see more details in Appendix B). The solution U(t 0 + ∆t) is then approximated by U 1 , computed as

U 1 = U 0 + s i=1 d i z i , d T := (d 1 , . . . , d s ) = (b 1 , . . . , b s )A -1 . (9)
If all the elements of the matrix of coefficients A are non-zero, then we say we are considering a fully IRK scheme [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF]. Moreover, if

a sj = b j , j = 1, . . . , s (10)
then the last stage corresponds to the solution U 1 and thus d T = (0, 0, . . . , 0, 1) in [START_REF] Bijl | Implicit time integration schemes for the unsteady compressible Navier-Stokes equations: Laminar flow[END_REF]. Methods satisfying [START_REF] Brix | Adaptive multiresolution methods: Practical issues on data structures, implementation and parallelization[END_REF] are called stiffly accurate [START_REF] Prothero | On the stability and accuracy of one-step methods for solving stiff systems of ordinary differential equations[END_REF] and are particularly appropriate for the solution of (stiff) singular perturbation problems and for differential-algebraic equations [START_REF] Hairer | Error of Runge-Kutta methods for stiff problems studied via differential algebraic equations[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF].

An IRK approximation amounts then to solving a nonlinear system of equations of size s×m×N . The latter can be achieved by considering a simplified Newton method for system [START_REF] Bijl | Iterative solution techniques for unsteady flow computations using higher order time integration schemes[END_REF], recast as

G(Z) := Z -∆tAF (Z) = 0, (11)
where Z := (z 1 , . . . , z s ) T and F (Z) := (F(t 0 + c 1 ∆t, U 0 + z 1 ), . . . , F(t 0 + c s ∆t, U 0 + z s )) T . The (k + 1)-th approximation of the solution Z is thus computed in two steps. First, we solve the following linear system for δZ k ∈ R s×m×N :

(Id s×m×N -∆tAJ)δZ k = -Z k + ∆tAF (Z k ), ( 12
)
where J is a a square block-matrix of size s × m × N consisting of s rows and s columns of size m × N given by the Jacobians J 0 := ∂ U F(t 0 , U 0 ). Then, the previous solution Z k is corrected according to

Z k+1 = Z k + δZ k . (13)
A standard way of initializing the iterative algorithm considers simply

z 0 i = U 0 , i = 1, . . . , s. (14)

SDIRK schemes

An alternative to solving a large nonlinear system of size s × m × N is to consider an SDIRK scheme where a ij = 0 for j > i, and with equal diagonal coefficients, i.e., a ii = γ, i = 1, . . . , s. In general, A-or L-stable SDIRK schemes can be built of order p ≤ s + 1 or p ≤ s, respectively (see [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF] § IV.6

for more details). However, the stage order q of these schemes, that is, the order achieved by one single stage, is limited to 1.

The main idea is thus to successively solve the s stages by considering an m × N -dimensional system at each stage, that is, for i = 1, . . . , s,

z i = ∆tγF(t 0 + c i ∆t, U 0 + z i ) + ∆t i-1 j=1 a ij F(t 0 + c j ∆t, U 0 + z j ), ( 15
)
where the second term in the right-hand side is already known at the current stage. Adopting the same simplified Newton technique, this time stage-wise, we have to solve the following linear system for δz k i ∈ R m×N :

(Id m×N -∆tγJ 0 )δz k i = -z k i + ∆tγF(t 0 + c i ∆t, U 0 + z k i ) + ∆t i-1 j=1 a ij F(t 0 + c j ∆t, U 0 + z j ). (16)
A simple way of initializing the iterative algorithm at each stage considers

z 0 1 = U 0 , z 0 i = U 0 + z i-1 , i = 2, . . . , s. (17)
We will denote as SDIRK2, SDIRK3, and SDIRK4, respectively, the second, third, and fourth order SDIRK schemes here considered (Butcher tableaux in Appendix C).

RadauIIA schemes

Fully IRK schemes with a number of stages below its approximation order can be built based on collocation methods [START_REF] Guillon | La résolution numérique des problèmes différentiels aux conditions initiales par des méthodes de collocation[END_REF][START_REF] Wright | Some relationships between implicit Runge-Kutta, collocation and Lanczos τ methods, and their stability properties[END_REF], together with the simplified order conditions introduced by Butcher [START_REF] Butcher | Implicit Runge-Kutta processes[END_REF]. In this case, the coefficients (b j , c j ) s j=1 correspond to the quadrature formula of order p such that 1 0 π(τ ) dτ = s j=1 b j π(c j ) for polynomials π(τ ) of degree ≤ p-1. Moreover, the coefficients in c and A, together with conditions for the stage order q, imply that at every stage i the quadrature formula ci 0 π(τ ) dτ = s j=1 a ij π(c j ) holds for polynomials π(τ ) of degree ≤ q -1. Depending on the quadrature formula considered, such as Gauss, Radau or Lobatto, different families of implicit Runge-Kutta methods can be constructed (for more details, see [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF] § IV.5).

In this work we consider the family of RadauIIA methods introduced by Ehle [START_REF] Ehle | On Padé approximations to the exponential function and A-stable methods for the numerical solution of initial value problems[END_REF], based on [START_REF] Butcher | Integration processes based on Radau quadrature formulas[END_REF], that consider Radau quadrature formulas [START_REF] Radau | Étude sur les formules d'approximation qui servent à calculer la valeur numérique d'une intégrale definie[END_REF] such that p = 2s -1 and q = s. These are Aand L-stable schemes that are stiffly accurate methods according to [START_REF] Brix | Adaptive multiresolution methods: Practical issues on data structures, implementation and parallelization[END_REF]. In particular we consider the third and fifth order RadauIIA schemes referred as Radau3 and Radau5 (Butcher tableaux in Appendix C). Note that, even though Gauss methods attain a maximum order of p = 2s [START_REF] Butcher | Implicit Runge-Kutta processes[END_REF][START_REF] Ehle | High order A-stable methods for the numerical solution of systems of DEs[END_REF], they are neither stiffly accurate nor L-stable schemes, which are both important properties for stiff problems. Approximations of lower order are obtained with Lobatto methods satisfying p = 2s -2 [START_REF] Butcher | Implicit Runge-Kutta processes[END_REF][START_REF] Ehle | High order A-stable methods for the numerical solution of systems of DEs[END_REF][START_REF] Chipman | A-stable Runge-Kutta processes[END_REF][START_REF] Axelsson | A note on a class of strongly A-stable methods[END_REF]. In particular the collocation methods with p = 2s -2 and q = s, known as the LobattoIIIA methods, yield stiffly accurate schemes, but these are only A-stable.

Numerical implementation

We now discuss some particular aspects concerning the numerical implementation. We consider the multiresolution procedure presented in [START_REF] Duarte | New resolution strategy for multi-scale reaction waves using time operator splitting, space adaptive multiresolution and dedicated high order implicit/explicit time integrators[END_REF]. For the sake of completeness some key details of this particular implementation will be first recalled, while more details and references can be found in [START_REF] Duarte | Méthodes numériques adaptatives pour la simulation de la dynamique de fronts de réaction multi-échelles en temps et en espace[END_REF].

Construction of multiresolution grids

The adapted grid is composed of a set of nested dyadic grids: S j , j = 0, 1, . . . , J, from the coarsest to the finest. They are generated by recursively refining a given cell depending on the local regularity of the time-dependent variables, measured by the details at a given time. These grids are implemented in a multi-dimensional Cartesian finite volume framework. Following [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF] a centered polynomial interpolation of accuracy order β = 2r + 1 is implemented for the prediction operator, computed with the r nearest neighboring cells in each direction; the procedure is exact for polynomials of degree 2r. Here we will only consider the case β = 3 with one neighboring cell per direction (r = 1) including the diagonals in multidimensional configurations. For instance, in the one-dimensional case [START_REF] Harten | Multiresolution algorithms for the numerical solution of hyperbolic conservation laws[END_REF] the latter is given by

u j+1,2k = u j,k + 1 8 (u j,k-1 -u j,k+1 ), u j+1,2k+1 = uf j,k + 1 8 (u j,k+1 -u j,k-1 ).
Higher order formulae can be found in [START_REF] Müller | Adaptive Multiscale Schemes for Conservation Laws[END_REF], while extension to multi-dimensional Cartesian grids is easily obtained by a tensorial product of the one-dimensional operator [START_REF] Bihari | Multiresolution schemes for the numerical solution of 2-D conservation laws I[END_REF][START_REF] Roussel | A conservative fully adaptive multiresolution algorithm for parabolic PDEs[END_REF]. In general the interpolation stencil is given by (2r + 1) d cells.

Data compression is achieved by means of thresholding by discarding the cells whose details are below a given tolerance (see [START_REF] Noskov | An implicit compact scheme solver for twodimensional multicomponent flows[END_REF]) and thus defining a compressed set Λ with the remaining cells. However, not all cells can be eliminated as this would prevent one from performing the multiresolution inter-grid operations. In particular all cells in the prediction interpolation stencils must be always available. Consequently, cells are gathered in a graded tree Λ , instead of Λ, that is, a data structure that satisfies the aforementioned conditions (see [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF] for more details). Notice that Λ ⊂ Λ and error estimates like (6) follow straightforwardly. Nevertheless, for the ease of reading we will keep the notation Λ in the following to refer to a graded tree. A graded tree-structure is hence used to represent data in the computer memory (see also [START_REF] Roussel | A conservative fully adaptive multiresolution algorithm for parabolic PDEs[END_REF]). Recalling the standard tree-structure terminology: if Ω µ ⊂ Ω λ belonging to consecutive grids, we say that Ω µ is a child of Ω λ and that Ω λ is the parent of Ω µ . We thus define the leaves L(Λ) of a tree Λ as the set of cells Ω λ , λ ∈ L(Λ), such that Ω λ has no children in Λ. Depending on the size of the computational domain, more graded trees may be needed. Therefore, cells are distributed in N R graded trees Λ r , r = 1, . . . , N R , where N R := N Rx N Ry N Rz , and N Rx , N Ry , and N Rz stand for the number of graded trees or roots per direction. The adapted grid is thus given by sets L(Λ r ), r = 1, . . . , N R , with a total number of cells: N L = NR r=1 #(L(Λ r )). If no adaptation is required, then the maximum number of cells will be N L = #(S J ) = N Rx N Ry N Rz 2 dJ , i.e., the size of the finest grid.

Input parameters for the multiresolution implementation are: the maximum grid-level J corresponding to the finest spatial discretization; the number of roots per direction N Rx , N Ry , and N Rz ; and the threshold parameter η MR , which defines the numerical accuracy of the compressed representations following (6).

Numerical function evaluations

Introducing the set I n L := {1, 2, . . . , N n L }, where N n L stands for the number of leaves during time t ∈ [t n , t n + ∆t n ], we define a bijective function h n : D(h n ) → I n L , with D(h n ) := NR r=1 L(Λ n r ).
The cells (Ω λ ) hn(λ)∈I n L correspond then to the adapted grid during the current timestep ∆t n , defined by the leaves of the tree representation. Considering again m = 1, the solution of the semi-discrete problem [START_REF] Bihari | Multiresolution schemes for the numerical solution of 2-D conservation laws I[END_REF] for t ∈ [t n , t n + ∆t n ] is similarly defined as U(t) = (u λ ) hn(λ)∈I n L , where u λ stands for the cell-average of variable u(t, x) in Ω λ according to [START_REF] Bendahmane | Adaptive multiresolution schemes with local time stepping for two-dimensional degenerate reaction-diffusion systems[END_REF]. The discrete function F(U) in ( 7) can be thus defined as

F(U) = (F λ (U)) hn(λ)∈I n L , where F λ (U) can be further decomposed into Φ λ (U)
and ω λ (U), coming from the discretization of differential operators and source terms, respectively.

In particular for a second order spatial discretization, considered in this work, the local source term ω λ (U) becomes ω(u λ ), that is, it is computed using the local cell-average values. During timestep ∆t n , the time-dependent problem (7) can be thus written at each cell Ω λ of the adapted grid as

d t u λ = F λ (U) = Φ λ (U) + ω(u λ ), t ∈ [t n , t n + ∆t n ], h n (λ) ∈ I n L , (18)
Φ λ (U) := |Ω λ | -1 µ |Γ λ,µ |Φ λ,µ , (19)
where the latter sum is made over all µ = λ such that the interface Γ λ,µ := Ω λ ∩ Ω µ is not trivial, i.e., over all the neighboring cells of Ω λ ; and Φ λ,µ accounts for the flux across each interface.

In the simplest (low order in space) schemes, the flux Φ λ,µ is typically a function of u λ and u µ only, while higher order schemes require considering additional cells. Without loss of generality, let us denote by R Φ [λ] the stencil required to compute fluxes associated with cell Ω λ . Here, we consider flux computation schemes such that all cells in R Φ [λ] belong to the same grid, that is, fluxes are computed on a locally uniform mesh. Problem ( 18) can be thus rewritten as

d t u λ = F λ (u λ ) λ∈R Φ [λ] , t ∈ [t n , t n + ∆t n ], h n (λ) ∈ I n L . (20)
The numerical integration of problem [START_REF] Bihari | Multiresolution schemes for the numerical solution of 2-D conservation laws I[END_REF] then involves evaluating function F λ in [START_REF] Duarte | Méthodes numériques adaptatives pour la simulation de la dynamique de fronts de réaction multi-échelles en temps et en espace[END_REF] for the N n L current cells. Moreover, for a given interface Γ λ,µ the following conservation property holds in a finite volume flux representation: Φ λ,µ + Φ µ,λ = 0. Computing Φ λ,µ for Ω λ amounts to evaluating also Φ µ,λ for the neighboring cell Ω µ . Let us denote Φ + λ,µ as the right flux for Ω λ and Φ - µ,λ as the left flux for Ω µ , along the direction normal to Γ λ,µ . Similarly, R + Φ [λ] stands for the stencil required to compute Φ + λ,µ and, naturally,

R - Φ [µ] ≡ R + Φ [λ]; we thus have that Φ - µ,λ = -Φ + λ,µ .
This property is thus exploited to save computations as fluxes are computed only once at each interface. The locally uniform grids are then defined by the stencil R + Φ [λ] enclosing the current leaf Ω λ . Ghost cells, computed according to the inter-grid prediction operation, are used whenever one cell in the current stencil is missing. These ghost cells are also added to the adapted grid at interfaces between cells of different sizes in order to compute numerical fluxes at the highest grid-level between two neighboring cells [START_REF] Roussel | A conservative fully adaptive multiresolution algorithm for parabolic PDEs[END_REF].

Notice that function h n is in practice used for indexation of leaves, identifying them regardless of their geometric layout. This is particularly useful to organize the computation of the entries of the Jacobian, as shown in Appendix D, and the linear system. All matrices are stored using a standard CSR (Compressed Sparse Row) format for sparse matrices.

Newton method and linear solver

The simplified Newton method to solve the nonlinear system (8) considers the linear system [START_REF] Butcher | Implicit Runge-Kutta processes[END_REF] for fully IRK schemes like the RadauIIA methods. System ( 12) is recast as

(∆t -1 Id s×m×N -AJ)δZ k = -∆t -1 Z k + AF (Z k ), (21)
mainly to avoid updating ∆tAJ when time step changes are required. Defining an accuracy tolerance η Newt , we consider the following stopping criterion for the iterative process:

δZ k 2 ≤ η Newt . (22)
Additionally, we define a convergence rate for the Newton solver as

Θ k = δZ k 2 δZ k-1 2 , k ≥ 1. ( 23
)
For the first iteration we set Θ 0 = δZ 0 2 /(2 max U 0 ). We also define a maximum number of Newton iterations k Newt,max , and inspired by [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF], computations are interrupted and restarted with a halved timestep in [START_REF] Duarte | A numerical strategy to discretize and solve the Poisson equation on dynamically adapted multiresolution grids for time-dependent streamer discharge simulations[END_REF], if any of the following happens:

• there is a k such that Θ k ≥ 1;

• for some k, we have that

(Θ k ) kNewt,max-k-1 δZ k 2 ≥ η Newt , (24)
where the left-hand side in ( 24) is a rough estimate of δZ kNewt,max-1 2 ;

• k Newt,max iterations have been performed and δZ kNewt,max-1 2 > η Newt .

Notice that if the timestep is halved, only diagonal entries in [∆t -1 Id s×m×N -AJ] need to be modified; however, the resulting new matrix must be again factorized.

In this work we have implemented the iterative GMRES method [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] to solve the linear system (21), with right-preconditioning based on an ILUT factorization [START_REF] Saad | ILUT: A dual threshold incomplete LU factorization[END_REF]. Considering a fixed Jacobian has the advantage that the factorization and preconditioning of matrix [∆t -1 Id s×m×N -AJ] needs to be performed only once, unless computations are restarted with a halved timestep. Notice that this is a purely algebraic problem that has completely lost any reminiscence of its original geometric layout, meaning that it is independent of the adapted grid generation or any other grid-related data structure or geometric consideration. Consequently, any linear solver could be used as a black box solver provided that it only needs the matrix entries and the right-hand side array as inputs. For an iterative linear solver like GMRES we define another accuracy tolerance, η LS , as stopping criterion. This tolerance is chosen such that η LS = κη Newt , with κ ≤ 1. In this work we consider, for instance, κ = 10 -2 , unless noted otherwise. If the linear solver is taking too many iterations to converge, noted henceforth as k LS,J iterations, we update the Jacobians in the Newton method. Re-factorization and preconditioning would be necessary in this case.

The same ideas apply for the numerical implementation of the SDIRK schemes, considering at each stage system

((∆tγ) -1 Id m×N -J 0 )δz k i = -(∆tγ) -1 z k i + F(t 0 + c i ∆t, U 0 + z k i ) + i-1 j=1 a ij γ F(t 0 + c j ∆t, U 0 + z j ), (25)
instead of [START_REF] Cohen | Wavelet Methods in Numerical Analysis[END_REF]. In principle the same block-matrix needs to be factorized for all Newton iterations at all stages. The stopping criterion [START_REF] Duarte | A new numerical strategy with space-time adaptivity and error control for multi-scale streamer discharge simulations[END_REF], the convergence rate (23), as well as the conditions for halving the time step are all applied stage-wise, that is, with δz k i 2 instead of δZ k 2 . Similarly, the Jacobian is updated in [START_REF] Dumont | Simulation of human ischemic stroke in realistic 3D geometry[END_REF] after k LS,J iterations.

Time-stepping strategy

Since we consider only A-stable IRK schemes in this work, we are uninhibited by stability issues in the choice of the timestep, which can be based solely on accuracy requirements. For some kinds of problems, a constant time step might be sufficient to capture the problem dynamics. However, more generally, an adaptive time-stepping could be considered in order to enhance the computational efficiency. In either case, the main goal is to define a time step ∆t such that the local error satisfies

U(t 0 + ∆t) -U 1 2 = C∆t p+1 ≤ η RK , (26)
where η RK is the desired accuracy tolerance for the p-th order IRK scheme. The advantage of higher order methods is that they can satisfy [START_REF] Ehle | High order A-stable methods for the numerical solution of systems of DEs[END_REF] with larger time steps than those achievable with conventional low order methods. A standard approach to time step control is based on numerically approximating the exact local error in [START_REF] Ehle | High order A-stable methods for the numerical solution of systems of DEs[END_REF], by considering a solution Û1 computed by a lower order method of order p < p (see, for instance, [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF]). In this way we use the computations at the n-th step to predict the local error at the next step,

err = U n -Ûn 2 ≈ Cn ∆t p+1 n , (27)
which defines a new time step,

∆t new = ∆t n η RK err 1/ p+1 , (28)
by assuming that η RK ≈ Cn+1 ∆t p+1 new with Cn+1 ≈ Cn . The next time step ∆t n+1 will be based on ∆t new if the current approximation error satisfies err ≤ η RK . Otherwise, the current n-th solution will be disregarded, and the same n-th step will be integrated again with ∆t new instead of ∆t n . The lower order approximations are defined in Appendix E.

Inspired by [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF], we define a safety factor ν k that depends on the current Newton iteration k, the current linear solver iteration k LS , and the maximum number of Newton iterations k Newt,max , as follows

ν k = ν × 2k Newt,max + 1 2k Newt,max + max(k, 0.5k LS ) , (29)
where ν > 0 is a standard safety factor close to 1. Here we typically consider ν = 0.9. For the SDIRK schemes, where more than one Newton solve is required per time step, k and k LS in (29) stand, respectively, for the maximum number of Newton and linear iterations performed within a given time step. The time step ∆t n+1 is thus defined as

∆t n+1 = min(ν k ∆t new , α∆t n ), (30)
where α > 1 limits the variation of successive time steps. Here we consider in general α = 1.5. However, if the computations were to be performed with a constant time step ∆t, we would consider the following time-stepping procedure ∆t n+1 = min(αν k ∆t n , ∆t), [START_REF] Guillon | La résolution numérique des problèmes différentiels aux conditions initiales par des méthodes de collocation[END_REF] which allows modifications on the chosen time step based on the performance of the Newton and linear solvers. In general the initial time step ∆t 0 should be set sufficiently small to account for potentially fast transients. The numerical accuracy of the time integration is defined by the user-provided tolerance parameter, η RK . The tolerance parameter for the Newton solver is set to a lower value: η Newt = κη RK , with κ < 1. In this way errors coming from both the Newton and linear solvers should remain smaller than those caused by the IRK scheme.

Numerical illustrations

We investigate the computational performance of the numerical strategy for three problems modeled by time-dependent stiff PDEs. In this work all the simulations were run on a standard laptop with an Intel Core i3 @ 2.27 GHz processor and a memory capacity of 1.8 GB.

The Belousov-Zhabotinski reaction

Let us consider the numerical approximation of a model for the Belousov-Zhabotinski (BZ) reaction, a catalyzed oxidation of an organic species by acid bromated ion (see [START_REF] Epstein | An Introduction to Nonlinear Chemical Dynamics[END_REF] for more details and illustrations). The present mathematical formulation [START_REF] Field | Oscillations in chemical systems. II. Thorough analysis of temporal oscillation in the bromate-cerium-malonic acid system[END_REF][START_REF] Gray | Chemical Oscillations and Instabilites[END_REF] takes into account three species: hypobromous acid HBrO 2 , bromide ions Br -, and cerium (IV). Denoting by a = [Ce(IV)], b = [HBrO 2 ], and c = [Br -], we obtain a very stiff system of three PDEs given by

∂ t a -D a ∂ 2 x a = 1 µ (-qa -ab + f c), ∂ t b -D b ∂ 2 x b = 1 ε (qa -ab + b(1 -b)) , ∂ t c -D c ∂ 2 x c = b -c,            (32)
where x ∈ R d , with real, positive parameters: f , small q, and small ε and µ, such that µ ε 1. In this study: ε = 10 -2 , µ = 10 -5 , f = 1.6, q = 2 × 10 -3 ; with diffusion coefficients: D a = 2.5 × 10 -3 , D b = 2.5 × 10 -3 , and D c = 1.5 × 10 -3 . The dynamical system associated with this problem models reactive, excitable media with a large time scale spectrum (see [START_REF] Gray | Chemical Oscillations and Instabilites[END_REF] for more details). The spatial configuration with the addition of diffusion involves propagating wavefronts with steep spatial gradients; in particular, two-dimensional spiral waves and three-dimensional scroll waves [START_REF] Duarte | New resolution strategy for multi-scale reaction waves using time operator splitting, space adaptive multiresolution and dedicated high order implicit/explicit time integrators[END_REF].

Numerical time integration errors

We consider problem (32) in a one-dimensional configuration with Neumann homogeneous boundary conditions, discretized on a uniform grid of 1024 cells over a space region of [0, 1]. A standard, second order, centered finite volumes scheme is employed for the diffusion term. No grid adaptation is considered here in order to assess only the numerical errors related to the time integration schemes. To obtain an initial condition, we initialize the problem with a discontinuous profile close to the left boundary; we then integrate in time until the BZ wavefronts are fully developed. Figure 1 shows the time evolution of the propagating waves for a time window of [0, 1]. In order to compute the local errors associated with the implicit solvers here considered, we define a reference solution for the resulting semi-discrete problem. The latter is chosen here as the solution obtained using the Radau5 scheme (51), computed with a fine tolerance: η RK = 10 -14 .

Starting from the solution at t = 0.5, Figure 2 shows the local errors associated with each IRK scheme for different time steps. Both tolerances for the Newton and the linear solver are set to η Newt = η LS = 10 -14 in these computations. Notice that the stiffest variable, a, is directly subject to a time scale given by the small parameter µ = 10 -5 . We are thus in practice interested in time steps larger than 10 -5 . Considering the stiff and non-stiff components of (32), a and c, respectively, we see the following numerical behavior. For the stiff variable (see Figure 2 (top)), local errors of O(∆t p+1 ) tend to O(∆t q+1 ) for relatively large time steps. For the non-stiff variable (see Figure 2 (bottom)), the order reduction goes from O(∆t p+1 ) to O(∆t q+2 ). Local errors for variable b (not shown), stiffer than c, also behave as the ones for the stiffest component, variable a. These results are consistent with the classical, theoretical bounds derived in [START_REF] Hairer | Error of Runge-Kutta methods for stiff problems studied via differential algebraic equations[END_REF] for stiff ODEs in singular perturbation form, that is, containing a small stiffness parameter given by µ in our case. These results highlight the importance of the stage order for IRK schemes and stiff problems. In this respect RadauIIA schemes perform better than SDIRK methods. The same can be said with respect to stiffly accurate schemes when comparing, for instance, SDIRK3 with Radau3; more accurate results are obtained with the latter. As a matter of fact, a well-known conclusion is that stiffly accurate schemes guarantee better accuracies for stiff problems [START_REF] Prothero | On the stability and accuracy of one-step methods for solving stiff systems of ordinary differential equations[END_REF][START_REF] Alexander | Diagonally implicit Runge-Kutta methods for stiff O.D.E.'s[END_REF][START_REF] Hairer | Error of Runge-Kutta methods for stiff problems studied via differential algebraic equations[END_REF]. Figure 2 also shows the error estimates err given by ( 27) for both SDIRK4 and Radau5 schemes. Notice that the actual local errors are bounded by err, which in particular overestimates them since err is computed using a third order, embedded scheme in both cases. Finally, it has to be remarked that higher order schemes perform better in terms of numerical accuracy than low order ones like the first order Euler method, even when order reduction appears and all methods show the same low order convergence.

Performance comparison

We now consider problem [START_REF] Hairer | Error of Runge-Kutta methods for stiff problems studied via differential algebraic equations[END_REF] in a two-dimensional configuration with Neumann homogeneous boundary conditions, using multiresolution analysis to adapt dynamically the spatial discretization grid. For the multiresolution analysis the following input parameters are considered: number of roots per direction, N Rx = N Ry = 1; maximum grid-level, J = 10; and accuracy tolerance, η MR = 10 -3 . The finest grid has a spatial resolution of 1024 × 1024 over a computational domain of [0, 1] × [0, 1]. We consider both SDIRK4 and Radau5 schemes with the following parameters: k Newt,max = 30, k LS,J = k Newt,max , and κ = 10 -1 , recalling that η LS = κη Newt = κ 2 η RK . The initial solution is taken at t = 2, when the spiral waves are fully developed (see Figure 3), and the PDEs are then integrated until t = 2.01. (See [START_REF] Duarte | New resolution strategy for multi-scale reaction waves using time operator splitting, space adaptive multiresolution and dedicated high order implicit/explicit time integrators[END_REF] for details on the initialization of this two-dimensional configuration.) The data compression, defined as the ratio in percentage between the active and the finest grids, is of about 15%. Figure 4 shows the evolution of time steps according to [START_REF] Gray | Chemical Oscillations and Instabilites[END_REF], considering α = 1.5 and ∆t 0 = 10 -4 at t = 2 for various accuracy tolerances: η RK between 10 -3 and 10 -6 . For this particular problem a roughly constant time step is attained, consistent with the quasi constant propagation speed of the wavefronts. Tables 1 and2 gather information on the performance of both solvers during the time window [2, 2.01]. As also seen in Figure 4, larger time steps for a given accuracy tolerance are used with the Radau5 scheme, even though both schemes consider a third order, embedded method to compute dynamically the integration time steps [START_REF] Epstein | An Introduction to Nonlinear Chemical Dynamics[END_REF]. For SDIRK4, increasing the accuracy of the Newton solver involves more iterations even when smaller time steps are considered, showing a rather low dependence on the time step size and a Newton solver piloted mainly by its accuracy tolerance. A different behavior is observed for Radau5 where smaller time steps involve roughly the same number of Newton iterations, regardless of the Newton accuracy tolerance, meaning that smaller time steps effectively improve the Newton solver. All this is a direct consequence of the initialization of the Newton solver; while for SDIRK4 the Newton solver at each stage is initialized using the previous stage solution and thus at some time within the current time step, this is not the case for the present Radau5 solver for which the larger the time step the worse the initial approximation.

In terms of the iterative linear solver, the number of iterations decreases considerably with smaller time steps even when tighter convergence tolerances are considered. This is a direct consequence of the better preconditioning of the more diagonal-dominant matrices in ( 12) and ( 16) for relatively small time steps. In terms of CPU time, following Tables 1 and2 we can see that SDIRK4 is approximately 4 to 5 times faster than Radau5. Updating the grid together with the multiresolution operations takes approximately 10 to 13 % for SDIRK4, whereas the time integration, 84-90 %. These numbers are within the range of values found in the literature for adaptive grid techniques (see, e.g., [START_REF] Duarte | Méthodes numériques adaptatives pour la simulation de la dynamique de fronts de réaction multi-échelles en temps et en espace[END_REF]). For Radau5 the multiresolution load goes down to 2 % with a roughly 98 % of the CPU time allocated to the time integration, showing a clear problem of performance. There are two main reasons why a straightforward implementation of Radau5 is not fully satisfactory. First of all, a better initialization of the Newton solver is required to improve the convergence rate of the linear solver, regardless of the solver considered. For example, in [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF] ( § IV.8) all stages are initialized by extrapolating from the previous time step and using an interpolation polynomial based on the quadrature order conditions. Even if an adaptive grid technique can considerably reduce the increase of data storage that the latter procedure involves for multi-dimensional PDEs, we still have to introduce additional operations to initialize grid points that were not present during the previous time step. The second problem is related to the size of the algebraic systems which are basically tripled in the case of Radau5. The latter heavily impacts the performance of the linear solver. In this particular implementation, most of the overload is related to the preconditioning ILUT solver which was implemented as a black box, contrary to the GMRES solver. A tailored ILUT solver implemented specifically for this data structure may have already improved its performance before considering parallel computing implementations.

Ignition model of diffusion flames

We now consider the mathematical model derived in [START_REF] Thévenin | Ignition dynamics of a diffusion flame rolled up in a vortex[END_REF] to investigate the ignition dynamics of a diffusion flame, formed while a reactive layer is being rolled-up in a vortex. The hydrodynamics is decoupled from species and energy transport equations by adopting a standard thermo-diffusive approximation, leading to a reaction-diffusion-convection model. A two-dimensional computational domain is considered where pure and fresh hydrogen at temperature T F,0 initially occupies the upper half part, while the remaining lower part of the domain is occupied by hot air at T O,0 . By defining a Schvab-Zeldo'vich variable Z and a reduced temperature θ given by

θ = T -T O,0 T F,0 -T O,0 , (33)
the mathematical model is given by a system of equations of the form [START_REF] Thévenin | Ignition dynamics of a diffusion flame rolled up in a vortex[END_REF]:

∂ t Z + v x ∂ x Z + v y ∂ y Z -∂ 2 x Z + ∂ 2 y Z = 0, ∂ t θ + v x ∂ x θ + v y ∂ y θ -∂ 2 x θ + ∂ 2 y θ = F (Z, θ),    (34)
F (Z, θ) = Da φχY O,0 1 -Z φτ + 1 χ (Z -θ) Z + τ χ (Z -θ) e (-τa/(1+τ θ)) ,
with physical constant parameters: Da = 1.65 × 10 7 , φ = 34.782608696, χ = 50, Y O,0 = 0.23, τ = -0.7, and τ a = 8, corresponding to T F,0 = 300 K and T O,0 = 1000 K. The velocity field (v x , v y ) is given by a single vortex centered on the planar interface between the two media, which varies strongly in time and space. Its tangential velocity is given by

v θ (r, t) = Re Sc r 1 -e -r 2 /(4 Sc t) , (35)
where r(x, y) stands for the distance to the vortex center (x 0 , y 0 ) = (0, 0), and with Reynolds and Schmidt numbers of Re = 1000 and Sc = 1, respectively. In Cartesian coordinates, the velocity of a counter-clockwise rotating vortex is thus given by

v x = y -y 0 r v θ , v y = - x -x 0 r v θ , r = (x -x 0 ) 2 + (y -y 0 ) 2 1/2 .
The physics of the phenomenon can be briefly described as follows. A rotating vortex is introduced immediately at t = 0. The resulting forced convection superposes to the diffusive mechanisms and accelerates the mixture of the gases. A diffusion flame then ignites along the contact surface of both media, taking into account the important difference of temperatures in those regions. Once the flame is completely ignited, it propagates outwards from the center of the computational domain. The complete phenomenon encompasses thus very different physical regimes like mixing, ignition, propagation, which can be characterized depending on the initial reactants configuration and on the imposed velocity field, as studied in detail in [START_REF] Thévenin | Ignition dynamics of a diffusion flame rolled up in a vortex[END_REF].

High order temporal approximations

We consider problem [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF] in a two-dimensional configuration with Neumann homogeneous boundary conditions, using multiresolution analysis to adapt dynamically the spatial discretization grid. The convective term is discretized in space using a standard first-order upwind scheme. As before the multiresolution analysis is parametrized as follows: number of roots per direction, N Rx = N Ry = 1; maximum grid-level, J = 10; and accuracy tolerance, η MR = 10 -3 . The finest grid has thus a spatial resolution of 1024 × 1024 over a computational domain of [-1, 1] × [-1, 1]. The model is simulated for a time window of [0, 1.5 × 10 -4 ], using all the time integration solvers previously described. In all cases the following parameters were chosen: k Newt,max = 30, k LS,J = k Newt,max , and κ = 10 -2 . For this highly unsteady problem the number of active grid cells increases from approximately 3% of 1024 2 for the initial inert configuration, up to 13% at the final time when the diffusion flame is fully ignited along the contact surface (see Figure 5). (See [START_REF] Duarte | Time-space adaptive numerical methods for the simulation of combustion fronts[END_REF] for further details on the initialization of this problem.) Figure 6 shows the time evolution of the maximum temperature T max throughout the computational domain. Notice that the initial T max corresponds to T O,0 = 1000 K, the hot air; however, the fuel is initially at a much lower temperature of T F,0 = 300 K and thus the local temperature changes are in fact even more dramatic. First, we consider a constant time step of ∆t = 10 -5 , but with the time-stepping procedure given by [START_REF] Guillon | La résolution numérique des problèmes différentiels aux conditions initiales par des méthodes de collocation[END_REF]. In all cases an initial time step of ∆t 0 = 10 -8 was considered, taking into account that the velocity field [START_REF] Harten | Adaptive multiresolution schemes for shock computations[END_REF] ] t [10 -4 ] step. A tolerance of η LS = 10 -5 is considered for the Newton solver. Figure 6 (top) clearly shows the difference between the approximations obtained with different discretization orders. As a reference solution we consider the one obtained with Radau5 and ∆t = 10 -6 . The first order Euler solver introduces an ignition delay of the order of one time step. This delay is subsequently corrected by increasing the order of the time discretization with the same time step. Considering the third order approximations, SDIRK3 and Radau3, the latter performs much better for this particular problem because of its L-stability capabilities, given the strong transients that the stiff system (34) models. As a matter of fact, the L-stable SDIRK2 performs also better than SDIRK3. The difference of quality of the approximations can be further assessed by considering the final temperature at t = 1.5 × 10 -4 , once the ignition process is achieved and the strongest transients resolved. Final T max goes from 2254.09 K for the Euler solver to 2227.96 K and 2224.74 K for SDIRK4 and Radau5, respectively. For comparison, the reference Radau5 solution yields 2224.95 K. With a smaller time step of ∆t = 10 -6 , the Euler solution still shows a considerable difference (see Figure 6 (bottom left)) with a final T max of 2219.54 K. All the other solvers with ∆t = 10 -6 yield solutions with less than 0.3 K of difference with respect to the Radau5 solution, except for SDIRK3 with about 1 K of difference.

SDIRK2 ∆t=10 -5 SDIRK3 ∆t=10 -5 SDIRK4 ∆t=10 -5 SDIRK4 η RK =10 -3
As an illustration, Figure 6 (bottom right) also depicts the time steps considered for the various SDIRK solvers. Recalling that the time-stepping strategy is actually influenced by the performance of the Newton and linear solvers, we can see the impact of strong physical changes during the ignition process as all solvers need to use at some point smaller time steps. Similar behaviors are observed for the Euler and Radau solvers (not shown). In particular it can be seen again how the non L-stable SDIRK3 is the most affected solver. Notice that for these numerical experiments we have chosen a relatively large k Newt,max as this allows for larger time steps according to (29); consequently, k LS,J is also large and the Jacobians are never recomputed. Time steps are thus reduced due to a bad convergence rate of the Newton solver; a more conservative lower k Newt,max prevents this bad convergence rates since time steps would not even attain 10 -5 . In general a careful tuning of parameters should be conducted with this "constant" time step strategy in order to get the best possible performance. This tuning can be highly problem-dependent. That is why a time-stepping strategy based on an accuracy tolerance is very convenient, as the time steps can be effectively adapted to the various physical scenarios within a prescribed accuracy, while reducing the importance of the many parameters related to the Newton and linear solvers. Time-adaptive solutions are shown in Figure 6 (bottom left) for SDIRK4 and Radau5 with η RK = 10 -3 . In terms of CPU time, with ∆t = 10 -5 the Euler solver takes approximately 3.6 minutes, compared to 5 and 14.9 for SDIRK4 and Radau5, respectively; but the physics simulated with the first order method diverges considerably from the right one. The time-adaptive SDIRK4 with η RK = 10 -3 takes approximately 4 minutes, becoming a very promising alternative to the cheaper but less accurate Euler scheme, especially if one takes into account that a more accurate Euler solver with ∆t = 10 -6 takes about 6 minutes.

Concluding remarks

We have considered high order, implicit integration schemes to solve stiff multi-dimensional PDEs on adaptive multiresolution grids. Such an adaptive technique yields highly compressed representations within a user-prescribed accuracy tolerance, considerably reducing the computational requirements of implicit Runge-Kutta schemes. In particular a competitive time-space adaptive strategy was introduced to simulate models involving different physical scenarios with a broad spectrum of time and space scales within a user-specified level of accuracy. By designing an appropriate procedure to evaluate functions and represent linear systems within the multiresolution data structure, we have implemented several implicit Runge-Kutta schemes of SDIRK-and RadauIIA-type. The resulting linear systems are completely independent of the grid generation or any other gridrelated data structure or geometric consideration. Solving the algebraic problems constitute then a separate aspect from the multiresolution analysis itself, while the same procedure remains perfectly valid for other space adaptive techniques. Three stiff models have been investigated to assess the computational performance of the numerical strategy in terms of accuracy and CPU time. The computational analyses have thus proved that stiff PDEs can be effectively approximated with high order time discretization schemes with very limited computational resources. In particular SDIRK schemes require roughly the same amount of memory than a standard, low order Euler method. More memory-demanding RadauIIA schemes can also be employed in conjunction with adapted grids; however, as previously discussed, further enhancements are required to achieve better computational performances. It was also shown that even in the presence of order reduction, high order schemes yield more accurate solutions than low order ones. The advantages of high order discretizations have been especially highlighted when dealing with highly unsteady problems. However, for problems of even larger size parallel computing capabilities must be developed within the current context to achieve overall satisfactory results. Additionally, high order space discretization schemes, well-suited for implicit schemes [START_REF] Noskov | An implicit compact scheme solver with application to chemically reacting flows[END_REF][START_REF] Noskov | An implicit compact scheme solver for twodimensional multicomponent flows[END_REF][START_REF] Dobbins | A fully implicit, compact finite difference method for the numerical solution of unsteady laminar flames[END_REF], could be also considered in conjunction with grid adaptation to further enhance the computational performance. These issues constitute particular topics of our current research.

A Details on multiresolution analysis

Defining Ω λ := Ω j,k , we denote |λ| := j if λ ∈ S j , while subscript k ∈ ∆ j ⊂ Z d corresponds to the position of the cell within S j . For instance, in Cartesian coordinates the univariate dyadic intervals in R are given by Ω λ = Ω j,k :=]2 -j k, 2 -j (k + 1)[, λ ∈ S j := {(j, k) s.t. j ∈ (0, 1, . . . , J), k ∈ Z}, [START_REF] Harten | Multiresolution algorithms for the numerical solution of hyperbolic conservation laws[END_REF] and the same follows for higher dimensions. Following [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF], the projection operator P j j-1 maps U j to U j-1 . It is obtained through exact averages computed at the finer level by

u λ = |Ω λ | -1 |µ|=|λ|+1,Ωµ⊂Ω λ |Ω µ |u µ . ( 37
)
As far as grids are nested, this projection operator is exact and unique [START_REF] Cohen | Wavelet Methods in Numerical Analysis[END_REF]. On the other hand, the prediction operator P j-1 j maps U j-1 to an approximation U j of U j . There are several choices to define P j-1 j , but two basic constraints are usually imposed:

1. The prediction is local, i.e., u λ depends on the values u µ on a finite stencil R λ surrounding Ω λ , where |λ| = |µ| + 1.

2. The prediction is consistent with the projection in the sense that

u λ = |Ω λ | -1 |µ|=|λ|+1,Ωµ⊂Ω λ |Ω µ | u µ ; (38)
i.e., one can retrieve the coarse cell averages from the predicted values:

P j j-1 • P j-1 j = Id. (39)
With these operators, for each cell Ω λ the prediction error or detail is defined as the difference between the exact and predicted values,

d λ := u λ -u λ , (40)
We can then construct as shown in [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF] a detail vector defined as D j := (d λ ) λ∈∇j , where the set ∇ j ⊂ S j is obtained by removing for each µ ∈ S j-1 one λ ∈ S j (Ω λ ⊂ Ω µ ) in order to avoid redundancy (considering ( 41)) and to get the one-to-one correspondence (4). Given a set of indices Λ ⊂ ∇ J , where ∇ J := J j=0 ∇ j with ∇ 0 := S 0 , the thresholding operator T Λ is such that leaves unchanged the components d λ of the multi-scale representation M J in (5) if λ ∈ Λ, and replaces it by 0 otherwise. Defining the level-dependent threshold values ( 0 , 1 , . . . , J ), the set Λ is given by λ

∈ Λ if d λ L p ≥ |λ| . (42)
Applying T Λ on the multi-scale decomposition M J of U J amounts then to building the multiresolution approximation U J := A Λ U J to U J , where the operator A Λ is given by

A Λ := M -1 T Λ M,
in which all details below a certain level of regularity have been discarded. The bound ( 6) is thus verified with the level-dependent threshold values: j = 2 d(j-J)/2 η MR , j = 0, 1, . . . , J.

B Details on implicit Runge-Kutta schemes

The solution U(t 0 + ∆t) of problem ( 7) is approximated by U 1 , computed as

g i = U 0 + ∆t s j=1
a ij F (t 0 + c j ∆t, g j ) , i = 1, . . . , s;

U 1 = U 0 + ∆t s j=1 b j F (t 0 + c j ∆t, g j ) , (43)
where the time dependence for F(U(t)) was added for the sake of clarity. Following the approach established in [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF], the set of arrays z 1 , . . . , z s , are defined such that z i = g i -U 0 , i = 1, . . . , s, and hence,

z i = ∆t s j=1
a ij F(t 0 + c j ∆t, U 0 + z j ), i = 1, . . . , s.

Therefore, knowing the solution z 1 , . . . , z s implies an explicit formula for U 1 in ( 44), for which s additional function evaluations are required. These extra computation can nevertheless be avoided if the matrix A is nonsingular, as seen in [START_REF] Bijl | Implicit time integration schemes for the unsteady compressible Navier-Stokes equations: Laminar flow[END_REF]. A standard iterative Newton solver for system [START_REF] Bijl | Iterative solution techniques for unsteady flow computations using higher order time integration schemes[END_REF] amounts to solving at each iteration a linear system that requires the inversion of the block-matrix:

  

Id m×N -∆ta 11 J 1 (z 1 ) . . .

-∆ta 1s J s (z s ) . . . . . . . . .

-∆ta s1 J 1 (z 1 ) . . . Id m×N -∆ta ss J s (z s )

   , (46)
where J i (z i ) := ∂ U F(t 0 + c i ∆t, U 0 + z i ), i = 1, . . . , s, stands for the Jacobian ∂ U F := ∂F(U)/∂U of size m × N , evaluated at the various inner stages. Following [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF], the simplified Newton solver [START_REF] Butcher | Implicit Runge-Kutta processes[END_REF] approximates all Jacobians in (46) as J i (z i ) ≈ J 0 := ∂ U F(t 0 , U 0 ), i = 1, . . . , s.

C Butcher tableau of IRK schemes

We consider a two-stage SDIRK scheme [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF] (Table II.7.2) given by

γ γ 1 -γ 1 -2γ γ 1/2 1/2 d 1 = 3γ -1 2γ 2 , d 2 = 1 2γ . ( 48
)
With γ = (2 ± √ 2)/2 the method ( 48) is of second order (p = 2, s = p) and L-stable. For γ = (3 + √ 3)/6 it becomes a third order scheme (p = 3, s = p -1), but is only A-stable. The

Figure 1 :

 1 Figure 1: One-dimensional BZ propagating waves for variables a (top left), b (top right), and c (bottom), at time intervals of 0.2 within [0, 1] from left to right.

Figure 2 :Figure 3 :

 23 Figure 2: Local L 2 -errors for stiff and non-stiff components, respectively, a (top) and c (bottom) using the SDIRK2, SDIRK3 and SDIRK4 schemes (left), and the Euler, Radau3 and Radau5 ones (right). Dashed lines of slopes 2 to 5 (top), 3 to 5 (bottom left), and 2, 4 and 5 (bottom right) are also depicted. Error estimates err given by (27) are indicated with red bullets (•) (top) for the SDIRK4 (left) and Radau5 (right) schemes.

Figure 5 :

 5 Figure 5: Two-dimensional ignition model. Temperature T deduced from (33) at t = 1.5 × 10 -4 (left) and the corresponding adapted grid for the zoomed region [-1, 0] × [-1, 0] (right).

Figure 6 :

 6 Figure 6: Evolution of the maximum temperature T max using a time step of ∆t = 10 -5 for the SDIRK (top left) and Radau (top right) solvers. The solution computed with Radau5 and ∆t = 10 -6 is depicted with a solid black line. Similarly, time-adaptive solutions based on a tolerance of η RK = 10 -3 are shown (bottom left). The various time steps considered for the SDIRK solvers are also illustrated (bottom right); in all cases, ∆t 0 = 10 -8 .

  or in terms of inter-level operations: d λ = u λ -P |λ|-1 |λ| • P |λ| |λ|-1 u λ . The consistency property (38) and the definitions of the projection operator (37) and of the detail (40) imply |λ|=|µ|+1,Ω λ ⊂Ωµ |Ω λ |d λ = 0.

Table 1 :

 1 Time integration with SDIRK4 for t ∈ [2, 2.01]: number of time steps, n; maximum time step used, max ∆t n ; maximum number of Newton iterations, max k; maximum number of GMRES iterations, max k LS ; CPU time in seconds.

	10 -2	η RK = 10 -3					10 -2	η RK = 10 -3
		η RK = 10 -4		SDIRK4			η RK = 10 -4		Radau5
		η RK = 10 -5						η RK = 10 -5
		η RK = 10 -6						η RK = 10 -6
	∆t					∆t	10 -3
							10 -4
	2	2.002	2.004	2.006	2.008	2.01	2	2.002	2.004	2.006	2.008	2.01
			t						t
	Figure 4: Time-stepping with different accuracy tolerances, η RK , for the SDIRK4 (left) and Radau5
	(right) schemes.
		η RK	n	max ∆t n	SDIRK4 max k (per stage) max k LS CPU time (s)
		10 -3 11 1.95 × 10 -3	13 (3)	15		171.93
		10 -4 15 1.01 × 10 -3	17 (4)	13		269.28
		10 -5 26 4.65 × 10 -4	17 (4)	11		472.14
		10 -6 45 2.51 × 10 -4	19 (4)	10		837.62

Table 2 :

 2

	η RK	n	max ∆t n	Radau5 max k max k LS CPU time (s)
	10 -3 10 2.56 × 10 -3	5	60	892.98
	10 -4 12 1.64 × 10 -3	6	56	1268.11
	10 -5 19 6.97 × 10 -4	6	35	2311.30
	10 -6 37 2.98 × 10 -4	5	21	3416.22

Time integration with Radau5 for t ∈ [2, 2.01]: number of time steps, n; maximum time step used, max ∆t n ; maximum number of Newton iterations, max k; maximum number of GMRES iterations, max k LS ; CPU time in seconds.

  radically changes during the first time

		2200	Euler SDIRK2						2200	Euler Radau3
			SDIRK3							Radau5
		2000	SDIRK4						2000
	max [K]	1600 1800	∆t =10 -5				max [K]	1600 1800	∆t =10 -5
	T							T
		1400							1400
		1200							1200
		1000							1000
			0.2	0.4	0.6	0.8	1	1.2		0.2	0.4	0.6	0.8		1	1.2
					t [10 -4 ]								t [10 -4 ]
		2000 2200	Radau5 η RK =10 -3 SDIRK4 η RK =10 -3 Euler ∆t=10 -6 Euler ∆t=10 -5					7 8 9 10
	T max [K]	1600 1800						∆t [10 -6	4 5 6
		1400							3
									2
		1200
									1
		1000							0
			0.2	0.4	0.6	0.8	1	1.2	0	0.2	0.4	0.6	0.8	1	1.2	1.4
					t [10 -4 ]

Bound (6) was similarly shown in[START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF] for both uniform and 1 -norm.

fourth order SDIRK4 scheme proposed in [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF] (Table IV

D Numerical computation of Jacobians

The Jacobian at a given stage s:

L , needs to be numerically approximated at each leaf of the adapted grid. The latter is done by considering for each i = 1, 2, . . . , N n L , λ = h -1 n (i), the following expression

with the one-dimensional array

and 1 µ (i) = 0, otherwise; δ µ is a small perturbation taken here as δ µ = (10 -16 × max(10 -5 , |u s µ |)) 0.5 , following [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF]. Recalling [START_REF] Duarte | Méthodes numériques adaptatives pour la simulation de la dynamique de fronts de réaction multi-échelles en temps et en espace[END_REF] and the conservation property, the latter perturbation is limited in practice to the stencil

) is used to compute the Jacobian entry at cell µ due to perturbed cell λ: J s hn(µ),hn(λ) . Notice that only cells contained in the adapted grid, i.e., the leaves, are perturbed to compute the Jacobians. Close to level interfaces where a ghost-cell Ω ν might be contained in the flux stencil, we perturb the latter with the same δ µ corresponding to its parent Ω µ ⊃ Ω ν , which is of course a leaf. Geometric proportions must then be taken into account when transferring fluxes between cells of different resolution.

E Embedded Runge-Kutta schemes

The error measures used in this work rely upon the lower order embedded IRK schemes introduced in [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF] for the SDIRK4 and Radau5 schemes. For SDIRK4, the error approximation err is computed from

with coefficients e T := (e 1 , . . . , e 5 ) = (23/6, 17/12, -125/4, 85/3, 1), which is accurate to third order (p = 3). In the case of Radau5, we consider

with K = 10 -1 , e T = K(-13 -7 √ 6, -13 + 7 √ 6, -1), which is also accurate to third order (p = 3). As pointed out in [START_REF] Kvaernø | Singly diagonally implicit Runge-Kutta methods with an explicit first stage[END_REF] neither of these error estimates has A-stability properties. A remedy was proposed in [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF] for their Radau5 solver, that is currently not implemented in our code.