Latent-space Physics: Towards Learning the Temporal Evolution of Fluid Flow

Introduction

The variables we use to describe real world physical systems typically take the form of complex functions with high dimensionality. Especially for transient numerical simulations, we usually employ continuous models to describe how these functions evolve over time. For such models, the field of computational methods has been highly successful at developing powerful numerical algorithms that accurately and efficiently predict how the natural phenomena under consideration will behave. In the following, we take a different view on this problem: instead of relying on analytic expressions, we use a deep learning approach to infer physical functions based on data. More specifically, we will focus on the temporal evolution of complex functions that arise in the context of fluid flows. Fluids encompass a large and important class of materials in human environments, and as such they’re particularly interesting candidates for learning models.

While other works have demonstrated that machine learning methods are highly competitive alternatives to traditional methods, e.g., for computing local interactions of particle based liquids [LJS∗15], to perform divergence free projections for a single point in time [TSSP16], or for adversarial training of high resolution flows [XFCT18], few works exist that target temporal evolutions of physical systems. While first works have considered predictions of Lagrangian objects such as rigid bodies [WZW∗17], and control of two dimensional interactions [MTP∗18], the question whether neural networks (NNs) can predict the evolution of complex three-dimensional functions such as pressure fields of fluids has not previously been addressed. We believe that this is a particularly interesting challenge, as it not only can lead to faster forward simulations, as we will demonstrate below, but also could be useful for giving NNs predictive capabilities for complex inverse problems.

The complexity of nature at human scales makes it necessary to finely discretize both space and time for traditional numerical methods, in turn leading to a large number of degrees of freedom. Key to our method is reducing the dimensionality of the problem using convolutional neural networks (CNNs) with respect to both time and space. Our method first learns to map the original, three-dimensional problem into a much smaller spatial latent space, at the same time learning the inverse mapping. We then train a second network that maps a collection of reduced representations into an encoding of the temporal evolution. This reduced temporal state is then used to output a sequence of spatial latent space representations, which are decoded to yield the full spatial data set for a point in time. A key advantage of CNNs in this context is that they give us a natural way to compute accurate and highly efficient non-linear representations. We will later on demonstrate that the setup for computing this reduced representation strongly influences how well the time network can predict changes over time, and we will demonstrate the generality of our approach with several liquid and single-phase problems. The specific contributions of this work are:

a first LSTM architecture to predict temporal evolutions of dense, physical 3D functions in learned latent spaces,

an efficient encoder and decoder architecture, which by means of a strong compression, yields a very fast simulation algorithm,

in addition to a detailed evaluation of training modalities.

Related Work and Background

Despite being a research topic for a long time [RHW88], the interest in neural network algorithms is a relatively new phenomenon, triggered by seminal works such as ImageNet [KSH12]. In computer graphics, such approaches have led to impressive results, e.g., for synthesizing novel viewpoints of natural scenes [FNPS16], to generate photorealistic face textures [SWH∗16], and to robustly transfer image styles between photographs [LPSB17], to name just a few examples. The underlying optimization approximates an unknown function $f^{*}(x)=y$, by minimizing an associated loss function $L$ such that $f(x,\theta)\approx y$. Here, $\theta$ denotes the degrees of freedom of the chosen representation for $f$. For our algorithm, we will consider deep neural networks. With the right choice of $L$, e.g., an $L_{2}$ norm in the simplest case, such a neural network will approximate the original function $f^{*}$ as closely as possible given its internal structure. A single layer $l$ of an NN can be written as $a^{l}=\sigma(W_{l}a^{l-1}+b_{l})$, where $a^{i}$ is the output of the i’th layer, $\sigma$ represents an activation function, and $W_{l},b_{l}$ denote weight matrix and bias, respectively. In the following, we collect the weights $W_{l},b_{l}$ for all layers $l$ in $\theta$.

The latent spaces of generative NNs were shown to be powerful tools in image processing and synthesis [RMC16, WZX∗16]. They provide a non-linear representation that is closely tied to the given data distribution, and our approach leverages such a latent space to predict the evolution of dense physical functions. While others have demonstrated that trained feature spaces likewise pose very suitable environments for high-level operations with impressive results [UGB∗16], we will focus on latent spaces of autoencoder networks in the following. The sequence-to-sequence methods which we will use for time prediction have so far predominantly found applications in the area of natural language processing, e.g., for tasks such as machine translation [SVL14]. These recurrent networks are especially popular for control tasks in reinforcement learning environments [MBM∗16]. Recently, impressive results were also achieved for tasks such as automatic video captioning [XYZM17].

