基于欧拉法的2D流体模拟

前言

为什么想要实现这个呢？其实以前看过些许视频，对流体的运动还挺感兴趣的，正好上个学期学了工程流体力学，所以我想来动手实现一下。

知识准备

上个学期学的工程流体力学没怎么认真学，现在回想起来只记得一点点了，就对欧拉法和拉格朗日法有点印象。欧拉法好像是追踪空间而拉格朗日法好像是追踪粒子，我想了想模拟数以万计的粒子之间的相互作用可能对我来说非常难，于是放弃了这个想法，所以只能用欧拉法了。其实欧拉法和texture这个东西非常契合，因为欧拉法研究的是流场，我们可以用texture来存储密度场和速度场以及其它物理量的场。初步想了下，最基本的肯定有密度场（颜色场）和速度场。如果想让流体比较真实地流动我们必须求解一个正确的速度场，这是我首先想到的，然后我就不知道该怎么实现了。去网上搜了搜基于欧拉法的流体模拟，搜到了对我帮助很大的两篇文章(Fast Fluid Dynamics Simulation on the GPU、Vorticity Confinement for Eulerian Fluid Simulations)还有一篇论文(Real-Time Fluid Dynamics for Games)。

流体运动的假设及术语

上面的文章和论文对要模拟的流体做出了不可压缩和均质这两个简单的假设，而这种流体随时间演化的规则可以归纳为下面两个等式。

\[\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{f} \]

\[\nabla \cdot \mathbf{u} = 0 \]

一个等式看着挺吓人的，反正我是看不懂，但是Fast Fluid Dynamics Simulation on the GPU这篇文章帮我概括了NS方程中的术语。第一流体会对流，流体的速度会导致流体传输物体、密度、速度还有其它量。这里有个很有趣的一点，速度让流体移动但是流体本身也携带着速度，因此速度让流体对流但是也会改变速度本身。第二还有压力，当外部的力作用于流体时它的影响并不会马上扩散到整个流体区域，离力近的分子会扩散力的影响。第三流体还会扩散，viscousity也就是粘度描述了流体对扩散的抗拒程度，但是这个现象一般不会模拟。第四流体会受到外部力的影响，而且主要体现在改变速度场，这是显而易见的一点。第二个等式很好理解，流体速度场中任意一点求得的散度都是0，这正是因为我们之前对要模拟的流体做了不可压缩和均质这两个简单的假设。

程序模拟流体运动的步骤

接下来我们了解下该如何用程序模拟这种流体的运动，Real-Time Fluid Dynamics for Games这篇论文概括了模拟这种流体运动的步骤：施加外部影响、扩散、对流，如下图所示。施加外部影响和对流这两个操作都会给我们一个新的速度场，但是新的速度场会不可避免地有散度，而我们需要的是一个不可压缩的速度场让流体对流，所以我们得用某种方法消除新的速度场的散度，这个过程就叫做project，下面的部分将说明该如何通过新的速度场来求解一个不可压缩速度场。

不可压缩速度场的计算

由亥姆霍兹分解定理和NS方程可知我们的速度场其实是由一个不可压缩的速度场（无散度）和一个压力的梯度场组成。

令新的速度场为w，它的不可压缩场为u，它的压力场为p那么可知

\[\mathbf{u} = \mathbf{w} - \nabla p \]

用这个公式我们可以求出不可压缩速度场，但前提是我们得先求出压力场，所以接下来我们求压力场，对新的速度场w求散度可知

\[\nabla \cdot \mathbf{w} = \nabla \cdot \mathbf{u} + \nabla^2 p \]

由于

\[\nabla \cdot \mathbf{u} = 0 \]

化简可得

\[\nabla \cdot \mathbf{w} = \nabla^2 p \]

Fast Fluid Dynamics Simulation on the GPU给出了相关的离散求解公式，我们令速度矢量为(u,v)

\[\nabla \cdot \mathbf{w}_{i,j} = \frac{u_{i+1,j}-u_{i-1,j}}{2 \Delta x} + \frac{v_{i,j+1}-v_{i,j-1}}{2 \Delta y} \]

\[\nabla^2 p_{i,j} = \frac{p_{i+1,j} - 2p_{i,j} + p_{i-1,j}}{(\Delta x)^2} + \frac{p_{i,j+1} - 2p_{i,j} + p_{i,j-1}}{(\Delta y)^2} \]

在这次模拟中Δx和Δy都等于1，于是公式最终化简变成了

\[\nabla \cdot \mathbf{w}_{i,j} = 0.5 \cdot \left( u_{i+1,j}-u_{i-1,j} + v_{i,j+1}-v_{i,j-1} \right) \]

