New Resolution Strategies for Multi-scale Reaction Waves: Optimal Time Operator Splitting and Space Adaptive Multiresolution (cid:63)

. We tackle the numerical simulation of reaction-diﬀusion equations modeling multi-scale reaction waves. This type of problems induces peculiar diﬃculties and potentially large stiﬀness which stem from the broad spectrum of temporal scales in the nonlinear chemical source term as well as from the presence of large spatial gradients in the reactive fronts, spatially very localized. In this paper, we introduce a new resolution strategy based on time operator splitting and space adaptive multiresolution in the context of very localized and stiﬀ reaction fronts. Thus, we aim at solving accurately complete models including all time and space scales of the phenomenon, considering large simulation domains with conventional computing resources. The eﬃciency is evaluated through the numerical simulation of conﬁgurations which were so far out of reach of standard methods in the ﬁeld of nonlinear chemical dynamics for 2D spiral waves and 3D scroll waves as an illustration


Introduction
Numerical simulations of multi-scale phenomena are commonly used for modeling purposes in many applications such as combustion, chemical vapor deposition, or air pollution modeling.In general, all these models raise several difficulties created by the large number of unknowns and the wide range of temporal scales due to large and detailed chemical kinetic mechanisms, as well as steep spatial gradients associated with very localized fronts of high chemical activity.Furthermore, a natural stumbling block to perform 3D simulations is the unreasonable memory requirements of standard numerical strategies for time dependent problems.
The most natural idea to overcome these difficulties is to use dedicated numerical methods and to solve the complete models where diffusion, reaction and eventually convection are coupled together.One aims at solving strongly coupled nonlinear systems with either a fully implicit method, or yet semi-implicit or linearized implicit methods instead (see [4] and references therein).However, the strong stability restrictions for the latter when dealing with very fast temporal scales, as well as the computing cost and the huge memory requirements of these methods, even if adaptive grids are used, make of these strategies difficult to be handled.
An alternative numerical strategy is then to combine implicit and explicit schemes to discretize nonlinear evolution problems in time.Further studies settled the appropriate numerical background for these methods called IMEX, which in particular might be conceived to solve stiff nonlinear problems [20,17].These methods are usually very efficient.Nevertheless, on the one hand, the feasibility of utilizing dedicated implicit solvers over a discretized domain becomes soon critical when treating large computational domains.And on the other hand, the time steps globally imposed over partial regions or the entire domain are strongly limited by either the stability restrictions of the explicit solver or by the fastest scales treated by the implicit scheme.
However, in many multi-scale problems, the fastest time scales do not play a leading role in the global physical phenomenon and thus, we might consider the possibility of using reduced models where these chemical scales have been previously relaxed.These simplified models provide reasonable predictions and the associated computing costs are significantly reduced in comparison with comprehensive chemical models.Nevertheless, these reduced models are usually accessible when the system is well-partitioned and the fast scales can be identified or isolated, a process that in realistic configurations, relies on sensitivity analysis which is most of the time difficult to conduct and justify.
It is then natural to envision a compromise, since the resolution of the fully coupled problem is most of the time out of reach and the appropriate definition of reduced models is normally difficult to establish.In this context, time operator splitting methods have been used for a long time and there exists a large literature showing the efficiency of such methods for evolution problems.In practice, it is firstly necessary to decouple numerically the reaction part from the rest of the physical phenomena like convection, diffusion or both, for which there also exist dedicated numerical methods.Hence, these techniques allow a completely independent optimization of the resolution of each subsystem which normally yields lower requirements of computing resources.
In the context of multi-scale waves, the dedicated methods chosen for each subsystem are then responsible for dealing with the fast scales associated with each one of them, in a separate manner; then, the composition of the global solution based on the splitting scheme should guarantee the good description of the global physical coupling.A rigorous numerical analysis is therefore required to better establish the conditions for which the latter fundamental constraint is verified.As a matter of fact, several works [21,18,4] proved that the standard numerical analysis of splitting schemes fails in presence of scales much faster than the splitting time step and motivated more rigorous studies for these stiff configurations [9,7] and in the case where spatial multi-scale phenomena arise as a consequence of steep spatial gradients [6].
We therefore introduce a new operator splitting strategy, based on these theoretical results, that considers on the one hand, a high order method like Radau5 [13], based on implicit Runge-Kutta schemes for stiff ODEs, to solve the reaction term; and on the other hand, another high order method like ROCK4 [1], based on explicit stabilized Runge-Kutta schemes, to solve the diffusion problem.Finally, the proposed numerical strategy is complemented by a mesh refinement technique based on Harten's pioneering work on adaptive multiresolution methods [14], being aware of the interest of adaptive mesh techniques for propagating waves exhibiting spatial multi-scale phenomena due to locally steep spatial gradients.The main goal is then to perform feasible and accurate in time and space simulations of the complete dynamics of multi-scale phenomena under study.
The paper is organized as follows: In Section 2, we first recall the time operator splitting schemes; then, in Section 3, we describe the construction of an optimal operator splitting scheme for multi-scale problems, and its coupling with a suitable grid adaptation strategy, the space adaptive multiresolution technique [2,15].2D and 3D simulations of a three species reaction-diffusion system modeling the Belousov-Zhabotinsky reaction are presented in Section 4 to illustrate the potential and performance of the method.We end in section 5 with some concluding remarks and prospects.

