Latent-space Physics: Towards Learning the Temporal Evolution of Fluid Flow
Steffen Wiewel Moritz Becher Nils Thuerey




    
    < g r a p h i c s >

    Our method encodes multiple steps of a simulation field, typically pressure, into a reduced latent representation with a convolutional neural network. A
    second neural network with LSTM units then predicts the latent space code for one or more future time steps, yielding large reductions in runtime compared
    to regular solvers.
    


Guessing models and the approachability ideal
    [
    December 30, 2023
=============================================


    We propose a method for the data-driven inference of temporal evolutions 
    of physical functions with deep learning. More specifically, we target fluid flows, i.e. Navier-Stokes problems,
    and we propose a novel LSTM-based approach to predict the changes of pressure fields over time.
    The central challenge in this context is the high dimensionality of Eulerian space-time data sets.
    We demonstrate for the first time that dense 3D+time functions of physics system can be predicted
    within the latent spaces of neural networks, and we arrive at a neural-network based simulation
    algorithm with significant practical speed-ups. We highlight the capabilities of our method with
    a series of complex liquid simulations, and with a set of single-phase buoyancy simulations.
    With a set of trained networks, our method is more than two orders of magnitudes faster than a
    traditional pressure solver. Additionally, we present and discuss a series of detailed evaluations
    for the different components of our algorithm.



<ccs2012>
<concept>
<concept_id>10010147.10010257.10010293.10010294</concept_id>
<concept_desc>Computing methodologies Neural networks</concept_desc>
<concept_significance>500</concept_significance>
</concept>
<ccs2012>
<concept_id>10010147.10010371.10010352.10010379</concept_id>
<concept_desc>Computing methodologies Physical simulation</concept_desc>
<concept_significance>500</concept_significance>
</concept>
<concept>


[500]Computing methodologies Neural networks
[500]Computing methodologies Physical simulation
    
  





§ 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 <cit.>,
to perform divergence free projections for a single point in time <cit.>,
or for adversarial training of high resolution flows <cit.>,
few works exist that target temporal evolutions of physical systems.
While first works have considered predictions of Lagrangian objects
such as rigid bodies <cit.>, and control of two dimensional 
interactions <cit.>,
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 <cit.>,
the interest in neural network algorithms is a relatively
new phenomenon, triggered by seminal works such as ImageNet <cit.>.
In computer graphics, such approaches have led to impressive results, e.g.,
for synthesizing novel viewpoints of natural scenes <cit.>,
to generate photorealistic face textures <cit.>,
and to robustly transfer image styles between photographs <cit.>,
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,θ) ≈ y.
Here, θ 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 = σ( W_l a^l-1 + b_l ), where a^i is the output of the i'th
layer, σ 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 θ.

The latent spaces of generative NNs were shown to be powerful tools in image processing and synthesis
<cit.>.
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 
<cit.>, 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 <cit.>.
These recurrent networks are especially popular for 
control tasks in reinforcement learning environments <cit.>.
Recently, impressive results were also achieved for tasks such as automatic 
video captioning <cit.>.

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. battaglia2016interaction predict
two-dimensional physics, a method that was subsequently extended to videos <cit.>.
Another line of work proposed a specialized architecture for two-dimensional
rigid bodies physics <cit.>,
while others have targeted predictions of liquid motions for robotic control <cit.>
Farimani et al. <cit.> 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 <cit.> 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) <cit.>.
In the field of weather prediction, short-term prediction architectures evolved that internally also make use of RNN layers <cit.>.
Other works create velocity field predictions by learning parameters for Reynolds Averaged Navier Stokes models from simulation data <cit.>.
In addition, learned Koopman operators <cit.> 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 <cit.>.
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

    ∂/∂t  + ·∇ =
    	- 1/ρ∇p + ν∇^2 + g⃗ ,   ∇· = 0 ,

