视频游戏的流体模拟（第 17 部分）

Michael J. Gourlay 博士

0/5 (0投票)

2014年5月28日

CPOL

17分钟阅读

18132

本文是本系列文章的第 17 篇，介绍了如何识别流体表面。您可以使用此信息来渲染表面或帮助计算表面张力。

流体表面识别

液体是具有自由表面（即形状不由其容器定义的表面）的流体。本文是本系列文章的第 17 篇，介绍了如何识别流体表面。您可以使用此信息来渲染表面或帮助计算表面张力。

第 1 部分总结了流体动力学；第 2 部分概述了流体模拟技术。第 3 部分和第 4 部分介绍了具有实时运行的双向流体-物体交互的涡粒流体模拟。第 5 部分对该模拟代码进行了剖析和优化。第 6 部分介绍了从涡量计算速度的微分方法，而第 7 部分展示了如何将流体模拟集成到典型的粒子系统中。第 8 部分解释了基于涡的方法的流体模拟如何处理流体中可变密度的问题；第 9 部分描述了如何近似浸入不同密度流体中的物体上的浮力和重力。第 10 部分介绍了密度如何随温度变化，热量如何在流体中传递，以及热量如何在物体和流体之间传递。第 11 部分增加了燃烧，一种产生热量的化学反应。第 12 部分解释了不正确的采样如何导致不必要的抖动运动，并描述了如何缓解它。第 13 部分增加了凸多面体和类升力。第 14 部分修改了这些多面体以表示基本容器。第 15 部分描述了一种基本的平滑粒子流体动力学 (SPH) 流体模拟，用于模拟容器中的流体。第 16 部分探讨了一种将 SPH 公式应用于涡粒方法的方法。

指示函数

流体模拟跟踪整个区域内的密度分布。如果您将密度分布量化为只有两个值（0 和 1），它将指示流体的位置和非位置。这样的函数称为指示函数，如图 1 所示。流体的表面位于指示函数为 0 和 1 之间的边界。如果您识别了这些过渡，您就识别了流体表面。

沿射线穿过物体的指示函数

图 1. 指示函数

一个问题：指示函数不利于跟踪移动的表面。假设您使用一些明确定义、无噪声的密度分布来识别流体表面。之后，您可以像平流被动示踪剂或涡子一样平流表面本身。然而，会产生一个问题：您将如何表示表面？您可以将其表示为一组多边形，然后将顶点平流为被动示踪粒子，但这会导致问题。例如，如果模拟开始时有两个流体团块，后来合并了（图 2）？或者一个团块分裂了？此类事件具有拓扑变化——流体表面连续性的变化。以这种方式跟踪表面将是极具问题的。

图 2. 流体表面的拓扑变化。两个团块合并或一个团块分裂。

尽管使用指示函数跟踪表面运动会产生更多问题而非解决问题，但至少在初始阶段，以这种方式识别表面仍然有用。稍后，我将提供一种更好的跟踪表面运动的解决方案，但目前让我们看看如何使用指示函数来识别表面。

使用指示函数查找表面

您正在跟踪流体表面这一事实意味着流体至少有两个组成部分——例如，空气和水。在理想的不可压缩流体中，每个组成部分在其内部都具有均匀的密度。我们将水内部的密度称为 ρ_inside，将水外部的密度称为 ρ_outside。

注意：请记住第 3 部分，模拟使用涡子来携带密度，但流体在没有涡子的地方也有密度——也就是说，在没有涡子的情况下，流体具有环境密度。此外，由于数值误差，涡子可能会收敛和发散，因此两个组成部分的密度都不是完全均匀的。这是一个小细节，但请记住，因为它会轻微影响算法。

在均匀网格上识别表面更容易，因此首先将密度传输到网格。第 3 部分描述了这样一个算法，而新的代码重用了相同的密度网格。

基本思想很简单，但其实现变得三维复杂。因此，首先考虑一维 (1D)，概念更容易理解，然后逐个案例向更高维度发展。每个案例都有一个单元格的概念，它是网格点的一个邻域。我们想在每个单元格内找到表面交叉点。