Using neural networks in the context of physics problems is a new area within the field of deep learning methods. Several works have targeted predictions of Lagrangian objects based on image data. E.g., Battaglia et al. \shortcitebattaglia2016interaction predict two-dimensional physics, a method that was subsequently extended to videos [WZW∗17]. Another line of work proposed a specialized architecture for two-dimensional rigid bodies physics [CUTT16], while others have targeted predictions of liquid motions for robotic control [SF17] Farimani et al. [FGP17] proposed an adversarial training approach to infer solutions for two-dimensional physics problems, such as heat diffusion and lid driven cavity flows. Other researchers have proposed networks to learn PDEs [LLMD17] by encoding the unknown differential operators with convolutions. To model the evolution of a subsurface multiphase flow others are using proper orthogonal decomposition in combination with recurrent neural networks (RNNs) [KE18]. In the field of weather prediction, short-term prediction architectures evolved that internally also make use of RNN layers [SCW∗15]. Other works create velocity field predictions by learning parameters for Reynolds Averaged Navier Stokes models from simulation data [LKT16]. In addition, learned Koopman operators [MWJK18] were proposed for representing temporal changes, whereas Lusch et al. are searching for representations of Koopman eigenfunctions to globally linearize dynamical systems using deep learning [LKB18]. While most of these works share our goal to infer Eulerian functions for physical models, they are limited to relatively simple, two dimensional problems. In contrast, we will demonstrate that our reduced latent space representation can work with complex functions with up to several million degrees of freedom.

We focus on flow physics, for which we employ the well established Navier-Stokes (NS) model. Its incompressible form is given by

where the most important quantities are flow velocity $\mathbf{u}$ and pressure $p$. The other parameters $\rho,\nu,\mathbf{g}$ denote density, kinematic viscosity and external forces, respectively. For liquids, we will assume that a signed-distance function $\phi$ is either advected with the flow, or reconstructed from a set of advected particles.

In the area of visual effects, Kass and Miller were the first to employ height field models [KM90], while Foster and Metaxas employed a first three-dimensional NS solver [FM96]. After Jos Stam proposed an unconditionally stable advection and time integration scheme \shortcitestam1999, it led to powerful liquid solvers based on the particle levelset [FF01], in conjunction with accurate free surface boundary conditions [ENGF03]. Since then, the fluid implicit particle (FLIP) method has been especially popular for detailed liquid simulations [ZB05], and we will use it to generate our training data. Solvers based on these algorithms have subsequently been extended with accurate boundary handling [BBB07] synthetic turbulence [KTJG08], or narrow band algorithms [FAW∗16], to name just a few examples. A good overview of fluid simulations for computer animation can be found in the book by R. Bridson \shortcitebridson2015. Aiming for more complex systems, other works have targeted coupled reduced order models, or sub-grid coupling effects [TLK16, FMB∗17]. While we do not target coupled fluid solvers in our work, these directions of course represent interesting future topics.

Beyond these primarily grid-based techniques, smoothed particle hydrodynamics (SPH) are a popular Lagrangian alternative [MCG03, MMCK14]. However, we will focus on Eulerian solvers in the following, as CNNs are particularly amenable to grid-based discretizations. The pressure function has received special attention, as the underlying iterative solver is often the most expensive component in an algorithm. E.g., techniques for dimensionality reduction [LZF10, ATW15], and fast solvers [MST10, ICS∗14] have been proposed to diminish its runtime impact.

In the context of fluid simulations and machine learning for animation, a regression forest based approach for learning SPH interactions has been proposed by Ladicky et al. \shortciteladicky2015data. Other graphics works have targeted learning flow descriptors with CNNs [CT17], or learning the statistics of splash formation [UHT17], and two-dimensional control problems [MTP∗18]. While pressure projection algorithms with CNNs [TSSP16, YYX16] shares similarities with our work on first sight, they are largely orthogonal. Instead of targeting divergence freeness for a single instance in time, our work aims for learning its temporal evolution over the course of many time steps. An important difference is also that the CNN-based projection so far has only been demonstrated for smoke simulations, similar to other model-reduced simulation approaches [TLP06]. For all of these methods, handling the strongly varying free surface boundary conditions remains an open challenge, and we demonstrate below that our approach works especially well for liquids.

Method

Given a set of parameters $\theta$ and a functional representation $f_{t}$ our goal is to predict the $o$ future states $\mathbf{x}(t+h)$ to $\mathbf{x}(t+oh)$ as closely as possible given a current state and a series of $n$ previous states, i.e.,

To provide better visual clarity the set of weights, i.e. learnable parameters, $\theta$ is omitted in the function definitions.

We will use CNNs for $f_{d}$ and $f_{e}$, and thus the space $\mathbf{c}$ is given by their learned latent space. We choose its dimensionality $m_{s}$ such that the temporal prediction problem above becomes feasible for dense three dimensional samples.

