Computational fluid dynamics (CFD) simulations (and other simulations that involve physicsbased models) often require fine spatial and temporal discretizations to achieve high levels of fidelity. These fine discretizations can make the corresponding simulations expensive to run due to the large number of degrees of freedom that must be modeled. In applications such as design, which can require running many simulations across multiple design configurations, this high computational cost can be prohibitive. Instead of running expensive simulations across the entire design space, a technique known as surrogate modeling is often used [1–3], wherein models are trained to approximate the mapping between design parameters and a scalar-valued objective that evaluates design quality. Surrogate models can be constructed from a limited amount of high-fidelity simulation data using techniques such as polynomial or Gaussian process regression, and allow for efficient evaluation of design objectives at new design points. However, while surrogate modeling approaches can enable faster analysis and design, they do not reduce the computational cost of running high-fidelity simulations.
In place of surrogate models, which only attempt to model the mapping from design parameters to high-level objectives, a great deal of recent attention has focused on using learning-based techniques to facilitate more efficient high-fidelity simulation. Many proposed techniques seek a low-dimensional, reduced representation of the modeled system that can be simulated more efficiently than the full system dynamics. This generally requires first learning a mapping from the high-dimensional dynamical system states to reduced states, and subsequently devising a method for integrating the reduced states forward in time. Projection-based approaches [4, 5] learn a mapping from high-dimensional states to a low-dimensional subspace, and then determine the reduced state dynamics through a projection onto known governing equations. The learned mappings are generally obtained through proper orthogonal decomposition (POD), which yields linear embeddings, although some recent work has explored the use of nonlinear embeddings [6, 7]. One drawback of projection-based approaches is that they must be integrated into existing simulation software, which can be challenging and time intensive [8].
An alternative approach for determining the reduced-state dynamics is to learn the dynamics directly from data. Provided the learned models are accurate, such a data-driven approach requires no integration with existing simulation software and no knowledge of the system governing equations. A number of recent works have applied two-stage (dimensionality reduction and dynamics learning) deep learning procedures to the task of modeling fluid flows [9–15]. These approaches generally require complex multi-stage training procedures, and often can only be applied to simulation at a single flow condition. In contrast, this paper proposes a new method, grounded in variational inference, for parameter-conditioned modeling of fluid flows, in which the state mapping and reduced state dynamics are learned jointly. The derived learning procedure allows the trained model to perform parameterized simulations, thus enabling the efficient simulation of modeled systems at a wide range of conditions. To motivate the proposed method, the discussion now proceeds by reviewing relevant background material related to generative modeling.
Generative modeling refers to modeling the process by which data is generated in the real world [16]. As such, physical models like the Navier–Stokes equations can be viewed as a form of generative model. In the absence of known governing equations, learned generative models can attempt to model real-world systems by approximating the distribution over data, ). Several powerful generative modeling paradigms have emerged in recent years, including variational autoencoders [17], generative adversarial networks [18], and flow-based models [19, 20]. The following section provides an introduction to variational autoencoders.
2.1 Variational Autoencoders
The goal in generative modeling is to learn ), a distribution over (possibly high-dimensional) data
, where
represents the parameters that govern the learned distribution. An accurate learned distribution will be close to the true data distribution,
), where the difference between the distributions can be measured through similarity measures such as the KL-divergence:
It can be shown that minimizing the KL-divergence between these distributions is equivalent to maximizing the log-likelihood the model assigns to the data [16]:
Success in generative modeling relies upon a model that (1) can be optimized efficiently (i.e. the gradients ) are known and easy to calculate) and (2) is sufficiently expressive in order to accurately approximate the true data distribution. For example, a joint Gaussian distribution over pixel color intensities could be used to represent a distribution over image data, but, if the true data distribution is non-Gaussian, then such a model would be incapable of representing
) even under the optimal parameter settings.
The expressiveness of modeled distributions can be enhanced by introducing latent variables, denoted by z. Latent variables are unobserved, but are assumed to summarize high-level information about the data, and are often assumed to be (much) lower-dimensional than the data. As an example, in image data the latent variable z may be taken to represent the class of an image (e.g. dog or human or car), and the data x would then represent a particular instantiation of that class. Under the presence of latent variables, the modeled distribution becomes:
When the distribution is modeled in this manner then, even if ) and
) are Gaussian distributions, the resulting marginal distribution
) can be quite complex.
Unfortunately, calculating (and in turn optimizing) the likelihood as defined in Eq. (3) is generally intractable due to the difficulty of evaluating the integral over z. Additionally, finding the posterior over z given x, defined as:
is also intractable due to the presence of ) in the denominator of the expression. In place of this intractable posterior distribution, variational inference [21] introduces an approximate posterior distribution
), which is defined by a set of parameters
. Multiplying and dividing the likelihood expression by this approximate posterior distribution yields:
Figure 1: Illustration of the training procedure for the variational autoencoder. Red lines indicate relationships relevant to the optimization objective.
Through Jensen’s inequality the following expression can be derived:
The right-hand side of the above expression is often referred to as the evidence lower bound (ELBO), and serves as a lower bound on the log-likelihood objective. It can be shown that the ELBO will be equal to log
i.e. the approximate posterior distribution matches the true posterior distribution over z. Thus, optimizing the evidence lower bound with respect to parameters and
serves the dual purpose of (1) approximately maximizing the true objective log
) and (2) driving the approximate posterior closer to the true posterior distribution [16].
The variational autoencoder (VAE) [17] presents a method for performing variational inference with neural network models. Much like standard autoencoders, VAEs consist of encoder and decoder components. The encoder (or recognition model) is denoted by ), and maps an input x to a distribution over the latent variable z. The decoder (or generative model) is denoted by
), and maps samples from
) to a distribution over x. The outputs of the neural networks are the distribution parameters, which are often assumed to be Gaussian with diagonal or constant covariance matrices.
Training a variational autoencoder consists of optimizing the set of parameters in order to maximize the evidence lower bound presented in Eq. (6). Figure 1 illustrates the VAE training procedure. Training inputs are fed into the encoder, which outputs parameters to a distribution over the latent variable z. Subsequently, samples from that distribution are mapped through the decoder to obtain a distribution over x. The optimization objective balances maximizing the likelihood assigned to x with minimizing the KL-divergence between the approximate posterior and the prior distribution. For simplicity, the prior is often assumed to be a standard Gaussian distribution. Evaluating the gradients with respect to the encoder parameters requires the use of a technique known as the reparameterization trick [22], which treats the encoder as a deterministic mapping subject to exogenous noise, and allows for the propagation of gradients through the sampling operation. A VAE that has been trained to maximize the evidence lower bound can subsequently be employed in a generative fashion,
Figure 2: Graphical models for VAE, CVAE, and sequential generative model.
where samples are drawn directly from the prior ) and passed through the decoder. This process can result in the generation of new data that appears as if it was drawn from
2.2 Conditional Variational Autoencoders
The conditional variational autoencoder (CVAE) [23] extends the VAE by attempting to learn ), a conditional model for the data distribution. Rather than modeling the distribution over the entire dataset, a CVAE models data distributions conditioned on context parameters c. The graphical model for the CVAE can be found in Fig. 2b.
As with the variational autoencoder, a lower bound to the log-likelihood objective can be derived for the CVAE:
Optimizing for this objective requires relatively few modifications to the VAE training procedure outlined in Section 2.1. Most notably, the encoder and decoder networks must be modified to accept the context parameters as inputs. Furthermore, the prior over z must be modeled as a conditional prior distribution, which is typically accomplished by introducing an auxiliary neural network that is trained to take c as an input and output the prior distribution parameters.
2.3 Generative Modeling of Sequential Data
Generative models for sequential data attempt to learn the joint distribution ), where
represent data with temporal relationships; here it is assumed that each
represents the state of a dynamical system at a given time step. In constructing sequential generative models, it is common to assume the existence of (lower-dimensional) latent variables
, also referred to as latent states, that govern the time evolution of the modeled system. Furthermore, the dynamics of the latent states are often assumed to be Markovian, meaning that a latent state
is conditionally independent of
given
(i.e.
summarizes all relevant past information). The graphical model that defines the problem structure can be found in Fig. 2c; a model of this form is commonly referred to as a state space model (SSM) [24].
Analogous to the material presented in Section 2.1, the likelihood expression can be written in a form that explicitly accounts for the presence of the latent variables:
Evaluating this integral is generally intractable due to the need to marginalize over the latent states. Likewise, evaluating the posterior distribution ) is generally intractable, once again motivating the introduction of an approximate posterior distribution
Given the conditional independence assumptions encoded by the graphical model, the conditional distribution ) can be rewritten as:
Likewise, the distribution ) can be simplified to:
Finally, the approximate posterior distribution can be expressed as:
Given these simplified expressions, a lower bound to the log-likelihood, log be derived by following the procedures outlined in Section 2.1. This lower bound is given by:
The ELBO thus consists of T likelihood components, representing the likelihood of states given latent states
KL-divergence components, representing the difference between
Figure 3: Graphical model for parameter-conditioned generative model. Parameters c are assumed to influence the time evolution of the latent states.
the approximate posterior and (conditional) prior distributions. If the latent states are assumed to follow non-Markovian dynamics, then the ELBO becomes:
where the elements that differentiate Eq. (14) from Eq. (13) are highlighted in blue.
As with the variational autoencoder, the prior distribution ) can be defined to be a standard Gaussian. Additionally, the parameters of the conditional prior distributions can be output by auxiliary neural networks. In contrast with the ELBOs for the VAE and CVAE, the ELBOs derived in this section contain distributions that require conditioning on sequences of variables. These distributions can be represented by standard feedforward neural networks, but such architectures are generally incapable of accommodating exceedingly long or variable-length sequences as inputs. This consideration motivates the use of recurrent neural networks, which allow for learning from arbitrarily long input sequences. Common recurrent neural network architectures include long short term memory (LSTM) [25] and gated recurrent unit (GRU) [26] networks.
2.4 Parameter-Conditioned Generative Modeling
This work proposes an extension to generative models for sequential data that is akin to the CVAE, wherein the latent state dynamics are assumed to be governed by parameters c. In the context of dynamical systems, c may represent physical parameters that can affect system behavior. Examples of such parameters include the mass and length of a pendulum, or the Reynolds number of a fluid flow. Thus, a parameter-conditioned generative model attempts to model the density
The graphical model for this problem can be found in Fig. 3. Following the same procedure as in Section 2.3, the evidence lower bound for the parameter-conditioned sequential generative model can be derived. For the case of non-Markovian latent state dynamics, this
lower bound is found to be:
Hence, all prior and approximate posterior distributions must now condition on the parameters c.
A trained parameter-conditioned generative model provides the opportunity to generate sequential data that evolves according to the prescribed parameters c. In the context of dynamical systems, this could enable the efficient study of how such parameters affect a system’s behavior. Given a parameter space of interest, training data can be collected by simulating (or experimentally evaluating) a system at a handful of parameter values. Subsequently, a generative model trained on the collected data could potentially provide insight into the nature of the studied system throughout the remainder of the parameter space. A generative model’s ability to accurately simulate a system throughout the parameter space is partly dependent on how well it learns to account for the manner in which the parameters affect the system’s dynamics. A description of a new objective that encourages the model to learn these relationships can be found in the next subsection.
2.5 Mutual Information Objective
Entropy, H(y), serves as a measure of uncertainty about a random variable y, and is defined as:
Intuitively, a random variable will have high entropy if its probability mass (or density) is spread evenly across its support. The mutual information I(w, y) between random variables w and y is defined as:
Thus, I(w, y) represents the reduction in uncertainty about w given knowledge of y; if w and y are independent, then I(w, y) = 0.
To ensure that parameter-conditioned generative models can be used to perform simulations throughout the parameter space, high mutual information is desired between the prescribed parameters c and the latent states . This can be explicitly enforced by training the generative model to maximize a mutual information objective,
addition to the evidence lower bound from Eq. (15). Unfortunately, maximizing the mutual information objective directly is challenging because it requires access to the generally intractable posterior distributions
In a procedure similar to variational inference, variational mutual information maximization [27–29] introduces a lower bound to the mutual information objective that can be tractably optimized. Let ) represent an approximate posterior distribution over c, which is defined by parameters
. Then the following lower bound to the mutual information objective can be derived:
It is assumed that p(c) represents the distribution over parameter values within the dataset, and thus H(c) is treated as a constant. Hence, raising the lower bound on the mutual information objective requires maximizing the log-likelihood assigned to c by the approximate posterior distribution. This can be accomplished in practice by defining a reconstruction model with parameters , which takes in values of
sampled from the prior distributions and outputs parameters to a distribution over c. All model parameters,
, and
, can then be jointly optimized to maximize the objective:
where represents the lower bound to the likelihood objective,
represents the lower bound to the mutual information objective, and
is a scalar value that trades off between the objectives.
The next section presents experimental results that evaluate the performance of parameter-conditioned generative models on simulating fluid flows.
This section applies parameter-conditioned generative modeling to the task of simulating fluid flows. Evaluations are performed in two problem domains: (1) two-dimensional airflow over a pair of counter-rotating cylinders and (2) three-dimensional airflow over a half-cylinder. The next section studies the ability of learned models to simulate two-dimensional airflow, evaluating the accuracy of generated solutions at a range of flow conditions. A discussion of generative modeling performance for three-dimensional flow follows thereafter.
3.1 Counter-Rotating Cylinders
The flow over non-streamlined, bluff bodies such as cylinders has been the subject of extensive experimental and numerical study. The scientific interest in such flows is motivated by a phenomenon known as vortex-induced vibration, where wake instabilities lead to unsteady vortex shedding over the surface of bluff bodies. Vortex-induced vibrations can arise in many engineered structures, such as buildings, bridges, and offshore piers. Vortex shedding gives rise to strong transverse forces, and is associated with higher drag and unsteady lift forces [30, 31]. The collapse of the Tacoma Narrows Bridge in 1940 serves as a particularly stark illustration of the destructive potential of the induced structural vibrations [32].
This section considers the test case of airflow over counter-rotating cylinders, as depicted in Fig. 4. The cylinder surfaces are assumed to be separated by a distance of 3d, where d is the cylinder diameter, and rotate at a fixed non-dimensional rotation speed given by:
Figure 4: Counter-rotating cylinder configuration.
where is the angular velocity of the cylinder and
is the velocity of the incident flow. As with flow over a single bluff body, flow over a pair of cylinders has been shown to exhibit unsteady vortex shedding [33–35]. Experimental and numerical results have demonstrated that there exists a critical rotation speed, Ω
, above which the wake instabilities disappear and the flow becomes steady, leaving a constant transverse load acting upon the cylinders [36– 38]. The thesis of Chan [32] performs a wide array of numerical experiments that characterize the behavior of the double cylinder system as the cylinder separation, Reynolds number, rotation speed, and rotation direction are varied. The experiments in this work only allow for variation in rotation speed and Reynolds number (i.e.
); the cylinder separation and rotation direction are assumed to be fixed.
All simulations are performed with the PyFR solver [39], which uses the flux reconstruction (FR) approach [40], a high-order numerical method for unstructured grids. The computational domain is defined as [50] in the stream-wise direction and [
30] in the cross-wise direction. The diameter of each cylinder is set to a value of one; the upper cylinder is centered at (0, 2), while the lower cylinder is centered at (0
2). A spatially varying velocity is applied to the surface of each cylinder to model cylinder rotation. The full mesh, depicted in the left figure of Fig. 5, consists of 3148 unstructured, quadratically curved quadrilateral elements. To align with previous work, the Mach number is set to M = 0.1, which implies a free-stream velocity of
Training a model to generate flow quantities across the entire region shown in the left figure of Fig. 5 would be challenging due to the limitations that memory constraints impose upon the size of neural network inputs. Furthermore, the non-uniform spacing between the solution points in the CFD simulations also presents a challenge, because many of the most powerful neural network architectures are designed to handle image data, and thus assume uniform spacing between points. To account for these challenges, the generative models are only trained to simulate fluid flow on a grid of equispaced points in the vicinity of the cylinders. Formatting the data in this manner allows for training on a set of image-like inputs of dimension 1754, where the first two input dimensions correspond to spatial dimensions and the final input dimension corresponds to different flow quantities. The
Figure 5: Left: View of mesh used for simulation of counter-rotating cylinders over entire computational domain. Right: Zoomed-in view of the mesh in the vicinity of the cylinders. Generative models are trained to only simulate flow in the region with fine mesh resolution.
Figure 6: Format of inputs to the generative model for Re = 150 and Ω = 0. Different physical quantities are treated as different channels in the input.
modeled region spans from in the stream-wise direction and
in the cross-wise direction, which is selected to coincide with the area of fine mesh resolution depicted in the right image of Fig. 5. Figure 6 shows an example training input at Reynolds number Re = 150 and rotation speed Ω = 0; the modeled flow quantities are density, x-velocity, y-velocity, and pressure. In all neural network inputs, flow quantities inside the cylinders are set to be free-stream values, and the network predictions for those points are overwritten with free-stream values during training and evaluation. The training data is generated by running CFD simulations at a range of flow conditions and saving image-like inputs every 500 solver steps.
Training data is generated at a range of flow conditions, corresponding to a variety of Reynolds numbers and rotation speeds. The set of Reynolds numbers considered is Re . This set is designed to span the range of Reynolds numbers over which vortex shedding occurs and the flow can still be considered two-dimensional. Because viscous forces are inversely proportional to Reynolds number for a fixed flow velocity, a change in Reynolds number has a more pronounced effect on flow behavior at Re = 75 than it does at e.g. Re = 150, and hence a finer spacing between simulation points is used at lower Reynolds numbers. At each Reynolds number, simulations are run at rotation speeds Ω
. This range of rotation speeds is selected such that the training data contains information about the system behavior with no cylinder rotation, as well as the behavior of the flow as Ω is increased, including rotation speeds in excess of Ω
. Figure 7 shows the points within the parameter space at which simulations are performed.
Figure 7: Illustration of parameter space for counter-rotating cylinder simulations. The points at which CFD simulations are run are represented by the black dots.
The hope is that a generative model could be trained on CFD data from each point in the plot, and subsequently be capable of simulating fluid flow anywhere in the shaded parameter domain in Fig. 7. By examining the training data, the effects of varying the rotation speed and Reynolds number can be observed. Figure 8 shows the time evolution of the lift coefficients for the top and bottom cylinders as a function of rotation speed at Reynolds number Re = 150. The lift coefficient is defined as:
where L is the lift force, is the free-stream density,
is the free-stream velocity, and d is the cylinder diameter. From the individual plots, two trends can be observed. First, the lift curves separate as the rotation speed is increased, with the mean lift coefficient becoming more positive on the top cylinder and more negative on the bottom cylinder. Second, the amplitude of oscillation in the lift curves is damped out at higher rotation speeds, with a constant transverse force acting on the cylinders at a rotation speed of Ω = 2. From examining these plots, it is apparent that the critical rotation speed, at which the flow transitions from unsteady to steady flow, is somewhere in the interval Ω
Figures 9 and 10 show the change in behavior of the cylinder system as the rotation speed and Reynolds number are varied. Figure 9 visualizes the x-velocity component of the training inputs for four different rotation speeds at a Reynolds number of Re = 150. It is clear from these images that the curved wake regions of low x-velocity begin to flatten out and emanate from different locations on the cylinders as the rotation speed increases. Figure 10 likewise illustrates the variation in flow behavior for four Reynolds numbers with a fixed rotation speed of Ω = 0. It can be seen that the extent of the region of low x-velocity decreases with increasing Reynolds number. Although it cannot be observed from these static images, an increase in Reynolds number also corresponds to a decreased period of vortex shedding.
Figure 8: Lift coefficients on top and bottom cylinder over time as a function of rotation speed.
3.2 Two-Dimensional Flow Training Details
This section provides details about training parameter-conditioned generative models on data from the counter-rotating cylinder system. All components of the generative model are jointly defined and trained in Tensorflow [41]. The parameters to the prior and approximate posterior distributions over , denoted by
) and
), are output by feedforward neural networks with a single hidden layer containing 128 units. The prior and approximate posterior distributions over all remaining latent states, denoted by
) and
), are represented by two-layer feedforward neural networks with 256 and 128 units in each hidden layer. Long short term memory recurrent neural networks with 256 hidden units are used to summarize information about the (latent) state sequences upon which these distributions are conditioned. As explained in Section 2.3, defining the distributions in this manner implicitly assumes that the dynamics of the latent states are non-Markovian; treating the latent states this way was empirically found to yield better performance than assuming Markovian dynamics in the latent space.
As defined above, the networks used to represent the approximate posterior distributions
Figure 9: Variation of x-velocity components of training inputs as a function of rotation speed at Reynolds number of Re = 150.
Figure 10: Variation of x-velocity components of training inputs as a function of Reynolds number at a rotation speed of Ω = 0.
may have an insufficient number of parameters to learn effectively from the very high-dimensional flow states. Therefore, a feature extractor network is introduced, which is trained to provide a mapping from the network inputs, , to a low-dimensional set of features. The approximate posterior distributions over the latent states are then conditioned upon these lower-dimensional features rather than the full flow state.
The feature extractor is comprised of a series of two-dimensional convolutional layers with ResNet [42] skip connections and 256, 128, 64, 32, 16, 8, and 4 filters in each layer. Convolutions are performed with a stride of two in the first five layers, and a striding of one is employed in the remaining layers. An affine transformation maps the output of the convolutional layers to a set of features with the same dimensionality as the latent states; based on performance in numerical experiments, the latent states are defined to be four-dimensional in all studies involving the counter-rotating cylinder system. The decoder network, denoted by ), is constructed to invert all operations performed by the feature extractor, which allows for the mapping of latent state samples back to the state space. Finally, the reconstruction network introduced in Section 2.5, denoted by
is defined to have a single 128-unit hidden layer. Rectified linear units [43] are used as the activation functions across all feedforward layers in the model.
During training, all network parameters are concurrently optimized to maximize the objective defined by Eq. (19). The likelihood components of the evidence lower bound are optimized by minimizing the error between the true flow states,
, and the output of the decoder network. This is equivalent to assuming that
) is a Laplacian distribution, where the mean of the distribution is output by the decoder network and the variance is constant.
error is used rather than
error because the
norm is more sensitive to small errors, which encourages more accurate simulation of the subtle temporal variations in the modeled flow quantities. To allow for the ELBO and mutual information objectives to be of the same scale, a weight of
is applied to the mutual information objective.
The generative model is trained on 20 distinct 32-step sequences from each flow condition
present within the training data. Training is initialized with a learning rate of 6 10
, which is decayed by a factor of 0.75 whenever the loss fails to decrease on a validation set. Training terminates once the learning rate falls below a value of 1
10
. The Adam optimizer [44] is used to perform updates to the model parameters. In total, training takes approximately 30 hours with the neural network parameters divided across two NVIDIA K40 GPUs. While this represents a significant upfront cost, flow simulations generated by the trained model are found to be approximately 120 times faster than the CFD simulations used to generate the training data. This speedup suggests that an immense amount of time could be saved by using generative models in place of CFD solvers, assuming that the generative models prove capable of performing accurate simulations at a wide range of flow conditions. The next section presents numerical experiments that evaluate the learned model’s ability to generate flow solutions for the counter-rotating cylinder system.
3.3 Two-Dimensional Flow Results
The quality of generated solutions can be evaluated by first qualitatively studying whether the trained model is capable of generating realistic solutions for multiple settings of the prescribed parameter c. Figure 11 presents a visual comparison of x-velocity values over time at four distinct flow conditions. For each flow condition, the top row of images corresponds to the ground-truth CFD solutions, while the bottom row contains solutions that were sampled from the generative model. To facilitate easier comparison, care was taken to approximately align the CFD and generated solutions in time such that both simulations show the system at similar moments within the vortex-shedding cycle. However, it should be noted that the initial conditions for the respective solutions are not the same, so discrepancies that arise should not necessarily be attributed to modeling errors by the generative model.
The results in Fig. 11 provide confidence that the generative model effectively captures the changes in flow behavior that occur as the Reynolds number and rotation speed are varied. Additionally, the generative model is shown to not just generate individual time snapshots that are visually similar to the corresponding fluid flow, it also generates solutions that appear to evolve in the same manner as the fluid flow. This temporal consistency is critical if the generative model is to be used to simulate the system for long time periods.
The most noteworthy takeaway from Fig. 11 is the visual similarity between the CFD and generated solutions at Reynolds number Re = 175 and rotation speed Ω = 1. These flow conditions are not contained within the training data, and yet the generative model is still able to generate solutions that are qualitatively similar to the ground-truth fluid flow. Thus, these results imply that it may be possible to use the trained generative model for simulation at many flow conditions, including flow conditions not contained in the training data.
A more quantitative manner of assessing the accuracy of the generated solutions is to look at how the lift forces acting upon the cylinder vary with prescribed Reynolds number and rotation speed. As noted previously in the discussion about Fig. 8, the lift curves for the top and bottom cylinders tend to separate as the rotation speed is increased, with the mean lift on the top cylinder becoming more positive and the mean lift on the bottom cylinder becoming more negative. To examine if the generated solutions are capable of recreating this trend, Fig. 12 presents the mean lift coefficient for the top and bottom cylinders as a
Figure 11: Visual comparison of x-velocity over time at four flow conditions. Solutions generated by the CFD solver are shown in the top row, while flows sampled from the generative model are shown in the bottom row.
Figure 12: Mean lift coefficient as a function of rotation speed for top and bottom cylinders with Reynolds number Re = 150.
Figure 13: Variation of dominant frequency in lift signal as a function of rotation speed and Reynolds number for CFD and generated solutions.
function of rotation speed for flow at Reynolds number Re = 150. It is clear from these plots that the mean lift coefficients calculated from the generated solutions are nearly identical to those determined from the CFD solutions.
In addition to trends in the mean lift coefficient, it can also be studied whether the generative model is capable of capturing variations in the frequency content of the lift curves. The time variation of the lift coefficient values presented in Fig. 8 is roughly sinusoidal, meaning that the corresponding power spectral density (PSD) should show a single peak at a dominant frequency denoted by . Figure 13 shows how the dominant frequency in the lift signal for the top cylinder varies with rotation speed and Reynolds number in both the CFD and generated solutions. Results are only shown up to a rotation speed of Ω = 1.5 because Ω = 1.75 is beyond the critical rotation speed for lower Reynolds numbers.
Figure 14: Depiction of wake measurement location.
Figure 15: Comparison of pressure values over time for CFD solutions (top) and generated solutions (bottom) as a function of rotation speed at Reynolds number Re = 150.
The curves associated with CFD solutions show an increase in the dominant frequency with Reynolds number and a slight decrease with rotation speed. The curves associated with the generated solutions appear to capture these trends quite well, with only minor discrepancies visible between the dominant frequencies identified from the CFD and generated solutions.
Lift is determined by integrating pressure values over the surface of an object. As such, lift values can be largely insensitive to local errors in the pressure distribution, as long as errors at one location on the object surface are balanced out by errors elsewhere. Therefore, the next set of numerical experiments evaluates how accurately the generative model is able to simulate the time evolution of flow quantities at a single point in space. The selected point, represented by the black dot in Fig. 14, is positioned directly in the wake of the top cylinder, located approximately two cylinder diameters downstream. The following discussion is restricted to the modeling of pressure values at this measurement location; similar results were obtained for density, x-velocity, and y-velocity.
Figure 15 shows the time evolution of pressure values at the measurement location as a function of rotation speed for flow at Reynolds number Re = 150. The top row shows the pressure signal obtained from the CFD simulations, while the bottom row shows the corresponding results for the generated solutions. These plots show that the nature of the
Figure 16: Variation of dominant frequency in pressure signal as a function of rotation speed and Reynolds number for CFD and generated solutions.
pressure signal changes dramatically as the rotation speed increases, and the generative model is able to track these changes quite well. Furthermore, note that the generative model correctly predicts a steady solution at Ω = 2.
Once again, the similarity between these signals can be evaluated quantitatively by considering the peaks in their power spectral density. In contrast with the PSD for the lift signals, the PSD for the pressure signals will generally have more than one peak. The comparison here thus considers only the frequency corresponding to the largest peak. Figure 16 shows the variation of with rotation speed and Reynolds number, and demonstrates strong agreement between the curves associated with the CFD and generated solutions. The largest relative error across all flow quantities, Reynolds numbers, and rotation speeds is 11%, suggesting that the generative model is quite effective at modeling the local properties of the flow. Furthermore, note that these curves include results from Reynolds number Re = 175, which is not present in the training data, implying that the model could be effective at modeling the flow properties even at unobserved flow conditions.
Having established confidence that the generative model can be used for accurate simulation at multiple flow conditions, the model is now used to carry out an experiment intended to mimic a design procedure. A design problem typically consists of a set of parameters and objectives, and the goal is to identify design parameter values that perform well according to the design objectives. Exhaustive searches of the design space can be impractical if evaluations of individual design parameters are expensive, as is the case with CFD simulations. However, given the relative low cost of sampling solutions from the generative model, it may be possible to rapidly simulate a given system at a variety of design parameters with the generative model in order to identify promising regions of the design space. In this experiment, the parameter-conditioned generative model is incorporated into a design-like procedure, where an exhaustive search across the parameter space is performed in order to estimate how the critical rotation speed, Ω, varies with Reynolds number.
Figure 17: Left: Norm of time derivative in generated solutions as a function of Ω at Re = 200. Right: Plot of parameter space with points representing identified critical rotation speeds and the plausible region represented by the shaded area.
The goal in this experiment is to identify, as a function of Reynolds number, a plausible range of rotation speeds within which the true critical rotation speed is likely to lie. This is accomplished by first using the generative model to generate solutions at a large number of Reynolds numbers and rotation speeds: Re This requires the generative model to simulate the cylinder system at nearly 13 000 different flow conditions. In total, running these simulations with the generative model takes a little under 10 hours; in contrast, running the same number of simulations with a CFD solver would take nearly two months.
The critical rotation speed corresponds to the transition from unsteady to steady flow. Thus, the critical rotation speed can be identified by considering the time variation between successive snapshots in generated solutions; this time variation should go to zero when steady flow is achieved. The blue line in the left plot of Fig. 17 shows the norm of the numerically estimated time derivative of generated solutions as a function of rotation speed for Re = 200. Note that the time derivatives do tend to decrease with rotation speed, with a sharp drop-off observed at Ω 9. Given the smooth variation in the time derivatives, it is somewhat difficult to identify definitively which rotation speed should correspond to Ω
. Instead, two time-derivative thresholds, represented by the black dotted lines in Fig. 17, are considered, which are assumed to yield a lower- and upper-bound on the true critical rotation speed. By considering these two thresholds, a plausible range for the critical rotation speed:
can be identified.
Furthermore, by identifying the plausible range at each Reynolds number, a plausible region of the parameter space can be determined, wherein all critical rotation speeds are assumed to lie. The shaded area in the right plot of Fig. 17 represents the plausible region identified from the simulations generated by the parameter-conditioned generative model. The size of this plausible region represents a significant reduction relative to the full
Figure 18: Visualization of the mesh for the half-cylinder domain.
parameter space considered in Fig. 7. To determine if this plausible region is reasonable, a bisection search was run at Reynolds numbers Re the true critical rotation speeds according to the CFD solver. These identified values are represented by the plotted points in Fig. 17, and, importantly, are all found to lie within the plausible region.
This section demonstrated that parameter-conditioned generative models can be used to efficiently and accurately simulate two-dimensional fluid flows. The next section studies whether the same techniques can be applied to modeling three-dimensional flows.
3.4 Three-Dimensional Flow over a Half-Cylinder
The test case under consideration is three-dimensional flow over an extruded half-cylinder at an effectively incompressible Mach number of M = 0.2. This flow exhibits several complex flow features, including separated shear layers, turbulent transition, and a fully turbulent wake. The dynamics of the half-cylinder are noteworthy relative to the circular cylinder because the flow separation point remains fixed at the boundary between the curved and flat surfaces. Flow over this configuration has been the subject of several numerical and experimental studies [13, 45–47].
In the given test case, the cylinder is taken to have a diameter d = 1 along its major axis. The domain is taken to be [9], and [0
] in the stream-, cross-, and span-wise directions, and the cylinder is positioned such that the back surface is centered at (0, 0, 0). A visualization of the mesh can be found in Fig. 18. The domain is divided into
non-overlapping, conforming hexahedral elements. As in Section 3.1, the system is simulated with the PyFR solver.
Because it is infeasible to train a model on data from the entire computational domain, the models are once again trained only on data sampled from equispaced points in the vicinity of the cylinder. The data at these points is obtained through multi-linear interpolation, where the ParaView software [48] is used to perform both the interpolation and flow visualizations. The sampling must be performed on a grid with sufficient fineness such that the vortical structures in the flow can be recovered. Thus, samples are drawn at a grid of points that is 128 64
32, spanning a domain of [
11], [
75], and [0
] in the stream-,
Figure 19: Visual comparison between flow generated by CFD simulation of entire computational domain and solution obtained by sampling flow quantities at a grid of equispaced points. The visualizations show iso-surfaces of vorticity magnitude, colored according to pressure.
cross-, and span-wise directions, respectively. This sampling procedure yields training inputs that are 128 64
32
5, where the first three input dimensions correspond to spatial dimensions and the final input dimension corresponds to the five modeled flow quantities: density, x-momentum, y-momentum, z-momentum, and pressure.
Figure 19 shows a qualitative comparison of the flow obtained by running CFD simulations across the entire computational domain against the flow obtained by sampling the flow quantities on the specified grid of points. The visualizations show iso-surfaces of vorticity magnitude, colored according to pressure. It can be seen that the sampled solution retains many of the key features of the original flow, but does not contain any information about the fluid flow far downstream of the cylinder. Hence, the sampled solution represents an approximation to the flow field obtained through CFD simulations, but such an approximation is necessary for training to be tractable.
3.5 Three-Dimensional Flow Training Details
In this experiment, only the Reynolds number is varied. Training data is generated by simulating the half-cylinder system at Reynolds numbers of Re flow at these Reynolds numbers exhibits a range of behaviors, with solutions at Re = 175 showing mostly laminar flow with some turbulence downstream, and solutions at Re = 300 exhibiting a significant amount of turbulence in the cylinder wake. The neural network architectures for the generative model are largely identical to those described for the two-dimensional flow in Section 3.2, with one notable exception: the two-dimensional convolutions in the feature extractor and decoder network are replaced with three-dimensional convolutions. Additionally, to account for the added complexities inherent to modeling three-dimensional, turbulent flow, 32-dimensional latent states are used in place of the four-dimensional latent states employed previously.
Parameter-conditioned generative models are trained on 165 distinct sequences from each flow condition. Training a single model takes approximately three days with the model parameters divided across two NVIDIA K40 GPUs and one NVIDIA GeForce GTX 1070 GPU. The longer training time can be attributed to the larger size of the training inputs and the additional computational cost associated with performing three-dimensional convolutions. Flow simulations generated by the trained model are once again found to be approximately 120 times faster than the corresponding CFD simulations. The next section studies how well the generative models are able to model the fluid flow.
3.6 Three-Dimensional Flow Results
The results presented here are solely meant to demonstrate that the generative modeling techniques can scale to larger problems; an in-depth quantitative evaluation of the generated results is reserved for future work. Therefore, only a qualitative evaluation of the generated solutions is presented in this section. Figures 20 and 21 present a visual comparison between the CFD and generated solutions at the four flow conditions contained within the training data. For each flow condition, the top row of images contains results from the CFD simulations, while the bottom row contains results sampled from the generative model. Figure 20 shows iso-surfaces of Q, the second invariant of the velocity gradient tensor, colored according to vorticity in the downstream direction. Figure 21 shows iso-surfaces of vorticity magnitude, colored according to pressure.
A strong visual similarity can be observed between the CFD solutions and the solutions sampled from the generative model. In particular, the generative model seems effective at capturing the variation in the quantity and nature of vortical structures with Reynolds number. It is worth emphasizing that the iso-surfaces in all images are functions of spatial derivatives of the flow velocities. These spatial derivatives are estimated through finite difference, and thus the presence of the iso-surfaces is strongly reliant upon accurately modeling the manner in which the velocity field varies in the cylinder wake. Hence, it is noteworthy that the trained model is capable of generating flow fields that not only contain these vortical structures, but also retain them across successive time steps while propagating them downstream.
This paper presented a method for learning parameter-conditioned sequential generative models. The central element of the proposed method is a variational inference procedure that enables the discovery of low-dimensional latent states capable of accurately capturing the dynamics of a modeled system. These inferred latent representations can be conditioned on parameters that are assumed to govern the modeled system’s behavior, thereby enabling the generative modeling of dynamical systems, where the modeled dynamics are a function of prescribed parameters.
The proposed method was first evaluated based on its ability to simulate two-dimensional airflow over a pair of counter-rotating cylinders. Extensive qualitative and quantitative experiments demonstrated that learned generative models were capable of effectively modeling both local and global properties of the flow field at a wide range of flow conditions. Furthermore, it was found that, relative to a CFD solver, a speedup of approximately 120can be obtained by running simulations with a trained generative model. A final set of
Figure 20: Comparison of iso-surfaces of Q-values, colored according to downstream vorticity, over time at four flow conditions. CFD solutions are shown in the top row, while generated solutions are shown in the bottom row.
Figure 21: Comparison of iso-surfaces of vorticity magnitude, colored according to pressure, over time at four flow conditions. CFD solutions are shown in the top row, while generated solutions are shown in the bottom row.
experiments with data from three-dimensional, turbulent flow simulations demonstrated that the generative modeling techniques can scale to complex and high-dimensional problems.
There are some limitations to the proposed approach that should be addressed in future work. First, in all experiments it was assumed that the flow geometries were fixed. However, many aerospace design problems, such as airfoil or wing design, focus predominantly on how changes in geometry affect the resulting flow. At the moment it is unclear how well generative models would perform if, for example, the separation between the cylinders were allowed to vary, thereby altering the location of the cylinders within the flow field. Hence, more work is needed to identify and address potential shortcomings in this regard. One further limitation is that the employed neural network architectures require the flow quantities to be sampled at uniformly spaced points. Future work should be dedicated to overcoming this restriction, which will make the techniques better suited to modeling fluid flows from simulations with unstructured grids.
This material is based upon work supported by the Stanford Strategic Energy Alliance. The authors would like to thank Antony Jameson for valuable guidance and feedback.
[1] J. Sacks, S. B. Schiller, and W. J. Welch, “Designs for computer experiments,” Technometrics, vol. 31, no. 1, pp. 41–47, 1989.
[2] N. V. Queipo, R. T. Haftka, W. Shyy, T. Goel, R. Vaidyanathan, and P. K. Tucker, “Surrogate-based analysis and optimization,” Progress in Aerospace Sciences, vol. 41, no. 1, pp. 1–28, 2005.
[3] A. I. Forrester and A. J. Keane, “Recent advances in surrogate-based optimization,” Progress in Aerospace Sciences, vol. 45, no. 1-3, pp. 50–79, 2009.
[4] K. Willcox and J. Peraire, “Balanced model reduction via the proper orthogonal decomposition,” AIAA Journal, vol. 40, no. 11, pp. 2323–2330, 2002.
[5] K. Carlberg, C. Bou-Mosleh, and C. Farhat, “Efficient non-linear model reduction via a least-squares petrov–galerkin projection and compressive tensor approximations,” International Journal for Numerical Methods in Engineering, vol. 86, no. 2, pp. 155– 181, 2011.
[6] K. Lee and K. Carlberg, “Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders,” ArXiv preprint arXiv:1812.08373, 2018.
[7] ——, “Deep conservation: A latent dynamics model for exact satisfaction of physical conservation laws,” ArXiv preprint arXiv:1909.09754, 2019.
[8] C. White, D. Ushizima, and C. Farhat, “Neural networks predict fluid dynamics solutions from tiny datasets,” ArXiv preprint arXiv:1902.00091, 2019.
[9] S. E. Otto and C. W. Rowley, “Linearly-recurrent autoencoder networks for learning dynamics,” ArXiv preprint arXiv:1712.01378, 2017.
[10] S. C. Puligilla and B. Jayaraman, “Deep multilayer convolution frameworks for data-driven learning of fluid flow dynamics,” in 2018 Fluid Dynamics Conference, 2018.
[11] J. Morton, F. D. Witherden, A. Jameson, and M. J. Kochenderfer, “Deep dynamical modeling and control of unsteady fluid flows,” in Advances in Neural Information Processing Systems (NeurIPS), 2018.
[12] B. Lusch, J. N. Kutz, and S. L. Brunton, “Deep learning for universal linear embeddings of nonlinear dynamics,” Nature Communications, vol. 9, no. 1, p. 4950, 2018.
[13] K. T. Carlberg, A. Jameson, M. J. Kochenderfer, J. Morton, L. Peng, and F. D. Witherden, “Recovering missing CFD data for high-order discretizations using deep neural networks and dynamics learning,” Journal of Computational Physics, vol. 395, pp. 105–124, 2019.
[14] S. Wiewel, M. Becher, and N. Thuerey, “Latent space physics: Towards learning the temporal evolution of fluid flow,” in Computer Graphics Forum, 2019.
[15] B. Kim, V. C. Azevedo, N. Thuerey, T. Kim, M. Gross, and B. Solenthaler, “Deep fluids: A generative network for parameterized fluid simulations,” in Computer Graphics Forum, 2019.
[16] D. P. Kingma and M. Welling, “An introduction to variational autoencoders,” ArXiv preprint arXiv:1906.02691, 2019.
[17] D. P Kingma and M. Welling, “Auto-encoding variational Bayes,” ArXiv preprint arXiv: 1312.6114, 2013.
[18] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial nets,” in Advances in Neural Information Processing Systems (NeurIPS), 2014.
[19] D. J. Rezende and S. Mohamed, “Variational inference with normalizing flows,” in International Conference on Machine Learning (ICML), 2015.
[20] L. Dinh, J. Sohl-Dickstein, and S. Bengio, “Density estimation using real NVP,” in International Conference on Learning Representations (ICLR), 2017.
[21] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe, “Variational inference: A review for statisticians,” Journal of the American Statistical Association, vol. 112, no. 518, pp. 859–877, 2017.
[22] D. P. Kingma, T. Salimans, and M. Welling, “Variational dropout and the local reparameterization trick,” in Advances in Neural Information Processing Systems (NeurIPS), 2015.
[23] K. Sohn, H. Lee, and X. Yan, “Learning structured output representation using deep conditional generative models,” in Advances in Neural Information Processing Systems (NeurIPS), 2015.
[24] J. Durbin and S. J. Koopman, Time series analysis by state space methods. Oxford University Press, 2012.
[25] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neural Computation, vol. 9, no. 8, pp. 1735–1780, 1997.
[26] J. Chung, C. Gulcehre, K. Cho, and Y. Bengio, “Empirical evaluation of gated recurrent neural networks on sequence modeling,” in Workshop on Deep Learning, Advances in Neural Information Processing Systems (NeurIPS), 2014.
[27] D. Barber and F. V. Agakov, “The IM algorithm: A variational approach to information maximization,” in Advances in Neural Information Processing Systems (NeurIPS), 2003.
[28] X. Chen, Y. Duan, R. Houthooft, J. Schulman, I. Sutskever, and P. Abbeel, “InfoGAN: Interpretable representation learning by information maximizing generative adversarial nets,” in Advances in Neural Information Processing Systems (NeurIPS), 2016.
[29] Y. Li, J. Song, and S. Ermon, “InfoGAIL: Interpretable imitation learning from visual demonstrations,” in Advances in Neural Information Processing Systems (NeurIPS), 2017.
[30] K. Roussopoulos, “Feedback control of vortex shedding at low Reynolds numbers,” Journal of Fluid Mechanics, vol. 248, pp. 267–296, 1993.
[31] S. J. Illingworth, “Model-based control of vortex shedding at low Reynolds numbers,” Theoretical and Computational Fluid Dynamics, vol. 30, no. 5, pp. 429–448, 2016.
[32] A. S. Chan, “Control and suppression of laminar vortex shedding off two-dimensional bluff bodies,” PhD thesis, Stanford University, 2012.
[33] C. Williamson, “Evolution of a single wake behind a pair of bluff bodies,” Journal of Fluid Mechanics, vol. 159, pp. 1–18, 1985.
[34] S. Kang, “Characteristics of flow over two circular cylinders in a side-by-side arrangement at low reynolds numbers,” Physics of Fluids, vol. 15, no. 9, pp. 2486–2498, 2003.
[35] I Peschard and P Le Gal, “Coupled wakes of cylinders,” Physical Review Letters, vol. 77, no. 15, p. 3122, 1996.
[36] H. S. Yoon, J. H. Kim, H. H. Chun, and H. J. Choi, “Laminar flow past two rotating circular cylinders in a side-by-side arrangement,” Physics of Fluids, vol. 19, no. 12, p. 128 103, 2007.
[37] A. S. Chan and A. Jameson, “Suppression of the unsteady vortex wakes of a circular cylinder pair by a doublet-like counter-rotation,” International Journal for Numerical Methods in Fluids, vol. 63, no. 1, pp. 22–39, 2010.
[38] A. S. Chan, P. A. Dewey, A. Jameson, C. Liang, and A. J. Smits, “Vortex suppression and drag reduction in the wake of counter-rotating cylinders,” Journal of Fluid Mechanics, vol. 679, pp. 343–382, 2011.
[39] F. D. Witherden, A. M. Farrington, and P. E. Vincent, “PyFR: An open source framework for solving advection–diffusion type problems on streaming architectures using the flux reconstruction approach,” Computer Physics Communications, vol. 185, no. 11, pp. 3028–3040, 2014.
[40] H. T. Huynh, “A flux reconstruction approach to high-order schemes including discontinuous galerkin methods,” in 18th AIAA Computational Fluid Dynamics Conference, 2007.
[41] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, et al., “Tensorflow: Large-scale machine learning on heterogeneous systems,” ArXiv preprint arXiv:1603.04467, 2015.
[42] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR), 2016.
[43] V. Nair and G. E. Hinton, “Rectified linear units improve restricted Boltzmann machines,” in International Conference on Machine Learning (ICML), 2010, pp. 807– 814.
[44] D. Kingma and J. Ba, “Adam: A method for stochastic optimization,” ArXiv preprint arXiv:1412.6980, 2014.
[45] Y Nakamura, “Vortex shedding from bluff bodies with splitter plates,” Journal of Fluids and Structures, vol. 10, no. 2, pp. 147–158, 1996.
[46] S. Kumarasamy and J. B. Barlow, “Computation of unsteady flow over a half-cylinder close to a moving wall,” Journal of Wind Engineering and Industrial Aerodynamics, vol. 69, pp. 239–248, 1997.
[47] A Santa Cruz, L David, J Pecheux, and A Texier, “Characterization by proper-orthogonal-decomposition of the passive controlled wake flow downstream of a half cylinder,” Experiments in Fluids, vol. 39, no. 4, pp. 730–742, 2005.
[48] J. Ahrens, B. Geveci, and C. Law, “Paraview: An end-user tool for large data visualization,” The Visualization Handbook, vol. 717, 2005.