SuViscoelasticityTest.cpp 12.7 KB
Newer Older
Cameron Stewart's avatar
Cameron Stewart committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
//======================================================================================================================
//
//  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 SuViscoelasticityTest.cpp
//! \ingroup lbm
//! \author Cameron Stewart <cstewart@icp.uni-stuttgart.de>
Cameron Stewart's avatar
Cameron Stewart committed
19
20
//! \brief D2Q9 Poiseuille flow simulation with Oldroyd-B viscoelasticity - Checks against analytical formula and
//!        reference data
Cameron Stewart's avatar
Cameron Stewart committed
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
//
//
//======================================================================================================================

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

#include "boundary/BoundaryHandling.h"

#include "core/Environment.h"
#include "core/SharedFunctor.h"

#include "domain_decomposition/SharedSweep.h"

#include "field/communication/PackInfo.h"

#include "lbm/boundary/all.h"
Cameron Stewart's avatar
Cameron Stewart committed
38
39
#include "lbm/field/AddToStorage.h"
#include "lbm/lattice_model/D2Q9.h"
Cameron Stewart's avatar
Cameron Stewart committed
40
41
42
43
44
45
#include "lbm/lattice_model/ForceModel.h"
#include "lbm/SuViscoelasticity.h"
#include "lbm/sweeps/CellwiseSweep.h"

#include "timeloop/SweepTimeloop.h"

cameron-stewart's avatar
cameron-stewart committed
46
namespace walberla {
Cameron Stewart's avatar
Cameron Stewart committed
47
48
49
50
51
52

//////////////
// TYPEDEFS //
//////////////

typedef GhostLayerField< Vector3<real_t>, 1> ForceField_T;
Cameron Stewart's avatar
Cameron Stewart committed
53
54
55
56
57
58

typedef lbm::collision_model::TRT CollisionModel_T;
typedef lbm::force_model::GuoField< ForceField_T > ForceModel_T;
typedef lbm::D2Q9< CollisionModel_T, false, ForceModel_T >  LatticeModel_T;
typedef LatticeModel_T::Stencil    Stencil_T;

Cameron Stewart's avatar
Cameron Stewart committed
59
60
61
62
63
typedef GhostLayerField<Matrix3<real_t>, 1> StressField_T;
typedef GhostLayerField< Vector3<real_t>, 1> VelocityField_T;
typedef lbm::PdfField< LatticeModel_T >  PdfField_T;
typedef walberla::uint8_t flag_t;
typedef FlagField< flag_t > FlagField_T;
Cameron Stewart's avatar
Cameron Stewart committed
64
const uint_t FieldGhostLayers = 2;
Cameron Stewart's avatar
Cameron Stewart committed
65

Cameron Stewart's avatar
Cameron Stewart committed
66
typedef lbm::NoSlip< LatticeModel_T, flag_t> NoSlip_T;
Michael Kuron's avatar
Michael Kuron committed
67
typedef BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T> BoundaryHandling_T;
Cameron Stewart's avatar
Cameron Stewart committed
68

Cameron Stewart's avatar
Cameron Stewart committed
69
70
71
///////////
// FLAGS //
///////////
Cameron Stewart's avatar
Cameron Stewart committed
72

Cameron Stewart's avatar
Cameron Stewart committed
73
74
const FlagUID Fluid_Flag( "fluid" );
const FlagUID NoSlip_Flag( "no slip" );
Cameron Stewart's avatar
Cameron Stewart committed
75

Cameron Stewart's avatar
Cameron Stewart committed
76
77
78
/////////////////////////////////////
// BOUNDARY HANDLING CUSTOMIZATION //
/////////////////////////////////////
Cameron Stewart's avatar
Cameron Stewart committed
79

Cameron Stewart's avatar
Cameron Stewart committed
80
81
82
class MyBoundaryHandling
{
public:
Cameron Stewart's avatar
Cameron Stewart committed
83

Cameron Stewart's avatar
Cameron Stewart committed
84
   MyBoundaryHandling( const BlockDataID & flagFieldID, const BlockDataID & pdfFieldID ) : flagFieldID_( flagFieldID ), pdfFieldID_( pdfFieldID ) {}
Cameron Stewart's avatar
Cameron Stewart committed
85

Cameron Stewart's avatar
Cameron Stewart committed
86
   BoundaryHandling_T * operator()( IBlock* const block, const StructuredBlockStorage* const storage ) const;
Cameron Stewart's avatar
Cameron Stewart committed
87

Cameron Stewart's avatar
Cameron Stewart committed
88
private:
Cameron Stewart's avatar
Cameron Stewart committed
89

Cameron Stewart's avatar
Cameron Stewart committed
90
91
   const BlockDataID flagFieldID_;
   const BlockDataID pdfFieldID_;
Cameron Stewart's avatar
Cameron Stewart committed
92

Cameron Stewart's avatar
Cameron Stewart committed
93
}; // class MyBoundaryHandling
Cameron Stewart's avatar
Cameron Stewart committed
94
95


Cameron Stewart's avatar
Cameron Stewart committed
96
97
98
99
BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const StructuredBlockStorage * const storage ) const
{
   WALBERLA_ASSERT_NOT_NULLPTR( block );
   WALBERLA_ASSERT_NOT_NULLPTR( storage );
Cameron Stewart's avatar
Cameron Stewart committed
100

Cameron Stewart's avatar
Cameron Stewart committed
101
102
   FlagField_T * flagField       = block->getData< FlagField_T >( flagFieldID_ );
   PdfField_T *  pdfField        = block->getData< PdfField_T > ( pdfFieldID_ );
Cameron Stewart's avatar
Cameron Stewart committed
103

Cameron Stewart's avatar
Cameron Stewart committed
104
   const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
Cameron Stewart's avatar
Cameron Stewart committed
105

Cameron Stewart's avatar
Cameron Stewart committed
106
   BoundaryHandling_T * handling = new BoundaryHandling_T( "moving obstacle boundary handling", flagField, fluid,
Michael Kuron's avatar
Michael Kuron committed
107
                                                           NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ) );
