diff --git a/apps/benchmarks/FreeSurfaceAdvection/DeformationField.cpp b/apps/benchmarks/FreeSurfaceAdvection/DeformationField.cpp
index fcbc4ca54fe71364c93c275000f9b210276a56c7..d5060093c973a3dbe7a3cbfe3e274459bc80e484 100644
--- a/apps/benchmarks/FreeSurfaceAdvection/DeformationField.cpp
+++ b/apps/benchmarks/FreeSurfaceAdvection/DeformationField.cpp
@@ -19,8 +19,9 @@
// This benchmark simulates the advection of a spherical bubble in a vortex field with very high deformation. The vortex
// field changes periodically so that the bubble returns to its initial position, where it should take its initial form.
// The relative geometrical error of the bubble's shape after one period is evaluated. There is no LBM flow simulation
-// performed here. It is a test case for the FSLBM's mass advection only. This benchmark is based on Viktor Haag's
-// master thesis (https://www10.cs.fau.de/publications/theses/2017/Haag_MT_2017.pdf).
+// performed here. It is a test case for the FSLBM's mass advection only. This benchmark is based on the work from
+// Liovic et al. (doi: 10.1016/j.compfluid.2005.09.003). The setup is identical as in Viktor Haag's master thesis
+// (https://www10.cs.fau.de/publications/theses/2017/Haag_MT_2017.pdf).
//======================================================================================================================
#include "core/Environment.h"
@@ -226,8 +227,7 @@ int main(int argc, char** argv)
}
// create the spherical bubble
- const geometry::Sphere sphereBubble(real_c(domainWidth) * Vector3< real_t >(real_c(0.5), real_c(0.75), real_c(0.25)),
- real_c(domainWidth) * real_c(0.15));
+ const geometry::Sphere sphereBubble(bubbleCenter, bubbleDiameter * real_c(0.5));
bubble_model::addBodyToFillLevelField< geometry::Sphere >(*blockForest, fillFieldID, sphereBubble, true);
// initialize domain boundary conditions from config file
diff --git a/apps/benchmarks/FreeSurfaceAdvection/SingleVortex.cpp b/apps/benchmarks/FreeSurfaceAdvection/SingleVortex.cpp
index d3c8e01ed776fac44d0b33f7cc7e94a6508c08ed..ab3894858e957339198dd4499caf03b8d5456fb6 100644
--- a/apps/benchmarks/FreeSurfaceAdvection/SingleVortex.cpp
+++ b/apps/benchmarks/FreeSurfaceAdvection/SingleVortex.cpp
@@ -80,14 +80,16 @@ inline Vector3< real_t > velocityProfile(Cell globalCell, real_t timePeriod, uin
const real_t xToDomainCenter = x - real_c(0.5);
const real_t yToDomainCenter = y - real_c(0.5);
const real_t r = real_c(std::sqrt(xToDomainCenter * xToDomainCenter + yToDomainCenter * yToDomainCenter));
- const real_t rTerm = (real_c(1) - real_c(2) * r) * (real_c(1) - real_c(2) * r);
+ const real_t tmp = real_c(1) - real_c(2) * r;
+ const real_t rTerm = tmp * tmp;
const real_t timeTerm = real_c(std::cos(math::pi * real_t(timestep) / timePeriod));
- const real_t velocityX = real_c(std::sin(real_c(2) * math::pi * y)) * real_c(std::sin(math::pi * x)) *
- real_c(std::sin(math::pi * x)) * timeTerm;
- const real_t velocityY = -real_c(std::sin(real_c(2) * math::pi * x)) * real_c(std::sin(math::pi * y)) *
- real_c(std::sin(math::pi * y)) * timeTerm;
+ const real_t sinpix = real_c(std::sin(math::pi * x));
+ const real_t sinpiy = real_c(std::sin(math::pi * y));
+
+ const real_t velocityX = real_c(std::sin(real_c(2) * math::pi * y)) * sinpix * sinpix * timeTerm;
+ const real_t velocityY = -real_c(std::sin(real_c(2) * math::pi * x)) * sinpiy * sinpiy * timeTerm;
const real_t velocityZ = rTerm * timeTerm;
return Vector3< real_t >(velocityX, velocityY, velocityZ);
@@ -111,7 +113,7 @@ int main(int argc, char** argv)
auto domainParameters = walberlaEnv.config()->getOneBlock("DomainParameters");
const uint_t domainWidth = domainParameters.getParameter< uint_t >("domainWidth");
- const real_t bubbleDiameter = real_c(domainWidth) * real_c(0.075);
+ const real_t bubbleDiameter = real_c(domainWidth) * real_c(0.3);
const Vector3< real_t > bubbleCenter = domainWidth * Vector3< real_t >(real_c(0.5), real_c(0.75), real_c(0.25));
// define domain size
@@ -226,8 +228,7 @@ int main(int argc, char** argv)
}
// create the spherical bubble
- const geometry::Sphere sphereBubble(real_c(domainWidth) * Vector3< real_t >(real_c(0.5), real_c(0.75), real_c(0.25)),
- real_c(domainWidth) * real_c(0.15));
+ const geometry::Sphere sphereBubble(bubbleCenter, bubbleDiameter * real_c(0.5));
bubble_model::addBodyToFillLevelField< geometry::Sphere >(*blockForest, fillFieldID, sphereBubble, true);
// initialize domain boundary conditions from config file
diff --git a/apps/benchmarks/FreeSurfaceAdvection/ZalesakDisk.cpp b/apps/benchmarks/FreeSurfaceAdvection/ZalesakDisk.cpp
index f2177d55950a4a2a57e113b7dcc12344aa4b459b..15f6b9885eaa55c12dee46a94d81e662f4840e42 100644
--- a/apps/benchmarks/FreeSurfaceAdvection/ZalesakDisk.cpp
+++ b/apps/benchmarks/FreeSurfaceAdvection/ZalesakDisk.cpp
@@ -72,8 +72,8 @@ using FreeSurfaceBoundaryHandling_T = FreeSurfaceBoundaryHandling< LatticeModel_
inline Vector3< real_t > velocityProfile(real_t angularVelocity, Cell globalCell, const Vector3< real_t >& domainCenter)
{
// add 0.5 to get Cell's center
- const real_t velocityX = -angularVelocity * ((real_c(globalCell.y()) + real_c(0.5)) - domainCenter[0]);
- const real_t velocityY = angularVelocity * ((real_c(globalCell.x()) + real_c(0.5)) - domainCenter[1]);
+ const real_t velocityX = -angularVelocity * ((real_c(globalCell.y()) + real_c(0.5)) - domainCenter[1]);
+ const real_t velocityY = angularVelocity * ((real_c(globalCell.x()) + real_c(0.5)) - domainCenter[0]);
return Vector3< real_t >(velocityX, velocityY, real_c(0));
}
diff --git a/apps/benchmarks/FreeSurfaceAdvection/functionality/GeometricalErrorEvaluator.h b/apps/benchmarks/FreeSurfaceAdvection/functionality/GeometricalErrorEvaluator.h
index b0a5d7bbae749142f0a25c108a35b274ef146038..0aee14538c64114a6a02979f71ff102bb0f003ec 100644
--- a/apps/benchmarks/FreeSurfaceAdvection/functionality/GeometricalErrorEvaluator.h
+++ b/apps/benchmarks/FreeSurfaceAdvection/functionality/GeometricalErrorEvaluator.h
@@ -15,7 +15,7 @@
//
//! \file GeometricalErrorEvaluator.h
//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
-//! \brief Compute the geometrical error in free-surface advection test cases.
+//! \brief Compute the relative geometrical error in free-surface advection test cases.
//======================================================================================================================
#include "blockforest/StructuredBlockForest.h"
@@ -85,7 +85,7 @@ class GeometricalErrorEvaluator
{
initialFillLevelSum_ += *initialfillFieldIt;
}
- }) // WALBERLA_FOR_ALL_CELLS
+ }) // WALBERLA_FOR_ALL_CELLS_OMP
}
mpi::allReduceInplace< real_t >(initialFillLevelSum_, mpi::SUM);
@@ -98,7 +98,6 @@ class GeometricalErrorEvaluator
const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo();
real_t geometricalError = real_c(0);
- real_t fillLevelSum = real_c(0);
for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
{
@@ -107,14 +106,12 @@ class GeometricalErrorEvaluator
const FlagField_T* const flagField = blockIt->getData< const FlagField_T >(flagFieldID);
WALBERLA_FOR_ALL_CELLS_OMP(initialfillFieldIt, initialfillField, fillFieldIt, fillField, flagFieldIt,
- flagField, omp parallel for schedule(static) reduction(+:geometricalError)
- reduction(+:fillLevelSum), {
+ flagField, omp parallel for schedule(static) reduction(+:geometricalError), {
if (flagInfo.isInterface(flagFieldIt) || flagInfo.isLiquid(flagFieldIt))
{
geometricalError += real_c(std::abs(*initialfillFieldIt - *fillFieldIt));
- fillLevelSum += *fillFieldIt;
}
- }) // WALBERLA_FOR_ALL_CELLS
+ }) // WALBERLA_FOR_ALL_CELLS_OMP
}
mpi::allReduceInplace< real_t >(geometricalError, mpi::SUM);