A second order directional split exponential integrator for systems of advection–diffusion–reaction equations

We propose a second order exponential scheme suitable for two-component coupled systems of stiff evolutionary advection–diffusion–reaction equations in two and three space dimensions. It is based on a directional splitting of the involved matrix functions, which allows for a simple yet efficient implementation through the computation of small-sized exponential-like functions and tensor-matrix products. The procedure straightforwardly extends to the case of an arbitrary number of components and to any space dimension. Several numerical examples in 2D and 3D with physically relevant (advective) Schnakenberg, FitzHugh–Nagumo, DIB, and advective Brusselator models clearly show the advantage of the approach against state-of-the-art techniques.


Introduction
Efficiently solving systems of Advection-Diffusion-Reaction (ADR) equations is of paramount importance in practical applications.In fact, these kinds of models effectively capture many physical, chemical, and biological phenomena, such as electrodeposition [9], biochemical reactions [27,33] and electric current flows [18], among the others.From a numerical point of view, the computation of approximated solutions gives rise to many challenges.For example, in this context there is often the intrinsic need of a fine spatial grid in order to correctly capture the dynamics of the system (the formation of the so-called Turing patterns [17], for instance).This, in turn, typically translates into working in a stiff regime [20], which requires a careful choice of the underlying numerical schemes.

˚Corresponding author
Email addresses: marco.caliari@univr.it(Marco Caliari), fabio.cassini@univr.it(Fabio Cassini) arXiv:2307.10684v1[math.NA] 20 Jul 2023 In particular, in this work we focus on the numerical integration of twocomponents systems of ADR equations in the form although the method presented later on in the manuscript easily extends to an arbitrary number of components.Here u, v : r0, T s ˆΩ Ă R d Ñ R represent the unknowns, K u , K v are linear advection-diffusion operators, while g u , g v are the nonlinear reaction terms.The choice of the latter basically determines the model and phenomenon under investigation.We also assume that the spatial domain Ω is the Cartesian product of one-dimensional intervals, that is Ω " ra 1 , b 1 s ˆ¨¨¨r a d , b d s.The system is finally completed with appropriate initial and boundary conditions.Throughout the manuscript, we will employ homogeneous Neumann boundary conditions, which appear to be widely used in the literature, but other types (such as Dirichlet or periodic) could be considered as well.We introduce a spatial grid of size n 1 ˆ¨¨¨ˆn d and apply the Method Of Lines (MOL) to model (1).We then obtain a system of stiff ordinary differential equations # u 1 ptq " K u uptq `gu puptq, vptqq, v 1 ptq " K v vptq `gv puptq, vptqq, where K u and K v are matrices of size N ˆN , with N " n 1 ¨¨¨n d , that discretize the linear operators K u and K v , respectively.The choice of a structured grid leads to matrices K u and K v with Kronecker sum structure.This means, for a matrix K P C N ˆN , that it can be decomposed as where Here, I µ is the identity matrix of size n µ ˆnµ and A µ is a matrix of small size n µ ˆnµ .In our context, the matrix A µ represents the one-dimensional discretization of the linear differential operator along the direction µ.The symbol b denotes the standard Kronecker product between matrices.Several methods can be employed for the time integration of system (2).Due to the presence of stiffness, an attractive choice are integrators of exponential or implicit type, useful to alleviate the stability restrictions on the time step size that come with explicit methods.Of course, in this case ad-hoc techniques have to be adopted in order to retain efficiency of computation.For instance, among schemes of exponential type, Lawson methods (also known as Integration Factor methods) are appealing [24], since they require just the action of the exponential matrix function and thus they can directly exploit the Kronecker sum structure of the system (see Reference [10] and Reference [11] for some insights on GPU scalability).On the other hand, it is well-known that they may suffer from a bad convergence behaviour, in particular when the boundary conditions are non-periodic (see Reference [21]).Different exponential integrators [22] rely on Krylov methods for the approximation of the matrix exponential or related functions with a fixed size of the Krylov space (see Reference [7]) or with incomplete orthogonalization (see References [19,26]).More recent techniques that exploit the underlying Kronecker sum structure are presented in References [12,16,17,28], while low degree rational approximations having real distinct poles, together with a dimensional splitting, are proposed in Reference [4].Concerning implicit methods, low-rank or other Model Order Reduction techniques together with matrix formulation (i.e, in the two-dimensional case d " 2) can be applied to classical schemes such as IMplicit EXplicit methods (see References [2,17]).We finally mention also the possibility to employ explicit schemes such as a Total Variation Diminishing explicit Runge-Kutta scheme (see Reference [32]), provided that a sufficiently small time step size is taken.
Following the preliminary yet valuable results obtained in Reference [10] for single equations, in this manuscript we propose an exponential scheme of second order for system (2).It is a directional split version of the well-known exponential Runge-Kutta method of order two ETD2RK, and it approximates the required four actions of matrix ϕ functions at the needed order by means of tensor-matrix products (using the so-called Tucker operator).This approach allows for a simple yet very efficient implementation, thanks to the high performance level 3 BLAS available for basically any kind of modern computer system.
Although the method can be straightforwardly implemented in any spatial dimension, we test it on several two-dimensional and three-dimensional systems of particular interest, namely The remaining part of paper is organized as follows.In Section 2 we briefly recall integrators of exponential type and the tensor framework needed to implement them efficiently.In Section 3 we present ETD2RKds, which is the proposed integrator with directional splitting.We proceed in Section 4 by testing it on a wide variety of two-dimensional and three-dimensional numerical examples with popular models in the literature (the DIB, the FitzHugh-Nagumo, the Schnakenberg, and the advective Brusselator models), highlighting the effectiveness and the efficiency of the approach.Finally we draw the conclusions in Section 5.

