Issue #354 analytical damper fix

elastica/dissipation.py

@@ -4,10 +4,11 @@
 Built in damper module implementations
 """
 
+import logging
 from abc import ABC, abstractmethod
-from typing import Any, Generic, TypeVar
+from typing import Any, Generic, TypeVar, TypeAlias, Callable
 
-from elastica.typing import RodType, SystemType
+from elastica.typing import RodType
 
 from numba import njit
 
@@ -74,6 +75,9 @@ def dampen_rates(self, system: T, time: np.float64) -> None:
         pass
 
 
+DampenType: TypeAlias = Callable[[RodType], None]
+
+
 class AnalyticalLinearDamper(DamperBase):
     """
     Analytical linear damper class. This class corresponds to the analytical version of
@@ -82,17 +86,44 @@ class AnalyticalLinearDamper(DamperBase):
 
     .. math::
 
-        \\mathbf{v}^{n+1} = \\mathbf{v}^n \\exp \\left( -  \\nu~dt  \\right)
+        m \\frac{\\partial \\mathbf{v}}{\\partial t} = -\\gamma_t \\mathbf{v}
 
-        \\pmb{\\omega}^{n+1} = \\pmb{\\omega}^n \\exp \\left( - \\frac{{\\nu}~m~dt } { \\mathbf{J}} \\right)
+        \\frac{\\mathbf{J}}{e} \\frac{\\partial \\pmb{\\omega}}{\\partial t} = -\\gamma_r \\pmb{\\omega}
 
     Examples
     --------