where the most important quantities are flow velocity u⃗ and pressure p. 
The other parameters ρ,ν,g⃗ denote density, kinematic viscosity and external forces, respectively.
For liquids, we will assume that a signed-distance function ϕ 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 <cit.>,
while Foster and Metaxas employed a first three-dimensional NS solver <cit.>.
After Jos Stam proposed an unconditionally stable advection and time integration scheme stam1999, 
it led to powerful liquid solvers based on the particle levelset <cit.>, 
in conjunction with accurate free surface boundary conditions <cit.>.
Since then, the fluid implicit particle (FLIP) method has been especially
popular for detailed liquid simulations <cit.>, and we will use it to generate our training data.
Solvers based on these algorithms have subsequently been extended 
with accurate boundary handling <cit.> synthetic turbulence <cit.>,
or narrow band algorithms <cit.>, to name just a few examples.
A good overview of fluid simulations for computer animation can be found in the 
book by R. Bridson bridson2015.
Aiming for more complex systems, other works have targeted
coupled reduced order models, or sub-grid coupling effects <cit.>.
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 <cit.>. 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 <cit.>, and fast solvers <cit.>
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. ladicky2015data.
Other graphics works have targeted learning flow descriptors with CNNs <cit.>, or
learning the statistics of splash formation <cit.>, and
two-dimensional control problems <cit.>.
While pressure projection algorithms with CNNs <cit.>
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 <cit.>. 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


The central goal of our method is to predict future states of a physical function x⃗. While previous works often consider
low dimensional Lagrangian states such as center of mass positions, x⃗ takes the form of a dense Eulerian function in 
our setting. E.g., it can represent a spatio-temporal pressure function, or a velocity field. Thus, we consider 
x⃗: ^3×^+→^d,
with d=1 for scalar functions such as pressure, and d=3 for vectorial functions such as velocity fields. As a consequence,
x⃗ has very high dimensionality when discretized, typically on the order of millions of spatial degrees of freedom.

Given a set of parameters θ and a functional representation f_t
our goal is to predict the o future states x⃗(t+h) to x⃗(t+oh) as closely as possible given a current state and a series of n previous states, i.e.,

    f_t ( x⃗(t-nh),...,x⃗(t-h), x⃗(t) ) ≈[x⃗(t+h), …, x⃗(t+oh)].

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

Due to the high dimensionality of x⃗, directly solving eq:ft_def would be very costly. Thus, we employ two additional functions f_d and f_e, that compute a low dimensional encoding. The encoder f_e maps into an m_s dimensional space ∈^m_s with f_e(x⃗(t))=^t, whereas the decoding function f_d reverses the mapping with f_d(^t)=x⃗(t).
Thus, f_d and f_e here represent spatial decoder and encoder models, respectively. 
Given such an en- and decoder, we rephrase the problem above as

    f̃_̃t̃( f_e( x⃗(t-nh) ),..., f_e( x⃗(t)) )    ≈[^t+h, …, ^t+oh] 
    
    		f_d( [^t+h, …, ^t+oh] )    ≈[x⃗(t+h), …, x⃗(t+oh)] 
    
    		f_d( f̃_̃t̃(f_e( x⃗(t-nh) ),..., f_e( x⃗(t))) )    ≈ f_t( x⃗(t-nh),..., x⃗(t) )



We will use CNNs for f_d and f_e, and thus the space  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.

Our prediction network that models the function f̃_̃t̃ will likewise employ an encoder-decoder structure, and turns the temporal stack of encoded data points f_e(x⃗) into a reduced representation , which we will refer to as temporal context below. 
The architectures of the spatial and temporal networks are described in the following sections sec:spacenet and sec:timenet, respectively.



 §.§ Reducing Dimensionality



In order to reduce the spatial dimensionality of our inference problem, we employ a fully convolutional autoencoder architecture <cit.>. 
Our autoencoder (AE) consists of the aforementioned encoding and decoding functions f_e,f_d and is trained to reconstruct the quantity x⃗ as accurately as possible w.r.t. an L_2 norm, i.e. 

    θ_d, θ_emin | f_d( f_e( x⃗(t) )) - x⃗(t) |_2,