Exponential integrators for systems of partial differential equations
The starting point for the exponential integrators under consideration is the variation-of-constants formula applied to system (2) with time step size τ " t n`1 ´tn , i.e., upt n`1 q " e τ K u upt n q `τ ż 1 0 e p1´θqτ K u g u pupt n `θτ q, vpt n `θτ qqdθ, vpt n`1 q " e τ K v vpt n q `τ ż 1 0 e p1´θqτ K v g v pupt n `θτ q, vpt n `θτ qqdθ. ( We then apply the trapezoidal quadrature rule, by introducing an explicit intermediate stage to approximate the evaluation at t n`1 (see also Reference [10]).
Putting everything together, we get the so-called Lawson2b scheme The method can be recovered also by multiplying the two equations of system (2) by the integration factors e ´tK u and e ´tK v , respectively, and applying the explicit Heun method (see for instance Reference [24]).Scheme (4) can be efficiently implemented by exploiting the Kronecker sum structure of the involved matrices.Indeed, for a matrix τ K with Kronecker sum structure, we can compute the action of its exponential on the vector w P C N by employing the representation e τ K w " vec `W ˆ1 e τ A1 ˆ2 ¨¨¨ˆd e τ A d ˘.
Here, W is an order-d tensor of size n 1 ˆ¨¨¨ˆn d , while vec is the operator which stacks the columns of the input tensor into a suitable column vector such that vecpW q " w.Finally, ˆµ denotes the µ-mode product (also known as tensormatrix product), and it can be realized in an efficient way by means of high performance level 3 BLAS operations at a computational cost Opn µ N q.The concatenation of d different µ-mode products, that we will denote as Ś d µ"1 , is called Tucker operator, and clearly has computational cost Oppn 1 `¨¨¨ǹ d qN q.The exponentials of the matrices τ A µ can be efficiently computed by a polynomial or rational approximation (see References [1,3,15,30]) together with a scaling and squaring algorithm, or even by a spectral decomposition (see Reference [17]).We refer the reader to References [11,14] for more details about formula (5), the Tucker operator and related operations (implemented in a free MATLAB package1 ).Here we just emphasize that the Tucker operator (5) in two dimensions (that is, K " A 2 ' A 1 and W P C n1ˆn2 ) coincides with the well-known formula e τ K w " vec ´eτA1 W `eτA2 ˘T¯.
If we instead interpolate the nonlinear functions g u and g v in formula (3) by a polynomial of degree one in θ at 0 and c ď 1, and we introduce an exponential Euler stage for the approximation at θ " c, we get a well-known class of exponential Runge-Kutta schemes of order two The matrix ϕ functions, for a generic matrix X P C N ˆN , are defined by Several techniques are available in the literature to compute the matrix ϕ functions (see References [6,10,25] for methods based on rational or polynomial approximations or quadrature formulas) or their action to vectors (see References [13,19,26] for Krylov methods with incomplete orthogonalization).Unfortunately, even if X " τ K, property (5) does not extend to the matrix ϕ functions.However, it is still possible exploit the Kronecker sum structure of K to compute the action of the matrix ϕ function over a vector w efficiently.In fact, if we introduce a quadrature formula for integral definition (7), we can approximate the desired quantity by computing the actions e p1´θiqτ K w using Tucker operators, where θ i are the quadrature nodes.In this way, the approximation via quadrature formula requires one Tucker operator for each quadrature node.The total number of quadrature nodes can be kept reasonably low by employing a modified scaling and squaring technique for the ϕ functions.Each iteration of the squaring procedure requires then the action of a matrix exponential, which can be realized again by a Tucker operator.We refer to References [12,16] for more details, also regarding different choices of the quadrature formula and of the scaling parameter.

