gLaSDI: Parametric Physics-informed Greedy Latent Space Dynamics Identification

Physical simulations have played an increasingly significant role in developments of engineering, science, and technology. The widespread applications of physical simulations in digital twins systems [1, 2] is one recent example. Many physical processes are mathematically modeled by time-dependent nonlinear partial differential equations (PDEs). As it is difficult or even impossible to obtain analytical solutions for many highly complicated problems, various numerical methods have been developed to approximate the analytical solutions. However, due to the complexity and the domain size of problems, high-fidelity forward physical simulations can be computationally intractable even with high performance computing, which prohibits their applications to problems that require a large number of forward simulations, such as design optimization [3, 4], optimal control [5], uncertainty quantification [6, 7], and inverse analysis [7, 8].

To achieve accurate and efficient physical simulations, data can play a key role. For example, various physics-constrained data-driven model reduction techniques have been developed, such as the projection-based reduced-order model (ROM), in which the state fields of the full-order model (FOM) are projected to a linear or nonlinear subspace so that the dimension of the state fields is significantly reduced. Popular linear projection techniques include the proper orthogonal decomposition (POD) [9], the reduced basis method [10], and the balanced truncation method [11], while autoencoders [12, 13] are often applied for nonlinear projection [14, 15, 16]. The linear-subspace ROM (LS-ROM) has been successfully applied to various problems, such as nonlinear heat conduction [17], Lagrangian hydrodynamics [18, 19, 20], nonlinear diffusion equations [17, 21], Burgers equations [20, 22, 23, 24], convection-diffusion equations [25, 26], Navier-Stokes equations [27, 28], Boltzmann transport problems [29, 30], fracture mechanics [31, 32], molecular dynamics [33, 34], fatigue analysis under cycling-induced plastic deformations [35], topology optimization [36, 37], structural design optimization [38, 39], etc. Despite successes of the classical LS-ROM in many applications, it is limited to the assumption that intrinsic solution space falls into a low-dimensional subspace, which means the solution space has a small Kolmogorov $n$-width. This assumption is not satisfied in advection-dominated systems with sharp gradients, moving shock fronts, and turbulence, which prohibits the applications of the LS-ROM approaches for these systems. On the other hand, it has been shown that nonlinear-subspace ROMs based on autoencoders outperforms the LS-ROM on advection-dominated systems [14, 40].

Several strategies have been developed to extend LS-ROM for addressing the challenge posed by advection-dominated systems, which can be mainly categorized into Lagrangian-based approaches and methods based on a transport-invariant coordinate frame. In the first category of approaches, Lagrangian coordinate grids are leveraged to build a ROM that propagates both the wave physics and the coordinate grid in time [41, 42, 43]. Although these methods work well, their applicability is limited by the requirement of full knowledge of the governing equations for obtaining the Lagrangian grid. The second category of strategies is based on transforming the system dynamics to a moving coordinate frame by adding a time-dependent shift to the spatial coordinates such that the system dynamics are absent of advection, such as the shifted POD method [44] and the implicit feature tracking algorithm based on a minimal-residual ROM [45]. Despite the effectiveness of these methods, the high computational costs prohibit their applications in practice.

Most aforementioned physics-constrained data-driven projection-based ROMs are intrusive, which require plugging the reduced-order solution representation into the discretized system of governing equations. Although the intrusive brings many benefits, such as extrapolation robustness, a requirement of less training data, and high accuracy, the implementation of the intrusive ROMs requires not only sufficient understanding of the numerical solver of the high-fidelity simulation but also access to the source code of the numerical solver.

In contrast, non-intrusive ROMs are purely data-driven. It requires neither access to the source code nor the knowledge of the high-fidelity solver. Many non-intrusive ROMs are constructed based on interpolation techniques that provide nonlinear mapping to relate inputs to outputs. Among various interpolation techniques, such as Gaussian processes [46, 47], radial basis functions [48, 49], Kriging [50, 51], neural networks (NNs) have been most popular due to their strong flexibility and capability supported by the universal approximation theorem [52]. NN-based surrogates have been applied to various physical simulations, such as fluid dynamics [53], particle simulations [54], bioinformatics [55], deep Koopman dynamical models [56], porous media flow [57, 58, 59, 60], etc. However, pure black-box NN-based surrogates lack interpretability and suffer from unstable and inaccurate generalization performance. For example, Swischuk, et al. [61] compared various ML models, including NNs, multivariate polynomial regression, $k$-nearest neighbors (k-NNs), and decision trees, used to learn nonlinear mapping between input parameters and low-dimensional representations of solution fields obtained by POD projection. The numerical examples of this study shows that the highly flexible NN-based model performs worst, which highlights the importance of choosing an appropriate ML strategy by considering the bias-variance trade off, especially when the training data coverage of the input space is sparse.

