FluidizedBedMEM.cpp 77.5 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/>.
//
16
//! \file FluidizedBedMEM.cpp
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
//! \author Dominik Schuster <dominik.schuster@fau.de>
//
//======================================================================================================================

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

#include "boundary/all.h"

#include "core/DataTypes.h"
#include "core/Environment.h"
#include "core/SharedFunctor.h"
#include "core/debug/Debug.h"
#include "core/debug/TestSubsystem.h"
#include "core/math/all.h"
#include "core/math/Random.h"
#include "core/timing/RemainingTimeLogger.h"

#include "domain_decomposition/SharedSweep.h"

#include "field/AddToStorage.h"
#include "field/communication/PackInfo.h"
39
#include "field/StabilityChecker.h"
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

#include "lbm/boundary/all.h"
#include "lbm/communication/PdfFieldPackInfo.h"
#include "lbm/field/AddToStorage.h"
#include "lbm/field/PdfField.h"
#include "lbm/lattice_model/D3Q19.h"
#include "lbm/sweeps/CellwiseSweep.h"
#include "lbm/sweeps/SplitPureSweep.h"
#include "lbm/lattice_model/SmagorinskyLES.h"
#include "lbm/PerformanceLogger.h"

#include "stencil/D3Q27.h"

#include "timeloop/SweepTimeloop.h"

#include "pe/basic.h"
#include "pe/cr/ICR.h"
#include "pe/Types.h"
#include "pe/vtk/BodyVtkOutput.h"
#include "pe/vtk/SphereVtkOutput.h"

#include "pe_coupling/mapping/all.h"
#include "pe_coupling/momentum_exchange_method/all.h"
#include "pe_coupling/utility/all.h"

#include <random>

#include "vtk/all.h"
#include "field/vtk/all.h"
#include "lbm/vtk/all.h"

#include "FluidizedBedMEMFunctionality.h"

namespace fluidizedBed {

///////////
// USING //
///////////

    using namespace walberla;
    using walberla::uint_t;

    const uint_t FieldGhostLayers = 1;

    typedef GhostLayerField<real_t, FieldGhostLayers> ScalarField_T;
    typedef GhostLayerField<pe::BodyID, FieldGhostLayers> BodyField_T;

    typedef walberla::uint8_t flag_t;
    typedef FlagField<flag_t> FlagField_T;

    typedef boost::tuple<pe::Sphere, pe::Plane> BodyTypeTuple;


/////////////////////////////////////
// BOUNDARY HANDLING CUSTOMIZATION //
/////////////////////////////////////

// class MyBoundaryHandling is responsible for creating the boundary handling and setting up the geometry at the outer domain boundaries

    template <typename BoundaryHandling_T, typename FlagField_T, typename PdfField_T, typename NoSlip_T, typename SimplePressure_T, typename UBB_T, typename MO_T>
    class MyBoundaryHandling {
    public:

        MyBoundaryHandling(const BlockDataID &flagField, const BlockDataID &pdfField, const BlockDataID &bodyField,
104
105
106
107
108
109
110
111
112
113
114
115
                //const shared_ptr<Vector3<real_t>> velocity,
                           const bool xPeriodic, const bool zPeriodic,
                           const bool spotInflow, const real_t spotDiameter, const real_t xCenter, const real_t  zCenter,
                           const real_t velocityMax, const real_t initialRamping, const real_t rampingTime,
                           const shared_ptr< lbm::TimeTracker > & timeTracker):
                flagField_(flagField), pdfField_(pdfField), bodyField_(bodyField)
                ,xPeriodic_(xPeriodic), zPeriodic_(zPeriodic), spotInflow_(spotInflow), spotDiameter_(spotDiameter)
                ,xCenter_(xCenter), zCenter_(zCenter), velocityMax_(velocityMax), initialRamping_(initialRamping)
                ,rampingTime_(rampingTime), timeTracker_(timeTracker)
                ,velocityFunctor_(FBfunc::VelocityFunctor_T (xCenter, zCenter, spotDiameter, velocityMax, initialRamping, rampingTime, spotInflow))
        {
        }
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

        BoundaryHandling_T *operator()(IBlock *const block, const StructuredBlockStorage *const storage) const;

    private:

        const BlockDataID flagField_;
        const BlockDataID pdfField_;
        const BlockDataID bodyField_;
        const bool xPeriodic_;
        const bool zPeriodic_;
        const bool spotInflow_;
        const real_t spotDiameter_;
        const real_t xCenter_;
        const real_t zCenter_;
        const real_t velocityMax_;
        const real_t initialRamping_;
        const real_t rampingTime_;
        const shared_ptr< lbm::TimeTracker > timeTracker_;
        FBfunc::VelocityFunctor_T velocityFunctor_;

    }; // class MyBoundaryHandling


