From dbd6d1616f58cf5e548d78aaeecc87b9dcbe8048 Mon Sep 17 00:00:00 2001 From: Chris Ohk Date: Tue, 11 Aug 2026 00:04:32 +0900 Subject: [PATCH 01/13] feat: add snow stress differential --- .../MPM/SnowConstitutiveModel-Impl.hpp | 73 ++++++++++++++++++- .../Particle/MPM/SnowConstitutiveModel.hpp | 9 +++ .../UnitTests/SnowConstitutiveModelTests.cpp | 56 ++++++++++++++ 3 files changed, 136 insertions(+), 2 deletions(-) diff --git a/Includes/Core/Particle/MPM/SnowConstitutiveModel-Impl.hpp b/Includes/Core/Particle/MPM/SnowConstitutiveModel-Impl.hpp index 416a1fa29d..e0773b0ce4 100644 --- a/Includes/Core/Particle/MPM/SnowConstitutiveModel-Impl.hpp +++ b/Includes/Core/Particle/MPM/SnowConstitutiveModel-Impl.hpp @@ -60,7 +60,7 @@ SnowConstitutiveModel::SnowConstitutiveModel(double youngsModulus, } template -typename SnowConstitutiveModel::State SnowConstitutiveModel::Update( +SnowConstitutiveModel::State SnowConstitutiveModel::Update( const MatrixType& deformationGradientIncrement, const State& state) const { ValidateDeformation(deformationGradientIncrement); @@ -99,7 +99,7 @@ typename SnowConstitutiveModel::State SnowConstitutiveModel::Update( } template -typename SnowConstitutiveModel::MatrixType +SnowConstitutiveModel::MatrixType SnowConstitutiveModel::ComputeKirchhoffStress(const State& state) const { ValidateDeformation(state.elastic); @@ -130,6 +130,75 @@ SnowConstitutiveModel::ComputeKirchhoffStress(const State& state) const return stress; } +template +SnowConstitutiveModel::MatrixType +SnowConstitutiveModel::ComputeFirstPiolaStressDifferential( + const State& state, const MatrixType& differential) const +{ + ValidateDeformation(state.elastic); + ValidateDeformation(state.plastic); + + if (!IsFinite(differential)) + { + throw std::invalid_argument{ "Invalid snow deformation differential." }; + } + + MatrixType u; + Vector singularValues; + MatrixType v; + + SVD(state.elastic, u, singularValues, v); + + const MatrixType principalDifferential = u.Transposed() * differential * v; + MatrixType omega; + + for (size_t i = 0; i < N; ++i) + { + for (size_t j = i + 1; j < N; ++j) + { + const double value = + (principalDifferential(i, j) - principalDifferential(j, i)) / + (singularValues[i] + singularValues[j]); + omega(i, j) = value; + omega(j, i) = -value; + } + } + const MatrixType rotationDifferential = u * omega * v.Transposed(); + const MatrixType& f = state.elastic; + const double determinant = f.Determinant(); + const MatrixType inverseTranspose = f.Inverse().Transposed(); + const MatrixType cofactor = determinant * inverseTranspose; + double determinantDifferential = 0.0; + + for (size_t row = 0; row < N; ++row) + { + for (size_t column = 0; column < N; ++column) + { + determinantDifferential += + cofactor(row, column) * differential(row, column); + } + } + + const MatrixType cofactorDifferential = + determinantDifferential * inverseTranspose - + determinant * inverseTranspose * differential.Transposed() * + inverseTranspose; + const double hardening = ComputeHardening(state); + const double mu = m_mu0 * hardening; + const double lambda = m_lambda0 * hardening; + const MatrixType result = + 2.0 * mu * (differential - rotationDifferential) + + lambda * (determinantDifferential * cofactor + + (determinant - 1.0) * cofactorDifferential); + + if (!IsFinite(result)) + { + throw std::invalid_argument{ "Non-finite snow stress differential." }; + } + + return result; +} + template double SnowConstitutiveModel::ComputeWaveSpeed(const State& state, double referenceDensity) const diff --git a/Includes/Core/Particle/MPM/SnowConstitutiveModel.hpp b/Includes/Core/Particle/MPM/SnowConstitutiveModel.hpp index 936f6f042a..4fde042ed5 100644 --- a/Includes/Core/Particle/MPM/SnowConstitutiveModel.hpp +++ b/Includes/Core/Particle/MPM/SnowConstitutiveModel.hpp @@ -90,6 +90,15 @@ class SnowConstitutiveModel final //! \return Fixed-corotated Kirchhoff stress. [[nodiscard]] MatrixType ComputeKirchhoffStress(const State& state) const; + //! + //! \brief Computes the first Piola stress differential. + //! + //! Evaluates `(d^2 Psi / d F_E d F_E) : differential` while holding the + //! plastic deformation fixed. + //! + [[nodiscard]] MatrixType ComputeFirstPiolaStressDifferential( + const State& state, const MatrixType& differential) const; + //! //! \brief Estimates the fastest elastic wave speed for the state. //! diff --git a/Tests/UnitTests/SnowConstitutiveModelTests.cpp b/Tests/UnitTests/SnowConstitutiveModelTests.cpp index f9342c3756..4225d2376d 100644 --- a/Tests/UnitTests/SnowConstitutiveModelTests.cpp +++ b/Tests/UnitTests/SnowConstitutiveModelTests.cpp @@ -34,6 +34,52 @@ MatrixD MakeQuarterTurn() return result; } +template +MatrixD ComputeFirstPiola(const SnowConstitutiveModel& model, + const SnowDeformationState& state) +{ + return model.ComputeKirchhoffStress(state) * + state.elastic.Inverse().Transposed(); +} + +template +void ExpectFirstPiolaDifferentialMatchesFiniteDifference() +{ + const SnowConstitutiveModel model{ 1000.0, 0.2, 0.025, 0.0075, 10.0 }; + SnowDeformationState state; + + state.elastic = MakeQuarterTurn() * MakeStretch(1.005); + state.plastic = MakeStretch(0.9); + + MatrixD differential; + differential(0, 0) = 0.2; + differential(0, 1) = -0.3; + differential(1, 0) = 0.4; + differential(1, 1) = -0.1; + + if constexpr (N == 3) + { + differential(2, 0) = 0.15; + differential(1, 2) = -0.2; + differential(2, 2) = 0.25; + } + + constexpr double eps = 1e-6; + auto plus = state; + auto minus = state; + + plus.elastic += eps * differential; + minus.elastic -= eps * differential; + + const MatrixD expected = + (ComputeFirstPiola(model, plus) - ComputeFirstPiola(model, minus)) / + (2.0 * eps); + const MatrixD actual = + model.ComputeFirstPiolaStressDifferential(state, differential); + + EXPECT_TRUE(actual.IsSimilar(expected, 1e-6)); +} + template void ExpectIdentityStateAndZeroStress() { @@ -257,6 +303,11 @@ void ExpectInvalidStatesRejected() EXPECT_THROW( (void)SnowConstitutiveModel{}.ComputeKirchhoffStress(inverted), std::invalid_argument); + + EXPECT_THROW( + (void)SnowConstitutiveModel{}.ComputeFirstPiolaStressDifferential( + {}, nonFinite), + std::invalid_argument); } template @@ -342,6 +393,11 @@ TEST(SnowConstitutiveModel, ElasticStress) RUN_FOR_2D_AND_3D(ExpectKnownElasticStress); } +TEST(SnowConstitutiveModel, FirstPiolaDifferential) +{ + RUN_FOR_2D_AND_3D(ExpectFirstPiolaDifferentialMatchesFiniteDifference); +} + TEST(SnowConstitutiveModel, Compression) { RUN_FOR_2D_AND_3D(ExpectCompressionClamped); From ef029c9875430d68a183aa8cab93c78d9e00fc1b Mon Sep 17 00:00:00 2001 From: Chris Ohk Date: Tue, 11 Aug 2026 00:04:49 +0900 Subject: [PATCH 02/13] feat: add conjugate residual solver --- Includes/Core/Math/CG-Impl.hpp | 73 +++++++++++++++++++++++++++++++++- Includes/Core/Math/CG.hpp | 14 ++++++- Tests/UnitTests/CGTests.cpp | 58 ++++++++++++++++++++++++++- 3 files changed, 142 insertions(+), 3 deletions(-) diff --git a/Includes/Core/Math/CG-Impl.hpp b/Includes/Core/Math/CG-Impl.hpp index a0065ab782..4ae7402387 100644 --- a/Includes/Core/Math/CG-Impl.hpp +++ b/Includes/Core/Math/CG-Impl.hpp @@ -13,6 +13,8 @@ #include +#include + namespace CubbyFlow { template @@ -32,6 +34,75 @@ void CG(const typename BLASType::MatrixType& A, lastResidualNorm); } +template +void CR(const typename BLASType::MatrixType& A, + const typename BLASType::VectorType& b, + unsigned int maxNumberOfIterations, double tolerance, + typename BLASType::VectorType* x, typename BLASType::VectorType* r, + typename BLASType::VectorType* d, typename BLASType::VectorType* q, + typename BLASType::VectorType* s, unsigned int* lastNumberOfIterations, + double* lastResidualNorm) +{ + BLASType::Residual(A, *x, b, r); + BLASType::Set(*r, d); + BLASType::MVM(A, *r, s); + BLASType::Set(*s, q); + + double rho = BLASType::Dot(*r, *s); + double residualNorm = BLASType::L2Norm(*r); + unsigned int iter = 0; + + while (residualNorm > tolerance && iter < maxNumberOfIterations) + { + const double denominator = BLASType::Dot(*q, *q); + + if (!std::isfinite(rho) || !std::isfinite(denominator) || + denominator <= 0.0 || rho == 0.0) + { + residualNorm = std::numeric_limits::infinity(); + break; + } + + const double alpha = rho / denominator; + + BLASType::AXPlusY(alpha, *d, *x, x); + BLASType::AXPlusY(-alpha, *q, *r, r); + + residualNorm = BLASType::L2Norm(*r); + ++iter; + + if (!std::isfinite(residualNorm)) + { + residualNorm = std::numeric_limits::infinity(); + break; + } + + if (residualNorm <= tolerance) + { + break; + } + + BLASType::MVM(A, *r, s); + + const double rhoNew = BLASType::Dot(*r, *s); + + if (!std::isfinite(rhoNew)) + { + residualNorm = std::numeric_limits::infinity(); + break; + } + + const double beta = rhoNew / rho; + + BLASType::AXPlusY(beta, *d, *r, d); + BLASType::AXPlusY(beta, *q, *s, q); + rho = rhoNew; + } + + *lastNumberOfIterations = iter; + *lastResidualNorm = residualNorm; +} + template void PCG(const typename BLASType::MatrixType& A, const typename BLASType::VectorType& b, @@ -113,4 +184,4 @@ void PCG(const typename BLASType::MatrixType& A, } } // namespace CubbyFlow -#endif \ No newline at end of file +#endif diff --git a/Includes/Core/Math/CG.hpp b/Includes/Core/Math/CG.hpp index 85aa9ac165..2dc295836c 100644 --- a/Includes/Core/Math/CG.hpp +++ b/Includes/Core/Math/CG.hpp @@ -48,6 +48,18 @@ void CG(const typename BLASType::MatrixType& A, typename BLASType::VectorType* s, unsigned int* lastNumberOfIterations, double* lastResidualNorm); +//! +//! \brief Solves a symmetric linear system with conjugate residual. +//! +template +void CR(const typename BLASType::MatrixType& A, + const typename BLASType::VectorType& b, + unsigned int maxNumberOfIterations, double tolerance, + typename BLASType::VectorType* x, typename BLASType::VectorType* r, + typename BLASType::VectorType* d, typename BLASType::VectorType* q, + typename BLASType::VectorType* s, unsigned int* lastNumberOfIterations, + double* lastResidualNorm); + //! //! \brief Solves pre-conditioned conjugate gradient. //! @@ -63,4 +75,4 @@ void PCG(const typename BLASType::MatrixType& A, #include -#endif \ No newline at end of file +#endif diff --git a/Tests/UnitTests/CGTests.cpp b/Tests/UnitTests/CGTests.cpp index 27c52840fe..a199adc891 100644 --- a/Tests/UnitTests/CGTests.cpp +++ b/Tests/UnitTests/CGTests.cpp @@ -92,4 +92,60 @@ TEST(PCG, Solve) EXPECT_LE(lastResidualNorm, std::numeric_limits::epsilon()); EXPECT_LE(lastNumIter, 2u); } -} \ No newline at end of file +} + +TEST(CR, SolveSymmetricSystems) +{ + using BLASType = BLAS; + + for (const auto& matrix : + { Matrix2x2D(4.0, 1.0, 1.0, 3.0), Matrix2x2D(2.0, 0.0, 0.0, -1.0) }) + { + const Vector2D expected(1.0, 1.0); + const Vector2D rhs = matrix * expected; + Vector2D x, r, d, q, s; + unsigned int iterations = 0; + double residual = 0.0; + + CR(matrix, rhs, 10, 1e-14, &x, &r, &d, &q, &s, &iterations, + &residual); + + EXPECT_TRUE(x.IsSimilar(expected, 1e-12)); + EXPECT_LE(residual, 1e-12); + EXPECT_LE(iterations, 2u); + } +} + +TEST(CR, ZeroIterationsReportsInitialResidual) +{ + using BLASType = BLAS; + const Matrix2x2D matrix(4.0, 1.0, 1.0, 3.0); + const Vector2D rhs(1.0, 2.0); + Vector2D x, r, d, q, s; + unsigned int iterations = 1; + double residual = 0.0; + + CR(matrix, rhs, 0, 0.0, &x, &r, &d, &q, &s, &iterations, + &residual); + + EXPECT_EQ(iterations, 0u); + EXPECT_DOUBLE_EQ(residual, std::sqrt(5.0)); + EXPECT_EQ(x, Vector2D{}); +} + +TEST(CR, AlreadyConverged) +{ + using BLASType = BLAS; + const Matrix2x2D matrix(4.0, 1.0, 1.0, 3.0); + Vector2D x(1.0, 1.0); + const Vector2D rhs = matrix * x; + Vector2D r, d, q, s; + unsigned int iterations = 1; + double residual = 1.0; + + CR(matrix, rhs, 10, 1e-14, &x, &r, &d, &q, &s, &iterations, + &residual); + + EXPECT_EQ(iterations, 0u); + EXPECT_LE(residual, 1e-14); +} From c33609b69c547585b30bec3c1892034ab318b1e8 Mon Sep 17 00:00:00 2001 From: Chris Ohk Date: Tue, 11 Aug 2026 00:05:16 +0900 Subject: [PATCH 03/13] feat: add semi-implicit snow MPM integration --- .../Particle/MPM/SnowMPMSolver-Impl.hpp | 436 +++++++++++++++++- .../Solver/Particle/MPM/SnowMPMSolver.hpp | 53 ++- Tests/UnitTests/SnowMPMSolverTests.cpp | 335 +++++++++++++- 3 files changed, 807 insertions(+), 17 deletions(-) diff --git a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp index 05e27ebe4b..03fc6a941d 100644 --- a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp +++ b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp @@ -21,6 +21,77 @@ namespace CubbyFlow { +template +struct SnowMPMSolver::LinearSystem +{ + const SnowMPMSolver* solver; + const Array1* activeNodes; + const Array1* nodeToActive; + const Array1* constrained; + double dtSquared; + + void Multiply(const VectorND& input, VectorND* output) const; +}; + +template +struct SnowMPMSolver::LinearSystemBLAS : BLAS +{ + using System = LinearSystem; + using Base = BLAS; + + static void MVM(const System& system, const VectorND& vector, + VectorND* result) + { + system.Multiply(vector, result); + } + + static void Residual(const System& system, const VectorND& x, + const VectorND& b, VectorND* result) + { + system.Multiply(x, result); + Base::AXPlusY(-1.0, *result, b, result); + } +}; + +template +void SnowMPMSolver::LinearSystem::Multiply(const VectorND& input, + VectorND* output) const +{ + VectorND projectedInput = input; + + for (size_t i = 0; i < projectedInput.GetRows(); ++i) + { + if ((*constrained)[i] != 0) + { + projectedInput[i] = 0.0; + } + } + + VectorND hessian(input.GetRows(), 0.0); + + solver->ApplyElasticHessian(*activeNodes, *nodeToActive, projectedInput, + &hessian); + output->Resize(input.GetRows(), 0.0); + + const auto& gridMass = solver->m_mpmSystemData->GridMass(); + + for (size_t slot = 0; slot < activeNodes->Length(); ++slot) + { + const double mass = gridMass((*activeNodes)[slot]); + + for (size_t axis = 0; axis < N; ++axis) + { + const size_t row = slot * N + axis; + + if ((*constrained)[row] == 0) + { + (*output)[row] = + mass * projectedInput[row] + dtSquared * hessian[row]; + } + } + } +} + template SnowMPMSolver::SnowMPMSolver(const SizeType& resolution, const VectorType& gridSpacing, @@ -67,6 +138,62 @@ void SnowMPMSolver::SetTimeStepLimitScale(double newScale) m_timeStepLimitScale = newScale; } +template +bool SnowMPMSolver::GetIsUsingSemiImplicit() const +{ + return m_isUsingSemiImplicit; +} + +template +void SnowMPMSolver::SetIsUsingSemiImplicit(bool isUsing) +{ + m_isUsingSemiImplicit = isUsing; +} + +template +unsigned int SnowMPMSolver::GetMaxNumberOfIterations() const +{ + return m_maxNumberOfIterations; +} + +template +void SnowMPMSolver::SetMaxNumberOfIterations( + unsigned int maxNumberOfIterations) +{ + m_maxNumberOfIterations = maxNumberOfIterations; +} + +template +double SnowMPMSolver::GetTolerance() const +{ + return m_tolerance; +} + +template +void SnowMPMSolver::SetTolerance(double tolerance) +{ + if (!std::isfinite(tolerance) || tolerance <= 0.0) + { + throw std::invalid_argument{ + "Semi-implicit tolerance must be positive and finite." + }; + } + + m_tolerance = tolerance; +} + +template +unsigned int SnowMPMSolver::GetLastNumberOfIterations() const +{ + return m_lastNumberOfIterations; +} + +template +double SnowMPMSolver::GetLastResidual() const +{ + return m_lastResidual; +} + template int SnowMPMSolver::GetClosedDomainBoundaryFlag() const { @@ -80,7 +207,7 @@ void SnowMPMSolver::SetClosedDomainBoundaryFlag(int flag) } template -typename SnowMPMSolver::Builder SnowMPMSolver::GetBuilder() +SnowMPMSolver::Builder SnowMPMSolver::GetBuilder() { return Builder{}; } @@ -141,7 +268,7 @@ unsigned int SnowMPMSolver::GetNumberOfSubTimeSteps( { desiredTimeStep = std::min(desiredTimeStep, minSpacing / maxVelocity); } - if (maxWaveSpeed > 0.0) + if (!m_isUsingSemiImplicit && maxWaveSpeed > 0.0) { desiredTimeStep = std::min(desiredTimeStep, minSpacing / maxWaveSpeed); } @@ -176,7 +303,28 @@ void SnowMPMSolver::OnBeginAdvanceTimeStep(double timeStepInSeconds) m_mpmSystemData->TransferFromParticlesToGrid(); InitializeReferenceVolumes(); UpdateGridVelocities(timeStepInSeconds); - ConstrainGridVelocities(); + + Array1 activeNodes; + Array1 nodeToActive; + + BuildActiveNodes(&activeNodes, &nodeToActive); + + Array1 constrained(activeNodes.Length() * N, uint8_t{ 0 }); + + ConstrainGridVelocities(activeNodes, nodeToActive, &constrained); + + if (m_isUsingSemiImplicit) + { + SolveGridVelocities(timeStepInSeconds, activeNodes, nodeToActive, + constrained); + ConstrainGridVelocities(activeNodes, nodeToActive, nullptr); + } + else + { + m_lastNumberOfIterations = 0; + m_lastResidual = 0.0; + } + UpdateDeformation(timeStepInSeconds); m_mpmSystemData->TransferFromGridToParticles(); } @@ -189,7 +337,7 @@ void SnowMPMSolver::OnEndAdvanceTimeStep(double timeStepInSeconds) } template -typename SnowMPMSolver::SizeType SnowMPMSolver::ClampIndex( +SnowMPMSolver::SizeType SnowMPMSolver::ClampIndex( const Vector& index, const SizeType& dataSize) { SizeType result; @@ -303,31 +451,81 @@ void SnowMPMSolver::UpdateGridVelocities(double timeStepInSeconds) } template -void SnowMPMSolver::ConstrainGridVelocities() +void SnowMPMSolver::BuildActiveNodes(Array1* activeNodes, + Array1* nodeToActive) const +{ + const auto& gridMass = m_mpmSystemData->GridMass(); + const auto dataView = gridMass.DataView(); + + activeNodes->Clear(); + nodeToActive->Resize(dataView.Length(), ssize_t{ -1 }); + nodeToActive->Fill(ssize_t{ -1 }); + + gridMass.ForEachDataPointIndex([&](const SizeType& index) { + if (gridMass(index) > 0.0) + { + (*nodeToActive)[dataView.Index(index)] = + static_cast(activeNodes->Length()); + activeNodes->Append(index); + } + }); +} + +template +void SnowMPMSolver::ConstrainGridVelocities( + const Array1& activeNodes, const Array1& nodeToActive, + Array1* constrained) { 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(); const auto dataPosition = gridVelocities.DataPosition(); + const auto dataView = gridMass.DataView(); const auto collider = this->GetCollider(); + if (constrained != nullptr) + { + constrained->Resize(activeNodes.Length() * N, uint8_t{ 0 }); + constrained->Fill(uint8_t{ 0 }); + } + gridVelocities.ParallelForEachDataPointIndex( - [&gridMass, &gridVelocities, &dataSize, &dataPosition, &collider, - this](const SizeType& index) { + [&gridMass, &gridVelocities, &dataSize, &dataPosition, &dataView, + &collider, &nodeToActive, constrained, this](const SizeType& index) { if (gridMass(index) <= 0.0) { return; } + const ssize_t active = nodeToActive[dataView.Index(index)]; + + if (active < 0) + { + return; + } + VectorType velocity = gridVelocities(index); + if (collider != nullptr) { + const VectorType incoming = velocity; VectorType position = dataPosition(index); + collider->ResolveCollision(0.0, 0.0, &position, &velocity); + + if (constrained != nullptr && velocity != incoming) + { + for (size_t axis = 0; axis < N; ++axis) + { + (*constrained)[static_cast(active) * N + axis] = + uint8_t{ 1 }; + } + } } for (size_t axis = 0; axis < N; ++axis) @@ -336,11 +534,23 @@ void SnowMPMSolver::ConstrainGridVelocities() index[axis] == 0 && velocity[axis] < 0.0) { velocity[axis] = 0.0; + + if (constrained != nullptr) + { + (*constrained)[static_cast(active) * N + axis] = + uint8_t{ 1 }; + } } if ((m_closedDomainBoundaryFlag & upperFlags[axis]) != 0 && index[axis] == dataSize[axis] - 1 && velocity[axis] > 0.0) { velocity[axis] = 0.0; + + if (constrained != nullptr) + { + (*constrained)[static_cast(active) * N + axis] = + uint8_t{ 1 }; + } } } @@ -349,7 +559,205 @@ void SnowMPMSolver::ConstrainGridVelocities() } template -typename SnowMPMSolver::MatrixType SnowMPMSolver::ComputeVelocityGradient( +void SnowMPMSolver::ApplyElasticHessian(const Array1& activeNodes, + const Array1& nodeToActive, + const VectorND& input, + VectorND* output) const +{ + const auto positions = m_mpmSystemData->Positions(); + const auto volumes = m_mpmSystemData->InitialVolumes(); + const auto states = m_mpmSystemData->DeformationStates(); + const auto& gridMass = m_mpmSystemData->GridMass(); + const auto dataSize = gridMass.DataSize(); + const auto spacing = gridMass.GridSpacing(); + const auto dataOrigin = gridMass.DataOrigin(); + const auto dataView = gridMass.DataView(); + + output->Resize(activeNodes.Length() * N, 0.0); + output->Fill(0.0); + + for (size_t p = 0; p < positions.Length(); ++p) + { + MatrixType differential; + const auto stencil = CubicBSplineKernel::GetStencil( + positions[p], spacing, dataOrigin); + + for (const auto& entry : stencil) + { + if (entry.weight == 0.0) + { + continue; + } + + const SizeType index = ClampIndex(entry.index, dataSize); + const ssize_t active = nodeToActive[dataView.Index(index)]; + + if (active < 0) + { + continue; + } + + VectorType velocityDifferential; + + for (size_t axis = 0; axis < N; ++axis) + { + velocityDifferential[axis] = + input[static_cast(active) * N + axis]; + } + + for (size_t row = 0; row < N; ++row) + { + for (size_t column = 0; column < N; ++column) + { + differential(row, column) += + velocityDifferential[row] * entry.gradient[column]; + } + } + } + + differential *= states[p].elastic; + + const MatrixType stressDifferential = + m_constitutiveModel.ComputeFirstPiolaStressDifferential( + states[p], differential); + + for (const auto& entry : stencil) + { + if (entry.weight == 0.0) + { + continue; + } + + const SizeType index = ClampIndex(entry.index, dataSize); + const ssize_t active = nodeToActive[dataView.Index(index)]; + + if (active < 0) + { + continue; + } + + const VectorType contribution = volumes[p] * stressDifferential * + states[p].elastic.Transposed() * + entry.gradient; + + for (size_t axis = 0; axis < N; ++axis) + { + (*output)[static_cast(active) * N + axis] += + contribution[axis]; + } + } + } +} + +template +void SnowMPMSolver::SolveGridVelocities(double timeStepInSeconds, + const Array1& activeNodes, + const Array1& nodeToActive, + const Array1& constrained) +{ + auto& gridVelocities = m_mpmSystemData->GridVelocities(); + const auto& gridVelocitiesBeforeUpdate = + m_mpmSystemData->GridVelocitiesBeforeUpdate(); + + m_lastNumberOfIterations = 0; + m_lastResidual = 0.0; + + try + { + if (activeNodes.IsEmpty()) + { + return; + } + + const size_t vectorSize = activeNodes.Length() * N; + const double dtSquared = timeStepInSeconds * timeStepInSeconds; + VectorND vStar(vectorSize, 0.0); + + for (size_t active = 0; active < activeNodes.Length(); ++active) + { + const VectorType velocity = gridVelocities(activeNodes[active]); + + for (size_t axis = 0; axis < N; ++axis) + { + vStar[active * N + axis] = velocity[axis]; + } + } + + VectorND hessian(vectorSize, 0.0); + + ApplyElasticHessian(activeNodes, nodeToActive, vStar, &hessian); + + VectorND rhs(vectorSize, 0.0); + + for (size_t i = 0; i < vectorSize; ++i) + { + if (constrained[i] == 0) + { + rhs[i] = -dtSquared * hessian[i]; + } + } + + const double initialResidual = LinearSystemBLAS::L2Norm(rhs); + + if (!std::isfinite(initialResidual)) + { + m_lastResidual = std::numeric_limits::infinity(); + throw std::runtime_error{ + "Semi-implicit snow solve failed to converge." + }; + } + if (initialResidual == 0.0) + { + return; + } + + const LinearSystem system{ this, &activeNodes, &nodeToActive, + &constrained, dtSquared }; + VectorND correction(vectorSize, 0.0); + VectorND residual(vectorSize, 0.0), direction(vectorSize, 0.0), + product(vectorSize, 0.0), image(vectorSize, 0.0); + double residualNorm = initialResidual; + + CR(system, rhs, m_maxNumberOfIterations, + m_tolerance * initialResidual, &correction, + &residual, &direction, &product, &image, + &m_lastNumberOfIterations, &residualNorm); + + m_lastResidual = residualNorm / initialResidual; + + VectorND nextVelocities(vStar + correction); + const bool isFinite = std::ranges::all_of( + nextVelocities, [](double value) { return std::isfinite(value); }); + + if (!isFinite || !std::isfinite(m_lastResidual) || + m_lastResidual > m_tolerance) + { + throw std::runtime_error{ + "Semi-implicit snow solve failed to converge." + }; + } + + for (size_t active = 0; active < activeNodes.Length(); ++active) + { + VectorType velocity; + + for (size_t axis = 0; axis < N; ++axis) + { + velocity[axis] = nextVelocities[active * N + axis]; + } + + gridVelocities(activeNodes[active]) = velocity; + } + } + catch (...) + { + gridVelocities.Set(gridVelocitiesBeforeUpdate); + throw; + } +} + +template +SnowMPMSolver::MatrixType SnowMPMSolver::ComputeVelocityGradient( size_t particleIndex) const { const auto positions = m_mpmSystemData->Positions(); @@ -443,7 +851,7 @@ void SnowMPMSolver::ConstrainParticlesToDomain() } template -typename SnowMPMSolver::Builder& SnowMPMSolver::Builder::WithResolution( +SnowMPMSolver::Builder& SnowMPMSolver::Builder::WithResolution( const SizeType& resolution) { m_resolution = resolution; @@ -451,7 +859,7 @@ typename SnowMPMSolver::Builder& SnowMPMSolver::Builder::WithResolution( } template -typename SnowMPMSolver::Builder& SnowMPMSolver::Builder::WithGridSpacing( +SnowMPMSolver::Builder& SnowMPMSolver::Builder::WithGridSpacing( const VectorType& gridSpacing) { m_gridSpacing = gridSpacing; @@ -459,7 +867,7 @@ typename SnowMPMSolver::Builder& SnowMPMSolver::Builder::WithGridSpacing( } template -typename SnowMPMSolver::Builder& SnowMPMSolver::Builder::WithOrigin( +SnowMPMSolver::Builder& SnowMPMSolver::Builder::WithOrigin( const VectorType& gridOrigin) { m_gridOrigin = gridOrigin; @@ -467,16 +875,14 @@ typename SnowMPMSolver::Builder& SnowMPMSolver::Builder::WithOrigin( } template -typename SnowMPMSolver::Builder& SnowMPMSolver::Builder::WithRadius( - double radius) +SnowMPMSolver::Builder& SnowMPMSolver::Builder::WithRadius(double radius) { m_radius = radius; return *this; } template -typename SnowMPMSolver::Builder& SnowMPMSolver::Builder::WithMass( - double mass) +SnowMPMSolver::Builder& SnowMPMSolver::Builder::WithMass(double mass) { m_mass = mass; return *this; diff --git a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp index cac661cc33..d47798b411 100644 --- a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp +++ b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp @@ -11,6 +11,8 @@ #ifndef CUBBYFLOW_SNOW_MPM_SOLVER_HPP #define CUBBYFLOW_SNOW_MPM_SOLVER_HPP +#include +#include #include #include #include @@ -27,6 +29,8 @@ namespace CubbyFlow //! //! Uses fixed-corotated elastoplastic snow response with the existing MPM //! particle-grid transfers and particle-solver collision lifecycle. +//! When semi-implicit integration is enabled, a non-convergent or non-finite +//! linear solve throws `std::runtime_error` before particle state is advanced. //! template class SnowMPMSolver : public std::conditional_t& index, const SizeType& dataSize); @@ -93,7 +125,21 @@ class SnowMPMSolver : public std::conditional_t* activeNodes, + Array1* nodeToActive) const; + + void ConstrainGridVelocities(const Array1& activeNodes, + const Array1& nodeToActive, + Array1* constrained); + + void SolveGridVelocities(double timeStepInSeconds, + const Array1& activeNodes, + const Array1& nodeToActive, + const Array1& constrained); + + void ApplyElasticHessian(const Array1& activeNodes, + const Array1& nodeToActive, + const VectorND& input, VectorND* output) const; [[nodiscard]] MatrixType ComputeVelocityGradient( size_t particleIndex) const; @@ -106,8 +152,13 @@ class SnowMPMSolver : public std::conditional_t> m_mpmSystemData; SnowConstitutiveModel m_constitutiveModel; + bool m_isUsingSemiImplicit = false; + unsigned int m_maxNumberOfIterations = 100; + unsigned int m_lastNumberOfIterations = 0; double m_timeStepLimitScale = 0.9; double m_maxVelocityGradient = 0.0; + double m_tolerance = 1e-6; + double m_lastResidual = 0.0; int m_closedDomainBoundaryFlag = DIRECTION_ALL; }; diff --git a/Tests/UnitTests/SnowMPMSolverTests.cpp b/Tests/UnitTests/SnowMPMSolverTests.cpp index a34960be8d..845c54a3ed 100644 --- a/Tests/UnitTests/SnowMPMSolverTests.cpp +++ b/Tests/UnitTests/SnowMPMSolverTests.cpp @@ -6,6 +6,7 @@ #include #include +#include #include #include #include @@ -90,6 +91,89 @@ void AddLinearParticleLattice(SnowMPMSolver* solver, double rate) data->AddParticles(positions, velocities); } +template +void AddCompressedParticleLattice(SnowMPMSolver* solver) +{ + auto data = solver->GetMPMSystemData(); + const auto spacing = data->GridMass().GridSpacing(); + const auto origin = data->GridMass().DataOrigin(); + Array1> positions; + + ForEachIndex(VectorUZ::MakeConstant(2), VectorUZ::MakeConstant(8), + [&](auto... rawIndices) { + const VectorUZ index{ rawIndices... }; + VectorD position = origin; + + for (size_t axis = 0; axis < N; ++axis) + { + position[axis] += + spacing[axis] * static_cast(index[axis]); + } + + positions.Append(position); + }); + + data->AddParticles(positions); + + const auto masses = data->ParticleMasses(); + auto volumes = data->InitialVolumes(); + auto states = data->DeformationStates(); + + for (size_t i = 0; i < positions.Length(); ++i) + { + volumes[i] = masses[i] / 400.0; + states[i].elastic(0, 0) = 0.98; + } +} + +template +std::shared_ptr> RunCompressedCase(bool semiImplicit) +{ + SnowMPMSolver solver{ VectorUZ::MakeConstant(10), + VectorD::MakeConstant(0.01) }; + + UseOneFixedStep(&solver); + solver.SetClosedDomainBoundaryFlag(DIRECTION_NONE); + solver.SetIsUsingSemiImplicit(semiImplicit); + + if (semiImplicit) + { + solver.SetTolerance(1e-4); + } + + solver.SetGravity({}); + solver.SetDragCoefficient(0.0); + AddCompressedParticleLattice(&solver); + + for (int frame = 0; frame < 3; ++frame) + { + solver.Update(Frame{ frame, 0.002 }); + } + + return solver.GetMPMSystemData(); +} + +template +double MaxParticleSpeed(const MPMSystemData& data) +{ + double result = 0.0; + + for (const auto& velocity : data.Velocities()) + { + for (double value : velocity) + { + if (!std::isfinite(value)) + { + return std::numeric_limits::infinity(); + } + } + + result = std::max(result, velocity.Length()); + } + + return result; +} + template size_t FindParticle(const MPMSystemData& data, const VectorD& position) { @@ -113,11 +197,150 @@ void ExpectEmptyUpdate() EXPECT_EQ(solver.GetMPMSystemData()->NumberOfParticles(), 0u); } +template +void ExpectEmptySemiImplicitUpdate() +{ + SnowMPMSolver solver; + + solver.SetIsUsingSemiImplicit(true); + solver.Update(Frame{ 0, 0.01 }); + + EXPECT_EQ(solver.GetMPMSystemData()->NumberOfParticles(), 0u); + EXPECT_EQ(solver.GetLastNumberOfIterations(), 0u); + EXPECT_DOUBLE_EQ(solver.GetLastResidual(), 0.0); +} + +template +void ExpectZeroResidualSemiImplicitUpdate() +{ + SnowMPMSolver solver{ VectorUZ::MakeConstant(8) }; + UseOneFixedStep(&solver); + solver.SetIsUsingSemiImplicit(true); + solver.SetGravity({}); + solver.SetDragCoefficient(0.0); + + auto data = solver.GetMPMSystemData(); + data->AddParticle(InteriorPosition(VectorD::MakeConstant(1.0))); + solver.Update(Frame{ 0, 0.01 }); + + EXPECT_TRUE(data->Velocities()[0].IsSimilar(VectorD{}, 0.0)); + EXPECT_EQ(solver.GetLastNumberOfIterations(), 0u); + EXPECT_DOUBLE_EQ(solver.GetLastResidual(), 0.0); +} + +template +void ExpectSmallStepMatchesExplicit() +{ + const auto run = [](bool semiImplicit) { + SnowMPMSolver solver{ VectorUZ::MakeConstant(8) }; + UseOneFixedStep(&solver); + solver.SetIsUsingSemiImplicit(semiImplicit); + solver.SetGravity({}); + solver.SetDragCoefficient(0.0); + + AddLinearParticleLattice(&solver, -0.1); + + solver.Update(Frame{ 0, 1e-8 }); + return solver.GetMPMSystemData(); + }; + + const auto explicitData = run(false); + const auto implicitData = run(true); + + ASSERT_EQ(explicitData->NumberOfParticles(), + implicitData->NumberOfParticles()); + + for (size_t i = 0; i < explicitData->NumberOfParticles(); ++i) + { + EXPECT_TRUE(explicitData->Velocities()[i].IsSimilar( + implicitData->Velocities()[i], 1e-8)); + EXPECT_TRUE(explicitData->DeformationStates()[i].elastic.IsSimilar( + implicitData->DeformationStates()[i].elastic, 1e-8)); + } +} + +template +void ExpectSemiImplicitStiffStability() +{ + const auto explicitData = RunCompressedCase(false); + const auto implicitData = RunCompressedCase(true); + + for (const auto& position : implicitData->Positions()) + { + for (double value : position) + { + EXPECT_TRUE(std::isfinite(value)); + } + } + for (const auto& velocity : implicitData->Velocities()) + { + for (double value : velocity) + { + EXPECT_TRUE(std::isfinite(value)); + } + } + for (const auto& state : implicitData->DeformationStates()) + { + for (double value : state.elastic) + { + EXPECT_TRUE(std::isfinite(value)); + } + for (double value : state.plastic) + { + EXPECT_TRUE(std::isfinite(value)); + } + } + + const double explicitSpeed = MaxParticleSpeed(*explicitData); + const double implicitSpeed = MaxParticleSpeed(*implicitData); + + EXPECT_GT(explicitSpeed, 1.0); + EXPECT_LT(implicitSpeed, explicitSpeed); + EXPECT_LT(implicitSpeed, 1.0); +} + +template +void ExpectSemiImplicitFailureRollback() +{ + SnowMPMSolver solver{ VectorUZ::MakeConstant(10), + VectorD::MakeConstant(0.01) }; + UseOneFixedStep(&solver); + solver.SetClosedDomainBoundaryFlag(DIRECTION_NONE); + solver.SetIsUsingSemiImplicit(true); + solver.SetMaxNumberOfIterations(0); + solver.SetTolerance(1e-12); + solver.SetGravity({}); + solver.SetDragCoefficient(0.0); + + AddCompressedParticleLattice(&solver); + + auto data = solver.GetMPMSystemData(); + const Array1> positionsBefore(data->Positions()); + const Array1> velocitiesBefore(data->Velocities()); + const Array1> statesBefore( + data->DeformationStates()); + + EXPECT_THROW(solver.Update(Frame{ 0, 0.002 }), std::runtime_error); + EXPECT_EQ(solver.GetLastNumberOfIterations(), 0u); + EXPECT_GT(solver.GetLastResidual(), solver.GetTolerance()); + + for (size_t i = 0; i < data->NumberOfParticles(); ++i) + { + EXPECT_EQ(data->Positions()[i], positionsBefore[i]); + EXPECT_EQ(data->Velocities()[i], velocitiesBefore[i]); + EXPECT_EQ(data->DeformationStates()[i].elastic, + statesBefore[i].elastic); + EXPECT_EQ(data->DeformationStates()[i].plastic, + statesBefore[i].plastic); + } +} + template void ExpectReferenceVolumeUsesCellVolume() { VectorD spacing = VectorD::MakeConstant(1.0); spacing[0] = 0.5; + if constexpr (N == 3) { spacing[2] = 2.0; @@ -126,6 +349,7 @@ void ExpectReferenceVolumeUsesCellVolume() SnowMPMSolver solver{ VectorUZ::MakeConstant(8), spacing, {}, 0.1, 2.0 }; + UseOneFixedStep(&solver); solver.SetGravity({}); @@ -134,6 +358,7 @@ void ExpectReferenceVolumeUsesCellVolume() solver.Update(Frame{ 0, 0.001 }); double cellVolume = 1.0; + for (size_t axis = 0; axis < N; ++axis) { cellVolume *= spacing[axis]; @@ -157,6 +382,7 @@ void ExpectUniformMotion() VectorD velocity; velocity[0] = 2.0; + const auto initialPosition = InteriorPosition(VectorD::MakeConstant(1.0)); auto data = solver.GetMPMSystemData(); @@ -177,15 +403,19 @@ void ExpectExternalForcesAppliedOnce() {}, 0.1, 2.0 }; + UseOneFixedStep(&gravitySolver); gravitySolver.SetDragCoefficient(0.0); VectorD gravity; gravity[1] = -2.0; + gravitySolver.SetGravity(gravity); + auto gravityData = gravitySolver.GetMPMSystemData(); gravityData->AddParticle(InteriorPosition(VectorD::MakeConstant(1.0))); gravitySolver.Update(Frame{ 0, 0.1 }); + EXPECT_TRUE(gravityData->Velocities()[0].IsSimilar(0.1 * gravity, 1e-12)); SnowMPMSolver dragSolver{ VectorUZ::MakeConstant(8), @@ -199,11 +429,14 @@ void ExpectExternalForcesAppliedOnce() VectorD velocity; velocity[0] = 1.0; + auto dragData = dragSolver.GetMPMSystemData(); dragData->AddParticle(InteriorPosition(VectorD::MakeConstant(1.0)), velocity); dragSolver.Update(Frame{ 0, 0.1 }); + velocity[0] = 0.9; + EXPECT_TRUE(dragData->Velocities()[0].IsSimilar(velocity, 1e-12)); } @@ -214,17 +447,20 @@ void ExpectLinearVelocityUpdatesDeformation() UseOneFixedStep(&solver); solver.SetGravity({}); solver.SetDragCoefficient(0.0); + AddLinearParticleLattice(&solver, 1.0); auto data = solver.GetMPMSystemData(); const auto target = InteriorPosition(VectorD::MakeConstant(1.0)); const size_t targetIndex = FindParticle(*data, target); + ASSERT_LT(targetIndex, data->NumberOfParticles()); solver.Update(Frame{ 0, 0.001 }); auto expected = Matrix::MakeIdentity(); expected(0, 0) = 1.001; + EXPECT_TRUE(data->DeformationStates()[targetIndex].elastic.IsSimilar( expected, 1e-10)); } @@ -242,7 +478,11 @@ void ExpectAdaptiveRestrictions() baseline.GetMPMSystemData()->InitialVolumes()[0] = baseline.GetMPMSystemData()->ParticleMasses()[0] / 400.0; baseline.Initialize(); + const unsigned int baselineSteps = baseline.NumberOfSubTimeSteps(0.1); + baseline.SetIsUsingSemiImplicit(true); + + EXPECT_EQ(baseline.NumberOfSubTimeSteps(0.1), 1u); TestableSnowMPMSolver finer{ resolution, VectorD::MakeConstant(0.5) }; finer.SetGravity({}); @@ -250,36 +490,54 @@ void ExpectAdaptiveRestrictions() finer.GetMPMSystemData()->InitialVolumes()[0] = finer.GetMPMSystemData()->ParticleMasses()[0] / 400.0; finer.Initialize(); + EXPECT_GT(finer.NumberOfSubTimeSteps(0.1), baselineSteps); TestableSnowMPMSolver fast{ resolution, unit }; fast.SetGravity({}); + VectorD fastVelocity; fastVelocity[0] = 100.0; + fast.GetMPMSystemData()->AddParticle(position, fastVelocity); fast.GetMPMSystemData()->InitialVolumes()[0] = fast.GetMPMSystemData()->ParticleMasses()[0] / 400.0; fast.Initialize(); + EXPECT_GT(fast.NumberOfSubTimeSteps(0.1), baselineSteps); + fast.SetIsUsingSemiImplicit(true); + + EXPECT_GT(fast.NumberOfSubTimeSteps(0.1), 1u); + TestableSnowMPMSolver softerDensity{ resolution, unit }; softerDensity.SetGravity({}); softerDensity.GetMPMSystemData()->AddParticle(position); softerDensity.GetMPMSystemData()->InitialVolumes()[0] = softerDensity.GetMPMSystemData()->ParticleMasses()[0] / 100.0; softerDensity.Initialize(); + EXPECT_GT(softerDensity.NumberOfSubTimeSteps(0.1), baselineSteps); TestableSnowMPMSolver deforming{ resolution, unit }; deforming.SetGravity({}); + AddLinearParticleLattice(&deforming, 100.0); + auto volumes = deforming.GetMPMSystemData()->InitialVolumes(); const auto masses = deforming.GetMPMSystemData()->ParticleMasses(); + for (size_t i = 0; i < volumes.Length(); ++i) { volumes[i] = masses[i] / 400.0; } + deforming.Initialize(); + + EXPECT_EQ(deforming.NumberOfSubTimeSteps(0.1), 56u); + + deforming.SetIsUsingSemiImplicit(true); + EXPECT_EQ(deforming.NumberOfSubTimeSteps(0.1), 56u); } @@ -287,11 +545,30 @@ template void ExpectParametersAndBuilder() { SnowMPMSolver solver; + EXPECT_FALSE(solver.GetIsUsingSemiImplicit()); + + solver.SetIsUsingSemiImplicit(true); + EXPECT_TRUE(solver.GetIsUsingSemiImplicit()); + EXPECT_EQ(solver.GetMaxNumberOfIterations(), 100u); + + solver.SetMaxNumberOfIterations(25); + EXPECT_EQ(solver.GetMaxNumberOfIterations(), 25u); + EXPECT_DOUBLE_EQ(solver.GetTolerance(), 1e-6); + + solver.SetTolerance(1e-8); + EXPECT_DOUBLE_EQ(solver.GetTolerance(), 1e-8); + EXPECT_THROW(solver.SetTolerance(0.0), std::invalid_argument); + EXPECT_THROW(solver.SetTolerance(std::numeric_limits::quiet_NaN()), + std::invalid_argument); + EXPECT_EQ(solver.GetLastNumberOfIterations(), 0u); + EXPECT_DOUBLE_EQ(solver.GetLastResidual(), 0.0); EXPECT_EQ(solver.GetClosedDomainBoundaryFlag(), DIRECTION_ALL); + solver.SetClosedDomainBoundaryFlag(DIRECTION_LEFT | DIRECTION_UP); EXPECT_EQ(solver.GetClosedDomainBoundaryFlag(), DIRECTION_LEFT | DIRECTION_UP); EXPECT_DOUBLE_EQ(solver.GetTimeStepLimitScale(), 0.9); + solver.SetTimeStepLimitScale(0.5); EXPECT_DOUBLE_EQ(solver.GetTimeStepLimitScale(), 0.5); EXPECT_THROW(solver.SetTimeStepLimitScale(0.0), std::invalid_argument); @@ -322,7 +599,7 @@ void ExpectParametersAndBuilder() } template -void ExpectClosedDomainWalls() +void ExpectClosedDomainWalls(bool semiImplicit = false) { constexpr std::array lowerFlags{ DIRECTION_LEFT, DIRECTION_DOWN, DIRECTION_BACK }; @@ -338,13 +615,16 @@ void ExpectClosedDomainWalls() TestableSnowMPMSolver solver{ resolution, spacing }; solver.SetClosedDomainBoundaryFlag(isUpper ? upperFlags[axis] : lowerFlags[axis]); + solver.SetIsUsingSemiImplicit(semiImplicit); solver.SetGravity({}); solver.SetDragCoefficient(0.0); VectorD position = VectorD::MakeConstant(2.0); VectorD velocity; + position[axis] = isUpper ? 3.75 : 0.25; velocity[axis] = isUpper ? 1.0 : -1.0; + auto data = solver.GetMPMSystemData(); data->AddParticle(position, velocity); @@ -353,6 +633,7 @@ void ExpectClosedDomainWalls() VectorUZ nodeIndex = VectorUZ::MakeConstant(2); nodeIndex[axis] = isUpper ? data->GridMass().DataSize()[axis] - 1 : 0; + ASSERT_GT(data->GridMass()(nodeIndex), 0.0); EXPECT_DOUBLE_EQ(data->GridVelocities()(nodeIndex)[axis], 0.0); } @@ -362,13 +643,16 @@ void ExpectClosedDomainWalls() { TestableSnowMPMSolver solver{ resolution, spacing }; solver.SetClosedDomainBoundaryFlag(boundaryFlag); + solver.SetIsUsingSemiImplicit(semiImplicit); solver.SetGravity({}); solver.SetDragCoefficient(0.0); VectorD position = VectorD::MakeConstant(2.0); VectorD velocity; + position[0] = 0.25; velocity[0] = boundaryFlag == DIRECTION_NONE ? -1.0 : 1.0; + auto data = solver.GetMPMSystemData(); data->AddParticle(position, velocity); @@ -376,11 +660,18 @@ void ExpectClosedDomainWalls() VectorUZ nodeIndex = VectorUZ::MakeConstant(2); nodeIndex[0] = 0; + ASSERT_GT(data->GridMass()(nodeIndex), 0.0); EXPECT_DOUBLE_EQ(data->GridVelocities()(nodeIndex)[0], velocity[0]); } } +template +void ExpectSemiImplicitClosedDomainWall() +{ + ExpectClosedDomainWalls(true); +} + template void ExpectMovingColliderAffectsGrid() { @@ -392,13 +683,16 @@ void ExpectMovingColliderAffectsGrid() VectorD normal; normal[0] = 1.0; + auto collider = std::make_shared>( std::make_shared>(normal, VectorD{})); collider->linearVelocity[0] = 1.0; + solver.SetCollider(collider); VectorD position = VectorD::MakeConstant(2.0); position[0] = 0.25; + auto data = solver.GetMPMSystemData(); data->AddParticle(position); @@ -406,6 +700,7 @@ void ExpectMovingColliderAffectsGrid() VectorUZ nodeIndex = VectorUZ::MakeConstant(2); nodeIndex[0] = 0; + ASSERT_GT(data->GridMass()(nodeIndex), 0.0); EXPECT_DOUBLE_EQ(data->GridVelocities()(nodeIndex)[0], 1.0); } @@ -422,17 +717,21 @@ void ExpectFrictionAffectsGrid() VectorD normal; normal[1] = 1.0; + VectorD point; point[1] = 2.0; + auto collider = std::make_shared>( std::make_shared>(normal, point)); collider->SetFrictionCoefficient(frictionCoefficient); + solver.SetCollider(collider); VectorD position = VectorD::MakeConstant(2.25); VectorD velocity; velocity[0] = 1.0; velocity[1] = -1.0; + auto data = solver.GetMPMSystemData(); data->AddParticle(position, velocity); @@ -440,6 +739,7 @@ void ExpectFrictionAffectsGrid() const VectorUZ nodeIndex = VectorUZ::MakeConstant(2); EXPECT_GT(data->GridMass()(nodeIndex), 0.0); + return data->GridVelocities()(nodeIndex); }; @@ -465,10 +765,13 @@ void ExpectParticleDomainProjection() VectorD position = VectorD::MakeConstant(2.0); VectorD velocity; + position[0] = -0.01; velocity[0] = -1.0; + auto data = solver.GetMPMSystemData(); data->AddParticle(position, velocity); + solver.Update(Frame{ 0, 1e-3 }); for (size_t axis = 0; axis < N; ++axis) @@ -508,6 +811,36 @@ TEST(SnowMPMSolver, EmptyUpdate) RUN_FOR_2D_AND_3D(ExpectEmptyUpdate); } +TEST(SnowMPMSolver, EmptySemiImplicitUpdate) +{ + RUN_FOR_2D_AND_3D(ExpectEmptySemiImplicitUpdate); +} + +TEST(SnowMPMSolver, ZeroResidualSemiImplicitUpdate) +{ + RUN_FOR_2D_AND_3D(ExpectZeroResidualSemiImplicitUpdate); +} + +TEST(SnowMPMSolver, SmallStepImplicitAgreement) +{ + RUN_FOR_2D_AND_3D(ExpectSmallStepMatchesExplicit); +} + +TEST(SnowMPMSolver, SemiImplicitStiffStability) +{ + RUN_FOR_2D_AND_3D(ExpectSemiImplicitStiffStability); +} + +TEST(SnowMPMSolver, SemiImplicitFailureRollback) +{ + RUN_FOR_2D_AND_3D(ExpectSemiImplicitFailureRollback); +} + +TEST(SnowMPMSolver, SemiImplicitClosedDomainWall) +{ + RUN_FOR_2D_AND_3D(ExpectSemiImplicitClosedDomainWall); +} + TEST(SnowMPMSolver, ReferenceVolume) { RUN_FOR_2D_AND_3D(ExpectReferenceVolumeUsesCellVolume); From 7e2b2ddc80ed7aa68661f88cdf18284ba6547dd7 Mon Sep 17 00:00:00 2001 From: Chris Ohk Date: Tue, 11 Aug 2026 00:05:36 +0900 Subject: [PATCH 04/13] feat: expose semi-implicit snow controls to Python --- .../Solver/Particle/MPM/SnowMPMSolver.cpp | 11 ++++++++++- Tests/PythonTests/test_snow_mpm_solver.py | 18 ++++++++++++++++++ 2 files changed, 28 insertions(+), 1 deletion(-) diff --git a/Sources/API/Python/Solver/Particle/MPM/SnowMPMSolver.cpp b/Sources/API/Python/Solver/Particle/MPM/SnowMPMSolver.cpp index 2df77d9249..522327c265 100644 --- a/Sources/API/Python/Solver/Particle/MPM/SnowMPMSolver.cpp +++ b/Sources/API/Python/Solver/Particle/MPM/SnowMPMSolver.cpp @@ -35,7 +35,16 @@ void AddSnowMPMSolver(pybind11::module& m, const char* name) &Solver::SetTimeStepLimitScale) .def_property("closedDomainBoundaryFlag", &Solver::GetClosedDomainBoundaryFlag, - &Solver::SetClosedDomainBoundaryFlag); + &Solver::SetClosedDomainBoundaryFlag) + .def_property("isUsingSemiImplicit", &Solver::GetIsUsingSemiImplicit, + &Solver::SetIsUsingSemiImplicit) + .def_property("maxNumberOfIterations", + &Solver::GetMaxNumberOfIterations, + &Solver::SetMaxNumberOfIterations) + .def_property("tolerance", &Solver::GetTolerance, &Solver::SetTolerance) + .def_property_readonly("lastNumberOfIterations", + &Solver::GetLastNumberOfIterations) + .def_property_readonly("lastResidual", &Solver::GetLastResidual); } } // namespace diff --git a/Tests/PythonTests/test_snow_mpm_solver.py b/Tests/PythonTests/test_snow_mpm_solver.py index 121c5c2919..2e420afb60 100644 --- a/Tests/PythonTests/test_snow_mpm_solver.py +++ b/Tests/PythonTests/test_snow_mpm_solver.py @@ -40,4 +40,22 @@ def test_snow_mpm_solver_api( solver.closedDomainBoundaryFlag = 5 assert solver.closedDomainBoundaryFlag == 5 + assert solver.isUsingSemiImplicit is False + solver.isUsingSemiImplicit = True + assert solver.isUsingSemiImplicit is True + + assert solver.maxNumberOfIterations == 100 + solver.maxNumberOfIterations = 25 + assert solver.maxNumberOfIterations == 25 + + assert solver.tolerance == pytest.approx(1e-6) + solver.tolerance = 1e-8 + assert solver.tolerance == pytest.approx(1e-8) + + with pytest.raises(ValueError): + solver.tolerance = 0.0 + + assert solver.lastNumberOfIterations == 0 + assert solver.lastResidual == pytest.approx(0.0) + solver.Update(pyCubbyFlow.Frame(0, 0.001)) From 22a3c06658d1b9d0f39cc02480a69270d748e64e Mon Sep 17 00:00:00 2001 From: Chris Ohk Date: Tue, 11 Aug 2026 09:20:48 +0900 Subject: [PATCH 05/13] fix: clear semi-implicit operator output --- Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp index 03fc6a941d..fb318c7b02 100644 --- a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp +++ b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp @@ -72,6 +72,7 @@ void SnowMPMSolver::LinearSystem::Multiply(const VectorND& input, solver->ApplyElasticHessian(*activeNodes, *nodeToActive, projectedInput, &hessian); output->Resize(input.GetRows(), 0.0); + output->Fill(0.0); const auto& gridMass = solver->m_mpmSystemData->GridMass(); From 26dcc5d14245271267e303a32f93df5a201c50ea Mon Sep 17 00:00:00 2001 From: Chris Ohk Date: Tue, 11 Aug 2026 09:22:25 +0900 Subject: [PATCH 06/13] test: cover conjugate residual breakdown --- Includes/Core/Math/CG.hpp | 17 +++++++++++++++++ Tests/UnitTests/CGTests.cpp | 21 +++++++++++++++++++++ 2 files changed, 38 insertions(+) diff --git a/Includes/Core/Math/CG.hpp b/Includes/Core/Math/CG.hpp index 2dc295836c..9087a1a1e6 100644 --- a/Includes/Core/Math/CG.hpp +++ b/Includes/Core/Math/CG.hpp @@ -51,6 +51,23 @@ void CG(const typename BLASType::MatrixType& A, //! //! \brief Solves a symmetric linear system with conjugate residual. //! +//! Iteration stops when the absolute L2 norm of the residual is at most +//! `tolerance`. On numerical breakdown, `lastResidualNorm` is set to infinity +//! and `x` contains the last valid partial solution. +//! +//! \param[in] A Symmetric system matrix. +//! \param[in] b Right-hand-side vector. +//! \param[in] maxNumberOfIterations Maximum number of iterations. +//! \param[in] tolerance Non-negative absolute residual L2 tolerance. +//! \param[in,out] x Initial guess and computed or partial solution. +//! \param[out] r Residual work vector. +//! \param[out] d Search-direction work vector. +//! \param[out] q Matrix-image work vector. +//! \param[out] s Residual-image work vector. +//! \param[out] lastNumberOfIterations Number of completed iterations. +//! \param[out] lastResidualNorm Final residual L2 norm, or infinity on +//! breakdown. +//! template void CR(const typename BLASType::MatrixType& A, const typename BLASType::VectorType& b, diff --git a/Tests/UnitTests/CGTests.cpp b/Tests/UnitTests/CGTests.cpp index a199adc891..dfc7d19e58 100644 --- a/Tests/UnitTests/CGTests.cpp +++ b/Tests/UnitTests/CGTests.cpp @@ -149,3 +149,24 @@ TEST(CR, AlreadyConverged) EXPECT_EQ(iterations, 0u); EXPECT_LE(residual, 1e-14); } + +TEST(CR, SingularSystemReportsBreakdown) +{ + using BLASType = BLAS; + const Matrix2x2D matrix; + const Vector2D rhs(1.0, 2.0); + Vector2D x; + Vector2D r; + Vector2D d; + Vector2D q; + Vector2D s; + unsigned int iterations = 1; + double residual = 0.0; + + CR(matrix, rhs, 10, 1e-14, &x, &r, &d, &q, &s, &iterations, + &residual); + + EXPECT_EQ(iterations, 0u); + EXPECT_EQ(x, Vector2D{}); + EXPECT_TRUE(std::isinf(residual)); +} From b475c060398b8920ed51509057ef721b17080201 Mon Sep 17 00:00:00 2001 From: Chris Ohk Date: Tue, 11 Aug 2026 09:22:59 +0900 Subject: [PATCH 07/13] docs: clarify snow solver failure contracts --- Includes/Core/Particle/MPM/SnowConstitutiveModel.hpp | 8 ++++++++ Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp | 9 +++++++-- 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/Includes/Core/Particle/MPM/SnowConstitutiveModel.hpp b/Includes/Core/Particle/MPM/SnowConstitutiveModel.hpp index 4fde042ed5..17d05065bb 100644 --- a/Includes/Core/Particle/MPM/SnowConstitutiveModel.hpp +++ b/Includes/Core/Particle/MPM/SnowConstitutiveModel.hpp @@ -96,6 +96,14 @@ class SnowConstitutiveModel final //! Evaluates `(d^2 Psi / d F_E d F_E) : differential` while holding the //! plastic deformation fixed. //! + //! \param[in] state Current elastic and plastic deformation state. + //! \param[in] differential Elastic deformation differential. + //! + //! \return First Piola stress differential. + //! + //! \throws std::invalid_argument If the state, differential, or computed + //! stress differential is invalid or non-finite. + //! [[nodiscard]] MatrixType ComputeFirstPiolaStressDifferential( const State& state, const MatrixType& differential) const; diff --git a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp index d47798b411..846b70fe28 100644 --- a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp +++ b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp @@ -29,8 +29,13 @@ namespace CubbyFlow //! //! Uses fixed-corotated elastoplastic snow response with the existing MPM //! particle-grid transfers and particle-solver collision lifecycle. -//! When semi-implicit integration is enabled, a non-convergent or non-finite -//! linear solve throws `std::runtime_error` before particle state is advanced. +//! When semi-implicit integration is enabled, failures are reported before +//! particle state is advanced. +//! +//! \throws std::runtime_error If the semi-implicit linear solve does not +//! converge or produces a non-finite result. +//! \throws std::invalid_argument If the constitutive state or its stress +//! differential is invalid or non-finite. //! template class SnowMPMSolver : public std::conditional_t Date: Tue, 11 Aug 2026 09:25:23 +0900 Subject: [PATCH 08/13] refactor: clarify semi-implicit solver data flow --- .../Solver/Particle/MPM/SnowMPMSolver-Impl.hpp | 17 ++++++++++------- Tests/UnitTests/CGTests.cpp | 17 ++++++++++++++--- Tests/UnitTests/SnowMPMSolverTests.cpp | 5 +++-- 3 files changed, 27 insertions(+), 12 deletions(-) diff --git a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp index fb318c7b02..e7f6f1d58f 100644 --- a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp +++ b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp @@ -462,7 +462,8 @@ void SnowMPMSolver::BuildActiveNodes(Array1* activeNodes, nodeToActive->Resize(dataView.Length(), ssize_t{ -1 }); nodeToActive->Fill(ssize_t{ -1 }); - gridMass.ForEachDataPointIndex([&](const SizeType& index) { + gridMass.ForEachDataPointIndex([&gridMass, activeNodes, dataView, + nodeToActive](const SizeType& index) { if (gridMass(index) > 0.0) { (*nodeToActive)[dataView.Index(index)] = @@ -715,8 +716,10 @@ void SnowMPMSolver::SolveGridVelocities(double timeStepInSeconds, const LinearSystem system{ this, &activeNodes, &nodeToActive, &constrained, dtSquared }; VectorND correction(vectorSize, 0.0); - VectorND residual(vectorSize, 0.0), direction(vectorSize, 0.0), - product(vectorSize, 0.0), image(vectorSize, 0.0); + VectorND residual(vectorSize, 0.0); + VectorND direction(vectorSize, 0.0); + VectorND product(vectorSize, 0.0); + VectorND image(vectorSize, 0.0); double residualNorm = initialResidual; CR(system, rhs, m_maxNumberOfIterations, @@ -727,10 +730,10 @@ void SnowMPMSolver::SolveGridVelocities(double timeStepInSeconds, m_lastResidual = residualNorm / initialResidual; VectorND nextVelocities(vStar + correction); - const bool isFinite = std::ranges::all_of( - nextVelocities, [](double value) { return std::isfinite(value); }); - - if (!isFinite || !std::isfinite(m_lastResidual) || + if (const bool isFinite = std::ranges::all_of( + nextVelocities, + [](double value) { return std::isfinite(value); }); + !isFinite || !std::isfinite(m_lastResidual) || m_lastResidual > m_tolerance) { throw std::runtime_error{ diff --git a/Tests/UnitTests/CGTests.cpp b/Tests/UnitTests/CGTests.cpp index dfc7d19e58..c419e90ab9 100644 --- a/Tests/UnitTests/CGTests.cpp +++ b/Tests/UnitTests/CGTests.cpp @@ -103,7 +103,11 @@ TEST(CR, SolveSymmetricSystems) { const Vector2D expected(1.0, 1.0); const Vector2D rhs = matrix * expected; - Vector2D x, r, d, q, s; + Vector2D x; + Vector2D r; + Vector2D d; + Vector2D q; + Vector2D s; unsigned int iterations = 0; double residual = 0.0; @@ -121,7 +125,11 @@ TEST(CR, ZeroIterationsReportsInitialResidual) using BLASType = BLAS; const Matrix2x2D matrix(4.0, 1.0, 1.0, 3.0); const Vector2D rhs(1.0, 2.0); - Vector2D x, r, d, q, s; + Vector2D x; + Vector2D r; + Vector2D d; + Vector2D q; + Vector2D s; unsigned int iterations = 1; double residual = 0.0; @@ -139,7 +147,10 @@ TEST(CR, AlreadyConverged) const Matrix2x2D matrix(4.0, 1.0, 1.0, 3.0); Vector2D x(1.0, 1.0); const Vector2D rhs = matrix * x; - Vector2D r, d, q, s; + Vector2D r; + Vector2D d; + Vector2D q; + Vector2D s; unsigned int iterations = 1; double residual = 1.0; diff --git a/Tests/UnitTests/SnowMPMSolverTests.cpp b/Tests/UnitTests/SnowMPMSolverTests.cpp index 845c54a3ed..38c3c2ff9e 100644 --- a/Tests/UnitTests/SnowMPMSolverTests.cpp +++ b/Tests/UnitTests/SnowMPMSolverTests.cpp @@ -74,7 +74,8 @@ void AddLinearParticleLattice(SnowMPMSolver* solver, double rate) Array1> positions; Array1> velocities; - ForEachIndex(dataSize, [&](auto... rawIndices) { + ForEachIndex(dataSize, [center, origin, &positions, rate, spacing, + &velocities](auto... rawIndices) { const VectorUZ index{ rawIndices... }; VectorD position = origin; for (size_t axis = 0; axis < N; ++axis) @@ -100,7 +101,7 @@ void AddCompressedParticleLattice(SnowMPMSolver* solver) Array1> positions; ForEachIndex(VectorUZ::MakeConstant(2), VectorUZ::MakeConstant(8), - [&](auto... rawIndices) { + [origin, &positions, spacing](auto... rawIndices) { const VectorUZ index{ rawIndices... }; VectorD position = origin; From 69ad89eb491484d2c1651bf7afb7e195c7ff78e3 Mon Sep 17 00:00:00 2001 From: Chris Ohk Date: Tue, 11 Aug 2026 09:27:27 +0900 Subject: [PATCH 09/13] refactor: split snow grid constraints --- .../Particle/MPM/SnowMPMSolver-Impl.hpp | 122 +++++++++++------- .../Solver/Particle/MPM/SnowMPMSolver.hpp | 11 ++ 2 files changed, 83 insertions(+), 50 deletions(-) diff --git a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp index e7f6f1d58f..a33451ff95 100644 --- a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp +++ b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp @@ -478,17 +478,9 @@ void SnowMPMSolver::ConstrainGridVelocities( const Array1& activeNodes, const Array1& nodeToActive, Array1* constrained) { - 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(); - const auto dataPosition = gridVelocities.DataPosition(); const auto dataView = gridMass.DataView(); - const auto collider = this->GetCollider(); if (constrained != nullptr) { @@ -497,8 +489,8 @@ void SnowMPMSolver::ConstrainGridVelocities( } gridVelocities.ParallelForEachDataPointIndex( - [&gridMass, &gridVelocities, &dataSize, &dataPosition, &dataView, - &collider, &nodeToActive, constrained, this](const SizeType& index) { + [this, constrained, &gridMass, &nodeToActive, + dataView](const SizeType& index) { if (gridMass(index) <= 0.0) { return; @@ -511,53 +503,83 @@ void SnowMPMSolver::ConstrainGridVelocities( return; } - VectorType velocity = gridVelocities(index); + ConstrainGridVelocityAtNode(index, static_cast(active), + constrained); + }); +} - if (collider != nullptr) - { - const VectorType incoming = velocity; - VectorType position = dataPosition(index); +template +void SnowMPMSolver::ConstrainGridVelocityAtNode(const SizeType& index, + size_t active, + Array1* constrained) +{ + auto& gridVelocities = m_mpmSystemData->GridVelocities(); + VectorType velocity = gridVelocities(index); - collider->ResolveCollision(0.0, 0.0, &position, &velocity); + ApplyGridColliderConstraint(index, active, constrained, &velocity); + ApplyGridDomainConstraint(index, active, constrained, &velocity); - if (constrained != nullptr && velocity != incoming) - { - for (size_t axis = 0; axis < N; ++axis) - { - (*constrained)[static_cast(active) * N + axis] = - uint8_t{ 1 }; - } - } - } + gridVelocities(index) = velocity; +} - for (size_t axis = 0; axis < N; ++axis) - { - if ((m_closedDomainBoundaryFlag & lowerFlags[axis]) != 0 && - index[axis] == 0 && velocity[axis] < 0.0) - { - velocity[axis] = 0.0; +template +void SnowMPMSolver::ApplyGridColliderConstraint(const SizeType& index, + size_t active, + Array1* constrained, + VectorType* velocity) const +{ + const auto collider = this->GetCollider(); + if (collider == nullptr) + { + return; + } - if (constrained != nullptr) - { - (*constrained)[static_cast(active) * N + axis] = - uint8_t{ 1 }; - } - } - if ((m_closedDomainBoundaryFlag & upperFlags[axis]) != 0 && - index[axis] == dataSize[axis] - 1 && velocity[axis] > 0.0) - { - velocity[axis] = 0.0; + const VectorType incoming = *velocity; + VectorType position = + m_mpmSystemData->GridVelocities().DataPosition()(index); - if (constrained != nullptr) - { - (*constrained)[static_cast(active) * N + axis] = - uint8_t{ 1 }; - } - } - } + collider->ResolveCollision(0.0, 0.0, &position, velocity); - gridVelocities(index) = velocity; - }); + if (constrained != nullptr && *velocity != incoming) + { + for (size_t axis = 0; axis < N; ++axis) + { + (*constrained)[active * N + axis] = uint8_t{ 1 }; + } + } +} + +template +void SnowMPMSolver::ApplyGridDomainConstraint(const SizeType& index, + size_t active, + Array1* constrained, + VectorType* velocity) const +{ + static constexpr std::array lowerFlags{ DIRECTION_LEFT, DIRECTION_DOWN, + DIRECTION_BACK }; + static constexpr std::array upperFlags{ DIRECTION_RIGHT, DIRECTION_UP, + DIRECTION_FRONT }; + const auto dataSize = m_mpmSystemData->GridVelocities().DataSize(); + + 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; + + if (constrained != nullptr) + { + (*constrained)[active * N + axis] = uint8_t{ 1 }; + } + } + } } template diff --git a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp index 846b70fe28..7b6f8820ba 100644 --- a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp +++ b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp @@ -137,6 +137,17 @@ class SnowMPMSolver : public std::conditional_t& nodeToActive, Array1* constrained); + void ConstrainGridVelocityAtNode(const SizeType& index, size_t active, + Array1* constrained); + + void ApplyGridColliderConstraint(const SizeType& index, size_t active, + Array1* constrained, + VectorType* velocity) const; + + void ApplyGridDomainConstraint(const SizeType& index, size_t active, + Array1* constrained, + VectorType* velocity) const; + void SolveGridVelocities(double timeStepInSeconds, const Array1& activeNodes, const Array1& nodeToActive, From eea4e6bf2da8ef174c35c9a33ff55091e5cacd7f Mon Sep 17 00:00:00 2001 From: Chris Ohk Date: Tue, 11 Aug 2026 09:31:19 +0900 Subject: [PATCH 10/13] refactor: split snow Hessian application --- .../Particle/MPM/SnowMPMSolver-Impl.hpp | 151 ++++++++++-------- .../Solver/Particle/MPM/SnowMPMSolver.hpp | 11 ++ 2 files changed, 99 insertions(+), 63 deletions(-) diff --git a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp index a33451ff95..afbffda214 100644 --- a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp +++ b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp @@ -583,96 +583,121 @@ void SnowMPMSolver::ApplyGridDomainConstraint(const SizeType& index, } template -void SnowMPMSolver::ApplyElasticHessian(const Array1& activeNodes, - const Array1& nodeToActive, - const VectorND& input, - VectorND* output) const +SnowMPMSolver::MatrixType +SnowMPMSolver::ComputeParticleDeformationDifferential( + const Stencil& stencil, const Array1& nodeToActive, + const VectorND& input, const MatrixType& elastic) const { - const auto positions = m_mpmSystemData->Positions(); - const auto volumes = m_mpmSystemData->InitialVolumes(); - const auto states = m_mpmSystemData->DeformationStates(); const auto& gridMass = m_mpmSystemData->GridMass(); const auto dataSize = gridMass.DataSize(); - const auto spacing = gridMass.GridSpacing(); - const auto dataOrigin = gridMass.DataOrigin(); const auto dataView = gridMass.DataView(); + MatrixType result; - output->Resize(activeNodes.Length() * N, 0.0); - output->Fill(0.0); - - for (size_t p = 0; p < positions.Length(); ++p) + for (const auto& entry : stencil) { - MatrixType differential; - const auto stencil = CubicBSplineKernel::GetStencil( - positions[p], spacing, dataOrigin); - - for (const auto& entry : stencil) + if (entry.weight == 0.0) { - if (entry.weight == 0.0) - { - continue; - } + continue; + } - const SizeType index = ClampIndex(entry.index, dataSize); - const ssize_t active = nodeToActive[dataView.Index(index)]; + const SizeType index = ClampIndex(entry.index, dataSize); + const ssize_t active = nodeToActive[dataView.Index(index)]; - if (active < 0) - { - continue; - } + if (active < 0) + { + continue; + } - VectorType velocityDifferential; + VectorType velocityDifferential; - for (size_t axis = 0; axis < N; ++axis) - { - velocityDifferential[axis] = - input[static_cast(active) * N + axis]; - } + for (size_t axis = 0; axis < N; ++axis) + { + velocityDifferential[axis] = + input[static_cast(active) * N + axis]; + } - for (size_t row = 0; row < N; ++row) + for (size_t row = 0; row < N; ++row) + { + for (size_t column = 0; column < N; ++column) { - for (size_t column = 0; column < N; ++column) - { - differential(row, column) += - velocityDifferential[row] * entry.gradient[column]; - } + result(row, column) += + velocityDifferential[row] * entry.gradient[column]; } } + } - differential *= states[p].elastic; + result *= elastic; + return result; +} - const MatrixType stressDifferential = - m_constitutiveModel.ComputeFirstPiolaStressDifferential( - states[p], differential); +template +void SnowMPMSolver::AccumulateParticleHessian( + const Stencil& stencil, const Array1& nodeToActive, double volume, + const MatrixType& elastic, const MatrixType& stressDifferential, + VectorND* output) const +{ + const auto& gridMass = m_mpmSystemData->GridMass(); + const auto dataSize = gridMass.DataSize(); + const auto dataView = gridMass.DataView(); - for (const auto& entry : stencil) + for (const auto& entry : stencil) + { + if (entry.weight == 0.0) { - if (entry.weight == 0.0) - { - continue; - } + continue; + } - const SizeType index = ClampIndex(entry.index, dataSize); - const ssize_t active = nodeToActive[dataView.Index(index)]; + const SizeType index = ClampIndex(entry.index, dataSize); + const ssize_t active = nodeToActive[dataView.Index(index)]; - if (active < 0) - { - continue; - } + if (active < 0) + { + continue; + } - const VectorType contribution = volumes[p] * stressDifferential * - states[p].elastic.Transposed() * - entry.gradient; + const VectorType contribution = + volume * stressDifferential * elastic.Transposed() * entry.gradient; - for (size_t axis = 0; axis < N; ++axis) - { - (*output)[static_cast(active) * N + axis] += - contribution[axis]; - } + for (size_t axis = 0; axis < N; ++axis) + { + (*output)[static_cast(active) * N + axis] += + contribution[axis]; } } } +template +void SnowMPMSolver::ApplyElasticHessian(const Array1& activeNodes, + const Array1& nodeToActive, + const VectorND& input, + VectorND* output) const +{ + const auto positions = m_mpmSystemData->Positions(); + const auto volumes = m_mpmSystemData->InitialVolumes(); + const auto states = m_mpmSystemData->DeformationStates(); + const auto& gridMass = m_mpmSystemData->GridMass(); + const auto spacing = gridMass.GridSpacing(); + const auto dataOrigin = gridMass.DataOrigin(); + + output->Resize(activeNodes.Length() * N, 0.0); + output->Fill(0.0); + + for (size_t p = 0; p < positions.Length(); ++p) + { + const auto stencil = CubicBSplineKernel::GetStencil( + positions[p], spacing, dataOrigin); + const MatrixType differential = ComputeParticleDeformationDifferential( + stencil, nodeToActive, input, states[p].elastic); + const MatrixType stressDifferential = + m_constitutiveModel.ComputeFirstPiolaStressDifferential( + states[p], differential); + + AccumulateParticleHessian(stencil, nodeToActive, volumes[p], + states[p].elastic, stressDifferential, + output); + } +} + template void SnowMPMSolver::SolveGridVelocities(double timeStepInSeconds, const Array1& activeNodes, diff --git a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp index 7b6f8820ba..1c11512850 100644 --- a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp +++ b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp @@ -122,6 +122,7 @@ class SnowMPMSolver : public std::conditional_t::Stencil; [[nodiscard]] static SizeType ClampIndex(const Vector& index, const SizeType& dataSize); @@ -157,6 +158,16 @@ class SnowMPMSolver : public std::conditional_t& nodeToActive, const VectorND& input, VectorND* output) const; + [[nodiscard]] MatrixType ComputeParticleDeformationDifferential( + const Stencil& stencil, const Array1& nodeToActive, + const VectorND& input, const MatrixType& elastic) const; + + void AccumulateParticleHessian(const Stencil& stencil, + const Array1& nodeToActive, + double volume, const MatrixType& elastic, + const MatrixType& stressDifferential, + VectorND* output) const; + [[nodiscard]] MatrixType ComputeVelocityGradient( size_t particleIndex) const; From 8c82c1e9946e99e8a7d8d2e9ca82f36db570b4fc Mon Sep 17 00:00:00 2001 From: Chris Ohk Date: Tue, 11 Aug 2026 09:33:46 +0900 Subject: [PATCH 11/13] refactor: split snow grid solve phases --- .../Particle/MPM/SnowMPMSolver-Impl.hpp | 201 +++++++++++------- .../Solver/Particle/MPM/SnowMPMSolver.hpp | 21 ++ 2 files changed, 144 insertions(+), 78 deletions(-) diff --git a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp index afbffda214..1769aec9f9 100644 --- a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp +++ b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp @@ -699,105 +699,150 @@ void SnowMPMSolver::ApplyElasticHessian(const Array1& activeNodes, } template -void SnowMPMSolver::SolveGridVelocities(double timeStepInSeconds, - const Array1& activeNodes, - const Array1& nodeToActive, - const Array1& constrained) +VectorND SnowMPMSolver::GatherActiveGridVelocities( + const Array1& activeNodes) const { - auto& gridVelocities = m_mpmSystemData->GridVelocities(); - const auto& gridVelocitiesBeforeUpdate = - m_mpmSystemData->GridVelocitiesBeforeUpdate(); - - m_lastNumberOfIterations = 0; - m_lastResidual = 0.0; + const auto& gridVelocities = m_mpmSystemData->GridVelocities(); + VectorND result(activeNodes.Length() * N, 0.0); - try + for (size_t active = 0; active < activeNodes.Length(); ++active) { - if (activeNodes.IsEmpty()) + const VectorType velocity = gridVelocities(activeNodes[active]); + + for (size_t axis = 0; axis < N; ++axis) { - return; + result[active * N + axis] = velocity[axis]; } + } - const size_t vectorSize = activeNodes.Length() * N; - const double dtSquared = timeStepInSeconds * timeStepInSeconds; - VectorND vStar(vectorSize, 0.0); + return result; +} - for (size_t active = 0; active < activeNodes.Length(); ++active) - { - const VectorType velocity = gridVelocities(activeNodes[active]); +template +VectorND SnowMPMSolver::BuildSemiImplicitRightHandSide( + double dtSquared, const Array1& activeNodes, + const Array1& nodeToActive, const Array1& constrained, + const VectorND& velocities) const +{ + VectorND hessian(velocities.GetRows(), 0.0); + ApplyElasticHessian(activeNodes, nodeToActive, velocities, &hessian); - for (size_t axis = 0; axis < N; ++axis) - { - vStar[active * N + axis] = velocity[axis]; - } + VectorND result(velocities.GetRows(), 0.0); + for (size_t i = 0; i < result.GetRows(); ++i) + { + if (constrained[i] == 0) + { + result[i] = -dtSquared * hessian[i]; } + } - VectorND hessian(vectorSize, 0.0); + return result; +} - ApplyElasticHessian(activeNodes, nodeToActive, vStar, &hessian); +template +VectorND SnowMPMSolver::SolveGridVelocityCorrection( + double dtSquared, const Array1& activeNodes, + const Array1& nodeToActive, const Array1& constrained, + const VectorND& rhs, double initialResidual) +{ + const LinearSystem system{ this, &activeNodes, &nodeToActive, &constrained, + dtSquared }; + VectorND correction(rhs.GetRows(), 0.0); + VectorND residual(rhs.GetRows(), 0.0); + VectorND direction(rhs.GetRows(), 0.0); + VectorND product(rhs.GetRows(), 0.0); + VectorND image(rhs.GetRows(), 0.0); + double residualNorm = initialResidual; - VectorND rhs(vectorSize, 0.0); + CR(system, rhs, m_maxNumberOfIterations, + m_tolerance * initialResidual, &correction, &residual, + &direction, &product, &image, + &m_lastNumberOfIterations, &residualNorm); - for (size_t i = 0; i < vectorSize; ++i) - { - if (constrained[i] == 0) - { - rhs[i] = -dtSquared * hessian[i]; - } - } + m_lastResidual = residualNorm / initialResidual; + return correction; +} - const double initialResidual = LinearSystemBLAS::L2Norm(rhs); +template +VectorND SnowMPMSolver::ComputeGridVelocityUpdate( + double timeStepInSeconds, const Array1& activeNodes, + const Array1& nodeToActive, const Array1& constrained) +{ + const double dtSquared = timeStepInSeconds * timeStepInSeconds; + const VectorND vStar = GatherActiveGridVelocities(activeNodes); + const VectorND rhs = BuildSemiImplicitRightHandSide( + dtSquared, activeNodes, nodeToActive, constrained, vStar); + const double initialResidual = LinearSystemBLAS::L2Norm(rhs); - if (!std::isfinite(initialResidual)) - { - m_lastResidual = std::numeric_limits::infinity(); - throw std::runtime_error{ - "Semi-implicit snow solve failed to converge." - }; - } - if (initialResidual == 0.0) - { - return; - } + if (!std::isfinite(initialResidual)) + { + m_lastResidual = std::numeric_limits::infinity(); + throw std::runtime_error{ + "Semi-implicit snow solve failed to converge." + }; + } + if (initialResidual == 0.0) + { + return vStar; + } + + const VectorND correction = + SolveGridVelocityCorrection(dtSquared, activeNodes, nodeToActive, + constrained, rhs, initialResidual); + VectorND result(vStar + correction); - const LinearSystem system{ this, &activeNodes, &nodeToActive, - &constrained, dtSquared }; - VectorND correction(vectorSize, 0.0); - VectorND residual(vectorSize, 0.0); - VectorND direction(vectorSize, 0.0); - VectorND product(vectorSize, 0.0); - VectorND image(vectorSize, 0.0); - double residualNorm = initialResidual; - - CR(system, rhs, m_maxNumberOfIterations, - m_tolerance * initialResidual, &correction, - &residual, &direction, &product, &image, - &m_lastNumberOfIterations, &residualNorm); - - m_lastResidual = residualNorm / initialResidual; - - VectorND nextVelocities(vStar + correction); - if (const bool isFinite = std::ranges::all_of( - nextVelocities, - [](double value) { return std::isfinite(value); }); - !isFinite || !std::isfinite(m_lastResidual) || - m_lastResidual > m_tolerance) + if (const bool isFinite = std::ranges::all_of( + result, [](double value) { return std::isfinite(value); }); + !isFinite || !std::isfinite(m_lastResidual) || + m_lastResidual > m_tolerance) + { + throw std::runtime_error{ + "Semi-implicit snow solve failed to converge." + }; + } + + return result; +} + +template +void SnowMPMSolver::StoreActiveGridVelocities( + const Array1& activeNodes, const VectorND& velocities) +{ + auto& gridVelocities = m_mpmSystemData->GridVelocities(); + + for (size_t active = 0; active < activeNodes.Length(); ++active) + { + VectorType velocity; + + for (size_t axis = 0; axis < N; ++axis) { - throw std::runtime_error{ - "Semi-implicit snow solve failed to converge." - }; + velocity[axis] = velocities[active * N + axis]; } - for (size_t active = 0; active < activeNodes.Length(); ++active) - { - VectorType velocity; + gridVelocities(activeNodes[active]) = velocity; + } +} - for (size_t axis = 0; axis < N; ++axis) - { - velocity[axis] = nextVelocities[active * N + axis]; - } +template +void SnowMPMSolver::SolveGridVelocities(double timeStepInSeconds, + const Array1& activeNodes, + const Array1& nodeToActive, + const Array1& constrained) +{ + auto& gridVelocities = m_mpmSystemData->GridVelocities(); + const auto& gridVelocitiesBeforeUpdate = + m_mpmSystemData->GridVelocitiesBeforeUpdate(); + + m_lastNumberOfIterations = 0; + m_lastResidual = 0.0; - gridVelocities(activeNodes[active]) = velocity; + try + { + if (!activeNodes.IsEmpty()) + { + const VectorND nextVelocities = ComputeGridVelocityUpdate( + timeStepInSeconds, activeNodes, nodeToActive, constrained); + StoreActiveGridVelocities(activeNodes, nextVelocities); } } catch (...) diff --git a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp index 1c11512850..2a106171a4 100644 --- a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp +++ b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp @@ -154,6 +154,27 @@ class SnowMPMSolver : public std::conditional_t& nodeToActive, const Array1& constrained); + [[nodiscard]] VectorND GatherActiveGridVelocities( + const Array1& activeNodes) const; + + [[nodiscard]] VectorND BuildSemiImplicitRightHandSide( + double dtSquared, const Array1& activeNodes, + const Array1& nodeToActive, const Array1& constrained, + const VectorND& velocities) const; + + [[nodiscard]] VectorND SolveGridVelocityCorrection( + double dtSquared, const Array1& activeNodes, + const Array1& nodeToActive, const Array1& constrained, + const VectorND& rhs, double initialResidual); + + [[nodiscard]] VectorND ComputeGridVelocityUpdate( + double timeStepInSeconds, const Array1& activeNodes, + const Array1& nodeToActive, + const Array1& constrained); + + void StoreActiveGridVelocities(const Array1& activeNodes, + const VectorND& velocities); + void ApplyElasticHessian(const Array1& activeNodes, const Array1& nodeToActive, const VectorND& input, VectorND* output) const; From 6e981e03cfc6ad518baee4ec052b102bf7f10840 Mon Sep 17 00:00:00 2001 From: Chris Ohk Date: Tue, 11 Aug 2026 09:35:08 +0900 Subject: [PATCH 12/13] test: split snow boundary scenarios --- Tests/UnitTests/SnowMPMSolverTests.cpp | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/Tests/UnitTests/SnowMPMSolverTests.cpp b/Tests/UnitTests/SnowMPMSolverTests.cpp index 38c3c2ff9e..ec94e1801b 100644 --- a/Tests/UnitTests/SnowMPMSolverTests.cpp +++ b/Tests/UnitTests/SnowMPMSolverTests.cpp @@ -600,7 +600,7 @@ void ExpectParametersAndBuilder() } template -void ExpectClosedDomainWalls(bool semiImplicit = false) +void ExpectClosedDomainWallStopsVelocity(bool semiImplicit) { constexpr std::array lowerFlags{ DIRECTION_LEFT, DIRECTION_DOWN, DIRECTION_BACK }; @@ -639,6 +639,13 @@ void ExpectClosedDomainWalls(bool semiImplicit = false) EXPECT_DOUBLE_EQ(data->GridVelocities()(nodeIndex)[axis], 0.0); } } +} + +template +void ExpectUnconstrainedDomainWallPreservesVelocity(bool semiImplicit) +{ + const auto resolution = VectorUZ::MakeConstant(4); + const auto spacing = VectorD::MakeConstant(1.0); for (int boundaryFlag : { DIRECTION_NONE, DIRECTION_LEFT }) { @@ -667,6 +674,13 @@ void ExpectClosedDomainWalls(bool semiImplicit = false) } } +template +void ExpectClosedDomainWalls(bool semiImplicit = false) +{ + ExpectClosedDomainWallStopsVelocity(semiImplicit); + ExpectUnconstrainedDomainWallPreservesVelocity(semiImplicit); +} + template void ExpectSemiImplicitClosedDomainWall() { From ccc4307c50b6f3c587628744d35a202120a9fa1b Mon Sep 17 00:00:00 2001 From: Chris Ohk Date: Tue, 11 Aug 2026 10:33:21 +0900 Subject: [PATCH 13/13] refactor: improve code readability --- Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp index 1769aec9f9..62a9f8a605 100644 --- a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp +++ b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp @@ -529,6 +529,7 @@ void SnowMPMSolver::ApplyGridColliderConstraint(const SizeType& index, VectorType* velocity) const { const auto collider = this->GetCollider(); + if (collider == nullptr) { return; @@ -725,9 +726,11 @@ VectorND SnowMPMSolver::BuildSemiImplicitRightHandSide( const VectorND& velocities) const { VectorND hessian(velocities.GetRows(), 0.0); + ApplyElasticHessian(activeNodes, nodeToActive, velocities, &hessian); VectorND result(velocities.GetRows(), 0.0); + for (size_t i = 0; i < result.GetRows(); ++i) { if (constrained[i] == 0)