Max Duarte
email: max.duarte@ecp.fr

Stéphane Descombes
email: stephane.descombes@unice.fr

Christian Tenaud
email: christian.tenaud@limsi.fr

Sébastien Candel
email: sebastien.candel@ecp.fr

Marc Massot
email: marc.massot@ecp.fr

Time-space adaptive numerical methods for the simulation of combustion fronts

Keywords: Laminar flames, time operator splitting, space adaptive multiresolution, time adaptive integration, error control

come

Introduction

Combustion simulation raises essentially multi-scale problems defined by a competition between a variety of processes evolving at different rates in time and featuring a broad range of length scales [1]. Such simulations are mathematically stiff, due to the broad spectrum of time scales in the nonlinear chemical terms and to the steep, spatially localized gradients in the reaction fronts. Direct numerical simulations (DNS) ensuring high fidelity numerical solutions with a sufficiently fine resolution of all temporal and spatial scales, necessarily require important computational resources. One approach to overcome computational restrictions in terms of CPU time and memory storage capacities is thus founded on the development of efficient techniques taking full advantage of massively parallel computing architectures and growth of computer power (see, e.g., [2,3]). Another strategy consists in reducing the modeling complexity and associated computational effort by more efficient handling of the time and spatial scales, while guaranteeing reliable predictive capabilities. This is exemplified by chemical kinetics reduction methods [4,5], tabulation techniques [6,7], or artificial flame thickening schemes [8] which all focus on an improved handling of time and spatial scales. The spatial filtering exploited in large eddy simulations (LES) relies on similar ideas. In LES the subgrid scales are modeled to effectively reduce the scale diversity and allow simulations of industrial configurations or scientific problems which could not have been handled by direct numerical simulations. Couplings of LES with other techniques such as tabulated chemistry models have also been investigated in order to enhance LES capabilities [9,10]. Additionally, modern LES flow solvers exploit parallel capabilities of modern computer architectures to simulate physically and geometrically complex configurations (see, e.g., [START_REF] Staffelbach | High Performance Computing for Computational Science -VECPAR 2008[END_REF]12]).

In general, the design and practical implementation of the numerical scheme constitute key elements for successful simulations. However, while many numerical methods carefully designed to ensure accurate solutions remain suitable for academic problems or small scale applications, realistic simulations of complex problems need to exploit methods of lower algorithmic complexity. As an illustration explicit time integration techniques are often implemented in industrial and scientific codes, for their inherent parallelism capabilities. Computer power can thus compensate severe restrictions of these schemes when they are applied to stiff problems. In this general context, it is important to develop alternative and more efficient strategies for such problems. The resulting algorithms aim at reducing computational requirements by thoroughly enhancing the numerical method of solution. This approach can then be coupled with other modeling considerations, like those presented previously, and with a suitable exploitation of modern computing capabilities. For instance, implicit time integration methods are investigated for complex applications [13,14,15,16] 2 , and efforts are also being made on computationally less demanding hybrid methods which combine implicit and explicit schemes such as IMEX [17,18] or operator splitting [19,20] techniques.

In this context, either for DNS or for other types of simulations on single or multi-processor architectures, research in terms of numerical methods is concentrated on accelerating schemes, ensuring adequate resolution properties and a proper handling of temporal and spatial scales of the problem. Nevertheless, due to the computational costs of realistic simulations, the validation of numerical results is restricted either to qualitative analysis for complex applications, or to precise comparisons with analytical or small scale academic configurations. The framework of the present investigation is therefore set by a new solution paradigm that complements and further extends the general approach, by introducing accuracy tracking capabilities in the numerical strategies. The objective is to control approximation errors associated to the numerical solution while capturing the dominant scales, and hence to estimate the quality of the results independent of the complexity or dimension of the computational model. We propose in this article a novel method, based on a dedicated operator splitting scheme with adaptive time-stepping features and space multiresolution analysis, relying on dynamic accuracy estimates.

A splitting approach offers a reduced degree of computational complexity and thus straightforward implementations with respect to IMEX methods, which require additional stability and order conditions that combine all inner implicit-explicit schemes [21,22,23,24]. However, appropriate criteria must be introduced to efficiently decouple the physical phenomena via splitting and to control the so-called splitting errors [25]. Splitting methods [26,27] have been used in the literature for decades, and were widely implemented and exploited for combustion problems to overcome classical restrictions on computational resources (see, e.g., [28,29,30,[START_REF] Oran | Numerical Simulation of Reacting Flows[END_REF][START_REF] Schwer | [END_REF]33,34,35]). A nice example is given by the numerical strategy developed by Knio et al. in [36,37] for reactive flows in a low Mach formulation with detailed chemical kinetics and transport parameters. The splitting scheme introduced by these authors combines the dedicated stiff implicit multi-step VODE solver [38] for the chemical reaction term with a second order explicit RKC scheme [39,40] for the diffusion problem. In this way, important gains of computational efficiency were achieved with a splitting time step not limited by the stiff reactive scales and set according to the extended stability domain of the RKC solver (convective stability constraints are less restrictive). Moreover, Day & Bell introduced in [41] another efficient low Mach solver with an operator splitting method coupled with an AMR (Adaptive Mesh Refinement) technique [42,43]. The reaction problem was solved with VODE, whereas the convection-diffusion term was explicitly integrated on the adapted mesh with local CFL time steps set by the corresponding grid size [44,45]. With these bases, further developments in terms of algorithm implementation and parallel computing techniques led to the effective simulation of three-dimensional turbulent premixed flames with detailed chemistry (see, e.g., [46,47]), a remarkable achievement for laboratory-scale turbulent flames (see, e.g., [START_REF] Bell | Proc. Nat. Acad. Sci[END_REF][START_REF] Bell | [END_REF]).

Considering the state of art and these recent advances, one may note that splitting schemes favor the use of dedicated numerical solvers of different nature as well as straightforward coupling with other techniques, with important gains in computational efficiency. Nevertheless, there are some open issues related to the construction of splitting schemes and the interaction of splitting errors with those originating from the inner implicit-explicit solvers (the influence of the latter ones on the global integration error was numerically illustrated, for instance, in [36,37]). Specifically, a critical matter underlined in the literature is the lack of precise criteria to properly choose the splitting time steps according to the physical decoupling capabilities of the problem and for a given accuracy. Another question is the extension of these strategies to highly dynamic problems for which neither a constant nor a stability-based variable splitting time step is adequate, taking into account that the explicit schemes are intended to handle the slow, non-stiff part of the equations.

To address these problems, we have first proposed in [50] a new approach in the design of splitting schemes for propagating waves modeled by stiff reaction-diffusion systems, in which the time integration errors were uniquely related to the splitting errors, even for large splitting time scales, based on mathematical analyses carried out mainly in [51,52]. The underlying idea is to decouple time integration errors by choosing one-step and high-order dedicated methods with time-stepping features for the split subproblems, to independently handle and solve the cor-responding physical-numerical time spectra, and such that the corresponding numerical errors remain negligible when compared with the splitting ones. The global error is then controlled by the splitting time step, defined according to the physical decoupling capabilities of the phenomenon and hence independently of standard stability constraints associated to mesh size or stiff source time scales. Additionally, the splitting scheme was coupled in [50] with a dynamic mesh refinement technique based on multiresolution (MR) analysis [53,54,[START_REF] Müller | Adaptive Multiscale Schemes for Conservation Laws[END_REF], previously restricted to non-stiff applications in the literature. For a given semi-discretized problem, the MR mathematical background allows a better monitoring of numerical errors introduced by the compressed spatial representations with respect to the original uniform grid discretization. Secondly, we have introduced in [56] a theoretical framework for stiff reaction-diffusion systems to dynamically select and adapt the splitting time step. The latter is based on local error estimates throughout the time integration in order to achieve a given controlled level of precision, even for very stiff problems.

The purpose of the present contribution is to carefully combine these two key ingredients and extend their applicability to a more complex level of modeling for convection-reactiondiffusion problems potentially modeling highly unsteady physical processes. The resulting solution scheme constitutes a fundamental building block for most combustion simulations. It provides an efficient algorithm in terms of both memory storage and computational efficiency, which allows multidimensional simulations assuring a given error tolerance, fixed in advance by the user. It is worth underlining that the purpose and challenge of the present contribution is to track and control the time integration errors as well as those originating from the compressed spatial representations for a given finest spatial resolution. The latter should be able to capture the whole set of length scales present in the physical problem one wishes to solve. In this way, this paper introduces a novel algorithm which will yield accurate solution of convectionreaction-diffusion problems with a large number of degrees of freedom related to required fine space resolution, still on standard computing platforms before fully exploiting massively parallel architectures.

To assess the efficiency of the method in terms of error control and computational cost, we consider the numerical simulation of premixed and diffusion flames of interest in combustion applications. This study specifically considers laminar flames interacting with vortex structures including propagation of flame fronts and self-ignition processes of reactive mixtures. Hydrodynamics are decoupled from species and energy transport equations by adopting the standard thermo-diffusive approximation, so that the time-space resolution of the multi-scale reactiondiffusion-convection problem can be thoroughly evaluated. The latter aspects constitute the major development investment in terms of numerical schemes, as shown in [36,41,37], and a fundamental and validating first stage for the present tools, before considering more complex models or a complete low Mach number formulation.

The paper is organized as follows. Section 2 defines the general flame/vortex interaction problem that will be solved numerically. Such flame vortex interactions were extensively investigated to describe fundamental combustion processes (see, e.g., [START_REF] Marble | Adv. in Aerospace Science[END_REF][START_REF] Karagozian | [END_REF]59,60,61,62]), and they also have been examined in some clever experiments (see, e.g., [63,64,65,66,67,68,69]). An AMR technique was also recently coupled with a level-set method in [70] to describe flame/vortex interaction problems. Section 3 depicts the main features of the space-time adaptive technique. The mathematical formulation as well as two-and three-dimensional numerical simulations are reported in Section 4, for the propagation of premixed flames interacting with a counter-rotating vortex pair. New issues arise in Section 5 which concerns the numerical simulation of an ignition process while a reactive layer is being rolled-up in a vortex.