    template <typename BoundaryHandling_T, typename FlagField_T, typename PdfField_T, typename NoSlip_T, typename SimplePressure_T, typename UBB_T, typename MO_T>
    BoundaryHandling_T *MyBoundaryHandling<BoundaryHandling_T, FlagField_T, PdfField_T, NoSlip_T, SimplePressure_T, UBB_T, MO_T>::operator()(IBlock *const block, const StructuredBlockStorage *const storage) const {
141
142
        WALBERLA_ASSERT_NOT_NULLPTR(block);
        WALBERLA_ASSERT_NOT_NULLPTR(storage);
143

144
145
146
        FlagField_T *flagField = block->getData<FlagField_T>(flagField_);
        PdfField_T *pdfField = block->getData<PdfField_T>(pdfField_);
        BodyField_T *bodyField = block->getData<BodyField_T>(bodyField_);
147

148
        const auto fluid = flagField->flagExists(FBfunc::Fluid_Flag) ? flagField->getFlag(FBfunc::Fluid_Flag) : flagField->registerFlag(FBfunc::Fluid_Flag);
149

150
151
152
153
154
        BoundaryHandling_T *handling = new BoundaryHandling_T("cf boundary handling", flagField, fluid,
                                                              boost::tuples::make_tuple(NoSlip_T("NoSlip", FBfunc::NoSlip_Flag, pdfField),
                                                                                        SimplePressure_T("Outlet", FBfunc::Outlet_Flag, pdfField, real_c(1.0)),
                                                                                        UBB_T("UBB", FBfunc::UBB_Flag, pdfField, timeTracker_, uint_t(0), velocityFunctor_, block->getAABB() ),
                                                                                        MO_T("MO", FBfunc::MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block)));
155

156
157
158
        const auto ubb = flagField->getFlag(FBfunc::UBB_Flag);
        const auto noSlip = flagField->getFlag(FBfunc::NoSlip_Flag);
        const auto outlet = flagField->getFlag(FBfunc::Outlet_Flag);
159

160
        CellInterval domainBB = storage->getDomainCellBB();
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
        domainBB.xMin() -= cell_idx_c(FieldGhostLayers);
        domainBB.xMax() += cell_idx_c(FieldGhostLayers);

        if (!xPeriodic_) {
// LEFT
            CellInterval left(domainBB.xMin(), domainBB.yMin(), domainBB.zMin(), domainBB.xMin(), domainBB.yMax(), domainBB.zMax());
            storage->transformGlobalToBlockLocalCellInterval(left, *block);
            handling->forceBoundary(noSlip, left);

// RIGHT
            CellInterval right(domainBB.xMax(), domainBB.yMin(), domainBB.zMin(), domainBB.xMax(), domainBB.yMax(), domainBB.zMax());
            storage->transformGlobalToBlockLocalCellInterval(right, *block);
            handling->forceBoundary(noSlip, right);
        }

        domainBB.zMin() -= cell_idx_c(FieldGhostLayers);
        domainBB.zMax() += cell_idx_c(FieldGhostLayers);

        if (!zPeriodic_) {
// FRONT
            CellInterval front(domainBB.xMin(), domainBB.yMin(), domainBB.zMin(), domainBB.xMax(), domainBB.yMax(), domainBB.zMin());
            storage->transformGlobalToBlockLocalCellInterval(front, *block);
            handling->forceBoundary(noSlip, front);

// BACK
            CellInterval back(domainBB.xMin(), domainBB.yMin(), domainBB.zMax(), domainBB.xMax(), domainBB.yMax(), domainBB.zMax());
            storage->transformGlobalToBlockLocalCellInterval(back, *block);
            handling->forceBoundary(noSlip, back);
        }


        domainBB.yMin() -= cell_idx_c(FieldGhostLayers);
        domainBB.yMax() += cell_idx_c(FieldGhostLayers);
195
196
197

// BOTTOM

198
        if (spotInflow_){
199

200
201
202
            CellInterval bottom(domainBB.xMin(), domainBB.yMin(), domainBB.zMin(), domainBB.xMax(), domainBB.yMin(), domainBB.zMax());
            storage->transformGlobalToBlockLocalCellInterval(bottom, *block);
            handling->forceBoundary(noSlip, bottom);
203

204
205
206
207
208
209
210
211
            for (cell_idx_t xIterator = domainBB.xMin(); xIterator <= domainBB.xMax(); ++xIterator ){
                for (cell_idx_t zIterator = domainBB.zMin(); zIterator <= domainBB.zMax(); ++zIterator ){
                    if ( sqrt ( (real_c(xIterator) - xCenter_)*(real_c(xIterator) - xCenter_) + (real_c(zIterator) - zCenter_)*(real_c(zIterator) - zCenter_) ) <= spotDiameter_ * 0.5 ){
                        CellInterval inflow (xIterator, domainBB.yMin(), zIterator, xIterator, domainBB.yMin(), zIterator);
                        storage->transformGlobalToBlockLocalCellInterval( inflow , *block);
                        handling->forceBoundary(ubb, inflow);
                    }
                }
212

213
            }
214
215


216
        } else {
217

218
219
220
            CellInterval bottom(domainBB.xMin(), domainBB.yMin(), domainBB.zMin(), domainBB.xMax(), domainBB.yMin(), domainBB.zMax());
            storage->transformGlobalToBlockLocalCellInterval(bottom, *block);
            handling->forceBoundary(ubb, bottom);
221

222
        }
223
224
225
226



// TOP
227
228
229
        CellInterval top(domainBB.xMin(), domainBB.yMax(), domainBB.zMin(), domainBB.xMax(), domainBB.yMax(), domainBB.zMax());
        storage->transformGlobalToBlockLocalCellInterval(top, *block);
        handling->forceBoundary(outlet, top);
230

231
        handling->fillWithDomain(FieldGhostLayers);
232

233
        return handling;
234
235
236
237
238
239
240
241
242
243
244
245
246
    }

    template <typename CollisionModel_T>
    struct LatticeModelCreator
    {

    };

    template<>
    struct LatticeModelCreator< lbm::collision_model::TRT >
    {
        static lbm::D3Q19<lbm::collision_model::TRT, false> makeLatticeModel(BlockDataID /*omegaFieldID*/, real_t omega)
        {
247
            return lbm::D3Q19<lbm::collision_model::TRT, false> (lbm::collision_model::TRT::constructWithMagicNumber(omega));
248
249
250
251
252
253
254
255
256
        }
    };


    template<>
    struct LatticeModelCreator< lbm::collision_model::SRT >
    {
        static lbm::D3Q19<lbm::collision_model::SRT, false> makeLatticeModel(BlockDataID /*omegaFieldID*/, real_t omega)
        {
257
            return lbm::D3Q19<lbm::collision_model::SRT, false> (lbm::collision_model::SRT(omega));
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
        }
    };

    template<>
    struct LatticeModelCreator< lbm::collision_model::SRTField<ScalarField_T> >
    {

        static lbm::D3Q19<lbm::collision_model::SRTField<ScalarField_T>, false> makeLatticeModel(BlockDataID omegaFieldID, real_t /*omega*/)
        {
            return lbm::D3Q19<lbm::collision_model::SRTField<ScalarField_T>, false> (lbm::collision_model::SRTField<ScalarField_T>(omegaFieldID));
        }
    };