Time Operator Splitting
Let us first set the general mathematical framework in this work.A class of multi-scale phenomena are modeled by general reaction-diffusion systems like: where f : R m → R m and u : R × R d → R m , with the diffusion matrix D(u), which is a tensor of order d × d × m.
In order to simplify the presentation, we consider problem (1) with linear diagonal diffusion, in which case the elements of the diffusion matrix are written as D ijk (u) = D k δ ij , so that the diffusion operator reduces to the heat operator with scalar diffusion coefficient D k for component u k of u, k = 1, . . ., m.In any case, the proposed numerical strategy normally deals with general problem (1).Performing a fine spatial discretization, we obtain the semi-discretized initial value problem: where B corresponds to the discretization of the Laplacian operator with the coefficients D k within; U and F (U) are arranged component-wise all over the discretized spatial domain.Considering a standard decoupling of the diffusion and reaction parts of (2), we denote X ∆t (U 0 ) as the numerical solution of the discretized diffusion equation: with initial data U D (0) = U 0 after an integration time step ∆t.We also denote by Y ∆t (U 0 ) the numerical solution of the reaction part: with initial data U R (0) = U 0 .The two Lie approximation formulae of the solution of system (2) are then defined by whereas the two Strang approximation formulae [19] are given by where ∆t is now the splitting time step.
It is well known that Lie formulae (5) (resp.Strang formulae ( 6)) are approximations of order 1 (resp.2) of the exact solution of (2) in the case where X ∆t and Y ∆t are the exact solutions X ∆t and Y ∆t of problems (3) and (4).Then, appropriate numerical approximations of X ∆t and Y ∆t are required in order to compute Lie and Strang formulae with the prescribed order.
Higher order splitting configurations are also possible.Nevertheless, the order conditions for such composition methods state that either negative time substeps or complex coefficients are necessary (see [13]).The formers imply normally important stability restrictions and more sophisticated numerical implementations.In the particular case of negative time steps, they are completely undesirable for PDEs that are ill-posed for negative time progression.

Construction of the Numerical Strategy
In this section, we introduce a new splitting strategy for multi-scale waves modeled by stiff reaction-diffusion systems.Once established the time integration strategy, we detail briefly the adaptive multiresolution method that we have implemented as mesh refinement technique for this new resolution technique.