Laminar flames coupled with vortex structures

We consider in this study the numerical simulation of laminar flames interacting with vortex structures. For the sake of simplicity the chemistry is modeled by a global, single step, irreversible reaction given by

ν F F + ν O O → ν P P, (1)
where ν k , k = F, O, P, stand for the stoichiometric coefficients for the fuel F, the oxidizer O, and the combustion products P. The reaction rate is controlled by an Arrhenius law with a relatively high activation energy giving rise to thin reaction layers. The following standard modeling assumptions are also introduced throughout this study: The thermo-diffusive assumption 7 essentially decouples the velocity field computation from the determination of species mass fractions and temperature. Known solutions of the incompressible Navier-Stokes equations may then be imposed, and the problem is reduced to solving the following species and energy balance equations:

∂ t Y k + v • ∂ x Y k -D ∂ 2 x Y k = - ν k W k ρ ẇ, ∂ t T + v • ∂ x T -D ∂ 2 x T = ν F W F Q ρc p ẇ,                (2)
with x ∈ R d , where W k is the molar mass and Y k , the corresponding mass fraction for k = F, O, P.

Variable T accounts for the temperature, and the reaction rate of progress ẇ is related to the rate of consumption of fuel ẇF , oxidizer ẇO , and products ẇP by ẇ = -( ẇF

/ν F ) = -( ẇO /ν O ) = ẇP /ν P .
The heat release per unit mass of fuel Q is defined by

Q = h F +(ν O W O /ν F W F )h O -(ν P W P /ν F W F )h P ,
where h k , k = F, O, P, is the corresponding species enthalpy. If one also considers a non-reacting diluent, noted by index N, the following equation must be included in (2):

∂ t Y N + v • ∂ x Y N -D ∂ 2 x Y N = 0. (3)
Notice that by definition, the mass fractions verify that

Y F + Y O + Y P + Y N = 1.
In the thermo-diffusive approximation the velocity field v(x, t) is usually defined analytically and imposed into (2). Throughout this study we consider a two-dimensional viscous core vortex configuration, which features an azimuthal velocity of the form:

v θ (r, t) = Γ 2πr 1 -e -r 2 /4νt , (4)
where Γ denotes the vortex circulation; r(x, y), the distance to the vortex center; and ν, the kinematic viscosity. This velocity field satisfies the incompressible Navier-Stokes equations; it is divergence-free, and it has a viscous core with a typical dimension of R ν ≈ (νt) 1/2 . Inside the core velocity increases in a quasi linear fashion as a function of the radial distance, and the fluid rotates like a solid body:

v θ (r, t) ≈ Γr 8πνt , r < R ν . (5)
A fast decay occurs immediately outside the core, and at large distances the flow tends to that of ideal line vortex:

v θ (r, t) ≈ Γ 2πr , r ≫ R ν . (6)

Time-space adaptive technique

For the reaction-diffusion-convection system defined by (2), we introduce a general timespace adaptive strategy: space adaptation is ensured by a multiresolution decomposition, whereas the time integration is performed by a dedicated splitting scheme with dynamic splitting time steps. The main idea of this splitting method is to apply high-order and one-step dedicated schemes to the reaction, diffusion, and convection subproblems, treated in an independent way [50]. Each solver is then intended to handle the fastest physical-numerical scales associated to each subsystem. The solution of the complete problem is finally reconstructed from the previous solutions according to the splitting scheme. Different multi-scale phenomena are thus considered separately during a splitting time step that is defined in order to achieve a prescribed numerical integration accuracy. This is enforced by using dynamic error estimates computed by an adaptive splitting scheme. These estimates measure the physical decoupling capabilities of the problem and settle the splitting time steps for a given accuracy [56]. Finally, the multiresolution analysis yields spatially adapted mesh representations with important gains in CPU time and memory usage.

A second order Strang scheme S is considered for the general problem (2) [27]. The solution U(t + ∆t) at time t + ∆t is computed from the previous solution U(t) at t by

S ∆t U(t) = R ∆t/2 D ∆t/2 C ∆t D ∆t/2 R ∆t/2 U(t), (7)
with U = (Y F , Y O , Y P , Y N , T ) T , and splitting time step ∆t. The operators R, D, C correspond, respectively, to the numerical integration of the reaction, diffusion, and convection problems, performed independently and successively in the order indicated in (7). Adaptive time stepping is considered for all three operators yielding variable reaction (∆t R1 and ∆t R2 ), diffusion (∆t D1 and ∆t D2 ), and convection (∆t C ) substeps, within the corresponding splitting time steps (∆t/2 or ∆t). Moreover, the reaction is described by spatially decoupled systems of ODEs, solved independently point by point with different time steps according to the local reactive intensity.

According to the adaptive splitting method introduced in [56], we also compute the embedded and lower order shifted-Strang splitting scheme:

S ∆t U(t) = R (1/2-δ)∆t D ∆t/2 C ∆t D ∆t/2 R (1/2+δ)∆t U(t), ( 8
)
where the δ parameter is defined such that the local error estimate err = S ∆t U(t) -S ∆t U(t) L 2 remains valid even for large splitting time steps [56]. Both solutions (7) and (8) have the first reactive half-step in common, and the following operators are applied simultaneously to both intermediate solutions: (R ∆t/2 U(t), R (1/2+δ)∆t U(t)) T , by putting them together as if one was solving a set of variables of twice the original number. The splitting time steps are then dynamically computed by

∆t new = υ ∆t η S ∆t U(t) -S ∆t U(t) L 2 , (9)
in order to ensure a prescribed accuracy η for each splitting solution (7). A safety factor is also considered: 0 < υ ≤ 1, close to 1 (υ = 0.9 in this study). The solution (7) at t + ∆t is accepted if err < η, and the integration proceeds with ∆t = ∆t new , according to (9). Otherwise it is rejected and computed once again with the new splitting time step ∆t new . The reaction and diffusion problems are solved, respectively, by the dedicated high-order one-step solvers3 : Radau5 and ROCK4, as in [50]. Radau5 [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF] is a fifth order implicit Runge-Kutta method exhibiting A-and L-stability properties to efficiently solve stiff systems of ODEs, whereas the ROCK4 scheme [START_REF] Abdulle | [END_REF] is formally a fourth order stabilized explicit Runge-Kutta method with extended stability domain along the negative real axis, well suited to numerically treat mildly stiff parabolic operators. Both methods implement adaptive time stepping techniques to guarantee computations within a prescribed accuracy tolerance, η Radau5 and η ROCK4 , set in this case smaller than the splitting tolerance η: (η Radau5 , η ROCK4 ) < η. Similarly, an explicit highorder in time and space, one-step monotonicity preserving scheme: OSMP, developed by [75], is implemented as the convective scheme. The latter employs TVD criteria to prevent spurious oscillations around discontinuities or sharp spatial gradients, combined with monotonicity preserving (MP) constraints for non-monotone data to avoid extrema clipping and local loss of accuracy common to TVD schemes near extrema. In this way monotone resolution of discontinuities and high order accuracy away from discontinuities (including smooth extrema) can be simultaneously attained [75]. Considering its explicit character, standard CFL stability restrictions are imposed to substeps ∆t C within each splitting time step ∆t.

Although the dynamic step size selection is made within a prescribed accuracy tolerance for Radau5 and ROCK4, this is not currently the case for the convective scheme for which time stepping is based only on stability constraints. As indicated in [50] and also in previous works (see, e.g., [37]), the attention was focused on the numerical solution of stiff reaction-diffusion systems. Hence we have extended the numerical strategy to reaction-diffusion-convection problems in which the main constraint for the convective term is given by the small time steps resulting from stability conditions related to fine spatial discretizations. In this context, the convective scheme must be at least of second order so that the numerical errors of the convection problem also remain negligible with respect to the splitting errors. The practical implementation of the convective scheme is summarized in Appendix A.

A fully adaptive multiresolution technique based on [54] is then coupled with the previous dedicated time adaptive operator splitting strategy. In this way, considering a finite volume discretization for problem (2) on a fine grid S J , the latter defines a set of dyadic nested meshes S j on which problem (2) can be represented, for j = 0, 1, • • • , J, from the coarsest to the finest grid 4 . We denote by U J split the numerical solution of the semi-discretized problem associated to system (2), computed at some time t by the time adaptive splitting scheme (7) on the uniform grid S J . Additionally, U MR split corresponds to the solution obtained with the proposed time-space adaptive scheme also using (7) this time on a dynamic adaptive mesh generated by multiresolution analysis. Based on the multiresolution mathematical background [54], the following holds

U MR split -U J split L 2 ≤ Cε, ( 10
)
for some positive C, where ε is a threshold parameter that sets the accuracy of the spatially adapted representations. Smaller ε implies more refined (less compressed) and hence more accurate solutions. Even though a rigorous mathematical proof of (10) is not yet available for parabolic problems, this property was already checked numerically (see, e.g., [START_REF] Roussel | [END_REF]79,50]).

The numerical accuracy of the simulations can then be described by

