MESA_PD_GranularGas.cpp 17.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
//======================================================================================================================
//
//  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   MESA_PD_GranularGas.cpp
//! \author Sebastian Eibl <sebastian.eibl@fau.de>
//
//======================================================================================================================

Sebastian Eibl's avatar
Sebastian Eibl committed
21
22
23
24
25
26
27
28
29
30
#include "Accessor.h"
#include "check.h"
#include "Contact.h"
#include "CreateParticles.h"
#include "NodeTimings.h"
#include "Parameters.h"
#include "SelectProperty.h"
#include "sortParticleStorage.h"
#include "SQLProperties.h"

31
32
33
34
35
36
37
38
#include <mesa_pd/vtk/ParticleVtkOutput.h>

#include <mesa_pd/collision_detection/AnalyticContactDetection.h>
#include <mesa_pd/data/LinkedCells.h>
#include <mesa_pd/data/ParticleAccessor.h>
#include <mesa_pd/data/ParticleStorage.h>
#include <mesa_pd/data/ShapeStorage.h>
#include <mesa_pd/domain/BlockForestDomain.h>
39
#include <mesa_pd/kernel/AssocToBlock.h>
40
41
42
43
44
45
46
47
#include <mesa_pd/kernel/DoubleCast.h>
#include <mesa_pd/kernel/ExplicitEulerWithShape.h>
#include <mesa_pd/kernel/InsertParticleIntoLinkedCells.h>
#include <mesa_pd/kernel/ParticleSelector.h>
#include <mesa_pd/kernel/SpringDashpot.h>
#include <mesa_pd/mpi/ContactFilter.h>
#include <mesa_pd/mpi/ReduceProperty.h>
#include <mesa_pd/mpi/SyncNextNeighbors.h>
48
#include <mesa_pd/mpi/SyncNextNeighborsBlockForest.h>
49
50
51
52
53
54
55
56
57
58
59
60
61
62

#include <mesa_pd/mpi/notifications/ForceTorqueNotification.h>

#include <blockforest/BlockForest.h>
#include <blockforest/Initialization.h>
#include <core/Abort.h>
#include <core/Environment.h>
#include <core/math/Random.h>
#include <core/mpi/Reduce.h>
#include <core/grid_generator/SCIterator.h>
#include <core/logging/Logging.h>
#include <core/OpenMP.h>
#include <core/timing/Timer.h>
#include <core/waLBerlaBuildInfo.h>
63
#include <sqlite/SQLite.h>
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
#include <vtk/VTKOutput.h>

#include <functional>
#include <memory>
#include <string>
#include <type_traits>

