初次接触流体模拟还是在游戏程序设计这门课上，这门课的大作业是一个小组合作项目，要求基于Unity或UE开发一款游戏，最好实现一个算法。一筹莫展之际，朋友提到他看到的一位前辈关于流体模拟的一系列文章SPH（光滑粒子流体动力学）流体模拟实现：算法总览_sph算法计算流体-CSDN博客……2019到2020年跨年的那一周，在我们在学院，三个人实现了最简单的流体模拟部分，做了一个非常简单的小游戏。这个项目以后有机会也会添加进来。在这个小游戏完成之际，我就产生了一些困惑，像是如何保证流体不可压缩、如何处理与其他物体的碰撞……很可惜，这些困惑搁置了许久，我也在不知不觉中渐渐远离了图形学和游戏开发。再到后来（此处省略许多经历QAQ），我跌跌撞撞地又回到了喜欢的方向上。在前辈的建议下，学习了Direct3D 12。在寻找几个练习项目时，又发现了一些非常优秀的前辈的项目，又看到了熟悉的SPH。一面是为了给自己找工作的时候添加一些项目经验，一面是弥补前些年留下的遗憾，我在D3D上用Computer Shader又实现了一次流体模拟。不同的是，这次实现的是基于PBD（Postion Based Dynamic）优化的SPH算法——PBF，也将部分计算从CPU搬到了GPU上。这篇博客将从PBD、PBF、D3D上PBF实现，这三个部分进行介绍。

PBD和PBF的许多公式网上已经介绍的非常多了，而且基本相似，这里仅仅做一些自己的理解和梳理，可能不是非常全面，最近时间太紧，具体的推导公式后续会补充上来。想要了解理论细节可以参看下面的链接。

本文参考了许多优秀的文章：

Unity物理引擎实战-基于PBD方法的水体模拟（一） - 知乎 (zhihu.com)

流体模拟Fluid Simulation：Position Based Fluid | YangWC's Blog

【深入浅出 Nvidia FleX】(1) Position Based Dynamics - 知乎 (zhihu.com)

【深入浅出 Nvidia FleX】(2) Position Based Fluids - 知乎 (zhihu.com)

效果如下：

流体模拟概述

计算机按照离散形式处理数据，而现实世界中的流体不是一帧一帧地运动，而是连续的。如何用计算机来处理现实世界中连续运动的流体是流体模拟首先要解决的问题之一。解决这个问题的方法主要有两个：一是拉格朗日方法，二是欧拉方法。

拉格朗日方法将流体视为一个个独立的粒子，每个粒子均携带足够多的物理量，通过在每一帧模拟每个粒子的运动实现对流体的模拟。

欧拉方法则是对空间进行划分，针对空间进行解算，在每一帧计算流体的运动。

多数现有的流体模拟算法使用在模拟水时使用拉格朗日方法，在模拟烟雾时则使用欧拉方法。

本文所提到的流体模拟算法是一种拉格朗日方法，在介绍PBF之前，先介绍与其密切相关的SPH算法。

光滑粒子动力学

光滑粒子动力学（Smoothed Particle Hydrodynamics，SPH）是经典的流体模拟方法。其核心思想是逐粒子计算。“对单个粒子附近区域（半径为r的球型区域）寻找若干邻居粒子，通过一套精妙的函数计算出每个邻居粒子对当前粒子的贡献，以求出当前粒子的部分物理量（例如密度、压强等）。得到了这些物理量以后我们再通过另一套精妙的函数来计算这些粒子在这帧得到的外力，从而推演出这一帧该粒子的受力情况，以及位置、速度、加速度的变化情况。”（这段话讲解的简单透彻，直接引用自这篇文章基于粒子表示的实时流体模拟与渲染（上） - 知乎 (zhihu.com)）。

SPH基于光滑模型（Smooth Model），所谓光滑模型可以在这样一个场景下理解：空间中有一组粒子，每个粒子都携带了一组物理量，为了估计空间中任意一点的物理量需要对该点周围粒子的物理量做插值。

上文提到的精妙的函数即对这些物理量做插值的函数。通常而言，该函数可写为： \[ A_i^{smooth}=\frac{1}{n}\sum_{j}V_jA_jW_{ij} \] 其中 $V_j$ 为第 $j$ 个粒子的体积，$A_j$ 为第 $j$个粒子的物理量，$W_{ij}$ 即Smoothing Kernel Function（光滑核函数），也可称为光滑核。