where θ_d, θ_e represent the parameters of the decoder and encoder, respectively.
We use a series of convolutional layers activated by leaky rectified linear units (LeakyReLU) <cit.> 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 ∈ [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. <cit.>, 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

    θ_e_0...k, θ_d_0...kmin | f_d_0∘ f_d_1∘ ... ∘ f_d_k( f_e_k∘ f_e_k-1∘ ... ∘ f_e_0( x⃗(t) )) - x⃗(t) |_2 ,

with θ_e_0...k, θ_d_0...k denoting the parameters of the sub-stack for pretraining stage k, and ∘ 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:autoencoder_architecture, 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 <cit.>.
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 tab:autoencoder_architecture.

In addition to this basic architecture, we will also evaluate a variational autoencoder <cit.> in sec:eval 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.



 §.§ Prediction of Future States






The prediction network, that models f̃_̃t̃, transforms a sequence of n+1 chronological, encoded input states X = (^t-nh, …, ^t-h, ^t) 
into a consecutive list of o predicted future states Y = (^t+h, …, ^t+oh).
The minimization problem solved during training thus aims for minimizing
the mean absolute error between the o predicted and ground truth states
with an L_1 norm:

    θ_tmin| f̃_̃t̃( ^t-nh, …, ^t-h, ^t) - [ ^t+h, …, ^t+oh] |_1  
    		 .

Here θ_t denotes the parameters of the prediction network, and [·,·] denotes concatenation of the  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  varies smoothly in time, the order of the entries is arbitrary and we cannot make any assumptions about local neighborhoods within .
As such, convolving  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.

The prediction network approximates the desired function f̃_̃t̃ with the help of an internal temporal context of dimension m_t, which we will denote as .
Thus, the first part of our prediction network represents a recurrent encoder module, transforming n+1 latent space points into a temporal context .
The time decoder module has a similar structure, and is likewise realized with layers of LSTM units.
This module takes a context  as input, and outputs a series of future, spatial latent space representations.
By means of its internal state, the time decoder is trained to predict o future states when receiving the same context  repeatedly.  
We use tanh activations for all LSTM layers, and hard sigmoid functions for efficient, internal activations. 

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:seq-to-seq that is used instead of the fully recurrent approach presented in fig:fullyrecurrent.



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:eval.








§ FLUID FLOW DATA


To generate fluid data sets for training we rely on a NS solver with operator splitting <cit.>
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 u⃗, evaluating external forces, and then computing
the harmonic pressure function p. 
In addition, a visible, passive quantity such as smoke density ρ, or a level-set representation ϕ
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. 

In addition, we consider a split pressure, which represents a residual quantity. 
Assuming a fluid at rest on a ground with height z_g,
a hydrostatic pressure p_s for cell at height z,
can be calculated as
p_s(z) = p(z_0) + 1/A∫_z^z_0∬_A g⃗ρ(h)  dx  dy  dh,
with z_0, p_0, A denoting
surface height, surface pressure, and cell area, respectively. 
As density and gravity can be treated as constant in our setting, this further simplifies to 
p_s = ρg⃗ (z - z_0).
While this form can be evaluated very efficiently, it has the drawback that it is only
valid for fluids in hydrostatic equilibrium, and typically cannot be used for dynamic simulations in 3D.
Given a data-driven method to predict pressure fields, we can incorporate the hydrostatic pressure 
into a 3D liquid simulation by decomposing the regular pressure field into hydrostatic and dynamic 
components p_t = p_s + p_d, such that our autoencoder separately receives and encodes the 
two fields p_s and p_d.
To differentiate between the classic pressure field created by the simulation and our extracted split pressure components p_s and p_d, the classic pressure field is called total pressure p_t in the following.
With this split pressure, the autoencoder could potentially put more emphasis on the small scale 
fluctuations p_d from the hydrostatic pressure gradient.

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:data_fields_raw, 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 u⃗ and the two pressure variants, total p_t and split (p_s and p_d) in the following. 
When targeting 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 <cit.>.
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.



 §.§ Interval Prediction

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=∞ 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.



 §.§ 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 [-1,1].



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:dataset_renderings, fig:dataset_liq128initial and fig:dataset_smoke 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 <cit.>.
More specifically, given two signed distance functions ϕ_r, ϕ_p representing reference and 
predicted surfaces,	we compute the surface error as	

    e_h = max( 1/|S_p| ∑_p⃗_1 ∈ S_pϕ_r(p⃗_1), 1/|S_r| ∑_p⃗_2 ∈ S_rϕ_p(p⃗_2)) / Δ x.

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




 §.§ Spatial Encoding

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 x⃗ under consideration with = f_e(x⃗),
and then restore it from its latent space representation  x⃗' = f_d(). 
In the following, we will compare flow velocity , 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:encdec_hausdorff and 
fig:Mencdec_PSNR 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:Mencdec_PSNR. 
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:encdec_renderings.







 §.§ Temporal Prediction


Next, we evaluate reconstruction quality when including the temporal prediction network.
Thus, now a quantity x⃗' 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:seq-to-seq, respectively.
Hence, 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·10^-2 with a recurrent dropout of 0.385,
and a learning rate of 1.26·10^-4 with a decay factor of 3.34·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:basenetworkcomparisoni1000 and fig:Mbasenetworkcomparisoni1000_psnr, with a surface visualization in  fig:prediction_renderings.
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:encdec_hausdorff and <ref>, 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:basenetworkcomparisoni15 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 .
A graph comparing accuracy for 1, 3 and 5 steps of output can be found in fig:intervals_comparison_levelset_multiout.
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:results).
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:Malgcomparisoni1000_psnr.
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:Malgcomparisoni1000_psnr). 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 sec:additional_results.
In fig:data_fields_lstm_noj 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:data_fields_lstm_noja 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:data_fields_prediction of sec:additional_results.