U J qe -U MR split L 2 ≤ U J qe -U J split L 2 splitting error + U J split -U MR split L 2 multiresolution error , (11)
where U J qe corresponds to a reference quasi-exact solution of the fully coupled reaction-diffusionconvection problem (2), discretized on the uniform grid S J . Most of the time, the latter solution is not available or too expensive to compute. The time stepping procedure (9) guarantees local time integration errors of order η; it then limits the numerical error introduced during each time integration step and monitors the splitting errors, i.e., the first term in the right-hand side of [START_REF] Staffelbach | High Performance Computing for Computational Science -VECPAR 2008[END_REF]. Similarly, the multiresolution decomposition for a given solution in time ensures a compressed representation of accuracy ε, and hence this parameter characterizes the multiresolution errors, i.e., the second term in the right-hand side of ( 11) through (10). Both η and ε are user defined parameters in order to achieve a desired level of accuracy with respect to a fully coupled solution of the problem on an equivalent uniform grid. The spatial discretization errors, which are in practice difficult to evaluate unless an analytical solution is known, are therefore not included in [START_REF] Staffelbach | High Performance Computing for Computational Science -VECPAR 2008[END_REF] and are settled by the order of discretization of the spatial operators and by the mesh size of S J , which in turn is limited by the available computational resources. In this application, second and third order spatial discretizations were considered for the diffusion and convection subsystems, respectively, in (2). More details on these space discretization errors are given in Appendix E as a complement to §5 for a planar strained flame for which an analytic solution can be derived.

Propagation of premixed flames

In the framework of problem (2) that models laminar flames interacting with vortices, one may first study the performance of the proposed numerical strategy for the simulation of premixed flames in two-and three-dimensional configurations. The model under consideration is borrowed from a configuration investigated by Laverdant & Candel in [80]. Some of these results were previously announced in [START_REF] Duarte | FVCA VI Problems & Perspectives[END_REF], in a more general context without any detailed analysis.

Model formulation

We consider a square computational domain where a mixture of fuel and oxidizer lies in the lower half-plane, while products occupy the upper half-plane. Fast but finite rate kinetics give rise to, a thin premixed laminar flame initially located at the mid-plane. By defining the progress variable c(x, y, t):

c = T -T o T b -T o = Y Fo -Y F Y Fo -Y Fb , (12)
where subscripts ( ) o and ( ) b denote, respectively, variables in the fresh mixture zone and the burnt product zone, we can derive the following time dependent PDE in a two-dimensional configuration [80]:

∂ t ⋆ c + v x,⋆ ∂ x ⋆ c + v y,⋆ ∂ y ⋆ c -∂ 2 x ⋆ c + ∂ 2 y ⋆ c = Da (1 -c)e -T a /(T o (1+τc)) , ( 13
)
where Da is a Damköhler number, τ = T b /T o -1, and subscript ( ) ⋆ stands for dimensionless variables; T a corresponds to the activation energy coming from an Arrhenius equation for the chemical reaction rate. Derivation of ( 13) from the original system ( 2) is briefly presented in Appendix B.

The dimensionless tangential velocity induced by the viscous core vortex (4) becomes

v θ,⋆ (r ⋆ , t ⋆ ) = Re Sc r ⋆ 1 -e -r 2 ⋆ /(4 Sc t ⋆ ) , (14)
where the Reynolds and Schmidt numbers are defined by

Re = Γ 2πν , Sc = ν D . (15)
In cartesian coordinates, the velocity is thus given by

v x,⋆ = y ⋆ -y 0,⋆ r ⋆ v θ,⋆ , v y,⋆ = - x ⋆ -x 0,⋆ r ⋆ v θ,⋆ , (16)
for a counter-clockwise rotating vortex with radius

r 2 ⋆ = (x ⋆ -x 0,⋆ ) 2 + (y ⋆ -y 0,⋆ ) 2 , (17)
where (x 0,⋆ , y 0,⋆ ) is the center of the vortex.

Analysis and performance of the method for several numerical simulations

In this application we consider two dynamic counter-rotating vortices, each one modeled by ( 14), interacting with a premixed flame governed by (13) in a two-dimensional computational domain. All the simulations presented in this work have been performed on an AMD Shanghai 2.7 GHz processor with a memory capacity of 32 GB.

Data initialization and simulation parameters

We solve problem (13) with Neumann homogeneous boundary conditions in a two-dimensional computational dimensionless domain of [-1, 1] 2 , i.e., L = 1 in (B.10). The initial condition corresponds to a premixed flame in the limit of large activation energy [80]:

c =        e (y ⋆ -y 0,⋆ )/∆ ⋆ , y ⋆ ≤ y 0,⋆ , 1, y ⋆ > y 0,⋆ , (18)
where ∆ ⋆ is the preheat zone thickness. The modeling parameters appearing in equations ( 13), (14), and ( 18) are as follows: Da = 2.5 × 10 9 , T a = 20000 K, T o = 300 K, τ = 6.72, Re = 1000, Sc = 1, and ∆ ⋆ = 0.02. In this part we consider a counter-rotating pair of vortices computed by superposition of two independent vortices, each modeled by (16) with opposite signs, and centered at (-0.25, -0.5) for the counter-clockwise vortex, and at (0.25, -0.5) for the clockwise one. We thus take y 0,⋆ = -0.5 into (18). Notice that the superposition of velocity fields of type ( 14) results in a new divergence-free field, but it verifies only approximately the Navier-Stokes equations. Velocities as well as the maximum time step ∆t max C,i for the stability of the convective scheme are computed before each time integration of the convection problem and updated after two time steps ∆t C,i according to (A.3).

The adaptive splitting accuracy tolerance in ( 9) is set to η = 10 -3 , unless noted otherwise, with η ROCK4 = 10 -5 and η Radau5 = 10 -7 for the ROCK4 and Radau5 solvers. The third order OSMP scheme is employed for the convection problem with a stability CFL condition equal to 1. In the computations, the shifting parameter δ for the lower order Strang scheme ( 8) is taken as a constant and sufficiently large value 5 of 0.05. Regarding (9), the following constraint was additionally considered for the dynamic computation of the splitting time steps:

∆t new = min          α 1 ∆t, υ ∆t η S ∆t U(t) -S ∆t U(t) L 2          , (19)
with α 1 equal to 1.5. This procedure is implemented to avoid large variations of time integration steps, taking also into account that the adapted grid is fixed during each time step. The time domain of integration is given by t ⋆ into [0, 4 × 10 -3 ]. For the multiresolution analysis, ε = 10 -2 and ε = 10 -3 are chosen as multiresolution threshold values in the following illustrations. (We shall discuss and analyze later on the choice of the tolerance parameters, in particular in §5.) The finest grid corresponds to a spatial discretization of 1024 2 points, i.e., J = 10 as finest grid level.

Characterization of the numerical performance of the method

According to the definition of the progress variable c in (12), the fresh mixture is given by c = 0, whereas c = 1 corresponds to the burnt gases. Starting from the planar premixed flame (18) at y ⋆ = -0.5, with fresh gases in the lower part (blue zone in the figure), Figure 1 (top) shows the time evolution of c and the interaction of the two imposed vortices with the flame front. The velocity field generated by the vortices is shown in Figure 2 (left), and is characterized by high values with localized strong gradients for the selected value of vortex Reynolds number of Re = 1000. As a consequence, the fresh mixture is drawn up towards the hot region (red zone in Figure 1 (top)) in the center region, at a rate which is much faster than the normal burning velocity of the reaction front, whereas hot gases propagate faster in the outer zones around the vortex cores.

The contour lines in Figure 1 (middle) account for the spatial thickness of the flame in which fresh gases react and burn, and where an important numerical effort is usually required to precisely describe the phenomenon. In this configuration the flame thickness is reduced from about 0.05 in the standard planar configuration, to approximatively 0.025 in regions where the flame surface is strained by the locally high velocity gradients, as shown in Figure 2 (left). A spatial mesh of 1024 2 points involves approximatively 10 discretization points throughout the flame front, and generates a reasonably good numerical representation of the problem. From a numerical point of view, the latter issue imposes a fine spatial discretization for a localized structure 100 times smaller than the global scale of the computational domain, and naturally justifies an adaptive mesh refinement technique. Figure 1 (bottom) shows the corresponding adapted grids. The representation involves 7 levels of different spatial discretization where the finest regions coincide with the propagating front. The data compressions DC illustrated in Figure 2 (right) are defined as 1 minus the ratio between the number of cells on the adapted grid AG and those on the finest uniform grid FG (1024 2 in this case), expressing the whole as a percentage:

DC = 1 - AG FG × 100. ( 20
)
For ε = 10 -2 and ε = 10 -3 , no more than, respectively, 10 % or 15 % of the 1024 2 points are necessary to represent the flame front within the prescribed tolerance.

In order to verify that the accuracy of the computations is settled by the accuracy tolerances, one may define for problem (13), discretized on a uniform mesh of 1024 2 : • A quasi-exact reference solution c J qe , obtained with the Strang scheme ( 7) with a small and constant splitting time step of ∆t = 10 -7 ;

• The splitting solution c J split is now obtained with the adaptive splitting scheme (7), with (19) and accuracy tolerance of η = 10 -3 , computed also on the uniform grid; and

• The time-space adaptive solution c MR split , composed of the adaptive splitting technique and the multiresolution representation with a spatial resolution equivalent to 1024 2 points in the finest grid level J = 10.

Defining also the numerical errors following [START_REF] Staffelbach | High Performance Computing for Computational Science -VECPAR 2008[END_REF]:

E J split = c J qe -c J split L 2 , E J MR = c J split -c MR split L 2 , E MR split = c J qe -c MR split L 2 , (21)
corresponding, respectively, to the time adaptive splitting, space adaptive multiresolution, and time-space adaptive approximations, the following Table 1 summarizes these errors where the solutions on adapted grids c MR split , were reconstructed on the finest grid only to compute the errors. With this choice of parameters, it is observed that the global accuracy of the numerical strategy E MR split is ruled by the time integration approximation error E J split , which is related to the local error tolerance η. The latter global error will remain practically independent of the multiresolution errors E J MR as seen in Table 1, for these and smaller threshold values. Notice however that a thresholding ε of 10 -2 yields errors of O(10 -3 ), and hence a better performance than expected. The latter is nevertheless a problem dependent feature and a safer choice generally would be to set ε equal to η at most, in order to guarantee sufficiently accurate spatial representations for the numerical time integration of the problem. Concerning the multiresolution errors, the proportionality with respect to ε is roughly verified, taking into account the improved accuracy for ε = 10 -2 , as previously noted. In this way we observe a good behavior of the method in terms of control of the numerical accuracy; a more detailed analysis is nevertheless postponed to §5.2.3 for a more complex configuration.