\[p_{i+1,j} + p_{i-1,j} + p_{i,j+1} + p_{i,j-1} - 4p_{i,j} = \nabla \cdot \mathbf{w}_{i,j} \]

在某一帧中新的速度场的散度是个定值而且可以求出来，所以这下一下子看出来了第二个等式其实是线性方程组的一部分，所以我们可以用GPU实现Jacobi迭代来求解我们的压力场，下面是求解公式

\[p_{i,j}^{\text{new}} = \frac{1}{4} \left( p_{i+1,j} + p_{i-1,j} + p_{i,j+1} + p_{i,j-1} - \nabla \cdot \mathbf{w}_{i,j} \right) \]

论文和文章中提到用迭代法计算压力场前应该先把压力场重置为0,迭代次数为20-50次即可。走到这里我们已经完成了最关键的消除散度部分，其实还可以在CPU这边用超松弛迭代法来求解我们的压力场。下面我们把流体模拟分为这几个详细的步骤：施加外部影响、求散度场、求压力场、求不可压缩速度场、流体对流、涡度限制（可选）、设置圆形障碍（可选），然后用hlsl把每个步骤实现。

实现过程

之前已经提到了会省略扩散部分，所以不实现这个。首先为了节省计算开销用于模拟的材质是颜色材质宽高的1/2，然后我准备在一个矩形区域内模拟流体运动，并且在这个矩形区域内放置一个圆形障碍物。

边界条件处理

由于是在一个放置了圆形障碍物的封闭矩形区域内模拟流体运动，因此速度、散度、压力在边界处的数值很值得注意，其实这涉及到边界条件的处理，如果你不考虑这个那么模拟肯定会失效。Real-Time Fluid Dynamics for Games这篇论文说明了该如何处理边界条件，此外它也说明了该如何模拟流体运动而且还贴出了具体的C语言代码来借鉴，于是我选择用这篇论文使用的边界处理方法，下面是这篇论文使用的边界处理函数。

void set_bnd(int N, int b, float* x) 
{ 
  int i; 
 
  for(i=1; i<=N; i++)
  { 
    x[IX(0  ,i)] = b==1 ? –x[IX(1,i)] : x[IX(1,i)]; 
    x[IX(N+1,i)] = b==1 ? –x[IX(N,i)] : x[IX(N,i)]; 
    x[IX(i,0  )] = b==2 ? –x[IX(i,1)] : x[IX(i,1)]; 
    x[IX(i,N+1)] = b==2 ? –x[IX(i,N)] : x[IX(i,N)]; 
  } 

  x[IX(0  ,0  )] = 0.5*(x[IX(1,0  )]+x[IX(0  ,1)]); 
  x[IX(0  ,N+1)] = 0.5*(x[IX(1,N+1)]+x[IX(0  ,N)]); 
  x[IX(N+1,0  )] = 0.5*(x[IX(N,0  )]+x[IX(N+1,1)]); 
  x[IX(N+1,N+1)] = 0.5*(x[IX(N,N+1)]+x[IX(N+1,N)]); 
}

对于速度边界，在竖直边界上的速度应该复制旁边Texel的速度然后翻转水平速度，而在水平边界上则做相反处理，最后角点处的速度为0。对于散度和压力的边界条件，采用的操作很简单直接复制旁边Texel的值即可。上述就是这篇论文用到的处理边界的方法，关于圆形边界的速度、压力处理部分会在后面介绍，我们会用到和set_bnd函数类似的技巧。

施加外部影响

施加的外部影响会作用到颜色场和速度场上，为了满足NS方程术语中压力所描述的规则，我们借助下方的两个公式来分别喷洒颜料和修改速度场。用鼠标的位置和随机HSV颜色来喷洒颜料，用鼠标的位置和鼠标位置的变化来修改速度场，修改速度场的时候记得注意边界条件。

\[exp \left( - \frac{ \left( x - x_{\mathbf{p}} \right) ^ 2 + \left( y - y_{\mathbf{p}} \right) ^2 }{r} \right) \cdot \mathbf{color} \]

\[exp \left( - \frac{ \left( x - x_{\mathbf{p}} \right) ^ 2 + \left( y - y_{\mathbf{p}} \right) ^2 }{r} \right) \cdot \Delta \mathbf{p} \]

float3 ColorAt(const uint2 loc)
{
    float2 texCoord = (float2(loc) + float2(0.5, 0.5)) * colorTexelSize;
        
    texCoord -= pos;
    
    const float aspectRatio = 16.0 / 9.0;
    
    texCoord.x *= aspectRatio;
    
    const float3 color = exp(-dot(texCoord, texCoord) / splatRadius) * splatColor.rgb * 0.15;
        
    const float3 curColor = colorReadTex[loc].rgb;
    
    return color + curColor;
}

