《Fluid Engine Development》 学习笔记4-预测校订不可压缩SPH-PCISPH
传统SPH方案的主要问题之一是时间步长限制。在原始的SPH中，咱们首先从当前设置计算密度，使用EOS计算压强，应用压力梯度，而后运行时间积分。这个过程意味着只须要必定的压缩量就能够触发内核半径内的压力，从而延迟计算。所以，咱们须要使用更小的时间步长(意味着更多的迭代)，这在计算上是昂贵的。或者，咱们可使用不那么严格的EOS，然而，这个解决方案可能会引入相似弹簧的振荡。微调参数如声速或粘度能够帮助避免此类问题。然而，这并非一个基本的解决方案，对用户来讲也是不切实际的。Solenthaler 和Pajarola经过在SPH模拟中引入预测-校订器概念来解决这个问题。这种又称为预测校订不可压缩SPH（PCISPH），它是一种偏差测量算法，假定测量值和指望密度的差值是偏差。本篇文章总结我在《Fluid Engine Development》学到关于PCISPH的知识。算法
算法原理
1.计算协力，预测速度位置，碰撞检测或碰撞处理数组
2.计算密度后测量密度偏差，利用密度偏差更新压强函数
3.计算新压强梯度力ui
4.屡次重复以上过程，获得使密度偏差最小的修正力（压强梯度力），而后使用累积力进行下一步spa
因为它是经过累积校订力使最终状态处于不可压缩状态，而不是经过多个SPH步骤保证不可压缩，因此它不须要在乎时间步长的问题code
算法代码实现结构以下，ci
void CalfFluidEngine::PCISPHSolver3::accumulatePressureForce(double timeStepInSeconds)
{
auto particles = GetSphData();
const size_t numberOfParticles = particles->GetNumberOfParticles();
const double delta = computeDelta(timeStepInSeconds);
const double targetDensity = particles->GetDensity();
const double mass = particles->GetParticleMass();
auto& p = particles->GetPressures();
auto& d = particles->GetDensities();
auto& x = particles->GetPositions();
auto& v = particles->GetVelocities();
auto& f = particles->GetForces();
//Initialize
for (unsigned int k = 0; k < _maxNumberOfIterations; ++k)
{
// Predict velocity and position
// Resolve collisions
// Compute pressure from density error
// Compute pressure gradient force
// Compute max density error
double maxDensityError = ......;
double densityErrorRatio = maxDensityError / targetDensity;
if (std::fabs(densityErrorRatio) < _maxDensityErrorRatio){
break;
}
}
//Accumlate pressure force
}
咱们能够看到预测-校订自己是一个循环体，循环在密度偏差小于指定阈值循环执行次数小于指定迭代次数时终止get
而循环体内所作的就是不断进行修正压力，直到密度偏差足够小数学
须要注意的是在执行这套算法前，粒子所受其余力包括粘度，重力已经所有计算完成，已所有累加到forces数组中it
void accumulateForces(double timeIntervalInSeconds)
{
ParticleSystemSolver3::accumulateForces(timeIntervalInSeconds);
accumulateViscosityForce();
accumulatePressureForce(timeIntervalInSeconds);
}
这边给出完整的算法实现代码
void CalfFluidEngine::PCISPHSolver3::accumulatePressureForce(double timeStepInSeconds)
{
auto particles = GetSphData();
const size_t numberOfParticles = particles->GetNumberOfParticles();
const double delta = computeDelta(timeStepInSeconds);
const double targetDensity = particles->GetDensity();
const double mass = particles->GetParticleMass();
auto& p = particles->GetPressures();
auto& d = particles->GetDensities();
auto& x = particles->GetPositions();
auto& v = particles->GetVelocities();
auto& f = particles->GetForces();
//Initialize
std::vector<double> predictedDensities(numberOfParticles, 0.0);
SphStandardKernel3 kernel(particles->GetKernelRadius());
tbb::parallel_for(
tbb::blocked_range<size_t>(0, numberOfParticles),
[&](const tbb::blocked_range<size_t> & b) {
for (size_t i = b.begin(); i != b.end(); ++i)
{
p[i] = 0.0;
_pressureForces[i] = Vector3D::zero;
_densityErrors[i] = 0.0;
predictedDensities[i] = d[i];
}
});
for (unsigned int k = 0; k < _maxNumberOfIterations; ++k)
{
// Predict velocity and position
tbb::parallel_for(
tbb::blocked_range<size_t>(0, numberOfParticles),
[&](const tbb::blocked_range<size_t> & b) {
for (size_t i = b.begin(); i != b.end(); ++i)
{
_tempVelocities[i] = v[i] +
(f[i] + _pressureForces[i]) / mass * timeStepInSeconds;
_tempPositions[i] = x[i] +
_tempVelocities[i] * timeStepInSeconds;
}
});
// Resolve collisions
ParticleSystemSolver3::resolveCollision(_tempPositions, _tempVelocities);
// Compute pressure from density error
tbb::parallel_for(
tbb::blocked_range<size_t>(0, numberOfParticles),
[&](const tbb::blocked_range<size_t> & b) {
for (size_t i = b.begin(); i != b.end(); ++i)
{
double weightSum = 0.0;
const auto& neighbors = particles->GetNeighborLists()[i];
for (size_t j : neighbors) {
double dist = Vector3D::Distance(_tempPositions[j], _tempPositions[i]);
weightSum += kernel(dist);
}
weightSum += kernel(0);
double density = mass * weightSum;
double densityError = (density - targetDensity);
double pressure = delta * densityError;
if (pressure < 0.0) {
pressure *= _negativePressureScale;
densityError *= _negativePressureScale;
}
p[i] += pressure;
predictedDensities[i] = density;
_densityErrors[i] = densityError;
}
});
// Compute pressure gradient force
tbb::parallel_for(
tbb::blocked_range<size_t>(0, numberOfParticles),
[&](const tbb::blocked_range<size_t> & b) {
for (size_t i = b.begin(); i != b.end(); ++i)
{
_pressureForces[i] = Vector3D::zero;
}
});
SphSystemSolver3::accumulatePressureForce(x, predictedDensities, p, _pressureForces);
// Compute max density error
double maxDensityError = 0.0;
for (size_t i = 0; i < numberOfParticles; ++i) {
maxDensityError = AbsMax(maxDensityError, _densityErrors[i]);
}
double densityErrorRatio = maxDensityError / targetDensity;
if (std::fabs(densityErrorRatio) < _maxDensityErrorRatio){
break;
}
}
//Accumlate pressure force
tbb::parallel_for(
tbb::blocked_range<size_t>(0, numberOfParticles),
[&](const tbb::blocked_range<size_t> & b) {
for (size_t i = b.begin(); i != b.end(); ++i)
{
f[i] += _pressureForces[i];
}
});
}
乍一看算法实现并不难，大部分原理都已经在以前的笔记提过，惟一须要额外关注的PCISPH利用密度偏差更新压强的数学公式
压强更新公式
压强经过 \(p_{i}(t) += \delta \rho_{error}\)更新，这个公式是怎么来的呢，咱们来推导一下。
首先回顾一下密度的计算公式 \(\rho = \sum_{j} m_{j}W(x_{ij},h)\) 其中 \(W(x_{ij},h)\)是光滑核函数,\(x_{ij}是x_{i},x_{j}间的距离\)，由于粒子质量一致，因此 \(\rho = m\sum_{j} W(x_{ij})\)
设当前时间为t，则 $\rho_{t+\Delta t} = m\sum_{j} W(x_{i}(t+\Delta t) - x_{j}(t+\Delta t)) = m\sum_{j} W(x_{i}(t) +\Delta x_{i}(t) - x_{j}(t) - \Delta x_{j}(t) ) $
设 $r_{ij}(t) = x_{i}(t) - x_{j}(t),\Delta r_{ij} = \Delta x_{i}(t) - \Delta x_{j}(t) $, 则 $\rho_{t+\Delta t} = m\sum_{j} W(r_{ij} (t) + \Delta r_{ij}(t)) $
而后泰勒展开可得 \(\rho_{t+\Delta t} = m\sum_{j} W(r_{ij} (t) )+ \nabla W(r_{ij}(t)) \cdot \Delta r_{ij}(t) = \rho_{i}(t) + \Delta \rho_{i}(t)\)
可求得密度增量为$\Delta \rho_{i}(t) =m \sum_{j}\nabla W(r_{ij}(t)) \cdot \Delta r_{ij}(t)= m \sum_{j} \nabla W(r_{ij}(t)) \cdot (\Delta x_{i}(t) - \Delta x_{j}(t)) $
\(= m (\Delta x_{i}(t)\sum _{j}\nabla W_{ij}\)
$ - \sum_{j}\nabla W_{ij}\Delta x_{j}(t)) $
由于 $ \Delta x_{i} = \Delta t ^{2} \frac{F_{i}^{p}}{m}$
而压力\(f_{p}= m^{2} \sum_{j}(\frac{p_{i}}{\rho _{i} ^{2}} + \frac{p_{j}}{\rho _{j} ^{2}}) \nabla W(|x - x_{j}|)\),可得 \(F_{j =i}^p = m^2 \sum _{j } \frac{p_{i} + p_{i}}{\rho_{0}^{2}}\nabla W_{ij}\)
代入一下可得\(\Delta x_{i} = - \Delta t^{2} m \frac{2p_{i}}{\rho_{0}^2}\sum _{j} \nabla W_{ij}\)，\(\Delta x_{j} = - \Delta t^{2} m \frac{2p_{i}}{\rho_{0}^2}\nabla W_{ij}\)
把上式代入密度增量公式，可得 $\Delta \rho_{i}(t) =\Delta t^{2} m^{2} \frac{2p_{i}}{\rho_{0}^2}(-\sum_{j}\nabla W_{ij} \cdot \sum_{j}\nabla W_{ij} - \sum_{j}(\nabla W_{ij} \cdot \nabla W_{ij})) $
把上式转换一下，可得压强计算公式 \(p_{i} = \frac{\Delta \rho_{i}(t) }{(-\sum_{j}\nabla W_{ij} \cdot \sum_{j}\nabla W_{ij} - \sum_{j}(\nabla W_{ij} \cdot \nabla W_{ij})\beta }\)
这串公式能够这么理解，当但愿增量密度 \(\Delta \rho _{i}(t)\),须要施加压强pi
其中 \(\beta =\Delta t^{2} m^{2} \frac{2}{\rho_{0}^2}\)
设当前密度为 \(\rho_{i}^{*}\) 目标密度为 \(\rho_{0}\) 则\(\Delta \rho_{i}(t) = \rho_{0} - \rho_{i}^{*} = - \rho_{error}\)
设 \(p_{i} = \delta \rho_{error}\) 可得 \(\delta = \frac{-1}{(-\sum_{j}\nabla W_{ij} \cdot \sum_{j}\nabla W_{ij} - \sum_{j}(\nabla W_{ij} \cdot \nabla W_{ij}))\beta }\)
\(p_{i}(t) += \delta \rho_{error}\)
经过压强更新公式，咱们能够不断的在矫正粒子密度使其接近不可压缩条件
代码实现以下
double CalfFluidEngine::PCISPHSolver3::computeDelta(double timeStepInSeconds)
{
auto particles = GetSphData();
const double kernelRadius = particles->GetKernelRadius();
std::vector<Vector3D> points;
BccLatticePointGenerator pointsGenerator;
Vector3D origin = Vector3D::zero;
BoundingBox3D sampleBound(origin, origin);
sampleBound.Expand(1.5 * kernelRadius);
pointsGenerator.Generate(sampleBound, particles->GetTargetSpacing(), &points);
SphSpikyKernel3 kernel(kernelRadius);
double denom = 0;
Vector3D denom1 = Vector3D::zero;
double denom2 = 0;
for (size_t i = 0; i < points.size(); ++i) {
const Vector3D& point = points[i];
double distanceSquared = point.SquareMagnitude();
if (distanceSquared < kernelRadius * kernelRadius) {
double distance = std::sqrt(distanceSquared);
Vector3D direction =
(distance > 0.0) ? point / distance : Vector3D::zero;
// grad(Wij)
Vector3D gradWij = kernel.Gradient(distance, direction);
denom1 += gradWij;
denom2 += Vector3D::Dot(gradWij,gradWij);
}
}
denom += -Vector3D::Dot(denom1,denom1) - denom2;
return (std::fabs(denom) > 0.0) ?
-1 / (computeBeta(timeStepInSeconds) * denom) : 0;
}
double CalfFluidEngine::PCISPHSolver3::computeBeta(double timeStepInSeconds)
{
auto particles = GetSphData();
double t = particles->GetParticleMass() * timeStepInSeconds
/ particles->GetDensity();
return 2.0 * t * t;
}
这边再一次列出计算密度偏差的代码
// Compute pressure from density error
tbb::parallel_for(
tbb::blocked_range<size_t>(0, numberOfParticles),
[&](const tbb::blocked_range<size_t> & b) {
for (size_t i = b.begin(); i != b.end(); ++i)
{
double weightSum = 0.0;
const auto& neighbors = particles->GetNeighborLists()[i];
for (size_t j : neighbors) {
double dist = Vector3D::Distance(_tempPositions[j], _tempPositions[i]);
weightSum += kernel(dist);
}
weightSum += kernel(0);
double density = mass * weightSum;
double densityError = (density - targetDensity);
double pressure = delta * densityError;
if (pressure < 0.0) {
pressure *= _negativePressureScale;
densityError *= _negativePressureScale;
}
p[i] += pressure;
predictedDensities[i] = density;
_densityErrors[i] = densityError;
}
});
- 1. 《Fluid Engine Development》 学习笔记2-基础
- 2. 《Fluid Engine Development》 学习笔记1-求解线性方程组
- 3. 《Fluid Engine Development》 学习笔记3-光滑粒子流体动力学
- 4. Unreal Engine 4 学习笔记1：初识
- 5. Unreal Engine 4 笔记
- 6. webpack学习笔记--压缩代码
- 7. Linux学习笔记压缩与打包
- 8. Redis学习笔记（五） 压缩列表
- 9. Unity 贴图压缩学习笔记
- 10. Struts2学习笔记（4）--校验构架
- 更多相关文章...
- • 您已经学习了 XML Schema，下一步学习什么呢？ - XML Schema 教程
- • 我们已经学习了 SQL，下一步学习什么呢？ - SQL 教程
- • Tomcat学习笔记（史上最全tomcat学习笔记）
- • 适用于PHP初学者的学习线路和建议
-
每一个你不满意的现在，都有一个你没有努力的曾经。