01_SolvingPDE.cpp 12.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
//======================================================================================================================
//
//  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/>.
//
Jean-Noël Grad's avatar
Jean-Noël Grad committed
16
//! \file 01_SolvingPDE.cpp
17
18
19
20
21
22
23
24
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================

#include "blockforest/Initialization.h"
#include "blockforest/communication/UniformBufferedScheme.h"

#include "core/Environment.h"
25
#include "core/math/Constants.h"
26
27
28
29
30
31
32
33
34
35
36
37
38
39

#include "field/Field.h"
#include "field/AddToStorage.h"
#include "field/communication/PackInfo.h"

#include "stencil/D2Q5.h"

#include "gui/Gui.h"
#include "timeloop/SweepTimeloop.h"

#include <cmath>
#include <vector>


40
namespace walberla {
41

42
using ScalarField = GhostLayerField<real_t, 1>;
43
using Stencil_T = stencil::D2Q5;
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59


// function to initialize the boundaries of the source and destination fields
// The values of the Dirichlet boundary conditions are stored in the ghost layers of the outer blocks.
void initBC( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDataID & srcID, const BlockDataID & dstID )
{
   // iterate all blocks with an iterator 'block'
   for( auto block = blocks->begin(); block != blocks->end(); ++block )
   {
      // The Dirichlet boundary condition is non-zero only at the top boundary.
      if( blocks->atDomainYMaxBorder( *block ) )
      {
         // get the field data out of the blocks
         auto src = block->getData< ScalarField >( srcID );
         auto dst = block->getData< ScalarField >( dstID );

60
         // obtain a CellInterval object that holds information about the number of cells in x,y,z direction of the field including ghost layers
61
62
63
64
65
66
67
68
69
70
71
72
         // Since src and dst have the same size, one object is enough.
         CellInterval xyz = src->xyzSizeWithGhostLayer();

         // To only loop over the top row of the cells, i.e., the ghost layer containing the domain boundary, restrict the range in y-direction.
         xyz.yMin() = xyz.yMax();

         // iterate all cells in the top boundary with the iterator 'cell'
         for( auto cell = xyz.begin(); cell != xyz.end(); ++cell )
         {
            // obtain the physical coordinate of the center of the current cell
            const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );

73
74
75
            // set the values of the Dirichlet boundary condition given by the function u(x,1) = sin(2*PI*x)*sinh(2*PI) in the source and destination field
            src->get( *cell ) = std::sin( real_c(2) * math::pi * p[0] ) * std::sinh( real_c(2) * math::pi );
            dst->get( *cell ) = std::sin( real_c(2) * math::pi * p[0] ) * std::sinh( real_c(2) * math::pi );
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
         }
      }
   }
}



// function to initialize the field holding the right-hand side function f
void initRHS( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDataID & rhsID )
{
   // iterate all blocks with an iterator 'block'
   for( auto block = blocks->begin(); block != blocks->end(); ++block )
   {
      // get the field data out of the block
      auto f = block->getData< ScalarField >( rhsID );

      // obtain a CellInterval object that holds information about the number of cells in x,y,z direction of the field excluding ghost layers
      CellInterval xyz = f->xyzSize();

      // iterate all (inner) cells in the field
      for( auto cell = xyz.begin(); cell != xyz.end(); ++cell )
      {
         // obtain the physical coordinate of the center of the current cell
         const Vector3< real_t > p = blocks->getBlockLocalCellCenter( *block, *cell );

101
102
103
         // set the right-hand side, given by the function f(x,y) = 4*PI*PI*sin(2*PI*x)*sinh(2*PI*y)
         f->get( *cell ) = real_c(4) * math::pi * math::pi * std::sin( real_c(2) * math::pi * p[0] ) *
                           std::sinh( real_c(2) * math::pi * p[1] );
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
      }
   }
}



// class containing the Jacobi sweep
class JacobiSweep
{
public:

   // constructor
   JacobiSweep( const BlockDataID & srcID, const BlockDataID & dstID, const BlockDataID & rhsID, const real_t & dx, const real_t & dy )
      : srcID_( srcID ), dstID_( dstID ), rhsID_( rhsID ), dx_( dx ), dy_( dy ) {}

