Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
No results found
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 5598 additions and 38 deletions
apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
View file @ 2dee476c
......@@ -80,11 +80,11 @@
#include "mesa_pd/domain/BlockForestDomain.h"
#include "mesa_pd/domain/BlockForestDataHandling.h"
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
#include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/LinearSpringDashpot.h"
#include "mesa_pd/kernel/NonLinearSpringDashpot.h"
#include "mesa_pd/kernel/ParticleSelector.h"
#include "mesa_pd/kernel/VelocityVerletWithShape.h"
#include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/mpi/ClearNextNeighborSync.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h"
......@@ -132,6 +132,10 @@ using ScalarField_T = GhostLayerField< real_t, 1>;
const uint_t FieldGhostLayers = 1;
#define USE_TRTLIKE
//#define USE_D3Q27TRTLIKE
//#define USE_CUMULANT
///////////
// FLAGS //
///////////
......@@ -260,9 +264,7 @@ public:
WALBERLA_MPI_SECTION()
{
mpi::allReduceInplace( pos[0], mpi::SUM );
mpi::allReduceInplace( pos[1], mpi::SUM );
mpi::allReduceInplace( pos[2], mpi::SUM );
mpi::allReduceInplace( pos, mpi::SUM );
mpi::allReduceInplace( transVel[2], mpi::SUM );
}
......@@ -363,6 +365,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
// create bounding planes
mesa_pd::data::Particle&& p0 = *ps->create(true);
p0.setPosition(simulationDomain.minCorner());
p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,1) ));
p0.setOwner(mpi::MPIManager::instance()->rank());
p0.setType(0);
......@@ -374,6 +377,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
//only create top plane when no outflow BC should be set there
mesa_pd::data::Particle&& p1 = *ps->create(true);
p1.setPosition(simulationDomain.maxCorner());
p1.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,-1) ));
p1.setOwner(mpi::MPIManager::instance()->rank());
p1.setType(0);
......@@ -758,10 +762,12 @@ int main( int argc, char **argv )
+ "_mn" + std::to_string(magicNumber);
if(useOmegaBulkAdaption) checkPointFileName +="_omegaBulkAdaption";
if(applyOutflowBCAtTop) checkPointFileName +="_outflowBC";
#ifdef USE_TRT_LIKE_LATTICE_MODEL
checkPointFileName +="_TRTlike";
#else
checkPointFileName += "_other";
#if defined(USE_TRTLIKE)
checkPointFileName += "_trtlike";
#elif defined(USE_CUMULANT)
checkPointFileName += "_cumulant";
#elif defined(USE_D3Q27TRTLIKE)
checkPointFileName += "_d3q27trtlike";
#endif
WALBERLA_LOG_INFO_ON_ROOT("Checkpointing:");
......@@ -852,10 +858,10 @@ int main( int argc, char **argv )
{
WALBERLA_LOG_INFO_ON_ROOT( "Initializing particles from checkpointing file!" );
particleStorageID = blocks->loadBlockData( checkPointFileName + "_mesa.txt", mesa_pd::domain::createBlockForestDataHandling(ps), "Particle Storage" );
} else {
particleStorageID = blocks->addBlockData(mesa_pd::domain::createBlockForestDataHandling(ps), "Particle Storage");
mesa_pd::mpi::ClearNextNeighborSync CNNS;
CNNS(*accessor);
} else {
particleStorageID = blocks->addBlockData(mesa_pd::domain::createBlockForestDataHandling(ps), "Particle Storage");
}
// bounding planes
......@@ -901,14 +907,19 @@ int main( int argc, char **argv )
////////////////////////
// add omega bulk field
BlockDataID omegaBulkFieldID = field::addToStorage<ScalarField_T>( blocks, "omega bulk field", omegaBulk, field::fzyx, uint_t(0) );
BlockDataID omegaBulkFieldID = field::addToStorage<ScalarField_T>( blocks, "omega bulk field", omegaBulk, field::fzyx);
// create the lattice model
real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega );
real_t lambda_d = lbm::collision_model::TRT::lambda_d( omega, magicNumber );
#ifdef WALBERLA_BUILD_WITH_CODEGEN
#ifdef USE_CUMULANT
WALBERLA_LOG_INFO_ON_ROOT("Using generated cumulant lattice model!");
LatticeModel_T latticeModel = LatticeModel_T(omega);
#elif defined(USE_TRTLIKE) || defined(USE_D3Q27TRTLIKE)
WALBERLA_LOG_INFO_ON_ROOT("Using generated TRT-like lattice model!");
LatticeModel_T latticeModel = LatticeModel_T(omegaBulkFieldID, lambda_d, lambda_e);
#endif
#else
WALBERLA_LOG_INFO_ON_ROOT("Using waLBerla built-in MRT lattice model and ignoring omega bulk field since not supported!");
LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::D3Q19MRT( omegaBulk, omegaBulk, lambda_d, lambda_e, lambda_e, lambda_d ));
......@@ -962,9 +973,9 @@ int main( int argc, char **argv )
syncCall();
real_t timeStepSizeRPD = real_t(1)/real_t(numRPDSubCycles);
mesa_pd::kernel::ExplicitEulerWithShape explEulerIntegrator(timeStepSizeRPD);
mesa_pd::kernel::VelocityVerletWithShapePreForceUpdate vvIntegratorPreForce(timeStepSizeRPD);
mesa_pd::kernel::VelocityVerletWithShapePostForceUpdate vvIntegratorPostForce(timeStepSizeRPD);
mesa_pd::kernel::ExplicitEuler explEulerIntegrator(timeStepSizeRPD);
mesa_pd::kernel::VelocityVerletPreForceUpdate vvIntegratorPreForce(timeStepSizeRPD);
mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(timeStepSizeRPD);
// linear model
mesa_pd::kernel::LinearSpringDashpot linearCollisionResponse(1);
......@@ -1117,11 +1128,7 @@ int main( int argc, char **argv )
}else if( reconstructorType == "Ext")
{
auto sphereNormalExtrapolationDirectionFinder = make_shared<lbm_mesapd_coupling::SphereNormalExtrapolationDirectionFinder>(blocks);
#ifdef USE_TRT_LIKE_LATTICE_MODEL
auto extrapolationReconstructor = lbm_mesapd_coupling::makeExtrapolationReconstructor<BoundaryHandling_T, lbm_mesapd_coupling::SphereNormalExtrapolationDirectionFinder, true>(blocks, boundaryHandlingID, sphereNormalExtrapolationDirectionFinder, uint_t(3), true);
#else
auto extrapolationReconstructor = lbm_mesapd_coupling::makeExtrapolationReconstructor<BoundaryHandling_T, lbm_mesapd_coupling::SphereNormalExtrapolationDirectionFinder, false>(blocks, boundaryHandlingID, sphereNormalExtrapolationDirectionFinder, uint_t(3), true);
#endif
timeloopAfterParticles.add() << BeforeFunction( fullPDFCommunicationScheme, "LBM Communication" )
<< Sweep( makeSharedSweep(lbm_mesapd_coupling::makePdfReconstructionManager<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, FormerMO_Flag, Fluid_Flag, extrapolationReconstructor, conserveMomentum)), "PDF Restore" );
} else
......@@ -1308,9 +1315,8 @@ int main( int argc, char **argv )
}
else
{
lbm_mesapd_coupling::RegularParticlesSelector sphereSelector;
lbm_mesapd_coupling::AddHydrodynamicInteractionKernel addHydrodynamicInteraction;
ps->forEachParticle(useOpenMP, sphereSelector, *accessor, addHydrodynamicInteraction, *accessor );
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addHydrodynamicInteraction, *accessor );
}
hydForce = getForce(sphereUid,*accessor) - lubForce - collisionForce;
......
apps/benchmarks/FluidParticleCouplingWithLoadBalancing/CMakeLists.txt 0 → 100644
View file @ 2dee476c
waLBerla_add_executable( NAME FluidParticleWorkloadEvaluation FILES FluidParticleWorkloadEvaluation.cpp DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk )
waLBerla_add_executable( NAME FluidParticleWorkloadDistribution FILES FluidParticleWorkloadDistribution.cpp DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk )
apps/benchmarks/FluidParticleCouplingWithLoadBalancing/FluidParticleWorkloadDistribution.cpp 0 → 100644
View file @ 2dee476c
//======================================================================================================================
//
// 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 FluidParticleWorkloadDistribution.cpp
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#include "blockforest/Initialization.h"
#include "blockforest/communication/UniformBufferedScheme.h"
#include "blockforest/loadbalancing/all.h"
#include "blockforest/AABBRefinementSelection.h"
#include "boundary/all.h"
#include "core/DataTypes.h"
#include "core/Environment.h"
#include "core/debug/Debug.h"
#include "core/debug/TestSubsystem.h"
#include "core/math/all.h"
#include "core/timing/RemainingTimeLogger.h"
#include "core/mpi/Broadcast.h"
#include "domain_decomposition/SharedSweep.h"
#include "domain_decomposition/BlockSweepWrapper.h"
#include "field/AddToStorage.h"
#include "field/StabilityChecker.h"
#include "field/communication/PackInfo.h"
#include "lbm/boundary/all.h"
#include "lbm/communication/PdfFieldPackInfo.h"
#include "lbm/field/AddToStorage.h"
#include "lbm/field/PdfField.h"
#include "lbm/field/VelocityFieldWriter.h"
#include "lbm/lattice_model/D3Q19.h"
#include "lbm/sweeps/CellwiseSweep.h"
#include "lbm_mesapd_coupling/amr/BlockInfo.h"
#include "lbm_mesapd_coupling/amr/InfoCollection.h"
#include "lbm_mesapd_coupling/amr/weight_assignment/WeightEvaluationFunctions.h"
#include "lbm_mesapd_coupling/amr/weight_assignment/WeightAssignmentFunctor.h"
#include "lbm_mesapd_coupling/amr/weight_assignment/MetisAssignmentFunctor.h"
#include "lbm_mesapd_coupling/mapping/ParticleMapping.h"
#include "lbm_mesapd_coupling/momentum_exchange_method/MovingParticleMapping.h"
#include "lbm_mesapd_coupling/momentum_exchange_method/boundary/CurvedLinear.h"
#include "lbm_mesapd_coupling/momentum_exchange_method/reconstruction/Reconstructor.h"
#include "lbm_mesapd_coupling/momentum_exchange_method/reconstruction/PdfReconstructionManager.h"
#include "lbm_mesapd_coupling/utility/AddForceOnParticlesKernel.h"
#include "lbm_mesapd_coupling/utility/ParticleSelector.h"
#include "lbm_mesapd_coupling/DataTypes.h"
#include "lbm_mesapd_coupling/utility/AverageHydrodynamicForceTorqueKernel.h"
#include "lbm_mesapd_coupling/utility/AddHydrodynamicInteractionKernel.h"
#include "lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
#include "lbm_mesapd_coupling/utility/LubricationCorrectionKernel.h"
#include "lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h"
#include "lbm_mesapd_coupling/utility/InspectionProbe.h"
#include "mesa_pd/collision_detection/AnalyticContactDetection.h"
#include "mesa_pd/data/ParticleAccessorWithShape.h"
#include "mesa_pd/data/ParticleStorage.h"
#include "mesa_pd/data/ShapeStorage.h"
#include "mesa_pd/data/DataTypes.h"
#include "mesa_pd/data/shape/HalfSpace.h"
#include "mesa_pd/data/shape/Sphere.h"
#include "mesa_pd/domain/BlockForestDomain.h"
#include "mesa_pd/domain/BlockForestDataHandling.h"
#include "mesa_pd/kernel/AssocToBlock.h"
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/LinearSpringDashpot.h"
#include "mesa_pd/kernel/ParticleSelector.h"
#include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/mpi/ClearNextNeighborSync.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/SyncNextNeighborsBlockForest.h"
#include "mesa_pd/mpi/ReduceProperty.h"
#include "mesa_pd/mpi/ReduceContactHistory.h"
#include "mesa_pd/mpi/ContactFilter.h"
#include "mesa_pd/mpi/notifications/ForceTorqueNotification.h"
#include "mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h"
#include "mesa_pd/vtk/ParticleVtkOutput.h"
#include "timeloop/SweepTimeloop.h"
#include "vtk/all.h"
#include "field/vtk/all.h"
#include "lbm/vtk/all.h"
#include "Utility.h"
namespace fluid_particle_workload_distribution
{
///////////
// USING //
///////////
using namespace walberla;
using walberla::uint_t;
using LatticeModel_T = lbm::D3Q19< lbm::collision_model::TRT, false >;
using Stencil_T = LatticeModel_T::Stencil;
using PdfField_T = lbm::PdfField<LatticeModel_T>;
using flag_t = walberla::uint8_t;
using FlagField_T = FlagField<flag_t>;
using ParticleField_T = lbm_mesapd_coupling::ParticleField_T;
using ParticleAccessor_T = mesa_pd::data::ParticleAccessorWithShape;
const uint_t FieldGhostLayers = 1;
using NoSlip_T = lbm::NoSlip<LatticeModel_T, flag_t> ;
using MO_T = lbm_mesapd_coupling::CurvedLinear< LatticeModel_T, FlagField_T, ParticleAccessor_T >;
using BoundaryHandling_T = BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T, MO_T >;
///////////
// FLAGS //
///////////
const FlagUID Fluid_Flag( "fluid" );
const FlagUID NoSlip_Flag( "no slip" );
const FlagUID MO_Flag( "moving obstacle" );
const FlagUID FormerMO_Flag( "former moving obstacle" );
///////////////////////////////////
// LOAD BALANCING FUNCTIONALITY //
//////////////////////////////////
/////////////////////
// BLOCK STRUCTURE //
/////////////////////
static void workloadAndMemoryAssignment( SetupBlockForest& forest )
{
for (auto &block : forest) {
block.setWorkload( numeric_cast< workload_t >( uint_t(1) << block.getLevel() ) );
block.setMemory( numeric_cast< memory_t >(1) );
}
}
static shared_ptr< StructuredBlockForest > createBlockStructure( const AABB & domainAABB, Vector3<uint_t> blockSizeInCells,
bool useBox, const std::string & loadDistributionStrategy,
bool keepGlobalBlockInformation = false )
{
SetupBlockForest sforest;
Vector3<uint_t> numberOfBlocksPerDirection( uint_c(domainAABB.size(0)) / blockSizeInCells[0],
uint_c(domainAABB.size(1)) / blockSizeInCells[1],
uint_c(domainAABB.size(2)) / blockSizeInCells[2] );
for(uint_t i = 0; i < 3; ++i )
{
WALBERLA_CHECK_EQUAL( numberOfBlocksPerDirection[i] * blockSizeInCells[i], uint_c(domainAABB.size(i)),
"Domain can not be decomposed in direction " << i << " into fine blocks of size " << blockSizeInCells[i] );
}
sforest.addWorkloadMemorySUIDAssignmentFunction( workloadAndMemoryAssignment );
Vector3<bool> periodicity( true, true, false);
if( useBox )
{
periodicity[0] = false;
periodicity[1] = false;
}
sforest.init( domainAABB,
numberOfBlocksPerDirection[0], numberOfBlocksPerDirection[1], numberOfBlocksPerDirection[2],
periodicity[0], periodicity[1], periodicity[2]);
// calculate process distribution
const memory_t memoryLimit = math::Limits< memory_t >::inf();
if( loadDistributionStrategy == "Hilbert" )
{
bool useHilbert = true;
sforest.balanceLoad( blockforest::StaticLevelwiseCurveBalance(useHilbert), uint_c( MPIManager::instance()->numProcesses() ), real_t(0), memoryLimit, true );
} else if ( loadDistributionStrategy == "Morton" )
{
bool useHilbert = false;
sforest.balanceLoad( blockforest::StaticLevelwiseCurveBalance(useHilbert), uint_c( MPIManager::instance()->numProcesses() ), real_t(0), memoryLimit, true );
} else if ( loadDistributionStrategy == "ParMetis" )
{
blockforest::StaticLevelwiseParMetis::Algorithm algorithm = blockforest::StaticLevelwiseParMetis::Algorithm::PARMETIS_PART_GEOM_KWAY;
blockforest::StaticLevelwiseParMetis staticParMetis(algorithm);
sforest.balanceLoad( staticParMetis, uint_c( MPIManager::instance()->numProcesses() ), real_t(0), memoryLimit, true );
} else if (loadDistributionStrategy == "Diffusive" )
{
// also use Hilbert curve here
bool useHilbert = true;
sforest.balanceLoad( blockforest::StaticLevelwiseCurveBalance(useHilbert), uint_c( MPIManager::instance()->numProcesses() ), real_t(0), memoryLimit, true );
} else
{
WALBERLA_ABORT("Load distribution strategy \"" << loadDistributionStrategy << "\t not implemented! - Aborting" );
}
WALBERLA_LOG_INFO_ON_ROOT( sforest );
// create StructuredBlockForest (encapsulates a newly created BlockForest)
shared_ptr< StructuredBlockForest > sbf =
make_shared< StructuredBlockForest >( make_shared< BlockForest >( uint_c( MPIManager::instance()->rank() ), sforest, keepGlobalBlockInformation ),
blockSizeInCells[0], blockSizeInCells[1], blockSizeInCells[2]);
sbf->createCellBoundingBoxes();
return sbf;
}
/////////////////////////////////////
// BOUNDARY HANDLING CUSTOMIZATION //
/////////////////////////////////////
class MyBoundaryHandling : public blockforest::AlwaysInitializeBlockDataHandling< BoundaryHandling_T >
{
public:
MyBoundaryHandling( const weak_ptr< StructuredBlockStorage > & blocks,
const BlockDataID & flagFieldID, const BlockDataID & pdfFieldID,
const BlockDataID & particleFieldID, const shared_ptr<ParticleAccessor_T>& ac ) :
blocks_( blocks ), flagFieldID_( flagFieldID ), pdfFieldID_( pdfFieldID ), particleFieldID_( particleFieldID ), ac_( ac )
{}
BoundaryHandling_T * initialize( IBlock * const block ) override;
private:
weak_ptr< StructuredBlockStorage > blocks_;
const BlockDataID flagFieldID_;
const BlockDataID pdfFieldID_;
const BlockDataID particleFieldID_;
shared_ptr<ParticleAccessor_T> ac_;
}; // class MyBoundaryHandling
BoundaryHandling_T * MyBoundaryHandling::initialize( IBlock * const block )
{
WALBERLA_ASSERT_NOT_NULLPTR( block );
auto * flagField = block->getData< FlagField_T >( flagFieldID_ );
auto * pdfField = block->getData< PdfField_T > ( pdfFieldID_ );
auto* particleField = block->getData<lbm_mesapd_coupling::ParticleField_T>(particleFieldID_);
const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
auto blocksPtr = blocks_.lock();
WALBERLA_CHECK_NOT_NULLPTR( blocksPtr );
BoundaryHandling_T * handling = new BoundaryHandling_T( "moving obstacle boundary handling", flagField, fluid,
NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
MO_T( "MO", MO_Flag, pdfField, flagField, particleField, ac_, fluid, *blocksPtr, *block ),
BoundaryHandling_T::Mode::ENTIRE_FIELD_TRAVERSAL);
handling->fillWithDomain( FieldGhostLayers );
return handling;
}
void clearBoundaryHandling( BlockForest & forest, const BlockDataID & boundaryHandlingID ) {
for( auto blockIt = forest.begin(); blockIt != forest.end(); ++blockIt )
{
BoundaryHandling_T * boundaryHandling = blockIt->getData<BoundaryHandling_T>(boundaryHandlingID);
boundaryHandling->clear( FieldGhostLayers );
}
}
void recreateBoundaryHandling( BlockForest & forest, const BlockDataID & boundaryHandlingID ) {
for( auto blockIt = forest.begin(); blockIt != forest.end(); ++blockIt )
{
BoundaryHandling_T * boundaryHandling = blockIt->getData<BoundaryHandling_T>(boundaryHandlingID);
boundaryHandling->fillWithDomain( FieldGhostLayers );
}
}
void clearParticleField( BlockForest & forest, const BlockDataID & particleFieldID, const ParticleAccessor_T & accessor ) {
for (auto blockIt = forest.begin(); blockIt != forest.end(); ++blockIt) {
ParticleField_T *particleField = blockIt->getData<ParticleField_T>(particleFieldID);
particleField->setWithGhostLayer(accessor.getInvalidUid());
}
}
//*******************************************************************************************************************
void evaluateTimers(WcTimingPool & timingPool,
const std::vector<std::vector<std::string> > & timerKeys,
std::vector<real_t> & timings )
{
for (auto & timingsIt : timings)
{
timingsIt = real_t(0);
}
timingPool.unifyRegisteredTimersAcrossProcesses();
for (auto i = uint_t(0); i < timerKeys.size(); ++i )
{
auto keys = timerKeys[i];
for (const auto &timerName : keys)
{
if(timingPool.timerExists(timerName))
{
timings[i] += real_c(timingPool[timerName].total());
}
}
}
}
uint_t evaluateEdgeCut(BlockForest & forest)
{
//note: only works for edges in uniform grids
auto edgecut = uint_t(0); // = edge weights between processes
for( auto blockIt = forest.begin(); blockIt != forest.end(); ++blockIt )
{
auto * block = static_cast<blockforest::Block*> (&(*blockIt));
real_t blockVolume = block->getAABB().volume();
real_t approximateEdgeLength = std::cbrt( blockVolume );
uint_t faceNeighborWeight = uint_c(approximateEdgeLength * approximateEdgeLength ); //common face
uint_t edgeNeighborWeight = uint_c(approximateEdgeLength); //common edge
uint_t cornerNeighborWeight = uint_c( 1 ); //common corner
for( const uint_t idx : blockforest::getFaceNeighborhoodSectionIndices() )
{
for (auto nb = uint_t(0); nb < block->getNeighborhoodSectionSize(idx); ++nb)
{
if( block->neighborExistsRemotely(idx,nb) ) edgecut += faceNeighborWeight;
}
}
for( const uint_t idx : blockforest::getEdgeNeighborhoodSectionIndices() )
{
for (auto nb = uint_t(0); nb < block->getNeighborhoodSectionSize(idx); ++nb)
{
if( block->neighborExistsRemotely(idx,nb) ) edgecut += edgeNeighborWeight;
}
}
for( const uint_t idx : blockforest::getCornerNeighborhoodSectionIndices() )
{
for (auto nb = uint_t(0); nb < block->getNeighborhoodSectionSize(idx); ++nb)
{
if( block->neighborExistsRemotely(idx,nb) ) edgecut += cornerNeighborWeight;
}
}
}
return edgecut;
}
void evaluateTotalSimulationTimePassed(WcTimingPool & timeloopTimingPool, const std::string & simulationString, const std::string & loadbalancingString,
double & totalSimTime, double & totalLBTime)
{
shared_ptr< WcTimingPool> reduced = timeloopTimingPool.getReduced(timing::REDUCE_TOTAL, 0);
double totalTime = 0.0;
WALBERLA_ROOT_SECTION(){
totalTime = (*reduced)[simulationString].total();
}
totalSimTime = totalTime;
double lbTime = 0.0;
WALBERLA_ROOT_SECTION(){
lbTime = (*reduced)[loadbalancingString].total();
}
totalLBTime = lbTime;
}
void createPlane( const shared_ptr<mesa_pd::data::ParticleStorage> & ps, const shared_ptr<mesa_pd::data::ShapeStorage> & ss,
const Vector3<real_t> position, const 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.getNormalized() ));
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);
WALBERLA_LOG_INFO_ON_ROOT("Created plane at position " << position << " with normal " << normal.getNormalized ());
}
void createBasicPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, const shared_ptr<mesa_pd::data::ShapeStorage> & ss, const math::AABB & simulationDomain)
{
createPlane(ps, ss, simulationDomain.minCorner(), Vector3<real_t>(0,0, 1));
//createPlane(ps, ss, Vector3<real_t>(simulationDomain.xMax() * 0.5, simulationDomain.yMax() * 0.5, simulationDomain.zMin() + std::max(simulationDomain.xMax(), simulationDomain.yMax()) * 0.5 ), Vector3<real_t>(0,1, 1));
createPlane(ps, ss, simulationDomain.maxCorner(), Vector3<real_t>(0,0,-1));
}
void addBoxPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, const shared_ptr<mesa_pd::data::ShapeStorage> & ss, const math::AABB & simulationDomain)
{
// add bounding planes (four horizontal directions)
createPlane(ps, ss, simulationDomain.minCorner(), Vector3<real_t>( 1, 0,0));
createPlane(ps, ss, simulationDomain.maxCorner(), Vector3<real_t>(-1, 0,0));
createPlane(ps, ss, simulationDomain.minCorner(), Vector3<real_t>( 0, 1,0));
createPlane(ps, ss, simulationDomain.maxCorner(), Vector3<real_t>( 0,-1,0));
}
void addHopperPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, const shared_ptr<mesa_pd::data::ShapeStorage> & ss, const math::AABB & simulationDomain,
real_t hopperRelHeight, real_t hopperRelOpening)
{
//hopper planes
real_t xMax = simulationDomain.xMax();
real_t yMax = simulationDomain.yMax();
real_t zMax = simulationDomain.zMax();
Vector3<real_t> p1(0,0,hopperRelHeight*zMax);
Vector3<real_t> n1(p1[2],0,hopperRelOpening*xMax-p1[0]);
Vector3<real_t> n2(0,p1[2],hopperRelOpening*yMax-p1[0]);
createPlane(ps, ss, p1, n1);
createPlane(ps, ss, p1, n2);
Vector3<real_t> p2(xMax,yMax,hopperRelHeight*zMax);
Vector3<real_t> n3(-p2[2],0,-((real_t(1)-hopperRelOpening)*xMax-p2[0]));
Vector3<real_t> n4(0,-p2[2],-((real_t(1)-hopperRelOpening)*yMax-p2[1]));
createPlane(ps, ss, p2, n3);
createPlane(ps, ss, p2, n4);
}
void evaluateParticleSimulation(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, const shared_ptr<mesa_pd::data::ShapeStorage> & ss,
uint_t numParticles, uint_t numGlobalParticles)
{
auto numShapes = ss->shapes.size();
uint_t numLocalParticles = 0;
//uint_t numGhostParticles = 0;
uint_t numGlobalParticlesOfRank = 0;
for(auto p = ps->begin(); p != ps->end(); ++p)
{
using namespace walberla::mesa_pd::data::particle_flags;
if (isSet(p->getFlags(), GHOST))
{
//++numGhostParticles;
} else if (isSet(p->getFlags(), GLOBAL))
{
++numGlobalParticlesOfRank;
} else
{
++numLocalParticles;
}
if(p->getShapeID() >= numShapes)
{
WALBERLA_LOG_INFO("Found invalid shape id " << *p);
}
}
//WALBERLA_LOG_INFO(numShapes << " " << numLocalParticles << " " << numGhostParticles);
if(numGlobalParticlesOfRank != numGlobalParticles)
{
WALBERLA_LOG_INFO("Number of global particles has changed to " << numGlobalParticlesOfRank);
}
mpi::reduceInplace(numLocalParticles, mpi::SUM);
if(numLocalParticles != numParticles)
{
WALBERLA_LOG_INFO_ON_ROOT("Number of particles has changed to " << numLocalParticles);
}
}
bool isnan(Vector3<real_t> vec)
{
return std::isnan(vec[0]) || std::isnan(vec[1]) || std::isnan(vec[2]);
}
void checkParticleProperties(const shared_ptr<mesa_pd::data::ParticleStorage> & ps)
{
for(auto p = ps->begin(); p != ps->end(); ++p)
{
if(isnan(p->getHydrodynamicForce())) WALBERLA_LOG_INFO("Found nan in hyd Force " << *p);
if(isnan(p->getOldHydrodynamicForce())) WALBERLA_LOG_INFO("Found nan in old hyd Force " << *p);
if(isnan(p->getForce())) WALBERLA_LOG_INFO("Found nan in Force " << *p);
if(isnan(p->getOldForce())) WALBERLA_LOG_INFO("Found nan in Force " << *p);
}
}
template< typename Probe_T>
void checkMapping(const shared_ptr< StructuredBlockStorage > & blocks, const BlockDataID pdfFieldID,
const BlockDataID boundaryHandlingID, Probe_T & probe )
{
for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
{
auto * pdfField = blockIt->getData< PdfField_T >( pdfFieldID );
auto * boundaryHandling = blockIt->getData< BoundaryHandling_T >( boundaryHandlingID );
WALBERLA_FOR_ALL_CELLS_XYZ(pdfField,
if (boundaryHandling->isDomain(x, y, z))
{
uint_t f = uint_t(0);
//for( uint_t f = uint_t(0); f < PdfField_T::F_SIZE; ++f )
//{
if( !walberla::field::internal::stabilityCheckerIsFinite( pdfField->get( x, y, z, cell_idx_c(f) ) ) )
{
Vector3< real_t > center;
blocks->getBlockLocalCellCenter( *blockIt, Cell(x,y,z), center );
Cell gCell(x,y,z);
blocks->transformBlockLocalToGlobalCell( gCell, *blockIt);
WALBERLA_LOG_INFO("Instability found in block local cell( " << x << ", " << y << ", " << z << " ) at index " << f
<< " = global cell ( " << gCell.x() << ", " << gCell.y() << ", " << gCell.z() << " ) with cell center ( " << center[0] << ", " << center[1] << ", " << center[2] << " )");
probe.setPosition(center);
real_t rho;
Vector3<real_t> velocity;
probe(rho, velocity);
}
//}
}
);
}
}
//*******************************************************************************************************************
/*!\brief Simulation of settling particles inside a rectangular column filled with viscous fluid
*
* This application is used in the paper
* Rettinger, Ruede - "Dynamic Load Balancing Techniques for Particulate Flow Simulations", 2019, Computation
* in Section 4 to apply the load estimator and to evaluate different load distribution strategies.
* It has here been adapted to the new mesapd and its coupling.
* More infos can be found in the PhD thesis by Christoph Rettinger.
*
* It, however, features several different command line arguments that can be used to tweak the simulation.
* The setup can be horizontally period, a box or a hopper geometry (configurable, as in the paper).
* The size, resolution and used blocks for the domain partitioning can be changed.
* Initially, all particles are pushed upwards to obtain a dense packing at the top plane.
*
* Most importantly, the load balancing can be modified:
* - load estimation strategies:
* - pure LBM = number of cells per block = constant workload per block
* - coupling based load estimator = use fitted function from Sec. 3 of paper
* - load distribution strategies:
* - space-filling curves: Hilbert and Morton
* - ParMETIS (and several algorithms and parameters)
* - diffusive (and options)
* - load balancing frequency
*/
//*******************************************************************************************************************
int main( int argc, char **argv )
{
debug::enterTestMode();
mpi::Environment env( argc, argv );
MPIManager::instance()->useWorldComm();
///////////////////
// Customization //
///////////////////
// simulation control
bool shortRun = false;
bool funcTest = false;
bool fileIO = true;
bool checkSimulation = false;
uint_t vtkWriteFreqDD = 0; //domain decomposition
uint_t vtkWriteFreqPa = 0; //particles
uint_t vtkWriteFreqFl = 0; //fluid
uint_t vtkWriteFreq = 0; //general
uint_t vtkWriteFreqInit = 0; //initial (particle-only) simulation
std::string baseFolder = "vtk_out_WorkloadDistribution"; // folder for vtk and file output
bool useProgressLogging = false;
// physical setup
auto GalileoNumber = real_t(50);
auto densityRatio = real_t(1.5);
auto diameter = real_t(15);
auto solidVolumeFraction = real_t(0.1);
auto blockSize = uint_t(32);
auto XBlocks = uint_t(12);
auto YBlocks = uint_t(12);
auto ZBlocks = uint_t(16);
bool useBox = false;
bool useHopper = false;
auto hopperRelHeight = real_t(0.5); // for hopper setup
auto hopperRelOpening = real_t(0.3); // for hopper setup
auto timesteps = uint_t(80000);
//numerical parameters
auto loadBalancingCheckFrequency = uint_t(100);
auto numRPDSubCycles = uint_t(10);
bool useBlockForestSync = false;
// load balancing
std::string loadEvaluationStrategy = "LBM"; //LBM, Fit
std::string loadDistributionStrategy = "Hilbert"; //Morton, Hilbert, ParMetis, Diffusive
real_t blockBaseWeight = real_t(1);
auto parMetis_ipc2redist = real_t(1000);
auto parMetisTolerance = real_t(-1);
std::string parMetisAlgorithmString = "ADAPTIVE_REPART";
auto diffusionFlowIterations = uint_t(15);
auto diffusionMaxIterations = uint_t(20);
bool useNoSlipForPlanes = false;
for( int i = 1; i < argc; ++i )
{
if( std::strcmp( argv[i], "--shortRun" ) == 0 ) { shortRun = true; continue; }
if( std::strcmp( argv[i], "--funcTest" ) == 0 ) { funcTest = true; continue; }
if( std::strcmp( argv[i], "--fileIO" ) == 0 ) { fileIO = true; continue; }
if( std::strcmp( argv[i], "--useProgressLogging" ) == 0 ) { useProgressLogging = true; continue; }
if( std::strcmp( argv[i], "--checkSimulation" ) == 0 ) { checkSimulation = true; continue; }
if( std::strcmp( argv[i], "--vtkWriteFreqDD" ) == 0 ) { vtkWriteFreqDD = uint_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--vtkWriteFreqPa" ) == 0 ) { vtkWriteFreqPa = uint_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--vtkWriteFreqFl" ) == 0 ) { vtkWriteFreqFl = uint_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--vtkWriteFreq" ) == 0 ) { vtkWriteFreq = uint_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--vtkWriteFreqInit" ) == 0 ) { vtkWriteFreqInit = uint_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--baseFolder" ) == 0 ) { baseFolder = argv[++i]; continue; }
if( std::strcmp( argv[i], "--densityRatio" ) == 0 ) { densityRatio = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--Ga" ) == 0 ) { GalileoNumber = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--diameter" ) == 0 ) { diameter = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--blockSize" ) == 0 ) { blockSize = uint_c(std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--XBlocks" ) == 0 ) { XBlocks = uint_c(std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--YBlocks" ) == 0 ) { YBlocks = uint_c(std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--ZBlocks" ) == 0 ) { ZBlocks = uint_c(std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--useBox" ) == 0 ) { useBox = true; continue; }
if( std::strcmp( argv[i], "--useHopper" ) == 0 ) { useHopper = true; continue; }
if( std::strcmp( argv[i], "--hopperHeight" ) == 0 ) { hopperRelHeight = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--hopperOpening" ) == 0 ) { hopperRelOpening = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--timesteps" ) == 0 ) { timesteps = uint_c(std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--numRPDSubCycles" ) == 0 ) { numRPDSubCycles = uint_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--useBlockForestSync" ) == 0 ) { useBlockForestSync = true; continue; }
if( std::strcmp( argv[i], "--useNoSlipForPlanes" ) == 0 ) { useNoSlipForPlanes = true; continue; }
if( std::strcmp( argv[i], "--blockBaseWeight" ) == 0 ) { blockBaseWeight = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--loadBalancingCheckFrequency" ) == 0 ) { loadBalancingCheckFrequency = uint_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--loadEvaluationStrategy" ) == 0 ) { loadEvaluationStrategy = argv[++i]; continue; }
if( std::strcmp( argv[i], "--loadDistributionStrategy" ) == 0 ) { loadDistributionStrategy = argv[++i]; continue; }
if( std::strcmp( argv[i], "--ipc2redist" ) == 0 ) { parMetis_ipc2redist = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--parMetisTolerance" ) == 0 ) { parMetisTolerance = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--parMetisAlgorithm" ) == 0 ) { parMetisAlgorithmString = argv[++i]; continue; }
if( std::strcmp( argv[i], "--diffusionFlowIterations" ) == 0 ) { diffusionFlowIterations = uint_c(std::atof(argv[++i])); continue; }
if( std::strcmp( argv[i], "--diffusionMaxIterations" ) == 0 ) { diffusionMaxIterations = uint_c(std::atof(argv[++i])); continue; }
WALBERLA_ABORT("Unrecognized command line argument found: " << argv[i]);
}
if( fileIO )
{
WALBERLA_ROOT_SECTION(){
// create base directory if it does not yet exist
filesystem::path tpath( baseFolder );
if( !filesystem::exists( tpath ) )
filesystem::create_directory( tpath );
}
WALBERLA_MPI_BARRIER();
}
if(useProgressLogging) logging::Logging::instance()->includeLoggingToFile(baseFolder + "/progress_logging.txt");
if( loadEvaluationStrategy != "LBM" && loadEvaluationStrategy != "Fit" && loadEvaluationStrategy != "FitMulti")
{
WALBERLA_ABORT("Invalid load evaluation strategy: " << loadEvaluationStrategy);
}
if( vtkWriteFreq != 0 )
{
vtkWriteFreqDD = vtkWriteFreq;
vtkWriteFreqPa = vtkWriteFreq;
vtkWriteFreqFl = vtkWriteFreq;
}
if( diameter > real_c(blockSize) )
{
WALBERLA_LOG_WARNING("Particle Synchronization might not work since bodies are large compared to block size!");
}
if( useHopper )
{
WALBERLA_CHECK(hopperRelHeight >= real_t(0) && hopperRelHeight <= real_t(1), "Invalid relative hopper height of " << hopperRelHeight);
WALBERLA_CHECK(hopperRelOpening >= real_t(0) && hopperRelOpening <= real_t(1), "Invalid relative hopper opening of " << hopperRelOpening);
}
//////////////////////////
// NUMERICAL PARAMETERS //
//////////////////////////
const Vector3<uint_t> domainSize( XBlocks * blockSize, YBlocks * blockSize, ZBlocks * blockSize );
const auto domainVolume = real_t(domainSize[0] * domainSize[1] * domainSize[2]);
const real_t sphereVolume = math::pi / real_t(6) * diameter * diameter * diameter;
const uint_t numberOfSediments = uint_c(std::ceil(solidVolumeFraction * domainVolume / sphereVolume));
real_t expectedSedimentVolumeFraction = (useBox||useHopper) ? real_t(0.45) : real_t(0.52);
const real_t expectedSedimentedVolume = real_t(1)/expectedSedimentVolumeFraction * real_c(numberOfSediments) * sphereVolume;
const real_t expectedSedimentedHeight = std::max(diameter, expectedSedimentedVolume / real_c(domainSize[0] * domainSize[1]));
const auto uRef = real_t(0.02);
const real_t xRef = diameter;
const real_t tRef = xRef / uRef;
const real_t gravitationalAcceleration = uRef * uRef / ( (densityRatio-real_t(1)) * diameter );
const real_t viscosity = uRef * diameter / GalileoNumber;
const real_t omega = lbm::collision_model::omegaFromViscosity(viscosity);
const real_t tau = real_t(1) / omega;
const auto dx = real_t(1);
const real_t overlap = real_t( 1.5 ) * dx;
timesteps = funcTest ? 1 : ( shortRun ? uint_t(100) : timesteps );
WALBERLA_LOG_INFO_ON_ROOT("Setup (in simulation, i.e. lattice, units):");
WALBERLA_LOG_INFO_ON_ROOT(" - domain size = " << domainSize);
WALBERLA_LOG_INFO_ON_ROOT(" - sediment diameter = " << diameter );
WALBERLA_LOG_INFO_ON_ROOT(" - Galileo number = " << GalileoNumber );
WALBERLA_LOG_INFO_ON_ROOT(" - number of sediments: " << numberOfSediments);
WALBERLA_LOG_INFO_ON_ROOT(" - densityRatio = " << densityRatio );
WALBERLA_LOG_INFO_ON_ROOT(" - fluid: relaxation time (tau) = " << tau << ", kin. visc = " << viscosity );
WALBERLA_LOG_INFO_ON_ROOT(" - gravitational acceleration = " << gravitationalAcceleration );
WALBERLA_LOG_INFO_ON_ROOT(" - reference values: x = " << xRef << ", t = " << tRef << ", vel = " << uRef);
WALBERLA_LOG_INFO_ON_ROOT(" - omega: " << omega);
if( vtkWriteFreqDD > 0 )
{
WALBERLA_LOG_INFO_ON_ROOT(" - writing vtk files of domain decomposition to folder \"" << baseFolder << "\" with frequency " << vtkWriteFreqDD);
}
if( vtkWriteFreqPa > 0 )
{
WALBERLA_LOG_INFO_ON_ROOT(" - writing vtk files of bodies data to folder \"" << baseFolder << "\" with frequency " << vtkWriteFreqPa);
}
if( vtkWriteFreqFl > 0 )
{
WALBERLA_LOG_INFO_ON_ROOT(" - writing vtk files of fluid data to folder \"" << baseFolder << "\" with frequency " << vtkWriteFreqFl);
}
if( useBox )
{
WALBERLA_LOG_INFO_ON_ROOT(" - using box setup");
}
else if ( useHopper )
{
WALBERLA_LOG_INFO_ON_ROOT(" - using hopper setup");
}
else
{
WALBERLA_LOG_INFO_ON_ROOT(" - using horizontally periodic domain");
}
WALBERLA_LOG_INFO_ON_ROOT(" - refinement / load balancing check frequency: " << loadBalancingCheckFrequency);
WALBERLA_LOG_INFO_ON_ROOT(" - load evaluation strategy: " << loadEvaluationStrategy);
WALBERLA_LOG_INFO_ON_ROOT(" - load distribution strategy: " << loadDistributionStrategy);
///////////////////////////
// BLOCK STRUCTURE SETUP //
///////////////////////////
Vector3<uint_t> blockSizeInCells( blockSize );
AABB simulationDomain( real_t(0), real_t(0), real_t(0), real_c(domainSize[0]), real_c(domainSize[1]), real_c(domainSize[2]) );
AABB sedimentDomain( real_t(0), real_t(0), real_c(domainSize[2]) - expectedSedimentedHeight, real_c(domainSize[0]), real_c(domainSize[1]), real_c(domainSize[2]) );
auto blocks = createBlockStructure( simulationDomain, blockSizeInCells, (useBox||useHopper), loadDistributionStrategy );
//write initial domain decomposition to file
if( vtkWriteFreqDD > 0 )
{
vtk::writeDomainDecomposition( blocks, "initial_domain_decomposition", baseFolder );
}
/////////
// RPD //
/////////
const real_t restitutionCoeff = real_t(0.97);
const real_t frictionCoeffStatic = real_t(0.8);
const real_t frictionCoeffDynamic = real_t(0.15);
const real_t collisionTime = real_t(4) * diameter; // from my paper
const real_t poissonsRatio = real_t(0.22);
const real_t kappa = real_t(2) * ( real_t(1) - poissonsRatio ) / ( real_t(2) - poissonsRatio ) ;
auto rpdDomain = std::make_shared<mesa_pd::domain::BlockForestDomain>(blocks->getBlockForestPointer());
//init data structures
auto ps = walberla::make_shared<mesa_pd::data::ParticleStorage>(1);
blocks->addBlockData(mesa_pd::domain::createBlockForestDataHandling(ps), "Particle Storage"); // returned ID is not used, but ps has to be known to blockforest
auto ss = walberla::make_shared<mesa_pd::data::ShapeStorage>();
auto accessor = walberla::make_shared<ParticleAccessor_T >(ps, ss);
real_t timeStepSizeRPD = real_t(1)/real_t(numRPDSubCycles);
mesa_pd::kernel::VelocityVerletPreForceUpdate vvIntegratorPreForce(timeStepSizeRPD);
mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(timeStepSizeRPD);
// types: 0 = wall, 1: sphere
mesa_pd::kernel::LinearSpringDashpot collisionResponse(2);
collisionResponse.setFrictionCoefficientDynamic(0,1,frictionCoeffDynamic);
collisionResponse.setFrictionCoefficientDynamic(1,1,frictionCoeffDynamic);
collisionResponse.setFrictionCoefficientStatic(0,1,frictionCoeffStatic);
collisionResponse.setFrictionCoefficientStatic(1,1,frictionCoeffStatic);
const real_t particleMass = densityRatio * sphereVolume;
const real_t effMass_sphereWall = particleMass;
const real_t effMass_sphereSphere = particleMass * particleMass / ( real_t(2) * particleMass );
collisionResponse.setStiffnessAndDamping(0,1,restitutionCoeff,collisionTime,kappa,effMass_sphereWall);
collisionResponse.setStiffnessAndDamping(1,1,restitutionCoeff,collisionTime,kappa,effMass_sphereSphere);
mesa_pd::mpi::ReduceProperty reduceProperty;
mesa_pd::mpi::ReduceContactHistory reduceAndSwapContactHistory;
mesa_pd::mpi::SyncNextNeighbors syncNextNeighborFunc;
mesa_pd::kernel::AssocToBlock associateToBlock(blocks->getBlockForestPointer());
mesa_pd::mpi::SyncNextNeighborsBlockForest syncNextNeighborBlockForestFunc;
// create bounding planes
createBasicPlaneSetup(ps, ss, simulationDomain);
if(useBox || useHopper) addBoxPlaneSetup(ps, ss, simulationDomain);
if(useHopper) addHopperPlaneSetup(ps, ss, simulationDomain, hopperRelHeight, hopperRelOpening);
auto numGlobalParticles = ps->size();
WALBERLA_LOG_INFO_ON_ROOT("Created " << numGlobalParticles << " global particles");
// add the sediments
WALBERLA_LOG_INFO_ON_ROOT("Starting creation of sediments");
AABB sedimentGenerationDomain( real_t(0), real_t(0), real_t(0.5)*real_c(domainSize[2]), real_c(domainSize[0]), real_c(domainSize[1]), real_c(domainSize[2]) );
auto xParticle = real_t(0);
auto yParticle = real_t(0);
auto zParticle = real_t(0);
auto rank = mpi::MPIManager::instance()->rank();
auto sphereShape = ss->create<mesa_pd::data::Sphere>( diameter * real_t(0.5) );
ss->shapes[sphereShape]->updateMassAndInertia(densityRatio);
std::mt19937 randomGenerator (static_cast<unsigned int>(2610)); // fixed seed: quasi-random and reproducible
for( uint_t nSed = 0; nSed < numberOfSediments; ++nSed )
{
WALBERLA_ROOT_SECTION()
{
xParticle = math::realRandom<real_t>(sedimentGenerationDomain.xMin(), sedimentGenerationDomain.xMax(),randomGenerator);
yParticle = math::realRandom<real_t>(sedimentGenerationDomain.yMin(), sedimentGenerationDomain.yMax(),randomGenerator);
zParticle = math::realRandom<real_t>(sedimentGenerationDomain.zMin(), sedimentGenerationDomain.zMax(),randomGenerator);
}
WALBERLA_MPI_SECTION()
{
mpi::broadcastObject( xParticle );
mpi::broadcastObject( yParticle );
mpi::broadcastObject( zParticle );
}
auto position = Vector3<real_t>( xParticle, yParticle, zParticle );
if (!rpdDomain->isContainedInProcessSubdomain(uint_c(rank), position)) continue;
auto p = ps->create();
p->setPosition(position);
p->setInteractionRadius(diameter * real_t(0.5));
p->setShapeID(sphereShape);
p->setType(1);
p->setOwner(rank);
}
if(useBlockForestSync)
{
ps->forEachParticle(false, mesa_pd::kernel::SelectLocal(), *accessor, associateToBlock, *accessor);
syncNextNeighborBlockForestFunc(*ps, blocks->getBlockForestPointer(), rpdDomain, overlap);
} else
{
syncNextNeighborFunc(*ps, *rpdDomain, overlap);
}
// Carry out particle-only simulation to obtain dense packing at top plane
// consists of three phases:
// 1: carefully resolve initial overlaps due to random generation, no gravity
// 2: apply low gravity in positive z-direction to create packing until convergence or targeted packing height reached
// 3: carry out a few time steps with gravity in negative direction to relay the system towards the real setup
const bool useOpenMP = false;
const real_t dt_RPD_Init = real_t(1);
const auto particleSimStepsPhase1 = uint_t(1000);
const auto maxParticleSimStepsPhase2 = (shortRun) ? uint_t(10) : uint_t(200000);
const auto particleSimStepsPhase3 = uint_t(std::sqrt(real_t(2)/std::fabs(gravitationalAcceleration)));
uint_t maxInitialParticleSimSteps = particleSimStepsPhase1 + maxParticleSimStepsPhase2 + particleSimStepsPhase3;
auto particleVtkOutput = make_shared<mesa_pd::vtk::ParticleVtkOutput>(ps);
particleVtkOutput->addOutput<mesa_pd::data::SelectParticleLinearVelocity>("velocity");
particleVtkOutput->addOutput<mesa_pd::data::SelectParticleOwner>("owner");
particleVtkOutput->setParticleSelector( [sphereShape](const mesa_pd::data::ParticleStorage::iterator& pIt) {return !mesa_pd::data::particle_flags::isSet(pIt->getFlags(), mesa_pd::data::particle_flags::GHOST) && pIt->getShapeID() == sphereShape;} ); //limit output to local sphere
auto particleVtkWriterInit = vtk::createVTKOutput_PointData(particleVtkOutput, "Particles_init", 1, baseFolder, "simulation_step");
real_t gravitationalAccelerationGeneration = gravitationalAcceleration;
auto oldMinParticlePosition = real_t(0);
real_t phase2ConvergenceLimit = std::fabs(gravitationalAccelerationGeneration);
real_t heightConvergenceThreshold = sedimentDomain.zMin();
uint_t beginOfPhase3SimStep = uint_t(0);
uint_t currentPhase = 1;
for(auto pet = uint_t(0); pet <= maxInitialParticleSimSteps; ++pet )
{
real_t maxPenetrationDepth = 0;
ps->forEachParticlePairHalf(useOpenMP, mesa_pd::kernel::ExcludeInfiniteInfinite(), *accessor,
[&collisionResponse, &rpdDomain, &maxPenetrationDepth, dt_RPD_Init]
(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))
{
maxPenetrationDepth = std::max(maxPenetrationDepth, std::abs(acd.getPenetrationDepth()));
collisionResponse(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(),
acd.getContactNormal(), acd.getPenetrationDepth(), dt_RPD_Init);
}
}
},
*accessor );
reduceAndSwapContactHistory(*ps);
mpi::allReduceInplace(maxPenetrationDepth, mpi::MAX);
reduceProperty.operator()<mesa_pd::ForceTorqueNotification>(*ps);
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, mesa_pd::kernel::ExplicitEuler(dt_RPD_Init), *accessor);
if(useBlockForestSync)
{
syncNextNeighborBlockForestFunc(*ps, blocks->getBlockForestPointer(), rpdDomain, overlap);
ps->forEachParticle(false, mesa_pd::kernel::SelectLocal(), *accessor, associateToBlock, *accessor);
} else
{
syncNextNeighborFunc(*ps, *rpdDomain, overlap);
}
if( vtkWriteFreqInit > uint_t(0) && pet % vtkWriteFreqInit == uint_t(0) )
{
particleVtkWriterInit->write();
}
if(currentPhase == 1)
{
// damp velocities to avoid too large ones
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor,
[](const size_t idx, ParticleAccessor_T& ac){
ac.setLinearVelocity(idx, ac.getLinearVelocity(idx) * real_t(0.5));
ac.setAngularVelocity(idx, ac.getAngularVelocity(idx) * real_t(0.5));
}, *accessor);
if(pet > particleSimStepsPhase1)
{
WALBERLA_LOG_INFO_ON_ROOT("Starting phase 2 of initial particle simulation, with height threshold = " << heightConvergenceThreshold);
currentPhase = 2;
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor,
[](const size_t idx, ParticleAccessor_T& ac){
ac.setLinearVelocity(idx, Vector3<real_t>(0.0));
ac.setAngularVelocity(idx, Vector3<real_t>(0.0));
}, *accessor);
}
} else if(currentPhase == 2)
{
Vector3<real_t> gravitationalForce( real_t(0), real_t(0), (densityRatio - real_t(1)) * gravitationalAccelerationGeneration * sphereVolume );
lbm_mesapd_coupling::AddForceOnParticlesKernel addGravitationalForce(gravitationalForce);
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addGravitationalForce, *accessor );
real_t minParticlePosition = sedimentGenerationDomain.zMax();
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor,
[&minParticlePosition](const size_t idx, ParticleAccessor_T& ac){
minParticlePosition = std::min(ac.getPosition(idx)[2], minParticlePosition);
}, *accessor);
WALBERLA_MPI_SECTION()
{
mpi::allReduceInplace(minParticlePosition, mpi::MIN);
}
WALBERLA_ROOT_SECTION()
{
if( pet % 100 == 0)
{
WALBERLA_LOG_INFO("[" << pet << "] Min position of all particles = " << minParticlePosition << " with goal height " << heightConvergenceThreshold);
}
}
if( minParticlePosition > heightConvergenceThreshold ) currentPhase = 3;
if( pet % 500 == 0)
{
if( std::fabs(minParticlePosition - oldMinParticlePosition) / minParticlePosition < phase2ConvergenceLimit ) currentPhase = 3;
oldMinParticlePosition = minParticlePosition;
}
if( currentPhase == 3)
{
WALBERLA_LOG_INFO_ON_ROOT("Starting phase 3 of initial particle simulation");
beginOfPhase3SimStep = pet;
}
} else if(currentPhase == 3)
{
Vector3<real_t> gravitationalForce( real_t(0), real_t(0), -(densityRatio - real_t(1)) * gravitationalAccelerationGeneration * sphereVolume );
lbm_mesapd_coupling::AddForceOnParticlesKernel addGravitationalForce(gravitationalForce);
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addGravitationalForce, *accessor );
if(pet - beginOfPhase3SimStep > particleSimStepsPhase3)
{
Vector3<real_t> initialParticleVelocity(real_t(0));
WALBERLA_LOG_INFO_ON_ROOT("Setting initial velocity " << initialParticleVelocity << " of all particles");
// reset velocities to avoid too large ones
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor,
[initialParticleVelocity](const size_t idx, ParticleAccessor_T& ac){
ac.setLinearVelocity(idx, initialParticleVelocity);
ac.setAngularVelocity(idx, Vector3<real_t>(0));
}, *accessor);
break;
}
}
}
WALBERLA_LOG_INFO_ON_ROOT("Sediment layer creation done!");
///////////////////////
// ADD DATA TO BLOCKS //
////////////////////////
// create the lattice model
LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( omega, lbm::collision_model::TRT::threeSixteenth ) );
// add PDF field
BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (fzyx)", latticeModel,
Vector3< real_t >( real_t(0) ), real_t(1),
FieldGhostLayers, field::fzyx );
// add flag field
BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>( blocks, "flag field", FieldGhostLayers );
// add particle field
BlockDataID particleFieldID = field::addToStorage<lbm_mesapd_coupling::ParticleField_T>( blocks, "particle field", accessor->getInvalidUid(), field::fzyx, FieldGhostLayers );
// add boundary handling & initialize outer domain boundaries
BlockDataID boundaryHandlingID = blocks->addBlockData( make_shared< MyBoundaryHandling >( blocks, flagFieldID, pdfFieldID, particleFieldID, accessor ), "boundary handling" );
Vector3<real_t> gravitationalForce( real_t(0), real_t(0), -(densityRatio - real_t(1)) * gravitationalAcceleration * sphereVolume );
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;});
lbm_mesapd_coupling::ParticleMappingKernel<BoundaryHandling_T> particleMappingKernel(blocks, boundaryHandlingID);
lbm_mesapd_coupling::MovingParticleMappingKernel<BoundaryHandling_T> movingParticleMappingKernel(blocks, boundaryHandlingID, particleFieldID);
lbm_mesapd_coupling::ParticleMappingKernel<BoundaryHandling_T> fixedParticleMappingKernel(blocks, boundaryHandlingID);
lbm_mesapd_coupling::InitializeHydrodynamicForceTorqueForAveragingKernel initializeHydrodynamicForceTorqueForAveragingKernel;
WALBERLA_LOG_INFO_ON_ROOT(" - Lubrication correction:");
WALBERLA_LOG_INFO_ON_ROOT(" - normal cut off distance = " << lubricationCorrectionKernel.getNormalCutOffDistance());
WALBERLA_LOG_INFO_ON_ROOT(" - tangential translational cut off distance = " << lubricationCorrectionKernel.getTangentialTranslationalCutOffDistance());
WALBERLA_LOG_INFO_ON_ROOT(" - tangential rotational cut off distance = " << lubricationCorrectionKernel.getTangentialRotationalCutOffDistance());
const real_t maximumLubricationCutOffDistance = std::max(lubricationCorrectionKernel.getNormalCutOffDistance(), std::max(lubricationCorrectionKernel.getTangentialRotationalCutOffDistance(), lubricationCorrectionKernel.getTangentialTranslationalCutOffDistance()));
std::function<void(void)> particleMappingCall = [ps, accessor, &movingParticleMappingKernel, &fixedParticleMappingKernel, useNoSlipForPlanes]()
{
// map planes into the LBM simulation -> act as no-slip boundaries
if(useNoSlipForPlanes)
{
ps->forEachParticle(false, lbm_mesapd_coupling::GlobalParticlesSelector(), *accessor, fixedParticleMappingKernel, *accessor, NoSlip_Flag);
}
else {
ps->forEachParticle(false, lbm_mesapd_coupling::GlobalParticlesSelector(), *accessor, movingParticleMappingKernel, *accessor, MO_Flag);
}
// map particles into the LBM simulation
ps->forEachParticle(false, lbm_mesapd_coupling::RegularParticlesSelector(), *accessor, movingParticleMappingKernel, *accessor, MO_Flag);
};
particleMappingCall();
//////////////////////////
// LOAD BALANCING UTILs //
//////////////////////////
auto & blockforest = blocks->getBlockForest();
blockforest.recalculateBlockLevelsInRefresh( false ); // = only load balancing, no refinement checking
blockforest.alwaysRebalanceInRefresh( true ); //load balancing every time refresh is triggered
blockforest.reevaluateMinTargetLevelsAfterForcedRefinement( false );
blockforest.allowRefreshChangingDepth( false );
blockforest.allowMultipleRefreshCycles( false ); // otherwise info collections are invalid
shared_ptr<lbm_mesapd_coupling::amr::InfoCollection> couplingInfoCollection = walberla::make_shared<lbm_mesapd_coupling::amr::InfoCollection>();
uint_t numberOfProcesses = uint_c(MPIManager::instance()->numProcesses());
if( loadDistributionStrategy == "Hilbert" || loadDistributionStrategy == "Morton")
{
bool curveAllGather = true;
bool balanceLevelwise = true;
if( loadDistributionStrategy == "Hilbert")
{
bool useHilbert = true;
blockforest.setRefreshPhantomBlockMigrationPreparationFunction( blockforest::DynamicCurveBalance< blockforest::PODPhantomWeight<real_t> >( useHilbert, curveAllGather, balanceLevelwise ) );
}
else if (loadDistributionStrategy == "Morton" )
{
bool useHilbert = false;
blockforest.setRefreshPhantomBlockMigrationPreparationFunction( blockforest::DynamicCurveBalance< blockforest::PODPhantomWeight<real_t> >( useHilbert, curveAllGather, balanceLevelwise ) );
}
if( loadEvaluationStrategy == "Fit" )
{
lbm_mesapd_coupling::amr::WeightEvaluationFunctor weightEvaluationFunctor(couplingInfoCollection, lbm_mesapd_coupling::amr::fittedTotalWeightEvaluationFunction);
lbm_mesapd_coupling::amr::WeightAssignmentFunctor weightAssignmentFunctor(weightEvaluationFunctor);
weightAssignmentFunctor.setBlockBaseWeight(blockBaseWeight);
blockforest.setRefreshPhantomBlockDataAssignmentFunction(weightAssignmentFunctor);
}
else if( loadEvaluationStrategy == "LBM" )
{
lbm_mesapd_coupling::amr::WeightAssignmentFunctor weightAssignmentFunctor(lbm_mesapd_coupling::amr::defaultWeightEvaluationFunction);
weightAssignmentFunctor.setBlockBaseWeight(blockBaseWeight);
blockforest.setRefreshPhantomBlockDataAssignmentFunction(weightAssignmentFunctor);
}
else
{
WALBERLA_ABORT("Invalid load evaluation strategy: " << loadEvaluationStrategy);
}
blockforest.setRefreshPhantomBlockDataPackFunction(blockforest::PODPhantomWeightPackUnpack<real_t>());
blockforest.setRefreshPhantomBlockDataUnpackFunction(blockforest::PODPhantomWeightPackUnpack<real_t>());
}
else if( loadDistributionStrategy == "ParMetis")
{
#ifndef WALBERLA_BUILD_WITH_PARMETIS
WALBERLA_ABORT( "You are trying to use ParMetis functionality but waLBerla is not configured to use it. Set 'WALBERLA_BUILD_WITH_PARMETIS' to 'ON' in your CMake cache to build against an installed version of ParMetis!" );
#endif
blockforest::DynamicParMetis::Algorithm parMetisAlgorithm = blockforest::DynamicParMetis::stringToAlgorithm(parMetisAlgorithmString);
blockforest::DynamicParMetis::WeightsToUse parMetisWeightsToUse = blockforest::DynamicParMetis::WeightsToUse::PARMETIS_BOTH_WEIGHTS;
blockforest::DynamicParMetis::EdgeSource parMetisEdgeSource = blockforest::DynamicParMetis::EdgeSource::PARMETIS_EDGES_FROM_EDGE_WEIGHTS;
blockforest::DynamicParMetis dynamicParMetis(parMetisAlgorithm, parMetisWeightsToUse, parMetisEdgeSource);
dynamicParMetis.setipc2redist(parMetis_ipc2redist);
auto numberOfBlocks = XBlocks * YBlocks * ZBlocks;
real_t loadImbalanceTolerance = (parMetisTolerance < real_t(1)) ? std::max(real_t(1.05), real_t(1) + real_t(1) / ( real_c(numberOfBlocks) / real_c(numberOfProcesses) ) ) : parMetisTolerance;
dynamicParMetis.setImbalanceTolerance(double(loadImbalanceTolerance), 0);
WALBERLA_LOG_INFO_ON_ROOT(" - ParMetis configuration: ");
WALBERLA_LOG_INFO_ON_ROOT(" - algorithm = " << dynamicParMetis.algorithmToString() );
WALBERLA_LOG_INFO_ON_ROOT(" - weights to use = " << dynamicParMetis.weightsToUseToString() );
WALBERLA_LOG_INFO_ON_ROOT(" - edge source = " << dynamicParMetis.edgeSourceToString() );
WALBERLA_LOG_INFO_ON_ROOT(" - ipc2redist parameter = " << dynamicParMetis.getipc2redist() );
WALBERLA_LOG_INFO_ON_ROOT(" - imbalance tolerance = " << dynamicParMetis.getImbalanceTolerance() );
blockforest.setRefreshPhantomBlockDataPackFunction(blockforest::DynamicParMetisBlockInfoPackUnpack());
blockforest.setRefreshPhantomBlockDataUnpackFunction(blockforest::DynamicParMetisBlockInfoPackUnpack());
blockforest.setRefreshPhantomBlockMigrationPreparationFunction( dynamicParMetis );
if( loadEvaluationStrategy == "Fit" )
{
lbm_mesapd_coupling::amr::WeightEvaluationFunctor weightEvaluationFunctor(couplingInfoCollection, lbm_mesapd_coupling::amr::fittedTotalWeightEvaluationFunction);
lbm_mesapd_coupling::amr::MetisAssignmentFunctor weightAssignmentFunctor(weightEvaluationFunctor); //attention: special METIS assignment functor!
weightAssignmentFunctor.setBlockBaseWeight(blockBaseWeight);
real_t weightMultiplicator = real_t(1000); // values from predictor are in range [0-5] which is too coarse when cast to int as done in parmetis
weightAssignmentFunctor.setWeightMultiplicator(weightMultiplicator);
blockforest.setRefreshPhantomBlockDataAssignmentFunction(weightAssignmentFunctor);
}
else if( loadEvaluationStrategy == "LBM" )
{
lbm_mesapd_coupling::amr::MetisAssignmentFunctor weightAssignmentFunctor(lbm_mesapd_coupling::amr::defaultWeightEvaluationFunction);
weightAssignmentFunctor.setBlockBaseWeight(blockBaseWeight);
blockforest.setRefreshPhantomBlockDataAssignmentFunction(weightAssignmentFunctor);
}
else
{
WALBERLA_ABORT("Invalid load evaluation strategy: " << loadEvaluationStrategy);
}
}
else if( loadDistributionStrategy == "Diffusive")
{
using DB_T = blockforest::DynamicDiffusionBalance< blockforest::PODPhantomWeight<real_t> >;
DB_T dynamicDiffusion(diffusionMaxIterations, diffusionFlowIterations );
dynamicDiffusion.setMode(DB_T::Mode::DIFFUSION_PUSH);
WALBERLA_LOG_INFO_ON_ROOT(" - Dynamic diffusion configuration: ");
WALBERLA_LOG_INFO_ON_ROOT(" - max iterations = " << dynamicDiffusion.getMaxIterations() );
WALBERLA_LOG_INFO_ON_ROOT(" - flow iterations = " << dynamicDiffusion.getFlowIterations());
blockforest.setRefreshPhantomBlockDataPackFunction(blockforest::PODPhantomWeightPackUnpack<real_t>());
blockforest.setRefreshPhantomBlockDataUnpackFunction(blockforest::PODPhantomWeightPackUnpack<real_t>());
blockforest.setRefreshPhantomBlockMigrationPreparationFunction( dynamicDiffusion );
if( loadEvaluationStrategy == "Fit" )
{
lbm_mesapd_coupling::amr::WeightEvaluationFunctor weightEvaluationFunctor(couplingInfoCollection, lbm_mesapd_coupling::amr::fittedTotalWeightEvaluationFunction);
lbm_mesapd_coupling::amr::WeightAssignmentFunctor weightAssignmentFunctor(weightEvaluationFunctor);
weightAssignmentFunctor.setBlockBaseWeight(blockBaseWeight);
blockforest.setRefreshPhantomBlockDataAssignmentFunction(weightAssignmentFunctor);
}
else if( loadEvaluationStrategy == "LBM" )
{
lbm_mesapd_coupling::amr::WeightAssignmentFunctor weightAssignmentFunctor(lbm_mesapd_coupling::amr::defaultWeightEvaluationFunction);
weightAssignmentFunctor.setBlockBaseWeight(blockBaseWeight);
blockforest.setRefreshPhantomBlockDataAssignmentFunction(weightAssignmentFunctor);
}
else
{
WALBERLA_ABORT("Invalid load evaluation strategy: " << loadEvaluationStrategy);
}
} else
{
WALBERLA_ABORT("Load distribution strategy \"" << loadDistributionStrategy << "\t not implemented! - Aborting" );
}
lbm_mesapd_coupling::amr::BlockInfo emptyExampleBlock(blockSize*blockSize*blockSize, 0, 0, 0, 0, numRPDSubCycles);
WALBERLA_LOG_INFO_ON_ROOT("An example empty block has the weight: " << lbm_mesapd_coupling::amr::fittedTotalWeightEvaluationFunction(emptyExampleBlock));
///////////////
// TIME LOOP //
///////////////
// create the timeloop
SweepTimeloop timeloop( blocks->getBlockStorage(), timesteps );
timeloop.addFuncBeforeTimeStep( RemainingTimeLogger( timeloop.getNrOfTimeSteps() ), "Remaining Time Logger" );
if( vtkWriteFreqPa != uint_t(0) ) {
auto particleVtkWriter = vtk::createVTKOutput_PointData(particleVtkOutput, "Particles", vtkWriteFreqPa, baseFolder, "simulation_step");
timeloop.addFuncBeforeTimeStep( vtk::writeFiles( particleVtkWriter ), "VTK (sphere data)" );
}
if( vtkWriteFreqFl != uint_t(0) ) {
// pdf field
auto pdfFieldVTK = vtk::createVTKOutput_BlockData(blocks, "fluid_field", vtkWriteFreqFl, 0, false, baseFolder);
field::FlagFieldCellFilter<FlagField_T> fluidFilter(flagFieldID);
fluidFilter.addFlag(Fluid_Flag);
pdfFieldVTK->addCellInclusionFilter(fluidFilter);
pdfFieldVTK->addCellDataWriter(
make_shared<lbm::VelocityVTKWriter<LatticeModel_T, float> >(pdfFieldID, "VelocityFromPDF"));
pdfFieldVTK->addCellDataWriter(
make_shared<lbm::DensityVTKWriter<LatticeModel_T, float> >(pdfFieldID, "DensityFromPDF"));
timeloop.addFuncBeforeTimeStep(vtk::writeFiles(pdfFieldVTK), "VTK (fluid field data)");
}
if( vtkWriteFreqDD != uint_t(0) ) {
auto domainDecompVTK = vtk::createVTKOutput_DomainDecomposition(blocks, "domain_decomposition", vtkWriteFreqDD, baseFolder );
timeloop.addFuncBeforeTimeStep( vtk::writeFiles(domainDecompVTK), "VTK (domain decomposition)");
}
// setup of the LBM communication for synchronizing the pdf field between neighboring blocks
blockforest::communication::UniformBufferedScheme< Stencil_T > optimizedPDFCommunicationScheme( blocks );
optimizedPDFCommunicationScheme.addPackInfo( make_shared< lbm::PdfFieldPackInfo< LatticeModel_T > >( pdfFieldID ) ); // optimized sync
auto sweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag );
// add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
timeloop.add() << BeforeFunction( optimizedPDFCommunicationScheme, "LBM Communication" )
<< Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingID ), "Boundary Handling" );
// streaming & collide
timeloop.add() << Sweep( makeSharedSweep(sweep), "Stream&Collide" );
SweepTimeloop timeloopAfterParticles( blocks->getBlockStorage(), timesteps );
// sweep for updating the particle mapping into the LBM simulation
bool conserveMomentum = false;
timeloopAfterParticles.add() << Sweep( lbm_mesapd_coupling::makeMovingParticleMapping<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, MO_Flag, FormerMO_Flag,
lbm_mesapd_coupling::RegularParticlesSelector(), conserveMomentum), "Particle Mapping" );
bool recomputeTargetDensity = false;
auto gradReconstructor = lbm_mesapd_coupling::makeGradsMomentApproximationReconstructor<BoundaryHandling_T>(blocks, boundaryHandlingID, omega, recomputeTargetDensity,true);
blockforest::communication::UniformBufferedScheme< Stencil_T > fullPDFCommunicationScheme( blocks );
fullPDFCommunicationScheme.addPackInfo( make_shared< field::communication::PackInfo< PdfField_T > >( pdfFieldID ) ); // full sync
timeloopAfterParticles.add() << BeforeFunction( fullPDFCommunicationScheme, "PDF Communication" )
<< Sweep( makeSharedSweep(lbm_mesapd_coupling::makePdfReconstructionManager<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, FormerMO_Flag, Fluid_Flag,
gradReconstructor, conserveMomentum) ), "PDF Restore" );
////////////////////////
// EXECUTE SIMULATION //
////////////////////////
uint_t loadEvaluationFrequency = loadBalancingCheckFrequency;
// file for simulation infos
std::string infoFileName( baseFolder + "/simulation_info.txt");
WALBERLA_ROOT_SECTION()
{
std::ofstream file;
file.open( infoFileName.c_str(), std::fstream::out | std::fstream::trunc );
file << "#i\t t\t tSim\t tLoadBal\t numProcs\t blocks (min/max/sum)\n";
file.close();
}
// process local timing measurements and predicted loads
std::string processLocalFiles(baseFolder + "/processLocalFiles");
WALBERLA_ROOT_SECTION()
{
filesystem::path tpath( processLocalFiles );
if( !filesystem::exists( tpath ) )
filesystem::create_directory( tpath );
}
std::string measurementFileProcessName(processLocalFiles + "/measurements_" + std::to_string(MPIManager::instance()->rank()) + ".txt");
{
std::ofstream file;
file.open( measurementFileProcessName.c_str(), std::fstream::out | std::fstream::trunc );
file << "#i\t t\t mTotSim\t mLB\t mLBM\t mBH\t mCoup1\t mCoup2\t mRPD\t cLBM\t cRB\t numBlocks\n";
file.close();
}
std::string predictionFileProcessName(processLocalFiles + "/predictions_" + std::to_string(MPIManager::instance()->rank()) + ".txt");
{
std::ofstream file;
file.open( predictionFileProcessName.c_str(), std::fstream::out | std::fstream::trunc );
file << "#i\t t\t wlLBM\t wlBH\t wlCoup1\t wlCoup2\t wlRPD\t edgecut\t numBlocks\n";
file.close();
}
std::vector< std::vector<std::string> > timerKeys;
std::vector<std::string> LBMTimer;
LBMTimer.emplace_back("Stream&Collide");
timerKeys.push_back(LBMTimer);
std::vector<std::string> bhTimer;
bhTimer.emplace_back("Boundary Handling");
timerKeys.push_back(bhTimer);
std::vector<std::string> couplingTimer1;
couplingTimer1.emplace_back("Particle Mapping");
std::vector<std::string> couplingTimer2;
couplingTimer2.emplace_back("PDF Restore");
timerKeys.push_back(couplingTimer1);
timerKeys.push_back(couplingTimer2);
std::vector<std::string> rpdTimer;
rpdTimer.emplace_back("RPD Force");
rpdTimer.emplace_back("RPD VV1");
rpdTimer.emplace_back("RPD VV2");
rpdTimer.emplace_back("RPD Lub");
rpdTimer.emplace_back("RPD Collision");
timerKeys.push_back(rpdTimer);
std::vector<std::string> LBMCommTimer;
LBMCommTimer.emplace_back("LBM Communication");
LBMCommTimer.emplace_back("PDF Communication");
timerKeys.push_back(LBMCommTimer);
std::vector<std::string> rpdCommTimer;
rpdCommTimer.emplace_back("Reduce Force Torque");
rpdCommTimer.emplace_back("Reduce Hyd Force Torque");
rpdCommTimer.emplace_back("Reduce Contact History");
rpdCommTimer.emplace_back("Sync");
timerKeys.push_back(rpdCommTimer);
std::vector<real_t> timings(timerKeys.size());
double oldmTotSim = 0.0;
double oldmLB = 0.0;
auto measurementFileCounter = uint_t(0);
auto predictionFileCounter = uint_t(0);
std::string loadEvaluationStep("load evaluation");
std::string loadBalancingStep("load balancing");
std::string simulationStep("simulation");
WcTimingPool timeloopTiming;
WcTimingPool simulationTiming;
lbm_mesapd_coupling::InspectionProbe<PdfField_T,BoundaryHandling_T,ParticleAccessor_T> probeForNaNs( Vector3<real_t>(0, 0, 0),
blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor,
true, true, "" );
// time loop
for (uint_t i = 0; i < timesteps; ++i )
{
// evaluate measurements (note: reflect simulation behavior BEFORE the evaluation)
if( loadEvaluationFrequency > 0 && i % loadEvaluationFrequency == 0 && i > 0 && fileIO)
{
simulationTiming[loadEvaluationStep].start();
// write process local timing measurements to files (per process, per load balancing step)
{
auto mTotSim = simulationTiming[simulationStep].total();
auto mLB = simulationTiming[loadBalancingStep].total();
evaluateTimers(timeloopTiming, timerKeys, timings);
auto & forest = blocks->getBlockForest();
uint_t numBlocks = forest.getNumberOfBlocks();
// write to process local file
std::ofstream file;
file.open( measurementFileProcessName.c_str(), std::ofstream::app );
file << measurementFileCounter << "\t " << real_c(i) / tRef << "\t"
<< mTotSim - oldmTotSim << "\t" << mLB - oldmLB << "\t";
for (real_t timing : timings) {
file << "\t" << timing;
}
file << "\t" << numBlocks << "\n";
file.close();
oldmTotSim = mTotSim;
oldmLB = mLB;
measurementFileCounter++;
// reset timer to have measurement from evaluation to evaluation point
timeloopTiming.clear();
}
// evaluate general simulation infos (on root)
{
double totalTimeToCurrentTimestep(0.0);
double totalLBTimeToCurrentTimestep(0.0);
evaluateTotalSimulationTimePassed(simulationTiming, simulationStep, loadBalancingStep,
totalTimeToCurrentTimestep, totalLBTimeToCurrentTimestep);
math::DistributedSample numberOfBlocks;
auto & forest = blocks->getBlockForest();
uint_t numBlocks = forest.getNumberOfBlocks();
numberOfBlocks.castToRealAndInsert(numBlocks);
numberOfBlocks.mpiGatherRoot();
WALBERLA_ROOT_SECTION()
{
std::ofstream file;
file.open( infoFileName.c_str(), std::ofstream::app );
file << i << "\t " << real_c(i) / tRef << "\t"
<< totalTimeToCurrentTimestep << "\t " << totalLBTimeToCurrentTimestep << "\t "
<< numberOfProcesses << "\t ";
file << uint_c(numberOfBlocks.min()) << "\t ";
file << uint_c(numberOfBlocks.max()) << "\t ";
file << uint_c(numberOfBlocks.sum()) << "\n ";
file.close();
}
}
simulationTiming[loadEvaluationStep].end();
}
if( loadBalancingCheckFrequency != 0 && i % loadBalancingCheckFrequency == 0)
{
if(useProgressLogging) walberla::logging::Logging::instance()->setFileLogLevel(logging::Logging::LogLevel::PROGRESS);
simulationTiming[loadBalancingStep].start();
if( loadEvaluationStrategy != "LBM" ) {
WALBERLA_LOG_INFO_ON_ROOT("Checking for load balancing...");
// update info collections for the particle presence based check and the load balancing:
auto &forest = blocks->getBlockForest();
lbm_mesapd_coupling::amr::updateAndSyncInfoCollection<BoundaryHandling_T,ParticleAccessor_T >(forest, boundaryHandlingID, *accessor, numRPDSubCycles, *couplingInfoCollection);
if(useProgressLogging) WALBERLA_LOG_INFO("Created info collection with size " << couplingInfoCollection->size());
auto numBlocksBefore = forest.getNumberOfBlocks();
if(useProgressLogging) WALBERLA_LOG_INFO("Number of blocks before refresh " << numBlocksBefore);
mesa_pd::mpi::ClearNextNeighborSync CNNS;
CNNS(*accessor);
if(useProgressLogging) WALBERLA_LOG_INFO_ON_ROOT("Cleared particle sync and starting refresh");
blocks->refresh();
auto numBlocksAfter = forest.getNumberOfBlocks();
if(useProgressLogging) WALBERLA_LOG_INFO("Number of blocks after refresh " << numBlocksAfter);
if(useProgressLogging) WALBERLA_LOG_INFO_ON_ROOT("Refresh finished, recreating all datastructures");
//WALBERLA_LOG_INFO(rank << ", " << numBlocksBefore << " -> " << numBlocksAfter);
rpdDomain->refresh();
if(useBlockForestSync)
{
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, associateToBlock, *accessor);
syncNextNeighborBlockForestFunc(*ps, blocks->getBlockForestPointer(), rpdDomain, overlap);
} else
{
syncNextNeighborFunc(*ps, *rpdDomain, overlap);
}
clearBoundaryHandling(forest, boundaryHandlingID);
clearParticleField(forest, particleFieldID, *accessor);
recreateBoundaryHandling(forest, boundaryHandlingID);
//NOTE: order in mapping is important: first all global/fixed particles,
// only then moving ones, which do not overwrite the mapping of the global/fixed ones
particleMappingCall();
}
simulationTiming[loadBalancingStep].end();
if(useProgressLogging) walberla::logging::Logging::instance()->setFileLogLevel(logging::Logging::LogLevel::INFO);
}
if(checkSimulation)
{
WALBERLA_LOG_INFO_ON_ROOT("Checking time step " << i );
evaluateParticleSimulation(ps, ss, numberOfSediments, numGlobalParticles);
checkMapping(blocks, pdfFieldID, boundaryHandlingID, probeForNaNs);
}
// evaluate predictions (note: reflect the predictions for all upcoming simulations, thus the corresponding measurements have to be taken afterwards)
if( loadEvaluationFrequency > 0 && i % loadEvaluationFrequency == 0 && fileIO)
{
simulationTiming[loadEvaluationStep].start();
// write process local load predictions to files (per process, per load balancing step)
{
auto wlLBM = real_t(0);
auto wlBH = real_t(0);
auto wlCoup1 = real_t(0);
auto wlCoup2 = real_t(0);
auto wlRPD = real_t(0);
auto & forest = blocks->getBlockForest();
lbm_mesapd_coupling::amr::updateAndSyncInfoCollection<BoundaryHandling_T,ParticleAccessor_T >(forest, boundaryHandlingID, *accessor, numRPDSubCycles, *couplingInfoCollection);
for( auto blockIt = forest.begin(); blockIt != forest.end(); ++blockIt ) {
auto * block = static_cast<blockforest::Block *> (&(*blockIt));
const auto &blockID = block->getId();
auto infoIt = couplingInfoCollection->find(blockID);
auto blockInfo = infoIt->second;
wlLBM += lbm_mesapd_coupling::amr::fittedLBMWeightEvaluationFunction(blockInfo);
wlBH += lbm_mesapd_coupling::amr::fittedBHWeightEvaluationFunction(blockInfo);
wlCoup1 += lbm_mesapd_coupling::amr::fittedCoup1WeightEvaluationFunction(blockInfo);
wlCoup2 += lbm_mesapd_coupling::amr::fittedCoup2WeightEvaluationFunction(blockInfo);
wlRPD += lbm_mesapd_coupling::amr::fittedRPDWeightEvaluationFunction(blockInfo);
}
// note: we count the edge weight doubled here in total (to and from the other process). ParMetis only counts one direction.
uint_t edgecut = evaluateEdgeCut(forest);
uint_t numBlocks = forest.getNumberOfBlocks();
std::ofstream file;
file.open( predictionFileProcessName.c_str(), std::ofstream::app );
file << predictionFileCounter << "\t " << real_c(i) / tRef << "\t"
<< wlLBM << "\t" << wlBH << "\t" << wlCoup1 << "\t" << wlCoup2 << "\t" << wlRPD << "\t"
<< edgecut << "\t" << numBlocks << "\n";
file.close();
predictionFileCounter++;;
}
simulationTiming[loadEvaluationStep].end();
}
simulationTiming[simulationStep].start();
// perform a single simulation step
timeloop.singleStep( timeloopTiming );
timeloopTiming["Reduce Hyd Force Torque"].start();
reduceProperty.operator()<mesa_pd::HydrodynamicForceTorqueNotification>(*ps);
timeloopTiming["Reduce Hyd Force Torque"].end();
timeloopTiming["RPD Force"].start();
if( i == 0 )
{
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
}
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, averageHydrodynamicForceTorque, *accessor );
timeloopTiming["RPD Force"].end();
if(checkSimulation)
{
checkParticleProperties(ps);
}
for(auto subCycle = uint_t(0); subCycle < numRPDSubCycles; ++subCycle )
{
timeloopTiming["RPD VV1"].start();
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, vvIntegratorPreForce, *accessor);
timeloopTiming["RPD VV1"].end();
timeloopTiming["Sync"].start();
if(useBlockForestSync)
{
syncNextNeighborBlockForestFunc(*ps, blocks->getBlockForestPointer(), rpdDomain, overlap);
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, associateToBlock, *accessor);
} else
{
syncNextNeighborFunc(*ps, *rpdDomain, overlap);
}
timeloopTiming["Sync"].end();
// lubrication correction
timeloopTiming["RPD Lub"].start();
ps->forEachParticlePairHalf(useOpenMP, mesa_pd::kernel::ExcludeInfiniteInfinite(), *accessor,
[&lubricationCorrectionKernel,maximumLubricationCutOffDistance, &rpdDomain]
(const size_t idx1, const size_t idx2, auto& ac)
{
mesa_pd::collision_detection::AnalyticContactDetection acd;
acd.getContactThreshold() = maximumLubricationCutOffDistance;
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 );
timeloopTiming["RPD Lub"].end();
// collision response
timeloopTiming["RPD Collision"].start();
ps->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 );
timeloopTiming["RPD Collision"].end();
timeloopTiming["Reduce Contact History"].start();
reduceAndSwapContactHistory(*ps);
timeloopTiming["Reduce Contact History"].end();
timeloopTiming["RPD Force"].start();
lbm_mesapd_coupling::AddHydrodynamicInteractionKernel addHydrodynamicInteraction;
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addHydrodynamicInteraction, *accessor );
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addGravitationalForce, *accessor );
timeloopTiming["RPD Force"].end();
timeloopTiming["Reduce Force Torque"].start();
reduceProperty.operator()<mesa_pd::ForceTorqueNotification>(*ps);
timeloopTiming["Reduce Force Torque"].end();
// integration
timeloopTiming["RPD VV2"].start();
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, vvIntegratorPostForce, *accessor);
timeloopTiming["RPD VV2"].end();
timeloopTiming["Sync"].start();
if(useBlockForestSync)
{
syncNextNeighborBlockForestFunc(*ps, blocks->getBlockForestPointer(), rpdDomain, overlap);
} else
{
syncNextNeighborFunc(*ps, *rpdDomain, overlap);
}
timeloopTiming["Sync"].end();
}
timeloopTiming["RPD Force"].start();
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, resetHydrodynamicForceTorque, *accessor );
timeloopTiming["RPD Force"].end();
// update particle mapping
timeloopAfterParticles.singleStep(timeloopTiming);
simulationTiming[simulationStep].end();
}
simulationTiming.logResultOnRoot();
return EXIT_SUCCESS;
}
} // namespace fluid_particle_workload_distribution
int main( int argc, char **argv ){
fluid_particle_workload_distribution::main(argc, argv);
}
apps/benchmarks/FluidParticleCouplingWithLoadBalancing/FluidParticleWorkloadEvaluation.cpp 0 → 100644
View file @ 2dee476c
//======================================================================================================================
//
// 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 FluidParticleWorkLoadEvaluation.cpp
//! \author Christoph Rettinger <christoph.rettinger@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/debug/Debug.h"
#include "core/debug/TestSubsystem.h"
#include "core/logging/Logging.h"
#include "core/math/all.h"
#include "core/timing/RemainingTimeLogger.h"
#include "core/mpi/MPIManager.h"
#include "core/mpi/Reduce.h"
#include "core/mpi/Broadcast.h"
#include "domain_decomposition/SharedSweep.h"
#include "field/AddToStorage.h"
#include "field/communication/PackInfo.h"
#include "lbm/communication/PdfFieldPackInfo.h"
#include "lbm/field/AddToStorage.h"
#include "lbm/field/PdfField.h"
#include "lbm/lattice_model/D3Q19.h"
#include "lbm/lattice_model/ForceModel.h"
#include "lbm/sweeps/CellwiseSweep.h"
#include "lbm/sweeps/SweepWrappers.h"
#include "lbm/BlockForestEvaluation.h"
#include "lbm_mesapd_coupling/mapping/ParticleMapping.h"
#include "lbm_mesapd_coupling/momentum_exchange_method/MovingParticleMapping.h"
#include "lbm_mesapd_coupling/momentum_exchange_method/boundary/CurvedLinear.h"
#include "lbm_mesapd_coupling/momentum_exchange_method/reconstruction/Reconstructor.h"
#include "lbm_mesapd_coupling/momentum_exchange_method/reconstruction/PdfReconstructionManager.h"
#include "lbm_mesapd_coupling/utility/AddForceOnParticlesKernel.h"
#include "lbm_mesapd_coupling/utility/ParticleSelector.h"
#include "lbm_mesapd_coupling/DataTypes.h"
#include "lbm_mesapd_coupling/utility/AverageHydrodynamicForceTorqueKernel.h"
#include "lbm_mesapd_coupling/utility/AddHydrodynamicInteractionKernel.h"
#include "lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
#include "lbm_mesapd_coupling/utility/LubricationCorrectionKernel.h"
#include "lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h"
#include "mesa_pd/collision_detection/AnalyticContactDetection.h"
#include "mesa_pd/data/ParticleAccessorWithShape.h"
#include "mesa_pd/data/ParticleStorage.h"
#include "mesa_pd/data/ShapeStorage.h"
#include "mesa_pd/data/DataTypes.h"
#include "mesa_pd/data/shape/HalfSpace.h"
#include "mesa_pd/data/shape/Sphere.h"
#include "mesa_pd/domain/BlockForestDomain.h"
#include "mesa_pd/domain/BlockForestDataHandling.h"
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/LinearSpringDashpot.h"
#include "mesa_pd/kernel/ParticleSelector.h"
#include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h"
#include "mesa_pd/mpi/ReduceContactHistory.h"
#include "mesa_pd/mpi/ContactFilter.h"
#include "mesa_pd/mpi/notifications/ForceTorqueNotification.h"
#include "mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h"
#include "mesa_pd/vtk/ParticleVtkOutput.h"
#include "timeloop/SweepTimeloop.h"
#include "vtk/all.h"
#include "field/vtk/all.h"
#include "lbm/vtk/all.h"
#include <vector>
#include <iostream>
namespace fluid_particle_workload_evaluation
{
///////////
// USING //
///////////
using namespace walberla;
using walberla::uint_t;
using LatticeModel_T = lbm::D3Q19< lbm::collision_model::TRT, false >;
using Stencil_T = LatticeModel_T::Stencil;
using PdfField_T = lbm::PdfField<LatticeModel_T>;
using flag_t = walberla::uint8_t;
using FlagField_T = FlagField<flag_t>;
const uint_t FieldGhostLayers = 1;
///////////
// FLAGS //
///////////
const FlagUID Fluid_Flag ( "fluid" );
const FlagUID MO_Flag ( "moving obstacle" );
const FlagUID FormerMO_Flag( "former moving obstacle" );
/////////////////////////////////////
// BOUNDARY HANDLING CUSTOMIZATION //
/////////////////////////////////////
template <typename ParticleAccessor_T>
class MyBoundaryHandling
{
public:
using MO_T = lbm_mesapd_coupling::CurvedLinear< LatticeModel_T, FlagField_T, ParticleAccessor_T >;
using Type = BoundaryHandling< FlagField_T, Stencil_T, MO_T >;
MyBoundaryHandling( const BlockDataID & flagFieldID, const BlockDataID & pdfFieldID,
const BlockDataID & particleFieldID, const shared_ptr<ParticleAccessor_T>& ac,
bool useEntireFieldTraversal) :
flagFieldID_( flagFieldID ), pdfFieldID_( pdfFieldID ), particleFieldID_( particleFieldID ), ac_( ac ), useEntireFieldTraversal_(useEntireFieldTraversal) {}
Type * operator()( IBlock* const block, const StructuredBlockStorage* const storage ) const
{
WALBERLA_ASSERT_NOT_NULLPTR( block );
WALBERLA_ASSERT_NOT_NULLPTR( storage );
auto * flagField = block->getData< FlagField_T >( flagFieldID_ );
auto * pdfField = block->getData< PdfField_T > ( pdfFieldID_ );
auto * particleField = block->getData< lbm_mesapd_coupling::ParticleField_T > ( particleFieldID_ );
const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
typename Type::Mode mode = (useEntireFieldTraversal_) ? Type::Mode::ENTIRE_FIELD_TRAVERSAL : Type::Mode::OPTIMIZED_SPARSE_TRAVERSAL ;
Type * handling = new Type( "moving obstacle boundary handling", flagField, fluid,
MO_T( "MO", MO_Flag, pdfField, flagField, particleField, ac_, fluid, *storage, *block ),
mode);
handling->fillWithDomain( FieldGhostLayers );
return handling;
}
private:
const BlockDataID flagFieldID_;
const BlockDataID pdfFieldID_;
const BlockDataID particleFieldID_;
shared_ptr<ParticleAccessor_T> ac_;
bool useEntireFieldTraversal_;
};
template< typename BoundaryHandling_T >
void evaluateFluidQuantities(const shared_ptr< StructuredBlockStorage > & blocks, const BlockDataID boundaryHandlingID,
uint_t & numCells, uint_t & numFluidCells, uint_t & numNBCells )
{
for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
{
auto * boundaryHandling = blockIt->getData< BoundaryHandling_T >( boundaryHandlingID );
auto xyzSize = boundaryHandling->getFlagField()->xyzSize();
numCells += xyzSize.numCells();
for(auto z = cell_idx_t(xyzSize.zMin()); z <= cell_idx_t(xyzSize.zMax()); ++z ){
for(auto y = cell_idx_t(xyzSize.yMin()); y <= cell_idx_t(xyzSize.yMax()); ++y ){
for(auto x = cell_idx_t(xyzSize.xMin()); x <= cell_idx_t(xyzSize.xMax()); ++x ) {
if (boundaryHandling->isDomain(x, y, z)) {
++numFluidCells;
}
if (boundaryHandling->isNearBoundary(x, y, z)) {
++numNBCells;
}
}
}
}
}
}
void evaluateRPDQuantities( const shared_ptr< mesa_pd::data::ParticleStorage > & ps,
uint_t & numLocalParticles, uint_t & numGhostParticles)
{
for (auto pIt = ps->begin(); pIt != ps->end(); ++pIt)
{
using namespace walberla::mesa_pd::data::particle_flags;
if (isSet(pIt->getFlags(), GHOST))
{
++numGhostParticles;
} else
{
//note: global particles are included here
// use if (!isSet(pIt->getFlags(), GLOBAL)) if should be excluded
++numLocalParticles;
}
}
}
void evaluateTimers(WcTimingPool & timingPool,
const std::vector<std::vector<std::string> > & timerKeys,
std::vector<double> & timings )
{
for (auto & timingsIt : timings)
{
timingsIt = 0.0;
}
timingPool.unifyRegisteredTimersAcrossProcesses();
double scalingFactor = 1000.0; // milliseconds
for (auto i = uint_t(0); i < timerKeys.size(); ++i )
{
auto keys = timerKeys[i];
for (const auto &timerName : keys)
{
if(timingPool.timerExists(timerName))
{
timings[i] += timingPool[timerName].total() * scalingFactor;
}
}
}
}
//*******************************************************************************************************************
/*! Application to evaluate the workload (time measurements) for a fluid-particle simulation
*
* This application is used in the paper
* Rettinger, Ruede - "Dynamic Load Balancing Techniques for Particulate Flow Simulations", 2019, Computation
* in Section 3 to develop and calibrate the workload estimator.
* The setup features settling particle inside a horizontally periodic box.
* A comprehensive description is given in Sec. 3.3 of the paper.
* It uses 4 x 4 x 5 blocks for domain partitioning.
* For each block ( = each process), the block local quantities are evaluated as well as the timing infos of
* the fluid-particle interaction algorithm. Those infos are then written to files that can be used later on
* for function fitting to obtain a workload estimator.
*
* NOTE: Since this estimator relies on timing measurements, this evaluation procedure should be carried out everytime
* a different implementation, hardware or algorithm is used.
*
*/
//*******************************************************************************************************************
int main( int argc, char **argv )
{
debug::enterTestMode();
mpi::Environment env( argc, argv );
auto solidVolumeFraction = real_t(0.2);
// LBM / numerical parameters
auto blockSize = uint_t(32);
auto uSettling = real_t(0.1); // characteristic settling velocity
auto diameter = real_t(10);
auto Ga = real_t(30); //Galileo number
auto numRPDSubCycles = uint_t(10);
auto vtkIOFreq = uint_t(0);
auto timestepsNonDim = real_t(2.5);
auto numSamples = uint_t(2000);
std::string baseFolder = "workload_files"; // folder for vtk and file output
bool noFileOutput = false;
bool useEntireFieldTraversal = true;
bool useFusedStreamCollide = false;
auto XBlocks = uint_t(4);
auto YBlocks = uint_t(4);
auto ZBlocks = uint_t(5);
real_t topWallOffsetFactor = real_t(1.05);
bool vtkDuringInit = false;
for( int i = 1; i < argc; ++i )
{
if( std::strcmp( argv[i], "--vtkIOFreq" ) == 0 ) { vtkIOFreq = uint_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--noFileOutput" ) == 0 ) { noFileOutput = true; continue; }
if( std::strcmp( argv[i], "--vtkDuringInit" ) == 0 ) { vtkDuringInit = true; continue; }
if( std::strcmp( argv[i], "--basefolder" ) == 0 ) { baseFolder = argv[++i]; continue; }
if( std::strcmp( argv[i], "--solidVolumeFraction" ) == 0 ) { solidVolumeFraction = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--diameter" ) == 0 ) { diameter = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--blockSize" ) == 0 ) { blockSize = uint_c(std::atof( argv[++i]) ); continue; }
if( std::strcmp( argv[i], "--XBlocks" ) == 0 ) { XBlocks = uint_c(std::atof( argv[++i]) ); continue; }
if( std::strcmp( argv[i], "--YBlocks" ) == 0 ) { YBlocks = uint_c(std::atof( argv[++i]) ); continue; }
if( std::strcmp( argv[i], "--ZBlocks" ) == 0 ) { ZBlocks = uint_c(std::atof( argv[++i]) ); continue; }
if( std::strcmp( argv[i], "--uSettling" ) == 0 ) { uSettling = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--topWallOffsetFactor" ) == 0 ) { topWallOffsetFactor = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--Ga" ) == 0 ) { Ga = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--timestepsNonDim" ) == 0 ) { timestepsNonDim = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--numRPDSubCycles" ) == 0 ) { numRPDSubCycles = uint_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--useEntireFieldTraversal" ) == 0 ) { useEntireFieldTraversal = true; continue; }
if( std::strcmp( argv[i], "--numSamples" ) == 0 ) { numSamples = uint_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--useFusedStreamCollide" ) == 0 ) { useFusedStreamCollide = true; continue; }
WALBERLA_ABORT("Unrecognized command line argument found: " << argv[i]);
}
WALBERLA_CHECK(diameter > real_t(1));
WALBERLA_CHECK(uSettling > real_t(0));
WALBERLA_CHECK(Ga > real_t(0));
WALBERLA_CHECK(solidVolumeFraction > real_t(0));
WALBERLA_CHECK(solidVolumeFraction < real_t(0.65));
///////////////////////////
// SIMULATION PROPERTIES //
///////////////////////////
if( MPIManager::instance()->numProcesses() != int(XBlocks * YBlocks * ZBlocks) )
{
WALBERLA_LOG_WARNING_ON_ROOT("WARNING! You have specified less or more processes than number of blocks -> the time measurements are no longer blockwise!")
}
if( diameter > real_c(blockSize) )
{
WALBERLA_LOG_WARNING_ON_ROOT("The bodies might be too large to work with the currently used synchronization!");
}
WALBERLA_LOG_INFO_ON_ROOT("Using setup with sedimenting particles -> creating two planes and applying gravitational force")
const uint_t XCells = blockSize * XBlocks;
const uint_t YCells = blockSize * YBlocks;
const uint_t ZCells = blockSize * ZBlocks;
const real_t topWallOffset = topWallOffsetFactor * real_t(blockSize); // move the top wall downwards to take away a certain portion of the overall domain
// determine number of spheres to generate, if necessary scale diameter a bit to reach desired solid volume fraction
real_t domainHeight = real_c(ZCells) - topWallOffset;
real_t fluidVolume = real_c( XCells * YCells ) * domainHeight;
real_t solidVolume = solidVolumeFraction * fluidVolume;
uint_t numberOfParticles = uint_c(std::ceil(solidVolume / ( math::pi / real_t(6) * diameter * diameter * diameter )));
diameter = std::cbrt( solidVolume / ( real_c(numberOfParticles) * math::pi / real_t(6) ) );
auto densityRatio = real_t(2.5);
real_t viscosity = uSettling * diameter / Ga;
const real_t omega = lbm::collision_model::omegaFromViscosity(viscosity);
const real_t gravitationalAcceleration = uSettling * uSettling / ( (densityRatio-real_t(1)) * diameter );
real_t tref = diameter / uSettling;
real_t Tref = domainHeight / uSettling;
uint_t timesteps = uint_c(timestepsNonDim * Tref);
const real_t dx = real_c(1);
WALBERLA_LOG_INFO_ON_ROOT("viscosity = " << viscosity);
WALBERLA_LOG_INFO_ON_ROOT("tau = " << real_t(1)/omega);
WALBERLA_LOG_INFO_ON_ROOT("diameter = " << diameter);
WALBERLA_LOG_INFO_ON_ROOT("solid volume fraction = " << solidVolumeFraction);
WALBERLA_LOG_INFO_ON_ROOT("domain size (in cells) = " << XCells << " x " << YCells << " x " << ZCells);
WALBERLA_LOG_INFO_ON_ROOT("number of bodies = " << numberOfParticles);
WALBERLA_LOG_INFO_ON_ROOT("gravitational acceleration = " << gravitationalAcceleration);
WALBERLA_LOG_INFO_ON_ROOT("Ga = " << Ga);
WALBERLA_LOG_INFO_ON_ROOT("uSettling = " << uSettling);
WALBERLA_LOG_INFO_ON_ROOT("tref = " << tref);
WALBERLA_LOG_INFO_ON_ROOT("Tref = " << Tref);
WALBERLA_LOG_INFO_ON_ROOT("timesteps = " << timesteps);
WALBERLA_LOG_INFO_ON_ROOT("number of workload samples = " << numSamples);
// create folder to store logging files
WALBERLA_ROOT_SECTION()
{
walberla::filesystem::path path1( baseFolder );
if( !walberla::filesystem::exists( path1 ) )
walberla::filesystem::create_directory( path1 );
}
///////////////////////////
// BLOCK STRUCTURE SETUP //
///////////////////////////
Vector3<bool> periodicity( true );
periodicity[2] = false;
// create domain
shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid( XBlocks, YBlocks, ZBlocks, blockSize, blockSize, blockSize, dx,
0, false, false, //one block per process!
periodicity[0], periodicity[1], periodicity[2], //periodicity
false );
/////////
// RPD //
/////////
const real_t restitutionCoeff = real_t(0.97);
const real_t frictionCoeffStatic = real_t(0.8);
const real_t frictionCoeffDynamic = real_t(0.15);
const real_t collisionTime = real_t(4) * diameter; // from my paper
const real_t poissonsRatio = real_t(0.22);
const real_t kappa = real_t(2) * ( real_t(1) - poissonsRatio ) / ( real_t(2) - poissonsRatio ) ;
auto rpdDomain = std::make_shared<mesa_pd::domain::BlockForestDomain>(blocks->getBlockForestPointer());
//init data structures
auto ps = walberla::make_shared<mesa_pd::data::ParticleStorage>(1);
auto ss = walberla::make_shared<mesa_pd::data::ShapeStorage>();
using ParticleAccessor_T = mesa_pd::data::ParticleAccessorWithShape;
auto accessor = walberla::make_shared<ParticleAccessor_T >(ps, ss);
real_t timeStepSizeRPD = real_t(1)/real_t(numRPDSubCycles);
mesa_pd::kernel::VelocityVerletPreForceUpdate vvIntegratorPreForce(timeStepSizeRPD);
mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(timeStepSizeRPD);
// types: 0 = wall, 1: sphere
mesa_pd::kernel::LinearSpringDashpot collisionResponse(2);
collisionResponse.setFrictionCoefficientDynamic(0,1,frictionCoeffDynamic);
collisionResponse.setFrictionCoefficientDynamic(1,1,frictionCoeffDynamic);
collisionResponse.setFrictionCoefficientStatic(0,1,frictionCoeffStatic);
collisionResponse.setFrictionCoefficientStatic(1,1,frictionCoeffStatic);
const real_t sphereVolume = math::pi / real_t(6) * diameter * diameter * diameter;
const real_t particleMass = densityRatio * sphereVolume;
const real_t effMass_sphereWall = particleMass;
const real_t effMass_sphereSphere = particleMass * particleMass / ( real_t(2) * particleMass );
collisionResponse.setStiffnessAndDamping(0,1,restitutionCoeff,collisionTime,kappa,effMass_sphereWall);
collisionResponse.setStiffnessAndDamping(1,1,restitutionCoeff,collisionTime,kappa,effMass_sphereSphere);
mesa_pd::mpi::ReduceProperty reduceProperty;
mesa_pd::mpi::ReduceContactHistory reduceAndSwapContactHistory;
//////////////
// COUPLING //
//////////////
// connect to pe
const real_t overlap = real_c( 1.5 ) * dx;
std::function<void(void)> syncCall = [&ps,&rpdDomain,overlap](){
mesa_pd::mpi::SyncNextNeighbors syncNextNeighborFunc;
syncNextNeighborFunc(*ps, *rpdDomain, overlap);
};
auto generationDomain = AABB( real_t(0), real_t(0), real_t(0), real_c(XCells), real_c(YCells), real_c(ZCells) - topWallOffset);
// create plane at top and bottom
mesa_pd::data::Particle&& p0 = *ps->create(true);
p0.setPosition(generationDomain.minCorner());
p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,1) ));
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);
mesa_pd::data::Particle&& p1 = *ps->create(true);
p1.setPosition(generationDomain.maxCorner());
p1.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,-1) ));
p1.setOwner(mpi::MPIManager::instance()->rank());
p1.setType(0);
mesa_pd::data::particle_flags::set(p1.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
mesa_pd::data::particle_flags::set(p1.getFlagsRef(), mesa_pd::data::particle_flags::FIXED);
auto sphereShape = ss->create<mesa_pd::data::Sphere>( diameter * real_t(0.5) );
ss->shapes[sphereShape]->updateMassAndInertia(densityRatio);
auto xParticle = real_t(0);
auto yParticle = real_t(0);
auto zParticle = real_t(0);
auto rank = mpi::MPIManager::instance()->rank();
for( uint_t nPart = 0; nPart < numberOfParticles; ++nPart )
{
WALBERLA_ROOT_SECTION()
{
xParticle = math::realRandom<real_t>(generationDomain.xMin(), generationDomain.xMax());
yParticle = math::realRandom<real_t>(generationDomain.yMin(), generationDomain.yMax());
zParticle = math::realRandom<real_t>(generationDomain.zMin(), generationDomain.zMax());
}
WALBERLA_MPI_SECTION()
{
mpi::broadcastObject( xParticle );
mpi::broadcastObject( yParticle );
mpi::broadcastObject( zParticle );
}
auto position = Vector3<real_t>( xParticle, yParticle, zParticle );
//WALBERLA_LOG_INFO_ON_ROOT(position);
if (!rpdDomain->isContainedInProcessSubdomain(uint_c(rank), position)) continue;
auto p = ps->create();
p->setPosition(position);
p->setInteractionRadius(diameter * real_t(0.5));
p->setShapeID(sphereShape);
p->setType(1);
p->setOwner(rank);
}
syncCall();
// resolve possible overlaps of the particles due to the random initialization via a particle-only simulation
const bool useOpenMP = false;
const real_t dt_RPD_Init = real_t(1);
const auto initialParticleSimSteps = uint_t(20000);
const real_t overlapLimit = real_t(0.001) * diameter;
auto particleVtkOutput = make_shared<mesa_pd::vtk::ParticleVtkOutput>(ps);
particleVtkOutput->addOutput<mesa_pd::data::SelectParticleLinearVelocity>("velocity");
particleVtkOutput->setParticleSelector( [sphereShape](const mesa_pd::data::ParticleStorage::iterator& pIt) {return !mesa_pd::data::particle_flags::isSet(pIt->getFlags(), mesa_pd::data::particle_flags::GHOST) && pIt->getShapeID() == sphereShape;} ); //limit output to local sphere
auto particleVtkWriterInit = vtk::createVTKOutput_PointData(particleVtkOutput, "Particles_init", 1, baseFolder, "simulation_step");
for(auto pet = uint_t(0); pet <= initialParticleSimSteps; ++pet )
{
real_t maxPenetrationDepth = 0;
ps->forEachParticlePairHalf(useOpenMP, mesa_pd::kernel::ExcludeInfiniteInfinite(), *accessor,
[&collisionResponse, &rpdDomain, &maxPenetrationDepth, dt_RPD_Init]
(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))
{
maxPenetrationDepth = std::max(maxPenetrationDepth, std::abs(acd.getPenetrationDepth()));
collisionResponse(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(),
acd.getContactNormal(), acd.getPenetrationDepth(), dt_RPD_Init);
}
}
},
*accessor );
reduceAndSwapContactHistory(*ps);
mpi::allReduceInplace(maxPenetrationDepth, mpi::MAX);
reduceProperty.operator()<mesa_pd::ForceTorqueNotification>(*ps);
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, mesa_pd::kernel::ExplicitEuler(dt_RPD_Init), *accessor);
syncCall();
if( pet % uint_t(20) == uint_t(0) )
{
if(vtkDuringInit)
{
particleVtkWriterInit->write();
}
WALBERLA_LOG_INFO_ON_ROOT(pet << " - current max overlap = " << maxPenetrationDepth / diameter * real_t(100) << "%");
}
if(maxPenetrationDepth < overlapLimit) break;
// reset velocities to avoid too large ones
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor,
[](const size_t idx, ParticleAccessor_T& ac){
ac.setLinearVelocity(idx, ac.getLinearVelocity(idx) * real_t(0.5));
ac.setAngularVelocity(idx, ac.getAngularVelocity(idx) * real_t(0.5));
}, *accessor);
}
// reset all velocities to zero
Vector3<real_t> initialSphereVelocity(real_t(0));
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, [](const size_t idx, ParticleAccessor_T& ac){
ac.getNewContactHistoryRef(idx).clear();
ac.getOldContactHistoryRef(idx).clear();
}, *accessor);
WALBERLA_LOG_INFO_ON_ROOT("Setting initial velocity " << initialSphereVelocity << " of all spheres");
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor,
[initialSphereVelocity](const size_t idx, ParticleAccessor_T& ac){
ac.setLinearVelocity(idx, initialSphereVelocity);
ac.setAngularVelocity(idx, Vector3<real_t>(real_t(0)));
}, *accessor);
syncCall();
///////////////////////
// ADD DATA TO BLOCKS //
////////////////////////
// create the lattice model
LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( omega ) );
// add PDF field
BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (fzyx)", latticeModel,
Vector3< real_t >( real_t(0) ), real_t(1),
uint_t(1), field::fzyx );
// add flag field
BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>( blocks, "flag field" );
// add particle field
BlockDataID particleFieldID = field::addToStorage<lbm_mesapd_coupling::ParticleField_T>( blocks, "particle field", accessor->getInvalidUid(), field::fzyx, FieldGhostLayers );
// add boundary handling & initialize outer domain boundaries
using BoundaryHandling_T = MyBoundaryHandling<ParticleAccessor_T>::Type;
BlockDataID boundaryHandlingID = blocks->addStructuredBlockData< BoundaryHandling_T >(MyBoundaryHandling<ParticleAccessor_T>( flagFieldID, pdfFieldID, particleFieldID, accessor, useEntireFieldTraversal), "boundary handling" );
Vector3<real_t> gravitationalForce( real_t(0), real_t(0), -(densityRatio - real_t(1)) * gravitationalAcceleration * sphereVolume );
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;});
lbm_mesapd_coupling::ParticleMappingKernel<BoundaryHandling_T> particleMappingKernel(blocks, boundaryHandlingID);
lbm_mesapd_coupling::MovingParticleMappingKernel<BoundaryHandling_T> movingParticleMappingKernel(blocks, boundaryHandlingID, particleFieldID);
WALBERLA_LOG_INFO_ON_ROOT(" - Lubrication correction:");
WALBERLA_LOG_INFO_ON_ROOT(" - normal cut off distance = " << lubricationCorrectionKernel.getNormalCutOffDistance());
WALBERLA_LOG_INFO_ON_ROOT(" - tangential translational cut off distance = " << lubricationCorrectionKernel.getTangentialTranslationalCutOffDistance());
WALBERLA_LOG_INFO_ON_ROOT(" - tangential rotational cut off distance = " << lubricationCorrectionKernel.getTangentialRotationalCutOffDistance());
const real_t maximumLubricationCutOffDistance = std::max(lubricationCorrectionKernel.getNormalCutOffDistance(), std::max(lubricationCorrectionKernel.getTangentialRotationalCutOffDistance(), lubricationCorrectionKernel.getTangentialTranslationalCutOffDistance()));
// map planes into the LBM simulation -> act as no-slip boundaries
ps->forEachParticle(false, lbm_mesapd_coupling::GlobalParticlesSelector(), *accessor, movingParticleMappingKernel, *accessor, MO_Flag);
// map particles into the LBM simulation
ps->forEachParticle(false, lbm_mesapd_coupling::RegularParticlesSelector(), *accessor, movingParticleMappingKernel, *accessor, MO_Flag);
lbm::BlockForestEvaluation<FlagField_T> bfEval(blocks, flagFieldID, Fluid_Flag);
WALBERLA_LOG_INFO_ON_ROOT(bfEval.loggingString());
///////////////
// TIME LOOP //
///////////////
// create the timeloop
SweepTimeloop timeloop( blocks->getBlockStorage(), timesteps );
timeloop.addFuncBeforeTimeStep( RemainingTimeLogger( timeloop.getNrOfTimeSteps() ), "Remaining Time Logger" );
if( vtkIOFreq != uint_t(0) )
{
auto particleVtkWriter = vtk::createVTKOutput_PointData(particleVtkOutput, "Particles", vtkIOFreq, baseFolder, "simulation_step");
timeloop.addFuncBeforeTimeStep( vtk::writeFiles( particleVtkWriter ), "VTK (sphere data)" );
// flag field
auto flagFieldVTK = vtk::createVTKOutput_BlockData( blocks, "flag_field", vtkIOFreq, 1, false, baseFolder );
flagFieldVTK->addCellDataWriter( make_shared< field::VTKWriter< FlagField_T > >( flagFieldID, "FlagField" ) );
timeloop.addFuncBeforeTimeStep( vtk::writeFiles( flagFieldVTK ), "VTK (flag field data)" );
// pdf field
auto pdfFieldVTK = vtk::createVTKOutput_BlockData( blocks, "fluid_field", vtkIOFreq, 0, false, baseFolder );
field::FlagFieldCellFilter< FlagField_T > fluidFilter( flagFieldID );
fluidFilter.addFlag( Fluid_Flag );
pdfFieldVTK->addCellInclusionFilter( fluidFilter );
pdfFieldVTK->addCellDataWriter( make_shared< lbm::VelocityVTKWriter< LatticeModel_T, float > >( pdfFieldID, "VelocityFromPDF" ) );
pdfFieldVTK->addCellDataWriter( make_shared< lbm::DensityVTKWriter < LatticeModel_T, float > >( pdfFieldID, "DensityFromPDF" ) );
timeloop.addFuncBeforeTimeStep( vtk::writeFiles( pdfFieldVTK ), "VTK (fluid field data)" );
auto domainDecompVTK = vtk::createVTKOutput_DomainDecomposition(blocks, "domain_decomposition", vtkIOFreq, baseFolder );
timeloop.addFuncBeforeTimeStep( vtk::writeFiles(domainDecompVTK), "VTK (domain decomposition)");
}
// setup of the LBM communication for synchronizing the pdf field between neighboring blocks
blockforest::communication::UniformBufferedScheme< Stencil_T > optimizedPDFCommunicationScheme( blocks );
optimizedPDFCommunicationScheme.addPackInfo( make_shared< lbm::PdfFieldPackInfo< LatticeModel_T > >( pdfFieldID ) ); // optimized sync
auto sweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag );
if( !useFusedStreamCollide )
{
// Collide
timeloop.add() << Sweep( makeCollideSweep(sweep), "Collide" );
}
// add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
timeloop.add() << BeforeFunction( optimizedPDFCommunicationScheme, "LBM Communication" )
<< Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingID ), "Boundary Handling" );
if( useFusedStreamCollide )
{
// streaming & collide
timeloop.add() << Sweep( makeSharedSweep(sweep), "Stream&Collide" );
} else
{
// streaming
timeloop.add() << Sweep( makeStreamSweep(sweep), "Stream" );
}
SweepTimeloop timeloopAfterParticles( blocks->getBlockStorage(), timesteps );
bool conserveMomentum = false;
// sweep for updating the particle mapping into the LBM simulation
timeloopAfterParticles.add() << Sweep( lbm_mesapd_coupling::makeMovingParticleMapping<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID,boundaryHandlingID, particleFieldID, accessor, MO_Flag, FormerMO_Flag,
lbm_mesapd_coupling::RegularParticlesSelector(), conserveMomentum), "Particle Mapping" );
bool recomputeTargetDensity = false;
auto gradReconstructor = lbm_mesapd_coupling::makeGradsMomentApproximationReconstructor<BoundaryHandling_T>(blocks, boundaryHandlingID, omega, recomputeTargetDensity,true);
blockforest::communication::UniformBufferedScheme< Stencil_T > fullPDFCommunicationScheme( blocks );
fullPDFCommunicationScheme.addPackInfo( make_shared< field::communication::PackInfo< PdfField_T > >( pdfFieldID ) ); // full sync
timeloopAfterParticles.add() << BeforeFunction( fullPDFCommunicationScheme, "PDF Communication" )
<< Sweep( makeSharedSweep(lbm_mesapd_coupling::makePdfReconstructionManager<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, FormerMO_Flag, Fluid_Flag,
gradReconstructor, conserveMomentum) ), "PDF Restore" );
////////////////////////
// EXECUTE SIMULATION //
////////////////////////
WcTimingPool timeloopTiming;
std::vector< std::vector<std::string> > timerKeys;
std::vector<std::string> LBMTimer;
LBMTimer.emplace_back("Stream&Collide");
LBMTimer.emplace_back("Stream");
LBMTimer.emplace_back("Collide");
timerKeys.push_back(LBMTimer);
std::vector<std::string> bhTimer;
bhTimer.emplace_back("Boundary Handling");
timerKeys.push_back(bhTimer);
std::vector<std::string> couplingTimer1;
couplingTimer1.emplace_back("Particle Mapping");
std::vector<std::string> couplingTimer2;
couplingTimer2.emplace_back("PDF Restore");
timerKeys.push_back(couplingTimer1);
timerKeys.push_back(couplingTimer2);
std::vector<std::string> rpdTimer;
rpdTimer.emplace_back("RPD Force");
rpdTimer.emplace_back("RPD VV1");
rpdTimer.emplace_back("RPD VV2");
rpdTimer.emplace_back("RPD Lub");
rpdTimer.emplace_back("RPD Collision");
timerKeys.push_back(rpdTimer);
uint_t numCells = uint_t(0);
uint_t numFluidCells = uint_t(0);
uint_t numNBCells = uint_t(0);
uint_t numLocalParticles = uint_t(0);
uint_t numGhostParticles = uint_t(0);
uint_t numContacts = uint_t(0);
std::vector<double> timings(timerKeys.size());
// every rank writes its own file -> numProcesses number of samples!
int myRank = MPIManager::instance()->rank();
std::string logFileName = baseFolder + "/load";
logFileName += "_settling";
logFileName += "_spheres";
logFileName += "_d" + std::to_string(int_c(std::ceil(diameter)));
logFileName += "_bs" + std::to_string(blockSize);
logFileName += "_" + std::to_string(myRank) + ".txt";
std::ofstream file;
if(!noFileOutput)
{
WALBERLA_LOG_INFO_ON_ROOT("Writing load info to file " << logFileName);
file.open( logFileName.c_str());
file << "# svf = " << solidVolumeFraction << ", d = " << diameter << ", domain = " << XCells << "x" << YCells << "x" << ZCells << "\n";
}
auto timeStepOfFirstTiming = uint_t(50);
if( timesteps - timeStepOfFirstTiming < numSamples )
{
WALBERLA_LOG_WARNING_ON_ROOT("Less actual time steps than number of required samples!");
}
uint_t nSample( 0 ); // number of current sample
real_t samplingFrequency = real_c(timesteps - timeStepOfFirstTiming) / real_c(numSamples);
// time loop
for (uint_t i = 0; i < timesteps; ++i )
{
// perform a single simulation step
timeloop.singleStep( timeloopTiming );
reduceProperty.operator()<mesa_pd::HydrodynamicForceTorqueNotification>(*ps);
timeloopTiming["RPD Force"].start();
if( i == 0 )
{
lbm_mesapd_coupling::InitializeHydrodynamicForceTorqueForAveragingKernel initializeHydrodynamicForceTorqueForAveragingKernel;
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
}
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, averageHydrodynamicForceTorque, *accessor );
timeloopTiming["RPD Force"].end();
for(auto subCycle = uint_t(0); subCycle < numRPDSubCycles; ++subCycle )
{
timeloopTiming["RPD VV1"].start();
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, vvIntegratorPreForce, *accessor);
timeloopTiming["RPD VV1"].end();
syncCall();
// lubrication correction
timeloopTiming["RPD Lub"].start();
ps->forEachParticlePairHalf(useOpenMP, mesa_pd::kernel::ExcludeInfiniteInfinite(), *accessor,
[&lubricationCorrectionKernel,maximumLubricationCutOffDistance, &rpdDomain,&numContacts]
(const size_t idx1, const size_t idx2, auto& ac)
{
mesa_pd::collision_detection::AnalyticContactDetection acd;
acd.getContactThreshold() = maximumLubricationCutOffDistance;
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());
++numContacts; // this then also includes all actual contacts since there the contact threshold is smaller
}
}
},
*accessor );
timeloopTiming["RPD Lub"].end();
// collision response
timeloopTiming["RPD Collision"].start();
ps->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 );
timeloopTiming["RPD Collision"].end();
reduceAndSwapContactHistory(*ps);
timeloopTiming["RPD Force"].start();
lbm_mesapd_coupling::AddHydrodynamicInteractionKernel addHydrodynamicInteraction;
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addHydrodynamicInteraction, *accessor );
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addGravitationalForce, *accessor );
timeloopTiming["RPD Force"].end();
reduceProperty.operator()<mesa_pd::ForceTorqueNotification>(*ps);
// integration
timeloopTiming["RPD VV2"].start();
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, vvIntegratorPostForce, *accessor);
timeloopTiming["RPD VV2"].end();
syncCall();
}
timeloopTiming["RPD Force"].start();
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, resetHydrodynamicForceTorque, *accessor );
timeloopTiming["RPD Force"].end();
// update particle mapping
timeloopAfterParticles.singleStep(timeloopTiming);
//WALBERLA_LOG_INFO_ON_ROOT(timeloopTiming);
// check if current time step should be included in sample
if( i >= uint_c( samplingFrequency * real_c(nSample) ) + timeStepOfFirstTiming )
{
// include -> evaluate all timers and quantities
evaluateFluidQuantities<BoundaryHandling_T>(blocks, boundaryHandlingID, numCells, numFluidCells, numNBCells);
evaluateRPDQuantities(ps, numLocalParticles, numGhostParticles);
evaluateTimers(timeloopTiming, timerKeys, timings);
if(!noFileOutput)
{
auto totalTime = std::accumulate(timings.begin(), timings.end(), 0.0 );
file << timeloop.getCurrentTimeStep() << " " << real_c(timeloop.getCurrentTimeStep()) / Tref << " "
<< numCells << " " << numFluidCells << " " << numNBCells << " "
<< numLocalParticles << " " << numGhostParticles << " "
<< real_c(numContacts) / real_c(numRPDSubCycles) << " " << numRPDSubCycles;
for (auto timing : timings) {
file << " " << timing;
}
file << " " << totalTime << "\n";
}
++nSample;
}
numCells = uint_t(0);
numFluidCells = uint_t(0);
numNBCells = uint_t(0);
numLocalParticles = uint_t(0);
numGhostParticles = uint_t(0);
numContacts = uint_t(0);
// reset timers to always include only a single time step in them
timeloopTiming.clear();
}
if(!noFileOutput) {
file.close();
}
WALBERLA_LOG_INFO_ON_ROOT("Simulation finished!");
return 0;
}
}
int main( int argc, char **argv ){
fluid_particle_workload_evaluation::main(argc, argv);
}
apps/benchmarks/FluidParticleCouplingWithLoadBalancing/Utility.h 0 → 100644
View file @ 2dee476c
#pragma once
#include "lbm_mesapd_coupling/amr/BlockInfo.h"
namespace walberla {
namespace lbm_mesapd_coupling {
namespace amr {
/*
* Result from the workload evaluation as described in
* Rettinger, Ruede - "Dynamic Load Balancing Techniques for Particulate Flow Simulations", 2019, Computation
*/
inline real_t fittedLBMWeightEvaluationFunction(const BlockInfo& blockInfo)
{
uint_t Ce = blockInfo.numberOfCells;
uint_t F = blockInfo.numberOfFluidCells;
real_t weight = real_t(7.597476065046571e-06) * real_c(Ce) + real_t(8.95723566283202e-05) * real_c(F) + real_t(-0.1526111388616016);
return std::max(weight,real_t(0));
}
inline real_t fittedBHWeightEvaluationFunction(const BlockInfo& blockInfo)
{
uint_t Ce = blockInfo.numberOfCells;
uint_t NB = blockInfo.numberOfNearBoundaryCells;
real_t weight = real_t(1.3067711379655123e-07) * real_c(Ce) + real_t(0.0007289549127205142) * real_c(NB) + real_t(-0.1575698071795788);
return std::max(weight,real_t(0));
}
inline real_t fittedRPDWeightEvaluationFunction(const BlockInfo& blockInfo)
{
uint_t Pl = blockInfo.numberOfLocalParticles;
uint_t Pg = blockInfo.numberOfGhostParticles;
uint_t Sc = blockInfo.numberOfRPDSubCycles;
real_t cPlPg2 = real_t(2.402288635599054e-05);
real_t cPl = real_t(0.00040932622363097144);
real_t cPg = real_t(0.0007268941363125683);
real_t c = real_t(2.01883028312316e-19);
real_t weight = real_c(Sc) * ( cPlPg2 * real_c(Pl+Pg) * real_c(Pl+Pg) + cPl * real_c(Pl) + cPg * real_c(Pg) + c );
return std::max(weight,real_t(0));
}
inline real_t fittedCoup1WeightEvaluationFunction(const BlockInfo& blockInfo)
{
uint_t Ce = blockInfo.numberOfCells;
uint_t F = blockInfo.numberOfFluidCells;
uint_t Pl = blockInfo.numberOfLocalParticles;
uint_t Pg = blockInfo.numberOfGhostParticles;
real_t weight = real_t(5.610203409278647e-06) * real_c(Ce) + real_t(-7.257635845636656e-07) * real_c(F) + real_t(0.02049703546054693) * real_c(Pl) + real_t(0.04248208493809902) * real_c(Pg) + real_t(-0.26609470510074784);
return std::max(weight,real_t(0));
}
inline real_t fittedCoup2WeightEvaluationFunction(const BlockInfo& blockInfo)
{
uint_t Ce = blockInfo.numberOfCells;
uint_t F = blockInfo.numberOfFluidCells;
uint_t Pl = blockInfo.numberOfLocalParticles;
uint_t Pg = blockInfo.numberOfGhostParticles;
real_t weight = real_t(7.198479654682179e-06) * real_c(Ce) + real_t(1.178247475854302e-06) * real_c(F) + real_t(-0.0026401549115124628) * real_c(Pl) + real_t(0.008459646786179298) * real_c(Pg) + real_t(-0.001077320113275954);
return std::max(weight,real_t(0));
}
inline real_t fittedTotalWeightEvaluationFunction(const BlockInfo& blockInfo)
{
return fittedLBMWeightEvaluationFunction(blockInfo) + fittedBHWeightEvaluationFunction(blockInfo) +
fittedRPDWeightEvaluationFunction(blockInfo) + fittedCoup1WeightEvaluationFunction(blockInfo) +
fittedCoup2WeightEvaluationFunction(blockInfo);
}
} //namespace amr
} //namespace lbm_mesapd_coupling
} //namespace walberla
apps/benchmarks/FluidizedBed/CMakeLists.txt 0 → 100644
View file @ 2dee476c
waLBerla_link_files_to_builddir("*.prm")
if (WALBERLA_BUILD_WITH_GPU_SUPPORT AND WALBERLA_BUILD_WITH_CODEGEN AND (CMAKE_CUDA_ARCHITECTURES GREATER_EQUAL 60 OR WALBERLA_BUILD_WITH_HIP))
waLBerla_add_executable(NAME FluidizedBed_PSM_GPU FILES FluidizedBedGPU.cpp
DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::gpu walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd walberla::timeloop walberla::vtk PSMCodegenPython_srt_sc1 )
endif ()
apps/benchmarks/FluidizedBed/FluidizedBedGPU.cpp 0 → 100644
View file @ 2dee476c
//======================================================================================================================
//
// 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 FluidizedBedGPU.cpp
//! \ingroup lbm_mesapd_coupling
//! \author Samuel Kemmler <samuel.kemmler@fau.de>
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//! \brief Modification of showcases/FluidizedBed/FluidizedBedPSM.cpp
//
//======================================================================================================================
#include "blockforest/Initialization.h"
#include "core/DataTypes.h"
#include "core/Environment.h"
#include "core/debug/Debug.h"
#include "core/grid_generator/SCIterator.h"
#include "core/logging/all.h"
#include "core/math/all.h"
#include "core/mpi/Broadcast.h"
#include "core/timing/RemainingTimeLogger.h"
#include "field/AddToStorage.h"
#include "field/vtk/all.h"
#include "geometry/InitBoundaryHandling.h"
#include "gpu/AddGPUFieldToStorage.h"
#include "gpu/DeviceSelectMPI.h"
#include "gpu/communication/UniformGPUScheme.h"
#include "lbm/PerformanceLogger.h"
#include "lbm/vtk/all.h"
#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 "vtk/all.h"
#include "InitializeDomainForPSM.h"
#include "PSMPackInfo.h"
#include "PSMSweepSplit.h"
#include "PSM_Density.h"
#include "PSM_InfoHeader.h"
#include "PSM_MacroGetter.h"
#include "PSM_NoSlip.h"
#include "PSM_UBB.h"
namespace fluidized_bed
{
///////////
// USING //
///////////
using namespace walberla;
using walberla::uint_t;
using flag_t = walberla::uint8_t;
using FlagField_T = FlagField< flag_t >;
using namespace lbm_mesapd_coupling::psm::gpu;
using PackInfo_T = pystencils::PSMPackInfo;
///////////
// FLAGS //
///////////
const FlagUID Fluid_Flag("Fluid");
const FlagUID NoSlip_Flag("NoSlip");
const FlagUID Inflow_Flag("Inflow");
const FlagUID Outflow_Flag("Outflow");
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;
}
struct FluidInfo
{
uint_t numFluidCells = 0;
real_t averageVelocity = 0_r;
real_t maximumVelocity = 0_r;
real_t averageDensity = 0_r;
real_t maximumDensity = 0_r;
void allReduce()
{
walberla::mpi::allReduceInplace(numFluidCells, walberla::mpi::SUM);
walberla::mpi::allReduceInplace(averageVelocity, walberla::mpi::SUM);
walberla::mpi::allReduceInplace(maximumVelocity, walberla::mpi::MAX);
;
walberla::mpi::allReduceInplace(averageDensity, walberla::mpi::SUM);
walberla::mpi::allReduceInplace(maximumDensity, walberla::mpi::MAX);
averageVelocity /= real_c(numFluidCells);
averageDensity /= real_c(numFluidCells);
}
};
std::ostream& operator<<(std::ostream& os, FluidInfo const& m)
{
return os << "Fluid Info: numFluidCells = " << m.numFluidCells << ", uAvg = " << m.averageVelocity
<< ", uMax = " << m.maximumVelocity << ", densityAvg = " << m.averageDensity
<< ", densityMax = " << m.maximumDensity;
}
FluidInfo evaluateFluidInfo(const shared_ptr< StructuredBlockStorage >& blocks, const BlockDataID& densityFieldID,
const BlockDataID& velocityFieldID)
{
FluidInfo info;
for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
{
auto densityField = blockIt->getData< DensityField_T >(densityFieldID);
auto velocityField = blockIt->getData< VelocityField_T >(velocityFieldID);
WALBERLA_FOR_ALL_CELLS_XYZ(
densityField, ++info.numFluidCells; Vector3< real_t > velocity(
velocityField->get(x, y, z, 0), velocityField->get(x, y, z, 1), velocityField->get(x, y, z, 2));
real_t density = densityField->get(x, y, z); real_t velMagnitude = velocity.length();
info.averageVelocity += velMagnitude; info.maximumVelocity = std::max(info.maximumVelocity, velMagnitude);
info.averageDensity += density; info.maximumDensity = std::max(info.maximumDensity, density);)
}
info.allReduce();
return info;
}
//////////
// MAIN //
//////////
//*******************************************************************************************************************
/*!\brief Basic simulation of a fluidization setup
*
* Initially, the mono-sized sphere are created on a structured grid inside the domain.
* The domain is either periodic or bounded by walls in the horizontal directions (x and y).
* In z-direction, a constant inflow from below is provided
* and a pressure boundary condition is set at the top, resembling an outflow boundary.
*
* The simulation is run for the given number of seconds (runtime).
*
* All parameters should be set via the input file.
*
* For the overall algorithm and the different model parameters, see
* Rettinger, Rüde - An efficient four-way coupled lattice Boltzmann - discrete element method for
* fully resolved simulations of particle-laden flows (2020, preprint: https://arxiv.org/abs/2003.01490)
*
*/
//*******************************************************************************************************************
int main(int argc, char** argv)
{
Environment env(argc, argv);
gpu::selectDeviceBasedOnMpiRank();
auto cfgFile = env.config();
if (!cfgFile) { WALBERLA_ABORT("Usage: " << argv[0] << " path-to-configuration-file \n"); }
WALBERLA_LOG_INFO_ON_ROOT("waLBerla revision: " << std::string(WALBERLA_GIT_SHA1).substr(0, 8));
WALBERLA_LOG_INFO_ON_ROOT("compiler flags: " << std::string(WALBERLA_COMPILER_FLAGS));
WALBERLA_LOG_INFO_ON_ROOT("build machine: " << std::string(WALBERLA_BUILD_MACHINE));
WALBERLA_LOG_INFO_ON_ROOT(*cfgFile);
// read all parameters from the config file
Config::BlockHandle physicalSetup = cfgFile->getBlock("PhysicalSetup");
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");
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 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()) << ")");
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 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");
Config::BlockHandle outputSetup = cfgFile->getBlock("Output");
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(
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
const real_t dt_SI = uInflow / uInflow_SI * dx_SI;
const real_t diameter = 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 * diameter * diameter * diameter;
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 * diameter;
WALBERLA_LOG_INFO_ON_ROOT("Simulation setup:");
WALBERLA_LOG_INFO_ON_ROOT(" - particles: diameter = " << diameter << ", densityRatio = " << densityRatio);
WALBERLA_LOG_INFO_ON_ROOT(" - fluid: kin. visc = " << viscosity << ", relaxation rate = " << omega);
WALBERLA_LOG_INFO_ON_ROOT(" - grav. acceleration = " << gravitationalAcceleration);
WALBERLA_LOG_INFO_ON_ROOT(" - Galileo number = " << GalileiNumber);
WALBERLA_LOG_INFO_ON_ROOT(" - particle Reynolds number = " << ReynoldsNumberParticle);
WALBERLA_LOG_INFO_ON_ROOT(" - domain size = " << domainSize);
WALBERLA_LOG_INFO_ON_ROOT(" - cells per blocks per direction = " << cellsPerBlockPerDirection);
WALBERLA_LOG_INFO_ON_ROOT(" - dx = " << dx_SI << " m");
WALBERLA_LOG_INFO_ON_ROOT(" - dt = " << dt_SI << " s");
WALBERLA_LOG_INFO_ON_ROOT(" - total time steps = " << numTimeSteps);
WALBERLA_LOG_INFO_ON_ROOT(" - particle generation spacing = " << particleGenerationSpacing);
WALBERLA_LOG_INFO_ON_ROOT(" - info spacing = " << infoSpacing);
WALBERLA_LOG_INFO_ON_ROOT(" - vtk spacing particles = " << vtkSpacingParticles
<< ", fluid slice = " << vtkSpacingFluid);
///////////////////////////
// BLOCK STRUCTURE SETUP //
///////////////////////////
const bool periodicInZ = false;
shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid(
numXBlocks, numYBlocks, numZBlocks, cellsPerBlockPerDirection[0], cellsPerBlockPerDirection[1],
cellsPerBlockPerDirection[2], dx, 0, false, false, periodicInX, periodicInY, periodicInZ, // periodicity
false);
auto simulationDomain = blocks->getDomain();
//////////////////
// RPD COUPLING //
//////////////////
auto rpdDomain = std::make_shared< mesa_pd::domain::BlockForestDomain >(blocks->getBlockForestPointer());
// init data structures
auto ps = walberla::make_shared< mesa_pd::data::ParticleStorage >(1);
auto ss = walberla::make_shared< mesa_pd::data::ShapeStorage >();
using ParticleAccessor_T = mesa_pd::data::ParticleAccessorWithShape;
auto accessor = walberla::make_shared< ParticleAccessor_T >(ps, ss);
// prevent particles from interfering with inflow and outflow by putting the bounding planes slightly in front
const real_t planeOffsetFromInflow = dx;
const real_t planeOffsetFromOutflow = dx;
createPlaneSetup(ps, ss, simulationDomain, periodicInX, periodicInY, planeOffsetFromInflow, planeOffsetFromOutflow);
auto sphereShape = ss->create< mesa_pd::data::Sphere >(diameter * real_t(0.5));
ss->shapes[sphereShape]->updateMassAndInertia(densityParticle);
// create spheres
auto generationDomain = simulationDomain.getExtended(-particleGenerationSpacing * 0.5_r);
for (auto pt : grid_generator::SCGrid(generationDomain, generationDomain.center(), particleGenerationSpacing))
{
if (rpdDomain->isContainedInProcessSubdomain(uint_c(mpi::MPIManager::instance()->rank()), pt))
{
mesa_pd::data::Particle&& p = *ps->create();
p.setPosition(pt);
p.setInteractionRadius(diameter * real_t(0.5));
p.setOwner(mpi::MPIManager::instance()->rank());
p.setShapeID(sphereShape);
p.setType(1);
p.setLinearVelocity(0.1_r * Vector3< real_t >(math::realRandom(
-uInflow, uInflow))); // set small initial velocity to break symmetries
}
}
///////////////////////
// ADD DATA TO BLOCKS //
////////////////////////
// add PDF field
BlockDataID pdfFieldID =
field::addToStorage< PdfField_T >(blocks, "pdf field (fzyx)", real_c(std::nan("")), field::fzyx);
BlockDataID pdfFieldGPUID = gpu::addGPUFieldToStorage< PdfField_T >(blocks, pdfFieldID, "pdf field GPU");
BlockDataID densityFieldID = field::addToStorage< DensityField_T >(blocks, "Density", real_t(0), field::fzyx);
BlockDataID velFieldID = field::addToStorage< VelocityField_T >(blocks, "Velocity", real_t(0), field::fzyx);
BlockDataID BFieldID =
field::addToStorage< lbm_mesapd_coupling::psm::gpu::BField_T >(blocks, "B field", 0, field::fzyx, 1);
// add flag field
BlockDataID flagFieldID = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
// 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 = " Boundaries"
"{"
"Border { direction T; walldistance -1; flag Outflow; }"
"Border { direction B; walldistance -1; flag Inflow; }";
if (!periodicInX)
{
boundariesBlockString += "Border { direction W; walldistance -1; flag NoSlip; }"
"Border { direction E; walldistance -1; flag NoSlip; }";
}
if (!periodicInY)
{
boundariesBlockString += "Border { direction S; walldistance -1; flag NoSlip; }"
"Border { direction N; walldistance -1; flag NoSlip; }";
}
boundariesBlockString += "}";
WALBERLA_ROOT_SECTION()
{
std::ofstream boundariesFile("boundaries.prm");
boundariesFile << boundariesBlockString;
boundariesFile.close();
}
WALBERLA_MPI_BARRIER()
auto boundariesCfgFile = Config();
boundariesCfgFile.readParameterFile("boundaries.prm");
auto boundariesConfig = boundariesCfgFile.getBlock("Boundaries");
// map boundaries into the LBM simulation
geometry::initBoundaryHandling< FlagField_T >(*blocks, flagFieldID, boundariesConfig);
geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldID, Fluid_Flag);
lbm::PSM_NoSlip noSlip(blocks, pdfFieldGPUID);
noSlip.fillFromFlagField< FlagField_T >(blocks, flagFieldID, NoSlip_Flag, Fluid_Flag);
lbm::PSM_UBB ubb(blocks, pdfFieldGPUID, inflowVec[0], inflowVec[1], inflowVec[2]);
ubb.fillFromFlagField< FlagField_T >(blocks, flagFieldID, Inflow_Flag, Fluid_Flag);
lbm::PSM_Density density_bc(blocks, pdfFieldGPUID, real_t(1));
density_bc.fillFromFlagField< FlagField_T >(blocks, flagFieldID, Outflow_Flag, Fluid_Flag);
///////////////
// TIME LOOP //
///////////////
// map particles into the LBM simulation
// note: planes are not mapped and are thus only visible to the particles, not to the fluid
// instead, the respective boundary conditions for the fluid are explicitly set, see the boundary handling
ParticleAndVolumeFractionSoA_T< 1 > particleAndVolumeFractionSoA(blocks, omega);
PSMSweepCollection psmSweepCollection(blocks, accessor, lbm_mesapd_coupling::RegularParticlesSelector(),
particleAndVolumeFractionSoA, particleSubBlockSize);
for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
{
psmSweepCollection.particleMappingSweep(&(*blockIt));
}
pystencils::InitializeDomainForPSM pdfSetter(
particleAndVolumeFractionSoA.BsFieldID, particleAndVolumeFractionSoA.BFieldID,
particleAndVolumeFractionSoA.particleVelocitiesFieldID, pdfFieldGPUID, real_t(0), real_t(0), real_t(0),
real_t(1.0), real_t(0), real_t(0), real_t(0));
for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
{
// pdfSetter requires particle velocities at cell centers
psmSweepCollection.setParticleVelocitiesSweep(&(*blockIt));
pdfSetter(&(*blockIt));
}
// setup of the LBM communication for synchronizing the pdf field between neighboring blocks
gpu::communication::UniformGPUScheme< Stencil_T > com(blocks, true, false);
com.addPackInfo(make_shared< PackInfo_T >(pdfFieldGPUID));
auto communication = std::function< void() >([&]() { com.communicate(); });
// create the timeloop
SweepTimeloop commTimeloop(blocks->getBlockStorage(), numTimeSteps);
SweepTimeloop timeloop(blocks->getBlockStorage(), numTimeSteps);
timeloop.addFuncBeforeTimeStep(RemainingTimeLogger(timeloop.getNrOfTimeSteps()), "Remaining Time Logger");
pystencils::PSM_MacroGetter getterSweep(BFieldID, densityFieldID, pdfFieldID, velFieldID, real_t(0.0), real_t(0.0),
real_t(0.0));
// vtk output
if (vtkSpacingParticles != uint_t(0))
{
// sphere
auto particleVtkOutput = make_shared< mesa_pd::vtk::ParticleVtkOutput >(ps);
particleVtkOutput->addOutput< mesa_pd::data::SelectParticleUid >("uid");
particleVtkOutput->addOutput< mesa_pd::data::SelectParticleLinearVelocity >("velocity");
particleVtkOutput->addOutput< mesa_pd::data::SelectParticleInteractionRadius >("radius");
// limit output to process-local spheres
particleVtkOutput->setParticleSelector([sphereShape](const mesa_pd::data::ParticleStorage::iterator& pIt) {
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", vtkSpacingParticles, vtkFolder);
timeloop.addFuncBeforeTimeStep(vtk::writeFiles(particleVtkWriter), "VTK (sphere data)");
}
if (vtkSpacingFluid != uint_t(0))
{
// velocity field, only a slice
auto pdfFieldVTK = vtk::createVTKOutput_BlockData(blocks, "fluid", vtkSpacingFluid, 0, false, vtkFolder);
pdfFieldVTK->addBeforeFunction(communication);
pdfFieldVTK->addBeforeFunction([&]() {
gpu::fieldCpy< PdfField_T, gpu::GPUField< real_t > >(blocks, pdfFieldID, pdfFieldGPUID);
gpu::fieldCpy< GhostLayerField< real_t, 1 >, BFieldGPU_T >(blocks, BFieldID,
particleAndVolumeFractionSoA.BFieldID);
for (auto& block : *blocks)
getterSweep(&block);
});
AABB sliceAABB(real_t(0), real_c(domainSize[1]) * real_t(0.5) - real_t(1), real_t(0), real_c(domainSize[0]),
real_c(domainSize[1]) * real_t(0.5) + real_t(1), real_c(domainSize[2]));
vtk::AABBCellFilter aabbSliceFilter(sliceAABB);
field::FlagFieldCellFilter< FlagField_T > fluidFilter(flagFieldID);
fluidFilter.addFlag(Fluid_Flag);
vtk::ChainedFilter combinedSliceFilter;
combinedSliceFilter.addFilter(fluidFilter);
combinedSliceFilter.addFilter(aabbSliceFilter);
pdfFieldVTK->addCellInclusionFilter(combinedSliceFilter);
pdfFieldVTK->addCellDataWriter(make_shared< field::VTKWriter< VelocityField_T > >(velFieldID, "Velocity"));
pdfFieldVTK->addCellDataWriter(make_shared< field::VTKWriter< DensityField_T > >(densityFieldID, "Density"));
pdfFieldVTK->addCellDataWriter(make_shared< field::VTKWriter< BField_T > >(BFieldID, "Fraction mapping field B"));
timeloop.addFuncBeforeTimeStep(vtk::writeFiles(pdfFieldVTK), "VTK (fluid field data)");
}
if (vtkSpacingFluid != uint_t(0) || vtkSpacingParticles != uint_t(0))
{
vtk::writeDomainDecomposition(blocks, "domain_decomposition", vtkFolder);
}
// add performance logging
const lbm::PerformanceLogger< FlagField_T > performanceLogger(blocks, flagFieldID, Fluid_Flag,
performanceLogFrequency);
timeloop.addFuncAfterTimeStep(performanceLogger, "Evaluate performance logging");
// add LBM communication function and boundary handling sweep
timeloop.add() << Sweep(deviceSyncWrapper(ubb.getSweep()), "Boundary Handling (UBB)");
timeloop.add() << Sweep(deviceSyncWrapper(density_bc.getSweep()), "Boundary Handling (Density)");
timeloop.add() << Sweep(deviceSyncWrapper(noSlip.getSweep()), "Boundary Handling (NoSlip)");
// stream + collide LBM step
pystencils::PSMSweepSplit PSMSweep(particleAndVolumeFractionSoA.BsFieldID, particleAndVolumeFractionSoA.BFieldID,
particleAndVolumeFractionSoA.particleForcesFieldID,
particleAndVolumeFractionSoA.particleVelocitiesFieldID, pdfFieldGPUID,
real_t(0.0), real_t(0.0), real_t(0.0), omega);
addPSMSweepsToTimeloops(commTimeloop, timeloop, com, psmSweepCollection, PSMSweep);
////////////////////////
// EXECUTE SIMULATION //
////////////////////////
WcTimingPool timeloopTiming;
const bool useOpenMP = true;
real_t linkedCellWidth = linkedCellWidthRation * diameter;
mesa_pd::data::LinkedCells linkedCells(rpdDomain->getUnionOfLocalAABBs().getExtended(linkedCellWidth),
linkedCellWidth);
mesa_pd::kernel::InsertParticleIntoLinkedCells ipilc;
// time loop
for (uint_t timeStep = 0; timeStep < numTimeSteps; ++timeStep)
{
// perform a single simulation step -> this contains LBM and setting of the hydrodynamic interactions
commTimeloop.singleStep(timeloopTiming);
timeloop.singleStep(timeloopTiming);
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)
{
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();
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)
{
timeloopTiming["Evaluate infos"].start();
auto particleInfo = evaluateParticleInfo(*accessor);
WALBERLA_LOG_INFO_ON_ROOT(particleInfo);
auto fluidInfo = evaluateFluidInfo(blocks, densityFieldID, velFieldID);
WALBERLA_LOG_INFO_ON_ROOT(fluidInfo);
if (particleBarriers) WALBERLA_MPI_BARRIER();
timeloopTiming["Evaluate infos"].end();
}
}
timeloopTiming.logResultOnRoot();
return EXIT_SUCCESS;
}
} // namespace fluidized_bed
int main(int argc, char** argv) { fluidized_bed::main(argc, argv); }
apps/benchmarks/FluidizedBed/input.prm 0 → 100644
View file @ 2dee476c
PhysicalSetup // all to be specified in SI units!
{
xSize 0.05; // = width
ySize 0.02; // = depth
zSize 0.08; // = height
periodicInX false;
periodicInY false;
runtime 0.1;
uInflow 0.005;
gravitationalAcceleration 9.81;
kinematicViscosityFluid 1e-5;
densityFluid 1000.;
particleDiameter 0.002;
densityParticle 1100.;
dynamicFrictionCoefficient 0.15;
coefficientOfRestitution 0.6;
collisionTimeFactor 1.0;
particleGenerationSpacing 0.00401; // 0.00401 or 0.00201
}
NumericalSetup
{
dx 0.0001; // in m
uInflow 0.01; // in LBM units, should be smaller than 0.1, this then determines dt
// product of number of blocks should be equal to number of used processes
numXBlocks 1;
numYBlocks 1;
numZBlocks 1;
useLubricationForces true;
numberOfParticleSubCycles 10;
particleSubBlockSize <10, 10, 10>;
linkedCellWidthRation 1.01;
particleBarriers true;
}
Output
{
infoSpacing 0.0; // in s
vtkSpacingParticles 0.0; // in s
vtkSpacingFluid 0.0; // in s
vtkFolder vtk_out;
performanceLogFrequency 500;
}
apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/CMakeLists.txt
View file @ 2dee476c
waLBerla_add_executable ( NAME ForcesOnSphereNearPlaneInShearFlow
DEPENDS blockforest boundary core domain_decomposition field lbm pe pe_coupling postprocessing timeloop vtk )
\ No newline at end of file
DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::pe walberla::pe_coupling walberla::postprocessing walberla::timeloop walberla::vtk )
\ No newline at end of file
apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/ForcesOnSphereNearPlaneInShearFlow.cpp
View file @ 2dee476c
......@@ -81,23 +81,23 @@ using walberla::uint_t;
//////////////
// PDF field, flag field & body field
typedef lbm::D3Q19< lbm::collision_model::TRT, false > LatticeModel_T;
using LatticeModel_T = lbm::D3Q19<lbm::collision_model::TRT, false>;
using Stencil_T = LatticeModel_T::Stencil;
using PdfField_T = lbm::PdfField<LatticeModel_T>;
using flag_t = walberla::uint8_t;
using FlagField_T = FlagField<flag_t>;
typedef GhostLayerField< pe::BodyID, 1 > BodyField_T;
using BodyField_T = GhostLayerField<pe::BodyID, 1>;
const uint_t FieldGhostLayers = 4;
// boundary handling
typedef pe_coupling::SimpleBB< LatticeModel_T, FlagField_T > MO_SBB_T;
typedef pe_coupling::CurvedLinear< LatticeModel_T, FlagField_T > MO_CLI_T;
using MO_SBB_T = pe_coupling::SimpleBB<LatticeModel_T, FlagField_T>;
using MO_CLI_T = pe_coupling::CurvedLinear<LatticeModel_T, FlagField_T>;
typedef BoundaryHandling< FlagField_T, Stencil_T, MO_SBB_T, MO_CLI_T > BoundaryHandling_T;
using BoundaryHandling_T = BoundaryHandling<FlagField_T, Stencil_T, MO_SBB_T, MO_CLI_T>;
typedef std::tuple< pe::Sphere, pe::Plane > BodyTypeTuple;
using BodyTypeTuple = std::tuple<pe::Sphere, pe::Plane>;
///////////
// FLAGS //
......@@ -291,13 +291,8 @@ public:
WALBERLA_MPI_SECTION()
{
mpi::allReduceInplace( force[0], mpi::SUM );
mpi::allReduceInplace( force[1], mpi::SUM );
mpi::allReduceInplace( force[2], mpi::SUM );
mpi::allReduceInplace( torque[0], mpi::SUM );
mpi::allReduceInplace( torque[1], mpi::SUM );
mpi::allReduceInplace( torque[2], mpi::SUM );
mpi::allReduceInplace( force, mpi::SUM );
mpi::allReduceInplace( torque, mpi::SUM );
}
if( fileIO_ )
......@@ -609,14 +604,14 @@ int main( int argc, char **argv )
LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( omega, lbm::collision_model::TRT::threeSixteenth, finestLevel ) );
// add PDF field
BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (zyxf)", latticeModel,
BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (fzyx)", latticeModel,
Vector3< real_t >( real_t(0) ), real_t(1),
FieldGhostLayers, field::zyxf );
FieldGhostLayers, field::fzyx );
// add flag field
BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>( blocks, "flag field", FieldGhostLayers );
// add body field
BlockDataID bodyFieldID = field::addToStorage<BodyField_T>( blocks, "body field", nullptr, field::zyxf, FieldGhostLayers );
BlockDataID bodyFieldID = field::addToStorage<BodyField_T>( blocks, "body field", nullptr, field::fzyx, FieldGhostLayers );
// add boundary handling
BlockDataID boundaryHandlingID = blocks->addStructuredBlockData< BoundaryHandling_T >( MyBoundaryHandling( flagFieldID, pdfFieldID, bodyFieldID ), "boundary handling" );
......
apps/benchmarks/FreeSurfaceAdvection/CMakeLists.txt 0 → 100644
View file @ 2dee476c
waLBerla_link_files_to_builddir( *.prm )
waLBerla_add_executable(NAME DeformationField
FILES DeformationField.cpp
DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::postprocessing walberla::timeloop walberla::vtk )
waLBerla_add_executable(NAME SingleVortex
FILES SingleVortex.cpp
DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::postprocessing walberla::timeloop walberla::vtk )
waLBerla_add_executable(NAME ZalesakDisk
FILES ZalesakDisk.cpp
DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::postprocessing walberla::timeloop walberla::vtk )
\ No newline at end of file
apps/benchmarks/FreeSurfaceAdvection/DeformationField.cpp 0 → 100644
View file @ 2dee476c
//======================================================================================================================
//
// 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 DeformationField.cpp
//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
//
// This benchmark simulates the advection of a spherical bubble in a vortex field with very high deformation. The vortex
// field changes periodically so that the bubble returns to its initial position, where it should take its initial form.
// The relative geometrical error of the bubble's shape after one period is evaluated. There is no LBM flow simulation
// performed here. It is a test case for the FSLBM's mass advection only. This benchmark is based on the work from
// Liovic et al. (doi: 10.1016/j.compfluid.2005.09.003). The setup is identical as in Viktor Haag's master thesis
// (https://www10.cs.fau.de/publications/theses/2017/Haag_MT_2017.pdf).
//======================================================================================================================
#include "core/Environment.h"
#include "field/Gather.h"
#include "lbm/PerformanceLogger.h"
#include "lbm/field/AddToStorage.h"
#include "lbm/free_surface/LoadBalancing.h"
#include "lbm/free_surface/SurfaceMeshWriter.h"
#include "lbm/free_surface/TotalMassComputer.h"
#include "lbm/free_surface/VtkWriter.h"
#include "lbm/free_surface/bubble_model/Geometry.h"
#include "lbm/free_surface/surface_geometry/SurfaceGeometryHandler.h"
#include "lbm/free_surface/surface_geometry/Utility.h"
#include "lbm/lattice_model/D3Q19.h"
#include "functionality/AdvectionDynamicsHandler.h"
#include "functionality/GeometricalErrorEvaluator.h"
namespace walberla
{
namespace free_surface
{
namespace DeformationField
{
using ScalarField_T = GhostLayerField< real_t, 1 >;
using VectorField_T = GhostLayerField< Vector3< real_t >, 1 >;
// Lattice model is only created for dummy purposes; no LBM simulation is performed
using CollisionModel_T = lbm::collision_model::SRT;
using LatticeModel_T = lbm::D3Q19< CollisionModel_T, true >;
using LatticeModelStencil_T = LatticeModel_T::Stencil;
using PdfField_T = lbm::PdfField< LatticeModel_T >;
using PdfCommunication_T = blockforest::SimpleCommunication< LatticeModelStencil_T >;
// the geometry computations in SurfaceGeometryHandler require meaningful values in the ghost layers in corner
// directions (flag field and fill level field); this holds, even if the lattice model uses a D3Q19 stencil
using CommunicationStencil_T =
typename std::conditional< LatticeModel_T::Stencil::D == uint_t(2), stencil::D2Q9, stencil::D3Q27 >::type;
using Communication_T = blockforest::SimpleCommunication< CommunicationStencil_T >;
using flag_t = uint32_t;
using FlagField_T = FlagField< flag_t >;
using FreeSurfaceBoundaryHandling_T = FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >;
// function describing the initialization velocity profile (in global cell coordinates)
inline Vector3< real_t > velocityProfile(Cell globalCell, real_t timePeriod, uint_t timestep,
const Vector3< real_t >& domainSize)
{
// add 0.5 to get the cell's center
const real_t x = (real_c(globalCell.x()) + real_c(0.5)) / domainSize[0];
const real_t y = (real_c(globalCell.y()) + real_c(0.5)) / domainSize[1];
const real_t z = (real_c(globalCell.z()) + real_c(0.5)) / domainSize[2];
const real_t timeTerm = real_c(std::cos(math::pi * real_t(timestep) / timePeriod));
const real_t sinpix = real_c(std::sin(math::pi * x));
const real_t sinpiy = real_c(std::sin(math::pi * y));
const real_t sinpiz = real_c(std::sin(math::pi * z));
const real_t sin2pix = real_c(std::sin(real_c(2) * math::pi * x));
const real_t sin2piy = real_c(std::sin(real_c(2) * math::pi * y));
const real_t sin2piz = real_c(std::sin(real_c(2) * math::pi * z));
const real_t velocityX = real_c(2) * sinpix * sinpix * sin2piy * sin2piz * timeTerm;
const real_t velocityY = -sin2pix * sinpiy * sinpiy * sin2piz * timeTerm;
const real_t velocityZ = -sin2pix * sin2piy * sinpiz * sinpiz * timeTerm;
return Vector3< real_t >(velocityX, velocityY, velocityZ);
}
int main(int argc, char** argv)
{
Environment walberlaEnv(argc, argv);
if (argc < 2) { WALBERLA_ABORT("Please specify a parameter file as input argument.") }
// print content of parameter file
WALBERLA_LOG_INFO_ON_ROOT(*walberlaEnv.config());
// get block forest parameters from parameter file
auto blockForestParameters = walberlaEnv.config()->getOneBlock("BlockForestParameters");
const Vector3< uint_t > cellsPerBlock = blockForestParameters.getParameter< Vector3< uint_t > >("cellsPerBlock");
const Vector3< bool > periodicity = blockForestParameters.getParameter< Vector3< bool > >("periodicity");
// get domain parameters from parameter file
auto domainParameters = walberlaEnv.config()->getOneBlock("DomainParameters");
const uint_t domainWidth = domainParameters.getParameter< uint_t >("domainWidth");
const real_t bubbleDiameter = real_c(domainWidth) * real_c(0.3);
const Vector3< real_t > bubbleCenter = domainWidth * Vector3< real_t >(real_c(0.35), real_c(0.35), real_c(0.35));
// define domain size
Vector3< uint_t > domainSize;
domainSize[0] = domainWidth;
domainSize[1] = domainWidth;
domainSize[2] = domainWidth;
// compute number of blocks as defined by domainSize and cellsPerBlock
Vector3< uint_t > numBlocks;
numBlocks[0] = uint_c(std::ceil(real_c(domainSize[0]) / real_c(cellsPerBlock[0])));
numBlocks[1] = uint_c(std::ceil(real_c(domainSize[1]) / real_c(cellsPerBlock[1])));
numBlocks[2] = uint_c(std::ceil(real_c(domainSize[2]) / real_c(cellsPerBlock[2])));
// get number of (MPI) processes
const uint_t numProcesses = uint_c(MPIManager::instance()->numProcesses());
WALBERLA_CHECK_LESS_EQUAL(numProcesses, numBlocks[0] * numBlocks[1] * numBlocks[2],
"The number of MPI processes is greater than the number of blocks as defined by "
"\"domainSize/cellsPerBlock\". This would result in unused MPI processes. Either decrease "
"the number of MPI processes or increase \"cellsPerBlock\".")
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numProcesses);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellsPerBlock);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(domainSize);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numBlocks);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(domainWidth);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(periodicity);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(bubbleDiameter);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(bubbleCenter);
// get physics parameters from parameter file
auto physicsParameters = walberlaEnv.config()->getOneBlock("PhysicsParameters");
const uint_t timesteps = physicsParameters.getParameter< uint_t >("timesteps");
const uint_t timestepsToInitialPosition = physicsParameters.getParameter< uint_t >("timestepsToInitialPosition");
// compute CFL number
const real_t dx_SI = real_c(1) / real_c(domainWidth);
const real_t dt_SI = real_c(3) / real_c(timestepsToInitialPosition);
const real_t CFL = dt_SI / dx_SI; // with velocity_SI = 1
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(CFL);
// dummy collision model (LBM not simulated in this benchmark)
const CollisionModel_T collisionModel = CollisionModel_T(real_c(1));
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(timestepsToInitialPosition);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(timesteps);
// read model parameters from parameter file
const auto modelParameters = walberlaEnv.config()->getOneBlock("ModelParameters");
const std::string pdfReconstructionModel = modelParameters.getParameter< std::string >("pdfReconstructionModel");
const std::string excessMassDistributionModel =
modelParameters.getParameter< std::string >("excessMassDistributionModel");
const std::string curvatureModel = modelParameters.getParameter< std::string >("curvatureModel");
const bool useSimpleMassExchange = modelParameters.getParameter< bool >("useSimpleMassExchange");
const real_t cellConversionThreshold = modelParameters.getParameter< real_t >("cellConversionThreshold");
const real_t cellConversionForceThreshold = modelParameters.getParameter< real_t >("cellConversionForceThreshold");
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfReconstructionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionThreshold);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionForceThreshold);
// read evaluation parameters from parameter file
const auto evaluationParameters = walberlaEnv.config()->getOneBlock("EvaluationParameters");
const uint_t performanceLogFrequency = evaluationParameters.getParameter< uint_t >("performanceLogFrequency");
const uint_t evaluationFrequency = evaluationParameters.getParameter< uint_t >("evaluationFrequency");
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(performanceLogFrequency);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(evaluationFrequency);
// create non-uniform block forest (non-uniformity required for load balancing)
const std::shared_ptr< StructuredBlockForest > blockForest =
createNonUniformBlockForest(domainSize, cellsPerBlock, numBlocks, periodicity);
// create lattice model
const LatticeModel_T latticeModel = LatticeModel_T(collisionModel);
// add pdf field
const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel, field::fzyx);
// add fill level field (initialized with 1, i.e., liquid everywhere)
const BlockDataID fillFieldID =
field::addToStorage< ScalarField_T >(blockForest, "Fill level field", real_c(1.0), field::fzyx, uint_c(2));
// add boundary handling
const std::shared_ptr< FreeSurfaceBoundaryHandling_T > freeSurfaceBoundaryHandling =
std::make_shared< FreeSurfaceBoundaryHandling_T >(blockForest, pdfFieldID, fillFieldID);
const BlockDataID flagFieldID = freeSurfaceBoundaryHandling->getFlagFieldID();
const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo();
// initialize the velocity profile
for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
{
PdfField_T* const pdfField = blockIt->getData< PdfField_T >(pdfFieldID);
FlagField_T* const flagField = blockIt->getData< FlagField_T >(flagFieldID);
WALBERLA_FOR_ALL_CELLS(pdfFieldIt, pdfField, flagFieldIt, flagField, {
// cell in block-local coordinates
const Cell localCell = pdfFieldIt.cell();
// get cell in global coordinates
Cell globalCell = pdfFieldIt.cell();
blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt, localCell);
// set velocity profile (CFL_SI = CFL_LBM = velocity_LBM * dt_LBM / dx_LBM = velocity_LBM)
const Vector3< real_t > initialVelocity =
CFL * velocityProfile(globalCell, real_c(timestepsToInitialPosition), uint_c(0), domainSize);
pdfField->setDensityAndVelocity(localCell, initialVelocity, real_c(1));
}) // WALBERLA_FOR_ALL_CELLS
}
// create the spherical bubble
const geometry::Sphere sphereBubble(bubbleCenter, bubbleDiameter * real_c(0.5));
bubble_model::addBodyToFillLevelField< geometry::Sphere >(*blockForest, fillFieldID, sphereBubble, true);
// initialize domain boundary conditions from config file
const auto boundaryParameters = walberlaEnv.config()->getOneBlock("BoundaryParameters");
freeSurfaceBoundaryHandling->initFromConfig(boundaryParameters);
// IMPORTANT REMARK: this must be called only after every solid flag has been set; otherwise, the boundary handling
// might not detect solid flags correctly
freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
// communication after initialization
Communication_T communication(blockForest, flagFieldID, fillFieldID);
communication();
PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
pdfCommunication();
const ConstBlockDataID initialFillFieldID =
field::addCloneToStorage< ScalarField_T >(blockForest, fillFieldID, "Initial fill level field");
// add bubble model
const std::shared_ptr< bubble_model::BubbleModelBase > bubbleModel =
std::make_shared< bubble_model::BubbleModelConstantPressure >(real_c(1));
// create timeloop
SweepTimeloop timeloop(blockForest, timesteps);
// add surface geometry handler
const SurfaceGeometryHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > geometryHandler(
blockForest, freeSurfaceBoundaryHandling, fillFieldID, curvatureModel, false, false, real_c(0));
geometryHandler.addSweeps(timeloop);
// get fields created by surface geometry handler
const ConstBlockDataID normalFieldID = geometryHandler.getConstNormalFieldID();
// add boundary handling for standard boundaries and free surface boundaries
const AdvectionDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
blockForest, pdfFieldID, flagFieldID, fillFieldID, normalFieldID, freeSurfaceBoundaryHandling, bubbleModel,
pdfReconstructionModel, excessMassDistributionModel, useSimpleMassExchange, cellConversionThreshold,
cellConversionForceThreshold);
dynamicsHandler.addSweeps(timeloop);
// add evaluator for geometrical
const std::shared_ptr< real_t > geometricalError = std::make_shared< real_t >(real_c(0));
const GeometricalErrorEvaluator< FreeSurfaceBoundaryHandling_T, FlagField_T, ScalarField_T >
geometricalErrorEvaluator(blockForest, freeSurfaceBoundaryHandling, initialFillFieldID, fillFieldID,
evaluationFrequency, geometricalError);
timeloop.addFuncAfterTimeStep(geometricalErrorEvaluator, "Evaluator: geometrical error");
// add evaluator for total and excessive mass (mass that is currently undistributed)
const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
const std::shared_ptr< real_t > excessMass = std::make_shared< real_t >(real_c(0));
const TotalMassComputer< FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T > totalMassComputer(
blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, dynamicsHandler.getConstExcessMassFieldID(),
evaluationFrequency, totalMass, excessMass);
timeloop.addFuncAfterTimeStep(totalMassComputer, "Evaluator: total mass");
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, BlockDataID(),
geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
geometryHandler.getObstNormalFieldID());
// add triangle mesh output of free surface
SurfaceMeshWriter< ScalarField_T, FlagField_T > surfaceMeshWriter(
blockForest, fillFieldID, flagFieldID, flagIDs::liquidInterfaceGasFlagIDs, real_c(0), walberlaEnv.config());
surfaceMeshWriter(); // write initial mesh
timeloop.addFuncAfterTimeStep(surfaceMeshWriter, "Writer: surface mesh");
// add logging for computational performance
const lbm::PerformanceLogger< FlagField_T > performanceLogger(
blockForest, flagFieldID, flagIDs::liquidInterfaceFlagIDs, performanceLogFrequency);
timeloop.addFuncAfterTimeStep(performanceLogger, "Evaluator: performance logging");
WcTimingPool timingPool;
for (uint_t t = uint_c(0); t != timesteps; ++t)
{
timeloop.singleStep(timingPool, true);
if (t % evaluationFrequency == uint_c(0))
{
WALBERLA_LOG_DEVEL_ON_ROOT("time step = " << t << "\n\t\ttotal mass = " << *totalMass << "\n\t\texcess mass = "
<< *excessMass << "\n\t\tgeometrical error = " << *geometricalError);
}
// set the constant velocity profile
for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
{
PdfField_T* const pdfField = blockIt->getData< PdfField_T >(pdfFieldID);
FlagField_T* const flagField = blockIt->getData< FlagField_T >(flagFieldID);
WALBERLA_FOR_ALL_CELLS(pdfFieldIt, pdfField, flagFieldIt, flagField, {
const Cell localCell = pdfFieldIt.cell();
// get cell in global coordinates
Cell globalCell = pdfFieldIt.cell();
blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt, localCell);
// set velocity profile (CFL_SI = CFL_LBM = velocity_LBM * dt_LBM / dx_LBM = velocity_LBM)
const Vector3< real_t > velocity =
CFL * velocityProfile(globalCell, real_c(timestepsToInitialPosition), t, domainSize);
pdfField->setDensityAndVelocity(localCell, velocity, real_c(1));
}) // WALBERLA_FOR_ALL_CELLS
}
pdfCommunication();
if (t % performanceLogFrequency == uint_c(0) && t > uint_c(0)) { timingPool.logResultOnRoot(); }
}
return EXIT_SUCCESS;
}
} // namespace DeformationField
} // namespace free_surface
} // namespace walberla
int main(int argc, char** argv) { return walberla::free_surface::DeformationField::main(argc, argv); }
\ No newline at end of file
apps/benchmarks/FreeSurfaceAdvection/DeformationField.prm 0 → 100644
View file @ 2dee476c
BlockForestParameters
{
cellsPerBlock < 25, 25, 25 >;
periodicity < 1, 1, 1 >;
}
DomainParameters
{
domainWidth 50;
}
PhysicsParameters
{
timestepsToInitialPosition 3000;
timesteps 3001;
}
ModelParameters
{
pdfReconstructionModel OnlyMissing;
excessMassDistributionModel EvenlyAllInterface;
curvatureModel FiniteDifferenceMethod;
useSimpleMassExchange false;
cellConversionThreshold 1e-2;
cellConversionForceThreshold 1e-1;
}
EvaluationParameters
{
evaluationFrequency 300;
performanceLogFrequency 10000;
}
BoundaryParameters
{
// X
//Border { direction W; walldistance -1; FreeSlip{} }
//Border { direction E; walldistance -1; FreeSlip{} }
// Y
//Border { direction N; walldistance -1; FreeSlip{} }
//Border { direction S; walldistance -1; FreeSlip{} }
// Z
//Border { direction T; walldistance -1; FreeSlip{} }
//Border { direction B; walldistance -1; FreeSlip{} }
}
MeshOutputParameters
{
writeFrequency 300;
baseFolder mesh-out;
}
VTK
{
fluid_field
{
writeFrequency 300;
ghostLayers 0;
baseFolder vtk-out;
samplingResolution 1;
writers
{
fill_level;
mapped_flag;
velocity;
density;
//curvature;
//normal;
//obstacle_normal;
//pdf;
//flag;
}
inclusion_filters
{
//liquidInterfaceFilter; // only include liquid and interface cells in VTK output
}
before_functions
{
//ghost_layer_synchronization; // only needed if writing the ghost layer
gas_cell_zero_setter; // sets velocity=0 and density=1 all gas cells before writing VTK output
}
}
domain_decomposition
{
writeFrequency 0;
baseFolder vtk-out;
outputDomainDecomposition true;
}
}
\ No newline at end of file
apps/benchmarks/FreeSurfaceAdvection/SingleVortex.cpp 0 → 100644
View file @ 2dee476c
//======================================================================================================================
//
// 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 SingleVortex.cpp
//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
//
// This benchmark simulates the advection of a spherical bubble in a vortex field.
// The vortex field changes periodically so that the bubble returns to its initial position, where it should take its
// initial form. The relative geometrical error of the bubble's shape after one period is evaluated. There is no LBM
// flow simulation performed here, it is a test case for the FSLBM's mass advection. This benchmark is based on the
// two-dimensional test case from Rider and Kothe (doi: 10.1006/jcph.1998.5906). The extension to three-dimensional
// space is based on the Viktor Haag's master thesis
// (https://www10.cs.fau.de/publications/theses/2017/Haag_MT_2017.pdf).
//======================================================================================================================
#include "core/Environment.h"
#include "field/Gather.h"
#include "lbm/PerformanceLogger.h"
#include "lbm/field/AddToStorage.h"
#include "lbm/free_surface/LoadBalancing.h"
#include "lbm/free_surface/SurfaceMeshWriter.h"
#include "lbm/free_surface/TotalMassComputer.h"
#include "lbm/free_surface/VtkWriter.h"
#include "lbm/free_surface/bubble_model/Geometry.h"
#include "lbm/free_surface/surface_geometry/SurfaceGeometryHandler.h"
#include "lbm/free_surface/surface_geometry/Utility.h"
#include "lbm/lattice_model/D3Q19.h"
#include "functionality/AdvectionDynamicsHandler.h"
#include "functionality/GeometricalErrorEvaluator.h"
namespace walberla
{
namespace free_surface
{
namespace SingleVortex
{
using ScalarField_T = GhostLayerField< real_t, 1 >;
using VectorField_T = GhostLayerField< Vector3< real_t >, 1 >;
// Lattice model is only created for dummy purposes; no LBM simulation is performed
using CollisionModel_T = lbm::collision_model::SRT;
using LatticeModel_T = lbm::D3Q19< CollisionModel_T, true >;
using LatticeModelStencil_T = LatticeModel_T::Stencil;
using PdfField_T = lbm::PdfField< LatticeModel_T >;
using PdfCommunication_T = blockforest::SimpleCommunication< LatticeModelStencil_T >;
// the geometry computations in SurfaceGeometryHandler require meaningful values in the ghost layers in corner
// directions (flag field and fill level field); this holds, even if the lattice model uses a D3Q19 stencil
using CommunicationStencil_T =
typename std::conditional< LatticeModel_T::Stencil::D == uint_t(2), stencil::D2Q9, stencil::D3Q27 >::type;
using Communication_T = blockforest::SimpleCommunication< CommunicationStencil_T >;
using flag_t = uint32_t;
using FlagField_T = FlagField< flag_t >;
using FreeSurfaceBoundaryHandling_T = FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >;
// function describing the initialization velocity profile (in global cell coordinates)
inline Vector3< real_t > velocityProfile(Cell globalCell, real_t timePeriod, uint_t timestep,
const Vector3< real_t >& domainSize)
{
// add 0.5 to get the cell's center
const real_t x = (real_c(globalCell.x()) + real_c(0.5)) / domainSize[0];
const real_t y = (real_c(globalCell.y()) + real_c(0.5)) / domainSize[1];
const real_t xToDomainCenter = x - real_c(0.5);
const real_t yToDomainCenter = y - real_c(0.5);
const real_t r = real_c(std::sqrt(xToDomainCenter * xToDomainCenter + yToDomainCenter * yToDomainCenter));
const real_t tmp = real_c(1) - real_c(2) * r;
const real_t rTerm = tmp * tmp;
const real_t timeTerm = real_c(std::cos(math::pi * real_t(timestep) / timePeriod));
const real_t sinpix = real_c(std::sin(math::pi * x));
const real_t sinpiy = real_c(std::sin(math::pi * y));
const real_t velocityX = real_c(std::sin(real_c(2) * math::pi * y)) * sinpix * sinpix * timeTerm;
const real_t velocityY = -real_c(std::sin(real_c(2) * math::pi * x)) * sinpiy * sinpiy * timeTerm;
const real_t velocityZ = rTerm * timeTerm;
return Vector3< real_t >(velocityX, velocityY, velocityZ);
}
int main(int argc, char** argv)
{
Environment walberlaEnv(argc, argv);
if (argc < 2) { WALBERLA_ABORT("Please specify a parameter file as input argument.") }
// print content of parameter file
WALBERLA_LOG_INFO_ON_ROOT(*walberlaEnv.config());
// get block forest parameters from parameter file
auto blockForestParameters = walberlaEnv.config()->getOneBlock("BlockForestParameters");
const Vector3< uint_t > cellsPerBlock = blockForestParameters.getParameter< Vector3< uint_t > >("cellsPerBlock");
const Vector3< bool > periodicity = blockForestParameters.getParameter< Vector3< bool > >("periodicity");
// get domain parameters from parameter file
auto domainParameters = walberlaEnv.config()->getOneBlock("DomainParameters");
const uint_t domainWidth = domainParameters.getParameter< uint_t >("domainWidth");
const real_t bubbleDiameter = real_c(domainWidth) * real_c(0.3);
const Vector3< real_t > bubbleCenter = domainWidth * Vector3< real_t >(real_c(0.5), real_c(0.75), real_c(0.25));
// define domain size
Vector3< uint_t > domainSize;
domainSize[0] = domainWidth;
domainSize[1] = domainWidth;
domainSize[2] = domainWidth;
// compute number of blocks as defined by domainSize and cellsPerBlock
Vector3< uint_t > numBlocks;
numBlocks[0] = uint_c(std::ceil(real_c(domainSize[0]) / real_c(cellsPerBlock[0])));
numBlocks[1] = uint_c(std::ceil(real_c(domainSize[1]) / real_c(cellsPerBlock[1])));
numBlocks[2] = uint_c(std::ceil(real_c(domainSize[2]) / real_c(cellsPerBlock[2])));
// get number of (MPI) processes
const uint_t numProcesses = uint_c(MPIManager::instance()->numProcesses());
WALBERLA_CHECK_LESS_EQUAL(numProcesses, numBlocks[0] * numBlocks[1] * numBlocks[2],
"The number of MPI processes is greater than the number of blocks as defined by "
"\"domainSize/cellsPerBlock\". This would result in unused MPI processes. Either decrease "
"the number of MPI processes or increase \"cellsPerBlock\".")
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numProcesses);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellsPerBlock);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(domainSize);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numBlocks);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(domainWidth);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(periodicity);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(bubbleDiameter);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(bubbleCenter);
// get physics parameters from parameter file
auto physicsParameters = walberlaEnv.config()->getOneBlock("PhysicsParameters");
const uint_t timesteps = physicsParameters.getParameter< uint_t >("timesteps");
const uint_t timestepsToInitialPosition = physicsParameters.getParameter< uint_t >("timestepsToInitialPosition");
// compute CFL number
const real_t dx_SI = real_c(1) / real_c(domainWidth);
const real_t dt_SI = real_c(3) / real_c(timestepsToInitialPosition);
const real_t CFL = dt_SI / dx_SI; // with velocity_SI = 1
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(CFL);
// dummy collision model (LBM not simulated in this benchmark)
const CollisionModel_T collisionModel = CollisionModel_T(real_c(1));
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(timestepsToInitialPosition);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(timesteps);
// read model parameters from parameter file
const auto modelParameters = walberlaEnv.config()->getOneBlock("ModelParameters");
const std::string pdfReconstructionModel = modelParameters.getParameter< std::string >("pdfReconstructionModel");
const std::string excessMassDistributionModel =
modelParameters.getParameter< std::string >("excessMassDistributionModel");
const std::string curvatureModel = modelParameters.getParameter< std::string >("curvatureModel");
const bool useSimpleMassExchange = modelParameters.getParameter< bool >("useSimpleMassExchange");
const real_t cellConversionThreshold = modelParameters.getParameter< real_t >("cellConversionThreshold");
const real_t cellConversionForceThreshold = modelParameters.getParameter< real_t >("cellConversionForceThreshold");
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfReconstructionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionThreshold);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionForceThreshold);
// read evaluation parameters from parameter file
const auto evaluationParameters = walberlaEnv.config()->getOneBlock("EvaluationParameters");
const uint_t performanceLogFrequency = evaluationParameters.getParameter< uint_t >("performanceLogFrequency");
const uint_t evaluationFrequency = evaluationParameters.getParameter< uint_t >("evaluationFrequency");
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(performanceLogFrequency);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(evaluationFrequency);
// create non-uniform block forest (non-uniformity required for load balancing)
const std::shared_ptr< StructuredBlockForest > blockForest =
createNonUniformBlockForest(domainSize, cellsPerBlock, numBlocks, periodicity);
// create lattice model
const LatticeModel_T latticeModel = LatticeModel_T(collisionModel);
// add pdf field
const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel, field::fzyx);
// add fill level field (initialized with 1, i.e., liquid everywhere)
const BlockDataID fillFieldID =
field::addToStorage< ScalarField_T >(blockForest, "Fill level field", real_c(1.0), field::fzyx, uint_c(2));
// add boundary handling
const std::shared_ptr< FreeSurfaceBoundaryHandling_T > freeSurfaceBoundaryHandling =
std::make_shared< FreeSurfaceBoundaryHandling_T >(blockForest, pdfFieldID, fillFieldID);
const BlockDataID flagFieldID = freeSurfaceBoundaryHandling->getFlagFieldID();
const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo();
// initialize the velocity profile
for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
{
PdfField_T* const pdfField = blockIt->getData< PdfField_T >(pdfFieldID);
FlagField_T* const flagField = blockIt->getData< FlagField_T >(flagFieldID);
WALBERLA_FOR_ALL_CELLS(pdfFieldIt, pdfField, flagFieldIt, flagField, {
// cell in block-local coordinates
const Cell localCell = pdfFieldIt.cell();
// get cell in global coordinates
Cell globalCell = pdfFieldIt.cell();
blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt, localCell);
// set velocity profile (CFL_SI = CFL_LBM = velocity_LBM * dt_LBM / dx_LBM = velocity_LBM)
const Vector3< real_t > initialVelocity =
CFL * velocityProfile(globalCell, real_c(timestepsToInitialPosition), uint_c(0), domainSize);
pdfField->setDensityAndVelocity(localCell, initialVelocity, real_c(1));
}) // WALBERLA_FOR_ALL_CELLS
}
// create the spherical bubble
const geometry::Sphere sphereBubble(bubbleCenter, bubbleDiameter * real_c(0.5));
bubble_model::addBodyToFillLevelField< geometry::Sphere >(*blockForest, fillFieldID, sphereBubble, true);
// initialize domain boundary conditions from config file
const auto boundaryParameters = walberlaEnv.config()->getOneBlock("BoundaryParameters");
freeSurfaceBoundaryHandling->initFromConfig(boundaryParameters);
// IMPORTANT REMARK: this must be called only after every solid flag has been set; otherwise, the boundary handling
// might not detect solid flags correctly
freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
// communication after initialization
Communication_T communication(blockForest, flagFieldID, fillFieldID);
communication();
PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
pdfCommunication();
const ConstBlockDataID initialFillFieldID =
field::addCloneToStorage< ScalarField_T >(blockForest, fillFieldID, "Initial fill level field");
// add bubble model
const std::shared_ptr< bubble_model::BubbleModelBase > bubbleModel =
std::make_shared< bubble_model::BubbleModelConstantPressure >(real_c(1));
// create timeloop
SweepTimeloop timeloop(blockForest, timesteps);
// add surface geometry handler
const SurfaceGeometryHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > geometryHandler(
blockForest, freeSurfaceBoundaryHandling, fillFieldID, curvatureModel, false, false, real_c(0));
geometryHandler.addSweeps(timeloop);
// get fields created by surface geometry handler
const ConstBlockDataID normalFieldID = geometryHandler.getConstNormalFieldID();
// add boundary handling for standard boundaries and free surface boundaries
const AdvectionDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
blockForest, pdfFieldID, flagFieldID, fillFieldID, normalFieldID, freeSurfaceBoundaryHandling, bubbleModel,
pdfReconstructionModel, excessMassDistributionModel, useSimpleMassExchange, cellConversionThreshold,
cellConversionForceThreshold);
dynamicsHandler.addSweeps(timeloop);
// add evaluator for geometrical
const std::shared_ptr< real_t > geometricalError = std::make_shared< real_t >(real_c(0));
const GeometricalErrorEvaluator< FreeSurfaceBoundaryHandling_T, FlagField_T, ScalarField_T >
geometricalErrorEvaluator(blockForest, freeSurfaceBoundaryHandling, initialFillFieldID, fillFieldID,
evaluationFrequency, geometricalError);
timeloop.addFuncAfterTimeStep(geometricalErrorEvaluator, "Evaluator: geometrical error");
// add evaluator for total and excessive mass (mass that is currently undistributed)
const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
const std::shared_ptr< real_t > excessMass = std::make_shared< real_t >(real_c(0));
const TotalMassComputer< FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T > totalMassComputer(
blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, dynamicsHandler.getConstExcessMassFieldID(),
evaluationFrequency, totalMass, excessMass);
timeloop.addFuncAfterTimeStep(totalMassComputer, "Evaluator: total mass");
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, BlockDataID(),
geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
geometryHandler.getObstNormalFieldID());
// add triangle mesh output of free surface
SurfaceMeshWriter< ScalarField_T, FlagField_T > surfaceMeshWriter(
blockForest, fillFieldID, flagFieldID, flagIDs::liquidInterfaceGasFlagIDs, real_c(0), walberlaEnv.config());
surfaceMeshWriter(); // write initial mesh
timeloop.addFuncAfterTimeStep(surfaceMeshWriter, "Writer: surface mesh");
// add logging for computational performance
const lbm::PerformanceLogger< FlagField_T > performanceLogger(
blockForest, flagFieldID, flagIDs::liquidInterfaceFlagIDs, performanceLogFrequency);
timeloop.addFuncAfterTimeStep(performanceLogger, "Evaluator: performance logging");
WcTimingPool timingPool;
for (uint_t t = uint_c(0); t != timesteps; ++t)
{
timeloop.singleStep(timingPool, true);
if (t % evaluationFrequency == uint_c(0))
{
WALBERLA_LOG_DEVEL_ON_ROOT("time step = " << t << "\n\t\ttotal mass = " << *totalMass << "\n\t\texcess mass = "
<< *excessMass << "\n\t\tgeometrical error = " << *geometricalError);
}
// set the constant velocity profile
for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
{
PdfField_T* const pdfField = blockIt->getData< PdfField_T >(pdfFieldID);
FlagField_T* const flagField = blockIt->getData< FlagField_T >(flagFieldID);
WALBERLA_FOR_ALL_CELLS(pdfFieldIt, pdfField, flagFieldIt, flagField, {
const Cell localCell = pdfFieldIt.cell();
// get cell in global coordinates
Cell globalCell = pdfFieldIt.cell();
blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt, localCell);
// set velocity profile (CFL_SI = CFL_LBM = velocity_LBM * dt_LBM / dx_LBM = velocity_LBM)
const Vector3< real_t > velocity =
CFL * velocityProfile(globalCell, real_c(timestepsToInitialPosition), t, domainSize);
pdfField->setDensityAndVelocity(localCell, velocity, real_c(1));
}) // WALBERLA_FOR_ALL_CELLS
}
pdfCommunication();
if (t % performanceLogFrequency == uint_c(0) && t > uint_c(0)) { timingPool.logResultOnRoot(); }
}
return EXIT_SUCCESS;
}
} // namespace SingleVortex
} // namespace free_surface
} // namespace walberla
int main(int argc, char** argv) { return walberla::free_surface::SingleVortex::main(argc, argv); }
\ No newline at end of file
apps/benchmarks/FreeSurfaceAdvection/SingleVortex.prm 0 → 100644
View file @ 2dee476c
BlockForestParameters
{
cellsPerBlock < 25, 25, 25 >;
periodicity < 1, 1, 1 >;
}
DomainParameters
{
domainWidth 50;
}
PhysicsParameters
{
timestepsToInitialPosition 3000;
timesteps 3001;
}
ModelParameters
{
pdfReconstructionModel OnlyMissing;
excessMassDistributionModel EvenlyAllInterface;
curvatureModel FiniteDifferenceMethod;
useSimpleMassExchange false;
cellConversionThreshold 1e-2;
cellConversionForceThreshold 1e-1;
}
EvaluationParameters
{
evaluationFrequency 300;
performanceLogFrequency 10000;
}
BoundaryParameters
{
// X
//Border { direction W; walldistance -1; FreeSlip{} }
//Border { direction E; walldistance -1; FreeSlip{} }
// Y
//Border { direction N; walldistance -1; FreeSlip{} }
//Border { direction S; walldistance -1; FreeSlip{} }
// Z
//Border { direction T; walldistance -1; FreeSlip{} }
//Border { direction B; walldistance -1; FreeSlip{} }
}
MeshOutputParameters
{
writeFrequency 300;
baseFolder mesh-out;
}
VTK
{
fluid_field
{
writeFrequency 300;
ghostLayers 0;
baseFolder vtk-out;
samplingResolution 1;
writers
{
fill_level;
mapped_flag;
velocity;
density;
//curvature;
//normal;
//obstacle_normal;
//pdf;
//flag;
}
inclusion_filters
{
//liquidInterfaceFilter; // only include liquid and interface cells in VTK output
}
before_functions
{
//ghost_layer_synchronization; // only needed if writing the ghost layer
gas_cell_zero_setter; // sets velocity=0 and density=1 all gas cells before writing VTK output
}
}
domain_decomposition
{
writeFrequency 0;
baseFolder vtk-out;
outputDomainDecomposition true;
}
}
\ No newline at end of file
apps/benchmarks/FreeSurfaceAdvection/ZalesakDisk.cpp 0 → 100644
View file @ 2dee476c
//======================================================================================================================
//
// 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 ZalesakDisk.cpp
//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
//
// This benchmark simulates the advection of a slotted disk of gas in a constant rotating velocity field. The disk
// returns to its initial position, where it should take its initial form. The relative geometrical error of the
// bubble's shape after one rotation is evaluated. There is no LBM flow simulation performed here, it is a test case for
// the FSLBM's mass advection. This benchmark is commonly referred to as Zalesak's rotating disk (see
// doi: 10.1016/0021-9991(79)90051-2). The setup chosen here is identical to the one used by Janssen (see
// doi: 10.1016/j.camwa.2009.08.064).
//======================================================================================================================
#include "core/Environment.h"
#include "field/Gather.h"
#include "lbm/PerformanceLogger.h"
#include "lbm/field/AddToStorage.h"
#include "lbm/free_surface/LoadBalancing.h"
#include "lbm/free_surface/SurfaceMeshWriter.h"
#include "lbm/free_surface/TotalMassComputer.h"
#include "lbm/free_surface/VtkWriter.h"
#include "lbm/free_surface/bubble_model/Geometry.h"
#include "lbm/free_surface/surface_geometry/SurfaceGeometryHandler.h"
#include "lbm/free_surface/surface_geometry/Utility.h"
#include "lbm/lattice_model/D2Q9.h"
#include "functionality/AdvectionDynamicsHandler.h"
#include "functionality/GeometricalErrorEvaluator.h"
namespace walberla
{
namespace free_surface
{
namespace ZalesakDisk
{
using ScalarField_T = GhostLayerField< real_t, 1 >;
using VectorField_T = GhostLayerField< Vector3< real_t >, 1 >;
// Lattice model is only created for dummy purposes; no LBM simulation is performed
using CollisionModel_T = lbm::collision_model::SRT;
using LatticeModel_T = lbm::D2Q9< CollisionModel_T, true >;
using LatticeModelStencil_T = LatticeModel_T::Stencil;
using PdfField_T = lbm::PdfField< LatticeModel_T >;
using PdfCommunication_T = blockforest::SimpleCommunication< LatticeModelStencil_T >;
// the geometry computations in SurfaceGeometryHandler require meaningful values in the ghost layers in corner
// directions (flag field and fill level field); this holds, even if the lattice model uses a D3Q19 stencil
using CommunicationStencil_T =
typename std::conditional< LatticeModel_T::Stencil::D == uint_t(2), stencil::D2Q9, stencil::D3Q27 >::type;
using Communication_T = blockforest::SimpleCommunication< CommunicationStencil_T >;
using flag_t = uint32_t;
using FlagField_T = FlagField< flag_t >;
using FreeSurfaceBoundaryHandling_T = FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >;
// function describing the initialization velocity profile (in global cell coordinates)
inline Vector3< real_t > velocityProfile(real_t angularVelocity, Cell globalCell, const Vector3< real_t >& domainCenter)
{
// add 0.5 to get Cell's center
const real_t velocityX = -angularVelocity * ((real_c(globalCell.y()) + real_c(0.5)) - domainCenter[1]);
const real_t velocityY = angularVelocity * ((real_c(globalCell.x()) + real_c(0.5)) - domainCenter[0]);
return Vector3< real_t >(velocityX, velocityY, real_c(0));
}
int main(int argc, char** argv)
{
Environment walberlaEnv(argc, argv);
if (argc < 2) { WALBERLA_ABORT("Please specify a parameter file as input argument.") }
// print content of parameter file
WALBERLA_LOG_INFO_ON_ROOT(*walberlaEnv.config());
// get block forest parameters from parameter file
auto blockForestParameters = walberlaEnv.config()->getOneBlock("BlockForestParameters");
const Vector3< uint_t > cellsPerBlock = blockForestParameters.getParameter< Vector3< uint_t > >("cellsPerBlock");
const Vector3< bool > periodicity = blockForestParameters.getParameter< Vector3< bool > >("periodicity");
// get domain parameters from parameter file
auto domainParameters = walberlaEnv.config()->getOneBlock("DomainParameters");
const uint_t domainWidth = domainParameters.getParameter< uint_t >("domainWidth");
const real_t diskRadius = real_c(domainWidth) * real_c(0.125);
const Vector3< real_t > diskCenter =
Vector3< real_t >(real_c(domainWidth) * real_c(0.5), real_c(domainWidth) * real_c(0.75), real_c(0.5));
const real_t diskSlotLength = real_c(2) * diskRadius - real_c(0.1) * real_c(domainWidth);
const real_t diskSlotWidth = real_c(0.06) * real_c(domainWidth);
// define domain size
Vector3< uint_t > domainSize;
domainSize[0] = domainWidth;
domainSize[1] = domainWidth;
domainSize[2] = uint_c(1);
const Vector3< real_t > domainCenter =
real_c(0.5) * Vector3< real_t >(real_c(domainSize[0]), real_c(domainSize[1]), real_c(domainSize[2]));
// compute number of blocks as defined by domainSize and cellsPerBlock
Vector3< uint_t > numBlocks;
numBlocks[0] = uint_c(std::ceil(real_c(domainSize[0]) / real_c(cellsPerBlock[0])));
numBlocks[1] = uint_c(std::ceil(real_c(domainSize[1]) / real_c(cellsPerBlock[1])));
numBlocks[2] = uint_c(std::ceil(real_c(domainSize[2]) / real_c(cellsPerBlock[2])));
// get number of (MPI) processes
const uint_t numProcesses = uint_c(MPIManager::instance()->numProcesses());
WALBERLA_CHECK_LESS_EQUAL(numProcesses, numBlocks[0] * numBlocks[1] * numBlocks[2],
"The number of MPI processes is greater than the number of blocks as defined by "
"\"domainSize/cellsPerBlock\". This would result in unused MPI processes. Either decrease "
"the number of MPI processes or increase \"cellsPerBlock\".")
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numProcesses);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellsPerBlock);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(domainSize);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numBlocks);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(domainWidth);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(periodicity);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(diskRadius);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(diskCenter);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(diskSlotLength);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(diskSlotWidth);
// get physics parameters from parameter file
auto physicsParameters = walberlaEnv.config()->getOneBlock("PhysicsParameters");
const uint_t timesteps = physicsParameters.getParameter< uint_t >("timesteps");
const uint_t timestepsFullRotation = physicsParameters.getParameter< uint_t >("timestepsFullRotation");
// compute CFL number
const real_t dx_SI = real_c(4) / real_c(domainWidth);
const real_t dt_SI = real_c(12.59652) / real_c(timestepsFullRotation);
const real_t CFL = dt_SI / dx_SI; // with velocity_SI = 1
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(CFL);
// dummy collision model (LBM not simulated in this benchmark)
const CollisionModel_T collisionModel = CollisionModel_T(real_c(1));
const real_t angularVelocity = real_c(2) * math::pi / real_c(timestepsFullRotation);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(timestepsFullRotation);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(timesteps);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(angularVelocity);
// read model parameters from parameter file
const auto modelParameters = walberlaEnv.config()->getOneBlock("ModelParameters");
const std::string pdfReconstructionModel = modelParameters.getParameter< std::string >("pdfReconstructionModel");
const std::string excessMassDistributionModel =
modelParameters.getParameter< std::string >("excessMassDistributionModel");
const std::string curvatureModel = modelParameters.getParameter< std::string >("curvatureModel");
const bool useSimpleMassExchange = modelParameters.getParameter< bool >("useSimpleMassExchange");
const real_t cellConversionThreshold = modelParameters.getParameter< real_t >("cellConversionThreshold");
const real_t cellConversionForceThreshold = modelParameters.getParameter< real_t >("cellConversionForceThreshold");
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfReconstructionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionThreshold);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionForceThreshold);
// read evaluation parameters from parameter file
const auto evaluationParameters = walberlaEnv.config()->getOneBlock("EvaluationParameters");
const uint_t performanceLogFrequency = evaluationParameters.getParameter< uint_t >("performanceLogFrequency");
const uint_t evaluationFrequency = evaluationParameters.getParameter< uint_t >("evaluationFrequency");
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(performanceLogFrequency);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(evaluationFrequency);
// create non-uniform block forest (non-uniformity required for load balancing)
const std::shared_ptr< StructuredBlockForest > blockForest =
createNonUniformBlockForest(domainSize, cellsPerBlock, numBlocks, periodicity);
// create lattice model
const LatticeModel_T latticeModel = LatticeModel_T(collisionModel);
// add pdf field
const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel, field::fzyx);
// add fill level field (initialized with 1, i.e., liquid everywhere)
const BlockDataID fillFieldID =
field::addToStorage< ScalarField_T >(blockForest, "Fill level field", real_c(1.0), field::fzyx, uint_c(2));
// add boundary handling
const std::shared_ptr< FreeSurfaceBoundaryHandling_T > freeSurfaceBoundaryHandling =
std::make_shared< FreeSurfaceBoundaryHandling_T >(blockForest, pdfFieldID, fillFieldID);
const BlockDataID flagFieldID = freeSurfaceBoundaryHandling->getFlagFieldID();
const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo();
// initialize the velocity profile
for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
{
PdfField_T* const pdfField = blockIt->getData< PdfField_T >(pdfFieldID);
WALBERLA_FOR_ALL_CELLS(pdfFieldIt, pdfField, {
// cell in block-local coordinates
const Cell localCell = pdfFieldIt.cell();
// get cell in global coordinates
Cell globalCell = pdfFieldIt.cell();
blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt, localCell);
// set velocity profile
const Vector3< real_t > initialVelocity = velocityProfile(angularVelocity, globalCell, domainCenter);
pdfField->setDensityAndVelocity(localCell, initialVelocity, real_c(1));
}) // WALBERLA_FOR_ALL_CELLS
}
// create the disk
const Vector3< real_t > diskBottomEnd = diskCenter - Vector3< real_t >(real_c(0), real_c(0), -real_c(domainSize[2]));
const Vector3< real_t > diskTopEnd = diskCenter - Vector3< real_t >(real_c(0), real_c(0), real_c(domainSize[2]));
const geometry::Cylinder disk(diskBottomEnd, diskTopEnd, diskRadius);
bubble_model::addBodyToFillLevelField< geometry::Cylinder >(*blockForest, fillFieldID, disk, true);
// create the disk's slot
const Vector3< real_t > slotMinCorner = Vector3< real_t >(diskCenter[0] - diskSlotWidth * real_c(0.5),
diskCenter[1] - diskRadius, -real_c(domainSize[2]));
const Vector3< real_t > slotMaxCorner = Vector3< real_t >(
diskCenter[0] + diskSlotWidth * real_c(0.5), diskCenter[1] - diskRadius + diskSlotLength, real_c(domainSize[2]));
AABB slotAABB(slotMinCorner, slotMaxCorner);
bubble_model::addBodyToFillLevelField< AABB >(*blockForest, fillFieldID, slotAABB, false);
// initialize domain boundary conditions from config file
const auto boundaryParameters = walberlaEnv.config()->getOneBlock("BoundaryParameters");
freeSurfaceBoundaryHandling->initFromConfig(boundaryParameters);
// IMPORTANT REMARK: this must be called only after every solid flag has been set; otherwise, the boundary handling
// might not detect solid flags correctly
freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
// communication after initialization
Communication_T communication(blockForest, flagFieldID, fillFieldID);
communication();
PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
pdfCommunication();
const ConstBlockDataID initialFillFieldID =
field::addCloneToStorage< ScalarField_T >(blockForest, fillFieldID, "Initial fill level field");
// add bubble model
const std::shared_ptr< bubble_model::BubbleModelBase > bubbleModel =
std::make_shared< bubble_model::BubbleModelConstantPressure >(real_c(1));
// create timeloop
SweepTimeloop timeloop(blockForest, timesteps);
// add surface geometry handler
const SurfaceGeometryHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > geometryHandler(
blockForest, freeSurfaceBoundaryHandling, fillFieldID, curvatureModel, false, false, real_c(0));
geometryHandler.addSweeps(timeloop);
// get fields created by surface geometry handler
const ConstBlockDataID normalFieldID = geometryHandler.getConstNormalFieldID();
// add boundary handling for standard boundaries and free surface boundaries
const AdvectionDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
blockForest, pdfFieldID, flagFieldID, fillFieldID, normalFieldID, freeSurfaceBoundaryHandling, bubbleModel,
pdfReconstructionModel, excessMassDistributionModel, useSimpleMassExchange, cellConversionThreshold,
cellConversionForceThreshold);
dynamicsHandler.addSweeps(timeloop);
// add evaluator for geometrical
const std::shared_ptr< real_t > geometricalError = std::make_shared< real_t >(real_c(0));
const GeometricalErrorEvaluator< FreeSurfaceBoundaryHandling_T, FlagField_T, ScalarField_T >
geometricalErrorEvaluator(blockForest, freeSurfaceBoundaryHandling, initialFillFieldID, fillFieldID,
evaluationFrequency, geometricalError);
timeloop.addFuncAfterTimeStep(geometricalErrorEvaluator, "Evaluator: geometrical errors");
// add evaluator for total and excessive mass (mass that is currently undistributed)
const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
const std::shared_ptr< real_t > excessMass = std::make_shared< real_t >(real_c(0));
const TotalMassComputer< FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T > totalMassComputer(
blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, dynamicsHandler.getConstExcessMassFieldID(),
evaluationFrequency, totalMass, excessMass);
timeloop.addFuncAfterTimeStep(totalMassComputer, "Evaluator: total mass");
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, BlockDataID(),
geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
geometryHandler.getObstNormalFieldID());
// add triangle mesh output of free surface
SurfaceMeshWriter< ScalarField_T, FlagField_T > surfaceMeshWriter(
blockForest, fillFieldID, flagFieldID, flagIDs::liquidInterfaceGasFlagIDs, real_c(0), walberlaEnv.config());
surfaceMeshWriter(); // write initial mesh
timeloop.addFuncAfterTimeStep(surfaceMeshWriter, "Writer: surface mesh");
// add logging for computational performance
const lbm::PerformanceLogger< FlagField_T > performanceLogger(
blockForest, flagFieldID, flagIDs::liquidInterfaceFlagIDs, performanceLogFrequency);
timeloop.addFuncAfterTimeStep(performanceLogger, "Evaluator: performance logging");
WcTimingPool timingPool;
for (uint_t t = uint_c(0); t != timesteps; ++t)
{
timeloop.singleStep(timingPool, true);
if (t % evaluationFrequency == uint_c(0))
{
WALBERLA_LOG_DEVEL_ON_ROOT("time step = " << t << "\n\t\ttotal mass = " << *totalMass << "\n\t\texcess mass = "
<< *excessMass << "\n\t\tgeometrical error = " << *geometricalError);
}
// set the constant velocity profile
for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
{
PdfField_T* const pdfField = blockIt->getData< PdfField_T >(pdfFieldID);
FlagField_T* const flagField = blockIt->getData< FlagField_T >(flagFieldID);
WALBERLA_FOR_ALL_CELLS(pdfFieldIt, pdfField, flagFieldIt, flagField, {
const Cell localCell = pdfFieldIt.cell();
// get cell in global coordinates
Cell globalCell = pdfFieldIt.cell();
blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt, localCell);
// set velocity profile
const Vector3< real_t > initialVelocity = velocityProfile(angularVelocity, globalCell, domainCenter);
pdfField->setDensityAndVelocity(localCell, initialVelocity, real_c(1));
}) // WALBERLA_FOR_ALL_CELLS
}
pdfCommunication();
if (t % performanceLogFrequency == uint_c(0) && t > uint_c(0)) { timingPool.logResultOnRoot(); }
}
return EXIT_SUCCESS;
}
} // namespace ZalesakDisk
} // namespace free_surface
} // namespace walberla
int main(int argc, char** argv) { return walberla::free_surface::ZalesakDisk::main(argc, argv); }
\ No newline at end of file
apps/benchmarks/FreeSurfaceAdvection/ZalesakDisk.prm 0 → 100644
View file @ 2dee476c
BlockForestParameters
{
cellsPerBlock < 100, 100, 1 >;
periodicity < 0, 0, 0 >;
}
DomainParameters
{
domainWidth 200;
}
PhysicsParameters
{
timestepsFullRotation 12570; // angularVelocity = 2 * pi / timestepsFullRotation
timesteps 12571;
}
ModelParameters
{
pdfReconstructionModel OnlyMissing;
excessMassDistributionModel EvenlyAllInterface;
curvatureModel FiniteDifferenceMethod;
useSimpleMassExchange false;
cellConversionThreshold 1e-2;
cellConversionForceThreshold 1e-1;
}
EvaluationParameters
{
evaluationFrequency 12570;
performanceLogFrequency 10000;
}
BoundaryParameters
{
// X
//Border { direction W; walldistance -1; FreeSlip{} }
//Border { direction E; walldistance -1; FreeSlip{} }
// Y
//Border { direction N; walldistance -1; FreeSlip{} }
//Border { direction S; walldistance -1; FreeSlip{} }
// Z
//Border { direction T; walldistance -1; FreeSlip{} }
//Border { direction B; walldistance -1; FreeSlip{} }
}
MeshOutputParameters
{
writeFrequency 0;
baseFolder mesh-out;
}
VTK
{
fluid_field
{
writeFrequency 12570;
ghostLayers 0;
baseFolder vtk-out;
samplingResolution 1;
writers
{
fill_level;
mapped_flag;
velocity;
density;
//curvature;
//normal;
//obstacle_normal;
//pdf;
//flag;
}
inclusion_filters
{
//liquidInterfaceFilter; // only include liquid and interface cells in VTK output
}
before_functions
{
//ghost_layer_synchronization; // only needed if writing the ghost layer
gas_cell_zero_setter; // sets velocity=0 and density=1 all gas cells before writing VTK output
}
}
domain_decomposition
{
writeFrequency 0;
baseFolder vtk-out;
outputDomainDecomposition true;
}
}
\ No newline at end of file
apps/benchmarks/FreeSurfaceAdvection/functionality/AdvectSweep.h 0 → 100644
View file @ 2dee476c
//======================================================================================================================
//
// 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 AdvectSweep.h
//! \ingroup free_surface
//! \author Martin Bauer
//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
//! \brief Sweep for mass advection and interface cell conversion marking (simplified StreamReconstructAdvectSweep).
//
//======================================================================================================================
#pragma once
#include "core/DataTypes.h"
#include "core/math/Vector3.h"
#include "core/mpi/Reduce.h"
#include "core/timing/TimingPool.h"
#include "field/FieldClone.h"
#include "field/FlagField.h"
#include "lbm/free_surface/FlagInfo.h"
#include "lbm/free_surface/bubble_model/BubbleModel.h"
#include "lbm/free_surface/dynamics/PdfReconstructionModel.h"
#include "lbm/free_surface/dynamics/functionality/AdvectMass.h"
#include "lbm/free_surface/dynamics/functionality/FindInterfaceCellConversion.h"
#include "lbm/free_surface/dynamics/functionality/GetOredNeighborhood.h"
#include "lbm/free_surface/dynamics/functionality/ReconstructInterfaceCellABB.h"
#include "lbm/sweeps/StreamPull.h"
namespace walberla
{
namespace free_surface
{
template< typename LatticeModel_T, typename BoundaryHandling_T, typename FlagField_T, typename FlagInfo_T,
typename ScalarField_T, typename VectorField_T >
class AdvectSweep
{
public:
using flag_t = typename FlagInfo_T::flag_t;
using PdfField_T = lbm::PdfField< LatticeModel_T >;
AdvectSweep(BlockDataID handlingID, BlockDataID fillFieldID, BlockDataID flagFieldID, BlockDataID pdfField,
const FlagInfo_T& flagInfo, const PdfReconstructionModel& pdfReconstructionModel,
bool useSimpleMassExchange, real_t cellConversionThreshold, real_t cellConversionForceThreshold)
: handlingID_(handlingID), fillFieldID_(fillFieldID), flagFieldID_(flagFieldID), pdfFieldID_(pdfField),
flagInfo_(flagInfo), neighborhoodFlagFieldClone_(flagFieldID), fillFieldClone_(fillFieldID),
pdfFieldClone_(pdfField), pdfReconstructionModel_(pdfReconstructionModel),
useSimpleMassExchange_(useSimpleMassExchange), cellConversionThreshold_(cellConversionThreshold),
cellConversionForceThreshold_(cellConversionForceThreshold)
{}
void operator()(IBlock* const block);
protected:
BlockDataID handlingID_;
BlockDataID fillFieldID_;
BlockDataID flagFieldID_;
BlockDataID pdfFieldID_;
FlagInfo_T flagInfo_;
// efficient clones of fields to provide temporary fields (for writing to)
field::FieldClone< FlagField_T, true > neighborhoodFlagFieldClone_;
field::FieldClone< ScalarField_T, true > fillFieldClone_;
field::FieldClone< PdfField_T, true > pdfFieldClone_;
PdfReconstructionModel pdfReconstructionModel_;
bool useSimpleMassExchange_;
real_t cellConversionThreshold_;
real_t cellConversionForceThreshold_;
}; // class AdvectSweep
template< typename LatticeModel_T, typename BoundaryHandling_T, typename FlagField_T, typename FlagInfo_T,
typename ScalarField_T, typename VectorField_T >
void AdvectSweep< LatticeModel_T, BoundaryHandling_T, FlagField_T, FlagInfo_T, ScalarField_T,
VectorField_T >::operator()(IBlock* const block)
{
// fetch data
ScalarField_T* const fillSrcField = block->getData< ScalarField_T >(fillFieldID_);
PdfField_T* const pdfSrcField = block->getData< PdfField_T >(pdfFieldID_);
FlagField_T* const flagField = block->getData< FlagField_T >(flagFieldID_);
// temporary fields that act as destination fields
FlagField_T* const neighborhoodFlagField = neighborhoodFlagFieldClone_.get(block);
ScalarField_T* const fillDstField = fillFieldClone_.get(block);
// combine all neighbor flags using bitwise OR and write them to the neighborhood field
// this is simply a pre-computation of often required values
// IMPORTANT REMARK: the "OredNeighborhood" is also required for each cell in the first ghost layer; this requires
// access to all first ghost layer cell's neighbors (i.e., to the second ghost layer)
WALBERLA_CHECK_GREATER_EQUAL(flagField->nrOfGhostLayers(), uint_c(2));
getOredNeighborhood< typename LatticeModel_T::Stencil >(flagField, neighborhoodFlagField);
// explicitly avoid OpenMP due to bubble model update (reportFillLevelChange)
WALBERLA_FOR_ALL_CELLS_OMP(
pdfSrcFieldIt, pdfSrcField, fillSrcFieldIt, fillSrcField, fillDstFieldIt, fillDstField, flagFieldIt, flagField,
neighborhoodFlagFieldIt, neighborhoodFlagField, omp critical, {
if (flagInfo_.isInterface(flagFieldIt))
{
const real_t rho = lbm::getDensity(pdfSrcField->latticeModel(), pdfSrcFieldIt);
// compute mass advection using post-collision PDFs (explicitly not PDFs updated by stream above)
const real_t deltaMass =
(advectMass< LatticeModel_T, FlagField_T, typename ScalarField_T::iterator,
typename PdfField_T::iterator, typename FlagField_T::iterator,
typename FlagField_T::iterator, FlagInfo_T >) (flagField, fillSrcFieldIt, pdfSrcFieldIt,
flagFieldIt, neighborhoodFlagFieldIt,
flagInfo_, useSimpleMassExchange_);
// update fill level after LBM stream and mass exchange
*fillDstFieldIt = *fillSrcFieldIt + deltaMass / rho;
}
else // treat non-interface cells
{
// manually adjust the fill level to avoid outdated fill levels being copied from fillSrcField
if (flagInfo_.isGas(flagFieldIt)) { *fillDstFieldIt = real_c(0); }
else
{
if (flagInfo_.isLiquid(flagFieldIt)) { *fillDstFieldIt = real_c(1); }
else // flag is e.g. obstacle or outflow
{
*fillDstFieldIt = *fillSrcFieldIt;
}
}
}
}) // WALBERLA_FOR_ALL_CELLS_XYZ_OMP
fillSrcField->swapDataPointers(fillDstField);
BoundaryHandling_T* const handling = block->getData< BoundaryHandling_T >(handlingID_);
// mark interface cell for conversion
findInterfaceCellConversions< LatticeModel_T >(handling, fillSrcField, flagField, neighborhoodFlagField, flagInfo_,
cellConversionThreshold_, cellConversionForceThreshold_);
}
} // namespace free_surface
} // namespace walberla
apps/benchmarks/FreeSurfaceAdvection/functionality/AdvectionDynamicsHandler.h 0 → 100644
View file @ 2dee476c
//======================================================================================================================
//
// 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 AdvectionDynamicsHandler.h
//! \ingroup surface_dynamics
//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
//! \brief Handles free surface advection (without LBM flow simulation; this is a simplified SurfaceDynamicsHandler).
//
//======================================================================================================================
#pragma once
#include "core/DataTypes.h"
#include "domain_decomposition/StructuredBlockStorage.h"
#include "field/AddToStorage.h"
#include "field/FlagField.h"
#include "lbm/blockforest/communication/SimpleCommunication.h"
#include "lbm/blockforest/communication/UpdateSecondGhostLayer.h"
#include "lbm/free_surface/BlockStateDetectorSweep.h"
#include "lbm/free_surface/FlagInfo.h"
#include "lbm/free_surface/boundary/FreeSurfaceBoundaryHandling.h"
#include "lbm/free_surface/bubble_model/BubbleModel.h"
#include "lbm/free_surface/dynamics/CellConversionSweep.h"
#include "lbm/free_surface/dynamics/ConversionFlagsResetSweep.h"
#include "lbm/free_surface/dynamics/ExcessMassDistributionModel.h"
#include "lbm/free_surface/dynamics/ExcessMassDistributionSweep.h"
#include "lbm/free_surface/dynamics/PdfReconstructionModel.h"
#include "lbm/sweeps/CellwiseSweep.h"
#include "lbm/sweeps/SweepWrappers.h"
#include "stencil/D3Q27.h"
#include "timeloop/SweepTimeloop.h"
#include "AdvectSweep.h"
namespace walberla
{
namespace free_surface
{
template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename VectorField_T >
class AdvectionDynamicsHandler
{
protected:
using Communication_T = blockforest::SimpleCommunication< typename LatticeModel_T::Stencil >;
// communication in corner directions (D2Q9/D3Q27) is required for all fields but the PDF field
using CommunicationStencil_T =
typename std::conditional< LatticeModel_T::Stencil::D == uint_t(2), stencil::D2Q9, stencil::D3Q27 >::type;
using CommunicationCorner_T = blockforest::SimpleCommunication< CommunicationStencil_T >;
using FreeSurfaceBoundaryHandling_T = FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >;
public:
AdvectionDynamicsHandler(const std::shared_ptr< StructuredBlockForest >& blockForest, BlockDataID pdfFieldID,
BlockDataID flagFieldID, BlockDataID fillFieldID, ConstBlockDataID normalFieldID,
const std::shared_ptr< FreeSurfaceBoundaryHandling_T >& freeSurfaceBoundaryHandling,
const std::shared_ptr< BubbleModelBase >& bubbleModel,
const std::string& pdfReconstructionModel, const std::string& excessMassDistributionModel,
bool useSimpleMassExchange, real_t cellConversionThreshold,
real_t cellConversionForceThreshold)
: blockForest_(blockForest), pdfFieldID_(pdfFieldID), flagFieldID_(flagFieldID), fillFieldID_(fillFieldID),
normalFieldID_(normalFieldID), bubbleModel_(bubbleModel),
freeSurfaceBoundaryHandling_(freeSurfaceBoundaryHandling), pdfReconstructionModel_(pdfReconstructionModel),
excessMassDistributionModel_({ excessMassDistributionModel }), useSimpleMassExchange_(useSimpleMassExchange),
cellConversionThreshold_(cellConversionThreshold), cellConversionForceThreshold_(cellConversionForceThreshold)
{
WALBERLA_CHECK(LatticeModel_T::compressible,
"The free surface lattice Boltzmann extension works only with compressible LBM models.");
if (excessMassDistributionModel_.isEvenlyAllInterfaceFallbackLiquidType())
{
// add additional field for storing excess mass in liquid cells
excessMassFieldID_ =
field::addToStorage< ScalarField_T >(blockForest_, "Excess mass", real_c(0), field::fzyx, uint_c(1));
}
}
ConstBlockDataID getConstExcessMassFieldID() const { return excessMassFieldID_; }
/********************************************************************************************************************
* The order of these sweeps is similar to page 40 in the dissertation of T. Pohl, 2008.
*******************************************************************************************************************/
void addSweeps(SweepTimeloop& timeloop) const
{
using StateSweep = BlockStateDetectorSweep< FlagField_T >;
const auto& flagInfo = freeSurfaceBoundaryHandling_->getFlagInfo();
const auto blockStateUpdate = StateSweep(blockForest_, flagInfo, flagFieldID_);
// empty sweeps required for using selectors (e.g. StateSweep::onlyGasAndBoundary)
const auto emptySweep = [](IBlock*) {};
// sweep for
// - reconstruction of PDFs in interface cells
// - streaming of PDFs in interface cells (and liquid cells on the same block)
// - advection of mass
// - update bubble volumes
// - marking interface cells for conversion
const AdvectSweep< LatticeModel_T, typename FreeSurfaceBoundaryHandling_T::BoundaryHandling_T, FlagField_T,
typename FreeSurfaceBoundaryHandling_T::FlagInfo_T, ScalarField_T, VectorField_T >
advectSweep(freeSurfaceBoundaryHandling_->getHandlingID(), fillFieldID_, flagFieldID_, pdfFieldID_, flagInfo,
pdfReconstructionModel_, useSimpleMassExchange_, cellConversionThreshold_,
cellConversionForceThreshold_);
// sweep acts only on blocks with at least one interface cell (due to StateSweep::fullFreeSurface)
timeloop.add()
<< Sweep(advectSweep, "Sweep: Advect", StateSweep::fullFreeSurface)
<< Sweep(emptySweep, "Empty sweep: Advect")
// do not communicate PDFs here:
// - stream on blocks with "StateSweep::fullFreeSurface" was performed here using post-collision PDFs
// - stream on other blocks is performed below and should also use post-collision PDFs
// => if PDFs were communicated here, the ghost layer of other blocks would have post-stream PDFs
<< AfterFunction(CommunicationCorner_T(blockForest_, fillFieldID_, flagFieldID_),
"Communication: after Advect sweep")
<< AfterFunction(blockforest::UpdateSecondGhostLayer< ScalarField_T >(blockForest_, fillFieldID_),
"Second ghost layer update: after Advect sweep (fill level field)")
<< AfterFunction(blockforest::UpdateSecondGhostLayer< FlagField_T >(blockForest_, flagFieldID_),
"Second ghost layer update: after Advect sweep (flag field)");
// convert cells
// - according to the flags from StreamReconstructAdvectSweep (interface -> gas/liquid)
// - to ensure a closed layer of interface cells (gas/liquid -> interface)
// - detect and register bubble merges/splits (bubble volumes are already updated in StreamReconstructAdvectSweep)
// - convert cells and initialize PDFs near inflow boundaries
const CellConversionSweep< LatticeModel_T, typename FreeSurfaceBoundaryHandling_T::BoundaryHandling_T,
ScalarField_T >
cellConvSweep(freeSurfaceBoundaryHandling_->getHandlingID(), pdfFieldID_, flagInfo, bubbleModel_.get());
timeloop.add() << Sweep(cellConvSweep, "Sweep: cell conversion", StateSweep::fullFreeSurface)
<< Sweep(emptySweep, "Empty sweep: cell conversion")
<< AfterFunction(Communication_T(blockForest_, pdfFieldID_),
"Communication: after cell conversion sweep (PDF field)")
// communicate the flag field also in corner directions
<< AfterFunction(CommunicationCorner_T(blockForest_, flagFieldID_),
"Communication: after cell conversion sweep (flag field)")
<< AfterFunction(blockforest::UpdateSecondGhostLayer< FlagField_T >(blockForest_, flagFieldID_),
"Second ghost layer update: after cell conversion sweep (flag field)");
// distribute excess mass:
// - excess mass: mass that is free after conversion from interface to gas/liquid cells
// - update the bubble model
if (excessMassDistributionModel_.isEvenlyType())
{
const ExcessMassDistributionSweepInterfaceEvenly< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T >
distributeMassSweep(excessMassDistributionModel_, fillFieldID_, flagFieldID_, pdfFieldID_, flagInfo);
timeloop.add() << Sweep(distributeMassSweep, "Sweep: excess mass distribution", StateSweep::fullFreeSurface)
<< Sweep(emptySweep, "Empty sweep: distribute excess mass")
<< AfterFunction(CommunicationCorner_T(blockForest_, fillFieldID_),
"Communication: after excess mass distribution sweep")
<< AfterFunction(
blockforest::UpdateSecondGhostLayer< ScalarField_T >(blockForest_, fillFieldID_),
"Second ghost layer update: after excess mass distribution sweep (fill level field)")
// update bubble model, i.e., perform registered bubble merges/splits; bubble merges/splits are
// already detected and registered by CellConversionSweep
<< AfterFunction(std::bind(&bubble_model::BubbleModelBase::update, bubbleModel_),
"Sweep: bubble model update");
}
else
{
if (excessMassDistributionModel_.isWeightedType())
{
const ExcessMassDistributionSweepInterfaceWeighted< LatticeModel_T, FlagField_T, ScalarField_T,
VectorField_T >
distributeMassSweep(excessMassDistributionModel_, fillFieldID_, flagFieldID_, pdfFieldID_, flagInfo,
normalFieldID_);
timeloop.add() << Sweep(distributeMassSweep, "Sweep: excess mass distribution", StateSweep::fullFreeSurface)
<< Sweep(emptySweep, "Empty sweep: distribute excess mass")
<< AfterFunction(CommunicationCorner_T(blockForest_, fillFieldID_),
"Communication: after excess mass distribution sweep")
<< AfterFunction(
blockforest::UpdateSecondGhostLayer< ScalarField_T >(blockForest_, fillFieldID_),
"Second ghost layer update: after excess mass distribution sweep (fill level field)")
// update bubble model, i.e., perform registered bubble merges/splits; bubble merges/splits
// are already detected and registered by CellConversionSweep
<< AfterFunction(std::bind(&bubble_model::BubbleModelBase::update, bubbleModel_),
"Sweep: bubble model update");
}
else
{
if (excessMassDistributionModel_.isEvenlyAllInterfaceFallbackLiquidType())
{
const ExcessMassDistributionSweepInterfaceAndLiquid< LatticeModel_T, FlagField_T, ScalarField_T,
VectorField_T >
distributeMassSweep(excessMassDistributionModel_, fillFieldID_, flagFieldID_, pdfFieldID_, flagInfo,
excessMassFieldID_);
timeloop.add()
// perform this sweep also on "onlyLBM" blocks because liquid cells also exchange excess mass here
<< Sweep(distributeMassSweep, "Sweep: excess mass distribution", StateSweep::fullFreeSurface)
<< Sweep(distributeMassSweep, "Sweep: excess mass distribution", StateSweep::onlyLBM)
<< Sweep(emptySweep, "Empty sweep: distribute excess mass")
<< AfterFunction(CommunicationCorner_T(blockForest_, fillFieldID_, excessMassFieldID_),
"Communication: after excess mass distribution sweep")
<< AfterFunction(blockforest::UpdateSecondGhostLayer< ScalarField_T >(blockForest_, fillFieldID_),
"Second ghost layer update: after excess mass distribution sweep (fill level field)")
// update bubble model, i.e., perform registered bubble merges/splits; bubble
// merges/splits are already detected and registered by CellConversionSweep
<< AfterFunction(std::bind(&bubble_model::BubbleModelBase::update, bubbleModel_),
"Sweep: bubble model update");
}
}
}
// reset all flags that signal cell conversions (except "keepInterfaceForWettingFlag")
ConversionFlagsResetSweep< FlagField_T > resetConversionFlagsSweep(flagFieldID_, flagInfo);
timeloop.add() << Sweep(resetConversionFlagsSweep, "Sweep: conversion flag reset", StateSweep::fullFreeSurface)
<< Sweep(emptySweep, "Empty sweep: conversion flag reset")
<< AfterFunction(CommunicationCorner_T(blockForest_, flagFieldID_),
"Communication: after excess mass distribution sweep")
<< AfterFunction(blockforest::UpdateSecondGhostLayer< FlagField_T >(blockForest_, flagFieldID_),
"Second ghost layer update: after excess mass distribution sweep (flag field)");
// update block states
timeloop.add() << Sweep(blockStateUpdate, "Sweep: block state update");
}
private:
std::shared_ptr< StructuredBlockForest > blockForest_;
BlockDataID pdfFieldID_;
BlockDataID flagFieldID_;
BlockDataID fillFieldID_;
ConstBlockDataID normalFieldID_;
std::shared_ptr< BubbleModelBase > bubbleModel_;
std::shared_ptr< FreeSurfaceBoundaryHandling_T > freeSurfaceBoundaryHandling_;
PdfReconstructionModel pdfReconstructionModel_;
ExcessMassDistributionModel excessMassDistributionModel_;
bool useSimpleMassExchange_;
real_t cellConversionThreshold_;
real_t cellConversionForceThreshold_;
BlockDataID excessMassFieldID_ = BlockDataID();
}; // class AdvectionDynamicsHandler
} // namespace free_surface
} // namespace walberla
\ No newline at end of file
apps/benchmarks/FreeSurfaceAdvection/functionality/GeometricalErrorEvaluator.h 0 → 100644
View file @ 2dee476c
//======================================================================================================================
//
// 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 GeometricalErrorEvaluator.h
//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
//! \brief Compute the relative geometrical error in free-surface advection test cases.
//======================================================================================================================
#include "blockforest/StructuredBlockForest.h"
#include "core/DataTypes.h"
#include "domain_decomposition/BlockDataID.h"
#include "field/iterators/IteratorMacros.h"
namespace walberla
{
namespace free_surface
{
template< typename FreeSurfaceBoundaryHandling_T, typename FlagField_T, typename ScalarField_T >
class GeometricalErrorEvaluator
{
public:
GeometricalErrorEvaluator(const std::weak_ptr< StructuredBlockForest >& blockForest,
const std::weak_ptr< const FreeSurfaceBoundaryHandling_T >& freeSurfaceBoundaryHandling,
const ConstBlockDataID& initialfillFieldID, const ConstBlockDataID& fillFieldID,
uint_t frequency, const std::shared_ptr< real_t >& geometricalError)
: blockForest_(blockForest), freeSurfaceBoundaryHandling_(freeSurfaceBoundaryHandling),
initialFillFieldID_(initialfillFieldID), fillFieldID_(fillFieldID), frequency_(frequency),
geometricalError_(geometricalError), executionCounter_(uint_c(0))
{}
void operator()()
{
if (frequency_ == uint_c(0)) { return; }
auto blockForest = blockForest_.lock();
WALBERLA_CHECK_NOT_NULLPTR(blockForest);
auto freeSurfaceBoundaryHandling = freeSurfaceBoundaryHandling_.lock();
WALBERLA_CHECK_NOT_NULLPTR(freeSurfaceBoundaryHandling);
if (executionCounter_ == uint_c(0))
{
computeInitialFillLevelSum(blockForest, freeSurfaceBoundaryHandling);
computeError(blockForest, freeSurfaceBoundaryHandling);
}
else
{
// only evaluate in given frequencies
if (executionCounter_ % frequency_ == uint_c(0)) { computeError(blockForest, freeSurfaceBoundaryHandling); }
}
++executionCounter_;
}
void computeInitialFillLevelSum(
const std::shared_ptr< const StructuredBlockForest >& blockForest,
const std::shared_ptr< const FreeSurfaceBoundaryHandling_T >& freeSurfaceBoundaryHandling)
{
const BlockDataID flagFieldID = freeSurfaceBoundaryHandling->getFlagFieldID();
const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo();
for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
{
const ScalarField_T* const initialfillField = blockIt->getData< const ScalarField_T >(initialFillFieldID_);
const FlagField_T* const flagField = blockIt->getData< const FlagField_T >(flagFieldID);
// avoid OpenMP here because initialFillLevelSum_ is a class member and not a regular variable
WALBERLA_FOR_ALL_CELLS_OMP(initialfillFieldIt, initialfillField, flagFieldIt, flagField, omp critical, {
if (flagInfo.isInterface(flagFieldIt) || flagInfo.isLiquid(flagFieldIt))
{
initialFillLevelSum_ += *initialfillFieldIt;
}
}) // WALBERLA_FOR_ALL_CELLS_OMP
}
mpi::allReduceInplace< real_t >(initialFillLevelSum_, mpi::SUM);
}
void computeError(const std::shared_ptr< const StructuredBlockForest >& blockForest,
const std::shared_ptr< const FreeSurfaceBoundaryHandling_T >& freeSurfaceBoundaryHandling)
{
const BlockDataID flagFieldID = freeSurfaceBoundaryHandling->getFlagFieldID();
const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo();
real_t geometricalError = real_c(0);
for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
{
const ScalarField_T* const initialfillField = blockIt->getData< const ScalarField_T >(initialFillFieldID_);
const ScalarField_T* const fillField = blockIt->getData< const ScalarField_T >(fillFieldID_);
const FlagField_T* const flagField = blockIt->getData< const FlagField_T >(flagFieldID);
WALBERLA_FOR_ALL_CELLS_OMP(initialfillFieldIt, initialfillField, fillFieldIt, fillField, flagFieldIt,
flagField, omp parallel for schedule(static) reduction(+:geometricalError), {
if (flagInfo.isInterface(flagFieldIt) || flagInfo.isLiquid(flagFieldIt))
{
geometricalError += real_c(std::abs(*initialfillFieldIt - *fillFieldIt));
}
}) // WALBERLA_FOR_ALL_CELLS_OMP
}
mpi::allReduceInplace< real_t >(geometricalError, mpi::SUM);
// compute L1 norms
*geometricalError_ = geometricalError / initialFillLevelSum_;
}
private:
std::weak_ptr< StructuredBlockForest > blockForest_;
std::weak_ptr< const FreeSurfaceBoundaryHandling_T > freeSurfaceBoundaryHandling_;
ConstBlockDataID initialFillFieldID_;
ConstBlockDataID fillFieldID_;
uint_t frequency_;
std::shared_ptr< real_t > geometricalError_;
uint_t executionCounter_;
real_t initialFillLevelSum_ = real_c(0);
}; // class GeometricalErrorEvaluator
} // namespace free_surface
} // namespace walberla
\ No newline at end of file