Table 2 includes the CPU times for half the time domain of study: t ⋆ into [0, 2 × 10 -3 ], taking into account that the reference quasi-exact solution is expensive to compute. The splitting CPU time accounts for the cost reduction with respect to the quasi-exact solution that considers t a small time step of the order of the fastest numerical scale (the convective CFL constraint in this case). Notice that more effective strategies could be implemented to obtain the coupled reference solution, and the previous values should be taken as one possible numerical indicator. Additionally, the CPU times related to the time-space adaptive technique account for the gain issued from the compressed data representation, if one compares them with the splitting CPU time. In this case these gains are coherent with the corresponding data compressions achieved with each threshold value in Figure 2 (right), for instance, about 90 % and 88 % for a threshold value ε of, respectively, 10 -2 and 10 -3 . The total CPU time for the time-space adaptive technique for the whole time domain of study [0, 4 × 10 -3 ] was about 80.73 and 98.38 minutes for ε = 10 -2 and ε = 10 -3 , respectively. Finally, Figure 3 illustrates the adaptive time steps considered in this problem. The splitting time steps are practically the same for both multiresolution tolerances. This indicates an appropriate spatial representation of the propagating fronts on the corresponding adapted meshes for the given set of tolerances (η and ε), considering that the time integration is performed on an adapted but fixed grid. Otherwise, any deficiency would be reflected by the local error estimates, and thus by the splitting time steps issued from the time adaptive scheme. An initial value of ∆t = 10 -8 was chosen in order to cope with the sudden appearance of the velocity field with very high maximum values of about 4 × 10 5 , with Re = 1000 into (14). Additionally it can be observed that the introduction of (19) to limit the growth of the splitting time steps yields a smooth evolution of the time stepping, and furthermore ensures an appropriate spatial representation of this particular initialization with a highly varying velocity field during the transient phase. Taking into account that the global physics is controlled by the propagation of the flame, the splitting time steps evolve until a practically constant value of ∆t ≈ 10 -5 . Nevertheless, time adaptation is needed to handle fast variations in the beginning of the process, and for the final total combustion of the fresh gases at some unknown time. In this configuration, all gases are burnt by t ⋆ = 3.5 × 10 -3 . The disclosed behavior of the time stepping procedure thus follows the physics of this particular problem and moreover justifies the choice of a constant shift parameter δ. However in a general situation, the dynamic evaluation of δ, as described in [56], should be included. Figure 3 also shows that an important decoupling of time scales is possible, where the splitting time step ∆t is globally at least 10 times larger than the inner integration steps for the split reaction, diffusion, and convection problems. This naturally yields important gains of computational efficiency, always within a prescribed accuracy. The reaction and diffusion time steps are dynamically set based on the accuracy tolerances η Radau5 and η ROCK4 , and for each half splitting time step we represent in Figure 3 the averaged values of the inner reaction and diffusion substeps in order to obtain clearer representations. Reaction steps are of the order of ∆t R ≈ 7 × 10 -7 at the flame front (shown in Figure 3), and they progressively increase up to ∆t R = ∆t/2 away from the highly reacting area. Diffusion time steps are of the order of ∆t D ≈ 10 -6 . The reaction and diffusion substeps corresponding to the second splitting half-step (not represented in Figure 3) behave qualitatively similar to those during the first half-step. The convection time steps are computed by (A.5), based on the maximum stability time steps ∆t max C,i , which are illustrated in Figure 3. This convective step ranges from ∆t C ≈ 2 × 10 -9 in the beginning to ∆t C ≈ 2 × 10 -7 and then ∆t C ≈ 5 × 10 -6 , due mainly to the constraining high Reynolds number considered.

⋆ [10 -3 ] E J split [10 -2 ] E J MR [10 -3 ] E MR split [10 -2 ] ε = 10 -2 ε = 10 -3 ε = 10 -2 ε = 10 -3 0.
The method performance has been investigated in two-dimensional configurations, for which quasi-exact solutions can be computed and allow detailed analyses. However, the time-space adaptive technique is easily extended to three-dimensional configurations, as illustrated in Figure 4 for the interaction of a premixed flame with a toroidal vortex. Since similar conclusions can be drawn also from this case, the details are briefly presented in Appendix C for the sake of legibility.

Ignition of diffusion flames

In this section, we investigate the ignition dynamics of a diffusion flame interacting with a vortex. The mathematical model was envisaged by Thévenin & Candel in [START_REF] Thévenin | [END_REF]. In what follows we consider several physical configurations investigated in [START_REF] Thévenin | [END_REF] and conduct a both qualitative and quantitative study on the performance of the present time-space adaptive method. Some preliminary results on one particular configuration were recently described in [START_REF] Duarte | Ann. Res. Briefs[END_REF].

Model formulation

Let us consider a two-dimensional computational domain where pure and fresh hydrogen with mass fraction Y F,0 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 with an oxidizer mass fraction Y O,0 . A single vortex modeled by (4) and centered on the planar interface between the two media, is then introduced to transport and mix both reactants. By defining a Schvab-Zeldo'vich variable Z and a reduced temperature θ given by

θ = T -T O,0 T F,0 -T O,0 , (22)
the original system (2) can be written as a reduced system of equations of the form [START_REF] Thévenin | [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, θ),            (23)
with

F(Z, θ) = Da φχY O,0 1 -Z φτ + 1 χ (Z -θ) Z + τ χ (Z -θ) e (-τ a /(1+τθ)) , (24)
where physical constant parameters φ, χ, and τ in F(Z, θ) are explicitly given in Appendix D; and τ a = T a /T O,0 is the reduced activation temperature.

Numerical simulations

We consider a two-dimensional computational domain with initially separated fresh fuel and hot air. A single vortex modeled by ( 14) constitute the velocity field imposed to the system.

Data initialization and simulation parameters

We consider the numerical solution of ( 23)-( 24) with Neumann homogeneous boundary conditions in a dimensionless domain of [-1, 1] 2 , i.e., L = 1 in (B.10). The initial solution is given by

Z(x ⋆ , y ⋆ ) = θ(x ⋆ , y ⋆ ) = 1 2 1 + tanh(α y ⋆ -y 0,⋆ ) , (25)
where α = 200. Notice that according to (D.3) and (22), Z(x ⋆ , y ⋆ ) = θ(x ⋆ , y ⋆ ) = 1 in the upper half-plane for the fresh fuel, and Z(x ⋆ , y ⋆ ) = θ(x ⋆ , y ⋆ ) = 0 in the lower part for the hot air. Therefore, considering the time dependent equation for Z in (23), Z must be included into [0, 1] at any time. Additionally, θ ≤ 1 according to (22). The initial mass fractions and temperatures of the gases are given by Y F,0 = 1, Y O,0 = 0.23, T F,0 = 300 K, and T O,0 . The following values define the characteristic parameters (D.2) for the evaluation of the reaction rate (24) and the velocity field ( 14): Q/c p = 5 × 10 4 K, s = 8, Da = 1.65 × 10 7 , T a = 8000 K, Sc = 1, and Re. The velocity field is given by the vortex rotating in the counter-clockwise direction, defined by ( 16) and centered at (0, 0). We thus take y 0,⋆ = 0 into (25). In the following computations we will consider different values of the vortex Reynolds number Re and the air temperature T O,0 in order to characterize different physical scenarios according to [START_REF] Thévenin | [END_REF].

Regarding the numerical strategy and unless otherwise noted the adaptive splitting accuracy tolerance in ( 9) is set to η = 10 -3 , with η ROCK4 = 10 -5 , η Radau5 = 10 -7 , and the third order OSMP scheme with a stability CFL condition equal to 1. As in the previous computations and after some preliminary runs, we consider the same constant shift parameter δ, equal to 0.05, for the lower order Strang scheme (8) and the growth limiting procedure for splitting time steps given by (19). The grid adaptation was performed with ε = 10 -3 , unless otherwise noted, for a finest grid corresponding to a spatial discretization of 1024 2 points, i.e., J = 10 as finest grid level.

Three different ignition dynamics

