[测试记录用,非技术分享文]UE4-Plugin-Shallow Water

Posted 2021-05-29 wx-泡

[测试记录用,非技术分享文]

Asher说搞个shallow water作为测试题，这就来研究一下。未知结果如何（可能连效果都实现不出来.. but,always do it..）

 

看了下项目架构，好像内容并不算太多，无可演示的效果，so..造轮子？

没上cpp，有个别蓝图..hlsl不少，盲猜是搞算法为主..

 

先看SPH文件夹吧（量体较少，maybe比较简单容易入手）

（BasicSPHSystem内只有SPHEmitter一个发射器）

找个SPH的算法 https://blog.csdn.net/liuyunduo/article/details/84098884

分别是重力，压力梯度（高往低，哈密顿算子∇），粘度*粒子之间速度差

展开一下，得到对应每个粒子（ri点）的加速度（影响流体运动的计算方法）：

然后考虑光滑核函数W + 考虑各点受力不均匀（用算术平均值） 的情况，展开得：

步骤应该是：

1.计算密度

2.计算压力

3.计算加速度

4.计算表面张力影响的加速度变化

5.计算粒子位置

6.渲染

 

-------------------------------------------------------------------------备注块start

备注部分知识点：

1.在流体模拟中，常用假设：流体是不可压缩的。表明流速场是无散的。也就是流体并不收缩或者扩散，在每一点处流入速率率等于流出速率。

但注意SPH算法方法违反了此项假设。因此SPH模拟的流体会表现出现实中很少见到的流体的“弹性”。（盲猜对比Asher那只史莱姆）

缓解办法：通过增大下方第4点中的k还有γ，使得流体坚硬。但需要更短的模拟时间步长（估计在这步上需考虑性能瓶颈了）

对于这个压缩性问题，并需要模拟复杂细节飞溅效果时，已有“预估纠正方法”来解决，拓展部分，书本没描述。但提供了参考文献名字。

2.纳维-斯托克斯方程。目前研究加速度场，故用单位为 每单位体积的质量的密度场来代替质量，用单位为 每单位体积上的力的力密度场来代替力。

3.密度、黏度、运动黏度（黏度系数/密度）数值参照表

4. 关于压强梯度部分的书本描述《基于物理的建模与动画》

 

图中14.6公式比计算所用的公式更为准确，记录，以备一用。

5.对于速度差（液体扩散），根据动量 ρu （密度&速度向量）的拉普拉斯项（二阶偏导）拆分化简可得。

6.对于表面张力，通过创造有限支撑区域（也就是ri离点r的距离）的颜色场，有粒子的位置值为1、之外的位置为0，可以得到光滑变化的颜色场c（x）。

因为场临近表面的地方很大，在其他地方几乎为零（也就是表面内外与表面边界的值区别很大）。

所以通过对矢量S进行阈值处理就可以定位流体表面。（轮廓检测）

且

因为表面场从无流体区域指向有流体区域，所以对负梯度归一化可以得到表面法线场

7.如果想调整表面形状，通过对上述表面曲面矢量场进行偏导，得出曲面矢量场的散度（这是度量曲面曲率的方式）。

表面张力正比于局部曲率，方向指向流体内部。即是

8.有限差分算法暂时不研究，先看看niagara怎么做SPH流体算法的。

-------------------------------------------------------------------------备注块end

 

回来NiagaraSystem，sprite的RenderBinding并无更改，那所有模块运算的都是中间变量了：

EmitterSpawn：

Set(Emitter)NeighborGird。这个节点是用于设置邻近网格生成所依赖的数据。在GenerateNeighborGird的时候用。

其中，

MaxNeighborsPerCell就是一个Cell单元里面所能容纳最大的粒子数。

CellSize就是每一个Cell单元的大小（长宽高）

WorldBBoxSize就是邻近网格总大小，500*500*500.

SetResolutionMethod就是分辨率方法，也就是现在总共有（500/10，500/10,500/10）个分辨率大小。

 

AttributeReader没啥好说，Niagara读取上一帧粒子数据的方法。备注：EmitterName一定要与发射器名字相同。

 

Set：一堆数值的设置，目前并不知道所有的数值的含义。先把知道的写了，回头知道再补：

RestDensity：静态流体粒子密度

h：光滑核半径（粒子的属性“扩散”到周围，用来算密度、压力、加速度的）

h_2和h_3：和半径的平方和次方

 

GradW_Const：压力场梯度的光滑核常数，只跟温度有关

 

