Lees Edwads: Add optional exp. decay to osc shear protocol

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 */

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"

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<CoreClass>(*m_protocol).m_initial_pos_offset},
-                    {"amplitude", std::get<CoreClass>(*m_protocol).m_amplitude},
-                    {"omega", std::get<CoreClass>(*m_protocol).m_omega},
-                    {"time_0", std::get<CoreClass>(*m_protocol).m_time_0}});
+    add_parameters(
+        {{"initial_pos_offset",
+          std::get<CoreClass>(*m_protocol).m_initial_pos_offset},
+         {"amplitude", std::get<CoreClass>(*m_protocol).m_amplitude},
+         {"omega", std::get<CoreClass>(*m_protocol).m_omega},
+         {"time_0", std::get<CoreClass>(*m_protocol).m_time_0},
+         {"decay_rate", std::get<CoreClass>(*m_protocol).m_decay_rate}});
   }
   std::shared_ptr<::LeesEdwards::ActiveProtocol> protocol() override {
     return m_protocol;

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)