LeesEdwards.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 GeneratedOutflowBC.cpp
//! \author Markus Holzer <markus.holzer@fau.de>
//
//======================================================================================================================
#include "blockforest/Initialization.h"
#include "blockforest/communication/UniformBufferedScheme.h"

#include "core/Environment.h"
#include "core/timing/RemainingTimeLogger.h"

#include "field/AddToStorage.h"
#include "field/vtk/VTKWriter.h"

#include "timeloop/SweepTimeloop.h"

// Generated Files
#include "LeesEdwards_InfoHeader.h"

using namespace walberla;

using PackInfo_T  = lbm::LeesEdwards_PackInfo;
using flag_t      = walberla::uint8_t;
using FlagField_T = FlagField< flag_t >;

//////////
// MAIN //
//////////

class LeesEdwardsUpdate
{
 public:
   LeesEdwardsUpdate(const shared_ptr< StructuredBlockForest >& blocks, BlockDataID fieldID, BlockDataID tmpfieldID, real_t offset)
      : blocks_(blocks), fieldID_(fieldID), tmpfieldID_(tmpfieldID), offset_(offset)
   {}

   void operator()(IBlock* block)
   {
      // TODO should dimension_x contain the ghost layers or not. At the moment value is 64 with GL it is 66. In the
      // lbmpy Leed Edwards this is mixed. Probably not good


      // Top cells
      if (blocks_->atDomainYMaxBorder(*block))
      {
         uint_t dimension_x = blocks_->getNumberOfXCells(*block);
         real_t weight      = fmod(offset_ + real_c(dimension_x), 1.0);

         auto pdf_field = block->getData< PdfField_T >(fieldID_);
         auto pdf_tmp_field = block->getData< PdfField_T >(tmpfieldID_);

         CellInterval ci;
         pdf_field->getGhostRegion(stencil::N, ci, 1, true);

         // shift
         for (auto cell = ci.begin(); cell != ci.end(); ++cell)
         {
            cell_idx_t x = cell->x();

            uint_t ind1 = uint_c(floor(x - offset_ + real_c(dimension_x))) % dimension_x;
            uint_t ind2 = uint_c(ceil(x - offset_ + real_c(dimension_x))) % dimension_x;

            for (uint_t q = 0; q < Stencil_T::Q; ++q)
            {
               pdf_tmp_field->get(*cell, q) = (1 - weight) * pdf_field->get(cell_idx_c(ind1), cell->y(), cell->z(), q) +
                                              weight * pdf_field->get(cell_idx_c(ind2), cell->y(), cell->z(), q);
            }
         }

         // swap
         for (auto cell = ci.begin(); cell != ci.end(); ++cell)
         {
            for (uint_t q = 0; q < Stencil_T::Q; ++q)
            {
               pdf_field->get(*cell, q) = pdf_tmp_field->get(*cell, q);
            }
         }
      }
      // Bottom cells
      if (blocks_->atDomainYMinBorder(*block))
      {
         uint_t dimension_x = blocks_->getNumberOfXCells(*block);
         real_t weight      = fmod(offset_ + real_c(dimension_x), 1.0);

         auto pdf_field = block->getData< PdfField_T >(fieldID_);
         auto pdf_tmp_field = block->getData< PdfField_T >(tmpfieldID_);

         CellInterval ci;
         pdf_field->getGhostRegion(stencil::S, ci, 1, true);

         // shift
         for (auto cell = ci.begin(); cell != ci.end(); ++cell)
         {
            cell_idx_t x = cell->x();

            uint_t ind1 = uint_c(floor(x + offset_ + real_c(dimension_x))) % dimension_x;
            uint_t ind2 = uint_c(ceil(x + offset_ + real_c(dimension_x))) % dimension_x;

            for (uint_t q = 0; q < Stencil_T::Q; ++q)
            {
               pdf_tmp_field->get(*cell, q) = (1 - weight) * pdf_field->get(cell_idx_c(ind1), cell->y(), cell->z(), q) +
                                              weight * pdf_field->get(cell_idx_c(ind2), cell->y(), cell->z(), q);
            }
         }

         // swap
         for (auto cell = ci.begin(); cell != ci.end(); ++cell)
         {
            for (uint_t q = 0; q < Stencil_T::Q; ++q)
            {
               pdf_tmp_field->get(*cell, q) = pdf_field->get(*cell, q);
            }
         }
      }
   }

