diff --git a/python/mesa_pd.py b/python/mesa_pd.py index 0d94b5a92b101d927cabcd4e44c1a99302870142..d3e062a4f3603090763585531e537d990935e1bc 100755 --- a/python/mesa_pd.py +++ b/python/mesa_pd.py @@ -103,6 +103,7 @@ if __name__ == '__main__': mpd.add(kernel.InsertParticleIntoSparseLinkedCells()) mpd.add(kernel.LinearSpringDashpot()) mpd.add(kernel.NonLinearSpringDashpot()) + mpd.add(kernel.PFCDamping()) mpd.add(kernel.SemiImplicitEuler()) mpd.add(kernel.SingleCast(ps)) mpd.add(kernel.SpringDashpot()) diff --git a/python/mesa_pd/kernel/PFCDamping.py b/python/mesa_pd/kernel/PFCDamping.py new file mode 100644 index 0000000000000000000000000000000000000000..c3faf5e9829f80652bcf93e02cb7357ed51e5636 --- /dev/null +++ b/python/mesa_pd/kernel/PFCDamping.py @@ -0,0 +1,30 @@ +# -*- coding: utf-8 -*- + +from mesa_pd.accessor import create_access +from mesa_pd.utility import generate_file + + +class PFCDamping: + def __init__(self, integrate_rotation=True): + self.context = {'bIntegrateRotation': integrate_rotation, + 'interface': [ + create_access("linearVelocity", "walberla::mesa_pd::Vec3", access="gs"), + create_access("force", "walberla::mesa_pd::Vec3", access="gs"), + create_access("flags", "walberla::mesa_pd::data::particle_flags::FlagT", access="g") + ] + } + + if integrate_rotation: + self.context['interface'].append(create_access("angularVelocity", "walberla::mesa_pd::Vec3", access="gs")) + self.context['interface'].append(create_access("torque", "walberla::mesa_pd::Vec3", access="gs")) + + def generate(self, module): + ctx = {'module': module, **self.context} + generate_file(module['module_path'], 'kernel/PFCDamping.templ.h', ctx) + + ctx["InterfaceTestName"] = "PFCDampingInterfaceCheck" + ctx["KernelInclude"] = "kernel/PFCDamping.h" + ctx["ExplicitInstantiation"] = \ + "template void kernel::PFCDamping::operator()(const size_t p_idx1, Accessor& ac) const;" + generate_file(module['test_path'], 'tests/CheckInterface.templ.cpp', ctx, + 'kernel/interfaces/PFCDampingInterfaceCheck.cpp') diff --git a/python/mesa_pd/kernel/__init__.py b/python/mesa_pd/kernel/__init__.py index 4a639693bde73fbcb50cff624750a058e1519271..fda6d3c46645040607b4c5e14d0509967e4ccea6 100644 --- a/python/mesa_pd/kernel/__init__.py +++ b/python/mesa_pd/kernel/__init__.py @@ -12,6 +12,7 @@ from .InsertParticleIntoLinkedCells import InsertParticleIntoLinkedCells from .InsertParticleIntoSparseLinkedCells import InsertParticleIntoSparseLinkedCells from .LinearSpringDashpot import LinearSpringDashpot from .NonLinearSpringDashpot import NonLinearSpringDashpot +from .PFCDamping import PFCDamping from .SemiImplicitEuler import SemiImplicitEuler from .SingleCast import SingleCast from .SpringDashpot import SpringDashpot @@ -32,6 +33,7 @@ __all__ = ['DoubleCast', 'InsertParticleIntoSparseLinkedCells', 'LinearSpringDashpot', 'NonLinearSpringDashpot', + 'PFCDamping', 'SemiImplicitEuler', 'SingleCast', 'SpringDashpot', diff --git a/python/mesa_pd/templates/kernel/PFCDamping.templ.h b/python/mesa_pd/templates/kernel/PFCDamping.templ.h new file mode 100644 index 0000000000000000000000000000000000000000..7dfc8ae0c734ce06d978138f82b53bf4ad1e3742 --- /dev/null +++ b/python/mesa_pd/templates/kernel/PFCDamping.templ.h @@ -0,0 +1,90 @@ +//====================================================================================================================== +// +// 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 . +// +//! \file PFCDamping.h +//! \author Igor Ostanin +//! \author Grigorii Drozdov, +// +//====================================================================================================================== + +//====================================================================================================================== +// +// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!! +// +//====================================================================================================================== + +#pragma once + +#include +#include + +namespace walberla { +namespace mesa_pd { +namespace kernel { + +/** + * PFC style damping + * + * This kernel requires the following particle accessor interface + * \code + {%- for prop in interface %} + {%- if 'g' in prop.access %} + * const {{prop.type}}& get{{prop.name | capFirst}}(const size_t p_idx) const; + {%- endif %} + {%- if 's' in prop.access %} + * void set{{prop.name | capFirst}}(const size_t p_idx, const {{prop.type}}& v); + {%- endif %} + {%- if 'r' in prop.access %} + * {{prop.type}}& get{{prop.name | capFirst}}Ref(const size_t p_idx); + {%- endif %} + * + {%- endfor %} + * \endcode + * + * \ingroup mesa_pd_kernel + */ +class PFCDamping +{ +public: + PFCDamping(const real_t alpha) : alpha_(alpha) {} + + template + void operator()(const size_t p_idx, Accessor& ac) const; +private: + real_t alpha_ = 0_r; +}; + +template +inline void PFCDamping::operator()(const size_t p_idx, + Accessor& ac) const +{ + static_assert(std::is_base_of::value, "please provide a valid accessor"); + + Vec3 damp_F(0,0,0); + Vec3 damp_M(0,0,0); + + for (size_t i = 0; i < 3; i++) + { + damp_F[i] = - alpha_ * std::fabs( ac.getForce(p_idx)[i] ) * math::sign( ac.getLinearVelocity(p_idx)[i] ); + damp_M[i] = - alpha_ * std::fabs( ac.getTorque(p_idx)[i] ) * math::sign( ac.getAngularVelocity(p_idx)[i] ); + } + + ac.setForce (p_idx, ac.getForce(p_idx) + damp_F); + ac.setTorque(p_idx, ac.getTorque(p_idx) + damp_M); +} + +} //namespace kernel +} //namespace mesa_pd +} //namespace walberla diff --git a/src/mesa_pd/kernel/PFCDamping.h b/src/mesa_pd/kernel/PFCDamping.h new file mode 100644 index 0000000000000000000000000000000000000000..386df12296a8354e08a0b76c3026e200bc7ebccd --- /dev/null +++ b/src/mesa_pd/kernel/PFCDamping.h @@ -0,0 +1,92 @@ +//====================================================================================================================== +// +// 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 . +// +//! \file PFCDamping.h +//! \author Igor Ostanin +//! \author Grigorii Drozdov, +// +//====================================================================================================================== + +//====================================================================================================================== +// +// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!! +// +//====================================================================================================================== + +#pragma once + +#include +#include + +namespace walberla { +namespace mesa_pd { +namespace kernel { + +/** + * PFC style damping + * + * This kernel requires the following particle accessor interface + * \code + * const walberla::mesa_pd::Vec3& getLinearVelocity(const size_t p_idx) const; + * void setLinearVelocity(const size_t p_idx, const walberla::mesa_pd::Vec3& v); + * + * const walberla::mesa_pd::Vec3& getForce(const size_t p_idx) const; + * void setForce(const size_t p_idx, const walberla::mesa_pd::Vec3& v); + * + * const walberla::mesa_pd::data::particle_flags::FlagT& getFlags(const size_t p_idx) const; + * + * const walberla::mesa_pd::Vec3& getAngularVelocity(const size_t p_idx) const; + * void setAngularVelocity(const size_t p_idx, const walberla::mesa_pd::Vec3& v); + * + * const walberla::mesa_pd::Vec3& getTorque(const size_t p_idx) const; + * void setTorque(const size_t p_idx, const walberla::mesa_pd::Vec3& v); + * + * \endcode + * + * \ingroup mesa_pd_kernel + */ +class PFCDamping +{ +public: + PFCDamping(const real_t alpha) : alpha_(alpha) {} + + template + void operator()(const size_t p_idx, Accessor& ac) const; +private: + real_t alpha_ = 0_r; +}; + +template +inline void PFCDamping::operator()(const size_t p_idx, + Accessor& ac) const +{ + static_assert(std::is_base_of::value, "please provide a valid accessor"); + + Vec3 damp_F(0,0,0); + Vec3 damp_M(0,0,0); + + for (size_t i = 0; i < 3; i++) + { + damp_F[i] = - alpha_ * std::fabs( ac.getForce(p_idx)[i] ) * math::sign( ac.getLinearVelocity(p_idx)[i] ); + damp_M[i] = - alpha_ * std::fabs( ac.getTorque(p_idx)[i] ) * math::sign( ac.getAngularVelocity(p_idx)[i] ); + } + + ac.setForce (p_idx, ac.getForce(p_idx) + damp_F); + ac.setTorque(p_idx, ac.getTorque(p_idx) + damp_M); +} + +} //namespace kernel +} //namespace mesa_pd +} //namespace walberla \ No newline at end of file diff --git a/tests/mesa_pd/CMakeLists.txt b/tests/mesa_pd/CMakeLists.txt index caf3ba2752b2630203a42de911480bb4969eda00..ad034a098c25011fb9a54971dfa080e8a0647487 100644 --- a/tests/mesa_pd/CMakeLists.txt +++ b/tests/mesa_pd/CMakeLists.txt @@ -128,6 +128,9 @@ waLBerla_execute_test( NAME MESA_PD_Kernel_LinearSpringDashpot ) waLBerla_compile_test( NAME MESA_PD_Kernel_LinkedCellsVsBruteForce FILES kernel/LinkedCellsVsBruteForce.cpp DEPENDS core ) waLBerla_execute_test( NAME MESA_PD_Kernel_LinkedCellsVsBruteForce PROCESSES 27 ) +waLBerla_compile_test( NAME MESA_PD_Kernel_PFCDamping FILES kernel/PFCDamping.cpp DEPENDS core ) +waLBerla_execute_test( NAME MESA_PD_Kernel_PFCDamping ) + waLBerla_compile_test( NAME MESA_PD_Kernel_SemiImplicitEuler FILES kernel/SemiImplicitEuler.cpp DEPENDS core ) waLBerla_execute_test( NAME MESA_PD_Kernel_SemiImplicitEuler ) diff --git a/tests/mesa_pd/kernel/PFCDamping.cpp b/tests/mesa_pd/kernel/PFCDamping.cpp new file mode 100644 index 0000000000000000000000000000000000000000..0329ef91128fc04edbc188124b925f978d16ca5f --- /dev/null +++ b/tests/mesa_pd/kernel/PFCDamping.cpp @@ -0,0 +1,64 @@ +//====================================================================================================================== +// +// This file is part of waLBerla. waLBerla is free software: you can +// redistribute it and/or modify it under the terms of the GNU General Public +// License as published by the Free Software Foundation, either version 3 of +// the License, or (at your option) any later version. +// +// waLBerla is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +// for more details. +// +// You should have received a copy of the GNU General Public License along +// with waLBerla (see COPYING.txt). If not, see . +// +//! \file PFCDamping.cpp +//! \author Igor Ostanin +//! \author Grigorii Drozdov +// +//====================================================================================================================== + +#include +#include + +#include + +#include + +namespace walberla { + +using namespace walberla::mesa_pd; + +int main( int argc, char ** argv ) +{ + Environment env(argc, argv); + WALBERLA_UNUSED(env); + mpi::MPIManager::instance()->useWorldComm(); + + //init data structures + data::SingleParticleAccessor ac; + + ac.setLinearVelocity( 0, Vec3(+2_r,-2_r,+2_r) ); + ac.setForce ( 0, Vec3(+2_r,+3_r,-4_r) ); + + ac.setAngularVelocity( 0, Vec3(+2_r,-2_r,+2_r) ); + ac.setTorque ( 0, Vec3(+3_r,+5_r,-2_r) ); + + //init kernels + kernel::PFCDamping damping( 0.1_r ); + + damping(0, ac); + + WALBERLA_CHECK_FLOAT_EQUAL(ac.getForce(0), Vec3(1.8_r, 3.3_r, -4.4_r)); + WALBERLA_CHECK_FLOAT_EQUAL(ac.getTorque(0), Vec3(2.7_r, 5.5_r, -2.2_r)); + + return EXIT_SUCCESS; +} + +} //namespace walberla + +int main( int argc, char ** argv ) +{ + return walberla::main(argc, argv); +} diff --git a/tests/mesa_pd/kernel/interfaces/PFCDampingInterfaceCheck.cpp b/tests/mesa_pd/kernel/interfaces/PFCDampingInterfaceCheck.cpp new file mode 100644 index 0000000000000000000000000000000000000000..3f52ab357266f412f2ae362dea55fcc1cd000c13 --- /dev/null +++ b/tests/mesa_pd/kernel/interfaces/PFCDampingInterfaceCheck.cpp @@ -0,0 +1,76 @@ +//====================================================================================================================== +// +// 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 . +// +//! \file PFCDampingInterfaceCheck.cpp +//! \author Sebastian Eibl +// +//====================================================================================================================== + +//====================================================================================================================== +// +// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!! +// +//====================================================================================================================== + +#include +#include + +#include + +#include + +namespace walberla { +namespace mesa_pd { + +class Accessor : public data::IAccessor +{ +public: + virtual ~Accessor() = default; + const walberla::mesa_pd::Vec3& getLinearVelocity(const size_t /*p_idx*/) const {return linearVelocity_;} + void setLinearVelocity(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { linearVelocity_ = v;} + + const walberla::mesa_pd::Vec3& getForce(const size_t /*p_idx*/) const {return force_;} + void setForce(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { force_ = v;} + + const walberla::mesa_pd::data::particle_flags::FlagT& getFlags(const size_t /*p_idx*/) const {return flags_;} + + const walberla::mesa_pd::Vec3& getAngularVelocity(const size_t /*p_idx*/) const {return angularVelocity_;} + void setAngularVelocity(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { angularVelocity_ = v;} + + const walberla::mesa_pd::Vec3& getTorque(const size_t /*p_idx*/) const {return torque_;} + void setTorque(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { torque_ = v;} + + + id_t getInvalidUid() const {return UniqueID::invalidID();} + size_t getInvalidIdx() const {return std::numeric_limits::max();} + /** + * @brief Returns the index of particle specified by uid. + * @param uid unique id of the particle to be looked up + * @return the index of the particle or std::numeric_limits::max() if the particle is not found + */ + size_t uidToIdx(const id_t& /*uid*/) const {return 0;} + size_t size() const { return 1; } +private: + walberla::mesa_pd::Vec3 linearVelocity_; + walberla::mesa_pd::Vec3 force_; + walberla::mesa_pd::data::particle_flags::FlagT flags_; + walberla::mesa_pd::Vec3 angularVelocity_; + walberla::mesa_pd::Vec3 torque_; +}; + +template void kernel::PFCDamping::operator()(const size_t p_idx1, Accessor& ac) const; + +} //namespace mesa_pd +} //namespace walberla \ No newline at end of file