This article aims to give a broad overview of how neural networks, Fully Convolutional neural networks in specific, are able to learn fluid flow around an obstacle by learning from examples.

This series is divided into three parts.

Part 1: A data-driven approach to CFD

Part 2: Implementation details (this post)

Part 3: Results
 

In part 1, we gained a high-level overview of the data-driven approach to CFD and the steps that are needed to make it work. In this second part, we explore some technical details of essential steps of the data-driven approach to CFD. Concretely, we want to focus on two aspects:

  • The architectural details of the neural network
  • Accuracy measurements for the predictions of the Neural Network

Results

As mentioned in the previous post, the network has a Fully-Convolutional architecture. In order to understand how the network processes the input and learns from given samples, lets dive into the code:

import tensorflow as tf
import numpy as np

def int_shape(x):

  list = x.get_shape().as_list()
  return list

def _activation_summary(x):
  """Helper to create summaries for activations.
  Creates a summary that provides a histogram of activations.
  Creates a summary that measure the sparsity of activations.
  Args:
    x: Tensor
  Returns:
    nothing
  """
  tensor_name = x.op.name
  tf.summary.histogram(tensor_name + '/activations', x)
  tf.summary.scalar(tensor_name + '/sparsity', tf.nn.zero_fraction(x))

def _variable(name, shape, initializer):
  """Helper to create a Variable.
  Args:
    name: name of the variable
    shape: list of ints
    initializer: initializer for Variable
  Returns:
    Variable Tensor
  """
  # getting rid of stddev for xavier ## testing this for faster convergence
  var = tf.get_variable(name, shape, initializer=initializer)
  return var

def conv_layer(inputs, kernel_size, stride, num_features, idx, tf.nn.relu=None):
  with tf.variable_scope('{0}_conv'.format(idx)) as scope:
    input_channels = int(inputs.get_shape()[3])

    weights = _variable('weights', shape=[kernel_size,kernel_size,input_channels,num_features],initializer=tf.contrib.layers.xavier_initializer_conv2d())
    biases = _variable('biases',[num_features],initializer=tf.contrib.layers.xavier_initializer_conv2d())

    conv = tf.nn.conv2d(inputs, weights, strides=[1, stride, stride, 1], padding='SAME')
    conv_biased = tf.nn.bias_add(conv, biases)
    if tf.nn.relu is not None:
      conv_biased = tf.nn.relu(conv_biased)
    return conv_biased

def transpose_conv_layer(inputs, kernel_size, stride, num_features, idx, tf.nn.relu=None):
  with tf.variable_scope('{0}_trans_conv'.format(idx)) as scope:
    input_channels = int(inputs.get_shape()[3])
    
    weights = _variable('weights', shape=[kernel_size,kernel_size,num_features,input_channels],initializer=tf.contrib.layers.xavier_initializer_conv2d())
    biases = _variable('biases',[num_features],initializer=tf.contrib.layers.xavier_initializer_conv2d())
    batch_size = tf.shape(inputs)[0]
    output_shape = tf.stack([tf.shape(inputs)[0], tf.shape(inputs)[1]*stride, tf.shape(inputs)[2]*stride, num_features]) 
    conv = tf.nn.conv2d_transpose(inputs, weights, output_shape, strides=[1,stride,stride,1], padding='SAME')
    conv_biased = tf.nn.bias_add(conv, biases)
    if tf.nn.relu is not None:
      conv_biased = tf.nn.relu(conv_biased)
    shape = int_shape(inputs)
    conv_biased = tf.reshape(conv_biased, [-1, shape[1]*stride, shape[2]*stride, num_features]) # new
    return conv_biased

def fc_layer(inputs, hiddens, idx, tf.nn.relu=None, flat = False):

  with tf.variable_scope('{0}_fc'.format(idx)) as scope:
    input_shape = inputs.get_shape().as_list()
    if flat:
      dim = input_shape[1]*input_shape[2]*input_shape[3]
      inputs_processed = tf.reshape(inputs, [-1,dim])
    else:
      dim = input_shape[1]
      inputs_processed = inputs
    
    weights = _variable('weights', shape=[dim,hiddens],initializer=tf.contrib.layers.xavier_initializer())
    biases = _variable('biases', [hiddens], initializer=tf.contrib.layers.xavier_initializer())
    output_biased = tf.add(tf.matmul(inputs_processed,weights),biases,name=str(idx)+'_fc')
    if tf.nn.relu is not None:
      output_biased = tf.nn.relu(ouput_biased)

    return output_biased

def nin(x, num_units, idx):
    """ a network in network layer (1x1 CONV) """
    s = int_shape(x)
    s[0] = 5
    x = tf.reshape(x, [np.prod(s[:-1]), s[-1]])
    x = fc_layer(x, num_units, idx)
    return tf.reshape(x, s[:-1]+[num_units])