每个网格点都有一个相关的密度值。该算法查找密度在内部值和外部值之间转换的位置。在 1D 中，一个单元格只是两个相邻网格点之间的线段，表面是该线段上的一个点。图 3 提供了一个示例。

图 3. 查找外部值和内部值之间的转换。P 和 Q 是具有外部密度值的样本。A 是一个样本，其值为 Z < A < M。B 是一个样本，其值为 Z < B < M。S 值是指示函数从外部过渡到内部的位置的估计。M 值是最大密度值的样本。x 是样本的位置；h 是样本之间的间隔。表面位于每个 z 处，即指示函数过渡到外部值的点。假设表面发生的 Q 和 A 之间的 x 值与 A/M 成线性比例。如果 I(x_i)| < I(x_i+1)，则 z = x_i+ h_i(1 - A/M)。如果 I(x_i) > I(x_i+1)，则 z = x_i+ h_i(A/M)。

这个想法可以识别表面，但您如何以数据形式表示该信息？换句话说，考虑到您可以计算单元格内的表面交叉点，您应该如何存储该信息？下一节将解释一种方法。

符号距离函数

符号距离函数 (SDF) 描述了区域中每个点到最近表面的距离。图 4 显示了一个示例。

沿射线穿过物体的 SDF

图 4. 符号距离函数

对于包含表面交叉点的单元格周围的每个网格点，您可以存储该点到表面交叉点的距离。换句话说，您可以在每个相邻于表面交叉点的网格点上编码到最近表面的 SDF。在 1D 中，您可以使用密度作为指示符来查找表面交叉点，并为其到每个相邻网格点的距离进行编码。

在二维 (2D) 中，这个想法变得更加复杂，实现也更加复杂。在 2D 中，单元格是四边形。在 2D 中，您必须考虑四个网格点，而不是像 1D 中那样考虑每对网格点之间的线段。您仍然从查找相邻网格点之间的表面交叉点开始，就像在 1D 情况一样。每个单元格的边界上可以有零到四个表面交叉点，每种配置会导致零、一个或两个终止于表面交叉点的线段——也就是说，2D 单元格内的表面是一条线段或两条线段。图 5 提供了一个示例。

图 5. 查找 2D 中的符号距离。每个网格点的符号距离取决于到线段的最短距离。在某些情况下，这恰好是到单元格边界上最近表面交叉点的距离。在其他情况下，这恰好是到穿过单元格边界上两个表面交叉点的直线的距离。

在三维中，您必须考虑盒子状单元格而不是四边形，情况会变得更加复杂。每个盒子有 12 条边，每条边都可能有表面交叉点。每个表面交叉点都是一个与表面重合的多边形的顶点。此外，与 1D 或 2D 不同，即使您知道所有表面交叉点，如果不查看相邻单元格，您也无法确定多边形的拓扑。

所有这些问题都有解决方案，但为了简洁起见，暂时将细节搁置。给定一个密度值单元格（和相邻单元格），您可以找到单元格内的表面。

给定一个告诉我们表面在单元格内位置的算法，您应该如何表示该信息？您可以生成多边形网格，如上所述，但这在跟踪表面移动时会导致问题。如果您改用符号距离，跟踪会更容易。（稍后我将解释如何以及为什么。）

一个像下面这样的算法会将符号距离值分配给每个相邻于表面交叉点的网格点

鉴于
- 均匀的密度值网格；
- 感兴趣的流体外部的密度，ρ_outside；以及
- 感兴趣的流体内部的密度，ρ_inside。
Do
- 计算 ρ_denominator= ρ_inside- ρ_outside；
- 创建一个与密度网格具有相同几何形状的浮点 SDF 值均匀网格；
- 将 SDF 均匀网格初始化为所有 FLT_MAX 值；
- 创建一个与密度网格具有相同几何形状的布尔 SDF-pinned 值均匀网格；
- 将 SDF-pinned 值初始化为 false；以及
- 对于每个网格点…
- 对于 {x, y, z} 中的每个相邻网格点…
  - 找到其密度，ρ_adjacent。
  - 如果当前和相邻的网格点位于表面的两侧…
    - 找到到表面的距离。
    - 有条件地更新网格点的 SDF。