    template <typename  collisionModel_T>
273
    int runSimulation ( Environment & env) {
274

275
        uint_t processes = MPIManager::instance()->numProcesses();
276

277
        FBfunc::SetupFB config;
278

279
        auto model_parameters = env.config()->getBlock("Model");
280

281
        const bool turb = model_parameters.getParameter<bool>("turbulence_model", false);
282
283
284
285
286

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

287
288
289
        typedef lbm::D3Q19<collisionModel_T, false> LatticeModel_T;
        typedef typename LatticeModel_T::Stencil Stencil_T;
        typedef lbm::PdfField<LatticeModel_T> PdfField_T;
290
291
292

// boundary handling

293
294
295
296
        typedef lbm::NoSlip<LatticeModel_T, flag_t> NoSlip_T;
        typedef lbm::SimplePressure<LatticeModel_T, flag_t> SimplePressure_T;
        typedef lbm::DynamicUBB<LatticeModel_T, flag_t, FBfunc::VelocityFunctor_T> UBB_T;
        typedef pe_coupling::CurvedLinear<LatticeModel_T, FlagField_T> MO_T;
297

298
299
        typedef boost::tuples::tuple<NoSlip_T, SimplePressure_T, UBB_T, MO_T> BoundaryConditions_T;
        typedef BoundaryHandling<FlagField_T, Stencil_T, BoundaryConditions_T> BoundaryHandling_T;
300

301
        typedef MyBoundaryHandling<BoundaryHandling_T, FlagField_T, PdfField_T, NoSlip_T, SimplePressure_T, UBB_T, MO_T> MyBoundaryHandling_T;
302
303
304
305
306
307


///////////////////////////
// SIMULATION PROPERTIES //
///////////////////////////

308
309
        // GEOMETRY PARAMETERS
        auto geo_parameters = env.config()->getBlock("Geometry");
310

311
312
313
314
315
        const real_t depth_SI = geo_parameters.getParameter<real_t>("depth", real_c(0.02));                         //[m]
        const real_t width_SI = geo_parameters.getParameter<real_t>("width", real_c(0.10));                         //[m]
        const real_t height_SI = geo_parameters.getParameter<real_t>("height", real_c(0.40));                         //[m]
        const real_t offsetTop_SI = geo_parameters.getParameter<real_t>("offset_top", real_c(0.00));                         //[m]
        const real_t offsetBot_SI = geo_parameters.getParameter<real_t>("offset_bot", real_c(0.00));                         //[m]
316

317
318
319
320
321
322
323
324
325
        if (!geo_parameters.isDefined("depth")) {
            WALBERLA_ABORT("You need to specify \"depth\" in the configuration file");
        }
        if (!geo_parameters.isDefined("width")) {
            WALBERLA_ABORT("You need to specify \"width\" in the configuration file");
        }
        if (!geo_parameters.isDefined("height")) {
            WALBERLA_ABORT("You need to specify \"height\" in the configuration file");
        }
326

327
        const real_t crossSectionArea_SI = width_SI * depth_SI;
328

329
330
        const bool xPeriodic = geo_parameters.getParameter<bool>("width_periodic", false);
        const bool zPeriodic = geo_parameters.getParameter<bool>("depth_periodic", false);
331

332
333
        const bool twoSpecies = geo_parameters.getParameter<bool>("two_species", false);
        //const bool createSpeciesWithRandomGenerator = geo_parameters.getParameter<bool>("create_species_with_random_generator", false);
334

335
336
337
338
339
340
341
342
343
        const real_t diameterA_SI = geo_parameters.getParameter<real_t>("diameter_A", real_c(0.00215));   //[m]
        const real_t diameterB_SI = geo_parameters.getParameter<real_t>("diameter_B", real_c(0.00215));   //[m]
        const real_t radiusStandardDeviationA = 0.5 * geo_parameters.getParameter<real_t>("diameter_standard_distribution_A", real_c(0.00000));
        const real_t radiusStandardDeviationB = 0.5 * geo_parameters.getParameter<real_t>("diameter_standard_distribution_B", real_c(0.00000));
        config.dx_SI = geo_parameters.getParameter<real_t>("cell_size", real_c(0.0002));
        config.diameterA = diameterA_SI / config.dx_SI;
        config.diameterB = diameterB_SI / config.dx_SI;
        const real_t radiusA = real_c(config.diameterA * 0.5);
        const real_t radiusB = real_c(config.diameterB * 0.5);
344

345
346
347
        const uint_t xBlocks = geo_parameters.getParameter<uint_t>("blocks_width", uint_c(5));
        const uint_t yBlocks = geo_parameters.getParameter<uint_t>("blocks_height", uint_c(20));
        const uint_t zBlocks = geo_parameters.getParameter<uint_t>("blocks_depth", uint_c(1));
348

349
350
351
        if (processes != (xBlocks * yBlocks * zBlocks)) {
            WALBERLA_ABORT("number of processes must be equal to number of blocks ( here: " << (xBlocks * yBlocks * zBlocks) << " )");
        }
352

353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
        config.xlength = uint_c(std::round(width_SI / config.dx_SI));
        config.ylength = uint_c(std::round((height_SI + offsetTop_SI + offsetBot_SI) / config.dx_SI));
        config.zlength = uint_c(std::round(depth_SI / config.dx_SI));
        config.offsetTop = offsetBot_SI / config.dx_SI;
        config.offsetBot = offsetBot_SI / config.dx_SI;

        const bool spotInflow = geo_parameters.getParameter<bool>("spot_inflow", false);
        real_t spotDiameter = geo_parameters.getParameter<real_t>("spot_diameter", real_c(0.01));
        spotDiameter = spotDiameter / config.dx_SI;

        // FLUIDIZED BED PARAMETERS
        auto fb_parameters = env.config()->getBlock("FluidizedBed");

        config.nrParticlesA = fb_parameters.getParameter<uint_t>("#particles_A", uint_c(0));
        config.nrParticlesB = fb_parameters.getParameter<uint_t>("#particles_B", uint_c(0));
        const real_t densitySolidA_SI = fb_parameters.getParameter<real_t>("density_particle_A",
                                                                           real_c(22.5));         //[kg/m3] // alumina katalyst
        const real_t densitySolidB_SI = fb_parameters.getParameter<real_t>("density_particle_B", real_c(22.5));         //[kg/m3]
        config.densityFluid_SI = fb_parameters.getParameter<real_t>("density_fluid", real_c(1.2));          //[kg/m3]      //water
        const real_t dynViscosity_SI = fb_parameters.getParameter<real_t>("dynamic_viscosity",
                                                                          real_c(1.84e-5));       //[Pa*s]  // viscosity of water
        const real_t kinViscosity_SI = dynViscosity_SI / config.densityFluid_SI;                                       //[m2/s]
        const real_t gravity_SI = fb_parameters.getParameter<real_t>("gravity", real_c(9.81));

        config.densityRatioA = densitySolidA_SI / config.densityFluid_SI;
        config.densityRatioB = densitySolidB_SI / config.densityFluid_SI;

        real_t velocity_SI;
        real_t particleReynoldsNumber_SI;
        if (fb_parameters.isDefined("velocity")) {
            velocity_SI = fb_parameters.getParameter<real_t>("velocity", real_c(0.35));
            particleReynoldsNumber_SI = velocity_SI * diameterA_SI / kinViscosity_SI;
        } else if (fb_parameters.isDefined("particle_RE")) {
            particleReynoldsNumber_SI = fb_parameters.getParameter<real_t>("particle_RE", real_c(48));
            velocity_SI = particleReynoldsNumber_SI * kinViscosity_SI / diameterA_SI;
        } else if (fb_parameters.isDefined("volume_flow_rate")) {
            real_t volumeFlowRate_SI = fb_parameters.getParameter<real_t>("volume_flow_rate", real_c(2.5));
            velocity_SI = real_c(volumeFlowRate_SI / crossSectionArea_SI / 3600);
            particleReynoldsNumber_SI = velocity_SI * diameterA_SI / kinViscosity_SI;
        } else {
            WALBERLA_ABORT("You need to specify either \"volume_flow_rate\" or \"velocity\" or \"particle_RE\" in the configuration file");
        }
395

396
397
        // SIMULATION PARAMETERS
        auto sim_parameters = env.config()->getBlock("Simulation");
398

399
        const bool useLubrication = sim_parameters.getParameter<bool>("use_lubrication", false);
400
401
402
403
404
405
406

        const bool createPackedBed = sim_parameters.getParameter<bool>("create_packed_bed", false);
        const bool checkMaximumParticleVelocity = sim_parameters.getParameter<bool>("check_maximum_particle_velocity", true);
        const real_t settleVelocity = sim_parameters.getParameter<real_t>("settle_velocity", real_c(0.002));
        const uint_t settleCheckFrequency = sim_parameters.getParameter<uint_t>("settle_check_frequency", uint_t(100));
        const real_t peAccelerator = sim_parameters.getParameter<real_t>("pe_accelerator", real_c(50));

407
        const bool runPreliminaryTimeloop = sim_parameters.getParameter<bool>("run_preliminary_timeloop", false);
408
409
        const real_t eps = sim_parameters.getParameter<real_t>("eps", real_c(0.05));

410
411
412
413
        const bool initializeFromCheckPointFile = sim_parameters.getParameter<bool>("initialize_from_checkpoint", false);
        const bool enableCheckpointing = sim_parameters.getParameter<bool>("enable_checkpointing", false);
        const std::string checkpointPathBodies = sim_parameters.getParameter<std::string>("path_to_checkpoint_bodies", "checkpoint_bodies.dat");
        const std::string checkpointPathPdf = sim_parameters.getParameter<std::string>("path_to_checkpoint_pdf", "checkpoint_pdf.dat");
414

415
416
417
418
        config.evaluationFreq = sim_parameters.getParameter<uint_t>("evaluation_frequency", uint_c(100));
        config.velCheckFreq = sim_parameters.getParameter<uint_t>("velocity_check_frequency", uint_c(10));
        config.checkpointFreq = sim_parameters.getParameter<uint_t>("checkpoint_frequency", uint_c(100000)); //has to be large
        config.stopVel = sim_parameters.getParameter<real_t>("abort_LBM_velocity", real_c(0.50));
419
420
421
        const bool useForceAveraging = sim_parameters.getParameter<bool>("use_force_averaging", true);
        config.lbmSubCycles = useForceAveraging ? real_t(2) : real_t(1);

422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
        if (sim_parameters.isDefined("time_step")) {
            config.dt_SI = sim_parameters.getParameter<real_t>("time_step", real_c(1e-6));
            config.velocity = velocity_SI * config.dt_SI / config.dx_SI;
        } else if (sim_parameters.isDefined("LBM_velocity")) {
            config.velocity = sim_parameters.getParameter("LBM_velocity", real_c(0.06));
            config.dt_SI = config.velocity * config.dx_SI / velocity_SI;
        } else if (sim_parameters.isDefined("omega")) {
            real_t omegaTemp = sim_parameters.getParameter<real_t>("omega", real_c(1.90));
            config.kinViscosity = lbm::collision_model::viscosityFromOmega(omegaTemp);
            config.dt_SI = config.kinViscosity / kinViscosity_SI * config.dx_SI * config.dx_SI;
            config.velocity = velocity_SI * config.dt_SI / config.dx_SI;
        } else if (sim_parameters.isDefined("tau")) {
            real_t tauTemp = sim_parameters.getParameter<real_t>("tau", real_c(0.52));
            config.kinViscosity = lbm::collision_model::viscosityFromOmega(real_c(1.0) / tauTemp);
            config.dt_SI = config.kinViscosity / kinViscosity_SI * config.dx_SI * config.dx_SI;
            config.velocity = velocity_SI * config.dt_SI / config.dx_SI;
        } else {
            WALBERLA_ABORT("You need to specify either \"time_step\" or \"LBM_velocity\" or \"omega\" or \"tau\" in the configuration file");
        }
441
442


443
        const real_t smagorinskyConstant = sim_parameters.getParameter<real_t>("smagorinsky_constant", real_c(0.1));
444

445
446
        uint_t timesteps;
        uint_t timestepsRamping;
447

448
449
450
451
452
453
454
455
        if (sim_parameters.isDefined("simulation_time")) {
            real_t simTime = sim_parameters.getParameter<real_t>("simulation_time", real_c(1.0));
            timesteps = uint_c(simTime / (config.lbmSubCycles * config.dt_SI));
        } else if (sim_parameters.isDefined("#time_steps")) {
            timesteps = sim_parameters.getParameter<uint_t>("#time_steps", uint_c(100000));
        } else {
            WALBERLA_ABORT("You need to specify either \"simulation_time\" or \"#time_steps\"  in the configuration file");
        }
456

457
458
459
460
461
462
        if (sim_parameters.isDefined("ramping_time")) {
            real_t rampingTime = sim_parameters.getParameter<real_t>("ramping_time", real_c(0.0));
            timestepsRamping = uint_c(rampingTime / (config.lbmSubCycles * config.dt_SI));
        } else {
            timestepsRamping = sim_parameters.getParameter<uint_t>("#ramping_steps", uint_c(0));
        }
463

464
        const real_t initialRampingVelocityFactor = sim_parameters.getParameter<real_t>("initial_ramping_scaling", real_c(0.50));
465

466
467
468
469
        const real_t initialSphereVel =
                geo_parameters.getParameter<real_t>("initial_sphere_random_velocity", real_c(0.0)) * (config.dt_SI / config.dx_SI);
        const real_t initialSphereDistance =
                geo_parameters.getParameter<real_t>("initial_sphere_distance", real_c(diameterA_SI * 1.5)) / config.dx_SI;
470
471


472
473
        config.kinViscosity = kinViscosity_SI * config.dt_SI / config.dx_SI / config.dx_SI;
        config.gravity = gravity_SI * config.dt_SI * config.dt_SI / config.dx_SI;
474

475
476
477
478
        const real_t FroudeNumber_SI = real_c(sqrt(velocity_SI * velocity_SI / (gravity_SI * diameterA_SI) / (config.densityRatioA - 1)));
        const real_t ArrheniusNumber_SI =
                (densitySolidA_SI - config.densityFluid_SI) * config.densityFluid_SI * diameterA_SI * diameterA_SI * diameterA_SI * gravity_SI /
                (dynViscosity_SI * dynViscosity_SI);
479

480
481
482
483
484
        const real_t particleReynoldsNumber_LBM = config.velocity * config.diameterA / config.kinViscosity;
        const real_t FroudeNumber_LBM = real_c(
                sqrt(config.velocity * config.velocity / (config.gravity * config.diameterA) / (config.densityRatioA - 1)));
        const real_t galileiNumber =
                config.gravity * config.diameterA * config.diameterA * config.diameterA / (config.kinViscosity * config.kinViscosity);
485

486
487
488
489
        const real_t UmfFroudenumber = real_c(33.7 * (sqrt(1 + 3.6e-5 * ArrheniusNumber_SI) - 1) / sqrt(ArrheniusNumber_SI));
        const real_t REpmf = real_c(sqrt(33.7 * 33.7 + 0.0494 * ArrheniusNumber_SI) - 33.7);            // Wen You 1966
        const real_t gravityPressure = real_c(config.nrParticlesA) * gravity_SI * (densitySolidA_SI - config.densityFluid_SI)
                                       * (1.0 / 6.0 * math::PI * diameterA_SI * diameterA_SI * diameterA_SI) / crossSectionArea_SI;
490

491
        const real_t omega = lbm::collision_model::omegaFromViscosity(config.kinViscosity);
492

493
494
        // MATERIAL PARAMETERS
        auto mat_parameters = env.config()->getBlock("Material");
495

496
497
        const bool useDEM = mat_parameters.getParameter<bool>("use_DEM", true);

498
        const uint_t HCSITSMaxIterations = mat_parameters.getParameter<uint_t>("HCSITS_max_iterations", uint_c(10));
499
500
        const real_t HCSITSRelaxationParameter = mat_parameters.getParameter<real_t>("HCSITS_relaxation_parameter", real_c(0.7));

501
502
503
504
        const real_t volumeA = real_c(4.0 / 3.0 * math::PI * radiusA * radiusA * radiusA);
        const real_t volumeB = real_c(4.0 / 3.0 * math::PI * radiusB * radiusB * radiusB);
        const real_t MijA = real_c(0.5 * config.densityRatioA * volumeA);
        const real_t MijB = real_c(0.5 * config.densityRatioB * volumeB);
505
506

        real_t MijWall = MijA;
507
508
        if (twoSpecies) {
            MijWall = 0.5 * (MijA + MijB);
509
510
        }

511
512
        const real_t cdsA = mat_parameters.getParameter<real_t>("static_friction_particles_A", real_c(0.5));           // static friction
        const real_t cdsB = mat_parameters.getParameter<real_t>("static_friction_particles_B", real_c(0.5));           // static friction
513
        const real_t cdsWall = mat_parameters.getParameter<real_t>("static_friction_particles_wall", real_c(0.5));        // static friction
514
515
        const real_t cdfA = mat_parameters.getParameter<real_t>("dynamic_friction_particles_A", real_c(0.5));          // dynamic friction
        const real_t cdfB = mat_parameters.getParameter<real_t>("dynamic_friction_particles_B", real_c(0.5));          // dynamic friction
516
517
518
        const real_t cdfWall = mat_parameters.getParameter<real_t>("dynamic_friction_wall", real_c(0.5));                 // dynamic friction

        const real_t poisson = real_c(0.2);                  // not used in DEM
519
520
        const real_t youngsModulus = 100000;                 // not used in DEM
        //const real_t restitution = real_c(0.8);              // not used in DEM
521

522
523
524
        // found in: Interface-resolved direct numerical simulation of the erosion of a sediment bed sheared by laminar channel flow
        //           Kidanemariam and Uhlmann
        //           DOI: 10.1016/j.ijmultiphaseflow.2014.08.008
525

526
527
        //const real_t normalizedstiffness = real_c(13000);
        //const real_t stiffness = normalizedstiffness * ((densityRatio-1)*gravity*volume/diameter);
528

529
530
        // found in: Comparison between finite volume and lattice-Boltzmann method simulations of gas-fluidised beds: bed expansion and particle–fluid interaction force
        // R. Third, Y. Chen, C. R. Müller
531

532
533
534
535
536
537
538
539
540
541
542
543
        const real_t restitutionParticlesA = mat_parameters.getParameter<real_t>("restitution_particles_A", real_c(0.2));
        const real_t restitutionParticlesB = mat_parameters.getParameter<real_t>("restitution_particles_B", real_c(0.2));
        const real_t restitutionWall = mat_parameters.getParameter<real_t>("restitution_wall", real_c(0.2));
        const real_t collisionTimeParticlesA = mat_parameters.getParameter<real_t>("collision_time_particles_A", real_c(10));
        const real_t collisionTimeParticlesB = mat_parameters.getParameter<real_t>("collision_time_particles_B", real_c(10));
        const real_t collisionTimeWall = mat_parameters.getParameter<real_t>("collision_time_wall", real_c(10));
        const real_t tangentialDampingFactorParticlesA = mat_parameters.getParameter<real_t>("tangential_damping_factor_particles_A",
                                                                                             real_c(1.0));
        const real_t tangentialDampingFactorParticlesB = mat_parameters.getParameter<real_t>("tangential_damping_factor_particles_B",
                                                                                             real_c(1.0));
        const real_t tangentialDampingFactorWall = mat_parameters.getParameter<real_t>("tangential_damping_factor_wall", real_c(1.0));
        const real_t timestepsPerCollision = mat_parameters.getParameter<real_t>("timesteps_per_collision", real_c(100));
544

545
        const real_t minCollisionTime = std::min(collisionTimeParticlesA, std::min(collisionTimeParticlesB, collisionTimeWall));
546

547
        config.peSubCycles = uint_c(std::round(timestepsPerCollision / minCollisionTime));
548
549
550
551

        // put eq(14) in eq(16) and solve for k_n

        // PARTICLES A
552
553
554
555
556
557
558
        const real_t stiffnessNparticlesA = (MijA * math::PI * math::PI) /
                                            (collisionTimeParticlesA * collisionTimeParticlesA *
                                             (1.0 - (log(restitutionParticlesA) * log(restitutionParticlesA))
                                                    / (math::PI * math::PI + log(restitutionParticlesA) * log(restitutionParticlesA))));
        const real_t dampingNparticlesA = -2.0 * sqrt(MijA * stiffnessNparticlesA) * log(restitutionParticlesA) /
                                          sqrt(math::PI * math::PI + log(restitutionParticlesA) * log(restitutionParticlesA));
        const real_t dampingTparticlesA = dampingNparticlesA * tangentialDampingFactorParticlesA;
559
560

        // PARTICLES B
561
562
563
564
565
566
567
        const real_t stiffnessNparticlesB = (MijB * math::PI * math::PI) /
                                            (collisionTimeParticlesB * collisionTimeParticlesB *
                                             (1.0 - (log(restitutionParticlesB) * log(restitutionParticlesB))
                                                    / (math::PI * math::PI + log(restitutionParticlesB) * log(restitutionParticlesB))));
        const real_t dampingNparticlesB = -2.0 * sqrt(MijB * stiffnessNparticlesB) * log(restitutionParticlesB) /
                                          sqrt(math::PI * math::PI + log(restitutionParticlesB) * log(restitutionParticlesB));
        const real_t dampingTparticlesB = dampingNparticlesB * tangentialDampingFactorParticlesB;
568
569

        // WALL
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
        const real_t stiffnessNwall = (MijWall * math::PI * math::PI) /
                                      (collisionTimeWall * collisionTimeWall * (1.0 - (log(restitutionWall) * log(restitutionWall))
                                                                                      / (math::PI * math::PI +
                                                                                         log(restitutionWall) * log(restitutionWall))));
        const real_t dampingNwall = -2.0 * sqrt(MijWall * stiffnessNwall) * log(restitutionWall) /
                                    sqrt(math::PI * math::PI + log(restitutionWall) * log(restitutionWall));
        const real_t dampingTwall = dampingNwall * tangentialDampingFactorWall;

        // cdn=cdt

        // EVALUATION PARAMETERS
        auto eval_parameters = env.config()->getBlock("Evaluation");

        const bool calculateCenterOfMass = eval_parameters.getParameter<bool>("mass_center", false);
        const bool calculatePressureDrop = eval_parameters.getParameter<bool>("pressure_drop", false);
        const bool calculateErgunPressure = eval_parameters.getParameter<bool>("ergun_pressure_drop", false);
        const bool calculatePressureDifference = eval_parameters.getParameter<bool>("pressure_difference", false);
        const bool calculateParticleFlux = eval_parameters.getParameter<bool>("particle_flux", false);
        const bool calculateGranularTempGlobal = eval_parameters.getParameter<bool>("granular_temperature_global", false);
        const bool calculateGranularTempLocal = eval_parameters.getParameter<bool>("granular_temperature_local", false);
        const bool calculateSolidVolumeFraction = eval_parameters.getParameter<bool>("solid_volume_fraction", false);
        //const bool calculateCollisions = eval_parameters.getParameter<bool>("collision_count", false);

        const bool box1 = eval_parameters.getParameter<bool>("box1", false);

        real_t box1xSize = eval_parameters.getParameter<real_t>("box1_Size_x", real_t(0.02));
        box1xSize = box1xSize / config.dx_SI;
        real_t box1ySize = eval_parameters.getParameter<real_t>("box1_Size_y", real_t(0.02));
        box1ySize = box1ySize / config.dx_SI;
        real_t box1zSize = eval_parameters.getParameter<real_t>("box1_Size_z", real_t(0.02));
        box1zSize = box1zSize / config.dx_SI;

        real_t box1xCenter = eval_parameters.getParameter<real_t>("box1_Center_x", real_t(0.05));
        box1xCenter = box1xCenter / config.dx_SI;
        real_t box1yCenter = eval_parameters.getParameter<real_t>("box1_Center_y", real_t(0.05));
        box1yCenter = box1yCenter / config.dx_SI;
        real_t box1zCenter = eval_parameters.getParameter<real_t>("box1_Center_z", real_t(0.01));
        box1zCenter = box1zCenter / config.dx_SI;

        const bool box2 = eval_parameters.getParameter<bool>("box2", false);

        real_t box2xSize = eval_parameters.getParameter<real_t>("box2_Size_x", real_t(0.02));
        box2xSize = box2xSize / config.dx_SI;
        real_t box2ySize = eval_parameters.getParameter<real_t>("box2_Size_y", real_t(0.02));
        box2ySize = box2ySize / config.dx_SI;
        real_t box2zSize = eval_parameters.getParameter<real_t>("box2_Size_z", real_t(0.02));
        box2zSize = box2zSize / config.dx_SI;

        real_t box2xCenter = eval_parameters.getParameter<real_t>("box2_Center_x", real_t(0.05));
        box2xCenter = box2xCenter / config.dx_SI;
        real_t box2yCenter = eval_parameters.getParameter<real_t>("box2_Center_y", real_t(0.10));
        box2yCenter = box2yCenter / config.dx_SI;
        real_t box2zCenter = eval_parameters.getParameter<real_t>("box2_Center_z", real_t(0.01));
        box2zCenter = box2zCenter / config.dx_SI;

        const bool box3 = eval_parameters.getParameter<bool>("box3", false);

        real_t box3xSize = eval_parameters.getParameter<real_t>("box3_Size_x", real_t(0.02));
        box3xSize = box3xSize / config.dx_SI;
        real_t box3ySize = eval_parameters.getParameter<real_t>("box3_Size_y", real_t(0.02));
        box3ySize = box3ySize / config.dx_SI;
        real_t box3zSize = eval_parameters.getParameter<real_t>("box3_Size_z", real_t(0.02));
        box3zSize = box3zSize / config.dx_SI;

        real_t box3xCenter = eval_parameters.getParameter<real_t>("box3_Center_x", real_t(0.05));
        box3xCenter = box3xCenter / config.dx_SI;
        real_t box3yCenter = eval_parameters.getParameter<real_t>("box3_Center_y", real_t(0.15));
        box3yCenter = box3yCenter / config.dx_SI;
        real_t box3zCenter = eval_parameters.getParameter<real_t>("box3_Center_z", real_t(0.01));
        box3zCenter = box3zCenter / config.dx_SI;

        // VTK PARAMETERS
        auto vtk_parameters = env.config()->getBlock("VTK");

        const bool writePackedBed = vtk_parameters.getParameter<bool>("write_packed_bed", false);
        const uint_t frequencyPackedBed = vtk_parameters.getParameter<uint_t>("frequency_packed_bed", uint_c(500));

        const bool writeSpheres = vtk_parameters.getParameter<bool>("write_spheres", false);
        const uint_t frequencySpheres = vtk_parameters.getParameter<uint_t>("frequency_spheres", uint_c(200));

        const bool writeFlagField = vtk_parameters.getParameter<bool>("write_flag_field", false);
        const uint_t frequencyFlagField = vtk_parameters.getParameter<uint_t>("frequency_flag_field", uint_c(200000));

        const bool writeFluidField = vtk_parameters.getParameter<bool>("write_fluid_field", false);
        const uint_t frequencyFluidField = vtk_parameters.getParameter<uint_t>("frequency_fluid_field", uint_c(50000));
        const real_t samplingResolutionFluidField = vtk_parameters.getParameter<real_t>("sampling_resolution_fluid_field", real_c(2.0));

        const bool writeOmega = vtk_parameters.getParameter<bool>("write_omega", false);
        const uint_t frequencyOmega = vtk_parameters.getParameter<uint_t>("frequency_omega", uint_c(50000));
        const real_t samplingResolutionOmega = vtk_parameters.getParameter<real_t>("sampling_resolution_omega", real_c(2.0));


        ///////////////////////////
        // BLOCK STRUCTURE SETUP //
        ///////////////////////////

        const uint_t xCells = config.xlength / xBlocks;
        const uint_t yCells = config.ylength / yBlocks;
        const uint_t zCells = config.zlength / zBlocks;

        // ASSERTS
        WALBERLA_CHECK_EQUAL(config.xlength, xCells * xBlocks);
        WALBERLA_CHECK_EQUAL(config.ylength, yCells * yBlocks);
        WALBERLA_CHECK_EQUAL(config.zlength, zCells * zBlocks);

        const real_t dx = real_t(1);
        const real_t overlap = real_c(1.5) * dx;

        auto blocks = blockforest::createUniformBlockGrid(xBlocks, yBlocks, zBlocks, xCells, yCells, zCells, dx, true,
                                                          xPeriodic, false, zPeriodic);

        // LOCAL VOLUME FOR SOLID FRACTION CALCULATION AND GRANULAR TEMPERATURE LOCALLY
        const AABB probe1(
                (blocks->getDomain().xMin() < (box1xCenter - box1xSize / 2.0)) ? (box1xCenter - box1xSize / 2.0) : blocks->getDomain().xMin(),
                (blocks->getDomain().yMin() < (box1yCenter - box1ySize / 2.0)) ? (box1yCenter - box1ySize / 2.0) : blocks->getDomain().yMin(),
                (blocks->getDomain().zMin() < (box1zCenter - box1zSize / 2.0)) ? (box1zCenter - box1zSize / 2.0) : blocks->getDomain().zMin(),
                (blocks->getDomain().xMax() > (box1xCenter + box1xSize / 2.0)) ? (box1xCenter + box1xSize / 2.0) : blocks->getDomain().xMax(),
                (blocks->getDomain().yMax() > (box1yCenter + box1ySize / 2.0)) ? (box1yCenter + box1ySize / 2.0) : blocks->getDomain().yMax(),
                (blocks->getDomain().zMax() > (box1zCenter + box1zSize / 2.0)) ? (box1zCenter + box1zSize / 2.0) : blocks->getDomain().zMax());

        const AABB probe2(
                (blocks->getDomain().xMin() < (box2xCenter - box2xSize / 2.0)) ? (box2xCenter - box2xSize / 2.0) : blocks->getDomain().xMin(),
                (blocks->getDomain().yMin() < (box2yCenter - box2ySize / 2.0)) ? (box2yCenter - box2ySize / 2.0) : blocks->getDomain().yMin(),
                (blocks->getDomain().zMin() < (box2zCenter - box2zSize / 2.0)) ? (box2zCenter - box2zSize / 2.0) : blocks->getDomain().zMin(),
                (blocks->getDomain().xMax() > (box2xCenter + box2xSize / 2.0)) ? (box2xCenter + box2xSize / 2.0) : blocks->getDomain().xMax(),
                (blocks->getDomain().yMax() > (box2yCenter + box2ySize / 2.0)) ? (box2yCenter + box2ySize / 2.0) : blocks->getDomain().yMax(),
                (blocks->getDomain().zMax() > (box2zCenter + box2zSize / 2.0)) ? (box2zCenter + box2zSize / 2.0) : blocks->getDomain().zMax());

        const AABB probe3(
                (blocks->getDomain().xMin() < (box3xCenter - box3xSize / 2.0)) ? (box3xCenter - box3xSize / 2.0) : blocks->getDomain().xMin(),
                (blocks->getDomain().yMin() < (box3yCenter - box3ySize / 2.0)) ? (box3yCenter - box3ySize / 2.0) : blocks->getDomain().yMin(),
                (blocks->getDomain().zMin() < (box3zCenter - box3zSize / 2.0)) ? (box3zCenter - box3zSize / 2.0) : blocks->getDomain().zMin(),
                (blocks->getDomain().xMax() > (box3xCenter + box3xSize / 2.0)) ? (box3xCenter + box3xSize / 2.0) : blocks->getDomain().xMax(),
                (blocks->getDomain().yMax() > (box3yCenter + box3ySize / 2.0)) ? (box3yCenter + box3ySize / 2.0) : blocks->getDomain().yMax(),
                (blocks->getDomain().zMax() > (box3zCenter + box3zSize / 2.0)) ? (box3zCenter + box3zSize / 2.0) : blocks->getDomain().zMax());

        // Transformation of AABB to CellInterval for evaluation functions
        if (box1) {
            config.probesAABB.push_back(probe1);
            config.probesCell.push_back(blocks->getCellBBFromAABB(probe1));
        }
        if (box2) {
            config.probesAABB.push_back(probe2);
            config.probesCell.push_back(blocks->getCellBBFromAABB(probe2));
        }
        if (box3) {
            config.probesAABB.push_back(probe3);
            config.probesCell.push_back(blocks->getCellBBFromAABB(probe3));
        }
719

720
721
        const real_t yLowerPlane = config.offsetBot <= 0.0 ? 0.5 : config.offsetBot - 0.5;
        const real_t yUpperPlane = config.offsetTop <= 0.0 ? real_c(config.ylength) - 0.5 : real_c(config.ylength) - config.offsetTop + 0.5;
722

723
724
725
726
        const AABB lowerPlaneAABB(blocks->getDomain().xMin(), yLowerPlane, blocks->getDomain().zMin(),
                                  blocks->getDomain().xMax(), yLowerPlane, blocks->getDomain().zMax());
        const AABB upperPlaneAABB(blocks->getDomain().xMin(), yUpperPlane, blocks->getDomain().zMin(),
                                  blocks->getDomain().xMax(), yUpperPlane, blocks->getDomain().zMax());
727

728
729
        const CellInterval lowerCellPlane = blocks->getCellBBFromAABB(lowerPlaneAABB);
        const CellInterval upperCellPlane = blocks->getCellBBFromAABB(upperPlaneAABB);
730

731
732
733
        /////////////////
        // PE COUPLING //
        /////////////////
734

735
736
        shared_ptr<pe::BodyStorage> globalBodyStorage = make_shared<pe::BodyStorage>();
        pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
737

738
        auto bodyStorageID = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "pe Body Storage");
739

740
741
742
743
744
745
746
747
748
        const auto materialA = pe::createMaterial("myMatA", config.densityRatioA,
                                                  restitutionParticlesA, cdsA, cdfA, poisson, youngsModulus, stiffnessNparticlesA,
                                                  dampingNparticlesA, dampingTparticlesA);
        const auto materialB = pe::createMaterial("myMatB", config.densityRatioB,
                                                  restitutionParticlesB, cdsB, cdfB, poisson, youngsModulus, stiffnessNparticlesB,
                                                  dampingNparticlesA, dampingTparticlesB);
        const auto materialWall = pe::createMaterial("myMatWall", 0.5 * (config.densityRatioA + config.densityRatioB),
                                                     restitutionWall, cdsWall, cdfWall, poisson, youngsModulus, stiffnessNwall, dampingNwall,
                                                     dampingTwall);
749
750