In order to reduce the spatial dimensionality of our inference problem, we employ a fully convolutional autoencoder architecture [MMCS11]. Our autoencoder (AE) consists of the aforementioned encoding and decoding functions $f_{e},f_{d}$ and is trained to reconstruct the quantity $\mathbf{x}$ as accurately as possible w.r.t. an $L_{2}$ norm, i.e.

where $\theta_{d},\theta_{e}$ represent the parameters of the decoder and encoder, respectively. We use a series of convolutional layers activated by leaky rectified linear units (LeakyReLU) [MHN13] for encoder and decoder, with a bottleneck layer of dimensionality $m_{s}$. This layer yields the latent space encoding that we use to predict the temporal evolution. Both encoder and decoder consist of 6 convolutional layers that increase / decrease the number of features by a factor of 2. In total, this yields a reduction factor of 256.

In the following, we will explain additional details of the autoencoder pre-training, and layer setup. We denote layers in the encoder and decoder stack as $f_{e_{i}}$ and $f_{d_{j}}$, where $i,j\in[0,l]$, with $i,j$ being integers, denote the depth from the input and output layers, $l$ being the depth of the latent space layer. In our network architecture, encoder and decoder layers with $i=j$ have to match, i.e., the output shape of $f_{e_{i}}$ has to be identical to that of $f_{d_{j}}$ and vice versa. This setup allows for a greedy, layer-wise pretraining of the autoencoder, as proposed by Bengio et al. [BLPL07], where beginning from a shallow single layer deep autoencoder, additional layers are added to the model forming a series of deeper models for each stage. The optimization problem of such a stacked autoencoder in pretraining is therefore formulated as

with $\theta_{e_{0...k}},\theta_{d_{0...k}}$ denoting the parameters of the sub-stack for pretraining stage $k$, and $\circ$ denoting composition of functions. For our final models, we typically use a depth $l=5$, and thus perform 6 runs of pretraining before training the complete model. This is illustrated in Fig. 1, where the paths from $f_{e_{k}}$ over $c_{k}$ to $f_{d_{k}}$ are only active in pretraining stage $k$. After pretraining only the path $f_{e_{5}}$ over $c_{5}$ to $f_{d_{5}}$ remains active.

In addition, our autoencoder does not use any pooling layers, but instead only relies on strided convolutional layers. This means we apply convolutions with a stride of $s$, skipping $s-1$ entries when applying the convolutional kernel. We assume the input is padded, and hence for $s=1$ the output size matches the input, while choosing a stride $s>1$ results in a downsampled output [ODO16]. Equivalently the decoder network employs strided transposed convolutions, where strided application increases the output dimensions by a factor of $s$. The details of the network architecture, with corresponding strides and kernel sizes can be found in Table 1.

In addition to this basic architecture, we will also evaluate a variational autoencoder [RMW14] in Sec. 5 that enforces a normalization on the latent space while keeping the presented AE layout identical. As no pre-trained models for physics problems are available, we found greedy pre-training of the autoencoder stack to be crucial for a stable and feasible training process.

2 Prediction of Future States

Here $\theta_{t}$ denotes the parameters of the prediction network, and $[\cdot,\cdot]$ denotes concatenation of the $\mathbf{c}$ vectors.

In contrast to the spatial reduction network above, which receives the full spatial input at once and infers a latent space coordinate without any data internal to the network, the prediction network uses a recurrent architecture for predicting the evolution over time. It receives a series of inputs one by one, and computes its output iteratively with the help of an internal network state. In contrast to the spatial reduction, which very heavily relies on convolutions, we cannot employ similar convolutions for the time data sets. While it is a valid assumption that each entry of a latent space vector $\mathbf{c}$ varies smoothly in time, the order of the entries is arbitrary and we cannot make any assumptions about local neighborhoods within $\mathbf{c}$. As such, convolving $\mathbf{c}$ with a filter along the latent space entries typically does not give meaningful results. Instead, our prediction network will use convolutions to translate the LSTM state into the latent space, in addition to fully connected layers of LSTM units.

Note that the iterative nature is shared by encoder and decoder module of the prediction network, i.e., the encoder actually internally produces $n+1$ contexts, the first $n$ of which are intermediate contexts. These intermediate contexts are only required for the feedback loop internal to the corresponding LSTM layer, and are discarded afterwards. We only keep the very last context in order to pass it to the decoder part of the network. This context is repeated $o$ times, in order for the decoder LSTM to infer the desired future states.

Despite the spatial dimensionality reduction with an autoencoder, the number of weights in LSTM layers can quickly grow due to their inherent internal feedback loops (typically equivalent to four fully connected layers). To prevent overfitting from exploding weight numbers in the LSTM layers, we propose a hybrid structure of LSTM units and convolutions as shown in Fig. 2 that is used instead of the fully recurrent approach presented in Fig. 3.