有关详细信息，请参阅附带代码中的 ComputeImmediateSDFFromDensity_Nearest。

SDF-pinned 网格仅指示哪些网格点的 SDF 具有值。这可能看起来是多余的，但稍后会派上用场。目前请忽略它。

运行此算法将得到一个 SDF 值均匀网格，但仅限于表面交叉点附近。其余区域有效具有未分配的值。

在跟踪表面运动时，了解离表面更远的 SDF 值很有用。下一节将介绍在给定一些已建立的 SDF 值（固定值）的情况下如何填充剩余的 SDF 值。

Eikonal 方程

您想要填充 SDF 值网格，前提是只有其中一些已被分配。这个问题恰好与求解受边界条件约束的偏微分方程相同。边界条件当然是 SDF 在某些（固定的）位置是已知的（来自上述算法）。感兴趣的微分方程表示，任意两点之间的 SDF 差值 u，其大小等于这两点之间的距离

这是 **Eikonal 方程** 的一个特例。Eikonal 这个词源自与 icon 相同的希腊词根，意为图像，因为它描述了光线的传播。求解 Eikonal 方程背后的思想是沿表面法线投射光线，并在沿每条光线行进时分配 SDF 值。该过程可行，但不是最有效的解决方案。

您还可以通过访问网格中的每个点并从相邻网格点传播 SDF 值来求解此方程（参见图 6）。

图 6. 在相邻网格点之间传播 SDF 值。粗黑线描绘了表面。箭头描绘了 SDF 值从表面传播到其余区域。请注意，此图中的 SDF 值被量化为线段。实际的 SDF 会更平滑——例如，如果网格具有无限分辨率，它将具有圆角。

粗略地说，这种方法的伪代码看起来很简单

对于每个网格点
- 使用 Eikonal 方程求解有限差分。

该算法类似于第 6 部分中用于求解泊松方程的 Gauss-Seidel 算法。但此算法应用于 Eikonal 方程时存在一些问题。与求解泊松方程一样，在每个网格点求解 Eikonal 方程需要已经拥有有限差分方案使用的所有相邻网格点的解。但与求解泊松方程不同，Eikonal 方程的解仅沿一个方向传播：顺风。信息从源传播，源是具有已求解、固定值的区域。对于 Eikonal 方程，有限差分模板仅使用上风值——即来自信息传播相反方向的值。

请记住，网格开始时有一些网格点的 SDF 是已知的（并且是固定的）。例如，假设初始已知值位于域的内部。如果求解器算法从第零个（例如，左上角）网格点开始，并沿索引值增加的方向（例如，向下和向右）前进，那么它将没有上风邻居，因此无法为该网格点分配值。如果相反，求解器从最后一个网格点开始，并沿索引值减小的方向前进，那么最后一个网格点将没有上风邻居，因此无法为该网格点分配值。因此，解决方案需要多次遍历网格，一次沿正增量，一次沿负增量，针对 {x,y,z} 中的每个轴。这总共需要对网格进行八次遍历。

{-,-,-}
{+,-,-}
{-,+,-}
{+,+,-}
{-,-,+}
{+,-,+}
{-,+,+}
{+,+,+}

ComputeRemainingSDFFromImmediateSDF_Nearest_Block 例程在附带的代码中实现了此算法，但请先阅读下一节以了解如何使用该例程。

使用波前并行化 Eikonal 方程求解器

如何并行化 Eikonal 求解器？因为它类似于 Gauss-Seidel 泊松求解器，您可能会想采用第 6 部分描述的相同方法：分割域并将线程分配给每个切片。与泊松求解器一样，每个网格点的解取决于相邻网格点的解，因此可能会想到红黑排序。但 Eikonal 方程的解严格地从源向外传播，因此解需要尊重该传播。

将 SDF 想象成由微风携带的一种物质，微风从源（表面）流出。在计算网格点的 SDF 之前，必须先计算其前驱的 SDF——那些上风的网格点。要并行化此算法，求解器必须有一种方法将作业添加到队列中，其中这些作业对应于其上风前驱已具有解的网格点。只有在解决了某个区域的所有输入后，您才能将解决该区域的作业添加到队列中。

