【splishsplash】splishsplash源码算法剖析：粘性力算法

前面三个文章讲的是模拟的主流程。而splishsplash作为一个图形学算法开发人员的工具，最大的特点是方便植入和测试新算法，尤其是非压力梯度力算法。

入口

在TimeStepDFSPH.cpp的第125行：

	sim->computeNonPressureForces();

这就是非压力梯度力的入口。

F11步进，到Simulation.cpp的496行

fm->computeViscosity();

这里计算粘性力，其中fm是FluidModel类指针

F11再次步进，到FluidModel.cpp的439行，这里判断一下是不是选择了粘性力。如果没选，m_viscosity是个空指针。

    if (m_viscosity)
        m_viscosity->step();

再次步进
利用多态性，step函数就会跳转到你具体选择的粘性力模型的step函数处。

控制流

从main函数开始，到Viscosity_Weiler2018::step()函数

1. main()中调用base->run()
2. SimulatorBase::run()中调用runSimulation()
3. SimulatorBase::runSimulation()中调用m_gui->run()
4. Simulator_GUI_imgui::run()中调用MiniGL::mainLoop()
5. MiniGL::mainLoop()中调用idlefunc()
6. SimulatorBase::timeStep()中调用Simulation::getCurrent()->getTimeStep()->step()
7. TimeStepDFSPH::step()中调用sim->computeNonPressureForces();
8. Simulation::computeNonPressureForces()中调用fm->computeViscosity();
9. FluidModel::computeViscosity()中调用m_viscosity->step()
10. 调用Viscosity_Weiler2018::step()

粘性力的step函数

我们这里选择的粘性力模型是：Viscosity_Weiler2018

进入具体的某个粘性力的step函数以后，一般是先取出一些数据（如粒子数目，支持半径，时间步长等）存储到局部变量里。而取出的方法靠FluidModel类型的指针，一般叫做m_model。

m_model定义于NonPressureForceBase.h的第15行

	protected:
		FluidModel *m_model;

由于是protected类型，所有派生类都会继承这个成员变量。粘性力也属于NonPressureForceBase类的派生类。

粘性算法

所有的流体仿真问题最后都会变成Ax=b问题

算法分为三部分：

1. 准备矩阵A
2. 准备矩阵b
3. 利用Eigen求解Ax=b

为此我们首先要知道A和b分别代表什么。我们首先看A代表什么？

我们去看Weiler那篇文献： [WKBB18] Marcel Weiler, Dan Koschier, Magnus Brand,
and Jan Bender. A physically consistent implicit viscosity solver for
SPH fluids. Computer Graphics Forum (Eurographics), 2018. URL:
https://doi.org/10.1111/cgf.13349
首先是方程7，计算了拉普拉斯项。这一项乘以粘度mu就是NS方程中的粘性力。

把它时间积分得到方程8

方程8是计算下一时刻的速度的公式。这是要求解的核心公式，下面把它变形成Ax=b的形式。

于是方程9就是变形以后的方程8。

根据方程9，得到了A的表达式，即方程10。

其中下表ij表示粒子i和它的邻居j。mij就是mi和mj的简单平均。

下面是作者的具体讲解，这里懒得翻译了。

大意就是利用CG来求解。并且利用了Jacobi的Preconditioner来加速。还利用了initial
guess来猜测合理初值以加速求解。CG采用matrix free，即一个数一个数地算而不是一下子算一个向量。

我们回到代码

下面这一部分就是准备矩阵A，准备向量b，以及准备preconditioner和initial guess的

MatrixReplacement A(3*m_model->numActiveParticles(), matrixVecProd, (void*) this);
m_solver.preconditioner().init(m_model->numActiveParticles(), diagonalMatrixElement, (void*)this);

m_solver.setTolerance(m_maxError);
m_solver.setMaxIterations(m_maxIter);
m_solver.compute(A);