正如上文所述，粒子的物理量（如密度、压强）是根据粒子领域的邻居粒子计算得到的。而在计算这些物理量时，邻居粒子对中心粒子物理量的贡献随距离增加而减小。这种随距离而衰减的函数被称为“光滑核函数”，最大影响半径为“光滑核半径”。

光滑核函数主要具有以下几个性质：

归一性：光滑核函数在支持域上的积分为1；
紧支域：光滑核函数在支持域外的值为0；
收敛性：光滑核函数在支持域内任意一点的值不小于0；
衰减性：光滑核半径增大时，光滑核函数的值应当是单调递减的；
对称性：光滑核函数应为偶函数；
光滑性：光滑核函数应当连续。

流体模拟中主要涉及的3个力为重力、压力和粘滞力。重力加速度是个常量，所以重力易于计算。而压力和粘滞力则需要借助上文提到的光滑模型。

压力即由压强差产生的力。密度高的位置压强大，若一个粒子四周的压强不均等，则会产生压强差，此时粒子便会受到压力的作用。因此，为了计算压强差，需要对粒子的压强求梯度。 \[ F_i^{Pressure}=-V_i\nabla P^{smooth} \] 其中 $V_i$ 是粒子体积，体积越大，粒子收到的压力越大，$\nabla P^{smooth}$ 是对压强求梯度，即计算压强差。

粘滞力产生的效果是让粒子朝同一个方向移动，即减小粒子之间的速度差。 \[ F_i^{Viscosity}=-vm_i\triangle_iV^{smooth} \] 其中，$v$ 是粘滞系数，$m_i$ 是粒子质量，$\triangle_iV^{smooth}$ 是拉普拉斯算子，该项实际上是在计算一个点与周围邻居的差值。

有了上述基础，SPH算法中，更新每个粒子位置的基本流程可以总结为：

计算粒子周围邻居；
计算重力、压力、粘滞力；
更新粒子速度；
更新粒子位置。

Position Based Dynamic

传统的物理模拟常常使用基于力的方法（Force Based），这类方法通过计算物体受力，求解物体加速度，进一步更新物体速度、位置等。对于基于力的方法而言，时间步长越小，数值越稳定，但带来的计算量也就更大。基于力的方法针对每一种动态物体，会有一个独立的解算器（Solver），各种Solver按照一定顺序进行计算，从而得到模拟结果，这样会带来大量冗余工作。基于位置动力学（Position Based Dynamic，PBD）则使用同一种方式表达动态物体，这种方法只需要一个Solver便可以模拟所有物体，并支持它们之间的双向交互效果。PBD直接处理位移，通过求解满足约束函数，来调整每个物体的位置。

常见的约束有：

距离约束，用于模拟布料，包括布料的拉伸和弯曲等行为；
形状约束，用于模拟刚体和软体；
密度约束，用于模拟流体；
体积约束，用于模拟充气体；
碰撞约束，计算粒子之间的碰撞反应，比如把粒子从穿透状态解开，计算摩擦力等等

概述

动力学物体由 $N$ 个顶点和 $M$ 个约束组成。

顶点 $i\in{[1,2,...,N]}$ 的质量为 $m_i$，位置为 $x_i$，速度为 $v_i$。

约束 $j\in{[1,2,...,N]}$ 有五种特性：

约束的基数 $n_j$ ，即第 $j$ 个约束所影响的顶点数；
约束是实值函数；
约束索引值为 $\{i_1, i_2, ..., i_{n_j}\}, i_k\in[1, 2, ..,N]$，即影响到的顶点；
每个约束都有对应的刚度参数 $k_j\in[0, ..., 1]$ ，可以理解为约束强度；
约束分为：
1. 等式约束$C_j(x_{i_1},...,x_{i_{n_j}})=0$；
2. 不等式约束$C_j(x_{i_1},...,x_{i_{n_j}})\geq0$：

PBD基本流程

PBD的基本流程如下：

初始化顶点位置、速度和质量倒数；
循环求解：
1. 用转换为位置约束的力（如重力）更新顶点的速度；
2. 对速度添加阻尼；
3. 根据速度计算位置，作为位置的预测值；
4. 为每个顶点生成碰撞约束；
5. 对约束进行迭代求解，校正粒子位置；
6. 更新粒子的位置和速度信息；
7. 根据摩擦（friction）和恢复（restitution）系数修正速度。

简单概括就是根据当前速度预测顶点的位置，然后在满足约束求解的条件下，修正顶点的位置和速度。

Position Based Fluids

PBD一开始是为模拟固体而设计的，Macklin等人[1]结合SPH插值方法，提出了基于PBD框架的流体模拟方法Position Based Fluids（PBF）。PBF算法依然沿用了SPH的光滑模型，但是在求解每一帧粒子的位置和速度时采用了PBD的方案。