Intel® Threading Building Blocks (Intel® TBB) 提供了 parallel_do，这有助于解决此类问题。parallel_do 的魔力在于 **feeders**（供给器）。工作例程会收到一块工作和一个供给器。供给器允许工作者添加工作块。因此，当工作者完成一块工作时，它可以确定该完成块是否会解除对某些未来块的阻塞，如果解除，则将该块供给 Intel TBB。之后，Intel TBB 将调度新块；当作业队列为空时，parallel_do 停止运行。

您可以将工作块等于网格点，但 parallel_do 机制会产生一些开销。因此，最好将多个网格点聚合到更大的块中。假设，最小的块是网格点，最大的块是整个域。较小的块导致更精细的粒度和更好的并行性，但开销更高。较大的块具有更粗的粒度和更差的并行性，但开销较低。要选择最佳块大小，请分析不同块大小的问题，以查看哪个大小导致最快的运行时间。

区块

本文随附的代码使用 Block 结构来表示一块工作。一个块代表均匀网格内的区域。每个块都有一个索引三元组：其 x、y 和 z 索引。

struct Block
{
    int mBlockIndices[ 3 ] ;    // Indices of this block.

Block 具有辅助例程，可在给定块大小时查找域开始和结束块，该域已被细分。

    /** Compute index of first voxel in block, for given axis.
    */
    int BeginR( size_t axis , const int voxelsPerBlock[ 3 ] , const int numGridPoints[ 3 ] , const int inc[ 3 ] ) const
    {
        // When running backwards, beginning index is at maximal corner of block.
        const int shift         = ( inc[ axis ] > 0 ) ? 0 : 1 ;
        // First compute unconstrained index -- as if all blocks are fully occupied.
        const int unconstrained = ( mBlockIndices[ axis ] + shift ) * voxelsPerBlock[ axis ] ;
        // Finally, constrain index to account for blocks that are not fully occupied by voxels.
        // That can happen when numGridPoints is not an integer multiple of voxelsPerBlock.
        return Min2( unconstrained , numGridPoints[ axis ] - 1 ) ;
    }

    /** Compute index of one-past-last voxel in block, for given axis.
    */
    int EndR( size_t axis , const int voxelsPerBlock[ 3 ] , const int numGridPoints[ 3 ] , const int inc[ 3 ] ) const
    {
        const int begin = BeginR( axis , voxelsPerBlock , numGridPoints , inc ) ;
        if( inc[ axis ] > 0 )
        {   // Running forward.  End is (one past) maximal corner of block.  Here, "past" means "after".
            const int unconstrained = begin + voxelsPerBlock[ axis ] ;
            return Min2( unconstrained , numGridPoints[ axis ] ) ;
        }
        else
        {   // Running backward.  End is (one past) minimal corner of block.  Here, "past" means "before".
            const int unconstrained = begin - voxelsPerBlock[ axis ] ;
            return Max2( unconstrained , -1 ) ;
        }
    }
} ;

Worker Functor（工作函数对象）

按照 Intel TBB 的惯例，创建一个函数对象来执行工作

/** Functor (function object) to compute SDF for remainder of domain, using Threading Building Blocks.
*/
class ComputeRemainingSDFFromImmediateSDF_Nearest_TBB
{
    UniformGrid< float > &      mSignedDistanceGrid         ;   /// Grid of SDF values.
    const UniformGrid< int > &  mSdfPinnedGrid              ;   /// Grid of flags indicating which SDF values are "pinned"
    int                         mInc[ 3 ]                   ;   /// Processing loop increment for each axis.  +1 fwd, -1 bkwd
    int                         mVoxelsPerBlock[ 3 ]        ;   /// Number of voxels per block, along each axis.
    tbb::atomic< int > *        mPendingPredecessorCount    ;   /// Number of predecessors each block depends on
    mutable tbb::atomic< bool > mutexLock                   ;   /// Mutex lock, used to synchronize debug output.

函数调用运算符描述了工作流程