VectorXr b(3*numParticles);
VectorXr x(3*numParticles);
VectorXr g(3*numParticles);

computeRHS(b, g);

值得指出的是，这些部分都是利用Eigen算矩阵求解的八股文，所以有必要好好研究一下。

其中

MatrixReplacement A(3*m_model->numActiveParticles(), matrixVecProd, (void*) this);

是实例化了一个类MatrixReplacement ，这个类是定义在Eigen中的。它需要三个参数：

1. 粒子数*3
2. 一个函数指针matrixVecProd
3. this指针，指向的是Viscosity类的实例

传递matrixVecProd的这种方法叫函数回调。粗略地讲，函数回调就是做填空题。
我一开始不知道具体要调用的函数是谁。也许是模型A对应的函数，也许是模型B对应的函数。所以我挖个空。等我真正调用的时候，把函数名字填进去。这样就实现了多态性。

总之，这步是准备矩阵A。

其中，

m_solver.preconditioner().init(m_model->numActiveParticles(), diagonalMatrixElement, (void*)this);

是准备preconditioner，与前面类似diagonalMatrixElement也是个函数指针。

其中

m_solver.compute(A);

就是把preconditioner乘到A上，从而完成预处理。

其中

computeRHS(b, g);

就是计算b。

g就是那个initial guess的x0，所以这一步也做了initial guess。

那么，b代表什么呢？我们还是要回到论文。

这里指出来，b必须是和边界条件有关系的。

由于边界处于粒子缺失（deficientcy，SPH的固有缺陷之一）的特殊状态，粘性也要特殊处理。

他们采用的是Akinci2012的方法。 AKINCI N., IHMSEN M., AKINCI G., SOLENTHALER
B.,TESCHNER M.: Versatile rigid-fluid coupling for incompressible SPH.
ACM Trans. Graph. 31, 4 (2012), 1–8.

方程12中的这个b正是我们需要的b

我们现在讲完了准备步骤的1和2.下面就是套八股文用Eigen来计算Ax=b

	START_TIMING("CG solve");
	x = m_solver.solveWithGuess(b, g);
	m_iterations = (int)m_solver.iterations();
	STOP_TIMING_AVG;
	INCREASE_COUNTER("Visco iterations", static_cast<Real>(m_iterations));

其中

m_iterations = (int)m_solver.iterations();

这一句正是最核心的计算步骤。所有的实际计算在这一步中发生。

当算完以后，最后还差一步将计算好的粘性力加到加速度里。

applyForces(x);

总结一下，大部分都是套用Eigen的八股文，而我们实际关心的部分在于：

• matrixVecProd这个函数准备矩阵A
• computeRHS这个函数准备向量b
• applyForces这个函数将计算好的粘性力加到加速度里

我们现在分别讨论

matrixVecProd这个函数准备矩阵A

值得注意的是，matrixVecProd写了两种版本，一种是有AVX并行加速的（默认开启），一种是没有的。利用#ifdef USE_AVX这个宏来控制开关。
两种代码大同小异。

他们的核心都是这一行
（AVX版本）

delta_ai_avx += V_gradW * ((d_mu_rho0 / density_j_avx) * (vi_avx - vj_avx).dot(xixj) / (xixj.squaredNorm() + h2_001));

（无AVX版本）

ai += d * mu * (model->getMass(neighborIndex) / density_j) * (vi - vj).dot(xixj) / (xixj.squaredNorm() + 0.01*h2) * gradW;

将公式翻译过来就是
$$\nabla W\frac{(d\mu {\rho }_0/{\rho }_j)({\mathbf{v}}_i-{\mathbf{v}}_j)\cdot ({\mathbf{x}}_i-{\mathbf{x}}_j)}{||({\mathbf{x}}_i-{\mathbf{x}}_j)|{|}^2+0.01h^2}$$
对应论文公式10

computeRHS这个函数准备向量b

applyForces这个函数将计算好的粘性力加到加速度里