DropInPool.cpp

//======================================================================================================================
//
//  This file is part of waLBerla. waLBerla is free software: you can
//  redistribute it and/or modify it under the terms of the GNU General Public
//  License as published by the Free Software Foundation, either version 3 of
//  the License, or (at your option) any later version.
//
//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
//  for more details.
//
//  You should have received a copy of the GNU General Public License along
//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
//! \file DropInPool.cpp
//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
//
//======================================================================================================================

#include "blockforest/Initialization.h"

#include "core/Environment.h"

#include "field/adaptors/AdaptorCreators.h"
#include "field/vtk/FlagFieldCellFilter.h"
#include "field/vtk/FlagFieldMapping.h"
#include "field/vtk/VTKWriter.h"

#include "lbm/PerformanceLogger.h"
#include "lbm/field/Adaptors.h"
#include "lbm/field/AddToStorage.h"
#include "lbm/free_surface/bubble_model/Geometry.h"
#include "lbm/free_surface/dynamics/SurfaceDynamicsHandler.h"
#include "lbm/free_surface/surface_geometry/SurfaceGeometryHandler.h"
#include "lbm/lattice_model/D3Q19.h"

#include "vtk/Initialization.h"

using namespace walberla;

using CollisionModel_T = lbm::collision_model::SRT;
using LatticeModel_T   = lbm::D3Q19< CollisionModel_T, true, lbm::force_model::GuoConstant, 2 >;
using Stencil_T        = LatticeModel_T::Stencil;
using PdfField_T       = lbm::PdfField< LatticeModel_T >;

using Comm = blockforest::SimpleCommunication< LatticeModel_T::CommunicationStencil >;

using flag_t      = uint16_t;
using FlagField_T = FlagField< flag_t >;
using FsBoundary  = free_surface::FreeSurfaceBoundary< LatticeModel_T, FlagField_T >;

using ScalarField_T = GhostLayerField< real_t, 1 >;
using VectorField_T = GhostLayerField< Vector3< real_t >, 1 >;

