From 4249d29ccf9fe60bdc6051541c4839669e0073b5 Mon Sep 17 00:00:00 2001
From: Lukas Werner <lks.werner@fau.de>
Date: Thu, 26 Aug 2021 08:56:21 +0000
Subject: [PATCH] Add virtual mass functionality to lbm-mesapd coupling
---
apps/CMakeLists.txt | 2 +-
apps/showcases/CMakeLists.txt | 1 +
.../CMakeLists.txt | 3 +
.../LightRisingParticleInFluidAMR.cpp | 1410 +++++++++++++++++
python/mesa_pd.py | 16 +
.../ParticlePresenceLevelDetermination.h | 57 +-
.../virtualmass/AddVirtualForceTorqueKernel.h | 88 +
.../virtualmass/InitializeVirtualMassKernel.h | 64 +
...ticleAccessorWithShapeVirtualMassWrapper.h | 58 +
src/mesa_pd/data/ParticleAccessor.h | 63 +
src/mesa_pd/data/ParticleAccessorWithShape.h | 3 +
src/mesa_pd/data/ParticleStorage.h | 203 +++
src/mesa_pd/mpi/notifications/ParseMessage.h | 7 +
.../notifications/ParticleCopyNotification.h | 28 +
.../ParticleMigrationNotification.h | 21 +
tests/lbm_mesapd_coupling/CMakeLists.txt | 3 +
.../utility/VirtualMass.cpp | 117 ++
17 files changed, 2118 insertions(+), 26 deletions(-)
create mode 100644 apps/showcases/LightRisingParticleInFluidAMR/CMakeLists.txt
create mode 100644 apps/showcases/LightRisingParticleInFluidAMR/LightRisingParticleInFluidAMR.cpp
create mode 100644 src/lbm_mesapd_coupling/utility/virtualmass/AddVirtualForceTorqueKernel.h
create mode 100644 src/lbm_mesapd_coupling/utility/virtualmass/InitializeVirtualMassKernel.h
create mode 100644 src/lbm_mesapd_coupling/utility/virtualmass/ParticleAccessorWithShapeVirtualMassWrapper.h
create mode 100644 tests/lbm_mesapd_coupling/utility/VirtualMass.cpp
diff --git a/apps/CMakeLists.txt b/apps/CMakeLists.txt
index 178380b57..965675d4d 100644
--- a/apps/CMakeLists.txt
+++ b/apps/CMakeLists.txt
@@ -32,4 +32,4 @@ endif()
# Python module
if ( WALBERLA_BUILD_WITH_PYTHON )
add_subdirectory( pythonmodule )
-endif()
\ No newline at end of file
+endif()
diff --git a/apps/showcases/CMakeLists.txt b/apps/showcases/CMakeLists.txt
index 1b23d4954..6b34ee252 100644
--- a/apps/showcases/CMakeLists.txt
+++ b/apps/showcases/CMakeLists.txt
@@ -2,6 +2,7 @@ add_subdirectory( BidisperseFluidizedBed )
add_subdirectory( CombinedResolvedUnresolved )
add_subdirectory( FluidizedBed )
add_subdirectory( HeatConduction )
+add_subdirectory( LightRisingParticleInFluidAMR )
add_subdirectory( Mixer )
add_subdirectory( PegIntoSphereBed )
if ( WALBERLA_BUILD_WITH_CODEGEN AND WALBERLA_BUILD_WITH_PYTHON )
diff --git a/apps/showcases/LightRisingParticleInFluidAMR/CMakeLists.txt b/apps/showcases/LightRisingParticleInFluidAMR/CMakeLists.txt
new file mode 100644
index 000000000..96fd54bdb
--- /dev/null
+++ b/apps/showcases/LightRisingParticleInFluidAMR/CMakeLists.txt
@@ -0,0 +1,3 @@
+waLBerla_add_executable ( NAME LightRisingParticleInFluidAMR
+ FILES LightRisingParticleInFluidAMR.cpp
+ DEPENDS core mesa_pd lbm lbm_mesapd_coupling domain_decomposition field vtk geometry postprocessing )
diff --git a/apps/showcases/LightRisingParticleInFluidAMR/LightRisingParticleInFluidAMR.cpp b/apps/showcases/LightRisingParticleInFluidAMR/LightRisingParticleInFluidAMR.cpp
new file mode 100644
index 000000000..6dcafd552
--- /dev/null
+++ b/apps/showcases/LightRisingParticleInFluidAMR/LightRisingParticleInFluidAMR.cpp
@@ -0,0 +1,1410 @@
+//======================================================================================================================
+//
+// 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 LightRisingParticleInFluidAMR.cpp
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+#include "blockforest/SetupBlockForest.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+#include "blockforest/communication/UniformToNonUniformPackInfoAdapter.h"
+#include "blockforest/loadbalancing/DynamicCurve.h"
+#include "blockforest/loadbalancing/StaticCurve.h"
+
+#include "boundary/all.h"
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/SharedFunctor.h"
+#include "core/debug/Debug.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/logging/all.h"
+#include "core/math/Random.h"
+#include "core/math/all.h"
+#include "core/mpi/Broadcast.h"
+#include "core/timing/RemainingTimeLogger.h"
+
+#include "domain_decomposition/BlockSweepWrapper.h"
+
+#include "field/AddToStorage.h"
+#include "field/StabilityChecker.h"
+#include "field/adaptors/AdaptorCreators.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/field/QCriterionFieldWriter.h"
+#include "lbm/field/VelocityFieldWriter.h"
+#include "lbm/lattice_model/D3Q19.h"
+#include "lbm/refinement/all.h"
+#include "lbm/sweeps/CellwiseSweep.h"
+#include "lbm/vtk/all.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+#include "vtk/all.h"
+
+#include <functional>
+
+#include "lbm_mesapd_coupling/DataTypes.h"
+#include "lbm_mesapd_coupling/amr/InfoCollection.h"
+#include "lbm_mesapd_coupling/amr/level_determination/ParticlePresenceLevelDetermination.h"
+#include "lbm_mesapd_coupling/mapping/ParticleMapping.h"
+#include "lbm_mesapd_coupling/momentum_exchange_method/MovingParticleMapping.h"
+#include "lbm_mesapd_coupling/momentum_exchange_method/boundary/CurvedLinear.h"
+#include "lbm_mesapd_coupling/momentum_exchange_method/reconstruction/ExtrapolationDirectionFinder.h"
+#include "lbm_mesapd_coupling/momentum_exchange_method/reconstruction/PdfReconstructionManager.h"
+#include "lbm_mesapd_coupling/momentum_exchange_method/reconstruction/Reconstructor.h"
+#include "lbm_mesapd_coupling/utility/AddForceOnParticlesKernel.h"
+#include "lbm_mesapd_coupling/utility/AddHydrodynamicInteractionKernel.h"
+#include "lbm_mesapd_coupling/utility/AverageHydrodynamicForceTorqueKernel.h"
+#include "lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h"
+#include "lbm_mesapd_coupling/utility/ParticleSelector.h"
+#include "lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
+#include "lbm_mesapd_coupling/utility/SubCyclingManager.h"
+#include "lbm_mesapd_coupling/utility/virtualmass/AddVirtualForceTorqueKernel.h"
+#include "lbm_mesapd_coupling/utility/virtualmass/InitializeVirtualMassKernel.h"
+#include "lbm_mesapd_coupling/utility/virtualmass/ParticleAccessorWithShapeVirtualMassWrapper.h"
+
+#include "mesa_pd/data/DataTypes.h"
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/ShapeStorage.h"
+#include "mesa_pd/data/shape/Sphere.h"
+#include "mesa_pd/domain/BlockForestDataHandling.h"
+#include "mesa_pd/domain/BlockForestDomain.h"
+#include "mesa_pd/kernel/ParticleSelector.h"
+#include "mesa_pd/kernel/VelocityVerlet.h"
+#include "mesa_pd/mpi/ClearGhostOwnerSync.h"
+#include "mesa_pd/mpi/ClearNextNeighborSync.h"
+#include "mesa_pd/mpi/ReduceProperty.h"
+#include "mesa_pd/mpi/SyncGhostOwners.h"
+#include "mesa_pd/mpi/SyncNextNeighbors.h"
+#include "mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h"
+#include "mesa_pd/vtk/ParticleVtkOutput.h"
+
+/**
+ * \brief Showcase for a stabilized fluid - very light particle simulation using virtual mass
+ * A light sphere, situated at the bottom of a cuboid domain, rises.
+ * Depending on the density of the sphere (--sphereDensity) and galileo number (--galileoNumber), the sphere
+ * then undergoes various movement patterns along its way to the top.
+ * Employs adaptive mesh refinement and checkpointing.
+ * Can output VTK files based on the q criterion to visualize vortices behind the rising sphere.
+ * Most default parameters should be sufficient and do not need explicit values.
+ * As used in the master thesis of Lukas Werner, 2020, and follow-up paper at http://doi.org/10.1002/fld.5034
+ */
+
+namespace light_rising_particle_amr {
+using namespace walberla;
+using walberla::uint_t;
+
+#ifdef BUILD_WITH_MRT
+using LatticeModel_T = lbm::D3Q19< lbm::collision_model::D3Q19MRT>;
+#else
+using LatticeModel_T = lbm::D3Q19< lbm::collision_model::TRT, false >;
+#endif
+
+using Stencil_T = LatticeModel_T::Stencil;
+using PdfField_T = lbm::PdfField<LatticeModel_T>;
+using ParticleField_T = lbm_mesapd_coupling::ParticleField_T;
+using VelocityField_T = GhostLayerField<Vector3<real_t>, 1>;
+using QCriterionField_T = GhostLayerField<real_t, 1>;
+
+using flag_t = walberla::uint8_t;
+using FlagField_T = FlagField<flag_t>;
+
+using ScalarField_T = GhostLayerField<real_t, 1>;
+
+const uint_t FieldGhostLayers = 4;
+
+using ParticleAccessor_T = mesa_pd::data::ParticleAccessorWithShape;
+
+using NoSlip_T = lbm::NoSlip<LatticeModel_T, flag_t> ;
+using MO_T = lbm_mesapd_coupling::CurvedLinear< LatticeModel_T, FlagField_T, ParticleAccessor_T >;
+using BoundaryHandling_T = BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T, MO_T >;
+
+
+//// flags
+
+const FlagUID Fluid_Flag("fluid");
+const FlagUID NoSlip_Flag("no slip");
+const FlagUID MO_Flag("moving obstacle");
+const FlagUID FormerMO_Flag("former moving obstacle");
+
+
+//// block structure
+
+static void refinementSelection(SetupBlockForest &forest, uint_t levels, AABB refinementBox) {
+ real_t dx = real_t(1); // dx on finest level
+ const uint_t finestLevel = levels - uint_t(1);
+ for (auto block = forest.begin(); block != forest.end(); ++block) {
+ uint_t blockLevel = block->getLevel();
+ uint_t levelScalingFactor = (uint_t(1) << (finestLevel - blockLevel));
+ real_t dxOnLevel = dx * real_c(levelScalingFactor);
+ AABB blockAABB = block->getAABB();
+ // extend block AABB by ghostlayers
+ AABB extendedBlockAABB = blockAABB.getExtended(dxOnLevel * real_c(FieldGhostLayers));
+ if (extendedBlockAABB.intersects(refinementBox))
+ if (blockLevel < finestLevel) // only if the block is not on the finest level
+ block->setMarker(true);
+ }
+}
+
+static void workloadAndMemoryAssignment(SetupBlockForest &forest) {
+ for (auto block = forest.begin(); block != forest.end(); ++block) {
+ block->setWorkload(numeric_cast<workload_t>(uint_t(1) << block->getLevel()));
+ block->setMemory(numeric_cast<memory_t>(1));
+ }
+}
+static shared_ptr<StructuredBlockForest>
+createBlockStructure(const AABB &domainAABB, Vector3<uint_t> blockSizeInCells, uint_t numberOfLevels, real_t diameter,
+ Vector3<real_t> spherePosition, Vector3<bool> periodicity,
+ bool readFromCheckPointFile, const std::string & forestCheckPointFileName,
+ bool keepGlobalBlockInformation = false) {
+ if (readFromCheckPointFile) {
+ WALBERLA_LOG_INFO_ON_ROOT("Reading block forest from file!");
+ return blockforest::createUniformBlockGrid(forestCheckPointFileName,
+ blockSizeInCells[0], blockSizeInCells[1], blockSizeInCells[2],
+ false);
+ } else {
+ SetupBlockForest sforest;
+ Vector3<uint_t> numberOfFineBlocksPerDirection(uint_c(domainAABB.size(0)) / blockSizeInCells[0],
+ uint_c(domainAABB.size(1)) / blockSizeInCells[1],
+ uint_c(domainAABB.size(2)) / blockSizeInCells[2]);
+ for (uint_t i = 0; i < 3; ++i) {
+ WALBERLA_CHECK_EQUAL(numberOfFineBlocksPerDirection[i] * blockSizeInCells[i], uint_c(domainAABB.size(i)),
+ "Domain can not be decomposed in direction " << i << " into fine blocks of size "
+ << blockSizeInCells[i]);
+ }
+
+ uint_t levelScalingFactor = (uint_t(1) << (numberOfLevels - uint_t(1)));
+ Vector3<uint_t> numberOfCoarseBlocksPerDirection(numberOfFineBlocksPerDirection / levelScalingFactor);
+ for (uint_t i = 0; i < 3; ++i) {
+ WALBERLA_CHECK_EQUAL(numberOfCoarseBlocksPerDirection[i] * levelScalingFactor,
+ numberOfFineBlocksPerDirection[i],
+ "Domain can not be refined in direction " << i
+ << " according to the specified number of levels!");
+ }
+ AABB refinementBox( std::floor(spherePosition[0] - real_t(0.5) * diameter),
+ std::floor(spherePosition[1] - real_t(0.5) * diameter),
+ std::floor(spherePosition[2] - real_t(0.5) * diameter),
+ std::ceil( spherePosition[0] + real_t(0.5) * diameter),
+ std::ceil( spherePosition[1] + real_t(0.5) * diameter),
+ std::ceil( spherePosition[2] + real_t(0.5) * diameter) );
+ WALBERLA_LOG_INFO_ON_ROOT(" - refinement box: " << refinementBox);
+
+ sforest.addRefinementSelectionFunction(
+ std::bind(refinementSelection, std::placeholders::_1, numberOfLevels, refinementBox));
+ sforest.addWorkloadMemorySUIDAssignmentFunction(workloadAndMemoryAssignment);
+
+ sforest.init(domainAABB, numberOfCoarseBlocksPerDirection[0], numberOfCoarseBlocksPerDirection[1],
+ numberOfCoarseBlocksPerDirection[2], periodicity[0], periodicity[1], periodicity[2]);
+
+ // calculate process distribution
+ const memory_t memoryLimit = math::Limits<memory_t>::inf();
+ sforest.balanceLoad(blockforest::StaticLevelwiseCurveBalance(true),
+ uint_c(MPIManager::instance()->numProcesses()), real_t(0), memoryLimit, true);
+
+ WALBERLA_LOG_INFO_ON_ROOT(sforest);
+
+ MPIManager::instance()->useWorldComm();
+
+ auto blockForest = make_shared<BlockForest>(uint_c(MPIManager::instance()->rank()), sforest,
+ keepGlobalBlockInformation);
+
+ // create StructuredBlockForest (encapsulates a newly created BlockForest)
+ shared_ptr<StructuredBlockForest> sbf =
+ make_shared<StructuredBlockForest>(
+ blockForest,
+ blockSizeInCells[0], blockSizeInCells[1], blockSizeInCells[2]);
+ sbf->createCellBoundingBoxes();
+
+ return sbf;
+ }
+}
+
+//// boundary handling
+
+class MyBoundaryHandling : public blockforest::AlwaysInitializeBlockDataHandling< BoundaryHandling_T > {
+public:
+ MyBoundaryHandling(const weak_ptr< StructuredBlockStorage > & storage, const BlockDataID& flagFieldID,
+ const BlockDataID& pdfFieldID, const BlockDataID& particleFieldID, const shared_ptr<ParticleAccessor_T>& ac) :
+ storage_(storage), flagFieldID_(flagFieldID), pdfFieldID_(pdfFieldID), particleFieldID_(particleFieldID),
+ ac_(ac) {
+
+ }
+
+ BoundaryHandling_T * initialize( IBlock * const block ) override;
+
+private:
+ weak_ptr< StructuredBlockStorage > storage_;
+
+ const BlockDataID flagFieldID_;
+ const BlockDataID pdfFieldID_;
+ const BlockDataID particleFieldID_;
+
+ shared_ptr<ParticleAccessor_T> ac_;
+};
+
+BoundaryHandling_T * MyBoundaryHandling::initialize( IBlock * const block )
+{
+ WALBERLA_ASSERT_NOT_NULLPTR( block );
+
+ auto* flagField = block->getData<FlagField_T>(flagFieldID_);
+ auto* pdfField = block->getData<PdfField_T>(pdfFieldID_);
+ auto* particleField = block->getData<lbm_mesapd_coupling::ParticleField_T>(particleFieldID_);
+
+ const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
+
+ auto storagePtr = storage_.lock();
+ WALBERLA_CHECK_NOT_NULLPTR( storagePtr );
+
+ BoundaryHandling_T * handling = new BoundaryHandling_T( "moving obstacle boundary handling", flagField, fluid,
+ NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
+ MO_T( "MO", MO_Flag, pdfField, flagField, particleField, ac_, fluid, *storagePtr, *block ) );
+
+ handling->fillWithDomain( FieldGhostLayers );
+
+ return handling;
+}
+
+void clearBoundaryHandling( BlockForest & forest, const BlockDataID & boundaryHandlingID ) {
+ for( auto blockIt = forest.begin(); blockIt != forest.end(); ++blockIt )
+ {
+ BoundaryHandling_T * boundaryHandling = blockIt->getData<BoundaryHandling_T>(boundaryHandlingID);
+ boundaryHandling->clear( FieldGhostLayers );
+ }
+}
+
+void recreateBoundaryHandling( BlockForest & forest, const BlockDataID & boundaryHandlingID ) {
+ for( auto blockIt = forest.begin(); blockIt != forest.end(); ++blockIt )
+ {
+ BoundaryHandling_T * boundaryHandling = blockIt->getData<BoundaryHandling_T>(boundaryHandlingID);
+ boundaryHandling->fillWithDomain( FieldGhostLayers );
+ }
+}
+
+void clearParticleField( BlockForest & forest, const BlockDataID & particleFieldID, const ParticleAccessor_T & accessor ) {
+ for (auto blockIt = forest.begin(); blockIt != forest.end(); ++blockIt) {
+ ParticleField_T *particleField = blockIt->getData<ParticleField_T>(particleFieldID);
+ particleField->setWithGhostLayer(accessor.getInvalidUid());
+ }
+}
+
+class SpherePropertyLogger {
+public:
+ SpherePropertyLogger(lbm_mesapd_coupling::SubCyclingManager* mesapdTimeStep, const shared_ptr<ParticleAccessor_T>& ac, walberla::id_t sphereUid,
+ const std::string& fileName, bool fileIO, bool writeHeader, real_t t_g, real_t u_g, real_t diameter) :
+ mesapdTimeStep_(mesapdTimeStep), ac_(ac), sphereUid_(sphereUid), fileName_(fileName), fileIO_(fileIO), t_g_(t_g),
+ u_g_(u_g), diameter_(diameter), position_(real_t(0)) {
+ if (fileIO_ && writeHeader) {
+ WALBERLA_ROOT_SECTION() {
+ std::ofstream file;
+ file.open(fileName_.c_str());
+ file << "#\t posZ\t uZ\t u\t rotX\t rotY\t rotZ\t hydFZ\t hydT\t tN\t posXN\t posYN\t posZN\t uXN\t uYN\t uZN\t uN\t wX\t wY\t wZ\n";
+ file.close();
+ }
+ }
+ }
+
+ void operator()() {
+ const uint_t timestep (mesapdTimeStep_->getCurrentTimeStep());
+ Vector3<real_t> transVel;
+ Vector3<real_t> hydForce;
+ Vector3<real_t> hydTorque;
+ Vector3<real_t> angVel;
+ Vector3<real_t> eulerRotation;
+
+ 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)) {
+ transVel = ac_->getLinearVelocity(idx);
+ hydForce = ac_->getHydrodynamicForce(idx);
+ hydTorque = ac_->getHydrodynamicTorque(idx);
+ angVel = ac_->getAngularVelocity(idx);
+ eulerRotation = ac_->getRotation(idx).getMatrix().getEulerAnglesXYZ();
+ }
+ }
+
+ WALBERLA_MPI_SECTION() {
+ 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);
+
+ mpi::allReduceInplace(hydForce[0], mpi::SUM);
+ mpi::allReduceInplace(hydForce[1], mpi::SUM);
+ mpi::allReduceInplace(hydForce[2], mpi::SUM);
+
+ mpi::allReduceInplace(hydTorque[0], mpi::SUM);
+ mpi::allReduceInplace(hydTorque[1], mpi::SUM);
+ mpi::allReduceInplace(hydTorque[2], mpi::SUM);
+
+ mpi::allReduceInplace(eulerRotation[0], mpi::SUM);
+ mpi::allReduceInplace(eulerRotation[1], mpi::SUM);
+ mpi::allReduceInplace(eulerRotation[2], mpi::SUM);
+ }
+
+ syncPosition();
+
+ if (fileIO_) {
+ writeToFile(timestep, position_, eulerRotation, transVel, angVel, hydForce, hydTorque);
+ }
+ }
+
+ Vector3<real_t> getPosition() const {
+ return position_;
+ }
+
+ void syncPosition() {
+ Vector3<real_t> pos(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)) {
+ pos = ac_->getPosition(idx);
+ }
+ }
+
+ WALBERLA_MPI_SECTION() {
+ mpi::allReduceInplace(pos[0], mpi::SUM);
+ mpi::allReduceInplace(pos[1], mpi::SUM);
+ mpi::allReduceInplace(pos[2], mpi::SUM);
+ }
+
+ position_ = pos;
+ }
+
+private:
+ void writeToFile(const uint_t timestep, const Vector3<real_t>& position, const Vector3<real_t>& eulerRotation,
+ const Vector3<real_t>& velocity, const Vector3<real_t>& angularVelocity,
+ const Vector3<real_t>& hydForce, const Vector3<real_t>& hydTorque) {
+ WALBERLA_ROOT_SECTION() {
+ std::ofstream file;
+ file.open(fileName_.c_str(), std::ofstream::app);
+
+ file << timestep << "\t"
+ << position[2] << "\t"
+ << velocity[2] << "\t"
+ << velocity.length() << "\t"
+ << eulerRotation[0] << "\t"
+ << eulerRotation[1] << "\t"
+ << eulerRotation[2] << "\t"
+ << hydForce[2] << "\t"
+ << hydTorque.length() << "\t"
+
+ << real_t(timestep) / t_g_ << "\t"
+ << position[0] / diameter_ << "\t"
+ << position[1] / diameter_ << "\t"
+ << position[2] / diameter_ << "\t"
+ << velocity[0] / u_g_ << "\t"
+ << velocity[1] / u_g_ << "\t"
+ << velocity[2] / u_g_ << "\t"
+ << velocity.length() / u_g_ << "\t"
+ << angularVelocity[0] * diameter_ / u_g_ << "\t"
+ << angularVelocity[1] * diameter_ / u_g_ << "\t"
+ << angularVelocity[2] * diameter_ / u_g_ << "\n";
+
+ file.close();
+ }
+ }
+
+ lbm_mesapd_coupling::SubCyclingManager* mesapdTimeStep_;
+ shared_ptr<ParticleAccessor_T> ac_;
+ const walberla::id_t sphereUid_;
+ std::string fileName_;
+ bool fileIO_;
+
+ real_t t_g_;
+ real_t u_g_;
+ real_t diameter_;
+
+ Vector3<real_t> position_;
+};
+
+void createCheckpointFiles(const std::string& checkPointFileName, std::vector<std::string>& oldCheckpointFiles,
+ uint_t lbmTimeStep, uint_t mesaTimeStep, const shared_ptr<BlockForest>& forest, const BlockDataID pdfFieldID,
+ const BlockDataID particleStorageID) {
+ WALBERLA_LOG_INFO_ON_ROOT("Writing checkpointing files in lbm time step " << lbmTimeStep << " (mesapd: " << mesaTimeStep << ").");
+
+ const std::string timeStepTag = std::to_string(lbmTimeStep) + "_" + std::to_string(mesaTimeStep);
+
+ const std::string lbmFileName = checkPointFileName + "_" + timeStepTag + "_lbm.txt";
+ const std::string mesaFileName = checkPointFileName + "_" + timeStepTag + "_mesa.txt";
+ const std::string forestFileName = checkPointFileName + "_" + timeStepTag + "_forest.txt";
+
+ forest->saveBlockData(lbmFileName, pdfFieldID);
+ forest->saveBlockData(mesaFileName, particleStorageID);
+ forest->saveToFile(forestFileName);
+
+ WALBERLA_ROOT_SECTION() {
+ if (!oldCheckpointFiles.empty()) {
+ for (const std::string& file: oldCheckpointFiles) {
+ filesystem::path path(file);
+ if (filesystem::exists(path)) {
+ filesystem::remove(path);
+ }
+ }
+ oldCheckpointFiles.clear();
+ WALBERLA_LOG_INFO("Deleted old checkpoint files.");
+ }
+
+ oldCheckpointFiles.push_back(lbmFileName);
+ oldCheckpointFiles.push_back(mesaFileName);
+ oldCheckpointFiles.push_back(forestFileName);
+ }
+}
+
+// main
+
+int main(int argc, char** argv) {
+ debug::enterTestMode();
+
+ mpi::Environment env(argc, argv);
+
+ logging::Logging::instance()->setStreamLogLevel(logging::Logging::INFO);
+ logging::Logging::instance()->setFileLogLevel(logging::Logging::INFO);
+
+ math::seedRandomGenerator( static_cast<unsigned int>(1337 * mpi::MPIManager::instance()->worldRank()) );
+
+ // general params
+ bool fileIO = true;
+ uint_t vtkIOFreq = 0;
+ uint_t fullVtkIOFreq = 0;
+ uint_t qCriterionVtkIOFreq = 0;
+ std::string vtkIOSelection = "sliced"; // "none", "full", "thicksliced" (half the sphere still in), "trackingsliced" (sphere is followed in y direction)
+ std::string baseFolder = "vtk_out_UnsettlingSphereDR";
+ auto vtkLowerVelocityLimit = real_t(0.003);
+ auto vtkLowerQCriterionLimit = real_t(1e-10);
+ auto propertyLoggingFreq = uint_t(1);
+
+ // numerical params
+ auto blockSize = uint_t(32);
+ auto XBlocks = uint_t(32);
+ auto YBlocks = uint_t(32);
+ auto ZBlocks = uint_t(64);
+ auto refinementCheckFrequency = uint_t(0);
+ auto initialZ = real_t(-1);
+ auto diameter = real_t(16);
+ bool offsetInitialPosition = true;
+ bool averageForceTorqueOverTwoTimeSteps = true;
+ auto numMESAPDSubCycles = uint_t(10);
+ auto u_ref = real_t(0.01);
+ auto densitySphere = real_t(0.1);
+ auto timesteps = uint_t(2000);
+ auto galileoNumber = real_t(300);
+ auto Lambda_Bulk = real_t(100); // if = 1, normal TRT, else not. Up to 100.
+ bool conserveMomentum = false;
+
+ // virtual mass
+ bool useVirtualMass = true;
+ auto C_v = real_t(0.5);
+ auto C_v_omega = real_t(0.5);
+
+ // dynamic refinement
+ auto numberOfLevels = uint_t(3);
+
+ bool useVorticityCriterion = false;
+ auto lowerFluidRefinementLimit = real_t(0);
+ auto upperFluidRefinementLimit = std::numeric_limits<real_t>::infinity();
+
+ bool useCurlCriterion = true;
+ auto upperCurlLimit = real_t(0.001);
+ auto curlCoarsenFactor = real_t(0.2);
+ auto curlLengthScaleWeight = real_t(2);
+
+ // checkpointing
+ bool writeCheckPointFile = false;
+ bool readFromCheckPointFile = false;
+ uint_t checkPointingFreq = 0;
+ uint_t checkPointLBMTimeStep = 0;
+ uint_t checkPointMESAPDTimeStep = 0;
+
+ for (int i = 1; i < argc; i++) {
+ if (std::strcmp(argv[i], "--shortRun") == 0) { fileIO = false; densitySphere = real_t(0.02); galileoNumber = 100; timesteps = 2; continue; }
+
+ if (std::strcmp(argv[i], "--noLogging") == 0) { fileIO = false; continue; }
+ if (std::strcmp(argv[i], "--vtkIOFreq") == 0) { vtkIOFreq = uint_c( std::stod( argv[++i] ) ); continue; }
+ if (std::strcmp(argv[i], "--fullVtkIOFreq") == 0) { fullVtkIOFreq = uint_c( std::stod( argv[++i] ) ); continue; }
+ if (std::strcmp(argv[i], "--qCriterionVtkIOFreq") == 0) { qCriterionVtkIOFreq = uint_c( std::stod( argv[++i] ) ); continue; }
+ if (std::strcmp(argv[i], "--vtkIOSelection") == 0) { vtkIOSelection = argv[++i]; continue; }
+ if (std::strcmp(argv[i], "--numMESAPDSubCycles") == 0) { numMESAPDSubCycles = uint_c( std::stod( argv[++i] ) ); continue; }
+ if (std::strcmp(argv[i], "--noForceAveraging") == 0) { averageForceTorqueOverTwoTimeSteps = false; continue; }
+ if (std::strcmp(argv[i], "--baseFolder") == 0) { baseFolder = argv[++i]; continue; }
+ if (std::strcmp(argv[i], "--vtkLowerVelocityLimit") == 0) { vtkLowerVelocityLimit = real_c(std::stod(argv[++i])); continue; }
+ if (std::strcmp(argv[i], "--vtkLowerQCriterionLimit") == 0) { vtkLowerQCriterionLimit = real_c(std::stod(argv[++i])); continue; }
+ if (std::strcmp(argv[i], "--propertyLoggingFreq") == 0) { propertyLoggingFreq = uint_c( std::stod( argv[++i] ) ); continue; }
+
+ if (std::strcmp(argv[i], "--dontOffsetInitialPos") == 0) { offsetInitialPosition = false; continue; }
+ if (std::strcmp(argv[i], "--sphereDensity") == 0) { densitySphere = real_c(std::stod(argv[++i])); continue; }
+ if (std::strcmp(argv[i], "--u_ref") == 0) { u_ref = real_c(std::stod(argv[++i])); continue; }
+ if (std::strcmp(argv[i], "--timesteps") == 0) { timesteps = uint_c(std::stod(argv[++i])); continue; }
+ if (std::strcmp(argv[i], "--blockSize") == 0) { blockSize = uint_c(std::stod(argv[++i])); continue; }
+ if (std::strcmp(argv[i], "--XBlocks") == 0) { XBlocks = uint_c(std::stod(argv[++i])); continue; }
+ if (std::strcmp(argv[i], "--YBlocks") == 0) { YBlocks = uint_c(std::stod(argv[++i])); continue; }
+ if (std::strcmp(argv[i], "--ZBlocks") == 0) { ZBlocks = uint_c(std::stod(argv[++i])); continue; }
+ if (std::strcmp(argv[i], "--refinementCheckFrequency") == 0) { refinementCheckFrequency = uint_c(std::stod(argv[++i])); continue; }
+ if (std::strcmp(argv[i], "--initialZ") == 0) { initialZ = real_c(std::stod(argv[++i])); continue; }
+ if (std::strcmp(argv[i], "--diameter") == 0) { diameter = real_c(std::stod(argv[++i])); continue; }
+ if (std::strcmp(argv[i], "--galileoNumber") == 0) { galileoNumber = real_c(std::stod(argv[++i])); continue; }
+ if (std::strcmp(argv[i], "--lambdaBulk") == 0) { Lambda_Bulk = real_c(std::stod(argv[++i])); continue; }
+
+ if (std::strcmp(argv[i], "--noVirtualMass") == 0) { useVirtualMass = false; continue; }
+ if (std::strcmp(argv[i], "--c_v") == 0) { C_v = real_c(std::stod(argv[++i])); continue; }
+ if (std::strcmp(argv[i], "--c_v_omega") == 0) { C_v_omega = real_c(std::stod(argv[++i])); continue; }
+
+ // dynamic refinement
+ if (std::strcmp(argv[i], "--numberOfLevels") == 0) { numberOfLevels = uint_c(std::stod(argv[++i])); continue; }
+ if (std::strcmp(argv[i], "--useVorticityCriterion") == 0) { useVorticityCriterion = true; useCurlCriterion = false; continue; }
+ if (std::strcmp(argv[i], "--lowerLimit") == 0) { lowerFluidRefinementLimit = real_c(std::stod(argv[++i])); continue; }
+ if (std::strcmp(argv[i], "--upperLimit") == 0) { upperFluidRefinementLimit = real_c(std::stod(argv[++i])); continue; }
+
+ if (std::strcmp(argv[i], "--upperCurlLimit") == 0) { upperCurlLimit = real_c(std::stod(argv[++i])); continue; }
+ if (std::strcmp(argv[i], "--curlCoarsenFactor") == 0) { curlCoarsenFactor = real_c(std::stod(argv[++i])); continue; }
+ if (std::strcmp(argv[i], "--curlLengthScaleWeight") == 0) { curlLengthScaleWeight = real_c(std::stod(argv[++i])); continue; }
+
+ if (std::strcmp(argv[i], "--writeCheckPointFile") == 0) { writeCheckPointFile = true; continue; }
+ if (std::strcmp(argv[i], "--readFromCheckPointFile") == 0) { readFromCheckPointFile = true; continue; }
+ if (std::strcmp(argv[i], "--checkPointingFreq") == 0) { checkPointingFreq = uint_c(std::stod( argv[++i]) ); continue; }
+ if (std::strcmp(argv[i], "--checkPointLBMTimeStep") == 0) { checkPointLBMTimeStep = uint_c(std::stod( argv[++i]) ); continue; }
+ if (std::strcmp(argv[i], "--checkPointMESAPDTimeStep") == 0) { checkPointMESAPDTimeStep = uint_c(std::stod( argv[++i]) ); continue; }
+
+ WALBERLA_ABORT("Unrecognized command line argument found: " << argv[i]);
+ }
+
+ WALBERLA_CHECK(!realIsEqual(u_ref, real_t(0)), "u_ref has to be non-zero.")
+ WALBERLA_CHECK(useCurlCriterion || useVorticityCriterion,
+ "Use curl or vorticity (legacy) criterion for refinement.");
+ WALBERLA_CHECK(!(useCurlCriterion && useVorticityCriterion),
+ "Using curl and vorticity criterion together makes no sense.");
+
+ // create base dir if it doesnt already exist
+ filesystem::path bpath(baseFolder);
+ if (!filesystem::exists(bpath)) {
+ filesystem::create_directory(bpath);
+ }
+ filesystem::path qpath(baseFolder + "/q_criterion_mesh");
+ if (filesystem::exists(qpath)) {
+ filesystem::remove_all(qpath);
+ }
+
+ //// numerical parameters
+
+ const auto densityFluid = real_t(1);
+
+ // calculate relaxation time from u_ref
+ real_t viscosity = u_ref * diameter / galileoNumber;
+ real_t omega = lbm::collision_model::omegaFromViscosity(viscosity);
+ real_t relaxationTime = real_t(1) / omega;
+
+ const auto d3 = diameter*diameter*diameter;
+ const auto G2 = galileoNumber*galileoNumber;
+ const auto v2 = viscosity*viscosity;
+ const auto gravitationalAcceleration = G2*v2 / (std::fabs(densitySphere-real_t(1))*d3);
+
+ const auto dx = real_t(1);
+ const auto overlap = real_t(4.5) * dx;
+
+ const Vector3<uint_t> domainSize( XBlocks * blockSize, YBlocks * blockSize, ZBlocks * blockSize );
+ if (realIsIdentical(initialZ, real_t(-1))) {
+ initialZ = real_c(blockSize) + real_t(0.5) * diameter;
+ }
+ Vector3<real_t> initialPosition(real_t(0.5) * real_c(domainSize[0]),
+ real_t(0.5) * real_c(domainSize[1]),
+ initialZ);
+ if (offsetInitialPosition) {
+ // offset sphere's initial position half a cell to introduce a slight disturbance.
+ initialPosition += Vector3<real_t>(real_t(0.5), real_t(0.5), real_t(0));
+ }
+
+ if (useVorticityCriterion && floatIsEqual(lowerFluidRefinementLimit, real_t(0)) && std::isinf(upperFluidRefinementLimit)) {
+ // use computed criterion instead of user input
+ lowerFluidRefinementLimit = real_t(0.05) * u_ref;
+ upperFluidRefinementLimit = real_t(0.1) * u_ref;
+ }
+
+ //const auto dx = real_t(1);
+
+ const auto u_g = std::sqrt(std::fabs(densitySphere - real_t(1)) * gravitationalAcceleration * diameter);
+ const auto t_g = std::sqrt(diameter / (std::fabs(densitySphere - real_t(1)) * gravitationalAcceleration));
+
+ const bool rising = densitySphere < real_t(1);
+
+ WALBERLA_LOG_INFO_ON_ROOT("Setup (in simulation, i.e. lattice, units):");
+ WALBERLA_LOG_INFO_ON_ROOT(" - number of processes = " << mpi::MPIManager::instance()->numProcesses());
+ WALBERLA_LOG_INFO_ON_ROOT(" - domain size = " << domainSize);
+ WALBERLA_LOG_INFO_ON_ROOT(" - domain size / D = " << (Vector3<real_t>(domainSize) / diameter));
+ WALBERLA_LOG_INFO_ON_ROOT(" - sphere: diameter = " << diameter << ", density = " << densitySphere << ", initial position: " << initialPosition);
+ WALBERLA_LOG_INFO_ON_ROOT(" - fluid: density = " << densityFluid << ", relaxation time (tau) = " << relaxationTime << ", kin. visc = " << viscosity );
+ WALBERLA_LOG_INFO_ON_ROOT(" - lambda bulk = " << Lambda_Bulk );
+ WALBERLA_LOG_INFO_ON_ROOT(" - galileo number = " << galileoNumber );
+ WALBERLA_LOG_INFO_ON_ROOT(" - gravitational velocity (u_g) = " << u_g);
+ WALBERLA_LOG_INFO_ON_ROOT(" - gravitational acceleration = " << gravitationalAcceleration );
+#ifdef BUILD_WITH_MRT
+ WALBERLA_LOG_INFO_ON_ROOT(" - lattice model: MRT");
+#else
+ WALBERLA_LOG_INFO_ON_ROOT(" - lattice model: TRT");
+#endif
+ WALBERLA_LOG_INFO_ON_ROOT(" - integrator = Velocity Verlet" );
+ WALBERLA_LOG_INFO_ON_ROOT(" - use virtual mass = " << useVirtualMass );
+ //WALBERLA_LOG_INFO_ON_ROOT(" - virtual mass acceleration estimation = " << vmAccelerationOption );
+ WALBERLA_LOG_INFO_ON_ROOT(" - C_v = " << C_v );
+ WALBERLA_LOG_INFO_ON_ROOT(" - C_v_omega = " << C_v_omega );
+ WALBERLA_LOG_INFO_ON_ROOT(" - timesteps = " << timesteps);
+ WALBERLA_LOG_INFO_ON_ROOT(" - averageForceTorqueOverTwoTimeSteps = " << averageForceTorqueOverTwoTimeSteps);
+ WALBERLA_LOG_INFO_ON_ROOT(" - conserve momentum = " << conserveMomentum);
+ if (useCurlCriterion) {
+ WALBERLA_LOG_INFO_ON_ROOT(" - using curl criterion with upper limit = " << upperCurlLimit << ", coarsen factor = " << curlCoarsenFactor << " and length scale weight = " << curlLengthScaleWeight);
+ }
+ if (useVorticityCriterion) {
+ WALBERLA_LOG_INFO_ON_ROOT(" - using vorticity criterion with lower limit = " << lowerFluidRefinementLimit << " and upper limit = " << upperFluidRefinementLimit );
+ }
+ if (vtkIOFreq > 0) {
+ WALBERLA_LOG_INFO_ON_ROOT(" - writing vtk files to folder \"" << baseFolder << "\" with frequency " << vtkIOFreq);
+ }
+
+ if(refinementCheckFrequency == 0 && numberOfLevels != 1){
+ // determine check frequency automatically based on maximum admissible velocity and block sizes
+ auto uMax = real_t(0.02); // this should never ever be reached
+ const uint_t finestLevel = numberOfLevels - uint_t(1);
+ const uint_t levelScalingFactor = ( uint_t(1) << finestLevel );
+ const uint_t lbmTimeStepsPerTimeLoopIteration = levelScalingFactor;
+ refinementCheckFrequency = uint_c((overlap + (real_c(blockSize) - real_t(2) * real_t(FieldGhostLayers)) * dx) / uMax) / lbmTimeStepsPerTimeLoopIteration;
+ }
+ WALBERLA_LOG_INFO_ON_ROOT(" - refinement / load balancing check frequency (coarse time steps): " << refinementCheckFrequency);
+
+ //// checkpointing setup
+
+ WALBERLA_LOG_INFO_ON_ROOT(" - read from checkpoint file = " << readFromCheckPointFile);
+ WALBERLA_LOG_INFO_ON_ROOT(" - write checkpoint file = " << writeCheckPointFile);
+
+ std::string checkPointFileName = "checkpointing";
+ if (readFromCheckPointFile || writeCheckPointFile) {
+ WALBERLA_LOG_INFO_ON_ROOT(" - checkPointingFileName = " << checkPointFileName);
+ }
+
+ std::string readCheckPointFileName = checkPointFileName + "_" + std::to_string(checkPointLBMTimeStep) + "_" +
+ std::to_string(checkPointMESAPDTimeStep);
+
+ if (readFromCheckPointFile && writeCheckPointFile) {
+ // decide which option to choose
+ if (filesystem::exists(readCheckPointFileName+"_lbm.txt")) {
+ WALBERLA_LOG_INFO_ON_ROOT("Checkpoint file already exists! Will skip writing check point file and start from this check point!");
+ writeCheckPointFile = false;
+ } else {
+ WALBERLA_LOG_INFO_ON_ROOT("Checkpoint file does not exists yet! Will skip reading check point file and just regularly start the simulation from the beginning!");
+ readFromCheckPointFile = false;
+ }
+ }
+
+ // write parameters to file
+ if (fileIO) {
+ WALBERLA_ROOT_SECTION() {
+ std::string paramsFileName(baseFolder + "/Params.txt");
+ std::ofstream file;
+ file.open(paramsFileName.c_str());
+ file
+ << "D\t "
+ << "Ga\t "
+ << "LBM\t "
+ << "pi_p\t "
+ << "tau\t "
+ << "visc\t "
+ << "u_g\t "
+ << "t_g\t "
+ << "acc_g\t "
+ << "vm\t "
+ << "C_v\t "
+ << "C_vw\t "
+ << "fa\t "
+ << "curlCrit\t "
+ << "curlMax\t "
+ << "curlCoaF\t "
+ << "curlLenS\t "
+ << "refChk\t "
+ << "numLvl\t "
+ << "sX\t "
+ << "sY\t "
+ << "sZ\t "
+ << "nProc\t "
+ << "\n";
+ file
+ << diameter << "\t "
+ << galileoNumber << "\t "
+ #ifdef BUILD_WITH_MRT
+ << "MRT\t "
+ #else
+ << "TRT\t "
+ #endif
+ << densitySphere << "\t "
+ << relaxationTime << "\t "
+ << viscosity << "\t "
+ << u_g << "\t "
+ << t_g << "\t "
+ << gravitationalAcceleration << "\t "
+ << (useVirtualMass ? "lt" : "no") << "\t "
+ << C_v << "\t "
+ << C_v_omega << "\t "
+ << averageForceTorqueOverTwoTimeSteps << "\t "
+ << useCurlCriterion << "\t "
+ << upperCurlLimit << "\t "
+ << curlCoarsenFactor << "\t "
+ << curlLengthScaleWeight << "\t "
+ << refinementCheckFrequency << "\t "
+ << numberOfLevels << "\t "
+ << domainSize[0] << "\t "
+ << domainSize[1] << "\t "
+ << domainSize[2] << "\t "
+ << mpi::MPIManager::instance()->numProcesses() << "\t "
+ << "\n";
+ file.close();
+ }
+ }
+
+ //
+
+ uint_t currentLBMTimeStep = 0;
+ uint_t currentMESAPDTimeStep = 0;
+
+ if (readFromCheckPointFile) {
+ currentLBMTimeStep = checkPointLBMTimeStep;
+ currentMESAPDTimeStep = checkPointMESAPDTimeStep;
+ }
+
+ //// block structure setup
+
+ const uint_t finestLevel = numberOfLevels - uint_t(1);
+
+ Vector3<uint_t> blockSizeInCells(blockSize);
+
+ AABB simulationDomain( real_t(0), real_t(0), real_t(0), real_c(domainSize[0]), real_c(domainSize[1]), real_c(domainSize[2]) );
+ Vector3<bool> periodicity(true, true, true);
+ auto blocks = createBlockStructure(simulationDomain, blockSizeInCells, numberOfLevels, diameter, initialPosition,
+ periodicity, readFromCheckPointFile, readCheckPointFileName+"_forest.txt");
+ auto forest = blocks->getBlockForestPointer();
+
+ if (vtkIOFreq > 0) {
+ vtk::writeDomainDecomposition(blocks, "initial_domain_decomposition", baseFolder);
+ }
+
+ //// dynamic refinement I
+
+ forest->recalculateBlockLevelsInRefresh(true);
+ forest->alwaysRebalanceInRefresh(false);
+ forest->reevaluateMinTargetLevelsAfterForcedRefinement(false);
+ forest->allowRefreshChangingDepth(false);
+
+ //// mesapd coupling
+
+ auto mesapdDomain = 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>();
+
+ auto accessor = walberla::make_shared<ParticleAccessor_T>(ps, ss);
+
+ BlockDataID particleStorageID;
+ if (readFromCheckPointFile) {
+ WALBERLA_LOG_INFO_ON_ROOT("Initializing particles from checkpointing file.");
+ particleStorageID = forest->loadBlockData( readCheckPointFileName + "_mesa.txt", mesa_pd::domain::createBlockForestDataHandling(ps), "Particle Storage" );
+ } else {
+ particleStorageID = forest->addBlockData(mesa_pd::domain::createBlockForestDataHandling(ps), "Particle Storage");
+ }
+
+
+ // bounding planes can be added here
+
+ auto sphereShape = ss->create<mesa_pd::data::Sphere>(diameter * real_t(0.5));
+ ss->shapes[sphereShape]->updateMassAndInertia(densitySphere);
+
+ walberla::id_t sphereUid = 0;
+ if (readFromCheckPointFile) {
+ for (auto pIt = ps->begin(); pIt != ps->end(); ++pIt) {
+ // find sphere in loaded data structure and store uid for later reference
+ if (pIt->getShapeID() == sphereShape) {
+ sphereUid = pIt->getUid();
+ pIt->getNeighborStateRef().clear();
+ pIt->getGhostOwnersRef().clear();
+ }
+ }
+ } else {
+ if (mesapdDomain->isContainedInProcessSubdomain(uint_c(mpi::MPIManager::instance()->rank()), initialPosition)) {
+ mesa_pd::data::Particle &&p = *ps->create();
+ p.setPosition(initialPosition);
+ p.setInteractionRadius(diameter * real_t(0.5));
+ p.setOwner(mpi::MPIManager::instance()->rank());
+ p.setShapeID(sphereShape);
+ sphereUid = p.getUid();
+ }
+ }
+ mpi::allReduceInplace(sphereUid, mpi::SUM);
+
+ if(sphereUid == 0) {
+ WALBERLA_ABORT("No sphere present - aborting!");
+ }
+
+ auto ghostOwnerSyncIterations = uint_t(real_t(0.5) * diameter/real_t(blockSize) + real_t(1));
+ WALBERLA_LOG_INFO_ON_ROOT(" - Number of req'd ghost owner sync iterations: " << ghostOwnerSyncIterations << ", thus using " << (ghostOwnerSyncIterations == 1 ? "SNN" : "SGO") << " for synchronization.");
+
+ bool useSyncNextNeighbors = (ghostOwnerSyncIterations == 1);
+
+ std::function<void(void)> syncCall = [ps, mesapdDomain, useSyncNextNeighbors, overlap](){
+ if (useSyncNextNeighbors) {
+ mesa_pd::mpi::SyncNextNeighbors SNN; // alt. SyncGhostOwners
+ SNN(*ps, *mesapdDomain, overlap);
+ } else {
+ // need to use syncghostowners since the particle is bigger than a block and thus, syncnextneighbors doesn't suffice.
+ mesa_pd::mpi::SyncGhostOwners SGO;
+ SGO(*ps, *mesapdDomain, overlap);
+ }
+ };
+
+ std::function<void(void)> initializationSyncCall = [&syncCall, ghostOwnerSyncIterations](){
+ for (uint_t j = 0; j < ghostOwnerSyncIterations; ++j) {
+ syncCall();
+ }
+ };
+
+ //initial sync
+ initializationSyncCall();
+
+ //// add data to blocks
+
+ // create the lattice model
+#ifdef BUILD_WITH_MRT
+ const real_t Lambda_Magic = real_t(3) / real_t(16);
+ const real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega );
+ const real_t lambda_d = lbm::collision_model::TRT::lambda_d( omega, Lambda_Magic );
+ const real_t omegaBulk = real_t(1) / (Lambda_Bulk * ( real_t(1) / omega - real_t(1)/ real_t(2) ) + real_t(1)/ real_t(2) );
+
+ LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::D3Q19MRT( omegaBulk, lambda_e, lambda_d, lambda_e, lambda_e, lambda_d, finestLevel ));
+#else
+ LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( omega, lbm::collision_model::TRT::threeSixteenth, finestLevel ) );
+#endif
+
+ // add PDF field
+ BlockDataID pdfFieldID;
+ if (readFromCheckPointFile) {
+ // add PDF field
+ WALBERLA_LOG_INFO_ON_ROOT( "Initializing PDF Field from checkpointing file!" );
+ shared_ptr< lbm::internal::PdfFieldHandling< LatticeModel_T > > dataHandling =
+ make_shared< lbm::internal::PdfFieldHandling< LatticeModel_T > >(blocks, latticeModel, false,
+ Vector3<real_t>(real_t(0)), real_t(1),
+ FieldGhostLayers, field::zyxf );
+
+ pdfFieldID = blocks->loadBlockData( readCheckPointFileName+"_lbm.txt", dataHandling, "pdf field" );
+
+ } else {
+ // add PDF field
+ pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >(blocks, "pdf field (zyxf)", latticeModel,
+ Vector3<real_t>(real_t(0)), real_t(1),
+ FieldGhostLayers, field::zyxf);
+ }
+
+ // add flag field
+ BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>( blocks, "flag field", FieldGhostLayers );
+
+ using FluidFilter_T = typename field::FlagFieldEvaluationFilter<FlagField_T>;
+ FluidFilter_T fluidFlagFieldEvaluationFilter(flagFieldID, Fluid_Flag);
+
+ // add particle field
+ BlockDataID particleFieldID = field::addToStorage<ParticleField_T>(blocks, "particle field",
+ accessor->getInvalidUid(),
+ field::zyxf,
+ FieldGhostLayers);
+
+ // add velocity field and utility
+ BlockDataID velocityFieldID = field::addToStorage<VelocityField_T>( blocks, "velocity field", Vector3<real_t>(real_t(0)), field::zyxf, uint_t(2) );
+
+ typedef lbm::VelocityFieldWriter< PdfField_T, VelocityField_T > VelocityFieldWriter_T;
+ BlockSweepWrapper< VelocityFieldWriter_T > velocityFieldWriter( blocks, VelocityFieldWriter_T( pdfFieldID, velocityFieldID ) );
+
+ shared_ptr<blockforest::communication::NonUniformBufferedScheme<stencil::D3Q27> > velocityCommunicationScheme = make_shared<blockforest::communication::NonUniformBufferedScheme<stencil::D3Q27> >( blocks );
+ velocityCommunicationScheme->addPackInfo( make_shared< field::refinement::PackInfo<VelocityField_T, stencil::D3Q27> >( velocityFieldID ) );
+
+ // add q criterion field (only needed for mesh output)
+ BlockDataID qCriterionFieldID = field::addToStorage<QCriterionField_T>(blocks, "q criterion field", real_t(0), field::zyxf, uint_t(1));
+
+ typedef lbm::QCriterionFieldWriter<VelocityField_T, QCriterionField_T, FluidFilter_T> QCriterionFieldWriter_T;
+ BlockSweepWrapper<QCriterionFieldWriter_T> qCriterionFieldWriter(blocks, QCriterionFieldWriter_T(blocks, velocityFieldID,
+ qCriterionFieldID, fluidFlagFieldEvaluationFilter));
+
+ // add boundary handling
+ BlockDataID boundaryHandlingID = blocks->addBlockData( make_shared< MyBoundaryHandling >( blocks, flagFieldID, pdfFieldID, particleFieldID, accessor ),
+ "boundary handling" );
+
+ // map particles into simulation
+ lbm_mesapd_coupling::RegularParticlesSelector sphereSelector;
+ lbm_mesapd_coupling::ParticleMappingKernel<BoundaryHandling_T> particleMappingKernel(blocks, boundaryHandlingID);
+ lbm_mesapd_coupling::MovingParticleMappingKernel<BoundaryHandling_T> movingParticleMappingKernel(blocks,
+ boundaryHandlingID,
+ particleFieldID);
+
+ std::function<void(void)> mappingCall = [ps, sphereSelector, accessor, &movingParticleMappingKernel](){
+ ps->forEachParticle(false, sphereSelector, *accessor, movingParticleMappingKernel, *accessor, MO_Flag);
+ };
+
+ mappingCall();
+
+ syncCall();
+
+ //// dynamic refinement II
+
+ blockforest::CombinedMinTargetLevelDeterminationFunctions minTargetLevelDeterminationFunctions;
+
+ auto couplingInfoCollection = walberla::make_shared<lbm_mesapd_coupling::amr::InfoCollection>();
+ lbm_mesapd_coupling::amr::ParticlePresenceLevelDetermination particlePresenceRefinement(couplingInfoCollection, finestLevel);
+
+ minTargetLevelDeterminationFunctions.add(particlePresenceRefinement);
+
+ if (useVorticityCriterion) {
+ // add refinement criterion based on vorticity magnitude
+ // "legacy"
+ lbm::refinement::VorticityBasedLevelDetermination<field::FlagFieldEvaluationFilter<FlagField_T>> vorticityRefinement(
+ velocityFieldID, fluidFlagFieldEvaluationFilter, upperFluidRefinementLimit, lowerFluidRefinementLimit, finestLevel);
+
+ minTargetLevelDeterminationFunctions.add(vorticityRefinement);
+ }
+
+ if (useCurlCriterion) {
+ // add refinement criterion based on curl (=vorticity) magnitude
+ // -> computes the same as vorticity, but has other params and uses correct length scale weighting
+ const real_t lowerCurlLimit = curlCoarsenFactor*upperCurlLimit;
+
+ lbm::refinement::CurlBasedLevelDetermination<field::FlagFieldEvaluationFilter<FlagField_T>> curlRefinement(
+ velocityFieldID, *blocks, fluidFlagFieldEvaluationFilter, finestLevel, upperCurlLimit, lowerCurlLimit,
+ curlLengthScaleWeight);
+
+ minTargetLevelDeterminationFunctions.add(curlRefinement);
+ }
+
+ forest->setRefreshMinTargetLevelDeterminationFunction(minTargetLevelDeterminationFunctions);
+
+ bool curveHilbert = false;
+ bool curveAllGather = true;
+
+ forest->setRefreshPhantomBlockMigrationPreparationFunction( blockforest::DynamicCurveBalance< blockforest::NoPhantomData >( curveHilbert, curveAllGather ) );
+
+
+ //// time loop
+
+ // setup of the LBM communication for synchronizing the pdf field between neighboring blocks
+ std::function< void () > commFunction;
+ blockforest::communication::UniformBufferedScheme< Stencil_T > scheme( blocks );
+ scheme.addPackInfo( make_shared< lbm::PdfFieldPackInfo< LatticeModel_T > >( pdfFieldID ) );
+ commFunction = scheme;
+
+ // create the timeloop
+ SweepTimeloop timeloop( blocks->getBlockStorage(), timesteps );
+ timeloop.setCurrentTimeStep(currentLBMTimeStep);
+
+ shared_ptr<WcTimingPool> timeloopTiming = make_shared<WcTimingPool>();
+
+ auto sweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag );
+ auto refinementTimestep = lbm::refinement::makeTimeStep< LatticeModel_T, BoundaryHandling_T >( blocks, sweep, pdfFieldID, boundaryHandlingID );
+
+ //// MESAPD Time Step
+
+ const real_t sphereVolume = math::pi / real_t(6) * diameter * diameter * diameter;
+ Vector3<real_t> gravitationalForce(real_t(0), real_t(0),
+ -(densitySphere - densityFluid) * gravitationalAcceleration * sphereVolume);
+
+ lbm_mesapd_coupling::SubCyclingManager mesapdTimeStep(numMESAPDSubCycles, timeloopTiming);
+ mesapdTimeStep.setCurrentTimeStep(currentMESAPDTimeStep);
+ bool useOpenMP = false;
+
+ // before subcycling
+ mesa_pd::mpi::ReduceProperty reduceProperty;
+ lbm_mesapd_coupling::InitializeHydrodynamicForceTorqueForAveragingKernel initializeHydrodynamicForceTorqueAveraging;
+ lbm_mesapd_coupling::AverageHydrodynamicForceTorqueKernel averageHydrodynamicForceTorque;
+
+ mesapdTimeStep.addFuncBeforeSubCycles([ps, accessor, useOpenMP, &reduceProperty, averageForceTorqueOverTwoTimeSteps,
+ &initializeHydrodynamicForceTorqueAveraging, &averageHydrodynamicForceTorque,
+ &mesapdTimeStep](){
+ reduceProperty.operator()<mesa_pd::HydrodynamicForceTorqueNotification>(*ps);
+ if (averageForceTorqueOverTwoTimeSteps) {
+ if (mesapdTimeStep.getCurrentTimeStep() == 0) {
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, initializeHydrodynamicForceTorqueAveraging, *accessor);
+ } else {
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, averageHydrodynamicForceTorque, *accessor);
+ }
+ }
+ }, "RPD Before Subcycling");
+
+ // during subcycling
+ using VirtualMass_ParticleAccessor_T = lbm_mesapd_coupling::ParticleAccessorWithShapeVirtualMassWrapper<ParticleAccessor_T>;
+ auto virtualMassAccessor = walberla::make_shared<VirtualMass_ParticleAccessor_T>(ps, ss);
+
+ mesa_pd::kernel::VelocityVerletPreForceUpdate vvIntegratorPreForce(real_t(1) / real_t(numMESAPDSubCycles));
+ mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(real_t(1) / real_t(numMESAPDSubCycles));
+ lbm_mesapd_coupling::AddForceOnParticlesKernel addGravitationalForce(gravitationalForce);
+ lbm_mesapd_coupling::AddHydrodynamicInteractionKernel addHydrodynamicInteraction;
+
+ mesapdTimeStep.addFuncDuringSubCycles([ps, accessor, virtualMassAccessor, useOpenMP, syncCall,
+ useVirtualMass,
+ vvIntegratorPreForce, addHydrodynamicInteraction, addGravitationalForce,
+ vvIntegratorPostForce, mesapdDomain,
+ C_v, C_v_omega, densityFluid](){
+ if (useVirtualMass) {
+ lbm_mesapd_coupling::InitializeVirtualMassKernel virtualMass;
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, virtualMass, *accessor, C_v, C_v_omega, densityFluid);
+
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *virtualMassAccessor, vvIntegratorPreForce, *virtualMassAccessor);
+ } else {
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, vvIntegratorPreForce, *accessor);
+ }
+ syncCall();
+
+ // Hydrodynamic force and torque got reduced onto particle local process before subcycling,
+ // thus only add them there.
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addHydrodynamicInteraction, *accessor);
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addGravitationalForce, *accessor);
+
+ if (useVirtualMass) {
+ lbm_mesapd_coupling::InitializeVirtualMassKernel virtualMass;
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, virtualMass, *accessor, C_v, C_v_omega, densityFluid);
+
+ lbm_mesapd_coupling::AddVirtualForceTorqueKernel addVirtualForceAndTorque(ps);
+
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addVirtualForceAndTorque, *accessor);
+ }
+
+ if (useVirtualMass) {
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *virtualMassAccessor, vvIntegratorPostForce, *virtualMassAccessor);
+ } else {
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, vvIntegratorPostForce, *accessor);
+ }
+
+ syncCall();
+ }, "RPD Subcycling");
+
+ // after subcycling
+ // evaluation
+ std::string loggingFileName(baseFolder + "/LoggingUnsettlingSphereDynamicRefinement");
+ if (useVirtualMass) {
+ loggingFileName += "_vm";
+ }
+ loggingFileName += ".txt";
+ if (fileIO) {
+ WALBERLA_LOG_INFO_ON_ROOT(" - writing logging output to file \"" << loggingFileName << "\"");
+ }
+ SpherePropertyLogger logger(&mesapdTimeStep, accessor, sphereUid, loggingFileName, fileIO,
+ !readFromCheckPointFile, t_g, u_g, diameter);
+
+ auto periodicLogger = [&mesapdTimeStep, &logger, propertyLoggingFreq](){
+ if (propertyLoggingFreq != 0 && mesapdTimeStep.getCurrentTimeStep() % propertyLoggingFreq == 0) {
+ logger();
+ }
+ };
+
+ mesapdTimeStep.addFuncAfterSubCycles(periodicLogger, "Sphere property logger" );
+
+ lbm_mesapd_coupling::ResetHydrodynamicForceTorqueKernel resetHydrodynamicForceTorque;
+ mesapdTimeStep.addFuncAfterSubCycles([ps, accessor, useOpenMP, resetHydrodynamicForceTorque](){
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, resetHydrodynamicForceTorque, *accessor);
+ }, "RPD After Subcycling");
+
+ // mesapdTimeStep seems to get copied or something by the functorwrapper
+ auto wrapper = lbm::refinement::FunctorWrapper([&mesapdTimeStep](){
+ mesapdTimeStep();
+ });
+ refinementTimestep->addPostStreamVoidFunction(wrapper, "rpd Time Step", finestLevel);
+
+
+ //// Particle Mapping
+
+ // sweep for updating the particle mapping into the LBM simulation
+ auto movingParticleMapping = lbm_mesapd_coupling::makeMovingParticleMapping<PdfField_T, BoundaryHandling_T>(blocks,
+ pdfFieldID,
+ boundaryHandlingID,
+ particleFieldID,
+ accessor,
+ MO_Flag,
+ FormerMO_Flag,
+ sphereSelector,
+ conserveMomentum);
+ refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper(movingParticleMapping, blocks),
+ "Particle Mapping", finestLevel);
+
+ //// PDF Restore
+
+ blockforest::communication::NonUniformBufferedScheme<Stencil_T> fullPDFCommunicationScheme(blocks);
+ auto uniformPackInfo = make_shared<field::communication::PackInfo<PdfField_T>>(pdfFieldID);
+ auto packInfoAdapter = make_shared<blockforest::communication::UniformToNonUniformPackInfoAdapter>(uniformPackInfo);
+ // using the pack info adapter is required, else the ghost layers are not correctly synced, resulting in invalid reconstructed cells (with NaN as value).
+ fullPDFCommunicationScheme.addPackInfo(packInfoAdapter);
+ refinementTimestep->addPostStreamVoidFunction(lbm::refinement::FunctorWrapper(fullPDFCommunicationScheme),
+ "PDF Communication", finestLevel);
+
+ // add sweep for restoring PDFs in cells previously occupied by particles
+ 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);
+ //auto equilibriumRecon = lbm_mesapd_coupling::makeEquilibriumReconstructor<BoundaryHandling_T>(
+ // blocks, boundaryHandlingID);
+ auto reconstructionManager = lbm_mesapd_coupling::makePdfReconstructionManager<PdfField_T, BoundaryHandling_T>(
+ blocks,
+ pdfFieldID,
+ boundaryHandlingID,
+ particleFieldID,
+ accessor,
+ FormerMO_Flag,
+ Fluid_Flag,
+ extrapolationSphereNormalReconstructor,
+ conserveMomentum);
+ refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper(*reconstructionManager, blocks),
+ "PDF Restore", finestLevel);
+
+ // add LBM sweep with refinement
+ timeloop.addFuncBeforeTimeStep( makeSharedFunctor( refinementTimestep ), "LBM refinement time step" );
+
+
+ //// VTK output
+
+ auto vtkInitialExecutionCount = uint_t(0);
+ std::string vtkFileTag = "";
+ if (readFromCheckPointFile) {
+ vtkInitialExecutionCount = checkPointLBMTimeStep;
+ vtkFileTag = "_" + std::to_string(vtkInitialExecutionCount);
+ }
+
+ // sphere(s)
+ auto particleVtkOutput = make_shared<mesa_pd::vtk::ParticleVtkOutput>(ps);
+ particleVtkOutput->addOutput<mesa_pd::data::SelectParticleOwner>("owner");
+ particleVtkOutput->addOutput<mesa_pd::data::SelectParticleLinearVelocity>("velocity");
+ particleVtkOutput->addOutput<mesa_pd::data::SelectParticleInteractionRadius>("interactionRadius");
+ auto particleVtkWriter = vtk::createVTKOutput_PointData(particleVtkOutput, "Particles"+vtkFileTag, vtkIOFreq,
+ baseFolder, "simulation_step", false, true, true, true,
+ vtkInitialExecutionCount);
+
+ // flag field left out
+
+ // pdf field
+ auto pdfFieldVTK = vtk::createVTKOutput_BlockData(blocks, "fluid_field"+vtkFileTag, vtkIOFreq, 0,
+ false, baseFolder, "simulation_step", false, true, true, true,
+ vtkInitialExecutionCount);
+
+ vtk::ChainedFilter combinedFilter;
+ if (vtkIOSelection == "sliced") {
+ AABB sliceAABB( real_t(0), real_c(domainSize[1])*real_t(0.5), real_t(0),
+ real_c(domainSize[0]), real_c(domainSize[1])*real_t(0.5)+real_t(1), real_c(domainSize[2]) );
+ vtk::AABBCellFilter aabbSliceFilter( sliceAABB );
+ combinedFilter.addFilter( aabbSliceFilter );
+ } else if (vtkIOSelection == "thicksliced") {
+ AABB sliceAABB( real_t(0), real_c(domainSize[1])*real_t(0.5)-real_t(diameter)*real_t(0.5)-real_t(1), real_t(0),
+ real_c(domainSize[0]), real_c(domainSize[1])*real_t(0.5), real_c(domainSize[2]) );
+ vtk::AABBCellFilter aabbSliceFilter( sliceAABB );
+ combinedFilter.addFilter( aabbSliceFilter );
+ } else if (vtkIOSelection == "trackingsliced") {
+ auto filter = [logger, domainSize](CellSet& filteredCells, const IBlock& block, const StructuredBlockStorage& storage, const uint_t ghostLayers = uint_t(0)){
+ auto yPosition = logger.getPosition()[1];
+ real_t minY = std::floor(yPosition);
+ real_t maxY = std::floor(yPosition)+real_t(1);
+ AABB sliceAABB( real_t(0), minY, real_t(0),
+ real_c(domainSize[0]), maxY, real_c(domainSize[2]) );
+ vtk::AABBCellFilter aabbSliceFilter( sliceAABB );
+ aabbSliceFilter(filteredCells, block, storage, ghostLayers);
+ };
+ combinedFilter.addFilter( filter );
+ logger.syncPosition(); // synchronize the position once to correctly supply the vtk output before the first timestep
+ } else if (vtkIOSelection == "velocity") {
+ field::VelocityCellFilter<PdfField_T, field::FlagFieldEvaluationFilter<FlagField_T>> filter(pdfFieldID,
+ fluidFlagFieldEvaluationFilter, vtkLowerVelocityLimit);
+ combinedFilter.addFilter(filter);
+ } else if (vtkIOSelection != "none" && vtkIOSelection != "full") {
+ WALBERLA_ABORT("Invalid vtk output selection.");
+ }
+
+ field::FlagFieldCellFilter< FlagField_T > fluidFilter(flagFieldID);
+ fluidFilter.addFlag(Fluid_Flag);
+ combinedFilter.addFilter(fluidFilter);
+
+ pdfFieldVTK->addCellInclusionFilter(combinedFilter);
+
+ pdfFieldVTK->addBeforeFunction(velocityFieldWriter);
+ pdfFieldVTK->addBeforeFunction(*velocityCommunicationScheme);
+
+ pdfFieldVTK->addCellDataWriter(make_shared<lbm::VelocityVTKWriter<LatticeModel_T, float>>(pdfFieldID, "VelocityFromPDF"));
+ pdfFieldVTK->addCellDataWriter(make_shared<lbm::DensityVTKWriter<LatticeModel_T, float>>(pdfFieldID, "DensityFromPDF"));
+ pdfFieldVTK->addCellDataWriter(make_shared<lbm::QCriterionVTKWriter<VelocityField_T, FluidFilter_T, float>>(blocks,
+ fluidFlagFieldEvaluationFilter, velocityFieldID, "QCriterionFromPDF"));
+ pdfFieldVTK->addCellDataWriter(make_shared<lbm::CurlMagnitudeVTKWriter<VelocityField_T, FluidFilter_T, float>>(blocks,
+ fluidFlagFieldEvaluationFilter, velocityFieldID, "CurlSensor", real_t(2)));
+
+ // full size vtk
+ auto fullPdfFieldVTK = vtk::createVTKOutput_BlockData(blocks, "fluid_field_full"+vtkFileTag, fullVtkIOFreq,
+ 0, false, baseFolder, "simulation_step", false, true, true, true,
+ vtkInitialExecutionCount);
+ fullPdfFieldVTK->addBeforeFunction(velocityFieldWriter);
+ fullPdfFieldVTK->addBeforeFunction(*velocityCommunicationScheme);
+ fullPdfFieldVTK->addCellInclusionFilter(fluidFilter);
+ fullPdfFieldVTK->addCellDataWriter(make_shared<lbm::VelocityVTKWriter<LatticeModel_T, float>>(pdfFieldID, "VelocityFromPDF"));
+ fullPdfFieldVTK->addCellDataWriter(make_shared<lbm::QCriterionVTKWriter<VelocityField_T , FluidFilter_T, float>>(blocks,
+ fluidFlagFieldEvaluationFilter, velocityFieldID, "QCriterionFromPDF"));
+ fullPdfFieldVTK->addCellDataWriter(make_shared<lbm::CurlMagnitudeVTKWriter<VelocityField_T, FluidFilter_T, float>>(blocks,
+ fluidFlagFieldEvaluationFilter, velocityFieldID, "CurlSensor", curlLengthScaleWeight));
+
+ // q value vtk
+ auto qCriterionVTK = vtk::createVTKOutput_BlockData(blocks, "q_value"+vtkFileTag, qCriterionVtkIOFreq, 0,
+ false, baseFolder, "simulation_step", false, true, true, true,
+ vtkInitialExecutionCount);
+
+ field::QCriterionCellFilter<VelocityField_T, field::FlagFieldEvaluationFilter<FlagField_T>> qFilter(velocityFieldID,
+ fluidFlagFieldEvaluationFilter, vtkLowerQCriterionLimit);
+ qCriterionVTK->addCellInclusionFilter(qFilter);
+ qCriterionVTK->addBeforeFunction(velocityFieldWriter);
+ qCriterionVTK->addBeforeFunction(*velocityCommunicationScheme);
+ qCriterionVTK->addCellDataWriter(make_shared<lbm::QCriterionVTKWriter<VelocityField_T, FluidFilter_T, float>>(blocks,
+ fluidFlagFieldEvaluationFilter, velocityFieldID, "QCriterionFromPDF"));
+ qCriterionVTK->addCellDataWriter(make_shared<lbm::VorticityComponentVTKWriter<VelocityField_T, FluidFilter_T>>(blocks,
+ fluidFlagFieldEvaluationFilter, velocityFieldID, uint_t(2), "VorticityZFromPDF", u_g/diameter));
+
+ // domain decomposition
+ auto domainDecompVTK = vtk::createVTKOutput_DomainDecomposition(blocks, "domain_decomposition"+vtkFileTag,
+ vtkIOFreq, baseFolder, "simulation_step", false, true, true, true,
+ vtkInitialExecutionCount);
+
+ timeloop.addFuncBeforeTimeStep(vtk::writeFiles(particleVtkWriter), "VTK (sparse data)");
+ if (vtkIOSelection != "none") {
+ timeloop.addFuncBeforeTimeStep(vtk::writeFiles(pdfFieldVTK), "VTK (fluid field data)");
+ }
+ timeloop.addFuncAfterTimeStep(vtk::writeFiles(domainDecompVTK), "VTK (domain decomposition)");
+
+ if (fullVtkIOFreq != uint_t(0)) {
+ timeloop.addFuncBeforeTimeStep(vtk::writeFiles(fullPdfFieldVTK), "VTK (full fluid field data)");
+ }
+
+ timeloop.addFuncBeforeTimeStep(vtk::writeFiles(qCriterionVTK), "VTK (q value data)");
+
+ timeloop.addFuncAfterTimeStep( RemainingTimeLogger(timeloop.getNrOfTimeSteps()), "Remaining Time Logger" );
+
+ real_t terminationPosition;
+ // when the sphere is only a diameter away from the initial position again
+ if (rising) {
+ terminationPosition = real_t(initialPosition[2]) - diameter;
+ } else {
+ terminationPosition = real_t(initialPosition[2]) + diameter;
+ }
+ if (terminationPosition < real_t(0)) {
+ terminationPosition = real_t(domainSize[2]) + terminationPosition;
+ } else if (terminationPosition > real_t(domainSize[2])) {
+ terminationPosition = terminationPosition - real_t(domainSize[2]);
+ }
+ WALBERLA_LOG_INFO_ON_ROOT(" - termination position = " << terminationPosition );
+
+ std::vector<std::string> oldCheckpointFiles;
+
+ //// execute simulation
+ bool hasCrossedBoundary = false;
+
+ for (uint_t i = 0; i < timesteps; ++i) {
+ if (refinementCheckFrequency != 0 && i % refinementCheckFrequency == 0) {
+ WALBERLA_LOG_INFO_ON_ROOT("Refinement: Refreshing...");
+
+ lbm_mesapd_coupling::amr::updateAndSyncInfoCollection<BoundaryHandling_T, ParticleAccessor_T>(*forest, boundaryHandlingID,
+ *accessor, numMESAPDSubCycles, *couplingInfoCollection);
+
+ // for the fluid property based check
+ if (useCurlCriterion || useVorticityCriterion) {
+ velocityFieldWriter();
+ (*velocityCommunicationScheme)();
+ }
+
+ // clear mesapd
+ if (useSyncNextNeighbors) {
+ mesa_pd::mpi::ClearNextNeighborSync CNNS;
+ CNNS(*accessor);
+ } else {
+ mesa_pd::mpi::ClearGhostOwnerSync CGOS;
+ CGOS(*accessor);
+ }
+
+ forest->refresh();
+ mesapdDomain->refresh();
+ initializationSyncCall();
+
+ clearBoundaryHandling(*forest, boundaryHandlingID);
+ clearParticleField(*forest, particleFieldID, *accessor);
+ recreateBoundaryHandling(*forest, boundaryHandlingID);
+
+ mappingCall();
+
+ std::vector<uint_t> blocksPerLevel(numberOfLevels);
+ for (uint_t l = 0; l < numberOfLevels; l++) {
+ blocksPerLevel[l] = blocks->getNumberOfBlocks(l);
+ }
+ mpi::reduceInplace(blocksPerLevel, mpi::SUM);
+ WALBERLA_ROOT_SECTION() {
+ WALBERLA_LOG_INFO("Refinement: Number of blocks per level: ");
+ for (uint_t l = 0; l < numberOfLevels; l++) {
+ WALBERLA_LOG_INFO(" - Level " << l << ": " << blocksPerLevel[l] << " blocks");
+ }
+ }
+ }
+
+ // perform a single simulation step
+ timeloop.singleStep( *timeloopTiming );
+
+ // check for termination
+ // check if the sphere crossed a boundary
+ auto zPosition = logger.getPosition()[2];
+ if (!hasCrossedBoundary && ((rising && zPosition < terminationPosition) ||
+ (!rising && zPosition > terminationPosition))) {
+ hasCrossedBoundary = true;
+ }
+ if (hasCrossedBoundary && ((rising && zPosition > terminationPosition)
+ || (!rising && zPosition < terminationPosition))) {
+ // end of simulation reached after crossing a boundary and then passing the termination position
+ WALBERLA_LOG_INFO_ON_ROOT(
+ "Sphere reached terminal position " << logger.getPosition() << " after " << i << " timesteps!");
+ break;
+ }
+ if (checkPointingFreq != 0 && i % checkPointingFreq == 0) {
+ createCheckpointFiles(checkPointFileName, oldCheckpointFiles, timeloop.getCurrentTimeStep(),
+ mesapdTimeStep.getCurrentTimeStep(), forest, pdfFieldID, particleStorageID);
+ }
+ }
+
+ size_t particleIdx = accessor->uidToIdx(sphereUid);
+ if (particleIdx != accessor->getInvalidIdx()) {
+ if (!mesa_pd::data::particle_flags::isSet(accessor->getFlags(particleIdx), mesa_pd::data::particle_flags::GHOST)) {
+ Vector3<real_t> terminalVelocity = accessor->getLinearVelocity(particleIdx);
+ real_t reynoldsNumber = terminalVelocity[2] * diameter / viscosity;
+ WALBERLA_LOG_INFO("Terminal velocity: " << terminalVelocity);
+ WALBERLA_LOG_INFO("Terminal Reynolds number: " << reynoldsNumber);
+ }
+ }
+
+ timeloopTiming->logResultOnRoot();
+
+ return EXIT_SUCCESS;
+}
+
+}
+
+int main(int argc, char **argv) {
+ light_rising_particle_amr::main(argc, argv);
+}
\ No newline at end of file
diff --git a/python/mesa_pd.py b/python/mesa_pd.py
index 4e748a1d7..3161a1a27 100755
--- a/python/mesa_pd.py
+++ b/python/mesa_pd.py
@@ -63,6 +63,22 @@ if __name__ == '__main__':
ps.add_property("oldHydrodynamicTorque", "walberla::mesa_pd::Vec3", defValue="real_t(0)",
syncMode="ON_OWNERSHIP_CHANGE")
+ # Properties for virtual mass:
+ ps.add_property("virtualMass", "walberla::real_t", defValue="real_t(0)",
+ syncMode="ON_OWNERSHIP_CHANGE")
+ ps.add_property("invMassIncludingVirtual", "walberla::real_t", defValue="real_t(0)",
+ syncMode="ON_OWNERSHIP_CHANGE")
+ ps.add_property("oldLinearAcceleration", "walberla::mesa_pd::Vec3", defValue="real_t(0)",
+ syncMode="ON_OWNERSHIP_CHANGE")
+ ps.add_property("invInertiaBF", "walberla::mesa_pd::Mat3", defValue="real_t(0)",
+ syncMode="ON_OWNERSHIP_CHANGE")
+ ps.add_property("virtualInertiaBF", "walberla::mesa_pd::Mat3", defValue="real_t(0)",
+ syncMode="ON_OWNERSHIP_CHANGE")
+ ps.add_property("invInertiaBFIncludingVirtual", "walberla::mesa_pd::Mat3", defValue="real_t(0)",
+ syncMode="ON_OWNERSHIP_CHANGE")
+ ps.add_property("oldAngularAcceleration", "walberla::mesa_pd::Vec3", defValue="real_t(0)",
+ syncMode="ON_OWNERSHIP_CHANGE")
+
# properties for VBond model
ps.add_property("clusterID", "int64_t", defValue="-1", syncMode="ON_GHOST_CREATION")
ps.add_property("segmentID", "int64_t", defValue="-1", syncMode="ON_GHOST_CREATION")
diff --git a/src/lbm_mesapd_coupling/amr/level_determination/ParticlePresenceLevelDetermination.h b/src/lbm_mesapd_coupling/amr/level_determination/ParticlePresenceLevelDetermination.h
index 9ad15259e..deb9b2a7d 100644
--- a/src/lbm_mesapd_coupling/amr/level_determination/ParticlePresenceLevelDetermination.h
+++ b/src/lbm_mesapd_coupling/amr/level_determination/ParticlePresenceLevelDetermination.h
@@ -25,7 +25,8 @@
namespace walberla {
namespace lbm_mesapd_coupling {
-
+namespace amr
+{
/*
* Class to determine the minimum level a block can be.
* For coupled LBM-PE simulations the following rules apply:
@@ -39,21 +40,24 @@ namespace lbm_mesapd_coupling {
*/
class ParticlePresenceLevelDetermination
{
-public:
-
- ParticlePresenceLevelDetermination( const shared_ptr<lbm_mesapd_coupling::InfoCollection> & infoCollection, uint_t finestLevel) :
- infoCollection_( infoCollection ), finestLevel_( finestLevel)
+ public:
+ ParticlePresenceLevelDetermination(const shared_ptr< InfoCollection > infoCollection,
+ uint_t finestLevel)
+ : infoCollection_(infoCollection), finestLevel_(finestLevel)
{}
- void operator()( std::vector< std::pair< const Block *, uint_t > > & minTargetLevels,
- std::vector< const Block * > &, const BlockForest & /*forest*/ ) {
- for (auto &minTargetLevel : minTargetLevels) {
+ void operator()(std::vector< std::pair< const Block*, uint_t > >& minTargetLevels, std::vector< const Block* >&,
+ const BlockForest& /*forest*/)
+ {
+ for (auto& minTargetLevel : minTargetLevels)
+ {
uint_t currentLevelOfBlock = minTargetLevel.first->getLevel();
- const uint_t numberOfParticlesInDirectNeighborhood = getNumberOfLocalAndShadowParticlesInNeighborhood(minTargetLevel.first);
+ const uint_t numberOfParticlesInDirectNeighborhood =
+ getNumberOfLocalAndShadowParticlesInNeighborhood(minTargetLevel.first);
- uint_t targetLevelOfBlock = currentLevelOfBlock; //keep everything as it is
- if ( numberOfParticlesInDirectNeighborhood > uint_t(0) )
+ uint_t targetLevelOfBlock = currentLevelOfBlock; // keep everything as it is
+ if (numberOfParticlesInDirectNeighborhood > uint_t(0))
{
// set block to finest level if there are particles nearby
targetLevelOfBlock = finestLevel_;
@@ -61,41 +65,44 @@ public:
else
{
// block could coarsen since there are no particles nearby
- if( currentLevelOfBlock > uint_t(0) )
- targetLevelOfBlock = currentLevelOfBlock - uint_t(1);
+ if (currentLevelOfBlock > uint_t(0)) targetLevelOfBlock = currentLevelOfBlock - uint_t(1);
}
- WALBERLA_CHECK_LESS_EQUAL(std::abs(int_c(targetLevelOfBlock) - int_c(currentLevelOfBlock)), uint_t(1), "Only level difference of maximum 1 allowed!");
+ WALBERLA_CHECK_LESS_EQUAL(std::abs(int_c(targetLevelOfBlock) - int_c(currentLevelOfBlock)), uint_t(1),
+ "Only level difference of maximum 1 allowed!");
minTargetLevel.second = targetLevelOfBlock;
}
}
-private:
-
- uint_t getNumberOfLocalAndShadowParticlesInNeighborhood(const Block * block) {
+ private:
+ uint_t getNumberOfLocalAndShadowParticlesInNeighborhood(const Block* block)
+ {
auto numParticles = uint_t(0);
// add particles of current block
const auto infoIt = infoCollection_->find(block->getId());
- WALBERLA_CHECK_UNEQUAL(infoIt, infoCollection_->end(), "Block with ID " << block->getId() << " not found in info collection!");
+ WALBERLA_CHECK_UNEQUAL(infoIt, infoCollection_->end(),
+ "Block with ID " << block->getId() << " not found in info collection!");
- numParticles += infoIt->second.numberOfLocalParticles + infoIt->second.numberOfShadowParticles;
+ numParticles += infoIt->second.numberOfLocalParticles + infoIt->second.numberOfGhostParticles;
// add particles of all neighboring blocks
- for(uint_t i = 0; i < block->getNeighborhoodSize(); ++i)
+ for (uint_t i = 0; i < block->getNeighborhoodSize(); ++i)
{
- const BlockID &neighborBlockID = block->getNeighborId(i);
- const auto infoItNeighbor = infoCollection_->find(neighborBlockID);
- WALBERLA_CHECK_UNEQUAL(infoItNeighbor, infoCollection_->end(), "Neighbor block with ID " << neighborBlockID << " not found in info collection!");
+ const BlockID& neighborBlockID = block->getNeighborId(i);
+ const auto infoItNeighbor = infoCollection_->find(neighborBlockID);
+ WALBERLA_CHECK_UNEQUAL(infoItNeighbor, infoCollection_->end(),
+ "Neighbor block with ID " << neighborBlockID << " not found in info collection!");
- numParticles += infoItNeighbor->second.numberOfLocalParticles + infoItNeighbor->second.numberOfShadowParticles;
+ numParticles += infoItNeighbor->second.numberOfLocalParticles + infoItNeighbor->second.numberOfGhostParticles;
}
return numParticles;
}
- shared_ptr<lbm_mesapd_coupling::InfoCollection> infoCollection_;
+ shared_ptr< InfoCollection > infoCollection_;
uint_t finestLevel_;
};
+} // namespace amr
} // namespace lbm_mesapd_coupling
} // namespace walberla
diff --git a/src/lbm_mesapd_coupling/utility/virtualmass/AddVirtualForceTorqueKernel.h b/src/lbm_mesapd_coupling/utility/virtualmass/AddVirtualForceTorqueKernel.h
new file mode 100644
index 000000000..6f4a7ea35
--- /dev/null
+++ b/src/lbm_mesapd_coupling/utility/virtualmass/AddVirtualForceTorqueKernel.h
@@ -0,0 +1,88 @@
+//======================================================================================================================
+//
+// 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 AddVirtualForceTorqueKernel.h
+//! \ingroup lbm_mesapd_coupling
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <utility>
+
+#include "mesa_pd/data/DataTypes.h"
+#include "mesa_pd/data/IAccessor.h"
+#include "mesa_pd/common/ParticleFunctions.h"
+#include "mesa_pd/data/ParticleStorage.h"
+
+namespace walberla {
+namespace lbm_mesapd_coupling {
+
+/**
+ * Kernel that sets a virtual force and torque on particles. It accesses a virtual mass and inertia, which have to be
+ * set beforehand using InitializeVirtualMassKernel.
+ * During the calculation of virtual force and torque linear and angular accelerations are used, for which acceleration
+ * estimators need to be supplied.
+ *
+ * This kernel requires the following particle attributes:
+ * ps.addProperty("virtualMass", "walberla::real_t", defValue="real_t(0)", syncMode="ON_OWNERSHIP_CHANGE")
+ * ps.addProperty("virtualInertiaBF", "walberla::mesa_pd::Mat3", defValue="real_t(0)", syncMode="ON_OWNERSHIP_CHANGE")
+ * ps.add_property("oldLinearAcceleration", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ON_OWNERSHIP_CHANGE")
+ * ps.add_property("oldAngularAcceleration", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ON_OWNERSHIP_CHANGE")
+ */
+class AddVirtualForceTorqueKernel
+{
+public:
+ explicit AddVirtualForceTorqueKernel(shared_ptr<mesa_pd::data::ParticleStorage> ps) : ps_(std::move(ps)){};
+
+ template <typename Accessor_T>
+ void operator()(const size_t idx, Accessor_T& ac) {
+ static_assert(std::is_base_of<mesa_pd::data::IAccessor, Accessor_T>::value, "please provide a valid accessor");
+
+ WALBERLA_CHECK_FLOAT_UNEQUAL(ac.getVirtualMass(idx), real_t(0.),
+ "No virtual mass set on body. Was the virtualMass kernel not called before?");
+
+ //// force
+ // obtain the acceleration estimation
+ Vector3<real_t> approximatedCurrentLinearAcceleration = ac.getOldLinearAcceleration(idx);
+
+ // prepare the next acceleration estimation for the following time step
+ const Vector3<real_t> virtualForce = ac.getVirtualMass(idx) * approximatedCurrentLinearAcceleration;
+ const Vector3<real_t> currentLinearAcceleration = (ac.getForce(idx) + virtualForce) /
+ (ac.getVirtualMass(idx) + ac.getMass(idx));
+ ac.setOldLinearAcceleration(idx, currentLinearAcceleration);
+
+ mesa_pd::addForceAtomic(idx, ac, virtualForce);
+
+ //// torque
+ // obtain the acceleration estimation
+ Vector3<real_t> approximatedCurrentAngularAcceleration = ac.getOldAngularAcceleration(idx);
+
+ // prepare the next acceleration estimation for the following time step
+ const Vector3<real_t> virtualTorque = ac.getVirtualInertiaBF(idx) * approximatedCurrentAngularAcceleration;
+ const Vector3<real_t> angularAcceleration = math::transformMatrixRART(ac.getRotation(idx).getMatrix(),
+ ac.getInvInertiaBFIncludingVirtual(idx)) * (ac.getTorque(idx) + virtualTorque);
+ ac.setOldAngularAcceleration(idx, angularAcceleration);
+
+ mesa_pd::addTorqueAtomic(idx, ac, virtualTorque);
+ }
+
+private:
+ const shared_ptr<mesa_pd::data::ParticleStorage> ps_;
+};
+
+}
+}
\ No newline at end of file
diff --git a/src/lbm_mesapd_coupling/utility/virtualmass/InitializeVirtualMassKernel.h b/src/lbm_mesapd_coupling/utility/virtualmass/InitializeVirtualMassKernel.h
new file mode 100644
index 000000000..f03cf4274
--- /dev/null
+++ b/src/lbm_mesapd_coupling/utility/virtualmass/InitializeVirtualMassKernel.h
@@ -0,0 +1,64 @@
+//======================================================================================================================
+//
+// 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 VirtualMass.h
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/IAccessor.h>
+
+namespace walberla {
+namespace lbm_mesapd_coupling {
+
+/**
+ * This kernel calculates virtual mass and inertia and sets them as attributes on a particle.
+ *
+ * It requires the following particle attributes:
+ * ps.addProperty("virtualMass", "walberla::real_t", defValue="real_t(0)", syncMode="ON_OWNERSHIP_CHANGE")
+ * ps.addProperty("invMassIncludingVirtual", "walberla::real_t", defValue="real_t(0)", syncMode="ON_OWNERSHIP_CHANGE")
+ * ps.addProperty("virtualInertiaBF", "walberla::mesa_pd::Mat3", defValue="real_t(0)", syncMode="ON_OWNERSHIP_CHANGE")
+ * ps.addProperty("invInertiaBFIncludingVirtual", "walberla::mesa_pd::Mat3", defValue="real_t(0)", syncMode="ON_OWNERSHIP_CHANGE")
+ */
+class InitializeVirtualMassKernel {
+public:
+ InitializeVirtualMassKernel() = default;
+
+ template <typename Accessor>
+ void operator()(size_t i, Accessor& ac, real_t C_v, real_t C_v_omega, real_t fluidDensity) const;
+};
+
+
+template <typename Accessor>
+inline void InitializeVirtualMassKernel::operator()(const size_t i, Accessor& ac, const real_t C_v, const real_t C_v_omega,
+ const real_t fluidDensity) const {
+ static_assert(std::is_base_of<mesa_pd::data::IAccessor, Accessor>::value, "please provide a valid accessor");
+
+ const real_t virtualMass = C_v * fluidDensity * ac.getVolume(i);
+ ac.setVirtualMass(i, virtualMass);
+ ac.setInvMassIncludingVirtual(i, real_t(1.) / (ac.getMass(i) + virtualMass));
+
+ const real_t angularVirtualMass = C_v_omega * fluidDensity * ac.getVolume(i);
+ const mesa_pd::Mat3 virtualInertiaBF = ac.getInertiaBF(i) * ac.getInvMass(i) * angularVirtualMass;
+ const mesa_pd::Mat3 inertiaBFIncludingVirtual = ac.getInertiaBF(i) + virtualInertiaBF;
+ ac.setVirtualInertiaBF(i, virtualInertiaBF);
+ ac.setInvInertiaBFIncludingVirtual(i, inertiaBFIncludingVirtual.getInverse());
+}
+
+} //namespace lbm_mesapd_coupling
+} //namespace walberla
\ No newline at end of file
diff --git a/src/lbm_mesapd_coupling/utility/virtualmass/ParticleAccessorWithShapeVirtualMassWrapper.h b/src/lbm_mesapd_coupling/utility/virtualmass/ParticleAccessorWithShapeVirtualMassWrapper.h
new file mode 100644
index 000000000..744a614fb
--- /dev/null
+++ b/src/lbm_mesapd_coupling/utility/virtualmass/ParticleAccessorWithShapeVirtualMassWrapper.h
@@ -0,0 +1,58 @@
+//======================================================================================================================
+//
+// 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 ParticleAccessorWithShapeAndVirtualMassWrapper.h
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/ParticleStorage.h>
+#include <mesa_pd/data/ShapeStorage.h>
+
+namespace walberla {
+namespace lbm_mesapd_coupling {
+
+/**
+ * !\brief Wraps a ParticleAccessor to change its getInvMass and getInvInertiaBF methods to include virtual mass.
+ * @tparam T A ParticleAccessor providing functions getInvMass and getInvInertiaBF, and their respective getInvMassIncludingVirtual and getInvInertiaBFIncludingVirtual.
+ */
+template <typename T>
+class ParticleAccessorWithShapeVirtualMassWrapper : public T {
+public:
+ ParticleAccessorWithShapeVirtualMassWrapper(std::shared_ptr<mesa_pd::data::ParticleStorage>& ps, std::shared_ptr<mesa_pd::data::ShapeStorage>& ss)
+ : T(ps, ss) {
+ }
+
+ auto getInvMass(const size_t p_idx) const {
+ return T::getInvMassIncludingVirtual(p_idx);
+ }
+
+ const auto getMass(const size_t p_idx) const {
+ return T::getMass(p_idx) + T::getVirtualMass(p_idx);
+ }
+
+ auto getInvInertiaBF(const size_t p_idx) const {
+ return T::getInvInertiaBFIncludingVirtual(p_idx);
+ }
+
+ const auto getInertiaBF(const size_t p_idx) const {
+ return T::getInertiaBF(p_idx) + T::getVirtualInertiaBF(p_idx);
+ }
+};
+
+}
+}
\ No newline at end of file
diff --git a/src/mesa_pd/data/ParticleAccessor.h b/src/mesa_pd/data/ParticleAccessor.h
index 62c61f6f7..14ed03dc2 100644
--- a/src/mesa_pd/data/ParticleAccessor.h
+++ b/src/mesa_pd/data/ParticleAccessor.h
@@ -164,6 +164,34 @@ public:
walberla::mesa_pd::Vec3& getOldHydrodynamicTorqueRef(const size_t p_idx) {return ps_->getOldHydrodynamicTorqueRef(p_idx);}
void setOldHydrodynamicTorque(const size_t p_idx, walberla::mesa_pd::Vec3 const & v) { ps_->setOldHydrodynamicTorque(p_idx, v);}
+ walberla::real_t const & getVirtualMass(const size_t p_idx) const {return ps_->getVirtualMass(p_idx);}
+ walberla::real_t& getVirtualMassRef(const size_t p_idx) {return ps_->getVirtualMassRef(p_idx);}
+ void setVirtualMass(const size_t p_idx, walberla::real_t const & v) { ps_->setVirtualMass(p_idx, v);}
+
+ walberla::real_t const & getInvMassIncludingVirtual(const size_t p_idx) const {return ps_->getInvMassIncludingVirtual(p_idx);}
+ walberla::real_t& getInvMassIncludingVirtualRef(const size_t p_idx) {return ps_->getInvMassIncludingVirtualRef(p_idx);}
+ void setInvMassIncludingVirtual(const size_t p_idx, walberla::real_t const & v) { ps_->setInvMassIncludingVirtual(p_idx, v);}
+
+ walberla::mesa_pd::Vec3 const & getOldLinearAcceleration(const size_t p_idx) const {return ps_->getOldLinearAcceleration(p_idx);}
+ walberla::mesa_pd::Vec3& getOldLinearAccelerationRef(const size_t p_idx) {return ps_->getOldLinearAccelerationRef(p_idx);}
+ void setOldLinearAcceleration(const size_t p_idx, walberla::mesa_pd::Vec3 const & v) { ps_->setOldLinearAcceleration(p_idx, v);}
+
+ walberla::mesa_pd::Mat3 const & getInvInertiaBF(const size_t p_idx) const {return ps_->getInvInertiaBF(p_idx);}
+ walberla::mesa_pd::Mat3& getInvInertiaBFRef(const size_t p_idx) {return ps_->getInvInertiaBFRef(p_idx);}
+ void setInvInertiaBF(const size_t p_idx, walberla::mesa_pd::Mat3 const & v) { ps_->setInvInertiaBF(p_idx, v);}
+
+ walberla::mesa_pd::Mat3 const & getVirtualInertiaBF(const size_t p_idx) const {return ps_->getVirtualInertiaBF(p_idx);}
+ walberla::mesa_pd::Mat3& getVirtualInertiaBFRef(const size_t p_idx) {return ps_->getVirtualInertiaBFRef(p_idx);}
+ void setVirtualInertiaBF(const size_t p_idx, walberla::mesa_pd::Mat3 const & v) { ps_->setVirtualInertiaBF(p_idx, v);}
+
+ walberla::mesa_pd::Mat3 const & getInvInertiaBFIncludingVirtual(const size_t p_idx) const {return ps_->getInvInertiaBFIncludingVirtual(p_idx);}
+ walberla::mesa_pd::Mat3& getInvInertiaBFIncludingVirtualRef(const size_t p_idx) {return ps_->getInvInertiaBFIncludingVirtualRef(p_idx);}
+ void setInvInertiaBFIncludingVirtual(const size_t p_idx, walberla::mesa_pd::Mat3 const & v) { ps_->setInvInertiaBFIncludingVirtual(p_idx, v);}
+
+ walberla::mesa_pd::Vec3 const & getOldAngularAcceleration(const size_t p_idx) const {return ps_->getOldAngularAcceleration(p_idx);}
+ walberla::mesa_pd::Vec3& getOldAngularAccelerationRef(const size_t p_idx) {return ps_->getOldAngularAccelerationRef(p_idx);}
+ void setOldAngularAcceleration(const size_t p_idx, walberla::mesa_pd::Vec3 const & v) { ps_->setOldAngularAcceleration(p_idx, v);}
+
int64_t const & getClusterID(const size_t p_idx) const {return ps_->getClusterID(p_idx);}
int64_t& getClusterIDRef(const size_t p_idx) {return ps_->getClusterIDRef(p_idx);}
void setClusterID(const size_t p_idx, int64_t const & v) { ps_->setClusterID(p_idx, v);}
@@ -337,6 +365,34 @@ public:
void setOldHydrodynamicTorque(const size_t /*p_idx*/, walberla::mesa_pd::Vec3 const & v) { oldHydrodynamicTorque_ = v;}
walberla::mesa_pd::Vec3& getOldHydrodynamicTorqueRef(const size_t /*p_idx*/) {return oldHydrodynamicTorque_;}
+ walberla::real_t const & getVirtualMass(const size_t /*p_idx*/) const {return virtualMass_;}
+ void setVirtualMass(const size_t /*p_idx*/, walberla::real_t const & v) { virtualMass_ = v;}
+ walberla::real_t& getVirtualMassRef(const size_t /*p_idx*/) {return virtualMass_;}
+
+ walberla::real_t const & getInvMassIncludingVirtual(const size_t /*p_idx*/) const {return invMassIncludingVirtual_;}
+ void setInvMassIncludingVirtual(const size_t /*p_idx*/, walberla::real_t const & v) { invMassIncludingVirtual_ = v;}
+ walberla::real_t& getInvMassIncludingVirtualRef(const size_t /*p_idx*/) {return invMassIncludingVirtual_;}
+
+ walberla::mesa_pd::Vec3 const & getOldLinearAcceleration(const size_t /*p_idx*/) const {return oldLinearAcceleration_;}
+ void setOldLinearAcceleration(const size_t /*p_idx*/, walberla::mesa_pd::Vec3 const & v) { oldLinearAcceleration_ = v;}
+ walberla::mesa_pd::Vec3& getOldLinearAccelerationRef(const size_t /*p_idx*/) {return oldLinearAcceleration_;}
+
+ walberla::mesa_pd::Mat3 const & getInvInertiaBF(const size_t /*p_idx*/) const {return invInertiaBF_;}
+ void setInvInertiaBF(const size_t /*p_idx*/, walberla::mesa_pd::Mat3 const & v) { invInertiaBF_ = v;}
+ walberla::mesa_pd::Mat3& getInvInertiaBFRef(const size_t /*p_idx*/) {return invInertiaBF_;}
+
+ walberla::mesa_pd::Mat3 const & getVirtualInertiaBF(const size_t /*p_idx*/) const {return virtualInertiaBF_;}
+ void setVirtualInertiaBF(const size_t /*p_idx*/, walberla::mesa_pd::Mat3 const & v) { virtualInertiaBF_ = v;}
+ walberla::mesa_pd::Mat3& getVirtualInertiaBFRef(const size_t /*p_idx*/) {return virtualInertiaBF_;}
+
+ walberla::mesa_pd::Mat3 const & getInvInertiaBFIncludingVirtual(const size_t /*p_idx*/) const {return invInertiaBFIncludingVirtual_;}
+ void setInvInertiaBFIncludingVirtual(const size_t /*p_idx*/, walberla::mesa_pd::Mat3 const & v) { invInertiaBFIncludingVirtual_ = v;}
+ walberla::mesa_pd::Mat3& getInvInertiaBFIncludingVirtualRef(const size_t /*p_idx*/) {return invInertiaBFIncludingVirtual_;}
+
+ walberla::mesa_pd::Vec3 const & getOldAngularAcceleration(const size_t /*p_idx*/) const {return oldAngularAcceleration_;}
+ void setOldAngularAcceleration(const size_t /*p_idx*/, walberla::mesa_pd::Vec3 const & v) { oldAngularAcceleration_ = v;}
+ walberla::mesa_pd::Vec3& getOldAngularAccelerationRef(const size_t /*p_idx*/) {return oldAngularAcceleration_;}
+
int64_t const & getClusterID(const size_t /*p_idx*/) const {return clusterID_;}
void setClusterID(const size_t /*p_idx*/, int64_t const & v) { clusterID_ = v;}
int64_t& getClusterIDRef(const size_t /*p_idx*/) {return clusterID_;}
@@ -389,6 +445,13 @@ private:
walberla::mesa_pd::Vec3 hydrodynamicTorque_;
walberla::mesa_pd::Vec3 oldHydrodynamicForce_;
walberla::mesa_pd::Vec3 oldHydrodynamicTorque_;
+ walberla::real_t virtualMass_;
+ walberla::real_t invMassIncludingVirtual_;
+ walberla::mesa_pd::Vec3 oldLinearAcceleration_;
+ walberla::mesa_pd::Mat3 invInertiaBF_;
+ walberla::mesa_pd::Mat3 virtualInertiaBF_;
+ walberla::mesa_pd::Mat3 invInertiaBFIncludingVirtual_;
+ walberla::mesa_pd::Vec3 oldAngularAcceleration_;
int64_t clusterID_;
int64_t segmentID_;
std::unordered_set<walberla::mpi::MPIRank> neighborState_;
diff --git a/src/mesa_pd/data/ParticleAccessorWithShape.h b/src/mesa_pd/data/ParticleAccessorWithShape.h
index 7f3413c12..8c884d370 100644
--- a/src/mesa_pd/data/ParticleAccessorWithShape.h
+++ b/src/mesa_pd/data/ParticleAccessorWithShape.h
@@ -37,7 +37,10 @@ public:
{}
const auto& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvMass();}
+ const auto& getMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getMass();}
const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvInertiaBF();}
+ const auto& getInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInertiaBF();}
+ const auto getVolume(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getVolume();}
data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)].get();}
private:
diff --git a/src/mesa_pd/data/ParticleStorage.h b/src/mesa_pd/data/ParticleStorage.h
index 81d628681..4d64d3405 100644
--- a/src/mesa_pd/data/ParticleStorage.h
+++ b/src/mesa_pd/data/ParticleStorage.h
@@ -100,6 +100,13 @@ public:
using hydrodynamicTorque_type = walberla::mesa_pd::Vec3;
using oldHydrodynamicForce_type = walberla::mesa_pd::Vec3;
using oldHydrodynamicTorque_type = walberla::mesa_pd::Vec3;
+ using virtualMass_type = walberla::real_t;
+ using invMassIncludingVirtual_type = walberla::real_t;
+ using oldLinearAcceleration_type = walberla::mesa_pd::Vec3;
+ using invInertiaBF_type = walberla::mesa_pd::Mat3;
+ using virtualInertiaBF_type = walberla::mesa_pd::Mat3;
+ using invInertiaBFIncludingVirtual_type = walberla::mesa_pd::Mat3;
+ using oldAngularAcceleration_type = walberla::mesa_pd::Vec3;
using clusterID_type = int64_t;
using segmentID_type = int64_t;
using neighborState_type = std::unordered_set<walberla::mpi::MPIRank>;
@@ -221,6 +228,34 @@ public:
oldHydrodynamicTorque_type& getOldHydrodynamicTorqueRef() {return storage_.getOldHydrodynamicTorqueRef(i_);}
void setOldHydrodynamicTorque(oldHydrodynamicTorque_type const & v) { storage_.setOldHydrodynamicTorque(i_, v);}
+ virtualMass_type const & getVirtualMass() const {return storage_.getVirtualMass(i_);}
+ virtualMass_type& getVirtualMassRef() {return storage_.getVirtualMassRef(i_);}
+ void setVirtualMass(virtualMass_type const & v) { storage_.setVirtualMass(i_, v);}
+
+ invMassIncludingVirtual_type const & getInvMassIncludingVirtual() const {return storage_.getInvMassIncludingVirtual(i_);}
+ invMassIncludingVirtual_type& getInvMassIncludingVirtualRef() {return storage_.getInvMassIncludingVirtualRef(i_);}
+ void setInvMassIncludingVirtual(invMassIncludingVirtual_type const & v) { storage_.setInvMassIncludingVirtual(i_, v);}
+
+ oldLinearAcceleration_type const & getOldLinearAcceleration() const {return storage_.getOldLinearAcceleration(i_);}
+ oldLinearAcceleration_type& getOldLinearAccelerationRef() {return storage_.getOldLinearAccelerationRef(i_);}
+ void setOldLinearAcceleration(oldLinearAcceleration_type const & v) { storage_.setOldLinearAcceleration(i_, v);}
+
+ invInertiaBF_type const & getInvInertiaBF() const {return storage_.getInvInertiaBF(i_);}
+ invInertiaBF_type& getInvInertiaBFRef() {return storage_.getInvInertiaBFRef(i_);}
+ void setInvInertiaBF(invInertiaBF_type const & v) { storage_.setInvInertiaBF(i_, v);}
+
+ virtualInertiaBF_type const & getVirtualInertiaBF() const {return storage_.getVirtualInertiaBF(i_);}
+ virtualInertiaBF_type& getVirtualInertiaBFRef() {return storage_.getVirtualInertiaBFRef(i_);}
+ void setVirtualInertiaBF(virtualInertiaBF_type const & v) { storage_.setVirtualInertiaBF(i_, v);}
+
+ invInertiaBFIncludingVirtual_type const & getInvInertiaBFIncludingVirtual() const {return storage_.getInvInertiaBFIncludingVirtual(i_);}
+ invInertiaBFIncludingVirtual_type& getInvInertiaBFIncludingVirtualRef() {return storage_.getInvInertiaBFIncludingVirtualRef(i_);}
+ void setInvInertiaBFIncludingVirtual(invInertiaBFIncludingVirtual_type const & v) { storage_.setInvInertiaBFIncludingVirtual(i_, v);}
+
+ oldAngularAcceleration_type const & getOldAngularAcceleration() const {return storage_.getOldAngularAcceleration(i_);}
+ oldAngularAcceleration_type& getOldAngularAccelerationRef() {return storage_.getOldAngularAccelerationRef(i_);}
+ void setOldAngularAcceleration(oldAngularAcceleration_type const & v) { storage_.setOldAngularAcceleration(i_, v);}
+
clusterID_type const & getClusterID() const {return storage_.getClusterID(i_);}
clusterID_type& getClusterIDRef() {return storage_.getClusterIDRef(i_);}
void setClusterID(clusterID_type const & v) { storage_.setClusterID(i_, v);}
@@ -322,6 +357,13 @@ public:
using hydrodynamicTorque_type = walberla::mesa_pd::Vec3;
using oldHydrodynamicForce_type = walberla::mesa_pd::Vec3;
using oldHydrodynamicTorque_type = walberla::mesa_pd::Vec3;
+ using virtualMass_type = walberla::real_t;
+ using invMassIncludingVirtual_type = walberla::real_t;
+ using oldLinearAcceleration_type = walberla::mesa_pd::Vec3;
+ using invInertiaBF_type = walberla::mesa_pd::Mat3;
+ using virtualInertiaBF_type = walberla::mesa_pd::Mat3;
+ using invInertiaBFIncludingVirtual_type = walberla::mesa_pd::Mat3;
+ using oldAngularAcceleration_type = walberla::mesa_pd::Vec3;
using clusterID_type = int64_t;
using segmentID_type = int64_t;
using neighborState_type = std::unordered_set<walberla::mpi::MPIRank>;
@@ -443,6 +485,34 @@ public:
oldHydrodynamicTorque_type& getOldHydrodynamicTorqueRef(const size_t idx) {return oldHydrodynamicTorque_[idx];}
void setOldHydrodynamicTorque(const size_t idx, oldHydrodynamicTorque_type const & v) { oldHydrodynamicTorque_[idx] = v; }
+ virtualMass_type const & getVirtualMass(const size_t idx) const {return virtualMass_[idx];}
+ virtualMass_type& getVirtualMassRef(const size_t idx) {return virtualMass_[idx];}
+ void setVirtualMass(const size_t idx, virtualMass_type const & v) { virtualMass_[idx] = v; }
+
+ invMassIncludingVirtual_type const & getInvMassIncludingVirtual(const size_t idx) const {return invMassIncludingVirtual_[idx];}
+ invMassIncludingVirtual_type& getInvMassIncludingVirtualRef(const size_t idx) {return invMassIncludingVirtual_[idx];}
+ void setInvMassIncludingVirtual(const size_t idx, invMassIncludingVirtual_type const & v) { invMassIncludingVirtual_[idx] = v; }
+
+ oldLinearAcceleration_type const & getOldLinearAcceleration(const size_t idx) const {return oldLinearAcceleration_[idx];}
+ oldLinearAcceleration_type& getOldLinearAccelerationRef(const size_t idx) {return oldLinearAcceleration_[idx];}
+ void setOldLinearAcceleration(const size_t idx, oldLinearAcceleration_type const & v) { oldLinearAcceleration_[idx] = v; }
+
+ invInertiaBF_type const & getInvInertiaBF(const size_t idx) const {return invInertiaBF_[idx];}
+ invInertiaBF_type& getInvInertiaBFRef(const size_t idx) {return invInertiaBF_[idx];}
+ void setInvInertiaBF(const size_t idx, invInertiaBF_type const & v) { invInertiaBF_[idx] = v; }
+
+ virtualInertiaBF_type const & getVirtualInertiaBF(const size_t idx) const {return virtualInertiaBF_[idx];}
+ virtualInertiaBF_type& getVirtualInertiaBFRef(const size_t idx) {return virtualInertiaBF_[idx];}
+ void setVirtualInertiaBF(const size_t idx, virtualInertiaBF_type const & v) { virtualInertiaBF_[idx] = v; }
+
+ invInertiaBFIncludingVirtual_type const & getInvInertiaBFIncludingVirtual(const size_t idx) const {return invInertiaBFIncludingVirtual_[idx];}
+ invInertiaBFIncludingVirtual_type& getInvInertiaBFIncludingVirtualRef(const size_t idx) {return invInertiaBFIncludingVirtual_[idx];}
+ void setInvInertiaBFIncludingVirtual(const size_t idx, invInertiaBFIncludingVirtual_type const & v) { invInertiaBFIncludingVirtual_[idx] = v; }
+
+ oldAngularAcceleration_type const & getOldAngularAcceleration(const size_t idx) const {return oldAngularAcceleration_[idx];}
+ oldAngularAcceleration_type& getOldAngularAccelerationRef(const size_t idx) {return oldAngularAcceleration_[idx];}
+ void setOldAngularAcceleration(const size_t idx, oldAngularAcceleration_type const & v) { oldAngularAcceleration_[idx] = v; }
+
clusterID_type const & getClusterID(const size_t idx) const {return clusterID_[idx];}
clusterID_type& getClusterIDRef(const size_t idx) {return clusterID_[idx];}
void setClusterID(const size_t idx, clusterID_type const & v) { clusterID_[idx] = v; }
@@ -575,6 +645,13 @@ public:
std::vector<hydrodynamicTorque_type> hydrodynamicTorque_ {};
std::vector<oldHydrodynamicForce_type> oldHydrodynamicForce_ {};
std::vector<oldHydrodynamicTorque_type> oldHydrodynamicTorque_ {};
+ std::vector<virtualMass_type> virtualMass_ {};
+ std::vector<invMassIncludingVirtual_type> invMassIncludingVirtual_ {};
+ std::vector<oldLinearAcceleration_type> oldLinearAcceleration_ {};
+ std::vector<invInertiaBF_type> invInertiaBF_ {};
+ std::vector<virtualInertiaBF_type> virtualInertiaBF_ {};
+ std::vector<invInertiaBFIncludingVirtual_type> invInertiaBFIncludingVirtual_ {};
+ std::vector<oldAngularAcceleration_type> oldAngularAcceleration_ {};
std::vector<clusterID_type> clusterID_ {};
std::vector<segmentID_type> segmentID_ {};
std::vector<neighborState_type> neighborState_ {};
@@ -616,6 +693,13 @@ ParticleStorage::Particle& ParticleStorage::Particle::operator=(const ParticleSt
getHydrodynamicTorqueRef() = rhs.getHydrodynamicTorque();
getOldHydrodynamicForceRef() = rhs.getOldHydrodynamicForce();
getOldHydrodynamicTorqueRef() = rhs.getOldHydrodynamicTorque();
+ getVirtualMassRef() = rhs.getVirtualMass();
+ getInvMassIncludingVirtualRef() = rhs.getInvMassIncludingVirtual();
+ getOldLinearAccelerationRef() = rhs.getOldLinearAcceleration();
+ getInvInertiaBFRef() = rhs.getInvInertiaBF();
+ getVirtualInertiaBFRef() = rhs.getVirtualInertiaBF();
+ getInvInertiaBFIncludingVirtualRef() = rhs.getInvInertiaBFIncludingVirtual();
+ getOldAngularAccelerationRef() = rhs.getOldAngularAcceleration();
getClusterIDRef() = rhs.getClusterID();
getSegmentIDRef() = rhs.getSegmentID();
getNeighborStateRef() = rhs.getNeighborState();
@@ -654,6 +738,13 @@ ParticleStorage::Particle& ParticleStorage::Particle::operator=(ParticleStorage:
getHydrodynamicTorqueRef() = std::move(rhs.getHydrodynamicTorqueRef());
getOldHydrodynamicForceRef() = std::move(rhs.getOldHydrodynamicForceRef());
getOldHydrodynamicTorqueRef() = std::move(rhs.getOldHydrodynamicTorqueRef());
+ getVirtualMassRef() = std::move(rhs.getVirtualMassRef());
+ getInvMassIncludingVirtualRef() = std::move(rhs.getInvMassIncludingVirtualRef());
+ getOldLinearAccelerationRef() = std::move(rhs.getOldLinearAccelerationRef());
+ getInvInertiaBFRef() = std::move(rhs.getInvInertiaBFRef());
+ getVirtualInertiaBFRef() = std::move(rhs.getVirtualInertiaBFRef());
+ getInvInertiaBFIncludingVirtualRef() = std::move(rhs.getInvInertiaBFIncludingVirtualRef());
+ getOldAngularAccelerationRef() = std::move(rhs.getOldAngularAccelerationRef());
getClusterIDRef() = std::move(rhs.getClusterIDRef());
getSegmentIDRef() = std::move(rhs.getSegmentIDRef());
getNeighborStateRef() = std::move(rhs.getNeighborStateRef());
@@ -693,6 +784,13 @@ void swap(ParticleStorage::Particle lhs, ParticleStorage::Particle rhs)
std::swap(lhs.getHydrodynamicTorqueRef(), rhs.getHydrodynamicTorqueRef());
std::swap(lhs.getOldHydrodynamicForceRef(), rhs.getOldHydrodynamicForceRef());
std::swap(lhs.getOldHydrodynamicTorqueRef(), rhs.getOldHydrodynamicTorqueRef());
+ std::swap(lhs.getVirtualMassRef(), rhs.getVirtualMassRef());
+ std::swap(lhs.getInvMassIncludingVirtualRef(), rhs.getInvMassIncludingVirtualRef());
+ std::swap(lhs.getOldLinearAccelerationRef(), rhs.getOldLinearAccelerationRef());
+ std::swap(lhs.getInvInertiaBFRef(), rhs.getInvInertiaBFRef());
+ std::swap(lhs.getVirtualInertiaBFRef(), rhs.getVirtualInertiaBFRef());
+ std::swap(lhs.getInvInertiaBFIncludingVirtualRef(), rhs.getInvInertiaBFIncludingVirtualRef());
+ std::swap(lhs.getOldAngularAccelerationRef(), rhs.getOldAngularAccelerationRef());
std::swap(lhs.getClusterIDRef(), rhs.getClusterIDRef());
std::swap(lhs.getSegmentIDRef(), rhs.getSegmentIDRef());
std::swap(lhs.getNeighborStateRef(), rhs.getNeighborStateRef());
@@ -732,6 +830,13 @@ std::ostream& operator<<( std::ostream& os, const ParticleStorage::Particle& p )
"hydrodynamicTorque : " << p.getHydrodynamicTorque() << "\n" <<
"oldHydrodynamicForce: " << p.getOldHydrodynamicForce() << "\n" <<
"oldHydrodynamicTorque: " << p.getOldHydrodynamicTorque() << "\n" <<
+ "virtualMass : " << p.getVirtualMass() << "\n" <<
+ "invMassIncludingVirtual: " << p.getInvMassIncludingVirtual() << "\n" <<
+ "oldLinearAcceleration: " << p.getOldLinearAcceleration() << "\n" <<
+ "invInertiaBF : " << p.getInvInertiaBF() << "\n" <<
+ "virtualInertiaBF : " << p.getVirtualInertiaBF() << "\n" <<
+ "invInertiaBFIncludingVirtual: " << p.getInvInertiaBFIncludingVirtual() << "\n" <<
+ "oldAngularAcceleration: " << p.getOldAngularAcceleration() << "\n" <<
"clusterID : " << p.getClusterID() << "\n" <<
"segmentID : " << p.getSegmentID() << "\n" <<
"neighborState : " << p.getNeighborState() << "\n" <<
@@ -841,6 +946,13 @@ inline ParticleStorage::iterator ParticleStorage::create(const id_t& uid)
hydrodynamicTorque_.emplace_back(real_t(0));
oldHydrodynamicForce_.emplace_back(real_t(0));
oldHydrodynamicTorque_.emplace_back(real_t(0));
+ virtualMass_.emplace_back(real_t(0));
+ invMassIncludingVirtual_.emplace_back(real_t(0));
+ oldLinearAcceleration_.emplace_back(real_t(0));
+ invInertiaBF_.emplace_back(real_t(0));
+ virtualInertiaBF_.emplace_back(real_t(0));
+ invInertiaBFIncludingVirtual_.emplace_back(real_t(0));
+ oldAngularAcceleration_.emplace_back(real_t(0));
clusterID_.emplace_back(-1);
segmentID_.emplace_back(-1);
neighborState_.emplace_back();
@@ -905,6 +1017,13 @@ inline ParticleStorage::iterator ParticleStorage::erase(iterator& it)
hydrodynamicTorque_.pop_back();
oldHydrodynamicForce_.pop_back();
oldHydrodynamicTorque_.pop_back();
+ virtualMass_.pop_back();
+ invMassIncludingVirtual_.pop_back();
+ oldLinearAcceleration_.pop_back();
+ invInertiaBF_.pop_back();
+ virtualInertiaBF_.pop_back();
+ invInertiaBFIncludingVirtual_.pop_back();
+ oldAngularAcceleration_.pop_back();
clusterID_.pop_back();
segmentID_.pop_back();
neighborState_.pop_back();
@@ -956,6 +1075,13 @@ inline void ParticleStorage::reserve(const size_t size)
hydrodynamicTorque_.reserve(size);
oldHydrodynamicForce_.reserve(size);
oldHydrodynamicTorque_.reserve(size);
+ virtualMass_.reserve(size);
+ invMassIncludingVirtual_.reserve(size);
+ oldLinearAcceleration_.reserve(size);
+ invInertiaBF_.reserve(size);
+ virtualInertiaBF_.reserve(size);
+ invInertiaBFIncludingVirtual_.reserve(size);
+ oldAngularAcceleration_.reserve(size);
clusterID_.reserve(size);
segmentID_.reserve(size);
neighborState_.reserve(size);
@@ -992,6 +1118,13 @@ inline void ParticleStorage::clear()
hydrodynamicTorque_.clear();
oldHydrodynamicForce_.clear();
oldHydrodynamicTorque_.clear();
+ virtualMass_.clear();
+ invMassIncludingVirtual_.clear();
+ oldLinearAcceleration_.clear();
+ invInertiaBF_.clear();
+ virtualInertiaBF_.clear();
+ invInertiaBFIncludingVirtual_.clear();
+ oldAngularAcceleration_.clear();
clusterID_.clear();
segmentID_.clear();
neighborState_.clear();
@@ -1029,6 +1162,13 @@ inline size_t ParticleStorage::size() const
//WALBERLA_ASSERT_EQUAL( uid_.size(), hydrodynamicTorque.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), oldHydrodynamicForce.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), oldHydrodynamicTorque.size() );
+ //WALBERLA_ASSERT_EQUAL( uid_.size(), virtualMass.size() );
+ //WALBERLA_ASSERT_EQUAL( uid_.size(), invMassIncludingVirtual.size() );
+ //WALBERLA_ASSERT_EQUAL( uid_.size(), oldLinearAcceleration.size() );
+ //WALBERLA_ASSERT_EQUAL( uid_.size(), invInertiaBF.size() );
+ //WALBERLA_ASSERT_EQUAL( uid_.size(), virtualInertiaBF.size() );
+ //WALBERLA_ASSERT_EQUAL( uid_.size(), invInertiaBFIncludingVirtual.size() );
+ //WALBERLA_ASSERT_EQUAL( uid_.size(), oldAngularAcceleration.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), clusterID.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), segmentID.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), neighborState.size() );
@@ -1459,6 +1599,69 @@ public:
walberla::mesa_pd::Vec3 const & operator()(const data::Particle& p) const {return p.getOldHydrodynamicTorque();}
};
///Predicate that selects a certain property from a Particle
+class SelectParticleVirtualMass
+{
+public:
+ using return_type = walberla::real_t;
+ walberla::real_t& operator()(data::Particle& p) const {return p.getVirtualMassRef();}
+ walberla::real_t& operator()(data::Particle&& p) const {return p.getVirtualMassRef();}
+ walberla::real_t const & operator()(const data::Particle& p) const {return p.getVirtualMass();}
+};
+///Predicate that selects a certain property from a Particle
+class SelectParticleInvMassIncludingVirtual
+{
+public:
+ using return_type = walberla::real_t;
+ walberla::real_t& operator()(data::Particle& p) const {return p.getInvMassIncludingVirtualRef();}
+ walberla::real_t& operator()(data::Particle&& p) const {return p.getInvMassIncludingVirtualRef();}
+ walberla::real_t const & operator()(const data::Particle& p) const {return p.getInvMassIncludingVirtual();}
+};
+///Predicate that selects a certain property from a Particle
+class SelectParticleOldLinearAcceleration
+{
+public:
+ using return_type = walberla::mesa_pd::Vec3;
+ walberla::mesa_pd::Vec3& operator()(data::Particle& p) const {return p.getOldLinearAccelerationRef();}
+ walberla::mesa_pd::Vec3& operator()(data::Particle&& p) const {return p.getOldLinearAccelerationRef();}
+ walberla::mesa_pd::Vec3 const & operator()(const data::Particle& p) const {return p.getOldLinearAcceleration();}
+};
+///Predicate that selects a certain property from a Particle
+class SelectParticleInvInertiaBF
+{
+public:
+ using return_type = walberla::mesa_pd::Mat3;
+ walberla::mesa_pd::Mat3& operator()(data::Particle& p) const {return p.getInvInertiaBFRef();}
+ walberla::mesa_pd::Mat3& operator()(data::Particle&& p) const {return p.getInvInertiaBFRef();}
+ walberla::mesa_pd::Mat3 const & operator()(const data::Particle& p) const {return p.getInvInertiaBF();}
+};
+///Predicate that selects a certain property from a Particle
+class SelectParticleVirtualInertiaBF
+{
+public:
+ using return_type = walberla::mesa_pd::Mat3;
+ walberla::mesa_pd::Mat3& operator()(data::Particle& p) const {return p.getVirtualInertiaBFRef();}
+ walberla::mesa_pd::Mat3& operator()(data::Particle&& p) const {return p.getVirtualInertiaBFRef();}
+ walberla::mesa_pd::Mat3 const & operator()(const data::Particle& p) const {return p.getVirtualInertiaBF();}
+};
+///Predicate that selects a certain property from a Particle
+class SelectParticleInvInertiaBFIncludingVirtual
+{
+public:
+ using return_type = walberla::mesa_pd::Mat3;
+ walberla::mesa_pd::Mat3& operator()(data::Particle& p) const {return p.getInvInertiaBFIncludingVirtualRef();}
+ walberla::mesa_pd::Mat3& operator()(data::Particle&& p) const {return p.getInvInertiaBFIncludingVirtualRef();}
+ walberla::mesa_pd::Mat3 const & operator()(const data::Particle& p) const {return p.getInvInertiaBFIncludingVirtual();}
+};
+///Predicate that selects a certain property from a Particle
+class SelectParticleOldAngularAcceleration
+{
+public:
+ using return_type = walberla::mesa_pd::Vec3;
+ walberla::mesa_pd::Vec3& operator()(data::Particle& p) const {return p.getOldAngularAccelerationRef();}
+ walberla::mesa_pd::Vec3& operator()(data::Particle&& p) const {return p.getOldAngularAccelerationRef();}
+ walberla::mesa_pd::Vec3 const & operator()(const data::Particle& p) const {return p.getOldAngularAcceleration();}
+};
+///Predicate that selects a certain property from a Particle
class SelectParticleClusterID
{
public:
diff --git a/src/mesa_pd/mpi/notifications/ParseMessage.h b/src/mesa_pd/mpi/notifications/ParseMessage.h
index 32ae118a5..ed3124bf2 100644
--- a/src/mesa_pd/mpi/notifications/ParseMessage.h
+++ b/src/mesa_pd/mpi/notifications/ParseMessage.h
@@ -153,6 +153,13 @@ void ParseMessage::operator()(int sender,
pIt->setHydrodynamicTorque(objparam.hydrodynamicTorque_);
pIt->setOldHydrodynamicForce(objparam.oldHydrodynamicForce_);
pIt->setOldHydrodynamicTorque(objparam.oldHydrodynamicTorque_);
+ pIt->setVirtualMass(objparam.virtualMass_);
+ pIt->setInvMassIncludingVirtual(objparam.invMassIncludingVirtual_);
+ pIt->setOldLinearAcceleration(objparam.oldLinearAcceleration_);
+ pIt->setInvInertiaBF(objparam.invInertiaBF_);
+ pIt->setVirtualInertiaBF(objparam.virtualInertiaBF_);
+ pIt->setInvInertiaBFIncludingVirtual(objparam.invInertiaBFIncludingVirtual_);
+ pIt->setOldAngularAcceleration(objparam.oldAngularAcceleration_);
WALBERLA_LOG_DETAIL( "Processed PARTICLE_MIGRATION_NOTIFICATION." );
diff --git a/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h b/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h
index f0944ef6d..f0f23574c 100644
--- a/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h
+++ b/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h
@@ -69,6 +69,13 @@ public:
walberla::mesa_pd::Vec3 hydrodynamicTorque {real_t(0)};
walberla::mesa_pd::Vec3 oldHydrodynamicForce {real_t(0)};
walberla::mesa_pd::Vec3 oldHydrodynamicTorque {real_t(0)};
+ walberla::real_t virtualMass {real_t(0)};
+ walberla::real_t invMassIncludingVirtual {real_t(0)};
+ walberla::mesa_pd::Vec3 oldLinearAcceleration {real_t(0)};
+ walberla::mesa_pd::Mat3 invInertiaBF {real_t(0)};
+ walberla::mesa_pd::Mat3 virtualInertiaBF {real_t(0)};
+ walberla::mesa_pd::Mat3 invInertiaBFIncludingVirtual {real_t(0)};
+ walberla::mesa_pd::Vec3 oldAngularAcceleration {real_t(0)};
int64_t clusterID {-1};
int64_t segmentID {-1};
};
@@ -103,6 +110,13 @@ inline data::ParticleStorage::iterator createNewParticle(data::ParticleStorage&
pIt->setHydrodynamicTorque(data.hydrodynamicTorque);
pIt->setOldHydrodynamicForce(data.oldHydrodynamicForce);
pIt->setOldHydrodynamicTorque(data.oldHydrodynamicTorque);
+ pIt->setVirtualMass(data.virtualMass);
+ pIt->setInvMassIncludingVirtual(data.invMassIncludingVirtual);
+ pIt->setOldLinearAcceleration(data.oldLinearAcceleration);
+ pIt->setInvInertiaBF(data.invInertiaBF);
+ pIt->setVirtualInertiaBF(data.virtualInertiaBF);
+ pIt->setInvInertiaBFIncludingVirtual(data.invInertiaBFIncludingVirtual);
+ pIt->setOldAngularAcceleration(data.oldAngularAcceleration);
pIt->setClusterID(data.clusterID);
pIt->setSegmentID(data.segmentID);
return pIt;
@@ -152,6 +166,13 @@ mpi::GenericSendBuffer<T,G>& operator<<( mpi::GenericSendBuffer<T,G> & buf, cons
buf << obj.particle_.getHydrodynamicTorque();
buf << obj.particle_.getOldHydrodynamicForce();
buf << obj.particle_.getOldHydrodynamicTorque();
+ buf << obj.particle_.getVirtualMass();
+ buf << obj.particle_.getInvMassIncludingVirtual();
+ buf << obj.particle_.getOldLinearAcceleration();
+ buf << obj.particle_.getInvInertiaBF();
+ buf << obj.particle_.getVirtualInertiaBF();
+ buf << obj.particle_.getInvInertiaBFIncludingVirtual();
+ buf << obj.particle_.getOldAngularAcceleration();
buf << obj.particle_.getClusterID();
buf << obj.particle_.getSegmentID();
return buf;
@@ -182,6 +203,13 @@ mpi::GenericRecvBuffer<T>& operator>>( mpi::GenericRecvBuffer<T> & buf, mesa_pd:
buf >> objparam.hydrodynamicTorque;
buf >> objparam.oldHydrodynamicForce;
buf >> objparam.oldHydrodynamicTorque;
+ buf >> objparam.virtualMass;
+ buf >> objparam.invMassIncludingVirtual;
+ buf >> objparam.oldLinearAcceleration;
+ buf >> objparam.invInertiaBF;
+ buf >> objparam.virtualInertiaBF;
+ buf >> objparam.invInertiaBFIncludingVirtual;
+ buf >> objparam.oldAngularAcceleration;
buf >> objparam.clusterID;
buf >> objparam.segmentID;
return buf;
diff --git a/src/mesa_pd/mpi/notifications/ParticleMigrationNotification.h b/src/mesa_pd/mpi/notifications/ParticleMigrationNotification.h
index 847368ed4..d7cf2b54d 100644
--- a/src/mesa_pd/mpi/notifications/ParticleMigrationNotification.h
+++ b/src/mesa_pd/mpi/notifications/ParticleMigrationNotification.h
@@ -52,6 +52,13 @@ public:
walberla::mesa_pd::Vec3 hydrodynamicTorque_ {real_t(0)};
walberla::mesa_pd::Vec3 oldHydrodynamicForce_ {real_t(0)};
walberla::mesa_pd::Vec3 oldHydrodynamicTorque_ {real_t(0)};
+ walberla::real_t virtualMass_ {real_t(0)};
+ walberla::real_t invMassIncludingVirtual_ {real_t(0)};
+ walberla::mesa_pd::Vec3 oldLinearAcceleration_ {real_t(0)};
+ walberla::mesa_pd::Mat3 invInertiaBF_ {real_t(0)};
+ walberla::mesa_pd::Mat3 virtualInertiaBF_ {real_t(0)};
+ walberla::mesa_pd::Mat3 invInertiaBFIncludingVirtual_ {real_t(0)};
+ walberla::mesa_pd::Vec3 oldAngularAcceleration_ {real_t(0)};
};
inline explicit ParticleMigrationNotification( const data::Particle& particle ) : particle_(particle) {}
@@ -89,6 +96,13 @@ mpi::GenericSendBuffer<T,G>& operator<<( mpi::GenericSendBuffer<T,G> & buf, cons
buf << obj.particle_.getHydrodynamicTorque();
buf << obj.particle_.getOldHydrodynamicForce();
buf << obj.particle_.getOldHydrodynamicTorque();
+ buf << obj.particle_.getVirtualMass();
+ buf << obj.particle_.getInvMassIncludingVirtual();
+ buf << obj.particle_.getOldLinearAcceleration();
+ buf << obj.particle_.getInvInertiaBF();
+ buf << obj.particle_.getVirtualInertiaBF();
+ buf << obj.particle_.getInvInertiaBFIncludingVirtual();
+ buf << obj.particle_.getOldAngularAcceleration();
return buf;
}
@@ -104,6 +118,13 @@ mpi::GenericRecvBuffer<T>& operator>>( mpi::GenericRecvBuffer<T> & buf, mesa_pd:
buf >> objparam.hydrodynamicTorque_;
buf >> objparam.oldHydrodynamicForce_;
buf >> objparam.oldHydrodynamicTorque_;
+ buf >> objparam.virtualMass_;
+ buf >> objparam.invMassIncludingVirtual_;
+ buf >> objparam.oldLinearAcceleration_;
+ buf >> objparam.invInertiaBF_;
+ buf >> objparam.virtualInertiaBF_;
+ buf >> objparam.invInertiaBFIncludingVirtual_;
+ buf >> objparam.oldAngularAcceleration_;
return buf;
}
diff --git a/tests/lbm_mesapd_coupling/CMakeLists.txt b/tests/lbm_mesapd_coupling/CMakeLists.txt
index 90cf13471..eeb208530 100644
--- a/tests/lbm_mesapd_coupling/CMakeLists.txt
+++ b/tests/lbm_mesapd_coupling/CMakeLists.txt
@@ -35,6 +35,9 @@ waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_MEM_UpdateParticleMapping PROCES
waLBerla_compile_test( NAME LBM_MESAPD_COUPLING_UTIL_LubricationCorrection FILES utility/LubricationCorrection.cpp DEPENDS mesa_pd lbm_mesapd_coupling )
waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_UTIL_LubricationCorrection PROCESSES 1 )
+waLBerla_compile_test( NAME LBM_MESAPD_COUPLING_UTIL_VirtualMass FILES utility/VirtualMass.cpp DEPENDS mesa_pd lbm_mesapd_coupling )
+waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_UTIL_VirtualMass PROCESSES 1 )
+
waLBerla_compile_test( NAME LBM_MESAPD_COUPLING_UTIL_InspectionProbe FILES utility/InspectionProbe.cpp DEPENDS core mesa_pd lbm lbm_mesapd_coupling domain_decomposition field )
waLBerla_execute_test( NAME LBM_MESAPD_COUPLING_UTIL_InspectionProbe PROCESSES 1 )
diff --git a/tests/lbm_mesapd_coupling/utility/VirtualMass.cpp b/tests/lbm_mesapd_coupling/utility/VirtualMass.cpp
new file mode 100644
index 000000000..4337d0be7
--- /dev/null
+++ b/tests/lbm_mesapd_coupling/utility/VirtualMass.cpp
@@ -0,0 +1,117 @@
+//======================================================================================================================
+//
+// 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 VirtualMass.cpp
+//! \ingroup lbm_mesapd_coupling
+//! \author Lukas Werner <lks.werner@fau.de>
+//! \brief Evaluates the algorithms behind the virtual mass kernels by comparing with reference values.
+//
+//======================================================================================================================
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/math/all.h"
+
+#include "lbm_mesapd_coupling/utility/virtualmass/AddVirtualForceTorqueKernel.h"
+#include "lbm_mesapd_coupling/utility/virtualmass/InitializeVirtualMassKernel.h"
+#include "lbm_mesapd_coupling/utility/virtualmass/ParticleAccessorWithShapeVirtualMassWrapper.h"
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/ShapeStorage.h"
+
+namespace virtual_mass_test
+{
+using namespace walberla;
+
+int main( int argc, char **argv )
+{
+ debug::enterTestMode();
+
+ mpi::Environment env( argc, argv );
+
+ auto ps = std::make_shared<mesa_pd::data::ParticleStorage>(1);
+ auto ss = std::make_shared<mesa_pd::data::ShapeStorage>();
+ using ParticleAccessor = mesa_pd::data::ParticleAccessorWithShape;
+ ParticleAccessor accessor(ps, ss);
+
+ real_t sphereRadius = real_t(1);
+ real_t sphereDensity = real_t(1);
+ auto sphereShape = ss->create<mesa_pd::data::Sphere>(sphereRadius);
+ ss->shapes[sphereShape]->updateMassAndInertia(sphereDensity);
+
+ mesa_pd::data::Particle&& p1 = *ps->create();
+ p1.setShapeID(sphereShape);
+ auto idx = p1.getIdx();
+
+ using VirtualMass_ParticleAccessor_T = lbm_mesapd_coupling::ParticleAccessorWithShapeVirtualMassWrapper<ParticleAccessor>;
+ auto virtualMassAccessor = walberla::make_shared<VirtualMass_ParticleAccessor_T>(ps, ss);
+
+ auto C_v = real_t(0.5);
+ auto C_v_omega = real_t(0.5);
+ auto fluidDensity = real_t(0.5);
+
+ lbm_mesapd_coupling::InitializeVirtualMassKernel virtualMassInit;
+ virtualMassInit(idx, accessor, C_v, C_v_omega, fluidDensity);
+
+ const real_t sphereVirtualMass = C_v * fluidDensity * accessor.getVolume(idx);
+
+ WALBERLA_CHECK_FLOAT_EQUAL(sphereVirtualMass, accessor.getVirtualMass(idx));
+ WALBERLA_CHECK_FLOAT_EQUAL(real_t(1.) / (accessor.getMass(idx) + sphereVirtualMass),
+ virtualMassAccessor->getInvMass(idx));
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getMass(idx) + sphereVirtualMass,
+ virtualMassAccessor->getMass(idx));
+
+ const real_t sphereAngularVirtualMass = C_v_omega * fluidDensity * accessor.getVolume(idx);
+ const mesa_pd::Mat3 sphereVirtualInertiaBF = accessor.getInertiaBF(idx) * accessor.getInvMass(idx) * sphereAngularVirtualMass;
+
+ WALBERLA_CHECK_FLOAT_EQUAL(sphereVirtualInertiaBF, accessor.getVirtualInertiaBF(idx));
+ WALBERLA_CHECK_FLOAT_EQUAL((accessor.getInertiaBF(idx) + sphereVirtualInertiaBF).getInverse(),
+ virtualMassAccessor->getInvInertiaBF(idx));
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getInertiaBF(idx) + sphereVirtualInertiaBF,
+ virtualMassAccessor->getInertiaBF(idx));
+
+ lbm_mesapd_coupling::AddVirtualForceTorqueKernel addVirtualForceTorque(ps);
+
+ accessor.setForce(idx, Vector3(real_t(1)));
+ auto oldForce = accessor.getForce(idx);
+ addVirtualForceTorque(idx, accessor);
+ WALBERLA_CHECK_FLOAT_EQUAL(oldForce, accessor.getForce(idx));
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getOldLinearAcceleration(idx),
+ Vector3<real_t>(real_t(1) / (accessor.getMass(idx) + accessor.getVirtualMass(idx))));
+
+ Vector3 oldAcceleration(real_t(2));
+ accessor.setOldLinearAcceleration(idx, oldAcceleration);
+ accessor.setOldAngularAcceleration(idx, oldAcceleration);
+ oldForce = Vector3(real_t(2));
+ accessor.setForce(idx, oldForce);
+ Vector3 oldTorque(real_t(2));
+ accessor.setTorque(idx, oldTorque);
+ addVirtualForceTorque(idx, accessor);
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getForce(idx), oldForce + sphereVirtualMass*oldAcceleration);
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getForce(idx),
+ accessor.getOldLinearAcceleration(idx)*(accessor.getMass(idx)+accessor.getVirtualMass(idx)));
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getTorque(idx), oldTorque + sphereVirtualInertiaBF*oldAcceleration);
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getTorque(idx),
+ (accessor.getInertiaBF(idx)+accessor.getVirtualInertiaBF(idx))*accessor.getOldAngularAcceleration(idx));
+
+ return 0;
+}
+
+} // namespace virtual_mass_test
+
+int main( int argc, char **argv ){
+ virtual_mass_test::main(argc, argv);
+}
--
GitLab