From 602c22e0aaf534899b9701fe3eedbee4c4143eba Mon Sep 17 00:00:00 2001
From: Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
Date: Wed, 5 Oct 2022 15:56:27 +0200
Subject: [PATCH] Call zero-setter in FSLBM only when writing VTK output

---
 .../FreeSurface/BubblyPoiseuille.prm          | 22 ++---
 apps/showcases/FreeSurface/CapillaryWave.prm  |  4 +-
 .../FreeSurface/DamBreakCylindrical.prm       |  4 +-
 .../FreeSurface/DamBreakRectangular.prm       | 23 ++---
 apps/showcases/FreeSurface/DropImpact.prm     | 23 ++---
 apps/showcases/FreeSurface/DropWetting.prm    | 23 ++---
 apps/showcases/FreeSurface/GravityWave.prm    |  4 +-
 apps/showcases/FreeSurface/MovingDrop.prm     | 23 ++---
 apps/showcases/FreeSurface/RisingBubble.prm   |  4 +-
 apps/showcases/FreeSurface/TaylorBubble.prm   |  3 +-
 src/lbm/free_surface/VtkWriter.h              | 90 +++++++++++--------
 11 files changed, 128 insertions(+), 95 deletions(-)

diff --git a/apps/showcases/FreeSurface/BubblyPoiseuille.prm b/apps/showcases/FreeSurface/BubblyPoiseuille.prm
index c6e7d413c..cfd43e279 100644
--- a/apps/showcases/FreeSurface/BubblyPoiseuille.prm
+++ b/apps/showcases/FreeSurface/BubblyPoiseuille.prm
@@ -85,20 +85,22 @@ VTK
          normal;
          obstacle_normal;
          mapped_flag;
-        }
+      }
 
-        inclusion_filters
-        {
-            // only include liquid and interface cells in VTK output
-            //liquidInterfaceFilter;
-        }
+      inclusion_filters
+      {
+         // only include liquid and interface cells in VTK output
+         //liquidInterfaceFilter;
+      }
 
-        before_functions
-        {
-            //ghost_layer_synchronization; // only needed if writing the ghost layer
-        }
+      before_functions
+      {
+         //ghost_layer_synchronization;   // only needed if writing the ghost layer
+         gas_cell_zero_setter;            // sets velocity=0 and density=1 all gas cells before writing VTK output
+      }
 
    }
+
    domain_decomposition
    {
       writeFrequency 0;
diff --git a/apps/showcases/FreeSurface/CapillaryWave.prm b/apps/showcases/FreeSurface/CapillaryWave.prm
index 8d172cadd..b73208639 100644
--- a/apps/showcases/FreeSurface/CapillaryWave.prm
+++ b/apps/showcases/FreeSurface/CapillaryWave.prm
@@ -95,9 +95,11 @@ VTK
 
       before_functions
       {
-         //ghost_layer_synchronization; // only needed if writing the ghost layer
+         //ghost_layer_synchronization;   // only needed if writing the ghost layer
+         gas_cell_zero_setter;            // sets velocity=0 and density=1 all gas cells before writing VTK output
       }
    }
+
    domain_decomposition
    {
       writeFrequency             100;
diff --git a/apps/showcases/FreeSurface/DamBreakCylindrical.prm b/apps/showcases/FreeSurface/DamBreakCylindrical.prm
index de81d688e..aee842982 100644
--- a/apps/showcases/FreeSurface/DamBreakCylindrical.prm
+++ b/apps/showcases/FreeSurface/DamBreakCylindrical.prm
@@ -98,10 +98,12 @@ VTK
 
         before_functions
         {
-            //ghost_layer_synchronization; // only needed if writing the ghost layer
+            //ghost_layer_synchronization;   // only needed if writing the ghost layer
+            gas_cell_zero_setter;            // sets velocity=0 and density=1 all gas cells before writing VTK output
         }
 
    }