首先介绍流体的不可压缩性。上文提到，PBD方法通过密度约束来模拟流体，即通过密度约束来保证流体不可压缩。在SPH算法中，粒子的密度为： \[ \rho_i=\sum_j{V_j\rho_jW_{ij}}=\sum_j{\frac{m_j}{\rho_j}\rho_jW_{ij}}=\sum_jm_jW_{ij} \] 为保证流体不可压缩，需要让粒子 $i$ 的密度 $\rho_i$ 与静态密度 $\rho_0$ 尽量保持一致。因此，需要对每个粒子施加一个常量密度约束： \[ C_i(p_1, p_2, ..., p_n)=\frac{\rho_i}{\rho_0}-1 \] 将密度代入可得： \[ C_i=\frac{\sum_j m_j W_{ij}}{\rho_0}-1 \] PBF针对粒子质量相同时的情况，所以可以省去 \[m_j\] ，即： \[ C_i=\frac{\sum_jW_{ij}}{\rho_0}-1 \] 对 $C_i$ 求梯度可得： \[ \nabla_{pk}C_i=\frac{1}{\rho_0}\sum_j\nabla_{p_k}W_{ij} \] 这里要分两种情况： \[ \nabla_{pk}C_i=\frac{1}{\rho_0}\begin{cases} \sum_j\nabla_{p_k}W_{ij}& \text{k=i}\\ -\nabla_{p_k}W_{ij}& \text{k=j} \end{cases} \] 我的理解是，这里分情况没有什么特别的，到了下面求位移变化量的时候就看得明白了。

接下来，计算PBD所需要的拉格朗日乘子 $\lambda$ ： \[ \lambda_i = -\frac{C_i}{\sum_k|\nabla_{p_k}C_i|^2} \]

软约束

为了避免计算 $\lambda$ 时出现除 $0$ 错误，在分母上添加一个常量 $\epsilon$ 。这种在一定程度上允许违反约束条件的约束称为软约束（soft constraint），上式变为： \[ \lambda_i = -\frac{C_i}{\sum_k|\nabla_{p_k}C_i|^2+\epsilon} \] 经过上述约束投影后，粒子 $i$ 对应的位移向量 $\triangle p_i$ 为： \[ \triangle p_i=\lambda_i\nabla_{p_i}C_i+\sum_j\lambda_j\nabla_{p_j}C_j \\ =\frac{1}{\rho_0}\sum_j\lambda_i\nabla_{p_i}W_{ij}+(-\frac{1}{\rho_0}\sum_j\lambda_j\nabla_{p_j}W_{ij}) \\ =\frac{1}{\rho_0}\sum_j\lambda_i\nabla_{p_i}W_{ij}+ \frac{1}{\rho_0}\sum_j\lambda_j\nabla_{p_i}W_{ij} \\ =\frac{1}{\rho_0}\sum_j(\lambda_i+\lambda_j)\nabla_{p_i}W_{ij} \] 第二行的符号变化是因为把 $j$ 换成了 $i$ ，相当于反向了，所以符号取反。

Tensile Instability

另外，为了防止在邻居粒子不足的情况下，导致求出的流体密度低于静态密度，由此造成压强为负数，原本粒子间的压力变为吸引力，出现凝聚现象。使用了人工排斥力计算模型，当粒子距离过近时，添加排斥力使其分开，避免凝聚现象。这样，$\triangle p_i$ 的计算公式变为： \[ \triangle p_i=\frac{1}{\rho_0}\sum_j(\lambda_i+\lambda_j+S_{corr})\nabla_{p_i}W_{ij} \] 其中，$S_{corr}$ 由下式计算得来： \[ S_{corr}=-k(\frac{W_{ij}}{Const})^n \] 上式中的 $ Const $ 是根据光滑核函数 $ W_{poly6}$ 计算得到的一个固定值，距离距离球心的距离一般取 $0.1$ 到 $0.3$ ；$k=0.1$ 可以视为表面张力参数；$n=4$ 。

Vorticity Confinement

PBD方法引入了阻尼，会导致额外的能量损失，导致本应存在的涡流消失。为了补充这部分能量，PBF引入Vorticity Confinement向系统注入能量。 \[ f_i^{vorticity}=\epsilon(N\times \omega_i)\\ \] 其中，$N=\frac{\eta}{|\eta|}$ ，$\eta=\nabla|\omega_i|$，粒子 $i$ 的旋度 $\omega_i$ 为： \[ \omega_i=\nabla \times v=\sum_j(v_j-v_i) \times \nabla_{p_j} W_{ij} \]

