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

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

src/script_interface/lees_edwards/OscillatoryShear.hpp

@@ -41,7 +41,8 @@ class OscillatoryShear : public Protocol {
                      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}});
+                    {"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

@@ -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