In recent years, several ROM methods have been integrated with latent-space learning algorithms. Kim, et al. [62] proposed a DeepFluids framework in which the autoencoder was applied for nonlinear projection and a latent-space time integrator was used to approximate the evolution of the solutions in the latent space. Xie, et al. [63] applied the POD for linear projection and a multi-step NN to propagate the latent-space dynamical solutions. Hoang, et al. [64] applied the POD to compress space-time solution space to obtain space-time reduced-order basis and examined several surrogate models to map input parameters to space-time basis coefficients, including multivariate polynomial regression, k-NNs, random forest, and NNs. Kadeethum, et al. [59] compared performance of the POD and autoencoder compression along with various latent space interpolation techniques, such as radial basis function and artificial neural networks. However, the latent-space dynamics models of these methods are complex and lack interpretability.

To improve the interpretability and generalization capability, it is critical to identify the underlying equations governing the latent-space dynamics. Many methods have been developed for the identification of interpretable governing laws from data, including symbolic regression that searches both parameters and the governing equations simultaneously [65, 66], parametric models that fit parameters to equations of a given form, such as the sparse identification of nonlinear dynamics (SINDy) [67], and operator inference [68, 69, 70]. Cranmer et al. [71] applied graph neural networks to learn sparse latent representations and symbolic regression with a genetic algorithm to discover explicit analytical relations of the learned latent representations, which enhances efficiency of symbolic regression to high-dimensional data [72]. Instead of genetic algorithms, other techniques have been applied to guide the equation search in symbolic regression, such as gradient descent [73, 74] and Monte Carlo Tree Search with asymptotic constraints in NNs [75].

Champion, et al. [76] applied an autoencoder for nonlinear projection and SINDy to identify simple ordinary differential equations (ODEs) that govern the latent-space dynamics. The autoencoder and the SINDy model were trained interactively to achieve simple latent-space dynamics. However, the proposed SINDy-autoencoder method is not parameterized and generalizable. Bai and Peng [77] proposed parametric non-intrusive ROMs that combine the POD for linear projection and regression surrogates to approximate dynamical systems of latent variables, including support vector machines with kernel functions, tree-based methods, k-NNs, vectorial kernel orthogonal greedy algorithm (VKOGA), and SINDy. The ROMs integrated with VKOGA and SINDy deliver superior cost versus error trade-off. Additionally, various non-intrusive ROMs have been developed based on POD-based linear projection with latent space dynamics captured by polynomials through operator inference [68, 82, 83, 84, 85, 69, 86, 70, 87, 78, 79, 80, 81]. For example, Qian, et al. [69] introduced a lifting map to transform non-polynomial physical dynamics to quadratic polynomial dynamics and then combined POD-based linear projection with operator inference to identify quadratic reduced models for dynamical systems. Due to the limitation of the POD-based linear projection, these non-intrusive ROMs have difficulties with advection-dominated problems.

To address this challenge, Issan and Kramer [88] recently proposed a non-intrusive ROM based on shifted operator inference by transforming the original coordinate frame of dynamical systems to a moving coordinate frame in which the dynamics are absent of translation and rotation. Fries, et al. [40] proposed a parametric latent space dynamics identification (LaSDI) framework in which an autoencoder was applied for nonlinear projection and a set of parametric dynamics identification (DI) models were introduced in the parameter space to identify local latent-space dynamics, as illustrated in Fig. 1. The LaSDI framework can be viewed as a generalization of aforementioned non-intrusive ROMs built upon latent-space dynamics identification, since it allows linear or nonlinear projection and enables latent-space dynamics to be captured by flexible DI models based on general nonlinear functions. However, since a sequential training procedure was adopted for the autoencoder and the DI models, the lack of interaction between the autoencoder and the DI models leads to a strong dependency on the complexity and quality of the latent-space dynamics on the autoencoder architecture, which could pose challenges to the subsequent training of the DI models and thus affect the model performances. Most importantly, all the above-mentioned approaches rely on predefined training samples, such as uniform or Latin hypercube sampling that may not be optimal in terms of the number of samples for achieving the best model performance in the prescribed parameter space. As the generation of the simulation data can be computationally expensive, it is important to minimize the number of samples.