For our prediction network we use two LSTM layers that infer an internal temporal representation of the data, followed by a final linear, convolutional network that translates the temporal representation into the corresponding latent space point. This convolution effectively represents a translation of the context information from the LSTM layer into latent space points that is constant for all output steps. This architecture effectively prevents overfitting, and ensures a high quality temporal prediction, as we will demonstrate below. In particular, we will show that this hybrid network outperforms networks purely based on LSTM layers, and significantly reduces the weight footprint. Additionally the prediction network architecture can be extended by applying multiple stacked convolution layers after the final LSTM layer. We found this hybrid architecture crucial for inferring the high-dimensional outputs of physical simulations.

While the autoencoder, thanks to its fully convolutional architecture, could be applied to inputs of varying size, the prediction network is trained for fixed latent space inputs, and internal context sizes. Correspondingly, when the latent space size $m_{s}$ changes, it influences the size of the prediction network’s layers. Hence, the prediction network has to be re-trained from scratch when the latent space size is changed. We have not found this critical in practice, because the prediction network takes significantly less time to train than the autoencoder, as we will discuss in Sec. 5.

Fluid Flow Data

To generate fluid data sets for training we rely on a NS solver with operator splitting [Bri15] to calculate the training data at discrete points in space and time. On a high level, the solver contains the following steps: computing motion of the fluid with the help of transport, i.e. advection steps, for the velocity $\mathbf{u}$, evaluating external forces, and then computing the harmonic pressure function $p$. In addition, a visible, passive quantity such as smoke density $\rho$, or a level-set representation $\phi$ for free surface flows is often advected in parallel to the velocity itself. Calculating the pressure typically involves solving an elliptic second-order PDE, and the gradient of the resulting pressure is used to make the flow divergence free.

Overall, these physical data sets differ significantly from data sets such as natural images that are targeted with other learning approaches. They are typically well structured, and less ambiguous due to a lack of projections, which motivates our goal to use learned models. At the same time they exhibit strong temporal changes, as is visible in Fig. 4, which make the temporal inference problem a non-trivial task.

Depending on the choice of physical quantity to infer with our framework, different simulation algorithms emerge. We will focus on velocity $\mathbf{u}$ and the two pressure variants, total $p_{t}$ and split ($p_{s}$ and $p_{d}$) in the following. When targeting $\mathbf{u}$ with our method, this means that we can omit velocity advection as well as pressure solve, while the inference of pressure means that we only omit the pressure solve, but still need to perform advection and velocity correction with the pressure gradient. While this pressure inference requires more computations, the pressure solve is typically the most time consuming part with a super-linear complexity, and as such both options have comparable runtimes. When predicting the pressure field, we also use a boundary condition alignment step for the free surface [ATW15]. It takes the form of three Jacobi iterations in a narrow band at the liquid surface in order to align the Dirichlet boundary conditions with the current position of the interface. This step is important for liquids, as it incorporates small scale dynamics, leaving the large-scale dynamics to a learned model.

A variant for both of these simulation algorithm classes is to only rely on the network prediction for a limited time interval of $i_{p}$ time steps, and then perform a single full simulation step without any network calculations, i.e., for $i_{p}=0$ the network is not used at all, while $i_{p}=\infty$ is identical to the full network prediction described in the previous paragraph. We will investigate prediction intervals on the order of $4$ to $14$ steps. This simulation variant represents a joint numerical time integration and network prediction, that can have advantages to prevent drift from the learned predictions. We will denote such versions as interval predictions below.

2 Data Sets

To demonstrate that our approach is applicable to a wide range of physics phenomena, we will show results with three different 3D data sets in the following. To ensure a sufficient amount of variance with respect to physical motions and dynamics, we use randomized simulation setups. We target scenes with high complexity, i.e., strong visible splashes and vortices, and large CFL (Courant-Friedrichs-Lewy) numbers (typically around 2-3), that measure how fast information travels from cell to cell in a complex simulation domain. For each of our data sets, we generate $n_{s}$ scenes of different initial conditions, for which we discard the first $n_{w}$ time steps, as these typically contain small and regular, and hence less representative dynamics. Afterwards, we store a fixed number of $n_{t}$ time steps as training data, resulting in a final size of $n_{s}n_{t}$ spatial data sets. Each data set content is normalized to the range of .

