Commit 95121137 authored by Sebastian Eibl's avatar Sebastian Eibl
Browse files

[ADD] heat conduction showcase

parent 6368d63b
add_subdirectory( BidisperseFluidizedBed )
add_subdirectory( CombinedResolvedUnresolved )
add_subdirectory( HeatConduction )
add_subdirectory( Mixer )
waLBerla_link_files_to_builddir( *.cfg )
waLBerla_add_executable( NAME HeatConduction
FILES HeatConduction.cpp
DEPENDS core mesa_pd vtk )
//======================================================================================================================
//
// 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
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;
}
//======================================================================================================================
//
// 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 );
}
//======================================================================================================================
//
// 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
......@@ -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")
......
......@@ -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_;
......
......@@ -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;