Commit a34d9928 authored by Nghia Truong's avatar Nghia Truong

BUG: Fix SPH boundary condition

This bug is due to invalid computation of boundary friction, causing incorrect damping (negative damping) that accelerate the colliding particles instead of slowing them down.
parent e48e52a3
......@@ -28,15 +28,16 @@ using namespace imstk;
///
/// \brief Generate a sphere-shape fluid object
///
StdVectorOfVec3d generateSphereShapeFluid(const double particleRadius)
StdVectorOfVec3d
generateSphereShapeFluid(const double particleRadius)
{
const double sphereRadius = 2.0;
const Vec3d sphereCenter(0, 1, 0);
const Vec3d sphereCenter(0, 1, 0);
const auto sphereRadiusSqr = sphereRadius * sphereRadius;
const auto spacing = 2.0 * particleRadius;
const auto N = static_cast<size_t>(2.0 * sphereRadius / spacing); // Maximum number of particles in each dimension
const Vec3d lcorner = sphereCenter - Vec3d(sphereRadius, sphereRadius, sphereRadius); // Cannot use auto here, due to Eigen bug
const auto sphereRadiusSqr = sphereRadius * sphereRadius;
const auto spacing = 2.0 * particleRadius;
const auto N = static_cast<size_t>(2.0 * sphereRadius / spacing); // Maximum number of particles in each dimension
const Vec3d lcorner = sphereCenter - Vec3d(sphereRadius, sphereRadius, sphereRadius); // Cannot use auto here, due to Eigen bug
StdVectorOfVec3d particles;
particles.reserve(N * N * N);
......@@ -63,10 +64,11 @@ StdVectorOfVec3d generateSphereShapeFluid(const double particleRadius)
///
/// \brief Generate a box-shape fluid object
///
StdVectorOfVec3d generateBoxShapeFluid(const double particleRadius)
StdVectorOfVec3d
generateBoxShapeFluid(const double particleRadius)
{
const double boxWidth = 4.0;
const Vec3d boxLowerCorner(-2, -3, -2);
const Vec3d boxLowerCorner(-2, -3, -2);
const auto spacing = 2.0 * particleRadius;
const auto N = static_cast<size_t>(boxWidth / spacing);
......@@ -95,7 +97,8 @@ StdVectorOfVec3d getBunny(); // Defined in Bunny.cpp
///
/// \brief Generate a bunny-shape fluid object
///
StdVectorOfVec3d generateBunnyShapeFluid(const double particleRadius)
StdVectorOfVec3d
generateBunnyShapeFluid(const double particleRadius)
{
LOG_IF(FATAL, (std::abs(particleRadius - 0.08) > 1e-6)) << "Particle radius for this scene must be 0.08";
StdVectorOfVec3d particles;
......@@ -105,7 +108,8 @@ StdVectorOfVec3d generateBunnyShapeFluid(const double particleRadius)
return particles;
}
std::shared_ptr<SPHObject> generateFluid(const std::shared_ptr<Scene>&scene, const double particleRadius)
std::shared_ptr<SPHObject>
generateFluid(const std::shared_ptr<Scene>& scene, const double particleRadius)
{
StdVectorOfVec3d particles;
switch (SCENE_ID)
......@@ -134,7 +138,7 @@ std::shared_ptr<SPHObject> generateFluid(const std::shared_ptr<Scene>&scene, con
// Create a visiual model
auto fluidVisualModel = std::make_shared<VisualModel>(fluidGeometry);
auto fluidMaterial = std::make_shared<RenderMaterial>();
auto fluidMaterial = std::make_shared<RenderMaterial>();
fluidMaterial->setColor(Color::Blue);
fluidMaterial->setSphereGlyphSize(particleRadius);
fluidVisualModel->setRenderMaterial(fluidMaterial);
......@@ -149,10 +153,15 @@ std::shared_ptr<SPHObject> generateFluid(const std::shared_ptr<Scene>&scene, con
if (SCENE_ID == 2) // highly viscous fluid
{
sphParams->m_RatioKernelOverParticleRadius = 6.0;
sphParams->m_ViscosityFluid = 0.5;
sphParams->m_ViscosityFluid = 0.5;
sphParams->m_SurfaceTensionStiffness = 5.0;
}
if (SCENE_ID == 3) // bunny-shaped fluid
{
sphParams->m_BoundaryFriction = 0.3;
}
sphModel->configure(sphParams);
sphModel->setTimeStepSizeType(TimeSteppingType::realTime);
......
......@@ -26,75 +26,70 @@
namespace imstk
{
SPHCollisionHandling::SPHCollisionHandling(const CollisionHandling::Side &side,
const std::shared_ptr<CollisionData>& colData,
const std::shared_ptr<CollidingObject> &obj) :
SPHCollisionHandling::SPHCollisionHandling(const CollisionHandling::Side& side,
const std::shared_ptr<CollisionData>& colData,
const std::shared_ptr<CollidingObject>& obj) :
CollisionHandling(Type::SPH, side, colData), m_SPHObject(std::dynamic_pointer_cast<SPHObject>(obj))
{
LOG_IF(FATAL, (!m_SPHObject)) << "Invalid SPH object";
}
void SPHCollisionHandling::setBoundaryFriction(const Real friction)
{
m_BoundaryFriction = friction;
if (m_BoundaryFriction < 0.0 || m_BoundaryFriction > 1.0)
{
LOG(WARNING) << "Invalid frictionLength coefficient (it must be in [0, 1])";
if (m_BoundaryFriction < 0)
{
m_BoundaryFriction = 0;
}
else if (m_BoundaryFriction > static_cast<Real>(1.0))
{
m_BoundaryFriction = static_cast<Real>(1.0);
}
}
}
void SPHCollisionHandling::processCollisionData()
void
SPHCollisionHandling::processCollisionData()
{
const auto& SPHModel = m_SPHObject->getSPHModel();
#if defined(DEBUG) || defined(_DEBUG) || !defined(NDEBUG)
LOG_IF(FATAL, !SPHModel) << "SPH model was not initialized";
LOG_IF(FATAL, (!SPHModel)) << "SPH model was not initialized";
#endif
auto& state = SPHModel->getState();
const auto boundaryFriction = m_SPHObject->getSPHModel()->getParameters()->m_BoundaryFriction;
#if defined(DEBUG) || defined(_DEBUG) || !defined(NDEBUG)
LOG_IF(FATAL, (boundaryFriction<0.0 || boundaryFriction>1.0))
<< "Invalid boundary friction coefficient (value must be in [0, 1])";
#endif
auto& state = SPHModel->getState();
ParallelUtils::parallelFor(m_colData->MAColData.size(),
[&](const size_t idx)
{
const auto& cd = m_colData->MAColData[idx];
const auto pidx = cd.nodeId; // Fluid particle index
const auto pidx = cd.nodeId; // Fluid particle index
auto n = cd.penetrationVector; // This vector should point into solid object
// Correct particle position:
state.getPositions()[pidx] -= n; // This should not cause race condition, since pidx is distinct
// Correct particle position
state.getPositions()[pidx] -= n;
// Correct particle velocity: slip boundary condition with friction
auto& vel = state.getVelocities()[pidx];
const auto nLengthSqr = n.squaredNorm();
if (nLengthSqr < Real(1e-20)) // Normalize n
{
return; // Too little penetration: ignore
return; // Too little penetration: ignore
}
n /= std::sqrt(nLengthSqr);
const auto vn = vel.dot(n);
vel -= n * vn; // From now, vel is parallel with the solid surface
if (m_BoundaryFriction > Real(1e-20))
// Correct particle velocity: slip boundary condition with friction
const auto oldVel = state.getVelocities()[pidx];
const auto vn = oldVel.dot(n);
// If particle is escaping the boundary, ignore it
if (vn > 0)
{
const auto velLength = vel.norm();
const auto frictionLength = -vn * m_BoundaryFriction;
if (frictionLength < velLength && velLength > Real(1e-10))
{
vel -= (vel / velLength) * frictionLength; // Subtract a friction from velocity, which is proportional to the amount of penetration
}
else
Vec3r correctedVel = oldVel - vn * n; // From now, vel is parallel with the solid surface
if (boundaryFriction > Real(1e-20))
{
vel = Vec3r::Zero();
const auto velLength = correctedVel.norm();
const auto frictionLength = vn * boundaryFriction; // This is always positive
if (frictionLength < velLength && velLength > Real(1e-10))
{
correctedVel -= (correctedVel / velLength) * frictionLength; // Subtract a friction from velocity, which is proportional to the amount of penetration
}
else
{
correctedVel = Vec3r::Zero();
}
}
state.getVelocities()[pidx] = correctedVel;
}
});
}
......
......@@ -33,25 +33,12 @@ public:
SPHCollisionHandling() = delete;
virtual ~SPHCollisionHandling() override = default;
///
/// \brief Compute forces based on collision data
///
virtual void processCollisionData() override;
///
/// \brief Set the friction coefficient
///
void setBoundaryFriction(const Real friction);
///
/// \brief Get the friction coefficient
///
Real getBoundaryFriction() const { return m_BoundaryFriction; }
private:
std::shared_ptr<SPHObject> m_SPHObject;
Real m_BoundaryFriction = Real(0.1);
};
} // end namespace imstk
......@@ -31,7 +31,8 @@ SPHModelConfig::SPHModelConfig(const Real particleRadius)
initialize();
}
void SPHModelConfig::initialize()
void
SPHModelConfig::initialize()
{
LOG_IF(FATAL, (std::abs(m_ParticleRadius) < Real(1e-6))) << "Particle radius was not set properly";
......@@ -45,7 +46,8 @@ void SPHModelConfig::initialize()
m_KernelRadiusSqr = m_KernelRadius * m_KernelRadius;
}
bool SPHModel::initialize()
bool
SPHModel::initialize()
{
LOG_IF(FATAL, (!m_Geometry)) << "Model geometry is not yet set! Cannot initialize without model geometry.";
......@@ -74,7 +76,8 @@ bool SPHModel::initialize()
return true;
}
void SPHModel::advanceTimeStep()
void
SPHModel::advanceTimeStep()
{
findParticleNeighbors();
computeNeighborRelativePositions();
......@@ -89,12 +92,14 @@ void SPHModel::advanceTimeStep()
moveParticles(getTimeStep());
}
void SPHModel::computeTimeStepSize()
void
SPHModel::computeTimeStepSize()
{
m_dt = (this->m_timeStepSizeType == TimeSteppingType::fixed) ? m_DefaultDt : computeCFLTimeStepSize();
}
Real SPHModel::computeCFLTimeStepSize()
Real
SPHModel::computeCFLTimeStepSize()
{
auto maxVel = ParallelUtils::ParallelReduce::findMaxL2Norm(getState().getVelocities());
......@@ -114,7 +119,8 @@ Real SPHModel::computeCFLTimeStepSize()
return timestep;
}
void SPHModel::findParticleNeighbors()
void
SPHModel::findParticleNeighbors()
{
m_NeighborSearcher->getNeighbors(getState().getFluidNeighborLists(), getState().getPositions());
if (m_Parameters->m_bDensityWithBoundary) // if considering boundary particles for computing fluid density
......@@ -124,14 +130,15 @@ void SPHModel::findParticleNeighbors()
}
}
void SPHModel::computeNeighborRelativePositions()
void
SPHModel::computeNeighborRelativePositions()
{
auto computeRelativePositions = [&](const Vec3r& ppos, const std::vector<size_t>& neighborList,
const StdVectorOfVec3r& allPositions, std::vector<NeighborInfo>& neighborInfo) {
for (const size_t q : neighborList)
{
const Vec3r& qpos = allPositions[q];
const Vec3r r = ppos - qpos;
const Vec3r r = ppos - qpos;
neighborInfo.push_back({ r, m_Parameters->m_RestDensity });
}
};
......@@ -152,7 +159,8 @@ void SPHModel::computeNeighborRelativePositions()
});
}
void SPHModel::collectNeighborDensity()
void
SPHModel::collectNeighborDensity()
{
// after computing particle densities, cache them into neighborInfo variable, next to relative positions
// this is usefull because relative positions and densities are accessed together multiple times
......@@ -174,7 +182,8 @@ void SPHModel::collectNeighborDensity()
});
}
void SPHModel::computeDensity()
void
SPHModel::computeDensity()
{
ParallelUtils::parallelFor(getState().getNumParticles(),
[&](const size_t p) {
......@@ -194,7 +203,8 @@ void SPHModel::computeDensity()
});
}
void SPHModel::normalizeDensity()
void
SPHModel::normalizeDensity()
{
if (!m_Parameters->m_bNormalizeDensity)
{
......@@ -244,7 +254,8 @@ void SPHModel::normalizeDensity()
std::swap(getState().getDensities(), getState().getNormalizedDensities());
}
void SPHModel::computePressureAcceleration()
void
SPHModel::computePressureAcceleration()
{
auto particlePressure = [&](const Real density) {
const Real error = std::pow(density / m_Parameters->m_RestDensity, 7) - Real(1);
......@@ -281,7 +292,8 @@ void SPHModel::computePressureAcceleration()
});
}
void SPHModel::updateVelocity(Real timestep)
void
SPHModel::updateVelocity(Real timestep)
{
ParallelUtils::parallelFor(getState().getNumParticles(),
[&](const size_t p) {
......@@ -289,7 +301,8 @@ void SPHModel::updateVelocity(Real timestep)
});
}
void SPHModel::computeViscosity()
void
SPHModel::computeViscosity()
{
ParallelUtils::parallelFor(getState().getNumParticles(),
[&](const size_t p) {
......@@ -338,7 +351,8 @@ void SPHModel::computeViscosity()
}
// Compute surface tension Akinci et at. 2013 model (Versatile Surface Tension and Adhesion for SPH Fluids)
void SPHModel::computeSurfaceTension()
void
SPHModel::computeSurfaceTension()
{
// Firstly compute surface normal for all particles
ParallelUtils::parallelFor(getState().getNumParticles(),
......@@ -408,7 +422,8 @@ void SPHModel::computeSurfaceTension()
});
}
void SPHModel::moveParticles(Real timestep)
void
SPHModel::moveParticles(Real timestep)
{
ParallelUtils::parallelFor(getState().getNumParticles(),
[&](const size_t p) {
......
......@@ -56,9 +56,9 @@ public:
////////////////////////////////////////////////////////////////////////////////
// material parameters
Real m_RestDensity = Real(1000.0);
Real m_RestDensitySqr = Real(1000000.0);
Real m_RestDensityInv = Real(1.0 / 1000.0);
Real m_RestDensity = Real(1000.0);
Real m_RestDensitySqr = Real(1000000.0);
Real m_RestDensityInv = Real(1.0 / 1000.0);
Real m_ParticleMass = Real(1);
Real m_ParticleMassScale = Real(0.95); // scale particle mass to a smaller value to maintain stability
......@@ -76,6 +76,7 @@ public:
Real m_ViscosityFluid = Real(1e-2);
Real m_ViscosityBoundary = Real(1e-5);
Real m_SurfaceTensionStiffness = Real(1);
Real m_BoundaryFriction = Real(0.1);
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
......@@ -242,10 +243,10 @@ private:
void moveParticles(const Real timestep);
std::shared_ptr<PointSet> m_Geometry;
SPHSimulationState m_SimulationState;
SPHSimulationState m_SimulationState;
Real m_dt; ///> time step size
Real m_DefaultDt; ///> default time step size
Real m_dt; ///> time step size
Real m_DefaultDt; ///> default time step size
SPHSimulationKernels m_Kernels; // must be initialized during model initialization
std::shared_ptr<SPHModelConfig> m_Parameters; // must be set before simulation
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment