diff --git a/src/core/lees_edwards/protocols.hpp b/src/core/lees_edwards/protocols.hpp index b583b13c5e..615558b1f8 100644 --- a/src/core/lees_edwards/protocols.hpp +++ b/src/core/lees_edwards/protocols.hpp @@ -53,22 +53,28 @@ struct LinearShear { /** Lees-Edwards protocol for oscillatory shearing */ struct OscillatoryShear { OscillatoryShear() - : m_initial_pos_offset{0.}, m_amplitude{0.}, m_omega{0.}, m_time_0{0.} {} + : m_initial_pos_offset{0.}, m_amplitude{0.}, m_omega{0.}, m_time_0{0.}, + m_decay_rate{0.} {} OscillatoryShear(double initial_offset, double amplitude, double omega, double time_0) : m_initial_pos_offset{initial_offset}, m_amplitude{amplitude}, m_omega{omega}, m_time_0{time_0} {} double pos_offset(double time) const { - return m_initial_pos_offset + - m_amplitude * std::sin(m_omega * (time - m_time_0)); + return m_initial_pos_offset + m_amplitude * + exp(-(time - m_time_0) * m_decay_rate) * + std::sin(m_omega * (time - m_time_0)); } double shear_velocity(double time) const { - return m_omega * m_amplitude * std::cos(m_omega * (time - m_time_0)); + return -m_decay_rate * exp(-(time - m_time_0) * m_decay_rate) * + m_amplitude * std::sin(m_omega * (time - m_time_0)) + + exp(-(time - m_time_0) * m_decay_rate) * m_omega * m_amplitude * + std::cos(m_omega * (time - m_time_0)); } double m_initial_pos_offset; double m_amplitude; double m_omega; double m_time_0; + double m_decay_rate; }; /** Type which holds the currently active protocol */ diff --git a/src/python/espressomd/lees_edwards.py b/src/python/espressomd/lees_edwards.py index 826fc1d729..71230f9a3a 100644 --- a/src/python/espressomd/lees_edwards.py +++ b/src/python/espressomd/lees_edwards.py @@ -96,6 +96,8 @@ class OscillatoryShear(ScriptInterfaceHelper): Radian frequency of the oscillation. time_0 : :obj:`float` Time offset of the oscillation. + decay_rate : :obj:`float` + Apply an exponential decay with the given rate to the oscillation's amplitude. Defaults to 0, i.e., no decay. """ _so_name = "LeesEdwards::OscillatoryShear" diff --git a/src/script_interface/lees_edwards/OscillatoryShear.hpp b/src/script_interface/lees_edwards/OscillatoryShear.hpp index d860c6c77c..81520f9e05 100644 --- a/src/script_interface/lees_edwards/OscillatoryShear.hpp +++ b/src/script_interface/lees_edwards/OscillatoryShear.hpp @@ -37,11 +37,13 @@ class OscillatoryShear : public Protocol { OscillatoryShear() : m_protocol{ std::make_shared<::LeesEdwards::ActiveProtocol>(CoreClass())} { - add_parameters({{"initial_pos_offset", - std::get(*m_protocol).m_initial_pos_offset}, - {"amplitude", std::get(*m_protocol).m_amplitude}, - {"omega", std::get(*m_protocol).m_omega}, - {"time_0", std::get(*m_protocol).m_time_0}}); + add_parameters( + {{"initial_pos_offset", + std::get(*m_protocol).m_initial_pos_offset}, + {"amplitude", std::get(*m_protocol).m_amplitude}, + {"omega", std::get(*m_protocol).m_omega}, + {"time_0", std::get(*m_protocol).m_time_0}, + {"decay_rate", std::get(*m_protocol).m_decay_rate}}); } std::shared_ptr<::LeesEdwards::ActiveProtocol> protocol() override { return m_protocol; diff --git a/testsuite/python/lees_edwards.py b/testsuite/python/lees_edwards.py index 4ebdb2dafb..e0af2956a5 100644 --- a/testsuite/python/lees_edwards.py +++ b/testsuite/python/lees_edwards.py @@ -28,10 +28,18 @@ import itertools +def deriv(f, x, h=1E-5): + # central difference quotient + return 1 / (2 * h) * (f(x + h) - f(x - h)) + + np.random.seed(42) params_lin = {'initial_pos_offset': 0.1, 'time_0': 0.1, 'shear_velocity': 1.2} params_osc = {'initial_pos_offset': 0.1, 'time_0': -2.1, 'amplitude': 2.3, - 'omega': 2.51} + 'omega': 2.51, "decay_rate": 0} +params_osc_decay = {'initial_pos_offset': 0.1, 'time_0': -2.1, 'amplitude': 2.3, + 'omega': 2.51, "decay_rate": 0.1} + lin_protocol = espressomd.lees_edwards.LinearShear(**params_lin) @@ -41,6 +49,8 @@ def get_lin_pos_offset(time, initial_pos_offset=None, osc_protocol = espressomd.lees_edwards.OscillatoryShear(**params_osc) +osc_decay_protocol = espressomd.lees_edwards.OscillatoryShear( + **params_osc_decay) off_protocol = espressomd.lees_edwards.Off() const_offset_protocol = espressomd.lees_edwards.LinearShear( initial_pos_offset=2.2, shear_velocity=0) @@ -117,10 +127,6 @@ def test_protocols(self): system.lees_edwards.pos_offset, expected_pos) # Check if the offset is determined correctly - - system.time = 0.0 - - # Oscillatory shear system.lees_edwards.protocol = osc_protocol # check parameter setter/getter consistency @@ -139,6 +145,8 @@ def test_protocols(self): self.assertAlmostEqual( system.lees_edwards.pos_offset, expected_pos) + system.time = 0.0 + # Check that time change during integration updates offsets system.integrator.run(1) time = system.time @@ -160,6 +168,30 @@ def test_protocols(self): self.assertEqual(system.lees_edwards.shear_direction, "y") self.assertEqual(system.lees_edwards.shear_plane_normal, "z") + # Oscillatory shear + # Oscillatory shear with exponential decay + system.lees_edwards.protocol = osc_decay_protocol + + # check parameter setter/getter consistency + self.assertEqual( + system.lees_edwards.protocol.get_params(), params_osc_decay) + + def osc_decay_pos(t): return params_osc_decay["initial_pos_offset"] +\ + params_osc_decay["amplitude"] * np.exp(-(t - params_osc_decay["time_0"]) * params_osc_decay["decay_rate"]) *\ + np.sin(params_osc_decay["omega"] * + (t - params_osc_decay["time_0"])) + + # check pos offset and shear velocity at different times, + # check that LE offsets are recalculated on simulation time change + for time in [0., 2.3]: + system.time = time + expected_pos = osc_decay_pos(time) + expected_vel = deriv(osc_decay_pos, time) + self.assertAlmostEqual( + system.lees_edwards.pos_offset, expected_pos) + self.assertAlmostEqual( + system.lees_edwards.shear_velocity, expected_vel) + # Check that LE is disabled correctly via parameter system.lees_edwards.protocol = None self.assertIsNone(system.lees_edwards.shear_direction)