LapW_Const：黏度（速度差）的光滑核常数，只跟温度有关（也就是K）

 

W_Const：计算密度的光滑核常数。

 

Emitter Update & Particle Spawn：

没啥特别的，就是无限发射粒子，在一个球形范围随机位置往某方向锥形发射。

 

Particle Update：

三个Simulation Stage：Generate Neighbor Grid 、Compute Density、Integrate。

每个stage理解成有顺序执行的就行了，先执行前一个stage再执行后一个stage。

只有在我们的逻辑对当前帧其他粒子的属性或某些中间阶段性的结果有强的顺序依赖时，才需要SimulationStage。--https://zhuanlan.zhihu.com/p/339468114

Generate Neighbor Gird：生成邻近网格（神奇的直译）

主要目的就是将当前粒子逐个放进对应的邻近网格中，排到一个一维哈希表里面，保存对应粒子的ExecutionIndex。

P.S：如果不构造这个邻近网格，下面计算密度和SPH的时候就要遍历所有粒子，粒子数量一旦多起来就会爆炸。

下面这块代码片段是对应上面这个图的注释，对应hlsl函数在引擎源码下Source/Plugins/FX/Niagara/Source/Niagara/Private/NiagaraDataInterfaceNeighborGrid3D.cpp

OutPosition = Position;         //传进来什么位置就是什么位置


#if GPU_SIMULATION          //只在GPU模拟下进行
		// query input floats
float3 UnitPos;
//坐标变换.Simulation到单位坐标，对应SimulationToUnit的矩阵已经在前面计算
//大概是这么样的一个矩阵：原向量右乘一下这个矩阵就得出对应目标向量了
//
//  scale.x    0        0       0
//     0    scale.y     0       0
//     0       0     scale.z    0
// offset.x offset.y offset.z   0
//
NeighborGrid.SimulationToUnit(Position, SimulationToUnit, UnitPos);         

int3 Index;
//得出目标向量（坐标）后，可以得出粒子所处的Index位置
//Out：Index.x,Index.y,Index.z
NeighborGrid.UnitToIndex(UnitPos, Index.x,Index.y,Index.z);

float3 TotalDelta = {0,0,0};
uint ConstraintCount = 0;

int3 NumCells;
//取出邻近网格体的三维单元数量。比如这个例子是（500/10,500/10,500/10），那就是一个维度50个
//Out：NumCells.x,NumCells.y,NumCells.z
/*
    for (int32 InstanceIdx = 0; InstanceIdx < Context.NumInstances; ++InstanceIdx)
	{
		NumCellsX.SetAndAdvance(NumCells.X);
		NumCellsY.SetAndAdvance(NumCells.Y);
		NumCellsZ.SetAndAdvance(NumCells.Z);
	}
*/
NeighborGrid.GetNumCells(NumCells.x, NumCells.y, NumCells.z);

int MaxNeighborsPerCell;
//得到每个单元拥有的最大粒子数，这个值在一开始Niagara SetEmitter设置NeighborGird的时候就设置了
//Out：MaxNeighborsPerCell
/*
	for (int32 InstanceIdx = 0; InstanceIdx < Context.NumInstances; ++InstanceIdx)
	{
		OutMaxNeighborsPerCell.SetAndAdvance(InstData->MaxNeighborsPerCell);
	}
*/
NeighborGrid.MaxNeighborsPerCell(MaxNeighborsPerCell);
		