[numthreads(16, 9, 1)]
void main(const uint2 DTid : SV_DispatchThreadID)
{
    const float scale = simTexelSize.x / colorTexelSize.x;
    
    const float scaledRadius = obstacleRadius * scale;
    
    const float2 scaledPosition = obstaclePosition * scale;
    
    const float2 dir = (float2(DTid) + float2(0.5, 0.5)) - scaledPosition;
    
    if (DTid.x == 0 || DTid.x == colorTextureSize.x - 1 || DTid.y == 0 || DTid.y == colorTextureSize.y - 1 || length(dir) < scaledRadius)
    {
        colorWriteTex[DTid] = float4(0.0, 0.0, 0.0, 1.0);
    }
    else
    {
        colorWriteTex[DTid] = float4(ColorAt(DTid), 1.0);
    }
}

float2 VelocityAt(const uint2 loc)
{
    float2 texCoord = (float2(loc) + float2(0.5, 0.5)) * simTexelSize;
    
    texCoord -= pos;
    
    const float aspectRatio = 16.0 / 9.0;
    
    texCoord.x *= aspectRatio;
    
    const float2 velocity = exp(-dot(texCoord, texCoord) / splatRadius) * posDelta;
    
    const float2 curVelocity = velocityReadTex[loc];
    
    return curVelocity + velocity;
}

[numthreads(16, 9, 1)]
void main(const uint2 DTid : SV_DispatchThreadID)
{
    //interior texel
    if (DTid.x > 0 && DTid.x < simTextureSize.x - 1 && DTid.y > 0 && DTid.y < simTextureSize.y - 1)
    {
        velocityWriteTex[DTid] = VelocityAt(DTid);
    }
    //row texel
    else if (DTid.x > 0 && DTid.x < simTextureSize.x - 1)
    {
        if (DTid.y == 0)
        {
            float2 velocity = VelocityAt(uint2(DTid.x, 1));
                
            velocity.y = -velocity.y;
                
            velocityWriteTex[DTid] = velocity;
        }
        else if (DTid.y == simTextureSize.y - 1)
        {
            float2 velocity = VelocityAt(uint2(DTid.x, simTextureSize.y - 2));
                
            velocity.y = -velocity.y;
                
            velocityWriteTex[DTid] = velocity;
        }
    }
    //column texel
    else if (DTid.y > 0 && DTid.y < simTextureSize.y - 1)
    {
        if (DTid.x == 0)
        {
            float2 velocity = VelocityAt(uint2(1, DTid.y));
                
            velocity.x = -velocity.x;
                
            velocityWriteTex[DTid] = velocity;
        }
        else if (DTid.x == simTextureSize.x - 1)
        {
            float2 velocity = VelocityAt(uint2(simTextureSize.x - 2, DTid.y));
                
            velocity.x = -velocity.x;
                
            velocityWriteTex[DTid] = velocity;
        }
    }
    //corner texel
    else if ((DTid.x == 0 || DTid.x == simTextureSize.x - 1) && (DTid.y == 0 || DTid.y == simTextureSize.y - 1))
    {
        velocityWriteTex[DTid] = float2(0.0, 0.0);
    }
}

求散度场

施加外部影响后我们先求新的速度场的散度场，计算直接按公式来即可，记得注意边界条件。

float DivergenceAt(const uint2 loc)
{
    const float uR = velocityTex[loc + uint2(1, 0)].x;
        
    const float uL = velocityTex[loc - uint2(1, 0)].x;
        
    const float vT = velocityTex[loc + uint2(0, 1)].y;
        
    const float vB = velocityTex[loc - uint2(0, 1)].y;
        
    const float divergence = 0.5 * (uR - uL + vT - vB);
    
    return divergence;
}