        // fixed geometry
751
752
753
754
755
756
        if (!zPeriodic) {
            pe::createPlane(*globalBodyStorage, 0, Vector3<real_t>(0, 0, 1), Vector3<real_t>(0, 0, 0),
                            materialWall);                                                           // front
            pe::createPlane(*globalBodyStorage, 0, Vector3<real_t>(0, 0, -1), Vector3<real_t>(0, 0, real_c(config.zlength)),
                            materialWall);                                     // back
        }
757

758
759
760
761
762
763
        if (!xPeriodic) {
            pe::createPlane(*globalBodyStorage, 0, Vector3<real_t>(1, 0, 0), Vector3<real_t>(0, 0, 0),
                            materialWall);                                                           // left
            pe::createPlane(*globalBodyStorage, 0, Vector3<real_t>(-1, 0, 0), Vector3<real_t>(real_c(config.xlength), 0, 0),
                            materialWall);                                     // right
        }
764

765
766
767
768
        pe::createPlane(*globalBodyStorage, 0, Vector3<real_t>(0, 1, 0), Vector3<real_t>(0, config.offsetBot, 0),
                        materialWall);                            // bottom
        pe::createPlane(*globalBodyStorage, 0, Vector3<real_t>(0, -1, 0), Vector3<real_t>(0, real_c(config.ylength) - config.offsetTop, 0),
                        materialWall);                            // top
769

770
        vtk::writeDomainDecomposition(blocks, "domain_decomposition");
771

772
773
774
        ///////////////////
        // PE SIMULATION //
        ///////////////////
775

776
777
        //create empty pdf field
        BlockDataID pdfFieldID;
778

779
780
        // create omega  field for turbulence model
        BlockDataID omegaFieldID = field::addToStorage<ScalarField_T>(blocks, "omega field", omega, field::fzyx, FieldGhostLayers);
781

782
783
        // create appropriate lattice model
        LatticeModel_T latticeModel = LatticeModelCreator<collisionModel_T>::makeLatticeModel(omegaFieldID, omega);
784

785
786
        uint_t numberCreatedParticlesA = uint_c(0);
        uint_t numberCreatedParticlesB = uint_c(0);
787

788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
        auto ccdID = blocks->addBlockData(pe::ccd::createHashGridsDataHandling(globalBodyStorage, bodyStorageID), "CCD");

