diff --git a/apps/benchmarks/MaterialTransport/thermal_PSM/thermal_PSM_particles/thermalPSMParticles.cpp b/apps/benchmarks/MaterialTransport/thermal_PSM/thermal_PSM_particles/thermalPSMParticles.cpp
index 9cb01d6330a41d8a1df69c5503b83c496d6de6e7..cdd13348b24b5dfdee8f01cf36ab235f0c8a9d71 100644
--- a/apps/benchmarks/MaterialTransport/thermal_PSM/thermal_PSM_particles/thermalPSMParticles.cpp
+++ b/apps/benchmarks/MaterialTransport/thermal_PSM/thermal_PSM_particles/thermalPSMParticles.cpp
@@ -46,15 +46,36 @@
#include "lbm_mesapd_coupling/DataTypesCodegen.h"
#include "lbm_mesapd_coupling/partially_saturated_cells_method/codegen/PSMSweepCollection.h"
+#include "lbm_mesapd_coupling/utility/AddForceOnParticlesKernel.h"
+#include "lbm_mesapd_coupling/utility/AddHydrodynamicInteractionKernel.h"
+#include "lbm_mesapd_coupling/utility/AverageHydrodynamicForceTorqueKernel.h"
+#include "lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h"
+#include "lbm_mesapd_coupling/utility/LubricationCorrectionKernel.h"
#include "lbm_mesapd_coupling/utility/ParticleSelector.h"
+#include "lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
+#include "mesa_pd/collision_detection/AnalyticContactDetection.h"
#include "mesa_pd/data/DataTypes.h"
+#include "mesa_pd/data/LinkedCells.h"
#include "mesa_pd/data/ParticleAccessorWithShape.h"
#include "mesa_pd/data/ParticleStorage.h"
#include "mesa_pd/data/ShapeStorage.h"
+#include "mesa_pd/data/shape/HalfSpace.h"
#include "mesa_pd/data/shape/Sphere.h"
+#include "mesa_pd/domain/BlockForestDataHandling.h"
#include "mesa_pd/domain/BlockForestDomain.h"
+#include "mesa_pd/kernel/AssocToBlock.h"
+#include "mesa_pd/kernel/DoubleCast.h"
+#include "mesa_pd/kernel/InsertParticleIntoLinkedCells.h"
+#include "mesa_pd/kernel/LinearSpringDashpot.h"
+#include "mesa_pd/kernel/ParticleSelector.h"
+#include "mesa_pd/kernel/VelocityVerlet.h"
+#include "mesa_pd/mpi/ContactFilter.h"
+#include "mesa_pd/mpi/ReduceContactHistory.h"
+#include "mesa_pd/mpi/ReduceProperty.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
+#include "mesa_pd/mpi/notifications/ForceTorqueNotification.h"
+#include "mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h"
#include "mesa_pd/vtk/ParticleVtkOutput.h"
#include "sqlite/SQLite.h"
@@ -102,6 +123,100 @@ const FlagUID Inflow_Concentration_Flag("Inflow_Concentration");
// const FlagUID Dirichlet_Concentration_Flag("Dirichlet_Concentration");
const FlagUID Neumann_Concentration_Flag("Neumann_Concentration");
+void createPlane(const shared_ptr< mesa_pd::data::ParticleStorage >& ps,
+ const shared_ptr< mesa_pd::data::ShapeStorage >& ss, Vector3< real_t > position,
+ Vector3< real_t > normal)
+{
+ mesa_pd::data::Particle&& p0 = *ps->create(true);
+ p0.setPosition(position);
+ p0.setInteractionRadius(std::numeric_limits< real_t >::infinity());
+ p0.setShapeID(ss->create< mesa_pd::data::HalfSpace >(normal));
+ p0.setOwner(mpi::MPIManager::instance()->rank());
+ p0.setType(0);
+ mesa_pd::data::particle_flags::set(p0.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+ mesa_pd::data::particle_flags::set(p0.getFlagsRef(), mesa_pd::data::particle_flags::FIXED);
+}
+
+void createPlaneSetup(const shared_ptr< mesa_pd::data::ParticleStorage >& ps,
+ const shared_ptr< mesa_pd::data::ShapeStorage >& ss, const math::AABB& simulationDomain,
+ bool periodicInX, bool periodicInY, real_t offsetAtInflow, real_t offsetAtOutflow)
+{
+ createPlane(ps, ss, simulationDomain.minCorner() + Vector3< real_t >(0, 0, offsetAtInflow),
+ Vector3< real_t >(0, 0, 1));
+ createPlane(ps, ss, simulationDomain.maxCorner() + Vector3< real_t >(0, 0, offsetAtOutflow),
+ Vector3< real_t >(0, 0, -1));
+
+ if (!periodicInX)
+ {
+ createPlane(ps, ss, simulationDomain.minCorner(), Vector3< real_t >(1, 0, 0));
+ createPlane(ps, ss, simulationDomain.maxCorner(), Vector3< real_t >(-1, 0, 0));
+ }
+
+ if (!periodicInY)
+ {
+ createPlane(ps, ss, simulationDomain.minCorner(), Vector3< real_t >(0, 1, 0));
+ createPlane(ps, ss, simulationDomain.maxCorner(), Vector3< real_t >(0, -1, 0));
+ }
+}
+
+struct ParticleInfo
+{
+ real_t averageVelocity = 0_r;
+ real_t maximumVelocity = 0_r;
+ uint_t numParticles = 0;
+ real_t maximumHeight = 0_r;
+ real_t particleVolume = 0_r;
+ real_t heightOfMass = 0_r;
+
+ void allReduce()
+ {
+ walberla::mpi::allReduceInplace(numParticles, walberla::mpi::SUM);
+ walberla::mpi::allReduceInplace(averageVelocity, walberla::mpi::SUM);
+ walberla::mpi::allReduceInplace(maximumVelocity, walberla::mpi::MAX);
+ walberla::mpi::allReduceInplace(maximumHeight, walberla::mpi::MAX);
+ walberla::mpi::allReduceInplace(particleVolume, walberla::mpi::SUM);
+ walberla::mpi::allReduceInplace(heightOfMass, walberla::mpi::SUM);
+
+ averageVelocity /= real_c(numParticles);
+ heightOfMass /= particleVolume;
+ }
+};
+
+std::ostream& operator<<(std::ostream& os, ParticleInfo const& m)
+{
+ return os << "Particle Info: uAvg = " << m.averageVelocity << ", uMax = " << m.maximumVelocity
+ << ", numParticles = " << m.numParticles << ", zMax = " << m.maximumHeight << ", Vp = " << m.particleVolume
+ << ", zMass = " << m.heightOfMass;
+}
+
+template< typename Accessor_T >
+ParticleInfo evaluateParticleInfo(const Accessor_T& ac)
+{
+ static_assert(std::is_base_of_v< mesa_pd::data::IAccessor, Accessor_T >, "Provide a valid accessor");
+
+ ParticleInfo info;
+ for (uint_t i = 0; i < ac.size(); ++i)
+ {
+ if (isSet(ac.getFlags(i), mesa_pd::data::particle_flags::GHOST)) continue;
+ if (isSet(ac.getFlags(i), mesa_pd::data::particle_flags::GLOBAL)) continue;
+
+ ++info.numParticles;
+ real_t velMagnitude = ac.getLinearVelocity(i).length();
+ real_t particleVolume = ac.getShape(i)->getVolume();
+ real_t height = ac.getPosition(i)[2];
+ info.averageVelocity += velMagnitude;
+ info.maximumVelocity = std::max(info.maximumVelocity, velMagnitude);
+ info.maximumHeight = std::max(info.maximumHeight, height);
+ info.particleVolume += particleVolume;
+ info.heightOfMass += particleVolume * height;
+ }
+
+ info.allReduce();
+
+ return info;
+}
+
+
//////////
// MAIN //
//////////
@@ -121,57 +236,134 @@ int main(int argc, char** argv)
WALBERLA_LOG_INFO_ON_ROOT("build machine: " << std::string(WALBERLA_BUILD_MACHINE));
WALBERLA_LOG_INFO_ON_ROOT(*cfgFile);
- // Read config file
- Config::BlockHandle numericalSetup = cfgFile->getBlock("NumericalSetup");
+ // read all parameters from the config file
+
+ Config::BlockHandle physicalSetup = cfgFile->getBlock("PhysicalSetup");
+ const bool showcase = physicalSetup.getParameter< bool >("showcase");
+ const real_t xSize_SI = physicalSetup.getParameter< real_t >("xSize");
+ const real_t ySize_SI = physicalSetup.getParameter< real_t >("ySize");
+ const real_t zSize_SI = physicalSetup.getParameter< real_t >("zSize");
+ const bool periodicInX = physicalSetup.getParameter< bool >("periodicInX");
+ const bool periodicInY = physicalSetup.getParameter< bool >("periodicInY");
+ const real_t runtime_SI = physicalSetup.getParameter< real_t >("runtime");
+ const real_t uInflow_SI = physicalSetup.getParameter< real_t >("uInflow");
+ const real_t gravitationalAcceleration_SI = physicalSetup.getParameter< real_t >("gravitationalAcceleration");
+ const real_t kinematicViscosityFluid_SI = physicalSetup.getParameter< real_t >("kinematicViscosityFluid");
+ const real_t densityFluid_SI = physicalSetup.getParameter< real_t >("densityFluid");
+ const real_t particleDiameter_SI = physicalSetup.getParameter< real_t >("particleDiameter");
+ const real_t densityParticle_SI = physicalSetup.getParameter< real_t >("densityParticle");
+ const real_t dynamicFrictionCoefficient = physicalSetup.getParameter< real_t >("dynamicFrictionCoefficient");
+ const real_t coefficientOfRestitution = physicalSetup.getParameter< real_t >("coefficientOfRestitution");
+ const real_t collisionTimeFactor = physicalSetup.getParameter< real_t >("collisionTimeFactor");
+ const real_t particleGenerationSpacing_SI = physicalSetup.getParameter< real_t >("particleGenerationSpacing");
+
- const uint_t numXBlocks = numericalSetup.getParameter< uint_t >("numXBlocks");
- const uint_t numYBlocks = numericalSetup.getParameter< uint_t >("numYBlocks");
- const uint_t numZBlocks = numericalSetup.getParameter< uint_t >("numZBlocks");
- const uint_t cellsPerBlockPerDirection =
- numericalSetup.getParameter< uint_t >("cellsPerBlockPerDirection"); // same in all directions
+ Config::BlockHandle numericalSetup = cfgFile->getBlock("NumericalSetup");
+ const real_t dx_SI = numericalSetup.getParameter< real_t >("dx");
+ const real_t uInflow = numericalSetup.getParameter< real_t >("uInflow");
+ const real_t timeStepSize = numericalSetup.getParameter< real_t >("dt");
+ const uint_t numXBlocks = numericalSetup.getParameter< uint_t >("numXBlocks");
+ const uint_t numYBlocks = numericalSetup.getParameter< uint_t >("numYBlocks");
+ const uint_t numZBlocks = numericalSetup.getParameter< uint_t >("numZBlocks");
WALBERLA_CHECK_EQUAL(numXBlocks * numYBlocks * numZBlocks, uint_t(MPIManager::instance()->numProcesses()),
"When using GPUs, the number of blocks ("
<< numXBlocks * numYBlocks * numZBlocks << ") has to match the number of MPI processes ("
<< uint_t(MPIManager::instance()->numProcesses()) << ")");
- const bool periodicInY = numericalSetup.getParameter< bool >("periodicInY");
- const bool periodicInZ = numericalSetup.getParameter< bool >("periodicInZ");
- const bool sendDirectlyFromGPU = numericalSetup.getParameter< bool >("sendDirectlyFromGPU");
- const bool useCommunicationHiding = numericalSetup.getParameter< bool >("useCommunicationHiding");
- const Vector3< uint_t > frameWidth = numericalSetup.getParameter< Vector3< uint_t > >("frameWidth");
- const real_t gravity = numericalSetup.getParameter< real_t >("gravity");
-
- const real_t Thot = numericalSetup.getParameter< real_t >("Thot");
- const real_t Tcold = numericalSetup.getParameter< real_t >("Tcold");
- const real_t Pr = numericalSetup.getParameter< real_t >("PrandtlNumber");
- const real_t Ra = numericalSetup.getParameter< real_t >("RayleighNumber");
- const real_t Ma = numericalSetup.getParameter< real_t >("Mach");
- const real_t relaxationRate = numericalSetup.getParameter< real_t >("relaxationRate");
- const uint_t timeSteps = numericalSetup.getParameter< uint_t >("timeSteps");
-
- const bool useParticles = numericalSetup.getParameter< bool >("useParticles");
- const real_t particleDiameter = numericalSetup.getParameter< real_t >("particleDiameter");
- const real_t particleGenerationSpacing = numericalSetup.getParameter< real_t >("particleGenerationSpacing");
- const Vector3< real_t > generationDomainFraction =
- numericalSetup.getParameter< Vector3< real_t > >("generationDomainFraction");
- const Vector3< uint_t > generationPointOfReferenceOffset =
- numericalSetup.getParameter< Vector3< uint_t > >("generationPointOfReferenceOffset");
- const bool useParticleOffset = numericalSetup.getParameter< bool >("useParticleOffset");
+ if ((periodicInX && numXBlocks == 1) || (periodicInY && numYBlocks == 1))
+ {
+ WALBERLA_ABORT("The number of blocks must be greater than 1 in periodic dimensions.")
+ }
+
+ const bool useLubricationForces = numericalSetup.getParameter< bool >("useLubricationForces");
+ const uint_t numberOfParticleSubCycles = numericalSetup.getParameter< uint_t >("numberOfParticleSubCycles");
+ const bool useParticles = numericalSetup.getParameter< bool >("useParticles");
+ const bool useIntegrators = numericalSetup.getParameter< bool >("useIntegrators");
const Vector3< uint_t > particleSubBlockSize =
numericalSetup.getParameter< Vector3< uint_t > >("particleSubBlockSize");
+ const real_t linkedCellWidthRation = numericalSetup.getParameter< real_t >("linkedCellWidthRation");
+ const bool particleBarriers = numericalSetup.getParameter< bool >("particleBarriers");
+ const Vector3< real_t > generationDomainFraction =
+ numericalSetup.getParameter< Vector3< real_t > >("generationDomainFraction");
- if ((periodicInY && numYBlocks == 1) || (periodicInZ && numZBlocks == 1))
- {
- WALBERLA_LOG_WARNING_ON_ROOT("Using only 1 block in periodic dimensions can lead to unexpected behavior.")
- }
+
+ Config::BlockHandle TemperatureSetup = cfgFile->getBlock("TemperatureSetup");
+ const real_t Thot = TemperatureSetup.getParameter< real_t >("Thot");
+ const real_t Tcold = TemperatureSetup.getParameter< real_t >("Tcold");
+ const real_t Pr = TemperatureSetup.getParameter< real_t >("PrandtlNumber");
+ const real_t Ra = TemperatureSetup.getParameter< real_t >("RayleighNumber");
+ const real_t Ma = TemperatureSetup.getParameter< real_t >("Mach");
+ const real_t particleTemperature = TemperatureSetup.getParameter< real_t >("particleTemperature");
Config::BlockHandle outputSetup = cfgFile->getBlock("Output");
- const uint_t vtkSpacing = outputSetup.getParameter< uint_t >("vtkSpacing");
+ const real_t infoSpacing_SI = outputSetup.getParameter< real_t >("infoSpacing");
+ const real_t vtkSpacingParticles_SI = outputSetup.getParameter< real_t >("vtkSpacingParticles");
+ const real_t vtkSpacingFluid_SI = outputSetup.getParameter< real_t >("vtkSpacingFluid");
const std::string vtkFolder = outputSetup.getParameter< std::string >("vtkFolder");
const uint_t performanceLogFrequency = outputSetup.getParameter< uint_t >("performanceLogFrequency");
+
+
+
+ // convert SI units to simulation (LBM) units and check setup
+
+ Vector3< uint_t > domainSize(uint_c(std::ceil(xSize_SI / dx_SI)), uint_c(std::ceil(ySize_SI / dx_SI)),
+ uint_c(std::ceil(zSize_SI / dx_SI)));
+ WALBERLA_CHECK_FLOAT_EQUAL(real_t(domainSize[0]) * dx_SI, xSize_SI, "domain size in x is not divisible by given dx");
+ WALBERLA_CHECK_FLOAT_EQUAL(real_t(domainSize[1]) * dx_SI, ySize_SI, "domain size in y is not divisible by given dx");
+ WALBERLA_CHECK_FLOAT_EQUAL(real_t(domainSize[2]) * dx_SI, zSize_SI, "domain size in z is not divisible by given dx");
+
+ Vector3< uint_t > cellsPerBlockPerDirection(domainSize[0] / numXBlocks, domainSize[1] / numYBlocks,
+ domainSize[2] / numZBlocks);
+
+ WALBERLA_CHECK_EQUAL(domainSize[0], cellsPerBlockPerDirection[0] * numXBlocks,
+ "number of cells in x of " << domainSize[0]
+ << " is not divisible by given number of blocks in x direction");
+ WALBERLA_CHECK_EQUAL(domainSize[1], cellsPerBlockPerDirection[1] * numYBlocks,
+ "number of cells in y of " << domainSize[1]
+ << " is not divisible by given number of blocks in y direction");
+ WALBERLA_CHECK_EQUAL(domainSize[2], cellsPerBlockPerDirection[2] * numZBlocks,
+ "number of cells in z of " << domainSize[2]
+ << " is not divisible by given number of blocks in z direction");
+
+ WALBERLA_CHECK_GREATER_EQUAL(
+ particleDiameter_SI / dx_SI, 5_r,
+ "Your numerical resolution is below 5 cells per diameter and thus too small for such simulations!");
+
+ const real_t densityRatio = densityParticle_SI / densityFluid_SI;
+ const real_t ReynoldsNumberParticle = uInflow_SI * particleDiameter_SI / kinematicViscosityFluid_SI;
+ const real_t GalileiNumber = std::sqrt((densityRatio - 1_r) * particleDiameter_SI * gravitationalAcceleration_SI) *
+ particleDiameter_SI / kinematicViscosityFluid_SI;
+
+ // in simulation units: dt = 1, dx = 1, densityFluid = 1
+
+ real_t dt_SI = timeStepSize;
+ if(uInflow != real_c(0)){dt_SI = uInflow / uInflow_SI * dx_SI;}
+ const real_t particleDiameter = particleDiameter_SI / dx_SI;
+ const real_t particleGenerationSpacing = particleGenerationSpacing_SI / dx_SI;
+ const real_t viscosity = kinematicViscosityFluid_SI * dt_SI / (dx_SI * dx_SI);
+ const real_t omega = lbm::collision_model::omegaFromViscosity(viscosity);
+ const real_t gravitationalAcceleration = gravitationalAcceleration_SI * dt_SI * dt_SI / dx_SI;
+ const real_t particleVolume = math::pi / 6_r * particleDiameter * particleDiameter * particleDiameter;
+
+ const real_t densityFluid = real_t(1);
+ const real_t densityParticle = densityRatio;
+ const real_t dx = real_t(1);
+
+ const uint_t numTimeSteps = uint_c(std::ceil(runtime_SI / dt_SI));
+ const uint_t infoSpacing = uint_c(std::ceil(infoSpacing_SI / dt_SI));
+ const uint_t vtkSpacingParticles = uint_c(std::ceil(vtkSpacingParticles_SI / dt_SI));
+ const uint_t vtkSpacingFluid = uint_c(std::ceil(vtkSpacingFluid_SI / dt_SI));
+
+ const Vector3< real_t > inflowVec(0_r, 0_r, uInflow);
+
+ const real_t poissonsRatio = real_t(0.22);
+ const real_t kappa = real_t(2) * (real_t(1) - poissonsRatio) / (real_t(2) - poissonsRatio);
+ const real_t particleCollisionTime = collisionTimeFactor * particleDiameter;
+
+
// Necessary Calculations
- const real_t length_conversion = (numXBlocks * real_c(cellsPerBlockPerDirection));
+ const real_t length_conversion = (numXBlocks * real_c(cellsPerBlockPerDirection[0]));
Vector3< uint_t > domainSizeLB(uint_c(length_conversion), uint_c(length_conversion), uint_c(length_conversion));
const real_t uCharacteristicLB = Ma * (std::sqrt(real_c(1 / 3.0)));
@@ -180,11 +372,10 @@ int main(int argc, char** argv)
const real_t delta_T = Thot - Tcold;
const real_t T0 = (Thot + Tcold) / (2 * delta_T);
const real_t kinematicViscosityLB = (Ma / std::sqrt(3)) * std::sqrt(Pr / Ra) * length_conversion;
- // std::min((Ma / std::sqrt(3)) * std::sqrt(Pr / Ra) * length_conversion, std::sqrt(3.6 / Ra));
const real_t thermalDiffusivityLB = kinematicViscosityLB / Pr;
const real_t time_conversion = (length_conversion * length_conversion) / (thermalDiffusivityLB);
- const real_t gravityLB = gravity * time_conversion * time_conversion / length_conversion;
+ const real_t gravityLB = gravitationalAcceleration_SI * time_conversion * time_conversion / length_conversion;
const real_t alphaLB = (uCharacteristicLB * uCharacteristicLB) / (gravityLB * delta_T * length_conversion);
const real_t rho_0 = 1.0;
@@ -202,6 +393,8 @@ int main(int argc, char** argv)
const real_t machnumber = uchar * std::sqrt(3);
const real_t thermal_diffusivity_2 = (Ma / std::sqrt(3)) / (std::sqrt(Ra * Pr)) * length_conversion;
const real_t ratio = length_conversion / thermalDiffusivityLB;
+
+ WALBERLA_LOG_INFO_ON_ROOT("dx_SI is " << dx_SI << " " << "dt_SI is " << dt_SI);
WALBERLA_LOG_INFO_ON_ROOT("length conversion is " << length_conversion << " " << "time conversion is "
<< time_conversion);
WALBERLA_LOG_INFO_ON_ROOT("kinematic viscosity is " << kinematicViscosityLB);
@@ -215,16 +408,16 @@ int main(int argc, char** argv)
WALBERLA_LOG_INFO_ON_ROOT("alpha LB is " << alphaLB);
WALBERLA_LOG_INFO_ON_ROOT("gravity LB is " << gravityLB);
WALBERLA_LOG_INFO_ON_ROOT("Temperature difference delta_t is " << delta_T);
- //WALBERLA_LOG_INFO_ON_ROOT("Number of time steps is " << timeSteps);
+ WALBERLA_LOG_INFO_ON_ROOT("Number of time steps is " << numTimeSteps);
///////////////////////////
// BLOCK STRUCTURE SETUP //
///////////////////////////
- const bool periodicInX = false;
+ const bool periodicInZ = false;
shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid(
- numXBlocks, numYBlocks, numZBlocks, cellsPerBlockPerDirection, cellsPerBlockPerDirection, cellsPerBlockPerDirection, real_t(1), uint_t(0),
+ numXBlocks, numYBlocks, numZBlocks, cellsPerBlockPerDirection[0], cellsPerBlockPerDirection[1], cellsPerBlockPerDirection[2], real_t(1), uint_t(0),
false, false, periodicInX, periodicInY, periodicInZ, // periodicity
false);
@@ -250,19 +443,36 @@ int main(int argc, char** argv)
WALBERLA_CHECK_FLOAT_EQUAL(simulationDomain.yMin(), real_t(0));
WALBERLA_CHECK_FLOAT_EQUAL(simulationDomain.zMin(), real_t(0));
- auto generationDomain = math::AABB::createFromMinMaxCorner(
- math::Vector3< real_t >(simulationDomain.xMax() * (real_t(1) - generationDomainFraction[0]) / real_t(2),
- simulationDomain.yMax() * (real_t(1) - generationDomainFraction[1]) / real_t(2),
- simulationDomain.zMax() * (real_t(1) - generationDomainFraction[2]) / real_t(2)),
- math::Vector3< real_t >(simulationDomain.xMax() * (real_t(1) + generationDomainFraction[0]) / real_t(2),
- simulationDomain.yMax() * (real_t(1) + generationDomainFraction[1]) / real_t(2),
- simulationDomain.zMax() * (real_t(1) + generationDomainFraction[2]) / real_t(2)));
- real_t particleOffset = particleGenerationSpacing / real_t(2);
- /*for (auto pt :
- grid_generator::SCGrid(generationDomain, generationDomain.center(),
- particleGenerationSpacing))*/
- //{
- auto pt = Vector3< real_t >(simulationDomain.xMax()/2,simulationDomain.yMax()/2,simulationDomain.zMax()/2);
+ if(showcase)
+ {
+ auto generationDomain = math::AABB::createFromMinMaxCorner(
+ math::Vector3< real_t >(simulationDomain.xMax() * (real_t(1) - generationDomainFraction[0]) / real_t(2),
+ simulationDomain.yMax() * (real_t(1) - generationDomainFraction[1]) / real_t(2),
+ simulationDomain.zMax() * (real_t(1) - generationDomainFraction[2]) / real_t(2)),
+ math::Vector3< real_t >(simulationDomain.xMax() * (real_t(1) + generationDomainFraction[0]) / real_t(2),
+ simulationDomain.yMax() * (real_t(1) + generationDomainFraction[1]) / real_t(2),
+ simulationDomain.zMax() * (real_t(1) + generationDomainFraction[2]) / real_t(2)));
+ real_t particleOffset = particleGenerationSpacing / real_t(2);
+ uint_t numparticles = 0;
+ for (auto pt : grid_generator::SCGrid(generationDomain, generationDomain.center(), particleGenerationSpacing))
+ {
+ // auto pt = Vector3< real_t >(simulationDomain.xMax()/2,simulationDomain.yMax()/2,simulationDomain.zMax()/2);
+ if (rpdDomain->isContainedInProcessSubdomain(uint_c(mpi::MPIManager::instance()->rank()), pt))
+ {
+ mesa_pd::data::Particle&& p = *ps->create();
+ p.setPosition(pt);
+ p.setInteractionRadius(particleDiameter * real_t(0.5));
+ p.setOwner(mpi::MPIManager::instance()->rank());
+ p.setShapeID(sphereShape);
+ p.setType(0);
+ p.setTemperature(particleTemperature);
+ }
+ numparticles += 1;
+ if (numparticles == 1000) { break; }
+ }
+ }
+ else{
+ auto pt = Vector3< real_t >(simulationDomain.xMax()/2,simulationDomain.yMax()/2,simulationDomain.zMax()/2);
if (rpdDomain->isContainedInProcessSubdomain(uint_c(mpi::MPIManager::instance()->rank()), pt))
{
mesa_pd::data::Particle&& p = *ps->create();
@@ -271,16 +481,12 @@ int main(int argc, char** argv)
p.setOwner(mpi::MPIManager::instance()->rank());
p.setShapeID(sphereShape);
p.setType(0);
- p.setTemperature(math::realRandom(real_t(0.2),real_t(1)));
-
- //p.setTemperature(0.308382);
- //p.setLinearVelocity(Vector3<real_t>(0.1,0.1,0.1));
-
-
+ p.setTemperature(particleTemperature);
}
- //}
+ }
}
+
////////////////////////
// ADD DATA TO BLOCKS //
///////////////////////
@@ -342,8 +548,41 @@ int main(int argc, char** argv)
field::addFlagFieldToStorage< FlagField_T >(blocks, "concentration flag field");
// Synchronize particles between the blocks for the correct mapping of ghost particles
- mesa_pd::mpi::SyncNextNeighbors syncNextNeighborFunc;
- syncNextNeighborFunc(*ps, *rpdDomain);
+ // set up RPD functionality
+ std::function< void(void) > syncCall = [&ps, &rpdDomain]() {
+ // keep overlap for lubrication
+ const real_t overlap = real_t(1.5);
+ mesa_pd::mpi::SyncNextNeighbors syncNextNeighborFunc;
+ syncNextNeighborFunc(*ps, *rpdDomain, overlap);
+ };
+
+ syncCall();
+
+ real_t timeStepSizeRPD = real_t(1) / real_t(numberOfParticleSubCycles);
+ mesa_pd::kernel::VelocityVerletPreForceUpdate vvIntegratorPreForce(timeStepSizeRPD);
+ mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(timeStepSizeRPD);
+ mesa_pd::kernel::LinearSpringDashpot collisionResponse(2);
+ collisionResponse.setFrictionCoefficientDynamic(0, 1, dynamicFrictionCoefficient);
+ collisionResponse.setFrictionCoefficientDynamic(1, 1, dynamicFrictionCoefficient);
+ real_t massSphere = densityParticle * particleVolume;
+ real_t meffSpherePlane = massSphere;
+ real_t meffSphereSphere = massSphere * massSphere / (real_t(2) * massSphere);
+ collisionResponse.setStiffnessAndDamping(0, 1, coefficientOfRestitution, particleCollisionTime, kappa,
+ meffSpherePlane);
+ collisionResponse.setStiffnessAndDamping(1, 1, coefficientOfRestitution, particleCollisionTime, kappa,
+ meffSphereSphere);
+ mesa_pd::kernel::AssocToBlock assoc(blocks->getBlockForestPointer());
+ mesa_pd::mpi::ReduceProperty reduceProperty;
+ mesa_pd::mpi::ReduceContactHistory reduceAndSwapContactHistory;
+
+ // set up coupling functionality
+ Vector3< real_t > gravitationalForce(real_t(0), real_t(0),
+ -(densityParticle - densityFluid) * gravitationalAcceleration * particleVolume);
+ lbm_mesapd_coupling::AddForceOnParticlesKernel addGravitationalForce(gravitationalForce);
+ lbm_mesapd_coupling::ResetHydrodynamicForceTorqueKernel resetHydrodynamicForceTorque;
+ lbm_mesapd_coupling::AverageHydrodynamicForceTorqueKernel averageHydrodynamicForceTorque;
+ lbm_mesapd_coupling::LubricationCorrectionKernel lubricationCorrectionKernel(
+ viscosity, [](real_t r) { return (real_t(0.001 + real_t(0.00007) * r)) * r; });
// Assemble boundary block string
std::string boundariesBlockString = " BoundariesFluid"
@@ -368,8 +607,8 @@ int main(int argc, char** argv)
boundariesBlockString += "\n BoundariesConcentration";
boundariesBlockString += "{"
- "Border { direction W; walldistance -1; flag Density_Concentration_east; }"
- "Border { direction E; walldistance -1; flag Density_Concentration_east; }";
+ "Border { direction W; walldistance -1; flag Neumann_Concentration; }"
+ "Border { direction E; walldistance -1; flag Neumann_Concentration; }";
if (!periodicInY)
{
@@ -519,8 +758,8 @@ int main(int argc, char** argv)
// time loop objects for communication with and without hiding
- SweepTimeloop commTimeloop(blocks->getBlockStorage(), timeSteps);
- SweepTimeloop timeloop(blocks->getBlockStorage(), timeSteps);
+ SweepTimeloop commTimeloop(blocks->getBlockStorage(), numTimeSteps);
+ SweepTimeloop timeloop(blocks->getBlockStorage(), numTimeSteps);
timeloop.addFuncBeforeTimeStep(RemainingTimeLogger(timeloop.getNrOfTimeSteps()), "Remaining Time Logger");
@@ -546,7 +785,7 @@ int main(int argc, char** argv)
#endif
// vtk output
- if (vtkSpacing != uint_t(0))
+ if (vtkSpacingParticles != uint_t(0))
{
// particles
auto particleVtkOutput = make_shared< mesa_pd::vtk::ParticleVtkOutput >(ps);
@@ -559,11 +798,13 @@ int main(int argc, char** argv)
return pIt->getShapeID() == sphereShape &&
!(mesa_pd::data::particle_flags::isSet(pIt->getFlags(), mesa_pd::data::particle_flags::GHOST));
});
- auto particleVtkWriter = vtk::createVTKOutput_PointData(particleVtkOutput, "particles", vtkSpacing, vtkFolder);
+ auto particleVtkWriter = vtk::createVTKOutput_PointData(particleVtkOutput, "particles", vtkSpacingParticles, vtkFolder);
timeloop.addFuncBeforeTimeStep(vtk::writeFiles(particleVtkWriter), "VTK (sphere data)");
+ }
+ if (vtkSpacingFluid != uint_t(0)){
// Fields
- auto vtkOutput_Fluid = vtk::createVTKOutput_BlockData(blocks, "vtk files fluid", vtkSpacing, 0, false, vtkFolder);
+ auto vtkOutput_Fluid = vtk::createVTKOutput_BlockData(blocks, "vtk files fluid", vtkSpacingFluid, 0, false, vtkFolder);
vtkOutput_Fluid->addBeforeFunction(communication_fluid);
@@ -586,7 +827,7 @@ int main(int argc, char** argv)
});
auto vtkOutput_Concentration =
- vtk::createVTKOutput_BlockData(blocks, "vtk files concentration", vtkSpacing, 0, false, vtkFolder);
+ vtk::createVTKOutput_BlockData(blocks, "vtk files concentration", vtkSpacingFluid, 0, false, vtkFolder);
vtkOutput_Concentration->addBeforeFunction(communication_concentration);
@@ -628,8 +869,7 @@ int main(int argc, char** argv)
timeloop.addFuncBeforeTimeStep(vtk::writeFiles(vtkOutput_Fluid), "VTK output Fluid");
timeloop.addFuncBeforeTimeStep(vtk::writeFiles(vtkOutput_Concentration), "VTK output Concentration");
}
-
- if (vtkSpacing != uint_t(0)) { vtk::writeDomainDecomposition(blocks, "domain_decomposition", vtkFolder); }
+ if (vtkSpacingFluid != uint_t(0)) { vtk::writeDomainDecomposition(blocks, "domain_decomposition", vtkFolder); }
////////////////////////////////////////////////////////////////////////////////////////////////
// add LBM communication, boundary handling and the LBM sweeps to the time loop for codegen //
@@ -639,10 +879,6 @@ int main(int argc, char** argv)
particleAndVolumeFractionSoA.particleForcesFieldID, particleAndVolumeFractionSoA.particleVelocitiesFieldID,
pdfFieldFluidCPUGPUID, velFieldFluidCPUGPUID, T0, alphaLB, gravityLB, omega_f, rho_0);
- pystencils::PSMFluidSweepSplit psmFluidSplitSweep(
- particleAndVolumeFractionSoA.BsFieldID, particleAndVolumeFractionSoA.BFieldID, densityConcentrationFieldCPUGPUID,particleAndVolumeFractionSoA.particleForcesFieldID,
- particleAndVolumeFractionSoA.particleVelocitiesFieldID,
- pdfFieldFluidCPUGPUID,velFieldFluidCPUGPUID,T0,alphaLB,gravityLB,omega_f,rho_0,frameWidth);
pystencils::PSMConcentrationSweep psmConcentrationSweep(
particleAndVolumeFractionSoA.BsFieldID, particleAndVolumeFractionSoA.BFieldID, densityConcentrationFieldCPUGPUID,
@@ -655,7 +891,6 @@ int main(int argc, char** argv)
pystencils::LBMFluidSweep lbmFluidSweep(densityConcentrationFieldCPUGPUID,pdfFieldFluidCPUGPUID,velFieldFluidCPUGPUID,T0,alphaLB,gravityLB,omega_f,rho_0);
- pystencils::LBMFluidSplitSweep lbmFluidSplitSweep(densityConcentrationFieldCPUGPUID,pdfFieldFluidCPUGPUID,velFieldFluidCPUGPUID,T0,alphaLB,gravityLB,omega_f,rho_0,frameWidth);
pystencils::LBMConcentrationSweep lbmConcentrationSweep(densityConcentrationFieldCPUGPUID, pdfFieldConcentrationCPUGPUID,
velFieldFluidCPUGPUID,qe,qk);
@@ -663,19 +898,6 @@ int main(int argc, char** argv)
- // Add LBM (fluid and concentration) communication function and boundary handling sweep
- if (useCommunicationHiding)
- {
- timeloop.add() << Sweep(deviceSyncWrapper(noSlip_fluid_bc.getSweep()), "Boundary Handling (No slip fluid)");
- timeloop.add() << Sweep(deviceSyncWrapper(neumann_concentration_bc.getSweep()),
- "Boundary Handling (Concentration Neumann)");
- timeloop.add() << Sweep(deviceSyncWrapper(density_concentration_bc_west.getSweep()),
- "Boundary Handling (Concentration Density west)");
- timeloop.add() << Sweep(deviceSyncWrapper(density_concentration_bc_east.getSweep()),
- "Boundary Handling (Concentration Density east)");
- }
- else
- {
timeloop.add() << BeforeFunction(communication_fluid, "LBM fluid Communication")
<< Sweep(deviceSyncWrapper(noSlip_fluid_bc.getSweep()), "Boundary Handling (No slip fluid)");
timeloop.add() << BeforeFunction(communication_concentration, "LBM concentration Communication")
@@ -685,33 +907,19 @@ int main(int argc, char** argv)
"Boundary Handling (Concentration Density west)");
timeloop.add() << Sweep(deviceSyncWrapper(density_concentration_bc_east.getSweep()),
"Boundary Handling (Concentration Density east)");
- }
+
if (useParticles)
{
-
- timeloop.add() << Sweep(deviceSyncWrapper(psmSweepCollection.particleMappingSweep), "Particle mapping");
- timeloop.add() << Sweep(deviceSyncWrapper(psmSweepCollection.setParticleVelocitiesSweep),
- "Set particle velocities");
- timeloop.add() << Sweep(deviceSyncWrapper(psmFluidSweep), "PSM Fluid sweep");
- //timeloop.add() << Sweep(deviceSyncWrapper(lbmConcentrationSweep), "LBM Concentration sweep");
- timeloop.add() << Sweep(deviceSyncWrapper(psmSweepCollection.setParticleTemperaturesSweep),
- "Set particle temperatures");
- timeloop.add() << Sweep(deviceSyncWrapper(psmConcentrationSweep), "PSM Concentration sweep");
-
- // after both the sweeps, reduce the particle forces.
- timeloop.add() << Sweep(deviceSyncWrapper(psmSweepCollection.reduceParticleForcesSweep),
- "Reduce particle forces");
-
+ addThermalPSMSweepToTimeloop(timeloop, psmSweepCollection, psmFluidSweep,psmConcentrationSweep);
}
else{
- /*WALBERLA_LOG_INFO_ON_ROOT("No particles used so executing normal fluid and concentration sweeps");
- timeloop.add() << Sweep(deviceSyncWrapper(lbmFluidSweep), "LBM Fluid sweep");
- timeloop.add() << Sweep(deviceSyncWrapper(lbmConcentrationSweep), "LBM Concentration sweep");*/
timeloop.add() << Sweep(deviceSyncWrapper(psmSweepCollection.particleMappingSweep), "Particle mapping");
timeloop.add() << Sweep(deviceSyncWrapper(psmSweepCollection.setParticleVelocitiesSweep),
"Set particle velocities");
+ timeloop.add() << Sweep(deviceSyncWrapper(psmSweepCollection.setParticleTemperaturesSweep),
+ "Set particle temperatures");
timeloop.add() << Sweep(deviceSyncWrapper(psmFluidSweep), "PSM Fluid sweep");
//timeloop.add() << Sweep(deviceSyncWrapper(lbmConcentrationSweep), "LBM Concentration sweep");
timeloop.add() << Sweep(deviceSyncWrapper(psmConcentrationSweep), "PSM Concentration sweep");
@@ -730,17 +938,177 @@ int main(int argc, char** argv)
timeloop.addFuncAfterTimeStep(performanceLogger, "Evaluate performance logging");
}
+ ////////////////////////
+ // EXECUTE SIMULATION //
+ ////////////////////////
+
WcTimingPool timeloopTiming;
- // TODO: maybe add warmup phase
+ const bool useOpenMP = true;
+ real_t linkedCellWidth = linkedCellWidthRation * particleDiameter;
+ mesa_pd::data::LinkedCells linkedCells(rpdDomain->getUnionOfLocalAABBs().getExtended(linkedCellWidth),
+ linkedCellWidth);
+ mesa_pd::kernel::InsertParticleIntoLinkedCells ipilc;
- for (uint_t timeStep = 0; timeStep < timeSteps; ++timeStep)
+ // time loop
+ for (uint_t timeStep = 0; timeStep < numTimeSteps; ++timeStep)
{
- if (useCommunicationHiding) { commTimeloop.singleStep(timeloopTiming); }
+ // perform a single simulation step -> this contains LBM and setting of the hydrodynamic interactions
+
timeloop.singleStep(timeloopTiming);
+
+ if(useParticles)
+ {
+ if (particleBarriers) WALBERLA_MPI_BARRIER();
+ timeloopTiming["RPD forEachParticle assoc"].start();
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, assoc, *accessor);
+ if (particleBarriers) WALBERLA_MPI_BARRIER();
+ timeloopTiming["RPD forEachParticle assoc"].end();
+ timeloopTiming["RPD reduceProperty HydrodynamicForceTorqueNotification"].start();
+ reduceProperty.operator()< mesa_pd::HydrodynamicForceTorqueNotification >(*ps);
+ if (particleBarriers) WALBERLA_MPI_BARRIER();
+ timeloopTiming["RPD reduceProperty HydrodynamicForceTorqueNotification"].end();
+
+ if (timeStep == 0)
+ {
+ lbm_mesapd_coupling::InitializeHydrodynamicForceTorqueForAveragingKernel
+ initializeHydrodynamicForceTorqueForAveragingKernel;
+ timeloopTiming["RPD forEachParticle initializeHydrodynamicForceTorqueForAveragingKernel"].start();
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor,
+ initializeHydrodynamicForceTorqueForAveragingKernel, *accessor);
+ if (particleBarriers) WALBERLA_MPI_BARRIER();
+ timeloopTiming["RPD forEachParticle initializeHydrodynamicForceTorqueForAveragingKernel"].end();
+ }
+ timeloopTiming["RPD forEachParticle averageHydrodynamicForceTorque"].start();
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, averageHydrodynamicForceTorque,
+ *accessor);
+ if (particleBarriers) WALBERLA_MPI_BARRIER();
+ timeloopTiming["RPD forEachParticle averageHydrodynamicForceTorque"].end();
+
+ for (auto subCycle = uint_t(0); subCycle < numberOfParticleSubCycles; ++subCycle)
+ {
+ if(useIntegrators)
+ {
+ timeloopTiming["RPD forEachParticle vvIntegratorPreForce"].start();
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, vvIntegratorPreForce,
+ *accessor);
+ if (particleBarriers) WALBERLA_MPI_BARRIER();
+ timeloopTiming["RPD forEachParticle vvIntegratorPreForce"].end();
+ }
+ timeloopTiming["RPD syncCall"].start();
+ syncCall();
+ if (particleBarriers) WALBERLA_MPI_BARRIER();
+ timeloopTiming["RPD syncCall"].end();
+
+ timeloopTiming["RPD linkedCells.clear"].start();
+ linkedCells.clear();
+ if (particleBarriers) WALBERLA_MPI_BARRIER();
+ timeloopTiming["RPD linkedCells.clear"].end();
+ timeloopTiming["RPD forEachParticle ipilc"].start();
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, ipilc, *accessor, linkedCells);
+ if (particleBarriers) WALBERLA_MPI_BARRIER();
+ timeloopTiming["RPD forEachParticle ipilc"].end();
+
+ if (useLubricationForces)
+ {
+ // lubrication correction
+ timeloopTiming["RPD forEachParticlePairHalf lubricationCorrectionKernel"].start();
+ linkedCells.forEachParticlePairHalf(
+ useOpenMP, mesa_pd::kernel::ExcludeInfiniteInfinite(), *accessor,
+ [&lubricationCorrectionKernel, &rpdDomain](const size_t idx1, const size_t idx2, auto& ac) {
+ mesa_pd::collision_detection::AnalyticContactDetection acd;
+ acd.getContactThreshold() = lubricationCorrectionKernel.getNormalCutOffDistance();
+ mesa_pd::kernel::DoubleCast double_cast;
+ mesa_pd::mpi::ContactFilter contact_filter;
+ if (double_cast(idx1, idx2, ac, acd, ac))
+ {
+ if (contact_filter(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), *rpdDomain))
+ {
+ double_cast(acd.getIdx1(), acd.getIdx2(), ac, lubricationCorrectionKernel, ac,
+ acd.getContactNormal(), acd.getPenetrationDepth());
+ }
+ }
+ },
+ *accessor);
+ if (particleBarriers) WALBERLA_MPI_BARRIER();
+ timeloopTiming["RPD forEachParticlePairHalf lubricationCorrectionKernel"].end();
+ }
+
+ // collision response
+ timeloopTiming["RPD forEachParticlePairHalf collisionResponse"].start();
+ linkedCells.forEachParticlePairHalf(
+ useOpenMP, mesa_pd::kernel::ExcludeInfiniteInfinite(), *accessor,
+ [&collisionResponse, &rpdDomain, timeStepSizeRPD](const size_t idx1, const size_t idx2, auto& ac) {
+ mesa_pd::collision_detection::AnalyticContactDetection acd;
+ mesa_pd::kernel::DoubleCast double_cast;
+ mesa_pd::mpi::ContactFilter contact_filter;
+ if (double_cast(idx1, idx2, ac, acd, ac))
+ {
+ if (contact_filter(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), *rpdDomain))
+ {
+ collisionResponse(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(),
+ acd.getContactNormal(), acd.getPenetrationDepth(), timeStepSizeRPD);
+ }
+ }
+ },
+ *accessor);
+ if (particleBarriers) WALBERLA_MPI_BARRIER();
+ timeloopTiming["RPD forEachParticlePairHalf collisionResponse"].end();
+
+ timeloopTiming["RPD reduceProperty reduceAndSwapContactHistory"].start();
+ reduceAndSwapContactHistory(*ps);
+ if (particleBarriers) WALBERLA_MPI_BARRIER();
+ timeloopTiming["RPD reduceProperty reduceAndSwapContactHistory"].end();
+
+ // add hydrodynamic force
+ lbm_mesapd_coupling::AddHydrodynamicInteractionKernel addHydrodynamicInteraction;
+ timeloopTiming["RPD forEachParticle addHydrodynamicInteraction + addGravitationalForce"].start();
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addHydrodynamicInteraction,
+ *accessor);
+
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addGravitationalForce, *accessor);
+ if (particleBarriers) WALBERLA_MPI_BARRIER();
+ timeloopTiming["RPD forEachParticle addHydrodynamicInteraction + addGravitationalForce"].end();
+
+ timeloopTiming["RPD reduceProperty ForceTorqueNotification"].start();
+ reduceProperty.operator()< mesa_pd::ForceTorqueNotification >(*ps);
+ if (particleBarriers) WALBERLA_MPI_BARRIER();
+ timeloopTiming["RPD reduceProperty ForceTorqueNotification"].end();
+
+ if(useIntegrators)
+ {
+ timeloopTiming["RPD forEachParticle vvIntegratorPostForce"].start();
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, vvIntegratorPostForce,
+ *accessor);
+ if (particleBarriers) WALBERLA_MPI_BARRIER();
+ timeloopTiming["RPD forEachParticle vvIntegratorPostForce"].end();
+ }
+ }
+
+ timeloopTiming["RPD syncCall"].start();
+ syncCall();
+ if (particleBarriers) WALBERLA_MPI_BARRIER();
+ timeloopTiming["RPD syncCall"].end();
+
+ timeloopTiming["RPD forEachParticle resetHydrodynamicForceTorque"].start();
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, resetHydrodynamicForceTorque,
+ *accessor);
+ if (particleBarriers) WALBERLA_MPI_BARRIER();
+ timeloopTiming["RPD forEachParticle resetHydrodynamicForceTorque"].end();
+ }
+ if (infoSpacing != 0 && timeStep % infoSpacing == 0 && useParticles == true)
+ {
+ timeloopTiming["Evaluate infos"].start();
+
+ auto particleInfo = evaluateParticleInfo(*accessor);
+ WALBERLA_LOG_INFO_ON_ROOT(particleInfo);
+
+ if (particleBarriers) WALBERLA_MPI_BARRIER();
+ timeloopTiming["Evaluate infos"].end();
+ }
}
+
timeloopTiming.logResultOnRoot();
- auto timeloopTimingReduced = timeloopTiming.getReduced();
return EXIT_SUCCESS;
}