Two of the three data sets contain liquids, while the additional one targets smoke simulations. The liquid data sets with spatial resolutions of $64^{3}$ and $128^{3}$ contain randomized sloshing waves and colliding bodies of liquid. The scene setup consists of a low basin, represented by a large volume of liquid at the bottom of the domain, and a tall but narrow pillar of liquid, that drops into it. Additionally a random amount, ranging from zero to three smaller liquid drops are placed randomly in the domain. For the $128^{3}$ data set we additionally include complex geometries for the initial liquid bodies, yielding a larger range of behavior. These data sets will be denoted as liquid64 and liquid128, respectively. In addition, we consider a data set containing single-phase flows with buoyant smoke which we will denote as smoke128. We place 4 to 10 inflow regions into an empty domain at rest, and then simulate the resulting plumes of hot smoke. As all setups are invariant w.r.t. rotations around the axis of gravity (Y in our setups), we augment the data sets by mirroring along XY and YZ. This leads to sizes of the data sets from 80k to 400k entries, and the $128^{3}$ data sets have a total size of 671GB. Rendered examples from all data sets can be found in Fig. 16, Fig. 17 and Fig. 18 in the supplemental document, as well as further information about the initial conditions and physical parameters of the fluids.

Evaluation and Training

In the following we will evaluate the different options discussed in the previous section with respect to their prediction accuracies. In terms of evaluation metrics, we will use PSNR (peak signal-to-noise ratio) as a baseline metric, in addition to a surface-based Hausdorff distance in order to more accurately compare the position of the liquid interface [HKR93, LDGN15]. More specifically, given two signed distance functions $\phi_{r},\phi_{p}$ representing reference and predicted surfaces, we compute the surface error as

Unless otherwise noted, the error measurements start after 50 steps of simulation, and are averaged for ten test scenes from the liquid64 setup.

We first evaluate the accuracy of only the spatial encoding, i.e., the autoencoder network in conjunction with a numerical time integration scheme. At the end of a fluid solving time step, we encode the physical variable $\mathbf{x}$ under consideration with $\mathbf{c}=f_{e}(\mathbf{x})$, and then restore it from its latent space representation $\mathbf{x}^{\prime}=f_{d}(\mathbf{c})$. In the following, we will compare flow velocity $\mathbf{u}$, total pressure $p_{t}$, and split pressure ($p_{s}$, $p_{d}$), all with a latent space size of $m_{s}=1024$. We train a new autoencoder for each quantity, and we additionally consider a variational autoencoder for the split pressure. Training times for the autoencoders were two days on average, including pre-training. To train the different autoencoders, we use 6 epochs of pretraining and 25 epochs of training using an Adam optimizer, with a learning rate of 0.001 and a decay factor of 0.005. For training we used 80% of the data set, 10% for validation during training, and another 10% for testing.

Fig. 6(a) and Fig. 6(e) show error measurements averaged for 10 simulations from the test data set. Given the complexity of the data, especially the total pressure variant exhibits very good representational capabilities with an average PSNR value of $69.14$. On the other hand, the velocity encoding introduces significantly larger errors in Fig. 6(e). Interestingly, neither the latent space normalization of the VAE, nor the split pressure data increase the reconstruction accuracy, i.e., the CNN does not benefit from the reduced data range of the pressure splitting approach. A visual comparison of the results can be found in Fig. 5.

2 Temporal Prediction

Next, we evaluate reconstruction quality when including the temporal prediction network. Thus, now a quantity $\mathbf{x}^{\prime}$ is inferred based on a series of previous latent space points. For the following tests, our prediction model uses a history of $6$, and infers the next time step, thus $o=1$, with a latent space size $m_{s}=1024$. For a resolution of $64^{3}$ the fully recurrent network contains 700, and 1500 units for the first and second LSTM layer of Fig. 2, respectively. Hence, $\mathbf{d}$ has a dimensionality of 700 for this setup. The two LSTM layers are followed by a convolutional layer targeting the $m_{s}$ latent space dimensions for our hybrid architecture, or alternatively another LSTM layer of size $m_{s}$ for the fully recurrent version. A dropout rate of $1.32\cdot 10^{-2}$ with a recurrent dropout of $0.385$, and a learning rate of $1.26\cdot 10^{-4}$ with a decay factor of $3.34\cdot 10^{-4}$ were used for all trainings of the prediction network. Training was run for $50$ epochs with RMSProp, with $319600$ training samples in each epoch, taking 2 hours, on average. Hyperparameters as well as the length of the time history used as input for the prediction network and the generated output time steps were chosen by utilizing a hyper parameter search, i.e. training multiple configurations of the same network with differing input-/output counts or hyperparameter settings.

The error measurements for simulations predicted by the combination of autoencoder and prediction network are shown in Fig. 6(b) and Fig. 6(f), with a surface visualization in Fig. 8. Given the autoencoder baseline, the prediction network does very well at predicting future states for the simulation variables. The accuracy only slightly decreases compared to Fig. 6(a) and 6(e), with an average PSNR value of $64.80$ (a decrease of only 6.2% w.r.t. the AE baseline). Here, it is also worth noting that the LSTM does not benefit from the normalized latent space of the VAE. On the contrary, the predictions without the regular AE exhibit a lower error.

