From 33d21edf1fde3f1414421bcbc391d48f3449d781 Mon Sep 17 00:00:00 2001
From: Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
Date: Tue, 14 Mar 2023 13:24:02 +0100
Subject: [PATCH] Fix possible integer overflow in center of mass computation
in FSLBM showcases
---
apps/showcases/FreeSurface/RisingBubble.cpp | 40 +++++++++---------
apps/showcases/FreeSurface/TaylorBubble.cpp | 46 ++++++++++-----------
2 files changed, 41 insertions(+), 45 deletions(-)
diff --git a/apps/showcases/FreeSurface/RisingBubble.cpp b/apps/showcases/FreeSurface/RisingBubble.cpp
index 6699fdf1e..27e714c31 100644
--- a/apps/showcases/FreeSurface/RisingBubble.cpp
+++ b/apps/showcases/FreeSurface/RisingBubble.cpp
@@ -90,15 +90,17 @@ class CenterOfMassComputer
++executionCounter_;
- // only evaluate in given frequencies
+ // only evaluate in given intervals
if (executionCounter_ % frequency_ != uint_c(0)) { return; }
const BlockDataID flagFieldID = freeSurfaceBoundaryHandling->getFlagFieldID();
const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo();
- *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);
+ *centerOfMass_ = Vector3< real_t >(real_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)
{
@@ -110,21 +112,18 @@ class CenterOfMassComputer
Cell globalCell = flagFieldIt.cell();
// transform local cell to global coordinates
blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt);
- cellSum += globalCell;
- ++cellCount;
+ cellSumVector += Vector3< real_t >(real_c(globalCell[0]), real_c(globalCell[1]), real_c(globalCell[2]));
+ cellCount += real_c(1);
}
}) // WALBERLA_FOR_ALL_CELLS
}
- mpi::allReduceInplace< uint_t >(cellCount, 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);
+ mpi::allReduceInplace< real_t >(cellCount, mpi::SUM);
+ mpi::allReduceInplace< real_t >(cellSumVector, mpi::SUM);
- (*centerOfMass_)[0] = real_c(cellSumVector[0]) / real_c(cellCount);
- (*centerOfMass_)[1] = real_c(cellSumVector[1]) / real_c(cellCount);
- (*centerOfMass_)[2] = real_c(cellSumVector[2]) / real_c(cellCount);
+ (*centerOfMass_)[0] = cellSumVector[0] / cellCount;
+ (*centerOfMass_)[1] = cellSumVector[1] / cellCount;
+ (*centerOfMass_)[2] = cellSumVector[2] / cellCount;
}
private:
@@ -412,13 +411,12 @@ int main(int argc, char** argv)
{
timeloop.singleStep(timingPool, true);
- // only evaluate if center of mass has moved "significantly"
- if ((*centerOfMass)[2] > (centerOfMassOld[2] + real_c(1)) && t > uint_c(0))
+ if (t % evaluationFrequency == uint_c(0) && t > uint_c(0))
{
WALBERLA_ROOT_SECTION()
{
- 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 riseVelocity = ((*centerOfMass)[2] - centerOfMassOld[2]) / real_c(t - timestepOld);
+ const real_t dragForce = real_c(4) / real_c(3) * gravitationalAccelerationZ * real_c(bubbleDiameter) /
(riseVelocity * riseVelocity);
WALBERLA_LOG_DEVEL("time step = " << t << "\n\t\tcenterOfMass = " << *centerOfMass
@@ -429,10 +427,10 @@ int main(int argc, char** argv)
const std::vector< real_t > resultVector{ (*centerOfMass)[2], riseVelocity, dragForce };
writeVectorToFile(resultVector, t, filename);
- }
- timestepOld = t;
- centerOfMassOld = *centerOfMass;
+ timestepOld = t;
+ centerOfMassOld = *centerOfMass;
+ }
}
// stop simulation before bubble hits the top wall
diff --git a/apps/showcases/FreeSurface/TaylorBubble.cpp b/apps/showcases/FreeSurface/TaylorBubble.cpp
index ce2a45302..9685f8fad 100644
--- a/apps/showcases/FreeSurface/TaylorBubble.cpp
+++ b/apps/showcases/FreeSurface/TaylorBubble.cpp
@@ -76,9 +76,9 @@ class CenterOfMassComputer
public:
CenterOfMassComputer(const std::weak_ptr< const StructuredBlockForest >& blockForest,
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),
- centerOfMass_(centerOfMass), interval_(interval), executionCounter_(uint_c(0))
+ centerOfMass_(centerOfMass), frequency_(frequency), executionCounter_(uint_c(0))
{}
void operator()()
@@ -92,14 +92,16 @@ class CenterOfMassComputer
++executionCounter_;
// 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 typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo();
- *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);
+ *centerOfMass_ = Vector3< real_t >(real_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)
{
@@ -111,21 +113,18 @@ class CenterOfMassComputer
Cell globalCell = flagFieldIt.cell();
// transform local cell to global coordinates
blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt);
- cellSum += globalCell;
- ++cellCount;
+ cellSumVector += Vector3< real_t >(real_c(globalCell[0]), real_c(globalCell[1]), real_c(globalCell[2]));
+ cellCount += real_c(1);
}
}) // WALBERLA_FOR_ALL_CELLS
}
- mpi::allReduceInplace< uint_t >(cellCount, 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);
+ mpi::allReduceInplace< real_t >(cellCount, mpi::SUM);
+ mpi::allReduceInplace< real_t >(cellSumVector, mpi::SUM);
- (*centerOfMass_)[0] = real_c(cellSumVector[0]) / real_c(cellCount);
- (*centerOfMass_)[1] = real_c(cellSumVector[1]) / real_c(cellCount);
- (*centerOfMass_)[2] = real_c(cellSumVector[2]) / real_c(cellCount);
+ (*centerOfMass_)[0] = cellSumVector[0] / cellCount;
+ (*centerOfMass_)[1] = cellSumVector[1] / cellCount;
+ (*centerOfMass_)[2] = cellSumVector[2] / cellCount;
}
private:
@@ -133,7 +132,7 @@ class CenterOfMassComputer
std::weak_ptr< const FreeSurfaceBoundaryHandling_T > freeSurfaceBoundaryHandling_;
std::shared_ptr< Vector3< real_t > > centerOfMass_;
- uint_t interval_;
+ uint_t frequency_;
uint_t executionCounter_;
}; // class CenterOfMassComputer
@@ -441,13 +440,12 @@ int main(int argc, char** argv)
{
timeloop.singleStep(timingPool, true);
- // only evaluate if center of mass has moved "significantly"
- if ((*centerOfMass)[2] > (centerOfMassOld[2] + real_c(1)) && t > uint_c(0))
+ if (t % evaluationFrequency == uint_c(0) && t > uint_c(0))
{
WALBERLA_ROOT_SECTION()
{
- 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 riseVelocity = ((*centerOfMass)[2] - centerOfMassOld[2]) / real_c(t - timestepOld);
+ const real_t dragForce = real_c(4) / real_c(3) * gravitationalAccelerationZ * real_c(bubbleDiameter) /
(riseVelocity * riseVelocity);
WALBERLA_LOG_DEVEL("time step = " << t << "\n\t\tcenterOfMass = " << *centerOfMass
@@ -458,12 +456,12 @@ int main(int argc, char** argv)
const std::vector< real_t > resultVector{ (*centerOfMass)[2], riseVelocity, dragForce };
writeVectorToFile(resultVector, t, filename);
+
+ timestepOld = t;
+ centerOfMassOld = *centerOfMass;
}
}
- timestepOld = t;
- centerOfMassOld = *centerOfMass;
-
// stop simulation before bubble hits the top wall
if ((*centerOfMass)[2] > stoppingHeight) { break; }
--
GitLab