diff --git a/apps/showcases/CMakeLists.txt b/apps/showcases/CMakeLists.txt
index 536f56ff6acdbdb102faf394361bb2fd661af6e8..266582c2283dca2e7d8dd44973a1e6517027fb46 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 0000000000000000000000000000000000000000..f9465cda0679d5d0708213892a509bcc211c5bfc
--- /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 0000000000000000000000000000000000000000..c5ed90565ea564ef350e1ff5a1931010bb22641a
--- /dev/null
+++ b/apps/showcases/HeatConduction/ContactDetection.h
@@ -0,0 +1,106 @@
+//======================================================================================================================
+//
+// 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 collision_detection::AnalyticContactDetection
+{
+public:
+ ///import all collision detection functions from AnalyticContactDetection
+ using AnalyticContactDetection::operator();
+
+ ///overwrite functions that should be different
+ 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>
+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.getRadiusAtTemperature(getIdx1()),
+ ac.getPosition(getIdx2()),
+ ac.getRadiusAtTemperature(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.getRadiusAtTemperature(getIdx1()),
+ ac.getPosition(getIdx2()),
+ p.getNormal(),
+ getContactPoint(),
+ getContactNormal(),
+ getPenetrationDepth(),
+ getContactThreshold());
+}
+
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/apps/showcases/HeatConduction/HeatConduction.cfg b/apps/showcases/HeatConduction/HeatConduction.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..7d1d729826e3f7e72297c9b009663122afd8ff9b
--- /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 0000000000000000000000000000000000000000..eea1b2e726dca851debcbea3da698ff843ca8f82
--- /dev/null
+++ b/apps/showcases/HeatConduction/HeatConduction.cpp
@@ -0,0 +1,302 @@
+//======================================================================================================================
+//
+// 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>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+// This showcase combines a classic spring-dashpot interaction model with heat conduction and linear thermal expansion.
+//
+// First, a set of spheres is dropped onto a plate. The setup is periodic in horizontal directions. The plate as well
+// as the spheres shares the same initial temperature. After the spheres have settled, the plate is heated and cooled
+// down in a sinus wave. Through thermal conduction the heat is transported through the pile of spheres. At the same
+// time a linear thermal expansion model makes the spheres grow and shrink. The parameters are chosen for a visual
+// effect and are far to large for a realistic physical simulation. The evolution of the simulation is written to disk
+// as vtk files which can be visualized with paraview.
+//======================================================================================================================
+
+#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->setRadiusAtTemperature(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::SelectParticleRadiusAtTemperature>("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 0000000000000000000000000000000000000000..c56ee0732258220a0b8e67edbaaeb519843d2db3
--- /dev/null
+++ b/apps/showcases/HeatConduction/ThermalExpansion.h
@@ -0,0 +1,84 @@
+//======================================================================================================================
+//
+// 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(const uint_t numParticleTypes);
+ ThermalExpansion(const ThermalExpansion& other) = default;
+ ThermalExpansion(ThermalExpansion&& other) = default;
+ ThermalExpansion& operator=(const ThermalExpansion& other) = default;
+ ThermalExpansion& operator=(ThermalExpansion&& other) = default;
+
+ void setLinearExpansionCoefficient(const size_t type, const real_t& val);
+ real_t getLinearExpansionCoefficient(const size_t type) const;
+
+ template <typename Accessor>
+ void operator()(const size_t p_idx, Accessor& ac) const;
+private:
+ uint_t numParticleTypes_;
+
+ std::vector<real_t> linearExpansionCoefficient_ {};
+};
+
+ThermalExpansion::ThermalExpansion(const uint_t numParticleTypes)
+{
+ numParticleTypes_ = numParticleTypes;
+
+ linearExpansionCoefficient_.resize(numParticleTypes, real_t(0));
+}
+
+
+inline void ThermalExpansion::setLinearExpansionCoefficient(const size_t type, const real_t& val)
+{
+ WALBERLA_ASSERT_LESS( type, numParticleTypes_ );
+ linearExpansionCoefficient_[type] = val;
+}
+
+
+inline real_t ThermalExpansion::getLinearExpansionCoefficient(const size_t type) const
+{
+ WALBERLA_ASSERT_LESS( type, numParticleTypes_ );
+ return linearExpansionCoefficient_[type];
+}
+
+template <typename Accessor>
+inline void ThermalExpansion::operator()(const size_t p_idx,
+ Accessor& ac) const
+{
+ static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
+
+ const auto Tc = ac.getTemperature(p_idx)-real_t(273);
+ ac.setRadiusAtTemperature(p_idx, ac.getRadius273K(p_idx) * (real_t(1) + Tc * getLinearExpansionCoefficient(ac.getType(p_idx))));
+ ac.setInteractionRadius(p_idx, ac.getRadiusAtTemperature(p_idx));
+}
+
+} //namespace kernel
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/python/mesa_pd.py b/python/mesa_pd.py
index 785d9a4ff67129ae687343d8dc00a54a86796b19..a6d7fa151dbde6177ae28c2324fff3e2d6ae93ef 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("radiusAtTemperature", "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/python/mesa_pd/kernel/TemperatureIntegration.py b/python/mesa_pd/kernel/TemperatureIntegration.py
index 426a9b3874805db5abac7b686baf791d0a46263a..006eb37232df257eb42c5ac2c7d8be1c44ce72ab 100644
--- a/python/mesa_pd/kernel/TemperatureIntegration.py
+++ b/python/mesa_pd/kernel/TemperatureIntegration.py
@@ -9,11 +9,12 @@ class TemperatureIntegration:
self.context = {'interface': []}
self.context['interface'].append(create_access("temperature", "walberla::real_t", access="gs"))
self.context['interface'].append(create_access("heatFlux", "walberla::real_t", access="gs"))
+ self.context['interface'].append(create_access("invMass", "walberla::real_t", access="g"))
self.context['interface'].append(create_access("type", "uint_t", access="g"))
def generate(self, module):
ctx = {'module': module, **self.context}
- ctx["parameters"] = ["invHeatCapacity"]
+ ctx["parameters"] = ["invSpecificHeat"]
generate_file(module['module_path'], 'kernel/TemperatureIntegration.templ.h', ctx)
ctx["InterfaceTestName"] = "TemperatureIntegrationInterfaceCheck"
diff --git a/python/mesa_pd/templates/kernel/HeatConduction.templ.h b/python/mesa_pd/templates/kernel/HeatConduction.templ.h
index 8fe94e6da8e0681556ca441617909f44a6c47d5e..ef12443d9e070444ef050ec77df4c8bc5b2a95ba 100644
--- a/python/mesa_pd/templates/kernel/HeatConduction.templ.h
+++ b/python/mesa_pd/templates/kernel/HeatConduction.templ.h
@@ -38,9 +38,9 @@ namespace mesa_pd {
namespace kernel {
/**
- * Basic DEM kernel
+ * Heat conduction interaction kernel
*
- * This DEM kernel supports spring&dashpot in normal direction as well as friction in tangential direction.
+ * This kernel implements a simple heat conduction: \frac{dQ}{dt} = - \alpha (T_2 - T_1)
*
* \code
{%- for prop in interface %}
diff --git a/python/mesa_pd/templates/kernel/TemperatureIntegration.templ.h b/python/mesa_pd/templates/kernel/TemperatureIntegration.templ.h
index 92aa9f58d9759cf41d4bb708c38b1f4114f025fd..2bdd7f362e6412830c55f78b2f233505ddfc472e 100644
--- a/python/mesa_pd/templates/kernel/TemperatureIntegration.templ.h
+++ b/python/mesa_pd/templates/kernel/TemperatureIntegration.templ.h
@@ -34,8 +34,8 @@ namespace mesa_pd {
namespace kernel {
/**
- * Kernel which explicitly integrates all particles in time.
- * This integrator integrates velocity and position.
+ * Kernel which explicitly integrates a particle in time.
+ * The heat flux is converted into a temperature change.
*
* This kernel requires the following particle accessor interface
* \code
@@ -53,8 +53,8 @@ namespace kernel {
{%- endfor %}
* \endcode
*
- * \pre All forces acting on the particles have to be set.
- * \post All forces are reset to 0.
+ * \pre Heat flux has to be set/reduced.
+ * \post Heat flux is reset to 0.
* \ingroup mesa_pd_kernel
*/
class TemperatureIntegration
@@ -112,14 +112,14 @@ inline real_t TemperatureIntegration::get{{param | capFirst}}(const size_t type)
{%- endfor %}
template <typename Accessor>
-inline void TemperatureIntegration::operator()(const size_t idx,
+inline void TemperatureIntegration::operator()(const size_t p_idx,
Accessor& ac) const
{
static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
//formula for heat capacity
- ac.setTemperature(idx, getInvHeatCapacity(ac.getType(idx)) * ac.getHeatFlux(idx) * dt_ + ac.getTemperature(idx));
- ac.setHeatFlux (idx, real_t(0));
+ ac.setTemperature(p_idx, getInvSpecificHeat(ac.getType(p_idx)) * ac.getInvMass(p_idx) * ac.getHeatFlux(p_idx) * dt_ + ac.getTemperature(p_idx));
+ ac.setHeatFlux (p_idx, real_t(0));
}
} //namespace kernel
diff --git a/src/mesa_pd/data/ParticleAccessor.h b/src/mesa_pd/data/ParticleAccessor.h
index 59e5753306256ae967be958c297696804fd3a0c7..a00eb26c3bc6770288e31a978abaebdcc7f73f2e 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 & getRadiusAtTemperature(const size_t p_idx) const {return ps_->getRadiusAtTemperature(p_idx);}
+ walberla::real_t& getRadiusAtTemperatureRef(const size_t p_idx) {return ps_->getRadiusAtTemperatureRef(p_idx);}
+ void setRadiusAtTemperature(const size_t p_idx, walberla::real_t const & v) { ps_->setRadiusAtTemperature(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 & getRadiusAtTemperature(const size_t /*p_idx*/) const {return radiusAtTemperature_;}
+ void setRadiusAtTemperature(const size_t /*p_idx*/, walberla::real_t const & v) { radiusAtTemperature_ = v;}
+ walberla::real_t& getRadiusAtTemperatureRef(const size_t /*p_idx*/) {return radiusAtTemperature_;}
+
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 radiusAtTemperature_;
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 a22470d675486c78bfa154765b33bd197d0a97b2..8e248aca12c376fbf214481a7f9c552ffd5ccfcc 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 radiusAtTemperature_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);}
+ radiusAtTemperature_type const & getRadiusAtTemperature() const {return storage_.getRadiusAtTemperature(i_);}
+ radiusAtTemperature_type& getRadiusAtTemperatureRef() {return storage_.getRadiusAtTemperatureRef(i_);}
+ void setRadiusAtTemperature(radiusAtTemperature_type const & v) { storage_.setRadiusAtTemperature(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 radiusAtTemperature_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; }
+ radiusAtTemperature_type const & getRadiusAtTemperature(const size_t idx) const {return radiusAtTemperature_[idx];}
+ radiusAtTemperature_type& getRadiusAtTemperatureRef(const size_t idx) {return radiusAtTemperature_[idx];}
+ void setRadiusAtTemperature(const size_t idx, radiusAtTemperature_type const & v) { radiusAtTemperature_[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<radiusAtTemperature_type> radiusAtTemperature_ {};
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();
+ getRadiusAtTemperatureRef() = rhs.getRadiusAtTemperature();
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());
+ getRadiusAtTemperatureRef() = std::move(rhs.getRadiusAtTemperatureRef());
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.getRadiusAtTemperatureRef(), rhs.getRadiusAtTemperatureRef());
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" <<
+ "radiusAtTemperature : " << p.getRadiusAtTemperature() << "\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));
+ radiusAtTemperature_.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();
+ radiusAtTemperature_.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);
+ radiusAtTemperature_.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();
+ radiusAtTemperature_.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(), radiusAtTemperature.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 SelectParticleRadiusAtTemperature
+{
+public:
+ using return_type = walberla::real_t;
+ walberla::real_t& operator()(data::Particle& p) const {return p.getRadiusAtTemperatureRef();}
+ walberla::real_t& operator()(data::Particle&& p) const {return p.getRadiusAtTemperatureRef();}
+ walberla::real_t const & operator()(const data::Particle& p) const {return p.getRadiusAtTemperature();}
+};
+///Predicate that selects a certain property from a Particle
class SelectParticleCurrentBlock
{
public:
diff --git a/src/mesa_pd/kernel/HeatConduction.h b/src/mesa_pd/kernel/HeatConduction.h
index 12826c523265f3fd53be7b0ab96661a6ae3cfd0c..57e60cbccb272452e51490cbc5cfd0a35c66840f 100644
--- a/src/mesa_pd/kernel/HeatConduction.h
+++ b/src/mesa_pd/kernel/HeatConduction.h
@@ -38,9 +38,9 @@ namespace mesa_pd {
namespace kernel {
/**
- * Basic DEM kernel
+ * Heat conduction interaction kernel
*
- * This DEM kernel supports spring&dashpot in normal direction as well as friction in tangential direction.
+ * This kernel implements a simple heat conduction: \frac{dQ}{dt} = - \alpha (T_2 - T_1)
*
* \code
* const walberla::real_t& getTemperature(const size_t p_idx) const;
diff --git a/src/mesa_pd/kernel/TemperatureIntegration.h b/src/mesa_pd/kernel/TemperatureIntegration.h
index bdfc5732d622164fc05658cb1d30c35285fdd5e8..c90eb3ffcf3b65bb97e77a31d127eb43a02830a9 100644
--- a/src/mesa_pd/kernel/TemperatureIntegration.h
+++ b/src/mesa_pd/kernel/TemperatureIntegration.h
@@ -34,8 +34,8 @@ namespace mesa_pd {
namespace kernel {
/**
- * Kernel which explicitly integrates all particles in time.
- * This integrator integrates velocity and position.
+ * Kernel which explicitly integrates a particle in time.
+ * The heat flux is converted into a temperature change.
*
* This kernel requires the following particle accessor interface
* \code
@@ -45,12 +45,14 @@ namespace kernel {
* const walberla::real_t& getHeatFlux(const size_t p_idx) const;
* void setHeatFlux(const size_t p_idx, const walberla::real_t& v);
*
+ * const walberla::real_t& getInvMass(const size_t p_idx) const;
+ *
* const uint_t& getType(const size_t p_idx) const;
*
* \endcode
*
- * \pre All forces acting on the particles have to be set.
- * \post All forces are reset to 0.
+ * \pre Heat flux has to be set/reduced.
+ * \post Heat flux is reset to 0.
* \ingroup mesa_pd_kernel
*/
class TemperatureIntegration
@@ -67,16 +69,16 @@ public:
/// assumes this parameter is symmetric
- void setInvHeatCapacity(const size_t type, const real_t& val);
+ void setInvSpecificHeat(const size_t type, const real_t& val);
- real_t getInvHeatCapacity(const size_t type) const;
+ real_t getInvSpecificHeat(const size_t type) const;
private:
real_t dt_ = real_t(0.0);
uint_t numParticleTypes_;
- std::vector<real_t> invHeatCapacity_ {};
+ std::vector<real_t> invSpecificHeat_ {};
};
TemperatureIntegration::TemperatureIntegration(const real_t dt, const uint_t numParticleTypes)
@@ -84,32 +86,32 @@ TemperatureIntegration::TemperatureIntegration(const real_t dt, const uint_t num
{
numParticleTypes_ = numParticleTypes;
- invHeatCapacity_.resize(numParticleTypes, real_t(0));
+ invSpecificHeat_.resize(numParticleTypes, real_t(0));
}
-inline void TemperatureIntegration::setInvHeatCapacity(const size_t type, const real_t& val)
+inline void TemperatureIntegration::setInvSpecificHeat(const size_t type, const real_t& val)
{
WALBERLA_ASSERT_LESS( type, numParticleTypes_ );
- invHeatCapacity_[type] = val;
+ invSpecificHeat_[type] = val;
}
-inline real_t TemperatureIntegration::getInvHeatCapacity(const size_t type) const
+inline real_t TemperatureIntegration::getInvSpecificHeat(const size_t type) const
{
WALBERLA_ASSERT_LESS( type, numParticleTypes_ );
- return invHeatCapacity_[type];
+ return invSpecificHeat_[type];
}
template <typename Accessor>
-inline void TemperatureIntegration::operator()(const size_t idx,
+inline void TemperatureIntegration::operator()(const size_t p_idx,
Accessor& ac) const
{
static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
//formula for heat capacity
- ac.setTemperature(idx, getInvHeatCapacity(ac.getType(idx)) * ac.getHeatFlux(idx) * dt_ + ac.getTemperature(idx));
- ac.setHeatFlux (idx, real_t(0));
+ ac.setTemperature(p_idx, getInvSpecificHeat(ac.getType(p_idx)) * ac.getInvMass(p_idx) * ac.getHeatFlux(p_idx) * dt_ + ac.getTemperature(p_idx));
+ ac.setHeatFlux (p_idx, real_t(0));
}
} //namespace kernel
diff --git a/src/mesa_pd/mpi/notifications/ParseMessage.h b/src/mesa_pd/mpi/notifications/ParseMessage.h
index 63aa0d542f4d969e6659e64507fb3119335c5473..caf0e0f890312b7566ffd92cd3194728fb09f095 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->setRadiusAtTemperature(objparam.radiusAtTemperature);
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 2ac5383c17fdda618019c338ef36f47437da1cbf..b01a4dbeb92514fbed8eebcbbc54c6ce92b2f89e 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 radiusAtTemperature {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->setRadiusAtTemperature(data.radiusAtTemperature);
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_.getRadiusAtTemperature();
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.radiusAtTemperature;
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 ca769b85c986a66b3193b74570a80f1de32ce1c9..c4e1c2417b0ca1f02d1b7cbaabd193234085a605 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 radiusAtTemperature {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->setRadiusAtTemperature(data.radiusAtTemperature);
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_.getRadiusAtTemperature();
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.radiusAtTemperature;
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 62d0c6da38b6d57cee8279b2297d4979824b8d41..d845082000b48affcd17cbfecb77247a6d00662d 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 radiusAtTemperature {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_.getRadiusAtTemperature();
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.radiusAtTemperature;
buf >> objparam.oldContactHistory;
buf >> objparam.temperature;
return buf;
diff --git a/tests/mesa_pd/kernel/TemperatureIntegration.cpp b/tests/mesa_pd/kernel/TemperatureIntegration.cpp
index fefa30cadb52420b7714d39a663213b071e9827c..b84fc9cd285400852b837084603c7185f2fbc9ac 100644
--- a/tests/mesa_pd/kernel/TemperatureIntegration.cpp
+++ b/tests/mesa_pd/kernel/TemperatureIntegration.cpp
@@ -18,13 +18,11 @@
//
//======================================================================================================================
-#include <mesa_pd/data/DataTypes.h>
#include <mesa_pd/data/ParticleAccessor.h>
#include <mesa_pd/kernel/TemperatureIntegration.h>
#include <core/Environment.h>
-#include <core/logging/Logging.h>
#include <iostream>
@@ -45,11 +43,12 @@ int main( int argc, char ** argv )
accessor.setType( 0, 0 );
accessor.setHeatFlux( 0, real_t(8) );
accessor.setTemperature( 0, real_t(5) );
+ accessor.setInvMass( 0, real_t(1) );
//init kernels
const real_t dt = real_t(1);
kernel::TemperatureIntegration integrator( dt, 1 );
- integrator.setInvHeatCapacity( 0, real_t(2) );
+ integrator.setInvSpecificHeat( 0, real_t(2) );
integrator(0, accessor);
diff --git a/tests/mesa_pd/kernel/interfaces/TemperatureIntegrationInterfaceCheck.cpp b/tests/mesa_pd/kernel/interfaces/TemperatureIntegrationInterfaceCheck.cpp
index e764dffa1b0a1a88b597108d2b64fa6adc4060df..dfc209c39f2d3f6890d1017ab687064d1f7d9deb 100644
--- a/tests/mesa_pd/kernel/interfaces/TemperatureIntegrationInterfaceCheck.cpp
+++ b/tests/mesa_pd/kernel/interfaces/TemperatureIntegrationInterfaceCheck.cpp
@@ -44,6 +44,8 @@ public:
const walberla::real_t& getHeatFlux(const size_t /*p_idx*/) const {return heatFlux_;}
void setHeatFlux(const size_t /*p_idx*/, const walberla::real_t& v) { heatFlux_ = v;}
+ const walberla::real_t& getInvMass(const size_t /*p_idx*/) const {return invMass_;}
+
const uint_t& getType(const size_t /*p_idx*/) const {return type_;}
@@ -59,6 +61,7 @@ public:
private:
walberla::real_t temperature_;
walberla::real_t heatFlux_;
+ walberla::real_t invMass_;
uint_t type_;
};