ChannelFlow.cpp 5.41 KB
Newer Older
Camille Castells's avatar
Camille Castells committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
//======================================================================================================================
//
//  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 ChannelFlow.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 "geometry/InitBoundaryHandling.h"
#include "timeloop/SweepTimeloop.h"

// Generated Files
#include "ChannelFlow.h"

using namespace walberla;

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

auto pdfFieldAdder = [](IBlock* const block, StructuredBlockStorage* const storage) {
   return new PdfField_T(storage->getNumberOfXCells(*block), storage->getNumberOfYCells(*block),
                         storage->getNumberOfZCells(*block), uint_t(1), field::fzyx,
                         make_shared< field::AllocateAligned< real_t, 64 > >());
};

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 real_t omega     = parameters.getParameter< real_t >("omega");
   const uint_t timesteps = parameters.getParameter< uint_t >("timesteps");
   real_t ForceX = parameters.getParameter< real_t >("ForceX");

   WALBERLA_LOG_INFO_ON_ROOT("Using a relaxation rate of " << omega)
   WALBERLA_LOG_INFO_ON_ROOT("Set force " << ForceX << " in x direction.")

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

   // create fields
   BlockDataID pdfFieldID     = blocks->addStructuredBlockData< PdfField_T >(pdfFieldAdder, "PDFs");
   BlockDataID velFieldID     = field::addToStorage< VelocityField_T >(blocks, "velocity", real_t(0.0), field::fzyx);
   BlockDataID densityFieldID = field::addToStorage< ScalarField_T >(blocks, "density", real_t(1.0), field::fzyx);

   BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");

   pystencils::ChannelFlow_MacroSetter setterSweep(densityFieldID, pdfFieldID, velFieldID, ForceX);
   for (auto& block : *blocks)
      setterSweep(&block);

   // create and initialize boundary handling
   const FlagUID fluidFlagUID("Fluid");

   auto boundariesConfig = walberlaEnv.config()->getOneBlock("Boundaries");

   lbm::ChannelFlow_NoSlip noSlip(blocks, pdfFieldID);

   geometry::initBoundaryHandling< FlagField_T >(*blocks, flagFieldId, boundariesConfig);
   geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldId, fluidFlagUID);

   noSlip.fillFromFlagField< FlagField_T >(blocks, flagFieldId, FlagUID("NoSlip"), fluidFlagUID);

   // 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::ChannelFlow_Sweep UpdateSweep(densityFieldID, pdfFieldID, velFieldID, ForceX, omega);

   // add LBM sweep and communication to time loop
   timeloop.add() << Sweep(noSlip, "noSlip boundary") << AfterFunction(communication, "communication");
   timeloop.add() << Sweep(UpdateSweep, "LB stream & collide");

   // 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, "ChannelFlow_VTK", vtkWriteFrequency, 0, false,
                                                      "vtk_out", "simulation_step", false, true, true, false, 0);

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

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

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

   WALBERLA_LOG_INFO_ON_ROOT("Starting Simulation with " << timesteps << " timesteps")
   timeloop.run();
   WALBERLA_LOG_INFO_ON_ROOT("Finished Simulation")

   return EXIT_SUCCESS;
}