   // functor that carries out the Jacobi sweep
   void operator()(IBlock * const block);

private:

   const BlockDataID srcID_;
   const BlockDataID dstID_;
   const BlockDataID rhsID_;

   const real_t dx_;
   const real_t dy_;
};

void JacobiSweep::operator()( IBlock * const block )
{
   // get the source, destination, and rhs field data from the block
   auto src = block->getData< ScalarField >( srcID_ );
   auto dst = block->getData< ScalarField >( dstID_ );
   auto rhs = block->getData< ScalarField >( rhsID_ );

   // iterate all cells (excluding ghost layers) of the fields
   // This macro is given a field (here src) with the matching size.
   // Inside the macro, x, y and z can be used as indices to access the data.
   // Since the three fields have the same size, these indices are valid for all of them.
   WALBERLA_FOR_ALL_CELLS_XYZ(src,
      // carries out the sweep for the current cell (x,y,z) with the respective prefactors
      dst->get(x,y,z) =  rhs->get(x,y,z);
      dst->get(x,y,z) += ( real_c(1) / (dx_ * dx_) ) * src->get( x+1,  y , z );
      dst->get(x,y,z) += ( real_c(1) / (dx_ * dx_) ) * src->get( x-1,  y , z );
      dst->get(x,y,z) += ( real_c(1) / (dy_ * dy_) ) * src->get(  x , y+1, z );
      dst->get(x,y,z) += ( real_c(1) / (dy_ * dy_) ) * src->get(  x , y-1, z );
150
      dst->get(x,y,z) /= ( real_c(2) / (dx_ * dx_) + real_c(2)/(dy_ * dy_) + real_c(4) * math::pi * math::pi );
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
   )

   // swap source and destination fields
   src->swapDataPointers( dst );
}



// class containing the Jacobi sweep using the stencil concept
class JacobiSweepStencil
{
public:

   // constructor
   JacobiSweepStencil( const BlockDataID & srcID, const BlockDataID & dstID, const BlockDataID & rhsID, const std::vector< real_t > & weights )
      : srcID_( srcID ), dstID_( dstID ), rhsID_( rhsID ), weights_( weights )
   {
      // store the center weight in inverse form to avoid divisions in the sweep
      // weights_[ Stencil_T::idx[ stencil::C ] ] = real_c(1) / weights_[ Stencil_T::idx[ stencil::C ] ];
   }

   // functor that carries out the Jacobi sweep
   void operator()( IBlock * block );

private:

   const BlockDataID srcID_;
   const BlockDataID dstID_;
   const BlockDataID rhsID_;

   std::vector< real_t > weights_;
};