[numthreads(16, 9, 1)]
void main(const uint2 DTid : SV_DispatchThreadID)
{
    //interior texel
    if (DTid.x > 0 && DTid.x < simTextureSize.x - 1 && DTid.y > 0 && DTid.y < simTextureSize.y - 1)
    {
        divergenceTex[DTid] = DivergenceAt(DTid);
    }
    //row texel
    else if (DTid.x > 0 && DTid.x < simTextureSize.x - 1)
    {
        if (DTid.y == 0)
        {
            divergenceTex[DTid] = DivergenceAt(uint2(DTid.x, 1));
        }
        else if (DTid.y == simTextureSize.y - 1)
        {
            divergenceTex[DTid] = DivergenceAt(uint2(DTid.x, simTextureSize.y - 2));
        }
    }
    //column texel
    else if (DTid.y > 0 && DTid.y < simTextureSize.y - 1)
    {
        if (DTid.x == 0)
        {
            divergenceTex[DTid] = DivergenceAt(uint2(1, DTid.y));
        }
        else if (DTid.x == simTextureSize.x - 1)
        {
            divergenceTex[DTid] = DivergenceAt(uint2(simTextureSize.x - 2, DTid.y));
        }
    }
    //corner texel
    else if (DTid.y == 0)
    {
        if (DTid.x == 0)
        {
            divergenceTex[DTid] = DivergenceAt(uint2(1, 1));
        }
        else if (DTid.x == simTextureSize.x - 1)
        {
            divergenceTex[DTid] = DivergenceAt(uint2(simTextureSize.x - 2, 1));
        }
    }
    else if (DTid.y == simTextureSize.y - 1)
    {
        if (DTid.x == 0)
        {
            divergenceTex[DTid] = DivergenceAt(uint2(1, simTextureSize.y - 2));
        }
        else if (DTid.x == simTextureSize.x - 1)
        {
            divergenceTex[DTid] = DivergenceAt(uint2(simTextureSize.x - 2, simTextureSize.y - 2));
        }
    }
}

求压力场

现在我们已求得散度场，可以用散度场和Jacobi迭代来求解压力场，也是按公式来即可，记得在求解压力场前先将其重置为0然后还要注意边界条件。

float PressureAt(const uint2 loc)
{
    const float R = pressureReadTex[loc + uint2(1, 0)];
        
    const float L = pressureReadTex[loc - uint2(1, 0)];
        
    const float T = pressureReadTex[loc + uint2(0, 1)];
        
    const float B = pressureReadTex[loc - uint2(0, 1)];
        
    const float C = divergenceTex[loc];
    
    const float pressure = (R + L + T + B - C) * 0.25;
    
    return pressure;
}

[numthreads(16, 9, 1)]
void main(const uint2 DTid : SV_DispatchThreadID)
{
    //interior texel
    if (DTid.x > 0 && DTid.x < simTextureSize.x - 1 && DTid.y > 0 && DTid.y < simTextureSize.y - 1)
    {
        pressureWriteTex[DTid] = PressureAt(DTid);
    }
    //row texel
    else if (DTid.x > 0 && DTid.x < simTextureSize.x - 1)
    {
        if (DTid.y == 0)
        {
            pressureWriteTex[DTid] = PressureAt(uint2(DTid.x, 1));
        }
        else if (DTid.y == simTextureSize.y - 1)
        {
            pressureWriteTex[DTid] = PressureAt(uint2(DTid.x, simTextureSize.y - 2));
        }
    }
    //column texel
    else if (DTid.y > 0 && DTid.y < simTextureSize.y - 1)
    {
        if (DTid.x == 0)
        {
            pressureWriteTex[DTid] = PressureAt(uint2(1, DTid.y));
        }
        else if (DTid.x == simTextureSize.x - 1)
        {
            pressureWriteTex[DTid] = PressureAt(uint2(simTextureSize.x - 2, DTid.y));
        }
    }
    //corner texel
    else if (DTid.y == 0)
    {
        if (DTid.x == 0)
        {
            pressureWriteTex[DTid] = PressureAt(uint2(1, 1));
        }
        else if (DTid.x == simTextureSize.x - 1)
        {
            pressureWriteTex[DTid] = PressureAt(uint2(simTextureSize.x - 2, 1));
        }
    }
    else if (DTid.y == simTextureSize.y - 1)
    {
        if (DTid.x == 0)
        {
            pressureWriteTex[DTid] = PressureAt(uint2(1, simTextureSize.y - 2));
        }
        else if (DTid.x == simTextureSize.x - 1)
        {
            pressureWriteTex[DTid] = PressureAt(uint2(simTextureSize.x - 2, simTextureSize.y - 2));
        }
    }
}

求不可压缩速度场

现在已求得压力场，我们可以按公式求出一个不可压缩的速度场，按公式来即可，记得注意边界条件。

float2 VelocityAt(const uint2 loc)
{
    const float R = pressureTex[loc + uint2(1, 0)];
    
    const float L = pressureTex[loc - uint2(1, 0)];
    
    const float T = pressureTex[loc + uint2(0, 1)];
    
    const float B = pressureTex[loc - uint2(0, 1)];
    
    const float2 velocity = velocityReadTex[loc] - 0.5 * float2(R - L, T - B);
    
    return velocity;
}