+
    domain_decomposition
    {
       writeFrequency 248;
diff --git a/apps/showcases/FreeSurface/DamBreakRectangular.prm b/apps/showcases/FreeSurface/DamBreakRectangular.prm
index 34fdca4b0..8c4bc96e7 100644
--- a/apps/showcases/FreeSurface/DamBreakRectangular.prm
+++ b/apps/showcases/FreeSurface/DamBreakRectangular.prm
@@ -88,20 +88,21 @@ VTK
          normal;
          obstacle_normal;
          mapped_flag;
-        }
+      }
 
-        inclusion_filters
-        {
-            // only include liquid and interface cells in VTK output
-            //liquidInterfaceFilter;
-        }
-
-        before_functions
-        {
-            //ghost_layer_synchronization; // only needed if writing the ghost layer
-        }
+      inclusion_filters
+      {
+         // only include liquid and interface cells in VTK output
+         //liquidInterfaceFilter;
+      }
 
+      before_functions
+      {
+         //ghost_layer_synchronization;   // only needed if writing the ghost layer
+         gas_cell_zero_setter;            // sets velocity=0 and density=1 all gas cells before writing VTK output
+      }
    }
+
    domain_decomposition
    {
       writeFrequency 991;
diff --git a/apps/showcases/FreeSurface/DropImpact.prm b/apps/showcases/FreeSurface/DropImpact.prm
index 07afacce7..76580d98f 100644
--- a/apps/showcases/FreeSurface/DropImpact.prm
+++ b/apps/showcases/FreeSurface/DropImpact.prm
@@ -87,20 +87,21 @@ VTK
          normal;
          obstacle_normal;
          mapped_flag;
-        }
+      }
 
-        inclusion_filters
-        {
-            // only include liquid and interface cells in VTK output
-            //liquidInterfaceFilter;
-        }
-
-        before_functions
-        {
-            //ghost_layer_synchronization; // only needed if writing the ghost layer
-        }
+      inclusion_filters
+      {
+         // only include liquid and interface cells in VTK output
+         //liquidInterfaceFilter;
+      }
 
+      before_functions
+      {
+         //ghost_layer_synchronization;   // only needed if writing the ghost layer
+         gas_cell_zero_setter;            // sets velocity=0 and density=1 all gas cells before writing VTK output
+      }
    }
+
    domain_decomposition
    {
       writeFrequency 372;
diff --git a/apps/showcases/FreeSurface/DropWetting.prm b/apps/showcases/FreeSurface/DropWetting.prm
index bdfbe4e1e..e10ae9273 100644
--- a/apps/showcases/FreeSurface/DropWetting.prm
+++ b/apps/showcases/FreeSurface/DropWetting.prm
@@ -85,20 +85,21 @@ VTK
          normal;
          obstacle_normal;
          mapped_flag;
-        }
+      }
 
-        inclusion_filters
-        {
-            // only include liquid and interface cells in VTK output
-            //liquidInterfaceFilter;
-        }
-
-        before_functions
-        {
-            //ghost_layer_synchronization; // only needed if writing the ghost layer
-        }
+      inclusion_filters
+      {
+         // only include liquid and interface cells in VTK output
+         //liquidInterfaceFilter;
+      }
 
+      before_functions
+      {
+         //ghost_layer_synchronization;   // only needed if writing the ghost layer
+         gas_cell_zero_setter;            // sets velocity=0 and density=1 all gas cells before writing VTK output
+      }
    }
+
    domain_decomposition
    {
       writeFrequency 100;
diff --git a/apps/showcases/FreeSurface/GravityWave.prm b/apps/showcases/FreeSurface/GravityWave.prm
index 363a59c84..0400d62da 100644
--- a/apps/showcases/FreeSurface/GravityWave.prm
+++ b/apps/showcases/FreeSurface/GravityWave.prm
@@ -95,9 +95,11 @@ VTK
 
       before_functions
       {
-         //ghost_layer_synchronization; // only needed if writing the ghost layer
+         //ghost_layer_synchronization;   // only needed if writing the ghost layer
+         gas_cell_zero_setter;            // sets velocity=0 and density=1 all gas cells before writing VTK output
       }
    }