void JacobiSweepStencil::operator()( IBlock * block )
{
   // get the source, destination and rhs field data from the block
   auto src = block->getData< ScalarField >( srcID_ );
   auto dst = block->getData< ScalarField >( dstID_ );
   auto rhs = block->getData< ScalarField >( rhsID_ );

   // iterate all cells (excluding ghost layers) of the fields and carry out the Jacobi sweep
   WALBERLA_FOR_ALL_CELLS_XYZ(src,

      dst->get(x,y,z) =  rhs->get(x,y,z);

      // iterate the neighboring cells and multiply their value with the respective weight
      for( auto dir = Stencil_T::beginNoCenter(); dir != Stencil_T::end(); ++dir )
         dst->get(x,y,z) -= weights_[ dir.toIdx() ] * src->getNeighbor(x,y,z,*dir);

      dst->get(x,y,z) /= weights_[ Stencil_T::idx[ stencil::C ] ];

      // When the center weight is stored in inverse form, a multiplication instead of a division can be applied here.
      // dst->get(x,y,z) *= weights_[ Stencil_T::idx[ stencil::C ] ];
   )

   // swap source and destination fields
   src->swapDataPointers( dst );
}



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

   // number of blocks in x,y,z direction
   const uint_t xBlocks = uint_c(1);
   const uint_t yBlocks = uint_c(1);
   const uint_t zBlocks = uint_c(1);

   // number of cells per block in x,y,z direction
   // two dimensional layout, so only one cell in z direction
   const uint_t xCells = uint_c(50);
   const uint_t yCells = uint_c(25);
   const uint_t zCells = uint_c(1);

   // domain size is [0,2]x[0,1]
   const real_t xMin = real_c(0);
   const real_t xMax = real_c(2);
   const real_t yMin = real_c(0);
   const real_t yMax = real_c(1);

   // mesh spacings in x and y direction
   const real_t dx = (xMax - xMin)/real_c( xBlocks*xCells + uint_c(1) );
   const real_t dy = (yMax - yMin)/real_c( yBlocks*yCells + uint_c(1) );

   // create an axis-aligned bounding box to define domain
   // defines a rectangular domain by specifying (xMin, yMin, zMin, xMax, yMax, zMax)
   // care has to be taken regarding the cell-centered data layout in waLBerla
   // in z-direction only one cell layer with arbitrary thickness, here dx, is present
   auto aabb = math::AABB( xMin + real_c(0.5)*dx, yMin + real_c(0.5)*dy, real_c(0),
                           xMax - real_c(0.5)*dx, yMax - real_c(0.5)*dy, dx );

   // create blocks
   shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid (
             aabb,                               // axis-aligned bounding box
             xBlocks, yBlocks, zBlocks,          // number of blocks in x,y,z direction
             xCells,  yCells,  zCells,           // how many cells per block (x,y,z)
             false,                              // one block per process - "false" means all blocks to one process
             false, false, false );              // no periodicity


   // add fields with ghost layers to all blocks
   // source and destination fields for the unknowns u
   BlockDataID srcID = field::addToStorage< ScalarField >( blocks, "src", real_c(0), field::zyxf, uint_c(1));
   BlockDataID dstID = field::addToStorage< ScalarField >( blocks, "dst", real_c(0), field::zyxf, uint_c(1));
   // field to store the function f
   BlockDataID rhsID = field::addToStorage< ScalarField >( blocks, "rhs", real_c(0), field::zyxf, uint_c(1));

   // initialize the field
   initRHS( blocks, rhsID );

   // set the Dirichlet boundary conditions
   initBC( blocks, srcID, dstID );

   // create the communication scheme for the D2Q5 stencil
   blockforest::communication::UniformBufferedScheme< stencil::D2Q5 > myCommScheme( blocks );
   // add a PackInfo that packs/unpacks the source field
   // Since the RHS field is unchanged, there is no need to communicate it.
   myCommScheme.addPackInfo( make_shared< field::communication::PackInfo<ScalarField> >( srcID ) );

   // set up the timeloop object
   SweepTimeloop timeloop ( blocks, uint_c(1) );

   // register the communication function and the Jacobi sweep in the timeloop
   // either the hard coded version...
   /*
   timeloop.add() << BeforeFunction( myCommScheme, "Communication" )
                  << Sweep( JacobiSweep( srcID, dstID, rhsID, dx, dy ), "JacobiSweep" );
   */

   // ...or the variant using the stencil concept
   // set up the stencil weights
   std::vector< real_t > weights( Stencil_T::Size );
285
   weights[ Stencil_T::idx[ stencil::C ] ] = real_c(2) / ( dx * dx ) + real_c(2) / ( dy * dy ) + real_c(4) * math::pi * math::pi;
286
287
288
289
290
291
292
293
294
295
296
297
298
299
   weights[ Stencil_T::idx[ stencil::E ] ] = real_c(-1) / ( dx * dx );
   weights[ Stencil_T::idx[ stencil::W ] ] = real_c(-1) / ( dx * dx );
   weights[ Stencil_T::idx[ stencil::N ] ] = real_c(-1) / ( dy * dy );
   weights[ Stencil_T::idx[ stencil::S ] ] = real_c(-1) / ( dy * dy );

   timeloop.add() << BeforeFunction( myCommScheme, "Communication" )
                  << Sweep( JacobiSweepStencil( srcID, dstID, rhsID, weights ), "JacobiSweepStencil" );

   // start the GUI and run the simulation
   GUI gui ( timeloop, blocks, argc, argv );
   gui.run();

   return 0;
}
300
301
302
303
304
305
}

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