//如果粒子所处的Index位置在邻近网格单元最大范围内，则执行逻辑
//比如这个粒子在x=14,y=30,z=22这个index中,那是在(50,50,50)中的,所以执行
//备注：一个Index里就是有最大MaxNeighborsPerCell个粒子数
if (Index.x >= 0 && Index.x < NumCells.x && 
    Index.y >= 0 && Index.y < NumCells.y && 
	Index.z >= 0 && Index.z < NumCells.z)
{
    int LinearIndex;
    //将三维的Index转换成一维的，就是转换成哈希表的意思
    //Out：LinearIndex
    NeighborGrid.IndexToLinear(Index.x, Index.y, Index.z, LinearIndex);

    int PreviousNeighborCount;
    //读取并且设置这个Index里的粒子数+1，返回+1之前的粒子数量（记住是这个Index里面的粒子数，并不是整个Neighbor的粒子总数量）
    //Out：PreviousNeighborCount
    /*
        void {FunctionName}(int In_Index, int In_Increment, out int PreviousNeighborCount)
        {				
            InterlockedAdd(RW{OutputParticleNeighborCount}[In_Index], In_Increment, PreviousNeighborCount);				
        }
    */
    NeighborGrid.SetParticleNeighborCount(LinearIndex, 1, PreviousNeighborCount);

    //如果+1后这个Index容纳的粒子数并无超过MaxNeighborsPerCell,那就执行逻辑
    if (PreviousNeighborCount < MaxNeighborsPerCell)
    {
        int NeighborGridLinear;
        //根据三维Index数量,求出当前粒子在一维哈希表中的下标
        /*
            void {FunctionName}(int In_IndexX, int In_IndexY, int In_IndexZ, int In_Neighbor, out int Out_Linear)
			{
				Out_Linear = In_Neighbor + In_IndexX * {MaxNeighborsPerCellName} + In_IndexY * {MaxNeighborsPerCellName}*{NumCellsName}.x + In_IndexZ * {MaxNeighborsPerCellName}*{NumCellsName}.x*{NumCellsName}.y;
			}
        */
        NeighborGrid.NeighborGridIndexToLinear(Index.x, Index.y, Index.z, PreviousNeighborCount, NeighborGridLinear);

        int IGNORE;
        //根据刚求出的下标，在一维哈希表中找到对应下标的元素，并将它的值改成粒子的ExecutionIndex
        /*
            void {FunctionName}(int In_Index, int In_ParticleNeighborIndex, out int Out_Ignore)
			{
				RW{OutputParticleNeighbors}[In_Index] = In_ParticleNeighborIndex;				
				Out_Ignore = 0;
			}
        */
        NeighborGrid.SetParticleNeighbor(NeighborGridLinear, InstanceIdx, IGNORE);
    }		
}
#endif

 

然后，看计算密度的NiagaraModuleScrip（ComputeDensity），也是SimulationStage的第二项：

这项主要是计算当前粒子的密度，需要满足两个条件：

1.遍历一个正方形范围内的所有粒子，以免遗漏一些距离较远但在邻近网格内粒子（那三层for）

2.邻近粒子需要在当前计算的粒子  所在的邻近网格中，如果不在则不计算（避免粒子数量太多计算爆炸）

Density = 0;

#if GPU_SIMULATION

bool Valid;

float3 UnitPos;
//坐标变换.Simulation到单位坐标，对应SimulationToUnit的矩阵已经在前面计算
//大概是这么样的一个矩阵：原向量右乘一下这个矩阵就得出对应目标向量了
//
//  scale.x    0        0       0
//     0    scale.y     0       0
//     0       0     scale.z    0
// offset.x offset.y offset.z   0
//
NeighborGrid.SimulationToUnit(Position, SimulationToUnit, UnitPos);

int3 Index;
//得出目标向量（坐标）后，可以得出粒子所处的Index位置
//Out：Index.x,Index.y,Index.z
NeighborGrid.UnitToIndex(UnitPos, Index.x,Index.y,Index.z);

int3 NumCells;
//取出邻近网格体的三维单元数量。比如这个例子是（500/10,500/10,500/10），那就是一个维度50个
//Out：NumCells.x,NumCells.y,NumCells.z
/*
    for (int32 InstanceIdx = 0; InstanceIdx < Context.NumInstances; ++InstanceIdx)
	{
		NumCellsX.SetAndAdvance(NumCells.X);
		NumCellsY.SetAndAdvance(NumCells.Y);
		NumCellsZ.SetAndAdvance(NumCells.Z);
	}
*/
NeighborGrid.GetNumCells(NumCells.x, NumCells.y, NumCells.z);

int MaxNeighborsPerCell;
//得到每个单元拥有的最大粒子数，这个值在一开始Niagara SetEmitter设置NeighborGird的时候就设置了
//Out：MaxNeighborsPerCell
/*
	for (int32 InstanceIdx = 0; InstanceIdx < Context.NumInstances; ++InstanceIdx)
	{
		OutMaxNeighborsPerCell.SetAndAdvance(InstData->MaxNeighborsPerCell);
	}
*/
NeighborGrid.MaxNeighborsPerCell(MaxNeighborsPerCell);

//定义相互影响的粒子半径范围
int HalfWidth = 1;

