Add collisions and validation for weakly compressible MPM fluids

Includes/Core/Solver/Particle/MPM/MPMFluidSolver-Impl.hpp

@@ -12,6 +12,7 @@
 #define CUBBYFLOW_MPM_FLUID_SOLVER_IMPL_HPP
 
 #include <algorithm>
+#include <array>
 #include <cmath>
 #include <limits>
 #include <stdexcept>
@@ -76,6 +77,18 @@
     m_timeStepLimitScale = newScale;
 }
 
+template <size_t N>
+int MPMFluidSolver<N>::GetClosedDomainBoundaryFlag() const
+{
+    return m_closedDomainBoundaryFlag;
+}
+
+template <size_t N>
+void MPMFluidSolver<N>::SetClosedDomainBoundaryFlag(int flag)
+{
+    m_closedDomainBoundaryFlag = flag;
+}
+
 template <size_t N>
 MPMFluidSolver<N>::Builder MPMFluidSolver<N>::GetBuilder()
 {
@@ -154,11 +167,21 @@
 void MPMFluidSolver<N>::OnBeginAdvanceTimeStep(double timeStepInSeconds)
 {
     m_mpmSystemData->TransferFromParticlesToGrid();
+
     InitializeReferenceVolumes();
     UpdateGridVelocities(timeStepInSeconds);
+    ConstrainGridVelocities();
+
     m_mpmSystemData->TransferFromGridToParticles(timeStepInSeconds);
 }
 
+template <size_t N>
+void MPMFluidSolver<N>::OnEndAdvanceTimeStep(double timeStepInSeconds)
+{
+    Base::OnEndAdvanceTimeStep(timeStepInSeconds);
+    ConstrainParticlesToDomain();
+}
+
 template <size_t N>
 MPMFluidSolver<N>::SizeType MPMFluidSolver<N>::ClampIndex(
     const Vector<ssize_t, N>& index, const SizeType& dataSize)
@@ -252,6 +275,86 @@
     }
 }
 
+template <size_t N>
+void MPMFluidSolver<N>::ConstrainGridVelocities()
+{
+    static constexpr std::array lowerFlags{ DIRECTION_LEFT, DIRECTION_DOWN,
+                                            DIRECTION_BACK };
+    static constexpr std::array upperFlags{ DIRECTION_RIGHT, DIRECTION_UP,
+                                            DIRECTION_FRONT };
+
+    const auto& gridMass = m_mpmSystemData->GridMass();
+    auto& gridVelocities = m_mpmSystemData->GridVelocities();
+    const auto dataSize = gridVelocities.DataSize();
+
+    gridVelocities.ParallelForEachDataPointIndex(
+        [this, &gridMass, &gridVelocities, dataSize](const SizeType& index) {
+            if (gridMass(index) <= 0.0)
+            {
+                return;
+            }
+
+            VectorType velocity = gridVelocities(index);
+
+            if (const auto& collider = this->GetCollider(); collider != nullptr)
+            {
+                VectorType position = gridVelocities.DataPosition()(index);
+                collider->ResolveCollision(0.0, 0.0, &position, &velocity);
+            }
+
+            for (size_t axis = 0; axis < N; ++axis)
+            {
+                const bool exceedsLower =
+                    (m_closedDomainBoundaryFlag & lowerFlags[axis]) != 0 &&
+                    index[axis] == 0 && velocity[axis] < 0.0;
+                const bool exceedsUpper =
+                    (m_closedDomainBoundaryFlag & upperFlags[axis]) != 0 &&
+                    index[axis] == dataSize[axis] - 1 && velocity[axis] > 0.0;
+
+                if (exceedsLower || exceedsUpper)
+                {
+                    velocity[axis] = 0.0;
+                }
+            }
+
+            gridVelocities(index) = velocity;
+        });
+}
+
+template <size_t N>
+void MPMFluidSolver<N>::ConstrainParticlesToDomain()
+{
+    static constexpr std::array lowerFlags{ DIRECTION_LEFT, DIRECTION_DOWN,
+                                            DIRECTION_BACK };
+    static constexpr std::array upperFlags{ DIRECTION_RIGHT, DIRECTION_UP,
+                                            DIRECTION_FRONT };
+
+    const auto domain = m_mpmSystemData->GridMass().GetBoundingBox();
+    auto positions = m_mpmSystemData->Positions();
+    auto velocities = m_mpmSystemData->Velocities();
+
+    ParallelFor(
+        ZERO_SIZE, positions.Length(),
+        [&domain, &positions, &velocities, this](size_t i) {
+            for (size_t axis = 0; axis < N; ++axis)
+            {
+                if ((m_closedDomainBoundaryFlag & lowerFlags[axis]) != 0 &&
+                    positions[i][axis] <= domain.lowerCorner[axis])
+                {
+                    positions[i][axis] = domain.lowerCorner[axis];
+                    velocities[i][axis] = std::max(velocities[i][axis], 0.0);
+                }
+
+                if ((m_closedDomainBoundaryFlag & upperFlags[axis]) != 0 &&
+                    positions[i][axis] >= domain.upperCorner[axis])
+                {
+                    positions[i][axis] = domain.upperCorner[axis];
+                    velocities[i][axis] = std::min(velocities[i][axis], 0.0);
+                }
+            }
+        });
+}
+
 template <size_t N>
 MPMFluidSolver<N>::Builder& MPMFluidSolver<N>::Builder::WithResolution(
     const SizeType& resolution)