        // set up collision response, here DEM solver
        auto fcdID = blocks->addBlockData(pe::fcd::createGenericFCDDataHandling<BodyTypeTuple, pe::fcd::AnalyticCollideFunctor>(), "FCD");

        std::unique_ptr<pe::cr::ICR> cr;
        if (useDEM) {
            cr = std::make_unique<pe::cr::DEM>(globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID);
            WALBERLA_LOG_INFO_ON_ROOT("Using DEM!");

        } else {
            cr = std::make_unique<pe::cr::HCSITS>(globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID);
            pe::cr::HCSITS *hcsits = static_cast<pe::cr::HCSITS *>(cr.get());
            hcsits->setMaxIterations(HCSITSMaxIterations);
            hcsits->setRelaxationParameter(HCSITSRelaxationParameter);
            WALBERLA_LOG_INFO_ON_ROOT("Using HCSITS!");
        }

        // set up synchronization procedure
        std::function<void(void)> syncCall = std::bind(pe::syncNextNeighbors<BodyTypeTuple>, std::ref(blocks->getBlockForest()), bodyStorageID,
                                                       static_cast<WcTimingTree *>(NULL), overlap, false);

810
        if (initializeFromCheckPointFile) {
811

812
            WALBERLA_LOG_INFO_ON_ROOT("Initializing data from checkpoint file!");
813

814
            WALBERLA_LOG_RESULT_ON_ROOT("Reading pdf from checkpoint file");
815

816
817
818
819
820
            shared_ptr<lbm::internal::PdfFieldHandling<LatticeModel_T> > pdfDataHandling =
                    make_shared<lbm::internal::PdfFieldHandling<LatticeModel_T> >(blocks, latticeModel, false, Vector3<real_t>(0), real_t(1),
                                                                                  FieldGhostLayers, field::fzyx);
            // add pdf field
            pdfFieldID = (blocks->getBlockStorage()).loadBlockData(checkpointPathPdf, pdfDataHandling, "pdf field (fzyx)");
821

822
823
824
825
            WALBERLA_LOG_RESULT_ON_ROOT("Reading body from checkpoint file");

            bodyStorageID = blocks->loadBlockData(checkpointPathBodies, pe::createStorageDataHandling<BodyTypeTuple>());

826
        } else {
827

828
            std::mt19937 rnd(MPIManager::instance()->rank());
829

830
831
832
            std::default_random_engine gen(MPIManager::instance()->rank());
            std::normal_distribution<real_t> distributionA(radiusA, radiusStandardDeviationA);
            std::normal_distribution<real_t> distributionB(radiusB, radiusStandardDeviationB);
833

834
            pe::SphereID sphere;
835

836
            real_t ratioAB = real_c(config.nrParticlesA) / real_c(config.nrParticlesB);
837

838
839
840
            real_t xInterval = initialSphereDistance;
            real_t yInterval = initialSphereDistance;
            real_t zInterval = initialSphereDistance;
841

842
843
            const uint_t rows = uint_c(real_c(config.xlength) / xInterval);
            const uint_t columns = uint_c(real_c(config.zlength) / zInterval);
844

845
846
            xInterval = real_c(config.xlength) / real_c(rows);
            zInterval = real_c(config.zlength) / real_c(columns);
847

848
849
850
            if (twoSpecies) {
                for (unsigned int i = 0; i < (config.nrParticlesA + config.nrParticlesB); ++i) {
                    if (real_c(numberCreatedParticlesA + 1) < ratioAB * real_c(numberCreatedParticlesB + 1)) {
851

852
853
854
855
856
                        sphere = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, i + 1,
                                                  Vector3<real_t>(0.5 * xInterval + xInterval * real_c((i % (rows * columns)) / columns),
                                                                  config.offsetBot + 0.5 * yInterval + yInterval * real_c(i / (rows * columns)),
                                                                  0.5 * zInterval + zInterval * real_c((i % columns))),
                                                  distributionA(gen), materialA);
857

858
                        if (sphere != NULL) {
859

860
                            ++numberCreatedParticlesA;
861

862
863
864
865
866
867
868
                            if (initialSphereVel > 0.0) {
                                sphere->setLinearVel(math::realRandom(-initialSphereVel, initialSphereVel, rnd),
                                                     math::realRandom(-initialSphereVel, initialSphereVel, rnd),
                                                     math::realRandom(-initialSphereVel, initialSphereVel, rnd));
                            }
                        }
                    } else {
869

870
871
872
873
874
                        sphere = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, i + 1,
                                                  Vector3<real_t>(0.5 * xInterval + xInterval * real_c((i % (rows * columns)) / columns),
                                                                  config.offsetBot + 0.5 * yInterval + yInterval * real_c(i / (rows * columns)),
                                                                  0.5 * zInterval + zInterval * real_c((i % columns))),
                                                  distributionB(gen), materialB);
875

876
                        if (sphere != NULL) {
877

878
                            ++numberCreatedParticlesB;
879

880
881
882
883
884
885
                            if (initialSphereVel > 0.0) {
                                sphere->setLinearVel(math::realRandom(-initialSphereVel, initialSphereVel, rnd),
                                                     math::realRandom(-initialSphereVel, initialSphereVel, rnd),
                                                     math::realRandom(-initialSphereVel, initialSphereVel, rnd));
                            }
                        }
886

887
888
889
                    }
                }
            } else {
890
                for (unsigned int i = 0; i < (config.nrParticlesA); ++i) {
891

892
893
894
895
896
                    sphere = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, i + 1,
                                              Vector3<real_t>(0.5 * xInterval + xInterval * real_c((i % (rows * columns)) / columns),
                                                              config.offsetBot + 0.5 * yInterval + yInterval * real_c(i / (rows * columns)),
                                                              0.5 * zInterval + zInterval * real_c((i % columns))),
                                              distributionA(gen), materialA);
897

898
                    if (sphere != NULL) {
899

900
                        ++numberCreatedParticlesA;
901

902
903
904
905
906
907
908
909
                        if (initialSphereVel > 0.0) {
                            sphere->setLinearVel(math::realRandom(-initialSphereVel, initialSphereVel, rnd),
                                                 math::realRandom(-initialSphereVel, initialSphereVel, rnd),
                                                 math::realRandom(-initialSphereVel, initialSphereVel, rnd));
                        }
                    }
                }
            }