§ 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:renderPreviewa) 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:renderPreviewb.
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× speedup compared to a parallelized state-of-the-art iterative MIC-CG pressure solver <cit.>, 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 subsec:temp_pred and can therefore be set according to the required accuracy.
Even when taking into account a factor of ca. 10× for GPUs due to their better memory bandwidth, this leaves a speedup by a factor of more than 15×, 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 tab:performance and tab:performance64, while tab:datasets 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.








 §.§ Limitations

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 <cit.>.
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.



eg-alpha-doi










 




 
  Supplemental Document for 
 










§ 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.

Therefore we briefly summarize the central equations for layers of LSTM units to gain understanding about the required weights.
Below, f,i,o,g,s will denote forget, input, output, update and result connections, respectively.
h denotes the result of the LSTM layer.
Subscripts denote time steps, while θ and b denote weight and bias.
Below, we assume tanh as output activation function.
The new state for time step t of an LSTM layer is then given by:

    f_t    = σ( θ_xf x_t + θ_hf h_t-1 + b_f ) 
    
    		i_t    = σ( θ_xi x_t + θ_hi h_t-1 + b_i ) 
    
    		o_t    = σ( θ_xo x_t + θ_ho h_t-1 + b_o ) 
    
    		g_t    = tanh( θ_xg x_t + θ_hg h_t-1 + b_g ) 
    
    		s_t    = f_t ⊙ s_t-1 + i_t ⊙ g_t 
    
    		h_t    = o_t ⊙tanh(s_t)


The local feedback loops for the gates of an LSTM unit all have trainable weights, and as such induce an n × 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 × 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 4 n^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_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_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. 






§ ADDITIONAL RESULTS




In fig:prediction_renderingsOld additional time-steps of the comparison from fig:prediction_renderings are shown. 
Here, different inferred simulation
quantities can be compared over the course of a simulation for different models.
In addition, fig:transparent_64,
<ref>, and <ref> 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·10^-3 on average. In contrast,
the pressure field predicted by our LSTM on average introduced a 2.1·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:data_fields_lstm_noj, 
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:data_fields_prediction.









§ FLUID SIMULATION SETUP



 

FLIP Simulation







In addition to sec:fluiddata, we provide more information on the simulation setup in the following.
To generate our liquid datasets we use a classic NS solver <cit.>. The timestep is fixed to 0.1, and pressure is computed with a conjugent gradient solver accuracy 
of 5·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 u⃗, pressure p,
and potentially visible quantities such as the levelset ϕ, we utilize 
the Fluid Implicit Particle (FLIP) <cit.> 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 ϕ and the velocity 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 ϕ.
The merged levelset is extrapolated within a narrow band region of 3, as described in the main text. 
After the levelset transformations 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 u⃗ as right hand side.
After completing the pressure solve, the gradient of the result ∇ p is subtracted from u⃗ yielding 
an approximation of a divergence free version of u⃗.
The PS velocities are updated based on the difference between post-advection version of u⃗ and the latest divergence free one.


 

Prediction Integration

The presented LSTM prediction framework supports predictions of different simulation fields from the FLIP simulation presented above.
The supported fields are the final 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 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 u⃗ or the p variants, have comparable runtimes.















§ 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:data_gen.



In fig:hyperparam (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:hyperparamreglstm 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:lstmtrainhist.
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.