[numthreads(16, 9, 1)]
void main(const uint2 DTid : SV_DispatchThreadID)
{
    //interior texel
    if (DTid.x > 0 && DTid.x < simTextureSize.x - 1 && DTid.y > 0 && DTid.y < simTextureSize.y - 1)
    {
        velocityWriteTex[DTid] = VelocityAt(DTid);
    }
    //row texel
    else if (DTid.x > 0 && DTid.x < simTextureSize.x - 1)
    {
        if (DTid.y == 0)
        {
            float2 velocity = VelocityAt(uint2(DTid.x, 1));
            
            velocity.y = -velocity.y;
            
            velocityWriteTex[DTid] = velocity;
        }
        else if (DTid.y == simTextureSize.y - 1)
        {
            float2 velocity = VelocityAt(uint2(DTid.x, simTextureSize.y - 2));
             
            velocity.y = -velocity.y;
                
            velocityWriteTex[DTid] = velocity;
        }
    }
    //column texel
    else if (DTid.y > 0 && DTid.y < simTextureSize.y - 1)
    {
        if (DTid.x == 0)
        {
            float2 velocity = VelocityAt(uint2(1, DTid.y));
            
            velocity.x = -velocity.x;
            
            velocityWriteTex[DTid] = velocity;
        }
        else if (DTid.x == simTextureSize.x - 1)
        {
            float2 velocity = VelocityAt(uint2(simTextureSize.x - 2, DTid.y));
            
            velocity.x = -velocity.x;
            
            velocityWriteTex[DTid] = velocity;
        }
    }
    //corner texel
    else if ((DTid.x == 0 || DTid.x == simTextureSize.x - 1) && (DTid.y == 0 || DTid.y == simTextureSize.y - 1))
    {
        velocityWriteTex[DTid] = float2(0.0, 0.0);
    }
}

流体对流

现在已求得一个不可压缩的速度场，我们可以用这个速度场来让流体对流。实现流体对流也很简单，通过当前的位置和当前位置的速度我们可以反向推出当前位置的下一帧的流体是哪里贡献的，用这张图可以轻易理解。

接下来我们分别实现颜色场对流以及速度场对流。

float3 ColorAt(const uint2 loc)
{
    float2 texCoord = (float2(loc) + float2(0.5, 0.5)) * colorTexelSize;
    
    texCoord -= perframeResource.deltaTime * velocityTex.SampleLevel(linearClampSampler, texCoord, 0.0) * simTexelSize;
    
    const float3 color = colorReadTex.SampleLevel(linearClampSampler, texCoord, 0.0).rgb;
    
    const float decay = 1.0 + colorDissipationSpeed * perframeResource.deltaTime;
    
    return color / decay;
}

[numthreads(16, 9, 1)]
void main(const uint2 DTid : SV_DispatchThreadID)
{
    if (DTid.x == 0 || DTid.x == colorTextureSize.x - 1 || DTid.y == 0 || DTid.y == colorTextureSize.y - 1)
    {
        colorWriteTex[DTid] = float4(0.0, 0.0, 0.0, 1.0);
    }
    else
    {
        colorWriteTex[DTid] = float4(ColorAt(DTid), 1.0);
    }
}

float2 VelocityAt(const uint2 loc)
{
    float2 texCoord = (float2(loc) + float2(0.5, 0.5)) * simTexelSize;
    
    texCoord -= perframeResource.deltaTime * velocityReadTex[loc] * simTexelSize;

    const float2 velocity = velocityReadTex.SampleLevel(linearClampSampler, texCoord, 0.0);
    
    const float decay = 1.0 + velocityDissipationSpeed * perframeResource.deltaTime;
    
    return velocity / decay;
}

[numthreads(16, 9, 1)]
void main(const uint2 DTid : SV_DispatchThreadID)
{
    //interior texel
    if (DTid.x > 0 && DTid.x < simTextureSize.x - 1 && DTid.y > 0 && DTid.y < simTextureSize.y - 1)
    {
        velocityWriteTex[DTid] = VelocityAt(DTid);
    }
    //row texel
    else if (DTid.x > 0 && DTid.x < simTextureSize.x - 1)
    {
        if (DTid.y == 0)
        {
            float2 velocity = VelocityAt(uint2(DTid.x, 1));
                
            velocity.y = -velocity.y;
                
            velocityWriteTex[DTid] = velocity;
        }
        else if (DTid.y == simTextureSize.y - 1)
        {
            float2 velocity = VelocityAt(uint2(DTid.x, simTextureSize.y - 2));
                
            velocity.y = -velocity.y;
                
            velocityWriteTex[DTid] = velocity;
        }
    }
    //column texel
    else if (DTid.y > 0 && DTid.y < simTextureSize.y - 1)
    {
        if (DTid.x == 0)
        {
            float2 velocity = VelocityAt(uint2(1, DTid.y));
                
            velocity.x = -velocity.x;
                
            velocityWriteTex[DTid] = velocity;
        }
        else if (DTid.x == simTextureSize.x - 1)
        {
            float2 velocity = VelocityAt(uint2(simTextureSize.x - 2, DTid.y));
                
            velocity.x = -velocity.x;
                
            velocityWriteTex[DTid] = velocity;
        }
    }
    //corner texel
    else if ((DTid.x == 0 || DTid.x == simTextureSize.x - 1) && (DTid.y == 0 || DTid.y == simTextureSize.y - 1))
    {
        velocityWriteTex[DTid] = float2(0.0, 0.0);
    }
}

