《Fluid Engine Development》 學習筆記4-預測校訂不可壓縮SPH-PCISPH

時間  2019-12-11
標籤 fluid engine development 學習 筆記 預測 校訂 不可 壓縮 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;
            }
        });
相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程