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 2427 additions and 484 deletions
apps/benchmarks/ComplexGeometry/ComplexGeometry.cpp
View file @ 7a6d3b57
...@@ -16,6 +16,7 @@ ...@@ -16,6 +16,7 @@
//! \file ComplexGeometry.cpp //! \file ComplexGeometry.cpp
//! \ingroup mesh //! \ingroup mesh
//! \author Christian Godenschwager <christian.godenschwager@fau.de> //! \author Christian Godenschwager <christian.godenschwager@fau.de>
//! \author Markus Holzer <markus.holzer@fau.de>
// //
//====================================================================================================================== //======================================================================================================================
...@@ -24,241 +25,268 @@ ...@@ -24,241 +25,268 @@
#include "blockforest/communication/UniformBufferedScheme.h" #include "blockforest/communication/UniformBufferedScheme.h"
#include "blockforest/loadbalancing/StaticParMetis.h" #include "blockforest/loadbalancing/StaticParMetis.h"
#include "core/SharedFunctor.h"
#include "core/Environment.h" #include "core/Environment.h"
#include "core/SharedFunctor.h"
#include "core/logging/Logging.h" #include "core/logging/Logging.h"
#include "core/math/IntegerFactorization.h" #include "core/math/IntegerFactorization.h"
#include "core/timing/RemainingTimeLogger.h"
#include "domain_decomposition/SharedSweep.h" #include "domain_decomposition/SharedSweep.h"
#include "field/AddToStorage.h" #include "field/AddToStorage.h"
#include "field/StabilityChecker.h" #include "field/StabilityChecker.h"
#include "geometry/InitBoundaryHandling.h"
#include "geometry/mesh/TriangleMesh.h"
#include "geometry/mesh/TriangleMeshIO.h"
#include "lbm/BlockForestEvaluation.h"
#include "lbm/PerformanceEvaluation.h"
#include "lbm/PerformanceLogger.h"
#include "lbm/boundary/factories/DefaultBoundaryHandling.h" #include "lbm/boundary/factories/DefaultBoundaryHandling.h"
#include "lbm/communication/PdfFieldPackInfo.h" #include "lbm/communication/PdfFieldPackInfo.h"
#include "lbm/communication/SparsePdfFieldPackInfo.h" #include "lbm/communication/SparsePdfFieldPackInfo.h"
#include "lbm/field/AddToStorage.h" #include "lbm/field/AddToStorage.h"
#include "lbm/lattice_model/CollisionModel.h"
#include "lbm/lattice_model/D3Q19.h" #include "lbm/lattice_model/D3Q19.h"
#include "lbm/lattice_model/D3Q27.h" #include "lbm/lattice_model/D3Q27.h"
#include "lbm/lattice_model/CollisionModel.h"
#include "lbm/lattice_model/ForceModel.h" #include "lbm/lattice_model/ForceModel.h"
#include "lbm/refinement/TimeStep.h" #include "lbm/refinement/TimeStep.h"
#include "lbm/sweeps/CellwiseSweep.h" #include "lbm/sweeps/CellwiseSweep.h"
#include "lbm/sweeps/SplitPureSweep.h" #include "lbm/sweeps/SplitPureSweep.h"
#include "lbm/vtk/VTKOutput.h" #include "lbm/vtk/VTKOutput.h"
#include "lbm/BlockForestEvaluation.h"
#include "lbm/PerformanceEvaluation.h"
#include "lbm/PerformanceLogger.h"
#include "geometry/mesh/TriangleMesh.h" #include "mesh/blockforest/BlockExclusion.h"
#include "geometry/mesh/TriangleMeshIO.h"
#include "geometry/InitBoundaryHandling.h"
#include "mesh/TriangleMeshes.h"
#include "mesh/MeshOperations.h"
#include "mesh/DistanceComputations.h"
#include "mesh/DistanceFunction.h"
#include "mesh/MeshIO.h"
#include "mesh/MatrixVectorOperations.h"
#include "mesh/blockforest/BlockForestInitialization.h" #include "mesh/blockforest/BlockForestInitialization.h"
#include "mesh/blockforest/BlockWorkloadMemory.h" #include "mesh/blockforest/BlockWorkloadMemory.h"
#include "mesh/blockforest/BlockExclusion.h"
#include "mesh/blockforest/RefinementSelection.h" #include "mesh/blockforest/RefinementSelection.h"
#include "mesh/distance_octree/DistanceOctree.h"
#include "mesh/boundary/BoundarySetup.h"
#include "mesh/boundary/BoundaryInfo.h" #include "mesh/boundary/BoundaryInfo.h"
#include "mesh/boundary/BoundaryLocation.h" #include "mesh/boundary/BoundaryLocation.h"
#include "mesh/boundary/BoundaryLocationFunction.h" #include "mesh/boundary/BoundaryLocationFunction.h"
#include "mesh/boundary/BoundarySetup.h"
#include "mesh/boundary/BoundaryUIDFaceDataSource.h" #include "mesh/boundary/BoundaryUIDFaceDataSource.h"
#include "mesh/boundary/ColorToBoundaryMapper.h" #include "mesh/boundary/ColorToBoundaryMapper.h"
#include "mesh/vtk/VTKMeshWriter.h"
#include "mesh/vtk/CommonDataSources.h"
#include "stencil/D3Q19.h" #include "stencil/D3Q19.h"
#include "timeloop/SweepTimeloop.h" #include "timeloop/SweepTimeloop.h"
#include "core/timing/RemainingTimeLogger.h"
#include <cmath> #include <cmath>
#include <vector>
#include <string> #include <string>
#include <vector>
namespace walberla { #include "mesh_common/DistanceComputations.h"
#include "mesh_common/DistanceFunction.h"
#include "mesh_common/MatrixVectorOperations.h"
#include "mesh_common/MeshIO.h"
#include "mesh_common/MeshOperations.h"
#include "mesh_common/TriangleMeshes.h"
#include "mesh_common/distance_octree/DistanceOctree.h"
#include "mesh_common/vtk/CommonDataSources.h"
#include "mesh_common/vtk/VTKMeshWriter.h"
namespace walberla
{
template< typename MeshType > template< typename MeshType >
void vertexToFaceColor( MeshType & mesh, const typename MeshType::Color & defaultColor ) void vertexToFaceColor(MeshType& mesh, const typename MeshType::Color& defaultColor)
{ {
WALBERLA_CHECK( mesh.has_vertex_colors() ); WALBERLA_CHECK(mesh.has_vertex_colors())
mesh.request_face_colors(); mesh.request_face_colors();
for( auto faceIt = mesh.faces_begin(); faceIt != mesh.faces_end(); ++faceIt ) for (auto faceIt = mesh.faces_begin(); faceIt != mesh.faces_end(); ++faceIt)
{ {
typename MeshType::Color vertexColor; typename MeshType::Color vertexColor;
bool useVertexColor = true; bool useVertexColor = true;
auto vertexIt = mesh.fv_iter( *faceIt ); auto vertexIt = mesh.fv_iter(*faceIt);
WALBERLA_ASSERT( vertexIt.is_valid() ); WALBERLA_ASSERT(vertexIt.is_valid())
vertexColor = mesh.color( *vertexIt ); vertexColor = mesh.color(*vertexIt);
++vertexIt; ++vertexIt;
while( vertexIt.is_valid() && useVertexColor ) while (vertexIt.is_valid() && useVertexColor)
{ {
if( vertexColor != mesh.color( *vertexIt ) ) if (vertexColor != mesh.color(*vertexIt)) useVertexColor = false;
useVertexColor = false;
++vertexIt; ++vertexIt;
} }
mesh.set_color( *faceIt, useVertexColor ? vertexColor : defaultColor ); mesh.set_color(*faceIt, useVertexColor ? vertexColor : defaultColor);
} }
} }
int main(int argc, char* argv[])
int main( int argc, char * argv[] )
{ {
Environment env( argc, argv ); Environment env(argc, argv);
if( !env.config() ) if (!env.config()) { WALBERLA_ABORT_NO_DEBUG_INFO("USAGE: " << argv[0] << " INPUT_FILE") }
{
WALBERLA_ABORT_NO_DEBUG_INFO( "USAGE: " << argv[0] << " INPUT_FILE" );
}
mpi::MPIManager::instance()->useWorldComm(); mpi::MPIManager::instance()->useWorldComm();
const auto & config = *( env.config() ); const auto& config = *(env.config());
Config::BlockHandle configBlock = config.getOneBlock("ComplexGeometry");
Config::BlockHandle configBlock = config.getOneBlock( "ComplexGeometry" ); const std::string meshFile = configBlock.getParameter< std::string >("meshFile");
const real_t dx = configBlock.getParameter< real_t >("coarseDx");
const real_t omega = configBlock.getParameter< real_t >("coarseOmega");
const std::string meshFile = configBlock.getParameter< std::string >( "meshFile" ); if (configBlock.getParameter< bool >("logLevelDetail"))
const real_t dx = configBlock.getParameter< real_t >( "coarseDx" ); walberla::logging::Logging::instance()->setLogLevel(walberla::logging::Logging::DETAIL);
const real_t omega = configBlock.getParameter< real_t >( "coarseOmega" );
//const uint_t blockPerProcess = configBlock.getParameter< uint_t >( "blocksPerProcess", uint_t(6) );
const uint_t timeSteps = configBlock.getParameter< uint_t >( "coarseTimeSteps" );
const Vector3<real_t> bodyForce = configBlock.getParameter< Vector3<real_t> >( "bodyForce" );
//const bool sparseCommunication = configBlock.getParameter< bool >( "sparseCommunication", true );
const Vector3<real_t> domainBlowUp = configBlock.getParameter< Vector3<real_t> >( "domainBlowUp", Vector3<real_t>(6) );
const Vector3<uint_t> blockSize = configBlock.getParameter< Vector3<uint_t> >( "blockSize", Vector3<uint_t>(16) );
uint_t numLevels = configBlock.getParameter< uint_t >( "numLevels", uint_t(2) );
numLevels = std::max( numLevels, uint_t(1) ); const uint_t timeSteps = configBlock.getParameter< uint_t >("coarseTimeSteps");
const Vector3< real_t > bodyForce = configBlock.getParameter< Vector3< real_t > >("bodyForce");
const Vector3< real_t > domainBlowUp =
configBlock.getParameter< Vector3< real_t > >("domainBlowUp", Vector3< real_t >(6));
const Vector3< uint_t > blockSize =
configBlock.getParameter< Vector3< uint_t > >("blockSize", Vector3< uint_t >(16));
uint_t numLevels = configBlock.getParameter< uint_t >("numLevels", uint_t(2));
const bool WriteDistanceOctree = configBlock.getParameter< bool >("WriteDistanceOctree", false);
const bool WriteSetupForestAndReturn = configBlock.getParameter< bool >("WriteSetupForestAndReturn", false);
//uint_t numProcesses = uint_c( MPIManager::instance()->numProcesses() ); numLevels = std::max(numLevels, uint_t(1));
const real_t fineDX = dx / real_c(std::pow(2, numLevels));
WALBERLA_LOG_DEVEL_VAR_ON_ROOT( meshFile ); // uint_t numProcesses = uint_c( MPIManager::instance()->numProcesses() );
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(meshFile)
auto mesh = make_shared< mesh::TriangleMesh >(); auto mesh = make_shared< mesh::TriangleMesh >();
mesh->request_vertex_colors(); mesh->request_vertex_colors();
WALBERLA_LOG_DEVEL_ON_ROOT( "Loading mesh" ); WALBERLA_LOG_DEVEL_ON_ROOT("Loading mesh")
mesh::readAndBroadcast( meshFile, *mesh); mesh::readAndBroadcast(meshFile, *mesh);
vertexToFaceColor( *mesh, mesh::TriangleMesh::Color(255,255,255) ); vertexToFaceColor(*mesh, mesh::TriangleMesh::Color(255, 255, 255));
WALBERLA_LOG_DEVEL_ON_ROOT("Adding distance info to mesh")
auto triDist = make_shared< mesh::TriangleDistance< mesh::TriangleMesh > >(mesh);
WALBERLA_LOG_DEVEL_ON_ROOT("Building distance octree")
auto distanceOctree = make_shared< mesh::DistanceOctree< mesh::TriangleMesh > >(triDist);
WALBERLA_LOG_DEVEL_ON_ROOT("done. Octree has height " << distanceOctree->height())
// write distance octree to file
if (WriteDistanceOctree) {
distanceOctree->writeVTKOutput("distanceOctree");
}
WALBERLA_LOG_DEVEL_ON_ROOT( "Adding distance info to mesh" ); auto aabb = computeAABB(*mesh);
auto triDist = make_shared< mesh::TriangleDistance<mesh::TriangleMesh> >( mesh ); aabb.scale(domainBlowUp);
WALBERLA_LOG_DEVEL_ON_ROOT( "Building distance octree" );
auto distanceOctree = make_shared< mesh::DistanceOctree<mesh::TriangleMesh> >( triDist );
WALBERLA_LOG_DEVEL_ON_ROOT( "done. Octree has height " << distanceOctree->height() );
distanceOctree->writeVTKOutput("distanceOctree"); mesh::ComplexGeometryStructuredBlockforestCreator bfc(aabb, Vector3< real_t >(dx));
auto aabb = computeAABB( *mesh ); bfc.setRootBlockExclusionFunction(mesh::makeExcludeMeshInterior(distanceOctree, dx));
aabb.scale( domainBlowUp ); bfc.setBlockExclusionFunction(mesh::makeExcludeMeshInteriorRefinement(distanceOctree, fineDX));
mesh::ComplexGeometryStructuredBlockforestCreator bfc( aabb, Vector3<real_t>( dx ), mesh::makeExcludeMeshInterior( distanceOctree, dx ) ); auto meshWorkloadMemory = mesh::makeMeshWorkloadMemory(distanceOctree, dx);
auto meshWorkloadMemory = mesh::makeMeshWorkloadMemory( distanceOctree, dx );
meshWorkloadMemory.setInsideCellWorkload(1); meshWorkloadMemory.setInsideCellWorkload(1);
meshWorkloadMemory.setOutsideCellWorkload(1); meshWorkloadMemory.setOutsideCellWorkload(1);
bfc.setWorkloadMemorySUIDAssignmentFunction( meshWorkloadMemory ); bfc.setWorkloadMemorySUIDAssignmentFunction(meshWorkloadMemory);
bfc.setPeriodicity( Vector3<bool>(true) ); bfc.setPeriodicity(Vector3< bool >(true));
bfc.setRefinementSelectionFunction( makeRefinementSelection( distanceOctree, numLevels - uint_t(1), dx, dx * real_t(1) ) ); bfc.setRefinementSelectionFunction(
makeRefinementSelection(distanceOctree, numLevels - uint_t(1), dx, dx * real_t(1)));
auto structuredBlockforest = bfc.createStructuredBlockForest( blockSize ); if (WriteSetupForestAndReturn)
{
WALBERLA_LOG_INFO_ON_ROOT("Setting up SetupBlockForest")
auto setupForest = bfc.createSetupBlockForest(blockSize);
WALBERLA_LOG_INFO_ON_ROOT("Writing SetupBlockForest to VTK file")
WALBERLA_ROOT_SECTION()
{
setupForest->writeVTKOutput("SetupBlockForest");
}
WALBERLA_LOG_INFO_ON_ROOT("Stopping program")
return EXIT_SUCCESS;
}
typedef lbm::D3Q19<lbm::collision_model::SRT, false, lbm::force_model::SimpleConstant> LatticeModel_T; auto structuredBlockforest = bfc.createStructuredBlockForest(blockSize);
using flag_t = walberla::uint8_t; typedef lbm::D3Q19< lbm::collision_model::SRT, false, lbm::force_model::SimpleConstant > LatticeModel_T;
using FlagField_T = FlagField<flag_t>;
using PdfField_T = lbm::PdfField<LatticeModel_T>;
LatticeModel_T latticeModel{ lbm::collision_model::SRT( omega ), lbm::force_model::SimpleConstant( bodyForce ) }; using flag_t = walberla::uint8_t;
using FlagField_T = FlagField< flag_t >;
using PdfField_T = lbm::PdfField< LatticeModel_T >;
LatticeModel_T latticeModel{ lbm::collision_model::SRT(omega), lbm::force_model::SimpleConstant(bodyForce) };
static const uint_t NUM_GHOSTLAYERS = 4; static const uint_t NUM_GHOSTLAYERS = 4;
BlockDataID pdfFieldId = lbm::addPdfFieldToStorage( structuredBlockforest, "pdf field", latticeModel, Vector3<real_t>(0), real_t(1), NUM_GHOSTLAYERS, field::fzyx ); BlockDataID pdfFieldId = lbm::addPdfFieldToStorage(structuredBlockforest, "pdf field", latticeModel,
BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >( structuredBlockforest, "flag field", NUM_GHOSTLAYERS ); Vector3< real_t >(0), real_t(1), NUM_GHOSTLAYERS, field::fzyx);
BlockDataID flagFieldId =
field::addFlagFieldToStorage< FlagField_T >(structuredBlockforest, "flag field", NUM_GHOSTLAYERS);
const FlagUID fluidFlagUID( "Fluid" ); const FlagUID fluidFlagUID("Fluid");
typedef lbm::DefaultBoundaryHandlingFactory< LatticeModel_T, FlagField_T > BHFactory; typedef lbm::DefaultBoundaryHandlingFactory< LatticeModel_T, FlagField_T > BHFactory;
auto boundariesConfig = configBlock.getOneBlock( "Boundaries" ); auto boundariesConfig = configBlock.getOneBlock("Boundaries");
BlockDataID boundaryHandlingId = BHFactory::addBoundaryHandlingToStorage( structuredBlockforest, "boundary handling", flagFieldId, pdfFieldId, fluidFlagUID, BlockDataID boundaryHandlingId = BHFactory::addBoundaryHandlingToStorage(
boundariesConfig.getParameter< Vector3<real_t> >( "velocity0", Vector3<real_t>() ), structuredBlockforest, "boundary handling", flagFieldId, pdfFieldId, fluidFlagUID,
boundariesConfig.getParameter< Vector3<real_t> >( "velocity1", Vector3<real_t>() ), boundariesConfig.getParameter< Vector3< real_t > >("velocity0", Vector3< real_t >()),
boundariesConfig.getParameter< real_t > ( "pressure0", real_c( 1.0 ) ), boundariesConfig.getParameter< Vector3< real_t > >("velocity1", Vector3< real_t >()),
boundariesConfig.getParameter< real_t > ( "pressure1", real_c( 1.001 ) ) ); boundariesConfig.getParameter< real_t >("pressure0", real_c(1.0)),
boundariesConfig.getParameter< real_t >("pressure1", real_c(1.001)));
mesh::ColorToBoundaryMapper<mesh::TriangleMesh> colorToBoundryMapper(( mesh::BoundaryInfo( BHFactory::getNoSlipBoundaryUID() ) )); mesh::ColorToBoundaryMapper< mesh::TriangleMesh > colorToBoundryMapper(
(mesh::BoundaryInfo(BHFactory::getNoSlipBoundaryUID())));
// colorToBoundryMapper.set( mesh::TriangleMesh::Color(255,0,0), mesh::BoundaryInfo( BHFactory::getPressure0BoundaryUID() ) ); // colorToBoundryMapper.set( mesh::TriangleMesh::Color(255,0,0), mesh::BoundaryInfo(
// colorToBoundryMapper.set( mesh::TriangleMesh::Color(0,0,255), mesh::BoundaryInfo( BHFactory::getPressure1BoundaryUID() ) ); // BHFactory::getPressure0BoundaryUID() ) ); colorToBoundryMapper.set( mesh::TriangleMesh::Color(0,0,255),
// colorToBoundryMapper.set( mesh::TriangleMesh::Color(255,255,255), mesh::BoundaryInfo( BHFactory::getNoSlipBoundaryUID() ) ); // mesh::BoundaryInfo( BHFactory::getPressure1BoundaryUID() ) ); colorToBoundryMapper.set(
// mesh::TriangleMesh::Color(255,255,255), mesh::BoundaryInfo( BHFactory::getNoSlipBoundaryUID() ) );
auto boundaryLocations = colorToBoundryMapper.addBoundaryInfoToMesh( *mesh ); auto boundaryLocations = colorToBoundryMapper.addBoundaryInfoToMesh(*mesh);
mesh::VTKMeshWriter< mesh::TriangleMesh > meshWriter( mesh, "meshBoundaries", 1 ); mesh::VTKMeshWriter< mesh::TriangleMesh > meshWriter(mesh, "meshBoundaries", 1);
meshWriter.addDataSource( make_shared< mesh::BoundaryUIDFaceDataSource< mesh::TriangleMesh > >( boundaryLocations ) ); meshWriter.addDataSource(make_shared< mesh::BoundaryUIDFaceDataSource< mesh::TriangleMesh > >(boundaryLocations));
meshWriter.addDataSource( make_shared< mesh::ColorFaceDataSource< mesh::TriangleMesh > >() ); meshWriter.addDataSource(make_shared< mesh::ColorFaceDataSource< mesh::TriangleMesh > >());
meshWriter.addDataSource( make_shared< mesh::ColorVertexDataSource< mesh::TriangleMesh > >() ); meshWriter.addDataSource(make_shared< mesh::ColorVertexDataSource< mesh::TriangleMesh > >());
meshWriter(); meshWriter();
WALBERLA_LOG_DEVEL_ON_ROOT( "Voxelizing mesh" ); WALBERLA_LOG_DEVEL_ON_ROOT("Voxelizing mesh")
mesh::BoundarySetup boundarySetup( structuredBlockforest, makeMeshDistanceFunction( distanceOctree ), NUM_GHOSTLAYERS ); mesh::BoundarySetup boundarySetup(structuredBlockforest, makeMeshDistanceFunction(distanceOctree), NUM_GHOSTLAYERS);
//WALBERLA_LOG_DEVEL( "Writing Voxelisation" ); // WALBERLA_LOG_DEVEL( "Writing Voxelisation" );
//boundarySetup.writeVTKVoxelfile(); // boundarySetup.writeVTKVoxelfile();
WALBERLA_LOG_DEVEL_ON_ROOT( "Setting up fluid cells" ); WALBERLA_LOG_DEVEL_ON_ROOT("Setting up fluid cells")
boundarySetup.setDomainCells<BHFactory::BoundaryHandling>( boundaryHandlingId, mesh::BoundarySetup::OUTSIDE ); boundarySetup.setDomainCells< BHFactory::BoundaryHandling >(boundaryHandlingId, mesh::BoundarySetup::OUTSIDE);
WALBERLA_LOG_DEVEL_ON_ROOT( "Setting up boundaries" ); WALBERLA_LOG_DEVEL_ON_ROOT("Setting up boundaries")
boundarySetup.setBoundaries<BHFactory::BoundaryHandling>( boundaryHandlingId, makeBoundaryLocationFunction( distanceOctree, boundaryLocations ), mesh::BoundarySetup::INSIDE ); boundarySetup.setBoundaries< BHFactory::BoundaryHandling >(
WALBERLA_LOG_DEVEL_ON_ROOT( "done" ); boundaryHandlingId, makeBoundaryLocationFunction(distanceOctree, boundaryLocations), mesh::BoundarySetup::INSIDE);
WALBERLA_LOG_DEVEL_ON_ROOT("done")
lbm::BlockForestEvaluation<FlagField_T>( structuredBlockforest, flagFieldId, fluidFlagUID ).logInfoOnRoot(); lbm::BlockForestEvaluation< FlagField_T >(structuredBlockforest, flagFieldId, fluidFlagUID).logInfoOnRoot();
lbm::PerformanceLogger<FlagField_T> perfLogger( structuredBlockforest, flagFieldId, fluidFlagUID, 100 ); lbm::PerformanceLogger< FlagField_T > perfLogger(structuredBlockforest, flagFieldId, fluidFlagUID, 100);
SweepTimeloop timeloop( structuredBlockforest->getBlockStorage(), timeSteps ); SweepTimeloop timeloop(structuredBlockforest->getBlockStorage(), timeSteps);
auto sweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldId, flagFieldId, fluidFlagUID ); auto sweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >(pdfFieldId, flagFieldId, fluidFlagUID);
auto refinementTimeStep = lbm::refinement::makeTimeStep<LatticeModel_T, BHFactory::BoundaryHandling > ( structuredBlockforest, sweep, pdfFieldId, boundaryHandlingId ); auto refinementTimeStep = lbm::refinement::makeTimeStep< LatticeModel_T, BHFactory::BoundaryHandling >(
timeloop.addFuncBeforeTimeStep( makeSharedFunctor( refinementTimeStep ), "Refinement time step" ); structuredBlockforest, sweep, pdfFieldId, boundaryHandlingId);
timeloop.addFuncBeforeTimeStep(makeSharedFunctor(refinementTimeStep), "Refinement time step");
// log remaining time // log remaining time
timeloop.addFuncAfterTimeStep( timing::RemainingTimeLogger( timeloop.getNrOfTimeSteps() ), "remaining time logger" ); timeloop.addFuncAfterTimeStep(timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps()), "remaining time logger");
// LBM stability check // LBM stability check
timeloop.addFuncAfterTimeStep( makeSharedFunctor( field::makeStabilityChecker< PdfField_T, FlagField_T >( env.config(), structuredBlockforest, pdfFieldId, timeloop.addFuncAfterTimeStep(makeSharedFunctor(field::makeStabilityChecker< PdfField_T, FlagField_T >(
flagFieldId, fluidFlagUID ) ), env.config(), structuredBlockforest, pdfFieldId, flagFieldId, fluidFlagUID)),
"LBM stability check" ); "LBM stability check");
timeloop.addFuncAfterTimeStep( perfLogger, "PerformanceLogger" );
timeloop.addFuncAfterTimeStep(perfLogger, "Evaluator: performance logging");
// add VTK output to time loop // add VTK output to time loop
lbm::VTKOutput< LatticeModel_T, FlagField_T >::addToTimeloop( timeloop, structuredBlockforest, env.config(), pdfFieldId, flagFieldId, fluidFlagUID ); lbm::VTKOutput< LatticeModel_T, FlagField_T >::addToTimeloop(timeloop, structuredBlockforest, env.config(),
pdfFieldId, flagFieldId, fluidFlagUID);
WcTimingPool timingPool; WcTimingPool timingPool;
WALBERLA_LOG_INFO_ON_ROOT( "Starting timeloop" ); WALBERLA_LOG_INFO_ON_ROOT("Starting timeloop")
timeloop.run( timingPool ); timeloop.run(timingPool);
WALBERLA_LOG_INFO_ON_ROOT( "Timeloop done" ); WALBERLA_LOG_INFO_ON_ROOT("Timeloop done")
timingPool.unifyRegisteredTimersAcrossProcesses(); timingPool.unifyRegisteredTimersAcrossProcesses();
timingPool.logResultOnRoot( timing::REDUCE_TOTAL, true ); timingPool.logResultOnRoot(timing::REDUCE_TOTAL, true);
return EXIT_SUCCESS; return EXIT_SUCCESS;
} }
} // namespace walberla } // namespace walberla
int main( int argc, char * argv[] ) int main(int argc, char* argv[]) { return walberla::main(argc, argv); }
{
return walberla::main( argc, argv );
}
apps/benchmarks/ComplexGeometry/test.conf
View file @ 7a6d3b57
ComplexGeometry ComplexGeometry
{ {
meshFile cube.obj; meshFile bunny.obj;
coarseDx 0.1; coarseDx 4;
coarseOmega 1.6; coarseOmega 1.6;
coarseTimeSteps 1; coarseTimeSteps 1;
numLevels 2; numLevels 4;
bodyForce <0.0001, 0, 0>; bodyForce <0.0001, 0, 0>;
blockSize <16,16,16>; blockSize <16,16,16>;
domainBlowUp <5,5,5>; // simulation domain is blow up factor times mesh size per dimension domainBlowUp <5,5,5>; // simulation domain is blow up factor times mesh size per dimension
WriteDistanceOctree false;
WriteSetupForestAndReturn true;
logLevelDetail false;
Boundaries { Boundaries {
} }
......
apps/benchmarks/CouetteFlow/CMakeLists.txt
View file @ 7a6d3b57
waLBerla_link_files_to_builddir( "*.dat" ) waLBerla_link_files_to_builddir( "*.dat" )
waLBerla_add_executable( NAME CouetteFlow DEPENDS blockforest boundary core field lbm postprocessing stencil timeloop vtk sqlite ) waLBerla_add_executable( NAME CouetteFlow DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::field walberla::lbm walberla::postprocessing walberla::stencil walberla::timeloop walberla::vtk walberla::sqlite )
############## ##############
# Some tests # # Some tests #
############## ##############
waLBerla_execute_test( NO_MODULE_LABEL NAME CouetteFlowTestNoCheckRelease COMMAND $<TARGET_FILE:CouetteFlow> ${CMAKE_CURRENT_SOURCE_DIR}/TestNoCheck.dat --trt --linear-exp PROCESSES 4 CONFIGURATIONS Release RelWithDbgInfo ) waLBerla_execute_test( NO_MODULE_LABEL NAME CouetteFlowTestNoCheckRelease COMMAND $<TARGET_FILE:CouetteFlow> ${CMAKE_CURRENT_SOURCE_DIR}/TestNoCheck.dat --trt --linear-exp PROCESSES 4 CONFIGURATIONS Release RelWithDbgInfo DEPENDS_ON_TARGETS CouetteFlow )
waLBerla_execute_test( NO_MODULE_LABEL NAME CouetteFlowTestNoCheckDebug COMMAND $<TARGET_FILE:CouetteFlow> ${CMAKE_CURRENT_SOURCE_DIR}/TestNoCheck.dat --trt --linear-exp PROCESSES 4 LABELS longrun CONFIGURATIONS Debug DebugOptimized ) waLBerla_execute_test( NO_MODULE_LABEL NAME CouetteFlowTestNoCheckDebug COMMAND $<TARGET_FILE:CouetteFlow> ${CMAKE_CURRENT_SOURCE_DIR}/TestNoCheck.dat --trt --linear-exp PROCESSES 4 LABELS longrun CONFIGURATIONS Debug DebugOptimized DEPENDS_ON_TARGETS CouetteFlow )
waLBerla_execute_test( NO_MODULE_LABEL NAME CouetteFlowTest0 COMMAND $<TARGET_FILE:CouetteFlow> ${CMAKE_CURRENT_SOURCE_DIR}/Test0.dat --trt --linear-exp LABELS longrun CONFIGURATIONS Release RelWithDbgInfo ) waLBerla_execute_test( NO_MODULE_LABEL NAME CouetteFlowTest0 COMMAND $<TARGET_FILE:CouetteFlow> ${CMAKE_CURRENT_SOURCE_DIR}/Test0.dat --trt --linear-exp LABELS longrun CONFIGURATIONS Release RelWithDbgInfo DEPENDS_ON_TARGETS CouetteFlow )
waLBerla_execute_test( NO_MODULE_LABEL NAME CouetteFlowTest2 COMMAND $<TARGET_FILE:CouetteFlow> ${CMAKE_CURRENT_SOURCE_DIR}/Test2.dat --trt --linear-exp LABELS longrun verylongrun PROCESSES 4 CONFIGURATIONS Release RelWithDbgInfo ) waLBerla_execute_test( NO_MODULE_LABEL NAME CouetteFlowTest2 COMMAND $<TARGET_FILE:CouetteFlow> ${CMAKE_CURRENT_SOURCE_DIR}/Test2.dat --trt --linear-exp LABELS longrun verylongrun PROCESSES 4 CONFIGURATIONS Release RelWithDbgInfo DEPENDS_ON_TARGETS CouetteFlow )
apps/benchmarks/CouetteFlow/CouetteFlow.cpp
View file @ 7a6d3b57
...@@ -118,21 +118,21 @@ using walberla::real_t; ...@@ -118,21 +118,21 @@ using walberla::real_t;
// TYPEDEFS // // TYPEDEFS //
////////////// //////////////
typedef lbm::D3Q15< lbm::collision_model::SRT, false > D3Q15_SRT_INCOMP; using D3Q15_SRT_INCOMP = lbm::D3Q15<lbm::collision_model::SRT, false>;
typedef lbm::D3Q15< lbm::collision_model::SRT, true > D3Q15_SRT_COMP; using D3Q15_SRT_COMP = lbm::D3Q15<lbm::collision_model::SRT, true>;
typedef lbm::D3Q15< lbm::collision_model::TRT, false > D3Q15_TRT_INCOMP; using D3Q15_TRT_INCOMP = lbm::D3Q15<lbm::collision_model::TRT, false>;
typedef lbm::D3Q15< lbm::collision_model::TRT, true > D3Q15_TRT_COMP; using D3Q15_TRT_COMP = lbm::D3Q15<lbm::collision_model::TRT, true>;
typedef lbm::D3Q19< lbm::collision_model::SRT, false > D3Q19_SRT_INCOMP; using D3Q19_SRT_INCOMP = lbm::D3Q19<lbm::collision_model::SRT, false>;
typedef lbm::D3Q19< lbm::collision_model::SRT, true > D3Q19_SRT_COMP; using D3Q19_SRT_COMP = lbm::D3Q19<lbm::collision_model::SRT, true>;
typedef lbm::D3Q19< lbm::collision_model::TRT, false > D3Q19_TRT_INCOMP; using D3Q19_TRT_INCOMP = lbm::D3Q19<lbm::collision_model::TRT, false>;
typedef lbm::D3Q19< lbm::collision_model::TRT, true > D3Q19_TRT_COMP; using D3Q19_TRT_COMP = lbm::D3Q19<lbm::collision_model::TRT, true>;
typedef lbm::D3Q19< lbm::collision_model::D3Q19MRT, false > D3Q19_MRT_INCOMP; using D3Q19_MRT_INCOMP = lbm::D3Q19<lbm::collision_model::D3Q19MRT, false>;
typedef lbm::D3Q27< lbm::collision_model::SRT, false > D3Q27_SRT_INCOMP; using D3Q27_SRT_INCOMP = lbm::D3Q27<lbm::collision_model::SRT, false>;
typedef lbm::D3Q27< lbm::collision_model::SRT, true > D3Q27_SRT_COMP; using D3Q27_SRT_COMP = lbm::D3Q27<lbm::collision_model::SRT, true>;
typedef lbm::D3Q27< lbm::collision_model::TRT, false > D3Q27_TRT_INCOMP; using D3Q27_TRT_INCOMP = lbm::D3Q27<lbm::collision_model::TRT, false>;
typedef lbm::D3Q27< lbm::collision_model::TRT, true > D3Q27_TRT_COMP; using D3Q27_TRT_COMP = lbm::D3Q27<lbm::collision_model::TRT, true>;
template< typename LatticeModel_T > template< typename LatticeModel_T >
struct Types struct Types
...@@ -162,19 +162,19 @@ template< typename LatticeModel_T, class Enable = void > ...@@ -162,19 +162,19 @@ template< typename LatticeModel_T, class Enable = void >
struct StencilString; struct StencilString;
template< typename LatticeModel_T > template< typename LatticeModel_T >
struct StencilString< LatticeModel_T, typename std::enable_if< std::is_same< typename LatticeModel_T::Stencil, stencil::D3Q15 >::value >::type > struct StencilString< LatticeModel_T, std::enable_if_t< std::is_same_v< typename LatticeModel_T::Stencil, stencil::D3Q15 > > >
{ {
static const char * str() { return "D3Q15"; } static const char * str() { return "D3Q15"; }
}; };
template< typename LatticeModel_T > template< typename LatticeModel_T >
struct StencilString< LatticeModel_T, typename std::enable_if< std::is_same< typename LatticeModel_T::Stencil, stencil::D3Q19 >::value >::type > struct StencilString< LatticeModel_T, std::enable_if_t< std::is_same_v< typename LatticeModel_T::Stencil, stencil::D3Q19 > > >
{ {
static const char * str() { return "D3Q19"; } static const char * str() { return "D3Q19"; }
}; };
template< typename LatticeModel_T > template< typename LatticeModel_T >
struct StencilString< LatticeModel_T, typename std::enable_if< std::is_same< typename LatticeModel_T::Stencil, stencil::D3Q27 >::value >::type > struct StencilString< LatticeModel_T, std::enable_if_t< std::is_same_v< typename LatticeModel_T::Stencil, stencil::D3Q27 > > >
{ {
static const char * str() { return "D3Q27"; } static const char * str() { return "D3Q27"; }
}; };
...@@ -184,22 +184,22 @@ template< typename LatticeModel_T, class Enable = void > ...@@ -184,22 +184,22 @@ template< typename LatticeModel_T, class Enable = void >
struct CollisionModelString; struct CollisionModelString;
template< typename LatticeModel_T > template< typename LatticeModel_T >
struct CollisionModelString< LatticeModel_T, typename std::enable_if< std::is_same< typename LatticeModel_T::CollisionModel::tag, struct CollisionModelString< LatticeModel_T, std::enable_if_t< std::is_same_v< typename LatticeModel_T::CollisionModel::tag,
lbm::collision_model::SRT_tag >::value >::type > lbm::collision_model::SRT_tag > > >
{ {
static const char * str() { return "SRT"; } static const char * str() { return "SRT"; }
}; };
template< typename LatticeModel_T > template< typename LatticeModel_T >
struct CollisionModelString< LatticeModel_T, typename std::enable_if< std::is_same< typename LatticeModel_T::CollisionModel::tag, struct CollisionModelString< LatticeModel_T, std::enable_if_t< std::is_same_v< typename LatticeModel_T::CollisionModel::tag,
lbm::collision_model::TRT_tag >::value >::type > lbm::collision_model::TRT_tag > > >
{ {
static const char * str() { return "TRT"; } static const char * str() { return "TRT"; }
}; };
template< typename LatticeModel_T > template< typename LatticeModel_T >
struct CollisionModelString< LatticeModel_T, typename std::enable_if< std::is_same< typename LatticeModel_T::CollisionModel::tag, struct CollisionModelString< LatticeModel_T, std::enable_if_t< std::is_same_v< typename LatticeModel_T::CollisionModel::tag,
lbm::collision_model::MRT_tag >::value >::type > lbm::collision_model::MRT_tag > > >
{ {
static const char * str() { return "MRT"; } static const char * str() { return "MRT"; }
}; };
...@@ -271,10 +271,10 @@ private: ...@@ -271,10 +271,10 @@ private:
static void workloadAndMemoryAssignment( SetupBlockForest& forest, const memory_t memoryPerBlock ) static void workloadAndMemoryAssignment( SetupBlockForest& forest, const memory_t memoryPerBlock )
{ {
for( auto block = forest.begin(); block != forest.end(); ++block ) for(auto & block : forest)
{ {
block->setWorkload( numeric_cast< workload_t >( uint_t(1) << block->getLevel() ) ); block.setWorkload( numeric_cast< workload_t >( uint_t(1) << block.getLevel() ) );
block->setMemory( memoryPerBlock ); block.setMemory( memoryPerBlock );
} }
} }
...@@ -293,7 +293,7 @@ static shared_ptr< SetupBlockForest > createSetupBlockForest( const blockforest: ...@@ -293,7 +293,7 @@ static shared_ptr< SetupBlockForest > createSetupBlockForest( const blockforest:
( setup.zCells + uint_t(2) * FieldGhostLayers ) ) * memoryPerCell; ( setup.zCells + uint_t(2) * FieldGhostLayers ) ) * memoryPerCell;
forest->addRefinementSelectionFunction( refinementSelectionFunctions ); forest->addRefinementSelectionFunction( refinementSelectionFunctions );
forest->addWorkloadMemorySUIDAssignmentFunction( std::bind( workloadAndMemoryAssignment, std::placeholders::_1, memoryPerBlock ) ); forest->addWorkloadMemorySUIDAssignmentFunction([memoryPerBlock](auto & sbf) { workloadAndMemoryAssignment(sbf, memoryPerBlock); });
forest->init( AABB( real_c(0), real_c(0), real_c(0), real_c( setup.xBlocks * setup.xCells ), forest->init( AABB( real_c(0), real_c(0), real_c(0), real_c( setup.xBlocks * setup.xCells ),
real_c( setup.yBlocks * setup.yCells ), real_c( setup.yBlocks * setup.yCells ),
...@@ -370,10 +370,10 @@ class MyBoundaryHandling ...@@ -370,10 +370,10 @@ class MyBoundaryHandling
{ {
public: public:
typedef lbm::NoSlip< LatticeModel_T, flag_t > NoSlip_T; using NoSlip_T = lbm::NoSlip<LatticeModel_T, flag_t>;
typedef lbm::SimpleUBB< LatticeModel_T, flag_t > UBB_T; using UBB_T = lbm::SimpleUBB<LatticeModel_T, flag_t>;
typedef BoundaryHandling< FlagField_T, typename Types<LatticeModel_T>::Stencil_T, NoSlip_T, UBB_T > BoundaryHandling_T; using BoundaryHandling_T = BoundaryHandling<FlagField_T, typename Types<LatticeModel_T>::Stencil_T, NoSlip_T, UBB_T>;
...@@ -616,7 +616,7 @@ struct AddRefinementTimeStep ...@@ -616,7 +616,7 @@ struct AddRefinementTimeStep
} }
else else
{ {
typedef lbm::SplitSweep< LatticeModel_T, FlagField_T > Sweep_T; using Sweep_T = lbm::SplitSweep<LatticeModel_T, FlagField_T>;
auto mySweep = make_shared< Sweep_T >( pdfFieldId, flagFieldId, Fluid_Flag ); auto mySweep = make_shared< Sweep_T >( pdfFieldId, flagFieldId, Fluid_Flag );
addRefinementTimeStep< LatticeModel_T, Sweep_T >( timeloop, blocks, pdfFieldId, boundaryHandlingId, timingPool, levelwiseTimingPool, addRefinementTimeStep< LatticeModel_T, Sweep_T >( timeloop, blocks, pdfFieldId, boundaryHandlingId, timingPool, levelwiseTimingPool,
...@@ -634,10 +634,10 @@ struct AddRefinementTimeStep ...@@ -634,10 +634,10 @@ struct AddRefinementTimeStep
}; };
template< typename LatticeModel_T > template< typename LatticeModel_T >
struct AddRefinementTimeStep< LatticeModel_T, typename std::enable_if< std::is_same< typename LatticeModel_T::CollisionModel::tag, lbm::collision_model::MRT_tag >::value || struct AddRefinementTimeStep< LatticeModel_T, std::enable_if_t< std::is_same_v< typename LatticeModel_T::CollisionModel::tag, lbm::collision_model::MRT_tag > ||
std::is_same< typename LatticeModel_T::Stencil, stencil::D3Q15 >::value || std::is_same_v< typename LatticeModel_T::Stencil, stencil::D3Q15 > ||
std::is_same< typename LatticeModel_T::Stencil, stencil::D3Q27 >::value std::is_same_v< typename LatticeModel_T::Stencil, stencil::D3Q27 >
>::type > > >
{ {
static void add( SweepTimeloop & timeloop, shared_ptr< blockforest::StructuredBlockForest > & blocks, static void add( SweepTimeloop & timeloop, shared_ptr< blockforest::StructuredBlockForest > & blocks,
const BlockDataID & pdfFieldId, const BlockDataID & flagFieldId, const BlockDataID & boundaryHandlingId, const BlockDataID & pdfFieldId, const BlockDataID & flagFieldId, const BlockDataID & boundaryHandlingId,
...@@ -718,9 +718,9 @@ void run( const shared_ptr< Config > & config, const LatticeModel_T & latticeMod ...@@ -718,9 +718,9 @@ void run( const shared_ptr< Config > & config, const LatticeModel_T & latticeMod
if( vtkBeforeTimeStep ) if( vtkBeforeTimeStep )
{ {
for( auto output = vtkOutputFunctions.begin(); output != vtkOutputFunctions.end(); ++output ) for(auto & output : vtkOutputFunctions)
timeloop.addFuncBeforeTimeStep( output->second.outputFunction, std::string("VTK: ") + output->first, timeloop.addFuncBeforeTimeStep( output.second.outputFunction, std::string("VTK: ") + output.first,
output->second.requiredGlobalStates, output->second.incompatibleGlobalStates ); output.second.requiredGlobalStates, output.second.incompatibleGlobalStates );
} }
// add 'refinement' LB time step to time loop // add 'refinement' LB time step to time loop
...@@ -733,11 +733,11 @@ void run( const shared_ptr< Config > & config, const LatticeModel_T & latticeMod ...@@ -733,11 +733,11 @@ void run( const shared_ptr< Config > & config, const LatticeModel_T & latticeMod
// evaluation // evaluation
const auto exactSolutionFunction = std::bind( exactVelocity, std::placeholders::_1, blocks->getDomain(), setup.maxVelocity_L ); const auto exactSolutionFunction = [domain = blocks->getDomain(), maxVel = setup.maxVelocity_L](auto & p) { return exactVelocity(p, domain, maxVel); };
auto volumetricFlowRate = field::makeVolumetricFlowRateEvaluation< VelocityAdaptor_T, FlagField_T >( configBlock, blocks, velocityAdaptorId, auto volumetricFlowRate = field::makeVolumetricFlowRateEvaluation< VelocityAdaptor_T, FlagField_T >( configBlock, blocks, velocityAdaptorId,
flagFieldId, Fluid_Flag, flagFieldId, Fluid_Flag,
std::bind( exactFlowRate, setup.flowRate_L ), [flowRate = setup.flowRate_L] { return exactFlowRate(flowRate); },
exactSolutionFunction ); exactSolutionFunction );
volumetricFlowRate->setNormalizationFactor( real_t(1) / setup.maxVelocity_L ); volumetricFlowRate->setNormalizationFactor( real_t(1) / setup.maxVelocity_L );
volumetricFlowRate->setDomainNormalization( Vector3<real_t>( real_t(1) ) ); volumetricFlowRate->setDomainNormalization( Vector3<real_t>( real_t(1) ) );
...@@ -766,14 +766,14 @@ void run( const shared_ptr< Config > & config, const LatticeModel_T & latticeMod ...@@ -766,14 +766,14 @@ void run( const shared_ptr< Config > & config, const LatticeModel_T & latticeMod
if( !vtkBeforeTimeStep ) if( !vtkBeforeTimeStep )
{ {
for( auto output = vtkOutputFunctions.begin(); output != vtkOutputFunctions.end(); ++output ) for(auto & output : vtkOutputFunctions)
timeloop.addFuncAfterTimeStep( output->second.outputFunction, std::string("VTK: ") + output->first, timeloop.addFuncAfterTimeStep( output.second.outputFunction, std::string("VTK: ") + output.first,
output->second.requiredGlobalStates, output->second.incompatibleGlobalStates ); output.second.requiredGlobalStates, output.second.incompatibleGlobalStates );
} }
// remaining time logger // remaining time logger
const double remainingTimeLoggerFrequency = configBlock.getParameter< double >( "remainingTimeLoggerFrequency", 3.0 ); const real_t remainingTimeLoggerFrequency = configBlock.getParameter< real_t >( "remainingTimeLoggerFrequency", real_c(3.0) );
timeloop.addFuncAfterTimeStep( timing::RemainingTimeLogger( timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency ), "Remaining time logger" ); timeloop.addFuncAfterTimeStep( timing::RemainingTimeLogger( timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency ), "Remaining time logger" );
// logging right before the simulation starts // logging right before the simulation starts
......
apps/benchmarks/DEM/CMakeLists.txt
View file @ 7a6d3b57
waLBerla_add_executable( NAME DEM FILES DEM.cpp DEPENDS blockforest core pe ) waLBerla_add_executable( NAME DEM FILES DEM.cpp DEPENDS walberla::blockforest walberla::core walberla::pe )
apps/benchmarks/DEM/DEM.cpp
View file @ 7a6d3b57
...@@ -48,7 +48,7 @@ int main( int argc, char** argv ) ...@@ -48,7 +48,7 @@ int main( int argc, char** argv )
using namespace walberla; using namespace walberla;
using namespace walberla::pe; using namespace walberla::pe;
typedef std::tuple<Sphere, Plane> BodyTuple ; using BodyTuple = std::tuple<Sphere, Plane> ;
walberla::MPIManager::instance()->initializeMPI( &argc, &argv ); walberla::MPIManager::instance()->initializeMPI( &argc, &argv );
......
apps/benchmarks/FieldCommunication/CMakeLists.txt
View file @ 7a6d3b57
...@@ -4,4 +4,4 @@ waLBerla_link_files_to_builddir( "*.py" ) ...@@ -4,4 +4,4 @@ waLBerla_link_files_to_builddir( "*.py" )
waLBerla_add_executable ( NAME FieldCommunication waLBerla_add_executable ( NAME FieldCommunication
DEPENDS blockforest core domain_decomposition field postprocessing sqlite ) DEPENDS walberla::blockforest walberla::core walberla::domain_decomposition walberla::field walberla::postprocessing walberla::sqlite walberla::python_coupling )
apps/benchmarks/FieldCommunication/FieldCommunication.cpp
View file @ 7a6d3b57
...@@ -34,9 +34,9 @@ template<typename Stencil_T> ...@@ -34,9 +34,9 @@ template<typename Stencil_T>
class SingleMessageBufferedScheme class SingleMessageBufferedScheme
{ {
public: public:
typedef Stencil_T Stencil; using Stencil = Stencil_T;
SingleMessageBufferedScheme( const weak_ptr_wrapper< StructuredBlockForest > & bf, const int tag = 17953 ) SingleMessageBufferedScheme( const weak_ptr< StructuredBlockForest > & bf, const int tag = 17953 )
: blockForest_( bf ), tag_( tag ) {} : blockForest_( bf ), tag_( tag ) {}
inline void addDataToCommunicate( const shared_ptr< communication::UniformPackInfo > &packInfo ) inline void addDataToCommunicate( const shared_ptr< communication::UniformPackInfo > &packInfo )
...@@ -67,7 +67,7 @@ public: ...@@ -67,7 +67,7 @@ public:
private: private:
std::vector< shared_ptr< UniformBufferedScheme< Stencil>> > schemes_; std::vector< shared_ptr< UniformBufferedScheme< Stencil>> > schemes_;
weak_ptr_wrapper< StructuredBlockForest > blockForest_; weak_ptr< StructuredBlockForest > blockForest_;
int tag_; int tag_;
}; };
...@@ -354,8 +354,8 @@ int main( int argc, char **argv ) ...@@ -354,8 +354,8 @@ int main( int argc, char **argv )
auto databaseBlock = config->getBlock( "Database" ); auto databaseBlock = config->getBlock( "Database" );
if ( databaseBlock ) if ( databaseBlock )
{ {
for ( auto it = databaseBlock.begin(); it != databaseBlock.end(); ++it ) for (const auto & it : databaseBlock)
stringProperties[it->first] = it->second; stringProperties[it.first] = it.second;
} }
realProperties["total_min"] = real_c( timingPool["totalTime"].min()) / real_c( iterations ); realProperties["total_min"] = real_c( timingPool["totalTime"].min()) / real_c( iterations );
......
apps/benchmarks/FieldCommunication/config.py
View file @ 7a6d3b57
...@@ -34,7 +34,8 @@ sng_network = supermuc_network_spread() ...@@ -34,7 +34,8 @@ sng_network = supermuc_network_spread()
class AlreadySimulated: class AlreadySimulated:
def __init__(self, db_file, properties=('processes0*processes1*processes2', 'layout', 'ghostLayers', 'cartesianCommunicator', 'stencil', def __init__(self, db_file, properties=('processes0*processes1*processes2', 'layout', 'ghostLayers',
'cartesianCommunicator', 'stencil',
'cellsPerBlock0', 'cellsPerBlock1', 'cellsPerBlock2', 'cellsPerBlock0', 'cellsPerBlock1', 'cellsPerBlock2',
'blocksPerProcess', 'localCommunicationMode', 'singleMessage', 'blocksPerProcess', 'localCommunicationMode', 'singleMessage',
'fieldsPdf', 'fieldsPdfOpt', 'fieldsVector', 'fieldsScalar', 'fieldsPdf', 'fieldsPdfOpt', 'fieldsVector', 'fieldsScalar',
......
apps/benchmarks/FluidParticleCoupling/CMakeLists.txt
View file @ 7a6d3b57
waLBerla_link_files_to_builddir( "*.dat" ) waLBerla_link_files_to_builddir( "*.dat" )
waLBerla_python_file_generates(GeneratedLBM.py if( WALBERLA_BUILD_WITH_CODEGEN )
GeneratedLBM.cpp GeneratedLBM.h
)
waLBerla_python_file_generates(GeneratedLBMWithForce.py waLBerla_generate_target_from_python(NAME FluidParticleCouplingGeneratedLBM FILE GeneratedLBM.py
GeneratedLBMWithForce.cpp GeneratedLBMWithForce.h OUT_FILES GeneratedLBM.cpp GeneratedLBM.h
) )
waLBerla_generate_target_from_python(NAME FluidParticleCouplingGeneratedLBMWithForce FILE GeneratedLBMWithForce.py
OUT_FILES GeneratedLBMWithForce.cpp GeneratedLBMWithForce.h
)
if( WALBERLA_BUILD_WITH_CODEGEN ) waLBerla_add_executable(NAME SphereWallCollision FILES SphereWallCollision.cpp
DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling
walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk FluidParticleCouplingGeneratedLBM )
waLBerla_add_executable ( NAME SphereWallCollision FILES SphereWallCollision.cpp GeneratedLBM.py waLBerla_add_executable(NAME SettlingSphereInBox FILES SettlingSphereInBox.cpp
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk ) DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling
walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk FluidParticleCouplingGeneratedLBM )
waLBerla_add_executable ( NAME SettlingSphereInBox FILES SettlingSphereInBox.cpp GeneratedLBM.py waLBerla_add_executable(NAME SphereMovingWithPrescribedVelocity FILES SphereMovingWithPrescribedVelocity.cpp
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk ) DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::mesa_pd walberla::lbm_mesapd_coupling
walberla::postprocessing walberla::timeloop walberla::vtk FluidParticleCouplingGeneratedLBM )
waLBerla_add_executable ( NAME SphereMovingWithPrescribedVelocity FILES SphereMovingWithPrescribedVelocity.cpp GeneratedLBM.py waLBerla_add_executable(NAME LubricationForceEvaluation FILES LubricationForceEvaluation.cpp
DEPENDS blockforest boundary core domain_decomposition field lbm mesa_pd lbm_mesapd_coupling postprocessing timeloop vtk ) DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd
walberla::postprocessing walberla::timeloop walberla::vtk FluidParticleCouplingGeneratedLBM )
waLBerla_add_executable ( NAME LubricationForceEvaluation FILES LubricationForceEvaluation.cpp GeneratedLBM.py waLBerla_add_executable(NAME DragForceSphere FILES DragForceSphere.cpp
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk ) DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd
walberla::postprocessing walberla::timeloop walberla::vtk FluidParticleCouplingGeneratedLBMWithForce )
waLBerla_add_executable ( NAME DragForceSphere FILES DragForceSphere.cpp GeneratedLBMWithForce.py waLBerla_add_executable(NAME ForcesOnSphereNearPlane FILES ForcesOnSphereNearPlane.cpp
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk ) DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd
walberla::postprocessing walberla::timeloop walberla::vtk FluidParticleCouplingGeneratedLBM )
waLBerla_add_executable ( NAME ForcesOnSphereNearPlane FILES ForcesOnSphereNearPlane.cpp GeneratedLBM.py waLBerla_add_executable(NAME ObliqueWetCollision FILES ObliqueWetCollision.cpp
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk ) DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd
walberla::postprocessing walberla::timeloop walberla::vtk FluidParticleCouplingGeneratedLBM )
waLBerla_add_executable ( NAME ObliqueWetCollision FILES ObliqueWetCollision.cpp GeneratedLBM.py waLBerla_add_executable(NAME MotionSettlingSphere FILES MotionSettlingSphere.cpp
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk ) DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd
walberla::postprocessing walberla::timeloop walberla::vtk FluidParticleCouplingGeneratedLBM )
else() else()
waLBerla_add_executable ( NAME SphereWallCollision FILES SphereWallCollision.cpp waLBerla_add_executable ( NAME SphereWallCollision FILES SphereWallCollision.cpp
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk ) DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk )
waLBerla_add_executable ( NAME SettlingSphereInBox FILES SettlingSphereInBox.cpp waLBerla_add_executable ( NAME SettlingSphereInBox FILES SettlingSphereInBox.cpp
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk ) DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk )
waLBerla_add_executable ( NAME SphereMovingWithPrescribedVelocity FILES SphereMovingWithPrescribedVelocity.cpp waLBerla_add_executable ( NAME SphereMovingWithPrescribedVelocity FILES SphereMovingWithPrescribedVelocity.cpp
DEPENDS blockforest boundary core domain_decomposition field lbm mesa_pd lbm_mesapd_coupling postprocessing timeloop vtk ) DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::mesa_pd walberla::lbm_mesapd_coupling walberla::postprocessing walberla::timeloop walberla::vtk )
waLBerla_add_executable ( NAME LubricationForceEvaluation FILES LubricationForceEvaluation.cpp waLBerla_add_executable ( NAME LubricationForceEvaluation FILES LubricationForceEvaluation.cpp
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk ) DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk )
waLBerla_add_executable ( NAME DragForceSphere FILES DragForceSphere.cpp waLBerla_add_executable ( NAME DragForceSphere FILES DragForceSphere.cpp
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk ) DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk )
waLBerla_add_executable ( NAME ForcesOnSphereNearPlane FILES ForcesOnSphereNearPlane.cpp waLBerla_add_executable ( NAME ForcesOnSphereNearPlane FILES ForcesOnSphereNearPlane.cpp
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk ) DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk )
waLBerla_add_executable ( NAME ObliqueWetCollision FILES ObliqueWetCollision.cpp waLBerla_add_executable ( NAME ObliqueWetCollision FILES ObliqueWetCollision.cpp
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk ) DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk )
endif() endif()
waLBerla_add_executable ( NAME ObliqueDryCollision FILES ObliqueDryCollision.cpp waLBerla_add_executable ( NAME ObliqueDryCollision FILES ObliqueDryCollision.cpp
DEPENDS blockforest core mesa_pd postprocessing ) DEPENDS walberla::blockforest walberla::core walberla::mesa_pd walberla::postprocessing )
\ No newline at end of file
apps/benchmarks/FluidParticleCoupling/DragForceSphere.cpp
View file @ 7a6d3b57
...@@ -79,6 +79,11 @@ ...@@ -79,6 +79,11 @@
#ifdef WALBERLA_BUILD_WITH_CODEGEN #ifdef WALBERLA_BUILD_WITH_CODEGEN
#include "GeneratedLBMWithForce.h" #include "GeneratedLBMWithForce.h"
#define USE_TRTLIKE
//#define USE_D3Q27TRTLIKE
//#define USE_CUMULANTTRT
//#define USE_CUMULANT
#endif #endif
namespace drag_force_sphere_mem namespace drag_force_sphere_mem
...@@ -529,17 +534,27 @@ int main( int argc, char **argv ) ...@@ -529,17 +534,27 @@ int main( int argc, char **argv )
real_t omegaBulk = (useSRT) ? lambda_e : lbm_mesapd_coupling::omegaBulkFromOmega(omega, bulkViscRateFactor); real_t omegaBulk = (useSRT) ? lambda_e : lbm_mesapd_coupling::omegaBulkFromOmega(omega, bulkViscRateFactor);
// add omega bulk field // 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 // create the lattice model
#ifdef WALBERLA_BUILD_WITH_CODEGEN #ifdef WALBERLA_BUILD_WITH_CODEGEN
#if defined(USE_TRTLIKE) || defined(USE_D3Q27TRTLIKE)
WALBERLA_LOG_INFO_ON_ROOT("Using generated TRT-like lattice model!"); WALBERLA_LOG_INFO_ON_ROOT("Using generated TRT-like lattice model!");
WALBERLA_LOG_INFO_ON_ROOT(" - magic number " << magicNumber); WALBERLA_LOG_INFO_ON_ROOT(" - magic number " << magicNumber);
WALBERLA_LOG_INFO_ON_ROOT(" - omegaBulk = " << omegaBulk << ", bulk visc. = " << lbm_mesapd_coupling::bulkViscosityFromOmegaBulk(omegaBulk) << " (bvrf " << bulkViscRateFactor << ")"); WALBERLA_LOG_INFO_ON_ROOT(" - omegaBulk = " << omegaBulk << ", bulk visc. = " << lbm_mesapd_coupling::bulkViscosityFromOmegaBulk(omegaBulk) << " (bvrf " << bulkViscRateFactor << ")");
WALBERLA_LOG_INFO_ON_ROOT(" - lambda_e " << lambda_e << ", lambda_d " << lambda_d << ", omegaBulk " << omegaBulk ); WALBERLA_LOG_INFO_ON_ROOT(" - lambda_e " << lambda_e << ", lambda_d " << lambda_d << ", omegaBulk " << omegaBulk );
WALBERLA_LOG_INFO_ON_ROOT(" - use omega bulk adaption = " << useOmegaBulkAdaption << " (adaption layer size = " << adaptionLayerSize << ")"); WALBERLA_LOG_INFO_ON_ROOT(" - use omega bulk adaption = " << useOmegaBulkAdaption << " (adaption layer size = " << adaptionLayerSize << ")");
LatticeModel_T latticeModel = LatticeModel_T(omegaBulkFieldID, setup.extForce, lambda_d, lambda_e); LatticeModel_T latticeModel = LatticeModel_T(omegaBulkFieldID, setup.extForce, lambda_d, lambda_e);
#elif defined(USE_CUMULANTTRT)
WALBERLA_LOG_INFO_ON_ROOT("Using generated cumulant TRT lattice model!");
WALBERLA_LOG_INFO_ON_ROOT(" - magic number " << magicNumber);
WALBERLA_LOG_INFO_ON_ROOT(" - lambda_e " << lambda_e << ", lambda_d " << lambda_d );
LatticeModel_T latticeModel = LatticeModel_T(setup.extForce, lambda_d, lambda_e);
#elif defined(USE_CUMULANT)
LatticeModel_T latticeModel = LatticeModel_T(setup.extForce, omega);
#endif
#else #else
WALBERLA_LOG_INFO_ON_ROOT("Using waLBerla built-in MRT lattice model and ignoring omega bulk field since not supported!"); WALBERLA_LOG_INFO_ON_ROOT("Using waLBerla built-in MRT lattice model and ignoring omega bulk field since not supported!");
WALBERLA_LOG_INFO_ON_ROOT(" - magic number " << magicNumber); WALBERLA_LOG_INFO_ON_ROOT(" - magic number " << magicNumber);
......
apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp
View file @ 7a6d3b57
...@@ -216,13 +216,8 @@ public: ...@@ -216,13 +216,8 @@ public:
WALBERLA_MPI_SECTION() WALBERLA_MPI_SECTION()
{ {
mpi::allReduceInplace( force[0], mpi::SUM ); mpi::allReduceInplace( force, mpi::SUM );
mpi::allReduceInplace( force[1], mpi::SUM ); mpi::allReduceInplace( torque, 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 );
} }
if( fileIO_ ) if( fileIO_ )
...@@ -309,6 +304,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con ...@@ -309,6 +304,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
// create bounding planes // create bounding planes
mesa_pd::data::Particle p0 = *ps->create(true); mesa_pd::data::Particle p0 = *ps->create(true);
p0.setPosition(simulationDomain.minCorner()); 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.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,1) ));
p0.setOwner(mpi::MPIManager::instance()->rank()); p0.setOwner(mpi::MPIManager::instance()->rank());
p0.setType(0); p0.setType(0);
...@@ -317,6 +313,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con ...@@ -317,6 +313,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p1 = *ps->create(true); mesa_pd::data::Particle p1 = *ps->create(true);
p1.setPosition(simulationDomain.maxCorner()); 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.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,-1) ));
p1.setOwner(mpi::MPIManager::instance()->rank()); p1.setOwner(mpi::MPIManager::instance()->rank());
p1.setType(0); p1.setType(0);
...@@ -407,6 +404,7 @@ int main( int argc, char **argv ) ...@@ -407,6 +404,7 @@ int main( int argc, char **argv )
bool initializeVelocityProfile = false; bool initializeVelocityProfile = false;
bool useOmegaBulkAdaption = false; bool useOmegaBulkAdaption = false;
real_t adaptionLayerSize = real_t(2); real_t adaptionLayerSize = real_t(2);
std::string boundaryCondition = "CLI"; // SBB, CLI
real_t relativeChangeConvergenceEps = real_t(1e-5); real_t relativeChangeConvergenceEps = real_t(1e-5);
real_t physicalCheckingFrequency = real_t(0.1); real_t physicalCheckingFrequency = real_t(0.1);
...@@ -420,6 +418,7 @@ int main( int argc, char **argv ) ...@@ -420,6 +418,7 @@ int main( int argc, char **argv )
if( std::strcmp( argv[i], "--timesteps" ) == 0 ) { maximumNonDimTimesteps = real_c( std::atof( argv[++i] ) ); continue; } if( std::strcmp( argv[i], "--timesteps" ) == 0 ) { maximumNonDimTimesteps = real_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--wallDistance" ) == 0 ) { normalizedWallDistance = real_c( std::atof( argv[++i] ) ); continue; } if( std::strcmp( argv[i], "--wallDistance" ) == 0 ) { normalizedWallDistance = real_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--Re" ) == 0 ) { ReynoldsNumberShear = real_c( std::atof( argv[++i] ) ); continue; } if( std::strcmp( argv[i], "--Re" ) == 0 ) { ReynoldsNumberShear = real_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--boundaryCondition" ) == 0 ) { boundaryCondition = argv[++i]; continue; }
if( std::strcmp( argv[i], "--velocity" ) == 0 ) { wallVelocity = real_c( std::atof( argv[++i] ) ); continue; } if( std::strcmp( argv[i], "--velocity" ) == 0 ) { wallVelocity = real_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--xOffset" ) == 0 ) { xOffsetOfSpherePosition = real_c( std::atof( argv[++i] ) ); continue; } if( std::strcmp( argv[i], "--xOffset" ) == 0 ) { xOffsetOfSpherePosition = real_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--yOffset" ) == 0 ) { yOffsetOfSpherePosition = real_c( std::atof( argv[++i] ) ); continue; } if( std::strcmp( argv[i], "--yOffset" ) == 0 ) { yOffsetOfSpherePosition = real_c( std::atof( argv[++i] ) ); continue; }
...@@ -438,6 +437,7 @@ int main( int argc, char **argv ) ...@@ -438,6 +437,7 @@ int main( int argc, char **argv )
WALBERLA_CHECK_GREATER_EQUAL(normalizedWallDistance, real_t(0.5)); WALBERLA_CHECK_GREATER_EQUAL(normalizedWallDistance, real_t(0.5));
WALBERLA_CHECK_GREATER_EQUAL(ReynoldsNumberShear, real_t(0)); WALBERLA_CHECK_GREATER_EQUAL(ReynoldsNumberShear, real_t(0));
WALBERLA_CHECK_GREATER_EQUAL(diameter, real_t(0)); WALBERLA_CHECK_GREATER_EQUAL(diameter, real_t(0));
WALBERLA_CHECK(boundaryCondition == "SBB" || boundaryCondition == "CLI");
////////////////////////// //////////////////////////
// NUMERICAL PARAMETERS // // NUMERICAL PARAMETERS //
...@@ -484,6 +484,7 @@ int main( int argc, char **argv ) ...@@ -484,6 +484,7 @@ int main( int argc, char **argv )
WALBERLA_LOG_INFO_ON_ROOT(" - use omega bulk adaption = " << useOmegaBulkAdaption << " (adaption layer size = " << adaptionLayerSize << ")"); WALBERLA_LOG_INFO_ON_ROOT(" - use omega bulk adaption = " << useOmegaBulkAdaption << " (adaption layer size = " << adaptionLayerSize << ")");
WALBERLA_LOG_INFO_ON_ROOT(" - sphere diameter = " << diameter << ", position = " << initialPosition << " ( xOffset = " << xOffsetOfSpherePosition << ", yOffset = " << yOffsetOfSpherePosition << " )"); WALBERLA_LOG_INFO_ON_ROOT(" - sphere diameter = " << diameter << ", position = " << initialPosition << " ( xOffset = " << xOffsetOfSpherePosition << ", yOffset = " << yOffsetOfSpherePosition << " )");
WALBERLA_LOG_INFO_ON_ROOT(" - base folder VTK = " << baseFolderVTK << ", base folder logging = " << baseFolderLogging ); WALBERLA_LOG_INFO_ON_ROOT(" - base folder VTK = " << baseFolderVTK << ", base folder logging = " << baseFolderLogging );
WALBERLA_LOG_INFO_ON_ROOT(" - boundary condition = " << boundaryCondition );
/////////////////////////// ///////////////////////////
// BLOCK STRUCTURE SETUP // // BLOCK STRUCTURE SETUP //
...@@ -562,7 +563,7 @@ int main( int argc, char **argv ) ...@@ -562,7 +563,7 @@ int main( int argc, char **argv )
//////////////////////// ////////////////////////
// add omega bulk field // 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 // create the lattice model
real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega ); real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega );
...@@ -594,11 +595,21 @@ int main( int argc, char **argv ) ...@@ -594,11 +595,21 @@ int main( int argc, char **argv )
lbm_mesapd_coupling::MovingParticleMappingKernel<BoundaryHandling_T> movingParticleMappingKernel(blocks, boundaryHandlingID, particleFieldID); lbm_mesapd_coupling::MovingParticleMappingKernel<BoundaryHandling_T> movingParticleMappingKernel(blocks, boundaryHandlingID, particleFieldID);
lbm_mesapd_coupling::AverageHydrodynamicForceTorqueKernel averageHydrodynamicForceTorque; lbm_mesapd_coupling::AverageHydrodynamicForceTorqueKernel averageHydrodynamicForceTorque;
// map sphere into the LBM simulation if(boundaryCondition == "CLI")
ps->forEachParticle(false, lbm_mesapd_coupling::RegularParticlesSelector(), *accessor, movingParticleMappingKernel, *accessor, CLI_Flag); {
// map sphere into the LBM simulation
ps->forEachParticle(false, lbm_mesapd_coupling::RegularParticlesSelector(), *accessor, movingParticleMappingKernel, *accessor, CLI_Flag);
// map planes
ps->forEachParticle(false, lbm_mesapd_coupling::GlobalParticlesSelector(), *accessor, movingParticleMappingKernel, *accessor, CLI_Flag);
}
else
{
// map sphere into the LBM simulation
ps->forEachParticle(false, lbm_mesapd_coupling::RegularParticlesSelector(), *accessor, movingParticleMappingKernel, *accessor, SBB_Flag);
// map planes
ps->forEachParticle(false, lbm_mesapd_coupling::GlobalParticlesSelector(), *accessor, movingParticleMappingKernel, *accessor, SBB_Flag);
}
// map planes
ps->forEachParticle(false, lbm_mesapd_coupling::GlobalParticlesSelector(), *accessor, movingParticleMappingKernel, *accessor, CLI_Flag);
// initialize Couette velocity profile in whole domain // initialize Couette velocity profile in whole domain
if(initializeVelocityProfile) if(initializeVelocityProfile)
...@@ -661,13 +672,13 @@ int main( int argc, char **argv ) ...@@ -661,13 +672,13 @@ int main( int argc, char **argv )
timeloop.add() << Sweep( makeSharedSweep( lbmSweep ), "cell-wise LB sweep" ); timeloop.add() << Sweep( makeSharedSweep( lbmSweep ), "cell-wise LB sweep" );
#endif #endif
// add force evaluation and logging // add force evaluation and logging
real_t normalizationFactor = math::pi / real_t(8) * densityFluid * shearRate * shearRate * wallDistance * wallDistance * diameter * diameter ; real_t normalizationFactor = math::pi / real_t(8) * densityFluid * shearRate * shearRate * wallDistance * wallDistance * diameter * diameter ;
std::string loggingFileName( baseFolderLogging + "/LoggingForcesNearPlane"); std::string loggingFileName( baseFolderLogging + "/LoggingForcesNearPlane");
loggingFileName += "_D" + std::to_string(uint_c(diameter)); loggingFileName += "_D" + std::to_string(uint_c(diameter));
loggingFileName += "_Re" + std::to_string(uint_c(ReynoldsNumberShear)); loggingFileName += "_Re" + std::to_string(uint_c(ReynoldsNumberShear));
loggingFileName += "_WD" + std::to_string(uint_c(normalizedWallDistance * real_t(1000))); loggingFileName += "_WD" + std::to_string(uint_c(normalizedWallDistance * real_t(1000)));
loggingFileName += "_" + boundaryCondition;
loggingFileName += "_bvrf" + std::to_string(uint_c(bulkViscRateFactor)); loggingFileName += "_bvrf" + std::to_string(uint_c(bulkViscRateFactor));
if(useOmegaBulkAdaption) loggingFileName += "_uOBA" + std::to_string(uint_c(adaptionLayerSize)); if(useOmegaBulkAdaption) loggingFileName += "_uOBA" + std::to_string(uint_c(adaptionLayerSize));
loggingFileName += "_mn" + std::to_string(magicNumber); loggingFileName += "_mn" + std::to_string(magicNumber);
...@@ -796,6 +807,7 @@ int main( int argc, char **argv ) ...@@ -796,6 +807,7 @@ int main( int argc, char **argv )
resultFileName += "_D" + std::to_string(uint_c(diameter)); resultFileName += "_D" + std::to_string(uint_c(diameter));
resultFileName += "_Re" + std::to_string(uint_c(ReynoldsNumberShear)); resultFileName += "_Re" + std::to_string(uint_c(ReynoldsNumberShear));
resultFileName += "_WD" + std::to_string(uint_c(normalizedWallDistance * real_t(1000))); resultFileName += "_WD" + std::to_string(uint_c(normalizedWallDistance * real_t(1000)));
resultFileName += "_" + boundaryCondition;
resultFileName += "_bvrf" + std::to_string(uint_c(bulkViscRateFactor)); resultFileName += "_bvrf" + std::to_string(uint_c(bulkViscRateFactor));
if(useOmegaBulkAdaption) resultFileName += "_uOBA" + std::to_string(uint_c(adaptionLayerSize)); if(useOmegaBulkAdaption) resultFileName += "_uOBA" + std::to_string(uint_c(adaptionLayerSize));
resultFileName += "_mn" + std::to_string(magicNumber); resultFileName += "_mn" + std::to_string(magicNumber);
......
apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py
View file @ 7a6d3b57
import sympy as sp import sympy as sp
from sympy.core.cache import clear_cache
import pystencils as ps import pystencils as ps
from lbmpy.creationfunctions import create_lb_method_from_existing, create_lb_method from lbmpy.creationfunctions import LBMConfig, LBMOptimisation, Method, Stencil, create_lb_method
from lbmpy_walberla import generate_lattice_model from lbmpy_walberla import generate_lattice_model
from pystencils_walberla import CodeGeneration from pystencils_walberla import CodeGeneration
from pystencils_walberla import get_vectorize_instruction_set
from lbmpy.creationfunctions import create_lb_collision_rule from lbmpy.creationfunctions import create_lb_collision_rule
from lbmpy.moments import is_even, get_order, MOMENT_SYMBOLS from lbmpy.moments import is_even, get_order, MOMENT_SYMBOLS
from lbmpy.stencils import get_stencil from lbmpy.stencils import LBStencil
from collections import OrderedDict
with CodeGeneration() as ctx: with CodeGeneration() as ctx:
omegaVisc = sp.Symbol("omega_visc")
omegaBulk = ps.fields("omega_bulk: [3D]", layout='fzyx')
omegaMagic = sp.Symbol("omega_magic")
stencil = get_stencil("D3Q19", 'walberla')
x, y, z = MOMENT_SYMBOLS generatedMethod = 'TRTlike'
one = sp.Rational(1, 1) # generatedMethod = 'D3Q27TRTlike'
sq = x ** 2 + y ** 2 + z ** 2 # generatedMethod = 'cumulant'
moments = [
[one, x, y, z], # [0, 3, 5, 7]
[sq - 1], # [1]
[3 * sq ** 2 - 6 * sq + 1], # [2]
[(3 * sq - 5) * x, (3 * sq - 5) * y, (3 * sq - 5) * z], # [4, 6, 8]
[3 * x ** 2 - sq, y ** 2 - z ** 2, x * y, y * z, x * z], # [9, 11, 13, 14, 15]
[(2 * sq - 3) * (3 * x ** 2 - sq), (2 * sq - 3) * (y ** 2 - z ** 2)], # [10, 12]
[(y ** 2 - z ** 2) * x, (z ** 2 - x ** 2) * y, (x ** 2 - y ** 2) * z] # [16, 17, 18]
]
# relaxation rate for first group of moments (1,x,y,z) is set to zero internally clear_cache()
relaxation_rates=[omegaBulk.center_vector, omegaBulk.center_vector, omegaMagic, omegaVisc, omegaVisc, omegaMagic]
methodWithForce = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False, dtype_string = 'float64' if ctx.double_accuracy else 'float32'
nested_moments=moments, relaxation_rates=relaxation_rates)
#print(methodWithForce.relaxation_rates) cpu_vectorize_info = {'instruction_set': get_vectorize_instruction_set(ctx)}
#print(methodWithForce.moment_matrix)
collision_rule = create_lb_collision_rule(lb_method=methodWithForce) if generatedMethod == 'TRTlike':
generate_lattice_model(ctx, 'GeneratedLBM', collision_rule) omegaVisc = sp.Symbol('omega_visc')
omegaBulk = ps.fields(f'omega_bulk: {dtype_string}[3D]', layout='fzyx')
omegaMagic = sp.Symbol('omega_magic')
stencil = LBStencil(Stencil.D3Q19)
x, y, z = MOMENT_SYMBOLS
one = sp.Rational(1, 1)
sq = x ** 2 + y ** 2 + z ** 2
moments = [
[one, x, y, z], # [0, 3, 5, 7]
[sq - 1], # [1]
[3 * sq ** 2 - 6 * sq + 1], # [2]
[(3 * sq - 5) * x, (3 * sq - 5) * y, (3 * sq - 5) * z], # [4, 6, 8]
[3 * x ** 2 - sq, y ** 2 - z ** 2, x * y, y * z, x * z], # [9, 11, 13, 14, 15]
[(2 * sq - 3) * (3 * x ** 2 - sq), (2 * sq - 3) * (y ** 2 - z ** 2)], # [10, 12]
[(y ** 2 - z ** 2) * x, (z ** 2 - x ** 2) * y, (x ** 2 - y ** 2) * z] # [16, 17, 18]
]
# relaxation rate for first group of moments (1,x,y,z) is set to zero internally
relaxation_rates = [omegaBulk.center, omegaBulk.center,
omegaMagic, omegaVisc, omegaVisc, omegaMagic]
lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, continuous_equilibrium=False,
zero_centered=False, delta_equilibrium=False,
nested_moments=moments, relaxation_rates=relaxation_rates)
lbm_opt = LBMOptimisation(cse_global=True)
methodWithForce = create_lb_method(lbm_config=lbm_config)
# print(methodWithForce.relaxation_rates)
# print(methodWithForce.moment_matrix)
collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
generate_lattice_model(ctx, 'GeneratedLBM', collision_rule, field_layout='fzyx', cpu_vectorize_info=cpu_vectorize_info)
if generatedMethod == 'D3Q27TRTlike':
omegaVisc = sp.Symbol('omega_visc')
omegaBulk = ps.fields(f'omega_bulk: {dtype_string}[3D]', layout='fzyx')
omegaMagic = sp.Symbol('omega_magic')
stencil = LBStencil(Stencil.D3Q27)
relaxation_rates = [omegaVisc, omegaBulk.center_vector, omegaMagic, omegaVisc, omegaMagic, omegaVisc]
lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, maxwellian_moments=False, weighted=True,
compressible=False, relaxation_rates=relaxation_rates)
lbm_opt = LBMOptimisation(cse_global=True)
methodWithForce = create_lb_method(lbm_config=lbm_config)
collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
generate_lattice_model(ctx, 'GeneratedLBM', collision_rule, field_layout='fzyx',
cpu_vectorize_info=cpu_vectorize_info)
# using create_with_continuous_maxwellian_eq_moments won't work probably
if generatedMethod == 'cumulant':
x, y, z = MOMENT_SYMBOLS
cumulants = [0] * 27
cumulants[0] = sp.sympify(1) # 000
cumulants[1] = x # 100
cumulants[2] = y # 010
cumulants[3] = z # 001
cumulants[4] = x*y # 110
cumulants[5] = x*z # 101
cumulants[6] = y*z # 011
cumulants[7] = x**2 - y**2 # 200 - 020
cumulants[8] = x**2 - z**2 # 200 - 002
cumulants[9] = x**2 + y**2 + z**2 # 200 + 020 + 002
cumulants[10] = x*y**2 + x*z**2 # 120 + 102
cumulants[11] = x**2*y + y*z**2 # 210 + 012
cumulants[12] = x**2*z + y**2*z # 201 + 021
cumulants[13] = x*y**2 - x*z**2 # 120 - 102
cumulants[14] = x**2*y - y*z**2 # 210 - 012
cumulants[15] = x**2*z - y**2*z # 201 - 021
cumulants[16] = x*y*z # 111
cumulants[17] = x**2*y**2 - 2*x**2*z**2 + y**2*z**2 # 220- 2*202 +022
cumulants[18] = x**2*y**2 + x**2*z**2 - 2*y**2*z**2 # 220 + 202 - 2*002
cumulants[19] = x**2*y**2 + x**2*z**2 + y**2*z**2 # 220 + 202 + 022
cumulants[20] = x**2*y*z # 211
cumulants[21] = x*y**2*z # 121
cumulants[22] = x*y*z**2 # 112
cumulants[23] = x**2*y**2*z # 221
cumulants[24] = x**2*y*z**2 # 212
cumulants[25] = x*y**2*z**2 # 122
cumulants[26] = x**2*y**2*z**2 # 222
def get_relaxation_rate(cumulant, omega):
if get_order(cumulant) <= 1:
return 0
elif get_order(cumulant) == 2 and cumulant != x**2 + y**2 + z**2:
return omega
else:
return 1
stencil = LBStencil(Stencil.D3Q27)
omega = sp.Symbol('omega')
rr_dict = OrderedDict((c, get_relaxation_rate(c, omega))
for c in cumulants)
from lbmpy.methods import create_with_continuous_maxwellian_eq_moments
my_method = create_with_continuous_maxwellian_eq_moments(stencil, rr_dict, cumulant=True, compressible=True,)
collision_rule = create_lb_collision_rule(lb_method=my_method,
optimization={'cse_global': True,
'cse_pdfs': False})
print(my_method.relaxation_rates)
generate_lattice_model(ctx, 'GeneratedLBM', collision_rule, field_layout='fzyx',
cpu_vectorize_info=cpu_vectorize_info)
apps/benchmarks/FluidParticleCoupling/GeneratedLBMWithForce.py
View file @ 7a6d3b57
import sympy as sp import sympy as sp
from sympy.core.cache import clear_cache
import pystencils as ps import pystencils as ps
from lbmpy.creationfunctions import create_lb_method_from_existing, create_lb_method from lbmpy.creationfunctions import LBMConfig, LBMOptimisation, ForceModel, Method, Stencil, create_lb_method
from lbmpy_walberla import generate_lattice_model from lbmpy_walberla import generate_lattice_model
from pystencils_walberla import CodeGeneration from pystencils_walberla import CodeGeneration
from pystencils_walberla import get_vectorize_instruction_set
from lbmpy.creationfunctions import create_lb_collision_rule from lbmpy.creationfunctions import create_lb_collision_rule
from lbmpy.moments import is_even, get_order, MOMENT_SYMBOLS from lbmpy.moments import MOMENT_SYMBOLS, is_even, get_order
from lbmpy.stencils import get_stencil from lbmpy.stencils import LBStencil
from lbmpy.forcemodels import Luo
from collections import OrderedDict
with CodeGeneration() as ctx: with CodeGeneration() as ctx:
forcing=(sp.symbols("fx"),0,0) forcing = (sp.symbols('fx'), 0, 0)
omegaVisc = sp.Symbol("omega_visc")
omegaBulk = ps.fields("omega_bulk: [3D]", layout='fzyx') generatedMethod = 'TRTlike'
omegaMagic = sp.Symbol("omega_magic") # generatedMethod = 'D3Q27TRTlike'
stencil = get_stencil("D3Q19", 'walberla') # generatedMethod = 'cumulant'
# generatedMethod = 'cumulantTRT'
forcemodel=Luo(forcing)
print('Generating ' + generatedMethod + ' LBM method')
x, y, z = MOMENT_SYMBOLS
one = sp.Rational(1, 1) clear_cache()
sq = x ** 2 + y ** 2 + z ** 2
moments = [ dtype_string = 'float64' if ctx.double_accuracy else 'float32'
[one, x, y, z], # [0, 3, 5, 7]
[sq - 1], # [1] cpu_vectorize_info = {'instruction_set': get_vectorize_instruction_set(ctx)}
[3 * sq ** 2 - 6 * sq + 1], # [2]
[(3 * sq - 5) * x, (3 * sq - 5) * y, (3 * sq - 5) * z], # [4, 6, 8] if generatedMethod == 'TRTlike':
[3 * x ** 2 - sq, y ** 2 - z ** 2, x * y, y * z, x * z], # [9, 11, 13, 14, 15] omegaVisc = sp.Symbol('omega_visc')
[(2 * sq - 3) * (3 * x ** 2 - sq), (2 * sq - 3) * (y ** 2 - z ** 2)], # [10, 12] omegaBulk = ps.fields(f'omega_bulk: {dtype_string}[3D]', layout='fzyx')
[(y ** 2 - z ** 2) * x, (z ** 2 - x ** 2) * y, (x ** 2 - y ** 2) * z] # [16, 17, 18] omegaMagic = sp.Symbol('omega_magic')
] stencil = LBStencil(Stencil.D3Q19)
# relaxation rate for first group of moments (1,x,y,z) is set to zero internally x, y, z = MOMENT_SYMBOLS
relaxation_rates=[omegaBulk.center_vector, omegaBulk.center_vector, omegaMagic, omegaVisc, omegaVisc, omegaMagic] one = sp.Rational(1, 1)
sq = x ** 2 + y ** 2 + z ** 2
methodWithForce = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False, moments = [
force_model=forcemodel, nested_moments=moments, relaxation_rates=relaxation_rates) [one, x, y, z], # [0, 3, 5, 7]
[sq - 1], # [1]
#print(methodWithForce.relaxation_rates) [3 * sq ** 2 - 6 * sq + 1], # [2]
#print(methodWithForce.moment_matrix) [(3 * sq - 5) * x, (3 * sq - 5) * y, (3 * sq - 5) * z], # [4, 6, 8]
[3 * x ** 2 - sq, y ** 2 - z ** 2, x * y, y * z, x * z], # [9, 11, 13, 14, 15]
collision_rule = create_lb_collision_rule(lb_method=methodWithForce) [(2 * sq - 3) * (3 * x ** 2 - sq), (2 * sq - 3) * (y ** 2 - z ** 2)], # [10, 12]
generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule) [(y ** 2 - z ** 2) * x, (z ** 2 - x ** 2) * y, (x ** 2 - y ** 2) * z] # [16, 17, 18]
]
# relaxation rate for first group of moments (1,x,y,z) is set to zero internally
relaxation_rates = [omegaBulk.center, omegaBulk.center,
omegaMagic, omegaVisc, omegaVisc, omegaMagic]
lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, continuous_equilibrium=False,
zero_centered=False, delta_equilibrium=False,
force=forcing, force_model=ForceModel.LUO,
nested_moments=moments, relaxation_rates=relaxation_rates)
lbm_opt = LBMOptimisation(cse_global=True)
methodWithForce = create_lb_method(lbm_config=lbm_config)
# print(methodWithForce.relaxation_rates)
# print(methodWithForce.moment_matrix)
collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule, field_layout='fzyx',
cpu_vectorize_info=cpu_vectorize_info)
if generatedMethod == 'D3Q27TRTlike':
omegaVisc = sp.Symbol('omega_visc')
omegaBulk = ps.fields(f'omega_bulk: {dtype_string}[3D]', layout='fzyx')
omegaMagic = sp.Symbol('omega_magic')
stencil = LBStencil(Stencil.D3Q27)
relaxation_rates = [omegaVisc, omegaBulk.center_vector, omegaMagic, omegaVisc, omegaMagic, omegaVisc]
lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, maxwellian_moments=False, weighted=True,
force=forcing, force_model=ForceModel.LUO,
compressible=False, relaxation_rates=relaxation_rates)
lbm_opt = LBMOptimisation(cse_global=True)
methodWithForce = create_lb_method(lbm_config=lbm_config)
collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule, field_layout='fzyx',
cpu_vectorize_info=cpu_vectorize_info)
if generatedMethod == 'cumulant':
x, y, z = MOMENT_SYMBOLS
cumulants = [0] * 27
cumulants[0] = sp.sympify(1) # 000
cumulants[1] = x # 100
cumulants[2] = y # 010
cumulants[3] = z # 001
cumulants[4] = x*y # 110
cumulants[5] = x*z # 101
cumulants[6] = y*z # 011
cumulants[7] = x**2 - y**2 # 200 - 020
cumulants[8] = x**2 - z**2 # 200 - 002
cumulants[9] = x**2 + y**2 + z**2 # 200 + 020 + 002
cumulants[10] = x*y**2 + x*z**2 # 120 + 102
cumulants[11] = x**2*y + y*z**2 # 210 + 012
cumulants[12] = x**2*z + y**2*z # 201 + 021
cumulants[13] = x*y**2 - x*z**2 # 120 - 102
cumulants[14] = x**2*y - y*z**2 # 210 - 012
cumulants[15] = x**2*z - y**2*z # 201 - 021
cumulants[16] = x*y*z # 111
cumulants[17] = x**2*y**2 - 2*x**2*z**2 + y**2*z**2 # 220- 2*202 +022
cumulants[18] = x**2*y**2 + x**2*z**2 - 2*y**2*z**2 # 220 + 202 - 2*002
cumulants[19] = x**2*y**2 + x**2*z**2 + y**2*z**2 # 220 + 202 + 022
cumulants[20] = x**2*y*z # 211
cumulants[21] = x*y**2*z # 121
cumulants[22] = x*y*z**2 # 112
cumulants[23] = x**2*y**2*z # 221
cumulants[24] = x**2*y*z**2 # 212
cumulants[25] = x*y**2*z**2 # 122
cumulants[26] = x**2*y**2*z**2 # 222
def get_relaxation_rate(cumulant, omega):
if get_order(cumulant) <= 1:
return 0
elif get_order(cumulant) == 2 and cumulant != x**2 + y**2 + z**2:
return omega
else:
return 1
stencil = LBStencil(Stencil.D3Q27)
omega = sp.Symbol('omega')
rr_dict = OrderedDict((c, get_relaxation_rate(c, omega))
for c in cumulants)
from lbmpy.methods import create_with_continuous_maxwellian_eq_moments
my_method = create_with_continuous_maxwellian_eq_moments(stencil, rr_dict, cumulant=True, compressible=True,
force_model=forcemodel)
collision_rule = create_lb_collision_rule(lb_method=my_method,
optimization={'cse_global': True,
'cse_pdfs': False})
print(my_method.relaxation_rates)
generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule, field_layout='fzyx',
cpu_vectorize_info=cpu_vectorize_info)
if generatedMethod == 'cumulantTRT':
x, y, z = MOMENT_SYMBOLS
cumulants = [0] * 27
cumulants[0] = sp.sympify(1) # 000
cumulants[1] = x # 100
cumulants[2] = y # 010
cumulants[3] = z # 001
cumulants[4] = x*y # 110
cumulants[5] = x*z # 101
cumulants[6] = y*z # 011
cumulants[7] = x**2 - y**2 # 200 - 020
cumulants[8] = x**2 - z**2 # 200 - 002
cumulants[9] = x**2 + y**2 + z**2 # 200 + 020 + 002
cumulants[10] = x*y**2 + x*z**2 # 120 + 102
cumulants[11] = x**2*y + y*z**2 # 210 + 012
cumulants[12] = x**2*z + y**2*z # 201 + 021
cumulants[13] = x*y**2 - x*z**2 # 120 - 102
cumulants[14] = x**2*y - y*z**2 # 210 - 012
cumulants[15] = x**2*z - y**2*z # 201 - 021
cumulants[16] = x*y*z # 111
cumulants[17] = x**2*y**2 - 2*x**2*z**2 + y**2*z**2 # 220- 2*202 +022
cumulants[18] = x**2*y**2 + x**2*z**2 - 2*y**2*z**2 # 220 + 202 - 2*002
cumulants[19] = x**2*y**2 + x**2*z**2 + y**2*z**2 # 220 + 202 + 022
cumulants[20] = x**2*y*z # 211
cumulants[21] = x*y**2*z # 121
cumulants[22] = x*y*z**2 # 112
cumulants[23] = x**2*y**2*z # 221
cumulants[24] = x**2*y*z**2 # 212
cumulants[25] = x*y**2*z**2 # 122
cumulants[26] = x**2*y**2*z**2 # 222
def get_relaxation_rate(cumulant, omegaVisc, omegaMagic):
if get_order(cumulant) <= 1:
return 0
elif is_even(cumulant):
return omegaVisc
else:
return omegaMagic
stencil = LBStencil(Stencil.D3Q27)
omegaVisc = sp.Symbol('omega_visc')
omegaMagic = sp.Symbol('omega_magic')
rr_dict = OrderedDict((c, get_relaxation_rate(c, omegaVisc, omegaMagic))
for c in cumulants)
from lbmpy.methods import create_with_continuous_maxwellian_eq_moments
my_method = create_with_continuous_maxwellian_eq_moments(stencil, rr_dict, cumulant=True, compressible=True,
force_model=forcemodel)
collision_rule = create_lb_collision_rule(lb_method=my_method,
optimization={'cse_global': True,
'cse_pdfs': False})
print(my_method.relaxation_rates)
generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule, field_layout='fzyx',
cpu_vectorize_info=cpu_vectorize_info)
apps/benchmarks/FluidParticleCoupling/LubricationForceEvaluation.cpp
View file @ 7a6d3b57
...@@ -52,6 +52,7 @@ ...@@ -52,6 +52,7 @@
#include "lbm_mesapd_coupling/mapping/ParticleMapping.h" #include "lbm_mesapd_coupling/mapping/ParticleMapping.h"
#include "lbm_mesapd_coupling/momentum_exchange_method/MovingParticleMapping.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/boundary/CurvedLinear.h"
#include "lbm_mesapd_coupling/momentum_exchange_method/boundary/SimpleBB.h"
#include "lbm_mesapd_coupling/utility/ParticleSelector.h" #include "lbm_mesapd_coupling/utility/ParticleSelector.h"
#include "lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h" #include "lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
#include "lbm_mesapd_coupling/utility/LubricationCorrectionKernel.h" #include "lbm_mesapd_coupling/utility/LubricationCorrectionKernel.h"
...@@ -115,7 +116,8 @@ const uint_t FieldGhostLayers = 1; ...@@ -115,7 +116,8 @@ const uint_t FieldGhostLayers = 1;
const FlagUID Fluid_Flag( "fluid" ); const FlagUID Fluid_Flag( "fluid" );
const FlagUID FreeSlip_Flag( "free slip" ); const FlagUID FreeSlip_Flag( "free slip" );
const FlagUID MO_Flag( "moving obstacle" ); const FlagUID MO_SBB_Flag( "moving obstacle sbb" );
const FlagUID MO_CLI_Flag( "moving obstacle cli" );
const FlagUID FormerMO_Flag( "former moving obstacle" ); const FlagUID FormerMO_Flag( "former moving obstacle" );
///////////////////////////////////// /////////////////////////////////////
...@@ -127,8 +129,9 @@ class MyBoundaryHandling ...@@ -127,8 +129,9 @@ class MyBoundaryHandling
public: public:
using FreeSlip_T = lbm::FreeSlip< LatticeModel_T, FlagField_T>; using FreeSlip_T = lbm::FreeSlip< LatticeModel_T, FlagField_T>;
using MO_T = lbm_mesapd_coupling::CurvedLinear< LatticeModel_T, FlagField_T, ParticleAccessor_T >; using MO_SBB_T = lbm_mesapd_coupling::SimpleBB< LatticeModel_T, FlagField_T, ParticleAccessor_T >;
using Type = BoundaryHandling< FlagField_T, Stencil_T, FreeSlip_T, MO_T >; using MO_CLI_T = lbm_mesapd_coupling::CurvedLinear< LatticeModel_T, FlagField_T, ParticleAccessor_T >;
using Type = BoundaryHandling< FlagField_T, Stencil_T, FreeSlip_T, MO_SBB_T, MO_CLI_T >;
MyBoundaryHandling( const BlockDataID & flagFieldID, const BlockDataID & pdfFieldID, MyBoundaryHandling( const BlockDataID & flagFieldID, const BlockDataID & pdfFieldID,
const BlockDataID & particleFieldID, const shared_ptr<ParticleAccessor_T>& ac) : const BlockDataID & particleFieldID, const shared_ptr<ParticleAccessor_T>& ac) :
...@@ -147,7 +150,8 @@ public: ...@@ -147,7 +150,8 @@ public:
Type * handling = new Type( "moving obstacle boundary handling", flagField, fluid, Type * handling = new Type( "moving obstacle boundary handling", flagField, fluid,
FreeSlip_T( "FreeSlip", FreeSlip_Flag, pdfField, flagField, fluid ), FreeSlip_T( "FreeSlip", FreeSlip_Flag, pdfField, flagField, fluid ),
MO_T( "MO", MO_Flag, pdfField, flagField, particleField, ac_, fluid, *storage, *block ) ); MO_SBB_T( "SBB", MO_SBB_Flag, pdfField, flagField, particleField, ac_, fluid, *storage, *block ),
MO_CLI_T( "CLI", MO_CLI_Flag, pdfField, flagField, particleField, ac_, fluid, *storage, *block ) );
const auto freeslip = flagField->getFlag( FreeSlip_Flag ); const auto freeslip = flagField->getFlag( FreeSlip_Flag );
...@@ -237,6 +241,7 @@ int main( int argc, char **argv ) ...@@ -237,6 +241,7 @@ int main( int argc, char **argv )
real_t magicNumber = real_t(3)/real_t(16); real_t magicNumber = real_t(3)/real_t(16);
bool useOmegaBulkAdaption = false; bool useOmegaBulkAdaption = false;
real_t adaptionLayerSize = real_t(2); real_t adaptionLayerSize = real_t(2);
bool useSBB = false;
// 1: translation in normal direction -> normal Lubrication force // 1: translation in normal direction -> normal Lubrication force
// 2: translation in tangential direction -> tangential Lubrication force and torque // 2: translation in tangential direction -> tangential Lubrication force and torque
...@@ -246,19 +251,20 @@ int main( int argc, char **argv ) ...@@ -246,19 +251,20 @@ int main( int argc, char **argv )
for( int i = 1; i < argc; ++i ) for( int i = 1; i < argc; ++i )
{ {
if( std::strcmp( argv[i], "--sphWallTest" ) == 0 ) { sphSphTest = false; continue; } if( std::strcmp( argv[i], "--sphWallTest" ) == 0 ) { sphSphTest = false; continue; }
if( std::strcmp( argv[i], "--noLogging" ) == 0 ) { fileIO = false; continue; } if( std::strcmp( argv[i], "--noLogging" ) == 0 ) { fileIO = false; continue; }
if( std::strcmp( argv[i], "--vtkIOFreq" ) == 0 ) { vtkIOFreq = uint_c( std::atof( argv[++i] ) ); continue; } if( std::strcmp( argv[i], "--vtkIOFreq" ) == 0 ) { vtkIOFreq = uint_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--setup" ) == 0 ) { setup = uint_c( std::atof( argv[++i] ) ); continue; } if( std::strcmp( argv[i], "--setup" ) == 0 ) { setup = uint_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--baseFolder" ) == 0 ) { baseFolder = argv[++i]; continue; } if( std::strcmp( argv[i], "--baseFolder" ) == 0 ) { baseFolder = argv[++i]; continue; }
if( std::strcmp( argv[i], "--diameter" ) == 0 ) { radius = real_t(0.5)*real_c(std::atof(argv[++i])); continue; } if( std::strcmp( argv[i], "--diameter" ) == 0 ) { radius = real_t(0.5)*real_c(std::atof(argv[++i])); continue; }
if( std::strcmp( argv[i], "--gapSize" ) == 0 ) { gapSize = real_c(std::atof(argv[++i])); continue; } if( std::strcmp( argv[i], "--gapSize" ) == 0 ) { gapSize = real_c(std::atof(argv[++i])); continue; }
if( std::strcmp( argv[i], "--bulkViscRateFactor" ) == 0 ) { bulkViscRateFactor = real_c(std::atof( argv[++i] ) ); continue; } if( std::strcmp( argv[i], "--bulkViscRateFactor" ) == 0 ) { bulkViscRateFactor = real_c(std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--tau" ) == 0 ) { tau = real_c(std::atof( argv[++i] ) ); continue; } if( std::strcmp( argv[i], "--tau" ) == 0 ) { tau = real_c(std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--fileName" ) == 0 ) { fileNameEnding = argv[++i]; continue; } if( std::strcmp( argv[i], "--fileName" ) == 0 ) { fileNameEnding = argv[++i]; continue; }
if( std::strcmp( argv[i], "--magicNumber" ) == 0 ) { magicNumber = real_c(std::atof(argv[++i])); continue; } if( std::strcmp( argv[i], "--magicNumber" ) == 0 ) { magicNumber = real_c(std::atof(argv[++i])); continue; }
if( std::strcmp( argv[i], "--useOmegaBulkAdaption" ) == 0 ) { useOmegaBulkAdaption = true; continue; } if( std::strcmp( argv[i], "--useOmegaBulkAdaption" ) == 0 ) { useOmegaBulkAdaption = true; continue; }
if( std::strcmp( argv[i], "--adaptionLayerSize" ) == 0 ) { adaptionLayerSize = real_c(std::atof( argv[++i] )); continue; } if( std::strcmp( argv[i], "--adaptionLayerSize" ) == 0 ) { adaptionLayerSize = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--useSBB" ) == 0 ) { useSBB = true; continue; }
WALBERLA_ABORT("Unrecognized command line argument found: " << argv[i]); WALBERLA_ABORT("Unrecognized command line argument found: " << argv[i]);
} }
...@@ -315,6 +321,7 @@ int main( int argc, char **argv ) ...@@ -315,6 +321,7 @@ int main( int argc, char **argv )
<< " time steps = " << timesteps << "\n" << " time steps = " << timesteps << "\n"
<< " fStokes = " << fStokes << "\n" << " fStokes = " << fStokes << "\n"
<< " setup = " << setup << "\n" << " setup = " << setup << "\n"
<< " useSBB = " << useSBB << "\n"
<< "-------------------------------------------------------\n" << "-------------------------------------------------------\n"
<< " domainSize = " << xSize << " x " << ySize << " x " << zSize << "\n" << " domainSize = " << xSize << " x " << ySize << " x " << zSize << "\n"
<< " blocks = " << xBlocks << " x " << yBlocks << " x " << zBlocks << "\n" << " blocks = " << xBlocks << " x " << yBlocks << " x " << zBlocks << "\n"
...@@ -374,7 +381,7 @@ int main( int argc, char **argv ) ...@@ -374,7 +381,7 @@ int main( int argc, char **argv )
} else if (setup == 2) } else if (setup == 2)
{ {
// only tangential translational velocity // only tangential translational velocity
p.setLinearVelocity(Vector3<real_t>( real_t(0), velocity, real_t(0))); p.setLinearVelocity(Vector3<real_t>( real_t(0), real_t(0), velocity));
p.setAngularVelocity(Vector3<real_t>( real_t(0), real_t(0), real_t(0))); p.setAngularVelocity(Vector3<real_t>( real_t(0), real_t(0), real_t(0)));
} else if (setup == 3) } else if (setup == 3)
{ {
...@@ -406,7 +413,7 @@ int main( int argc, char **argv ) ...@@ -406,7 +413,7 @@ int main( int argc, char **argv )
} else if (setup == 2) } else if (setup == 2)
{ {
// only tangential translational velocity // only tangential translational velocity
p.setLinearVelocity(Vector3<real_t>( real_t(0), -velocity, real_t(0))); p.setLinearVelocity(Vector3<real_t>( real_t(0), real_t(0), -velocity));
p.setAngularVelocity(Vector3<real_t>( real_t(0), real_t(0), real_t(0))); p.setAngularVelocity(Vector3<real_t>( real_t(0), real_t(0), real_t(0)));
} else if (setup == 3) } else if (setup == 3)
{ {
...@@ -436,6 +443,7 @@ int main( int argc, char **argv ) ...@@ -436,6 +443,7 @@ int main( int argc, char **argv )
// create two planes // create two planes
mesa_pd::data::Particle&& p0 = *ps->create(true); mesa_pd::data::Particle&& p0 = *ps->create(true);
p0.setPosition(referenceVector); p0.setPosition(referenceVector);
p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(1,0,0) )); p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(1,0,0) ));
p0.setOwner(mpi::MPIManager::instance()->rank()); p0.setOwner(mpi::MPIManager::instance()->rank());
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::INFINITE);
...@@ -444,6 +452,7 @@ int main( int argc, char **argv ) ...@@ -444,6 +452,7 @@ int main( int argc, char **argv )
mesa_pd::data::Particle&& p1 = *ps->create(true); mesa_pd::data::Particle&& p1 = *ps->create(true);
p1.setPosition(Vector3<real_t>(real_c(xSize),0,0)); p1.setPosition(Vector3<real_t>(real_c(xSize),0,0));
p1.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(-1,0,0) )); p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(-1,0,0) ));
p1.setOwner(mpi::MPIManager::instance()->rank()); p1.setOwner(mpi::MPIManager::instance()->rank());
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::INFINITE);
...@@ -493,7 +502,7 @@ int main( int argc, char **argv ) ...@@ -493,7 +502,7 @@ int main( int argc, char **argv )
//////////////////////// ////////////////////////
// add omega bulk field // 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 // create the lattice model
real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega ); real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega );
...@@ -546,7 +555,13 @@ int main( int argc, char **argv ) ...@@ -546,7 +555,13 @@ int main( int argc, char **argv )
/////////////// ///////////////
// map particles into the LBM simulation // map particles into the LBM simulation
ps->forEachParticle(false, mesa_pd::kernel::SelectAll(), *accessor, movingParticleMappingKernel, *accessor, MO_Flag); if(useSBB)
{
ps->forEachParticle(false, mesa_pd::kernel::SelectAll(), *accessor, movingParticleMappingKernel, *accessor, MO_SBB_Flag);
} else
{
ps->forEachParticle(false, mesa_pd::kernel::SelectAll(), *accessor, movingParticleMappingKernel, *accessor, MO_CLI_Flag);
}
// create the timeloop // create the timeloop
SweepTimeloop timeloop( blocks->getBlockStorage(), timesteps ); SweepTimeloop timeloop( blocks->getBlockStorage(), timesteps );
...@@ -637,7 +652,7 @@ int main( int argc, char **argv ) ...@@ -637,7 +652,7 @@ int main( int argc, char **argv )
real_t curTorqueNorm = real_t(0); real_t curTorqueNorm = real_t(0);
real_t oldTorqueNorm = real_t(0); real_t oldTorqueNorm = real_t(0);
real_t convergenceLimit = real_t(1e-4); real_t convergenceLimit = real_t(1e-5);
// time loop // time loop
for (uint_t i = 1; i <= timesteps; ++i ) for (uint_t i = 1; i <= timesteps; ++i )
...@@ -690,18 +705,10 @@ int main( int argc, char **argv ) ...@@ -690,18 +705,10 @@ int main( int argc, char **argv )
WALBERLA_MPI_SECTION() WALBERLA_MPI_SECTION()
{ {
mpi::allReduceInplace( hydForce[0], mpi::SUM ); mpi::allReduceInplace( hydForce, mpi::SUM );
mpi::allReduceInplace( hydForce[1], mpi::SUM ); mpi::reduceInplace( lubForce, mpi::SUM );
mpi::allReduceInplace( hydForce[2], mpi::SUM ); mpi::allReduceInplace( hydTorque, mpi::SUM );
mpi::reduceInplace( lubForce[0], mpi::SUM ); mpi::reduceInplace( lubTorque, mpi::SUM );
mpi::reduceInplace( lubForce[1], mpi::SUM );
mpi::reduceInplace( lubForce[2], mpi::SUM );
mpi::allReduceInplace( hydTorque[0], mpi::SUM );
mpi::allReduceInplace( hydTorque[1], mpi::SUM );
mpi::allReduceInplace( hydTorque[2], mpi::SUM );
mpi::reduceInplace( lubTorque[0], mpi::SUM );
mpi::reduceInplace( lubTorque[1], mpi::SUM );
mpi::reduceInplace( lubTorque[2], mpi::SUM );
} }
curForceNorm = hydForce.length(); curForceNorm = hydForce.length();
...@@ -747,6 +754,7 @@ int main( int argc, char **argv ) ...@@ -747,6 +754,7 @@ int main( int argc, char **argv )
loggingFileName += "_bvrf" + std::to_string(uint_c(bulkViscRateFactor)); loggingFileName += "_bvrf" + std::to_string(uint_c(bulkViscRateFactor));
loggingFileName += "_mn" + std::to_string(float(magicNumber)); loggingFileName += "_mn" + std::to_string(float(magicNumber));
if( useOmegaBulkAdaption ) loggingFileName += "_uOBA" + std::to_string(uint_c(adaptionLayerSize)); if( useOmegaBulkAdaption ) loggingFileName += "_uOBA" + std::to_string(uint_c(adaptionLayerSize));
if( useSBB ) loggingFileName += "_SBB";
if( !fileNameEnding.empty()) loggingFileName += "_" + fileNameEnding; if( !fileNameEnding.empty()) loggingFileName += "_" + fileNameEnding;
loggingFileName += ".txt"; loggingFileName += ".txt";
......
apps/benchmarks/FluidParticleCoupling/MotionSettlingSphere.cpp 0 → 100644
View file @ 7a6d3b57
//======================================================================================================================
//
// 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 MotionSettlingSphere.cpp
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#include "boundary/all.h"
#include "blockforest/all.h"
#include "core/all.h"
#include "domain_decomposition/all.h"
#include "field/AddToStorage.h"
#include "field/communication/PackInfo.h"
#include "field/vtk/all.h"
#include "lbm/boundary/all.h"
#include "lbm/communication/PdfFieldPackInfo.h"
#include "lbm/field/AddToStorage.h"
#include "lbm/field/PdfField.h"
#include "lbm/lattice_model/D3Q19.h"
#include "lbm/sweeps/CellwiseSweep.h"
#include "lbm/sweeps/SweepWrappers.h"
#include "lbm/vtk/all.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/SimpleBB.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/ExtrapolationDirectionFinder.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/OmegaBulkAdaption.h"
#include "lbm_mesapd_coupling/utility/InspectionProbe.h"
#include "lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h"
#include "mesa_pd/data/ParticleStorage.h"
#include "mesa_pd/data/ParticleAccessorWithShape.h"
#include "mesa_pd/data/ShapeStorage.h"
#include "mesa_pd/data/DataTypes.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/ParticleSelector.h"
#include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/SyncGhostOwners.h"
#include "mesa_pd/mpi/ReduceProperty.h"
#include "mesa_pd/mpi/notifications/ForceTorqueNotification.h"
#include "mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h"
#include "mesa_pd/vtk/ParticleVtkOutput.h"
#include "stencil/D3Q27.h"
#include "timeloop/SweepTimeloop.h"
#include "vtk/BlockCellDataWriter.h"
#include "vtk/Initialization.h"
#include "vtk/VTKOutput.h"
#include <functional>
#include <memory>
#ifdef WALBERLA_BUILD_WITH_CODEGEN
#include "GeneratedLBM.h"
#endif
namespace motion_settling_sphere{
///////////
// USING //
///////////
using namespace walberla;
using walberla::uint_t;
//////////////
// TYPEDEFS //
//////////////
using namespace walberla;
using walberla::uint_t;
#ifdef WALBERLA_BUILD_WITH_CODEGEN
using LatticeModel_T = lbm::GeneratedLBM;
#else
using LatticeModel_T = lbm::D3Q19< lbm::collision_model::D3Q19MRT>;
#endif
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 ScalarField_T = GhostLayerField< real_t, 1>;
const uint_t FieldGhostLayers = 1;
///////////
// FLAGS //
///////////
const FlagUID Fluid_Flag ( "fluid" );
const FlagUID UBB_Flag ( "velocity bounce back" );
const FlagUID Outlet_Flag ( "outlet" );
const FlagUID MEM_BB_Flag ( "moving obstacle BB" );
const FlagUID MEM_CLI_Flag ( "moving obstacle CLI" );
const FlagUID FormerMEM_Flag( "former moving obstacle" );
//////////////////////////////
// Coupling algorithm enums //
//////////////////////////////
enum MEMVariant { BB, CLI, MR };
MEMVariant to_MEMVariant( const std::string& s )
{
if( s == "BB" ) return MEMVariant::BB;
if( s == "CLI" ) return MEMVariant::CLI;
throw std::runtime_error("invalid conversion from text to MEMVariant");
}
std::string MEMVariant_to_string ( const MEMVariant& m )
{
if( m == MEMVariant::BB ) return "BB";
if( m == MEMVariant::CLI ) return "CLI";
throw std::runtime_error("invalid conversion from MEMVariant to string");
}
/////////////////////////////////////
// BOUNDARY HANDLING CUSTOMIZATION //
/////////////////////////////////////
template <typename ParticleAccessor_T>
class MyBoundaryHandling
{
public:
using UBB_T = lbm::SimpleUBB< LatticeModel_T, flag_t >;
using Outlet_T = lbm::SimplePressure< LatticeModel_T, flag_t >;
using MO_SBB_T = lbm_mesapd_coupling::SimpleBB< LatticeModel_T, FlagField_T, ParticleAccessor_T >;
using MO_CLI_T = lbm_mesapd_coupling::CurvedLinear< LatticeModel_T, FlagField_T, ParticleAccessor_T >;
using Type = BoundaryHandling< FlagField_T, Stencil_T, UBB_T, Outlet_T, MO_SBB_T, MO_CLI_T >;
MyBoundaryHandling( const BlockDataID & flagFieldID, const BlockDataID & pdfFieldID,
const BlockDataID & particleFieldID, const shared_ptr<ParticleAccessor_T>& ac,
Vector3<real_t> uInfty) :
flagFieldID_( flagFieldID ), pdfFieldID_( pdfFieldID ), particleFieldID_( particleFieldID ),
ac_( ac ), velocity_( uInfty )
{}
Type * operator()( IBlock* const block, const StructuredBlockStorage* const storage ) const
{
WALBERLA_ASSERT_NOT_NULLPTR( block );
WALBERLA_ASSERT_NOT_NULLPTR( storage );
FlagField_T * flagField = block->getData< FlagField_T >( flagFieldID_ );
PdfField_T * 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 );
Type * handling = new Type( "moving obstacle boundary handling", flagField, fluid,
UBB_T( "UBB", UBB_Flag, pdfField, velocity_),
Outlet_T( "Outlet", Outlet_Flag, pdfField, real_t(1) ),
MO_SBB_T ( "MEM_BB", MEM_BB_Flag, pdfField, flagField, particleField, ac_, fluid, *storage, *block ),
MO_CLI_T( "MEM_CLI", MEM_CLI_Flag, pdfField, flagField, particleField, ac_, fluid, *storage, *block ) );
const auto ubb = flagField->getFlag( UBB_Flag );
const auto outlet = flagField->getFlag( Outlet_Flag );
CellInterval domainBB = storage->getDomainCellBB();
domainBB.xMin() -= cell_idx_c( FieldGhostLayers ); // -1
domainBB.xMax() += cell_idx_c( FieldGhostLayers ); // cellsX
domainBB.yMin() -= cell_idx_c( FieldGhostLayers ); // 0
domainBB.yMax() += cell_idx_c( FieldGhostLayers ); // cellsY+1
domainBB.zMin() -= cell_idx_c( FieldGhostLayers ); // 0
domainBB.zMax() += cell_idx_c( FieldGhostLayers ); // cellsZ+1
// BOTTOM
// -1........-1........-1
// cellsX+1..cellsY+1..-1
CellInterval bottom( domainBB.xMin(), domainBB.yMin(), domainBB.zMin(), domainBB.xMax(), domainBB.yMax(), domainBB.zMin() );
storage->transformGlobalToBlockLocalCellInterval( bottom, *block );
handling->forceBoundary( ubb, bottom );
// TOP
// -1........-1........cellsZ+1
// cellsX+1..cellsY+1..cellsZ+1
CellInterval top( domainBB.xMin(), domainBB.yMin(), domainBB.zMax(), domainBB.xMax(), domainBB.yMax(), domainBB.zMax() );
storage->transformGlobalToBlockLocalCellInterval( top, *block );
handling->forceBoundary( outlet, top );
handling->fillWithDomain( FieldGhostLayers );
return handling;
}
private:
const BlockDataID flagFieldID_;
const BlockDataID pdfFieldID_;
const BlockDataID particleFieldID_;
shared_ptr<ParticleAccessor_T> ac_;
const Vector3<real_t> velocity_;
}; // class MyBoundaryHandling
/*
* Functionality for continuous logging of sphere properties during initial (resting) simulation
*/
template< typename ParticleAccessor_T>
class RestingSphereForceEvaluator
{
public:
RestingSphereForceEvaluator( const shared_ptr< ParticleAccessor_T > & ac, walberla::id_t sphereUid,
uint_t averageFrequency, const std::string & basefolder ) :
ac_( ac ), sphereUid_( sphereUid ), averageFrequency_( averageFrequency )
{
WALBERLA_ROOT_SECTION() {
// initially write header in file
std::ofstream file;
filename_ = basefolder;
filename_ += "/log_init.txt";
file.open(filename_.c_str());
file << "# f_z_current f_z_average f_x f_y\n";
file.close();
}
}
// evaluate the acting drag force on a fixed sphere
void operator()(const uint_t timestep)
{
// average the force over averageFrequency consecutive timesteps since there are small fluctuations in the force
auto currentForce = getForce();
currentAverage_ += currentForce[2];
WALBERLA_ROOT_SECTION()
{
std::ofstream file;
file.open( filename_.c_str(), std::ofstream::app );
file.precision(8);
file << timestep << "\t" << currentForce[2] << "\t" << dragForceNew_
<< "\t" << currentForce[0] << "\t" << currentForce[1] << std::endl;
file.close();
}
if( timestep % averageFrequency_ != 0) return;
dragForceOld_ = dragForceNew_;
dragForceNew_ = currentAverage_ / real_c( averageFrequency_ );
currentAverage_ = real_t(0);
}
// Return the relative temporal change in the drag force
real_t getForceDiff() const
{
return std::fabs( ( dragForceNew_ - dragForceOld_ ) / dragForceNew_ );
}
real_t getDragForce() const
{
return dragForceNew_;
}
private:
// Obtain the drag force acting on the sphere
Vector3<real_t> getForce()
{
Vector3<real_t> force(real_t( 0 ));
size_t idx = ac_->uidToIdx(sphereUid_);
if( idx != ac_->getInvalidIdx())
{
if(!mesa_pd::data::particle_flags::isSet( ac_->getFlags(idx), mesa_pd::data::particle_flags::GHOST))
{
force = ac_->getHydrodynamicForce(idx);
}
}
WALBERLA_MPI_SECTION()
{
mpi::allReduceInplace( force, mpi::SUM );
}
return force;
}
shared_ptr< ParticleAccessor_T > ac_;
const walberla::id_t sphereUid_;
real_t currentAverage_ = real_t(0);
uint_t averageFrequency_;
std::string filename_;
real_t dragForceOld_ = real_t(0);
real_t dragForceNew_ = real_t(0);
};
/*
* Functionality for continuous logging of sphere properties during moving simulation
*/
template< typename ParticleAccessor_T>
class MovingSpherePropertyEvaluator
{
public:
MovingSpherePropertyEvaluator( const shared_ptr< ParticleAccessor_T > & ac, walberla::id_t sphereUid, const std::string & basefolder,
const Vector3<real_t> & u_infty, const real_t Galileo, const real_t GalileoSim,
const real_t gravity, const real_t viscosity, const real_t diameter,
const real_t densityRatio ) :
ac_( ac ), sphereUid_( sphereUid ),
u_infty_( u_infty ), gravity_( gravity ), viscosity_( viscosity ),
diameter_( diameter ), densityRatio_( densityRatio )
{
WALBERLA_ROOT_SECTION()
{
std::ofstream fileSetup;
std::string filenameSetup = basefolder;
filenameSetup += "/setup.txt";
fileSetup.open( filenameSetup.c_str() );
fileSetup << "# setup parameters: \n";
fileSetup << "processes = " << MPIManager::instance()->numProcesses() << "\n";
fileSetup << "Galileo number (targeted) = " << Galileo << "\n";
fileSetup << "Galileo number (simulated) = " << GalileoSim << "\n";
fileSetup << "gravity = " << gravity << "\n";
fileSetup << "viscosity = " << viscosity << "\n";
fileSetup << "diameter = " << diameter << "\n";
fileSetup << "uInfty = " << u_infty_[2] << "\n";
real_t u_ref = std::sqrt( std::fabs(densityRatio - real_t(1)) * gravity * diameter );
fileSetup << "u_{ref} = " << u_ref << "\n";
fileSetup.close();
}
filenameRes_ = basefolder;
filenameRes_ += "/log_results.txt";
filenameAdd_ = basefolder;
filenameAdd_ += "/log_additional.txt";
WALBERLA_ROOT_SECTION()
{
// initially write header in file
std::ofstream fileRes;
fileRes.open( filenameRes_.c_str() );
// raw data
fileRes << "#\t x_p_x\t x_p_y\t x_p_z\t u_p_x\t u_p_y\t u_p_z\t omega_p_x\t omega_p_y\t omega_p_z\t u_{pr}_x\t u_{pr}_y\t u_{pr}_z\t ";
// processed data, 14 - 18
fileRes << "Re\t u_{pV}\t u_{pH}\t omega_{pV}\t omega_{pH}\t alpha\n";
fileRes.close();
std::ofstream fileAdd;
fileAdd.open( filenameAdd_.c_str() );
fileAdd << "# forceX forceY forceZ torqueX torqueY torqueZ\n";
fileAdd.close();
}
}
// evaluate and write the sphere properties
void operator()(const uint_t timestep)
{
Vector3<real_t> particleTransVel( real_t(0) );
Vector3<real_t> particleAngularVel( real_t(0) );
Vector3<real_t> particlePos( real_t(0) );
Vector3<real_t> particleForce( real_t(0) );
Vector3<real_t> particleTorque( real_t(0) );
size_t idx = ac_->uidToIdx(sphereUid_);
if( idx != ac_->getInvalidIdx())
{
if(!mesa_pd::data::particle_flags::isSet( ac_->getFlags(idx), mesa_pd::data::particle_flags::GHOST))
{
particlePos = ac_->getPosition(idx);
particleTransVel = ac_->getLinearVelocity(idx);
particleAngularVel = ac_->getAngularVelocity(idx);
particleForce = ac_->getHydrodynamicForce(idx);
particleTorque = ac_->getHydrodynamicTorque(idx);
}
}
// reduce to root
WALBERLA_MPI_SECTION()
{
mpi::allReduceInplace( particlePos, mpi::SUM );
mpi::reduceInplace( particleTransVel, mpi::SUM );
mpi::reduceInplace( particleAngularVel, mpi::SUM );
mpi::reduceInplace( particleForce, mpi::SUM );
mpi::reduceInplace( particleTorque, mpi::SUM );
}
particleHeight_ = particlePos[2];
//only root process has all the data
WALBERLA_ROOT_SECTION()
{
Vector3<real_t> u_p_r( particleTransVel[0] - u_infty_[0], particleTransVel[1] - u_infty_[1], particleTransVel[2] - u_infty_[2]);
real_t u_p_H = std::sqrt( u_p_r[0] * u_p_r[0] + u_p_r[1] * u_p_r[1]);
real_t u_p_V = u_p_r[2];
real_t omega_p_H = std::sqrt( particleAngularVel[0] * particleAngularVel[0] + particleAngularVel[1] * particleAngularVel[1] );
real_t omega_p_V = particleAngularVel[2];
real_t u_ref = std::sqrt( std::fabs(densityRatio_ - real_t(1)) * gravity_ * diameter_ );
real_t Reynolds = u_p_r.length() * diameter_ / viscosity_;
// results
std::ofstream fileRes;
fileRes.open( filenameRes_.c_str(), std::ofstream::app );
fileRes.precision(8);
fileRes << timestep
<< "\t" << particlePos[0] << "\t" << particlePos[1] << "\t" << particlePos[2]
<< "\t" << particleTransVel[0] << "\t" << particleTransVel[1] << "\t" << particleTransVel[2]
<< "\t" << particleAngularVel[0] << "\t" << particleAngularVel[1] << "\t" << particleAngularVel[2]
<< "\t" << u_p_r[0] << "\t" << u_p_r[1] << "\t" << u_p_r[2]
<< "\t" << Reynolds << "\t" << u_p_V/u_ref << "\t" << u_p_H/u_ref
<< "\t" << omega_p_V*(diameter_/u_ref) << "\t" << omega_p_H*(diameter_/u_ref)
<< "\t" << std::atan( u_p_H/std::fabs(u_p_V) )
<< std::endl;
fileRes.close();
// additional
std::ofstream fileAdd;
fileAdd.open( filenameAdd_.c_str(), std::ofstream::app );
fileAdd.precision(8);
fileAdd << timestep
<< "\t" << particleForce[0] << "\t" << particleForce[1] << "\t" << particleForce[2]
<< "\t" << particleTorque[0] << "\t" << particleTorque[1] << "\t" << particleTorque[2]
<< std::endl;
fileAdd.close();
}
}
real_t getParticleHeight() const
{
return particleHeight_;
}
private:
shared_ptr< ParticleAccessor_T > ac_;
const walberla::id_t sphereUid_;
std::string filenameRes_;
std::string filenameAdd_;
Vector3<real_t> u_infty_;
real_t gravity_;
real_t viscosity_;
real_t diameter_;
real_t densityRatio_;
real_t particleHeight_;
};
/*
* Result evaluation for the written VTK files
*
* input: vel (fluid velocity), u_p (particle velocity), u_infty (inflow velocity)
*
* needed data:
* relative flow velocity: vel_r = vel - u_p
* relative particle velocity: u_p_r = u_p - u_infty
*
* magnitude of relative particle velocity in horizontal plane: u_p_H = sqrt( (u_p_r)_x^2 + (u_p_r)_y^2 )
* unit vector of particle motion in horizontal plane: e_p_H = ( (u_p_r)_x, (u_p_r)_y, 0) / u_p_H
* unit vector perpendicular to e_p_H and e_z: e_p_Hz_perp = ( -(u_p_r)_y, (u_p_r)_x, 0) / u_p_H
* unit vector of particle motion: e_p_parallel = u_p_r / ||u_p_r||
* unit vector perpendicular to e_p_parallel and e_p_Hz_perp: e_p_perp = e_p_Hz_perp x e_p_parallel
* projected fluid velocities: vel_r_parallel = vel_r * (-e_p_parallel)
* vel_r_perp = vel_r * e_p_perp
* vel_r_Hz_perp = vel_r * e_p_Hz_perp
*/
template< typename ParticleAccessor_T>
class VTKInfoLogger
{
public:
VTKInfoLogger( SweepTimeloop* timeloop, const shared_ptr< ParticleAccessor_T > & ac, walberla::id_t sphereUid,
const std::string & baseFolder, const Vector3<real_t>& u_infty ) :
timeloop_( timeloop ), ac_( ac ), sphereUid_( sphereUid ), baseFolder_( baseFolder ), u_infty_( u_infty )
{ }
void operator()()
{
Vector3<real_t> u_p( real_t(0) );
size_t idx = ac_->uidToIdx(sphereUid_);
if( idx != ac_->getInvalidIdx())
{
if(!mesa_pd::data::particle_flags::isSet( ac_->getFlags(idx), mesa_pd::data::particle_flags::GHOST))
{
u_p = ac_->getLinearVelocity(idx);
}
}
WALBERLA_MPI_SECTION()
{
mpi::allReduceInplace( u_p, mpi::SUM );
}
Vector3<real_t> u_p_r = u_p - u_infty_;
real_t u_p_r_sqr_mag = u_p_r.sqrLength();
real_t u_p_H = std::sqrt( u_p_r[0] * u_p_r[0] + u_p_r[1] * u_p_r[1]);
Vector3<real_t> e_p_H (u_p_r[0], u_p_r[1], real_t(0));
e_p_H /= u_p_H;
Vector3<real_t> e_p_Hz_perp (-u_p_r[1], u_p_r[0], real_t(0));
e_p_Hz_perp /= u_p_H;
Vector3<real_t> e_p_parallel = u_p_r / u_p_r.length();
Vector3<real_t> e_p_perp = e_p_Hz_perp%e_p_parallel;
WALBERLA_ROOT_SECTION()
{
std::ofstream file;
std::string filename(baseFolder_);
filename += "/log_vtk/log_vtk_";
filename += std::to_string( timeloop_->getCurrentTimeStep() );
filename += ".txt";
file.open( filename.c_str() );
file.precision(8);
file << "u_p = " << u_p << "\n";
file << "u_p_r_2 = " << u_p_r_sqr_mag << "\n\n";
file << "e_{pH} = " << e_p_H << "\n";
file << "e_{pparallel} = " << e_p_parallel << "\n";
file << "e_{pperp} = " << e_p_perp << "\n";
file << "e_{pHzperp} = " << e_p_Hz_perp << "\n";
file.close();
}
}
private:
SweepTimeloop* timeloop_;
shared_ptr< ParticleAccessor_T > ac_;
const walberla::id_t sphereUid_;
std::string baseFolder_;
Vector3<real_t> u_infty_;
};
/*
* This extensive, physical test case simulates a single, heavy sphere in ambient fluid flow.
* It is based on the benchmark proposed in
* Uhlmann, Dusek - "The motion of a single heavy sphere in ambient fluid: A benchmark for interface-resolved
* particulate flow simulations with significant relative velocities" IJMF (2014),
* doi: 10.1016/j.ijmultiphaseflow.2013.10.010
* Results for LBM done with waLBerla are published in
* Rettinger, Ruede - "A comparative study of fluid-particle coupling methods for fully resolved
* lattice Boltzmann simulations" CaF (2017),
* doi: 10.1016/j.compfluid.2017.05.033
* The setup and the benchmark itself are explained in detail in these two papers.
* Several logging files are written to evaluate the benchmark quantities.
*
* This case is the same as can be found in apps/benchmarks/MotionSingleHeavySphere
* but using the lbm_mesapd_coupling, instead of the pe_coupling.
* Compared to this previous work, significantly different outcomes can be observed when using the CLI boundary condition.
* This can be traced back to the distinct way of force averaging.
* The force averaging applied in the previous study (two LBM steps, then average) introduces a error in the Galilean invariance.
* Surprisingly, this error can result in better agreement to the reference solutions.
*
*/
int main( int argc, char **argv )
{
debug::enterTestMode();
mpi::Environment env( argc, argv );
bool noViscosityIterations = false;
bool vtkIOInit = false;
bool longDom = false;
bool broadDom = false;
bool useOmegaBulkAdaption = false;
real_t bulkViscRateFactor = real_t(1);
uint_t averageForceTorqueOverTwoTimeStepsMode = 1; // 0: no, 1: running, 2: 2LBM
bool conserveMomentum = false;
bool useDiffusiveScaling = true;
bool useVelocityVerlet = true;
bool initFluidVelocity = true;
bool vtkOutputGhostlayers = false;
MEMVariant memVariant = MEMVariant::BB;
std::string reconstructorType = "Grad"; // Eq, EAN, Ext, Grad
real_t diameter = real_t(18);
uint_t XBlocks = uint_t(4);
uint_t YBlocks = uint_t(4);
real_t viscosity = real_t(0.01); // used in diffusive scaling
real_t uIn = real_t(0.02); // used for acoustic scaling
/*
* 0: custom
* 1: Ga = 144
* 2: Ga = 178.46
* 3: Ga = 190
* 4: Ga = 250
*/
uint_t simulationCase = 1;
real_t Galileo = real_t(1);
real_t Re_target = real_t(1);
real_t timestepsNonDim = real_t(1);
real_t vtkWriteSpacingNonDim = real_t(3); // write every 3 non-dim timesteps
uint_t vtkNumOutputFiles = 10; // write only 10 vtk files: the 10 final ones
std::string folderEnding = "";
////////////////////////////
// COMMAND LINE ARGUMENTS //
////////////////////////////
for( int i = 1; i < argc; ++i )
{
if( std::strcmp( argv[i], "--case" ) == 0 ) {simulationCase = uint_c( std::atof( argv[++i] ) ); continue;}
if( std::strcmp( argv[i], "--Ga" ) == 0 ){Galileo = real_c( std::atof( argv[++i] ) ); continue;}
if( std::strcmp( argv[i], "--Re" ) == 0 ){Re_target = 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], "--diameter" ) == 0 ){diameter = real_c( std::atof( argv[++i] ) ); continue;}
if( std::strcmp( argv[i], "--noViscosityIterations" ) == 0 ){noViscosityIterations = true; continue;}
if( std::strcmp( argv[i], "--viscosity" ) == 0 ){viscosity = real_c( std::atof( argv[++i] ) ); continue;}
if( std::strcmp( argv[i], "--uIn" ) == 0 ){uIn = real_c( std::atof( argv[++i] ) ); continue;}
if( std::strcmp( argv[i], "--vtkIOInit" ) == 0 ){vtkIOInit = true; continue;}
if( std::strcmp( argv[i], "--longDom" ) == 0 ){longDom = true; continue;}
if( std::strcmp( argv[i], "--broadDom" ) == 0 ){broadDom = true; continue;}
if( std::strcmp( argv[i], "--variant" ) == 0 ){memVariant = to_MEMVariant( argv[++i] ); continue;}
if( std::strcmp( argv[i], "--useOmegaBulkAdaption" ) == 0 ){useOmegaBulkAdaption = true; continue;}
if( std::strcmp( argv[i], "--bvrf" ) == 0 ){bulkViscRateFactor = real_c( std::atof( argv[++i] ) ); continue;}
if( std::strcmp( argv[i], "--forceTorqueAverageMode" ) == 0 ){averageForceTorqueOverTwoTimeStepsMode = uint_c( std::atof( argv[++i] ) ); continue;}
if( std::strcmp( argv[i], "--conserveMomentum" ) == 0 ){conserveMomentum = true; continue;}
if( std::strcmp( argv[i], "--useDiffusiveScaling" ) == 0 ){useDiffusiveScaling = true; 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], "--folderEnding" ) == 0 ){folderEnding = argv[++i]; continue;}
if( std::strcmp( argv[i], "--reconstructorType" ) == 0 ) {reconstructorType = argv[++i]; continue;}
if( std::strcmp( argv[i], "--useEuler" ) == 0 ) {useVelocityVerlet = false; continue;}
if( std::strcmp( argv[i], "--noFluidVelocityInit" ) == 0 ) {initFluidVelocity = false; continue;}
if( std::strcmp( argv[i], "--vtkOutputGhostlayers" ) == 0 ) {vtkOutputGhostlayers = true; continue;}
if( std::strcmp( argv[i], "--vtkNumOutputFiles" ) == 0 ){vtkNumOutputFiles = uint_c( std::atof( argv[++i] ) ); continue;}
if( std::strcmp( argv[i], "--vtkWriteSpacingNonDim" ) == 0 ){vtkWriteSpacingNonDim = real_c( std::atof( argv[++i] ) ); continue;}
WALBERLA_ABORT("command line argument unknown: " << argv[i] );
}
///////////////////////////
// SIMULATION PROPERTIES //
///////////////////////////
const real_t radius = real_t(0.5) * diameter;
uint_t xlength = uint_c( diameter * real_t(5.34) );
uint_t ylength = xlength;
uint_t zlength = uint_c( diameter * real_t(16) );
if(broadDom)
{
xlength *= uint_t(2);
ylength *= uint_t(2);
}
if(longDom)
{
zlength *= uint_t(3);
}
// estimate Reynolds number (i.e. inflow velocity) based on Galileo number
// values are taken from the original simulation of Uhlmann, Dusek
switch( simulationCase )
{
case 0:
WALBERLA_LOG_INFO_ON_ROOT("Using custom simulation case -> you have to provide Ga, Re, and time steps via command line arguments!");
break;
case 1:
Galileo = real_t(144);
Re_target = real_t(185.08);
timestepsNonDim = real_t(100);
break;
case 2:
Galileo = real_t(178.46);
Re_target = real_t(243.01);
timestepsNonDim = real_t(250);
break;
case 3:
Galileo = real_t(190);
Re_target = real_t(262.71);
timestepsNonDim = real_t(250);
break;
case 4:
Galileo = real_t(250);
Re_target = real_t(365.10);
timestepsNonDim = real_t(510);
break;
default:
WALBERLA_ABORT("Simulation case is not supported!");
}
if( useDiffusiveScaling)
{
// estimate fitting inflow velocity (diffusive scaling, viscosity is fixed)
// attention: might lead to large velocities for small diameters
uIn = Re_target * viscosity / diameter;
} else {
// estimate fitting viscosity (acoustic scaling: velocity is fixed)
// note that this is different to the approach taking in the paper, where an diffusive scaling with visc = 0.01 is taken
// attention: might lead to small viscosities for small diameters
viscosity = uIn * diameter / Re_target;
}
real_t omega = lbm::collision_model::omegaFromViscosity(viscosity);
real_t omegaBulk = lbm_mesapd_coupling::omegaBulkFromOmega(omega, bulkViscRateFactor);
Vector3<real_t> uInfty = Vector3<real_t>( real_t(0), real_t(0), uIn );
const real_t densityRatio = real_t(1.5);
const uint_t averageFrequency = uint_c( ( ( uint_c(Galileo) >= 200) ? real_t(500) : real_t(2) ) * diameter / uIn ); // for initial simulation
const real_t convergenceLimit = real_t(1e-4);
const real_t convergenceLimitGalileo = real_t(1e-4);
const real_t dx = real_t(1);
const real_t magicNumberTRT = lbm::collision_model::TRT::threeSixteenth;
const uint_t timestepsInit = uint_c( ( ( uint_c(Galileo) >= 200) ? real_t(3000) : real_t(100) ) * diameter / uIn ); // maximum number of time steps for the initial simulation
const uint_t writeFrequencyInit = uint_t(1000); // vtk write frequency init
///////////////////////////
// BLOCK STRUCTURE SETUP //
///////////////////////////
const int numProcs = MPIManager::instance()->numProcesses();
uint_t ZBlocks = uint_c(numProcs) / ( XBlocks * YBlocks );
const uint_t XCells = xlength / XBlocks;
const uint_t YCells = ylength / YBlocks;
const uint_t ZCells = zlength / ZBlocks;
if( (xlength != XCells * XBlocks) || (ylength != YCells * YBlocks) || (zlength != ZCells * ZBlocks) )
{
WALBERLA_ABORT("Domain partitioning does not fit to total domain size!");
}
WALBERLA_LOG_INFO_ON_ROOT("Motion settling sphere simulation of case " << simulationCase << " -> Ga = " << Galileo);
WALBERLA_LOG_INFO_ON_ROOT("Diameter: " << diameter);
WALBERLA_LOG_INFO_ON_ROOT("Domain: " << xlength << " x " << ylength << " x " << zlength);
WALBERLA_LOG_INFO_ON_ROOT("Processes: " << XBlocks << " x " << YBlocks << " x " << ZBlocks);
WALBERLA_LOG_INFO_ON_ROOT("Subdomains: " << XCells << " x " << YCells << " x " << ZCells);
WALBERLA_LOG_INFO_ON_ROOT("Tau: " << real_t(1)/omega);
WALBERLA_LOG_INFO_ON_ROOT("Viscosity: " << viscosity);
WALBERLA_LOG_INFO_ON_ROOT("uIn: " << uIn);
WALBERLA_LOG_INFO_ON_ROOT("Re_infty: " << uIn * diameter / viscosity);
auto blocks = blockforest::createUniformBlockGrid( XBlocks, YBlocks, ZBlocks, XCells, YCells, ZCells, dx, true,
true, true, false );
//////////////////
// 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);
real_t xParticle = real_t(0);
real_t yParticle = real_t(0);
real_t zParticle = real_t(0);
// root determines particle position, then broadcasts it
WALBERLA_ROOT_SECTION()
{
if( simulationCase == 1 )
{
xParticle = real_c( xlength ) * real_t(0.5);
yParticle = real_c( ylength ) * real_t(0.5);
}
else if( simulationCase == 4 )
{
// add random perturbance for chaotic regime
walberla::math::seedRandomGenerator( std::mt19937::result_type(std::time(nullptr)) );
xParticle = real_c( xlength ) * real_t(0.5) + walberla::math::realRandom( real_t(-0.5), real_t(0.5) );
yParticle = real_c( ylength ) * real_t(0.5) + walberla::math::realRandom( real_t(-0.5), real_t(0.5) );
}
else
{
//add small perturbance to sphere position to break stability due to numerical symmetry
real_t perturbance = real_t(0.35);
xParticle = real_c( xlength ) * real_t(0.5) + perturbance;
yParticle = real_c( ylength ) * real_t(0.5);
}
zParticle = (longDom) ? ( diameter * real_t(16) + real_c( xlength ) ) : real_c( xlength );
}
// broadcast to other ranks
WALBERLA_MPI_SECTION()
{
mpi::broadcastObject( xParticle );
mpi::broadcastObject( yParticle );
mpi::broadcastObject( zParticle );
}
// add sphere
Vector3<real_t> position( xParticle, yParticle, zParticle );
auto sphereShape = ss->create<mesa_pd::data::Sphere>( radius );
ss->shapes[sphereShape]->updateMassAndInertia(densityRatio);
walberla::id_t sphereUid = 0;
if (rpdDomain->isContainedInProcessSubdomain( uint_c(mpi::MPIManager::instance()->rank()), position ))
{
mesa_pd::data::Particle&& p = *ps->create();
p.setPosition(position);
p.setInteractionRadius(radius);
p.setOwner(mpi::MPIManager::instance()->rank());
p.setShapeID(sphereShape);
sphereUid = p.getUid();
}
mpi::allReduceInplace(sphereUid, mpi::SUM);
// set up synchronization procedure
const real_t overlap = real_t( 1.5 ) * dx;
std::function<void(void)> syncCall;
if( XBlocks <= uint_t(4) )
{
WALBERLA_LOG_INFO_ON_ROOT("Using next neighbor sync!")
syncCall = [&ps,&rpdDomain,overlap](){
mesa_pd::mpi::SyncNextNeighbors syncNextNeighborFunc;
syncNextNeighborFunc(*ps, *rpdDomain, overlap);
};
syncCall();
}
else
{
WALBERLA_LOG_INFO_ON_ROOT("Using ghost owner sync!")
syncCall = [&ps,&rpdDomain,overlap](){
mesa_pd::mpi::SyncGhostOwners syncGhostOwnersFunc;
syncGhostOwnersFunc(*ps, *rpdDomain, overlap);
};
for(uint_t i = 0; i < uint_c(XBlocks/2); ++i) syncCall();
}
WALBERLA_LOG_INFO_ON_ROOT("Initial sphere position: " << position);
///////////////////////
// ADD DATA TO BLOCKS //
////////////////////////
// add omega bulk field
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, magicNumberTRT );
#ifdef WALBERLA_BUILD_WITH_CODEGEN
WALBERLA_LOG_INFO_ON_ROOT("Using generated TRT-like lattice model!");
LatticeModel_T latticeModel = LatticeModel_T(omegaBulkFieldID, lambda_d, lambda_e);
#else
WALBERLA_LOG_INFO_ON_ROOT("Using waLBerla built-in TRT lattice model and ignoring omega bulk!");
LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( omega, magicNumberTRT ) );
#endif
// add PDF field
// initial velocity in domain = inflow velocity
BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field", latticeModel, initFluidVelocity ? uInfty : Vector3<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, uInfty), "boundary handling" );
// kernels
mesa_pd::mpi::ReduceProperty reduceProperty;
lbm_mesapd_coupling::AddHydrodynamicInteractionKernel addHydrodynamicInteraction;
lbm_mesapd_coupling::ResetHydrodynamicForceTorqueKernel resetHydrodynamicForceTorque;
lbm_mesapd_coupling::AverageHydrodynamicForceTorqueKernel averageHydrodynamicForceTorque;
lbm_mesapd_coupling::InitializeHydrodynamicForceTorqueForAveragingKernel initializeHydrodynamicForceTorqueForAveragingKernel;
lbm_mesapd_coupling::RegularParticlesSelector sphereSelector;
lbm_mesapd_coupling::MovingParticleMappingKernel<BoundaryHandling_T> movingParticleMappingKernel(blocks, boundaryHandlingID, particleFieldID);
// mapping of sphere required by MEM variants
// sets the correct flags
if( memVariant == MEMVariant::CLI )
{
ps->forEachParticle(false, sphereSelector, *accessor, movingParticleMappingKernel, *accessor, MEM_CLI_Flag);
}else {
ps->forEachParticle(false, sphereSelector, *accessor, movingParticleMappingKernel, *accessor, MEM_BB_Flag);
}
// base folder to store all logs and vtk output
std::string basefolder ("MSHS_");
basefolder += std::to_string( uint_c( Galileo ) );
basefolder += "_";
basefolder += std::to_string( uint_c( diameter ) );
basefolder += "_MEM_";
basefolder += MEMVariant_to_string( memVariant );
basefolder += "_bvrf";
basefolder += std::to_string(int(bulkViscRateFactor));
if( useOmegaBulkAdaption ) basefolder += "Adapt";
basefolder += "_" + reconstructorType;
if( longDom ) basefolder += "_longDom";
if( broadDom ) basefolder += "_broadDom";
if( conserveMomentum ) basefolder += "_conserveMomentum";
if( !folderEnding.empty() ) basefolder += "_" + folderEnding;
WALBERLA_LOG_INFO_ON_ROOT("Basefolder for simulation results: " << basefolder);
// create base directory if it does not yet exist
filesystem::path tpath( basefolder );
if( !filesystem::exists( tpath ) )
filesystem::create_directory( tpath );
// 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
blockforest::communication::UniformBufferedScheme< Stencil_T > fullPDFCommunicationScheme( blocks );
fullPDFCommunicationScheme.addPackInfo( make_shared< field::communication::PackInfo< PdfField_T > >( pdfFieldID ) ); // full sync
// omega bulk adaption
using OmegaBulkAdapter_T = lbm_mesapd_coupling::OmegaBulkAdapter<ParticleAccessor_T, decltype(sphereSelector)>;
real_t defaultOmegaBulk = lbm_mesapd_coupling::omegaBulkFromOmega(omega, real_t(1));
real_t adaptionLayerSize = real_t(2);
shared_ptr<OmegaBulkAdapter_T> omegaBulkAdapter = make_shared<OmegaBulkAdapter_T>(blocks, omegaBulkFieldID, accessor, defaultOmegaBulk, omegaBulk, adaptionLayerSize, sphereSelector);
for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt) {
(*omegaBulkAdapter)(blockIt.get());
}
//////////////////////////////////////
// TIME LOOP FOR INITIAL SIMULATION //
//////////////////////////////////////
// Initialization simulation: fixed sphere and simulate until convergence (of drag force)
// create the timeloop
SweepTimeloop timeloopInit( blocks->getBlockStorage(), timestepsInit );
// add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
auto bhSweep = BoundaryHandling_T::getBlockSweep( boundaryHandlingID );
timeloopInit.add() << BeforeFunction( optimizedPDFCommunicationScheme, "LBM Communication" )
<< Sweep(bhSweep, "Boundary Handling" );
// stream + collide LBM step
#ifdef WALBERLA_BUILD_WITH_CODEGEN
auto lbmSweep = LatticeModel_T::Sweep( pdfFieldID );
timeloopInit.add() << Sweep( lbmSweep, "LB sweep" );
#else
auto lbmSweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag );
timeloopInit.add() << Sweep( makeSharedSweep( lbmSweep ), "cell-wise LB sweep" );
#endif
if( vtkIOInit )
{
auto pdfFieldVTKInit = vtk::createVTKOutput_BlockData( blocks, "fluid_field_init", writeFrequencyInit, 0, false, basefolder );
field::FlagFieldCellFilter< FlagField_T > fluidFilterInit( flagFieldID );
fluidFilterInit.addFlag( Fluid_Flag );
pdfFieldVTKInit->addCellInclusionFilter( fluidFilterInit );
pdfFieldVTKInit->addCellDataWriter( make_shared< lbm::VelocityVTKWriter< LatticeModel_T, float > >( pdfFieldID, "VelocityFromPDF" ) );
pdfFieldVTKInit->addCellDataWriter( make_shared< lbm::DensityVTKWriter < LatticeModel_T, float > >( pdfFieldID, "DensityFromPDF" ) );
timeloopInit.addFuncAfterTimeStep( vtk::writeFiles( pdfFieldVTKInit ), "VTK (fluid field data)" );
}
timeloopInit.addFuncAfterTimeStep( RemainingTimeLogger( timeloopInit.getNrOfTimeSteps(), real_t(30) ), "Remaining Time Logger" );
////////////////////////////////
// EXECUTE INITIAL SIMULATION //
////////////////////////////////
real_t gravity = real_t(1);
real_t GalileoSim = real_t(1);
real_t ReynoldsSim = real_t(1);
real_t u_ref = real_t(1);
WALBERLA_LOG_INFO_ON_ROOT("Starting initialization phase (sphere is kept fixed).");
WALBERLA_LOG_INFO_ON_ROOT("Iterating, and adapting the viscosity, until the targeted Galileo number is set.");
RestingSphereForceEvaluator<ParticleAccessor_T> forceEval( accessor, sphereUid, averageFrequency, basefolder );
while (true) {
WcTimingPool timeloopInitTiming;
WALBERLA_LOG_INFO_ON_ROOT("(Re-)Starting initial simulation.");
for( uint_t i = 1; i < timestepsInit; ++i ){
ps->forEachParticle(false, mesa_pd::kernel::SelectAll(), *accessor, resetHydrodynamicForceTorque, *accessor );
timeloopInit.singleStep( timeloopInitTiming );
reduceProperty.operator()<mesa_pd::HydrodynamicForceTorqueNotification>(*ps);
forceEval(timeloopInit.getCurrentTimeStep()+1);
// check if the relative change in the average drag force is below the specified convergence criterion
if (forceEval.getForceDiff() < convergenceLimit && i > 2 * std::max(averageFrequency, zlength) )
{
// conditions to break:
// - force diff sufficiently small
// - AND more time steps than twice the domain size in z direction to ensure new information due
// to possible viscosity change has traveled twice through domain,
// or twice the average frequency to have converged averages
WALBERLA_LOG_INFO_ON_ROOT("Drag force converged at time " << timeloopInit.getCurrentTimeStep());
break;
}
// if simulation gets unstable
if( std::isnan(forceEval.getDragForce()) )
{
WALBERLA_ABORT("NAN value detected in drag force during initial simulation, exiting....");
}
}
WALBERLA_LOG_INFO_ON_ROOT("Initial simulation has ended.")
//evaluate the gravitational force necessary to keep the sphere at a approximately fixed position
gravity = forceEval.getDragForce() / ( (densityRatio - real_t(1) ) * diameter * diameter * diameter * math::pi / real_t(6) );
GalileoSim = std::sqrt( ( densityRatio - real_t(1) ) * gravity * diameter * diameter * diameter ) / viscosity;
ReynoldsSim = uIn * diameter / viscosity;
u_ref = std::sqrt( std::fabs(densityRatio - real_t(1)) * gravity * diameter );
WALBERLA_LOG_INFO_ON_ROOT("Acting gravity = " << gravity );
WALBERLA_LOG_INFO_ON_ROOT("Simulated Galileo number = " << GalileoSim );
WALBERLA_LOG_INFO_ON_ROOT("Targeted Galileo number = " << Galileo );
WALBERLA_LOG_INFO_ON_ROOT("Reynolds number infty = " << ReynoldsSim );
if ( noViscosityIterations )
{
timeloopInitTiming.logResultOnRoot();
WALBERLA_LOG_INFO_ON_ROOT("Terminate iterations since viscosity should not be changed, flag \"--noViscosityIterations\"");
break;
}
// if simulated Galileo number is close enough to targeted Galileo number, stop the initial simulation
if( std::abs( GalileoSim - Galileo ) / Galileo < convergenceLimitGalileo )
{
timeloopInitTiming.logResultOnRoot();
WALBERLA_LOG_INFO_ON_ROOT("Iterations converged, simulated Galileo number is close enough to targeted one");
break;
}
// else update the simulation parameter accordingly and resume simulation
real_t diff = GalileoSim/Galileo;
viscosity = diff * viscosity;
omega = lbm::collision_model::omegaFromViscosity( viscosity );
lambda_e = lbm::collision_model::TRT::lambda_e( omega );
lambda_d = lbm::collision_model::TRT::lambda_d( omega, magicNumberTRT );
// also adapt omega bulk
omegaBulk = lbm_mesapd_coupling::omegaBulkFromOmega(omega, bulkViscRateFactor);
defaultOmegaBulk = lbm_mesapd_coupling::omegaBulkFromOmega(omega, real_t(1));
omegaBulkAdapter->setDefaultOmegaBulk(defaultOmegaBulk);
omegaBulkAdapter->setAdaptedOmegaBulk(omegaBulk);
// iterate all blocks with an iterator 'block' and change the collision model
for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
{
// get the field data out of the block
auto pdf = blockIt->getData< PdfField_T > ( pdfFieldID );
#ifdef WALBERLA_BUILD_WITH_CODEGEN
pdf->latticeModel().omega_magic_ = lambda_d;
pdf->latticeModel().omega_visc_ = lambda_e;
#else
pdf->latticeModel().collisionModel().resetWithMagicNumber( omega, magicNumberTRT );
#endif
(*omegaBulkAdapter)(blockIt.get());
}
WALBERLA_LOG_INFO_ON_ROOT("==> Adapting viscosity:");
WALBERLA_LOG_INFO_ON_ROOT("New viscosity = " << viscosity );
WALBERLA_LOG_INFO_ON_ROOT("New omega = " << omega );
WALBERLA_LOG_INFO_ON_ROOT("New omega bulk = " << omegaBulk );
}
if( averageForceTorqueOverTwoTimeStepsMode == 1 )
{
// maintain a good initial guess for averaging
ps->forEachParticle(false, mesa_pd::kernel::SelectLocal(), *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
}
ps->forEachParticle(false, mesa_pd::kernel::SelectAll(), *accessor, resetHydrodynamicForceTorque, *accessor );
///////////////
// TIME LOOP //
///////////////
// actual simulation: freely moving sphere with acting gravity
// calculate the number of timesteps
real_t dtSim = (averageForceTorqueOverTwoTimeStepsMode == 2) ? real_t(2) : real_t(1);
const real_t t_ref = ( diameter / u_ref );
const uint_t timesteps = uint_c( timestepsNonDim * t_ref / dtSim );
// set vtk write frequency accordingly
const uint_t writeFrequency = uint_c( vtkWriteSpacingNonDim * t_ref ); // vtk write frequency
const uint_t initWriteCallsToSkip = uint_c(std::max( 0, int( timesteps - ( vtkNumOutputFiles - uint_t(1) ) * writeFrequency - uint_t(1) ) ) ); // write calls to be skipped
WALBERLA_LOG_INFO_ON_ROOT("Starting simulation timeloop with " << timesteps << " timesteps!");
WALBERLA_LOG_INFO_ON_ROOT("Sphere is allowed to move freely under action of gravity");
SweepTimeloop timeloop( blocks->getBlockStorage(), uint_c(dtSim) * timesteps );
timeloop.addFuncBeforeTimeStep( RemainingTimeLogger( timeloop.getNrOfTimeSteps(), real_t(30) ), "Remaining Time Logger" );
// vtk output
auto pdfFieldVTK = vtk::createVTKOutput_BlockData( blocks, "fluid_field", writeFrequency, (vtkOutputGhostlayers) ? FieldGhostLayers : 0, false, basefolder );
pdfFieldVTK->setInitialWriteCallsToSkip( initWriteCallsToSkip );
pdfFieldVTK->addBeforeFunction( fullPDFCommunicationScheme );
// function to output plane infos for vtk output
pdfFieldVTK->addBeforeFunction(VTKInfoLogger<ParticleAccessor_T>( &timeloop, accessor, sphereUid, basefolder, uInfty ));
// create folder for log_vtk files to not pollute the basefolder
filesystem::path tpath2( basefolder+"/log_vtk" );
if( !filesystem::exists( tpath2 ) )
filesystem::create_directory( tpath2 );
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)" );
// 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(bhSweep, "Boundary Handling" );
// stream + collide LBM step
#ifdef WALBERLA_BUILD_WITH_CODEGEN
timeloop.add() << Sweep( lbmSweep, "LB sweep" );
#else
timeloop.add() << Sweep( makeSharedSweep( lbmSweep ), "cell-wise LB sweep" );
#endif
SweepTimeloop timeloopAfterParticles( blocks->getBlockStorage(), timesteps );
// sweep for updating the pe body mapping into the LBM simulation
if( memVariant == MEMVariant::CLI )
timeloopAfterParticles.add() << Sweep( lbm_mesapd_coupling::makeMovingParticleMapping<PdfField_T, BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, MEM_CLI_Flag, FormerMEM_Flag, sphereSelector, conserveMomentum), "Particle Mapping" );
else
timeloopAfterParticles.add() << Sweep( lbm_mesapd_coupling::makeMovingParticleMapping<PdfField_T, BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, MEM_BB_Flag, FormerMEM_Flag, sphereSelector, conserveMomentum), "Particle Mapping" );
// sweep for restoring PDFs in cells previously occupied by particles
if( reconstructorType == "EAN")
{
auto sphereNormalExtrapolationDirectionFinder = make_shared<lbm_mesapd_coupling::SphereNormalExtrapolationDirectionFinder>(blocks);
auto equilibriumAndNonEquilibriumSphereNormalReconstructor = lbm_mesapd_coupling::makeEquilibriumAndNonEquilibriumReconstructor<BoundaryHandling_T>(blocks, boundaryHandlingID, sphereNormalExtrapolationDirectionFinder, uint_t(3), true);
auto reconstructionManager = lbm_mesapd_coupling::makePdfReconstructionManager<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, FormerMEM_Flag, Fluid_Flag, equilibriumAndNonEquilibriumSphereNormalReconstructor, conserveMomentum);
timeloopAfterParticles.add() << BeforeFunction( fullPDFCommunicationScheme, "PDF Communication" )
<< Sweep( makeSharedSweep(reconstructionManager), "PDF Restore" );
} else if( reconstructorType == "Ext" )
{
auto sphereNormalExtrapolationDirectionFinder = make_shared<lbm_mesapd_coupling::SphereNormalExtrapolationDirectionFinder>(blocks);
auto extrapolationSphereNormalReconstructor = lbm_mesapd_coupling::makeExtrapolationReconstructor<BoundaryHandling_T,lbm_mesapd_coupling::SphereNormalExtrapolationDirectionFinder,true>(blocks, boundaryHandlingID, sphereNormalExtrapolationDirectionFinder, uint_t(3), true);
timeloopAfterParticles.add() << BeforeFunction( fullPDFCommunicationScheme, "PDF Communication" )
<< Sweep( makeSharedSweep(lbm_mesapd_coupling::makePdfReconstructionManager<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, FormerMEM_Flag, Fluid_Flag, extrapolationSphereNormalReconstructor, conserveMomentum) ), "PDF Restore" );
} else if( reconstructorType == "Grad")
{
auto gradReconstructor = lbm_mesapd_coupling::makeGradsMomentApproximationReconstructor<BoundaryHandling_T>(blocks, boundaryHandlingID, omega, false, true, true);
timeloopAfterParticles.add() << BeforeFunction( fullPDFCommunicationScheme, "PDF Communication" )
<< Sweep( makeSharedSweep(lbm_mesapd_coupling::makePdfReconstructionManager<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, FormerMEM_Flag, Fluid_Flag, gradReconstructor, conserveMomentum) ), "PDF Restore" );
} else if( reconstructorType == "Eq")
{
timeloopAfterParticles.add() << Sweep( makeSharedSweep(lbm_mesapd_coupling::makePdfReconstructionManager<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, FormerMEM_Flag, Fluid_Flag, conserveMomentum) ), "PDF Restore" );
} else {
WALBERLA_ABORT("Unknown reconstructor type " << reconstructorType);
}
// update bulk omega in all cells to adapt to changed particle position
if( useOmegaBulkAdaption )
{
timeloopAfterParticles.add() << Sweep( makeSharedSweep(omegaBulkAdapter), "Omega Bulk Adapter");
}
Vector3<real_t> extForcesOnSphere( real_t(0), real_t(0), - gravity * ( densityRatio - real_t(1) ) * diameter * diameter * diameter * math::pi / real_t(6));
lbm_mesapd_coupling::AddForceOnParticlesKernel addGravitationalForce(extForcesOnSphere);
mesa_pd::kernel::VelocityVerletPreForceUpdate vvIntegratorPreForce(dtSim);
mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(dtSim);
mesa_pd::kernel::ExplicitEuler explicitEulerIntegrator(dtSim);
MovingSpherePropertyEvaluator<ParticleAccessor_T> movingSpherePropertyEvaluator( accessor, sphereUid, basefolder, uInfty, Galileo, GalileoSim, gravity, viscosity, diameter, densityRatio );
////////////////////////
// EXECUTE SIMULATION //
////////////////////////
WcTimingPool timeloopTiming;
const bool useOpenMP = false;
// time loop
for (uint_t i = 0; i < timesteps; ++i )
{
// perform a single simulation step -> this contains LBM and setting of the hydrodynamic interactions
timeloop.singleStep( timeloopTiming );
if( averageForceTorqueOverTwoTimeStepsMode == 2)
{
// store current force and torque in fOld and tOld
ps->forEachParticle(useOpenMP, sphereSelector, *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
// reset f and t
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, resetHydrodynamicForceTorque, *accessor );
timeloop.singleStep( timeloopTiming );
// f = (f+fOld)/2
ps->forEachParticle(useOpenMP, sphereSelector, *accessor, averageHydrodynamicForceTorque, *accessor );
reduceProperty.operator()<mesa_pd::HydrodynamicForceTorqueNotification>(*ps);
}
else
{
reduceProperty.operator()<mesa_pd::HydrodynamicForceTorqueNotification>(*ps);
if( averageForceTorqueOverTwoTimeStepsMode == 1 )
{
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, averageHydrodynamicForceTorque, *accessor );
}
}
timeloopTiming["RPD"].start();
if( useVelocityVerlet)
{
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, vvIntegratorPreForce, *accessor);
syncCall();
}
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addHydrodynamicInteraction, *accessor );
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addGravitationalForce, *accessor );
reduceProperty.operator()<mesa_pd::ForceTorqueNotification>(*ps);
if( useVelocityVerlet ) {
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, vvIntegratorPostForce, *accessor);
}
else
{
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, explicitEulerIntegrator, *accessor);
}
syncCall();
timeloopTiming["RPD"].end();
// evaluation
timeloopTiming["Logging"].start();
movingSpherePropertyEvaluator((i+1)*uint_c(dtSim));
timeloopTiming["Logging"].end();
// reset after logging here
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, resetHydrodynamicForceTorque, *accessor );
timeloopAfterParticles.singleStep( timeloopTiming );
if(movingSpherePropertyEvaluator.getParticleHeight() > real_c(zlength) - real_t(3) * diameter)
{
pdfFieldVTK->forceWrite(timeloop.getCurrentTimeStep());
WALBERLA_LOG_WARNING_ON_ROOT("Sphere is too close to outflow, stopping simulation");
break;
}
}
timeloopTiming.logResultOnRoot();
return EXIT_SUCCESS;
}
} // namespace motion_settling_sphere
int main( int argc, char **argv ){
motion_settling_sphere::main(argc, argv);
}
apps/benchmarks/FluidParticleCoupling/ObliqueDryCollision.cpp
View file @ 7a6d3b57
...@@ -26,8 +26,8 @@ ...@@ -26,8 +26,8 @@
#include "mesa_pd/data/ShapeStorage.h" #include "mesa_pd/data/ShapeStorage.h"
#include "mesa_pd/kernel/DoubleCast.h" #include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/ExplicitEulerWithShape.h" #include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/VelocityVerletWithShape.h" #include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/kernel/LinearSpringDashpot.h" #include "mesa_pd/kernel/LinearSpringDashpot.h"
#include "mesa_pd/mpi/ReduceContactHistory.h" #include "mesa_pd/mpi/ReduceContactHistory.h"
...@@ -119,22 +119,24 @@ int main( int argc, char ** argv ) ...@@ -119,22 +119,24 @@ int main( int argc, char ** argv )
data::Particle&& p = *ps->create(); data::Particle&& p = *ps->create();
p.setPosition(Vec3(0,0,2*radius)); p.setPosition(Vec3(0,0,2*radius));
p.setLinearVelocity(Vec3(uTin, 0., -uNin)); p.setLinearVelocity(Vec3(uTin, 0., -uNin));
p.setInteractionRadius(radius);
p.setType(0); p.setType(0);
// create plane // create plane
data::Particle&& p0 = *ps->create(true); data::Particle&& p0 = *ps->create(true);
p0.setPosition(Vec3(0,0,0)); p0.setPosition(Vec3(0,0,0));
p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p0.setShapeID(ss->create<data::HalfSpace>(Vector3<real_t>(0,0,1))); p0.setShapeID(ss->create<data::HalfSpace>(Vector3<real_t>(0,0,1)));
p0.setType(0); p0.setType(0);
data::particle_flags::set(p0.getFlagsRef(), data::particle_flags::INFINITE); data::particle_flags::set(p0.getFlagsRef(), data::particle_flags::INFINITE);
data::particle_flags::set(p0.getFlagsRef(), data::particle_flags::FIXED); data::particle_flags::set(p0.getFlagsRef(), data::particle_flags::FIXED);
// velocity verlet // velocity verlet
kernel::VelocityVerletWithShapePreForceUpdate vvPreForce( dt ); kernel::VelocityVerletPreForceUpdate vvPreForce( dt );
kernel::VelocityVerletWithShapePostForceUpdate vvPostForce( dt ); kernel::VelocityVerletPostForceUpdate vvPostForce( dt );
// explicit euler // explicit euler
kernel::ExplicitEulerWithShape explEuler( dt ); kernel::ExplicitEuler explEuler( dt );
// collision response // collision response
collision_detection::AnalyticContactDetection acd; collision_detection::AnalyticContactDetection acd;
......
apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
View file @ 7a6d3b57
...@@ -80,10 +80,10 @@ ...@@ -80,10 +80,10 @@
#include "mesa_pd/domain/BlockForestDomain.h" #include "mesa_pd/domain/BlockForestDomain.h"
#include "mesa_pd/domain/BlockForestDataHandling.h" #include "mesa_pd/domain/BlockForestDataHandling.h"
#include "mesa_pd/kernel/DoubleCast.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/LinearSpringDashpot.h"
#include "mesa_pd/kernel/ParticleSelector.h" #include "mesa_pd/kernel/ParticleSelector.h"
#include "mesa_pd/kernel/VelocityVerletWithShape.h" #include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h" #include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h" #include "mesa_pd/mpi/ReduceProperty.h"
#include "mesa_pd/mpi/ReduceContactHistory.h" #include "mesa_pd/mpi/ReduceContactHistory.h"
...@@ -268,15 +268,9 @@ public: ...@@ -268,15 +268,9 @@ public:
WALBERLA_MPI_SECTION() WALBERLA_MPI_SECTION()
{ {
mpi::allReduceInplace( pos[0], mpi::SUM ); mpi::allReduceInplace( pos, mpi::SUM );
mpi::allReduceInplace( pos[1], mpi::SUM ); mpi::allReduceInplace( transVel, mpi::SUM );
mpi::allReduceInplace( pos[2], mpi::SUM ); mpi::allReduceInplace( angVel, mpi::SUM );
mpi::allReduceInplace( transVel[0], mpi::SUM );
mpi::allReduceInplace( transVel[1], mpi::SUM );
mpi::allReduceInplace( transVel[2], mpi::SUM );
mpi::allReduceInplace( angVel[0], mpi::SUM );
mpi::allReduceInplace( angVel[1], mpi::SUM );
mpi::allReduceInplace( angVel[2], mpi::SUM );
} }
position_ = pos; position_ = pos;
...@@ -370,6 +364,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con ...@@ -370,6 +364,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
// create bounding planes // create bounding planes
mesa_pd::data::Particle&& p0 = *ps->create(true); mesa_pd::data::Particle&& p0 = *ps->create(true);
p0.setPosition(simulationDomain.minCorner()); 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.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,1) ));
p0.setOwner(mpi::MPIManager::instance()->rank()); p0.setOwner(mpi::MPIManager::instance()->rank());
p0.setType(0); p0.setType(0);
...@@ -381,6 +376,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con ...@@ -381,6 +376,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
//only create top plane when no outflow BC should be set there //only create top plane when no outflow BC should be set there
mesa_pd::data::Particle&& p1 = *ps->create(true); mesa_pd::data::Particle&& p1 = *ps->create(true);
p1.setPosition(simulationDomain.maxCorner()); 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.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,-1) ));
p1.setOwner(mpi::MPIManager::instance()->rank()); p1.setOwner(mpi::MPIManager::instance()->rank());
p1.setType(0); p1.setType(0);
...@@ -449,9 +445,7 @@ Vector3<real_t> getForce(walberla::id_t uid, ParticleAccessor_T & ac) ...@@ -449,9 +445,7 @@ Vector3<real_t> getForce(walberla::id_t uid, ParticleAccessor_T & ac)
} }
WALBERLA_MPI_SECTION() WALBERLA_MPI_SECTION()
{ {
mpi::allReduceInplace( force[0], mpi::SUM ); mpi::allReduceInplace( force, mpi::SUM );
mpi::allReduceInplace( force[1], mpi::SUM );
mpi::allReduceInplace( force[2], mpi::SUM );
} }
return force; return force;
} }
...@@ -467,9 +461,7 @@ Vector3<real_t> getTorque(walberla::id_t uid, ParticleAccessor_T & ac) ...@@ -467,9 +461,7 @@ Vector3<real_t> getTorque(walberla::id_t uid, ParticleAccessor_T & ac)
} }
WALBERLA_MPI_SECTION() WALBERLA_MPI_SECTION()
{ {
mpi::allReduceInplace( torque[0], mpi::SUM ); mpi::allReduceInplace( torque, mpi::SUM );
mpi::allReduceInplace( torque[1], mpi::SUM );
mpi::allReduceInplace( torque[2], mpi::SUM );
} }
return torque; return torque;
} }
...@@ -815,7 +807,7 @@ int main( int argc, char **argv ) ...@@ -815,7 +807,7 @@ int main( int argc, char **argv )
//////////////////////// ////////////////////////
// add omega bulk field // 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 // create the lattice model
real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega ); real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega );
...@@ -853,9 +845,9 @@ int main( int argc, char **argv ) ...@@ -853,9 +845,9 @@ int main( int argc, char **argv )
syncCall(); syncCall();
real_t timeStepSizeRPD = real_t(1)/real_t(numRPDSubCycles); real_t timeStepSizeRPD = real_t(1)/real_t(numRPDSubCycles);
mesa_pd::kernel::ExplicitEulerWithShape explEulerIntegrator(timeStepSizeRPD); mesa_pd::kernel::ExplicitEuler explEulerIntegrator(timeStepSizeRPD);
mesa_pd::kernel::VelocityVerletWithShapePreForceUpdate vvIntegratorPreForce(timeStepSizeRPD); mesa_pd::kernel::VelocityVerletPreForceUpdate vvIntegratorPreForce(timeStepSizeRPD);
mesa_pd::kernel::VelocityVerletWithShapePostForceUpdate vvIntegratorPostForce(timeStepSizeRPD); mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(timeStepSizeRPD);
// linear model // linear model
mesa_pd::kernel::LinearSpringDashpot linearCollisionResponse(1); mesa_pd::kernel::LinearSpringDashpot linearCollisionResponse(1);
......
apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp
View file @ 7a6d3b57
...@@ -72,9 +72,9 @@ ...@@ -72,9 +72,9 @@
#include "mesa_pd/data/shape/Sphere.h" #include "mesa_pd/data/shape/Sphere.h"
#include "mesa_pd/domain/BlockForestDomain.h" #include "mesa_pd/domain/BlockForestDomain.h"
#include "mesa_pd/kernel/DoubleCast.h" #include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/ExplicitEulerWithShape.h" #include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/ParticleSelector.h" #include "mesa_pd/kernel/ParticleSelector.h"
#include "mesa_pd/kernel/VelocityVerletWithShape.h" #include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h" #include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h" #include "mesa_pd/mpi/ReduceProperty.h"
#include "mesa_pd/mpi/ContactFilter.h" #include "mesa_pd/mpi/ContactFilter.h"
...@@ -186,9 +186,10 @@ class SpherePropertyLogger ...@@ -186,9 +186,10 @@ class SpherePropertyLogger
public: public:
SpherePropertyLogger( const shared_ptr< ParticleAccessor_T > & ac, walberla::id_t sphereUid, SpherePropertyLogger( const shared_ptr< ParticleAccessor_T > & ac, walberla::id_t sphereUid,
const std::string & fileName, bool fileIO, const std::string & fileName, bool fileIO,
real_t dx_SI, real_t dt_SI, real_t diameter, real_t gravitationalForceMag) : real_t dx_SI, real_t dt_SI, real_t diameter, real_t gravitationalForceMag, real_t uref) :
ac_( ac ), sphereUid_( sphereUid ), fileName_( fileName ), fileIO_(fileIO), ac_( ac ), sphereUid_( sphereUid ), fileName_( fileName ), fileIO_(fileIO),
dx_SI_( dx_SI ), dt_SI_( dt_SI ), diameter_( diameter ), gravitationalForceMag_( gravitationalForceMag ), dx_SI_( dx_SI ), dt_SI_( dt_SI ), diameter_( diameter ),
gravitationalForceMag_( gravitationalForceMag ), uref_(uref), tref_(diameter / uref),
position_( real_t(0) ), maxVelocity_( real_t(0) ) position_( real_t(0) ), maxVelocity_( real_t(0) )
{ {
if ( fileIO_ ) if ( fileIO_ )
...@@ -197,7 +198,7 @@ public: ...@@ -197,7 +198,7 @@ public:
{ {
std::ofstream file; std::ofstream file;
file.open( fileName_.c_str() ); file.open( fileName_.c_str() );
file << "#\t t\t posZ\t gapZ/D\t velZ\t velZ_SI\t fZ\t fZ/fGravi\n"; file << "#\t t\t t/tref\t posZ\t gapZ/D\t velZ\t velZ_SI\t velZ/uref\t fZ\t fZ/fGravi\n";
file.close(); file.close();
} }
} }
...@@ -222,17 +223,9 @@ public: ...@@ -222,17 +223,9 @@ public:
WALBERLA_MPI_SECTION() WALBERLA_MPI_SECTION()
{ {
mpi::allReduceInplace( pos[0], mpi::SUM ); mpi::allReduceInplace( pos, mpi::SUM );
mpi::allReduceInplace( pos[1], mpi::SUM ); mpi::allReduceInplace( transVel, mpi::SUM );
mpi::allReduceInplace( pos[2], mpi::SUM ); mpi::allReduceInplace( hydForce, mpi::SUM );
mpi::allReduceInplace( transVel[0], mpi::SUM );
mpi::allReduceInplace( transVel[1], mpi::SUM );
mpi::allReduceInplace( transVel[2], mpi::SUM );
mpi::allReduceInplace( hydForce[0], mpi::SUM );
mpi::allReduceInplace( hydForce[1], mpi::SUM );
mpi::allReduceInplace( hydForce[2], mpi::SUM );
} }
position_ = pos[2]; position_ = pos[2];
...@@ -264,10 +257,9 @@ private: ...@@ -264,10 +257,9 @@ private:
auto velocity_SI = velocity * dx_SI_ / dt_SI_; auto velocity_SI = velocity * dx_SI_ / dt_SI_;
auto normalizedHydForce = hydForce / gravitationalForceMag_; auto normalizedHydForce = hydForce / gravitationalForceMag_;
file << timestep << "\t" << real_c(timestep) * dt_SI_ << "\t" << real_c(timestep) / tref_
file << timestep << "\t" << real_c(timestep) * dt_SI_ << "\t"
<< "\t" << position[2] << "\t" << scaledPosition[2] - real_t(0.5) << "\t" << position[2] << "\t" << scaledPosition[2] - real_t(0.5)
<< "\t" << velocity[2] << "\t" << velocity_SI[2] << "\t" << velocity[2] << "\t" << velocity_SI[2] << "\t" << velocity[2] / uref_
<< "\t" << hydForce[2] << "\t" << normalizedHydForce[2] << "\t" << hydForce[2] << "\t" << normalizedHydForce[2]
<< "\n"; << "\n";
file.close(); file.close();
...@@ -278,7 +270,7 @@ private: ...@@ -278,7 +270,7 @@ private:
const walberla::id_t sphereUid_; const walberla::id_t sphereUid_;
std::string fileName_; std::string fileName_;
bool fileIO_; bool fileIO_;
real_t dx_SI_, dt_SI_, diameter_, gravitationalForceMag_; real_t dx_SI_, dt_SI_, diameter_, gravitationalForceMag_, uref_, tref_;
real_t position_; real_t position_;
real_t maxVelocity_; real_t maxVelocity_;
...@@ -290,6 +282,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con ...@@ -290,6 +282,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
// create bounding planes // create bounding planes
mesa_pd::data::Particle p0 = *ps->create(true); mesa_pd::data::Particle p0 = *ps->create(true);
p0.setPosition(simulationDomain.minCorner()); 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.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,1) ));
p0.setOwner(mpi::MPIManager::instance()->rank()); p0.setOwner(mpi::MPIManager::instance()->rank());
p0.setType(0); p0.setType(0);
...@@ -298,6 +291,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con ...@@ -298,6 +291,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p1 = *ps->create(true); mesa_pd::data::Particle p1 = *ps->create(true);
p1.setPosition(simulationDomain.maxCorner()); 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.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,-1) ));
p1.setOwner(mpi::MPIManager::instance()->rank()); p1.setOwner(mpi::MPIManager::instance()->rank());
p1.setType(0); p1.setType(0);
...@@ -306,6 +300,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con ...@@ -306,6 +300,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p2 = *ps->create(true); mesa_pd::data::Particle p2 = *ps->create(true);
p2.setPosition(simulationDomain.minCorner()); p2.setPosition(simulationDomain.minCorner());
p2.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p2.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(1,0,0) )); p2.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(1,0,0) ));
p2.setOwner(mpi::MPIManager::instance()->rank()); p2.setOwner(mpi::MPIManager::instance()->rank());
p2.setType(0); p2.setType(0);
...@@ -314,6 +309,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con ...@@ -314,6 +309,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p3 = *ps->create(true); mesa_pd::data::Particle p3 = *ps->create(true);
p3.setPosition(simulationDomain.maxCorner()); p3.setPosition(simulationDomain.maxCorner());
p3.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p3.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(-1,0,0) )); p3.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(-1,0,0) ));
p3.setOwner(mpi::MPIManager::instance()->rank()); p3.setOwner(mpi::MPIManager::instance()->rank());
p3.setType(0); p3.setType(0);
...@@ -322,6 +318,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con ...@@ -322,6 +318,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p4 = *ps->create(true); mesa_pd::data::Particle p4 = *ps->create(true);
p4.setPosition(simulationDomain.minCorner()); p4.setPosition(simulationDomain.minCorner());
p4.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p4.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,1,0) )); p4.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,1,0) ));
p4.setOwner(mpi::MPIManager::instance()->rank()); p4.setOwner(mpi::MPIManager::instance()->rank());
p4.setType(0); p4.setType(0);
...@@ -330,6 +327,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con ...@@ -330,6 +327,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p5 = *ps->create(true); mesa_pd::data::Particle p5 = *ps->create(true);
p5.setPosition(simulationDomain.maxCorner()); p5.setPosition(simulationDomain.maxCorner());
p5.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p5.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,-1,0) )); p5.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,-1,0) ));
p5.setOwner(mpi::MPIManager::instance()->rank()); p5.setOwner(mpi::MPIManager::instance()->rank());
p5.setType(0); p5.setType(0);
...@@ -385,6 +383,8 @@ int main( int argc, char **argv ) ...@@ -385,6 +383,8 @@ int main( int argc, char **argv )
real_t adaptionLayerSize = real_t(2); real_t adaptionLayerSize = real_t(2);
bool useLubricationCorrection = true; bool useLubricationCorrection = true;
bool useGalileoParameterization = false;
for( int i = 1; i < argc; ++i ) for( int i = 1; i < argc; ++i )
{ {
if( std::strcmp( argv[i], "--shortrun" ) == 0 ) { shortrun = true; continue; } if( std::strcmp( argv[i], "--shortrun" ) == 0 ) { shortrun = true; continue; }
...@@ -406,6 +406,7 @@ int main( int argc, char **argv ) ...@@ -406,6 +406,7 @@ int main( int argc, char **argv )
if( std::strcmp( argv[i], "--useOmegaBulkAdaption" ) == 0 ) { useOmegaBulkAdaption = true; continue; } if( std::strcmp( argv[i], "--useOmegaBulkAdaption" ) == 0 ) { useOmegaBulkAdaption = true; continue; }
if( std::strcmp( argv[i], "--adaptionLayerSize" ) == 0 ) { adaptionLayerSize = real_c(std::atof( argv[++i] )); continue; } if( std::strcmp( argv[i], "--adaptionLayerSize" ) == 0 ) { adaptionLayerSize = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--noLubricationCorrection" ) == 0 ) { useLubricationCorrection = false; continue; } if( std::strcmp( argv[i], "--noLubricationCorrection" ) == 0 ) { useLubricationCorrection = false; continue; }
if( std::strcmp( argv[i], "--useGalileoParameterization" ) == 0 ) { useGalileoParameterization = true; continue; }
WALBERLA_ABORT("Unrecognized command line argument found: " << argv[i]); WALBERLA_ABORT("Unrecognized command line argument found: " << argv[i]);
} }
...@@ -433,6 +434,9 @@ int main( int argc, char **argv ) ...@@ -433,6 +434,9 @@ int main( int argc, char **argv )
real_t densityFluid_SI; real_t densityFluid_SI;
real_t dynamicViscosityFluid_SI; real_t dynamicViscosityFluid_SI;
real_t expectedSettlingVelocity_SI; real_t expectedSettlingVelocity_SI;
// expected velocity given as uMax in experiments of ten Cate (Table 2, E1-E4), with uInfty from Table 1
// > slightly different Re than in ten Cate's Table 1
switch( fluidType ) switch( fluidType )
{ {
case 1: case 1:
...@@ -494,8 +498,9 @@ int main( int argc, char **argv ) ...@@ -494,8 +498,9 @@ int main( int argc, char **argv )
const real_t diameter = diameter_SI / dx_SI; const real_t diameter = diameter_SI / dx_SI;
const real_t sphereVolume = math::pi / real_t(6) * diameter * diameter * diameter; const real_t sphereVolume = math::pi / real_t(6) * diameter * diameter * diameter;
//const real_t dt_SI = characteristicVelocity / ug_SI * dx_SI; // this uses Ga for parameterization
const real_t dt_SI = characteristicVelocity / expectedSettlingVelocity_SI * dx_SI; // this uses Re for parameterization (only possible since settling velocity is known) const real_t dt_SI = (useGalileoParameterization) ? characteristicVelocity / ug_SI * dx_SI : // this uses Ga for parameterization, where ug is the characteristic velocity
characteristicVelocity / expectedSettlingVelocity_SI * dx_SI; // this uses Re for parameterization (only possible since settling velocity is known)
const real_t viscosity = kinematicViscosityFluid_SI * dt_SI / ( dx_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 omega = lbm::collision_model::omegaFromViscosity(viscosity);
...@@ -602,7 +607,7 @@ int main( int argc, char **argv ) ...@@ -602,7 +607,7 @@ int main( int argc, char **argv )
//////////////////////// ////////////////////////
// add omega bulk field // 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 // create the lattice model
real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega ); real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega );
...@@ -639,9 +644,9 @@ int main( int argc, char **argv ) ...@@ -639,9 +644,9 @@ int main( int argc, char **argv )
syncCall(); syncCall();
mesa_pd::kernel::ExplicitEulerWithShape explEulerIntegrator(real_t(1)/real_t(numRPDSubCycles)); mesa_pd::kernel::ExplicitEuler explEulerIntegrator(real_t(1)/real_t(numRPDSubCycles));
mesa_pd::kernel::VelocityVerletWithShapePreForceUpdate vvIntegratorPreForce(real_t(1)/real_t(numRPDSubCycles)); mesa_pd::kernel::VelocityVerletPreForceUpdate vvIntegratorPreForce(real_t(1)/real_t(numRPDSubCycles));
mesa_pd::kernel::VelocityVerletWithShapePostForceUpdate vvIntegratorPostForce(real_t(1)/real_t(numRPDSubCycles)); mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(real_t(1)/real_t(numRPDSubCycles));
mesa_pd::mpi::ReduceProperty reduceProperty; mesa_pd::mpi::ReduceProperty reduceProperty;
...@@ -779,17 +784,19 @@ int main( int argc, char **argv ) ...@@ -779,17 +784,19 @@ int main( int argc, char **argv )
// evaluation functionality // evaluation functionality
std::string loggingFileName( baseFolder + "/LoggingSettlingSphere_"); std::string loggingFileName( baseFolder + "/LoggingSettlingSphere_");
loggingFileName += std::to_string(fluidType); loggingFileName += std::to_string(fluidType);
loggingFileName += "_res" + std::to_string(numberOfCellsInHorizontalDirection);
loggingFileName += "_recon" + reconstructorType; loggingFileName += "_recon" + reconstructorType;
loggingFileName += "_bvrf" + std::to_string(uint_c(bulkViscRateFactor)); loggingFileName += "_bvrf" + std::to_string(uint_c(bulkViscRateFactor));
loggingFileName += "_mn" + std::to_string(float(magicNumber)); loggingFileName += "_mn" + std::to_string(float(magicNumber));
if( useOmegaBulkAdaption ) loggingFileName += "_uOBA" + std::to_string(uint_c(adaptionLayerSize)); if( useOmegaBulkAdaption ) loggingFileName += "_uOBA" + std::to_string(uint_c(adaptionLayerSize));
if( useGalileoParameterization ) loggingFileName += "_Ga";
if( !fileNameEnding.empty()) loggingFileName += "_" + fileNameEnding; if( !fileNameEnding.empty()) loggingFileName += "_" + fileNameEnding;
loggingFileName += ".txt"; loggingFileName += ".txt";
if( fileIO ) if( fileIO )
{ {
WALBERLA_LOG_INFO_ON_ROOT(" - writing logging output to file \"" << loggingFileName << "\""); WALBERLA_LOG_INFO_ON_ROOT(" - writing logging output to file \"" << loggingFileName << "\"");
} }
SpherePropertyLogger<ParticleAccessor_T> logger( accessor, sphereUid, loggingFileName, fileIO, dx_SI, dt_SI, diameter, -gravitationalForce[2] ); SpherePropertyLogger<ParticleAccessor_T> logger( accessor, sphereUid, loggingFileName, fileIO, dx_SI, dt_SI, diameter, -gravitationalForce[2], characteristicVelocity );
//////////////////////// ////////////////////////
...@@ -832,7 +839,7 @@ int main( int argc, char **argv ) ...@@ -832,7 +839,7 @@ int main( int argc, char **argv )
syncCall(); syncCall();
} }
ps->forEachParticle(useOpenMP, sphereSelector, *accessor, addHydrodynamicInteraction, *accessor ); ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addHydrodynamicInteraction, *accessor );
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addGravitationalForce, *accessor ); ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addGravitationalForce, *accessor );
// lubrication correction // lubrication correction
......
apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp
View file @ 7a6d3b57
...@@ -76,7 +76,7 @@ ...@@ -76,7 +76,7 @@
#include "mesa_pd/data/shape/Sphere.h" #include "mesa_pd/data/shape/Sphere.h"
#include "mesa_pd/domain/BlockForestDomain.h" #include "mesa_pd/domain/BlockForestDomain.h"
#include "mesa_pd/kernel/DoubleCast.h" #include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/ExplicitEulerWithShape.h" #include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/ParticleSelector.h" #include "mesa_pd/kernel/ParticleSelector.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h" #include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h" #include "mesa_pd/mpi/ReduceProperty.h"
...@@ -129,7 +129,8 @@ const uint_t FieldGhostLayers = 1; ...@@ -129,7 +129,8 @@ const uint_t FieldGhostLayers = 1;
const FlagUID Fluid_Flag( "fluid" ); const FlagUID Fluid_Flag( "fluid" );
const FlagUID NoSlip_Flag( "no slip" ); const FlagUID NoSlip_Flag( "no slip" );
const FlagUID MO_Flag( "moving obstacle" ); const FlagUID MO_SBB_Flag( "moving obstacle sbb" );
const FlagUID MO_CLI_Flag( "moving obstacle cli" );
const FlagUID FormerMO_Flag( "former moving obstacle" ); const FlagUID FormerMO_Flag( "former moving obstacle" );
///////////////////////////////////// /////////////////////////////////////
...@@ -141,8 +142,9 @@ class MyBoundaryHandling ...@@ -141,8 +142,9 @@ class MyBoundaryHandling
public: public:
using NoSlip_T = lbm::NoSlip< LatticeModel_T, flag_t >; using NoSlip_T = lbm::NoSlip< LatticeModel_T, flag_t >;
using MO_T = lbm_mesapd_coupling::CurvedLinear< LatticeModel_T, FlagField_T, ParticleAccessor_T >; using MO_SBB_T = lbm_mesapd_coupling::SimpleBB< LatticeModel_T, FlagField_T, ParticleAccessor_T >;
using Type = BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T, MO_T >; using MO_CLI_T = lbm_mesapd_coupling::CurvedLinear< LatticeModel_T, FlagField_T, ParticleAccessor_T >;
using Type = BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T, MO_SBB_T, MO_CLI_T >;
MyBoundaryHandling( const BlockDataID & flagFieldID, const BlockDataID & pdfFieldID, MyBoundaryHandling( const BlockDataID & flagFieldID, const BlockDataID & pdfFieldID,
const BlockDataID & particleFieldID, const shared_ptr<ParticleAccessor_T>& ac) : const BlockDataID & particleFieldID, const shared_ptr<ParticleAccessor_T>& ac) :
...@@ -161,7 +163,8 @@ public: ...@@ -161,7 +163,8 @@ public:
Type * handling = new Type( "moving obstacle boundary handling", flagField, fluid, Type * handling = new Type( "moving obstacle boundary handling", flagField, fluid,
NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ), NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
MO_T( "MO", MO_Flag, pdfField, flagField, particleField, ac_, fluid, *storage, *block ) ); MO_SBB_T( "SBB", MO_SBB_Flag, pdfField, flagField, particleField, ac_, fluid, *storage, *block ),
MO_CLI_T( "CLI", MO_CLI_Flag, pdfField, flagField, particleField, ac_, fluid, *storage, *block ) );
handling->fillWithDomain( FieldGhostLayers ); handling->fillWithDomain( FieldGhostLayers );
...@@ -192,7 +195,7 @@ public: ...@@ -192,7 +195,7 @@ public:
const std::string & fileName, bool fileIO, const std::string & fileName, bool fileIO,
real_t diameter, real_t velocity, real_t forceMag) : real_t diameter, real_t velocity, real_t forceMag) :
ac_( ac ), sphereUid_( sphereUid ), fileName_( fileName ), fileIO_(fileIO), ac_( ac ), sphereUid_( sphereUid ), fileName_( fileName ), fileIO_(fileIO),
diameter_( diameter ), velocity_(velocity), forceMag_( forceMag ), diameter_( diameter ), velocityRef_(velocity), forceMag_( forceMag ),
position_( real_t(0) ) position_( real_t(0) )
{ {
if ( fileIO_ ) if ( fileIO_ )
...@@ -201,7 +204,7 @@ public: ...@@ -201,7 +204,7 @@ public:
{ {
std::ofstream file; std::ofstream file;
file.open( fileName_.c_str() ); file.open( fileName_.c_str() );
file << "#\t t\t posZ\t posZ/D\t fZ\t fZ/fSt\n"; file << "#\t t\t posZ\t posZ/D\t velZ\t velZ/velRef\t fZ\t fZ/fSt\n";
file.close(); file.close();
} }
} }
...@@ -210,6 +213,7 @@ public: ...@@ -210,6 +213,7 @@ public:
void operator()() void operator()()
{ {
real_t pos(real_t(0)); real_t pos(real_t(0));
real_t vel(real_t(0));
real_t hydForce(real_t(0)); real_t hydForce(real_t(0));
size_t idx = ac_->uidToIdx(sphereUid_); size_t idx = ac_->uidToIdx(sphereUid_);
...@@ -218,6 +222,7 @@ public: ...@@ -218,6 +222,7 @@ public:
if(!mesa_pd::data::particle_flags::isSet( ac_->getFlags(idx), mesa_pd::data::particle_flags::GHOST)) if(!mesa_pd::data::particle_flags::isSet( ac_->getFlags(idx), mesa_pd::data::particle_flags::GHOST))
{ {
pos = ac_->getPosition(idx)[2]; pos = ac_->getPosition(idx)[2];
vel = ac_->getLinearVelocity(idx)[2];
hydForce = ac_->getHydrodynamicForce(idx)[2]; hydForce = ac_->getHydrodynamicForce(idx)[2];
} }
} }
...@@ -225,10 +230,12 @@ public: ...@@ -225,10 +230,12 @@ public:
WALBERLA_MPI_SECTION() WALBERLA_MPI_SECTION()
{ {
mpi::allReduceInplace( pos, mpi::SUM ); mpi::allReduceInplace( pos, mpi::SUM );
mpi::allReduceInplace( vel, mpi::SUM );
mpi::allReduceInplace( hydForce, mpi::SUM ); mpi::allReduceInplace( hydForce, mpi::SUM );
} }
position_ = pos; position_ = pos;
velocity_ = vel;
hydForce_ = hydForce; hydForce_ = hydForce;
} }
...@@ -238,7 +245,12 @@ public: ...@@ -238,7 +245,12 @@ public:
return position_; return position_;
} }
void writeToFile( const uint_t timestep, real_t density1, real_t density2, real_t totalMass ) real_t getForce() const
{
return hydForce_;
}
void writeToFile( const uint_t timestep, real_t density1, real_t density2, real_t totalMass, uint_t countSolidCells, uint_t countFluidSolidLinks, Vector3<real_t> fluidVelocity )
{ {
WALBERLA_ROOT_SECTION() WALBERLA_ROOT_SECTION()
{ {
...@@ -249,10 +261,13 @@ public: ...@@ -249,10 +261,13 @@ public:
auto normalizedHydForce = hydForce_ / std::abs(forceMag_); auto normalizedHydForce = hydForce_ / std::abs(forceMag_);
file << std::setprecision(10) file << std::setprecision(10)
<< timestep << "\t" << real_c(timestep) / (diameter_ / velocity_) << "\t" << timestep << "\t" << real_c(timestep) / (diameter_ / velocityRef_) << "\t"
<< "\t" << position_ << "\t" << scaledPosition << "\t" << position_ << "\t" << scaledPosition
<< "\t" << velocity_ << "\t" << velocity_ / velocityRef_
<< "\t" << hydForce_ << "\t" << normalizedHydForce << "\t" << hydForce_ << "\t" << normalizedHydForce
<< "\t" << density1 << "\t" << density2 << "\t" << totalMass << "\t" << density1 << "\t" << density2 << "\t" << totalMass
<< "\t" << countSolidCells << "\t" << countFluidSolidLinks
<< "\t" << fluidVelocity[0] << "\t" << fluidVelocity[1] << "\t" << fluidVelocity[2]
<< "\n"; << "\n";
file.close(); file.close();
} }
...@@ -263,43 +278,51 @@ private: ...@@ -263,43 +278,51 @@ private:
const walberla::id_t sphereUid_; const walberla::id_t sphereUid_;
std::string fileName_; std::string fileName_;
bool fileIO_; bool fileIO_;
real_t diameter_, velocity_, forceMag_; real_t diameter_, velocityRef_, forceMag_;
real_t position_, hydForce_; real_t position_, velocity_, hydForce_;
}; };
void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, const shared_ptr<mesa_pd::data::ShapeStorage> & ss, const math::AABB & simulationDomain) void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, const shared_ptr<mesa_pd::data::ShapeStorage> & ss, const math::AABB & simulationDomain, real_t velocity = real_t(0))
{ {
mesa_pd::data::Particle p2 = *ps->create(true); mesa_pd::data::Particle p2 = *ps->create(true);
p2.setPosition(simulationDomain.minCorner()); p2.setPosition(simulationDomain.minCorner());
p2.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p2.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(1,0,0) )); p2.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(1,0,0) ));
p2.setOwner(mpi::MPIManager::instance()->rank()); p2.setOwner(mpi::MPIManager::instance()->rank());
p2.setType(0); p2.setType(0);
p2.setLinearVelocity(Vector3<real_t>(real_t(0),real_t(0),velocity));
mesa_pd::data::particle_flags::set(p2.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE); mesa_pd::data::particle_flags::set(p2.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
mesa_pd::data::particle_flags::set(p2.getFlagsRef(), mesa_pd::data::particle_flags::FIXED); mesa_pd::data::particle_flags::set(p2.getFlagsRef(), mesa_pd::data::particle_flags::FIXED);
mesa_pd::data::Particle p3 = *ps->create(true); mesa_pd::data::Particle p3 = *ps->create(true);
p3.setPosition(simulationDomain.maxCorner()); p3.setPosition(simulationDomain.maxCorner());
p3.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p3.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(-1,0,0) )); p3.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(-1,0,0) ));
p3.setOwner(mpi::MPIManager::instance()->rank()); p3.setOwner(mpi::MPIManager::instance()->rank());
p3.setType(0); p3.setType(0);
p3.setLinearVelocity(Vector3<real_t>(real_t(0),real_t(0),velocity));
mesa_pd::data::particle_flags::set(p3.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE); mesa_pd::data::particle_flags::set(p3.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
mesa_pd::data::particle_flags::set(p3.getFlagsRef(), mesa_pd::data::particle_flags::FIXED); mesa_pd::data::particle_flags::set(p3.getFlagsRef(), mesa_pd::data::particle_flags::FIXED);
mesa_pd::data::Particle p4 = *ps->create(true); mesa_pd::data::Particle p4 = *ps->create(true);
p4.setPosition(simulationDomain.minCorner()); p4.setPosition(simulationDomain.minCorner());
p4.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p4.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,1,0) )); p4.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,1,0) ));
p4.setOwner(mpi::MPIManager::instance()->rank()); p4.setOwner(mpi::MPIManager::instance()->rank());
p4.setType(0); p4.setType(0);
p4.setLinearVelocity(Vector3<real_t>(real_t(0),real_t(0),velocity));
mesa_pd::data::particle_flags::set(p4.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE); mesa_pd::data::particle_flags::set(p4.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
mesa_pd::data::particle_flags::set(p4.getFlagsRef(), mesa_pd::data::particle_flags::FIXED); mesa_pd::data::particle_flags::set(p4.getFlagsRef(), mesa_pd::data::particle_flags::FIXED);
mesa_pd::data::Particle p5 = *ps->create(true); mesa_pd::data::Particle p5 = *ps->create(true);
p5.setPosition(simulationDomain.maxCorner()); p5.setPosition(simulationDomain.maxCorner());
p5.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p5.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,-1,0) )); p5.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,-1,0) ));
p5.setOwner(mpi::MPIManager::instance()->rank()); p5.setOwner(mpi::MPIManager::instance()->rank());
p5.setType(0); p5.setType(0);
p5.setLinearVelocity(Vector3<real_t>(real_t(0),real_t(0),velocity));
mesa_pd::data::particle_flags::set(p5.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE); mesa_pd::data::particle_flags::set(p5.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
mesa_pd::data::particle_flags::set(p5.getFlagsRef(), mesa_pd::data::particle_flags::FIXED); mesa_pd::data::particle_flags::set(p5.getFlagsRef(), mesa_pd::data::particle_flags::FIXED);
...@@ -324,6 +347,25 @@ real_t getDensityAtPosition(const shared_ptr<StructuredBlockStorage> & blocks, B ...@@ -324,6 +347,25 @@ real_t getDensityAtPosition(const shared_ptr<StructuredBlockStorage> & blocks, B
return density; return density;
} }
template< typename VelocityInterpolator_T>
Vector3<real_t> getVelocityAtPosition(const shared_ptr<StructuredBlockStorage> & blocks, BlockDataID velocityInterpolatorID, Vector3<real_t> position)
{
Vector3<real_t> vel(real_t(0));
for( auto & block: *blocks)
{
if(block.getAABB().contains(position))
{
auto velocityInterpolator = block.getData<VelocityInterpolator_T>(velocityInterpolatorID);
velocityInterpolator->get(position, &vel);
}
}
mpi::reduceInplace(vel, mpi::SUM);
return vel;
}
template< typename BoundaryHandling_T> template< typename BoundaryHandling_T>
real_t getAverageDensityInSystem(const shared_ptr<StructuredBlockStorage> & blocks, BlockDataID pdfFieldID, BlockDataID boundaryHandlingID) real_t getAverageDensityInSystem(const shared_ptr<StructuredBlockStorage> & blocks, BlockDataID pdfFieldID, BlockDataID boundaryHandlingID)
{ {
...@@ -348,6 +390,41 @@ real_t getAverageDensityInSystem(const shared_ptr<StructuredBlockStorage> & bloc ...@@ -348,6 +390,41 @@ real_t getAverageDensityInSystem(const shared_ptr<StructuredBlockStorage> & bloc
return totalMass/real_c(count); return totalMass/real_c(count);
} }
template< typename BoundaryHandling_T>
void evaluateMapping(const shared_ptr<StructuredBlockStorage> & blocks, BlockDataID boundaryHandlingID,
uint_t & countSolidCells, uint_t & countFluidSolidLinks)
{
countSolidCells = uint_t(0);
countFluidSolidLinks = uint_t(0);
for( auto & block: *blocks)
{
auto boundaryHandling = block.getData<BoundaryHandling_T>(boundaryHandlingID);
auto flagField = boundaryHandling->getFlagField();
//auto innerDomain = flagField->xyzSize();
WALBERLA_FOR_ALL_CELLS_XYZ(flagField,
if(!boundaryHandling->isDomain(x,y,z) )
{
++countSolidCells;
} else
{
for( auto neighborDir = stencil::D3Q27::beginNoCenter(); neighborDir != stencil::D3Q27::end(); ++neighborDir )
{
Cell neighbor( x + neighborDir.cx(), y + neighborDir.cy(), z + neighborDir.cz() );
// check if neighbor cell is a solid cell
if( !boundaryHandling->isDomain( neighbor ) )
{
++countFluidSolidLinks;
}
}
}
);
}
mpi::reduceInplace(countSolidCells, mpi::SUM);
mpi::reduceInplace(countFluidSolidLinks, mpi::SUM);
}
////////// //////////
// MAIN // // MAIN //
////////// //////////
...@@ -389,13 +466,15 @@ int main( int argc, char **argv ) ...@@ -389,13 +466,15 @@ int main( int argc, char **argv )
std::string baseFolder = "vtk_out_MovingWithPrescribedVelocity"; std::string baseFolder = "vtk_out_MovingWithPrescribedVelocity";
std::string fileNameEnding = ""; std::string fileNameEnding = "";
bool logDensity = true; bool logDensity = true;
bool logMapping= true;
bool logFluidVelocity = true;
bool vtkOutputAtEnd = true; bool vtkOutputAtEnd = true;
//numerical parameters //numerical parameters
bool averageForceTorqueOverTwoTimeSteps = true;
bool conserveMomentum = false; bool conserveMomentum = false;
uint_t numRPDSubCycles = uint_t(1);
std::string reconstructorType = "Grad"; // Eq, EAN, Ext, Grad std::string reconstructorType = "Grad"; // Eq, EAN, Ext, Grad
std::string forceTorqueAveragingType = "runningAvg"; // noAvg, runningAvg, twoLBM,
std::string boundaryCondition = "CLI"; // SBB, CLI
real_t bulkViscRateFactor = real_t(1); real_t bulkViscRateFactor = real_t(1);
real_t magicNumber = real_t(0.1875); real_t magicNumber = real_t(0.1875);
bool useOmegaBulkAdaption = false; bool useOmegaBulkAdaption = false;
...@@ -409,21 +488,27 @@ int main( int argc, char **argv ) ...@@ -409,21 +488,27 @@ int main( int argc, char **argv )
real_t velocity = real_t(0.02); real_t velocity = real_t(0.02);
real_t Re = real_t(164); real_t Re = real_t(164);
real_t diameter = real_t(20); real_t diameter = real_t(20);
real_t numberOfPasses = real_t(3); real_t numberOfTimeStepsNonDim = real_t(20);
real_t domainWidthNonDim = real_t(4); real_t domainWidthNonDim = real_t(4);
real_t domainHeightNonDim = real_t(8); real_t domainHeightNonDim = real_t(8);
real_t accelerationFactor = real_t(3);
bool usePeriodicSetup = false; bool usePeriodicSetup = false;
bool artificiallyAccelerateSphere = false; bool artificiallyAccelerateSphere = false;
bool fixedSphere = false;
bool useGalileanInvariantSetup = false;
for( int i = 1; i < argc; ++i ) for( int i = 1; i < argc; ++i )
{ {
if( std::strcmp( argv[i], "--funcTest" ) == 0 ) { funcTest = true; continue; } if( std::strcmp( argv[i], "--funcTest" ) == 0 ) { funcTest = true; continue; }
if( std::strcmp( argv[i], "--noLogging" ) == 0 ) { fileIO = false; continue; } if( std::strcmp( argv[i], "--noLogging" ) == 0 ) { fileIO = false; continue; }
if( std::strcmp( argv[i], "--noVtkOutputAtEnd" ) == 0 ) { vtkOutputAtEnd = false; continue; } if( std::strcmp( argv[i], "--noVtkOutputAtEnd" ) == 0 ) { vtkOutputAtEnd = false; continue; }
if( std::strcmp( argv[i], "--logDensity" ) == 0 ) { logDensity = true; continue; } if( std::strcmp( argv[i], "--noLogDensity" ) == 0 ) { logDensity = false; continue; }
if( std::strcmp( argv[i], "--noLogMapping" ) == 0 ) { logMapping = false; continue; }
if( std::strcmp( argv[i], "--noLogFluidVelocity" ) == 0 ) { logFluidVelocity = false; continue; }
if( std::strcmp( argv[i], "--vtkIOFreq" ) == 0 ) { vtkIOFreq = uint_c( std::atof( argv[++i] ) ); continue; } if( std::strcmp( argv[i], "--vtkIOFreq" ) == 0 ) { vtkIOFreq = 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], "--averagingType" ) == 0 ) { forceTorqueAveragingType = argv[++i]; continue; }
if( std::strcmp( argv[i], "--noForceAveraging" ) == 0 ) { averageForceTorqueOverTwoTimeSteps = false; continue; } if( std::strcmp( argv[i], "--boundaryCondition" ) == 0 ) { boundaryCondition = argv[++i]; continue; }
if( std::strcmp( argv[i], "--conserveMomentum" ) == 0 ) { conserveMomentum = true; continue; } if( std::strcmp( argv[i], "--conserveMomentum" ) == 0 ) { conserveMomentum = true; continue; }
if( std::strcmp( argv[i], "--baseFolder" ) == 0 ) { baseFolder = argv[++i]; continue; } if( std::strcmp( argv[i], "--baseFolder" ) == 0 ) { baseFolder = argv[++i]; continue; }
if( std::strcmp( argv[i], "--reconstructorType" ) == 0 ) { reconstructorType = argv[++i]; continue; } if( std::strcmp( argv[i], "--reconstructorType" ) == 0 ) { reconstructorType = argv[++i]; continue; }
...@@ -437,16 +522,24 @@ int main( int argc, char **argv ) ...@@ -437,16 +522,24 @@ int main( int argc, char **argv )
if( std::strcmp( argv[i], "--fileName" ) == 0 ) { fileNameEnding = argv[++i]; continue; } if( std::strcmp( argv[i], "--fileName" ) == 0 ) { fileNameEnding = argv[++i]; continue; }
if( std::strcmp( argv[i], "--useOmegaBulkAdaption" ) == 0 ) { useOmegaBulkAdaption = true; continue; } if( std::strcmp( argv[i], "--useOmegaBulkAdaption" ) == 0 ) { useOmegaBulkAdaption = true; continue; }
if( std::strcmp( argv[i], "--adaptionLayerSize" ) == 0 ) { adaptionLayerSize = real_c(std::atof( argv[++i] )); continue; } if( std::strcmp( argv[i], "--adaptionLayerSize" ) == 0 ) { adaptionLayerSize = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--numberOfPasses" ) == 0 ) { numberOfPasses = real_c(std::atof( argv[++i] )); continue; } if( std::strcmp( argv[i], "--numberOfTimeStepsNonDim" ) == 0 ) { numberOfTimeStepsNonDim = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--initialSpherePosition" )== 0 ) { initialSpherePosition = real_c(std::atof( argv[++i] )); continue; } if( std::strcmp( argv[i], "--initialSpherePosition" )== 0 ) { initialSpherePosition = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--usePeriodicSetup" ) == 0 ) { usePeriodicSetup = true; continue; } if( std::strcmp( argv[i], "--usePeriodicSetup" ) == 0 ) { usePeriodicSetup = true; continue; }
if( std::strcmp( argv[i], "--artificiallyAccelerateSphere" ) == 0 ) { artificiallyAccelerateSphere = true; continue; } if( std::strcmp( argv[i], "--artificiallyAccelerateSphere" ) == 0 ) { artificiallyAccelerateSphere = true; continue; }
if( std::strcmp( argv[i], "--accelerationFactor" ) == 0 ) { accelerationFactor = real_c(std::atof( argv[++i] )); continue; }
if( std::strcmp( argv[i], "--blocksInX" ) == 0 ) { blocksInX = uint_c( std::atof( argv[++i] ) ); continue; } if( std::strcmp( argv[i], "--blocksInX" ) == 0 ) { blocksInX = uint_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--blocksInY" ) == 0 ) { blocksInY = uint_c( std::atof( argv[++i] ) ); continue; } if( std::strcmp( argv[i], "--blocksInY" ) == 0 ) { blocksInY = uint_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--blocksInZ" ) == 0 ) { blocksInZ = uint_c( std::atof( argv[++i] ) ); continue; } if( std::strcmp( argv[i], "--blocksInZ" ) == 0 ) { blocksInZ = uint_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--useGalileanInvariantSetup" ) == 0 ) { useGalileanInvariantSetup = true; continue; }
if( std::strcmp( argv[i], "--fixedSphere" ) == 0 ) { fixedSphere = true; continue; }
WALBERLA_ABORT("Unrecognized command line argument found: " << argv[i]); WALBERLA_ABORT("Unrecognized command line argument found: " << argv[i]);
} }
if(forceTorqueAveragingType != "noAvg" && forceTorqueAveragingType != "runningAvg" && forceTorqueAveragingType != "twoLBM")
{
WALBERLA_ABORT("Averaging type not implemented: " << forceTorqueAveragingType);
}
if( funcTest ) if( funcTest )
{ {
walberla::logging::Logging::instance()->setLogLevel(logging::Logging::LogLevel::WARNING); walberla::logging::Logging::instance()->setLogLevel(logging::Logging::LogLevel::WARNING);
...@@ -468,24 +561,33 @@ int main( int argc, char **argv ) ...@@ -468,24 +561,33 @@ int main( int argc, char **argv )
const real_t omega = lbm::collision_model::omegaFromViscosity(viscosity); const real_t omega = lbm::collision_model::omegaFromViscosity(viscosity);
const real_t relaxationTime = real_t(1) / omega; const real_t relaxationTime = real_t(1) / omega;
const real_t omegaBulk = lbm_mesapd_coupling::omegaBulkFromOmega(omega, bulkViscRateFactor); const real_t omegaBulk = lbm_mesapd_coupling::omegaBulkFromOmega(omega, bulkViscRateFactor);
const real_t tref = diameter / velocity;
const Vector3<uint_t> domainSize( uint_c(domainWidthNonDim * diameter), uint_c(domainWidthNonDim * diameter), uint_c(domainHeightNonDim * diameter) ); const Vector3<uint_t> domainSize( uint_c(domainWidthNonDim * diameter), uint_c(domainWidthNonDim * diameter), uint_c(domainHeightNonDim * diameter) );
WALBERLA_LOG_INFO_ON_ROOT("Setup (in simulation, i.e. lattice, units):"); 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(" - domain size = " << domainSize);
WALBERLA_LOG_INFO_ON_ROOT(" - Re = " << Re); WALBERLA_LOG_INFO_ON_ROOT(" - Re = " << Re);
WALBERLA_LOG_INFO_ON_ROOT(" - tref = " << tref);
WALBERLA_LOG_INFO_ON_ROOT(" - acceleration factor (a_acc) = " << accelerationFactor);
WALBERLA_LOG_INFO_ON_ROOT(" - sphere: diameter = " << diameter << ", constant velocity = " << velocity ); WALBERLA_LOG_INFO_ON_ROOT(" - sphere: diameter = " << diameter << ", constant velocity = " << velocity );
WALBERLA_LOG_INFO_ON_ROOT(" - relaxation time (tau) = " << relaxationTime << ", omega = " << omega << " kin. visc = " << viscosity ); WALBERLA_LOG_INFO_ON_ROOT(" - relaxation time (tau) = " << relaxationTime << ", omega = " << omega << ", kin. visc = " << viscosity );
WALBERLA_LOG_INFO_ON_ROOT(" - magic number " << magicNumber); WALBERLA_LOG_INFO_ON_ROOT(" - magic number " << magicNumber);
WALBERLA_LOG_INFO_ON_ROOT(" - omegaBulk = " << omegaBulk << ", bulk visc. = " << lbm_mesapd_coupling::bulkViscosityFromOmegaBulk(omegaBulk) << " (bvrf " << bulkViscRateFactor << ")"); WALBERLA_LOG_INFO_ON_ROOT(" - omegaBulk = " << omegaBulk << ", bulk visc. = " << lbm_mesapd_coupling::bulkViscosityFromOmegaBulk(omegaBulk) << " (bvrf " << bulkViscRateFactor << ")");
WALBERLA_LOG_INFO_ON_ROOT(" - use omega bulk adaption = " << useOmegaBulkAdaption << " (adaption layer size = " << adaptionLayerSize << ")"); WALBERLA_LOG_INFO_ON_ROOT(" - use omega bulk adaption = " << useOmegaBulkAdaption << " (adaption layer size = " << adaptionLayerSize << ")");
WALBERLA_LOG_INFO_ON_ROOT(" - reconstructor type = " << reconstructorType ); WALBERLA_LOG_INFO_ON_ROOT(" - reconstructor type = " << reconstructorType );
if( vtkIOFreq > 0 ) if( vtkIOFreq > 0 ) WALBERLA_LOG_INFO_ON_ROOT(" - writing vtk files to folder \"" << baseFolder << "\" with frequency " << vtkIOFreq);
{ if( useGalileanInvariantSetup ) {
WALBERLA_LOG_INFO_ON_ROOT(" - writing vtk files to folder \"" << baseFolder << "\" with frequency " << vtkIOFreq); WALBERLA_LOG_INFO_ON_ROOT("ATTENTION: Using Galilean Invariant Setup = Moving side walls, fixed sphere");
artificiallyAccelerateSphere = false;
usePeriodicSetup = false;
conserveMomentum = false;
fixedSphere = true;
} }
/////////////////////////// ///////////////////////////
// BLOCK STRUCTURE SETUP // // BLOCK STRUCTURE SETUP //
/////////////////////////// ///////////////////////////
...@@ -530,6 +632,7 @@ int main( int argc, char **argv ) ...@@ -530,6 +632,7 @@ int main( int argc, char **argv )
// bounding planes // bounding planes
if(!usePeriodicSetup) createPlaneSetup(ps,ss,blocks->getDomain()); if(!usePeriodicSetup) createPlaneSetup(ps,ss,blocks->getDomain());
if(useGalileanInvariantSetup) createPlaneSetup(ps,ss,blocks->getDomain(),-velocity);
// create sphere and store Uid // create sphere and store Uid
Vector3<real_t> initialPosition( real_t(0.5) * real_c(domainSize[0]), real_t(0.5) * real_c(domainSize[1]), initialSpherePosition * real_c(domainSize[2])); Vector3<real_t> initialPosition( real_t(0.5) * real_c(domainSize[0]), real_t(0.5) * real_c(domainSize[1]), initialSpherePosition * real_c(domainSize[2]));
...@@ -544,7 +647,7 @@ int main( int argc, char **argv ) ...@@ -544,7 +647,7 @@ int main( int argc, char **argv )
p.setInteractionRadius(diameter * real_t(0.5)); p.setInteractionRadius(diameter * real_t(0.5));
p.setOwner(mpi::MPIManager::instance()->rank()); p.setOwner(mpi::MPIManager::instance()->rank());
p.setShapeID(sphereShape); p.setShapeID(sphereShape);
p.setLinearVelocity(Vector3<real_t>(0.0,0.0,velocity)); if(!useGalileanInvariantSetup && !artificiallyAccelerateSphere) p.setLinearVelocity(Vector3<real_t>(real_t(0),real_t(0),velocity));
sphereUid = p.getUid(); sphereUid = p.getUid();
} }
mpi::allReduceInplace(sphereUid, mpi::SUM); mpi::allReduceInplace(sphereUid, mpi::SUM);
...@@ -554,7 +657,7 @@ int main( int argc, char **argv ) ...@@ -554,7 +657,7 @@ int main( int argc, char **argv )
//////////////////////// ////////////////////////
// add omega bulk field // 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 // create the lattice model
real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega ); real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega );
...@@ -569,7 +672,8 @@ int main( int argc, char **argv ) ...@@ -569,7 +672,8 @@ int main( int argc, char **argv )
// add PDF field // add PDF field
BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (fzyx)", latticeModel, BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (fzyx)", latticeModel,
Vector3< real_t >( real_t(0) ), real_t(1), useGalileanInvariantSetup ? Vector3<real_t>(real_t(0), real_t(0), -velocity) : Vector3< real_t >( real_t(0)),
real_t(1),
uint_t(1), field::fzyx ); uint_t(1), field::fzyx );
// add flag field // add flag field
BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>( blocks, "flag field" ); BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>( blocks, "flag field" );
...@@ -584,9 +688,13 @@ int main( int argc, char **argv ) ...@@ -584,9 +688,13 @@ int main( int argc, char **argv )
// interpolation functionality // interpolation functionality
using DensityAdaptor_T = typename lbm::Adaptor< LatticeModel_T >::Density; using DensityAdaptor_T = typename lbm::Adaptor< LatticeModel_T >::Density;
BlockDataID densityAdaptorID = field::addFieldAdaptor< DensityAdaptor_T >( blocks, pdfFieldID, "density adaptor" ); BlockDataID densityAdaptorID = field::addFieldAdaptor< DensityAdaptor_T >( blocks, pdfFieldID, "density adaptor" );
using VelocityAdaptor_T = typename lbm::Adaptor< LatticeModel_T >::VelocityVector;
BlockDataID velocityAdaptorID = field::addFieldAdaptor< VelocityAdaptor_T >( blocks, pdfFieldID, "velocity adaptor" );
using DensityInterpolator_T = typename field::TrilinearFieldInterpolator<DensityAdaptor_T, FlagField_T>; using DensityInterpolator_T = typename field::TrilinearFieldInterpolator<DensityAdaptor_T, FlagField_T>;
BlockDataID densityInterpolatorID = field::addFieldInterpolator< DensityInterpolator_T, FlagField_T >( blocks, densityAdaptorID, flagFieldID, Fluid_Flag ); BlockDataID densityInterpolatorID = field::addFieldInterpolator< DensityInterpolator_T, FlagField_T >( blocks, densityAdaptorID, flagFieldID, Fluid_Flag );
using VelocityInterpolator_T = typename field::TrilinearFieldInterpolator<VelocityAdaptor_T, FlagField_T>;
BlockDataID velocityInterpolatorID = field::addFieldInterpolator< VelocityInterpolator_T, FlagField_T >( blocks, velocityAdaptorID, flagFieldID, Fluid_Flag );
// set up RPD functionality // set up RPD functionality
std::function<void(void)> syncCall = [ps,rpdDomain, adaptionLayerSize ](){ std::function<void(void)> syncCall = [ps,rpdDomain, adaptionLayerSize ](){
...@@ -597,7 +705,6 @@ int main( int argc, char **argv ) ...@@ -597,7 +705,6 @@ int main( int argc, char **argv )
syncCall(); syncCall();
mesa_pd::kernel::ExplicitEulerWithShape explEulerIntegrator(real_t(1)/real_t(numRPDSubCycles));
mesa_pd::mpi::ReduceProperty reduceProperty; mesa_pd::mpi::ReduceProperty reduceProperty;
...@@ -615,10 +722,18 @@ int main( int argc, char **argv ) ...@@ -615,10 +722,18 @@ int main( int argc, char **argv )
/////////////// ///////////////
// map planes into the LBM simulation -> act as no-slip boundaries // map planes into the LBM simulation -> act as no-slip boundaries
ps->forEachParticle(false, lbm_mesapd_coupling::GlobalParticlesSelector(), *accessor, particleMappingKernel, *accessor, NoSlip_Flag); //ps->forEachParticle(false, lbm_mesapd_coupling::GlobalParticlesSelector(), *accessor, particleMappingKernel, *accessor, NoSlip_Flag);
// map particles into the LBM simulation // map particles into the LBM simulation
ps->forEachParticle(false, sphereSelector, *accessor, movingParticleMappingKernel, *accessor, MO_Flag); if(boundaryCondition == "SBB")
{
ps->forEachParticle(false, mesa_pd::kernel::SelectAll(), *accessor, movingParticleMappingKernel, *accessor, MO_SBB_Flag);
} else if(boundaryCondition == "CLI")
{
ps->forEachParticle(false, mesa_pd::kernel::SelectAll(), *accessor, movingParticleMappingKernel, *accessor, MO_CLI_Flag);
} else{
WALBERLA_ABORT("Unknown boundary condition " << boundaryCondition);
}
// setup of the LBM communication for synchronizing the pdf field between neighboring blocks // setup of the LBM communication for synchronizing the pdf field between neighboring blocks
...@@ -631,7 +746,10 @@ int main( int argc, char **argv ) ...@@ -631,7 +746,10 @@ int main( int argc, char **argv )
// create the timeloop // create the timeloop
const uint_t timesteps = uint_c(numberOfPasses * real_c(domainSize[2]) / velocity); real_t dtSim = real_t(1);
if (forceTorqueAveragingType == "twoLBM") dtSim = real_t(2);
const uint_t timesteps = uint_t( numberOfTimeStepsNonDim * tref);
WALBERLA_LOG_INFO_ON_ROOT("Running for " << timesteps << " time steps"); WALBERLA_LOG_INFO_ON_ROOT("Running for " << timesteps << " time steps");
SweepTimeloop timeloop( blocks->getBlockStorage(), timesteps ); SweepTimeloop timeloop( blocks->getBlockStorage(), timesteps );
...@@ -645,6 +763,7 @@ int main( int argc, char **argv ) ...@@ -645,6 +763,7 @@ int main( int argc, char **argv )
auto particleVtkOutput = make_shared<mesa_pd::vtk::ParticleVtkOutput>(ps); auto particleVtkOutput = make_shared<mesa_pd::vtk::ParticleVtkOutput>(ps);
particleVtkOutput->addOutput<mesa_pd::data::SelectParticleOwner>("owner"); particleVtkOutput->addOutput<mesa_pd::data::SelectParticleOwner>("owner");
particleVtkOutput->addOutput<mesa_pd::data::SelectParticleLinearVelocity>("velocity"); particleVtkOutput->addOutput<mesa_pd::data::SelectParticleLinearVelocity>("velocity");
particleVtkOutput->setParticleSelector( [sphereShape](const mesa_pd::data::ParticleStorage::iterator& pIt) {return pIt->getShapeID() == sphereShape;} ); //limit output to sphere
auto particleVtkWriter = vtk::createVTKOutput_PointData(particleVtkOutput, "Particles", vtkIOFreq, baseFolder, "simulation_step"); auto particleVtkWriter = vtk::createVTKOutput_PointData(particleVtkOutput, "Particles", vtkIOFreq, baseFolder, "simulation_step");
timeloop.addFuncBeforeTimeStep( vtk::writeFiles( particleVtkWriter ), "VTK (sphere data)" ); timeloop.addFuncBeforeTimeStep( vtk::writeFiles( particleVtkWriter ), "VTK (sphere data)" );
...@@ -667,37 +786,86 @@ int main( int argc, char **argv ) ...@@ -667,37 +786,86 @@ int main( int argc, char **argv )
} }
// add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
auto bhSweep = BoundaryHandling_T::getBlockSweep( boundaryHandlingID );
timeloop.add() << BeforeFunction( optimizedPDFCommunicationScheme, "LBM Communication" )
<< Sweep(bhSweep, "Boundary Handling" );
// stream + collide LBM step
#ifdef WALBERLA_BUILD_WITH_CODEGEN
auto lbmSweep = LatticeModel_T::Sweep( pdfFieldID );
//timeloop.add() << Sweep( [&lbmSweep](IBlock * const block){lbmSweep.streamCollide(block);}, "LB sweep" );
timeloop.add() << Sweep(lbmSweep, "LB sweep" );
#else
auto lbmSweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag );
timeloop.add() << Sweep( makeSharedSweep( lbmSweep ), "cell-wise LB sweep" );
#endif
SweepTimeloop timeloopAfterParticles( blocks->getBlockStorage(), timesteps );
// sweep for updating the particle mapping into the LBM simulation // sweep for updating the particle mapping into the LBM simulation
timeloop.add() << Sweep( lbm_mesapd_coupling::makeMovingParticleMapping<PdfField_T, BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, MO_Flag, FormerMO_Flag, sphereSelector, conserveMomentum), "Particle Mapping" ); if(boundaryCondition == "SBB")
{
timeloopAfterParticles.add() << Sweep( lbm_mesapd_coupling::makeMovingParticleMapping<PdfField_T, BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, MO_SBB_Flag, FormerMO_Flag, sphereSelector, conserveMomentum), "Particle Mapping" );
} else if(boundaryCondition == "CLI")
{
timeloopAfterParticles.add() << Sweep( lbm_mesapd_coupling::makeMovingParticleMapping<PdfField_T, BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, MO_CLI_Flag, FormerMO_Flag, sphereSelector, conserveMomentum), "Particle Mapping" );
} else{
WALBERLA_ABORT("Unknown boundary condition " << boundaryCondition);
}
bool useValuesFromGhostLayerForReconstruction = true; // true: full sync needed
// sweep for restoring PDFs in cells previously occupied by particles // sweep for restoring PDFs in cells previously occupied by particles
if( reconstructorType == "EAN") if( reconstructorType == "EAN")
{ {
auto sphereNormalExtrapolationDirectionFinder = make_shared<lbm_mesapd_coupling::SphereNormalExtrapolationDirectionFinder>(blocks); //auto sphereNormalExtrapolationDirectionFinder = make_shared<lbm_mesapd_coupling::SphereNormalExtrapolationDirectionFinder>(blocks);
auto equilibriumAndNonEquilibriumSphereNormalReconstructor = lbm_mesapd_coupling::makeEquilibriumAndNonEquilibriumReconstructor<BoundaryHandling_T>(blocks, boundaryHandlingID, sphereNormalExtrapolationDirectionFinder, uint_t(3), true); auto sphereNormalExtrapolationDirectionFinder = make_shared<lbm_mesapd_coupling::FlagFieldNormalExtrapolationDirectionFinder<BoundaryHandling_T> >(blocks,boundaryHandlingID);
auto equilibriumAndNonEquilibriumSphereNormalReconstructor = lbm_mesapd_coupling::makeEquilibriumAndNonEquilibriumReconstructor<BoundaryHandling_T>(blocks, boundaryHandlingID, sphereNormalExtrapolationDirectionFinder, uint_t(3), useValuesFromGhostLayerForReconstruction);
auto reconstructionManager = lbm_mesapd_coupling::makePdfReconstructionManager<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, FormerMO_Flag, Fluid_Flag, equilibriumAndNonEquilibriumSphereNormalReconstructor, conserveMomentum); auto reconstructionManager = lbm_mesapd_coupling::makePdfReconstructionManager<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, FormerMO_Flag, Fluid_Flag, equilibriumAndNonEquilibriumSphereNormalReconstructor, conserveMomentum);
timeloop.add() << BeforeFunction( fullPDFCommunicationScheme, "PDF Communication" ) timeloopAfterParticles.add() << BeforeFunction( fullPDFCommunicationScheme, "PDF Communication" )
<< Sweep( makeSharedSweep(reconstructionManager), "PDF Restore" ); << Sweep( makeSharedSweep(reconstructionManager), "PDF Restore" );
} else if( reconstructorType == "Grad") } else if( reconstructorType == "Grad")
{ {
bool recomputeDensity = false; bool recomputeDensity = false;
auto gradReconstructor = lbm_mesapd_coupling::makeGradsMomentApproximationReconstructor<BoundaryHandling_T>(blocks, boundaryHandlingID, omega, recomputeDensity, true, true); bool useCentralDifferences = true;
auto gradReconstructor = lbm_mesapd_coupling::makeGradsMomentApproximationReconstructor<BoundaryHandling_T>(blocks, boundaryHandlingID, omega, recomputeDensity, useCentralDifferences, useValuesFromGhostLayerForReconstruction);
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" );
} else if( reconstructorType == "GradDen")
{
bool recomputeDensity = true;
bool useCentralDifferences = true;
auto gradReconstructor = lbm_mesapd_coupling::makeGradsMomentApproximationReconstructor<BoundaryHandling_T>(blocks, boundaryHandlingID, omega, recomputeDensity, useCentralDifferences, useValuesFromGhostLayerForReconstruction);
timeloop.add() << BeforeFunction( fullPDFCommunicationScheme, "PDF Communication" ) 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" ); << Sweep( makeSharedSweep(lbm_mesapd_coupling::makePdfReconstructionManager<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, FormerMO_Flag, Fluid_Flag, gradReconstructor, conserveMomentum) ), "PDF Restore" );
} else if( reconstructorType == "Eq") } else if( reconstructorType == "Eq")
{ {
timeloop.add() << Sweep( makeSharedSweep(lbm_mesapd_coupling::makePdfReconstructionManager<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, FormerMO_Flag, Fluid_Flag, conserveMomentum) ), "PDF Restore" ); auto eqReconstructor = lbm_mesapd_coupling::makeEquilibriumReconstructor<BoundaryHandling_T>(blocks, boundaryHandlingID, useValuesFromGhostLayerForReconstruction);
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, eqReconstructor, conserveMomentum) ), "PDF Restore" );
} else if( reconstructorType == "Ext") } else if( reconstructorType == "Ext")
{ {
auto sphereNormalExtrapolationDirectionFinder = make_shared<lbm_mesapd_coupling::SphereNormalExtrapolationDirectionFinder>(blocks); // gets easily unstable
auto extrapolationReconstructor = lbm_mesapd_coupling::makeExtrapolationReconstructor<BoundaryHandling_T, lbm_mesapd_coupling::SphereNormalExtrapolationDirectionFinder, true>(blocks, boundaryHandlingID, sphereNormalExtrapolationDirectionFinder, uint_t(3), true); //auto sphereNormalExtrapolationDirectionFinder = make_shared<lbm_mesapd_coupling::SphereNormalExtrapolationDirectionFinder>(blocks);
timeloop.add() << BeforeFunction( fullPDFCommunicationScheme, "LBM Communication" ) auto sphereNormalExtrapolationDirectionFinder = make_shared<lbm_mesapd_coupling::FlagFieldNormalExtrapolationDirectionFinder<BoundaryHandling_T> >(blocks,boundaryHandlingID);
<< Sweep( makeSharedSweep(lbm_mesapd_coupling::makePdfReconstructionManager<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, FormerMO_Flag, Fluid_Flag, extrapolationReconstructor, conserveMomentum)), "PDF Restore" ); auto extrapolationReconstructor = lbm_mesapd_coupling::makeExtrapolationReconstructor<BoundaryHandling_T, lbm_mesapd_coupling::FlagFieldNormalExtrapolationDirectionFinder<BoundaryHandling_T>, false>(blocks, boundaryHandlingID, sphereNormalExtrapolationDirectionFinder, uint_t(3), useValuesFromGhostLayerForReconstruction);
} else { 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 if( reconstructorType == "ExtNS")
{
//auto sphereNormalExtrapolationDirectionFinder = make_shared<lbm_mesapd_coupling::SphereNormalExtrapolationDirectionFinder>(blocks);
auto sphereNormalExtrapolationDirectionFinder = make_shared<lbm_mesapd_coupling::FlagFieldNormalExtrapolationDirectionFinder<BoundaryHandling_T> >(blocks,boundaryHandlingID);
auto extrapolationReconstructor = lbm_mesapd_coupling::makeExtrapolationReconstructor<BoundaryHandling_T, lbm_mesapd_coupling::FlagFieldNormalExtrapolationDirectionFinder<BoundaryHandling_T>, true>(blocks, boundaryHandlingID, sphereNormalExtrapolationDirectionFinder, uint_t(3), useValuesFromGhostLayerForReconstruction);
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 {
WALBERLA_ABORT("Unknown reconstructor type " << reconstructorType); WALBERLA_ABORT("Unknown reconstructor type " << reconstructorType);
} }
...@@ -707,35 +875,37 @@ int main( int argc, char **argv ) ...@@ -707,35 +875,37 @@ int main( int argc, char **argv )
using OmegaBulkAdapter_T = lbm_mesapd_coupling::OmegaBulkAdapter<ParticleAccessor_T, decltype(sphereSelector)>; using OmegaBulkAdapter_T = lbm_mesapd_coupling::OmegaBulkAdapter<ParticleAccessor_T, decltype(sphereSelector)>;
real_t defaultOmegaBulk = lbm_mesapd_coupling::omegaBulkFromOmega(omega, real_t(1)); real_t defaultOmegaBulk = lbm_mesapd_coupling::omegaBulkFromOmega(omega, real_t(1));
shared_ptr<OmegaBulkAdapter_T> omegaBulkAdapter = make_shared<OmegaBulkAdapter_T>(blocks, omegaBulkFieldID, accessor, defaultOmegaBulk, omegaBulk, adaptionLayerSize, sphereSelector); shared_ptr<OmegaBulkAdapter_T> omegaBulkAdapter = make_shared<OmegaBulkAdapter_T>(blocks, omegaBulkFieldID, accessor, defaultOmegaBulk, omegaBulk, adaptionLayerSize, sphereSelector);
timeloop.add() << Sweep( makeSharedSweep(omegaBulkAdapter), "Omega Bulk Adapter"); timeloopAfterParticles.add() << Sweep( makeSharedSweep(omegaBulkAdapter), "Omega Bulk Adapter");
// initially adapt
for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt) {
(*omegaBulkAdapter)(blockIt.get());
}
} }
// add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
auto bhSweep = BoundaryHandling_T::getBlockSweep( boundaryHandlingID );
timeloop.add() << BeforeFunction( optimizedPDFCommunicationScheme, "LBM Communication" )
<< Sweep(bhSweep, "Boundary Handling" );
// stream + collide LBM step
#ifdef WALBERLA_BUILD_WITH_CODEGEN
auto lbmSweep = LatticeModel_T::Sweep( pdfFieldID );
timeloop.add() << Sweep( lbmSweep, "LB sweep" );
#else
auto lbmSweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag );
timeloop.add() << Sweep( makeSharedSweep( lbmSweep ), "cell-wise LB sweep" );
#endif
// evaluation functionality // evaluation functionality
std::string loggingFileName( baseFolder + "/LoggingPrescribedMovingSphere"); std::string fileNamePostFix("");
loggingFileName += "_Re" + std::to_string(uint_c(Re)); fileNamePostFix += "_Re" + std::to_string(uint_c(Re));
loggingFileName += "_D" + std::to_string(uint_c(diameter)); fileNamePostFix += "_D" + std::to_string(uint_c(diameter));
loggingFileName += "_vel" + std::to_string(float(velocity)); fileNamePostFix += "_vel" + std::to_string(float(velocity));
loggingFileName += "_recon" + reconstructorType; fileNamePostFix += "_" + boundaryCondition;
loggingFileName += "_bvrf" + std::to_string(uint_c(bulkViscRateFactor)); fileNamePostFix += "_recon" + reconstructorType;
loggingFileName += "_mn" + std::to_string(float(magicNumber)); fileNamePostFix += "_" + forceTorqueAveragingType;
if( useOmegaBulkAdaption ) loggingFileName += "_uOBA" + std::to_string(uint_c(adaptionLayerSize)); fileNamePostFix += "_bvrf" + std::to_string(uint_c(bulkViscRateFactor));
if( conserveMomentum ) loggingFileName += "_conserveMomentum"; //fileNamePostFix += "_mn" + std::to_string(float(magicNumber));
if( artificiallyAccelerateSphere ) loggingFileName += "_acc"; if( useOmegaBulkAdaption ) fileNamePostFix += "_uOBA" + std::to_string(uint_c(adaptionLayerSize));
if( !fileNameEnding.empty()) loggingFileName += "_" + fileNameEnding;
if( useGalileanInvariantSetup )
{
fileNamePostFix += "_galileanInvariant";
}
else{
if( conserveMomentum ) fileNamePostFix += "_conserveMomentum";
if( artificiallyAccelerateSphere ) fileNamePostFix += "_acc";
if( fixedSphere ) fileNamePostFix += "_fixed";
}
if( !fileNameEnding.empty()) fileNamePostFix += "_" + fileNameEnding;
std::string loggingFileName( baseFolder + "/Log" + fileNamePostFix);
loggingFileName += ".txt"; loggingFileName += ".txt";
if( fileIO ) if( fileIO )
...@@ -753,43 +923,41 @@ int main( int argc, char **argv ) ...@@ -753,43 +923,41 @@ int main( int argc, char **argv )
const bool useOpenMP = false; const bool useOpenMP = false;
const real_t densityRatio = real_t(7800) / real_t(935); mesa_pd::kernel::ExplicitEuler explEulerIntegrator(dtSim);
const real_t responseTime = densityRatio * diameter * diameter / ( real_t(18) * viscosity ); lbm_mesapd_coupling::InitializeHydrodynamicForceTorqueForAveragingKernel initializeHydrodynamicForceTorqueForAveragingKernel;
WALBERLA_LOG_INFO_ON_ROOT(" - response time " << responseTime);
const real_t accelerationFactor = real_t(1) / (real_t(0.1) * responseTime);
// time loop // time loop
for (uint_t i = 0; i < timesteps; ++i ) for (uint_t i = 0; i < uint_c(real_c(timesteps)/dtSim); ++i )
{ {
// perform a single simulation step -> this contains LBM and setting of the hydrodynamic interactions // perform a single simulation step -> this contains LBM and setting of the hydrodynamic interactions
timeloop.singleStep( timeloopTiming ); timeloop.singleStep( timeloopTiming );
timeloopTiming["RPD"].start(); if( forceTorqueAveragingType == "twoLBM")
reduceProperty.operator()<mesa_pd::HydrodynamicForceTorqueNotification>(*ps);
if( averageForceTorqueOverTwoTimeSteps )
{ {
if( i == 0 ) // store current force and torque in fOld and tOld
{ ps->forEachParticle(useOpenMP, sphereSelector, *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
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 );
}
reduceProperty.operator()<mesa_pd::ForceTorqueNotification>(*ps); // reset f and t
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, resetHydrodynamicForceTorque, *accessor );
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, explEulerIntegrator, *accessor);
syncCall();
timeloopTiming["RPD"].end(); timeloop.singleStep( timeloopTiming );
if( artificiallyAccelerateSphere ) // f = (f+fOld)/2
ps->forEachParticle(useOpenMP, sphereSelector, *accessor, averageHydrodynamicForceTorque, *accessor );
reduceProperty.operator()<mesa_pd::HydrodynamicForceTorqueNotification>(*ps);
} else
{ {
// accelerating after reading from a checkpoint file does not make sense, as current actual runtime is not known reduceProperty.operator()<mesa_pd::HydrodynamicForceTorqueNotification>(*ps);
real_t newSphereVel = velocity * (std::exp(-accelerationFactor * real_t(i) ) - real_t(1));
ps->forEachParticle(useOpenMP, sphereSelector, *accessor, [newSphereVel](const size_t idx, ParticleAccessor_T& ac){ ac.setLinearVelocity(idx, Vector3<real_t>(real_t(0), real_t(0), newSphereVel));}, *accessor); if( forceTorqueAveragingType == "runningAvg" )
{
if( i == 0 )
{
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
}
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, averageHydrodynamicForceTorque, *accessor );
}
} }
// evaluation // evaluation
...@@ -797,28 +965,68 @@ int main( int argc, char **argv ) ...@@ -797,28 +965,68 @@ int main( int argc, char **argv )
logger(); logger();
if( std::isnan(logger.getForce()) )
{
WALBERLA_ABORT("Nan found in force - aborting");
}
real_t density1 = real_t(0); real_t density1 = real_t(0);
real_t density2 = real_t(0); real_t density2 = real_t(0);
real_t averageDensity = real_t(0); real_t averageDensity = real_t(0);
auto movingProbePosition = Vector3<real_t>(real_c(domainSize[0])*real_t(0.75), real_c(domainSize[1])*real_t(0.5), logger.getPosition());
if(logDensity) if(logDensity)
{ {
density1 = getDensityAtPosition<DensityInterpolator_T >(blocks, densityInterpolatorID, Vector3<real_t>(real_c(domainSize[0])*real_t(0.25), real_c(domainSize[1])*real_t(0.5), real_c(domainSize[2])*real_t(0.5))); density1 = getDensityAtPosition<DensityInterpolator_T >(blocks, densityInterpolatorID, Vector3<real_t>(real_c(domainSize[0])*real_t(0.25), real_c(domainSize[1])*real_t(0.5), real_c(domainSize[2])*real_t(0.5)));
density2 = getDensityAtPosition<DensityInterpolator_T >(blocks, densityInterpolatorID, Vector3<real_t>(real_c(domainSize[0])*real_t(0.75), real_c(domainSize[1])*real_t(0.5), logger.getPosition())); density2 = getDensityAtPosition<DensityInterpolator_T >(blocks, densityInterpolatorID, movingProbePosition);
averageDensity = getAverageDensityInSystem<BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID); averageDensity = getAverageDensityInSystem<BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID);
} }
logger.writeToFile(i, density1, density2, averageDensity);
Vector3<real_t> velocityAtProbePos(real_t(0));
if(logFluidVelocity)
{
velocityAtProbePos = getVelocityAtPosition<VelocityInterpolator_T >(blocks, velocityInterpolatorID, movingProbePosition);
}
uint_t countSolidCells = 0;
uint_t countFluidSolidLinks = 0;
if(logMapping) evaluateMapping<BoundaryHandling_T>(blocks, boundaryHandlingID, countSolidCells, countFluidSolidLinks);
logger.writeToFile(i*uint_c(dtSim), density1, density2, averageDensity, countSolidCells, countFluidSolidLinks, velocityAtProbePos);
timeloopTiming["Logging"].end(); timeloopTiming["Logging"].end();
// reset after logging here // reset after logging here
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, resetHydrodynamicForceTorque, *accessor ); ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, resetHydrodynamicForceTorque, *accessor );
if(!fixedSphere)
{
timeloopTiming["RPD"].start();
if( artificiallyAccelerateSphere )
{
real_t newSphereVel = - velocity * (std::exp(- real_t(i) * dtSim * accelerationFactor/ tref ) - real_t(1));
ps->forEachParticle(useOpenMP, sphereSelector, *accessor, [newSphereVel](const size_t idx, ParticleAccessor_T& ac){ ac.setLinearVelocity(idx, Vector3<real_t>(real_t(0), real_t(0), newSphereVel));}, *accessor);
}
// there is no force acting on sphere
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, explEulerIntegrator, *accessor);
syncCall();
timeloopTiming["RPD"].end();
// update mapping, restore PDFs
timeloopAfterParticles.singleStep( timeloopTiming );
}
} }
timeloopTiming.logResultOnRoot(); timeloopTiming.logResultOnRoot();
if(vtkOutputAtEnd) if(vtkOutputAtEnd)
{ {
std::string vtkFileName = "fluid_field"; std::string vtkFileName = "fluid_field_final"+fileNamePostFix;
vtkFileName += "_bvrf" + std::to_string(uint_c(bulkViscRateFactor)); vtkFileName += "_bvrf" + std::to_string(uint_c(bulkViscRateFactor));
if( useOmegaBulkAdaption ) vtkFileName += "_uOBA" + std::to_string(uint_c(adaptionLayerSize)); if( useOmegaBulkAdaption ) vtkFileName += "_uOBA" + std::to_string(uint_c(adaptionLayerSize));
...@@ -843,6 +1051,13 @@ int main( int argc, char **argv ) ...@@ -843,6 +1051,13 @@ int main( int argc, char **argv )
pdfFieldVTK->addCellDataWriter( make_shared< lbm::DensityVTKWriter < LatticeModel_T, float > >( pdfFieldID, "DensityFromPDF" ) ); pdfFieldVTK->addCellDataWriter( make_shared< lbm::DensityVTKWriter < LatticeModel_T, float > >( pdfFieldID, "DensityFromPDF" ) );
vtk::writeFiles(pdfFieldVTK)(); vtk::writeFiles(pdfFieldVTK)();
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 pIt->getShapeID() == sphereShape;} ); //limit output to sphere
auto particleVtkWriter = vtk::createVTKOutput_PointData(particleVtkOutput, "Particles_final"+fileNamePostFix, uint_t(1), baseFolder);
vtk::writeFiles(particleVtkWriter)();
} }
......