涡度限制（可选）

到刚才那一步就完成了所有基础的步骤，但是我自己跑了下代码发现流体的运动稍微有点没意思。后面又看了下之前提到的的文章的扩展部分，于是发现了一个讲解Vorticity Confinement即涡度限制的一个部分，这个部分讲了使用涡度限制的话能让流体有那种炫酷的旋转特性，这篇Vorticity Confinement for Eulerian Fluid Simulations文章提供了涡度限制实现的具体代码。首先我们求得涡度场然后对速度场进行涡度限制。

我们令涡度为ω

\[\omega_{i,j} = 0.5 \cdot \left( u_{i,j+1} - u_{i,j-1} + v_{i-1,j} - v_{i+1,j} \right) \]

当前位置的加速度为

\[\mathbf{a_{i,j}} = \frac{0.5 \cdot \left( \left| \omega_{i,j-1} \right| - \left| \omega_{i,j+1} \right| , \left| \omega_{i+1,j} \right| - \left| \omega_{i-1,j} \right| \right)}{\left\| 0.5 \cdot \left( \left| \omega_{i,j-1} \right| - \left| \omega_{i,j+1} \right| , \left| \omega_{i+1,j} \right| - \left| \omega_{i-1,j} \right| \right) \right\|} \cdot vorticityIntensity \cdot \omega_{i,j} \]

接下来我们分别实现求涡度场和涡度限制

float VorticityAt(const uint2 loc)
{
    const float vR = velocityTex[loc + uint2(1, 0)].y;
    
    const float vL = velocityTex[loc - uint2(1, 0)].y;
    
    const float uT = velocityTex[loc + uint2(0, 1)].x;
    
    const float uB = velocityTex[loc - uint2(0, 1)].x;
    
    const float vorticity = 0.5 * (uT - uB + vL - vR);
    
    return vorticity;
}

[numthreads(16, 9, 1)]
void main(const uint2 DTid : SV_DispatchThreadID)
{
    //interior texel
    if (DTid.x > 0 && DTid.x < simTextureSize.x - 1 && DTid.y > 0 && DTid.y < simTextureSize.y - 1)
    {
        vorticityTex[DTid] = VorticityAt(DTid);
    }
    //row texel
    else if (DTid.x > 0 && DTid.x < simTextureSize.x - 1)
    {
        if (DTid.y == 0)
        {
            vorticityTex[DTid] = VorticityAt(uint2(DTid.x, 1));
        }
        else if (DTid.y == simTextureSize.y - 1)
        {
            vorticityTex[DTid] = VorticityAt(uint2(DTid.x, simTextureSize.y - 2));
        }
    }
    //column texel
    else if (DTid.y > 0 && DTid.y < simTextureSize.y - 1)
    {
        if (DTid.x == 0)
        {
            vorticityTex[DTid] = VorticityAt(uint2(1, DTid.y));
        }
        else if (DTid.x == simTextureSize.x - 1)
        {
            vorticityTex[DTid] = VorticityAt(uint2(simTextureSize.x - 2, DTid.y));
        }
    }
    //corner texel
    else if (DTid.y == 0)
    {
        if (DTid.x == 0)
        {
            vorticityTex[DTid] = VorticityAt(uint2(1, 1));
        }
        else if (DTid.x == simTextureSize.x - 1)
        {
            vorticityTex[DTid] = VorticityAt(uint2(simTextureSize.x - 2, 1));
        }
    }
    else if (DTid.y == simTextureSize.y - 1)
    {
        if (DTid.x == 0)
        {
            vorticityTex[DTid] = VorticityAt(uint2(1, simTextureSize.y - 2));
        }
        else if (DTid.x == simTextureSize.x - 1)
        {
            vorticityTex[DTid] = VorticityAt(uint2(simTextureSize.x - 2, simTextureSize.y - 2));
        }
    }
}

