LeesEdwards.cpp 8.3 KB
Newer Older
Markus Holzer's avatar
Markus Holzer 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
//======================================================================================================================
//
//  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 //
//////////

Markus Holzer's avatar
Markus Holzer committed
44
45
46
class LeesEdwardsUpdate
{
 public:
Jean-Noël Grad's avatar
Jean-Noël Grad committed
47
48
   LeesEdwardsUpdate(const shared_ptr< StructuredBlockForest >& blocks, BlockDataID fieldID, BlockDataID tmpfieldID, real_t offset)
      : blocks_(blocks), fieldID_(fieldID), tmpfieldID_(tmpfieldID), offset_(offset)
Markus Holzer's avatar
Markus Holzer committed
49
50
   {}

Markus Holzer's avatar
Markus Holzer committed
51
   void operator()(IBlock* block)
Markus Holzer's avatar
Markus Holzer committed
52
   {
Markus Holzer's avatar
Markus Holzer committed
53
54
55
56
57
58
59
60
61
62
63
      // 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_);
Jean-Noël Grad's avatar
Jean-Noël Grad committed
64
         auto pdf_tmp_field = block->getData< PdfField_T >(tmpfieldID_);
Markus Holzer's avatar
Markus Holzer committed
65
66
67
68

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

Jean-Noël Grad's avatar
Jean-Noël Grad committed
69
         // shift
Markus Holzer's avatar
Markus Holzer committed
70
71
72
73
         for (auto cell = ci.begin(); cell != ci.end(); ++cell)
         {
            cell_idx_t x = cell->x();

Jean-Noël Grad's avatar
Jean-Noël Grad committed
74
75
            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;
Markus Holzer's avatar
Markus Holzer committed
76
77
78

            for (uint_t q = 0; q < Stencil_T::Q; ++q)
            {
Jean-Noël Grad's avatar
Jean-Noël Grad committed
79
80
81
82
83
84
85
86
87
88
89
               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);
Markus Holzer's avatar
Markus Holzer committed
90
91
92
93
94
95
96
97
98
99
            }
         }
      }
      // 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_);
Jean-Noël Grad's avatar
Jean-Noël Grad committed
100
         auto pdf_tmp_field = block->getData< PdfField_T >(tmpfieldID_);
Markus Holzer's avatar
Markus Holzer committed
101
102
103
104

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

Jean-Noël Grad's avatar
Jean-Noël Grad committed
105
         // shift
Markus Holzer's avatar
Markus Holzer committed
106
107
108
109
         for (auto cell = ci.begin(); cell != ci.end(); ++cell)
         {
            cell_idx_t x = cell->x();

Jean-Noël Grad's avatar
Jean-Noël Grad committed
110
111
            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;
Markus Holzer's avatar
Markus Holzer committed
112
113
114

            for (uint_t q = 0; q < Stencil_T::Q; ++q)
            {
Jean-Noël Grad's avatar
Jean-Noël Grad committed
115
116
117
118
119
120
121
122
123
124
125
               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);
Markus Holzer's avatar
Markus Holzer committed
126
127
128
            }
         }
      }
Markus Holzer's avatar
Markus Holzer committed
129
130
131
   }

 private:
Markus Holzer's avatar
Markus Holzer committed
132
   const shared_ptr< StructuredBlockForest >& blocks_;
Markus Holzer's avatar
Markus Holzer committed
133
   BlockDataID fieldID_;
Jean-Noël Grad's avatar
Jean-Noël Grad committed
134
   BlockDataID tmpfieldID_;
Markus Holzer's avatar
Markus Holzer committed
135
136
137
   real_t offset_;
};

Markus Holzer's avatar
Markus Holzer committed
138
139
140
141
142
143
144
145
146
147
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));
Markus Holzer's avatar
Markus Holzer committed
148
   const real_t offset    = parameters.getParameter< real_t >("offset", real_c(10));
Markus Holzer's avatar
Markus Holzer committed
149
150
151
152
153
154

   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);
Jean-Noël Grad's avatar
Jean-Noël Grad committed
155
   BlockDataID pdfTmpFieldID  = field::addToStorage< PdfField_T >(blocks, "TmpPDFs", real_t(1.0), field::fzyx);
Markus Holzer's avatar
Markus Holzer committed
156
   BlockDataID velFieldID     = field::addToStorage< VectorField_T >(blocks, "velocity", real_t(0), field::fzyx);
Markus Holzer's avatar
Markus Holzer committed
157
   BlockDataID forceFieldID   = field::addToStorage< VectorField_T >(blocks, "force", real_t(0), field::fzyx);
Markus Holzer's avatar
Markus Holzer committed
158
159
160
161
162
163
   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);

164
165
166
167
168
169
170
171
172
173
174
175
   // 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);
   }

Markus Holzer's avatar
Markus Holzer committed
176
177
178
179
180
181
182
183
184
185
186
187
   // 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");
Markus Holzer's avatar
Markus Holzer committed
188
   timeloop.add() << BeforeFunction(communication, "communication")
Jean-Noël Grad's avatar
Jean-Noël Grad committed
189
                  << Sweep(LeesEdwardsUpdate(blocks, pdfFieldID, pdfTmpFieldID, offset), "Lees Edwards");
Markus Holzer's avatar
Markus Holzer committed
190
   timeloop.add() << Sweep(StreamSweep, "stream");
Markus Holzer's avatar
Markus Holzer committed
191
192
193
194
195
196
197
198
199

   // 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)
   {
Markus Holzer's avatar
Markus Holzer committed
200
201
      auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "Lees_Edwards", vtkWriteFrequency, 0, false, "vtk_out",
                                                      "simulation_step", false, true, true, false, 0);
Markus Holzer's avatar
Markus Holzer committed
202
203
204
205
206
207
208
209
210
211
212
213
214

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