 private:
   const shared_ptr< StructuredBlockForest >& blocks_;
   BlockDataID fieldID_;
   BlockDataID tmpfieldID_;
   real_t offset_;
};

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

   auto blocks = blockforest::createUniformBlockGridFromConfig(walberlaEnv.config());

   // read parameters
   auto parameters = walberlaEnv.config()->getOneBlock("Parameters");

   const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(10));
   const real_t offset    = parameters.getParameter< real_t >("offset", real_c(10));

   const double remainingTimeLoggerFrequency =
      parameters.getParameter< double >("remainingTimeLoggerFrequency", 3.0); // in seconds

   // create fields
   BlockDataID pdfFieldID     = field::addToStorage< PdfField_T >(blocks, "PDFs", real_t(1.0), field::fzyx);
   BlockDataID pdfTmpFieldID  = field::addToStorage< PdfField_T >(blocks, "TmpPDFs", real_t(1.0), field::fzyx);
   BlockDataID velFieldID     = field::addToStorage< VectorField_T >(blocks, "velocity", real_t(0), field::fzyx);
   BlockDataID forceFieldID   = field::addToStorage< VectorField_T >(blocks, "force", real_t(0), field::fzyx);
   BlockDataID densityFieldID = field::addToStorage< ScalarField_T >(blocks, "density", real_t(1.0), field::fzyx);

   pystencils::LeesEdwards_Setter setterSweep(densityFieldID, forceFieldID, pdfFieldID);
   for (auto& block : *blocks)
      setterSweep(&block);

   // create impulsion
   {
      auto &block = *(blocks->begin());
      auto pdf_field = block.template getData<PdfField_T>(pdfFieldID);
      auto const dim_x = blocks->getNumberOfXCells(block);
      auto const dim_y = blocks->getNumberOfYCells(block);
      Cell const cell1{dim_x / 2, dim_y - 4, 0u};
      Cell const cell2{dim_x / 2, dim_y - 5, 0u};
      pdf_field->get(cell1, cell_idx_t(2)) = real_t(1e8);
      pdf_field->get(cell2, cell_idx_t(2)) = real_t(1e8);
   }

   // create time loop
   SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);

   // create communication for PdfField
   blockforest::communication::UniformBufferedScheme< Stencil_T > communication(blocks);
   communication.addPackInfo(make_shared< PackInfo_T >(pdfFieldID));

   pystencils::LeesEdwards_Collision CollisionSweep(densityFieldID, forceFieldID, pdfFieldID, velFieldID);
   pystencils::LeesEdwards_Stream StreamSweep(densityFieldID, forceFieldID, pdfFieldID, velFieldID);

   // add LBM sweep and communication to time loop
   timeloop.add() << Sweep(CollisionSweep, "collision");
   timeloop.add() << BeforeFunction(communication, "communication")
                  << Sweep(LeesEdwardsUpdate(blocks, pdfFieldID, pdfTmpFieldID, offset), "Lees Edwards");
   timeloop.add() << Sweep(StreamSweep, "stream");

   // log remaining time
   timeloop.addFuncAfterTimeStep(timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency),
                                 "remaining time logger");

   // VTK Writer
   uint_t vtkWriteFrequency = parameters.getParameter< uint_t >("vtkWriteFrequency", 0);
   if (vtkWriteFrequency > 0)
   {
      auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "Lees_Edwards", vtkWriteFrequency, 0, false, "vtk_out",
                                                      "simulation_step", false, true, true, false, 0);

      auto velWriter     = make_shared< field::VTKWriter< VectorField_T > >(velFieldID, "velocity");
      auto densityWriter = make_shared< field::VTKWriter< ScalarField_T > >(densityFieldID, "density");

      vtkOutput->addCellDataWriter(velWriter);
      vtkOutput->addCellDataWriter(densityWriter);

      timeloop.addFuncBeforeTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
   }
   timeloop.run();

   return EXIT_SUCCESS;
}