def res_block(x, a=None, filter_size=16, keep_p=1.0, stride=1, gated=False, name="resnet"):
  orig_x = x
  orig_x_int_shape = int_shape(x)
  if orig_x_int_shape[3] == 1:
    x_1 = conv_layer(x, 3, stride, filter_size, name + '_conv_1')
  else:
    x_1 = conv_layer(tf.nn.relu(x), 3, stride, filter_size, name + '_conv_1')
  if a is not None:
    shape_a = int_shape(a)
    shape_x_1 = int_shape(x_1)
    a = tf.pad(
      a, [[0, 0], [0, shape_x_1[1]-shape_a[1]], [0, shape_x_1[2]-shape_a[2]],
      [0, 0]])
    x_1 += nin(tf.nn.relu(a), filter_size, name + '_nin')
  x_1 = tf.nn.relu(x_1)
  if keep_p < 1.0:
    x_1 = tf.nn.dropout(x_1, keep_prob=keep_p)
  if not gated:
    x_2 = conv_layer(x_1, 3, 1, filter_size, name + '_conv_2')
  else:
    x_2 = conv_layer(x_1, 3, 1, filter_size*2, name + '_conv_2')
    x_2_1, x_2_2 = tf.split(axis=3, num_or_size_splits=2, value=x_2)
    x_2 = x_2_1 * tf.nn.sigmoid(x_2_2)

  if int(orig_x.get_shape()[2]) > int(x_2.get_shape()[2]):
    assert(int(orig_x.get_shape()[2]) == 2*int(x_2.get_shape()[2]), "res net block only supports stride 2")
    orig_x = tf.nn.avg_pool(orig_x, [1,2,2,1], [1,2,2,1], padding='SAME')

  # pad it
  out_filter = filter_size
  in_filter = int(orig_x.get_shape()[3])
  if out_filter != in_filter:
    orig_x = tf.pad(
        orig_x, [[0, 0], [0, 0], [0, 0],
        [(out_filter-in_filter), 0]])
  return orig_x + x_2

def conv_res(inputs, tf.nn.relu_name, nr_res_blocks, filter_size, keep_prob=1.0, gated=True):
  """Builds conv part of net.
  Args:
    inputs: input images
    keep_prob: dropout layer
  """
  # store for as
  a = []
  # res_1
  x = inputs
  for i in range(nr_res_blocks):
    x = res_block(x, filter_size=filter_size, tf.nn.relu=tf.nn.relu, keep_p=keep_prob, gated=gated, name="resnet_1_" + str(i))
  
  # res_2
  a.append(x)
  filter_size = 2 * filter_size
  x = res_block(x, filter_size=filter_size, tf.nn.relu=tf.nn.relu, keep_p=keep_prob, stride=2, gated=gated, name="resnet_2_downsample")
  for i in range(nr_res_blocks):
    x = res_block(x, filter_size=filter_size, tf.nn.relu=tf.nn.relu, keep_p=keep_prob, gated=gated, name="resnet_2_" + str(i))
  
  # res_3
  a.append(x)
  filter_size = 2 * filter_size
  x = res_block(x, filter_size=filter_size, tf.nn.relu=tf.nn.relu, keep_p=keep_prob, stride=2, gated=gated, name="resnet_3_downsample")
  for i in range(nr_res_blocks):
    x = res_block(x, filter_size=filter_size, tf.nn.relu=tf.nn.relu, keep_p=keep_prob, gated=gated, name="resnet_3_" + str(i))
  
  # res_4
  a.append(x)
  filter_size = 2 * filter_size
  x = res_block(x, filter_size=filter_size, tf.nn.relu=tf.nn.relu, keep_p=keep_prob, stride=2, gated=gated, name="resnet_4_downsample")
  for i in range(nr_res_blocks):
    x = res_block(x, filter_size=filter_size, tf.nn.relu=tf.nn.relu, keep_p=keep_prob, gated=gated, name="resnet_4_" + str(i))
  
  # res_5
  a.append(x)
  filter_size = 2 * filter_size
  x = res_block(x, filter_size=filter_size, tf.nn.relu=tf.nn.relu, keep_p=keep_prob, stride=2, gated=gated, name="resnet_5_downsample")
  for i in range(nr_res_blocks):
    x = res_block(x, filter_size=filter_size, tf.nn.relu=tf.nn.relu, keep_p=keep_prob, gated=gated, name="resnet_5_" + str(i))
  

  # res_up_1
  filter_size = filter_size // 2
  x = transpose_conv_layer(x, 3, 2, filter_size, "up_conv_1")
  #x = PS(x,2,512)
  for i in range(nr_res_blocks):
    if i == 0:
      x = res_block(x, a=a[-1], filter_size=filter_size, tf.nn.relu=tf.nn.relu, keep_p=keep_prob, gated=gated, name="resnet_up_1_" + str(i))
    else:
      x = res_block(x, filter_size=filter_size, tf.nn.relu=tf.nn.relu, keep_p=keep_prob, gated=gated, name="resnet_up_1_" + str(i))
  

  # res_up_2
  filter_size = filter_size // 2
  x = transpose_conv_layer(x, 3, 2, filter_size, "up_conv_2")
  #x = PS(x,2,512)
  for i in range(nr_res_blocks):
    if i == 0:
      x = res_block(x, a=a[-2], filter_size=filter_size, tf.nn.relu=tf.nn.relu, keep_p=keep_prob, gated=gated, name="resnet_up_2_" + str(i))
    else:
      x = res_block(x, filter_size=filter_size, tf.nn.relu=tf.nn.relu, keep_p=keep_prob, gated=gated, name="resnet_up_2_" + str(i))

  # res_up_3
  filter_size = filter_size // 2
  x = transpose_conv_layer(x, 3, 2, filter_size, "up_conv_3")
  #x = PS(x,2,512)
  for i in range(nr_res_blocks):
    if i == 0:
      x = res_block(x, a=a[-3], filter_size=filter_size, tf.nn.relu=tf.nn.relu, keep_p=keep_prob, gated=gated, name="resnet_up_3_" + str(i))
    else:
      x = res_block(x, filter_size=filter_size, tf.nn.relu=tf.nn.relu, keep_p=keep_prob, gated=gated, name="resnet_up_3_" + str(i))
 
  # res_up_4
  filter_size = filter_size // 2
  x = transpose_conv_layer(x, 3, 2, filter_size, "up_conv_4")
  #x = PS(x,2,512)
  for i in range(nr_res_blocks):
    if i == 0:
      x = res_block(x, a=a[-4], filter_size=filter_size, tf.nn.relu=tf.nn.relu, keep_p=keep_prob, gated=gated, name="resnet_up_4_" + str(i))
    else:
      x = res_block(x, filter_size=filter_size, tf.nn.relu=tf.nn.relu, keep_p=keep_prob, gated=gated, name="resnet_up_4_" + str(i))
  
  # res_up_2
  x = conv_layer(x, 3, 1, 2, "last_conv")
  x = tf.nn.tanh(x) 
  return x

