diff --git a/apps/benchmarks/CMakeLists.txt b/apps/benchmarks/CMakeLists.txt
index 01d57440109ea02878e19231b10a36e012d93d8c..aabba8b459bbd6e7dbe05c5e42a34728824bbb4d 100644
--- a/apps/benchmarks/CMakeLists.txt
+++ b/apps/benchmarks/CMakeLists.txt
@@ -7,6 +7,7 @@ add_subdirectory( CouetteFlow )
add_subdirectory( FluidParticleCoupling )
add_subdirectory( ForcesOnSphereNearPlaneInShearFlow )
add_subdirectory( GranularGas )
+add_subdirectory( IntegratorAccuracy )
add_subdirectory( LennardJones )
add_subdirectory( NonUniformGrid )
add_subdirectory( MotionSingleHeavySphere )
diff --git a/apps/benchmarks/IntegratorAccuracy/CMakeLists.txt b/apps/benchmarks/IntegratorAccuracy/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8e95e831cc1fb275731b12bc14460cf64eedc729
--- /dev/null
+++ b/apps/benchmarks/IntegratorAccuracy/CMakeLists.txt
@@ -0,0 +1,4 @@
+waLBerla_link_files_to_builddir( "*.ipynb" )
+
+waLBerla_add_executable ( NAME IntegratorAccuracy
+ DEPENDS core mesa_pd )
diff --git a/apps/benchmarks/IntegratorAccuracy/IntegratorAccuracy.cpp b/apps/benchmarks/IntegratorAccuracy/IntegratorAccuracy.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..af2dc8a3aaa54f805822eabcbd4139fb49d3cdf7
--- /dev/null
+++ b/apps/benchmarks/IntegratorAccuracy/IntegratorAccuracy.cpp
@@ -0,0 +1,274 @@
+//======================================================================================================================
+//
+// 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 IntegratorAccuracy.cpp
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "mesa_pd/data/ParticleAccessor.h"
+
+#include "mesa_pd/kernel/ExplicitEuler.h"
+#include "mesa_pd/kernel/SemiImplicitEuler.h"
+#include "mesa_pd/kernel/VelocityVerlet.h"
+
+#include "core/Environment.h"
+#include "core/math/all.h"
+
+#include <iostream>
+
+namespace walberla {
+
+using namespace walberla::mesa_pd;
+
+class SingleParticleAccessorWithShape : public data::SingleParticleAccessor
+{
+public:
+ void setInvMass(const size_t /*p_idx*/, const real_t &val)
+ { invMass_ = val; }
+
+ const auto &getInvMass(const size_t /*p_idx*/) const
+ { return invMass_; }
+
+ void setInvInertiaBF(const size_t /*p_idx*/, const Mat3 &val)
+ { invInertiaBF_ = val; }
+
+ const auto &getInvInertiaBF(const size_t /*p_idx*/) const
+ { return invInertiaBF_; }
+
+private:
+ real_t invMass_;
+ Mat3 invInertiaBF_;
+};
+
+struct Oscillator
+{
+ real_t amplitude = 1.5_r;
+ real_t k = 0.1_r;
+ real_t damping = 0_r;
+ real_t mass = 0.9_r;
+ real_t dt = 0.2_r;
+ real_t phaseFraction = 0_r;
+ real_t periods = 10_r;
+ real_t phase = phaseFraction * math::pi;
+ real_t dampingRatio = damping / (2_r * std::sqrt(mass * k));
+ real_t omega = std::sqrt(k / mass) * std::sqrt(1_r - dampingRatio * dampingRatio);
+ real_t decay = std::sqrt(k / mass) * dampingRatio;
+ real_t durationOnePeriod = 2_r * math::pi / omega;
+ uint_t timeSteps = uint_c(periods * durationOnePeriod / dt);
+
+ Oscillator(int argc, char **argv)
+ {
+ for (int i = 1; i < argc; ++i)
+ {
+ if (std::strcmp(argv[i], "--dt") == 0)
+ {
+ dt = real_c(std::atof(argv[++i]));
+ continue;
+ }
+ if (std::strcmp(argv[i], "--amplitude") == 0)
+ {
+ amplitude = real_c(std::atof(argv[++i]));
+ continue;
+ }
+ if (std::strcmp(argv[i], "--k") == 0)
+ {
+ k = real_c(std::atof(argv[++i]));
+ continue;
+ }
+ if (std::strcmp(argv[i], "--damping") == 0)
+ {
+ damping = real_c(std::atof(argv[++i]));
+ continue;
+ }
+ if (std::strcmp(argv[i], "--mass") == 0)
+ {
+ mass = real_c(std::atof(argv[++i]));
+ continue;
+ }
+ if (std::strcmp(argv[i], "--phaseFraction") == 0)
+ {
+ phaseFraction = real_c(std::atof(argv[++i]));
+ continue;
+ }
+ if (std::strcmp(argv[i], "--periods") == 0)
+ {
+ periods = real_c(std::atof(argv[++i]));
+ continue;
+ }
+ WALBERLA_ABORT("Unrecognized command line argument found: " << argv[i]);
+ }
+
+ update();
+ }
+
+ void update()
+ {
+ phase = phaseFraction * math::pi;
+ dampingRatio = damping / (2_r * std::sqrt(mass * k));
+ omega = std::sqrt(k / mass) * std::sqrt(1_r - dampingRatio * dampingRatio);
+ decay = std::sqrt(k / mass) * dampingRatio;
+ durationOnePeriod = 2_r * math::pi / omega;
+ timeSteps = uint_c(periods * durationOnePeriod / dt);
+ }
+
+ Vec3 getForce(const Vec3 &pos, const Vec3 &vel) const
+ {
+ return -k * pos - damping * vel;
+ }
+
+ real_t getEnergy(const real_t &pos, const real_t &vel) const
+ {
+ return 0.5_r * mass * vel * vel + 0.5_r * k * pos * pos;
+ }
+
+ real_t analyticalPos(real_t t) const
+ {
+ return amplitude * std::exp(-decay * t) * std::cos(omega * t + phase);
+ }
+
+ real_t analyticalVel(real_t t) const
+ {
+ return -decay * amplitude * std::exp(-decay * t) * std::cos(omega * t + phase)
+ -amplitude * std::exp(-decay * t) * omega * std::sin(omega * t + phase);
+ }
+};
+
+struct ExplicitEuler
+{
+ ExplicitEuler(real_t dt) : integrator(dt) {}
+ void operator()(SingleParticleAccessorWithShape& particle,
+ const Oscillator& osc)
+ {
+ particle.setForce(0, osc.getForce(particle.getPosition(0),
+ particle.getLinearVelocity(0)));
+ integrator(0, particle);
+ }
+ kernel::ExplicitEuler integrator;
+};
+
+struct SemiImplicitEuler
+{
+ SemiImplicitEuler(real_t dt) : integrator(dt) {}
+ void operator()(SingleParticleAccessorWithShape& particle,
+ const Oscillator& osc)
+ {
+ particle.setForce(0, osc.getForce(particle.getPosition(0),
+ particle.getLinearVelocity(0)));
+ integrator(0, particle);
+ }
+ kernel::SemiImplicitEuler integrator;
+};
+
+struct VelocityVerlet
+{
+ VelocityVerlet(real_t dt) : preVV(dt), postVV(dt) {}
+ void operator()(SingleParticleAccessorWithShape& particle,
+ const Oscillator& osc)
+ {
+ preVV(0, particle);
+ particle.setForce(0, osc.getForce(particle.getPosition(0),
+ particle.getLinearVelocity(0)));
+ postVV(0, particle);
+ }
+ kernel::VelocityVerletPreForceUpdate preVV;
+ kernel::VelocityVerletPostForceUpdate postVV;
+};
+
+struct AccuracyResult
+{
+ real_t maxPosDeviation;
+ real_t maxVelDeviation;
+ real_t maxEneDeviation;
+};
+
+template <typename Integrator>
+AccuracyResult checkIntegrator(const Oscillator& osc)
+{
+ //init data structures
+ SingleParticleAccessorWithShape particle;
+
+ //first dummy argument is needed to fulfill accessor interface
+ particle.setPosition(0, Vec3(0, 0, osc.analyticalPos(0_r)));
+ particle.setLinearVelocity(0, Vec3(0, 0, osc.analyticalVel(0_r)));
+ particle.setInvMass(0, 1_r / osc.mass);
+ particle.setForce(0, osc.getForce(Vec3(0, 0, osc.analyticalPos(real_t(0))),
+ Vec3(0, 0, osc.analyticalVel(real_t(0)))));
+ particle.setOldForce(0, osc.getForce(Vec3(0, 0, osc.analyticalPos(-osc.dt)),
+ Vec3(0, 0, osc.analyticalVel(-osc.dt))));
+
+ // explicit euler
+ Integrator integrator(osc.dt);
+
+ real_t maxPosDeviation = 0_r;
+ real_t maxVelDeviation = 0_r;
+ real_t maxEneDeviation = 0_r;
+
+ for (auto i = uint_t(0); i <= osc.timeSteps; ++i)
+ {
+ real_t refPos = osc.analyticalPos(real_c(i) * osc.dt);
+ real_t refVel = osc.analyticalVel(real_c(i) * osc.dt);
+ real_t refEne = osc.getEnergy(refPos, refVel);
+
+ maxPosDeviation = std::max(maxPosDeviation, std::abs(particle.getPosition(0)[2] - refPos));
+ maxVelDeviation = std::max(maxVelDeviation, std::abs(particle.getLinearVelocity(0)[2] - refVel));
+ maxEneDeviation = std::max(maxEneDeviation, std::abs(osc.getEnergy(particle.getPosition(0)[2], particle.getLinearVelocity(0)[2]) - refEne));
+
+ std::cout << real_t(i) * osc.dt << " "
+ << refPos << " "
+ << refVel << " "
+ << refEne << " "
+ << particle.getPosition(0)[2] << " "
+ << particle.getLinearVelocity(0)[2] << " "
+ << osc.getEnergy(particle.getPosition(0)[2], particle.getLinearVelocity(0)[2]) << " "
+ << maxPosDeviation << " "
+ << maxVelDeviation << std::endl;
+
+ integrator(particle, osc);
+ }
+
+ return {maxPosDeviation, maxVelDeviation, maxEneDeviation};
+}
+
+int main(int argc, char **argv)
+{
+ mpi::Environment env(argc, argv);
+ WALBERLA_UNUSED(env);
+ mpi::MPIManager::instance()->useWorldComm();
+
+ Oscillator osc(argc, argv);
+
+ AccuracyResult res;
+// res = walberla::checkIntegrator<walberla::ExplicitEuler>(osc);
+// res = walberla::checkIntegrator<walberla::SemiImplicitEuler>(osc);
+ res = walberla::checkIntegrator<walberla::VelocityVerlet>(osc);
+
+ return EXIT_SUCCESS;
+}
+
+} //namespace walberla
+
+
+/*
+ * Simulates a harmonic oscillator to test the accuracy of the integrators.
+ * Playground for integrator analysis. The corresponding unit test is located at
+ * tests/mesa_pd/kernel/IntegratorAccuracy.cpp
+ */
+int main(int argc, char **argv)
+{
+ return walberla::main(argc, argv);
+}
+
diff --git a/apps/benchmarks/IntegratorAccuracy/IntegratorAccuracy.ipynb b/apps/benchmarks/IntegratorAccuracy/IntegratorAccuracy.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..f0cc8feef715c2af5e8e6192175d4cec00f09b44
--- /dev/null
+++ b/apps/benchmarks/IntegratorAccuracy/IntegratorAccuracy.ipynb
@@ -0,0 +1,156 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# error analysis of integration kernels"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import matplotlib.pyplot as plt\n",
+ "import numpy as np"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "integrator = 'euler'\n",
+ "data01 = np.loadtxt(f'{integrator}01.txt').transpose()\n",
+ "data02 = np.loadtxt(f'{integrator}02.txt').transpose()\n",
+ "data04 = np.loadtxt(f'{integrator}04.txt').transpose()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(data01[0], data01[7], label='dt = 0.1')\n",
+ "plt.plot(data02[0], data02[7], label='dt = 0.2')\n",
+ "plt.plot(data04[0], data04[7], label='dt = 0.4')\n",
+ "plt.xlabel('time')\n",
+ "plt.ylabel('max pos error')\n",
+ "plt.legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(data01[0], data01[8], label='dt = 0.1')\n",
+ "plt.plot(data02[0], data02[8], label='dt = 0.2')\n",
+ "plt.plot(data04[0], data04[8], label='dt = 0.4')\n",
+ "plt.xlabel('time')\n",
+ "plt.ylabel('max vel error')\n",
+ "plt.legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "final_errors = np.ones(3)\n",
+ "final_errors[0] = np.abs(data01[7][-1])\n",
+ "final_errors[1] = np.abs(data02[7][-1])\n",
+ "final_errors[2] = np.abs(data04[7][-1])\n",
+ "\n",
+ "plt.bar(np.arange(3), final_errors)\n",
+ "plt.show()\n",
+ "\n",
+ "print(final_errors[2] / final_errors[1])\n",
+ "print(final_errors[1] / final_errors[0])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(data01[0], data01[1], label='analytic')\n",
+ "plt.plot(data01[0], data01[4], label='dt = 0.1')\n",
+ "plt.plot(data02[0], data02[4], label='dt = 0.2')\n",
+ "plt.plot(data04[0], data04[4], label='dt = 0.4')\n",
+ "plt.xlabel('time')\n",
+ "plt.ylabel('pos')\n",
+ "plt.legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(data01[0], data01[2], label='analytic')\n",
+ "plt.plot(data01[0], data01[5], label='dt = 0.1')\n",
+ "plt.plot(data02[0], data02[5], label='dt = 0.2')\n",
+ "plt.plot(data04[0], data04[5], label='dt = 0.4')\n",
+ "plt.xlabel('time')\n",
+ "plt.ylabel('vel')\n",
+ "plt.legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(data01[0], data01[6] / data01[3][0], label='dt = 0.1')\n",
+ "plt.plot(data02[0], data02[6] / data01[3][0], label='dt = 0.2')\n",
+ "plt.plot(data04[0], data04[6] / data01[3][0], label='dt = 0.4')\n",
+ "plt.xlabel('time')\n",
+ "plt.ylabel('relative energy')\n",
+ "plt.legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.7.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/python/mesa_pd.py b/python/mesa_pd.py
index cc8b5309131ae2e90b3cd32a04bd77599fef9256..0d94b5a92b101d927cabcd4e44c1a99302870142 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.SemiImplicitEuler())
mpd.add(kernel.SingleCast(ps))
mpd.add(kernel.SpringDashpot())
mpd.add(kernel.SpringDashpotSpring())
diff --git a/python/mesa_pd/kernel/SemiImplicitEuler.py b/python/mesa_pd/kernel/SemiImplicitEuler.py
new file mode 100644
index 0000000000000000000000000000000000000000..f5b6632a1f57d7dba5631896a3bd26b5a5cce177
--- /dev/null
+++ b/python/mesa_pd/kernel/SemiImplicitEuler.py
@@ -0,0 +1,34 @@
+# -*- coding: utf-8 -*-
+
+from mesa_pd.accessor import create_access
+from mesa_pd.utility import generate_file
+
+
+class SemiImplicitEuler:
+ def __init__(self, integrate_rotation=True):
+ self.context = {'bIntegrateRotation': integrate_rotation,
+ 'interface': [
+ create_access("position", "walberla::mesa_pd::Vec3", access="gs"),
+ create_access("linearVelocity", "walberla::mesa_pd::Vec3", access="gs"),
+ create_access("invMass", "walberla::real_t", access="g"),
+ 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("rotation", "walberla::mesa_pd::Rot3", access="gs"))
+ self.context['interface'].append(create_access("angularVelocity", "walberla::mesa_pd::Vec3", access="gs"))
+ self.context['interface'].append(create_access("invInertiaBF", "walberla::mesa_pd::Mat3", access="g"))
+ 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/SemiImplicitEuler.templ.h', ctx)
+
+ ctx["InterfaceTestName"] = "SemiImplicitEulerInterfaceCheck"
+ ctx["KernelInclude"] = "kernel/SemiImplicitEuler.h"
+ ctx["SemiImplicitEulerInstantiation"] = \
+ "template void kernel::SemiImplicitEuler::operator()(const size_t p_idx1, Accessor& ac) const;"
+ generate_file(module['test_path'], 'tests/CheckInterface.templ.cpp', ctx,
+ 'kernel/interfaces/SemiImplicitEulerInterfaceCheck.cpp')
diff --git a/python/mesa_pd/kernel/__init__.py b/python/mesa_pd/kernel/__init__.py
index 317646d1e961b70c81c12e877161b353d9840bd8..4a639693bde73fbcb50cff624750a058e1519271 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 .SemiImplicitEuler import SemiImplicitEuler
from .SingleCast import SingleCast
from .SpringDashpot import SpringDashpot
from .SpringDashpotSpring import SpringDashpotSpring
@@ -31,6 +32,7 @@ __all__ = ['DoubleCast',
'InsertParticleIntoSparseLinkedCells',
'LinearSpringDashpot',
'NonLinearSpringDashpot',
+ 'SemiImplicitEuler',
'SingleCast',
'SpringDashpot',
'SpringDashpotSpring',
diff --git a/python/mesa_pd/templates/kernel/ExplicitEuler.templ.h b/python/mesa_pd/templates/kernel/ExplicitEuler.templ.h
index e3916291b679d64f6559e0b35d4b651690f75dc7..79b1518b3f32e6557b49463edf4a8f8f07ffe46e 100644
--- a/python/mesa_pd/templates/kernel/ExplicitEuler.templ.h
+++ b/python/mesa_pd/templates/kernel/ExplicitEuler.templ.h
@@ -34,7 +34,12 @@ namespace mesa_pd {
namespace kernel {
/**
- * Kernel which explicitly integrates all particles in time.
+ * Explicit Euler integration.
+ * Uses the 0.5at^2 extension for position integration.
+ * Boils down to using v_{t+0.5dt} for position integration.
+ * This version exhibits increased stability compared to standard explicit euler.
+ *
+ * Wachs, A. Particle-scale computational approaches to model dry and saturated granular flows of non-Brownian, non-cohesive, and non-spherical rigid bodies. Acta Mech 230, 1919–1980 (2019). https://doi.org/10.1007/s00707-019-02389-9
*
* This kernel requires the following particle accessor interface
* \code
@@ -75,15 +80,19 @@ inline void ExplicitEuler::operator()(const size_t idx,
if (!data::particle_flags::isSet( ac.getFlags(idx), data::particle_flags::FIXED))
{
- ac.setPosition (idx, ac.getInvMass(idx) * ac.getForce(idx) * dt_ * dt_ + ac.getLinearVelocity(idx) * dt_ + ac.getPosition(idx));
- ac.setLinearVelocity(idx, ac.getInvMass(idx) * ac.getForce(idx) * dt_ + ac.getLinearVelocity(idx));
+ ac.setPosition (idx, 0.5_r * ac.getInvMass(idx) * ac.getForce(idx) * dt_ * dt_ +
+ ac.getLinearVelocity(idx) * dt_ +
+ ac.getPosition(idx));
+ ac.setLinearVelocity(idx, ac.getInvMass(idx) * ac.getForce(idx) * dt_ +
+ ac.getLinearVelocity(idx));
{%- if bIntegrateRotation %}
const Vec3 wdot = math::transformMatrixRART(ac.getRotation(idx).getMatrix(),
ac.getInvInertiaBF(idx)) * ac.getTorque(idx);
// Calculating the rotation angle
- const Vec3 phi( ac.getAngularVelocity(idx) * dt_ + wdot * dt_ * dt_);
+ const Vec3 phi( 0.5_r * wdot * dt_ * dt_ +
+ ac.getAngularVelocity(idx) * dt_ );
// Calculating the new orientation
auto rotation = ac.getRotation(idx);
diff --git a/python/mesa_pd/templates/kernel/SemiImplicitEuler.templ.h b/python/mesa_pd/templates/kernel/SemiImplicitEuler.templ.h
new file mode 100644
index 0000000000000000000000000000000000000000..7045cbfe1ec8cc74f72d9554d16f10503745f01a
--- /dev/null
+++ b/python/mesa_pd/templates/kernel/SemiImplicitEuler.templ.h
@@ -0,0 +1,109 @@
+//======================================================================================================================
+//
+// 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 SemiImplicitEuler.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// 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 {
+
+/**
+ * Semi-implicit Euler integration for position and velocity.
+ *
+ * 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
+ *
+ * \pre All forces and torques acting on the particles have to be set.
+ * \post All forces and torques are reset to 0.
+ * \ingroup mesa_pd_kernel
+ */
+class SemiImplicitEuler
+{
+public:
+ explicit SemiImplicitEuler(const real_t dt) : dt_(dt) {}
+
+ template <typename Accessor>
+ void operator()(const size_t i, Accessor& ac) const;
+private:
+ real_t dt_ = real_t(0.0);
+};
+
+template <typename Accessor>
+inline void SemiImplicitEuler::operator()(const size_t idx,
+ Accessor& ac) const
+{
+ static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
+
+ if (!data::particle_flags::isSet( ac.getFlags(idx), data::particle_flags::FIXED))
+ {
+ ac.setLinearVelocity( idx, ac.getInvMass(idx) * ac.getForce(idx) * dt_ +
+ ac.getLinearVelocity(idx));
+ ac.setPosition ( idx, ac.getLinearVelocity(idx) * dt_ +
+ ac.getPosition(idx));
+
+ {%- if bIntegrateRotation %}
+ const Vec3 wdot = math::transformMatrixRART(ac.getRotation(idx).getMatrix(),
+ ac.getInvInertiaBF(idx)) * ac.getTorque(idx);
+
+ ac.setAngularVelocity(idx, wdot * dt_ +
+ ac.getAngularVelocity(idx));
+
+ // Calculating the rotation angle
+ const Vec3 phi( ac.getAngularVelocity(idx) * dt_ );
+
+ // Calculating the new orientation
+ auto rotation = ac.getRotation(idx);
+ rotation.rotate( phi );
+ ac.setRotation(idx, rotation);
+
+ {%- endif %}
+ }
+
+ ac.setForce (idx, Vec3(real_t(0), real_t(0), real_t(0)));
+ {%- if bIntegrateRotation %}
+ ac.setTorque(idx, Vec3(real_t(0), real_t(0), real_t(0)));
+ {%- endif %}
+}
+
+} //namespace kernel
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/python/mesa_pd/templates/kernel/VelocityVerlet.templ.h b/python/mesa_pd/templates/kernel/VelocityVerlet.templ.h
index 148ac8c8d96d71d4cb7a0af97805689236d46a35..98d97b628174fc04f082517336bc709766b62a19 100644
--- a/python/mesa_pd/templates/kernel/VelocityVerlet.templ.h
+++ b/python/mesa_pd/templates/kernel/VelocityVerlet.templ.h
@@ -41,6 +41,8 @@ namespace kernel {
* integration is only complete when both functions are called. The integration
* is symplectic.
*
+ * Wachs, A. Particle-scale computational approaches to model dry and saturated granular flows of non-Brownian, non-cohesive, and non-spherical rigid bodies. Acta Mech 230, 1919–1980 (2019). https://doi.org/10.1007/s00707-019-02389-9
+ *
* This kernel requires the following particle accessor interface
* \code
{%- for prop in interface %}
diff --git a/src/mesa_pd/kernel/ExplicitEuler.h b/src/mesa_pd/kernel/ExplicitEuler.h
index 8887219b19e206bf4e5b193b7090e3ae989886ce..639e3671a8712875796060e4cace0701ae6f74ff 100644
--- a/src/mesa_pd/kernel/ExplicitEuler.h
+++ b/src/mesa_pd/kernel/ExplicitEuler.h
@@ -34,7 +34,12 @@ namespace mesa_pd {
namespace kernel {
/**
- * Kernel which explicitly integrates all particles in time.
+ * Explicit Euler integration.
+ * Uses the 0.5at^2 extension for position integration.
+ * Boils down to using v_{t+0.5dt} for position integration.
+ * This version exhibits increased stability compared to standard explicit euler.
+ *
+ * Wachs, A. Particle-scale computational approaches to model dry and saturated granular flows of non-Brownian, non-cohesive, and non-spherical rigid bodies. Acta Mech 230, 1919–1980 (2019). https://doi.org/10.1007/s00707-019-02389-9
*
* This kernel requires the following particle accessor interface
* \code
@@ -87,13 +92,17 @@ inline void ExplicitEuler::operator()(const size_t idx,
if (!data::particle_flags::isSet( ac.getFlags(idx), data::particle_flags::FIXED))
{
- ac.setPosition (idx, ac.getInvMass(idx) * ac.getForce(idx) * dt_ * dt_ + ac.getLinearVelocity(idx) * dt_ + ac.getPosition(idx));
- ac.setLinearVelocity(idx, ac.getInvMass(idx) * ac.getForce(idx) * dt_ + ac.getLinearVelocity(idx));
+ ac.setPosition (idx, 0.5_r * ac.getInvMass(idx) * ac.getForce(idx) * dt_ * dt_ +
+ ac.getLinearVelocity(idx) * dt_ +
+ ac.getPosition(idx));
+ ac.setLinearVelocity(idx, ac.getInvMass(idx) * ac.getForce(idx) * dt_ +
+ ac.getLinearVelocity(idx));
const Vec3 wdot = math::transformMatrixRART(ac.getRotation(idx).getMatrix(),
ac.getInvInertiaBF(idx)) * ac.getTorque(idx);
// Calculating the rotation angle
- const Vec3 phi( ac.getAngularVelocity(idx) * dt_ + wdot * dt_ * dt_);
+ const Vec3 phi( 0.5_r * wdot * dt_ * dt_ +
+ ac.getAngularVelocity(idx) * dt_ );
// Calculating the new orientation
auto rotation = ac.getRotation(idx);
diff --git a/src/mesa_pd/kernel/SemiImplicitEuler.h b/src/mesa_pd/kernel/SemiImplicitEuler.h
new file mode 100644
index 0000000000000000000000000000000000000000..cf0df0734ea6bc7e5120809ef044c0e96f1e708b
--- /dev/null
+++ b/src/mesa_pd/kernel/SemiImplicitEuler.h
@@ -0,0 +1,115 @@
+//======================================================================================================================
+//
+// 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 SemiImplicitEuler.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// 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 {
+
+/**
+ * Semi-implicit Euler integration for position and velocity.
+ *
+ * This kernel requires the following particle accessor interface
+ * \code
+ * const walberla::mesa_pd::Vec3& getPosition(const size_t p_idx) const;
+ * void setPosition(const size_t p_idx, const walberla::mesa_pd::Vec3& v);
+ *
+ * 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::real_t& getInvMass(const size_t p_idx) const;
+ *
+ * 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::Rot3& getRotation(const size_t p_idx) const;
+ * void setRotation(const size_t p_idx, const walberla::mesa_pd::Rot3& v);
+ *
+ * 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::Mat3& getInvInertiaBF(const size_t p_idx) const;
+ *
+ * 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
+ *
+ * \pre All forces and torques acting on the particles have to be set.
+ * \post All forces and torques are reset to 0.
+ * \ingroup mesa_pd_kernel
+ */
+class SemiImplicitEuler
+{
+public:
+ explicit SemiImplicitEuler(const real_t dt) : dt_(dt) {}
+
+ template <typename Accessor>
+ void operator()(const size_t i, Accessor& ac) const;
+private:
+ real_t dt_ = real_t(0.0);
+};
+
+template <typename Accessor>
+inline void SemiImplicitEuler::operator()(const size_t idx,
+ Accessor& ac) const
+{
+ static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
+
+ if (!data::particle_flags::isSet( ac.getFlags(idx), data::particle_flags::FIXED))
+ {
+ ac.setLinearVelocity( idx, ac.getInvMass(idx) * ac.getForce(idx) * dt_ +
+ ac.getLinearVelocity(idx));
+ ac.setPosition ( idx, ac.getLinearVelocity(idx) * dt_ +
+ ac.getPosition(idx));
+ const Vec3 wdot = math::transformMatrixRART(ac.getRotation(idx).getMatrix(),
+ ac.getInvInertiaBF(idx)) * ac.getTorque(idx);
+
+ ac.setAngularVelocity(idx, wdot * dt_ +
+ ac.getAngularVelocity(idx));
+
+ // Calculating the rotation angle
+ const Vec3 phi( ac.getAngularVelocity(idx) * dt_ );
+
+ // Calculating the new orientation
+ auto rotation = ac.getRotation(idx);
+ rotation.rotate( phi );
+ ac.setRotation(idx, rotation);
+ }
+
+ ac.setForce (idx, Vec3(real_t(0), real_t(0), real_t(0)));
+ ac.setTorque(idx, Vec3(real_t(0), real_t(0), real_t(0)));
+}
+
+} //namespace kernel
+} //namespace mesa_pd
+} //namespace walberla
\ No newline at end of file
diff --git a/src/mesa_pd/kernel/VelocityVerlet.h b/src/mesa_pd/kernel/VelocityVerlet.h
index f029652c27e4c1e5e4ac6e81d27b083e59e3a7b6..69ed50dcb468111399c11a70534087af9255a69f 100644
--- a/src/mesa_pd/kernel/VelocityVerlet.h
+++ b/src/mesa_pd/kernel/VelocityVerlet.h
@@ -41,6 +41,8 @@ namespace kernel {
* integration is only complete when both functions are called. The integration
* is symplectic.
*
+ * Wachs, A. Particle-scale computational approaches to model dry and saturated granular flows of non-Brownian, non-cohesive, and non-spherical rigid bodies. Acta Mech 230, 1919–1980 (2019). https://doi.org/10.1007/s00707-019-02389-9
+ *
* This kernel requires the following particle accessor interface
* \code
* const walberla::mesa_pd::Vec3& getPosition(const size_t p_idx) const;
diff --git a/tests/mesa_pd/CMakeLists.txt b/tests/mesa_pd/CMakeLists.txt
index 3742a2ef2dcabdd68d50fe0c436fd0b21d360435..caf3ba2752b2630203a42de911480bb4969eda00 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_SemiImplicitEuler FILES kernel/SemiImplicitEuler.cpp DEPENDS core )
+waLBerla_execute_test( NAME MESA_PD_Kernel_SemiImplicitEuler )
+
waLBerla_compile_test( NAME MESA_PD_Kernel_SingleCast FILES kernel/SingleCast.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_Kernel_SingleCast )
@@ -197,4 +200,4 @@ if (WALBERLA_MESAPD_CONVEX_POLYHEDRON_AVAILABLE)
waLBerla_compile_test( NAME MESA_PD_COLLISIONDETECTION_ConvexPolyhedron_GJK_EPA FILES collision_detection/ConvexPolyhedron_GJK_EPA.cpp DEPENDS core mesa_pd mesh_common )
waLBerla_execute_test( NAME MESA_PD_COLLISIONDETECTION_ConvexPolyhedron_GJK_EPA )
-endif()
\ No newline at end of file
+endif()
diff --git a/tests/mesa_pd/DropTestAnalytic.cpp b/tests/mesa_pd/DropTestAnalytic.cpp
index 7c3c89ed0b3caec83e9bb699a1062a2e26ef44d9..3b257e23c0a3cac077300db2a1fe3bc79087b038 100644
--- a/tests/mesa_pd/DropTestAnalytic.cpp
+++ b/tests/mesa_pd/DropTestAnalytic.cpp
@@ -25,7 +25,7 @@
#include <mesa_pd/data/ShapeStorage.h>
#include <mesa_pd/kernel/DoubleCast.h>
-#include <mesa_pd/kernel/ExplicitEuler.h>
+#include <mesa_pd/kernel/SemiImplicitEuler.h>
#include <mesa_pd/kernel/ParticleSelector.h>
#include <mesa_pd/kernel/SpringDashpot.h>
@@ -102,7 +102,7 @@ int main( int argc, char ** argv )
real_t dt = real_t(0.00001);
// Init kernels
- kernel::ExplicitEuler explicitEuler( dt );
+ kernel::SemiImplicitEuler implEuler( dt );
kernel::SpringDashpot dem(1);
auto meff = real_t(1.0) / ss->shapes[sp->getShapeID()]->getInvMass();
dem.setParametersFromCOR(0,0,real_t(0.9), dt * real_t(20), meff);
@@ -152,7 +152,7 @@ int main( int argc, char ** argv )
ps->forEachParticle(false,
kernel::SelectLocal(),
accessor,
- explicitEuler,
+ implEuler,
accessor);
// if(i%1 == 0)
diff --git a/tests/mesa_pd/DropTestGeneral.cpp b/tests/mesa_pd/DropTestGeneral.cpp
index 3d668694c7f3429e12354e561a24389ca4e87443..0ab23c4fa42b9265200c11045c6b79d4773dd18d 100644
--- a/tests/mesa_pd/DropTestGeneral.cpp
+++ b/tests/mesa_pd/DropTestGeneral.cpp
@@ -25,16 +25,14 @@
#include <mesa_pd/data/ShapeStorage.h>
#include <mesa_pd/kernel/DoubleCast.h>
-#include <mesa_pd/kernel/ExplicitEuler.h>
+#include <mesa_pd/kernel/SemiImplicitEuler.h>
#include <mesa_pd/kernel/ParticleSelector.h>
#include <mesa_pd/kernel/SpringDashpot.h>
#include <core/Abort.h>
#include <core/Environment.h>
#include <core/logging/Logging.h>
-#include <core/waLBerlaBuildInfo.h>
-#include <functional>
#include <memory>
#include <string>
#include <type_traits>
@@ -122,7 +120,7 @@ int main( int argc, char ** argv )
real_t dt = real_t(0.00001);
// Init kernels
- kernel::ExplicitEuler explicitEuler( dt );
+ kernel::SemiImplicitEuler implEuler( dt );
kernel::SpringDashpot dem(1);
auto meff = real_t(1.0) / ss->shapes[sp->getShapeID()]->getInvMass();
dem.setParametersFromCOR(0,0,real_t(0.9), dt * real_t(20), meff);
@@ -172,7 +170,7 @@ int main( int argc, char ** argv )
ps->forEachParticle(false,
kernel::SelectLocal(),
accessor,
- explicitEuler,
+ implEuler,
accessor);
// if(i%1 == 0)
diff --git a/tests/mesa_pd/kernel/CoefficientOfRestitutionLSD.cpp b/tests/mesa_pd/kernel/CoefficientOfRestitutionLSD.cpp
index 7e316264f934e5088e33028a4dc9fe96714f440b..ac16d145e73e75ccfa02cfc151aa09b046766f32 100644
--- a/tests/mesa_pd/kernel/CoefficientOfRestitutionLSD.cpp
+++ b/tests/mesa_pd/kernel/CoefficientOfRestitutionLSD.cpp
@@ -13,7 +13,7 @@
// 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 DEMIntegratorAccuracy.cpp
+//! \file CoefficientOfRestitutionLSD.cpp
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
@@ -26,7 +26,7 @@
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/VelocityVerlet.h"
-#include "mesa_pd/kernel/ExplicitEuler.h"
+#include "mesa_pd/kernel/SemiImplicitEuler.h"
#include "mesa_pd/kernel/LinearSpringDashpot.h"
#include "core/Environment.h"
@@ -136,7 +136,7 @@ int main( int argc, char** argv )
collision_detection::AnalyticContactDetection acd;
kernel::DoubleCast double_cast;
- kernel::ExplicitEuler explEuler(dt);
+ kernel::SemiImplicitEuler implEuler(dt);
kernel::VelocityVerletPreForceUpdate vvPreForce( dt );
kernel::VelocityVerletPostForceUpdate vvPostForce( dt );
@@ -160,7 +160,7 @@ int main( int argc, char** argv )
auto force = accessor->getForce(0);
if(useVelocityVerlet) vvPostForce(0,*accessor);
- else explEuler(0, *accessor);
+ else implEuler(0, *accessor);
WALBERLA_LOG_INFO(steps << ": penetration = " << acd.getPenetrationDepth() << " || vel = " << accessor->getLinearVelocity(0)[2] << " || force = " << force[2]);
diff --git a/tests/mesa_pd/kernel/CoefficientOfRestitutionNLSD.cpp b/tests/mesa_pd/kernel/CoefficientOfRestitutionNLSD.cpp
index c422a70b497527a3ff6ef5aaf4181168a7ec9ddb..c2bfbcfdb7160555bc8945a8400b7a06fdf248cb 100644
--- a/tests/mesa_pd/kernel/CoefficientOfRestitutionNLSD.cpp
+++ b/tests/mesa_pd/kernel/CoefficientOfRestitutionNLSD.cpp
@@ -13,7 +13,7 @@
// 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 DEMIntegratorAccuracy.cpp
+//! \file CoefficientOfRestitutionNLSD.cpp
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
@@ -26,7 +26,7 @@
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/VelocityVerlet.h"
-#include "mesa_pd/kernel/ExplicitEuler.h"
+#include "mesa_pd/kernel/SemiImplicitEuler.h"
#include "mesa_pd/kernel/NonLinearSpringDashpot.h"
#include "mesa_pd/mpi/ReduceContactHistory.h"
@@ -137,7 +137,7 @@ int main( int argc, char** argv )
collision_detection::AnalyticContactDetection acd;
kernel::DoubleCast double_cast;
- kernel::ExplicitEuler explEuler(dt);
+ kernel::SemiImplicitEuler implEuler(dt);
kernel::VelocityVerletPreForceUpdate vvPreForce( dt );
kernel::VelocityVerletPostForceUpdate vvPostForce( dt );
@@ -162,7 +162,7 @@ int main( int argc, char** argv )
auto force = accessor->getForce(0);
if(useVelocityVerlet) vvPostForce(0,*accessor);
- else explEuler(0, *accessor);
+ else implEuler(0, *accessor);
WALBERLA_LOG_INFO(steps << ": penetration = " << acd.getPenetrationDepth() << " || vel = " << accessor->getLinearVelocity(0)[2] << " || force = " << force[2]);
diff --git a/tests/mesa_pd/kernel/CoefficientOfRestitutionSD.cpp b/tests/mesa_pd/kernel/CoefficientOfRestitutionSD.cpp
index 300e0eb916d6f19330ccc873733387fbd47981e3..910f5da3d3d68029778879db3c5bd8b641c8dbcb 100644
--- a/tests/mesa_pd/kernel/CoefficientOfRestitutionSD.cpp
+++ b/tests/mesa_pd/kernel/CoefficientOfRestitutionSD.cpp
@@ -13,7 +13,7 @@
// 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 DEMIntegratorAccuracy.cpp
+//! \file CoefficientOfRestitutionSD.cpp
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
@@ -26,7 +26,7 @@
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/VelocityVerlet.h"
-#include "mesa_pd/kernel/ExplicitEuler.h"
+#include "mesa_pd/kernel/SemiImplicitEuler.h"
#include "mesa_pd/kernel/SpringDashpot.h"
#include "core/Environment.h"
@@ -136,7 +136,7 @@ int main( int argc, char** argv )
collision_detection::AnalyticContactDetection acd;
kernel::DoubleCast double_cast;
- kernel::ExplicitEuler explEuler(dt);
+ kernel::SemiImplicitEuler implEuler(dt);
kernel::VelocityVerletPreForceUpdate vvPreForce( dt );
kernel::VelocityVerletPostForceUpdate vvPostForce( dt );
@@ -160,7 +160,7 @@ int main( int argc, char** argv )
auto force = accessor->getForce(0);
if(useVelocityVerlet) vvPostForce(0,*accessor);
- else explEuler(0, *accessor);
+ else implEuler(0, *accessor);
WALBERLA_LOG_INFO(steps << ": penetration = " << acd.getPenetrationDepth() << " || vel = " << accessor->getLinearVelocity(0)[2] << " || force = " << force[2]);
diff --git a/tests/mesa_pd/kernel/ExplicitEuler.cpp b/tests/mesa_pd/kernel/ExplicitEuler.cpp
index 1685ff1aed29b315bfd36be07ef45b0465b44908..17d55505b02fe16895825eca3fa2c1588eb98c71 100644
--- a/tests/mesa_pd/kernel/ExplicitEuler.cpp
+++ b/tests/mesa_pd/kernel/ExplicitEuler.cpp
@@ -66,8 +66,10 @@ int main( int argc, char ** argv )
accessor.setAngularVelocity( 0, angVel);
accessor.setForce( 0, force);
accessor.setTorque( 0, torque);
- accessor.setInvMass( 0, real_t(1.23456));
- accessor.setInvInertiaBF( 0, Mat3(real_t(1.23456), real_t(0), real_t(0), real_t(0), real_t(1.23456), real_t(0), real_t(0), real_t(0), real_t(1.23456)));
+ accessor.setInvMass( 0, 1.23456_r);
+ accessor.setInvInertiaBF( 0, Mat3(1.23456_r, 0_r, 0_r,
+ 0_r, 1.23456_r, 0_r,
+ 0_r, 0_r, 1.23456_r ));
//init kernels
const real_t dt = real_t(1);
@@ -83,12 +85,16 @@ int main( int argc, char ** argv )
WALBERLA_CHECK_FLOAT_EQUAL(accessor.getTorque(0), Vec3(0));
//check velocity
- WALBERLA_CHECK_FLOAT_EQUAL(accessor.getLinearVelocity(0), force * accessor.getInvMass(0) * dt + linVel);
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getLinearVelocity(0), force * accessor.getInvMass(0) * dt +
+ linVel);
WALBERLA_CHECK_FLOAT_EQUAL(accessor.getAngularVelocity(0), wdot * dt + angVel);
//check position
- WALBERLA_CHECK_FLOAT_EQUAL(accessor.getPosition(0), linVel * dt + force * accessor.getInvMass(0) * dt * dt);
- WALBERLA_CHECK_FLOAT_EQUAL(accessor.getRotation(0).getQuaternion(), Quat( (wdot * dt + angVel).getNormalized(), (wdot * dt + angVel).length() * dt ));
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getPosition(0), linVel * dt +
+ 0.5_r * force * accessor.getInvMass(0) * dt * dt);
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getRotation(0).getQuaternion(),
+ Quat( (wdot * dt + angVel).getNormalized(),
+ (0.5_r * wdot * dt + angVel).length() * dt ));
accessor.setPosition( 0, Vec3(0,0,0));
accessor.setRotation( 0, Rot3(Quat()));
@@ -96,8 +102,10 @@ int main( int argc, char ** argv )
accessor.setAngularVelocity( 0, angVel);
accessor.setForce( 0, force);
accessor.setTorque( 0, torque);
- accessor.setInvMass( 0, real_t(1.23456));
- accessor.setInvInertiaBF( 0, Mat3(real_t(1.23456), real_t(0), real_t(0), real_t(0), real_t(1.23456), real_t(0), real_t(0), real_t(0), real_t(1.23456)));
+ accessor.setInvMass( 0, 1.23456_r);
+ accessor.setInvInertiaBF( 0, Mat3(1.23456_r, 0_r, 0_r,
+ 0_r, 1.23456_r, 0_r,
+ 0_r, 0_r, 1.23456_r ));
data::particle_flags::set( accessor.getFlagsRef(0), data::particle_flags::FIXED );
integrator(0, accessor);
diff --git a/tests/mesa_pd/kernel/GenerateLinkedCells.cpp b/tests/mesa_pd/kernel/GenerateLinkedCells.cpp
index e7a60f36ed57cfc8a05fbda3c14d884b291f73c6..d4c2f30c678c7bb03ed4bc565c3a163fee6a531a 100644
--- a/tests/mesa_pd/kernel/GenerateLinkedCells.cpp
+++ b/tests/mesa_pd/kernel/GenerateLinkedCells.cpp
@@ -18,7 +18,7 @@
//
//======================================================================================================================
-#include <mesa_pd/kernel/ExplicitEuler.h>
+#include <mesa_pd/kernel/SemiImplicitEuler.h>
#include <mesa_pd/kernel/ForceLJ.h>
#include <mesa_pd/kernel/InsertParticleIntoLinkedCells.h>
#include <mesa_pd/kernel/ParticleSelector.h>
@@ -72,7 +72,7 @@ int main( int argc, char ** argv )
//init kernels
kernel::InsertParticleIntoLinkedCells ipilc;
kernel::ForceLJ lj(1);
- kernel::ExplicitEuler integrator( real_t(0.01) );
+ kernel::SemiImplicitEuler integrator( real_t(0.01) );
//timeloop
for (auto timestep = 0; timestep < 100; ++timestep)
diff --git a/tests/mesa_pd/kernel/IntegratorAccuracy.cpp b/tests/mesa_pd/kernel/IntegratorAccuracy.cpp
index 695bae5b9662299ceb77583d726fab3e012e4b62..e97c47106b1cd741e0652654601dd986f3d0cb1b 100644
--- a/tests/mesa_pd/kernel/IntegratorAccuracy.cpp
+++ b/tests/mesa_pd/kernel/IntegratorAccuracy.cpp
@@ -15,170 +15,252 @@
//
//! \file IntegratorAccuracy.cpp
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
//
//======================================================================================================================
-#include "mesa_pd/data/ParticleAccessorWithShape.h"
-#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/ParticleAccessor.h"
-#include "mesa_pd/kernel/ParticleSelector.h"
-#include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/kernel/ExplicitEuler.h"
+#include "mesa_pd/kernel/SemiImplicitEuler.h"
+#include "mesa_pd/kernel/VelocityVerlet.h"
#include "core/Environment.h"
-#include "core/logging/Logging.h"
#include "core/math/all.h"
#include <iostream>
-namespace integrator_accuracy
-{
-
-using namespace walberla;
+namespace walberla {
using namespace walberla::mesa_pd;
-
-Vec3 getForce(const Vec3& pos, real_t k)
+class SingleParticleAccessorWithShape : public data::SingleParticleAccessor
{
- return -k * pos;
-}
+public:
+ void setInvMass(const size_t /*p_idx*/, const real_t &val)
+ { invMass_ = val; }
-real_t analyticalTrajectory(real_t amplitude, real_t timeStep, real_t omega, real_t phase)
-{
- return amplitude * std::cos(omega * timeStep + phase);
-}
+ const auto &getInvMass(const size_t /*p_idx*/) const
+ { return invMass_; }
-real_t analyticalVelocity(real_t amplitude, real_t timeStep, real_t omega, real_t phase)
-{
- return -amplitude * omega * std::sin(omega * timeStep + phase);
-}
+ void setInvInertiaBF(const size_t /*p_idx*/, const Mat3 &val)
+ { invInertiaBF_ = val; }
+ const auto &getInvInertiaBF(const size_t /*p_idx*/) const
+ { return invInertiaBF_; }
-/*
- * Simulates a harmonic oscillator to test the accuracy of the integrators.
- * The error of the maximum position and the maximum velocity is compared against the analytical values.
- * Via command line arguments, the simulation can be adapted.
- * Currently tested integrators:
- * - explicit Euler (default)
- * - velocity verlet (-useVV)
- */
-int main( int argc, char ** argv )
-{
- mpi::Environment env(argc, argv);
- WALBERLA_UNUSED(env);
- mpi::MPIManager::instance()->useWorldComm();
+private:
+ real_t invMass_;
+ Mat3 invInertiaBF_;
+};
- real_t amplitude = real_t(1.5);
- real_t k = real_t(0.1);
- real_t mass = real_t(0.9);
- real_t dt = real_t(0.2);
- bool useVelocityVerlet = false;
- real_t phaseFraction = real_t(0);
- real_t periods = real_t(1.1);
+struct Oscillator
+{
+ real_t amplitude = 1.5_r;
+ real_t k = 0.1_r;
+ real_t damping = 0_r;
+ real_t mass = 0.9_r;
+ real_t dt = 0.2_r;
+ real_t phaseFraction = 0_r;
+ real_t periods = 10_r;
+ real_t phase = phaseFraction * math::pi;
+ real_t dampingRatio = damping / (2_r * std::sqrt(mass * k));
+ real_t omega = std::sqrt(k / mass) * std::sqrt(1_r - dampingRatio * dampingRatio);
+ real_t decay = std::sqrt(k / mass) * dampingRatio;
+ real_t durationOnePeriod = 2_r * math::pi / omega;
+ uint_t timeSteps = uint_c(periods * durationOnePeriod / dt);
- for( int i = 1; i < argc; ++i )
+ void update()
{
- if( std::strcmp( argv[i], "--useVV" ) == 0 ) { useVelocityVerlet = true; continue; }
- if( std::strcmp( argv[i], "--dt" ) == 0 ) { dt = real_c(std::atof( argv[++i] )); continue; }
- if( std::strcmp( argv[i], "--amplitude" ) == 0 ) { amplitude = real_c(std::atof( argv[++i] )); continue; }
- if( std::strcmp( argv[i], "--k" ) == 0 ) { k = real_c(std::atof( argv[++i] )); continue; }
- if( std::strcmp( argv[i], "--mass" ) == 0 ) { mass = real_c(std::atof( argv[++i] )); continue; }
- if( std::strcmp( argv[i], "--phaseFraction" ) == 0 ) { phaseFraction = real_c(std::atof( argv[++i] )); continue; }
- if( std::strcmp( argv[i], "--periods" ) == 0 ) { periods = real_c(std::atof( argv[++i] )); continue; }
- WALBERLA_ABORT("Unrecognized command line argument found: " << argv[i]);
+ phase = phaseFraction * math::pi;
+ dampingRatio = damping / (2_r * std::sqrt(mass * k));
+ omega = std::sqrt(k / mass) * std::sqrt(1_r - dampingRatio * dampingRatio);
+ decay = std::sqrt(k / mass) * dampingRatio;
+ durationOnePeriod = 2_r * math::pi / omega;
+ timeSteps = uint_c(periods * durationOnePeriod / dt);
}
+ Vec3 getForce(const Vec3 &pos, const Vec3 &vel) const
+ {
+ return -k * pos - damping * vel;
+ }
- real_t phase = phaseFraction * math::pi;
- real_t omega = std::sqrt(k / mass);
- real_t durationOnePeriod = real_t(2) * math::pi / omega;
- uint_t timeSteps = uint_c(periods * durationOnePeriod / dt);
+ real_t getEnergy(const real_t &pos, const real_t &vel) const
+ {
+ return 0.5_r * mass * vel * vel + 0.5_r * k * pos * pos;
+ }
- WALBERLA_LOG_INFO("omega = " << omega << ", T = " << real_t(2) * math::pi / omega << ", time steps = " << timeSteps << ", phase = " << phase << ", periods = " << periods);
+ real_t analyticalPos(real_t t) const
+ {
+ return amplitude * std::exp(-decay * t) * std::cos(omega * t + phase);
+ }
- //initialize particle
- const auto pos = Vec3(0,0,analyticalTrajectory(amplitude, real_t(0), omega, phase));
- const auto linVel = Vec3(0,0,analyticalVelocity(amplitude, real_t(0), omega, phase));
+ real_t analyticalVel(real_t t) const
+ {
+ return -decay * amplitude * std::exp(-decay * t) * std::cos(omega * t + phase)
+ -amplitude * std::exp(-decay * t) * omega * std::sin(omega * t + phase);
+ }
+};
- WALBERLA_LOG_INFO("Initial pos = " << pos << ", vel = " << linVel);
+struct ExplicitEuler
+{
+ ExplicitEuler(real_t dt) : integrator(dt) {}
+ void operator()(SingleParticleAccessorWithShape& particle,
+ const Oscillator& osc)
+ {
+ particle.setForce(0, osc.getForce(particle.getPosition(0),
+ particle.getLinearVelocity(0)));
+ integrator(0, particle);
+ }
+ kernel::ExplicitEuler integrator;
+};
- //init data structures
- auto ps = std::make_shared<data::ParticleStorage>(1);
- auto ss = std::make_shared<data::ShapeStorage>();
- auto smallSphere = ss->create<data::Sphere>( real_t(1) );
- ss->shapes[smallSphere]->updateMassAndInertia(real_t(2707));
- data::ParticleAccessorWithShape accessor(ps, ss);
-
- data::Particle&& p = *ps->create();
- p.setPosition(pos);
- p.setLinearVelocity(linVel);
- p.setForce(getForce(pos, k));
- p.setOldForce(getForce(pos, k));
- p.setInvMass(real_t(1) / mass);
- p.setShapeID(smallSphere);
-
- // velocity verlet
- kernel::VelocityVerletPreForceUpdate preForce( dt );
- kernel::VelocityVerletPostForceUpdate postForce( dt );
+struct SemiImplicitEuler
+{
+ SemiImplicitEuler(real_t dt) : integrator(dt) {}
+ void operator()(SingleParticleAccessorWithShape& particle,
+ const Oscillator& osc)
+ {
+ particle.setForce(0, osc.getForce(particle.getPosition(0),
+ particle.getLinearVelocity(0)));
+ integrator(0, particle);
+ }
+ kernel::SemiImplicitEuler integrator;
+};
- // explicit euler
- kernel::ExplicitEuler explEuler( dt );
-
- real_t maxVel = 0.;
- real_t maxRefVel = 0.;
- real_t maxPos = 0.;
- real_t maxRefPos = 0.;
-
- for (auto i = uint_t(1); i <= timeSteps; ++i)
+struct VelocityVerlet
+{
+ VelocityVerlet(real_t dt) : preVV(dt), postVV(dt) {}
+ void operator()(SingleParticleAccessorWithShape& particle,
+ const Oscillator& osc)
{
+ preVV(0, particle);
+ particle.setForce(0, osc.getForce(particle.getPosition(0),
+ particle.getLinearVelocity(0)));
+ postVV(0, particle);
+ }
+ kernel::VelocityVerletPreForceUpdate preVV;
+ kernel::VelocityVerletPostForceUpdate postVV;
+};
- if( useVelocityVerlet ) ps->forEachParticle(false, kernel::SelectAll(), accessor, preForce, accessor);
- p.setForce( getForce(p.getPosition(), k) );
- auto force = p.getForce();
+struct AccuracyResult
+{
+ real_t maxPosDeviation;
+ real_t maxVelDeviation;
+ real_t maxEneDeviation;
+};
- if( useVelocityVerlet ) ps->forEachParticle(false, kernel::SelectAll(), accessor, postForce, accessor);
- else ps->forEachParticle(false, kernel::SelectAll(), accessor, explEuler, accessor);
+template <typename Integrator>
+AccuracyResult checkIntegrator(const Oscillator& osc)
+{
+ //init data structures
+ SingleParticleAccessorWithShape particle;
- real_t refPos = analyticalTrajectory(amplitude, real_c(i) * dt, omega, phase);
- real_t refVel = analyticalVelocity(amplitude, real_c(i) * dt, omega, phase);
+ //first dummy argument is needed to fulfill accessor interface
+ particle.setPosition(0, Vec3(0, 0, osc.analyticalPos(0_r)));
+ particle.setLinearVelocity(0, Vec3(0, 0, osc.analyticalVel(0_r)));
+ particle.setInvMass(0, 1_r / osc.mass);
+ particle.setForce(0, osc.getForce(Vec3(0, 0, osc.analyticalPos(real_t(0))),
+ Vec3(0, 0, osc.analyticalVel(real_t(0)))));
+ particle.setOldForce(0, osc.getForce(Vec3(0, 0, osc.analyticalPos(-osc.dt)),
+ Vec3(0, 0, osc.analyticalVel(-osc.dt))));
- WALBERLA_LOG_INFO(i << ": pos = " << p.getPosition()[2] << " " << refPos
- << " || vel = " << p.getLinearVelocity()[2] << " " << refVel
- << " || force = " << force[2] );
+ // explicit euler
+ Integrator integrator(osc.dt);
- maxPos = std::max(maxPos, std::abs(p.getPosition()[2]));
- maxRefPos = std::max(maxRefPos, std::abs(refPos));
+ real_t maxPosDeviation = 0_r;
+ real_t maxVelDeviation = 0_r;
+ real_t maxEneDeviation = 0_r;
- maxVel = std::max(maxVel, std::abs(p.getLinearVelocity()[2]));
- maxRefVel = std::max(maxRefVel, std::abs(refVel));
+ for (auto i = uint_t(0); i <= osc.timeSteps; ++i)
+ {
+ real_t refPos = osc.analyticalPos(real_c(i) * osc.dt);
+ real_t refVel = osc.analyticalVel(real_c(i) * osc.dt);
+ real_t refEne = osc.getEnergy(refPos, refVel);
+
+ maxPosDeviation = std::max(maxPosDeviation, std::abs(particle.getPosition(0)[2] - refPos));
+ maxVelDeviation = std::max(maxVelDeviation, std::abs(particle.getLinearVelocity(0)[2] - refVel));
+ maxEneDeviation = std::max(maxEneDeviation, std::abs(osc.getEnergy(particle.getPosition(0)[2], particle.getLinearVelocity(0)[2]) - refEne));
+
+ integrator(particle, osc);
}
- real_t relativePositionError = ( maxPos - maxRefPos ) / maxRefPos;
- WALBERLA_LOG_INFO("error in position = " << relativePositionError * 100. << "%");
+ return {maxPosDeviation, maxVelDeviation, maxEneDeviation};
+}
- real_t relativeVelocityError = ( maxVel - maxRefVel ) / maxRefVel;
- WALBERLA_LOG_INFO("error in velocity = " << relativeVelocityError * 100. << "%");
+int main(int argc, char **argv)
+{
+ mpi::Environment env(argc, argv);
+ WALBERLA_UNUSED(env);
+ mpi::MPIManager::instance()->useWorldComm();
- if( useVelocityVerlet )
- {
- WALBERLA_CHECK_LESS(relativePositionError, real_t(0.01), "Error in position too large!");
- WALBERLA_CHECK_LESS(relativeVelocityError, real_t(0.01), "Error in velocity too large!");
- }
- else
+ if (std::is_same<real_t, float>::value)
{
- WALBERLA_CHECK_LESS(relativePositionError, real_t(0.11), "Error in position too large!");
- WALBERLA_CHECK_LESS(relativeVelocityError, real_t(0.10), "Error in velocity too large!");
+ WALBERLA_LOG_WARNING("waLBerla build in sp mode: skipping test due to low precision");
+ return EXIT_SUCCESS;
}
+ Oscillator osc;
+
+ AccuracyResult res;
+ osc.dt = 0.1_r;
+ osc.update();
+ res = walberla::checkIntegrator<walberla::ExplicitEuler>(osc);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxPosDeviation, 1.03068993874562120e+00_r);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxVelDeviation, 3.33225581358688350e-01_r);
+ res = walberla::checkIntegrator<walberla::SemiImplicitEuler>(osc);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxPosDeviation, 2.92576360173544339e-02_r);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxVelDeviation, 1.45120298258364505e-03_r);
+ res = walberla::checkIntegrator<walberla::VelocityVerlet>(osc);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxPosDeviation, 4.24116598691555435e-03_r);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxVelDeviation, 1.45059589844465098e-03_r);
+
+ osc.dt = 0.2_r;
+ osc.update();
+ res = walberla::checkIntegrator<walberla::ExplicitEuler>(osc);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxPosDeviation, 2.76387000209972467e+00_r);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxVelDeviation, 8.88433661185269896e-01_r);
+ res = walberla::checkIntegrator<walberla::SemiImplicitEuler>(osc);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxPosDeviation, 6.70464626869100577e-02_r);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxVelDeviation, 5.81009940233766925e-03_r);
+ res = walberla::checkIntegrator<walberla::VelocityVerlet>(osc);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxPosDeviation, 1.69147419522671719e-02_r);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxVelDeviation, 5.78432113979018836e-03_r);
+
+ osc.dt = 0.4_r;
+ osc.update();
+ res = walberla::checkIntegrator<walberla::ExplicitEuler>(osc);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxPosDeviation, 1.04680753378045406e+01_r);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxVelDeviation, 3.34470580215144420e+00_r);
+ res = walberla::checkIntegrator<walberla::SemiImplicitEuler>(osc);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxPosDeviation, 1.68291142727994780e-01_r);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxVelDeviation, 2.33193930919295134e-02_r);
+ res = walberla::checkIntegrator<walberla::VelocityVerlet>(osc);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxPosDeviation, 6.71909796751584687e-02_r);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxVelDeviation, 2.29906986019860066e-02_r);
+
+ //check energy conservation
+ osc.dt = 0.4_r;
+ osc.periods = 1000;
+ osc.update();
+ //res = walberla::checkIntegrator<walberla::ExplicitEuler>(osc);
+ //explicit euler is not symplectic
+ res = walberla::checkIntegrator<walberla::SemiImplicitEuler>(osc);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxEneDeviation, 8.03571427989904774e-03_r);
+ res = walberla::checkIntegrator<walberla::VelocityVerlet>(osc);
+ WALBERLA_CHECK_FLOAT_EQUAL( res.maxEneDeviation, 4.99960610503419334e-04_r);
return EXIT_SUCCESS;
}
-} //namespace integrator_accuracy
+} //namespace walberla
+
-int main( int argc, char ** argv )
+/*
+ * Simulates a harmonic oscillator to test the accuracy of the integrators.
+ */
+int main(int argc, char **argv)
{
- return integrator_accuracy::main(argc, argv);
+ return walberla::main(argc, argv);
}
diff --git a/tests/mesa_pd/kernel/LinearSpringDashpot.cpp b/tests/mesa_pd/kernel/LinearSpringDashpot.cpp
index 18beb0732157060883c716146297fa803e3b82cb..64f099e2679602570afa4f3b6c6922751038d0d2 100644
--- a/tests/mesa_pd/kernel/LinearSpringDashpot.cpp
+++ b/tests/mesa_pd/kernel/LinearSpringDashpot.cpp
@@ -26,7 +26,7 @@
#include "mesa_pd/data/ShapeStorage.h"
#include "mesa_pd/kernel/DoubleCast.h"
-#include "mesa_pd/kernel/ExplicitEuler.h"
+#include "mesa_pd/kernel/SemiImplicitEuler.h"
#include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/kernel/LinearSpringDashpot.h"
#include "mesa_pd/mpi/ReduceContactHistory.h"
@@ -152,7 +152,7 @@ int main( int argc, char ** argv )
kernel::VelocityVerletPostForceUpdate vvPostForce( dt );
// explicit euler
- kernel::ExplicitEuler explEuler( dt );
+ kernel::SemiImplicitEuler implEuler( dt );
// collision response
collision_detection::AnalyticContactDetection acd;
@@ -188,7 +188,7 @@ int main( int argc, char ** argv )
rch(*ps);
if(useVelocityVerlet) vvPostForce(0,*accessor);
- else explEuler(0, *accessor);
+ else implEuler(0, *accessor);
++steps;
} while (double_cast(0, 1, *accessor, acd, *accessor ) || p.getLinearVelocity()[2] < 0);
diff --git a/tests/mesa_pd/kernel/SemiImplicitEuler.cpp b/tests/mesa_pd/kernel/SemiImplicitEuler.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2d3a22f84fac0c4ecf720221dfb70c95466c3350
--- /dev/null
+++ b/tests/mesa_pd/kernel/SemiImplicitEuler.cpp
@@ -0,0 +1,133 @@
+//======================================================================================================================
+//
+// 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 SemiImplicitEuler.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/ParticleAccessor.h>
+
+#include <mesa_pd/kernel/SemiImplicitEuler.h>
+
+#include <core/Environment.h>
+#include <core/logging/Logging.h>
+
+#include <iostream>
+
+namespace walberla {
+
+using namespace walberla::mesa_pd;
+
+class SingleParticleAccessor : public data::SingleParticleAccessor
+{
+public:
+ const walberla::real_t& getInvMass(const size_t /*p_idx*/) const {return invMass_;}
+ void setInvMass(const size_t /*p_idx*/, const walberla::real_t& v) { invMass_ = v;}
+ const walberla::mesa_pd::Mat3& getInvInertiaBF(const size_t /*p_idx*/) const {return invInertiaBF_;}
+ void setInvInertiaBF(const size_t /*p_idx*/, const walberla::mesa_pd::Mat3& v) { invInertiaBF_ = v;}
+
+ walberla::real_t invMass_;
+ walberla::mesa_pd::Mat3 invInertiaBF_;
+};
+
+int main( int argc, char ** argv )
+{
+ Environment env(argc, argv);
+ WALBERLA_UNUSED(env);
+ mpi::MPIManager::instance()->useWorldComm();
+
+ //init data structures
+ SingleParticleAccessor accessor;
+
+ //initialize particle
+ const auto linVel = Vec3(1,2,3);
+ const auto angVel = Vec3(1,2,3);
+
+ const auto force = Vec3(1,2,3);
+ const auto torque = Vec3(1,2,3);
+
+ accessor.setPosition( 0, Vec3(0,0,0));
+ accessor.setRotation( 0, Rot3(Quat()));
+ accessor.setLinearVelocity( 0, linVel);
+ accessor.setAngularVelocity( 0, angVel);
+ accessor.setForce( 0, force);
+ accessor.setTorque( 0, torque);
+ accessor.setInvMass( 0, 1.23456_r);
+ accessor.setInvInertiaBF( 0, Mat3(1.23456_r, 0_r, 0_r,
+ 0_r, 1.23456_r, 0_r,
+ 0_r, 0_r, 1.23456_r ));
+
+ //init kernels
+ const real_t dt = real_t(1);
+ kernel::SemiImplicitEuler integrator( dt );
+
+ integrator(0, accessor);
+
+ const auto& R = accessor.getRotation(0).getMatrix();
+ const auto wdot = R * accessor.getInvInertiaBF(0) * R.getTranspose() * torque;
+
+ //check force
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getForce(0), Vec3(0));
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getTorque(0), Vec3(0));
+
+ //check velocity
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getLinearVelocity(0), force * accessor.getInvMass(0) * dt +
+ linVel);
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getAngularVelocity(0), wdot * dt + angVel);
+
+ //check position
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getPosition(0), linVel * dt +
+ force * accessor.getInvMass(0) * dt * dt);
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getRotation(0).getQuaternion(),
+ Quat( (wdot * dt + angVel).getNormalized(),
+ (wdot * dt + angVel).length() * dt ));
+
+ accessor.setPosition( 0, Vec3(0,0,0));
+ accessor.setRotation( 0, Rot3(Quat()));
+ accessor.setLinearVelocity( 0, linVel);
+ accessor.setAngularVelocity( 0, angVel);
+ accessor.setForce( 0, force);
+ accessor.setTorque( 0, torque);
+ accessor.setInvMass( 0, 1.23456_r);
+ accessor.setInvInertiaBF( 0, Mat3(1.23456_r, 0_r, 0_r,
+ 0_r, 1.23456_r, 0_r,
+ 0_r, 0_r, 1.23456_r ));
+ data::particle_flags::set( accessor.getFlagsRef(0), data::particle_flags::FIXED );
+
+ integrator(0, accessor);
+
+ //check force
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getForce(0), Vec3(0));
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getTorque(0), Vec3(0));
+
+ //check velocity
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getLinearVelocity(0), linVel);
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getAngularVelocity(0), angVel);
+
+ //check position
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getPosition(0), Vec3(0,0,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getRotation(0).getQuaternion(), Quat());
+
+ return EXIT_SUCCESS;
+}
+
+} //namespace walberla
+
+int main( int argc, char ** argv )
+{
+ return walberla::main(argc, argv);
+}
diff --git a/tests/mesa_pd/kernel/SpherePile.cpp b/tests/mesa_pd/kernel/SpherePile.cpp
index 67e1087a749cb50e0c29f77b40ac197a5598f0a4..3f71e929724fb98a3cdc6cb73095a45f6d79711a 100644
--- a/tests/mesa_pd/kernel/SpherePile.cpp
+++ b/tests/mesa_pd/kernel/SpherePile.cpp
@@ -26,7 +26,7 @@
#include "mesa_pd/data/ShapeStorage.h"
#include "mesa_pd/kernel/DoubleCast.h"
-#include "mesa_pd/kernel/ExplicitEuler.h"
+#include "mesa_pd/kernel/SemiImplicitEuler.h"
#include "mesa_pd/kernel/SpringDashpot.h"
#include "mesa_pd/kernel/SpringDashpotSpring.h"
#include "mesa_pd/mpi/ReduceContactHistory.h"
@@ -98,7 +98,7 @@ int main(int argc, char **argv)
false);
// explicit euler
- kernel::ExplicitEuler explEuler(dt);
+ kernel::SemiImplicitEuler implEuler(dt);
collision_detection::AnalyticContactDetection acd;
kernel::DoubleCast double_cast;
kernel::SpringDashpot sd(1);
@@ -147,7 +147,7 @@ int main(int argc, char **argv)
ps->forEachParticle(false,
kernel::SelectLocal(),
*ac,
- explEuler,
+ implEuler,
*ac);
}
WALBERLA_LOG_DEVEL_VAR(sp->getPosition());
@@ -193,7 +193,7 @@ int main(int argc, char **argv)
ps->forEachParticle(false,
kernel::SelectLocal(),
*ac,
- explEuler,
+ implEuler,
*ac);
}
WALBERLA_LOG_DEVEL_VAR(sp->getPosition());
diff --git a/tests/mesa_pd/kernel/interfaces/SemiImplicitEulerInterfaceCheck.cpp b/tests/mesa_pd/kernel/interfaces/SemiImplicitEulerInterfaceCheck.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ad498f1d909ff1cf404d41ad92808603b26669b6
--- /dev/null
+++ b/tests/mesa_pd/kernel/interfaces/SemiImplicitEulerInterfaceCheck.cpp
@@ -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 SemiImplicitEulerInterfaceCheck.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/SemiImplicitEuler.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& getPosition(const size_t /*p_idx*/) const {return position_;}
+ void setPosition(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { position_ = v;}
+
+ 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::real_t& getInvMass(const size_t /*p_idx*/) const {return invMass_;}
+
+ 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::Rot3& getRotation(const size_t /*p_idx*/) const {return rotation_;}
+ void setRotation(const size_t /*p_idx*/, const walberla::mesa_pd::Rot3& v) { rotation_ = v;}
+
+ 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::Mat3& getInvInertiaBF(const size_t /*p_idx*/) const {return invInertiaBF_;}
+
+ 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 position_;
+ walberla::mesa_pd::Vec3 linearVelocity_;
+ walberla::real_t invMass_;
+ walberla::mesa_pd::Vec3 force_;
+ walberla::mesa_pd::data::particle_flags::FlagT flags_;
+ walberla::mesa_pd::Rot3 rotation_;
+ walberla::mesa_pd::Vec3 angularVelocity_;
+ walberla::mesa_pd::Mat3 invInertiaBF_;
+ walberla::mesa_pd::Vec3 torque_;
+};
+
+
+
+} //namespace mesa_pd
+} //namespace walberla
\ No newline at end of file