namespace walberla {
namespace mesa_pd {

int main( int argc, char ** argv )
{
   using namespace walberla::timing;

   Environment env(argc, argv);
   auto mpiManager = walberla::mpi::MPIManager::instance();
   mpiManager->useWorldComm();

   logging::Logging::instance()->setStreamLogLevel(logging::Logging::INFO);
   logging::Logging::instance()->setFileLogLevel(logging::Logging::INFO);

   WALBERLA_LOG_INFO_ON_ROOT( "config file: " << argv[1] );
   WALBERLA_LOG_INFO_ON_ROOT( "waLBerla Revision: " << WALBERLA_GIT_SHA1 );

   math::seedRandomGenerator( static_cast<unsigned int>(1337 * mpiManager->worldRank()) );

   WALBERLA_LOG_INFO_ON_ROOT("*** READING CONFIG FILE ***");
   auto cfg = env.config();
   if (cfg == nullptr) WALBERLA_ABORT("No config specified!");
   const Config::BlockHandle mainConf  = cfg->getBlock( "GranularGas" );

Sebastian Eibl's avatar
Sebastian Eibl committed
95
96
   Parameters params;
   loadFromConfig(params, mainConf);
97
98
99
100
101
102
103
104
105

   WALBERLA_LOG_INFO_ON_ROOT("*** BLOCKFOREST ***");
   // create forest
   auto forest = blockforest::createBlockForestFromConfig( mainConf );
   if (!forest)
   {
      WALBERLA_LOG_INFO_ON_ROOT( "No BlockForest created ... exiting!");
      return EXIT_SUCCESS;
   }
106
   auto domain = std::make_shared<domain::BlockForestDomain> (forest);
107
108
109
110
111
112
113
114
115
116
117
118
119
120

   auto simulationDomain = forest->getDomain();
   auto localDomain = forest->begin()->getAABB();
   for (auto& blk : *forest)
   {
      localDomain.merge(blk.getAABB());
   }

   WALBERLA_LOG_INFO_ON_ROOT("*** SETUP - START ***");

   //init data structures
   auto ps = std::make_shared<data::ParticleStorage>(100);
   auto ss = std::make_shared<data::ShapeStorage>();
   ParticleAccessorWithShape accessor(ps, ss);
Sebastian Eibl's avatar
Sebastian Eibl committed
121
   data::LinkedCells     lc(localDomain.getExtended(params.spacing), params.spacing );
122

Sebastian Eibl's avatar
Sebastian Eibl committed
123
   auto  smallSphere = ss->create<data::Sphere>( params.radius );
124
125
126
   ss->shapes[smallSphere]->updateMassAndInertia(real_t(2707));
   for (auto& iBlk : *forest)
   {
Sebastian Eibl's avatar
Sebastian Eibl committed
127
      for (auto pt : grid_generator::SCGrid(iBlk.getAABB(),
128
                                            Vector3<real_t>(params.spacing) * real_c(0.5) + params.shift,
Sebastian Eibl's avatar
Sebastian Eibl committed
129
                                            params.spacing))
130
131
      {
         WALBERLA_CHECK(iBlk.getAABB().contains(pt));
Sebastian Eibl's avatar
Sebastian Eibl committed
132
         createSphere(*ps, pt, params.radius, smallSphere);
133
134
135
136
137
138
      }
   }
   int64_t numParticles = int64_c(ps->size());
   walberla::mpi::reduceInplace(numParticles, walberla::mpi::SUM);
   WALBERLA_LOG_INFO_ON_ROOT("#particles created: " << numParticles);

Sebastian Eibl's avatar
Sebastian Eibl committed
139
   auto shift = (params.spacing - params.radius - params.radius) * real_t(0.5);
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
   auto confiningDomain = simulationDomain.getExtended(shift);

   if (!forest->isPeriodic(0))
   {
      createPlane(*ps, *ss, confiningDomain.minCorner(), Vec3(+1,0,0));
      createPlane(*ps, *ss, confiningDomain.maxCorner(), Vec3(-1,0,0));
   }

   if (!forest->isPeriodic(1))
   {
      createPlane(*ps, *ss, confiningDomain.minCorner(), Vec3(0,+1,0));
      createPlane(*ps, *ss, confiningDomain.maxCorner(), Vec3(0,-1,0));
   }

   if (!forest->isPeriodic(2))
   {
      createPlane(*ps, *ss, confiningDomain.minCorner(), Vec3(0,0,+1));
      createPlane(*ps, *ss, confiningDomain.maxCorner(), Vec3(0,0,-1));
   }

   WALBERLA_LOG_INFO_ON_ROOT("*** SETUP - END ***");

   WALBERLA_LOG_INFO_ON_ROOT("*** VTK ***");
Sebastian Eibl's avatar
Sebastian Eibl committed
163
   auto vtkDomainOutput = walberla::vtk::createVTKOutput_DomainDecomposition( forest, "domain_decomposition", 1, params.vtk_out, "simulation_step" );
164
   auto vtkOutput       = make_shared<mesa_pd::vtk::ParticleVtkOutput>(ps) ;
Sebastian Eibl's avatar
Sebastian Eibl committed
165
   auto vtkWriter       = walberla::vtk::createVTKOutput_PointData(vtkOutput, "Bodies", 1, params.vtk_out, "simulation_step", false, false);
166
167
168
169
170
171
   vtkOutput->addOutput<SelectRank>("rank");
   vtkOutput->addOutput<data::SelectParticleOwner>("owner");
   //   vtkDomainOutput->write();

   WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - START ***");
   // Init kernels
Sebastian Eibl's avatar
Sebastian Eibl committed
172
   kernel::ExplicitEulerWithShape        explicitEulerWithShape( params.dt );
173
174
175
176
177
178
179
   kernel::InsertParticleIntoLinkedCells ipilc;
   kernel::SpringDashpot                 dem(1);
   dem.setStiffness(0, 0, real_t(0));
   dem.setDampingN (0, 0, real_t(0));
   dem.setDampingT (0, 0, real_t(0));
   dem.setFriction (0, 0, real_t(0));
   collision_detection::AnalyticContactDetection              acd;
180
   kernel::AssocToBlock                  assoc(forest);
181
182
183
   kernel::DoubleCast                    double_cast;
   mpi::ContactFilter                    contact_filter;
   mpi::ReduceProperty                   RP;
184
   mpi::SyncNextNeighborsBlockForest     SNN;
185
186

   // initial sync
187
188
   ps->forEachParticle(false, kernel::SelectLocal(), accessor, assoc, accessor);
   SNN(*ps, forest, domain);
Sebastian Eibl's avatar
Sebastian Eibl committed
189
190
   sortParticleStorage(*ps, params.sorting, lc.domain_, uint_c(lc.numCellsPerDim_[0]));
//   vtkWriter->write();
191

Sebastian Eibl's avatar
Sebastian Eibl committed
192
   for (int64_t outerIteration = 0; outerIteration < params.numOuterIterations; ++outerIteration)
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
   {
      WALBERLA_LOG_INFO_ON_ROOT("*** RUNNING OUTER ITERATION " << outerIteration << " ***");

      WcTimer      timer;
      WcTimingPool tp;
      auto    SNNBytesSent     = SNN.getBytesSent();
      auto    SNNBytesReceived = SNN.getBytesReceived();
      auto    SNNSends         = SNN.getNumberOfSends();
      auto    SNNReceives      = SNN.getNumberOfReceives();
      auto    RPBytesSent      = RP.getBytesSent();
      auto    RPBytesReceived  = RP.getBytesReceived();
      auto    RPSends          = RP.getNumberOfSends();
      auto    RPReceives       = RP.getNumberOfReceives();
      int64_t contactsChecked  = 0;
      int64_t contactsDetected = 0;
      int64_t contactsTreated  = 0;
Sebastian Eibl's avatar
Sebastian Eibl committed
209
      if (params.bBarrier) WALBERLA_MPI_BARRIER();
Sebastian Eibl's avatar
Sebastian Eibl committed
210
      timer.start();
Sebastian Eibl's avatar
Sebastian Eibl committed
211
      for (int64_t i=0; i < params.simulationSteps; ++i)
212
213
214
215
216
217
      {
         //      if (i % visSpacing == 0)
         //      {
         //         vtkWriter->write();
         //      }

218
219
220
221
222
         tp["AssocToBlock"].start();
         ps->forEachParticle(false, kernel::SelectLocal(), accessor, assoc, accessor);
         if (params.bBarrier) WALBERLA_MPI_BARRIER();
         tp["AssocToBlock"].end();

223
224
225
         tp["GenerateLinkedCells"].start();
         lc.clear();
         ps->forEachParticle(true, kernel::SelectAll(), accessor, ipilc, accessor, lc);
Sebastian Eibl's avatar
Sebastian Eibl committed
226
         if (params.bBarrier) WALBERLA_MPI_BARRIER();
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
         tp["GenerateLinkedCells"].end();

         tp["DEM"].start();
         contactsChecked  = 0;
         contactsDetected = 0;
         contactsTreated  = 0;
         lc.forEachParticlePairHalf(true,
                                    kernel::SelectAll(),
                                    accessor,
                                    [&](const size_t idx1, const size_t idx2, auto& ac)
         {
            ++contactsChecked;
            if (double_cast(idx1, idx2, ac, acd, ac ))
            {
               ++contactsDetected;
242
               if (contact_filter(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), *domain))
243
244
245
246
247
248
249
               {
                  ++contactsTreated;
                  dem(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), acd.getContactNormal(), acd.getPenetrationDepth());
               }
            }
         },
         accessor );
Sebastian Eibl's avatar
Sebastian Eibl committed
250
         if (params.bBarrier) WALBERLA_MPI_BARRIER();
251
252
253
254
         tp["DEM"].end();

         tp["ReduceForce"].start();
         RP.operator()<ForceTorqueNotification>(*ps);
Sebastian Eibl's avatar
Sebastian Eibl committed
255
         if (params.bBarrier) WALBERLA_MPI_BARRIER();
256
257
258
259
260
         tp["ReduceForce"].end();

         tp["Euler"].start();
         //ps->forEachParticle(false, [&](const size_t idx){WALBERLA_CHECK_EQUAL(ps->getForce(idx), Vec3(0,0,0), *(*ps)[idx] << "\n" << idx);});
         ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEulerWithShape, accessor);
Sebastian Eibl's avatar
Sebastian Eibl committed
261
         if (params.bBarrier) WALBERLA_MPI_BARRIER();
262
263
264
         tp["Euler"].end();

