Merge branch 'pfc_damping' into 'master'

[ADD] PFCDamping kernel to mesa_pd

See merge request walberla/walberla!369

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())

python/mesa_pd/kernel/PFCDamping.py

+# -*- 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')

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',

python/mesa_pd/templates/kernel/PFCDamping.templ.h

+//======================================================================================================================
+//
+//  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

src/mesa_pd/kernel/PFCDamping.h

+//======================================================================================================================
+//
+//  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

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 )
 

tests/mesa_pd/kernel/PFCDamping.cpp

+//======================================================================================================================
+//
+//  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);
+}

tests/mesa_pd/kernel/interfaces/PFCDampingInterfaceCheck.cpp

+//======================================================================================================================
+//
+//  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