The network is build in function conv_res in line 129. Here, first, the residual blocks are built on top of each other while scaling the image down and adding feature maps. Later, in the deconvolutional part of the network, the image is scaled up again, assigning two floating point values to each of the input pixels, representing the flow velocity in x- and y-directions.

Prediction accuracy estimation

In order to be able to make useful measurements of the validity of the prediction of a learning algorithm, one has to define metrics that describe the level of accuracy that the learning algorithms offer. The loss of a neural network is the variable that is being minimized by the optimizer during training and describes, how much of an error the neural net produces. Different measurements of how accurate the outputs of the neural network are needed to express the validity of the predictions.

CFD simulations are typically performed with a specific goal in mind. In many applications, a certain metric such as maximum flow velocity, pressure drop, or vorticity is being calculated and optimized. The output of the data-driven approach to CFD presented in this work does merely predicts the flow velocity field for each pixel or voxel for the pixel/voxel grid that is used as an input to the neural network.

For scalar values, a relative error between an estimated value v and the actual value \(v_0\) is defined as

\[\begin{align*} e_{rel} = \frac{|v-v_0|}{v} \end{align*}\]

Finally, it is necessary to find a suitable way to map the relative error \(e_rel\) of a prediction to the accuracy a of the neural network. If the relative error of the prediction is 0%, the accuracy should be 1.0. In contrast, if the relative error equal or greater than 100%, the accuracy should be close to zero. A Gaussian distribution for this mapping is proposed as follows:

\[\begin{align*} a(e_{rel}; \sigma^2) = e^{-\frac{e_{rel}^2}{2 \sigma^2}} \end{align*}\]

t is proposed that a relative error of 10% should lead to a prediction accuracy of 90%. Solving the equation for the variable sigma yields \(\sigma^2 = 0.048\).

In this work, four metrics describing the accuracy of the neural net output were examined:

Fraction of correct pixels

The number of correctly predicted pixels in an image gives an intuitive metric for how well the neural network predicts fluid flow behavior. As the network will never be able to predict the fluid flow velocity down to the last digit of a floating point number, the following approach is proposed: If the absolute error between the network prediction and the actual flow velocity is smaller than \(tv\), the respective pixel is declared as predicted correctly. The variable \(v\) denotes the average flow velocity in the simulation domain, and \(t\) denotes a threshold. A value of \(t = 0.01\) is suggested. The pixel-based accuracy of the network prediction can thus be defined as:

\[\begin{align*} a_{pixel} = \frac{N_{correctpixels}}{N_{totalpixels}} \end{align*}\]
  • Mass flow rate

Another important validation metric for an internal flow simulation is the mass flow rate through the cavity. Denoting \(\dot{m}\) as the mass flow rate, the relative error of the prediction of the mass flow rate can be defined as

\[\begin{align*} e_{rel, massflowrate} = \frac{\dot{m}_{simulation} - \dot{m}_{prediction}}{\dot{m}_{simulation}} \end{align*}\]

This leads to the following definition of the accuracy in terms of mass flow rate:

\[\begin{align*} a_{massflowrate} = a(e_{rel, massflowrate}; \sigma^2) \end{align*}\]

Drag value

In many external fluid flow simulation use-cases, the drag of a body in a fluid is an essential value to obtain, which is equal to the force the fluid acts on the body. The lower the drag, the lesser the resistance of the fluid. Therefore, the difference in drag was chosen as an additional metric for the success of the neural network predictions.

An acceptable approximation for the drag can be obtained with the following assumptions:

  • Steady flow
  • Incompressible flow
  • Uniform static pressure

Applying momentum conservation in the control volume \(\Omega\), the following momentum balance is obtained:

\[\begin{align*} F_{ext} + \int \int_{\partial \Omega} \rho v_{i}^2 dx dy = 0, i = {1,2} \end{align*}\]

while \(\partial \Omega\) denotes the boundary of the simulation domain and \(F_ext\) denotes an external force on the control volume identical to the drag force. Neglecting the velocity component in y-direction and assuming that the flow through the upper and lower boundaries due to the solid wall boundary conditions is negligible, the equation can be simplified to

\[\begin{align*} F_{ext} = \int_A^B \rho U_{infty}^2 dx - \int_C^D \rho v_x^2 dx \end{align*}\]

The velocity at the inlet is denoted as \(U_{infty}\). This derivation does not adhere to mass conservation as it neglects the flow through the upper and lower boundaries, though it holds a simple yet accurate estimate of the drag. Thus, the relative error of the prediction of the drag value is defined as

\[\begin{align*} e_{rel, Drag} = \frac{Drag_{simulation} - Drag_{prediction}}{Drag_{simulation}} \end{align*}\]

and the prediction accuracy may finally be defined as:

\[\begin{align*} a_{drag} = a(e_{rel, drag}; \sigma^2) \end{align*}\]

Conservation of mass

Numerical CFD solvers aim to find a solution to the continuity equation and the momentum equation. For an incompressible fluid, the continuity equation dictates that the divergence of the velocity vector field is zero for every point in the simulation domain. This follows the intuition that at no point in the simulation domain fluid springs into existence (divergence would be greater than zero) or ceases to exist (divergence would be smaller than zero).

By design, the Finite Volume Method, which is used for the SIMPLE algorithm, preserves this property of the fluid even in a discretized form. A data-driven approach should as well obey this rule. In order to check the physical consistency of the predicted flow field, a discretized form of the divergence operator on the predicted output images is proposed. If \(u_{i,j,k}\) denotes the \(j\)-th and \(k\)-th velocity component in \(i\)-direction, a discrete divergence operator \(div^h\) may be defined as:

\[\begin{align*} div^h (u_{i,j,k}) = \frac{1}{2} (u_{i,j+1,k} - u_{i,j-1,k} + u_{i,j,k+1} - u_{i,j,k-1} ) \end{align*}\]

This equation is obtained with a central differences approach to approximate the partial derivatives in all dimensions of space. Note that the so-defined divergence of a tensor of shape [2,M,N] has the shape [M-2,N-2] as the central difference for the outermost pixel is not defined.

In order to define an overall divergence accuracy, the discrete divergence at each grid point is accumulated over the whole grid with:

\[\begin{align*} div = \frac{1}{NM} \sum_i^N \sum_j^M div_{ij} \end{align*}\]

The prefactors in front of the summation term norm the accumulated divergence to the size of the grid in order to define an accuracy independent of grid sizes.

Finally, a mapping of the normed divergence accumulation to an accuracy estimation must be defined. In this work, a value of sigma is proposed for which a normed divergence of \(10^-3\), an accuracy of 90% is reached. Which leads to a value of \(\sigma = 0.00222\).

Thus, an estimate of the physical validity of the prediction can be obtained with the accuracy measurement:

\[\begin{align*} a_{div} = a(e_{div}; \sigma^2) \end{align*}\]

This marks the end of our dive into the mechanics of predicting steady-state fluid flows. I hope you enjoyed this visit!

The next chapter will illustrate the results obtained with this data-driven approach to CFD. Be prepared for over 9000 images!