From e562fd7979e5f41163be4faf1bea92776d58a06d Mon Sep 17 00:00:00 2001
From: Sebastian Eibl <sebastian.eibl@fau.de>
Date: Fri, 4 Dec 2020 16:29:26 +0100
Subject: [PATCH] [ADD] PFCDamping kernel to mesa_pd
---
python/mesa_pd.py | 1 +
python/mesa_pd/kernel/PFCDamping.py | 30 ++++++
python/mesa_pd/kernel/__init__.py | 2 +
.../templates/kernel/PFCDamping.templ.h | 90 ++++++++++++++++++
src/mesa_pd/kernel/PFCDamping.h | 92 +++++++++++++++++++
tests/mesa_pd/CMakeLists.txt | 3 +
tests/mesa_pd/kernel/PFCDamping.cpp | 64 +++++++++++++
.../interfaces/PFCDampingInterfaceCheck.cpp | 76 +++++++++++++++
8 files changed, 358 insertions(+)
create mode 100644 python/mesa_pd/kernel/PFCDamping.py
create mode 100644 python/mesa_pd/templates/kernel/PFCDamping.templ.h
create mode 100644 src/mesa_pd/kernel/PFCDamping.h
create mode 100644 tests/mesa_pd/kernel/PFCDamping.cpp
create mode 100644 tests/mesa_pd/kernel/interfaces/PFCDampingInterfaceCheck.cpp
diff --git a/python/mesa_pd.py b/python/mesa_pd.py
index 0d94b5a92..d3e062a4f 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 000000000..c3faf5e98
--- /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 4a639693b..fda6d3c46 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 000000000..7dfc8ae0c
--- /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 <http://www.gnu.org/licenses/>.
+//
+//! \file PFCDamping.h
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov, <drozd013@umn.edu>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/IAccessor.h>
+
+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 <typename Accessor>
+ void operator()(const size_t p_idx, Accessor& ac) const;
+private:
+ real_t alpha_ = 0_r;
+};
+
+template <typename Accessor>
+inline void PFCDamping::operator()(const size_t p_idx,
+ Accessor& ac) const
+{
+ static_assert(std::is_base_of<data::IAccessor, Accessor>::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 000000000..386df1229
--- /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 <http://www.gnu.org/licenses/>.
+//
+//! \file PFCDamping.h
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov, <drozd013@umn.edu>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/IAccessor.h>
+
+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 <typename Accessor>
+ void operator()(const size_t p_idx, Accessor& ac) const;
+private:
+ real_t alpha_ = 0_r;
+};
+
+template <typename Accessor>
+inline void PFCDamping::operator()(const size_t p_idx,
+ Accessor& ac) const
+{
+ static_assert(std::is_base_of<data::IAccessor, Accessor>::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 caf3ba275..ad034a098 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 000000000..0329ef911
--- /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 <http://www.gnu.org/licenses/>.
+//
+//! \file PFCDamping.cpp
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//
+//======================================================================================================================
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/ParticleAccessor.h>
+
+#include <mesa_pd/kernel/PFCDamping.h>
+
+#include <core/Environment.h>
+
+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 000000000..3f52ab357
--- /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 <http://www.gnu.org/licenses/>.
+//
+//! \file PFCDampingInterfaceCheck.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
+//
+//======================================================================================================================
+
+#include <mesa_pd/data/IAccessor.h>
+#include <mesa_pd/kernel/PFCDamping.h>
+
+#include <core/UniqueID.h>
+
+#include <map>
+
+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<int>::invalidID();}
+ size_t getInvalidIdx() const {return std::numeric_limits<size_t>::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<size_t>::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
--
GitLab