In this study, we propose a parametric adaptive greedy latent space dynamics identification (gLaSDI) framework for accurate, efficient, and robust physics-informed data-driven reduced-order modeling. To maximize and accelerate the exploration of the parameter space for optimal performance, an adaptive greedy sampling algorithm integrated with a physics-informed residual-based error indicator and random-subset evaluation is introduced to search for the optimal and minimal training samples on the fly. The proposed gLaSDI framework contains an autoencoder for nonlinear projection to discover intrinsic latent representations and a set of local DI models to capture local latent-space dynamics in the parameter space, which is further exploited by an efficient k-NN convex interpolation scheme. The concept of gLaSDI is similar to the active learning algorithms that are allowed to choose training data from which it learns [89, 90]. The DI models in gLaSDI can be viewed as the active learners in active learning, which are selected to maximize “diversity” of learners in the parameter space and minimize prediction errors. The autoencoder training and dynamics identification in the gLaSDI take place interactively to achieve an optimal identification of simple latent-space dynamics. The effectiveness and enhanced performance of the proposed gLaSDI framework is demonstrated by modeling various nonlinear dynamical problems with a comparison with the LaSDI framework [40].

The remainder of this paper is organized as follows. The governing equations of dynamical systems is introduced in Section 2. In Section 3, the ingredients of the proposed gLaSDI framework, the mathematical formulations, the training and testing algorithms are introduced. In Section 4, the effectiveness and capability of the proposed gLaSDI framework are examined by modeling various nonlinear dynamical problems, including Burgers equations, nonlinear heat conduction, and radial advection. The effects of various factors on model performance are investigated, including the number of nearest neighbors for convex interpolation, the latent-space dimension, the complexity of the DI models, and the size of the parameter space. A performance comparison between uniform sampling, i.e., LaSDI, and the physics-informed greedy sampling, i.e., gLaSDI, is also presented. Concluding remarks and discussions are summarized in Section 5.

A parameterized dynamical system characterized by a system of ordinary differential equations (ODEs) is considered

where $T\in\mathbb{R}_{+}$ is the final time; $\boldsymbol{\mu}\in\mathcal{D}\subseteq\mathbb{R}^{n_{\mu}}$ is the parameter in a parameter domain $\mathcal{D}$. $\mathbf{u}(t;\boldsymbol{\mu}):[0,T]\times\mathcal{D}\rightarrow\mathbb{R}^{N_{u}}$ is the parameterized time-dependent solution to the dynamical system; $\mathbf{f}:\mathbb{R}^{N_{u}}\times[0,T]\times\mathcal{D}\rightarrow\mathbb{R}^{N_{u}}$ denotes the velocity of $\mathbf{u}$; $\mathbf{u}_{0}$ is the initial state of $\mathbf{u}$. Eq. (1) can be considered as a semi-discretized equation of a system of partial differential equations (PDEs) with a spatial domain $\Omega\subseteq\mathbb{R}^{d},d\in\mathbb{N}(3)$, and $\mathbb{N}(N):=\{1,...,N\}$. Spatial discretization can be performed by numerical methods, such as the finite element method.

A uniform time discretization is considered in this study, with a time step size $\Delta t\in\mathbb{R}_{+}$ and $t_{n}=t_{n-1}+\Delta t$ for $n\in\mathbb{N}(N_{t})$ where $t_{0}=0$, $N_{t}\in\mathbb{N}$. Various explicit or implicit time integration schemes can be applied to solve Eq. (1). For example, with the implicit backward Euler time integrator, the approximate solutions to Eq. (1) can be obtained by solving the following nonlinear system of equations

where $\mathbf{u}_{n}:=\mathbf{u}(t^{n};\boldsymbol{\mu})$, and $\mathbf{f}_{n}:=\mathbf{f}(\mathbf{u}(t^{n};\boldsymbol{\mu}),t^{n};\boldsymbol{\mu})$. The residual function of Eq. (2) is expressed as