Time Integration Strategy
The standard orders achieved with a Lie or Strang scheme are no longer valid when we consider very stiff reactive or diffusive terms (see [9]).Furthermore, if the fastest time scales play a leading role in the global physics of the phenomenon, then the solution obtained by means of a splitting composition scheme will surely fail to capture the global dynamics of the phenomenon, unless we consider splitting time steps of the same order of such scales.
In the opposite case where these fast scales are not directly related to the physical evolution of the phenomenon, larger splitting time steps might be considered, but order reductions may then appear due to short-life transients associated to the fast variables.This is usually the case for propagating reaction waves where for instance, the speed of propagation is much slower than the chemical scales.In this context, it has been proved in [9] that better performances are expected while ending the splitting scheme by the time integration of the reaction part (4) or in a more general case, the part involving the fastest time scales of the phenomenon (see a numerical case in [7]).In particular, in the case of stiff reaction-diffusion systems with linear diagonal diffusion, no order loss is expected for the L ∆t 1 and S ∆t 2 schemes when faster scales are present in the reactive term.However, one must also take into account possible order reductions coming this time from space multi-scale phenomena due to steep spatial gradients whenever large splitting time steps are considered, as analyzed in [6].
All these theoretical considerations give us some insight into the numerical behavior of splitting techniques and thus, help us to select among the various splitting alternatives, depending on the nature of the problem.Nevertheless, the choice of suitable time integration methods for each subsystem is mandatory not only to guarantee such theoretical descriptions but also to take advantage of the particular features of each independent subproblem in order to solve them as accurate as possible with reasonable resources, as it is detailed in the following.Time Integration of the Reaction: Radau5.Radau5 [13] is not only an A-stable method, but also L-stable, so that very stiff systems of ODEs might be solved without any stability problem.It considers also an adapting time step strategy which guarantees a requested accuracy of the numerical integration and at the same time, allows to discriminate stiff zones from regular ones; hence, smaller time steps correspond to stiffer behaviors.It is a high order method (formally of order 5, which at worst might be reduced to 3) and thus, all error coming from the time integration will be bounded by the one due to the splitting procedure itself.
Nevertheless, this high order method is achieved thanks to an implicit Runge-Kutta scheme, this means that in a general case, nonlinear systems must be solved throughout the time integration process.Even if the solving system tools are highly optimized (which are based on modified Newton's methods), these procedures become very expensive for large systems and important memory requirements are needed in order to carry out these computations.As a consequence, the size of the system of equations to be solved is terribly limited by the computing resources.However, in a splitting scheme context, we easily overcome this difficulty because the reactive term of (2) is a system of ODEs without spatial coupling.Therefore, a local approach node by node is adopted where the memory requirements are only set by the number of local unknowns, which normally does not exceed conventional memory resources.Even more, this approach also allows a straightforward parallelization of computations where no data exchange is needed among nodes (see for example the numerical implementations achieved in [10]).
Another very important feature of this strategy is that precious computation time is saved because we adapt the time integration step only at nodes where the reaction phenomenon takes place.For multi-scale reaction waves, this happens in a very low percentage of the spatial domain, normally only in the neighborhood of the wavefront.Therefore, larger time steps are considered at nodes with a chemistry at (partial) equilibrium.This would not be possible if we integrated the entire reaction-diffusion system (2) at once.Time Integration of the Diffusion: ROCK4.If we now consider ROCK4 [1], we recall that one of the most important advantages of such method is its explicit character, hence the simplicity of its implementation.In fact, no sophisticated Linear Algebra tools are needed (no resolution of linear systems required) and thus, the resolution is based on simple matrix-vector products.Nevertheless, the computation cost relies directly on the requested quantity of such products, that is the number of internal stages s needed over one time integration step ∆t.The memory requirements are also reduced as a consequence of its explicit scheme; nevertheless we must keep in mind that these requirements increase proportionally with the number of nodes considered over the spatial domain.
ROCK4 is formally a stabilized explicit Runge-Kutta method and such methods feature extended stability domain along the negative real axis.Therefore, in order to guarantee stability for a fixed time step ∆t, the number of stages s needed is directly related to the spectral radius ρ(∂g/∂v) (considering a general problem such as v = g(v)), since it should verify: The method then is very appropriate for diffusion problems because of the usual predominance of negative real eigenvalues for which the method is efficiently stable.A very suitable example is the linear diagonal diffusion problem (3) with only negative real eigenvalues and constant spectral radius ρ(B).In our particular applications, the diffusive phenomenon has a leading role of propagator of perturbations over the (partial) equilibrium nodes that result on excitation of the reactive schemes and thus, the propagation of the reaction wave.The resulting self-similar character implies that the number of stages needed will remain practically constant throughout the evolution of the phenomenon.The spectral radius must be previously estimated (for example, using the Gershgoring theorem or even numerically, as proposed by the ROCK4 solver by means of a nonlinear power method).Once again, the implementation of this diffusion solver over the entire reactiondiffusion system (2) will not be appropriate under neither theoretical nor practical considerations, and highlights the inherited advantages of operator splitting.ROCK4 is also a high order method (order 4); therefore, the theoretical operator splitting analysis rest valid and the overall time integration errors are mainly due to the splitting scheme, where all the inner reaction and diffusion time scales are properly solved by these high order dedicated solvers.

Mesh Refinement Technique
We are concerned with the propagation of reacting wavefronts, hence important reactive activity as well as steep spatial gradients are localized phenomena.This implies that if we consider the resolution of reactive problem (4), a considerable amount of computing time is spent on nodes that are practically at (partial) equilibrium (see for example a precise computing time evaluation in [10]).Moreover, there is no need to represent these quasi-stationary regions with the same spatial discretization needed to describe the reacting wavefront, so that the diffusion problem (3) might also be solved over a smaller number of nodes.An adapted mesh obtained by a multiresolution process [2,15] then turn out to be a very convenient solution to overcome these difficulties.
Therefore, let us consider a set of nested spatial grids for j = 0, 1, • • • , J from the coarsest to the finest one.A multiresolution transformation allows to represent a discretized function as values on a coarser grid plus a series of local estimates at different levels of such nested grids.These estimates are related to the difference of values obtained by inter-level transformations.The theoretical background of such configuration (see [3]) states that wavelet coefficients can be then defined as prediction errors, and they will retain these estimates when going from a coarse to a finer grid.Hence, the main idea is to use the decay of the wavelet coefficients to obtain information on local regularity of the solution: lower wavelet coefficients are associated to local regular spatial configurations and vice-versa.
This representation yields to a thresholding process that builds dynamically the corresponding adapted grid on which the solutions are represented; then the error committed by the multiresolution transformation is proportional to ε, where ε is a threshold parameter [14,3].Thus, numerical experiments, based on this error estimate and on the splitting ones, allow one to properly choose the various simulation parameters in order to predict the level of accuracy of the simulation.As a consequence, a more precise control of error can be drawn out of this optimal combination of methods.
In this last section, we present some numerical illustrations of the proposed strategy.A problem coming from the nonlinear chemical dynamics is described and treated.The performance of the method is discussed in the context of 2D and 3D simulations.

Mathematical Model of Study
We are concerned with the numerical approximation of a model of the Belousov-Zhabotinski reaction, a catalyzed oxidation of an organic species by acid bromated ion (for more details and illustrations, see [11]).We thus consider the model introduced in [12] and coming from the classic of Field, Koros and Noyes (FKN) (1972), which takes into account three species: HBrO 2 (hypobromous acid), 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 partial differential equations: with diffusion coefficients D a , D b and D c , and some real positive parameters f , small q, and small , µ, such that µ .The dynamical system associated with this system models reactive excitable media with a large time scale spectrum (see [12] for more details).Moreover, the spatial configuration with addition of diffusion generates propagating wavefronts with steep spatial gradients.Hence, this model presents all the difficulties associated with a stiff multi-scale configuration.The advantages of applying a splitting strategy to these models have already been studied and presented in [8].In what follows, we will consider 2D and 3D configurations of problem (7).

2D BZ Equation
We first consider the 2D application of problem (7) with homogeneous Neumann boundary conditions and the following parameters, taken from a preliminary study [8]: = 10 −2 , µ = 10 −5 , f = 1.6 and 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 phenomenon is studied over a time domain of [0, 4] and a space region of [0, 1] × [0, 1].We define the reference or quasi-exact solution as the resolution of the coupled reactiondiffusion problem (7) on an uniform mesh of 256 × 256 performed by Radau5 with very fine tolerances.The main limitation to perform such computation on finer grids comes from the important memory requirements of Radau5.
Let us consider an application of the proposed MR/Splitting numerical strategy with 8 nested dyadic grids with N = 2 2×8 = 65536 = 256 × 256 cells on the finest grid j = J = 8.The time integration method uses the RDR Strang S ∆t 2 scheme with Radau5 for the time integration of the reaction term and ROCK4 for the diffusive part.The chosen splitting time step ∆t ≈ 3.9 × 10 −3 guarantees a fine resolution of the wave speed (relative error lower than 0.2%).The spiral waves simulated can be seen into Figure 1, where colors represent the grid levels at t = 2, when the wave is fully developed, and at t = 4, after a complete rotation period: the adapted grids are tightened around the stiff regions and clearly propagate along the waves.L 2 -error estimates show that for ε ≤ 10 −2 , the multiresolution errors become negligible compared to the splitting ones.Figure 2 shows these error estimates: proportionality to ε is verified whereas the splitting error depends only on the splitting time step.Nevertheless, this numerical strategy becomes really interesting on behalf of more accurate simulations of phenomena requiring very fine spatial resolution, for which standard strategies are simply unfeasible with conventional computing resources.Therefore, we study now the same problem but with finer spatial discretizations on the finest grid at j = J.Table 1 summarizes these results for the threshold value ε = 10 −2 .The data compression (DC) is defined as one minus the ratio between the number of cells on the adapted grid (AG) and those on the finest uniform grid (FG), expressing the whole as a percentage: Data compression increases with the number of levels as the space scales present in the problem are better discriminated by finer spatial resolutions.
In order to take into account the memory requirements of each resolution strategy for a fine spatial resolution of 1024 × 1024, we estimate the array size of the working space needed by Radau5 and ROCK4: where W 1 and W 2 are the number of unknowns solved by Radau5 and ROCK4.In the case of an uniform mesh, the total number of unknowns is W = 3×1024× 1024 ≈ 3.15 × 10 6 and thus, the global size L required for each solver is: 1. Quasi-exact: W 1 = W ≈ 3.15 × 10 6 and L = L 1 ≈ 4 × 10 13 .Therefore, on the one hand, it is hopeless trying to solve problem ( 7) with the quasi-exact strategy for these very fine discretizations, at least with standard computing resources.For instance, for these 2D simulations we have used an Intel(R) Core(TM)2 processor of 2 GHz with memory capacity of 16 Gb.And on the other hand, also a splitting strategy becomes more difficult to implement since the diffusion term is solved considering the entire spatial domain at once.A major advantage of the proposed numerical strategy is the possibility of representing results in a highly compressed way.As an example, Figure 3 shows the adapted grids obtained for the 1024 × 1024 configuration with ε = 10 −2 .

3D BZ Equation
We consider problem (7), now in a 3D configuration with the same parameters considered in the 2D case for a time domain of [0, 2] and in a space region of [0, 1]×[0, 1]×[0, 1].In order to explore the feasibility and potential advantages of the method, let us consider 9 nested dyadic grids with N = 2 3×9 = 134217728 = 512 × 512 × 512 cells on the finest grid j = J = 9.For these 3D implementations, we used an AMD-Shanghai processor of 2.7 GHz with memory capacity of 32 Gb.Then,with a threshold value of ε = 10 −1 and a splitting time step ∆t ≈ 7.8 × 10 −3 , the implementation of the proposed numerical strategy features data compressions that goes from 95.54% for the initial condition, 89.62% at t = 1 when the scroll waves are fully developed and 87.05% at final time t = 2.
For this configuration, a two times larger splitting time step is chosen in order to have splitting and multiresolution errors potentially of the same order.Smaller threshold values yield larger simulation domains which are not longer feasible with the considered computing resource and the current state of development of the code.Figure 4 shows the evolution of the finest grid of the adapted grid during the period of study.We distinguish clearly the development of the scroll wave where colors represent the values of variable a: the finest regions correspond to the neighborhood of the wavefront.
Performing the same comparison concerning memory requirements, the total number of unknowns for this case is W = 3 × 512 × 512 × 512 ≈ 4.03 × 10 8 and the global size of L required by each solver is:

Conclusions
The present work proposes then a simple and solid numerical strategy that couples adaptive multiresolution techniques with a new operator splitting strategy for multi-scale reactions waves modeled by stiff reaction-diffusion systems.It considers on the one hand, dedicated high order time integration methods to properly solve the entire spectrum of temporal scales of both the reaction and the diffusion part; and on the other hand, an adaptive multiresolution technique to represent and treat more accurately local spatial gradients associated to the wave front.The resulting highly compressed data representations as well as the accurate and feasible resolution of these stiff phenomena prove that large computational domains previously out of reach can be successfully simulated with conventional computing resources.
We have focused our attention on reaction-diffusion systems in order to settle the foundations for simulation of more complex phenomena with fully convection-reaction-diffusion systems or more detailed models.In particular, we search to extend this strategy to multi-scale problems in combustion domain, where other numerical techniques have already been proposed (see for instance [5,16]) and where there is a continuous research of new and more efficient methods.Therefore, an important amount of work is still in progress concerning on the one hand, programming features such as data structures, optimized routines and parallelization strategies.And on the other hand, numerical analysis of theoretical aspects, which may surely allow the introduction of dynamical error estimators and even adapting splitting time steps for more general multi-scale phenomena; these are particular topics of our current research.