Commit 0de68a60 authored by Sreekanth Arikatla's avatar Sreekanth Arikatla
Browse files

Merge branch 'Parallel-SPH' into 'master'

Parallel SPH

See merge request !318
parents cb157931 12a93d9d
......@@ -4,7 +4,7 @@
include(imstkAddExternalProject)
imstk_add_external_project(tbb
GIT_REPOSITORY https://gitlab.kitware.com/iMSTK/TBB-CMake.git
GIT_TAG tbb44u4
GIT_TAG b066defc0229a1e92d7a200eb3fe0f7e35945d95
INSTALL_COMMAND ${SKIP_STEP_COMMAND}
CMAKE_CACHE_ARGS
-DTBB_BUILD_TESTS:BOOL=OFF
......
......@@ -21,12 +21,16 @@ project(Example-SPHFluid)
#-----------------------------------------------------------------------------
# Create executable
#-----------------------------------------------------------------------------
add_executable(${PROJECT_NAME} Main.cpp Fluid.cpp Solid.cpp Bunny.cpp)
add_executable(${PROJECT_NAME}-BallDrop SPHFluid-BallDrop.cpp Fluid.cpp Solid.cpp Bunny.cpp)
add_executable(${PROJECT_NAME}-HighViscousity SPHFluid-HighViscousity.cpp Fluid.cpp Solid.cpp Bunny.cpp)
add_executable(${PROJECT_NAME}-BunnyShape SPHFluid-BunnyShape.cpp Fluid.cpp Solid.cpp Bunny.cpp)
#-----------------------------------------------------------------------------
# Link libraries to executable
#-----------------------------------------------------------------------------
target_link_libraries(${PROJECT_NAME} SimulationManager)
target_link_libraries(${PROJECT_NAME}-BallDrop SimulationManager)
target_link_libraries(${PROJECT_NAME}-HighViscousity SimulationManager)
target_link_libraries(${PROJECT_NAME}-BunnyShape SimulationManager)
#-----------------------------------------------------------------------------
# Add shaders
......
/*=========================================================================
Library: iMSTK
Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
& Imaging in Medicine, Rensselaer Polytechnic Institute.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0.txt
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=========================================================================*/
// Generate a fluid simulation example in with a sphere-shape fluid dropping onto the ground
#define SCENE_ID 1
#include "SPHFluidExample.hpp"
\ No newline at end of file
/*=========================================================================
Library: iMSTK
Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
& Imaging in Medicine, Rensselaer Polytechnic Institute.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0.txt
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=========================================================================*/
// Generate a fluid simulation example in with a bunny-shape fluid dropping onto a box
#define SCENE_ID 3
#include "SPHFluidExample.hpp"
\ No newline at end of file
/*=========================================================================
Library: iMSTK
Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
& Imaging in Medicine, Rensselaer Polytechnic Institute.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0.txt
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=========================================================================*/
// Generate a fluid simulation example in with a box-shape of highly viscous fluid dropping onto the ground
#define SCENE_ID 2
#include "SPHFluidExample.hpp"
\ No newline at end of file
......@@ -19,9 +19,13 @@
=========================================================================*/
#include <string>
#include "imstkSimulationManager.h"
#include "imstkSPHObject.h"
#include "imstkAPIUtilities.h"
#include "imstkPlane.h"
#include "imstkSphere.h"
using namespace imstk;
......@@ -29,60 +33,52 @@ std::shared_ptr<SPHObject> generateFluid(const std::shared_ptr<Scene>&scene, int
std::vector<std::shared_ptr<CollidingObject>> generateSolids(const std::shared_ptr<Scene>& scene, int sceneIdx);
///
/// \brief This example demonstrates the fluid simulation using SPH
/// \brief Usage: ./SPHFluid [scene=<scene_id>] [threads=<num_threads>] [radius=<particle_radius>], where scene_id is in [1..3]
/// \brief Example: ./SPHFluid scene=1 threads=8 radius=0.01
/// \brief Usage: ./SPHFluid [threads=<num_threads>] [radius=<particle_radius>]
/// \brief Example: ./SPHFluid threads=8 radius=0.01
///
int main(int argc, char* argv[])
{
// SimulationManager must be created first
auto sdk = std::make_shared<SimulationManager>(0);
auto sdk = std::make_shared<SimulationManager>(0);
int threads = -1;
int sceneIdx = 1;
double particleRadius = 0.1;
// Parse command line arguments
for (int i = 1; i < argc; ++i)
{
auto param = std::string(argv[i]);
if (param.find("threads") != std::string::npos &&
if (param.find("threads") == 0 &&
param.find_first_of("=") != std::string::npos)
{
threads = std::stoi(param.substr(param.find_first_of("=") + 1));
}
else if (param.find("scene") != std::string::npos &&
param.find_first_of("=") != std::string::npos)
{
sceneIdx = std::stoi(param.substr(param.find_first_of("=") + 1));
if (sceneIdx < 1)
{
sceneIdx = 1;
}
else if (sceneIdx > 3)
{
sceneIdx = 3;
}
LOG(INFO) << "Scene ID: " << sceneIdx;
}
else if (param.find("radius") != std::string::npos &&
else if (param.find("radius") == 0 &&
param.find_first_of("=") != std::string::npos)
{
particleRadius = std::stod(param.substr(param.find_first_of("=") + 1));
LOG(INFO) << "Particle radius: " << particleRadius;
}
else
{
LOG(FATAL) << "Invalid argument";
}
}
// Particle in this scene is pre-generated using particle radius 0.08
if (sceneIdx == 3)
// Override particle radius for scene3 because particles in this scene were pre-generated using particle radius 0.08
if (SCENE_ID == 3)
{
particleRadius = 0.08;
}
// Set thread pool size (nthreads <= 0 means using all logical cores)
sdk->setThreadPoolSize(threads);
auto scene = sdk->createNewScene("SPH Fluid");
// Generate fluid and solid objects
auto fluidObj = generateFluid(scene, sceneIdx, particleRadius);
auto solids = generateSolids(scene, sceneIdx);
auto fluidObj = generateFluid(scene, SCENE_ID, particleRadius);
auto solids = generateSolids(scene, SCENE_ID);
// Collision between fluid and solid objects
auto colGraph = scene->getCollisionGraph();
......
/*=========================================================================
Library: iMSTK
Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
& Imaging in Medicine, Rensselaer Polytechnic Institute.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0.txt
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=========================================================================*/
#pragma once
#include <atomic>
#include "imstkMath.h"
namespace imstk
{
namespace ParallelUtils
{
///
/// \brief Perform an atomic operation: target = f(target, operand)
///
template<class T, class Function>
void atomicOp(T& target, const T operand, Function&& f)
{
std::atomic<T>& tgt = *(static_cast<std::atomic<T>*>(&target));
T cur_val = target;
T new_val;
do
{
new_val = f(cur_val, operand);
}
while (!tgt.compare_exchange_weak(cur_val, new_val));
}
///
/// \brief Atomic addition for scalar numbers: target = target + operand
///
template<class T>
void atomicAdd(T& target, const T operand)
{
atomicOp(target, operand, [](T a, T b) { return a + b; });
}
///
/// \brief Atomic subtraction for scalar numbers: target = target - operand
///
template<class T>
void atomicSubtract(T& target, const T operand)
{
atomicOp(target, operand, [](T a, T b) { return a - b; });
}
///
/// \brief Atomic multiplication for scalar numbers: target = target * operand
///
template<class T>
void atomicMultiply(T& target, const T operand)
{
atomicOp(target, operand, [](T a, T b) { return a * b; });
}
///
/// \brief Atomic division for scalar numbers: target = target / operand
///
template<class T>
void atomicDivide(T& target, const T operand)
{
atomicOp(target, operand, [](T a, T b) { return a / b; });
}
///
/// \brief Atomic addition for two vectors: target = target + operand
///
template<class T, int N>
void atomicAdd(Eigen::Matrix<T, N, 1>& target, const Eigen::Matrix<T, N, 1>& operand)
{
for (int i = 0; i < N; ++i)
{
atomicAdd(target[i], operand[i]);
}
}
///
/// \brief Atomic subtraction for two vectors: target = target - operand
///
template<class T, int N>
void atomicSubtract(Eigen::Matrix<T, N, 1>& target, const Eigen::Matrix<T, N, 1>& operand)
{
for (int i = 0; i < N; ++i)
{
atomicSubtract(target[i], operand[i]);
}
}
///
/// \brief Atomic multiplication for a vector and a scalar number: target = target * operand
///
template<class T, int N>
void atomicMultiply(Eigen::Matrix<T, N, 1>& target, const T operand)
{
for (int i = 0; i < N; ++i)
{
atomicMultiply(target[i], operand);
}
}
///
/// \brief Atomic division for a vector and a scalar number: target = target / operand
///
template<class T, int N>
void atomicDivide(Eigen::Matrix<T, N, 1>& target, const T operand)
{
for (int i = 0; i < N; ++i)
{
atomicDivide(target[i], operand);
}
}
}// end namespace namespace ParallelUtils
}// end namespace imstk
/*=========================================================================
Library: iMSTK
Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
& Imaging in Medicine, Rensselaer Polytechnic Institute.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0.txt
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=========================================================================*/
#pragma once
#include <tbb/tbb.h>
namespace imstk
{
namespace ParallelUtils
{
///
/// \brief Execute a function in parallel over a range [beginIdx, endIdx) of indices
///
template<class IndexType, class Function>
void parallelFor(const IndexType beginIdx, const IndexType endIdx, Function&& function)
{
tbb::parallel_for(tbb::blocked_range<IndexType>(beginIdx, endIdx),
[&](const tbb::blocked_range<IndexType>& r) {
for (IndexType i = r.begin(), iEnd = r.end(); i < iEnd; ++i)
{
function(i);
}
});
}
///
/// \brief Execute a function in parallel over a range [0, endIdx) of indices
///
template<class IndexType, class Function>
void parallelFor(const IndexType endIdx, Function&& function)
{
parallelFor(IndexType(0), endIdx, std::forward<Function>(function));
}
///
/// \brief Execute a 2D function in parallel over a range of indices in the x dimension,
/// indices in the y dimension are scanned sequentially
///
template<class IndexType, class Function>
void parallelFor2Dx(const IndexType beginX, const IndexType endX,
const IndexType beginY, const IndexType endY,
Function&& function)
{
parallelFor(beginX, endX,
[&](IndexType i) {
for (IndexType j = beginY; j < endY; ++j)
{
function(i, j);
}
});
}
///
/// \brief Execute a 2D function in parallel over a range of indices in the y dimension,
/// indices in the x dimension are scanned sequentially
///
template<class IndexType, class Function>
void parallelFor2Dy(const IndexType beginX, const IndexType endX,
const IndexType beginY, const IndexType endY,
Function&& function)
{
parallelFor(beginY, endY,
[&](IndexType j) {
for (IndexType i = beginX; i < endX; ++i)
{
function(i, j);
}
});
}
///
/// \brief Execute a 3D function in parallel over a range of indices in the x dimension,
/// indices in the y and z dimensions are scanned sequentially
///
template<class IndexType, class Function>
void parallelFor3Dx(const IndexType beginX, const IndexType endX,
const IndexType beginY, const IndexType endY,
const IndexType beginZ, const IndexType endZ,
Function&& function)
{
parallelFor(beginX, endX,
[&](IndexType i) {
for (IndexType j = beginY; j < endY; ++j)
{
for (IndexType k = beginZ; k < endZ; ++k)
{
function(i, j, k);
}
}
});
}
///
/// \brief Execute a 3D function in parallel over a range of indices in the y dimension,
/// indices in the x and z dimensions are scanned sequentially
///
template<class IndexType, class Function>
void parallelFor3Dy(const IndexType beginX, const IndexType endX,
const IndexType beginY, const IndexType endY,
const IndexType beginZ, const IndexType endZ,
Function&& function)
{
parallelFor(beginY, endY,
[&](IndexType j) {
for (IndexType i = beginX; i < endX; ++i)
{
for (IndexType k = beginZ; k < endZ; ++k)
{
function(i, j, k);
}
}
});
}
///
/// \brief Execute a 3D function in parallel over a range of indices in the z dimension,
/// indices in the x and y dimensions are scanned sequentially
///
template<class IndexType, class Function>
void parallelFor3Dz(const IndexType beginX, const IndexType endX,
const IndexType beginY, const IndexType endY,
const IndexType beginZ, const IndexType endZ,
Function&& function)
{
parallelFor(beginX, endX,
[&](IndexType i) {
for (IndexType j = beginY; j < endY; ++j)
{
for (IndexType k = beginZ; k < endZ; ++k)
{
function(i, j, k);
}
}
});
}
} // end namespace ParallelUtils
} // end namespace imstk
/*=========================================================================
Library: iMSTK
Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
& Imaging in Medicine, Rensselaer Polytechnic Institute.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0.txt
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=========================================================================*/
#pragma once
#include "imstkMath.h"
#include <tbb/tbb.h>
#undef min
#undef max
namespace imstk
{
namespace ParallelUtils
{
///
/// \brief The ParallelReduce class
/// \brief A class for executing reduce operations in parallel
///
class ParallelReduce
{
///
/// \brief Private helper class, providing operator() using in std::parallel_reduce
/// for finding max L2 norm of an array of Vec3r
///
class MaxL2NormFunctor
{
public:
MaxL2NormFunctor(const StdVectorOfVec3r& data) : m_Data(data) {}
MaxL2NormFunctor(MaxL2NormFunctor& pObj, tbb::split) : m_Data(pObj.m_Data) {}
void operator()(const tbb::blocked_range<size_t>& r)
{
for (size_t i = r.begin(); i != r.end(); ++i)
{
Real mag2 = m_Data[i].squaredNorm();
m_Result = m_Result > mag2 ? m_Result : mag2;
}
}
void join(MaxL2NormFunctor& pObj) { m_Result = m_Result > pObj.m_Result ? m_Result : pObj.m_Result; }
Real getResult() const { return std::sqrt(m_Result); }
private:
Real m_Result = 0;
const StdVectorOfVec3r& m_Data;
};
///
/// \brief Private helper class, providing operator() using in std::parallel_reduce
/// for finding axis-aligned bounding box of a point set
///
class AABBFunctor
{
public:
AABBFunctor(const StdVectorOfVec3r& data) : m_Data(data) { if (data.size() > 0) { m_UpperCorner = data[0]; } }
AABBFunctor(AABBFunctor& pObj, tbb::split) : m_Data(pObj.m_Data) {}
void operator()(const tbb::blocked_range<size_t>& r)
{
for (size_t i = r.begin(); i != r.end(); ++i)
{
const auto& vec = m_Data[i];
for (int j = 0; j < 3; ++j)
{
m_LowerCorner[j] = (m_LowerCorner[j] < vec[j]) ? m_LowerCorner[j] : vec[j];
m_UpperCorner[j] = (m_UpperCorner[j] > vec[j]) ? m_UpperCorner[j] : vec[j];
}
}
}
void join(AABBFunctor& pObj)
{
for (int j = 0; j < 3; ++j)
{
m_LowerCorner[j] = (m_LowerCorner[j] < pObj.m_LowerCorner[j]) ? m_LowerCorner[j] : pObj.m_LowerCorner[j];
m_UpperCorner[j] = (m_UpperCorner[j] > pObj.m_UpperCorner[j]) ? m_UpperCorner[j] : pObj.m_UpperCorner[j];
}
}
const Vec3r& getLowerCorner() const { return m_LowerCorner; }
const Vec3r& getUpperCorner() const { return m_UpperCorner; }
private:
Vec3r m_LowerCorner = Vec3r(std::numeric_limits<Real>::max(),
std::numeric_limits<Real>::max(),
std::numeric_limits<Real>::max());
Vec3r m_UpperCorner = Vec3r(-std::numeric_limits<Real>::max(),
-std::numeric_limits<Real>::max(),
-std::numeric_limits<Real>::max());
const StdVectorOfVec3r& m_Data;
};
public:
///
/// \brief Find the maximum value of L2 norm from the input data array
///
static Real findMaxL2Norm(const StdVectorOfVec3r& data)
{
MaxL2NormFunctor pObj(data);