Includes/Core/Solver/Particle/MPM/MPMFluidSolver.hpp

@@ -28,6 +28,13 @@ namespace CubbyFlow
 //! Uses the Tait equation of state with existing MPM particle-grid transfers
 //! and particle-solver collision lifecycle.
 //!
+//! \note Choose the artificial speed of sound high enough to keep the maximum
+//! relative density error near or below 3%; ten times the expected peak flow
+//! speed is a common starting point. Adaptive stepping limits each step to
+//! `timeStepLimitScale * minGridSpacing / (speedOfSound + maxParticleSpeed)`.
+//! This explicit limit assumes weak compression; reduce the scale or increase
+//! the sound speed when density error or pressure oscillation is excessive.
+//!
 template <size_t N>
 class MPMFluidSolver : public std::conditional_t<N == 2, ParticleSystemSolver2,
                                                  ParticleSystemSolver3>
@@ -67,6 +74,12 @@ class MPMFluidSolver : public std::conditional_t<N == 2, ParticleSystemSolver2,
     //! Sets the adaptive time-step scale in `(0, 1]`.
     void SetTimeStepLimitScale(double newScale);
 
+    //! Returns the closed domain boundary flag.
+    [[nodiscard]] int GetClosedDomainBoundaryFlag() const;
+
+    //! Sets the closed domain boundary flag.
+    void SetClosedDomainBoundaryFlag(int flag);
+
     //! Returns a builder for MPMFluidSolver.
     [[nodiscard]] static Builder GetBuilder();
 
@@ -84,6 +97,9 @@ class MPMFluidSolver : public std::conditional_t<N == 2, ParticleSystemSolver2,
     //! Advances particle-grid fluid state before base particle integration.
     void OnBeginAdvanceTimeStep(double timeStepInSeconds) override;
 
+    //! Projects particles back into selected closed domain boundaries.
+    void OnEndAdvanceTimeStep(double timeStepInSeconds) override;
+
  private:
     [[nodiscard]] static SizeType ClampIndex(const Vector<ssize_t, N>& index,
                                              const SizeType& dataSize);
@@ -92,9 +108,14 @@ class MPMFluidSolver : public std::conditional_t<N == 2, ParticleSystemSolver2,
 
     void UpdateGridVelocities(double timeStepInSeconds);
 
+    void ConstrainGridVelocities();
+
+    void ConstrainParticlesToDomain();
+
     std::shared_ptr<MPMFluidSystemData<N>> m_mpmSystemData;
     MPMFluidConstitutiveModel<N> m_constitutiveModel;
     double m_timeStepLimitScale = 0.9;
+    int m_closedDomainBoundaryFlag = DIRECTION_ALL;
 };
 
 //! Front-end to create MPMFluidSolver objects step by step.

Sources/API/Python/Solver/Particle/MPM/MPMFluidSolver.cpp

@@ -38,6 +38,9 @@ void AddMPMFluidSolver(pybind11::module& m, const char* name)
         .def_property_readonly("mpmSystemData", &Solver::GetMPMSystemData)
         .def_property("timeStepLimitScale", &Solver::GetTimeStepLimitScale,
                       &Solver::SetTimeStepLimitScale)
+        .def_property("closedDomainBoundaryFlag",
+                      &Solver::GetClosedDomainBoundaryFlag,
+                      &Solver::SetClosedDomainBoundaryFlag)
         .def_property_readonly(
             "targetDensity",
             [](const Solver& solver) {

Tests/ManualTests/MPMFluidSolverTests.cpp

@@ -0,0 +1,79 @@
+// This code is based on Jet framework.
+// Copyright (c) 2018 Doyub Kim
+// CubbyFlow is voxel-based fluid simulation engine for computer games.
+// Copyright (c) 2020 CubbyFlow Team
+// Core Part: Chris Ohk, Junwoo Hwang, Jihong Sin, Seungwoo Yoo
+// AI Part: Dongheon Cho, Minseo Kim
+// We are making my contributions/submissions to this project solely in our
+// personal capacity and are not conveying any rights to any intellectual
+// property of any third parties.
+
+#include "gtest/gtest.h"
+
+#include <Core/Solver/Particle/MPM/MPMFluidSolver.hpp>
+
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+
+using namespace CubbyFlow;
+
+TEST(MPMFluidSolverManual, DamBreakReportsDensityErrorAndMomentum)
+{
+    constexpr double targetDensity = 1000.0;
+    constexpr double particleSpacing = 0.05;
+    constexpr double particleMass =
+        targetDensity * particleSpacing * particleSpacing;
+    MPMFluidSolver2 solver{ { 10, 10 },   { 0.1, 0.1 },  {},  0.02,
+                            particleMass, targetDensity, 30.0 };
+    solver.SetDragCoefficient(0.0);
+    solver.SetTimeStepLimitScale(0.5);
+
+    Array1<Vector2D> positions;
+    for (size_t j = 0; j < 6; ++j)
+    {
+        for (size_t i = 0; i < 3; ++i)
+        {
+            positions.Append(
+                { 0.05 + particleSpacing * static_cast<double>(i),
+                  0.05 + particleSpacing * static_cast<double>(j) });
+        }
+    }
+
+    auto data = solver.GetMPMSystemData();
+    data->AddParticles(positions);
+
+    for (Frame frame{ 0, 1.0 / 240.0 }; frame.index < 40; ++frame)
+    {
+        solver.Update(frame);
+    }
+
+    const auto masses = data->ParticleMasses();
+    const auto volumes = data->InitialVolumes();
+    const auto volumeRatios = data->VolumeRatios();
+    const auto velocities = data->Velocities();
+    double maxDensityError = 0.0;
+    Vector2D totalMomentum;
+
+    for (size_t i = 0; i < data->NumberOfParticles(); ++i)
+    {
+        const double density = solver.GetConstitutiveModel().ComputeDensity(
+            masses[i], volumes[i] * volumeRatios[i]);
+        maxDensityError =
+            std::max(maxDensityError, std::abs(density / targetDensity - 1.0));
+        totalMomentum += masses[i] * velocities[i];
+
+        EXPECT_TRUE(std::isfinite(density));
+        EXPECT_TRUE(std::isfinite(velocities[i].x));
+        EXPECT_TRUE(std::isfinite(velocities[i].y));
+    }
+
+    std::cout << "MPM dam-break max density error: " << 100.0 * maxDensityError
+              << "%\n"
+              << "MPM dam-break total momentum: (" << totalMomentum.x << ", "
+              << totalMomentum.y << ")\n";
+
+    EXPECT_TRUE(std::isfinite(maxDensityError));
+    EXPECT_TRUE(std::isfinite(totalMomentum.x));
+    EXPECT_TRUE(std::isfinite(totalMomentum.y));
+}

Tests/PythonTests/test_mpm_fluid_solver.py

@@ -47,9 +47,12 @@ def test_mpm_fluid_solver_api(
     assert solver.eosExponent == pytest.approx(5.0)
     assert solver.negativePressureScale == pytest.approx(0.25)
     assert solver.timeStepLimitScale == pytest.approx(0.9)
+    assert solver.closedDomainBoundaryFlag == pyCubbyFlow.DIRECTION_ALL
 
     solver.timeStepLimitScale = 0.5
     assert solver.timeStepLimitScale == pytest.approx(0.5)
+    solver.closedDomainBoundaryFlag = pyCubbyFlow.DIRECTION_LEFT
+    assert solver.closedDomainBoundaryFlag == pyCubbyFlow.DIRECTION_LEFT
     with pytest.raises(AttributeError):
         solver.targetDensity = 1000.0