float2 VelocityAt(const uint2 loc)
{
    const float L = vorticityTex[loc - uint2(1, 0)];
    
    const float R = vorticityTex[loc + uint2(1, 0)];
    
    const float T = vorticityTex[loc + uint2(0, 1)];
    
    const float B = vorticityTex[loc - uint2(0, 1)];
    
    const float C = vorticityTex[loc];
    
    float2 acc = 0.5 * float2(abs(B) - abs(T), abs(R) - abs(L));
    
    acc = acc / (length(acc) + 1e-5) * vorticityIntensity * C;
    
    const float2 curVelocity = velocityReadTex[loc];
    
    return curVelocity + acc * perframeResource.deltaTime;
}

[numthreads(16, 9, 1)]
void main(const uint2 DTid : SV_DispatchThreadID)
{
    //interior texel
    if (DTid.x > 0 && DTid.x < simTextureSize.x - 1 && DTid.y > 0 && DTid.y < simTextureSize.y - 1)
    {
        velocityWriteTex[DTid] = VelocityAt(DTid);
    }
    //row texel
    else if (DTid.x > 0 && DTid.x < simTextureSize.x - 1)
    {
        if (DTid.y == 0)
        {
            float2 velocity = VelocityAt(uint2(DTid.x, 1));
                
            velocity.y = -velocity.y;
                
            velocityWriteTex[DTid] = velocity;
        }
        else if (DTid.y == simTextureSize.y - 1)
        {
            float2 velocity = VelocityAt(uint2(DTid.x, simTextureSize.y - 2));
                
            velocity.y = -velocity.y;
                
            velocityWriteTex[DTid] = velocity;
        }
    }
    //column texel
    else if (DTid.y > 0 && DTid.y < simTextureSize.y - 1)
    {
        if (DTid.x == 0)
        {
            float2 velocity = VelocityAt(uint2(1, DTid.y));
                
            velocity.x = -velocity.x;
                
            velocityWriteTex[DTid] = velocity;
        }
        else if (DTid.x == simTextureSize.x - 1)
        {
            float2 velocity = VelocityAt(uint2(simTextureSize.x - 2, DTid.y));
                
            velocity.x = -velocity.x;
                
            velocityWriteTex[DTid] = velocity;
        }
    }
    //corner texel
    else if ((DTid.x == 0 || DTid.x == simTextureSize.x - 1) && (DTid.y == 0 || DTid.y == simTextureSize.y - 1))
    {
        velocityWriteTex[DTid] = float2(0.0, 0.0);
    }
}

设置圆形障碍（可选）

关于在圆形边界处速度和压力的处理其实很简单，我们可以借助线性采样器来帮我们采样周围的Texel然后复制翻转速度以及复制压力值。我选择采样以圆心到圆形边界处为方向且距离为2个Texel的位置，采样距离为1个Texel应该也没关系，反正只要保证连续性应该就行，下面是具体实现的代码。

[numthreads(16, 9, 1)]
void main(const uint2 DTid : SV_DispatchThreadID)
{
    const float2 position = float2(DTid) + float2(0.5, 0.5);
    
    const float2 dir = position - obstaclePosition;
    
    const float len = length(dir);
    
    if (len < obstacleRadius)
    {
        velocityWriteTex[DTid] = float2(0.0, 0.0);
    }
    else if (len < obstacleRadius + 1.0)
    {
        velocityWriteTex[DTid] = -velocityReadTex.SampleLevel(linearClampSampler, (position + normalize(dir) * 2.0) * simTexelSize, 0.0);
    }
    else
    {
        velocityWriteTex[DTid] = velocityReadTex[DTid];
    }
}

[numthreads(16, 9, 1)]
void main(const uint2 DTid : SV_DispatchThreadID)
{
    const float2 position = float2(DTid) + float2(0.5, 0.5);
    
    const float2 dir = position - obstaclePosition;
    
    const float len = length(dir);
    
    if (len < obstacleRadius)
    {
        pressureWriteTex[DTid] = 0.0;
    }
    else if (len < obstacleRadius + 1.0)
    {
        pressureWriteTex[DTid] = pressureReadTex.SampleLevel(linearClampSampler, (position + normalize(dir) * 2.0) * simTexelSize, 0.0);
    }
    else
    {
        pressureWriteTex[DTid] = pressureReadTex[DTid];
    }
}

后处理

上述就是流体模拟必要以及可选的操作，接下来的渲染环节我们可以做点后处理特效。关于对流体颜色材质进行后处理我只想到了泛光特效于是问了下GPT，AI给了如下意见