Fig. 6(c) shows an evaluation of the interval prediction scheme explained above. Here we employ the LSTM for $i_{p}=14$ consecutive steps, and then perform a single regular simulation step. This especially improves the pressure predictions, for which the average surface error after 100 steps is still below two cells.

We also evaluate how well our model can predict future states based on a single set of inputs. For multiple output steps, i.e. $o>1$, our model predicts several latent space points from a single time context $\mathbf{d}$. A graph comparing accuracy for 1, 3 and 5 steps of output can be found in Fig. 6(h). It is apparent that the accuracy barely degrades when multiple steps are predicted at once. However, this case is significantly more efficient for our model. E.g., the $o=3$ prediction only requires 30% more time to evaluate, despite generating three times as many predictions (details can be found in Sec. 6). Thus, the LSTM context successfully captures the state of the temporal latent space evolution, such that the model can predict future states almost as far as the given input history.

A comparison of a fully recurrent LSTM with our proposed hybrid alternative can be found in Fig. 6(g). In this scene, representative for our other test runs, the hybrid architecture outperforms the fully recurrent (FR) version in terms of accuracy, while using 8.9m fewer weights than the latter. The full network sizes are 19.5m weights for hybrid, and 28.4m for the FR network. We additionally evaluate a hybrid architecture with an additional conv. layer of size 4096 with tanh activation after the LSTM decoder layer (V2 in Fig. 6(g)). This variant yields similar error measurements to the original hybrid architecture. We found in general that additional layers did not significantly improve prediction quality in our setting. The FR version for the $128^{3}$ data set below requires 369.4m weights due to its increased latent space dimensionality, which turned out to be infeasible. Our hybrid variant has 64m weights, which is still a significant number, but yields accurate predictions and reasonable training times. Thus, in the following tests, a total pressure inference model with a hybrid LSTM architecture for $o=1$ will be used unless otherwise noted.

To clearly show the full data sets and their evolution over the course of a temporal prediction, we have trained a two-dimensional model, the details of which are given in App. B. In Fig. 7 sequences of the ground truth data are compared to the corresponding autoencoder baseline, and the outputs of our prediction network. Even though the autoencoder produces noise due to the strong compression, the temporal predictions closely match the autoencoder baseline, and the network is able to reproduce the complex behavior of the underlying simulations. E.g., the two waves forming on the right hand side of the domain in Fig. 7a indicate that the network successfully learned an abstraction of the temporal evolution of the flow. Further tests of the full prediction network with autoencoder models that utilize gradient losses to circumvent the visual noise, yielded no better prediction capabilities than the presented autoencoder with L2 loss. Additional 2D examples using the presented autoencoder can be found in Fig. 15 of App. B.

Results

We now apply our model to the additional data sets with higher spatial resolutions, and we will highlight the resulting performance in more detail. First, we demonstrate how our method performs on the liquid128 data set, with its eight times larger number of degrees of freedom per volume. Correspondingly, we use a latent space size of $m_{s}=8192$, and a prediction network with LSTM layers of size 1000 and 1500. Despite the additional complexity of this data set, our method successfully predicts the temporal evolution of the pressure fields, with an average PSNR of 44.8. The lower value compared to the $64^{3}$ case is most likely caused by the higher intricacy of the $128^{3}$ data. Fig. 9a) shows a more realistically rendered simulation for $i_{p}=4$. This setup contains a shape that was not part of any training data simulations. Our model successfully handles this new configuration, as well as other situations shown in the accompanying video. This indicates that our model generalizes to a broad class of physical behavior. To evaluate long term stability, we have additionally simulated a scene for 650 time steps which successfully comes to rest. This simulation and additional scenes can be found in the supplemental video.

A trained model for the smoke128 data set can be seen in Fig. 9b. Despite the significantly different physics, our approach successfully predicts the evolution and motion of the vortex structures. However, we noticed a tendency to underestimate pressure values, and to reduce small-scale motions. Thus, while our model successfully captures a significant part of the underlying physics, there is a clear room for improvement for this data set.

Our method also leads to significant speedups compared to regular pressure solvers, especially for larger volumes. For example the pressure inference by the prediction network for a $128^{3}$ volume takes 9.5ms, on average. Including the times to encode and decode the respective simulation fields of resolution $128^{3}$ (4.1ms and 3.3ms, respectively) this represents a 155$\times$ speedup compared to a parallelized state-of-the-art iterative MIC-CG pressure solver [Bri15], running with eight threads. While the latter yields a higher overall accuracy, and runs on a CPU instead of a GPU, it also represents a highly optimized numerical method. We believe the speedup of our LSTM version indicates a huge potential for very fast physics solvers with learned models.