int main(int argc, char** argv)
{
   Environment walberlaEnv(argc, argv);

   if (argc < 2) { WALBERLA_ABORT("Please specify a parameter file as input argument.") }

   // read domain parameters from parameter file
   auto domParam                            = walberlaEnv.config()->getOneBlock("DomainParameters");
   const Vector3< uint_t > numBlocks        = domParam.getParameter< Vector3< uint_t > >("blocks");
   const uint_t dropSize                    = domParam.getParameter< uint_t >("dropSize");
   const Vector3< uint_t > domainSizeFactor = domParam.getParameter< Vector3< uint_t > >("domainSizeFactor");
   const Vector3< uint_t > domainSize       = domainSizeFactor * dropSize;
   const Vector3< bool > periodic           = domParam.getParameter< Vector3< bool > >("periodic");
   const Vector3< uint_t > cellsPerBlock{ uint_c(ceil(real_c(domainSize[0]) / real_c(numBlocks[0]))),
                                          uint_c(ceil(real_c(domainSize[1]) / real_c(numBlocks[1]))),
                                          uint_c(ceil(real_c(domainSize[2]) / real_c(numBlocks[2]))) };

   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numBlocks);
   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(dropSize);
   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(domainSize);
   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(periodic);
   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellsPerBlock);

   // create uniform block grid
   auto blocks = blockforest::createUniformBlockGrid(numBlocks[0], numBlocks[1], numBlocks[2], cellsPerBlock[0],
                                                     cellsPerBlock[1], cellsPerBlock[2], real_c(1), true, periodic[0],
                                                     periodic[1], periodic[2], true);

   // read simulation parameters from parameter file
   auto simParam                 = walberlaEnv.config()->getOneBlock("SimulationParameters");
   const real_t omega            = simParam.getParameter< real_t >("omega");
   const Vector3< real_t > force = simParam.getParameter< Vector3< real_t > >("force");

   const uint_t timesteps       = simParam.getParameter< uint_t >("timesteps");
   const uint_t perfLogInterval = simParam.getParameter< uint_t >("perfLogInterval");

   const CollisionModel_T collisionModel = lbm::collision_model::SRT(omega);
   const real_t viscosity                = collisionModel.viscosity();

   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(omega);
   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(force);
   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(viscosity);

   // create lattice model
   LatticeModel_T latticeModel = LatticeModel_T(collisionModel, lbm::force_model::GuoConstant(force));

   // add pdf field
   const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blocks, "PDF field", latticeModel);

   // add fill level field (initialized with 0, i.e., gas everywhere)
   const BlockDataID fillFieldID = field::addToStorage< ScalarField_T >(blocks, "Fill level field", real_c(0.0));

   // add boundary handling
   FsBoundary boundary(blocks, pdfFieldID, fillFieldID);
   const BlockDataID flagFieldID           = boundary.getFlagFieldID();
   const typename FsBoundary::Flags& flags = boundary.getFlagInfo();

   // add various adaptors (for simplified access to macroscopic quantities)
   const BlockDataID densityAdaptorID =
      field::addFieldAdaptor< typename lbm::Adaptor< LatticeModel_T >::Density >(blocks, pdfFieldID, "DensityAdaptor");
   const BlockDataID velocityAdaptorID =
      field::addFieldAdaptor< typename lbm::Adaptor< LatticeModel_T >::VelocityVector >(blocks, pdfFieldID,
                                                                                        "VelocityVectorAdaptor");
   field::addFieldAdaptor< typename lbm::Adaptor< LatticeModel_T >::Velocity >(blocks, pdfFieldID, "VelocityAdaptor");

   // initialize a pool of liquid at the bottom of the domain in z-direction
   const AABB poolAABB(real_c(0), real_c(0), real_c(0), real_c(domainSize[0]), real_c(domainSize[1]),
                       real_c(domainSize[2]) * real_c(0.3));

   for (auto blockIterator = blocks->begin(); blockIterator != blocks->end(); ++blockIterator)
   {
      ScalarField_T* fillField = blockIterator->getData< ScalarField_T >(fillFieldID);

      // determine the cells that are relevant for a block (still in global coordinates)
      const CellInterval relevantCellBB = blocks->getCellBBFromAABB(poolAABB.getIntersection(blockIterator->getAABB()));

      // transform the global coordinates of relevant cells to block local coordinates
      CellInterval blockLocalCellBB;
      blocks->transformGlobalToBlockLocalCellInterval(blockLocalCellBB, *blockIterator, relevantCellBB);

      // clang-format off
         WALBERLA_FOR_ALL_CELLS_IN_INTERVAL_XYZ(blockLocalCellBB,
                                                fillField->get(x,y,z) = real_c(1);
         ) // WALBERLA_FOR_ALL_CELLS_IN_INTERVAL_XYZ
      // clang-format on
   }

   // initialize a cubical drop of liquid in the center of the domain
   const real_t dropSizeHalf = real_c(dropSize) * real_c(0.5);
   const AABB dropAABB(
      real_c(domainSize[0]) * real_c(0.5) - dropSizeHalf, real_c(domainSize[1]) * real_c(0.5) - dropSizeHalf,
      real_c(domainSize[2]) * real_c(0.75) - dropSizeHalf, real_c(domainSize[0]) * real_c(0.5) + dropSizeHalf,
      real_c(domainSize[1]) * real_c(0.5) + dropSizeHalf, real_c(domainSize[2]) * real_c(0.75) + dropSizeHalf);

   for (auto blockIterator = blocks->begin(); blockIterator != blocks->end(); ++blockIterator)
   {
      ScalarField_T* fillField = blockIterator->getData< ScalarField_T >(fillFieldID);

      // determine the cells that are relevant for a block (still in global coordinates)
      CellInterval relevantCellBB = blocks->getCellBBFromAABB(dropAABB.getIntersection(blockIterator->getAABB()));

      // transform the global coordinates of relevant cells to block local coordinates
      CellInterval blockLocalCellBB;
      blocks->transformGlobalToBlockLocalCellInterval(blockLocalCellBB, *blockIterator, relevantCellBB);

      // clang-format off
            WALBERLA_FOR_ALL_CELLS_IN_INTERVAL_XYZ(blockLocalCellBB,
                                                   fillField->get(x,y,z) = real_c(1);
            ) // WALBERLA_FOR_ALL_CELLS_IN_INTERVAL_XYZ
      // clang-format on
   }

   // initialize boundary conditions from config file
   auto boundaryParameters = walberlaEnv.config()->getOneBlock("BoundaryParameters");
   boundary.initFromConfig(boundaryParameters);
   boundary.initFlagsFromFillLevel();

   // initial communication
   Comm(blocks, pdfFieldID, fillFieldID, flagFieldID)();

   // add bubble model
   const bool disableSplits = true; // necessary if a gas bubble could split
   auto bubbleModel         = make_shared< bubble_model::BubbleModel >(blocks, disableSplits);
   bubbleModel->initFromFillLevelField(fillFieldID);

   // create timeloop
   SweepTimeloop tl(blocks, simParam.getParameter< uint_t >("timesteps"));

   // add surface geometry handler
   free_surface::SurfaceGeometryHandler< FlagField_T, ScalarField_T, VectorField_T > sgh(blocks, flagFieldID,
                                                                                         fillFieldID);
   sgh.setBoundaryFlagUIDs(FsBoundary::noSlipFlagID + FsBoundary::pressureFlag1ID + FsBoundary::pressureFlag2ID +
                           FsBoundary::outletFlagID + FsBoundary::ubbFlagID + FsBoundary::freeSlipFlagID);
   sgh.setContactAngle(simParam.getParameter< real_t >("contactAngle", 90));

   ConstBlockDataID curvatureField = sgh.getCurvatureFieldID();
   ConstBlockDataID normalField    = sgh.getNormalFieldID();

   sgh.addSweeps(tl);

   // add boundary handling for standard boundaries and free surface boundaries
   free_surface::SurfaceDynamicsHandler< LatticeModel_T, FlagField_T > dynamicsHandler(
      blocks, pdfFieldID, flagFieldID, fillFieldID, normalField, curvatureField, boundary, bubbleModel, real_c(0));

   dynamicsHandler.addSweeps(tl);

   // add VTK output
   auto vtkConfigFunc = [&](std::vector< shared_ptr< vtk::BlockCellDataWriterInterface > >& writers,
                            std::map< std::string, vtk::VTKOutput::CellFilter >& filters,
                            std::map< std::string, vtk::VTKOutput::BeforeFunction >& beforeFuncs) {
      using field::VTKWriter;
      writers.push_back(
         make_shared< VTKWriter< typename lbm::Adaptor< LatticeModel_T >::Density > >(densityAdaptorID, "density"));
      writers.push_back(make_shared< VTKWriter< typename lbm::Adaptor< LatticeModel_T >::VelocityVector > >(
         velocityAdaptorID, "velocity"));
      writers.push_back(make_shared< VTKWriter< PdfField_T, float > >(pdfFieldID, "PDFs"));
      writers.push_back(make_shared< VTKWriter< ScalarField_T, float > >(fillFieldID, "fill_levels"));
      writers.push_back(make_shared< VTKWriter< ScalarField_T, float > >(curvatureField, "curvature"));
      writers.push_back(make_shared< VTKWriter< VectorField_T, float > >(normalField, "normals"));
      writers.push_back(make_shared< VTKWriter< FlagField_T, float > >(flagFieldID, "flags"));

      auto flagMapper =
         make_shared< field::FlagFieldMapping< FlagField_T, walberla::uint8_t > >(flagFieldID, "mapped_flag_field");
      flagMapper->addMapping(FsBoundary::Flags::interfaceID, 1);
      flagMapper->addMapping(FsBoundary::Flags::liquidID, 2);
      flagMapper->addMapping(FsBoundary::Flags::gasID, 3);
      flagMapper->addMapping(FsBoundary::noSlipFlagID, 4);
      flagMapper->addMapping(FsBoundary::pressureFlag1ID, 5);
      flagMapper->addMapping(FsBoundary::pressureFlag2ID, 6);
      flagMapper->addMapping(FsBoundary::ubbFlagID, 7);

      writers.push_back(flagMapper);

      auto liquidInterfaceFilter = field::FlagFieldCellFilter< FlagField_T >(flagFieldID);
      liquidInterfaceFilter.addFlag(FsBoundary::Flags::liquidID);
      liquidInterfaceFilter.addFlag(FsBoundary::Flags::interfaceID);
      filters["liquidInterfaceFilter"] = liquidInterfaceFilter;

      auto preVTKComm = blockforest::communication::UniformBufferedScheme< stencil::D3Q27 >(blocks);
      preVTKComm.addPackInfo(make_shared< field::communication::PackInfo< PdfField_T > >(pdfFieldID));
      preVTKComm.addPackInfo(make_shared< field::communication::PackInfo< FlagField_T > >(flagFieldID));
      preVTKComm.addPackInfo(make_shared< field::communication::PackInfo< ScalarField_T > >(fillFieldID));

      beforeFuncs["sync"] = preVTKComm;
   };
   std::map< std::string, vtk::SelectableOutputFunction > vtkOutputFunctions;

   vtk::initializeVTKOutput(vtkOutputFunctions, vtkConfigFunc, blocks, walberlaEnv.config());

   for (auto output = vtkOutputFunctions.begin(); output != vtkOutputFunctions.end(); ++output)
   {
      tl.addFuncBeforeTimeStep(output->second.outputFunction, std::string("VTK: ") + output->first,
                               output->second.requiredGlobalStates, output->second.incompatibleGlobalStates);
   }

   // set velocity and density to zero in obstacle and gas cells (only for visualization purposes)
   // the PDF values in these fields are not important and thus not automatically set
   bool enableZeroSetter = true;
   if (enableZeroSetter)
   {
      auto function = [&](IBlock* block) {
         using namespace free_surface;
         PdfField_T* pdfFieldPtr   = block->getData< PdfField_T >(pdfFieldID);
         FlagField_T* flagFieldPtr = block->getData< FlagField_T >(flagFieldID);
         WALBERLA_FOR_ALL_CELLS(
            pdfIter, pdfFieldPtr, flagIter, flagFieldPtr,
            if (flags.isGas(flagIter) || flags.isObstacle(flagIter))
               pdfFieldPtr->setDensityAndVelocity(pdfIter.cell(), Vector3< real_t >(0), real_c(1.0));)
      };
      tl.add() << Sweep(function, "Zerosetter");
   }

   // add logging for computational performance
   lbm::PerformanceLogger< FlagField_T > perfLogger(blocks, flagFieldID, FsBoundary::Flags::liquidInterfaceIDs,
                                                    perfLogInterval);
   tl.addFuncAfterTimeStep(perfLogger, "PerformanceLogger");

   WcTimingPool timingPool;

   for (uint_t t = 0; t != timesteps; ++t)
   {
      if (t % uint_c(100) == uint_c(0)) { WALBERLA_LOG_DEVEL_ON_ROOT("Performing timestep=" << t); }
      tl.singleStep(timingPool, true);

      if (t % perfLogInterval == uint_c(0) && t > uint_c(0)) { timingPool.logResultOnRoot(); }
   }

   return EXIT_SUCCESS;
}