電影工業中的流體模擬（七）----PCISPH

C語言 · 發表 2018-11-04 19:58:37

摘要：上一節我們介紹瞭如何用標準形式的SPH方法求解納維斯托克斯方程。該方法可以保證系統的質量守恆，能夠捕捉浪花飛濺和泡沫等細節現象，在遊戲中得到大量使用，比如 Nvidia PhysX 物理引擎使用該方法來實時模擬水和煙霧(注：從PhysX 3.4 版本開始，SPH方法被棄用，水體模...

上一節我們介紹瞭如何用標準形式的SPH方法求解納維斯托克斯方程。該方法可以保證系統的質量守恆，能夠捕捉浪花飛濺和泡沫等細節現象，在遊戲中得到大量使用，比如 ofollow,noindex">Nvidia PhysX 物理引擎使用該方法來實時模擬水和煙霧(注：從PhysX 3.4 版本開始，SPH方法被棄用，水體模擬推薦使用更穩定也更高效的NVIDIA FleX，這個我們在後續文章中再給大家介紹）。

我司大神Matthias Müller [1]首次將SPH方法引入圖形學領域，以互動式幀率模擬了可壓縮流體，隨著硬體特別是GPU計算能力的快速提升，遊戲中的流體模擬效能得到了大幅提升[2]。但是標準SPH方法存在一個很大的問題，即由於使用氣態方程計算壓強，而導致流體可壓縮。而真實感流體模擬要求密度誤差小以避免流體出現明顯的體積變化，從而導致諸如水體的自由表面出現不自然振盪這樣的視覺失真。針對這一問題，已有若干解決方案，本文主要介紹預測-校正不可壓縮SPH（Predictive-Corrective Incompressible SPH, PCISPH）方法[3]。

至於為什麼挑 PCISPH 來介紹，額...因為...在這個實驗室待過，有感情^-^ 扯遠了... 其實是因為我自己有開原始碼框架可以結合起來介紹，這樣更容易理解。連結： CUDA" target="_blank" rel="nofollow,noindex">UnifiedParticleFrameworkCUDA（A unified particle framework similar to NVIDIA FleX.）

一、 演算法思路

PCISPH演算法如下：（不知道怎麼在知乎裡匯入Latex演算法，那就乾脆截圖吧）

演算法思路：

首先，遍歷所有粒子並查詢其光滑核函式支撐半徑範圍內的鄰居粒子資訊。其次，計算除壓力之外所受到的其他力包括粘性力、重力 Line"/> ，表面張力和其他外力。並且將每個粒子的壓強及壓力初始化為。然後，進入預測-校正迴圈。在該迴圈一開始，根據之前計算出的粒子所受合力預測其速度和位置，並根據預測出來的流體粒子位置資訊預測其密度和密度誤差，並根據第二部分將要介紹的公式更新粒子壓強值。在預測-校正迴圈最後計算壓力，並判斷迴圈繼續的條件是否滿足，即密度誤差大於指定閾值或者迴圈執行次數小於指定最小迭代次數。如果該條件滿足，則表示流體還沒有達到不可壓縮狀態，應該繼續迭代以更改流體壓強。反之當該條件不滿足時表示流體已具備不可壓縮性，可以退出預測-校正迴圈。在單個時間步的最後，即可根據粒子所受合外力計算下一個時間步的速度和位置資訊。

二、理論

PCISPH最重要的理論創新是壓強更新公式：

下面進行推導。如下圖所示，SPH方法是通過對粒子的物理量進行插值計算得到連續流體區域內任一點處的對應值。

流體粒子的密度的計算公式為

這裡，是粒子所處光滑半徑範圍內粒子的質量，，為光滑核函式。

我們只考慮粒子大小一樣，即粒子質量都為的情況，則流體密度計算公式可以簡寫成（這裡光滑半徑為常量，故可以省略）：

設當前時間步為，則下一時間步時流體粒子的密度為

其中，，

由一階泰勒展開，可以得到：

代入上面的密度公式：

而我們知道

故密度增量公式為：

這裡，

文獻[3]中採用的流體粒子壓力計算公式為：

這裡可以做個假設：相鄰流體粒子具有相同的壓強，即，而我們的目標是希望流體不可壓縮，即流體粒子的密度都為靜態密度，則壓力計算公式可以改寫為：

PCISPH演算法在每次迭代校正時只考慮流體粒子的壓力，通過迭代計算出的壓力讓粒子之間不至於靠的過近（粒子之間距離太近可以理解為流體可壓縮）。在只考慮壓力的情況下，根據時間積分可以計算出粒子在該壓力作用下的位移：

前面我們計算粒子的壓力是通過對所有鄰居粒子求和得到，這裡如果我們只考慮來自於粒子的壓力，則其值為：

根據牛頓第三定律，粒子受到的來自於粒子的壓力為：

同樣，由於粒子的壓力導致粒子產生的位移為：

將上面的公式代入密度增量公式：

由上式可得：

其中，

該公式可以這麼理解：在迭代過程中，我們需要不斷校正流體粒子的密度使其接近靜態密度從而滿足不可壓縮性條件，當密度變化為時需要的壓強值為。每次迭代需要預測流體粒子的密度值，這個值與靜態密度之間的誤差為。當前粒子的密度為，而我們的目標是讓粒子密度為，這需要的密度改變數為。因此，我們需要施加的壓強值為：

這個壓強計算公式在鄰居粒子數目不足的時候會導致計算錯誤，解決辦法是進行一次預計算，在流體粒子周圍充滿鄰居粒子的情況下計算係數，使得。

則每次迭代中壓強值按如下公式進行更新：

三、實現

UnifiedParticleFrameworkCUDA

程式碼分CPU & GPU版本，也分單相流體和流固耦合模擬，為確保簡潔，本文結合CPU單相不可壓縮流體程式碼介紹，對GPU實現感興趣的同學可以參考CudaParticlePhysics()。PCISPH核心演算法實現主要在下面這段程式碼中：

void UnifiedPhysics::CpuParticlePhysicsPureFluidsPCISPH()
{
NeighborSearchOMP();

CalculateExternalForcesPCISPHOMP();

PredictionCorrectionStep();

UpdateLiquidParticlePosVelCPUOMP();

BoundaryHandlingPureFluidOMP();

// resort particles array using radix sort using their z-index values
SortParticles();
}

首先我們要進行鄰居粒子查詢：

我們使用均勻網格對模擬區域進行劃分,如下圖所示：

每個流體粒子根據其空間位置被劃分為其中一個網格單元。因此，為了確定流體粒子的鄰居節點資訊，只需將搜尋範圍限制在該粒子所在網格單元及其相鄰26個網格單元內的全部流體粒子。

鄰居粒子查詢程式碼為：

void UnifiedPhysics::NeighborSearchOMP()
{
int i;
int chunk = 100;

#pragma omp parallel default(shared) private(i)
{
#pragma omp for schedule(dynamic, chunk) nowait
for (i = 0; i < particles_.size(); i++)
GetNeighbors(i);
}
}

void UnifiedPhysics::GetNeighbors(const int i)
{
UnifiedParticle &p = particles_[i];

// set query box
vmml::Vector3f min_query_box(fc_->scales * (p.position_ - fc_->globalSupportRadius)), max_query_box(fc_->scales * (p.position_ + fc_->globalSupportRadius));
// limit to bounding box of spatial domain
for (int i = 0; i < 3; i++) {
min_query_box[i] = MAX(min_query_box[i], 0);
max_query_box[i] = MIN(max_query_box[i], GRID_RESOLUTION);
}

// get index ranges covering the given query box
ZRanges ranges;
zIndex.BoxQuery(min_query_box, max_query_box, ranges);

// sort ranges
QuickSort(ranges.begin(), ranges.end());

// reset numbers of neighbors
std::vector<int> &neighbor_indices = neighbor_indices_[i];
neighbor_indices.clear();

// identify particles within ranges
std::vector<UnifiedParticle>::iterator q;
int counter = 0;
for (int k = 0; k < ranges.size(); k++) {
// the elements in the range shall already be sorted for lower_bound algorithm
for (q = std::lower_bound(particles_.begin(), particles_.end(), ranges[k].first); q < particles_.end() && q->index_ <= ranges[k].second; q++) {
int j = q - particles_.begin();
if (j == i) continue;

float dist = p.position_.distance(q->position_);
if (dist < fc_->globalSupportRadius) {
counter++;
neighbor_indices.push_back(j);
}
}
}
}

然後，我們要計算除壓力之外所受到的其他力，對應程式碼為：

void UnifiedPhysics::CalculateExternalForcesPCISPHOMP()
{
int i;
int chunk = 100;

#pragma omp parallel default(shared) private(i)
#pragma omp for schedule(dynamic, chunk) nowait
for (i = 0; i < particles_.size(); i++)
CalculateExternalForcesPCISPH(i);
}

void UnifiedPhysics::CalculateExternalForcesPCISPH(const int i)
{
// this method is called when using PCISPH
UnifiedParticle &p = particles_[i];
std::vector<int> &indices = neighbor_indices_[i];

p.force_ = (0.0, 0.0, 0.0);
const float restVolume = fc_->initialMass/ fc_->fluidRestDensity;

for (int j = 0; j < indices.size(); j++) {
float dist = p.position_.distance(particles_[indices[j]].position_);
if (dist < fc_->globalSupportRadius) {
float kernelVisc = my_kernel_->kernelViscosityLaplacianLut(dist); // symmetric

// compute artificial viscosity according to MCG03
// negative symmetry
vmml::Vector3f v_ij = p.velocity_ - particles_[indices[j]].velocity_;
p.force_ -= v_ij * restVolume * restVolume * fc_->fluidViscConst * kernelVisc;
}
}

// add gravity
p.force_.y -= fc_->initialMass * fc_->gravityConst; 

// add boundary forces 
if (fc_->addBoundaryForce)
p.force_ += BoxBoundaryForce(i);

// init some quantities which are going to be used in the prediction step
p.correction_pressure_ = 0.0f;
p.correction_pressure_force_ = (0.0, 0.0, 0.0);
}

接下來就是進入預測-校正迭代：

void UnifiedPhysics::PredictionCorrectionStep()
{
int i;
int chunk = 100;
int particlesSize = particles_.size();

density_error_too_large_ = true; // loop has to be executed at least once

int iteration = 0;
while( (iteration < fc_->minLoops) || ((density_error_too_large_) && (iteration < fc_->maxLoops)) )
{
#pragma omp parallel default(shared) private(i)
#pragma omp for schedule(dynamic, chunk) nowait
for(int i = 0; i < particles_.size(); i++)
PredictPositionAndVelocity(i);

max_predicted_density_ = 0.0; // loop termination criterion

#pragma omp parallel default(shared) private(i)
#pragma omp for schedule(dynamic, chunk) nowait
for(int i = 0; i < particles_.size(); i++)
ComputePredictedDensityAndPressure(i);

// check loop termination criterion
float densityErrorInPercent = MAX(0.1f * max_predicted_density_ - 100.0f, 0.0f); // 100/1000 * max_predicted_density_ - 100; 

if(fc_->printDebuggingInfo==1)
std::cout << "ERROR: " << densityErrorInPercent << "%" << std::endl;

// set flag to terminate loop if density error is smaller than requested
if(densityErrorInPercent < fc_->maxDensityErrorAllowed) 
density_error_too_large_ = false; // stop loop

#pragma omp parallel default(shared) private(i)
#pragma omp for schedule(dynamic, chunk) nowait
for(int i = 0; i < particles_.size(); i++)
ComputeCorrectivePressureForce(i);

iteration++;
}
}

迭代一開始，根據計算出的粒子所受合力預測其速度和位置：

void UnifiedPhysics::PredictPositionAndVelocity(const int i)
{
// v' = v + delta_t * a
// a = F / m
// v' = v + delta_t * F * V / m
// v' = v + delta_t * F * 1/density

// compute predicted position and velocity
UnifiedParticle &p = particles_[i];
if (p.type_ == LIQUID_PARTICLE)
{
p.predicted_velocity_ = p.velocity_ + (p.force_ + p.correction_pressure_force_) * fc_->deltaT / fc_->initialMass;
p.predicted_position_ = p.position_ + p.predicted_velocity_ * fc_->deltaT; // always use position at time t

CollisionHandling(p.predicted_position_, p.predicted_velocity_); 

// init some quantities which are going to be used in the next step
p.predicted_density_ = 0.0;

// check if particle has valid position (==not nan) to detect instabilities
if(!(p.predicted_position_[0] == p.predicted_position_[0]) ||
!(p.predicted_position_[1] == p.predicted_position_[1]) ||
!(p.predicted_position_[2] == p.predicted_position_[2]))
{
std::cout << "Particle has invalid predictedPosition!!" << std::endl;
abort();
}
}
}

然後根據預測出來的流體粒子位置資訊預測其密度，密度誤差並更新壓強：

void UnifiedPhysics::ComputePredictedDensityAndPressure(const int i)
{
// sph density
UnifiedParticle& p = particles_[i];
if (p.type_ == LIQUID_PARTICLE)
{
std::vector<int> &indices = neighbor_indices_[i];

p.predicted_density_ = kernel_self_ * fc_->initialMass;

// neighborhoods of current positions are used and not of predicted positions
// this might lead to problems for large velocities
for (int j = 0; j < indices.size(); j++) {
// we need to recompute the distance and kernel values using the predicted positions
float dist = p.predicted_position_.distance(particles_[indices[j]].predicted_position_);
if (dist < fc_->globalSupportRadius) {
float kernelValue = my_kernel_->kernelM4Lut(dist); // use LUT to make loop execution faster

// sph dens, symmetric version
p.predicted_density_ += kernelValue * fc_->initialMass;
}
}

// compute density error, correct only compression errors
p.density_error_ = MAX(p.predicted_density_ - fc_->fluidRestDensity, 0.0f);

// update pressure
// densityErrorFactor is precomputed and used as a constant (approximation)
// do not allow negative pressure corrections
p.correction_pressure_ += MAX(p.density_error_ * fc_->densityErrorFactor, 0.0f);

// max_predicted_density_ is needed for the loop termination criterion
// find maximal predicted density of all particles
max_predicted_density_ = MAX(max_predicted_density_, p.predicted_density_);
}
}

最後計算壓力：

void UnifiedPhysics::ComputeCorrectivePressureForce(const int i)
{
// init the pressure forces (after checking loop criterion / return above)
UnifiedParticle& p = particles_[i];
if (p.type_ == LIQUID_PARTICLE)
{
std::vector<int> &indices = neighbor_indices_[i];

p.correction_pressure_force_ = (0.0, 0.0, 0.0);

// compute pressure forces with correction pressure
float densSq = fc_->fluidRestDensity*fc_->fluidRestDensity;
float p_i = p.correction_pressure_;

for (int j = 0; j < indices.size(); j++) {
UnifiedParticle& neigh = particles_[indices[j]];
float dist = p.position_.distance(neigh.position_);
if (dist < fc_->globalSupportRadius) {
float p_j = neigh.correction_pressure_;

float kernelGradientValue = my_kernel_->kernelPressureGradLut(dist);
vmml::Vector3f kernelGradient = (p.position_ - neigh.position_) * kernelGradientValue;

// sum up pressure force according to Monaghan
float grad = p_i / densSq + p_j / (fc_->fluidRestDensity*fc_->fluidRestDensity);
p.correction_pressure_force_ -= kernelGradient * grad * fc_->initialMass * fc_->initialMass;
}
}
}
}

預測-校正迭代退出後，即可進行數值積分步驟更新粒子位置和速度資訊：

void UnifiedPhysics::UpdateLiquidParticlePosVelCPUOMP()
{
int i;
int chunk = 100;

#pragma omp parallel default(shared) private(i)
{
#pragma omp for schedule(dynamic, chunk) nowait
for (i = 0; i < particles_.size(); i++)
UpdateLiquidParticlePosVelSphCPU(i);
}
}

void UnifiedPhysics::UpdateLiquidParticlePosVelSphCPU(const int i)
{
UnifiedParticle &p = particles_[i];
if(p.type_ == LIQUID_PARTICLE)
{
// in case of PCISPH add the pressureForce to total force
if(fc_->physicsType == 'i') 
p.force_ += p.correction_pressure_force_;

static const float invMass = 1.0 / fc_->initialMass;
static const float deltaT = fc_->deltaT;
/************************************************************************/
/*Symplectic Euler Integration*/
/************************************************************************/
//*
p.velocity_ += p.force_ * invMass * deltaT;
p.position_ += p.velocity_ * deltaT;
//*/

/************************************************************************/
/*Leapfrog integration*/
/************************************************************************/
/*
vmml::Vector3f vel_half_previous = p.velocity_leapfrog;// read old velocity_leapfrog v(t-1/2)
vmml::Vector3f vel_half_next = vel_half_previous + p.force * invMass * deltaT;// calculate new velocity_leapfrog v(t+1/2) = v(t-1/2) + a(t) * dt
p.velocity_leapfrog = vel_half_next;// update new velocity_leapfrog 
p.velocity = (vel_half_previous + vel_half_next) * 0.5f;// update new velocity  v(t+1) = [v(t-1/2) + v(t+1/2)] * 0.5
p.position += vel_half_next * deltaT;// update new position p(t+1) = p(t) + v(t+1/2) * dt
//*/
}
}

當然，還需要進行與周圍物體的碰撞檢測與響應，下面這段程式碼只考慮簡單情況：

void UnifiedPhysics::BoundaryHandlingPureFluidOMP()
{
const vmml::Vector3f minBox = fc_->collisionBox.getMin();
const vmml::Vector3f maxBox = fc_->collisionBox.getMax();
static const float damping = -0.5; // 0.1
static const int numParticles = particles_.size(); 

int i;
int chunk = 100;
#pragma omp parallel default(shared) private(i)
{
#pragma omp for schedule(dynamic, chunk) nowait
for (i = 0; i < numParticles; i++)
{
UnifiedParticle* p = &particles_[i];

if(p->type_ == LIQUID_PARTICLE)
{
// non rigid & non boundary particles
BoundaryHandlingBoxPerNonRigidParticle(minBox, maxBox, damping, p->position_, p->velocity_);
}

// get new spatial index
p->index_ = zIndex.CalcIndex(fc_->scales * p->position_);
}
}
}

這樣，一個時間步的PCISPH模擬就算完成。

更多實現細節請參看原始碼。

總結：

PCISPH方法基於拆分思想，將作用於流體粒子上的力分為兩類：壓力和其他力。首先，使用其他力驅動流體粒子，預測粒子的中間位置和速度等物理資訊。然後進入預測校正迴圈，即利用壓力對預測的粒子物理屬性值進行迭代預測-校正，當密度誤差達到指定閾值即滿足指定不可壓縮性時退出預測-校正步驟。最後，使用合外力對所有粒子進行時間積分求解出下一個時間步長的位置和速度資訊。

有興趣的同學可以閱讀原始碼，框架裡面提供了標準SPH/WCSPH/PCISPH/DEM(Discrete Element Method)等演算法的CPU&GPU實現。

PS：寫這篇文章的時候，翻了下以前的程式碼，四五年前的了，記憶有些模糊，看到了下面這個：

enum ParticleObjectType
{
RIGID_PARTICLE= 0,// rigid body
SOFT_PARTICLE= 1,// soft/deformable body
LIQUID_PARTICLE= 2,// liquids
GRANULAR_PARTICLE= 3,// granular materials
CLOTH_PARTICLE= 4,// cloth, shell, thin object, etc...
SMOKE_PARTICLE= 5,// smoke
FROZEN_PARTICLE= 6// frozen boundaries (static or kinematic boundaries)
};

想起來當時做這個project時候的確是想用particle一統江湖，利用CUDA並行優化加速，然後模擬各種狂拽酷炫吊炸天的特效，這跟FleX構思如此一致，當然程式碼質量確實差了好大一截。人都是喜新厭舊的，有了新歡（FleX）就忘了舊愛（UnifiedParticleFrameworkCUDA）。有空再給大家介紹下FleX。

參考文獻：

[1] Müller, Matthias, David Charypar, and Markus Gross. "Particle-based fluid simulation for interactive applications." Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation . Eurographics Association, 2003.

[2] Schirm, Simon, Mark Harris. "Fast GPU Fluid Simulation in PhysX." Chapter 7.3 of Game Programming Gems 8, Adam Lake.

[3] Solenthaler, Barbara, and Renato Pajarola. "Predictive-corrective incompressible SPH." ACM transactions on graphics (TOG) . Vol. 28. No. 3. ACM, 2009.

電影工業中的流體模擬（七）----PCISPH

您可能也會喜歡…