Interface learning of multiphysics and multiscale systems
Shady E. Ahmed, Omer San, Kursat Kara, Rami Younis, Adil Rasheed
IInterface learning of multiphysics and multiscale systems
Shady E. Ahmed, Omer San, ∗ and Kursat Kara School of Mechanical & Aerospace Engineering, Oklahoma State University, Stillwater, OK 74078, USA.
Rami Younis
The McDougall School of Petroleum Engineering,The University of Tulsa, Tulsa, OK 74104, USA.
Adil Rasheed
Department of Engineering Cybernetics, Norwegian University of Science and Technology, N-7465, Trondheim, Norway. (Dated: June 19, 2020)Complex natural or engineered systems comprise multiple characteristic scales, multiple spa-tiotemporal domains, and even multiple physical closure laws. To address such challenges, weintroduce an interface learning paradigm and put forth a data-driven closure approach based onmemory embedding to provide physically correct boundary conditions at the interface. Our findingsin solving a bi-zonal nonlinear advection-diffusion equation show the proposed approach’s promisefor multiphysics and multiscale systems. We also highlight that high-performance computing envi-ronments can benefit from this methodology to reduce communication costs among processing unitsin emerging machine learning ready heterogeneous platforms toward exascale era.
Numerical simulations of multiscale, multicompo-nent, multiphysics and multidisciplinary systems requireproper treatment of interface boundary conditions amongsolvers. Otherwise, in a naive implementation, the stiffestpart of the domain dictates the overall spatial mesh reso-lution and time stepping requirements, making such sim-ulations computationally daunting. Most of such prob-lems that incorporate some sort of information share canbe put into the following six categories, explained withexamples as follows:
1. FOM-ROM coupling:
With the emergence of digi-tal twin like technologies, there is a demand for lightermodels that can run in real time [1, 2]. In the contextof weather prediction, the full order model (FOM) hasbeen in use for a long time; however, they are incapableof modelling phenomena associated with scales smallerthan what the coarse mesh can handle (like buildings andsmall terrain variations). These fine scale flow structurescan be modelled using a much refined mesh but thenthe simulations become computationally intractable. Totackle this problem a large variety of reduced order mod-els (ROMs) are being developed. In order to make theseROMs realistic there is a need to couple them to FOMmodel so that the FOM provides the interface conditions(both in space and time) to the ROM.
2. Multiphysics and multiscale coupling:
Various flowdichotomies with a multiphysics coupling of interactingsubsystems can be identified in a gas turbine flow. Therotating parts and wall turbulence largely govern the flowwithin the compressor and turbine sections. On the otherhand, for the flow within the combustor, chemical reac-tions, heat release, acoustics, and the presence of fuelspray come into play. Thus, a simulation of the flow ∗ [email protected] within the combustion chamber is significantly more ex-pensive and demanding than other sections. It would re-quire a finer and more sophisticated mesh, smaller timestep, and less numerical simplifications. Therefore, us-ing a unified global solver for the whole system would beeither too expensive (matching the level of the fidelityrequired for the expensive part), or unacceptably inaccu-rate (following the level of fidelity required for the inex-pensive part). Instead, multiple solvers are usually uti-lized to address different parts [3, 4], and information istransferred between solvers.
3. Geometric multiscale:
Blood flow in the whole circu-latory system is mathematically described by means ofheterogeneous problems featuring different degree of de-tail and different geometric dimension that interact to-gether through appropriate interface coupling conditions.Proper exchange of interface conditions between modelsoperating at different geometric approximations open al-together new vistas for biofluids simulations [5–7].
4. Model fusion:
Turbulence modelling generally requiresan apriori selection of the most suited model to handlea particular kind of flow. However, it is seldom thatone model is sufficient for different kind of zones in thecomputational domain. To alleviate this problem, hy-brid and blending models have been extensively utilizedto lift technical barriers in industrial applications, espe-cially in settings where the Reynolds-averaged Navier-Stokes (RANS) approach is not sufficient and large eddysimulation (LES) is expensive [8–10]. The approach canbe extended to blend any number of turbulence modelsprovided the exchange of information at the interface canbe accurately modelled [11].
5. Nested solvers:
To decrease the computational costrequired for an accurate representation of the numerousinterconnected physical systems, e.g., oceanic and atmo-spheric flows, several classes of nested models have beendeveloped and form the basis of highly successful appli- a r X i v : . [ phy s i c s . c o m p - ph ] J un Multiphysics& multiscaleGeometricmultiscale
INTERFACELEARNING
ROM-FOMcoupling DomaindecompositionNestedsolversModel fusion
Toward non-iterative domain decompositionFocus on coherent dynamics and patternsExchange physics at the interface Classification and blendingAn informed extrapolation Effective cyberinfrastructure
Generate a reduced order model (ROM)representation for zones with coherentstructures to supplement full order model (FOM)simulations, e.g., active flow control applicationsto delay separationIdentify the interactions betweenvarious physical phenomena at ascale and mitigate resourcesoptimally considering e.g., fast andslow dynamics, microscale andmacroscale descriptions, constitutivemodels, fluid-structure interactions Establish a mechanism to provide accurateprolongation operations from a low fidelity(averaged) space to a high fidelity space, e.g.,coupling between 0D lumped, 1D distributed, andfull 3D models in cardiovascular system modeling Utilize archival space-time backgroundinformation for reducing communicationcost, relaxing synchronization, andimproving load balancingDevelop a hierarchical spectralnudging approach for complexdynamical systems and architectmachine learning models to bridgethese nested solvers in anautomated fashionTrain machine learning (i.e., reinforcementlearning) tools for a continuous decision makingtoward scale-aware and physics-specific blendingmechanisms among numerous models, e.g.,smart hybridization criteria for RANS/LES models
FIG. 1. Big picture of the interface learning paradigm considering numerous scientific and engineering interpretations. cations and research at numerous weather and climatecenters. Enforcing consistent flow conditions betweensuccessive nesting levels is also considered one form ofinterface matching. For example, a spectral nudging ap-proach has been successfully implemented to force thelarge-scale atmospheric states from global climate mod-els onto a regional climate model [12–17].
6. Domain decomposition:
Since various zones in mul-tiscale systems are connected through interfaces, datasharing, and communicating consistent interface bound-ary conditions among respective solvers are inevitable.Likewise, multirate and locally adaptive stepping meth-ods can yield a mismatch at the space-time interface,and simple interpolation or extrapolation might lead tosolution divergence or instabilities [18]. An analogoussituation usually occurs in parallel computing environ-ments with domain decomposition and distribution overseparate processors with message passing interface tocommunicate information between processors. The het-erogeneity of different processing units creates an asyn-chronous computational environment, and the slowestprocessors will control the computational speed unlessload-balancing is performed [19, 20].In short we can conclude that developing novelmethodologies to model the information exchange at theinterface will have far reaching impact on a large varietyof problems as shown in Fig. 1. To this end, the currentletter puts forth an approach based on memory embed-ding via machine learning to provide physically correctinterfacial conditions. In particular, the proposed tech-nique relies on the time history of local information to es-timate consistent boundary conditions at the sub-domainboundaries without the need to resolve the neighboringregions (on the other side of the interface). It enables us to focus our computational resources on the regionor scales of interest. The proof-of-concept computationson a bi-zonal one-dimensional Burgers’ problem show-case the proposed approach’s promise for stiff multiscalesystems. We also highlight that high-performance com-puting environments can benefit from this methodologyto reduce communication costs among processing units.
Two-component system — As a demonstration of in-terface learning, we consider an application to the one-dimensional (1D) viscous Burgers problem. It combinesthe effects of viscous diffusion, friction, and nonlinearadvection, and thus serves as a prototypical test bed forseveral numerical simulations studies. In order to mimicmultiscale/multiphysics systems, we suppose the domainconsists of two distinguishable zones corresponding to dif-ferent physical parameters as follows, ∂u∂t + u ∂u∂x = ν ∂ u∂x − γu, (1)( ν, γ ) = (cid:40) (10 − ,
0) for 0 ≤ x ≤ x b (10 − ,
1) for x b < x ≤ , (2)where x b is the spatial location of the interface. In anaive implementation, a numerical solution of this prob-lem would imply the use of a grid resolution and timestep corresponding to the stiffest part (i.e., the left partin this case), all over the domain unless we adopt an im-plicit scheme which is unconditionally stable but requiresa nonlinear solver typically at each time step. Certainly,this puts an excessive and unnecessary computationalburden. For example, if we opt to using a spatial res-olution of 4096 grid spacings with a simple forward intime central in space (FTCS) scheme, the maximum timestep that can be used in the left zone is approximately2 . × − (i.e., δt ≤ δx / (2 ν ) based on von Neumannstability analysis). Instead, the right zone gives the flex-ibility of using a 100 times larger time step. However,resolving the whole domain simultaneously would dictatethe smaller time step, even if we are only interested in theright zone. A similar scenario would take place in multi-component systems with varying spatial grid resolutions,where a unified resolution all over the domain becomesunpractical. Thus, we explore the introduction of a mem-ory embedding architecture to enable resolving the zoneof interest independently of the rest of the domain. Memory embedding of interface boundaries — For ma-chine learning applicability, a pattern must exist andmost fluid flows are dominated by coherent structures.Thus, our underlying hypothesis is the existence of a dy-namical context or correlation between the past values offlow features at the interface and its one-sided neighbors(i.e., u ( x b , t n ), u ( x b + δx, t n ), u ( x b + 2 δx, t n ), . . . ), andthe future state at the interface (i.e., u ( x b , t n + δt )). Thiscorresponds to the Learn from Past (LP) model in Fig. 2.Since we incorporate a fully explicit time stepping schemein our simulation, the interface neighboring points mightbe evolved in time before the interface condition is up-dated. Thus, a variant of the LP model based on a combi-nation between old and updated values, namely the
Learnfrom Past and Present (LPP) model, can be utilized aswell. Furthermore, we extend this mapping to take intoaccount the time history dependence in a non-Markovianmanner through the adoption of recurrent neural net-works. Those exploit an internal state feature that re-serves information from past input to learn the context to improve and refine the output. For the neural networkarchitecture, we use a simplistic long short-term mem-ory (LSTM) of two layers, 20 neurons each. Althoughmore sophisticated ML architectures and/or numericalschemes (e.g., [21–24]) might be utilized, the main objec-tive of the present study is to emphasize the potential ofneural networks to advance computational fluid dynam-ics (CFD) simulations for multiscale and multicomponentsystems.
Proof-of-concept results — For the demonstration andassessment of the introduced methodology of interfacelearning, we consider two examples of varying complexity.In the first example, we address the problem of a trav-elling square wave, where the initial condition is definedwith an amplitude of 1 in the left zone (i.e., 0 ≤ x ≤ x b ),and zero in the right zone. In other words, the interface isexactly placed at the discontinuity location of the initialpropagating wave. So, the wave is guaranteed to instan-taneously enter the right zone once the flow is triggered.In a second example of increasing complexity, we studythe evolution of a pulse wave completely contained in aportion of the left region. The initiation of flow dynam-ics in the truncated domain is controlled by the interplaybetween advection, diffusion, and friction in different re-gions. For both examples, we solve the presented viscous1D Burgers problem for a time span of [0 ,
1] using a timestep of 2 . × − to resolve the whole domain x ∈ [0 , LP Model interfaceleft boundary right boundaryLPP Model From LSTMFrom macro-solver Flow variable at Flow variable at
FIG. 2. Different models to utilize LSTM mapping for learn-ing boundary conditions at interface over a spatial grid resolution of 4096. For external bound-ary conditions, we assume zero Dirichlet boundary con-ditions (i.e., u (0 , t ) = u (1 , t ) = 0). Data snapshots arestored every 100 time steps (corresponding to the coarsetime step of 2 . × − ).For interface learning, we consider two schemes/modelsfor the training as illustrated in Fig. 2 to learn the dy-namics at the internal boundary separating the two com-partments of different physical properties. For testing,we consider the truncated 1D domain, where x b ≤ x ≤ . × − , thus denoted the macro-solver here. We adopt the LSTM learning to update theleft boundary condition (i.e., at x = x b ). For the rightboundary (i.e., at x = 1), we keep the global zero Dirich-let conditions. Example 1: travelling square wave . In this example,we consider the following initial condition, u ( x,
0) = (cid:40) ≤ x ≤ x b x b < x ≤ . (3)In particular, we generate data for x b ∈{ / , / , / , / , / , / , / } , and we use fielddata at x b ∈ { / , / , / , / } for training and reservethe remaining cases for the out-of-sample testing. To en-hance the neural network performance, we augment theinput vector with the spatial and temporal informationas well. In other words, LP Model can be interpretedas the mapping u ( x b , t n + δt ) = G (cid:0) u ( x b , t n ) , u ( x b + δx, t n ) , u ( x b + 2 δx, t n ) , . . . , x b , x b + δx, x b + 2 δx, . . . , t n (cid:1) , x u ( x , t ) timetimetimetimetime x b = 0.25 x u ( x , t ) timetimetimetimetime x b = 0.50 x u ( x , t ) timetimetimetimetime x b = 0.75 x u ( x , t ) timetimetimetimetime x u ( x , t ) timetimetimetimetimeTrue LSTM BC Closure x u ( x , t ) timetimetimetimetime LP ModelLPP Model
FIG. 3. Results for LSTM boundary condition closure for different values of x b . Predicted velocity fields are shown at t ∈ { . , . , . , . , . } . while LPP Model learns the map of u ( x b , t n + δt ) = G (cid:0) u ( x b , t n ) , u ( x b + δx, t n ) , u ( x b + 2 δx, t n ) , . . . , u ( x b + δx, t n + δt ) , u ( x b + 2 δx, t n + δt ) , . . . , x b , x b + δx, x b +2 δx, . . . , t n , t n + δt (cid:1) .We compare the predicted velocity field within thetruncated domain using the proposed LSTM boundarycondition (BC) closure approach with respect to the truesolution obtained by solving the whole domain. Wenote here that the LSTM BC closure results are basedon utilizing the macro-solver (i.e., using a time step of2 . × − ), while the true solution is obtained by adopt-ing the micro-solver (i.e., using a time step of 2 . × − ).The spatio-temporal evolution of the velocity field for x b ∈ { / , / , / } is shown in Fig. 3, where x b is thelocation of the interface. We note that Fig. 3 corre-sponds to a three-point stencil for the LSTM mapping.In other words, the LP model uses values of u ( x b , t n ), u ( x b + δx, t n ), and u ( x b + 2 δx, t n ) for the predictionof u ( x b , t n + δt ), while the LPP model uses u ( x b , t n ), u ( x b + δx, t n ), u ( x b + 2 δx, t n ), u ( x b + δx, t n + δt ), and u ( x b +2 δx, t n + δt ). Visual results advocate the capabilityof the presented approach of predicting accurate valuesfor the interface boundary condition at different times.For more quantitative assessment, we also compute theresulting 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 (cid:0) u T ( x i , t n ) − u P ( x i , t n ) (cid:1) , (4) where u T is the true velocity field, and u P represents thepredictions by the LSTM BC closure approach. In theabove formula, N x stands for the number of grid pointsinvolved only in the truncated domain. In other words,it considers only the flow field values within [ x b , old values of the flow fields.Moreover, this might be attributed to the sub-optimalarchitecture we use for the LSTM. Although we foundthat results are not very sensitive to the given hyper-parameters , further tuning might be required to provide optimal performance. We also see from Table I that anincrease in the stencil size from 2 to 3 improves results.Nonetheless, a 2-point stencil mapping still provides ac-ceptable predictions, confirming the validity and robust-ness of the LSTM memory embedding skills to yield phys-ically consistent and accurate state estimates at the in-terface using local information, and may hold immensepotential for designing ML-ready predictive engines inphysical sciences. Example 2: pulse problem . In this case, we consideran initial condition of a pulse, completely contained inthe left region and study its propagation and travel fromthe left to right compartments. In particular, the initial
TABLE I. RMSE of LSTM boundary condition closure re-sults using different models with two-point and three-pointmapping. LP Model LPP Model x b .
125 1 . × − . × − . × − . × − .
250 4 . × − . × − . × − . × − .
375 8 . × − . × − . × − . × − .
500 2 . × − . × − . × − . × − .
625 6 . × − . × − . × − . × − .
750 1 . × − . × − . × − . × − .
875 2 . × − . × − . × − . × − pulse can be represented as u ( x,
0) = (cid:40) ≤ x ≤ w p w p < x ≤ , (5)where w p is the pulse width. For illustration, westore results corresponding to 7 varying pulse widthsas w p ∈ { . , . , . , . , . , . , . } . We use w p ∈ { . , . , . , . } for training and validation,while we assign w p ∈ { . , . , . } for out-of-sampletesting. For interface, we consider a fixed interface loca-tion at x b = 0 . before the pulse travelsinto the truncated zone. Thus the state at the interfaceis more dependent on the flow dynamics in both domainpartitions. Since the pulse width is a key factor in thisproblem setting, we augment our input vector with w p as well. For this particular example, we found that en-forcing higher memory embedding is crucial in providingaccurate results. Specifically, we adopt a sliding windowof a three-time step history (also called a lookback of 3)in our LSTM implementation. Results for the LP andLPP schemes are shown in Fig. 4 using 3-point mapping.We find that both schemes can sufficiently learn the in-terface dynamics and accurately predict its condition atout-of-sample settings.In Table II, we report the RMSE obtained from the in-terface learning approach using the LP an LPP schemesfor the 1D pulse problem described above. We can seethat memory embedding for interface learning yields veryaccurate results for truncated domain computations. Re-sults in this table include those obtained from both in-sample and out-of-sample data. However, we highlighthere that even the in-sample training and testing aresubstantially different. During the training phase, theLSTM is always fed with true inputs obtained from themicro-solver applied over the whole domain. In contrast,for testing, the LSTM only sees the true state at initialtime. After that, the output of the LSTM is provided tothe macro-solver to evolve in time, and the evolved statesare then given to the LSTM as new input. This recursivedeployment can introduce numerical noise as a result of the temporal advancement with larger time-step as wellas uncertainties at interface conditions. TABLE II. RMSE of LSTM boundary condition closure forthe pulse example using different models with two-point andthree-point stencils.LP Model LPP Model w p .
20 2 . × − . × − . × − . × − .
21 3 . × − . × − . × − . × − .
22 2 . × − . × − . × − . × − .
23 1 . × − . × − . × − . × − .
24 7 . × − . × − . × − . × − .
25 8 . × − . × − . × − . × − .
26 1 . × − . × − . × − . × − Conclusions — In this work, we have demonstratedthe potential of machine learning tools to advance andfacilitate CFD simulations of multiscale, multicompo-nent systems. In particular, we have shown the ca-pability of memory embedding to learn the dynamicsat the interface between different zones. This is espe-cially beneficial where the domain contains zones withstrong dynamics and components with complex config-uration that might dictate a very fine mesh resolutionsand time stepping. The proposed approach enables usto focus our efforts onto the domain portion of inter-est, while satisfying physically consistent interface con-ditions. It can serve as a non-iterative domain decompo-sition method. Toward model fusion technologies, suchan interface learning methodology might also hold signif-icant promise for the development of blending criteria inhybrid RANS/LES models. A proof-of-concept has beendemonstrated using the 1D viscous Burgers equation overa two-zone domain with different physical parameters.An LSTM is used to bypass the micro-solver correspond-ing the stiff region and provide valid interface boundaryconditions to enable the macro-solver to run indepen-dently. We have illustrated the success and robustness ofthe proposed methodology using different learning config-urations. Both LP and LLP models are the key conceptsintroduced in this letter, especially for designing intel-ligent boundary closure schemes, which may bear hugepotential in many scientific disciplines. Finally, we em-phasize that a similar interface closure technique can beadopted in high performance computing environments,to minimize the communication cost and delay betweendifferent asynchronous processors, a topic we would liketo pursue further in the future.This material is based upon work supported by theU.S. Department of Energy, Office of Science, Office ofAdvanced Scientific Computing Research under AwardNumber de-sc0019290. O.S. gratefully acknowledgestheir support.The data that support the findings of this study areavailable from the corresponding author upon request. x u ( x , t ) timetimetimetimetime w p = 0.21 x u ( x , t ) timetimetimetimetime w p = 0.23 x u ( x , t ) timetimetimetimetime w p = 0.25 x u ( x , t ) timetimetimetimetime x u ( x , t ) timetimetimetimetimeTrue LSTM BC Closure x u ( x , t ) timetimetimetimetime LP ModelLPP Model
FIG. 4. Results for LSTM boundary condition closure for the pulse problem using different values of w p . Predicted velocityfields are shown at t ∈ { . , . , . , . , . } .[1] B. Peherstorfer, K. Willcox, and M. Gunzburger, SIAMReview , 550 (2018).[2] F. Naets, D. De Gregoriis, and W. Desmet, InternationalJournal for Numerical Methods in Engineering , 765(2019).[3] S. Shankaran, J. Alonso, M.-F. Liou, N.-S. Liu, andR. Davis, in (2001) p. 974.[4] J. Xu, C. Huang, and K. Duraisamy, AIAA Journal ,618 (2020).[5] A. Quarteroni and A. Veneziani, Multiscale Modeling &Simulation , 173 (2003).[6] T. Passerini, M. De Luca, L. Formaggia, A. Quarteroni,and A. Veneziani, Journal of Engineering Mathematics , 319 (2009).[7] A. Quarteroni, A. Veneziani, and C. Vergara, ComputerMethods in Applied Mechanics and Engineering , 193(2016).[8] P. Sagaut, Multiscale and multiresolution approaches inturbulence: LES, DES and hybrid RANS/LES methods:applications and guidelines (Imperial College Press, Dan-vers, MA, USA, 2013).[9] A. Fadai-Ghotbi, C. Friess, R. Manceau, and J. Bor´ee,Physics of Fluids , 055104 (2010).[10] M. L. Shur, P. R. Spalart, M. K. Strelets, and A. K.Travin, International Journal of Heat and Fluid Flow ,1638 (2008).[11] R. Maulik, O. San, J. D. Jacob, and C. Crick, Journalof Fluid Mechanics , 784 (2019).[12] K. M. Waldron, J. Paegle, and J. D. Horel, MonthlyWeather Review , 529 (1996). [13] H. von Storch, H. Langenberg, and F. Feser, MonthlyWeather Review , 3664 (2000).[14] R. Radu, M. D´equ´e, and S. Somot, Tellus A: DynamicMeteorology and Oceanography , 898 (2008).[15] G. Miguez-Macho, G. L. Stenchikov, and A. Robock,Journal of Geophysical Research: Atmospheres (2004).[16] B. Rockel, C. L. Castro, R. A. Pielke Sr, H. von Storch,and G. Leoncini, Journal of Geophysical Research: At-mospheres (2008).[17] M. Schubert-Frisius, F. Feser, H. von Storch, andS. Rast, Monthly Weather Review , 909 (2017).[18] M. J. Gander and L. Halpern, in Domain DecompositionMethods in Science and Engineering XX (Springer, 2013)pp. 377–385.[19] D. A. Donzis and K. Aditya, Journal of ComputationalPhysics , 370 (2014).[20] A. Mittal and S. Girimaji, Physical Review E , 033304(2017).[21] A. Rasheed, O. San, and T. Kvamsdal, IEEE Access ,21980 (2020).[22] A. Bowers, L. Rebholz, A. Takhirov, and C. Trenchea,International Journal for Numerical Methods in Fluids , 805 (2012).[23] C. C. Manica, M. Neda, M. Olshanskii, and L. G. Reb-holz, ESAIM: Mathematical Modelling and NumericalAnalysis , 277 (2011).[24] J. Borggaard and T. Iliescu, Applied Mathematics Let-ters19