+
    domain_decomposition
    {
       writeFrequency             0;
diff --git a/apps/showcases/FreeSurface/MovingDrop.prm b/apps/showcases/FreeSurface/MovingDrop.prm
index 184963db1..dff95569f 100644
--- a/apps/showcases/FreeSurface/MovingDrop.prm
+++ b/apps/showcases/FreeSurface/MovingDrop.prm
@@ -85,20 +85,21 @@ VTK
          normal;
          obstacle_normal;
          mapped_flag;
-        }
+      }
 
-        inclusion_filters
-        {
-            // only include liquid and interface cells in VTK output
-            //liquidInterfaceFilter;
-        }
-
-        before_functions
-        {
-            //ghost_layer_synchronization; // only needed if writing the ghost layer
-        }
+      inclusion_filters
+      {
+         // only include liquid and interface cells in VTK output
+         //liquidInterfaceFilter;
+      }
 
+      before_functions
+      {
+         //ghost_layer_synchronization;   // only needed if writing the ghost layer
+         gas_cell_zero_setter;            // sets velocity=0 and density=1 all gas cells before writing VTK output
+      }
    }
+
    domain_decomposition
    {
       writeFrequency 0;
diff --git a/apps/showcases/FreeSurface/RisingBubble.prm b/apps/showcases/FreeSurface/RisingBubble.prm
index 5d483426b..62f1053d2 100644
--- a/apps/showcases/FreeSurface/RisingBubble.prm
+++ b/apps/showcases/FreeSurface/RisingBubble.prm
@@ -96,9 +96,11 @@ VTK
 
       before_functions
       {
-         //ghost_layer_synchronization; // only needed if writing the ghost layer
+         //ghost_layer_synchronization;   // only needed if writing the ghost layer
+         gas_cell_zero_setter;            // sets velocity=0 and density=1 all gas cells before writing VTK output
       }
    }
+
    domain_decomposition
    {
       writeFrequency             445;
diff --git a/apps/showcases/FreeSurface/TaylorBubble.prm b/apps/showcases/FreeSurface/TaylorBubble.prm
index ee3d8c9f5..e413afc7f 100644
--- a/apps/showcases/FreeSurface/TaylorBubble.prm
+++ b/apps/showcases/FreeSurface/TaylorBubble.prm
@@ -98,7 +98,8 @@ VTK
 
       before_functions
       {
-         //ghost_layer_synchronization; // only needed if writing the ghost layer
+         //ghost_layer_synchronization;   // only needed if writing the ghost layer
+         gas_cell_zero_setter;            // sets velocity=0 and density=1 all gas cells before writing VTK output
       }
    }
    domain_decomposition
diff --git a/src/lbm/free_surface/VtkWriter.h b/src/lbm/free_surface/VtkWriter.h
index f1ff2b93c..f98549a8f 100644
--- a/src/lbm/free_surface/VtkWriter.h
+++ b/src/lbm/free_surface/VtkWriter.h
@@ -120,6 +120,60 @@ void addVTKOutput(const std::weak_ptr< StructuredBlockForest >& blockForestPtr,
          std::make_shared< field::communication::PackInfo< VectorField_T > >(obstacleNormalFieldID));
 
       beforeFuncs["ghost_layer_synchronization"] = preVTKComm;
