From d135973120ede60b71c0faf91108f520086311c2 Mon Sep 17 00:00:00 2001 From: Chris Ohk Date: Sat, 8 Aug 2026 21:53:26 +0900 Subject: [PATCH 1/4] feat: add snow elastic wave speed estimate --- .../MPM/SnowConstitutiveModel-Impl.hpp | 82 +++++++++++++++++-- .../Particle/MPM/SnowConstitutiveModel.hpp | 15 ++++ .../UnitTests/SnowConstitutiveModelTests.cpp | 76 +++++++++++++++++ 3 files changed, 166 insertions(+), 7 deletions(-) diff --git a/Includes/Core/Particle/MPM/SnowConstitutiveModel-Impl.hpp b/Includes/Core/Particle/MPM/SnowConstitutiveModel-Impl.hpp index df0efb1bcf..416a1fa29d 100644 --- a/Includes/Core/Particle/MPM/SnowConstitutiveModel-Impl.hpp +++ b/Includes/Core/Particle/MPM/SnowConstitutiveModel-Impl.hpp @@ -113,13 +113,7 @@ SnowConstitutiveModel::ComputeKirchhoffStress(const State& state) const const MatrixType rotation = u * v.Transposed(); const double elasticDeterminant = state.elastic.Determinant(); - const double hardening = - std::exp(m_hardeningCoefficient * (1.0 - state.plastic.Determinant())); - - if (!std::isfinite(hardening)) - { - throw std::invalid_argument{ "Non-finite snow hardening." }; - } + const double hardening = ComputeHardening(state); const double mu = m_mu0 * hardening; const double lambda = m_lambda0 * hardening; @@ -136,6 +130,80 @@ SnowConstitutiveModel::ComputeKirchhoffStress(const State& state) const return stress; } +template +double SnowConstitutiveModel::ComputeWaveSpeed(const State& state, + double referenceDensity) const +{ + ValidateDeformation(state.elastic); + ValidateDeformation(state.plastic); + + if (!std::isfinite(referenceDensity) || referenceDensity <= 0.0) + { + throw std::invalid_argument{ "Invalid snow reference density." }; + } + + MatrixType u; + Vector singularValues; + MatrixType v; + SVD(state.elastic, u, singularValues, v); + + const double elasticDeterminant = state.elastic.Determinant(); + const double totalDeterminant = + elasticDeterminant * state.plastic.Determinant(); + const double currentDensity = referenceDensity / totalDeterminant; + const double hardening = ComputeHardening(state); + const double mu = m_mu0 * hardening; + const double lambda = m_lambda0 * hardening; + double maxCandidate = 0.0; + + for (size_t a = 0; a < N; ++a) + { + const double sigmaA = singularValues[a]; + maxCandidate = std::max( + maxCandidate, 2.0 * mu * sigmaA * sigmaA + + lambda * elasticDeterminant * elasticDeterminant); + + for (size_t b = 0; b < N; ++b) + { + if (a == b) + { + continue; + } + + const double sigmaB = singularValues[b]; + maxCandidate = std::max(maxCandidate, 2.0 * mu * sigmaB * sigmaB * + (sigmaA + sigmaB - 1.0) / + (sigmaA + sigmaB)); + } + } + + const double kappa = maxCandidate / totalDeterminant; + const double waveSpeed = std::sqrt(kappa / currentDensity); + + if (!std::isfinite(totalDeterminant) || totalDeterminant <= 0.0 || + !std::isfinite(currentDensity) || currentDensity <= 0.0 || + !std::isfinite(kappa) || kappa <= 0.0 || !std::isfinite(waveSpeed)) + { + throw std::invalid_argument{ "Invalid snow wave speed." }; + } + + return waveSpeed; +} + +template +double SnowConstitutiveModel::ComputeHardening(const State& state) const +{ + const double hardening = + std::exp(m_hardeningCoefficient * (1.0 - state.plastic.Determinant())); + + if (!std::isfinite(hardening)) + { + throw std::invalid_argument{ "Non-finite snow hardening." }; + } + + return hardening; +} + template bool SnowConstitutiveModel::IsFinite(const MatrixType& matrix) { diff --git a/Includes/Core/Particle/MPM/SnowConstitutiveModel.hpp b/Includes/Core/Particle/MPM/SnowConstitutiveModel.hpp index 51d8c714ca..936f6f042a 100644 --- a/Includes/Core/Particle/MPM/SnowConstitutiveModel.hpp +++ b/Includes/Core/Particle/MPM/SnowConstitutiveModel.hpp @@ -90,7 +90,22 @@ class SnowConstitutiveModel final //! \return Fixed-corotated Kirchhoff stress. [[nodiscard]] MatrixType ComputeKirchhoffStress(const State& state) const; + //! + //! \brief Estimates the fastest elastic wave speed for the state. + //! + //! Evaluates the principal-basis acoustic-tensor candidates for the + //! fixed-corotated snow energy. + //! + //! \param[in] state Current deformation state. + //! \param[in] referenceDensity Positive material density at rest. + //! + //! \return Estimated elastic wave speed. + [[nodiscard]] double ComputeWaveSpeed(const State& state, + double referenceDensity) const; + private: + [[nodiscard]] double ComputeHardening(const State& state) const; + [[nodiscard]] static bool IsFinite(const MatrixType& matrix); static void ValidateDeformation(const MatrixType& matrix); diff --git a/Tests/UnitTests/SnowConstitutiveModelTests.cpp b/Tests/UnitTests/SnowConstitutiveModelTests.cpp index 19056092fe..f9342c3756 100644 --- a/Tests/UnitTests/SnowConstitutiveModelTests.cpp +++ b/Tests/UnitTests/SnowConstitutiveModelTests.cpp @@ -161,6 +161,62 @@ void ExpectPlasticCompressionHardens() 1e-10); } +template +void ExpectAnalyticWaveSpeed() +{ + const SnowConstitutiveModel model{ 1000.0, 0.2, 0.025, 0.0075, 0.0 }; + + EXPECT_NEAR(model.ComputeWaveSpeed({}, 400.0), 5.0 / 3.0, 1e-12); + + SnowDeformationState state; + state.elastic = MakeStretch(1.005); + + EXPECT_NEAR(model.ComputeWaveSpeed(state, 400.0), 1.675, 1e-12); +} + +template +void ExpectWaveSpeedIncludesHardening() +{ + const SnowConstitutiveModel model{ 1000.0, 0.2, 0.025, 0.0075, 10.0 }; + + SnowDeformationState unpacked; + unpacked.elastic = MakeStretch(1.005); + + auto packed = unpacked; + packed.plastic = MakeStretch(0.8); + + EXPECT_NEAR(model.ComputeWaveSpeed(packed, 400.0) / + model.ComputeWaveSpeed(unpacked, 400.0), + std::exp(1.0), 1e-10); +} + +template +void ExpectSpecialWaveSpeedCandidate() +{ + const SnowConstitutiveModel model{ 1000.0, -0.5, 0.025, 0.0075, 0.0 }; + SnowDeformationState state; + state.elastic = 2.0 * MatrixD::MakeIdentity(); + + EXPECT_NEAR(model.ComputeWaveSpeed(state, 600.0), std::sqrt(10.0), 1e-12); +} + +template +void ExpectInvalidWaveSpeedInputRejected() +{ + const SnowConstitutiveModel model; + + EXPECT_THROW((void)model.ComputeWaveSpeed({}, 0.0), std::invalid_argument); + EXPECT_THROW((void)model.ComputeWaveSpeed( + {}, std::numeric_limits::quiet_NaN()), + std::invalid_argument); + + SnowDeformationState inverted; + inverted.elastic = MakeStretch(-1.0); + + EXPECT_THROW((void)model.ComputeWaveSpeed(inverted, 400.0), + std::invalid_argument); +} + template void ExpectRotatedStretchProjectsAndSoftens() { @@ -322,6 +378,26 @@ TEST(SnowConstitutiveModel, Hardening) RUN_FOR_2D_AND_3D(ExpectPlasticCompressionHardens); } +TEST(SnowConstitutiveModel, WaveSpeed) +{ + RUN_FOR_2D_AND_3D(ExpectAnalyticWaveSpeed); +} + +TEST(SnowConstitutiveModel, WaveSpeedHardening) +{ + RUN_FOR_2D_AND_3D(ExpectWaveSpeedIncludesHardening); +} + +TEST(SnowConstitutiveModel, SpecialWaveSpeedCandidate) +{ + RUN_FOR_2D_AND_3D(ExpectSpecialWaveSpeedCandidate); +} + +TEST(SnowConstitutiveModel, InvalidWaveSpeedInput) +{ + RUN_FOR_2D_AND_3D(ExpectInvalidWaveSpeedInputRejected); +} + TEST(SnowConstitutiveModel, RotatedStretchProjectionAndSoftening) { RUN_FOR_2D_AND_3D(ExpectRotatedStretchProjectsAndSoftens); From 55f01c14065e956fa580fd4c8efe7e5c018147fd Mon Sep 17 00:00:00 2001 From: Chris Ohk Date: Sat, 8 Aug 2026 21:53:50 +0900 Subject: [PATCH 2/4] feat: implement CPU snow MPM solvers --- .../Particle/MPM/SnowMPMSolver-Impl.hpp | 399 ++++++++++++++++++ .../Solver/Particle/MPM/SnowMPMSolver.hpp | 133 ++++++ Tests/UnitTests/SnowMPMSolverTests.cpp | 358 ++++++++++++++++ 3 files changed, 890 insertions(+) create mode 100644 Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp create mode 100644 Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp create mode 100644 Tests/UnitTests/SnowMPMSolverTests.cpp diff --git a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp new file mode 100644 index 0000000000..6c3a98f516 --- /dev/null +++ b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver-Impl.hpp @@ -0,0 +1,399 @@ +// 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. + +#ifndef CUBBYFLOW_SNOW_MPM_SOLVER_IMPL_HPP +#define CUBBYFLOW_SNOW_MPM_SOLVER_IMPL_HPP + +#include +#include +#include +#include + +namespace CubbyFlow +{ +template +SnowMPMSolver::SnowMPMSolver(const SizeType& resolution, + const VectorType& gridSpacing, + const VectorType& gridOrigin, double radius, + double mass) + : Base{ radius, mass }, + m_mpmSystemData{ std::make_shared>( + resolution, gridSpacing, gridOrigin) } +{ + if (!std::isfinite(radius) || radius < 0.0 || !std::isfinite(mass) || + mass <= 0.0) + { + throw std::invalid_argument{ "Invalid snow MPM particle parameters." }; + } + + m_mpmSystemData->SetRadius(radius); + m_mpmSystemData->SetMass(mass); + this->SetParticleSystemData(m_mpmSystemData); + this->SetIsUsingFixedSubTimeSteps(false); +} + +template +std::shared_ptr> SnowMPMSolver::GetMPMSystemData() const +{ + return m_mpmSystemData; +} + +template +double SnowMPMSolver::GetTimeStepLimitScale() const +{ + return m_timeStepLimitScale; +} + +template +void SnowMPMSolver::SetTimeStepLimitScale(double newScale) +{ + if (!std::isfinite(newScale) || newScale <= 0.0 || newScale > 1.0) + { + throw std::invalid_argument{ + "Time-step limit scale must be in (0, 1]." + }; + } + + m_timeStepLimitScale = newScale; +} + +template +typename SnowMPMSolver::Builder SnowMPMSolver::GetBuilder() +{ + return Builder{}; +} + +template +void SnowMPMSolver::OnInitialize() +{ + Base::OnInitialize(); + m_mpmSystemData->TransferFromParticlesToGrid(); + InitializeReferenceVolumes(); + m_maxVelocityGradient = ComputeMaxVelocityGradient(); +} + +template +unsigned int SnowMPMSolver::GetNumberOfSubTimeSteps( + double timeIntervalInSeconds) const +{ + if (!std::isfinite(timeIntervalInSeconds) || timeIntervalInSeconds <= 0.0) + { + return 1; + } + + const auto velocities = m_mpmSystemData->Velocities(); + const auto masses = m_mpmSystemData->ParticleMasses(); + const auto volumes = m_mpmSystemData->InitialVolumes(); + const auto states = m_mpmSystemData->DeformationStates(); + const auto spacing = m_mpmSystemData->GridMass().GridSpacing(); + double minSpacing = spacing[0]; + double maxVelocity = 0.0; + double maxWaveSpeed = 0.0; + + for (size_t axis = 1; axis < N; ++axis) + { + minSpacing = std::min(minSpacing, spacing[axis]); + } + + for (size_t i = 0; i < velocities.Length(); ++i) + { + maxVelocity = std::max(maxVelocity, velocities[i].Length()); + + if (volumes[i] == 0.0) + { + continue; + } + + if (!std::isfinite(volumes[i]) || volumes[i] < 0.0) + { + throw std::invalid_argument{ "Invalid snow reference volume." }; + } + + maxWaveSpeed = + std::max(maxWaveSpeed, m_constitutiveModel.ComputeWaveSpeed( + states[i], masses[i] / volumes[i])); + } + + double desiredTimeStep = std::numeric_limits::infinity(); + if (maxVelocity > 0.0) + { + desiredTimeStep = std::min(desiredTimeStep, minSpacing / maxVelocity); + } + if (maxWaveSpeed > 0.0) + { + desiredTimeStep = std::min(desiredTimeStep, minSpacing / maxWaveSpeed); + } + if (m_maxVelocityGradient > 0.0) + { + desiredTimeStep = + std::min(desiredTimeStep, 0.2 / m_maxVelocityGradient); + } + + if (!std::isfinite(desiredTimeStep)) + { + return 1; + } + + const double count = std::ceil(timeIntervalInSeconds / + (m_timeStepLimitScale * desiredTimeStep)); + const double maxCount = + static_cast(std::numeric_limits::max()); + + return static_cast(std::clamp(count, 1.0, maxCount)); +} + +template +void SnowMPMSolver::AccumulateForces(double timeStepInSeconds) +{ + (void)timeStepInSeconds; +} + +template +void SnowMPMSolver::OnBeginAdvanceTimeStep(double timeStepInSeconds) +{ + m_mpmSystemData->TransferFromParticlesToGrid(); + InitializeReferenceVolumes(); + UpdateGridVelocities(timeStepInSeconds); + UpdateDeformation(timeStepInSeconds); + m_mpmSystemData->TransferFromGridToParticles(); +} + +template +typename SnowMPMSolver::SizeType SnowMPMSolver::ClampIndex( + const Vector& index, const SizeType& dataSize) +{ + SizeType result; + + for (size_t axis = 0; axis < N; ++axis) + { + result[axis] = static_cast(std::clamp( + index[axis], 0, static_cast(dataSize[axis] - 1))); + } + + return result; +} + +template +void SnowMPMSolver::InitializeReferenceVolumes() +{ + const auto positions = m_mpmSystemData->Positions(); + const auto masses = m_mpmSystemData->ParticleMasses(); + auto volumes = m_mpmSystemData->InitialVolumes(); + const auto& gridMass = m_mpmSystemData->GridMass(); + const auto dataSize = gridMass.DataSize(); + const auto spacing = gridMass.GridSpacing(); + const auto dataOrigin = gridMass.DataOrigin(); + double cellVolume = 1.0; + + for (size_t axis = 0; axis < N; ++axis) + { + cellVolume *= spacing[axis]; + } + + if (!std::isfinite(cellVolume) || cellVolume <= 0.0) + { + throw std::invalid_argument{ "Invalid snow MPM cell volume." }; + } + + for (size_t i = 0; i < volumes.Length(); ++i) + { + if (volumes[i] != 0.0) + { + if (!std::isfinite(volumes[i]) || volumes[i] < 0.0) + { + throw std::invalid_argument{ "Invalid snow reference volume." }; + } + continue; + } + + double referenceDensity = 0.0; + const auto stencil = CubicBSplineKernel::GetStencil( + positions[i], spacing, dataOrigin); + + for (const auto& entry : stencil) + { + referenceDensity += gridMass(ClampIndex(entry.index, dataSize)) * + entry.weight / cellVolume; + } + + if (!std::isfinite(referenceDensity) || referenceDensity <= 0.0) + { + throw std::invalid_argument{ "Invalid snow reference density." }; + } + + volumes[i] = masses[i] / referenceDensity; + } +} + +template +void SnowMPMSolver::UpdateGridVelocities(double timeStepInSeconds) +{ + const auto positions = m_mpmSystemData->Positions(); + const auto velocities = m_mpmSystemData->Velocities(); + const auto masses = m_mpmSystemData->ParticleMasses(); + const auto volumes = m_mpmSystemData->InitialVolumes(); + const auto states = m_mpmSystemData->DeformationStates(); + const auto& gridMass = m_mpmSystemData->GridMass(); + auto& gridVelocities = m_mpmSystemData->GridVelocities(); + const auto dataSize = gridMass.DataSize(); + const auto spacing = gridMass.GridSpacing(); + const auto dataOrigin = gridMass.DataOrigin(); + + for (size_t i = 0; i < positions.Length(); ++i) + { + const VectorType relativeVelocity = + velocities[i] - this->GetWind()->Sample(positions[i]); + const VectorType externalForce = + masses[i] * this->GetGravity() - + this->GetDragCoefficient() * relativeVelocity; + const MatrixType stress = + m_constitutiveModel.ComputeKirchhoffStress(states[i]); + const auto stencil = CubicBSplineKernel::GetStencil( + positions[i], spacing, dataOrigin); + + for (const auto& entry : stencil) + { + if (entry.weight == 0.0) + { + continue; + } + + const auto index = ClampIndex(entry.index, dataSize); + const double nodeMass = gridMass(index); + if (nodeMass > 0.0) + { + gridVelocities(index) += + timeStepInSeconds * + (entry.weight * externalForce - + volumes[i] * stress * entry.gradient) / + nodeMass; + } + } + } +} + +template +typename SnowMPMSolver::MatrixType SnowMPMSolver::ComputeVelocityGradient( + size_t particleIndex) const +{ + const auto positions = m_mpmSystemData->Positions(); + const auto& gridVelocities = m_mpmSystemData->GridVelocities(); + const auto dataSize = gridVelocities.DataSize(); + const auto stencil = CubicBSplineKernel::GetStencil( + positions[particleIndex], gridVelocities.GridSpacing(), + gridVelocities.DataOrigin()); + MatrixType result; + + for (const auto& entry : stencil) + { + if (entry.weight != 0.0) + { + const VectorType velocity = + gridVelocities(ClampIndex(entry.index, dataSize)); + for (size_t row = 0; row < N; ++row) + { + for (size_t column = 0; column < N; ++column) + { + result(row, column) += + velocity[row] * entry.gradient[column]; + } + } + } + } + + return result; +} + +template +double SnowMPMSolver::ComputeMaxVelocityGradient() const +{ + double result = 0.0; + + for (size_t i = 0; i < m_mpmSystemData->NumberOfParticles(); ++i) + { + result = std::max(result, ComputeVelocityGradient(i).AbsMax()); + } + + return result; +} + +template +void SnowMPMSolver::UpdateDeformation(double timeStepInSeconds) +{ + auto states = m_mpmSystemData->DeformationStates(); + m_maxVelocityGradient = 0.0; + + for (size_t i = 0; i < states.Length(); ++i) + { + const MatrixType velocityGradient = ComputeVelocityGradient(i); + m_maxVelocityGradient = + std::max(m_maxVelocityGradient, velocityGradient.AbsMax()); + states[i] = m_constitutiveModel.Update( + MatrixType::MakeIdentity() + timeStepInSeconds * velocityGradient, + states[i]); + } +} + +template +typename SnowMPMSolver::Builder& SnowMPMSolver::Builder::WithResolution( + const SizeType& resolution) +{ + m_resolution = resolution; + return *this; +} + +template +typename SnowMPMSolver::Builder& SnowMPMSolver::Builder::WithGridSpacing( + const VectorType& gridSpacing) +{ + m_gridSpacing = gridSpacing; + return *this; +} + +template +typename SnowMPMSolver::Builder& SnowMPMSolver::Builder::WithOrigin( + const VectorType& gridOrigin) +{ + m_gridOrigin = gridOrigin; + return *this; +} + +template +typename SnowMPMSolver::Builder& SnowMPMSolver::Builder::WithRadius( + double radius) +{ + m_radius = radius; + return *this; +} + +template +typename SnowMPMSolver::Builder& SnowMPMSolver::Builder::WithMass( + double mass) +{ + m_mass = mass; + return *this; +} + +template +SnowMPMSolver SnowMPMSolver::Builder::Build() const +{ + return SnowMPMSolver{ m_resolution, m_gridSpacing, m_gridOrigin, m_radius, + m_mass }; +} + +template +std::shared_ptr> SnowMPMSolver::Builder::MakeShared() const +{ + return std::make_shared(m_resolution, m_gridSpacing, + m_gridOrigin, m_radius, m_mass); +} +} // namespace CubbyFlow + +#endif diff --git a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp new file mode 100644 index 0000000000..06acb5009c --- /dev/null +++ b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp @@ -0,0 +1,133 @@ +// 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. + +#ifndef CUBBYFLOW_SNOW_MPM_SOLVER_HPP +#define CUBBYFLOW_SNOW_MPM_SOLVER_HPP + +#include +#include +#include + +#include +#include + +namespace CubbyFlow +{ +//! +//! \brief N-D material point method solver for snow. +//! +//! Uses fixed-corotated elastoplastic snow response with the existing MPM +//! particle-grid transfers and particle-solver collision lifecycle. +//! +template +class SnowMPMSolver : public std::conditional_t +{ + public: + static_assert(N == 2 || N == 3, "Snow MPM supports only 2-D and 3-D."); + + using Base = std::conditional_t; + using MatrixType = Matrix; + using VectorType = Vector; + using SizeType = Vector; + + class Builder; + + //! Constructs an empty snow MPM solver. + explicit SnowMPMSolver( + const SizeType& resolution = SizeType::MakeConstant(32), + const VectorType& gridSpacing = VectorType::MakeConstant(1.0), + const VectorType& gridOrigin = VectorType{}, double radius = 1e-3, + double mass = 1e-3); + + //! Returns the owned MPM particle and grid state. + [[nodiscard]] std::shared_ptr> GetMPMSystemData() const; + + //! Returns the scale applied to adaptive time-step limits. + [[nodiscard]] double GetTimeStepLimitScale() const; + + //! Sets the adaptive time-step scale in `(0, 1]`. + void SetTimeStepLimitScale(double newScale); + + //! Returns a builder for SnowMPMSolver. + [[nodiscard]] static Builder GetBuilder(); + + protected: + //! Initializes particle-grid state before the first adaptive step query. + void OnInitialize() override; + + //! Returns the number of adaptive sub-time-steps. + [[nodiscard]] unsigned int GetNumberOfSubTimeSteps( + double timeIntervalInSeconds) const override; + + //! Snow forces are integrated on the background grid. + void AccumulateForces(double timeStepInSeconds) override; + + //! Advances particle-grid snow state before base particle integration. + void OnBeginAdvanceTimeStep(double timeStepInSeconds) override; + + private: + [[nodiscard]] static SizeType ClampIndex(const Vector& index, + const SizeType& dataSize); + + void InitializeReferenceVolumes(); + + void UpdateGridVelocities(double timeStepInSeconds); + + [[nodiscard]] MatrixType ComputeVelocityGradient( + size_t particleIndex) const; + + [[nodiscard]] double ComputeMaxVelocityGradient() const; + + void UpdateDeformation(double timeStepInSeconds); + + std::shared_ptr> m_mpmSystemData; + SnowConstitutiveModel m_constitutiveModel; + double m_timeStepLimitScale = 0.9; + double m_maxVelocityGradient = 0.0; +}; + +//! Front-end to create SnowMPMSolver objects step by step. +template +class SnowMPMSolver::Builder final +{ + public: + [[nodiscard]] Builder& WithResolution(const SizeType& resolution); + + [[nodiscard]] Builder& WithGridSpacing(const VectorType& gridSpacing); + + [[nodiscard]] Builder& WithOrigin(const VectorType& gridOrigin); + + [[nodiscard]] Builder& WithRadius(double radius); + + [[nodiscard]] Builder& WithMass(double mass); + + [[nodiscard]] SnowMPMSolver Build() const; + + [[nodiscard]] std::shared_ptr MakeShared() const; + + private: + SizeType m_resolution = SizeType::MakeConstant(32); + VectorType m_gridSpacing = VectorType::MakeConstant(1.0); + VectorType m_gridOrigin; + double m_radius = 1e-3; + double m_mass = 1e-3; +}; + +using SnowMPMSolver2 = SnowMPMSolver<2>; +using SnowMPMSolver3 = SnowMPMSolver<3>; +using SnowMPMSolver2Ptr = std::shared_ptr; +using SnowMPMSolver3Ptr = std::shared_ptr; +} // namespace CubbyFlow + +#include + +#endif diff --git a/Tests/UnitTests/SnowMPMSolverTests.cpp b/Tests/UnitTests/SnowMPMSolverTests.cpp new file mode 100644 index 0000000000..1942d15ccd --- /dev/null +++ b/Tests/UnitTests/SnowMPMSolverTests.cpp @@ -0,0 +1,358 @@ +#include "gtest/gtest.h" + +#include +#include + +#include +#include +#include + +using namespace CubbyFlow; + +namespace +{ +template +using VectorD = Vector; + +template +using VectorUZ = Vector; + +template +class TestableSnowMPMSolver final : public SnowMPMSolver +{ + public: + using SnowMPMSolver::SnowMPMSolver; + + void Initialize() + { + this->OnInitialize(); + } + + [[nodiscard]] unsigned int NumberOfSubTimeSteps(double interval) const + { + return this->GetNumberOfSubTimeSteps(interval); + } +}; + +template +VectorD InteriorPosition(const VectorD& spacing) +{ + VectorD result; + for (size_t axis = 0; axis < N; ++axis) + { + result[axis] = 3.0 * spacing[axis]; + } + return result; +} + +template +void UseOneFixedStep(SnowMPMSolver* solver) +{ + solver->SetIsUsingFixedSubTimeSteps(true); + solver->SetNumberOfFixedSubTimeSteps(1); +} + +template +void AddLinearParticleLattice(SnowMPMSolver* solver, double rate) +{ + auto data = solver->GetMPMSystemData(); + const auto dataSize = data->GridMass().DataSize(); + const auto spacing = data->GridMass().GridSpacing(); + const auto origin = data->GridMass().DataOrigin(); + const double center = + origin[0] + 0.5 * static_cast(dataSize[0] - 1) * spacing[0]; + Array1> positions; + Array1> velocities; + + ForEachIndex(dataSize, [&](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]); + } + + VectorD velocity; + velocity[0] = rate * (position[0] - center); + positions.Append(position); + velocities.Append(velocity); + }); + + data->AddParticles(positions, velocities); +} + +template +size_t FindParticle(const MPMSystemData& data, const VectorD& position) +{ + const auto positions = data.Positions(); + for (size_t i = 0; i < positions.Length(); ++i) + { + if (positions[i].IsSimilar(position, 1e-12)) + { + return i; + } + } + + return positions.Length(); +} + +template +void ExpectEmptyUpdate() +{ + SnowMPMSolver solver; + solver.Update(Frame{ 0, 0.01 }); + EXPECT_EQ(solver.GetMPMSystemData()->NumberOfParticles(), 0u); +} + +template +void ExpectReferenceVolumeUsesCellVolume() +{ + VectorD spacing = VectorD::MakeConstant(1.0); + spacing[0] = 0.5; + if constexpr (N == 3) + { + spacing[2] = 2.0; + } + + SnowMPMSolver solver{ + VectorUZ::MakeConstant(8), spacing, {}, 0.1, 2.0 + }; + UseOneFixedStep(&solver); + solver.SetGravity({}); + + auto data = solver.GetMPMSystemData(); + data->AddParticle(InteriorPosition(spacing)); + solver.Update(Frame{ 0, 0.001 }); + + double cellVolume = 1.0; + for (size_t axis = 0; axis < N; ++axis) + { + cellVolume *= spacing[axis]; + } + + EXPECT_NEAR(data->InitialVolumes()[0], + cellVolume * static_cast(size_t{ 1 } << N), 1e-12); + EXPECT_TRUE(data->DeformationStates()[0].elastic.IsSimilar( + Matrix::MakeIdentity(), 1e-12)); + EXPECT_TRUE(data->DeformationStates()[0].plastic.IsSimilar( + Matrix::MakeIdentity(), 1e-12)); +} + +template +void ExpectUniformMotion() +{ + SnowMPMSolver solver{ VectorUZ::MakeConstant(8) }; + UseOneFixedStep(&solver); + solver.SetGravity({}); + solver.SetDragCoefficient(0.0); + + VectorD velocity; + velocity[0] = 2.0; + const auto initialPosition = + InteriorPosition(VectorD::MakeConstant(1.0)); + auto data = solver.GetMPMSystemData(); + data->AddParticle(initialPosition, velocity); + + solver.Update(Frame{ 0, 0.01 }); + + EXPECT_TRUE(data->Velocities()[0].IsSimilar(velocity, 1e-12)); + EXPECT_TRUE(data->Positions()[0].IsSimilar( + initialPosition + 0.01 * velocity, 1e-12)); +} + +template +void ExpectExternalForcesAppliedOnce() +{ + SnowMPMSolver gravitySolver{ VectorUZ::MakeConstant(8), + VectorD::MakeConstant(1.0), + {}, + 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), + VectorD::MakeConstant(1.0), + {}, + 0.1, + 2.0 }; + UseOneFixedStep(&dragSolver); + dragSolver.SetGravity({}); + dragSolver.SetDragCoefficient(2.0); + + 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)); +} + +template +void ExpectLinearVelocityUpdatesDeformation() +{ + SnowMPMSolver solver{ VectorUZ::MakeConstant(8) }; + 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)); +} + +template +void ExpectAdaptiveRestrictions() +{ + const auto resolution = VectorUZ::MakeConstant(8); + const auto unit = VectorD::MakeConstant(1.0); + const auto position = InteriorPosition(unit); + + TestableSnowMPMSolver baseline{ resolution, unit }; + baseline.SetGravity({}); + baseline.GetMPMSystemData()->AddParticle(position); + baseline.GetMPMSystemData()->InitialVolumes()[0] = + baseline.GetMPMSystemData()->ParticleMasses()[0] / 400.0; + baseline.Initialize(); + const unsigned int baselineSteps = baseline.NumberOfSubTimeSteps(0.1); + + TestableSnowMPMSolver finer{ resolution, VectorD::MakeConstant(0.5) }; + finer.SetGravity({}); + finer.GetMPMSystemData()->AddParticle(0.5 * position); + 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); + + 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); +} + +template +void ExpectParametersAndBuilder() +{ + SnowMPMSolver solver; + 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); + EXPECT_THROW(solver.SetTimeStepLimitScale(1.1), std::invalid_argument); + EXPECT_THROW( + solver.SetTimeStepLimitScale(std::numeric_limits::quiet_NaN()), + std::invalid_argument); + + const auto resolution = VectorUZ::MakeConstant(7); + const auto spacing = VectorD::MakeConstant(0.25); + const auto origin = VectorD::MakeConstant(-1.0); + auto built = SnowMPMSolver::GetBuilder() + .WithResolution(resolution) + .WithGridSpacing(spacing) + .WithOrigin(origin) + .WithRadius(0.2) + .WithMass(3.0) + .Build(); + + const auto data = built.GetMPMSystemData(); + EXPECT_EQ(data->GridMass().Resolution(), resolution); + EXPECT_EQ(data->GridMass().GridSpacing(), spacing); + EXPECT_EQ(data->GridMass().Origin(), origin); + EXPECT_DOUBLE_EQ(data->Radius(), 0.2); + EXPECT_DOUBLE_EQ(data->Mass(), 3.0); + + EXPECT_NE(SnowMPMSolver::GetBuilder().MakeShared(), nullptr); +} +} // namespace + +#define RUN_FOR_2D_AND_3D(function) \ + do \ + { \ + { \ + SCOPED_TRACE("N = 2"); \ + function<2>(); \ + } \ + { \ + SCOPED_TRACE("N = 3"); \ + function<3>(); \ + } \ + } while (false) + +TEST(SnowMPMSolver, EmptyUpdate) +{ + RUN_FOR_2D_AND_3D(ExpectEmptyUpdate); +} + +TEST(SnowMPMSolver, ReferenceVolume) +{ + RUN_FOR_2D_AND_3D(ExpectReferenceVolumeUsesCellVolume); +} + +TEST(SnowMPMSolver, UniformMotion) +{ + RUN_FOR_2D_AND_3D(ExpectUniformMotion); +} + +TEST(SnowMPMSolver, ExternalForces) +{ + RUN_FOR_2D_AND_3D(ExpectExternalForcesAppliedOnce); +} + +TEST(SnowMPMSolver, LinearVelocityGradient) +{ + RUN_FOR_2D_AND_3D(ExpectLinearVelocityUpdatesDeformation); +} + +TEST(SnowMPMSolver, AdaptiveRestrictions) +{ + RUN_FOR_2D_AND_3D(ExpectAdaptiveRestrictions); +} + +TEST(SnowMPMSolver, ParametersAndBuilder) +{ + RUN_FOR_2D_AND_3D(ExpectParametersAndBuilder); +} From 3f3efef99b7747f4b81955e4bb3c09bb21ad50af Mon Sep 17 00:00:00 2001 From: Chris Ohk Date: Sat, 8 Aug 2026 21:54:13 +0900 Subject: [PATCH 3/4] feat: expose snow MPM solvers to Python --- .../Solver/Particle/MPM/SnowMPMSolver.hpp | 19 ++++++++ .../Solver/Particle/MPM/SnowMPMSolver.cpp | 47 +++++++++++++++++++ Sources/API/Python/main.cpp | 3 ++ Tests/PythonTests/test_snow_mpm_solver.py | 39 +++++++++++++++ 4 files changed, 108 insertions(+) create mode 100644 Includes/API/Python/Solver/Particle/MPM/SnowMPMSolver.hpp create mode 100644 Sources/API/Python/Solver/Particle/MPM/SnowMPMSolver.cpp create mode 100644 Tests/PythonTests/test_snow_mpm_solver.py diff --git a/Includes/API/Python/Solver/Particle/MPM/SnowMPMSolver.hpp b/Includes/API/Python/Solver/Particle/MPM/SnowMPMSolver.hpp new file mode 100644 index 0000000000..c2c2ad2486 --- /dev/null +++ b/Includes/API/Python/Solver/Particle/MPM/SnowMPMSolver.hpp @@ -0,0 +1,19 @@ +// 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. + +#ifndef CUBBYFLOW_PYTHON_SNOW_MPM_SOLVER_HPP +#define CUBBYFLOW_PYTHON_SNOW_MPM_SOLVER_HPP + +#include + +void AddSnowMPMSolver2(pybind11::module& m); +void AddSnowMPMSolver3(pybind11::module& m); + +#endif diff --git a/Sources/API/Python/Solver/Particle/MPM/SnowMPMSolver.cpp b/Sources/API/Python/Solver/Particle/MPM/SnowMPMSolver.cpp new file mode 100644 index 0000000000..12e9afb4b4 --- /dev/null +++ b/Sources/API/Python/Solver/Particle/MPM/SnowMPMSolver.cpp @@ -0,0 +1,47 @@ +// 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 +#include + +using namespace CubbyFlow; + +namespace +{ +template +void AddSnowMPMSolver(pybind11::module& m, const char* name) +{ + using Solver = SnowMPMSolver; + using SolverPtr = std::shared_ptr; + using VectorType = Vector; + using SizeType = Vector; + + pybind11::class_(m, name) + .def(pybind11::init(), + pybind11::arg("resolution") = SizeType::MakeConstant(32), + pybind11::arg("gridSpacing") = VectorType::MakeConstant(1.0), + pybind11::arg("gridOrigin") = VectorType{}, + pybind11::arg("radius") = 1e-3, pybind11::arg("mass") = 1e-3) + .def_property_readonly("mpmSystemData", &Solver::GetMPMSystemData) + .def_property("timeStepLimitScale", &Solver::GetTimeStepLimitScale, + &Solver::SetTimeStepLimitScale); +} +} // namespace + +void AddSnowMPMSolver2(pybind11::module& m) +{ + AddSnowMPMSolver<2, ParticleSystemSolver2>(m, "SnowMPMSolver2"); +} + +void AddSnowMPMSolver3(pybind11::module& m) +{ + AddSnowMPMSolver<3, ParticleSystemSolver3>(m, "SnowMPMSolver3"); +} diff --git a/Sources/API/Python/main.cpp b/Sources/API/Python/main.cpp index cce67f6902..27f6303537 100644 --- a/Sources/API/Python/main.cpp +++ b/Sources/API/Python/main.cpp @@ -92,6 +92,7 @@ #include #include #include +#include #include #include #include @@ -331,6 +332,8 @@ PYBIND11_MODULE(pyCubbyFlow, m) AddAPICSolver3(m); AddParticleSystemSolver2(m); AddParticleSystemSolver3(m); + AddSnowMPMSolver2(m); + AddSnowMPMSolver3(m); AddSPHSolver2(m); AddSPHSolver3(m); AddPCISPHSolver2(m); diff --git a/Tests/PythonTests/test_snow_mpm_solver.py b/Tests/PythonTests/test_snow_mpm_solver.py new file mode 100644 index 0000000000..3bb68f14b4 --- /dev/null +++ b/Tests/PythonTests/test_snow_mpm_solver.py @@ -0,0 +1,39 @@ +import pytest + +import pyCubbyFlow + + +@pytest.mark.parametrize( + "solver_type,data_type,resolution,spacing,origin", + [ + ( + pyCubbyFlow.SnowMPMSolver2, + pyCubbyFlow.MPMSystemData2, + pyCubbyFlow.Vector2UZ(4, 4), + pyCubbyFlow.Vector2D(0.5, 1.0), + pyCubbyFlow.Vector2D(-1.0, -2.0), + ), + ( + pyCubbyFlow.SnowMPMSolver3, + pyCubbyFlow.MPMSystemData3, + pyCubbyFlow.Vector3UZ(4, 4, 4), + pyCubbyFlow.Vector3D(0.5, 1.0, 2.0), + pyCubbyFlow.Vector3D(-1.0, -2.0, -3.0), + ), + ], +) +def test_snow_mpm_solver_api( + solver_type, data_type, resolution, spacing, origin +): + solver = solver_type(resolution, spacing, origin, 0.1, 2.0) + + assert isinstance(solver.mpmSystemData, data_type) + assert solver.mpmSystemData.gridMass.resolution == resolution + assert solver.mpmSystemData.gridMass.gridSpacing == spacing + assert solver.mpmSystemData.gridMass.gridOrigin == origin + + assert solver.timeStepLimitScale == pytest.approx(0.9) + solver.timeStepLimitScale = 0.5 + assert solver.timeStepLimitScale == pytest.approx(0.5) + + solver.Update(pyCubbyFlow.Frame(0, 0.001)) From be2cd1fa18255203dfc227ae0283878de83a4265 Mon Sep 17 00:00:00 2001 From: Chris Ohk Date: Sat, 8 Aug 2026 23:34:53 +0900 Subject: [PATCH 4/4] fix: include snow constitutive model directly --- Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp index 06acb5009c..a6060edcac 100644 --- a/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp +++ b/Includes/Core/Solver/Particle/MPM/SnowMPMSolver.hpp @@ -12,6 +12,7 @@ #define CUBBYFLOW_SNOW_MPM_SOLVER_HPP #include +#include #include #include