{'fraction_non_alphanumeric': 0.031155637081981092, 'fraction_numerical': 0.015577818540990546, 'mean_word_length': 4.43413586873179, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Assessment of a transient homogeneous reactor through in situ adaptive tabulation

The development of computational models for the numerical simulation of chemically reacting flows operating in the turbulent regime requires the solution of partial differential equations that represent the balance of mass, linear momentum, chemical species, and energy. The chemical reactions of the model may involve detailed reaction mechanisms for the description of the physicochemical phenomena. One of the biggest challenges is the stiffness of the numerical simulation of these models and the nonlinear nature of species rate of reaction. This work presents a study of in situ adaptive tabulation (ISAT) technique, focusing on the accuracy, efficiency, and memory usage in the simulation of homogeneous stirred reactor models using simple and complex reaction mechanisms. The combustion of carbon monoxide with oxygen and methane with air mixtures are considered, using detailed reaction mechanisms with 4 and 53 species, 3 and 325 reactions, respectively. The results of these simulations indicate that the developed implementation of ISAT technique has a absolute global error smaller than 1 %. Moreover, ISAT technique provides gains, in terms of computational time, of up to 80% when compared with the direct integration of the full chemical kinetics. However, in terms of memory usage the present implementation of ISAT technique is found to be excessively demanding.


Introduction
Combustion is a physicochemical phenomena characterized by a sequence of chemical reactions (most of then exothermic) which converts fuel/oxidizer mixtures (fresh gases) into combustion products (burned gases). Due to the large amount of energy released in these exothermic reactions, the combustion process is widely used in the operation of industrial devices such as gas turbines, process furnaces, etc.
Design and optimization of devices which operate based on combustion processes is an important engineering task, which demands the development of computational models as predictive tools. Such computational models may require the solution of partial differential equations that represent the balance of mass, linear momentum, chemical species, and energy. These models may include a detailed kinetic mechanism for the description of physicochemical phenomena involved [28,6,22,35]. Typically, such reaction mechanisms for mixtures of hydrocarbons with air involve tenths of species, hundreds of elementary reactions, and timescales that vary up to nine orders of magnitude [36,20,21].
Concerning the numerical simulation of these models, the challenge is related to reaction rate of the chemical species, which is nonlinear and presents a strong dependence with reaction mechanism dimension. Therefore, numerical simulation of a model with a detailed reaction mechanism is computationally demanding, which justifies the development of techniques to reduce the associated computational cost.
Several techniques to reduce the computational cost of these models are available in the literature. Essentially, they may be split into two classes. The first class includes techniques that seek to reduce the dimension of reaction mechanisms, such as Quasi-Steady State Ap-proximation [8], Rate-Controlled Constrained Equilibrium [13], Computational Singular Perturbation [15], the reduction to Skeleton Mechanisms [16], Intrinsic Low-Dimensional Manifold [19], Proper Orthogonal Decomposition [11], and Invariant Constrained Equilibrium Edge Pre-image Curve [26]. The second class includes techniques that seek efficient algorithms to solve complex models, such as Look-Up Table [3], Repro-Modeling [32], Piece-Wise Reusable Implementation of Solution Mapping [31], Artificial Neural Network [12], and In Situ Adaptive Tabulation (ISAT) [24,29,17,18].
A characterization of one type of transient homogeneous reactor, using simple and complex thermochemistry, is presented in this this work. The analysis explores the effects of mixing to residence time scale ratios as reactor controlling parameter. The ISAT technique is employed to reduce the computational time spent to approximate the solution of the equations that governs the evolution of reactive mixture in this transient reactor. This work also investigates accuracy, performance, and memory usage of the present ISAT technique implementation. The accuracy analysis uses local and global metrics to quantify the error incurred by ISAT compared with a reference solution and investigates the statistical seed influence on the results. The performance is investigated by comparing the computational time gain obtained by ISAT with respect to a standard numerical integration procedure. The memory usage analysis quantifies the amount of memory used by ISAT implementation.
This work is organized as follows: the transient reactor of interest, its governing equation as well the numerical scheme used are presented in the second section. The third section presents the basic theory of ISAT technique. The fourth section discusses about the reactor configurations studied and ISAT implementation issues. Finally, in the fifth section, the main conclusions are summarized and a path for a future work is suggested.
2 Transient homogeneous reactor model

The geometry of reactive systems
In the framework of the transported probability density function (PDF) models [23], a reactive flow may be described by an ensemble of stochastic particles, which mimic the behavior of the fluid system.
Consider a transient spatially homogeneous reactive mixture evolving adiabatically and at constant pressure in a continuous flow reactor. The thermodynamical state of a fluid particle in this reactor may be completely determined by the mass fraction Y i , where i = 1, · · · , n s , of the n s chemical species and the specific enthalpy h, which can be lumped in the composition vector defined as where the superscript T denotes the transposition operation. One should note that, due to the invariance of the system number of atoms, which ensures the total conservation of the mass, the components of vector φ are not linearly independent. The composition vector of each particle in this flow reactor evolves according to where t denotes the time, Γ is the rate of change due to mixing, and S is the rate of change associated with the chemical reactions. Integrating Eq.(2) from an initial time t 0 to a time t one obtain the reaction mapping which is the solution of Eq.(2) starting from the initial composition φ(t 0 ) = φ 0 . The reaction mapping corresponds to a trajectory in composition space, which, for large values of t, tends to the equilibrium composition for the given enthalpy and mass fractions on φ 0 . The composition space is the (n s +1)-dimensional Euclidean space where the first direction is associated with the enthalpy and the remaining n s are related to the chemical species.

Pairwise mixing stirred reactor
The classical Partially Stirred Reactor (PaSR), used in [4,27], describes Γ by the interaction by exchange with the mean (IEM) micromixing model [7], but, for the purpose of testing a thermochemistry reduction technique, it is desirable to employ a mixing model that leads to a composition region accessed during the solution process which is "wider" than that provided by the IEM model. A modified version of PaSR model called Pairwise Mixing Stirred Reactor (PMSR), [24], is designed to yield a much larger accessed region, and, hence, should provide a stringent test on the ability of ISAT technique to yield a reduction in computational time. In the PMSR model the reactor consists of an even number n p of particles, initially arranged in pairs (j 1 , j 2 ) such that the particles (1, 2), (3,4), · · · , (n p − 1, n p ) are partners. Given a time step, ∆t, for each discrete time k∆t, where k is an integer, the model is characterized by three types of events: inflow, outflow, and pairing. The inflow and outflow events consist of randomly selecting n in = ceil 1 2 n p (∆t/τ r ) pairs of particles, being τ r the residence time within the reactor, and exchanging their thermodynamical properties by the properties of a prescribed inflow. The pairing event consists of randomly selecting for pairing a number of pairs of particles, different from the inflow particles, equal to n pair = ceil 1 2 n p (∆t/τ p ) , being τ p the pairwise time. Then the chosen particles (inflow/outflow and paring) are randomly shuffled. Between these discrete times, the pairs of particles (j 1 , j 2 ) evolve according to dφ dt dφ dt where τ m is the mixing time.

Numerical Integration
An operator splitting technique [37] is employed to numerically integrate Eq.(2). The overall process of integration via operator splitting technique can be represented as where given an initial composition φ 0 and a time step ∆t, the first fractional step integrates the pure mixing system, to obtain φ (j) mix (t + ∆t). Then, the pure chemical reaction system, is integrated from an initial composition φ mix (t + ∆t) over a time step ∆t and gives φ(t + ∆t).
The operator splitting technique allows to solve each term in the evolution equation, Eq.(2), separately, using specific efficient numerical methods to treat the particular features inherent to the physical phenomenon modeled by each term [7].

Linearized reaction mapping
Consider a composition φ and an initial composition φ 0 , so that Taylor expansion of the reaction mapping of φ around φ 0 is where δφ = φ − φ 0 , A(φ 0 , t) is a n φ × n φ matrix, called mapping gradient matrix, with components given by · denotes the Euclidean norm of a vector, and O( δφ 2 ) denotes terms that have order δφ 2 .
The linear approximation R l (φ, t), obtained by neglecting the high order terms of Eq.(9), is second-order accurate at a connected region of composition space centered at φ 0 . The shape of this region is unknown before the calculations, but ISAT algorithm approximates this region by a hyper-ellipsoid, as will be shown in section 3.2. The local error of this linear approximation is defined as the Euclidean norm of the difference between the reaction mapping at φ and the linear approximation for it around φ 0 ,

Ellipsoid of accuracy
The accuracy of the linear approximation at φ 0 is controlled only if the local error is smaller than a positive error tolerance ε tol , which is heuristically chosen. The region of accuracy is defined as the connected region of the composition space centered at φ 0 where local error is not greater than ε tol . As shown in [24], this region is approximated by a hyper-ellipsoid centered at φ 0 which is dubbed ellipsoid of accuracy (EOA), and is mathematically represented by the following inequation where the EOA Cholesky matrix, denoted by L, is lower triangular [10]. The adaptive step of ISAT algorithm involves the solution of the following geometric problem: given a hyper-ellipsoid centered at φ 0 and a query composition, φ q , outside it, determine a new hyper-ellipsoid of minimum hyper-volume, centered at φ 0 , which encloses both the original hyper-ellipsoid and the point φ q . The solution of this problem is presented by [25] and is not shown here for sake of brevity.

Tabulation procedure
Initially ISAT algorithm receives the time step ∆t and the tolerance ε tol . Then, in every time step, ISAT algorithm receives a query composition φ q and returns an approximation for the corresponding reaction mapping R(φ q , t). This approximation is obtained via numerical integration of Eq. (8) or by the linear approximation R l (φ q , t).
During the reactive flow calculation, the computed values are sequentially stored in a table for future use. This process is known as in situ tabulation. The ISAT table, which is created by the tabulation process, includes the initial composition φ 0 , the reaction mapping R(φ 0 , t), and the mapping gradient matrix A(φ 0 , t). Using these elements it is possible to construct the linear approximation. As the calculation proceeds, a new query composition, φ q , is received by ISAT and the table is transversed until a φ 0 is found that is close to φ q . Depending on the accuracy, the linear approximation around φ 0 is returned or the reaction mapping of φ q is obtained by direct integration of Eq. (8).
The ISAT table is a binary search tree, since this data structure allows for searching an information in O(log 2 n tab ) operations, where n tab is the total number entries in the tree, if the tree is balanced [14]. The binary search tree is basically formed by two types of elements, nodes, and leaves. Each leaf of the tree stores the following data: Each node of the binary search tree has an associated cutting plane. This plane is defined by a normal vector and a scalar such that all compositions φ with v T φ > a are located to the "right" of the cutting plane, as sketched in Figure 1. The cutting plane construction defines a search criterion in the binary search tree. If during the calculation a query point φ q is encountered within the region of accuracy, i.e., ε ≤ ε tol , but outside the estimate of EOA, then the EOA growth proceeds as detailed in [25]. The first three items stored in the binary search tree leaf φ 0 , R(φ 0 , t) and A(φ 0 , t) are computed once, whereas L changes whenever the EOA is grown.
Once a query composition φ q is received by ISAT table, the binary search tree is initialized as a single leaf and the exact value of the reaction mapping is returned.
The subsequent steps are: 1. Given a query composition the tree is transversed until a leaf (φ 0 ) is found. (12) is used to determine if φ q is inside EOA or not. 3. If φ q is inside EOA, the reaction mapping is given by the linear approximation. This is the first of four outcomes, called retrieve. 4. If φ q is outside EOA, direct integration is used to compute the reaction mapping, and the local error is measured by Eq.(11). 5. If the local error is smaller than tolerance, ε tol , the EOA is grown according to the procedure presented in [25] and the reaction mapping is returned. This outcome is called growth. 6. If local error is greater than the tolerance ε tol and the maximum number of entries in the binary search tree is not reached, a new record is stored in the binary search tree based on φ q and the reaction mapping is returned. The original leaf is replaced by a node with the left leaf representing the old composition φ 0 and the right leaf the new one φ q as shown in Figure 2. This outcome is an addition. 7. If the local error is greater than the tolerance ε tol and the maximum number of entries in the binary search tree is reacted, the reaction mapping is returned. This outcome is called direct evaluation. Figure 2 Binary search tree before and after the addition of a new node.

Results and discussion
This section presents the results of numerical simulations performed in order to assess characteristics of PMSR reactor. This study investigates three configurations of a PMSR reactor, two mixtures of CO/O 2 , and one mixture of CH 4 /air, using different time scales to define their behavior. Accuracy, performance, and memory usage of ISAT technique implementation are also investigated. For this purpose, the approximation for solution of the model equation, Eq. (2), obtained by ISAT technique are compared with a reference solution from direct integration (DI) procedure, described in section 2.3.

CO/O 2 mixture
The first two cases studied consider a 1024 particles PMSR, initially filled with a fuel-lean mixture of CO/O 2 at 2948.5 K and 1 atm. At every time step, a fuel-lean mixture of CO/O 2 enters the reactor at T in = 300 K and p in = 1 atm. Both fuel-leam mixtures have equivalence ratio equal to 0.7. The constant pressure and enthalpy equilibrium state associated to the inflow mixture is reached at T eq = 2948.5 K. The reaction mechanism used to describe CO with O 2 kinetics involves 4 species and 3 reactions, [9]. For the first case a time scales configuration with τ m /τ r = τ p /τ r = 1/2 is used, so that the pairwise/mixing time scales are of the same order of magnitude as the residence time. In this situation, partially stirred reactor (PaSR) conditions result within the reactor. For the second case, time scales configuration adopted is τ m /τ r = τ p /τ r = 1/10, so that the pairwise/mixing time scales are small when compared to the residence time. Thus, the reactor should behave almost as a perfect stirred reactor (PSR), where the processes of mixing and pairing occur instantaneously. These cases use a binary search tree with a maximum of 50, 000 entries, ∆t = 10 µs, and ε tol = 10 −3 .

CH 4 /air mixture
The third case studied consists of a PMSR, with 1024 particles, initially filled with the combustion products of a stoichiometric mixture of CH 4 /air at 2100 K and 1 atm. At every time step, a stoichiometric mixture of CH 4 /air enters the reactor at T in = 300 K and p in = 1 atm. The constant pressure and enthalpy equilibrium state associated to this mixture is reached at T eq = 2225.5 K. The reaction of CH 4 with air kinetics is described by GRI mechanism version 3.0 [30], with 53 species and 325 reactions. For this PMSR one have τ m /τ r = τ p /τ r = 1/4, so that it is expected to behave like a PaSR. This simulation uses a binary search tree with a maximum of 60,000 entries, ∆t = 0.1 ms and ε tol = 10 −3 .

Analysis of reactor behavior
and the ensemble average and the ensemble variance operators are respectively defined as being ψ a generic property of the reactive system. For case 1 results, which span over a range of 500 residence times, one can observe an excellent qualitative agreement for T * and T 2 * . The ensemble average value rapidly drops from the initial value to reach the statistically steady-state regime around T * = 0.3.
The ensemble variance rapidly grows, then decreases until it reaches the statistically steady-state value around T 2 * = 0.13. The analysis of the two figures shows that the statistically steady-state regime is reached after τ * = 10.
In case 2, where a range of 250 residence times have been computed, one can also observe an excellent qualitative agreement for T * and T 2 * . Again, the overall history of the PMSR is the same for DI and ISAT. Similar results, not shown here for sake of brevity, have been obtained for the other thermochemical properties of the reactors. The difference among both cases is due to the behavior of each reactor at the statistically steady state     regime. Indeed, the case 2 reactor behavior is governed by a competition between the chemical and residence times mostly. Therefore, the thermodynamical properties steady-state PDF is concentrated over a smaller range than in case 1, where the mixing and pairing time scales are large when compared to the residence time. This behavior is illustrated in Figures 5 and 6, which present the comparison between DI and ISAT computations of mean histograms, averaged over the last 50 residence times, of T * and Y O for cases 1 and 2. These figures underscore the influence of PMSR controlling parameters, i.e., the time scales ratios τ m /τ r and τ p /τ r , on the thermochemical conditions prevailing within each reactor. Indeed, the temperature within the case 2 reactor is such that almost only burned gases are found. On the other hand, case 1 reactor is characterized by a bimodal temperature distribution with a large probability of finding T * = 0.04 and a broader temperature distribution leaning toward the burned gases. Such a distribution, as it could be expected, is reflected on Y O histogram, which also exhibits a bimodal distribution.
The comparison between DI and ISAT results for the first two statistical moments of T * and Y OH , for case 3, are presented in Figures 7 and 8. One can note reasonable and good agreements for the T * and the Y OH statistics, respectively. The results of case 3 show a large discrepancy for the T * than that obtained in case 1. The reaction mechanism of methane is much more complex than the mechanism used to model the carbon monoxide chemical kinetics in case 1, which would in principle lead to the PMSR with methane to assume a wider range of possible thermodynamic states. Thus, it could be expected that the present binary search tree, with 60k entries (almost similar to that used in case 1), would yield comparatively lower accuracy. Finally, Figure 9 presents the comparison between DI and ISAT computations of the mean histograms, averaged over the last 50 residence times, of the T * and Y HCO for case 3. The temperature histogram presents a bimodal distribution with large probability of finding T * ≥ 0.9 and a broader temperature distribution leaning to the fresh gases. On the other hand, Y HCO histogram shows a distribution essentially concentrated in the fresh gases region, and nearly homogeneous elsewhere. This behavior illustrates the fact that HCO, an       intermediate species, appears in low concentration in the burned gases.

Analysis of ISAT accuracy
The relative local error of ψ, which is a time dependent function, is defined as     where subscripts DI and ISAT denote DI and ISAT calculations of ψ, respectively. This metric represents a local measure of error incurred by ISAT.
The evolution of the relative local errors of T * and Y O statistical moments are presented in Figures 10 and 11. Concerning case 1 errors, one can observe a large statistical variation due to the stochastic nature of the PMSR model, with amplitudes reaching 13%. In case 2, relative errors of the order of 1% only can be observed.
All realizations of the stochastic processes in DI and ISAT calculations shown above use the same statistical seed for pseudorandom number generator. However, since the mixing model is statistical in nature, it could be expected that the seed value may influence the ISAT behavior. Therefore, if the DI calculation seed is kept fixed and ISAT seed is changed, the results for T * and ε r , T * are expected to be modified. This can be seen in Figure 12, where one may observe, by comparison with Figures 3 and 10, an increase in the ε r , T * for both cases.
The analysis of the graphs in Figure 12 indicates that 1024 particles are not sufficient to guarantee the statistical independence of the computed results. This hypothesis is also confirmed if one considers Figure 13, where it is possible to see discrepancies in DI and ISAT mean histograms of T * for both cases, even if, qualitatively these histograms are not much different than those shown in Figure 5. Thus, one can conclude that although the results not exhibit the statistical independence with a sample of 1024 particles, they do not vary much with a sample of this size. Possibly a sample of 4096 particles is sufficient to ensure the independence of the results.
The mean relative error of ψ over an interval ∆τ , defined as is a metric that gives a global measure of the error incurred by ISAT. The mean relative errors associated to the statistical moments of the T * and the Y O are presented in Table 1. One can observe that, when the same seed is used for ISAT and DI, the mean error of the average properties is smaller than 1.2% in case 1, and of 0.1% in case 2 only. For T * and Y O ensemble variances, the mean error is large in case 1. Concerning the simulations where the used seeds are different, the mean error of the properties average/variance are an order of magnitude greater than the previous one. In the early development of ISAT technique [24] it was noted that the choice of the tolerance could affect the accuracy of the problem solution. In order to investigate the effect of the tolerance on the present results, Figure 14 presents the absolute global error, ε g , as a function of ISAT error tolerance for cases 1 and 2. The absolute global error is defined, over a time interval ∆τ , as where φ denotes a vector which the components are the ensemble averages of φ components. As may be seen in Figure 14, when the same statistical seed is used, ε g decreases linearly as ε tol is reduced, for both studied cases. However, if different statistical seeds are used, one can note a limit behavior where ε g does not decrease if ε tol is reduced. This saturation in ε g value (∼ 10 −2 ) indicates that a decrease in ε tol value below 10 −3 is not effective in cases where the statistical seeds are different.
The evolution of the relative local error of the statistical moments of T * and Y OH , for case 3, are presented in Figure 15. Compared with the results of case 1, the errors of case 3 present more oscillations. Furthermore, for this case, ε g = 2.4 × 10 −3 is an order of magnitude larger than the values obtained in cases 1 and 2, also using ε tol = 10 −3 . Since different seeds are used for DI and ISAT calculations in case 3, large errors in the results are to be expected.

Analysis of ISAT performance
The comparison of time evolution of ISAT outputs (number of additions, growths, retrieves, and direct evaluations) and of binary search tree height, as well as the corresponding rates of change, for cases 1 and 2, are presented in Figures 16 and 17. A first important observation is that the number of additions in both cases reaches the maximum prescribed value in the binary search tree of 50k. As a consequence of binary search tree saturation, the additions curve reaches a steady state after 5.4 and 1.2 residence times (τ r ), in the first and second cases, respectively. This indicates that ISAT table is saturated earlier when mixing is faster (τ m /τ r small). Figures 16 and 17 also show the evolution of binary tree height, which reaches steady state 1.8 τ r in the first case and 0.4 τ r in the second case. It is also noteworthy that, in both cases, the tree height is an order of magnitude smaller than the total number of entries in the tree (∼ 17k in case 1 and ∼ 7k in case 2). This difference between height and total entries in data table ensures the efficiency of search process for a new query, which could be performed up to three and seven times faster in cases 1 and 2, respectively, than a vector search.
In the first case, Figures 16 and 17 show that the number of growths presents a sharp rate of change around τ r , whereas, in the second case, this occurs around 5 τ r . In both cases, growth steady state occurs after 10 τ r . During both simulations the number of growths is always smaller than additions. This indicates that the desirable massive increase of the ellipsoids of accuracy, to improve the estimate of the region of accuracy, is not observed. This behavior might circumstantial to the reaction system of the carbon monoxide considered, since, due to its simplicity (only three reactions), a small part of the realizable region should be assessed during the course of the calculation.
The Figure 16 shows that, after tree saturation occurs, the number of retrieves, and direct evaluations exceed the number of additions in both cases. In case 1 there is a higher occurrence of retrieve events, whereas in case 2 direct evaluations prevail. The number of retrieves exhibits a linear limit behavior in both cases. The ISAT behavior for the second case reflects the fact     that the binary tree of this case is poor, i.e., contains too few compositions in the region accessed by the calculation. As a consequence, the number of direct evaluations vastly outnumbers ISAT operations. This explains the error behavior observed in section 4.4.
When ISAT technique is more efficient than DI procedure, the computational time spent by ISAT is smaller than computational time spent by DI. The computational time spent by ISAT is the sum of the computa-tional time spent at each of its possible outputs. Therefore, the efficiency condition can be stated as where n A , n G , n R , n DE , n DI are the number of additions, growths, retrieves, direct evaluations, and direct integrations, respectively, whereas τ A , τ G , τ R , τ DE , and      τ DI are the corresponding average duration at each operation.
Assuming that it is possible to show [5] that a necessary (but not sufficient) condition for a calculation using ISAT to be faster than the same calculation using DI is For CO/O 2 mixtures estimates lead to τ A /τ DI ≈ 10, [5]. It is possible to see in Figure 16 that, for case 1, n R approximately exceeds n A by a factor of 130, whereas in      case 2 this factor is only 6. Therefore, as τ A /τ DI < 10 in case 2, ISAT calculations are not expected to be faster than DI procedure.
As can be seen in Table 2, where a comparison of computational time is shown, cases 1 and 2, for ε tol = 10 −3 are computed using DI in 4.0 ks and 2.0 ks respectively, whereas with the use of ISAT, the same cases spent 2.3 ks and 2.1 ks, respectively. Speed-up factors of 1.7 (case 1) and 1.0 (case 2) are obtained, where the speed-up factor is defined as the ratio between the computational time spent by DI and the computational time spent by ISAT. This table also allows to compare the computational time spent by DI and ISAT for different values of error tolerance. An increase in processing time is obtained as ISAT error tolerance is reduced, which is to be expected, given the fact that lower values of ε tol correspond to a smaller region of accuracy. Indeed, as ε tol is decreased, it is less likely that ISAT returns a retrieve, which is ISAT output with lower computational cost. In case 1, for all values of ε tol , ISAT technique offers an advantage, in terms of computational time, when compared to the process of di-rect integration, reducing on average the computational time in 42%. On the other hand, in case 2, no reduction in computational time is seen. As discussed above, this behavior is natural, once the necessary condition for efficiency of the algorithm is not reached.
Aiming to analyze the asymptotic behavior of ISAT speed-up, case 1 has also been simulated for 50,000 time steps. This computation with ISAT spends 204.7 ks, whereas if DI were to be used the same calculation would spend ∼ 1024.0 ks (speed-up factor of 5.0). Therefore, in this asymptotic case, ISAT spends approximately 80% less time than DI. The pioneer work of Pope (1997) [24] reports an asymptotic speed-up factor of 1000, but this impressive factor has not been observed in the study developed here.
The evolution of ISAT algorithm outputs, ISAT binary search tree height, and the corresponding rates of change, for case 3, are presented in Figure 18. As in cases 1 and 2 the number of additions reaches the maximum allowed value in the binary search tree, 60k for this case. The additions curve reaches a steady state after τ * = 3.2. The binary search tree maximum height for case 3 is ∼ 7k as in case 2. Here the steady state of the tree height is observed at τ * = 1.28, as can be seen in Figure 18. The behavior of the number of growths is quite similar to that of case 1, where the greater rate of change occurs near τ * = 2 but, now, the steady state is reached before τ * = 10.
For case 3, nevertheless, the number of growths is not always smaller than the number of additions. Initially the growths exceed additions by nearly an order of magnitude. After τ * = 0.2 the additions exceeded the number of growths, only to be overcome again once steady state is reached. This behavior indicates that the ellipsoids of accuracy show a considerable period of adaptation, which allows to access, with precision, a larger portion of the realizable region.
As in the previous cases, after tree saturation, the number of retrieves and direct evaluations exceed the number of additions. The number of direct evaluations overcomes total number of retrieves in approximately 60%, which is not negligible, but less than that observed in case 2. This indicates that the binary search tree covers a significant portion of the realizable region in the composition space. From Figure 18 one can estimate τ A /τ DI ≈ 50 in order that ISAT technique be more efficient than DI.
For case 3, the computational time spent by DI and ISAT are 689 ks and 455 ks, respectively. Therefore, a speed-up of 1.5 is observed. In this case, the effect of other values of ISAT error tolerance has not been investigated. This PMSR has 1024 particles and its evolution is computed during τ * = 150. Thus, the system of governing equations defined by Eq.(8) is solved 7, 680, 000 times. In this case ISAT technique allows to save 34%, in terms of computational time. For problems that require solving Eq.(8) several times, one could speculate that ISAT technique would provide an even better performance improvement, since more retrieves are expected when compared with the more costly direct integration.

Analysis of ISAT memory usage
The studied cases 1 and 2 are both modeled by a reaction mechanism with 4 species and use a binary search tree with 50,000 entries for ISAT simulations. These parameters lead to a memory consumption by ISAT algorithm of approximately 40 Mbytes, which is very small when compared with the available memory in the workstation used.
On the other hand, case 3 is modeled by a reaction mechanism with 53 species and use a binary search tree with 60,000 entries. This case uses approximately 3.3 Gbytes of memory. It is noteworthy that this value of memory usage is not negligible, when compared to the total memory available on the workstation. Thus, for practical purposes, this tree is the largest one that may be used to simulate a PMSR with this methane reaction mechanism and 1024 particles. This underscores what is perhaps the greatest shortcoming of ISAT technique, its huge expense of memory.

Concluding remarks
This work assessed characteristics of a PMSR reactor through the numerical simulation of reactor configurations with simple (CO/O 2 ) and complex (CH 4 /air) reaction mechanisms and also investigated ISAT technique as an option to evaluate a computational model with detailed thermochemistry. For this purpose, the accuracy, performance, and memory usage of the corresponding technique were analyzed through the comparison of ISAT results with a reference result, obtained from direct integration of PMSR model equation. The studied cases analysis showed that ISAT technique has a good accuracy from a global point of view. Also, it was possible to identify the statistical seed effect on the control of absolute global error, which decreases monotonically as ISAT error tolerance was reduced, when the same seed is used for ISAT and DI calculations. On the other hand, when different seeds were used, a limit value for ISAT error tolerance was observed. In terms of performance, ISAT technique allows to reduce the com-   putational time of the simulations in all cases studied. For CO/O 2 cases, a speed-up of up to 5.0 was achieved, whereas for CH 4 /air, the algorithm allowed to save 34% in terms of computational time. Moreover, ISAT technique presented the desired feature of speed-up factor increase with the complexity of analyzed system. Regarding the memory usage, ISAT technique showed to be very demanding. This work will be extended by coupling ISAT technique to the hybrid LES/PDF model by [2,34,1,33] for description of turbulent combustion, when a detailed reaction mechanism would allow a better description of combustion. In this context, ISAT could be a viable option that may be able to decrease, to an acceptable level, the simulation time.