  void operator() ( const Block & block , tbb::parallel_do_feeder< Block > & feeder ) const
  {
      SetFloatingPointControlWord( mMasterThreadFloatingPointControlWord ) ;
      SetMmxControlStatusRegister( mMasterThreadMmxControlStatusRegister ) ;
      ComputeRemainingSDFFromImmediateSDF_Nearest_Block( mSignedDistanceGrid , mSdfPinnedGrid , mVoxelsPerBlock , mInc , block ) ;
      int nextIndices[ 3 ] = { block.mBlockIndices[ 0 ] + mInc[ 0 ]
                             , block.mBlockIndices[ 1 ] + mInc[ 1 ]
                             , block.mBlockIndices[ 2 ] + mInc[ 2 ] } ;
      FeedBlockToLoop( feeder, nextIndices, 0,         nextIndices[ 0 ] , block.mBlockIndices[ 1 ] , block.mBlockIndices[ 2 ] ) ;
      FeedBlockToLoop( feeder, nextIndices, 1, block.mBlockIndices[ 0 ] ,         nextIndices[ 1 ] , block.mBlockIndices[ 2 ] ) ;
      FeedBlockToLoop( feeder, nextIndices, 2, block.mBlockIndices[ 0 ] , block.mBlockIndices[ 1 ] ,         nextIndices[ 2 ] ) ;
  }

辅助例程 FeedBlockToLoop 将新工作供给线程调度器

  void FeedBlockToLoop( tbb::parallel_do_feeder< Block > & feeder , const int nextIndices[ 3 ] , const int axis
                       , int ix , int iy , int iz ) const
  {
      if( nextIndices[ axis ] != EndBlock( axis ) )
      {   // Next block along given direction is still within domain.
          // Atomically decrement pending-predecessor count for that block.
          if( -- Predecessor( ix , iy , iz ) == 0 )
          {   // Pending-predecessor count for adjacent block reached zero, so it is ready to process.
              // Add that block to feeder so parallel_do will process it.
              feeder.add( Block( ix , iy , iz ) ) ;
          }
      }
  }

此例程执行魔幻部分：它跟踪每个块正在等待多少个待处理的前驱块。当该计数达到零时，FeedBlockToLoop 将块供给循环。

请注意，使用原子整数来跟踪待处理前驱块的数量。它必须是原子的，因为多个线程可能同时尝试递减或读取该值。Intel TBB 提供原子整数以使此操作变得容易。Predecessor 例程只是 atomic<int> 的一个方便的访问器。

    tbb::atomic< int > & Predecessor( int ix , int iy , int iz ) const
    {
        const size_t offset = ix + NumBlocks( 0 ) * ( iy + NumBlocks( 1 ) * iz ) ;
        ASSERT( mPendingPredecessorCount[ offset ] >= 0 ) ;
        return mPendingPredecessorCount[ offset ] ;
    }

在启动工作之前，您必须初始化函数对象

    void Init( int incX , int incY , int incZ )
    {
        // Initialize predecessor counts for blocks.
        const int numBlocksX = NumBlocks( 0 ) ;
        const int numBlocksY = NumBlocks( 1 ) ;
        const int numBlocksZ = NumBlocks( 2 ) ;
        const int numBlocks = numBlocksX * numBlocksY * numBlocksZ ;
        mInc[ 0 ] = incX ;
        mInc[ 1 ] = incY ;
        mInc[ 2 ] = incZ ;
        mPendingPredecessorCount = new tbb::atomic< int >[ numBlocks ] ;
        int beginBlock[ 3 ];
        BeginBlock( beginBlock ) ;
        for( int iz = 0 ; iz < numBlocksZ ; ++ iz )
        for( int iy = 0 ; iy < numBlocksY ; ++ iy )
        for( int ix = 0 ; ix < numBlocksX ; ++ ix )
        {   // For each block...
            const size_t offset = ix + numBlocksX * ( iy + numBlocksY * iz ) ;
            mPendingPredecessorCount[ offset ] = ( ix != beginBlock[0] ) + ( iy != beginBlock[1] ) + ( iz != beginBlock[2] ) ;
        }
    }

Init 分配并初始化待处理前驱计数数组。

请注意，开始块（第一个要处理的块）取决于传播方向……

