From c83095dacdd1057ecde44756272ab76882aad674 Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Fri, 8 Mar 2024 21:19:00 +0100 Subject: [PATCH 1/3] Lees Edwads: Add optional exp. decay to osc shear protocol --- src/core/lees_edwards/protocols.hpp | 12 +++--- .../lees_edwards/OscillatoryShear.hpp | 3 +- testsuite/python/lees_edwards.py | 41 ++++++++++++++++--- 3 files changed, 45 insertions(+), 11 deletions(-) diff --git a/src/core/lees_edwards/protocols.hpp b/src/core/lees_edwards/protocols.hpp index 8c4fc0f936..bb45fd8084 100644 --- a/src/core/lees_edwards/protocols.hpp +++ b/src/core/lees_edwards/protocols.hpp @@ -53,22 +53,24 @@ 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) + double time_0, double decay_rate) : m_initial_pos_offset{initial_offset}, - m_amplitude{amplitude}, m_omega{omega}, m_time_0{time_0} {} + m_amplitude{amplitude}, m_omega{omega}, m_time_0{time_0}, m_decay_rate{decay_rate} {} double pos_offset(double time) const { return m_initial_pos_offset + - m_amplitude * std::sin(m_omega * (time - m_time_0)); + 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/script_interface/lees_edwards/OscillatoryShear.hpp b/src/script_interface/lees_edwards/OscillatoryShear.hpp index d860c6c77c..807e590506 100644 --- a/src/script_interface/lees_edwards/OscillatoryShear.hpp +++ b/src/script_interface/lees_edwards/OscillatoryShear.hpp @@ -41,7 +41,8 @@ class OscillatoryShear : public Protocol { 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}}); + {"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 ee90f51439..60e72e87d5 100644 --- a/testsuite/python/lees_edwards.py +++ b/testsuite/python/lees_edwards.py @@ -27,11 +27,17 @@ import numpy as np 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 +47,7 @@ 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() @@ -109,10 +116,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 @@ -131,6 +134,9 @@ 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 @@ -151,6 +157,31 @@ def test_protocols(self): shear_direction="y", shear_plane_normal="z", protocol=lin_protocol) 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) + osc_decay_pos = lambda t: 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 From fa441bad4cf23c578071ac013094c47988766624 Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Wed, 14 Aug 2024 16:28:53 +0200 Subject: [PATCH 2/3] Formatting --- .../lees_edwards/OscillatoryShear.hpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/script_interface/lees_edwards/OscillatoryShear.hpp b/src/script_interface/lees_edwards/OscillatoryShear.hpp index 807e590506..81520f9e05 100644 --- a/src/script_interface/lees_edwards/OscillatoryShear.hpp +++ b/src/script_interface/lees_edwards/OscillatoryShear.hpp @@ -37,12 +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}, - {"decay_rate", std::get(*m_protocol).m_decay_rate}}); + 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; From 97ac1731ea45684c01751d0cefb24c66db8be3e2 Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Wed, 14 Aug 2024 16:31:46 +0200 Subject: [PATCH 3/3] doc --- src/python/espressomd/lees_edwards.py | 2 ++ 1 file changed, 2 insertions(+) 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"