《Fluid Engine Development》 學習筆記4-預測校正不可壓縮SPH-PCISPH
傳統SPH方案的主要問題之一是時間步長限制。在原始的SPH中，我們首先從當前設定計算密度，使用EOS計算壓強，應用壓力梯度，然後執行時間積分。這個過程意味著只需要一定的壓縮量就可以觸發核心半徑內的壓力，從而延遲計算。因此，我們需要使用更小的時間步長(意味著更多的迭代)，這在計算上是昂貴的。或者，我們可以使用不那麼嚴格的EOS，然而，這個解決方案可能會引入類似彈簧的振盪。微調引數如聲速或粘度可以幫助避免此類問題。然而，這並不是一個基本的解決方案，對使用者來說也是不切實際的。Solenthaler 和Pajarola通過在SPH模擬中引入預測-校正器概念來解決這個問題。這種又稱為預測校正不可壓縮SPH（PCISPH），它是一種誤差測量演算法，假定測量值和期望密度的差值是誤差。本篇文章總結我在《Fluid Engine Development》學到關於PCISPH的知識。
演算法原理
1.計算合力，預測速度位置，碰撞檢測或碰撞處理
2.計算密度後測量密度誤差，利用密度誤差更新壓強
3.計算新壓強梯度力
4.多次重複以上過程，得到使密度誤差最小的修正力（壓強梯度力），然後使用累積力進行下一步
由於它是通過累積校正力使最終狀態處於不可壓縮狀態，而不是通過多個SPH步驟保證不可壓縮，所以它不需要在意時間步長的問題
演算法程式碼實現結構如下，
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 }
我們可以看到預測-校正本身是一個迴圈體，迴圈在密度誤差小於指定閾值迴圈執行次數小於指定迭代次數時終止
而迴圈體內所做的就是不斷進行修正壓力，直到密度誤差足夠小
需要注意的是在執行這套演算法前，粒子所受其他力包括粘度，重力已經全部計算完成，已全部累加到forces陣列中
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;
}
});