// loop over all neighboring grid cells and find potential colliding particles
for (int x = -HalfWidth; x <= HalfWidth; ++x) {
for (int y = -HalfWidth; y <= HalfWidth; ++y) {
for (int z = -HalfWidth; z <= HalfWidth; ++z) {

for (int i = 0; i < MaxNeighborsPerCell; ++i)
{
    //计算出当前遍历的邻近粒子的index坐标
    const int3 IndexToUse =Index + int3(x,y,z);

	int NeighborLinearIndex;
    //将邻近粒子Index三维坐标转为一维哈希表Index位置
	NeighborGrid.NeighborGridIndexToLinear(IndexToUse.x, IndexToUse.y, IndexToUse.z, i, NeighborLinearIndex);

	int CurrNeighborIdx;
    //应该是读取上一个SimulationStage储存的ExecutionIndex
    /*
        void {FunctionName}(int In_Index, out int Out_ParticleNeighborIndex)
        {
            Out_ParticleNeighborIndex = {ParticleNeighbors}[In_Index];				
        }
    */
	NeighborGrid.GetParticleNeighbor(NeighborLinearIndex, CurrNeighborIdx);

	float3 OtherPos;
    //读取这个邻近粒子上一帧的坐标
    AttributeReader.GetVectorByIndex<Attribute="Position">(CurrNeighborIdx, Valid, OtherPos);

    //计算距离dist（也就是公式的r）
	const float3 delta = Position - OtherPos;
	const float dist = length(delta);

    //邻近粒子在邻近网格内且不等于自身，且在核函数半径范围内
	if (IndexToUse.x >= 0 && IndexToUse.x < NumCells.x && 
		IndexToUse.y >= 0 && IndexToUse.y < NumCells.y && 
		IndexToUse.z >= 0 && IndexToUse.z < NumCells.z && 
		CurrNeighborIdx != InstanceIdx &&
       CurrNeighborIdx != -1 &&
       dist > 0 && dist < h)
	{
          // increment density 
          //W_Const = 315.0/(64.0*3.1415 * pow(Emitter.h,9));光滑半径核常数
          //求和得最终密度（备注，粒子质量Mass已经在spawn particle的时候设成1了
          //所以这里乘上质量也就乘1，不写没区别）     
          float W = W_Const * pow(h_2 - dist*dist,3);
          Density += W;
	}
}       
}}}

#endif

 

然后，看ComputeSPHForce，两部分的计算，一部分算压力，一部分算黏度（速度差）

压力部分有个修正，暂时不知道是为了修正什么东西而存在的

黏度就正常照着公式算的

OutForce = float3(0,0,0);
#if GPU_SIMULATION

bool Valid;

float3 UnitPos;
//坐标变换.Simulation到单位坐标，对应SimulationToUnit的矩阵已经在前面计算
//大概是这么样的一个矩阵：原向量右乘一下这个矩阵就得出对应目标向量了
//
//  scale.x    0        0       0
//     0    scale.y     0       0
//     0       0     scale.z    0
// offset.x offset.y offset.z   0
//
NeighborGrid.SimulationToUnit(Position, SimulationToUnit, UnitPos);

int3 Index;
//得出目标向量（坐标）后，可以得出粒子所处的Index位置
//Out：Index.x,Index.y,Index.z
NeighborGrid.UnitToIndex(UnitPos, Index.x,Index.y,Index.z);

float3 TotalDelta = {0,0,0};

int3 NumCells;
//取出邻近网格体的三维单元数量。比如这个例子是（500/10,500/10,500/10），那就是一个维度50个
//Out：NumCells.x,NumCells.y,NumCells.z
/*
    for (int32 InstanceIdx = 0; InstanceIdx < Context.NumInstances; ++InstanceIdx)
	{
		NumCellsX.SetAndAdvance(NumCells.X);
		NumCellsY.SetAndAdvance(NumCells.Y);
		NumCellsZ.SetAndAdvance(NumCells.Z);
	}
*/
NeighborGrid.GetNumCells(NumCells.x, NumCells.y, NumCells.z);

int MaxNeighborsPerCell;
//得到每个单元拥有的最大粒子数，这个值在一开始Niagara SetEmitter设置NeighborGird的时候就设置了
//Out：MaxNeighborsPerCell
/*
	for (int32 InstanceIdx = 0; InstanceIdx < Context.NumInstances; ++InstanceIdx)
	{
		OutMaxNeighborsPerCell.SetAndAdvance(InstData->MaxNeighborsPerCell);
	}
*/
NeighborGrid.MaxNeighborsPerCell(MaxNeighborsPerCell);

int HalfWidth = 1;//Width / 2;

//压力
float3 PressureForce = float3(0,0,0);
//黏度
float3 ViscosityForce = float3(0,0,0);

