From 95121137f5e17264b084cdbc247107d6834c3dcd Mon Sep 17 00:00:00 2001
From: Sebastian Eibl <sebastian.eibl@fau.de>
Date: Mon, 25 May 2020 18:17:57 +0200
Subject: [PATCH] [ADD] heat conduction showcase
---
apps/showcases/CMakeLists.txt | 1 +
apps/showcases/HeatConduction/CMakeLists.txt | 5 +
.../HeatConduction/ContactDetection.h | 170 ++++++++++
.../HeatConduction/HeatConduction.cfg | 15 +
.../HeatConduction/HeatConduction.cpp | 291 ++++++++++++++++++
.../HeatConduction/ThermalExpansion.h | 51 +++
python/mesa_pd.py | 2 +
src/mesa_pd/data/ParticleAccessor.h | 9 +
src/mesa_pd/data/ParticleStorage.h | 29 ++
src/mesa_pd/mpi/notifications/ParseMessage.h | 1 +
.../notifications/ParticleCopyNotification.h | 4 +
.../ParticleGhostCopyNotification.h | 4 +
.../ParticleUpdateNotification.h | 3 +
13 files changed, 585 insertions(+)
create mode 100644 apps/showcases/HeatConduction/CMakeLists.txt
create mode 100644 apps/showcases/HeatConduction/ContactDetection.h
create mode 100644 apps/showcases/HeatConduction/HeatConduction.cfg
create mode 100644 apps/showcases/HeatConduction/HeatConduction.cpp
create mode 100644 apps/showcases/HeatConduction/ThermalExpansion.h
diff --git a/apps/showcases/CMakeLists.txt b/apps/showcases/CMakeLists.txt
index 536f56ff6..266582c22 100644
--- a/apps/showcases/CMakeLists.txt
+++ b/apps/showcases/CMakeLists.txt
@@ -1,3 +1,4 @@
add_subdirectory( BidisperseFluidizedBed )
add_subdirectory( CombinedResolvedUnresolved )
+add_subdirectory( HeatConduction )
add_subdirectory( Mixer )
diff --git a/apps/showcases/HeatConduction/CMakeLists.txt b/apps/showcases/HeatConduction/CMakeLists.txt
new file mode 100644
index 000000000..f9465cda0
--- /dev/null
+++ b/apps/showcases/HeatConduction/CMakeLists.txt
@@ -0,0 +1,5 @@
+waLBerla_link_files_to_builddir( *.cfg )
+
+waLBerla_add_executable( NAME HeatConduction
+ FILES HeatConduction.cpp
+ DEPENDS core mesa_pd vtk )
diff --git a/apps/showcases/HeatConduction/ContactDetection.h b/apps/showcases/HeatConduction/ContactDetection.h
new file mode 100644
index 000000000..f1c0fd5d1
--- /dev/null
+++ b/apps/showcases/HeatConduction/ContactDetection.h
@@ -0,0 +1,170 @@
+//======================================================================================================================
+//
+// 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 ContactDetection.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/collision_detection/AnalyticContactDetection.h>
+#include <mesa_pd/collision_detection/AnalyticCollisionFunctions.h>
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/shape/BaseShape.h>
+#include <mesa_pd/data/shape/HalfSpace.h>
+#include <mesa_pd/data/shape/Sphere.h>
+
+#include <core/logging/Logging.h>
+
+namespace walberla {
+namespace mesa_pd {
+
+class ContactDetection
+{
+public:
+ size_t& getIdx1() {return idx1_;}
+ size_t& getIdx2() {return idx2_;}
+ Vec3& getContactPoint() {return contactPoint_;}
+ Vec3& getContactNormal() {return contactNormal_;}
+ real_t& getPenetrationDepth() {return penetrationDepth_;}
+
+ const size_t& getIdx1() const {return idx1_;}
+ const size_t& getIdx2() const {return idx2_;}
+ const Vec3& getContactPoint() const {return contactPoint_;}
+ const Vec3& getContactNormal() const {return contactNormal_;}
+ const real_t& getPenetrationDepth() const {return penetrationDepth_;}
+
+ real_t& getContactThreshold() {return contactThreshold_;}
+
+ template <typename GEO1_T, typename GEO2_T, typename Accessor>
+ bool operator()( const size_t idx1,
+ const size_t idx2,
+ const GEO1_T& geo1,
+ const GEO2_T& geo2,
+ Accessor& ac);
+
+ template <typename Accessor>
+ bool operator()( const size_t idx1,
+ const size_t idx2,
+ const data::Sphere& geo1,
+ const data::Sphere& geo2,
+ Accessor& ac);
+
+ template <typename Accessor>
+ bool operator()( const size_t idx1,
+ const size_t idx2,
+ const data::Sphere& s,
+ const data::HalfSpace& p,
+ Accessor& ac );
+
+ template <typename Accessor>
+ bool operator()( const size_t idx1,
+ const size_t idx2,
+ const data::HalfSpace& p,
+ const data::Sphere& s,
+ Accessor& ac);
+
+private:
+ size_t idx1_;
+ size_t idx2_;
+ Vec3 contactPoint_;
+ Vec3 contactNormal_;
+ real_t penetrationDepth_;
+
+ real_t contactThreshold_ = real_t(0.0);
+};
+
+template <typename GEO1_T, typename GEO2_T, typename Accessor>
+inline bool ContactDetection::operator()( const size_t /*idx1*/,
+ const size_t /*idx2*/,
+ const GEO1_T& /*geo1*/,
+ const GEO2_T& /*geo2*/,
+ Accessor& /*ac*/)
+{
+ WALBERLA_ABORT("Collision not implemented!")
+}
+
+template <typename Accessor>
+inline bool ContactDetection::operator()( const size_t idx1,
+ const size_t idx2,
+ const data::Sphere& geo1,
+ const data::Sphere& geo2,
+ Accessor& ac)
+{
+ using namespace collision_detection::analytic;
+ WALBERLA_ASSERT_UNEQUAL(idx1, idx2, "colliding with itself!");
+
+ //ensure collision order idx2 has to be larger than idx1
+ if (ac.getUid(idx2) < ac.getUid(idx1))
+ return operator()(idx2, idx1, geo2, geo1, ac);
+
+ getIdx1() = idx1;
+ getIdx2() = idx2;
+ return detectSphereSphereCollision(ac.getPosition(getIdx1()),
+ ac.getRadius(getIdx1()),
+ ac.getPosition(getIdx2()),
+ ac.getRadius(getIdx2()),
+ getContactPoint(),
+ getContactNormal(),
+ getPenetrationDepth(),
+ getContactThreshold());
+}
+
+template <typename Accessor>
+inline bool ContactDetection::operator()( const size_t idx1,
+ const size_t idx2,
+ const data::Sphere& s,
+ const data::HalfSpace& p,
+ Accessor& ac )
+{
+ using namespace collision_detection::analytic;
+ WALBERLA_ASSERT_UNEQUAL(idx1, idx2, "colliding with itself!");
+
+ getIdx1() = idx1;
+ getIdx2() = idx2;
+ return detectSphereHalfSpaceCollision(ac.getPosition(getIdx1()),
+ ac.getRadius(getIdx1()),
+ ac.getPosition(getIdx2()),
+ p.getNormal(),
+ getContactPoint(),
+ getContactNormal(),
+ getPenetrationDepth(),
+ getContactThreshold());
+}
+
+template <typename Accessor>
+inline bool ContactDetection::operator()( const size_t idx1,
+ const size_t idx2,
+ const data::HalfSpace& p,
+ const data::Sphere& s,
+ Accessor& ac)
+{
+ return operator()(idx2, idx1, s, p, ac);
+}
+
+inline
+std::ostream& operator<<( std::ostream& os, const ContactDetection& ac )
+{
+ os << "idx1: " << ac.getIdx1() << "\n" <<
+ "idx2: " << ac.getIdx2() << "\n" <<
+ "contact point: " << ac.getContactPoint() << "\n" <<
+ "contact normal: " << ac.getContactNormal() << "\n" <<
+ "penetration depth: " << ac.getPenetrationDepth() << std::endl;
+ return os;
+}
+
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/apps/showcases/HeatConduction/HeatConduction.cfg b/apps/showcases/HeatConduction/HeatConduction.cfg
new file mode 100644
index 000000000..7d1d72982
--- /dev/null
+++ b/apps/showcases/HeatConduction/HeatConduction.cfg
@@ -0,0 +1,15 @@
+HeatConduction
+{
+ simulationCorner < 0, 0, 0 >;
+ simulationDomain < 0.3, 0.3, 0.3 >;
+ blocks < 3,3,1 >;
+ isPeriodic < 1, 1, 0 >;
+
+ radius 0.004;
+ spacing 0.010;
+ vMax 0.00000001;
+
+ dt 0.000001;
+ simulationSteps 30000;
+ visSpacing 100;
+}
diff --git a/apps/showcases/HeatConduction/HeatConduction.cpp b/apps/showcases/HeatConduction/HeatConduction.cpp
new file mode 100644
index 000000000..72bf5755d
--- /dev/null
+++ b/apps/showcases/HeatConduction/HeatConduction.cpp
@@ -0,0 +1,291 @@
+//======================================================================================================================
+//
+// 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 HeatConduction.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include <mesa_pd/vtk/ParticleVtkOutput.h>
+
+#include <mesa_pd/collision_detection/BroadPhase.h>
+
+#include <mesa_pd/data/LinkedCells.h>
+#include <mesa_pd/data/ParticleAccessorWithShape.h>
+#include <mesa_pd/data/ParticleStorage.h>
+#include <mesa_pd/data/ShapeStorage.h>
+
+#include <mesa_pd/domain/BlockForestDomain.h>
+
+#include <mesa_pd/kernel/DoubleCast.h>
+#include <mesa_pd/kernel/ExplicitEulerWithShape.h>
+#include <mesa_pd/kernel/HeatConduction.h>
+#include <mesa_pd/kernel/InsertParticleIntoLinkedCells.h>
+#include <mesa_pd/kernel/SpringDashpot.h>
+#include <mesa_pd/kernel/TemperatureIntegration.h>
+
+#include <mesa_pd/mpi/ContactFilter.h>
+#include <mesa_pd/mpi/SyncNextNeighbors.h>
+#include <mesa_pd/mpi/ReduceProperty.h>
+
+#include <mesa_pd/mpi/notifications/ForceTorqueNotification.h>
+#include <mesa_pd/mpi/notifications/HeatFluxNotification.h>
+
+#include "ContactDetection.h"
+#include "ThermalExpansion.h"
+
+#include <blockforest/Initialization.h>
+#include <core/Abort.h>
+#include <core/Environment.h>
+#include <core/math/Random.h>
+#include <core/mpi/Reduce.h>
+#include <core/grid_generator/SCIterator.h>
+#include <core/logging/Logging.h>
+#include <core/OpenMP.h>
+#include <core/timing/Timer.h>
+#include <core/timing/TimingPool.h>
+#include <core/waLBerlaBuildInfo.h>
+#include <sqlite/SQLite.h>
+#include <vtk/VTKOutput.h>
+
+#include <functional>
+#include <memory>
+#include <string>
+#include <mesa_pd/kernel/ParticleSelector.h>
+
+namespace walberla {
+namespace mesa_pd {
+
+class CustomParticleAccessor : public data::ParticleAccessorWithShape
+{
+public:
+ CustomParticleAccessor(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss, const real_t radius273K)
+ : ParticleAccessorWithShape(ps, ss)
+ , radius273K_(radius273K)
+ {}
+
+ real_t getRadius273K(const size_t /*p_idx*/) const {return radius273K_;}
+private:
+ real_t radius273K_;
+};
+
+auto createPlane( const std::shared_ptr<data::ParticleStorage>& ps,
+ const std::shared_ptr<data::ShapeStorage>& ss,
+ const Vec3& pos,
+ const Vec3& normal)
+{
+ auto p0 = ps->create(true);
+ p0->setPosition( pos );
+ p0->setShapeID( ss->create<data::HalfSpace>( normal ) );
+ p0->setOwner( walberla::mpi::MPIManager::instance()->rank() );
+ p0->setType( 0 );
+ p0->setTemperature( 273 );
+ data::particle_flags::set(p0->getFlagsRef(), data::particle_flags::INFINITE);
+ data::particle_flags::set(p0->getFlagsRef(), data::particle_flags::FIXED);
+ data::particle_flags::set(p0->getFlagsRef(), data::particle_flags::NON_COMMUNICATING);
+ return p0;
+}
+
+int main( int argc, char ** argv )
+{
+ using namespace walberla::timing;
+
+ Environment env(argc, argv);
+ walberla::mpi::MPIManager::instance()->useWorldComm();
+
+ WALBERLA_LOG_INFO_ON_ROOT( "config file: " << argv[1] );
+ WALBERLA_LOG_INFO_ON_ROOT( "waLBerla Revision: " << WALBERLA_GIT_SHA1 );
+
+ math::seedRandomGenerator( static_cast<unsigned int>(1337 * walberla::mpi::MPIManager::instance()->worldRank()) );
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** READING CONFIG FILE ***");
+ auto cfg = env.config();
+ if (cfg == nullptr) WALBERLA_ABORT("No config specified!");
+ const Config::BlockHandle mainConf = cfg->getBlock( "HeatConduction" );
+
+ const real_t spacing = mainConf.getParameter<real_t>("spacing", real_t(1.0) );
+ WALBERLA_LOG_INFO_ON_ROOT("spacing: " << spacing);
+
+ const real_t radius = mainConf.getParameter<real_t>("radius", real_t(0.5) );
+ WALBERLA_LOG_INFO_ON_ROOT("radius: " << radius);
+
+ const real_t vMax = mainConf.getParameter<real_t>("vMax", real_t(0.5) );
+ WALBERLA_LOG_INFO_ON_ROOT("vMax: " << vMax);
+
+ int64_t numOuterIterations = mainConf.getParameter<int64_t>("numOuterIterations", 10 );
+ WALBERLA_LOG_INFO_ON_ROOT("numOuterIterations: " << numOuterIterations);
+
+ int64_t simulationSteps = mainConf.getParameter<int64_t>("simulationSteps", 10 );
+ WALBERLA_LOG_INFO_ON_ROOT("simulationSteps: " << simulationSteps);
+
+ real_t dt = mainConf.getParameter<real_t>("dt", real_c(0.01) );
+ WALBERLA_LOG_INFO_ON_ROOT("dt: " << dt);
+
+ const int visSpacing = mainConf.getParameter<int>("visSpacing", 1000 );
+ WALBERLA_LOG_INFO_ON_ROOT("visSpacing: " << visSpacing);
+ const std::string path = mainConf.getParameter<std::string>("path", "vtk_out" );
+ WALBERLA_LOG_INFO_ON_ROOT("path: " << path);
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** BLOCKFOREST ***");
+ // create forest
+ auto forest = blockforest::createBlockForestFromConfig( mainConf );
+ if (!forest)
+ {
+ WALBERLA_LOG_INFO_ON_ROOT( "No BlockForest created ... exiting!");
+ return EXIT_SUCCESS;
+ }
+ domain::BlockForestDomain domain(forest);
+
+ auto simulationDomain = forest->getDomain();
+ auto localDomain = forest->begin()->getAABB();
+ for (auto& blk : *forest)
+ {
+ localDomain.merge(blk.getAABB());
+ }
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** SETUP - START ***");
+
+ //init data structures
+ auto ps = std::make_shared<data::ParticleStorage>(100);
+ auto ss = std::make_shared<data::ShapeStorage>();
+ CustomParticleAccessor accessor(ps, ss, radius);
+ data::LinkedCells lc(localDomain.getExtended(spacing), spacing+spacing );
+
+ auto smallSphere = ss->create<data::Sphere>( radius );
+ ss->shapes[smallSphere]->updateMassAndInertia(real_t(2707));
+ for (auto& iBlk : *forest)
+ {
+ for (auto pt : grid_generator::SCGrid(iBlk.getAABB(), Vector3<real_t>(spacing, spacing, spacing) * real_c(0.5), spacing))
+ {
+ WALBERLA_CHECK(iBlk.getAABB().contains(pt));
+
+ auto p = ps->create();
+ p->getPositionRef() = pt;
+ p->getInteractionRadiusRef() = radius;
+ p->setLinearVelocity( Vec3(math::realRandom(-vMax, vMax),math::realRandom(-vMax, vMax),math::realRandom(-vMax, vMax)) );
+ p->getShapeIDRef() = smallSphere;
+ p->getOwnerRef() = walberla::mpi::MPIManager::instance()->rank();
+ p->getTypeRef() = 0;
+ p->setTemperature( 273 );
+ p->setRadius(radius);
+ }
+ }
+ int64_t numParticles = int64_c(ps->size());
+ walberla::mpi::reduceInplace(numParticles, walberla::mpi::SUM);
+ WALBERLA_LOG_INFO_ON_ROOT("#particles created: " << numParticles);
+
+ auto plane = createPlane(ps, ss, simulationDomain.minCorner(), Vec3(0,0,1));
+ createPlane(ps, ss, simulationDomain.maxCorner(), Vec3(0,0,-1));
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** SETUP - END ***");
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** VTK ***");
+ auto vtkDomainOutput = walberla::vtk::createVTKOutput_DomainDecomposition( forest, "domain_partitioning", 1, "vtk", "simulation_step" );
+ auto vtkOutput = make_shared<mesa_pd::vtk::ParticleVtkOutput>(ps) ;
+ auto vtkWriter = walberla::vtk::createVTKOutput_PointData(vtkOutput, "particles", 1, "vtk", "simulation_step", false, false);
+ vtkOutput->addOutput<data::SelectParticleOwner>("owner");
+ vtkOutput->addOutput<data::SelectParticleTemperature>("temperature");
+ vtkOutput->addOutput<data::SelectParticleRadius>("radius");
+ vtkOutput->addOutput<data::SelectParticleLinearVelocity>("linVel");
+ vtkOutput->setParticleSelector([smallSphere](const data::ParticleStorage::iterator& pIt){ return pIt->getShapeID() == smallSphere;});
+ vtkDomainOutput->write();
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - START ***");
+ // Init kernels
+ kernel::ExplicitEuler explicitEuler( dt );
+ kernel::InsertParticleIntoLinkedCells ipilc;
+ kernel::HeatConduction heatConduction(1);
+ heatConduction.setConductance(0, 0, real_t(1));
+ kernel::SpringDashpot dem(1);
+ dem.setStiffness(0, 0, real_t(8.11e6));
+ dem.setDampingN (0, 0, real_t(6.86e1));
+ dem.setDampingT (0, 0, real_t(6.86e1));
+ dem.setFriction (0, 0, real_t(1.2));
+ kernel::TemperatureIntegration temperatureIntegration(dt, 1);
+ temperatureIntegration.setInvSpecificHeat(0, real_t(0.05) / dt / ss->shapes[smallSphere]->getInvMass());
+ kernel::ThermalExpansion thermalExpansion(1);
+ thermalExpansion.setLinearExpansionCoefficient(0, real_t(0.0005));
+
+ ContactDetection acd;
+ kernel::DoubleCast double_cast;
+ mpi::ContactFilter contact_filter;
+ mpi::ReduceProperty RP;
+ mpi::SyncNextNeighbors SNN;
+
+ // initial sync
+ SNN(*ps, domain);
+
+ for (int64_t i=0; i < simulationSteps; ++i)
+ {
+ if (i % visSpacing == 0)
+ {
+ vtkWriter->write();
+ }
+
+ if (i > visSpacing * 100)
+ {
+ auto rad = (real_t(i) - real_t(10000)) * real_t(5e-5) * math::pi;
+ plane->setTemperature(real_t(273) + real_t(300) * std::sin(rad));
+ }
+
+ ps->forEachParticle(false, kernel::SelectMaster(), accessor, [&](const size_t idx, auto& ac){ac.setForce(idx, Vec3(0,0,real_t(-9.81) ) );}, accessor);
+
+ lc.clear();
+ ps->forEachParticle(true, kernel::SelectAll(), accessor, ipilc, accessor, lc);
+
+ lc.forEachParticlePairHalf(true, kernel::SelectAll(), accessor,
+ [&](const size_t idx1, const size_t idx2, auto& ac)
+ {
+ if (collision_detection::isInInteractionDistance(idx1, idx2, ac))
+ {
+ if (double_cast(idx1, idx2, ac, acd, ac ))
+ {
+ if (contact_filter(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), domain))
+ {
+ dem(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), acd.getContactNormal(), acd.getPenetrationDepth());
+ heatConduction(acd.getIdx1(), acd.getIdx2(), ac);
+ }
+ }
+ }
+ },
+ accessor );
+
+ RP.operator()<ForceTorqueNotification>(*ps);
+
+ RP.operator()<HeatFluxNotification>(*ps);
+
+ ps->forEachParticle(true, kernel::SelectMaster(), accessor, explicitEuler, accessor);
+
+ ps->forEachParticle(true, kernel::SelectMaster(), accessor, temperatureIntegration, accessor);
+
+ if( i % 10 == 0 )
+ {
+ ps->forEachParticle(true, kernel::SelectMaster(), accessor, thermalExpansion, accessor);
+ }
+
+ SNN(*ps, domain);
+ }
+ WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - END ***");
+
+ return EXIT_SUCCESS;
+}
+
+} // namespace mesa_pd
+} // namespace walberla
+
+int main( int argc, char* argv[] )
+{
+ return walberla::mesa_pd::main( argc, argv );
+}
diff --git a/apps/showcases/HeatConduction/ThermalExpansion.h b/apps/showcases/HeatConduction/ThermalExpansion.h
new file mode 100644
index 000000000..e553c2692
--- /dev/null
+++ b/apps/showcases/HeatConduction/ThermalExpansion.h
@@ -0,0 +1,51 @@
+//======================================================================================================================
+//
+// 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 ThermalExpansion.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/IAccessor.h>
+
+namespace walberla {
+namespace mesa_pd {
+namespace kernel {
+
+class ThermalExpansion
+{
+public:
+ ThermalExpansion() = default;
+
+ template <typename Accessor>
+ void operator()(const size_t i, Accessor& ac) const;
+};
+
+template <typename Accessor>
+inline void ThermalExpansion::operator()(const size_t idx,
+ Accessor& ac) const
+{
+ static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
+
+ ac.setRadius(idx, real_t(0.004) + (ac.getTemperature(idx)-real_t(273)) * real_t(0.002) * real_t(0.001));
+ ac.setInteractionRadius(idx, ac.getRadius(idx));
+}
+
+} //namespace kernel
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/python/mesa_pd.py b/python/mesa_pd.py
index 785d9a4ff..dd53b2f4c 100755
--- a/python/mesa_pd.py
+++ b/python/mesa_pd.py
@@ -28,6 +28,8 @@ if __name__ == '__main__':
ps.add_property("torque", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
ps.add_property("oldTorque", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ON_OWNERSHIP_CHANGE")
+ ps.add_property("radius", "walberla::real_t", defValue="real_t(0)", syncMode="ALWAYS")
+
ps.add_include("blockforest/BlockForest.h")
ps.add_property("currentBlock", "blockforest::BlockID", defValue="", syncMode="NEVER")
diff --git a/src/mesa_pd/data/ParticleAccessor.h b/src/mesa_pd/data/ParticleAccessor.h
index 59e575330..753546900 100644
--- a/src/mesa_pd/data/ParticleAccessor.h
+++ b/src/mesa_pd/data/ParticleAccessor.h
@@ -108,6 +108,10 @@ public:
walberla::mesa_pd::Vec3& getOldTorqueRef(const size_t p_idx) {return ps_->getOldTorqueRef(p_idx);}
void setOldTorque(const size_t p_idx, walberla::mesa_pd::Vec3 const & v) { ps_->setOldTorque(p_idx, v);}
+ walberla::real_t const & getRadius(const size_t p_idx) const {return ps_->getRadius(p_idx);}
+ walberla::real_t& getRadiusRef(const size_t p_idx) {return ps_->getRadiusRef(p_idx);}
+ void setRadius(const size_t p_idx, walberla::real_t const & v) { ps_->setRadius(p_idx, v);}
+
blockforest::BlockID const & getCurrentBlock(const size_t p_idx) const {return ps_->getCurrentBlock(p_idx);}
blockforest::BlockID& getCurrentBlockRef(const size_t p_idx) {return ps_->getCurrentBlockRef(p_idx);}
void setCurrentBlock(const size_t p_idx, blockforest::BlockID const & v) { ps_->setCurrentBlock(p_idx, v);}
@@ -269,6 +273,10 @@ public:
void setOldTorque(const size_t /*p_idx*/, walberla::mesa_pd::Vec3 const & v) { oldTorque_ = v;}
walberla::mesa_pd::Vec3& getOldTorqueRef(const size_t /*p_idx*/) {return oldTorque_;}
+ walberla::real_t const & getRadius(const size_t /*p_idx*/) const {return radius_;}
+ void setRadius(const size_t /*p_idx*/, walberla::real_t const & v) { radius_ = v;}
+ walberla::real_t& getRadiusRef(const size_t /*p_idx*/) {return radius_;}
+
blockforest::BlockID const & getCurrentBlock(const size_t /*p_idx*/) const {return currentBlock_;}
void setCurrentBlock(const size_t /*p_idx*/, blockforest::BlockID const & v) { currentBlock_ = v;}
blockforest::BlockID& getCurrentBlockRef(const size_t /*p_idx*/) {return currentBlock_;}
@@ -351,6 +359,7 @@ private:
walberla::mesa_pd::Vec3 angularVelocity_;
walberla::mesa_pd::Vec3 torque_;
walberla::mesa_pd::Vec3 oldTorque_;
+ walberla::real_t radius_;
blockforest::BlockID currentBlock_;
uint_t type_;
int nextParticle_;
diff --git a/src/mesa_pd/data/ParticleStorage.h b/src/mesa_pd/data/ParticleStorage.h
index a22470d67..fd1c33cf8 100644
--- a/src/mesa_pd/data/ParticleStorage.h
+++ b/src/mesa_pd/data/ParticleStorage.h
@@ -86,6 +86,7 @@ public:
using angularVelocity_type = walberla::mesa_pd::Vec3;
using torque_type = walberla::mesa_pd::Vec3;
using oldTorque_type = walberla::mesa_pd::Vec3;
+ using radius_type = walberla::real_t;
using currentBlock_type = blockforest::BlockID;
using type_type = uint_t;
using nextParticle_type = int;
@@ -162,6 +163,10 @@ public:
oldTorque_type& getOldTorqueRef() {return storage_.getOldTorqueRef(i_);}
void setOldTorque(oldTorque_type const & v) { storage_.setOldTorque(i_, v);}
+ radius_type const & getRadius() const {return storage_.getRadius(i_);}
+ radius_type& getRadiusRef() {return storage_.getRadiusRef(i_);}
+ void setRadius(radius_type const & v) { storage_.setRadius(i_, v);}
+
currentBlock_type const & getCurrentBlock() const {return storage_.getCurrentBlock(i_);}
currentBlock_type& getCurrentBlockRef() {return storage_.getCurrentBlockRef(i_);}
void setCurrentBlock(currentBlock_type const & v) { storage_.setCurrentBlock(i_, v);}
@@ -293,6 +298,7 @@ public:
using angularVelocity_type = walberla::mesa_pd::Vec3;
using torque_type = walberla::mesa_pd::Vec3;
using oldTorque_type = walberla::mesa_pd::Vec3;
+ using radius_type = walberla::real_t;
using currentBlock_type = blockforest::BlockID;
using type_type = uint_t;
using nextParticle_type = int;
@@ -369,6 +375,10 @@ public:
oldTorque_type& getOldTorqueRef(const size_t idx) {return oldTorque_[idx];}
void setOldTorque(const size_t idx, oldTorque_type const & v) { oldTorque_[idx] = v; }
+ radius_type const & getRadius(const size_t idx) const {return radius_[idx];}
+ radius_type& getRadiusRef(const size_t idx) {return radius_[idx];}
+ void setRadius(const size_t idx, radius_type const & v) { radius_[idx] = v; }
+
currentBlock_type const & getCurrentBlock(const size_t idx) const {return currentBlock_[idx];}
currentBlock_type& getCurrentBlockRef(const size_t idx) {return currentBlock_[idx];}
void setCurrentBlock(const size_t idx, currentBlock_type const & v) { currentBlock_[idx] = v; }
@@ -531,6 +541,7 @@ public:
std::vector<angularVelocity_type> angularVelocity_ {};
std::vector<torque_type> torque_ {};
std::vector<oldTorque_type> oldTorque_ {};
+ std::vector<radius_type> radius_ {};
std::vector<currentBlock_type> currentBlock_ {};
std::vector<type_type> type_ {};
std::vector<nextParticle_type> nextParticle_ {};
@@ -569,6 +580,7 @@ ParticleStorage::Particle& ParticleStorage::Particle::operator=(const ParticleSt
getAngularVelocityRef() = rhs.getAngularVelocity();
getTorqueRef() = rhs.getTorque();
getOldTorqueRef() = rhs.getOldTorque();
+ getRadiusRef() = rhs.getRadius();
getCurrentBlockRef() = rhs.getCurrentBlock();
getTypeRef() = rhs.getType();
getNextParticleRef() = rhs.getNextParticle();
@@ -604,6 +616,7 @@ ParticleStorage::Particle& ParticleStorage::Particle::operator=(ParticleStorage:
getAngularVelocityRef() = std::move(rhs.getAngularVelocityRef());
getTorqueRef() = std::move(rhs.getTorqueRef());
getOldTorqueRef() = std::move(rhs.getOldTorqueRef());
+ getRadiusRef() = std::move(rhs.getRadiusRef());
getCurrentBlockRef() = std::move(rhs.getCurrentBlockRef());
getTypeRef() = std::move(rhs.getTypeRef());
getNextParticleRef() = std::move(rhs.getNextParticleRef());
@@ -640,6 +653,7 @@ void swap(ParticleStorage::Particle lhs, ParticleStorage::Particle rhs)
std::swap(lhs.getAngularVelocityRef(), rhs.getAngularVelocityRef());
std::swap(lhs.getTorqueRef(), rhs.getTorqueRef());
std::swap(lhs.getOldTorqueRef(), rhs.getOldTorqueRef());
+ std::swap(lhs.getRadiusRef(), rhs.getRadiusRef());
std::swap(lhs.getCurrentBlockRef(), rhs.getCurrentBlockRef());
std::swap(lhs.getTypeRef(), rhs.getTypeRef());
std::swap(lhs.getNextParticleRef(), rhs.getNextParticleRef());
@@ -676,6 +690,7 @@ std::ostream& operator<<( std::ostream& os, const ParticleStorage::Particle& p )
"angularVelocity : " << p.getAngularVelocity() << "\n" <<
"torque : " << p.getTorque() << "\n" <<
"oldTorque : " << p.getOldTorque() << "\n" <<
+ "radius : " << p.getRadius() << "\n" <<
"currentBlock : " << p.getCurrentBlock() << "\n" <<
"type : " << p.getType() << "\n" <<
"nextParticle : " << p.getNextParticle() << "\n" <<
@@ -782,6 +797,7 @@ inline ParticleStorage::iterator ParticleStorage::create(const id_t& uid)
angularVelocity_.emplace_back(real_t(0));
torque_.emplace_back(real_t(0));
oldTorque_.emplace_back(real_t(0));
+ radius_.emplace_back(real_t(0));
currentBlock_.emplace_back();
type_.emplace_back(0);
nextParticle_.emplace_back(-1);
@@ -843,6 +859,7 @@ inline ParticleStorage::iterator ParticleStorage::erase(iterator& it)
angularVelocity_.pop_back();
torque_.pop_back();
oldTorque_.pop_back();
+ radius_.pop_back();
currentBlock_.pop_back();
type_.pop_back();
nextParticle_.pop_back();
@@ -891,6 +908,7 @@ inline void ParticleStorage::reserve(const size_t size)
angularVelocity_.reserve(size);
torque_.reserve(size);
oldTorque_.reserve(size);
+ radius_.reserve(size);
currentBlock_.reserve(size);
type_.reserve(size);
nextParticle_.reserve(size);
@@ -924,6 +942,7 @@ inline void ParticleStorage::clear()
angularVelocity_.clear();
torque_.clear();
oldTorque_.clear();
+ radius_.clear();
currentBlock_.clear();
type_.clear();
nextParticle_.clear();
@@ -958,6 +977,7 @@ inline size_t ParticleStorage::size() const
//WALBERLA_ASSERT_EQUAL( uid_.size(), angularVelocity.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), torque.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), oldTorque.size() );
+ //WALBERLA_ASSERT_EQUAL( uid_.size(), radius.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), currentBlock.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), type.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), nextParticle.size() );
@@ -1291,6 +1311,15 @@ public:
walberla::mesa_pd::Vec3 const & operator()(const data::Particle& p) const {return p.getOldTorque();}
};
///Predicate that selects a certain property from a Particle
+class SelectParticleRadius
+{
+public:
+ using return_type = walberla::real_t;
+ walberla::real_t& operator()(data::Particle& p) const {return p.getRadiusRef();}
+ walberla::real_t& operator()(data::Particle&& p) const {return p.getRadiusRef();}
+ walberla::real_t const & operator()(const data::Particle& p) const {return p.getRadius();}
+};
+///Predicate that selects a certain property from a Particle
class SelectParticleCurrentBlock
{
public:
diff --git a/src/mesa_pd/mpi/notifications/ParseMessage.h b/src/mesa_pd/mpi/notifications/ParseMessage.h
index 63aa0d542..02b36290f 100644
--- a/src/mesa_pd/mpi/notifications/ParseMessage.h
+++ b/src/mesa_pd/mpi/notifications/ParseMessage.h
@@ -117,6 +117,7 @@ void ParseMessage::operator()(int sender,
pIt->setLinearVelocity(objparam.linearVelocity);
pIt->setRotation(objparam.rotation);
pIt->setAngularVelocity(objparam.angularVelocity);
+ pIt->setRadius(objparam.radius);
pIt->setOldContactHistory(objparam.oldContactHistory);
pIt->setTemperature(objparam.temperature);
diff --git a/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h b/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h
index 2ac5383c1..614f8bea3 100644
--- a/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h
+++ b/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h
@@ -61,6 +61,7 @@ public:
walberla::mesa_pd::Rot3 rotation {};
walberla::mesa_pd::Vec3 angularVelocity {real_t(0)};
walberla::mesa_pd::Vec3 oldTorque {real_t(0)};
+ walberla::real_t radius {real_t(0)};
uint_t type {0};
std::map<walberla::id_t, walberla::mesa_pd::data::ContactHistory> oldContactHistory {};
walberla::real_t temperature {real_t(0)};
@@ -92,6 +93,7 @@ inline data::ParticleStorage::iterator createNewParticle(data::ParticleStorage&
pIt->setRotation(data.rotation);
pIt->setAngularVelocity(data.angularVelocity);
pIt->setOldTorque(data.oldTorque);
+ pIt->setRadius(data.radius);
pIt->setType(data.type);
pIt->setOldContactHistory(data.oldContactHistory);
pIt->setTemperature(data.temperature);
@@ -138,6 +140,7 @@ mpi::GenericSendBuffer<T,G>& operator<<( mpi::GenericSendBuffer<T,G> & buf, cons
buf << obj.particle_.getRotation();
buf << obj.particle_.getAngularVelocity();
buf << obj.particle_.getOldTorque();
+ buf << obj.particle_.getRadius();
buf << obj.particle_.getType();
buf << obj.particle_.getOldContactHistory();
buf << obj.particle_.getTemperature();
@@ -165,6 +168,7 @@ mpi::GenericRecvBuffer<T>& operator>>( mpi::GenericRecvBuffer<T> & buf, mesa_pd:
buf >> objparam.rotation;
buf >> objparam.angularVelocity;
buf >> objparam.oldTorque;
+ buf >> objparam.radius;
buf >> objparam.type;
buf >> objparam.oldContactHistory;
buf >> objparam.temperature;
diff --git a/src/mesa_pd/mpi/notifications/ParticleGhostCopyNotification.h b/src/mesa_pd/mpi/notifications/ParticleGhostCopyNotification.h
index ca769b85c..aba424550 100644
--- a/src/mesa_pd/mpi/notifications/ParticleGhostCopyNotification.h
+++ b/src/mesa_pd/mpi/notifications/ParticleGhostCopyNotification.h
@@ -58,6 +58,7 @@ public:
size_t shapeID {};
walberla::mesa_pd::Rot3 rotation {};
walberla::mesa_pd::Vec3 angularVelocity {real_t(0)};
+ walberla::real_t radius {real_t(0)};
uint_t type {0};
std::map<walberla::id_t, walberla::mesa_pd::data::ContactHistory> oldContactHistory {};
walberla::real_t temperature {real_t(0)};
@@ -82,6 +83,7 @@ inline data::ParticleStorage::iterator createNewParticle(data::ParticleStorage&
pIt->setShapeID(data.shapeID);
pIt->setRotation(data.rotation);
pIt->setAngularVelocity(data.angularVelocity);
+ pIt->setRadius(data.radius);
pIt->setType(data.type);
pIt->setOldContactHistory(data.oldContactHistory);
pIt->setTemperature(data.temperature);
@@ -121,6 +123,7 @@ mpi::GenericSendBuffer<T,G>& operator<<( mpi::GenericSendBuffer<T,G> & buf, cons
buf << obj.particle_.getShapeID();
buf << obj.particle_.getRotation();
buf << obj.particle_.getAngularVelocity();
+ buf << obj.particle_.getRadius();
buf << obj.particle_.getType();
buf << obj.particle_.getOldContactHistory();
buf << obj.particle_.getTemperature();
@@ -141,6 +144,7 @@ mpi::GenericRecvBuffer<T>& operator>>( mpi::GenericRecvBuffer<T> & buf, mesa_pd:
buf >> objparam.shapeID;
buf >> objparam.rotation;
buf >> objparam.angularVelocity;
+ buf >> objparam.radius;
buf >> objparam.type;
buf >> objparam.oldContactHistory;
buf >> objparam.temperature;
diff --git a/src/mesa_pd/mpi/notifications/ParticleUpdateNotification.h b/src/mesa_pd/mpi/notifications/ParticleUpdateNotification.h
index 62d0c6da3..0214aa977 100644
--- a/src/mesa_pd/mpi/notifications/ParticleUpdateNotification.h
+++ b/src/mesa_pd/mpi/notifications/ParticleUpdateNotification.h
@@ -50,6 +50,7 @@ public:
walberla::mesa_pd::Vec3 linearVelocity {real_t(0)};
walberla::mesa_pd::Rot3 rotation {};
walberla::mesa_pd::Vec3 angularVelocity {real_t(0)};
+ walberla::real_t radius {real_t(0)};
std::map<walberla::id_t, walberla::mesa_pd::data::ContactHistory> oldContactHistory {};
walberla::real_t temperature {real_t(0)};
};
@@ -86,6 +87,7 @@ mpi::GenericSendBuffer<T,G>& operator<<( mpi::GenericSendBuffer<T,G> & buf, cons
buf << obj.particle_.getLinearVelocity();
buf << obj.particle_.getRotation();
buf << obj.particle_.getAngularVelocity();
+ buf << obj.particle_.getRadius();
buf << obj.particle_.getOldContactHistory();
buf << obj.particle_.getTemperature();
return buf;
@@ -100,6 +102,7 @@ mpi::GenericRecvBuffer<T>& operator>>( mpi::GenericRecvBuffer<T> & buf, mesa_pd:
buf >> objparam.linearVelocity;
buf >> objparam.rotation;
buf >> objparam.angularVelocity;
+ buf >> objparam.radius;
buf >> objparam.oldContactHistory;
buf >> objparam.temperature;
return buf;
--
GitLab