我觉得第6个Normal Mapping for Lighting和第9个Edge Detection & Highlighting非常好而且对我来说比较容易实现。下面我们先实现Normal Mapping for Lighting接着施加泛光特效，最后实现Edge Detection & Highlighting。为了实现Normal Mapping for Lighting我们得先把颜色转化为高度，有了高度我们可以计算法线，最后进行冯氏光照着色。

首先获取当前Texel以及周边Texel的颜色信息

const float3 center = colorTex[DTid].rgb;
    
const float3 right = colorTex[DTid + uint2(1, 0)].rgb;

const float3 left = colorTex[DTid - uint2(1, 0)].rgb;

const float3 top = colorTex[DTid + uint2(0, 1)].rgb;

const float3 bottom = colorTex[DTid - uint2(0, 1)].rgb;

如果要将颜色转化为高度我只想到了一个方法：那就是点乘(0.299, 0.587, 0.114)转化为亮度

const float leftHeight = dot(left, float3(0.299, 0.587, 0.114));

const float rightHeight = dot(right, float3(0.299, 0.587, 0.114));

const float topHeight = dot(top, float3(0.299, 0.587, 0.114));

const float bottomHeight = dot(bottom, float3(0.299, 0.587, 0.114));

接下来用下面的公式计算法线

\[\mathbf{n} = normalize \left( - \frac{\Delta h}{\Delta x} , - \frac{\Delta h}{\Delta y} , bumpScale \right) \]

转化为代码

const float dx = (rightHeight - leftHeight) * 0.5;

const float dy = (topHeight - bottomHeight) * 0.5;

const float3 normal = normalize(float3(-dx, -dy, bumpScale));

这里要记得把bumpScale调到合适的值，我用的数值是1/300，要不然会出现下图这个情况，出现了大量的持续时间比较长的圆形物体，这会让漫反射着色变得非常不好看

接下来引入光线方向然后进行冯氏光照着色

const float3 light = float3(0.0, 0.0, 1.0);
    
const float3 ambient = kA * center;

const float3 diffuse = kD * center * saturate(dot(normal, light));

const float3 color = ambient + diffuse;

originTex[DTid] = float4(color, 1.0);

记得把kA和kD这两个系数调整到合适的值，我用的是0.6和0.4，其实还想加入镜面高光来着，但是发现没带来质的改变，完成Normal Mapping for Lighting后我们加入泛光特效，代码我很久之前就写好了

TextureRenderView* outputTexture = nullptr;

if (config.phongShading)
{
	context->setPipelineState(phongShadeState.Get());
	context->setCSConstants({ colorTex->read()->getAllSRVIndex(),phongShadeTexture->getUAVMipIndex(0) }, 0);
	context->transitionResources();
	context->dispatch(phongShadeTexture->getTexture()->getWidth() / 16, phongShadeTexture->getTexture()->getHeight() / 9, 1);
	context->uavBarrier({ phongShadeTexture->getTexture() });

	outputTexture = effect->process(phongShadeTexture);
}
else
{
	outputTexture = effect->process(colorTex->read());
}

对比了下结果，使用了冯氏光照更清晰一些，细细品味了下有点像是在喷洒油彩一样，我感觉其实各有所好吧，下面是对比的截图

使用：

不使用：

加入泛光特效后我们实现Edge Detection & Highlighting，采用的步骤其实和Normal Mapping for Lighting非常相似，我们先把颜色转化为亮度，接着求亮度改变的幅值，最后用幅值强调边缘处的颜色。

const float3 center = colorTex[DTid].rgb;
    
const float3 right = colorTex[DTid + uint2(1, 0)].rgb;

const float3 left = colorTex[DTid - uint2(1, 0)].rgb;

const float3 top = colorTex[DTid + uint2(0, 1)].rgb;

const float3 bottom = colorTex[DTid - uint2(0, 1)].rgb;

const float leftLuma = dot(left, float3(0.299, 0.587, 0.114));

const float rightLuma = dot(right, float3(0.299, 0.587, 0.114));

const float topLuma = dot(top, float3(0.299, 0.587, 0.114));

const float bottomLuma = dot(bottom, float3(0.299, 0.587, 0.114));
   
const float dx = (rightLuma - leftLuma) * 0.5;

const float dy = (topLuma - bottomLuma) * 0.5;

const float edgeMagnitude = length(float2(dx, dy)) * edgeMagnitudeScale;

const float3 color = center * (1.0 + edgeMagnitude);

edgeHighlightTex[DTid] = float4(color, 1.0);

对比了下结果，使用了边缘检测及高光后边缘处的颜色确实更亮了一些，调高edgeMagnitudeScale会让流体有种果冻的感觉，我觉得这个也是看人喜好来说的，下面是对比的截图

使用：