910

911
912
913
914
            WALBERLA_MPI_SECTION() {
                mpi::allReduceInplace(numberCreatedParticlesA, mpi::SUM);
                mpi::allReduceInplace(numberCreatedParticlesB, mpi::SUM);
            }
915

916
917
918
919
            WALBERLA_ROOT_SECTION() {
                std::cout << "Number of created particles A is:     " << numberCreatedParticlesA << std::endl;
                std::cout << "Number of created particles B is:     " << numberCreatedParticlesB << std::endl << std::endl;
            }
920

921
922
923
924
925
926
927
928
929
930
931
932
            ///////////////////
            // PACKED BED ////
            //////////////////

            if (createPackedBed) {

                // spheres_init
                auto sphereVtkOutputInit = make_shared<pe::SphereVtkOutput>(bodyStorageID, blocks->getBlockStorage());
                auto sphereVtkWriterInit = vtk::createVTKOutput_PointData(sphereVtkOutputInit, "spheres_init");

                sphereVtkWriterInit->write();

933
                cr->setGlobalLinearAcceleration(pe::Vec3(0.0, -(config.gravity), 0.0));
934

935
936
937
938
939
940
941
942
                real_t avgVel = real_c(100.0);
                real_t maxVel = real_c(100.0);
                real_t currentSettleVelocity = real_c(100.0);
                uint_t itPE = uint_c(0);
                real_t velocityConversion = config.dx_SI / config.dt_SI;

                while ((currentSettleVelocity > settleVelocity) || (itPE < 4999)) {

943
944
                    cr->timestep(peAccelerator / real_c(config.peSubCycles));
                    syncCall();
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
                    ++itPE;

                    if (writePackedBed && itPE % frequencyPackedBed == 0) {
                        sphereVtkWriterInit->write();
                    }

                    if ((itPE) % settleCheckFrequency == 0) {

                        avgVel =
                                FBfunc::getAvgVel(blocks, bodyStorageID, numberCreatedParticlesA + numberCreatedParticlesB) * velocityConversion;
                        maxVel = FBfunc::getMaxVel(blocks, bodyStorageID) * velocityConversion;

                        if (checkMaximumParticleVelocity) {
                            currentSettleVelocity = maxVel;
                        } else {
                            currentSettleVelocity = avgVel;
                        }

                        WALBERLA_ROOT_SECTION() {
                            std::cout << "Time step " << itPE << std::endl;
                            std::cout << "Average Velocity:   " << avgVel << " m/s" << std::endl;
                            std::cout << "Maximum Velocity:   " << maxVel << " m/s" << std::endl << std::endl;
                        }
                    }
969
970
971
                }


972
                cr->setGlobalLinearAcceleration(pe::Vec3(0.0, 0.0, 0.0));
973
                FBfunc::restSphere(blocks, bodyStorageID);
974

975
976
977
                avgVel = FBfunc::getAvgVel(blocks, bodyStorageID, (numberCreatedParticlesA + numberCreatedParticlesB)) * velocityConversion;
                avgVel = FBfunc::getAvgVel(blocks, bodyStorageID, (numberCreatedParticlesA + numberCreatedParticlesB)) * velocityConversion;
                maxVel = FBfunc::getMaxVel(blocks, bodyStorageID) * velocityConversion;
978

979
980
981
982
                WALBERLA_ROOT_SECTION() {
                    std::cout << "||||||||||||| END OF PE SIM ||||||||||||||" << std::endl;
                    std::cout << "Average Velocity:   " << avgVel << " m/s" << std::endl;
                    std::cout << "Maximum Velocity:   " << maxVel << " m/s" << std::endl << std::endl;
983
                }
984
            }
985
986


987
988
989
990
991
992
993
994
995
996
997
998
999
1000
            ////////////////////
            // PACKED BED END //
            ////////////////////


            WALBERLA_ROOT_SECTION() {
                std::cout << "\nSetup Information\n" << std::endl;
                std::cout << "Physical time step is:                " << config.dt_SI << " s" << std::endl;
                std::cout << "Physical cell size is:                " << config.dx_SI << " m" << std::endl;
                std::cout << "Particle Reynoldsnumber(LBM) is:      " << particleReynoldsNumber_LBM << std::endl;
                std::cout << "Froudenumber(LBM) is:                 " << FroudeNumber_LBM << std::endl;
                std::cout << "Minimum Fluidization Re number is:    " << REpmf << std::endl;
                std::cout << "LBM velocity is:                      " << config.velocity << std::endl;
                std::cout << "LBM viscosity is:                     " << config.kinViscosity << std::endl;
For faster browsing, not all history is shown. View entire blame