Cameron Stewart's avatar
Cameron Stewart committed
108

Cameron Stewart's avatar
Cameron Stewart committed
109
   const auto noSlip = flagField->getFlag(NoSlip_Flag);
Cameron Stewart's avatar
Cameron Stewart committed
110

Cameron Stewart's avatar
Cameron Stewart committed
111
112
113
   CellInterval domainBB = storage->getDomainCellBB();
   domainBB.xMin() -= cell_idx_c( 1 );
   domainBB.xMax() += cell_idx_c( 1 );
Cameron Stewart's avatar
Cameron Stewart committed
114

Cameron Stewart's avatar
Cameron Stewart committed
115
116
   domainBB.yMin() -= cell_idx_c( 1 );
   domainBB.yMax() += cell_idx_c( 1 );
Cameron Stewart's avatar
Cameron Stewart committed
117

Cameron Stewart's avatar
Cameron Stewart committed
118
119
120
121
122
123
124
125
126
   // SOUTH
   CellInterval south( domainBB.xMin(), domainBB.yMin(), domainBB.zMin(), domainBB.xMax(), domainBB.yMin(), domainBB.zMax() );
   storage->transformGlobalToBlockLocalCellInterval( south, *block );
   handling->forceBoundary( noSlip, south );

   // NORTH
   CellInterval north( domainBB.xMin(), domainBB.yMax(), domainBB.zMin(), domainBB.xMax(), domainBB.yMax(), domainBB.zMax() );
   storage->transformGlobalToBlockLocalCellInterval( north, *block );
   handling->forceBoundary( noSlip, north );
Cameron Stewart's avatar
Cameron Stewart committed
127