实际计算过程中，需要先计算出每个粒子的 $\rho_i, \lambda_i$，然后一边计算与邻居粒子的在光滑核函数下的值，和$S_{corr}$ 项，一边累加$\triangle p_i$。然后，通过Vorticity Confinement向系统补充能量。最后更新还要粒子速度和位置。对于固体粒子，也就是容器，按照相同的方式处理。

D3D上的流体模拟

为了在D3D平台上实现流体模拟，这里借用了龙书给出的最基本的项目结构，在其中添加了几个类和Shader。主要涉及到这样几个部分：

SPH.cpp：C++类，构造流体、虚粒子容器、设置流体模拟所需的一些常量（如重力、流体密度、虚拟网格尺寸）
Grid.hlsl：Computer Shader，根据每个粒子的位置计算其所在的网格
Sort.hlsl：Computer Shader，用于实现双调排序，加速粒子索引速度
PBD：Computer Shader，位置预测、求解约束函数、位置更新等操作。

HLSL中涉及到这样几个结构体：

cbuffer cbSPHConstants : register(b0)
{
	float HashScale;
	float3 HashSize; 
	float3 HashTranslate; 
	float KernelRadius;
	float KernelRadius2;

	float POLY6ZERO;
	float POLY6;
	float SPIKY_GRAD;
	float VISC_LAP;

	float3 gravity;
	float damping; 
	float density0;
	float viscosity;
	float fluidmass;
	int NumParticles;
	int NumBoundaryParticles;

	int solverIterations;
	float deltaTime;
	float invDensity0;
	
	int TotalParticles;
};


struct Particle
{
	float density;
	float3 velocity;
	float3 position;
	float3 predictPosition;
	float lambda;
	float4 color;
};

通过设置根签名，规定数据传递的槽。渲染过程中，在不同的槽设置不同的数据：

void VecAddCSApp::BuildRootSignature()
{
	
	CD3DX12_ROOT_PARAMETER slotRootParameter[8];

	slotRootParameter[0].InitAsUnorderedAccessView(0);
	slotRootParameter[1].InitAsUnorderedAccessView(1);
	slotRootParameter[2].InitAsUnorderedAccessView(2);
	slotRootParameter[3].InitAsShaderResourceView(0);
	slotRootParameter[4].InitAsShaderResourceView(1);
	slotRootParameter[5].InitAsConstantBufferView(0);
	slotRootParameter[6].InitAsUnorderedAccessView(3);
	slotRootParameter[7].InitAsConstants(4, 1);

	CD3DX12_ROOT_SIGNATURE_DESC rootSigDesc(8, slotRootParameter,
		0, nullptr,
		D3D12_ROOT_SIGNATURE_FLAG_NONE);
	...
}

设置HLSL函数：