         tp["SNN"].start();
265
         SNN(*ps, forest, domain);
Sebastian Eibl's avatar
Sebastian Eibl committed
266
         if (params.bBarrier) WALBERLA_MPI_BARRIER();
267
268
269
         tp["SNN"].end();
      }
      timer.end();
Sebastian Eibl's avatar
Sebastian Eibl committed
270

271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
      SNNBytesSent     = SNN.getBytesSent();
      SNNBytesReceived = SNN.getBytesReceived();
      SNNSends         = SNN.getNumberOfSends();
      SNNReceives      = SNN.getNumberOfReceives();
      RPBytesSent      = RP.getBytesSent();
      RPBytesReceived  = RP.getBytesReceived();
      RPSends          = RP.getNumberOfSends();
      RPReceives       = RP.getNumberOfReceives();
      walberla::mpi::reduceInplace(SNNBytesSent, walberla::mpi::SUM);
      walberla::mpi::reduceInplace(SNNBytesReceived, walberla::mpi::SUM);
      walberla::mpi::reduceInplace(SNNSends, walberla::mpi::SUM);
      walberla::mpi::reduceInplace(SNNReceives, walberla::mpi::SUM);
      walberla::mpi::reduceInplace(RPBytesSent, walberla::mpi::SUM);
      walberla::mpi::reduceInplace(RPBytesReceived, walberla::mpi::SUM);
      walberla::mpi::reduceInplace(RPSends, walberla::mpi::SUM);
      walberla::mpi::reduceInplace(RPReceives, walberla::mpi::SUM);
      auto cC = walberla::mpi::reduce(contactsChecked, walberla::mpi::SUM);
      auto cD = walberla::mpi::reduce(contactsDetected, walberla::mpi::SUM);
      auto cT = walberla::mpi::reduce(contactsTreated, walberla::mpi::SUM);
      WALBERLA_LOG_DEVEL_ON_ROOT( "SNN bytes communicated:   " << SNNBytesSent << " / " << SNNBytesReceived );
      WALBERLA_LOG_DEVEL_ON_ROOT( "SNN communication partners: " << SNNSends << " / " << SNNReceives );
      WALBERLA_LOG_DEVEL_ON_ROOT( "RP bytes communicated:  " << RPBytesSent << " / " << RPBytesReceived );
      WALBERLA_LOG_DEVEL_ON_ROOT( "RP communication partners: " << RPSends << " / " << RPReceives );
      WALBERLA_LOG_DEVEL_ON_ROOT( "contacts checked/detected/treated: " << cC << " / " << cD << " / " << cT );

Sebastian Eibl's avatar
Sebastian Eibl committed
296
297
298
299
300
301
      auto timer_reduced = walberla::timing::getReduced(timer, REDUCE_TOTAL, 0);
      double PUpS = 0.0;
      WALBERLA_ROOT_SECTION()
      {
         WALBERLA_LOG_INFO_ON_ROOT(*timer_reduced);
         WALBERLA_LOG_INFO_ON_ROOT("runtime: " << timer_reduced->max());
Sebastian Eibl's avatar
Sebastian Eibl committed
302
         PUpS = double_c(numParticles) * double_c(params.simulationSteps) / double_c(timer_reduced->max());
Sebastian Eibl's avatar
Sebastian Eibl committed
303
304
305
         WALBERLA_LOG_INFO_ON_ROOT("PUpS: " << PUpS);
      }

306
307
308
309
      auto tp_reduced = tp.getReduced();
      WALBERLA_LOG_INFO_ON_ROOT(*tp_reduced);
      WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - END ***");

Sebastian Eibl's avatar
Sebastian Eibl committed
310
      if (params.checkSimulation)
311
      {
Sebastian Eibl's avatar
Sebastian Eibl committed
312
         check(*ps, *forest, params.spacing);
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
      }

