Merge branch 'fix-error-velocity-field-zalesak-disk' into 'master'

Fix errors in FSLBM advection test cases

See merge request walberla/walberla!596

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

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

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));
 }

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);