Cameron Stewart's avatar
Cameron Stewart committed
128
129
130
131
   handling->fillWithDomain( FieldGhostLayers );

   return handling;
}
Cameron Stewart's avatar
Cameron Stewart committed
132

Cameron Stewart's avatar
Cameron Stewart committed
133
134
135
//////////////////////
// Poiseuille Force //
//////////////////////
Cameron Stewart's avatar
Cameron Stewart committed
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165

template < typename BoundaryHandling_T >
class ConstantForce
{
public:

   ConstantForce( BlockDataID forceFieldId, BlockDataID boundaryHandlingId, real_t force)
         : forceFieldId_( forceFieldId ), boundaryHandlingId_(boundaryHandlingId), force_(force)
   {}

   void operator()( IBlock * block )
   {
      ForceField_T *forceField = block->getData< ForceField_T >(forceFieldId_);
      BoundaryHandling_T *boundaryHandling = block->getData< BoundaryHandling_T >( boundaryHandlingId_ );

      WALBERLA_FOR_ALL_CELLS_XYZ(forceField,
                                 {
                                    Cell cell(x,y,z);
                                    if (boundaryHandling->isDomain(cell)) {
                                       forceField->get(cell)[0] += force_;
                                    }
                                 })
   }

private:

   BlockDataID forceFieldId_, boundaryHandlingId_;
   real_t force_;
};

Cameron Stewart's avatar
Cameron Stewart committed
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
////////////////////
// Take Test Data //
////////////////////

class TestData
{
public:
   TestData(Timeloop & timeloop, shared_ptr< StructuredBlockForest > blocks, BlockDataID pdfFieldId, BlockDataID stressFieldId, uint_t timesteps, uint_t blockSize, real_t L, real_t H, real_t uExpected)
            : timeloop_(timeloop), blocks_(blocks), pdfFieldId_(pdfFieldId), stressFieldId_(stressFieldId), timesteps_(timesteps), blockSize_(blockSize), t0_(0), t1_(0), L_(L), H_(H), uMax_(0.0), uPrev_(0.0),
            uExpected_(uExpected), uSteady_(0.0) {}

   void operator()() {
      for (auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt) {
         PdfField_T *pdf = blockIt.get()->getData<PdfField_T>(pdfFieldId_);

181
         if (blockIt.get()->getAABB().contains(float_c(L_/2.0), float_c(H_/2.0), 0)) {
Cameron Stewart's avatar
Cameron Stewart committed
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
            uCurr_ = pdf->getVelocity(int32_c(blockSize_/2), int32_c(blockSize_/2), 0)[0]/uExpected_;
            tCurr_ = timeloop_.getCurrentTimeStep();
            if (tCurr_ == timesteps_ - 1){
               uSteady_ = uCurr_;
            }
            if (maxFlag_ == 0) {
               if (uCurr_ >= uPrev_) {
                  uMax_ = uCurr_;
               } else {
                  t0_ = tCurr_;
                  maxFlag_ = 1;
               }
               uPrev_ = uCurr_;
            }
            else if (maxFlag_ == 1) {
               if ((uCurr_ - 1.0) <= (uMax_ - 1.0) / std::exp(1)) {
                  t1_ = tCurr_ - t0_;
                  maxFlag_ = 2;
               }
            }
         }
      }
   }

   real_t getUSteady(){
      return uSteady_;
   }

   real_t getUMax(){
      return uMax_;
   }

   real_t getT0(){
      return real_c(t0_);
   }

   real_t getT1(){
      return real_c(t1_);
   }

private:
   Timeloop & timeloop_;
   shared_ptr< StructuredBlockForest > blocks_;
   BlockDataID pdfFieldId_, stressFieldId_;
   uint_t timesteps_, blockSize_, t0_, t1_, tCurr_;
   uint_t maxFlag_ = 0;
   real_t L_, H_, uMax_, uPrev_, uCurr_, uExpected_, uSteady_;
};