      WALBERLA_LOG_INFO_ON_ROOT("*** SQL OUTPUT - START ***");
      numParticles = 0;
      int64_t numGhostParticles = 0;
      ps->forEachParticle(false,
                          kernel::SelectAll(),
                          accessor,
                          [&numParticles, &numGhostParticles](const size_t idx, auto& ac)
      {
         if (data::particle_flags::isSet( ac.getFlagsRef(idx), data::particle_flags::GHOST))
         {
            ++numGhostParticles;
         } else
         {
            ++numParticles;
         }
      },
      accessor);
      walberla::mpi::reduceInplace(numParticles, walberla::mpi::SUM);
      walberla::mpi::reduceInplace(numGhostParticles, walberla::mpi::SUM);
      walberla::mpi::reduceInplace(contactsChecked, walberla::mpi::SUM);
      walberla::mpi::reduceInplace(contactsDetected, walberla::mpi::SUM);
      walberla::mpi::reduceInplace(contactsTreated, walberla::mpi::SUM);
      double linkedCellsVolume = lc.domain_.volume();
      walberla::mpi::reduceInplace(linkedCellsVolume, walberla::mpi::SUM);
      size_t numLinkedCells = lc.cells_.size();
      walberla::mpi::reduceInplace(numLinkedCells, walberla::mpi::SUM);
341
342
343
      size_t local_aabbs         = domain->getNumLocalAABBs();
      size_t neighbor_subdomains = domain->getNumNeighborSubdomains();
      size_t neighbor_processes  = domain->getNumNeighborProcesses();
344
345
346
      walberla::mpi::reduceInplace(local_aabbs, walberla::mpi::SUM);
      walberla::mpi::reduceInplace(neighbor_subdomains, walberla::mpi::SUM);
      walberla::mpi::reduceInplace(neighbor_processes, walberla::mpi::SUM);
Sebastian Eibl's avatar
Sebastian Eibl committed
347
348