In this part we reproduce some of the computations performed in [START_REF] Thévenin | [END_REF] with the proposed numerical strategy. The main goal of this example is to evaluate the capabilities of the time-space adaptive scheme to simulate different physical scenarios by means of a qualitative comparison with previous results in [START_REF] Thévenin | [END_REF]. Thévenin & Candel considered a standard alternate direction implicit technique [START_REF] Peyret | Computational Methods for Fluid Flow[END_REF] with dynamic time stepping ruled by both an advective and a chemical time step limitation such that the increase of the local product mass fraction is restricted to 1 % over one time step. Very small time steps were thus required during thermal runaway and the propagation phase for a space discretization typically of 500 2 points [START_REF] Thévenin | [END_REF]. In their paper these authors identified and described three different configurations with their corresponding ignition dynamics. Taking into account the various inputs of the model, two parameters were chosen and varied in the computations to switch from one regime to another: the vortex Reynolds number Re and the air temperature T O,0 . The upper half-plane is initially occupied by fresh fuel at T F,0 = 300 K, whereas the remaining lower half contains hot air at T O,0 . The counter-clockwise rotating vortex (16) centered on the planar interface is introduced immediately at t ⋆ = 0. The resulting forced convection superposes to the diffusive mechanisms and accelerates the mixture of the gases. As a consequence, a diffusion flame ignites along the contact surface of both media, taking into account the important difference of temperatures in those regions. Additionally if the velocity field is sufficiently strong, depending on Re, it will entrain fresh gases into the vortex core which will react with an intensity set by the mixing temperature of gases of about (T F,0 + T O,0 )/2. In general these locally lower temperatures result in a delayed ignition of the core unless air of sufficiently high temperature is initially considered. The complete behavior is clearly a function of the initial reactants configuration and of the imposed velocity field, as studied in detail in [START_REF] Thévenin | [END_REF]. The first configuration illustrated in Figure 5 corresponds to an ignition process with moderate air temperatures, that is T O,0 less than 1000 K. In particular T O,0 was set to 800 K in these computations with a vortex Reynolds number Re of 160. Instantaneous values of the heat release F given by ( 24) are shown in Figure 5. Note that the scale changes in each snapshot in order to get a better representation of the flame structure. The hot air temperature is enough to favor the ignition of a diffusion flame along the braids, after a long but finite time. Some time after, a diffusion flame is also ignited near the central vortex core and progressively consumes it in a completely independent way of the burning process in the braids. Two well-separated flames are hence generated. For this configuration the central core is characterized by well mixed gases and the resulting mixture temperature is low so that the ignition delay is quasi-infinite as a consequence of both relatively low initial temperatures and a small Reynolds flow value. This burning dynamics corresponds to a diffusion ignition mode. These results suitably reproduce the dynamics described in [START_REF] Thévenin | [END_REF], even for the corresponding maximum heat release. Figure 5 shows also the dynamic variation of the splitting time step ∆t for the given accuracy tolerance of η = 10 -3 according to the physics of the problem, illustrated in this example by the thermal runaway and the rejection of the corresponding splitting time steps. The maximum temperature (max T ) is derived out of θ into (22). More precise analyses on the performance of the method will be conducted in a particular configuration in the next section. This simulation takes approximately 10 minutes of CPU time for t ⋆ up to 10 -3 . Figure 6 illustrates a configuration with an air temperature T O,0 of 2000 K and a vortex Reynolds number Re of 700, following the second example in [START_REF] Thévenin | [END_REF]. This case describes an ignition process for a higher air temperature, characterized by a quick self-ignition of the vortex center resulting in a premixed flame that propagates outwards from the center. Actually, this is exactly the opposite configuration to the previous one in terms of flow and initial temperature parameters. A diffusion flame is independently developed in the braids together with a premixed flame on the fuel side of the domain, close to the diffusion flame layer. As studied in [START_REF] Thévenin | [END_REF], this configuration corresponds to the premixed ignition mode. Once again these computations qualitatively reproduce previous results but with higher values of heat release. This is a direct consequence of the finer spatial resolution used in the present computations that yields better discrete representations of the highly nonlinear function F. In the present illustration a shorter distance can be observed between both diffusion and premixed flames in the braids. Additional runs demonstrated that this distance is related to the initial solution (25) by means of the α parameter, which is very likely not the same in both computations 6 . Note that thermal runaway is almost instantaneous for high temperatures and takes place during the transient phase associated to the initial highly increasing velocity field, and consequently to sufficiently small splitting time steps. This simulation takes approximately 1 minute of CPU time for t ⋆ up to 10 -5 .

Finally, Figure 7 illustrates the third and last case corresponding to an intermediate air temperature T O,0 of 1200 K and vortex Reynolds number Re of 500. A mixed ignition regime is observed in this configuration in which the diffusion flame in the braids is ignited before the vortex core. In accordance to the intermediate flow and temperature conditions, the dynamics observed in the vortex center lies in between previous configurations; the mixture in the core eventually burns and is progressively consumed by both types of flames: from the outer edge by the diffusion flame originating from the braids as in the first case, and outwards by the selfignition of the well mixed core as in the second case. Results in [START_REF] Thévenin | [END_REF] are well retrieved but with an improved representation of the heat release distribution. Note that scales of heat release rate change in each snapshot for better visualization. Thermal runaway takes place after the initial vortex transient as in the first case, and the simulation takes approximately 4 minutes of CPU time for t ⋆ up to 10 -4 . Notice that all of these computations were performed in a very small amount of CPU time, thanks to the time-space adaptive capabilities of the method as we shall see in the following section, and hence exhaustive parametric studies could be easily conducted for different applications. Further interpretations of the various regimes of ignition and resulting configurations can be found in [START_REF] Thévenin | [END_REF].

In all three cases, calculations indicate that in regions where the flame is essentially of diffusion type, its thickness changes as a factor of the local strain rate: where the latter is low the thickness increases as observed in the braids at a distance from the vortex core. Conversely, the layer gets thinner closer to the vortex core as the local strain rates get considerably larger. This behavior takes place if the strain rate is below the critical extinction value and under these circumstances the effect of the chemical time on the thickness variation is relatively limited.

Characterization of numerical errors and performance of the method

In the previous section the ignition dynamics was briefly described from a purely physical point of view for different configurations in order to establish qualitative comparisons with previous results, and thus confirm the validity of the present results. In this part we are interested in evaluating the efficiency of the present method in a quantitative way and in particular throughout the fast transition phase, numerically the most difficult part of the problem, contrary to the previous situation in Figures 5,6, and 7 where the study was carried out practically after the thermal runaway process. Once again, as done in §4.2.2, we will define a sufficiently accurate reference solution in order to evaluate the numerical errors of the method, taking into account that no analytic solution exists for such vortex-flame configurations. A brief study is nevertheless presented in Appendix E to complement the following analysis for a much simpler problem with analytic solution: a planar strained flame in the limit of fast chemistry. We consider in this section fresh fuel initially at T F,0 equal to 300 K and hot air at a temperature T O,0 of 1000 K, corresponding to the third and last ignition mode for intermediate air temperatures, but with a higher vortex Reynolds number Re of 1000. By decreasing the air temperature we aim at shifting the sudden change in physics observed in Figure 7 and certainly not known at the start of the calculation, away from the initial transient behavior and corresponding small splitting time steps as seen in Figure 6. Beyond the inherent stiffness of the governing equations, a higher Reynolds number also imposes more severe convective conditions. In this way the chosen configuration qualitatively reproduces the physics encountered in the previous part but under numerically tougher constraints. Figure 8 (top) illustrates the thermal runaway during the ignition process. Temperatures are displayed at times t ⋆ of 5 × 10 -5 , 10 -4 , and 1.5 × 10 -4 , for which the maximum temperature is given, respectively, by approximately 1094 K, 2108 K, and 2209 K. The corresponding adapted grids are also shown in Figure 8 (bottom) for a more illustrative three-dimensional representation of the temperature. We have verified that the mixture transverse size and the corresponding flame thickness are of the order of 0.025, similar to the previous propagating case. A spatial discretization of 1024 2 points is thus reasonably accurate. Figure 8 (bottom) shows that the finest regions dynamically identify the local rise in temperature until ignition of the entire contact surface. In particular, the initial front does not require a full representation on the finest grid for a threshold tolerance of ε = 10 -3 (see Figure 8 (bottom left)). The complete refinement of the front then takes place after ignition and subsequent formation of steeper gradients. One finds that for smaller threshold values the initial front itself lies within the finest grid, but it is interesting to retain the current configuration to illustrate a limit case.

∆t ∆t R1 ∆t D1 ∆t C η=10 -3 -ε=10 -3
Figure 9: Two-dimensional flame ignition in a vortex field. Time adaptation is well apparent in the time evolution of splitting time steps ∆t with accuracy tolerance η = 10 -3 (left), as well as the reactive ∆t R1 and diffusive ∆t D1 substeps during the first splitting half-step, and the convective ones ∆t C in the Strang scheme (7) (right). Solutions were computed on a dynamically adapted grid corresponding to 1024 2 points at the finest level J = 10 with ε = 10 -3 .

From now on we focus on the time interval defined by t ⋆ and [0, 1.5 × 10 -4 ], which contains the physical transition from inert mixing to ignition of the reactants. Figure 9 illustrates the resulting time adaptation featured by the numerical strategy. As in the previous computations in §4, an initial splitting time step of ∆t = 10 -8 with the growth limiting procedure (19) were considered to properly handle the inclusion of the vortex and the fast variation of the velocity field. The splitting step increases until t ⋆ ≈ 6.5 × 10 -5 (∆t ≈ 1.89 × 10 -5 ) during the mixing phase, and one then finds a series of rejected steps for the given accuracy tolerance η of 10 -3 . The splitting time step is thus reduced down to the time scale needed to guarantee the prescribed accuracy: ∆t ≈ 1.12 × 10 -7 . This behavior naturally coincides with the sudden ignition of the flame and the subsequent fast propagation along the contact surface, once a certain temperature is locally reached after the initial mixing of reactants. A dynamic adaptation of the splitting time step is hence mandatory to identify these changes in the physical behavior of the phenomenon and to suitably describe this process. As in the previous case in §4, a constant and relatively large δ-shift in ( 8) is accurate enough, since strong variations in time stepping are associated to important time step reductions for which the theoretical framework of the adaptive splitting scheme is certainly valid [56]. Nevertheless, the dynamic evaluation of δ would be required for more general configurations.

The resulting reaction, diffusion, and convection time integration steps for η = 10 -3 are displayed in Figure 9 (right). Only the minimum reaction steps are represented which are of the order of ∆t R ≈ 10 -6 and ∆t R ≈ 5 × 10 -7 at the flame front, respectively, before and after complete ignition of the flame. For each splitting time step the local reaction time steps progressively increase from the depicted values up to ∆t R = ∆t/2 for the points lying away from the highly reactive area. Diffusion time steps are of the order of ∆t D ≈ 10 -6 , whereas the convective step ranges from ∆t C ≈ 10 -9 in the beginning to ∆t C ≈ 10 -7 . Once again, the convection time step is the most constraining step considering the high Reynolds number value adopted in the present calculations.