void VecAddCSApp::BuildShadersAndInputLayout()
{

	mShaders["standardVS"] = d3dUtil::CompileShader(L"Shaders\\color.hlsl", nullptr, "VS", "vs_5_0");
	mShaders["opaquePS"] = d3dUtil::CompileShader(L"Shaders\\color.hlsl", nullptr, "PS", "ps_5_0");

	mInputLayout =
	{
		{ "POSITION", 0, DXGI_FORMAT_R32G32B32_FLOAT, 0, 0, D3D12_INPUT_CLASSIFICATION_PER_VERTEX_DATA, 0 },
		{ "COLOR", 0, DXGI_FORMAT_R32G32B32A32_FLOAT, 0, 12, D3D12_INPUT_CLASSIFICATION_PER_VERTEX_DATA, 0 }
	};

	mShaders["UpdateIndexMapCS"] = d3dUtil::CompileShader(L"Shaders\\Grid.hlsl", nullptr, "UpdateIndexMap", "cs_5_0");
	mShaders["ClearTableCS"] = d3dUtil::CompileShader(L"Shaders\\Grid.hlsl", nullptr, "ClearTable", "cs_5_0");
	mShaders["UpdateTableCS"] = d3dUtil::CompileShader(L"Shaders\\Grid.hlsl", nullptr, "UpdateTable", "cs_5_0");

	// BitonicSort
	mShaders["BitonicSortCS"] = d3dUtil::CompileShader(L"Shaders\\Sort.hlsl", nullptr, "BitonicSort", "cs_5_0");
	mShaders["MatrixTransposeCS"] = d3dUtil::CompileShader(L"Shaders\\Sort.hlsl", nullptr, "MatrixTranspose", "cs_5_0");
	mShaders["FillCS"] = d3dUtil::CompileShader(L"Shaders\\Sort.hlsl", nullptr, "Fill", "cs_5_0");
	mShaders["CopyCS"] = d3dUtil::CompileShader(L"Shaders\\Sort.hlsl", nullptr, "Copy", "cs_5_0");

	// SPH
	mShaders["PredictPosCS"] = d3dUtil::CompileShader(L"Shaders\\PDB2.hlsl", nullptr, "PredictPos", "cs_5_0");
	mShaders["ComputeDensityCS"] = d3dUtil::CompileShader(L"Shaders\\PDB2.hlsl", nullptr, "ComputeDensity", "cs_5_0");
	mShaders["SolveConstraintCS"] = d3dUtil::CompileShader(L"Shaders\\PDB2.hlsl", nullptr, "SolveConstraint", "cs_5_0");
	mShaders["UpdateVelocitiesCS"] = d3dUtil::CompileShader(L"Shaders\\PDB2.hlsl", nullptr, "UpdateVelocities", "cs_5_0");
	mShaders["SolveViscosityCS"] = d3dUtil::CompileShader(L"Shaders\\PDB2.hlsl", nullptr, "SolveViscosity", "cs_5_0");
	mShaders["UpdatePositionsCS"] = d3dUtil::CompileShader(L"Shaders\\PDB2.hlsl", nullptr, "UpdatePositions", "cs_5_0");
	mShaders["ComputePsiCS"] = d3dUtil::CompileShader(L"Shaders\\PDB2.hlsl", nullptr, "ComputePsi", "cs_5_0");

}

当然，后续还要设置PSO，这里不再赘述。

最为重要的逻辑在Draw函数中：