cameron-stewart's avatar
cameron-stewart committed
231
} // namespace walberla
Cameron Stewart's avatar
Cameron Stewart committed
232
233
234
235
236

//////////
// Main //
//////////

Cameron Stewart's avatar
Cameron Stewart committed
237
int main(int argc, char ** argv ){
cameron-stewart's avatar
cameron-stewart committed
238
   using namespace walberla;
Cameron Stewart's avatar
Cameron Stewart committed
239

cameron-stewart's avatar
cameron-stewart committed
240
   Environment env( argc, argv );
Cameron Stewart's avatar
Cameron Stewart committed
241

Cameron Stewart's avatar
Cameron Stewart committed
242
   // read parameter
Cameron Stewart's avatar
Cameron Stewart committed
243
244
245
246
   shared_ptr<StructuredBlockForest> blocks = blockforest::createUniformBlockGridFromConfig( env.config() );
   auto parameters = env.config()->getOneBlock( "Parameters" );

   // extract some constants from the parameters
247
248
249
250
251
252
253
254
255
256
   const real_t eta_s     = parameters.getParameter< real_t > ("eta_s");
   const real_t force     = parameters.getParameter< real_t > ("force");
   const real_t eta_p     = parameters.getParameter< real_t > ("eta_p");
   const real_t lambda_p  = parameters.getParameter< real_t > ("lambda_p");
   const uint_t period    = parameters.getParameter< uint_t > ("period");
   const real_t L         = parameters.getParameter< real_t > ("L");
   const real_t H         = parameters.getParameter< real_t > ("H");
   const uint_t blockSize = parameters.getParameter< uint_t > ("blockSize");
   const bool   shortrun  = parameters.getParameter<   bool > ("shortrun");
   const uint_t timesteps = shortrun ? uint_t(2) : parameters.getParameter< uint_t > ("timesteps");
Cameron Stewart's avatar
Cameron Stewart committed
257

Cameron Stewart's avatar
Cameron Stewart committed
258
   // reference data
259
   const real_t          uExpected       = force*H*H/(real_c(8.0)*(eta_s + eta_p));
Cameron Stewart's avatar
Cameron Stewart committed
260
261
262
   const real_t          uMax            = parameters.getParameter< real_t > ("uMax");
   const real_t          t0              = parameters.getParameter< real_t > ("t0");
   const real_t          t1              = parameters.getParameter< real_t > ("t1");
Cameron Stewart's avatar
Cameron Stewart committed
263
264

   // create fields
Cameron Stewart's avatar
Cameron Stewart committed
265
266
   BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field", FieldGhostLayers);
   BlockDataID forceFieldId = field::addToStorage<ForceField_T>( blocks, "Force Field", Vector3<real_t>(0.0), field::zyxf, FieldGhostLayers);
Cameron Stewart's avatar
Cameron Stewart committed
267
   LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::TRT::constructWithMagicNumber( walberla::lbm::collision_model::omegaFromViscosity(eta_s)), lbm::force_model::GuoField<ForceField_T>( forceFieldId ) );
Cameron Stewart's avatar
Cameron Stewart committed
268
   BlockDataID pdfFieldId = lbm::addPdfFieldToStorage( blocks, "pdf field", latticeModel, Vector3<real_t>(), real_c(1.0), FieldGhostLayers );
Cameron Stewart's avatar
Cameron Stewart committed
269
270
271
   BlockDataID stressId = walberla::field::addToStorage<StressField_T>( blocks, "Stress Field", Matrix3<real_t>(0.0), field::zyxf, FieldGhostLayers);
   BlockDataID stressOldId = walberla::field::addToStorage<StressField_T>( blocks, "Old Stress Field", Matrix3<real_t>(0.0), field::zyxf, FieldGhostLayers);
   BlockDataID velocityId = walberla::field::addToStorage<VelocityField_T> (blocks, "Velocity Field", Vector3<real_t>(0.0), field::zyxf, FieldGhostLayers);
Cameron Stewart's avatar
Cameron Stewart committed
272