In what follows we investigate the accuracy and computational costs of the numerical approximations issued from the proposed adaptive method, as previously performed in §4. We thus consider again the quasi-exact and splitting solutions represented on a uniform grid of 1024 2 points, and computed with, respectively, a constant splitting time step ∆t of 10 -7 and a splitting time stepping according to the η accuracy tolerance. Let us underline that the quasi-exact approximation is indeed a splitting one but with a sufficiently small time step such that all time scales are practically coupled (see, e.g., Figure 9) and the splitting errors are negligible. Alternatively the adaptive method was implemented with different splitting accuracy tolerances η of 10 -3 , 10 -4 , 10 -5 , and 10 -6 . In each case the multiresolution threshold value ε was taken equal to the corresponding η parameter, i.e., η = ε, following previous results in §4 7 ; we shall further discuss this choice later on. Figure 10 (left) illustrates the splitting time stepping for different η tolerances.

No rejection of time steps is observed for η equal to 10 -5 or smaller, and several time steps in the transient phase coincide for different tolerances due to the limiting growth procedure (19). Time integration on a uniform grid yields essentially the same splitting time stepping observed in Figure 10. Additionally Figure 10 (right) illustrates the time evolution of data compressions DC.

Recalling that time integration is performed on a constant adapted grid and that the remeshing period is set by the splitting time step, the dynamic evaluation of the splitting time step allows an adequate updating of the spatial representation and consequently the necessary refinement of the spatial configuration corresponding to the new physics of ignition. The latter is particularly observed for the case ε = 10 -3 , a limit case for which the initial front is not fully represented on the finest grid, as previously remarked in Figure 8. Multiresolution, splitting, and combined time-space adaptive approximation errors deduced from (21) are gathered in Table 3 8 , computed for variable θ at two different times t ⋆ : 5 × 10 -5 and 10 -4 , approximately before and after the main thermal runaway mechanism. It can be seen that multiresolution errors E J MR reproduce the dependence on the threshold value ε into (10), where a relatively larger constant C of O(10) can be noticed for ε = 10 -4 and 10 -5 . Additionally, splitting errors E J split are effectively controlled by the local error accuracy η, except for t ⋆ = 5 × 10 -5 where an apparently constant error is maintained, very likely associated to the numerical initialization of the problem. Finally the global error of the method E MR split is shown to be composed of both multiresolution and splitting errors as established in (11) 9 . In general we expect that same values for both tolerances η and ε involve errors of the same order and an effective control of the global error, sum of all numerical approximations. Nevertheless, since local time integration errors, controlled by η, accumulate in time and on the other hand the multiresolution transform for a given ε is performed at each time iteration on the current solution, we could also expect that the former errors will eventually pilot the global one in [START_REF] Staffelbach | High Performance Computing for Computational Science -VECPAR 2008[END_REF], as seen previously in §4 for long time integration domains and in Table 3 for η = ε = 10 -6 . To conclude, one key point is that the compressed spatial representations must be accurate enough to guarantee a reliable numerical solution of the time dependent problem, taking into account that the time integration is performed on a fixed adapted grid during each time step. One simple way of enforcing this behavior considers threshold values ε smaller than η, and at most equal-valued tolerances, as illustrated in this study. Table 4 summarizes the CPU times in minutes for all previous numerical solutions considered in Table 3. Important gains in CPU time are achieved with the adaptive splitting technique, which are greatly improved by the time-space adaptive strategy. For instance, for a set of parameters η = 10 -3 and ε = 10 -3 , splitting adaptation implies a gain factor of about 3.25 with respect to the quasi-exact solution (from U J qe to U J split in Table 4), further increased to about 75 with both time and space adaptation (from U J qe to U MR split ). This global gain comes indeed from both adaptive procedures since a multiresolution solution with ε = 10 -3 and constant splitting time step ∆t of 10 -7 implies a factor of "only" 12 (this computation requires about 56.80 minutes of CPU (Notice that a quasiexact solution with a time stepping defined by the strongest numerical restriction, in this case set by the stability of the convective scheme, yields even more expensive computations as can be inferred from Figure 9 (right).) Time reductions related to the multiresolution representation (from U J split to U MR split ) are consistent with those achieved by data compressions in Figure 10 (right), for example, about 95 %, 91 %, and 84 % for a threshold value ε of, respectively, 10 -3 , 10 -4 , and 10 -5 . Conversely, gains related to splitting adaptation (from U J qe to U J split ) seem to be roughly the same in these computations for η equal to 10 -3 , 10 -4 , or 10 -5 , which underlines the impact of the inner solvers like Radau5 or ROCK4 on the global performance of the time integration method. In this particular study all computing parameters of these solvers were set and maintained to their standard default values which in some cases unnecessarily increase the number of function evaluations, for instance, for Jacobian or spectral radius computations in Radau5 or ROCK4.

Finally, it is important to point out that the global efficiency of the time-space method is certainly a problem dependent feature, that is nevertheless very likely to be high for strongly dynamic mechanisms and localized spatial structures, as illustrated in the present problem.

Conclusion and perspective

Numerical results obtained with the present time-space adaptive technique support the conclusions that different multi-scale physical configurations can be successfully simulated and that the error can be effectively controlled. Important gains in computational efficiency are achieved because of highly compressed data representations, as well as a dynamic splitting technique with adequate solvers, independent time stepping procedures, and splitting time steps not restricted by stability constraints. It is shown that the adaptive splitting scheme is critical in handling fast transients, and that it allows to deal with the difficult problem of the sudden ignition of a flame in a vortex field, which features a variety of spatial and time scales and different regimes of flame initiation. This kind of multi-scale problem would become extremely expensive, if it were envisaged without these adaptive capabilities, at least with standard computational resources.

Concerning error control features, the present numerical strategy is perfectly inscribed in a new solution scheme paradigm that aims at enhancing high fidelity numerical simulations with tools providing estimates "on the fly" of the quality of numerical results for general multi-dimensional configurations. The present scheme is such that for a given spatial discretization, the numerical accuracy of the simulations is set by two parameters:

• The threshold value ε of the multiresolution decomposition, which balances data compression and numerical errors related to compressed data representations;

• The accuracy tolerance η of the time splitting technique, which limits the degree of decoupling of the physical phenomena and hence, controls the numerical time integration errors.

This provides a solid basis for more detailed and complex numerical simulations. This is exemplified, for instance, with applications that include detailed chemical features, presented in [START_REF] Duarte | [END_REF], or more complex models such as gas discharge problems [86]. In the particular context of low Mach combustion problems, a numerical solver can be constructed by considering multiresolution decomposition and the present adaptive splitting technique embedded, for instance, in a classical pressure-projection technique. Notice than in the latter approach, species and temperature equations are evolved in time with an imposed flow during each time step [36,41] as in the illustrations considered in this work. Although the extension to systems involving large number of variables like in detailed chemical kinetics models is straightforward [START_REF] Duarte | [END_REF], the integration of the source term leads to important computational costs. This is so despite of the fact that the latter is embarrassingly parallel in the framework of operator splitting 10 and that computational requirements can be reduced by data compression resulting from the multiresolution analysis. This issue requires further studies in order to obtain high efficiency in terms of load balancing on parallel architectures and to better exploit current computational resources.

where the convection steps ∆t C are limited by the stability restrictions of the numerical scheme. This splitting procedure allows integrations with the theoretically right accuracy and easily ensures multi-dimensional MP/TVD constraints [75]. The stability constraint is given by a standard CFL condition (inversely proportional to the finest spatial grid among all dimensions), and thus a better solution considers rather

C 2∆t C U(t) = C ∆t C x C ∆t C y C ∆t C z C ∆t C z C ∆t C y C ∆t C x U(t), (A.2)
instead of (A.1), to better ensure the same numerical diffusion in all three directions and to preserve the isotropy of the computations. Furthermore, at each time step ∆t C we need to perform three steps in (A.2) to advance the solution, instead of five in (A.1). In the splitting scheme (7), the operator C is thus given by

C ∆t = I C i=1 C 2∆t C,i , (A.3)
such that 2I C convection steps ∆t C,i are performed within the global splitting step ∆t, and formula (A.2) is recast as

C 2∆t C,i = C ∆t C,i x C ∆t C,i y C ∆t C,i z C ∆t C,i z C ∆t C,i y C ∆t C,i x . (A.4)
The intermediate time step ∆t C,i is the same for all points over the computational domain, and it is given by

∆t C,i = min         ∆t max C,i , ∆t 2 - i-1 i ′ =1 ∆t C,i ′         , (A.5)
where ∆t max C,i is the current maximum convection time step within the stability domain of the numerical scheme. The previous procedure is general and remains valid for any convective scheme and for both linear and nonlinear transport problems, with time-or space-varying transport velocities.

verifies a time dependent equation of type (2) without source term like (3), whereas from a simple energy balance relation in an adiabatic framework, we get

c p (T b -T o ) = Q(Y Fo -Y Fb ). (B.4)
By evaluating θ Zo and θ Zb in (B.3) and from (B.4), it can be seen that θ Z is constant throughout the flame. Consequently, a progress variable c(x, y, t) can be introduced: taking into account that Y Fb = 0 into (12), with

c = T -T o T b -T o = Y Fo -Y F Y Fo -Y Fb . ( B
B ⋆ = B 1 W O W F ρ 2 o Y Oo T 2 o . (B.7)
Hence for the fuel mass fraction equation in (2), one now has

∂ t Y F + v x ∂ x Y F + v y ∂ y Y F -D ∂ 2 x Y F + ∂ 2 y Y F = - B ⋆ ρ o Y Fo (1 -c)e -T a /(T o (1+τc)) , (B.8)
which may be written as

∂ t c + v x ∂ x c + v y ∂ y c -D ∂ 2 x c + ∂ 2 y c = B ⋆ ρ o (1 -c)e -T a /(T o (1+τc)) . (B.9)
Considering a square computational domain of size 2L, a characteristic diffusion time τ d = L 2 /D, and a velocity V = D/L, one may define the following dimensionless variables:

x ⋆ = x L , y ⋆ = y L , v x,⋆ = v x V , v y,⋆ = v y V , t ⋆ = t τ d . (B.10)
We finally obtain [80]: The three-dimensional velocity field is defined by a toroidal vortex directly deduced from the previous two-dimensional velocity field. In each plane containing the z-axis, we consider thus a pair of counter-rotating vortices computed as usual with (14) where the radius is defined by r 2 ⋆ = (x ⋆x 0,⋆ ) 2 + (y ⋆y 0,⋆ ) 2 + (z ⋆z 0,⋆ ) 2 , (C.1) centered at x 2 0,⋆ + y 2 0,⋆ 1/2 = 0.25, and with z 0,⋆ = -0.5. Although the resulting field is not divergence free, it suffices to construct a three-dimensional configuration to illustrate the numerical capabilities of the method. The same simulations can be performed exactly in the same way with more physically consistent velocity fields.

∂ t ⋆ c + v x,⋆ ∂ x ⋆ c + v y,⋆ ∂ y ⋆ c -∂ 2 x ⋆ c + ∂ 2 y ⋆ c = Da (1 -c)e -T a /(T o (1+τc)) , (
Figure 4 shows the interaction of the initial premixed flame with the toroidal vortex, and the corresponding adapted grids on which the solutions are computed. Additionally, Figure C.11 shows the corresponding time steps of integration. We retrieve a qualitatively similar behavior with respect to the previous two-dimensional case, in terms of splitting time steps and the time stepping for each split subproblem. The fresh gases are completely burnt at a time t ⋆ = 3 × 10 -3 . Considering the obtained data compression, no more than 18 % of the 256 3 points are required. This simulation took approximatively 17.26 hours of CPU time. where, for a given dimensionless strain rate β(t ⋆ ), the velocity field is defined by v x,⋆ = βx ⋆ , v y,⋆ = -βy ⋆ .

(E.2)

The assumption of an infinitely fast chemistry implies also that Y F (t ⋆ )Y O (t ⋆ ) = 0 throughout the domain.

The symmetry of the resulting physical configuration with respect to x ⋆ = 0 yields a onedimensional formulation of the problem at the latter position, noticing that all time dependent variables in (E.1) will depend only on (y ⋆ , t ⋆ ). Using the following definitions and the analytic solution for (E.1) may be cast in the form (see, e.g., [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]):

Y F = Y F,0 1 -erf(κ f ) erf(κ) -erf(κ f ) , κ > κ f , (E.5)
where κ f corresponds to the actual position of the flame in the frame of reference of similarity variable κ, and it is related to the global equivalence ratio φ, defined in (D.2), through

erf(κ f ) = 1 -φ 1 + φ . (E.6)
For a given space representation with mesh size ∆y ⋆ corresponding to a grid level J, the steady solution (Y F ) J given by (E.5) will be then compared to the numerical approximation obtained with the multiresolution/splitting solution (Y F ) MR split , applied to problem (E.1) with (E.2). Contrary to the general framework of the present paper, an analytic solution is available for this particular problem, and the numerical errors characterized in [START_REF] Staffelbach | High Performance Computing for Computational Science -VECPAR 2008[END_REF] .

(E.7) Notice that while (Y F ) J is simply evaluated in (E.5) at grid points of a given uniform mesh J, the quasi-exact solution defined in §3: (Y F ) J qe in (E.7), represents a numerical approximation to (Y F ) J , resulting from a fine numerical integration of the coupled diffusion-convection problem (E.1) on the same uniform grid.

To simplify the expressions defined in (E.3), we consider the case with a constant strain rate β; we set L = 1 for the computational domain and the following physical parameters: Y F,0 = 1, Y O,0 = 0.23, and s = 8, as in §5, for a global equivalence ratio φ of 34.78 and κ f ≈ -1.35 according to (E.6). Then, for a given mesh J we initialize the problem by evaluating (E.5) with β = β 0 = 10 2 and t ⋆ = t⋆ = β -1 0 in (E.3). After the first time integration step ∆t, we switch β to β 1 = 10 4 and continue to numerically integrate (E.1) for t ⋆ into [0, 8 × β - 1 1 ]. The initial and final solutions are shown in Figure E.12 (left) for both analytic and numerical solutions with a multiresolution parameter of ε = 10 -3 and grid level J = 10 (equivalent to a space discretization of 1024 points) for the latter one. Figure E.12: Planar strained diffusion flame at x ⋆ = 0. Left: initial and final profiles of Y F , computed with analytic expression (E.5) (resp., solid and dashed lines) and multiresolution approximation (resp., bullets and circles) with ε = 10 -3 and a finest discretization ∆y ⋆ corresponding to 1024 points in y ⋆ (J = 10). Right: total numerical error (E.7) computed for different space discretizations ∆y ⋆ and multiresolution accuracies ε.

To illustrate the convergence of the numerical approximation (Y F ) MR split towards the analytic solution (Y F ) J with strain rate β 1 , we consider a small and constant splitting time step, ∆t = 10 -8 , such that the time integration (splitting) error becomes negligible in (E.7), i.e., (Y F ) J qe ≈ (Y F ) J split . In this way according to (E.7) the numerical error is mainly composed of the space discretization error, which should behave like ∆y 2 ⋆ , and the multiresolution one, proportional to ε following the bound in (10). Figure E.12 (right) shows the convergence of the numerical solution with ∆y ⋆ , measured by the total numerical error in (E.7) at x ⋆ = 0, for multiresolution approximations with different level of accuracy: from ε equal to 10 -2 to 10 -7 . The error behaves approximately like ∆y 2.77 ⋆ for any ε and relatively coarse grids (J < 9 or ∆y ⋆ > 3.9 × 10 -3 ); and like ∆y 1.48 ⋆ for finer grids and sufficiently small ε such that the total error is mainly due to the space discretization in (E.7) 11 . Alternatively, for sufficiently fine representations, the multiresolution error becomes dominant and hence the accuracy of the numerical approximation is mostly related to ε, regardless of the space discretization. An optimum choice of ∆y ⋆ and ε is then retrieved if both multiresolution and space discretization errors are of the same order. However, since analytic solutions are not available in general applications, it would be preferable to ensure that the multiresolution error will dominate the space discretization one in (E.7), and thus the total numerical error can be effectively controlled through the choice of ε. This is the approach considered throughout this study.

Figure 1 :

 1 Figure1: Two-dimensional propagating flame. Time evolution of progress variable c at t ⋆ = 5×10 -4 (left), 10 -3 (center), and 1.5 × 10 -3 (right). Top: red (resp., blue) zone corresponds to burnt (resp., fresh) gases, c = 1 (resp., c = 0). Middle: contour lines with c = 0 -0.99 and ∆c = 0.11. Bottom: dynamic adapted grid corresponding to 1024 2 points at the finest level J = 10 with ε = 10 -3 .

3 Figure 2 :

 32 Figure 2: Two-dimensional propagating flame. Left: velocity v y,⋆ at y ⋆ = -0.5, at t ⋆ = 10 -3 and t ⋆ = 3 × 10 -3 . Right: time evolution of data compressions DC for ε = 10 -2 and ε = 10 -3 .

3 Figure 3 :

 33 Figure 3: Two-dimensional propagating flame. Time evolution of the splitting time step ∆t, the reactive ∆t R1 and diffusive ∆t D1 substeps during the first splitting half-step, and the convective ones ∆t C according to the Strang scheme (7) with tolerances η = 10 -3 and ε = 10 -2 (left) or ε = 10 -3 (right).

Figure 4 :

 4 Figure 4: Three-dimensional propagating flame. Time evolution of progress variable c at t ⋆ = 5 × 10 -4 (left) and 1.5 × 10 -3 (right). Top: isosurfaces for c equal to 0.01 (blue) and 0.99 (red). Bottom: dynamic adapted grids (right) corresponding to 256 3 points at the finest level J = 8. Contour lines with c = 0.01 -0.91 and ∆c = 0.3.

Figure 5 :

 5 Figure 5: Two-dimensional flame ignition in a vortex field. Instantaneous heat release rate F at t ⋆ = 4 × 10 -4 (top left), 8 × 10 -4 (top right), and 10 -3 (bottom left). Diffusion ignition mode for a moderate air temperature T O,0 of 800 K and vortex Reynolds number Re of 160. Bottom right: time evolution of splitting time steps and maximum temperature T , deduced from θ. Rejected time steps are indicated with black bullets (•) while maximum temperatures for previous snapshots are marked with (•).

Figure 6 :

 6 Figure 6: Two-dimensional flame ignition in a vortex field. Instantaneous heat release rate F at t ⋆ = 4 × 10 -6 (top left), 7 × 10 -6 (top right), and 9 × 10 -6 (bottom left). Premixed ignition mode for a high air temperature T O,0 of 2000 K and vortex Reynolds number Re of 700. Bottom right: time evolution of splitting time steps and maximum temperature T , deduced from θ. Maximum temperatures for previous snapshots are marked with circles (•).

Figure 7 :

 7 Figure 7: Two-dimensional flame ignition in a vortex field. Instantaneous heat release rate F at t ⋆ = 4 × 10 -5 (top left), 7 × 10 -5 (top right), and 8.5 × 10 -5 (bottom left). Mixed ignition mode for an intermediate air temperature T O,0 of 1200 K and vortex Reynolds number Re of 500. Bottom right: time evolution of splitting time steps and maximum temperature T , deduced from θ. Rejected time steps are indicated with black bullets (•) while maximum temperatures for previous snapshots are marked with (•).

Figure 8 :

 8 Figure 8: Two-dimensional flame ignition in a vortex field. Time evolution of temperature T at t ⋆ = 5 × 10 -5 (right), 10 -4 (center), and 1.5 × 10 -4 (left). Initial temperature of the fresh fuel: T F,0 = 300 K, and of the hot air: T O,0 = 1000 K. Bottom: dynamic adapted grid corresponding to 1024 2 points at the finest level J = 10 with ε = 10 -3 .

5 Figure 10 :

 510 Figure 10: Two-dimensional flame ignition in a vortex field. Time and space adaptation given by the time evolution of splitting time steps ∆t (left) and of data compressions DC (right) for several accuracy tolerances η and ε.

. 5 )

 5 Defining τ = T b /T o -1, we obtain that T/T o = 1 + τc and the reaction rate (B.2) becomes ẇ = B ⋆ Y Fo (1c)e -T a /(T o (1+τc)) , (B.6)

  B.11) where Da = B ⋆ τ d /ρ o = τ d /τ ch is a Damköhler number and τ ch = B ⋆ /ρ o , a chemical time.

2 Figure C. 11 :

 211 Figure C.11: Three-dimensional propagating flame. Left: time evolution of the splitting time step ∆t, the reactive ∆t R1 and diffusive ∆t D1 substeps during the first splitting half-step, and the convective ones ∆t C with tolerances η = 10 -3 and ε = 10 -2 . Right: time evolution of data compressions DC with ε = 10 -2 .

  can be recast as(Y F ) J -(Y F ) MR split L 2 ≤ (Y F ) J -(Y F ) J qe L 2 space discretization error + (Y F ) J qe -(Y F ) J split L 2 splitting error + (Y F ) J split -(Y F ) MR split L 2 multiresolution error

  1. Mass diffusion velocities of chemical species are expressed by Fick's law. 2. Thermal diffusion of species (Soret-Dufour effect) is neglected. 3. Different species have constant and equal diffusion coefficients with respect to the mixture, noted D. 4. Constant pressure specific heats of all species are constant with the same value c p . 5. Lewis numbers corresponding to all species are equal to one. 6. The rate of pressure change in time is negligible. 7. Density variations associated to chemical heat release are neglected.

Table 1 :

 1 L 2 numerical errors for the time adaptive splitting (E J split ), space adaptive multiresolution (E J MR ), and time-space adaptive (E MR split ) strategies evaluated at different times. Finest grid: 1024 2 .

Table 2 :

 2 CPU time in minutes for the time-space adaptive, the time adaptive splitting, and the quasi-exact strategies for t ⋆ into [0, 2 × 10 -3 ]. Finest grid: 1024 2 .

	MR-splitting ε = splitting quasi-exact 10 -2 10 -3
	CPU time (min.) 56.27	71.05	589.00	6603.26

Table 3 :

 3 L 2 numerical errors for the space adaptive multiresolution (E J 10 -3 1.87 × 10 -3 3.26 × 10 -3 10 -4 10 -4 3.32 × 10 -3 3.26 × 10 -3 10 -5 10 -5 2.47 × 10 -4 1.82 × 10 -4 10 -6 10 -6 1.12 × 10 -4 2.41 × 10 -5

	MR ), time adaptive splitting (E J split ), and time-space

Table 4 :

 4 CPU time in minutes for the time-space adaptive, the time adaptive splitting, and the quasi-exact strategies for t ⋆ into [0, 1.5 × 10 -4 ] and several accuracy tolerances η and ε. Finest grid: 1024 2 . contrary to the factor of 23 observed in Table4from U J split to U MR split .

		η	ε	CPU time (min.)
		10 -3 10 -3	8.93
	U MR split	10 -4 10 -4 10 -5 10 -5	16.52 34.68
		10 -6 10 -6	131.91
	U J split U J qe	10 -3 10 -4 10 -5 10 -6 -	-----	207.52 196.47 225.95 480.22 674.69
	time),

