Commit 4cb13c54 authored by Sebastian Eibl's avatar Sebastian Eibl
Browse files

reworked accuracy test of mesa_pd integrators

- checking order of integrator
- checking energy conservation
- added semi implicit euler
- improved documentation
- playground for integrator analysis
parent a87980ad
Pipeline #28548 passed with stages
in 399 minutes and 36 seconds
......@@ -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 )
......
waLBerla_link_files_to_builddir( "*.ipynb" )
waLBerla_add_executable ( NAME IntegratorAccuracy
DEPENDS core mesa_pd )
//======================================================================================================================
//
// 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);
}
%% Cell type:markdown id: tags:
# error analysis of integration kernels
%% Cell type:code id: tags:
``` python
import matplotlib.pyplot as plt
import numpy as np
```
%% Cell type:code id: tags:
``` python
integrator = 'euler'
data01 = np.loadtxt(f'{integrator}01.txt').transpose()
data02 = np.loadtxt(f'{integrator}02.txt').transpose()
data04 = np.loadtxt(f'{integrator}04.txt').transpose()
```
%% Cell type:code id: tags:
``` python
plt.plot(data01[0], data01[7], label='dt = 0.1')
plt.plot(data02[0], data02[7], label='dt = 0.2')
plt.plot(data04[0], data04[7], label='dt = 0.4')
plt.xlabel('time')
plt.ylabel('max pos error')
plt.legend()
plt.show()
```
%% Cell type:code id: tags:
``` python
plt.plot(data01[0], data01[8], label='dt = 0.1')
plt.plot(data02[0], data02[8], label='dt = 0.2')
plt.plot(data04[0], data04[8], label='dt = 0.4')
plt.xlabel('time')
plt.ylabel('max vel error')
plt.legend()
plt.show()
```
%% Cell type:code id: tags:
``` python
final_errors = np.ones(3)
final_errors[0] = np.abs(data01[7][-1])
final_errors[1] = np.abs(data02[7][-1])
final_errors[2] = np.abs(data04[7][-1])
plt.bar(np.arange(3), final_errors)
plt.show()
print(final_errors[2] / final_errors[1])
print(final_errors[1] / final_errors[0])
```
%% Cell type:code id: tags:
``` python
plt.plot(data01[0], data01[1], label='analytic')
plt.plot(data01[0], data01[4], label='dt = 0.1')
plt.plot(data02[0], data02[4], label='dt = 0.2')
plt.plot(data04[0], data04[4], label='dt = 0.4')
plt.xlabel('time')
plt.ylabel('pos')
plt.legend()
plt.show()
```
%% Cell type:code id: tags:
``` python
plt.plot(data01[0], data01[2], label='analytic')
plt.plot(data01[0], data01[5], label='dt = 0.1')
plt.plot(data02[0], data02[5], label='dt = 0.2')
plt.plot(data04[0], data04[5], label='dt = 0.4')
plt.xlabel('time')
plt.ylabel('vel')
plt.legend()
plt.show()
```
%% Cell type:code id: tags:
``` python
plt.plot(data01[0], data01[6] / data01[3][0], label='dt = 0.1')
plt.plot(data02[0], data02[6] / data01[3][0], label='dt = 0.2')
plt.plot(data04[0], data04[6] / data01[3][0], label='dt = 0.4')
plt.xlabel('time')
plt.ylabel('relative energy')
plt.legend()
plt.show()
```
%% Cell type:code id: tags:
``` python
```
......@@ -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())
......
# -*- 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')
......@@ -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',
......
......@@ -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);
......
//======================================================================================================================
//
// 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));