+
+      // set velocity and density to zero in obstacle and gas cells (only for visualization purposes); the PDF values in
+      // these cells are not important and thus not set during the simulation;
+      // only enable this functionality if the non-liquid and non-interface cells are not excluded anyway
+      const auto vtkConfigBlock        = config->getOneBlock("VTK");
+      const auto fluidFieldConfigBlock = vtkConfigBlock.getBlock("fluid_field");
+      if (fluidFieldConfigBlock)
+      {
+         auto inclusionFiltersConfigBlock = fluidFieldConfigBlock.getBlock("inclusion_filters");
+
+         // liquidInterfaceFilter limits VTK-output to only liquid and interface cells
+         if (!inclusionFiltersConfigBlock.isDefined("liquidInterfaceFilter"))
+         {
+            class ZeroSetter
+            {
+             public:
+               ZeroSetter(const weak_ptr< StructuredBlockForest >& blockForest, const BlockDataID& pdfFieldID,
+                          const ConstBlockDataID& flagFieldID,
+                          const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo)
+                  : blockForest_(blockForest), pdfFieldID_(pdfFieldID), flagFieldID_(flagFieldID), flagInfo_(flagInfo)
+               {}
+
+               void operator()()
+               {
+                  auto blockForest = blockForest_.lock();
+                  WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+                  for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+                  {
+                     PdfField_T* const pdfField         = blockIt->template getData< PdfField_T >(pdfFieldID_);
+                     const FlagField_T* const flagField = blockIt->template getData< const FlagField_T >(flagFieldID_);
+                     WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(pdfField, uint_c(1), {
+                        const typename PdfField_T::Ptr pdfFieldPtr(*pdfField, x, y, z);
+                        const typename FlagField_T::ConstPtr flagFieldPtr(*flagField, x, y, z);
+
+                        if (flagInfo_.isGas(*flagFieldPtr) || flagInfo_.isObstacle(*flagFieldPtr))
+                        {
+                           pdfField->setDensityAndVelocity(pdfFieldPtr.cell(), Vector3< real_t >(real_c(0)),
+                                                           real_c(1.0));
+                        }
+                     }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
+                  }
+               }
+
+             private:
+               weak_ptr< StructuredBlockForest > blockForest_;
+               BlockDataID pdfFieldID_;
+               ConstBlockDataID flagFieldID_;
+               typename FreeSurfaceBoundaryHandling_T::FlagInfo_T flagInfo_;
+            };
+
+            beforeFuncs["gas_cell_zero_setter"] = ZeroSetter(blockForest, pdfFieldID, flagFieldID, flagInfo);
+         }
+      }
    };
 
    // add VTK output to timeloop
@@ -130,42 +184,6 @@ void addVTKOutput(const std::weak_ptr< StructuredBlockForest >& blockForestPtr,
       timeloop.addFuncBeforeTimeStep(output->second.outputFunction, std::string("VTK: ") + output->first,
                                      output->second.requiredGlobalStates, output->second.incompatibleGlobalStates);
    }
-
-   // only enable the zerosetter (see below) if the non-liquid and non-interface cells are not excluded anyway
-   bool enableZeroSetter            = true;
-   const auto vtkConfigBlock        = config->getOneBlock("VTK");
-   const auto fluidFieldConfigBlock = vtkConfigBlock.getBlock("fluid_field");
-   if (fluidFieldConfigBlock)
-   {
-      auto inclusionFiltersConfigBlock = fluidFieldConfigBlock.getBlock("inclusion_filters");
-
-      if (inclusionFiltersConfigBlock.isDefined("liquidInterfaceFilter"))
-      {
-         // liquidInterfaceFilter is defined which limits VTK-output to only liquid and interface cells
-         enableZeroSetter = false;
-      }
-   }
-
-   // set velocity and density to zero in obstacle and gas cells (only for visualization purposes); the PDF values in
-   // these cells are not important and thus not set during the simulation
-   if (enableZeroSetter)
-   {
-      const auto function = [&](IBlock* block) {
-         using namespace free_surface;
-         PdfField_T* const pdfField         = block->getData< PdfField_T >(pdfFieldID);
-         const FlagField_T* const flagField = block->getData< const FlagField_T >(flagFieldID);
-         WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(pdfField, uint_c(1), {
-            const typename PdfField_T::Ptr pdfFieldPtr(*pdfField, x, y, z);
-            const typename FlagField_T::ConstPtr flagFieldPtr(*flagField, x, y, z);
-
-            if (flagInfo.isGas(*flagFieldPtr) || flagInfo.isObstacle(*flagFieldPtr))
-            {
-               pdfField->setDensityAndVelocity(pdfFieldPtr.cell(), Vector3< real_t >(0), real_c(1.0));
-            }
-         }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
-      };
-      timeloop.add() << Sweep(function, "VTK: zero-setting");
-   }
 }
 
 } // namespace free_surface
-- 
GitLab