Solving the dynamical system of equation (Eq. (1)) could be computationally expensive, especially when the solution dimension ($N_{u}$) is large and the computational domain ($\Omega$) is geometrically complex. In this work, an efficient and accurate data-driven reduced-order modeling framework based on physics-informed greedy latent space dynamics identification is proposed, which will be discussed in details in the next section.

The ingredients of the proposed physics-informed parametric adaptive greedy latent space dynamics identification (gLaSDI) framework are introduced in this section, including autoencoders, dynamics identification (DI) models, $k$-nearest neighbors (k-NN) convex interpolation, and an adaptive greedy sampling procedure with a physics-informed error indicator. An autoencoder is trained to discover an intrinsic nonlinear latent representation of high-dimensional data from dynamical PDEs, while training cases are sampled on the fly and associated local DI models are trained simultaneously to capture localized latent-space dynamics. An interactive training procedure is employed for the autoencoder and local DI models, which is referred to as the interactive Autoencoder-DI training, as shown in Fig. 2. The interaction between the autoencoder and local (in the parameter space) DI models enables the identification of simple and smooth latent-space dynamics and therefore accurate and efficient data-driven reduced-order modeling. In the following, Sections 3.1 - 3.3 review the basics of autoencoders, the dynamics identification method, and the convex interpolation scheme, respectively, which are fundamental to the physics-informed adaptive greedy sampling algorithm introduced in Section 3.4. Section 3.5 summarizes the algorithms for gLaSDI training and testing.

The physics-informed adaptive greedy sampling can be performed on either a continuous parameter space ($\mathcal{D}$) or a discrete parameter space ($\mathcal{D}^{h}\subseteq\mathcal{D}$). In the following demonstration, a discrete parameter space ($\mathcal{D}^{h}$) is considered. $\mathbb{D}\subseteq\mathcal{D}^{h}$ denotes a set of $N_{\mu}$ selected training sample points.

Let us consider a training sample point $\boldsymbol{\mu}^{(i)}\in\mathcal{D}^{h}$, $i\in\mathbb{N}(N_{\mu})$ and $\mathbf{u}_{n}^{(i)}\in\mathbb{R}^{N_{u}}$ as the solution at the $n$-th time step of the dynamical system in Eq. (1) with the training sample point $\boldsymbol{\mu}^{(i)}$. The solutions at all time steps are arranged in a snapshot matrix denoted as $\mathbf{U}^{(i)}=[\mathbf{u}_{0}^{(i)},...,\mathbf{u}_{N_{t}}^{(i)}]\in\mathbb{R}^{N_{u}\times(N_{t}+1)}$. Concatenating snapshot matrices corresponding to all training sample points gives a full snapshot matrix $\mathbf{U}\in\mathbb{R}^{N_{u}\times(N_{t}+1)N_{\mu}}$

The performance of gLaSDI is demonstrated by solving four numerical problems: one-dimensional (1D) Burgers’ equation, two-dimensional (2D) Burgers’ equation, nonlinear time-dependent heat conduction, and radial advection. The effects of the number of nearest neighbors $k$ on the model performance are discussed. In each of the numerical examples, the gLaSDI’s performance is compared with that of LaSDI, that is the one without adaptive greedy sampling. Note that the physics-informed adaptive greedy sampling can be performed on either a continuous parameter space ($\mathcal{D}$) or a discrete parameter space ($\mathcal{D}^{h}\subseteq\mathcal{D}$). In the following examples, a discrete parameter space ($\mathcal{D}^{h}$) is considered. For all numerical experiments, we always start with the initial database $\mathcal{DB}_{0}$ and the associated parameter set $\mathbb{D}_{0}$ of four corner points of the parameter space.

The training of gLaSDI are performed on a NVIDIA V100 (Volta) GPU of the Livermore Computing Lassen system at the Lawrence Livermore National Laboratory, with 3,168 NVIDIA CUDA Cores and 64 GB GDDR5 GPU Memory. The open-source TensorFlow library [99] and the Adam optimizer [100] are employed for model training. The testing of gLaSDI and high-fidelity simulations are both performed on an IBM Power9 CPU with 128 cores and 3.5 GHz.