It however, also leads to a degradation of accuracy compared to a regular iterative solver. The degradation can be controlled by chosing an appropiate prediction interval as described in Sec. 5.2 and can therefore be set according to the required accuracy. Even when taking into account a factor of ca. 10$\times$ for GPUs due to their better memory bandwidth, this leaves a speedup by a factor of more than 15$\times$, pointing to a significant gain in efficiency for our LSTM-based prediction. In addition, we measured the speedup for our (not particularly optimized) implementation, where we include data transfer to and from the GPU for each simulation step. This very simple implementation already yields practical speedups of 10x for an interval prediction with $i_{p}=14$. Details can be found in Table 5 and Table 6, while Table 2 summarizes the sizes of our data sets. All measurements were created with the tensorflow timeline tools on Intel i7 6700k (4GHz) and Nvidia Geforce GTX970.

While we have shown that our approach leads to large speed-ups and robust simulations for a significant variety of fluid scenes, there are several areas with room for improvements and follow up work. First, our LSTM at the moment strongly relies on the AE, which primarily encodes large scale scale dynamics, while small scale dynamics are integrated by the alignment of free surface boundary conditions [ATW15]. Also, our current, relatively simple AE can introduce a certain amount of noise in the solutions, which, however, can potentially be alleviated by different network architectures.

Overall, improving the AE network is important in order to improve the quality of the temporal predictions. Our experiments also show that larger data sets should directly translate into improved predictions. This is especially important for the latent space data set, which cannot be easily augmented.

Conclusions

With this work we arrive at three important conclusions: first, deep neural network architectures can successfully predict the temporal evolution of dense physical functions, second, learned latent spaces in conjunction with LSTM-CNN hybrids are highly suitable for this task, and third, they can yield very significant increases in simulation performance.

In this way, we arrive at a data-driven solver that yields practical speed-ups, and at its core is more than 150x faster than a regular pressure solve. We believe that our work represents an important first step towards deep-learning powered simulation algorithms. On the other hand, given the complexity of the problem at hand, our approach represents only a first step. There are numerous, highly interesting avenues for future research, ranging from improving the accuracy of the predictions, over performance considerations, to using such physics predictions as priors for inverse problems.

References

Appendix A Long-short Term Memory Units and Dimensionality

A central challenge for deep learning problems involving fluid flow is the large number of degrees of freedom present in three-dimensional data sets. This quickly leads to layers with large numbers of nodes – from hundreds to thousands per layer. Here, a potentially unexpected side effect of using LSTM nodes is the number of weights they require.

The local feedback loops for the gates of an LSTM unit all have trainable weights, and as such induce an $n\times n$ weight matrix for $n$ LSTM units. E.g., even for a simple network with a one dimensional input and output, and a single hidden layer of 1000 LSTM units, with only $2\times 1000$ connections and weights between in-, output and middle layer, the LSTM layer internally stores $1000^{2}$ weights for its temporal feedback loop. In practice, LSTM units have input, forget and output gates in addition to the feedback connections, leading to $4n^{2}$ internal weights for an LSTM layer of size $n$. Correspondingly, the number of weights of such a layer with $n_{o}$ nodes, i.e., outputs, and $n_{i}$ inputs is given by $n_{\text{lstm}}=4(n_{o}^{2}+n_{o}(n_{i}+1))$. In contrast, the number of weights for the 1D convolutions we propose in the main document is $n_{\text{conv-1d}}=n_{o}k(n_{i}+1)$, with a kernel size $k=1$.

Keeping the number of weights at a minimum is in general extremely important to prevent overfitting, reduce execution times, and to arrive at networks which are able to generalize. To prevent the number of weights from exploding due to large LSTM layers, we propose the mixed use of LSTM units and convolutions for our final temporal network architecture. Here, we change the decoder part of the network to consist of a single dense LSTM layer that generates a sequence of $o$ vectors of size $m_{t_{d}}$. Instead of processing these vectors with another dense LSTM layer as before, we concatenate the outputs into a single tensor, and employ a single one-dimensional convolution translating the intermediate vector dimension into the required $m_{s}$ dimension for the latent space. Thus, the 1D convolution works along the vector content, and is applied in the same way to all $o$ outputs. Unlike the dense LSTM layers, the 1D convolution does not have a quadratic weight footprint, and purely depends on the size of input and output vectors.

Appendix B Additional Results

In Fig. 11 additional time-steps of the comparison from Fig. 8 are shown. Here, different inferred simulation quantities can be compared over the course of a simulation for different models. In addition, Fig. 12, 13, and 14 show more realistic renderings of our liquid64, liquid128, and smoke128 models, respectively.