void VecAddCSApp::Draw(const GameTimer& gt)
{
	
	// 添加计算着色器代码
	// Reuse the memory associated with command recording.
	// We can only reset when the associated command lists have finished execution on the GPU.
	ThrowIfFailed(mDirectCmdListAlloc->Reset());

	auto sphConstants = mCurrFrameResource->mSPHConstants->Resource();

	// A command list can be reset after it has been added to the command queue via ExecuteCommandList.
	// Reusing the command list reuses memory.
	//（1）填充IndexMap
	ThrowIfFailed(mCommandList->Reset(mDirectCmdListAlloc.Get(), mPSOs["UpdateIndexMapCS"].Get()));
	mCommandList->SetComputeRootSignature(mRootSignature.Get());

	mCommandList->SetComputeRootUnorderedAccessView(0, mIndexMap->GetGPUVirtualAddress());
	mCommandList->SetComputeRootUnorderedAccessView(1, mTable->GetGPUVirtualAddress());

	mCommandList->SetComputeRootShaderResourceView(3, mFluidParticles_Buffer->GetGPUVirtualAddress());
	mCommandList->SetComputeRootShaderResourceView(4, mBoundaryParticles_Buffer->GetGPUVirtualAddress());

	mCommandList->SetComputeRootConstantBufferView(5, sphConstants->GetGPUVirtualAddress());
	mCommandList->SetComputeRootUnorderedAccessView(6, mDebug->GetGPUVirtualAddress());

	int dispatchX = (mSPH.get()->mSPHConstants.TotalParticles / 128) + 1;
	mCommandList->Dispatch(dispatchX, 1, 1);

	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::Transition(mIndexMap.Get(),
		D3D12_RESOURCE_STATE_UNORDERED_ACCESS, D3D12_RESOURCE_STATE_COMMON));
	//mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mIndexMap.Get()));
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mTable.Get()));
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mDebug.Get()));
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::Transition(mFluidParticles_Buffer.Get(),
		D3D12_RESOURCE_STATE_COMMON, D3D12_RESOURCE_STATE_UNORDERED_ACCESS));
	
	
	//（2）排序
	int x, y, z;
	mCommandList->SetPipelineState(mPSOs["FillCS"].Get());
	int level = 0;
	int levelMask = 0;
	int width = mSPH.get()->mSort->Count;
	int height = 0;

	mCommandList->SetComputeRoot32BitConstants(7, 1, &level, 0);
	mCommandList->SetComputeRoot32BitConstants(7, 1, &levelMask, 1);
	mCommandList->SetComputeRoot32BitConstants(7, 1, &width, 2);
	mCommandList->SetComputeRoot32BitConstants(7, 1, &height, 3);
	// Data  : mGPUSortBuffer1
	// Input : mIndexMap
	mCommandList->SetComputeRootUnorderedAccessView(0, mGPUSortBuffer1->GetGPUVirtualAddress());
	mCommandList->SetComputeRootShaderResourceView(3, mIndexMap->GetGPUVirtualAddress());
	x = mSPH.get()->mSort->NumElements / mSPH.get()->mSort->THREADS;
	mCommandList->Dispatch(x, 1, 1);
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mGPUSortBuffer1.Get()));
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mDebug.Get()));
	
	height = mSPH.get()->mSort->MATRIX_HEIGHT;
	width = mSPH.get()->mSort->MATRIX_WIDTH;
	mCommandList->SetComputeRoot32BitConstants(7, 1, &width, 2);
	mCommandList->SetComputeRoot32BitConstants(7, 1, &height, 3);
	for (int l = 2; l <= mSPH.get()->mSort->BITONIC_BLOCK_SIZE; l = l * 2)
	{
		mCommandList->SetPipelineState(mPSOs["BitonicSortCS"].Get());
		// Sort the row data
		level = l;
		levelMask = level;
		mCommandList->SetComputeRoot32BitConstants(7, 1, &level, 0);
		mCommandList->SetComputeRoot32BitConstants(7, 1, &levelMask, 1);
		x = mSPH.get()->mSort->NumElements / mSPH.get()->mSort->BITONIC_BLOCK_SIZE;
		mCommandList->Dispatch(x, 1, 1);
		//添加资源屏障
		mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mGPUSortBuffer1.Get()));
		mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mDebug.Get()));
	}
	
	// Then sort the rows and columns for the levels > than the block size
	// Transpose. Sort the Columns. Transpose. Sort the Rows.
	for (int l = (mSPH.get()->mSort->BITONIC_BLOCK_SIZE * 2); l <= mSPH.get()->mSort->NumElements; l = l * 2)
	{
		//（a）Transpose the data from buffer 1 into buffer 2
		mCommandList->SetPipelineState(mPSOs["MatrixTransposeCS"].Get());
		level = l / mSPH.get()->mSort->BITONIC_BLOCK_SIZE;
		levelMask = (l & ~mSPH.get()->mSort->NumElements) / mSPH.get()->mSort->BITONIC_BLOCK_SIZE;
		mCommandList->SetComputeRoot32BitConstants(7, 1, &level, 0);
		mCommandList->SetComputeRoot32BitConstants(7, 1, &levelMask, 1);
		mCommandList->SetComputeRoot32BitConstants(7, 1, &width, 2);
		mCommandList->SetComputeRoot32BitConstants(7, 1, &height, 3);
		// 转换状态
		mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::Transition(mGPUSortBuffer1.Get(),
			D3D12_RESOURCE_STATE_UNORDERED_ACCESS, D3D12_RESOURCE_STATE_COMMON));
		// Data  : mGPUSortBuffer2：UAV
		// Input : mGPUSortBuffer1: SRV
		mCommandList->SetComputeRootUnorderedAccessView(0, mGPUSortBuffer2->GetGPUVirtualAddress());
		mCommandList->SetComputeRootShaderResourceView(3, mGPUSortBuffer1->GetGPUVirtualAddress());

		x = mSPH.get()->mSort->MATRIX_WIDTH / mSPH.get()->mSort->TRANSPOSE_BLOCK_SIZE;
		y = mSPH.get()->mSort->MATRIX_HEIGHT / mSPH.get()->mSort->TRANSPOSE_BLOCK_SIZE;
		mCommandList->Dispatch(x, y, 1);
		// 资源屏障
		mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mGPUSortBuffer2.Get()));
		mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mDebug.Get()));

		// break;
		//（b）Sort the transposed column data
		mCommandList->SetPipelineState(mPSOs["BitonicSortCS"].Get());
		mCommandList->SetComputeRootUnorderedAccessView(0, mGPUSortBuffer2->GetGPUVirtualAddress());
		mCommandList->Dispatch(mSPH.get()->mSort->NumElements / mSPH.get()->mSort->BITONIC_BLOCK_SIZE, 1, 1);
		// 资源屏障
		mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mGPUSortBuffer2.Get()));
		mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mDebug.Get()));

		// （c）Transpose the data from buffer 2 back into buffer 1
		mCommandList->SetPipelineState(mPSOs["MatrixTransposeCS"].Get());
		level = mSPH.get()->mSort->BITONIC_BLOCK_SIZE;
		levelMask = l;
		mCommandList->SetComputeRoot32BitConstants(7, 1, &level, 0);
		mCommandList->SetComputeRoot32BitConstants(7, 1, &levelMask, 1);
		mCommandList->SetComputeRoot32BitConstants(7, 1, &height, 2);
		mCommandList->SetComputeRoot32BitConstants(7, 1, &width, 3);
		// 转换状态
		mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::Transition(mGPUSortBuffer1.Get(),
			D3D12_RESOURCE_STATE_COMMON, D3D12_RESOURCE_STATE_UNORDERED_ACCESS));
		mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::Transition(mGPUSortBuffer2.Get(),
			D3D12_RESOURCE_STATE_UNORDERED_ACCESS, D3D12_RESOURCE_STATE_COMMON));
		// Data  : mGPUSortBuffer1: UAV
		// Input : mGPUSortBuffer2: SRV
		mCommandList->SetComputeRootUnorderedAccessView(0, mGPUSortBuffer1->GetGPUVirtualAddress());
		mCommandList->SetComputeRootShaderResourceView(3, mGPUSortBuffer2->GetGPUVirtualAddress());
		x = mSPH.get()->mSort->MATRIX_HEIGHT / mSPH.get()->mSort->TRANSPOSE_BLOCK_SIZE;
		y = mSPH.get()->mSort->MATRIX_WIDTH / mSPH.get()->mSort->TRANSPOSE_BLOCK_SIZE;
		mCommandList->Dispatch(x, y, 1);
		// 资源屏障
		mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mGPUSortBuffer1.Get()));
		mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mDebug.Get()));

		// （d）Sort the row data
		mCommandList->SetPipelineState(mPSOs["BitonicSortCS"].Get());
		mCommandList->SetComputeRootUnorderedAccessView(0, mGPUSortBuffer1->GetGPUVirtualAddress());
		mCommandList->Dispatch(mSPH.get()->mSort->NumElements / mSPH.get()->mSort->BITONIC_BLOCK_SIZE, 1, 1);
		// 资源屏障
		mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mGPUSortBuffer1.Get()));
		mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mDebug.Get()));

		// 收尾。转换状态，以便下一次循环正常
		mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::Transition(mGPUSortBuffer2.Get(),
			D3D12_RESOURCE_STATE_COMMON, D3D12_RESOURCE_STATE_UNORDERED_ACCESS));
	}
	mCommandList->SetPipelineState(mPSOs["CopyCS"].Get());
	width = mSPH.get()->mSort->Count;
	mCommandList->SetComputeRoot32BitConstants(7, 1, &width, 2);
	// 转换状态
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::Transition(mIndexMap.Get(),
		D3D12_RESOURCE_STATE_COMMON, D3D12_RESOURCE_STATE_UNORDERED_ACCESS));
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::Transition(mGPUSortBuffer1.Get(),
		D3D12_RESOURCE_STATE_UNORDERED_ACCESS, D3D12_RESOURCE_STATE_COMMON));
	// Data  : mIndexMap
	// Input : mGPUSortBuffer1
	mCommandList->SetComputeRootUnorderedAccessView(0, mIndexMap->GetGPUVirtualAddress());
	mCommandList->SetComputeRootShaderResourceView(3, mGPUSortBuffer1->GetGPUVirtualAddress());
	mCommandList->Dispatch(mSPH.get()->mSort->NumElements / mSPH.get()->mSort->THREADS, 1, 1);
	// 资源屏障
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mIndexMap.Get()));
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mDebug.Get()));
	//（3）
	mCommandList->SetPipelineState(mPSOs["ClearTableCS"].Get());
	mCommandList->SetComputeRootUnorderedAccessView(1, mTable->GetGPUVirtualAddress());
	mCommandList->Dispatch(dispatchX, 1, 1);
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mTable.Get()));
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mDebug.Get()));

	//（4）
	mCommandList->SetPipelineState(mPSOs["UpdateTableCS"].Get());
	mCommandList->SetComputeRootUnorderedAccessView(0, mIndexMap->GetGPUVirtualAddress());
	mCommandList->SetComputeRootUnorderedAccessView(1, mTable->GetGPUVirtualAddress());
	mCommandList->Dispatch(dispatchX, 1, 1);
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mIndexMap.Get()));
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mTable.Get()));
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mDebug.Get()));

	// (5) PDB
	// (5-0) ComputePsi
	mCommandList->SetPipelineState(mPSOs["ComputePsiCS"].Get());
	mCommandList->SetComputeRootUnorderedAccessView(0, mFluidParticles_Buffer->GetGPUVirtualAddress());
	mCommandList->SetComputeRootUnorderedAccessView(1, mBoundaryParticles_Buffer->GetGPUVirtualAddress());
	mCommandList->SetComputeRootUnorderedAccessView(6, mDebug->GetGPUVirtualAddress());
	mCommandList->SetComputeRootShaderResourceView(3, mIndexMap->GetGPUVirtualAddress());
	mCommandList->SetComputeRootShaderResourceView(4, mTable->GetGPUVirtualAddress());
	mCommandList->SetComputeRootConstantBufferView(5, sphConstants->GetGPUVirtualAddress());
	mCommandList->Dispatch(dispatchX, 1, 1);
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mBoundaryParticles_Buffer.Get()));
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mDebug.Get()));


	// (5-a) PredictPos
	mCommandList->SetPipelineState(mPSOs["PredictPosCS"].Get());
	mCommandList->Dispatch(dispatchX, 1, 1);
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mFluidParticles_Buffer.Get()));
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mDebug.Get()));

	// (5-b) ComputeDensity
	mCommandList->SetPipelineState(mPSOs["ComputeDensityCS"].Get());
	mCommandList->Dispatch(dispatchX, 1, 1);
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mFluidParticles_Buffer.Get()));
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mDebug.Get()));

	
	// (5-c) SolveConstraint
	mCommandList->SetPipelineState(mPSOs["SolveConstraintCS"].Get());
	mCommandList->Dispatch(dispatchX, 1, 1);
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mFluidParticles_Buffer.Get()));
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mDebug.Get()));

	// (3-4) UpdateVelocities
	mCommandList->SetPipelineState(mPSOs["UpdateVelocitiesCS"].Get());
	mCommandList->Dispatch(dispatchX, 1, 1);
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mFluidParticles_Buffer.Get()));
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mDebug.Get()));

	// (3-5) SolveViscosity
	mCommandList->SetPipelineState(mPSOs["SolveViscosityCS"].Get());
	mCommandList->Dispatch(dispatchX, 1, 1);
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mFluidParticles_Buffer.Get()));
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mDebug.Get()));

	// (3-6) UpdatePositions
	mCommandList->SetPipelineState(mPSOs["UpdatePositionsCS"].Get());
	mCommandList->Dispatch(dispatchX, 1, 1);
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mFluidParticles_Buffer.Get()));
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::UAV(mDebug.Get()));
	

	// Schedule to copy the data to the default buffer to the readback buffer.
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::Transition(mFluidParticles_Buffer.Get(),
		D3D12_RESOURCE_STATE_UNORDERED_ACCESS, D3D12_RESOURCE_STATE_COPY_SOURCE));
	mCommandList->CopyResource(mReadBackBuffer.Get(), mFluidParticles_Buffer.Get());
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::Transition(mFluidParticles_Buffer.Get(),
		D3D12_RESOURCE_STATE_COPY_SOURCE, D3D12_RESOURCE_STATE_UNORDERED_ACCESS));
	
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::Transition(mIndexMap.Get(),
		D3D12_RESOURCE_STATE_UNORDERED_ACCESS, D3D12_RESOURCE_STATE_COPY_SOURCE));
	mCommandList->CopyResource(mIndexMapBack.Get(), mIndexMap.Get());
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::Transition(mIndexMap.Get(),
		D3D12_RESOURCE_STATE_COPY_SOURCE, D3D12_RESOURCE_STATE_UNORDERED_ACCESS));

	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::Transition(mDebug.Get(),
		D3D12_RESOURCE_STATE_UNORDERED_ACCESS, D3D12_RESOURCE_STATE_COPY_SOURCE));
	mCommandList->CopyResource(mDebugBack.Get(), mDebug.Get());
	mCommandList->ResourceBarrier(1, &CD3DX12_RESOURCE_BARRIER::Transition(mDebug.Get(),
		D3D12_RESOURCE_STATE_COPY_SOURCE, D3D12_RESOURCE_STATE_UNORDERED_ACCESS));


	// Done recording commands.
	ThrowIfFailed(mCommandList->Close());
	// Add the command list to the queue for execution.
	ID3D12CommandList* cmdsLists1[] = { mCommandList.Get() };
	mCommandQueue->ExecuteCommandLists(_countof(cmdsLists1), cmdsLists1);
	// Wait for the work to finish.
	FlushCommandQueue();

	...

}

最近忙毕业，时间太紧张，QAQ，后续会将详细内容补充上来……

[1] Macklin, Miles, and Matthias Müller. "Position based fluids." ACM Transactions on Graphics (TOG) 32.4 (2013): 104.

Simplify

流体模拟——Position Based Fluids