Directional splitting of the matrix ϕ functions
As already mentioned, the splitting property e τ K w " e τ Ab1 ¨¨¨e τ A bd w, which is at the basis of equivalence (5), does not extend directly to ϕ functions.However, it has been observed in Reference [10] that a simple directional splitting of the matrix ϕ functions, that is gives an approximation compatible with second order integrators.The advantage of using the above approximation is that it can be efficiently computed by thanks to the properties of the Kronecker product.The approximation introduced by the splitting can affect in general the global error when applying the exponential integrator, but the preliminary numerical experiments reported in Reference [10] show that the overall computational cost required to reach the same level of accuracy is smaller.This is possible since only small sized matrix ϕ functions have to be computed, and the Tucker operator (which can be efficiently realized as mentioned in Section 2) has to be applied.
Concerning specifically class ( 6), an attractive choice is c " 1 for which the corresponding method, known as ETD2RK (see Reference [6]), requires the computations of four actions of the matrix ϕ functions instead of six.Then, embedding the splitting formula (8) in the ETD2RK integrator, we obtain a scheme that we call "directional split ETD2RK", and more concisely ETD2RKds.If we now write explicitly K u and K v in formula (6) as we finally obtain a complete expression of the proposed scheme, i.e., Notice that in the ETD2RKds scheme above, we are able to compute the actions of K u and K v in tensor form thanks to the equivalences Each of these two evaluations costs as a Tucker operator.Clearly, scheme (9) can be straightforwardly generalized to systems of semilinear equations with more than two components and Kronecker sum structure.Also, note that the case d " 2 can be written down explicitly without introducing tensor notation as T Remark that for ETD2RKds method, U n2 and V n2 turn out to be approximations of U pt n`1 q and V pt n`1 q, respectively, by the directional split version of the well-known exponential Euler method.Therefore, the scheme would allow for a construction of a pair of embedded methods, useful in variable step size integration.Finally, we also notice that the intermediate stages, as well the approximations of the discretizations of the components u and v of the system at time t n`1 , could be computed simultaneously in parallel.

Numerical experiments
In this section we show the effectiveness of the proposed directional split exponential integrator ETD2RKds by applying it to several two-dimensional and three-dimensional models of great interest.All the problems are discretized in space with standard second order finite differences on a grid of equispaced discretization points and n 1 " . . ." n d " n (for a total number of degrees of freedom equal to N " n d ).The computation of the small sized matrix exponentials is performed by the built-in MATLAB function expm, while the matrix ϕ functions are approximated by the phiquad function2 introduced in Reference [10].It performs the modified scaling and squaring algorithm and approximates integral definition (7) by the Gauss-Legendre-Lobatto formula.Whenever a random perturbation of an equilibrium is considered as initial condition, we fix the seed of the random number generator by the MATLAB command rng('default').
We compare our scheme with some other exponential or implicit integrators, developed both for the two-dimensional and three-dimensional case.They are all implemented in MATLAB language, and the arising linear systems or matrix functions are always computed by direct methods (that is, without using iterative procedures).In particular, we consider • an implementation of Lawson2b in which the action of the matrix exponential is realized by a Tucker operator (which does not introduce any directional splitting error) • an implementation of ETD2RK in which the action of the matrix ϕ functions is realized by the phiks routine3 (see Reference [13]) at a sufficiently high accuracy (depending on the input tolerance) • the Exponential Time Differencing Real Distinct Poles Integrating Factor (ETD-RDP-IF4 , see Reference [4]) method, in which the exponential functions are approximated by a low-degree Padé expansion (with a dimensional splitting for the solution of the arising linear systems) • the MATLAB solver ode23tb, which implements a variable step size diagonally implicit Runge-Kutta pair 2(3) (DIRK23, see Reference [23]), suggested for stiff problems with "crude error tolerances", fed with the exact Jacobian of the system.
Moreover, we test also the performance of the MATLAB solver ode23, which implements a variable time step explicit Runge-Kutta pair 3(2) (RK32, see Reference [8]).This method is considered in the comparisons in order to be sure that the examples, with the selected parameters and space discretizations, are in a stiff regime in which explicit methods suffer of a time step size restriction.Whenever calling MATLAB built-in integrators, since we are just interested in the solution at final times, we use a proper output function (through the option OutputFcn) in order not to save the solutions at intermediate steps and hence waste memory.In all the numerical experiments we first compare the performances of the methods in reaching a common range of accuracies with respect to a reference solution computed with ETD2RKds and a sufficiently large number of time steps.We therefore select some sequences of given time steps (for the fixed time step size implementations) or tolerances (for the variable time step size implementations) for the different methods and measure the overall needed wallclock time.The error is always measured as the 2-norm of the relative Frobenius norm of the solutions U and V at final time.Then, we run our proposed exponential scheme up to a larger final time, in order to observe the correct evolution dynamics (such as the pattern formation for the variable u).
All the experiments are performed on an Intel ® Core ™ i7-10750H CPU with six physical cores and 16GB of RAM.As a software, we use MathWorks MATLAB ® R2022a.All the codes to reproduce the examples can be found in a maintained GitHub repository5 .

Two-dimensional morpho-chemical DIB model
We consider the morpho-chemical DIB model (see References [2,17]) 1 p1 `av 2 uqp1 ´vqp1 ´av 3 p1 ´vqq ´av 4 vp1 `av 5 uqp1 `av 3 vqq, (10) in Ω " r0, 20s 2 .The model describes the electrodeposition process for metal growth.The unknowns u and v represent the morphology of the metal deposit and its surface chemical composition, respectively.The parameters, taken from Reference [2], are δ u " 1, ρ " 25{4, a u 1 " 10, a u 2 " 1, a u 3 " 66, a u 4 " 0.5, δ v " 20, , and a v 5 " 1.5.The particular choice of the parameter a v 4 yields the equilibrium pu e , v e q " p0, a u 4 q susceptible of Turing instability.The initial conditions are u 0 px 1 , x 2 q " u e `10 ´5 ¨Up0, 1q and v 0 px 1 , x 2 q " v e `10 ´5 ¨Up0, 1q, where Up0, 1q denotes the uniformly distributed random variable in p0, 1q.The spatial domain is discretized with a grid of N " 200 2 points.We first integrate the system up to the final time T " 2.5 with the methods described above, and a number of time steps chosen, for each fixed step size method, as reported in Table 1.Concerning the MATLAB variable step size Runge-Kutta integrators, Table 1 reports the values of the options AbsTol and RelTol of the functions ode23tb and ode23.The results are presented in Figure 1a The plot is clearly divided into two parts.On the right we have the more expensive methods, with an overall wall-clock time of roughly 100 seconds.In particular, the ETD-RDP-IF method (specifically designed for these type of models) and the ETD2RK method (with input tolerance set to 1e´3) do not perform much better than an explicit method.This means that for the considered parameters and spatial discretization the problem does not appear to be excessively stiff.In fact, the implicit method is roughly five times faster than the explicit one and with performances similar to the Lawson2b method.The ETD2RKds method is the fastest among all and, in particular, ten times faster than the corresponding version with "accurate" ϕ functions ETD2RK.Figure 1a: Results for the simulation of the DIB model (10) with N " 200 2 spatial discretization points.The number of time steps (or tolerance) for each integrator is reported in Table 1.
The final simulation time is T " 2.5.The simulation wall-clock time is 194 seconds.We observe the same pattern also by using a smaller time step τ " 1.25e´3.However, if a larger time step size τ " 2e´2 is used, a "wrong" Turing pattern appears (right plot).
Then, we integrate the system up to T " 100 with ETD2RKds and time step size equal to the largest one in the previous simulation, that is τ " 2e´3.The observed steady Turing pattern is shown in Figure 1b, left plot.The wall-clock time needed for this experiment is 194 seconds.The same pattern shows up if we use a smaller time step size (τ " 1.25e´3, corresponding to the smallest time step size in the previous simulation).It was observed in Reference [17] that a spatial domain not large enough or a number of discretization points not high enough can prevent from observing the correct Turing pattern.Here, where the domain and the number of discretization points have been properly chosen, we also notice that if we take a ten times larger time step size, a steady pattern still appears.However, due to accuracy issues it is clearly not the correct one (see Figure 1b, right plot).

Two-dimensional Schnakenberg model
We consider the Schnakenberg model (see References [2,5,17,31]) in the spatial domain Ω " r0, 1s 2 .The unknowns u and v represent two chemical concentrations in autocatalytic reactions.The parameters, taken from Reference [2], are δ u " 1, δ v " 10, ρ " 1000, a u " 0.1, and a v " 0.9.The equilibrium pu e , v e q " pa u `av , a v {pa u `av q 2 q is susceptible of Turing instability.The initial data are u 0 " u e `10 ´5 ¨Up0, 1q and v 0 " v e `10 ´5 ¨Up0, 1q.The spatial domain is discretized with a grid of N " 150 2 points.The final simulation time is T " 0.25.The results are resumed in Figure 2a.The numbers of time steps (or tolerances) for each integrator are collected in  In this case, we observe that the ETD2RK (with input tolerance 1e´4) and ETD-RDP-IF schemes basically show the same results, being however not competitive with the other methods.An overall neat advantage be obtained by applying the proposed methods ETD2RKds.
Finally, we perform a simulation just with ETD2RKds up to final time T " 2 with 2000 time steps.The overall wall-clock time is roughly 3 seconds.The observed pattern is in agreement with that reported in the literature (compare Figure 2b and Reference [2, Fig. 8 (a)]).2. The final simulation time is T " 0.25.

Two-dimensional FitzHugh-Nagumo model
We consider the FitzHugh-Nagumo model (see Reference [2]) in Ω " r0, πs 2 .The model describes the flow of an electric current through a nerve fiber.The unknowns u and v represent the electric potential and the recovery variable, respectively.The parameters, taken from Reference [2], are δ u " 1, δ v " 42.1887, ρ " 65.731, a v 1 " 11, and a v 2 " 0.1.With these choices, the equilibrium pu e , v e q " p0, 0q is susceptible of Turing instability.The initial solutions are u 0 px 1 , x 2 q " 10 ´3 ¨Up0, 1q and v 0 px 1 , x 2 q " 10 ´3 ¨Up0, 1q.The spatial domain is discretized by a grid of N " 100 2 points.In the first experiment we simulate up to T " 10.The numbers of time steps (or tolerances) for each integrator are collected in Table 3.The results are presented in Figure 3a.
First of all, we notice that the Lawson2b and ETD-RDP-IF methods require many more time steps than the exponential Runge-Kutta methods, and still they cannot reach the same accuracies.Their wall-clock time is in fact at least one order of magnitude larger than the most efficient methods.On the   other hand, as in the previous examples, the ETD2RKds method is not heavily affected by the direction splitting error.Indeed, it reaches almost the same accuracies obtained by the ETD2RK method (with input tolerance 1e´3), but with a much smaller computational time.
In the second experiment, we set the final simulation time to T " 50 and integrate the problem with the ETD2RKds method with 20000 time steps, for an overall computational time of 12 seconds.Again, the obtained pattern agrees with that reported in the literature (compare Figure 3b and Reference [2, Fig. 4 (a)]), i.e., a square pattern corresponding to the cosine modes p4, 4q.

Three-dimensional advective Schnakenberg system
The diffusion-reaction Schnakenberg system models the limit cycle behaviors of two-component chemical reactions.An advection term was introduced Figure 3a: Results for the simulation of 2D FitzHugh-Nagumo model (12) with N " 100 2 spatial discretization points.The number of time steps (or tolerance) for each integrator is reported in Table 3.The final simulation time is T " 10.
in Reference [29] in order to study its effect on Turing patterns (see also Reference [7]).We therefore study the system The computational domain is Ω " r0, 1s 3 , the advection parameters are α u " α v " 0.01, the diffusion parameters are δ u " 0.05 and δ v " 1, the reaction parameter is ρ " 100, and the concentration parameters are a u " 0.1305 and a v " 0.7695, respectively.As initial solution we take up0, x 1 , x 2 , x 3 q " a u `av `10 ´5e ´100ppx1´1{3q 2 `px2´1{2q 2 `px3´1{3q 2 q , vp0, x 1 , x 2 , x 3 q " a v pa u `av q 2 , which corresponds to a small deviation from the steady state solution pu e , v e q " pa u `av , a v {pa u `av q 2 q.We integrate the system up to final time T " 0.4, with the number of time steps (or tolerance) given in Table 4.The degrees of freedom in space are N " 80 3 .The results are given in Figure 4a.First of all, despite the aforementioned cares in the usage of the DIRK23 method from the MATLAB ODE Suite, the ode23tb function stopped to run after some seconds due to too large memory requirements.Concerning the other methods, we clearly see that the exponential Runge-Kutta (with input tolerance 1e´5) and the Lawson schemes perform way better then the other ones.Overall, the best method turns out to be again ETD2RKds.
Then, we repeat the experiment up to final time T " 0.8, with the same number of degrees of freedom in space and using as time integrator ETD2RKds with τ " 8e´3 (the largest one in the previous simulation).As already observed in the literature (see Reference [7]), the Turing instability leads to a spot-like pattern, depicted in Figure 4b.The overall simulation time is roughly 3 seconds.

Three-dimensional FitzHugh-Nagumo model
We consider again the FitzHugh-Nagumo model but now in a three-dimensional domain Ω " r0, πs 3 .The parameters δ u " 1, δ v " 42.1887, ρ " 24.649, a v 1 " 11, and a v 2 " 0.1 were obtained by following the analysis performed in Reference [18] in order to achieve a square pattern with modes p2, 2, 2q.With these choices, the equilibrium pu e , v e q " p0, 0q is susceptible of Turing instability.The initial solutions are u 0 px 1 , x 2 , x 3 q " 10 ´3 Üp0, 1q and v 0 px 1 , x 2 , x 3 q " 10 ´3 ¨Up0, 1q.The spatial domain is discretized by a grid of N " 64 3 points.We simulate up to the final time T " 10 with the number of time steps (or tolerance) given in  As in the previous example, the implicit Runge-Kutta method DIRK23 could not terminate due to excessive memory requirements.Among the other methods, ETD2RKds is again the best one, with a speedup of almost one order of magnitude with respect to the second one, that is the exponential Runge-Kutta method with the ϕ functions approximated by phiks (with input tolerance 5e´4).The RK32 and ETD-RDP-IF methods did not reach comparable accuracies within 10 4 seconds.
We then simulate the system up to T " 100 with 15000 time steps and the ETD2RKds method.We report in Figure 5b the obtained pattern which agrees with the theoretical expectation.The overall wall-clock time was about 123 seconds.

Three-dimensional advective Brusselator system
We consider the Brusselator system presented in Reference [7].The unknowns u and v represent the concentration of the activator and inhibitor in the chemical reaction, respectively.The spatial domain is Ω " r0, 1s 3 , the diffusion coefficients are δ u " 0.01, δ v " 0.02, the  14) with N " 64 3 spatial discretization points.The number of time steps for each integrator is reported in Table 5.The final simulation time is T " 10. parameter coefficients are a u 1 " 1 and a u 2 " 2, α u " α v " 0.1 and the final simulation time is T " 5. We first consider the case of a space discretization with N " 64 3 discretization points.The initial datum is upx 1 , x 2 , x 3 q " 1 sinp2πx 1 q sinp2πx 2 q sinp2πx 3 q and vpx 1 , x 2 , x 3 q " 3. The results are collected in    (15) with N " 64 3 spatial discretization points.The number of time steps (or tolerance) for each integrator is reported in Table 6.The final simulation time is T " 1.
Again the plot is split into two parts, with the Lawson2b method for the first time slightly faster then ETD2RKds.The other methods are almost one order of magnitude slower (ETD2RK was run with input tolerance 5e´6).
We then simulate the system up to T " 5 with 100 time steps with ETD2RKds, with an overall computational time of 1 second.As already observed in Reference [7], the solution approaches the equilibrium state with pu e , v e q " pa u 2 , a u 1 {a u 2 q, see Figure 6b.

Conclusions
In this paper we showed how it is possible to effectively exploit the Kronecker sum structure which arises by semidiscretizing in space two-component systems of coupled advection-diffusion-reaction equations.The proposed second order exponential-type time marching scheme, which is based on a directional splitting of the involved matrix functions and named ETD2RKds, is shown to outperform well-established techniques on a variety of physically relevant models from the literature, such as two-dimensional DIB, FitzHugh-Nagumo and Schnakenberg and three-dimensional FitzHugh-Nagumo, advective Schnakenberg and advective Brusselator models.The procedure easily generalizes to any number of components and any number of spatial dimensions.As future work we plan to perform simulations at HPC level (i.e., using server multi-core CPUs and GPUs) to further enhance the performances of the directional splitting procedure.

Figure 1b :
Figure 1b: Turing pattern for DIB model (10) obtained at final time T " 100 with N " 200 2 space discretization points with ETD2RKds and with time step size τ " 2e´3 (left plot).The simulation wall-clock time is 194 seconds.We observe the same pattern also by using a smaller time step τ " 1.25e´3.However, if a larger time step size τ " 2e´2 is used, a "wrong" Turing pattern appears (right plot).

Figure 2a :
Figure2a: Results for the simulation of the 2D Schnakenberg model(11) with N " 150 2 spatial discretization points.The number of time steps (or tolerance) for each integrator is reported in Table2.The final simulation time is T " 0.25.
pattern for 2D Schnakenberg model (11) obtained at final time T " 2 with N " 150 2 space discretization points with ETD2RKds.The time step size employed is τ " 1e´3.The simulation wall-clock time is 3 seconds.

Figure 3b :
Figure 3b: Turing pattern for 2D FitzHugh-Nagumo model (12) obtained at final time T " 50 with N " 100 2 space discretization points with ETD2RKds.The time step size employed is τ " 2.5e´3.The simulation wall-clock time is 12 seconds.

Figure 4a : 3 Figure 4b :
Figure 4a: Results for the simulation of the 3D Schnakenberg model (13) with N " 80 3 spatial discretization points.The number of time steps (or tolerance) for each integrator is reported in Table4.The final simulation time is T " 0.4.

Figure 5a :
Figure 5a: Results for the simulation of the 3D FitzHugh-Nagumo model (14) with N " 64 3 spatial discretization points.The number of time steps for each integrator is reported in Table5.The final simulation time is T " 10.

3 Figure 5b :
Figure 5b: Turing pattern for 3D FitzHugh-Nagumo model (14) obtained at final time T " 100 with N " 64 3 space discretization points and ETD2RKds.The time step size employed is τ " 6.67e´3.The wall-clock simulation time is 123 seconds.The reported slice (left plot) is the 32nd along the x 3 direction and the isosurface value (right plot) is 0.08.

Figure 6a :
Figure6a: Results for the simulation of the 3D Brusselator model(15) with N " 64 3 spatial discretization points.The number of time steps (or tolerance) for each integrator is reported in Table6.The final simulation time is T " 1.

Figure 6b :
Figure 6b: Steady state of 3D Brusselator model (15) obtained at final time T " 5 with N " 64 3 space discretization points and ETD2RKds.The time step size employed is τ " 5e´2.The simulation wall-clock time is 1 second.The reported slices are the last along the x 3 direction for the unknown u (left plot) and the unknown v (right plot).

Table 1 :
(10)mbers of time steps and tolerances used in the methods tested for 2D DIB model(10)up to final time T " 2.5.

Table 2 :
Numbers of time steps and tolerances used in the methods tested for 2D Schnakenberg model (11) up to final time T " 0.25.

Table 4 :
Numbers of time steps and tolerances used in the methods tested for 3D Schnakenberg model (13) up to final time T " 0.4.The DIRK23 method interrupted due to excessive memory requirements with tolerance 1e´1.

Table 5 .
The results are reported in Figure5a.

Table 5 :
Numbers of time steps used in the methods tested for 3D FitzHugh-Nagumo model (14) up to final time T " 10.The DIRK23 method interrupted due to excessive memory requirements with tolerance 1e´1.The RK32 and the ETD-RDP-IF methods required too much computational time.

Table 6 :
(15)ers of time steps and tolerances used in the methods tested for 3D Brusselator model(15)up to final time T " 1.The DIRK23 method interrupted due to excessive memory requirements with tolerance 1e´1.