Cameron Stewart's avatar
Cameron Stewart committed
273
274
   // add boundary handling
   BlockDataID boundaryHandlingId = blocks->addStructuredBlockData< BoundaryHandling_T >( MyBoundaryHandling( flagFieldId, pdfFieldId ), "boundary handling" );
Cameron Stewart's avatar
Cameron Stewart committed
275
276
277
278
279

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

   // create communication for PdfField
Cameron Stewart's avatar
Cameron Stewart committed
280
   blockforest::communication::UniformBufferedScheme< Stencil_T > communication( blocks );
Cameron Stewart's avatar
Cameron Stewart committed
281
   communication.addPackInfo( make_shared< field::communication::PackInfo< PdfField_T > >( pdfFieldId ) );
Cameron Stewart's avatar
Cameron Stewart committed
282
   auto testData = make_shared< TestData >(TestData(timeloop, blocks, pdfFieldId, stressId, timesteps, blockSize, L, H, uExpected));
Cameron Stewart's avatar
Cameron Stewart committed
283

Cameron Stewart's avatar
Cameron Stewart committed
284
   // structure timeloop
Cameron Stewart's avatar
Cameron Stewart committed
285
   timeloop.add() << BeforeFunction( communication, "communication" )
Cameron Stewart's avatar
Cameron Stewart committed
286
287
                  << Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingId ), "boundary handling" );
   timeloop.add() << BeforeFunction( lbm::viscoelastic::Su<LatticeModel_T, BoundaryHandling_T>(blocks, forceFieldId, pdfFieldId, boundaryHandlingId, stressId, stressOldId, velocityId,
Cameron Stewart's avatar
Cameron Stewart committed
288
                                                                                               lambda_p, eta_p, period, true), "viscoelasticity")
Cameron Stewart's avatar
Cameron Stewart committed
289
290
291
292
                  << Sweep( ConstantForce<BoundaryHandling_T>(forceFieldId, boundaryHandlingId, force),"Poiseuille Force");
   timeloop.add() << Sweep( makeSharedSweep( lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldId, flagFieldId, Fluid_Flag ) ),
                            "LB stream & collide" )
                  << AfterFunction(makeSharedFunctor(testData), "test data");
Cameron Stewart's avatar
Cameron Stewart committed
293

Cameron Stewart's avatar
Cameron Stewart committed
294
   timeloop.run();
Cameron Stewart's avatar
Cameron Stewart committed
295
296


297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
   if(!shortrun)
   {
      // compare to reference data
      real_t errSteady = real_c(fabs(testData->getUSteady() - real_c(1.0))/real_c(1.0));
      real_t errMax = real_c(fabs(testData->getUMax() - uMax)/uMax);
      real_t errt0 = real_c(fabs(testData->getT0() - t0)/t0);
      real_t errt1 = real_c(fabs(testData->getT1() - t1)/t1);

      WALBERLA_LOG_RESULT("Steady State Velocity Error: " << errSteady );
      WALBERLA_LOG_RESULT("Maximum Velocity Error: " << errMax );
      WALBERLA_LOG_RESULT("Time of Maximum Error: " << errt0 );
      WALBERLA_LOG_RESULT("Decay Time Error: " << errt1 );

      // check that errors < 1%
      if (errSteady < 0.01 && errMax < 0.01 && errt0 < 0.01 && errt1 < 0.01){
         WALBERLA_LOG_RESULT("Success" );
         return EXIT_SUCCESS;
      }
      else {
         WALBERLA_LOG_RESULT("Failure" );
         return EXIT_FAILURE;
      }
   } else{
Cameron Stewart's avatar
Cameron Stewart committed
320
321
      return EXIT_SUCCESS;
   }
Cameron Stewart's avatar
Cameron Stewart committed
322
323

}
cameron-stewart's avatar
cameron-stewart committed
324

cameron-stewart's avatar
cameron-stewart committed
325