      uint_t runId = uint_c(-1);
349
350
351
352
353
354
355
356
357
358
      WALBERLA_ROOT_SECTION()
      {
         std::map< std::string, walberla::int64_t > integerProperties;
         std::map< std::string, double >            realProperties;
         std::map< std::string, std::string >       stringProperties;

         stringProperties["walberla_git"]         = WALBERLA_GIT_SHA1;
         stringProperties["tag"]                  = "mesa_pd";
         integerProperties["mpi_num_processes"]   = mpiManager->numProcesses();
         integerProperties["omp_max_threads"]     = omp_get_max_threads();
359
360
361
362
363
         realProperties["PUpS"]                   = double_c(PUpS);
         realProperties["timer_min"]              = timer_reduced->min();
         realProperties["timer_max"]              = timer_reduced->max();
         realProperties["timer_average"]          = timer_reduced->average();
         realProperties["timer_total"]            = timer_reduced->total();
Sebastian Eibl's avatar
Sebastian Eibl committed
364
         integerProperties["outerIteration"]      = int64_c(outerIteration);
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
         integerProperties["num_particles"]       = numParticles;
         integerProperties["num_ghost_particles"] = numGhostParticles;
         integerProperties["contacts_checked"]    = contactsChecked;
         integerProperties["contacts_detected"]   = contactsDetected;
         integerProperties["contacts_treated"]    = contactsTreated;
         integerProperties["local_aabbs"]         = int64_c(local_aabbs);
         integerProperties["neighbor_subdomains"] = int64_c(neighbor_subdomains);
         integerProperties["neighbor_processes"]  = int64_c(neighbor_processes);
         integerProperties["SNNBytesSent"]        = SNNBytesSent;
         integerProperties["SNNBytesReceived"]    = SNNBytesReceived;
         integerProperties["SNNSends"]            = SNNSends;
         integerProperties["SNNReceives"]         = SNNReceives;
         integerProperties["RPBytesSent"]         = RPBytesSent;
         integerProperties["RPBytesReceived"]     = RPBytesReceived;
         integerProperties["RPSends"]             = RPSends;
         integerProperties["RPReceives"]          = RPReceives;
         realProperties["linkedCellsVolume"]      = linkedCellsVolume;
         integerProperties["numLinkedCells"]      = int64_c(numLinkedCells);
Sebastian Eibl's avatar
Sebastian Eibl committed
383
384
385
386
387
         realProperties["PUpS"]                   = double_c(PUpS);
         realProperties["timer_min"]              = timer_reduced->min();
         realProperties["timer_max"]              = timer_reduced->max();
         realProperties["timer_average"]          = timer_reduced->average();
         realProperties["timer_total"]            = timer_reduced->total();
Sebastian Eibl's avatar
Sebastian Eibl committed
388

389
         addBuildInfoToSQL( integerProperties, realProperties, stringProperties );
Sebastian Eibl's avatar
Sebastian Eibl committed
390
391
392
393
         saveToSQL(params, integerProperties, realProperties, stringProperties );
         addDomainPropertiesToSQL(*forest, integerProperties, realProperties, stringProperties);
         addSlurmPropertiesToSQL(integerProperties, realProperties, stringProperties);

394
395
         runId = sqlite::storeRunInSqliteDB( params.sqlFile, integerProperties, stringProperties, realProperties );
         sqlite::storeTimingPoolInSqliteDB( params.sqlFile, runId, *tp_reduced, "Timeloop" );
Sebastian Eibl's avatar
Sebastian Eibl committed
396
397
398
399
400
      }

      if (params.storeNodeTimings)
      {
         storeNodeTimings(runId, params.sqlFile, "NodeTiming", tp);
401
402
403
404
405
406
407
408
409
410
411
412
413
414
      }
      WALBERLA_LOG_INFO_ON_ROOT("*** SQL OUTPUT - END ***");
   }

   return EXIT_SUCCESS;
}

} // namespace mesa_pd
} // namespace walberla

int main( int argc, char* argv[] )
{
   return walberla::mesa_pd::main( argc, argv );
}