-    How to set analytical linear damper for rod or rigid body:
+    The current AnalyticalLinearDamper class supports three types of protocols:
+
+    1.  Uniform damping constant: the user provides the keyword argument `uniform_damping_constant`
+        of dimension (1/T). This leads to an exponential damping constant that is used for both
+        translation and rotational velocities.
+
+    >>> simulator.dampen(rod).using(
+    ...     AnalyticalLinearDamper,
+    ...     uniform_damping_constant=0.1,
+    ...     time_step = 1E-4,   # Simulation time-step
+    ... )
+
+    2.  Physical damping constant: separate exponential coefficients are computed for the
+        translational and rotational velocities, based on user-defined
+        `translational_damping_constant` and `rotational_damping_constant`.
+
+    >>> simulator.dampen(rod).using(
+    ...     AnalyticalLinearDamper,
+    ...     translational_damping_constant=0.1,
+    ...     rotational_damping_constant=0.05,
+    ...     time_step = 1E-4,   # Simulation time-step
+    ... )
+
+    3.  Damping constant: this protocol follows the original algorithm where the damping
+        constants for translational and rotational velocities are assumed to be numerically
+        identical. This leads to dimensional inconsistencies (see
+        https://github.com/GazzolaLab/PyElastica/issues/354). Hence, this option will be deprecated
+        in version 0.4.0.
 
     >>> simulator.dampen(rod).using(
     ...     AnalyticalLinearDamper,
-    ...     damping_constant=0.1,
+    ...     damping_constant=0.1,   # To be deprecated in 0.4.0
     ...     time_step = 1E-4,   # Simulation time-step
     ... )
 
@@ -108,58 +139,147 @@ class AnalyticalLinearDamper(DamperBase):
     1. Set a high value for `damping_constant` to first acheive a stable simulation.
     2. If you feel the simulation is overdamped, reduce `damping_constant` until you
        feel the simulation is underdamped, and expected dynamics are recovered.
-
-    Attributes
-    ----------
-    translational_damping_coefficient: numpy.ndarray
-        1D array containing data with 'float' type.
-        Damping coefficient acting on translational velocity.
-    rotational_damping_coefficient : numpy.ndarray
-        1D array containing data with 'float' type.
-        Damping coefficient acting on rotational velocity.
     """
 
-    def __init__(
-        self, damping_constant: float, time_step: float, **kwargs: Any
-    ) -> None:
-        """
-        Analytical linear damper initializer
-
-        Parameters
-        ----------
-        damping_constant : float
-            Damping constant for the analytical linear damper.
-        time_step : float
-            Time-step of simulation
-        """
+    def __init__(self, time_step: np.float64, **kwargs: Any) -> None:
         super().__init__(**kwargs)
-        # Compute the damping coefficient for translational velocity
+
+        damping_constant = kwargs.get("damping_constant", None)
+        uniform_damping_constant = kwargs.get("uniform_damping_constant", None)
+        translational_damping_constant = kwargs.get(
+            "translational_damping_constant", None
+        )
+        rotational_damping_constant = kwargs.get("rotational_damping_constant", None)
+
+        self._dampen_rates_protocol: DampenType
+
+        if (
+            (damping_constant is not None)
+            and (uniform_damping_constant is None)
+            and (translational_damping_constant is None)
+            and (rotational_damping_constant is None)
+        ):
+            logging.warning(
+                "Analytical linear damping using generic damping constant "
+                "will be deprecated in 0.4.0"
+            )
+            self._dampen_rates_protocol = self._deprecated_damping_protocol(
+                damping_constant=damping_constant, time_step=time_step
+            )
+
+        elif (
+            (damping_constant is None)
+            and (uniform_damping_constant is not None)
+            and (translational_damping_constant is None)
+            and (rotational_damping_constant is None)
+        ):
+            self._dampen_rates_protocol = self._uniform_damping_protocol(
+                uniform_damping_constant=uniform_damping_constant, time_step=time_step
+            )
+
+        elif (
+            (damping_constant is None)
+            and (uniform_damping_constant is None)
+            and (translational_damping_constant is not None)
+            and (rotational_damping_constant is not None)
+        ):
+            self._dampen_rates_protocol = self._physical_damping_protocol(
+                translational_damping_constant=translational_damping_constant,
+                rotational_damping_constant=rotational_damping_constant,
+                time_step=time_step,
+            )
+
+        else:
+            message = (
+                "AnalyticalLinearDamper usage:\n"
+                "\tsimulator.dampen(rod).using(\n"
+                "\t\tAnalyticalLinearDamper,\n"
+                "\t\ttranslational_damping_constant=...,\n"
+                "\t\trotational_damping_constant=...,\n"
+                "\t\ttime_step=...,\n"
+                "\t)\n"
+                "\tor\n"
+                "\tsimulator.dampen(rod).using(\n"
+                "\t\tAnalyticalLinearDamper,\n"
+                "\t\tuniform_damping_constant=...,\n"
+                "\t\ttime_step=...,\n"
+                "\t)\n"
+                "\tor (deprecated in 0.4.0)\n"
+                "\tsimulator.dampen(rod).using(\n"
+                "\t\tAnalyticalLinearDamper,\n"
+                "\t\tdamping_constant=...,\n"
+                "\t\ttime_step=...,\n"
+                "\t)\n"
+            )
+            raise ValueError(message)
+
+    def _deprecated_damping_protocol(
+        self, damping_constant: np.float64, time_step: np.float64
+    ) -> DampenType:
         nodal_mass = self._system.mass
-        self.translational_damping_coefficient = np.exp(-damping_constant * time_step)
+        self._translational_damping_coefficient = np.exp(-damping_constant * time_step)
 
-        # Compute the damping coefficient for exponential velocity
         if self._system.ring_rod_flag:
             element_mass = nodal_mass
         else:
             element_mass = 0.5 * (nodal_mass[1:] + nodal_mass[:-1])
             element_mass[0] += 0.5 * nodal_mass[0]
             element_mass[-1] += 0.5 * nodal_mass[-1]
-        self.rotational_damping_coefficient = np.exp(
+        self._rotational_damping_coefficient = np.exp(
             -damping_constant
             * time_step
             * element_mass
             * np.diagonal(self._system.inv_mass_second_moment_of_inertia).T
         )
 
-    def dampen_rates(self, system: RodType, time: np.float64) -> None:
-        system.velocity_collection[:] = (
-            system.velocity_collection * self.translational_damping_coefficient
+        def dampen_rates_protocol(rod: RodType) -> None:
+            rod.velocity_collection *= self._translational_damping_coefficient
+            rod.omega_collection *= np.power(
+                self._rotational_damping_coefficient, rod.dilatation
+            )
+
+        return dampen_rates_protocol
+
+    def _uniform_damping_protocol(
+        self, uniform_damping_constant: np.float64, time_step: np.float64
+    ) -> DampenType:
+        self._translational_damping_coefficient = (
+            self._rotational_damping_coefficient
+        ) = np.exp(-uniform_damping_constant * time_step)
+
+        def dampen_rates_protocol(rod: RodType) -> None:
+            rod.velocity_collection *= self._translational_damping_coefficient
+            rod.omega_collection *= self._rotational_damping_coefficient
+
+        return dampen_rates_protocol
+
+    def _physical_damping_protocol(
+        self,
+        translational_damping_constant: np.float64,
+        rotational_damping_constant: np.float64,
+        time_step: np.float64,
+    ) -> DampenType:
+        nodal_mass = self._system.mass
+        self._translational_damping_coefficient = np.exp(
+            -translational_damping_constant / nodal_mass * time_step
         )
 
-        system.omega_collection[:] = system.omega_collection * np.power(
-            self.rotational_damping_coefficient, system.dilatation
+        inv_moi = np.diagonal(self._system.inv_mass_second_moment_of_inertia).T
+        self._rotational_damping_coefficient = np.exp(
+            -rotational_damping_constant * inv_moi * time_step
         )
 
+        def dampen_rates_protocol(rod: RodType) -> None:
+            rod.velocity_collection *= self._translational_damping_coefficient
+            rod.omega_collection *= np.power(
+                self._rotational_damping_coefficient, rod.dilatation
+            )
+
+        return dampen_rates_protocol
+
+    def dampen_rates(self, rod: RodType, time: np.float64) -> None:
+        self._dampen_rates_protocol(rod)
+
 
 class LaplaceDissipationFilter(DamperBase):
     """