To enhance spatial resolution for reacting flows problems, compact schemes of high order were considered in[14,15,16], as well as locally refined grids for steady-state solutions of low Mach problems[13].

One may note that for stiff PDEs in splitting configurations, an important loss in efficiency is expected with stiff multi-step solvers like VODE or LSODE[71], mainly because of the expensive and less accurate starting procedure of multi-step schemes at each splitting time step, as demonstrated, for instance, in[72].

In this work, as in most multiresolution applications, we consider Cartesian meshes. Nevertheless, further developments in data structure conception as the ones introduced in[76] extend these ideas to more general geometries, as described, for instance, in[START_REF] Ballmann | Summary of Flow Modulation and Fluid-Structure Interaction Findings[END_REF] for compressible flows and fluid-structure interactions problems.

It was shown in[56] that large values of δ (noted ε in[56]) will in general extend the valid working region of the adaptive splitting technique, i.e., local error estimates err should remain valid for larger splitting time steps.

Smaller values of α imply larger distances between flames without significant alteration of the dynamics of the vortex center. Computations in[START_REF] Thévenin | [END_REF] seem to consider α equal to 100 for this case.

The configuration given by η = 10 -3 and ε = 10 -2 was briefly investigated in[START_REF] Duarte | Ann. Res. Briefs[END_REF].

