Interface learning in fluid dynamics: statistical inference of closures within micro-macro coupling models
IInterface learning in fluid dynamics: statistical inference ofclosures within micro-macro coupling models
Suraj Pawar, Shady E. Ahmed, and Omer San a) School of Mechanical & Aerospace Engineering, Oklahoma State University, Stillwater,OK 74078, USA. (Dated: 12 August 2020)
Many complex multiphysics systems in fluid dynamics involve using solvers withvaried levels of approximations in different regions of the computational domainto resolve multiple spatiotemporal scales present in the flow. The accuracy of thesolution is governed by how the information is exchanged between these solvers atthe interface and several methods have been devised for such coupling problems.In this article, we construct a data-driven model by spatially coupling a microscalelattice Boltzmann method (LBM) solver and macroscale finite difference method(FDM) solver for reaction-diffusion systems. The coupling between the micro-macrosolvers has one to many mapping at the interface leading to the interface closureproblem, and we propose a statistical inference method based on neural networksto learn this closure relation. The performance of the proposed framework in abifidelity setting partitioned between the FDM and LBM domain shows its promisefor complex systems where analytical relations between micro-macro solvers are notavailable.Keywords: Interface closure, Machine learning, Micro-macro coupling, Multiscalesystems
Introduction — Many problems in science and engineering are multiscale in nature con-sisting of different spatial and temporal scales that interact with each other . For example,simulation of turbulent and separated flows requires resolving a wide range of interactingscales . The traditional numerical setup usually focuses on the macroscale evolution ofthe system and either eliminates meso/microscale processes or consider them through someform of closure models . For many complex systems, modeling coarse-scale, macroscopicbehavior might be impractical without simplified assumptions, and one might have to resortto modeling fine-scale, microscale processes to get accurate solution . However, using themicroscale solver over the entire computational domain is computationally prohibitive evenin the exascale computing era. This difficulty can be overcome by utilizing a micro-macrocoupled solver to recover accurate solutions in a computationally efficient manner.There have been several studies that deal with coupling micro-macro solvers, such asfinite difference method-lattice Boltzmann method (FDM-LBM) solver , finite vol-ume method-molecular dynamics (FVM-MD) solver , finite volume method-lattice Boltz-mann method (FVM-LBM) solver , molecular dynamics-lattice Boltzmann (MD-LBM)solver . These approaches can also be interpreted in a broad sense as domain decompo-sition methods, where two solvers simultaneously advance the multiscale and multiphysicsproblems and the information is exchanged across the interface between solvers . One ofthe major challenges in these micro-macro solvers is a mismatch in the kind and number ofvariables used by different models. For example, microscopic solvers such as LBM describethe evolution of the system through particle distribution functions that are restricted tomove on a grid with certain velocities only. On the other hand, the macroscale solverslike partial differential equations (PDEs) model the system in terms of observables likeflow velocity, pressure, and density. The coupling between two solvers should not pro-duce any discontinuity in the interface region. Also, the global quality of the numericalsolution will be dependent on the treatment of interface boundary condition. Therefore,the interface boundary closure is a major challenge in modeling and computation of these a) Electronic mail: [email protected] a r X i v : . [ phy s i c s . c o m p - ph ] A ug micro-macro solvers. Some of the methods to determine interface boundary closure arebased on physical quantities of another solver or from an analytical relation between microand macro solvers . However, the derivation of such analytical expressions for complexgeometries or systems with complex local interactions might be impractical and for thesesystems iterative numerical algorithms have been proposed to close interface boundaryconditions .The challenges in the treatment of interface boundary conditions in complex systemsoffer opportunities for researchers to develop statistical inference approaches for couplingmicro-macro solvers. In fact, several data-driven approaches have been proposed in theliterature to learn the correlation between coarse-scale PDEs and fine-scale microscopicprocesses . The machine learning (ML) algorithms are also proven to be successful indiscovering hidden PDEs from macroscopic observations data . The recent advance-ment in the efficient implementation of ML algorithms along with the huge amount ofarchival data generated from high-fidelity numerical simulations and experiments has madestatistical inference approaches an attractive choice for multiscale systems. To this end,we propose a data-driven ML-based interface closure framework for the spatial couplingof microscale LBM and macroscale FDM solvers. Our framework is solver agnostic in thesense that it can be extended to any type of micro-macro coupled models. We illustrate ourframework for a reaction-diffusion system (the FitzHugh-Nagumo model) and also compareour results against two methods for solving interface closure problems based on numericalapproximations of the analytical relations. FitzHugh-Nagumo model — We illustrate the statistical interface closure for micro-macrocoupled solvers using the FitzHugh-Nagumo model. This model consists of two reactive-diffusion partial differential equations, whose dynamics is governed by ∂u∂t = D u ∂ u∂x + u − u − v, (1) ∂v∂t = D v ∂ v∂x + (cid:15) ( u − a v − a ) , (2)where D u and D v are the diffusion coefficients of u and v , a and a are model pa-rameters, and (cid:15) represents a kinetic bifurcation parameter. We set the parameters to a = 2 , a = − . , (cid:15) = 0 . , D u = 1 , and D v = 4 as suggested by Theodoropoulos, Qian,and Kevrekidis . Our spatial domain extends between [0 ,
20] and the model is integratedfrom time t = 0 to t = 450.We employ FDM as a macro solver and LBM as the micro solver. For the FDM solver,Equation 1 and Equation 2 are discretized with explicit forward in time and the centraldifference in space as follow u ( x, t + ∆ t ) = u ( x, t ) + D u ∆ t ∆ x ( u ( x + ∆ x, t ) − u ( x, t ) + u ( x − ∆ x, t )) + ∆ t ( u − u − v ) , (3) v ( x, t + ∆ t ) = v ( x, t ) + D v ∆ t ∆ x ( v ( x + ∆ x, t ) − v ( x, t ) + v ( x − ∆ x, t )) + (cid:15) ∆ t ( u − a v − a ) , (4)where ∆ x and ∆ t are the spatial and temporal discretization steps, respectively. Since theabove discretization is explicit in time, the time step in the above scheme is restricted by thestability condition ∆ t < C ∆ x . In the FitzHugh-Nagumo model, D v > D u , and thereforethe stability condition will be governed by the diffusion constant of the v equation, i.e.,∆ t < ∆ x / (2 D v ).The LBM method used as a micro solver describes the evolution of particle distributionfunctions f i ( x, t ) discretized in space x , and time t along the i th direction with velocity c i . We utilize the D1Q3 model, where D1 denotes the one-dimensional domain, and Q3stands for three velocities. The lattice Boltzmann equation (LBE) describing the evolutionof particle distribution is given as f li ( x + c i ∆ x, t + ∆ t ) = f li ( x, t ) − ω l ( f li ( x, t ) − f l, eq i ( x, t )) + R li ( x, t ) , l ∈ { u, v } (5)for c i = i ∈ {− , , } . The first step of the LBM is the collision step, where this diffusivecollisions are modeled by the BhatnagarGrossKrook (BGK) model as a relaxation to theequilibrium particle distribution f l, eq i ( x, t ) with a relaxation coefficient ω l (denoted by thesecond term in Equation 5), and the reactions are modeled by the term R li ( x, t ) . Thesecond step of the LBM is the propagation of the particle distribution functions to a neigh-boring grid based on the i th direction as shown on the left hand side of Equation 5. Theequilibrium particle distribution function, relaxation coefficients, and the reaction term arecomputed as below f u, eq i ( x, t ) = 13 u ( x, t ) , f v, eq i ( x, t ) = 13 v ( x, t ) , (6) ω u = 21 + 3 D u ∆ t ∆ x , ω v = 21 + 3 D v ∆ t ∆ x , (7) R ui ( x, t ) = ( u − u − v )3 , R vi ( x, t ) = (cid:15) ( u − a v − a )3 . (8)The observable u and v are defined as the zeroth order moment of the distribution functions u = i =1 (cid:88) i = − f ui ( x, t ) , v = i =1 (cid:88) i = − f vi ( x, t ) . (9) Coupling between FDM and LBM — In order demonstrate the coupling betweenmacroscale FDM solver and microscale LBM solver, we divide our computational do-main into two non-overlapping subdomains as shown in Figure 1. The left domain is solvedusing FDM and the right domain is resolved with LBM. The interface lie at the centerof the domain, i.e., at L/
2, where L is the length of the whole domain. Therefore, thelast grid point of the FDM is located at L/ − ∆ x/ L/ x/
2. For our analysis, we assume that the spatial discretization ∆ x and timestep ∆ t for both subdomains are the same. When the hybrid solver is used to model thedynamics of reactive-diffusion system, the information has to be exchanged carefully acrossthe interface as the FDM and LBM use different sets of variables, namely ( u, v ) versus( f ui , f vi ). In order to solve the FDM equation at the last point of the left domain, weneed information of ( u, v ) at x FDM + ∆ x location which can be obtained using the relationprovided in Equation 9 as follow u ( x FDM + ∆ x, t ) = i =1 (cid:88) i = − f ui ( x LBM , t ) , (10) v ( x FDM + ∆ x, t ) = i =1 (cid:88) i = − f vi ( x LBM , t ) . (11)The inverse mappings from FDM to LBM, i.e., u ( x FDM , t ) → f ui ( x FDM , t ) and v ( x FDM , t ) → f vi ( x FDM , t ) for i ∈ {− , , } are not straightforward. We need f l ( x FDM , t ) to evolve parti-cle distribution at x LBM from time t to t + ∆ t using Equation 5. The particle distributionin LBM can be reduced in terms of the FDM variables using a multiscale ChapmanEnskogexpansion as follow f li ( x, t ) = f l, [0] i ( x, t ) + f l, [1] i ( x, t ) + f l, [2] i ( x, t ) + . . . , (12) FDM LBMMacroscale PDE Mesoscale
FIG. 1. The exchange of information across the interface between two domains. The left domainis solved using FDM and the right domain is solved with LBM and the information is exchangedwith an appropriate boundary condition. where the superscript inside a square bracket denotes the order of approximations.We can reconstruct f l ( x FDM , t ) using the zeroth order ChapmanEnskog relations(CE-0)as shown below f u ( x FDM , t ) = f u, [0]1 ( x FDM , t ) = u ( x FDM , t )3 , (13) f v ( x FDM , t ) = f v, [0]1 ( x FDM , t ) = v ( x FDM , t )3 . (14)The above relation is the simplest approximation one can use and the relation given inEquation 12 describes the particle distribution in LBM up to the second order. Once theparticle distribution function is reconstructed, it has to be propagated from time t to t + ∆ t in two steps. First, the collision and reaction go from time t to t ∗ , and then the post-collision distribution f l, ∗ i is propagated from time t ∗ to t + ∆ t . The post-collision particledistributions are computed as follow f u, ∗ ( x FDM , t ) = (1 − ω u ) f u ( x FDM , t )+ ω u u ( x FDM , t ) + ∆ t R ui ( x, t ) , (15) f v, ∗ ( x FDM , t ) = (1 − ω v ) f v ( x FDM , t )+ ω v v ( x FDM , t ) + ∆ t R vi ( x, t ) . (16)Finally, the post-collision particle distributions are propagated to x LBM as given below f u ( x LBM , t + ∆ t ) = f u, ∗ ( x FDM , t ) , (17) f v ( x LBM , t + ∆ t ) = f v, ∗ ( x FDM , t ) (18)The analytical relation provided in Equation 12 might not always be available or might bedifficult to derive numerically due to complex geometry or complex local interaction forceterms. There are iterative numerical algorithms like constrained runs (CR) scheme thatcan be applied to approximate particle distribution numerically. The correlation betweenmicroscale and macroscale variables can also be learned with data-driven methods and weexplore the feasibility of feedforward neural network to approximate these relations. Wetrain the neural network to learn the particle distribution f l ( x LBM , t ) at time t based on thelocal information at the previous time step. More concretely, our neural network is trainedto learn the below mapping { u ( x FDM , t − ∆ t ) , v ( x FDM , t − ∆ t ) , f u − ( x LBM , t − ∆ t ) ,f v − ( x LBM , t − ∆ t ) } → { f u ( x LBM , t ) , f v ( x LBM , t ) } . (19)The training data for the neural network is generated by simulating the whole computationaldomain using FDM and LBM separately. We employ bounce-back boundary condition onleft and right side of the domain for LBM simulation. The initial condition for the trainingdata is generated with α = 0 . u ( x,
0) = α tanh( x − L/ , (20) v ( x,
0) = 0 . (cid:18) (cid:18) x − L/ (cid:19)(cid:19) . (21)During testing, we use α = 1 . Results and discussion — Here, we present the results of coupled FDM-LBM solver for theFitzHugh-Nagumo model and compare it against the LBM simulation. The computationaldomain is discretized with a spatial grid of size ∆ x = 0 . t = 0 .
001 from time t = 0 to t = 450. The results are provided for α = 1 . u as given in Equation 20. Figure 2 shows the plot of u and v for the LBM simulation, the FDM-LBM solver with CE-0 coupling at the interface,and the difference between two simulations. It can be seen that the CE-0 coupling is notable to predict the accurate dynamics after time t ≈
80. The local error for CE-0 couplingat the interface depends upon the dominant term left out in the ChapmanEnskog expansiongiven in Equation 12. Therefore, the CE-0 introduces an error which is first order in ∆ x .This error seems to have contaminated the solution and the coupled solver based on CE-0leads to the poor prediction of the u and v fields.In Figure 3, we depict the results for coupled FDM-LBM solver using ML-based cou-pling at the interface. The ML-based coupling can predict the correct dynamics of both u and v fields with a substantial improvement compared to the CE-0 coupling. In otherwords, the ML-based interface closure is able to learn the correlation between microscaleand macroscale variables solely from the data and produces more accurate prediction thanthe CE-0 coupling. This finding clearly illustrates the potential of ML for learning interfaceclosure in micro-macro solvers. Although the ML-based closure yields a sufficiently accuratesolution, we reiterate here that for the model problem investigated in this study, there arehigher-order approximations that can be utilized to reconstruct particle distribution func-tions. For example, the first-order approximation of f l ( x FDM , t ) using the ChapmanEnskogrelations (CE-1) is given below f u ( x FDM , t ) = u ( x FDM , t )3 − ∆ x ω u u ( x LBM , t ) − u ( x FDM − ∆ x, t )2∆ x , (22) f v ( x FDM , t ) = v ( x FDM , t )3 − ∆ x ω v v ( x LBM , t ) − v ( x FDM − ∆ x, t )2∆ x , (23)where the first term is the zeroth order approximation, and the second term is the firstorder derivative computed using the central difference scheme.For a fair comparison, we illustrate the results of FDM-LBM solver equipped with CE-1coupling at the interface. Figure 4 displays the u and v fields for the LBM simulation,the FDM-LBM solver with CE-1 coupling at the interface, and the difference betweentwo simulations. The agreement between LBM and hybrid solver is excellent and evenbetter than the ML-based closure at the interface. This clearly shows limitations of data-driven methods compared to higher-order numerical approximation of analytical relationsbetween microscale and macroscale solvers. However, obtaining such higher-order analyticalapproximation for complex geometries and complex local interactions might not be possible.Furthermore, for many complex systems, the evolution equations may not be available inclosed forms and for such systems, the ML-based interface closure will be attractive,especially, in the age of data. x u (LBM) 0 200 40005101520 v (LBM)0 200 40005101520 x u (CE-0) 0 200 40005101520 v (CE-0)0 200 400 t x u (Difference) 0 200 400 t v (Difference) − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . FIG. 2. Spatio-temporal prediction of the u and v field of the FitzHugh-Nagumo model by hybridsolver with CE-0 coupling at the interface. We assess the quantitative performance of different coupling algorithms investigated inthis study using the root mean-squares error (RMSE) defined asRMSE = (cid:118)(cid:117)(cid:117)(cid:116) N t N x N t (cid:88) n =1 N x (cid:88) i =1 ( φ T ( x i , t n ) − φ P ( x i , t n )) , (24)where N t is the total number of time snapshots, N x is the total number of grid points, φ is the variable of interest. The superscripts T and P stand for the truth solution (i.e.,LBM simulation) and predicted solution (i.e., FDM-LBM solver simulation). We store450 temporal snapshots between time t = 0 to t = 450 for the calculation of the RMSE.Table I reports the RMSE for different coupling algorithms presented in this study. The x u (LBM) 0 200 40005101520 v (LBM)0 200 40005101520 x u (ML) 0 200 40005101520 v (ML)0 200 400 t x u (Difference) 0 200 400 t v (Difference) − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . FIG. 3. Spatio-temporal prediction of the u and v field of the FitzHugh-Nagumo model by hybridsolver with ML based coupling at the interface. quantitative assessment implies that the error is minimum with CE-1 coupling and the errorfor the ML-based coupling is substantially less than the CE-0 coupling. We believe thatthe error for the ML-based coupling can be further reduced by including more informationin the input space of the neural network or using a deeper neural network architecture.However, our experiments with hyperparameters for the neural network suggest that thedeeper network leads to overfitting for this relatively simple problem and might not lead toan improvement in the prediction. TABLE I. The root mean squared error between the true and predicted solution for u and v fordifferent coupling algorithms.Couplingalgorithm u v CE-0 6 . × − . × − CE-1 9 . × − . × − ML 1 . × − . × − x u (LBM) 0 200 40005101520 v (LBM)0 200 40005101520 x u (CE-1) 0 200 40005101520 v (CE-1)0 200 400 t x u (Difference) 0 200 400 t v (Difference) − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . FIG. 4. Spatio-temporal prediction of the u and v field of the FitzHugh-Nagumo model by hybridsolver with CE-1 coupling at the interface. Conclusion — In this work, we have introduced a framework based on machine learning(ML) to learn the interface closure between microscale lattice Boltzmann method (LBM)and macroscale finite difference method (FDM) solver. This framework offers a great poten-tial for coupled micro-macro solver where there is one to many mapping across the interfaceand the true relation between microscale and macroscale variables is unknown or difficultto derive, as in the case of complex geometries and multiphysics systems. The frameworkis modular enough that it can be implemented to different coupling models such as RANS-LES, FDM-FVM, FVM-LBM, and others. Our experiments with the FitzHugh-Nagumomodel suggest that the ML-based interface closure framework is able to produce sufficientlyaccurate dynamics over a longer time. Also, the interface boundary closure based on thefirst-order approximation of ChapmanEnskog expansion (CE-1) produced even better pre-diction than the ML-based closure. It is no surprise that the interface closure based on thenumerical approximation of the analytical relations can recover the dynamics almost ex-actly. However, for many multiphysics systems like viscoelastic fluids or porous media, suchrelations might not exist and for such systems, we envision that our framework will be moreadvantageous. This is the first step towards building data-driven interface closure modelsfor micro-macro solvers, and our future efforts will be directed to extend this framework tomore complex higher-dimensional multiphysics problems in fluid dynamics.This material is based upon work supported by the U.S. Department of Energy, Officeof Science, Office of Advanced Scientific Computing Research under Award Number DE-SC0019290. O.S. gratefully acknowledges their support.The data that support the findings of this study are available from the correspondingauthor upon request. E. Weinan, B. Engquist, X. Li, W. Ren, and E. Vanden-Eijnden, “Heterogeneous multiscale methods: areview,” Communications in computational physics , 367–450 (2007). P. Koumoutsakos, “Multiscale flow simulations using particles,” Annu. Rev. Fluid Mech. , 457–487(2005). J. Fish,
Multiscale methods: bridging the scales in science and engineering (Oxford University Press onDemand, 2010). K. Matouˇs, M. G. Geers, V. G. Kouznetsova, and A. Gillman, “A review of predictive nonlinear theoriesfor multiscale modeling of heterogeneous materials,” Journal of Computational Physics , 192–220(2017). P. Sagaut,
Multiscale and multiresolution approaches in turbulence: LES, DES and hybrid RANS/LESmethods: applications and guidelines (World Scientific, 2013). N. G. Van Kampen,
Stochastic processes in physics and chemistry , Vol. 1 (Elsevier, 1992). R. Keunings, “Micro-macro methods for the multiscale simulation of viscoelastic flow using molecularmodels of kinetic theory,” Rheology reviews , 67–98 (2004). S. Attinger, “Generalized coarse graining procedures for flow in porous media,” Computational Geo-sciences , 253–273 (2003). I. G. Kevrekidis and G. Samaey, “Equation-free multiscale computation: Algorithms and applications,”Annual review of physical chemistry (2009). P. Albuquerque, D. Alemani, B. Chopard, and P. Leone, “Coupling a lattice boltzmann and a finitedifference scheme,” in
International Conference on Computational Science (Springer, 2004) pp. 540–547. P. Van Leemput, C. Vandekerckhove, W. Vanroose, and D. Roose, “Accuracy of hybrid lattice boltz-mann/finite difference schemes for reaction-diffusion systems,” Multiscale Modeling & Simulation , 838–857 (2007). X. Nie, S. Chen, M. Robbins, et al. , “A continuum and molecular dynamics hybrid method for micro-andnano-fluid flow,” Journal of Fluid Mechanics , 55 (2004). S. C. Mishra and H. K. Roy, “Solving transient conduction and radiation heat transfer problems usingthe lattice boltzmann method and the finite volume method,” Journal of Computational Physics ,89–107 (2007). H. Joshi, A. Agarwal, B. Puranik, C. Shu, and A. Agrawal, “A hybrid fvm–lbm method for single andmulti-fluid compressible flow problems,” International journal for numerical methods in fluids , 403–427(2010). H.-B. Luan, H. Xu, L. Chen, D.-L. Sun, Y.-L. He, and W.-Q. Tao, “Evaluation of the coupling scheme offvm and lbm for fluid flows around complex geometries,” International journal of heat and mass transfer , 1975–1985 (2011). A. Dupuis, E. Kotsalis, and P. Koumoutsakos, “Coupling lattice boltzmann and molecular dynamicsmodels for dense fluids,” Physical Review E , 046704 (2007). T. F. Chan, R. Glowinski, J. P´eriaux, and O. B. Wildlund, “Domain decomposition methods,” in
Domaindecomposition methods: proceedings of the Second International Symposium on Domain DecompositionMethods, Los Angeles, California, January 14-16, 1988 (1989). H. Tang, R. Haynes, and G. Houzeaux, “A review of domain decomposition methods for simulation offluid flows: Concepts, algorithms, and applications,” Archives of Computational Methods in Engineering, 1–33 (2020). P. Van Leemput, W. Vanroose, and D. Roose, “Numerical and analytical spatial coupling of a latticeboltzmann model and a partial differential equation,” in
Model reduction and coarse-graining approachesfor multiscale phenomena (Springer, 2006) pp. 423–441. C. W. Gear, T. J. Kaper, I. G. Kevrekidis, and A. Zagaris, “Projecting to a slow manifold: Singularlyperturbed systems and legacy codes,” SIAM Journal on Applied Dynamical Systems , 711–732 (2005). S. Sirisup, G. E. Karniadakis, D. Xiu, and I. G. Kevrekidis, “Equation-free/galerkin-free pod-assistedcomputation of incompressible flows,” Journal of Computational Physics , 568–587 (2005). I. G. Kevrekidis, C. W. Gear, J. M. Hyman, P. G. Kevrekidid, O. Runborg, C. Theodoropoulos, et al. ,“Equation-free, coarse-grained multiscale computation: Enabling mocroscopic simulators to performsystem-level analysis,” Communications in Mathematical Sciences , 715–762 (2003). C. Siettos, C. Gear, and I. Kevrekidis, “An equation-free approach to agent-based computation: Bifur-cation analysis and control of stationary states,” EPL (Europhysics Letters) , 48007 (2012). S. Lee, M. Kooshkbaghi, K. Spiliotis, C. I. Siettos, and I. G. Kevrekidis, “Coarse-scale pdes from fine-scale observations via machine learning,” Chaos: An Interdisciplinary Journal of Nonlinear Science ,013141 (2020). M. Raissi and G. E. Karniadakis, “Hidden physics models: Machine learning of nonlinear partial differ-ential equations,” Journal of Computational Physics , 125–141 (2018). S. H. Rudy, S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Data-driven discovery of partial differentialequations,” Science Advances , e1602614 (2017). Z. Long, Y. Lu, X. Ma, and B. Dong, “Pde-net: Learning pdes from data,” in
International Conferenceon Machine Learning (2018) pp. 3208–3216. C. Theodoropoulos, Y.-H. Qian, and I. G. Kevrekidis, “coarse stability and bifurcation analysis usingtime-steppers: A reaction-diffusion example,” Proceedings of the National Academy of Sciences , 9840–9843 (2000). S. Chen and G. D. Doolen, “Lattice boltzmann method for fluid flows,” Annual review of fluid mechanics , 329–364 (1998). S. Ponce Dawson, S. Chen, and G. D. Doolen, “Lattice boltzmann computations for reaction-diffusionequations,” The Journal of chemical physics , 1514–1523 (1993). B. Chopard, A. Dupuis, A. Masselot, and P. Luthi, “Cellular automata and lattice boltzmann techniques:An approach to model and simulate complex systems,” Advances in complex systems , 103–246 (2002). C. W. Gear and I. G. Kevrekidis, “Constraint-defined manifolds: a legacy code approach to low-dimensional computation,” Journal of Scientific Computing25