// loop over all neighboring grid cells and find potential colliding particles
for (int x = -HalfWidth; x <= HalfWidth; ++x) {
for (int y = -HalfWidth; y <= HalfWidth; ++y) {
for (int z = -HalfWidth; z <= HalfWidth; ++z) {
for (int i = 0; i < MaxNeighborsPerCell; ++i)
{
    const int3 IndexToUse =Index + int3(x,y,z);

	int NeighborLinearIndex;
    //将邻近粒子Index三维坐标转为一维哈希表Index位置
	NeighborGrid.NeighborGridIndexToLinear(IndexToUse.x, IndexToUse.y, IndexToUse.z, i, NeighborLinearIndex);

	int CurrNeighborIdx;
    //按照一维Index位置读取上一个SimulationStage储存的ExecutionIndex
    /*
        void {FunctionName}(int In_Index, out int Out_ParticleNeighborIndex)
        {
            Out_ParticleNeighborIndex = {ParticleNeighbors}[In_Index];				
        }
    */
	NeighborGrid.GetParticleNeighbor(NeighborLinearIndex, CurrNeighborIdx);

	float3 OtherPos;
    //读取这个邻近粒子上一帧的坐标
    AttributeReader.GetVectorByIndex<Attribute="Position">(CurrNeighborIdx, Valid, OtherPos);

    //计算距离dist（也就是公式的r）
	const float3 delta = Position - OtherPos;
	const float dist = length(delta);

	if (IndexToUse.x >= 0 && IndexToUse.x < NumCells.x && 
		IndexToUse.y >= 0 && IndexToUse.y < NumCells.y && 
		IndexToUse.z >= 0 && IndexToUse.z < NumCells.z && 
		CurrNeighborIdx != InstanceIdx && 
       CurrNeighborIdx != -1 &&
       dist > 0 && dist < h)
	{
     	float OtherDensity = 0;
        //读取上一帧此邻近粒子的密度
        AttributeReader.GetFloatByIndex<Attribute="Density">(CurrNeighborIdx, Valid, OtherDensity);	

        // mass is always 1 for now
     	//float OtherMass = 0;
        //AttributeReader.GetFloatByIndex<Attribute="Mass">(CurrNeighborIdx, Valid, OtherMass);	

     	float3 OtherVelocity = float3(0,0,0);
        //获取上一帧此邻近粒子的速度
        AttributeReader.GetVectorByIndex<Attribute="Velocity">(CurrNeighborIdx, Valid, OtherVelocity);	

        //GradW_Const压力场梯度的光滑核常数: -45.0 / (3.1415 * pow(Emitter.h,6))
        float3 GradW = GradW_Const * pow(h - dist,2);

        //理想气体状态方程，根据粒子密度和静态流体密度RestDensity算粒子压力的
       float Pressure = k * (Density - RestDensity);
       float PressureOther = k * (OtherDensity - RestDensity);

        if (ClampPressure)
        {
            Pressure = max(Pressure, 0);
            PressureOther = max(PressureOther, 0);
        }
    
        //W_Const计算密度的光滑核常数：315.0/(64.0*3.1415 * pow(Emitter.h,9))
        float W_Density = (W_Const * pow(h_2 - dist*dist,3));

        //这块修正的暂时不知道用来处理什么情况的
        float h_Corr = h*PressureCorrectionSize;
        float W_Corr = (W_Const * pow(h_2 - h_Corr*h_Corr,3));

       float3 GradWSCorr = GradW_Const * pow(h - h*PressureCorrectionSize,2);
       float sCorr = -PressureCorrectionMult * pow(W_Density/W_Corr, PressureCorrectionExponent);

        PressureForce -= 1. * (Pressure + PressureOther + sCorr) / (2. * OtherDensity) * GradW * delta / dist;

        // viscosity force
        //LapW_Const：15.0/(2.*3.1415 * pow(Emitter.h,3))
        //对应方程3.17
        float LapW = LapW_Const * ((-dist*dist*dist) / (2. * h_3) + dist*dist/(h*h) + h / (2. * dist) - 1);

        //对应方程3.18
        ViscosityForce += Viscosity  / OtherDensity * (OtherVelocity - Velocity) * LapW;

    }
}       
}}}

OutForce = PressureForce + ViscosityForce;


// old integration
//OutVelocity = Velocity + TotalForce * Dt / 1.;
//float VelMag = length(OutVelocity);
//OutVelocity = OutVelocity / VelMag * min(MaxVelocity, VelMag);
//OutPosition += OutVelocity * Dt;
#endif