As our solve indirectly targets divergence, we also measured how well the predicted pressure fields enforce divergence freeness over time. As a baseline, the numerical solver led to a residual divergence of $3.1\cdot 10^{-3}$ on average. In contrast, the pressure field predicted by our LSTM on average introduced a $2.1\cdot 10^{-4}$ increase of divergence per time step. Thus, the per time step error is well below the accuracy of our reference solver, and especially in combination with the interval predictions, we did not notice any significant changes in mass conservation compared to the reference simulations.

To visualize the temporal prediction capabilities as depicted in Fig. 7, the spatial encoding task of the total pressure $p_{t}$ approach was reduced to two spatial dimensions. For this purpose a 2D autoencoder network was trained on a dataset of resolution $64^{2}$. The temporal prediction network was trained as described in the main document. Additional sequences of the ground truth, the autoencoder baseline, and the temporal prediction by the LSTM network are shown in Fig. 15.

Appendix C Fluid Simulation Setup

In addition to Sec. 4, we provide more information on the simulation setup in the following. To generate our liquid datasets we use a classic NS solver [Bri15]. The timestep is fixed to $0.1$, and pressure is computed with a conjugent gradient solver accuracy of $5\cdot 10^{-5}$. The external forces in our setup only consist of a gravity vector of $(0.0,-0.01,0.0)$ that is applied after every velocity advection step. No additional viscosity or surface tension forces are included.

In addition to the central quantities of a fluid solve, flow velocity $\mathbf{u}$, pressure $p$, and potentially visible quantities such as the levelset $\phi$, we utilize the Fluid Implicit Particle (FLIP) [ZB05] method, which represents a grid-particle hybrid. It is used in this work on the one hand to generate the liquid datasets and on the other to be the base of our neural network driven interval prediction simulation.

To give a general overview of how the simulation proceeds, we shortly describe the computations executed for every time step in the following. In each simulation step we first advect the FLIP particle system $PS$, the levelset $\phi$ and the velocity $\mathbf{u}$ itself with the current velocity grid. Afterwards a second levelset containing the particle surface is created based on the current $PS$ configuration and is merged with $\phi$. The merged levelset is extrapolated within a narrow band region of $3$, as described in the main text. After the levelset transformations $\mathbf{u}$ is updated with the $PS$ velocities and the external forces like gravity are applied on the result, followed by the enforcement of the static wall boundary conditions. Next, a pressure field $p$ is computed via a Poisson solve using the divergence of $\mathbf{u}$ as right hand side. After completing the pressure solve, the gradient of the result $\nabla p$ is subtracted from $\mathbf{u}$ yielding an approximation of a divergence free version of $\mathbf{u}$. The $PS$ velocities are updated based on the difference between post-advection version of $\mathbf{u}$ and the latest divergence free one.

The presented LSTM prediction framework supports predictions of different simulation fields from the FLIP simulation presented above. The supported fields are the final $\mathbf{u}$ at the end of the simulation loop, the solved pressure $p$ and the decomposed version of $p$ with $p_{s}$ and $p_{d}$, i.e. the hydrostatic and dynamic components of the regular pressure field, respectively. For the prediction of these fields we supply multiple architectures that are compared in the main text. Those are the total pressure, variational split pressure, split pressure and velocity versions. The difference between the variational split pressure and the default split pressure approach is the architecture of the autoencoder, whereas the temporal prediction network stays the same.

Depending on the prediction architecture, the inferred, decoded predicted field is used instead of executing the corresponding numerical approximation step. When targeting $\mathbf{u}$ with our method, this means that we can omit velocity advection as well as pressure solve, while the inference of $p$ by the split or total pressure architecture means that we only omit the pressure solve, but still need to perform advection and velocity correction with the pressure gradient. While the latter requires more computations, the pressure solve is typically the most time consuming part, with a super-linear complexity, and as such both options, using either $\mathbf{u}$ or the $p$ variants, have comparable runtimes.

Appendix D Hyperparameters

To find appropriate hyperparameters for the prediction network, a large number of training runs with varying parameters were executed on a subset of the total training data domain. The subset consisted of $100$ scenes of the training data set discussed in Sec. 4.2.

In Fig. 19 (a-d), examples for those searches are shown. Each circle shown in the graphs represents the final result of one complete training run with the parameters given on the axes. The color represents the mean absolute error of the training error, ranging from purple (the best) to yellow (the worst). The size of the circle corresponds to the validation error, i.e., the most important quantity we are interested in. The best two runs are highlighted with a dark coloring. These searches yield interesting results, e.g. Fig. 19(a) shows that the network performed best without any weight decay regularization applied.

Choosing good parameters leads to robust learning behavior in the training process, an example is shown in Fig. 20. Note that it is possible for the validation error to be lower than the training error as dropout is not used for computing the validation loss. The mean absolute error of the prediction on a test set of $40$ scenes, which was generated independently from the training and validation data, was $0.0201$ for this case. These results suggest that the network generalizes well with the given hyperparameters.