The values corresponding to the case η = ε = 10 -3 are slightly different in[START_REF] Duarte | Ann. Res. Briefs[END_REF], because in this study all L 2 -norm errors are normalized by (max θmin θ), as it is done in the numerical code, taking into account that θ ≤ 1 and can be negative according to(22).

Although this expression is demonstrated in practice, further studies are required to gain an insight into the interaction of these space and time approximation errors.

A numerical strategy is considered embarrassingly parallel when no effort is required to separate the problem into a number of parallel tasks. In the present study, due to operator splitting, the solution of the reaction term at each spatial node represents an independent system of ODEs with size equal to the number of physical variables of the problem.

Even smaller splitting time steps are needed for finer discretizations in order to assure negligible time integration errors in (E.7) and a behavior closer to ∆y 2 ⋆ .

Acknowledgements

This research was supported by a fundamental project grant from ANR (French National Research Agency -ANR Blancs): Séchelles (PI S. Descombes -2009Descombes - -2013) ) and by a DIGITEO RTRA project: MUSE (PI M. Massot -2010-2014). The support of the France-Stanford Center for Interdisciplinary Studies through a collaborative project grant entitled "Multi-scale mathematical modeling and numerical methods for multiphase and reactive flows" (PIs: P. Moin and M. Massot) has been very helpful. We also acknowledge the computational resources of the Mesocentre of Ecole Centrale Paris where some of the simulations were performed.

Appendix A. Numerical implementation of the convective scheme A Strang dimensional splitting [27] is implemented for the convective scheme to handle multi-dimensional configurations:

x U(t), (A.1)

Appendix B. Derivation of the modeling equations for the propagation of premixed flames

The reaction rate into (2) is modeled by the following Arrhenius equation [80]:

where B 1 is a pre-exponential factor and T a , the activation energy. For premixed laminar flames, the mixture may be assumed to be fuel lean with a high diluent concentration. Hence the reaction rate is controlled by the fuel mass concentration, whereas the oxidizer mass fraction is nearly constant and equal to its upstream value,

Let us recall that subscripts ( ) o , ( ) b , and ( ) ⋆ respectively indicate, fresh mixture zone, burnt product zone, and dimensionless variables. The variation of the ratios ρ o /ρ b is assumed negligible because the pressure is essentially constant according to hypothesis 6. With these hypotheses, the composite Schvab-Zeldo'vich variable:

Appendix C. Three-dimensional interaction of a premixed flame with a toroidal vortex

We consider the solution of ( 13) over a dimensionless computational domain of [-1, 1] 3 . The same data initialization as in the two-dimensional case is considered with, respectively, z ⋆ and z 0,⋆ , instead of y ⋆ and y 0,⋆ in (18), as well as the same previous modeling parameters. The adaptive splitting accuracy tolerance in ( 9) is also set to η = 10 -3 , with η ROCK4 = 10 -5 , η Radau5 = 10 -7 , and ε = 10 -2 for the multiresolution threshold value, following the previous twodimensional results. The time domain of integration is given by t ⋆ into [0, 3.5 × 10 -3 ], whereas the finest grid corresponds to a spatial discretization of 256 3 points, i.e., J = 8 as finest grid level.

Appendix D. Details on the modeling equations for the ignition of diffusion flames

The dynamics of the models studied in §5 are described by the general set of equations ( 2) with (3), where the reaction rate is modeled by

in which A is a pre-exponential factor and T a , the activation energy.

One may construct Schvab-Zeldo'vich variables by combining the reacting species mass fractions with proper coefficients, to obtain a balance equation without source term, analogous to the equation (3) governing Y N . Thus, introducing the reduced total heat released χ; the normalized 29 temperature difference between reactants τ; the product to fuel stoichiometric ratio σ; the global equivalence ratio corresponding to a complete mixing between reactants in their initial state φ; and the stoichiometric factor s, defined respectively by

one may define the following variables [START_REF] Thévenin | [END_REF]:

recalling that θ is the reduced temperature defined as

This set of variables (Z 1 , Z 2 , Z 3 ) are initially equal and follow the same balance equation without reaction term (like (3)) and with the same boundary conditions. They are therefore equal at each point and for all times to the same value Z. Introducing the same dimensionless variables previously defined in (B.10), we can derive the following system of equations [START_REF] Thévenin | [END_REF]:

where τ a = T a /T O,0 is the reduced activation temperature and with the Damköhler number defined by Da = ρAτ d .

Appendix E. Planar strained diffusion flame in the fast chemistry limit

In order to present a simple and complementary validation of the numerical strategy established in §3, we briefly investigate in this part a planar strained diffusion flame with an infinitely fast chemistry such that an analytic solution to the problem can be derived. The theoretical study of this kind of configurations can be found, for instance, in the book of Poinsot & Veynante [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF].

Let us consider a square two-dimensional domain of size 2L in which fuel with mass fraction Y F,0 at temperature T F,0 is injected from the top of the domain at y = L, while air at T O,0 with an oxidizer mass fraction of Y O,0 is also introduced, this time from the opposite side at y = -L. The latter counterflow configuration results in a steady planar diffusion flame placed at some intermediate ordinate y between L and -L. Using the modeling assumptions established in §2 in combination with the fast chemistry limit, and using the dimensionless variables previously defined in (B.10), leads to the following time dependent equations for k = F, O: