Skip to content
Snippets Groups Projects
Commit 33d21edf authored by Christoph Schwarzmeier's avatar Christoph Schwarzmeier
Browse files

Fix possible integer overflow in center of mass computation in FSLBM showcases

parent 96634fea
No related merge requests found
Hide whitespace changes
Inline Side-by-side
apps/showcases/FreeSurface/RisingBubble.cpp
+ 19
− 21
View file @ 33d21edf
...@@ -90,15 +90,17 @@ class CenterOfMassComputer ...@@ -90,15 +90,17 @@ class CenterOfMassComputer
++executionCounter_; ++executionCounter_;
// only evaluate in given frequencies // only evaluate in given intervals
if (executionCounter_ % frequency_ != uint_c(0)) { return; } if (executionCounter_ % frequency_ != uint_c(0)) { return; }
const BlockDataID flagFieldID = freeSurfaceBoundaryHandling->getFlagFieldID(); const BlockDataID flagFieldID = freeSurfaceBoundaryHandling->getFlagFieldID();
const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo(); const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo();
*centerOfMass_ = Vector3< real_t >(real_c(0)); *centerOfMass_ = Vector3< real_t >(real_c(0));
Cell cellSum = Cell(cell_idx_c(0), cell_idx_c(0), cell_idx_c(0));
uint_t cellCount = uint_c(0); // use real_t to avoid integer overflow in sum below when using high grid resolutions
Vector3< real_t > cellSumVector = Vector3< real_t >(real_c(0));
real_t cellCount = real_c(0);
for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt) for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
{ {
...@@ -110,21 +112,18 @@ class CenterOfMassComputer ...@@ -110,21 +112,18 @@ class CenterOfMassComputer
Cell globalCell = flagFieldIt.cell(); Cell globalCell = flagFieldIt.cell();
// transform local cell to global coordinates // transform local cell to global coordinates
blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt); blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt);
cellSum += globalCell; cellSumVector += Vector3< real_t >(real_c(globalCell[0]), real_c(globalCell[1]), real_c(globalCell[2]));
++cellCount; cellCount += real_c(1);
} }
}) // WALBERLA_FOR_ALL_CELLS }) // WALBERLA_FOR_ALL_CELLS
} }
mpi::allReduceInplace< uint_t >(cellCount, mpi::SUM); mpi::allReduceInplace< real_t >(cellCount, mpi::SUM);
mpi::allReduceInplace< real_t >(cellSumVector, mpi::SUM);
Vector3< cell_idx_t > cellSumVector =
Vector3< cell_idx_t >(cellSum[0], cellSum[1], cellSum[2]); // use Vector3 to limit calls to MPI-reduce to one
mpi::allReduceInplace< cell_idx_t >(cellSumVector, mpi::SUM);
(*centerOfMass_)[0] = real_c(cellSumVector[0]) / real_c(cellCount); (*centerOfMass_)[0] = cellSumVector[0] / cellCount;
(*centerOfMass_)[1] = real_c(cellSumVector[1]) / real_c(cellCount); (*centerOfMass_)[1] = cellSumVector[1] / cellCount;
(*centerOfMass_)[2] = real_c(cellSumVector[2]) / real_c(cellCount); (*centerOfMass_)[2] = cellSumVector[2] / cellCount;
} }
private: private:
...@@ -412,13 +411,12 @@ int main(int argc, char** argv) ...@@ -412,13 +411,12 @@ int main(int argc, char** argv)
{ {
timeloop.singleStep(timingPool, true); timeloop.singleStep(timingPool, true);
// only evaluate if center of mass has moved "significantly" if (t % evaluationFrequency == uint_c(0) && t > uint_c(0))
if ((*centerOfMass)[2] > (centerOfMassOld[2] + real_c(1)) && t > uint_c(0))
{ {
WALBERLA_ROOT_SECTION() WALBERLA_ROOT_SECTION()
{ {
real_t riseVelocity = ((*centerOfMass)[2] - centerOfMassOld[2]) / real_c(t - timestepOld); const real_t riseVelocity = ((*centerOfMass)[2] - centerOfMassOld[2]) / real_c(t - timestepOld);
const real_t dragForce = real_c(4) / real_c(3) * gravitationalAccelerationZ * real_c(bubbleDiameter) / const real_t dragForce = real_c(4) / real_c(3) * gravitationalAccelerationZ * real_c(bubbleDiameter) /
(riseVelocity * riseVelocity); (riseVelocity * riseVelocity);
WALBERLA_LOG_DEVEL("time step = " << t << "\n\t\tcenterOfMass = " << *centerOfMass WALBERLA_LOG_DEVEL("time step = " << t << "\n\t\tcenterOfMass = " << *centerOfMass
...@@ -429,10 +427,10 @@ int main(int argc, char** argv) ...@@ -429,10 +427,10 @@ int main(int argc, char** argv)
const std::vector< real_t > resultVector{ (*centerOfMass)[2], riseVelocity, dragForce }; const std::vector< real_t > resultVector{ (*centerOfMass)[2], riseVelocity, dragForce };
writeVectorToFile(resultVector, t, filename); writeVectorToFile(resultVector, t, filename);
}
timestepOld = t; timestepOld = t;
centerOfMassOld = *centerOfMass; centerOfMassOld = *centerOfMass;
}
} }
// stop simulation before bubble hits the top wall // stop simulation before bubble hits the top wall
......
apps/showcases/FreeSurface/TaylorBubble.cpp
+ 22
− 24
View file @ 33d21edf
...@@ -76,9 +76,9 @@ class CenterOfMassComputer ...@@ -76,9 +76,9 @@ class CenterOfMassComputer
public: public:
CenterOfMassComputer(const std::weak_ptr< const StructuredBlockForest >& blockForest, CenterOfMassComputer(const std::weak_ptr< const StructuredBlockForest >& blockForest,
const std::weak_ptr< const FreeSurfaceBoundaryHandling_T >& freeSurfaceBoundaryHandling, const std::weak_ptr< const FreeSurfaceBoundaryHandling_T >& freeSurfaceBoundaryHandling,
uint_t interval, const std::shared_ptr< Vector3< real_t > >& centerOfMass) uint_t frequency, const std::shared_ptr< Vector3< real_t > >& centerOfMass)
: blockForest_(blockForest), freeSurfaceBoundaryHandling_(freeSurfaceBoundaryHandling), : blockForest_(blockForest), freeSurfaceBoundaryHandling_(freeSurfaceBoundaryHandling),
centerOfMass_(centerOfMass), interval_(interval), executionCounter_(uint_c(0)) centerOfMass_(centerOfMass), frequency_(frequency), executionCounter_(uint_c(0))
{} {}
void operator()() void operator()()
...@@ -92,14 +92,16 @@ class CenterOfMassComputer ...@@ -92,14 +92,16 @@ class CenterOfMassComputer
++executionCounter_; ++executionCounter_;
// only evaluate in given intervals // only evaluate in given intervals
if (executionCounter_ % interval_ != uint_c(0)) { return; } if (executionCounter_ % frequency_ != uint_c(0)) { return; }
const BlockDataID flagFieldID = freeSurfaceBoundaryHandling->getFlagFieldID(); const BlockDataID flagFieldID = freeSurfaceBoundaryHandling->getFlagFieldID();
const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo(); const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo();
*centerOfMass_ = Vector3< real_t >(real_c(0)); *centerOfMass_ = Vector3< real_t >(real_c(0));
Cell cellSum = Cell(cell_idx_c(0), cell_idx_c(0), cell_idx_c(0));
uint_t cellCount = uint_c(0); // use real_t to avoid integer overflow in sum below when using high grid resolutions
Vector3< real_t > cellSumVector = Vector3< real_t >(real_c(0));
real_t cellCount = real_c(0);
for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt) for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
{ {
...@@ -111,21 +113,18 @@ class CenterOfMassComputer ...@@ -111,21 +113,18 @@ class CenterOfMassComputer
Cell globalCell = flagFieldIt.cell(); Cell globalCell = flagFieldIt.cell();
// transform local cell to global coordinates // transform local cell to global coordinates
blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt); blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt);
cellSum += globalCell; cellSumVector += Vector3< real_t >(real_c(globalCell[0]), real_c(globalCell[1]), real_c(globalCell[2]));
++cellCount; cellCount += real_c(1);
} }
}) // WALBERLA_FOR_ALL_CELLS }) // WALBERLA_FOR_ALL_CELLS
} }
mpi::allReduceInplace< uint_t >(cellCount, mpi::SUM); mpi::allReduceInplace< real_t >(cellCount, mpi::SUM);
mpi::allReduceInplace< real_t >(cellSumVector, mpi::SUM);
Vector3< cell_idx_t > cellSumVector =
Vector3< cell_idx_t >(cellSum[0], cellSum[1], cellSum[2]); // use Vector3 to limit calls to MPI-reduce to one
mpi::allReduceInplace< cell_idx_t >(cellSumVector, mpi::SUM);
(*centerOfMass_)[0] = real_c(cellSumVector[0]) / real_c(cellCount); (*centerOfMass_)[0] = cellSumVector[0] / cellCount;
(*centerOfMass_)[1] = real_c(cellSumVector[1]) / real_c(cellCount); (*centerOfMass_)[1] = cellSumVector[1] / cellCount;
(*centerOfMass_)[2] = real_c(cellSumVector[2]) / real_c(cellCount); (*centerOfMass_)[2] = cellSumVector[2] / cellCount;
} }
private: private:
...@@ -133,7 +132,7 @@ class CenterOfMassComputer ...@@ -133,7 +132,7 @@ class CenterOfMassComputer
std::weak_ptr< const FreeSurfaceBoundaryHandling_T > freeSurfaceBoundaryHandling_; std::weak_ptr< const FreeSurfaceBoundaryHandling_T > freeSurfaceBoundaryHandling_;
std::shared_ptr< Vector3< real_t > > centerOfMass_; std::shared_ptr< Vector3< real_t > > centerOfMass_;
uint_t interval_; uint_t frequency_;
uint_t executionCounter_; uint_t executionCounter_;
}; // class CenterOfMassComputer }; // class CenterOfMassComputer
...@@ -441,13 +440,12 @@ int main(int argc, char** argv) ...@@ -441,13 +440,12 @@ int main(int argc, char** argv)
{ {
timeloop.singleStep(timingPool, true); timeloop.singleStep(timingPool, true);
// only evaluate if center of mass has moved "significantly" if (t % evaluationFrequency == uint_c(0) && t > uint_c(0))
if ((*centerOfMass)[2] > (centerOfMassOld[2] + real_c(1)) && t > uint_c(0))
{ {
WALBERLA_ROOT_SECTION() WALBERLA_ROOT_SECTION()
{ {
real_t riseVelocity = ((*centerOfMass)[2] - centerOfMassOld[2]) / real_c(t - timestepOld); const real_t riseVelocity = ((*centerOfMass)[2] - centerOfMassOld[2]) / real_c(t - timestepOld);
const real_t dragForce = real_c(4) / real_c(3) * gravitationalAccelerationZ * real_c(bubbleDiameter) / const real_t dragForce = real_c(4) / real_c(3) * gravitationalAccelerationZ * real_c(bubbleDiameter) /
(riseVelocity * riseVelocity); (riseVelocity * riseVelocity);
WALBERLA_LOG_DEVEL("time step = " << t << "\n\t\tcenterOfMass = " << *centerOfMass WALBERLA_LOG_DEVEL("time step = " << t << "\n\t\tcenterOfMass = " << *centerOfMass
...@@ -458,12 +456,12 @@ int main(int argc, char** argv) ...@@ -458,12 +456,12 @@ int main(int argc, char** argv)
const std::vector< real_t > resultVector{ (*centerOfMass)[2], riseVelocity, dragForce }; const std::vector< real_t > resultVector{ (*centerOfMass)[2], riseVelocity, dragForce };
writeVectorToFile(resultVector, t, filename); writeVectorToFile(resultVector, t, filename);
timestepOld = t;
centerOfMassOld = *centerOfMass;
} }
} }
timestepOld = t;
centerOfMassOld = *centerOfMass;
// stop simulation before bubble hits the top wall // stop simulation before bubble hits the top wall
if ((*centerOfMass)[2] > stoppingHeight) { break; } if ((*centerOfMass)[2] > stoppingHeight) { break; }
......
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