    void BeginBlock( int indices[ 3 ] ) const
    {
        for( int axis = 0 ; axis < 3 ; ++ axis )
        {
            if( mInc[ axis ] > 0 )
            {
                indices[ axis ] = 0 ;
            }
            else
            {
                indices[ axis ] = NumBlocks( axis ) - 1 ;
            }
        }
    }

……原因是上面提到的：为了求解 Eikonal 方程，求解器必须进行多次遍历，一次可以处理解可以传播的每个八分量的特征方向。

Driver 例程启动作业

static void Driver( UniformGrid< float > & signedDistanceGrid , const UniformGrid< int > & sdfPinnedGrid )
{
    static const int voxelsPerBlock[ 3 ] = { 8 , 8 , 8 } ;   // Tune for optimal performance

    ComputeRemainingSDFFromImmediateSDF_Nearest_TBB functor( signedDistanceGrid , sdfPinnedGrid , voxelsPerBlock ) ;

    static const int numPasses = 8 ;
    for( int pass = 0 ; pass < numPasses ; ++ pass )
    {   // For each pass...
        // Assign begin, end and increment for each axis.
        // Start at extremal corner.
        // Note that although this algorithm looks backwards from the direction of "travel",
        // it has to start at the first gridpoint even though it has no backward neighbors
        // *along that direction* because the other directions usually do have backward neighbors.
        // The backward-indexer below thus checks each index before accessing its backward neighbor.
        if( 0 == ( pass & 1 ))
        {
            functor.Init( -1 , -1 , -1 ) ;
        } // … (See accompanying code for details)
        else
        {
            functor.Init( 1 , 1 , 1 ) ;
        }

        Block rootBlock ;

        // Compute indices of root block based on inc.
        functor.BeginBlock( rootBlock.mBlockIndices ) ;

        // Invoke parallel_do on rootBlock, using functor as the loop body.
        // The end argument passed in here is somewhat bogus;
        // parallel_do expects an InputIterator, but there is only one root block.
        // This call therefore passes in the address of the root block (for begin), and one element past it (for end).
        // parallel_do will treat this as an STL InputIterator, because it obeys all properties of one:
        // copy-constructible, copy-assignable, destructible, comparable using == and !=, dereferenceable, incrementable.
        // The only thing that matters here for the "end" value is that, once parallel_do "increments" the "iterator"
        // and compares it to "end", the loop will terminate.
        tbb::parallel_do( & rootBlock , & rootBlock + 1 , functor ) ;
    }
}

请注意，在这种情况下，parallel_do 例程获得了一个要处理的块。那不是很多工作量。除了第一个块之外的所有工作都通过工作者本身供给 Intel TBB。

结果是工作沿着波前对角线传播到整个域，如图 7 所示。因此，这种并行化模式称为波前模式。

图 7. 波前模式

摘要

本文展示了如何使用指示函数初始化符号距离值来识别流体表面，如何通过求解 Eikonal 方程填充域的其余部分为符号距离值，以及如何使用 Intel TBB 并行化 Eikonal 求解器，使用波前并行编程模式。

未来文章

随附代码中使用的粒子渲染对于液体来说效果不佳。最好只渲染液体表面。本文介绍了如何识别表面。稍后的文章将介绍如何从中提取多边形网格，以便您可以渲染表面。您可以使用表面信息来模拟表面张力，这也是未来文章的主题。

粒子方法，包括 SPH 和涡方法，依赖于空间划分来加速邻域搜索。本系列中使用的均匀网格过于简单，效率不高，并且填充它比应有的花费更多时间。研究各种空间划分算法以确定哪种运行速度最快，尤其是借助多线程，将是有价值的。未来的文章将探讨这些问题。

延伸阅读

赵 (2004), A fast sweeping method for Eikonal equations, Mathematics of Computation, 74(250), pages 603–627.
- 描述了本文遵循的 Eikonal 方程的解。
Mark Snir, Parallel Programming Patterns.
- 描述了波前并行编程模式。
Intel® Threading Building Blocks Documentation (http://software.intel.com/sites/products/documentation/doclib/tbb_sa/help/index.htm#tbb_userguide/Design_Patterns/Wavefront.htm)
- 描述了如何使用 Intel TBB 实现波前模式。

本文档中重印的任何软件源代码均根据软件许可证提供，并且只能根据该许可证的条款使用或复制。