From 9b67a1ea082c588422a3fb269825178e761f8ff8 Mon Sep 17 00:00:00 2001 From: Songyuan Cui Date: Fri, 27 Sep 2024 15:46:24 -0500 Subject: [PATCH 01/12] fix: replace analytical linear damper with dimensionally consistent versions, and move dampers into dedicated directory --- elastica/dissipation/__init__.py | 9 + .../dissipation/analytical_linear_damper.py | 194 ++++++++++++++++++ elastica/dissipation/damper_base.py | 62 ++++++ .../laplacian_dissipation_filter.py} | 161 +-------------- 4 files changed, 269 insertions(+), 157 deletions(-) create mode 100644 elastica/dissipation/__init__.py create mode 100644 elastica/dissipation/analytical_linear_damper.py create mode 100644 elastica/dissipation/damper_base.py rename elastica/{dissipation.py => dissipation/laplacian_dissipation_filter.py} (61%) diff --git a/elastica/dissipation/__init__.py b/elastica/dissipation/__init__.py new file mode 100644 index 00000000..a9ddd44d --- /dev/null +++ b/elastica/dissipation/__init__.py @@ -0,0 +1,9 @@ +__doc__ = """ +(added in version 0.3.0) + +Built in damper module implementations +""" + +from .damper_base import DamperBase +from .analytical_linear_damper import AnalyticalLinearDamper +from .laplacian_dissipation_filter import LaplaceDissipationFilter diff --git a/elastica/dissipation/analytical_linear_damper.py b/elastica/dissipation/analytical_linear_damper.py new file mode 100644 index 00000000..584096e9 --- /dev/null +++ b/elastica/dissipation/analytical_linear_damper.py @@ -0,0 +1,194 @@ +import logging +from typing import Any, Callable, TypeAlias + +import numpy as np + +from elastica.typing import RodType +from elastica.dissipation.damper_base import DamperBase + + +DampenType: TypeAlias = Callable[[RodType], None] + + +class AnalyticalLinearDamper(DamperBase): + """ + Analytical linear damper class. This class corresponds to the analytical version of + a linear damper, and uses the following equations to damp translational and + rotational velocities: + + .. math:: + + \\mathbf{v}^{n+1} = \\mathbf{v}^n \\exp \\left( - \\nu~dt \\right) + + \\pmb{\\omega}^{n+1} = \\pmb{\\omega}^n \\exp \\left( - \\frac{{\\nu}~m~dt } { \\mathbf{J}} \\right) + + Examples + -------- + How to set analytical linear damper for rod or rigid body: + + >>> simulator.dampen(rod).using( + ... AnalyticalLinearDamper, + ... damping_constant=0.1, + ... time_step = 1E-4, # Simulation time-step + ... ) + + Notes + ----- + Since this class analytically treats the damping term, it is unconditionally stable + from a timestep perspective, i.e. the presence of damping does not impose any additional + restriction on the simulation timestep size. This implies that when using + AnalyticalLinearDamper, one can set `damping_constant` as high as possible, without worrying + about the simulation becoming unstable. This now leads to a streamlined procedure + for tuning the `damping_constant`: + + 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, time_step: np.float64, **kwargs: Any) -> None: + super().__init__(**kwargs) + + 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: + raise ValueError( + "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" + ) + + 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) + + 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( + -damping_constant + * time_step + * element_mass + * np.diagonal(self._system.inv_mass_second_moment_of_inertia).T + ) + + 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.view() + self._translational_damping_coefficient = np.exp( + -translational_damping_constant / nodal_mass * time_step + ) + + 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) diff --git a/elastica/dissipation/damper_base.py b/elastica/dissipation/damper_base.py new file mode 100644 index 00000000..264a1d59 --- /dev/null +++ b/elastica/dissipation/damper_base.py @@ -0,0 +1,62 @@ +from abc import ABC, abstractmethod +from typing import Any, Generic, TypeVar +import numpy as np + + +T = TypeVar("T") + + +class DamperBase(Generic[T], ABC): + """Base class for damping module implementations. + + Notes + ----- + All damper classes must inherit DamperBase class. + + + Attributes + ---------- + system : RodBase + + """ + + _system: T + + # TODO typing can be made better + def __init__(self, *args: Any, **kwargs: Any) -> None: + """Initialize damping module""" + try: + self._system = kwargs["_system"] + except KeyError: + raise KeyError( + "Please use simulator.dampen(...).using(...) syntax to establish " + "damping." + ) + + @property + def system(self) -> T: + """ + get system (rod or rigid body) reference + + Returns + ------- + SystemType + + """ + return self._system + + @abstractmethod + def dampen_rates(self, system: T, time: np.float64) -> None: + # TODO: In the future, we can remove rod and use self.system + """ + Dampen rates (velocity and/or omega) of a rod object. + + Parameters + ---------- + system : SystemType + System (rod or rigid-body) object. + time : float + The time of simulation. + + """ + pass diff --git a/elastica/dissipation.py b/elastica/dissipation/laplacian_dissipation_filter.py similarity index 61% rename from elastica/dissipation.py rename to elastica/dissipation/laplacian_dissipation_filter.py index 32ff4866..66c6e913 100644 --- a/elastica/dissipation.py +++ b/elastica/dissipation/laplacian_dissipation_filter.py @@ -1,164 +1,11 @@ -__doc__ = """ -(added in version 0.3.0) - -Built in damper module implementations -""" - -from abc import ABC, abstractmethod -from typing import Any, Generic, TypeVar - -from elastica.typing import RodType, SystemType - -from numba import njit +from typing import Any import numpy as np from numpy.typing import NDArray +from numba import njit - -T = TypeVar("T") - - -class DamperBase(Generic[T], ABC): - """Base class for damping module implementations. - - Notes - ----- - All damper classes must inherit DamperBase class. - - - Attributes - ---------- - system : RodBase - - """ - - _system: T - - # TODO typing can be made better - def __init__(self, *args: Any, **kwargs: Any) -> None: - """Initialize damping module""" - try: - self._system = kwargs["_system"] - except KeyError: - raise KeyError( - "Please use simulator.dampen(...).using(...) syntax to establish " - "damping." - ) - - @property - def system(self) -> T: - """ - get system (rod or rigid body) reference - - Returns - ------- - SystemType - - """ - return self._system - - @abstractmethod - def dampen_rates(self, system: T, time: np.float64) -> None: - # TODO: In the future, we can remove rod and use self.system - """ - Dampen rates (velocity and/or omega) of a rod object. - - Parameters - ---------- - system : SystemType - System (rod or rigid-body) object. - time : float - The time of simulation. - - """ - pass - - -class AnalyticalLinearDamper(DamperBase): - """ - Analytical linear damper class. This class corresponds to the analytical version of - a linear damper, and uses the following equations to damp translational and - rotational velocities: - - .. math:: - - \\mathbf{v}^{n+1} = \\mathbf{v}^n \\exp \\left( - \\nu~dt \\right) - - \\pmb{\\omega}^{n+1} = \\pmb{\\omega}^n \\exp \\left( - \\frac{{\\nu}~m~dt } { \\mathbf{J}} \\right) - - Examples - -------- - How to set analytical linear damper for rod or rigid body: - - >>> simulator.dampen(rod).using( - ... AnalyticalLinearDamper, - ... damping_constant=0.1, - ... time_step = 1E-4, # Simulation time-step - ... ) - - Notes - ----- - Since this class analytically treats the damping term, it is unconditionally stable - from a timestep perspective, i.e. the presence of damping does not impose any additional - restriction on the simulation timestep size. This implies that when using - AnalyticalLinearDamper, one can set `damping_constant` as high as possible, without worrying - about the simulation becoming unstable. This now leads to a streamlined procedure - for tuning the `damping_constant`: - - 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 - """ - super().__init__(**kwargs) - # Compute the damping coefficient for translational velocity - nodal_mass = self._system.mass - 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( - -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 - ) - - system.omega_collection[:] = system.omega_collection * np.power( - self.rotational_damping_coefficient, system.dilatation - ) +from elastica.typing import RodType +from elastica.dissipation.damper_base import DamperBase class LaplaceDissipationFilter(DamperBase): From 9828699ee71dc8fd0d956eca66bf8371849b4714 Mon Sep 17 00:00:00 2001 From: Songyuan Cui Date: Thu, 10 Oct 2024 14:47:41 -0500 Subject: [PATCH 02/12] refactor: refine error message formatting --- elastica/dissipation/analytical_linear_damper.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/elastica/dissipation/analytical_linear_damper.py b/elastica/dissipation/analytical_linear_damper.py index 584096e9..660b987f 100644 --- a/elastica/dissipation/analytical_linear_damper.py +++ b/elastica/dissipation/analytical_linear_damper.py @@ -104,7 +104,7 @@ def __init__(self, time_step: np.float64, **kwargs: Any) -> None: ) else: - raise ValueError( + message = ( "AnalyticalLinearDamper usage:\n" "\tsimulator.dampen(rod).using(\n" "\t\tAnalyticalLinearDamper,\n" @@ -125,6 +125,7 @@ def __init__(self, time_step: np.float64, **kwargs: Any) -> None: "\t\ttime_step=...,\n" "\t)\n" ) + raise ValueError(message) def _deprecated_damping_protocol( self, damping_constant: np.float64, time_step: np.float64 From e420603670438bcbf8eb7eaa05bd888a2da26e01 Mon Sep 17 00:00:00 2001 From: Songyuan Cui Date: Thu, 10 Oct 2024 14:48:42 -0500 Subject: [PATCH 03/12] refactor: move test_dissipation.py to dedicated directory --- .../test_analytical_linear_damper.py | 193 ++++++++++++++++++ tests/test_dissipation/test_damper_base.py | 43 ++++ .../test_laplacian_dissipation_filter.py} | 120 +---------- 3 files changed, 244 insertions(+), 112 deletions(-) create mode 100644 tests/test_dissipation/test_analytical_linear_damper.py create mode 100644 tests/test_dissipation/test_damper_base.py rename tests/{test_dissipation.py => test_dissipation/test_laplacian_dissipation_filter.py} (61%) diff --git a/tests/test_dissipation/test_analytical_linear_damper.py b/tests/test_dissipation/test_analytical_linear_damper.py new file mode 100644 index 00000000..a87a0a9d --- /dev/null +++ b/tests/test_dissipation/test_analytical_linear_damper.py @@ -0,0 +1,193 @@ +__doc__ = """ Test implementation of the analytical linear damper""" + +from itertools import combinations + +import pytest +import numpy as np +from numpy.testing import assert_allclose + +from elastica.dissipation import AnalyticalLinearDamper +from elastica.utils import Tolerance +from tests.test_rod.mock_rod import MockTestRod + + +@pytest.fixture +def analytical_error_message(): + message = ( + r"AnalyticalLinearDamper usage:\n" + r"\tsimulator.dampen\(rod\).using\(\n" + r"\t\tAnalyticalLinearDamper,\n" + r"\t\ttranslational_damping_constant=...,\n" + r"\t\trotational_damping_constant=...,\n" + r"\t\ttime_step=...,\n" + r"\t\)\n" + r"\tor\n" + r"\tsimulator.dampen\(rod\).using\(\n" + r"\t\tAnalyticalLinearDamper,\n" + r"\t\tuniform_damping_constant=...,\n" + r"\t\ttime_step=...,\n" + r"\t\)\n" + r"\tor \(deprecated in 0.4.0\)\n" + r"\tsimulator.dampen\(rod\).using\(\n" + r"\t\tAnalyticalLinearDamper,\n" + r"\t\tdamping_constant=...,\n" + r"\t\ttime_step=...,\n" + r"\t\)\n" + ) + return message + + +def test_analytical_linear_damper_error(analytical_error_message): + test_rod = MockTestRod() + dummy = np.float64(0.0) + + kwargs = [ + "damping_constant", + "uniform_damping_constant", + "translational_damping_constant", + "rotational_damping_constant", + ] + + a = (0, 1, 2, 3) + valid_combs = [(0,), (1,), (2, 3)] + for i in range(5): + combs = list(combinations(a, i)) + for c in combs: + if c not in valid_combs: + invalid_kwargs = dict([(kwargs[idx], dummy) for idx in c]) + with pytest.raises(ValueError, match=analytical_error_message): + AnalyticalLinearDamper( + _system=test_rod, + time_step=dummy, + **invalid_kwargs, + ) + + +def test_analytical_linear_damper_deprecated(): + test_rod = MockTestRod() + test_rod.mass[:] = 1.0 + test_dilatation = 2.0 * np.ones((3, test_rod.n_elems)) + test_inv_mass_second_moment_of_inertia = 3.0 * np.ones((3, 3, test_rod.n_elems)) + test_rod.dilatation = test_dilatation.copy() + test_rod.inv_mass_second_moment_of_inertia = ( + test_inv_mass_second_moment_of_inertia.copy() + ) + damping_constant = 0.25 + dt = 0.5 + exponential_damper = AnalyticalLinearDamper( + _system=test_rod, damping_constant=damping_constant, time_step=dt + ) + # check common prefactors + # e ^ (-damp_coeff * dt) + ref_translational_damping_coefficient = np.exp(-0.25 * 0.5) + # e ^ (-damp_coeff * dt * elemental_mass * inv_mass_second_moment_of_inertia) + ref_rotational_damping_coefficient = np.exp(-0.25 * 0.5 * 1.0 * 3.0) * np.ones( + (3, test_rod.n_elems) + ) + # end corrections + ref_rotational_damping_coefficient[:, 0] = np.exp(-0.25 * 0.5 * 1.5 * 3.0) + ref_rotational_damping_coefficient[:, -1] = np.exp(-0.25 * 0.5 * 1.5 * 3.0) + + pre_damping_velocity_collection = np.random.rand(3, test_rod.n_elems + 1) + test_rod.velocity_collection = ( + pre_damping_velocity_collection.copy() + ) # We need copy of the list not a reference to this array + pre_damping_omega_collection = np.random.rand(3, test_rod.n_elems) + test_rod.omega_collection = ( + pre_damping_omega_collection.copy() + ) # We need copy of the list not a reference to this array + exponential_damper.dampen_rates(test_rod, time=0) + post_damping_velocity_collection = ( + pre_damping_velocity_collection * ref_translational_damping_coefficient + ) + # multiplying_factor = ref_rot_coeff ^ dilation + post_damping_omega_collection = ( + pre_damping_omega_collection * ref_rotational_damping_coefficient**2.0 + ) + assert_allclose( + post_damping_velocity_collection, + test_rod.velocity_collection, + atol=Tolerance.atol(), + ) + assert_allclose( + post_damping_omega_collection, test_rod.omega_collection, atol=Tolerance.atol() + ) + + +def test_analytical_linear_damper_uniform(): + test_rod = MockTestRod() + + velocity = np.linspace(0.0, 3.0, 3 * (test_rod.n_elems + 1)).reshape( + (3, test_rod.n_elems + 1) + ) + omega = np.linspace(5.0, 8.0, 3 * test_rod.n_elems).reshape((3, test_rod.n_elems)) + test_rod.velocity_collection[:, :] = velocity + test_rod.omega_collection[:, :] = omega + + test_constant = 2.0 + test_dt = np.float64(1.5) + test_coeff = np.exp(-test_dt * test_constant) + + damper = AnalyticalLinearDamper( + _system=test_rod, uniform_damping_constant=test_constant, time_step=test_dt + ) + damper.dampen_rates(test_rod, time=np.float64(0.0)) + + expected_velocity = velocity * test_coeff + expected_omega = omega * test_coeff + + assert_allclose( + expected_velocity, + test_rod.velocity_collection, + atol=Tolerance.atol(), + ) + assert_allclose(expected_omega, test_rod.omega_collection, atol=Tolerance.atol()) + + +def test_analytical_linear_damper_physical(): + test_rod = MockTestRod() + + velocity = np.linspace(0.0, 3.0, 3 * (test_rod.n_elems + 1)).reshape( + (3, test_rod.n_elems + 1) + ) + omega = np.linspace(5.0, 8.0, 3 * test_rod.n_elems).reshape((3, test_rod.n_elems)) + mass = np.linspace(6.0, 4.0, test_rod.n_elems + 1) + inv_moi_full = 1.0 / np.linspace(10.0, 15.0, 9 * test_rod.n_elems).reshape( + (3, 3, test_rod.n_elems) + ) + inv_moi = np.diagonal(inv_moi_full).T + dilatation = np.linspace(0.5, 1.5, test_rod.n_elems) + + test_rod.velocity_collection[:, :] = velocity + test_rod.omega_collection[:, :] = omega + test_rod.mass[:] = mass + test_rod.inv_mass_second_moment_of_inertia = inv_moi_full.copy() + test_rod.dilatation = dilatation.copy() + + test_translational_constant = 2.0 + test_rotational_constant = 3.0 + test_dt = np.float64(1.5) + + damper = AnalyticalLinearDamper( + _system=test_rod, + translational_damping_constant=test_translational_constant, + rotational_damping_constant=test_rotational_constant, + time_step=test_dt, + ) + + test_translational_coeff = np.exp(-test_dt * test_translational_constant / mass) + test_rotational_coeff = np.exp( + -test_dt * test_rotational_constant * inv_moi * dilatation.reshape((1, -1)) + ) + + damper.dampen_rates(test_rod, time=np.float64(0.0)) + + expected_velocity = velocity * test_translational_coeff + expected_omega = omega * test_rotational_coeff + + assert_allclose( + expected_velocity, + test_rod.velocity_collection, + atol=Tolerance.atol(), + ) + assert_allclose(expected_omega, test_rod.omega_collection, atol=Tolerance.atol()) diff --git a/tests/test_dissipation/test_damper_base.py b/tests/test_dissipation/test_damper_base.py new file mode 100644 index 00000000..fe2e25c6 --- /dev/null +++ b/tests/test_dissipation/test_damper_base.py @@ -0,0 +1,43 @@ +__doc__ = """ Test implementation of the damper base class""" + +import numpy as np +from numpy.testing import assert_allclose + +from elastica.dissipation import DamperBase +from elastica.utils import Tolerance +from tests.test_rod.mock_rod import MockTestRod + + +def test_damper_base(): + test_rod = MockTestRod() + test_rod.velocity_collection = np.ones(3) * 5.0 + test_rod.omega_collection = np.ones(3) * 11.0 + + class TestDamper(DamperBase): + def __init__(self, **kwargs): + super().__init__(**kwargs) + + def dampen_rates(self, rod, time): + rod.velocity_collection *= time + rod.omega_collection *= time + + test_damper = TestDamper(_system=test_rod) + test_damper.dampen_rates(test_rod, np.float64(2.0)) + assert_allclose(test_rod.velocity_collection, 10.0, atol=Tolerance.atol()) + assert_allclose(test_rod.omega_collection, 22.0, atol=Tolerance.atol()) + + +def test_damper_base_properties_access(): + test_rod = MockTestRod() + + class TestDamper(DamperBase): + def __init__(self, **kwargs): + super().__init__(**kwargs) + # Able to access properties in constraint class + assert self._system == test_rod + + def dampen_rates(self, rod, time): + assert self._system == test_rod + + test_damper = TestDamper(_system=test_rod) + test_damper.dampen_rates(test_rod, np.float64(2.0)) diff --git a/tests/test_dissipation.py b/tests/test_dissipation/test_laplacian_dissipation_filter.py similarity index 61% rename from tests/test_dissipation.py rename to tests/test_dissipation/test_laplacian_dissipation_filter.py index 8c48c680..76d92c56 100644 --- a/tests/test_dissipation.py +++ b/tests/test_dissipation/test_laplacian_dissipation_filter.py @@ -1,118 +1,14 @@ -__doc__ = """ Test Dissipation module for in Elastica implementation""" - -# System imports -from elastica.dissipation import ( - DamperBase, - AnalyticalLinearDamper, - LaplaceDissipationFilter, -) -from elastica.utils import Tolerance +__doc__ = """ Test implementation of the Laplacian dissipation filter""" +import pytest import numpy as np from numpy.testing import assert_allclose -import pytest - +from elastica.dissipation import LaplaceDissipationFilter +from elastica.utils import Tolerance from tests.test_rod.mock_rod import MockTestRod, MockTestRingRod -def test_damper_base(): - test_rod = MockTestRod() - test_rod.velocity_collection = np.ones(3) * 5.0 - test_rod.omega_collection = np.ones(3) * 11.0 - - class TestDamper(DamperBase): - def __init__(self, **kwargs): - super().__init__(**kwargs) - - def dampen_rates(self, rod, time): - rod.velocity_collection *= time - rod.omega_collection *= time - - test_damper = TestDamper(_system=test_rod) - test_damper.dampen_rates(test_rod, 2) - assert_allclose(test_rod.velocity_collection, 10.0, atol=Tolerance.atol()) - assert_allclose(test_rod.omega_collection, 22.0, atol=Tolerance.atol()) - - -def test_damper_base_properties_access(): - test_rod = MockTestRod() - - class TestDamper(DamperBase): - def __init__(self, **kwargs): - super().__init__(**kwargs) - # Able to access properties in constraint class - assert self._system == test_rod - - def dampen_rates(self, rod, time): - assert self._system == test_rod - - test_damper = TestDamper(_system=test_rod) - test_damper.dampen_rates(test_rod, 2) - - -def test_analytical_linear_damper(): - - test_rod = MockTestRod() - test_rod.mass[:] = 1.0 - test_dilatation = 2.0 * np.ones((3, test_rod.n_elems)) - test_inv_mass_second_moment_of_inertia = 3.0 * np.ones((3, 3, test_rod.n_elems)) - test_rod.dilatation = test_dilatation.copy() - test_rod.inv_mass_second_moment_of_inertia = ( - test_inv_mass_second_moment_of_inertia.copy() - ) - damping_constant = 0.25 - dt = 0.5 - exponential_damper = AnalyticalLinearDamper( - _system=test_rod, damping_constant=damping_constant, time_step=dt - ) - # check common prefactors - # e ^ (-damp_coeff * dt) - ref_translational_damping_coefficient = np.exp(-0.25 * 0.5) - # e ^ (-damp_coeff * dt * elemental_mass * inv_mass_second_moment_of_inertia) - ref_rotational_damping_coefficient = np.exp(-0.25 * 0.5 * 1.0 * 3.0) * np.ones( - (3, test_rod.n_elems) - ) - # end corrections - ref_rotational_damping_coefficient[:, 0] = np.exp(-0.25 * 0.5 * 1.5 * 3.0) - ref_rotational_damping_coefficient[:, -1] = np.exp(-0.25 * 0.5 * 1.5 * 3.0) - assert_allclose( - exponential_damper.translational_damping_coefficient, - ref_translational_damping_coefficient, - atol=Tolerance.atol(), - ) - assert_allclose( - exponential_damper.rotational_damping_coefficient, - ref_rotational_damping_coefficient, - atol=Tolerance.atol(), - ) - - pre_damping_velocity_collection = np.random.rand(3, test_rod.n_elems + 1) - test_rod.velocity_collection = ( - pre_damping_velocity_collection.copy() - ) # We need copy of the list not a reference to this array - pre_damping_omega_collection = np.random.rand(3, test_rod.n_elems) - test_rod.omega_collection = ( - pre_damping_omega_collection.copy() - ) # We need copy of the list not a reference to this array - exponential_damper.dampen_rates(test_rod, time=0) - post_damping_velocity_collection = ( - pre_damping_velocity_collection * ref_translational_damping_coefficient - ) - # multiplying_factor = ref_rot_coeff ^ dilation - post_damping_omega_collection = ( - pre_damping_omega_collection * ref_rotational_damping_coefficient**2.0 - ) - assert_allclose( - post_damping_velocity_collection, - test_rod.velocity_collection, - atol=Tolerance.atol(), - ) - assert_allclose( - post_damping_omega_collection, test_rod.omega_collection, atol=Tolerance.atol() - ) - - @pytest.mark.parametrize("filter_order", [-1, 0, 3.2]) def test_laplace_dissipation_filter_init_invalid_filter_order(filter_order): test_rod = MockTestRod() @@ -151,7 +47,7 @@ def test_laplace_dissipation_filter_for_constant_field(filter_order): ) test_rod.velocity_collection[...] = 2.0 test_rod.omega_collection[...] = 3.0 - filter_damper.dampen_rates(system=test_rod, time=0) + filter_damper.dampen_rates(system=test_rod, time=np.float64(0.0)) # filter should keep a spatially invariant field unaffected post_damping_velocity_collection = 2.0 * np.ones_like(test_rod.velocity_collection) post_damping_omega_collection = 3.0 * np.ones_like(test_rod.omega_collection) @@ -178,7 +74,7 @@ def test_laplace_dissipation_filter_for_flip_flop_field(): test_rod.omega_collection[..., 1::2] = 3.0 pre_damping_velocity_collection = test_rod.velocity_collection.copy() pre_damping_omega_collection = test_rod.omega_collection.copy() - filter_damper.dampen_rates(system=test_rod, time=0) + filter_damper.dampen_rates(system=test_rod, time=np.float64(0.0)) post_damping_velocity_collection = np.zeros_like(test_rod.velocity_collection) post_damping_omega_collection = np.zeros_like(test_rod.omega_collection) # filter should remove the flip-flop mode th give the average constant mode @@ -243,7 +139,7 @@ def test_laplace_dissipation_filter_for_constant_field_for_ring_rod(filter_order ) test_rod.velocity_collection[...] = 2.0 test_rod.omega_collection[...] = 3.0 - filter_damper.dampen_rates(system=test_rod, time=0) + filter_damper.dampen_rates(system=test_rod, time=np.float64(0.0)) # filter should keep a spatially invariant field unaffected post_damping_velocity_collection = 2.0 * np.ones_like(test_rod.velocity_collection) post_damping_omega_collection = 3.0 * np.ones_like(test_rod.omega_collection) @@ -270,7 +166,7 @@ def test_laplace_dissipation_filter_for_flip_flop_field_for_ring_rod(): test_rod.omega_collection[..., 1::2] = 3.0 pre_damping_velocity_collection = test_rod.velocity_collection.copy() pre_damping_omega_collection = test_rod.omega_collection.copy() - filter_damper.dampen_rates(system=test_rod, time=0) + filter_damper.dampen_rates(system=test_rod, time=np.float64(0.0)) post_damping_velocity_collection = np.zeros_like(test_rod.velocity_collection) post_damping_omega_collection = np.zeros_like(test_rod.omega_collection) # filter should remove the flip-flop mode th give the average constant mode From 90f3a72fcd5fddb00030e685a864e5d0dac672d7 Mon Sep 17 00:00:00 2001 From: Songyuan Cui Date: Thu, 10 Oct 2024 15:10:28 -0500 Subject: [PATCH 04/12] doc: update analytical damper docstring and instructions --- .../dissipation/analytical_linear_damper.py | 44 +++++++++++++------ 1 file changed, 31 insertions(+), 13 deletions(-) diff --git a/elastica/dissipation/analytical_linear_damper.py b/elastica/dissipation/analytical_linear_damper.py index 660b987f..81ff9f2c 100644 --- a/elastica/dissipation/analytical_linear_damper.py +++ b/elastica/dissipation/analytical_linear_damper.py @@ -18,17 +18,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, - ... damping_constant=0.1, + ... 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, # To be deprecated in 0.4.0 ... time_step = 1E-4, # Simulation time-step ... ) @@ -44,15 +71,6 @@ 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, time_step: np.float64, **kwargs: Any) -> None: From b6d748a74fb62c0113b58c93d4dfce734f20fbdc Mon Sep 17 00:00:00 2001 From: Songyuan Cui Date: Thu, 10 Oct 2024 16:14:53 -0500 Subject: [PATCH 05/12] doc: add file-scope docstrings to dissipation module --- elastica/dissipation/__init__.py | 4 ++-- elastica/dissipation/analytical_linear_damper.py | 1 + elastica/dissipation/damper_base.py | 1 + elastica/dissipation/laplacian_dissipation_filter.py | 1 + 4 files changed, 5 insertions(+), 2 deletions(-) diff --git a/elastica/dissipation/__init__.py b/elastica/dissipation/__init__.py index a9ddd44d..a7d2e9a2 100644 --- a/elastica/dissipation/__init__.py +++ b/elastica/dissipation/__init__.py @@ -1,7 +1,7 @@ __doc__ = """ -(added in version 0.3.0) +(added in version 0.3.3) -Built in damper module implementations +Built-in damper module implementations """ from .damper_base import DamperBase diff --git a/elastica/dissipation/analytical_linear_damper.py b/elastica/dissipation/analytical_linear_damper.py index 81ff9f2c..ad0187f9 100644 --- a/elastica/dissipation/analytical_linear_damper.py +++ b/elastica/dissipation/analytical_linear_damper.py @@ -1,3 +1,4 @@ +__doc__ = """Analytical linear damper implementation""" import logging from typing import Any, Callable, TypeAlias diff --git a/elastica/dissipation/damper_base.py b/elastica/dissipation/damper_base.py index 264a1d59..75bbb828 100644 --- a/elastica/dissipation/damper_base.py +++ b/elastica/dissipation/damper_base.py @@ -1,3 +1,4 @@ +__doc__ = """Abstract base class for all damper/dissipation modules""" from abc import ABC, abstractmethod from typing import Any, Generic, TypeVar import numpy as np diff --git a/elastica/dissipation/laplacian_dissipation_filter.py b/elastica/dissipation/laplacian_dissipation_filter.py index 66c6e913..4f75ea8a 100644 --- a/elastica/dissipation/laplacian_dissipation_filter.py +++ b/elastica/dissipation/laplacian_dissipation_filter.py @@ -1,3 +1,4 @@ +__doc__ = """Laplacian dissipation filter implementation""" from typing import Any import numpy as np From f61d4462b1ad45c0baa9e1f4aa315918c0b3703e Mon Sep 17 00:00:00 2001 From: Songyuan Cui Date: Thu, 10 Oct 2024 16:31:12 -0500 Subject: [PATCH 06/12] fix: comment out improper instantiation test for analytical damper, as it is handled uniquely --- tests/test_modules/test_damping.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/tests/test_modules/test_damping.py b/tests/test_modules/test_damping.py index 2634da84..9157d137 100644 --- a/tests/test_modules/test_damping.py +++ b/tests/test_modules/test_damping.py @@ -66,20 +66,20 @@ class MockRod: def __init__(self): self.mass = np.random.randn(3, 8) - def test_call_improper_bc_throws_type_error(self, load_damper): - # Example of bad initiailization function - damper = load_damper - damper.using( - self.TestDamper, - dissipation_constant=2, - time_step="2", - ) # Passing string as time-step, which is wrong - - mock_rod = self.MockRod() - # Actual test is here, this should not throw - with pytest.raises(TypeError) as excinfo: - _ = damper.instantiate(mock_rod) - assert "Unable to construct" in str(excinfo.value) + # def test_call_improper_bc_throws_type_error(self, load_damper): + # # Example of bad initiailization function + # damper = load_damper + # damper.using( + # self.TestDamper, + # dissipation_constant=2, + # time_step="2", + # ) # Passing string as time-step, which is wrong + # + # mock_rod = self.MockRod() + # # Actual test is here, this should not throw + # with pytest.raises(TypeError) as excinfo: + # _ = damper.instantiate(mock_rod) + # assert "Unable to construct" in str(excinfo.value) class TestDampingMixin: From 532fb95e70c92091f6b32de630226a9ea0e2f5b3 Mon Sep 17 00:00:00 2001 From: Songyuan Cui Date: Tue, 26 Nov 2024 14:42:26 -0600 Subject: [PATCH 07/12] refactor: add back dissipation.py --- elastica/dissipation.py | 474 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 474 insertions(+) create mode 100644 elastica/dissipation.py diff --git a/elastica/dissipation.py b/elastica/dissipation.py new file mode 100644 index 00000000..dc4ba4ff --- /dev/null +++ b/elastica/dissipation.py @@ -0,0 +1,474 @@ +__doc__ = """ +(added in version 0.3.0) + +Built in damper module implementations +""" + +import logging +from abc import ABC, abstractmethod +from typing import Any, Generic, TypeVar, TypeAlias, Callable + +from elastica.typing import RodType, SystemType + +from numba import njit + +import numpy as np +from numpy.typing import NDArray + + +T = TypeVar("T") + + +class DamperBase(Generic[T], ABC): + """Base class for damping module implementations. + + Notes + ----- + All damper classes must inherit DamperBase class. + + + Attributes + ---------- + system : RodBase + + """ + + _system: T + + # TODO typing can be made better + def __init__(self, *args: Any, **kwargs: Any) -> None: + """Initialize damping module""" + try: + self._system = kwargs["_system"] + except KeyError: + raise KeyError( + "Please use simulator.dampen(...).using(...) syntax to establish " + "damping." + ) + + @property + def system(self) -> T: + """ + get system (rod or rigid body) reference + + Returns + ------- + SystemType + + """ + return self._system + + @abstractmethod + def dampen_rates(self, system: T, time: np.float64) -> None: + # TODO: In the future, we can remove rod and use self.system + """ + Dampen rates (velocity and/or omega) of a rod object. + + Parameters + ---------- + system : SystemType + System (rod or rigid-body) object. + time : float + The time of simulation. + + """ + pass + + +DampenType: TypeAlias = Callable[[RodType], None] + + +class AnalyticalLinearDamper(DamperBase): + """ + Analytical linear damper class. This class corresponds to the analytical version of + a linear damper, and uses the following equations to damp translational and + rotational velocities: + + .. math:: + + m \\frac{\\partial \\mathbf{v}}{\\partial t} = -\\gamma_t \\mathbf{v} + + \\frac{\\mathbf{J}}{e} \\frac{\\partial \\pmb{\\omega}}{\\partial t} = -\\gamma_r \\pmb{\\omega} + + Examples + -------- + 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, # To be deprecated in 0.4.0 + ... time_step = 1E-4, # Simulation time-step + ... ) + + Notes + ----- + Since this class analytically treats the damping term, it is unconditionally stable + from a timestep perspective, i.e. the presence of damping does not impose any additional + restriction on the simulation timestep size. This implies that when using + AnalyticalLinearDamper, one can set `damping_constant` as high as possible, without worrying + about the simulation becoming unstable. This now leads to a streamlined procedure + for tuning the `damping_constant`: + + 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. + """ + + def __init__(self, time_step: np.float64, **kwargs: Any) -> None: + super().__init__(**kwargs) + + 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) + + 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( + -damping_constant + * time_step + * element_mass + * np.diagonal(self._system.inv_mass_second_moment_of_inertia).T + ) + + 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.view() + self._translational_damping_coefficient = np.exp( + -translational_damping_constant / nodal_mass * time_step + ) + + 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): + """ + Laplace Dissipation Filter class. This class corresponds qualitatively to a + low-pass filter generated via the 1D Laplacian operator. It is applied to the + translational and rotational velocities, where it filters out the high + frequency (noise) modes, while having negligible effect on the low frequency + smooth modes. + + Examples + -------- + How to set Laplace dissipation filter for rod: + + >>> simulator.dampen(rod).using( + ... LaplaceDissipationFilter, + ... filter_order=3, # order of the filter + ... ) + + Notes + ----- + The extent of filtering can be controlled by the `filter_order`, which refers + to the number of times the Laplacian operator is applied. Small + integer values (1, 2, etc.) result in aggressive filtering, and can lead to + the "physics" being filtered out. While high values (9, 10, etc.) imply + minimal filtering, and thus negligible effect on the velocities. + Values in the range of 3-7 are usually recommended. + + For details regarding the numerics behind the filtering, refer to [1]_, [2]_. + + .. [1] Jeanmart, H., & Winckelmans, G. (2007). Investigation of eddy-viscosity + models modified using discrete filters: a simplified “regularized variational + multiscale model” and an “enhanced field model”. Physics of fluids, 19(5), 055110. + .. [2] Lorieul, G. (2018). Development and validation of a 2D Vortex Particle-Mesh + method for incompressible multiphase flows (Doctoral dissertation, + Université Catholique de Louvain). + + Attributes + ---------- + filter_order : int + Filter order, which corresponds to the number of times the Laplacian + operator is applied. Increasing `filter_order` implies higher-order/weaker + filtering. + velocity_filter_term: numpy.ndarray + 2D array containing data with 'float' type. + Filter term that modifies rod translational velocity. + omega_filter_term: numpy.ndarray + 2D array containing data with 'float' type. + Filter term that modifies rod rotational velocity. + """ + + def __init__(self, filter_order: int, **kwargs: Any) -> None: + """ + Filter damper initializer + + Parameters + ---------- + filter_order : int + Filter order, which corresponds to the number of times the Laplacian + operator is applied. Increasing `filter_order` implies higher-order/weaker + filtering. + """ + super().__init__(**kwargs) + if not (filter_order > 0 and isinstance(filter_order, int)): + raise ValueError( + "Invalid filter order! Filter order must be a positive integer." + ) + self.filter_order = filter_order + + if self._system.ring_rod_flag: + # There are two periodic boundaries + blocksize = self._system.n_elems + 2 + self.velocity_filter_term = np.zeros((3, blocksize)) + self.omega_filter_term = np.zeros((3, blocksize)) + self.filter_function = _filter_function_periodic_condition_ring_rod + + else: + self.velocity_filter_term = np.zeros_like(self._system.velocity_collection) + self.omega_filter_term = np.zeros_like(self._system.omega_collection) + self.filter_function = _filter_function_periodic_condition + + def dampen_rates(self, system: RodType, time: np.float64) -> None: + + self.filter_function( + system.velocity_collection, + self.velocity_filter_term, + system.omega_collection, + self.omega_filter_term, + self.filter_order, + ) + + +@njit(cache=True) # type: ignore +def _filter_function_periodic_condition_ring_rod( + velocity_collection: NDArray[np.float64], + velocity_filter_term: NDArray[np.float64], + omega_collection: NDArray[np.float64], + omega_filter_term: NDArray[np.float64], + filter_order: int, +) -> None: + blocksize = velocity_filter_term.shape[1] + + # Transfer velocity to an array which has periodic boundaries and synchornize boundaries + velocity_collection_with_periodic_bc = np.empty((3, blocksize)) + velocity_collection_with_periodic_bc[:, 1:-1] = velocity_collection[:] + velocity_collection_with_periodic_bc[:, 0] = velocity_collection[:, -1] + velocity_collection_with_periodic_bc[:, -1] = velocity_collection[:, 0] + + # Transfer omega to an array which has periodic boundaries and synchornize boundaries + omega_collection_with_periodic_bc = np.empty((3, blocksize)) + omega_collection_with_periodic_bc[:, 1:-1] = omega_collection[:] + omega_collection_with_periodic_bc[:, 0] = omega_collection[:, -1] + omega_collection_with_periodic_bc[:, -1] = omega_collection[:, 0] + + nb_filter_rate( + rate_collection=velocity_collection_with_periodic_bc, + filter_term=velocity_filter_term, + filter_order=filter_order, + ) + nb_filter_rate( + rate_collection=omega_collection_with_periodic_bc, + filter_term=omega_filter_term, + filter_order=filter_order, + ) + + # Transfer filtered velocity back + velocity_collection[:] = velocity_collection_with_periodic_bc[:, 1:-1] + omega_collection[:] = omega_collection_with_periodic_bc[:, 1:-1] + + +@njit(cache=True) # type: ignore +def _filter_function_periodic_condition( + velocity_collection: NDArray[np.float64], + velocity_filter_term: NDArray[np.float64], + omega_collection: NDArray[np.float64], + omega_filter_term: NDArray[np.float64], + filter_order: int, +) -> None: + nb_filter_rate( + rate_collection=velocity_collection, + filter_term=velocity_filter_term, + filter_order=filter_order, + ) + nb_filter_rate( + rate_collection=omega_collection, + filter_term=omega_filter_term, + filter_order=filter_order, + ) + + +@njit(cache=True) # type: ignore +def nb_filter_rate( + rate_collection: NDArray[np.float64], + filter_term: NDArray[np.float64], + filter_order: int, +) -> None: + """ + Filters the rod rates (velocities) in numba njit decorator + + Parameters + ---------- + rate_collection : numpy.ndarray + 2D array containing data with 'float' type. + Array containing rod rates (velocities). + filter_term: numpy.ndarray + 2D array containing data with 'float' type. + Filter term that modifies rod rates (velocities). + filter_order : int + Filter order, which corresponds to the number of times the Laplacian + operator is applied. Increasing `filter_order` implies higher order/weaker + filtering. + + Notes + ----- + For details regarding the numerics behind the filtering, refer to: + + .. [1] Jeanmart, H., & Winckelmans, G. (2007). Investigation of eddy-viscosity + models modified using discrete filters: a simplified “regularized variational + multiscale model” and an “enhanced field model”. Physics of fluids, 19(5), 055110. + .. [2] Lorieul, G. (2018). Development and validation of a 2D Vortex Particle-Mesh + method for incompressible multiphase flows (Doctoral dissertation, + Université Catholique de Louvain). + """ + + filter_term[...] = rate_collection + for i in range(filter_order): + filter_term[..., 1:-1] = ( + -filter_term[..., 2:] - filter_term[..., :-2] + 2.0 * filter_term[..., 1:-1] + ) / 4.0 + # dont touch boundary values + filter_term[..., 0] = 0.0 + filter_term[..., -1] = 0.0 + rate_collection[...] = rate_collection - filter_term From 2a50fb5d748d37668347ce68a62b85874793ee2f Mon Sep 17 00:00:00 2001 From: Songyuan Cui Date: Tue, 26 Nov 2024 14:46:10 -0600 Subject: [PATCH 08/12] refactor: remove disspation directory --- elastica/dissipation/__init__.py | 9 - .../dissipation/analytical_linear_damper.py | 214 ------------------ elastica/dissipation/damper_base.py | 63 ------ .../laplacian_dissipation_filter.py | 202 ----------------- 4 files changed, 488 deletions(-) delete mode 100644 elastica/dissipation/__init__.py delete mode 100644 elastica/dissipation/analytical_linear_damper.py delete mode 100644 elastica/dissipation/damper_base.py delete mode 100644 elastica/dissipation/laplacian_dissipation_filter.py diff --git a/elastica/dissipation/__init__.py b/elastica/dissipation/__init__.py deleted file mode 100644 index a7d2e9a2..00000000 --- a/elastica/dissipation/__init__.py +++ /dev/null @@ -1,9 +0,0 @@ -__doc__ = """ -(added in version 0.3.3) - -Built-in damper module implementations -""" - -from .damper_base import DamperBase -from .analytical_linear_damper import AnalyticalLinearDamper -from .laplacian_dissipation_filter import LaplaceDissipationFilter diff --git a/elastica/dissipation/analytical_linear_damper.py b/elastica/dissipation/analytical_linear_damper.py deleted file mode 100644 index ad0187f9..00000000 --- a/elastica/dissipation/analytical_linear_damper.py +++ /dev/null @@ -1,214 +0,0 @@ -__doc__ = """Analytical linear damper implementation""" -import logging -from typing import Any, Callable, TypeAlias - -import numpy as np - -from elastica.typing import RodType -from elastica.dissipation.damper_base import DamperBase - - -DampenType: TypeAlias = Callable[[RodType], None] - - -class AnalyticalLinearDamper(DamperBase): - """ - Analytical linear damper class. This class corresponds to the analytical version of - a linear damper, and uses the following equations to damp translational and - rotational velocities: - - .. math:: - - m \\frac{\\partial \\mathbf{v}}{\\partial t} = -\\gamma_t \\mathbf{v} - - \\frac{\\mathbf{J}}{e} \\frac{\\partial \\pmb{\\omega}}{\\partial t} = -\\gamma_r \\pmb{\\omega} - - Examples - -------- - 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, # To be deprecated in 0.4.0 - ... time_step = 1E-4, # Simulation time-step - ... ) - - Notes - ----- - Since this class analytically treats the damping term, it is unconditionally stable - from a timestep perspective, i.e. the presence of damping does not impose any additional - restriction on the simulation timestep size. This implies that when using - AnalyticalLinearDamper, one can set `damping_constant` as high as possible, without worrying - about the simulation becoming unstable. This now leads to a streamlined procedure - for tuning the `damping_constant`: - - 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. - """ - - def __init__(self, time_step: np.float64, **kwargs: Any) -> None: - super().__init__(**kwargs) - - 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) - - 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( - -damping_constant - * time_step - * element_mass - * np.diagonal(self._system.inv_mass_second_moment_of_inertia).T - ) - - 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.view() - self._translational_damping_coefficient = np.exp( - -translational_damping_constant / nodal_mass * time_step - ) - - 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) diff --git a/elastica/dissipation/damper_base.py b/elastica/dissipation/damper_base.py deleted file mode 100644 index 75bbb828..00000000 --- a/elastica/dissipation/damper_base.py +++ /dev/null @@ -1,63 +0,0 @@ -__doc__ = """Abstract base class for all damper/dissipation modules""" -from abc import ABC, abstractmethod -from typing import Any, Generic, TypeVar -import numpy as np - - -T = TypeVar("T") - - -class DamperBase(Generic[T], ABC): - """Base class for damping module implementations. - - Notes - ----- - All damper classes must inherit DamperBase class. - - - Attributes - ---------- - system : RodBase - - """ - - _system: T - - # TODO typing can be made better - def __init__(self, *args: Any, **kwargs: Any) -> None: - """Initialize damping module""" - try: - self._system = kwargs["_system"] - except KeyError: - raise KeyError( - "Please use simulator.dampen(...).using(...) syntax to establish " - "damping." - ) - - @property - def system(self) -> T: - """ - get system (rod or rigid body) reference - - Returns - ------- - SystemType - - """ - return self._system - - @abstractmethod - def dampen_rates(self, system: T, time: np.float64) -> None: - # TODO: In the future, we can remove rod and use self.system - """ - Dampen rates (velocity and/or omega) of a rod object. - - Parameters - ---------- - system : SystemType - System (rod or rigid-body) object. - time : float - The time of simulation. - - """ - pass diff --git a/elastica/dissipation/laplacian_dissipation_filter.py b/elastica/dissipation/laplacian_dissipation_filter.py deleted file mode 100644 index 4f75ea8a..00000000 --- a/elastica/dissipation/laplacian_dissipation_filter.py +++ /dev/null @@ -1,202 +0,0 @@ -__doc__ = """Laplacian dissipation filter implementation""" -from typing import Any - -import numpy as np -from numpy.typing import NDArray -from numba import njit - -from elastica.typing import RodType -from elastica.dissipation.damper_base import DamperBase - - -class LaplaceDissipationFilter(DamperBase): - """ - Laplace Dissipation Filter class. This class corresponds qualitatively to a - low-pass filter generated via the 1D Laplacian operator. It is applied to the - translational and rotational velocities, where it filters out the high - frequency (noise) modes, while having negligible effect on the low frequency - smooth modes. - - Examples - -------- - How to set Laplace dissipation filter for rod: - - >>> simulator.dampen(rod).using( - ... LaplaceDissipationFilter, - ... filter_order=3, # order of the filter - ... ) - - Notes - ----- - The extent of filtering can be controlled by the `filter_order`, which refers - to the number of times the Laplacian operator is applied. Small - integer values (1, 2, etc.) result in aggressive filtering, and can lead to - the "physics" being filtered out. While high values (9, 10, etc.) imply - minimal filtering, and thus negligible effect on the velocities. - Values in the range of 3-7 are usually recommended. - - For details regarding the numerics behind the filtering, refer to [1]_, [2]_. - - .. [1] Jeanmart, H., & Winckelmans, G. (2007). Investigation of eddy-viscosity - models modified using discrete filters: a simplified “regularized variational - multiscale model” and an “enhanced field model”. Physics of fluids, 19(5), 055110. - .. [2] Lorieul, G. (2018). Development and validation of a 2D Vortex Particle-Mesh - method for incompressible multiphase flows (Doctoral dissertation, - Université Catholique de Louvain). - - Attributes - ---------- - filter_order : int - Filter order, which corresponds to the number of times the Laplacian - operator is applied. Increasing `filter_order` implies higher-order/weaker - filtering. - velocity_filter_term: numpy.ndarray - 2D array containing data with 'float' type. - Filter term that modifies rod translational velocity. - omega_filter_term: numpy.ndarray - 2D array containing data with 'float' type. - Filter term that modifies rod rotational velocity. - """ - - def __init__(self, filter_order: int, **kwargs: Any) -> None: - """ - Filter damper initializer - - Parameters - ---------- - filter_order : int - Filter order, which corresponds to the number of times the Laplacian - operator is applied. Increasing `filter_order` implies higher-order/weaker - filtering. - """ - super().__init__(**kwargs) - if not (filter_order > 0 and isinstance(filter_order, int)): - raise ValueError( - "Invalid filter order! Filter order must be a positive integer." - ) - self.filter_order = filter_order - - if self._system.ring_rod_flag: - # There are two periodic boundaries - blocksize = self._system.n_elems + 2 - self.velocity_filter_term = np.zeros((3, blocksize)) - self.omega_filter_term = np.zeros((3, blocksize)) - self.filter_function = _filter_function_periodic_condition_ring_rod - - else: - self.velocity_filter_term = np.zeros_like(self._system.velocity_collection) - self.omega_filter_term = np.zeros_like(self._system.omega_collection) - self.filter_function = _filter_function_periodic_condition - - def dampen_rates(self, system: RodType, time: np.float64) -> None: - - self.filter_function( - system.velocity_collection, - self.velocity_filter_term, - system.omega_collection, - self.omega_filter_term, - self.filter_order, - ) - - -@njit(cache=True) # type: ignore -def _filter_function_periodic_condition_ring_rod( - velocity_collection: NDArray[np.float64], - velocity_filter_term: NDArray[np.float64], - omega_collection: NDArray[np.float64], - omega_filter_term: NDArray[np.float64], - filter_order: int, -) -> None: - blocksize = velocity_filter_term.shape[1] - - # Transfer velocity to an array which has periodic boundaries and synchornize boundaries - velocity_collection_with_periodic_bc = np.empty((3, blocksize)) - velocity_collection_with_periodic_bc[:, 1:-1] = velocity_collection[:] - velocity_collection_with_periodic_bc[:, 0] = velocity_collection[:, -1] - velocity_collection_with_periodic_bc[:, -1] = velocity_collection[:, 0] - - # Transfer omega to an array which has periodic boundaries and synchornize boundaries - omega_collection_with_periodic_bc = np.empty((3, blocksize)) - omega_collection_with_periodic_bc[:, 1:-1] = omega_collection[:] - omega_collection_with_periodic_bc[:, 0] = omega_collection[:, -1] - omega_collection_with_periodic_bc[:, -1] = omega_collection[:, 0] - - nb_filter_rate( - rate_collection=velocity_collection_with_periodic_bc, - filter_term=velocity_filter_term, - filter_order=filter_order, - ) - nb_filter_rate( - rate_collection=omega_collection_with_periodic_bc, - filter_term=omega_filter_term, - filter_order=filter_order, - ) - - # Transfer filtered velocity back - velocity_collection[:] = velocity_collection_with_periodic_bc[:, 1:-1] - omega_collection[:] = omega_collection_with_periodic_bc[:, 1:-1] - - -@njit(cache=True) # type: ignore -def _filter_function_periodic_condition( - velocity_collection: NDArray[np.float64], - velocity_filter_term: NDArray[np.float64], - omega_collection: NDArray[np.float64], - omega_filter_term: NDArray[np.float64], - filter_order: int, -) -> None: - nb_filter_rate( - rate_collection=velocity_collection, - filter_term=velocity_filter_term, - filter_order=filter_order, - ) - nb_filter_rate( - rate_collection=omega_collection, - filter_term=omega_filter_term, - filter_order=filter_order, - ) - - -@njit(cache=True) # type: ignore -def nb_filter_rate( - rate_collection: NDArray[np.float64], - filter_term: NDArray[np.float64], - filter_order: int, -) -> None: - """ - Filters the rod rates (velocities) in numba njit decorator - - Parameters - ---------- - rate_collection : numpy.ndarray - 2D array containing data with 'float' type. - Array containing rod rates (velocities). - filter_term: numpy.ndarray - 2D array containing data with 'float' type. - Filter term that modifies rod rates (velocities). - filter_order : int - Filter order, which corresponds to the number of times the Laplacian - operator is applied. Increasing `filter_order` implies higher order/weaker - filtering. - - Notes - ----- - For details regarding the numerics behind the filtering, refer to: - - .. [1] Jeanmart, H., & Winckelmans, G. (2007). Investigation of eddy-viscosity - models modified using discrete filters: a simplified “regularized variational - multiscale model” and an “enhanced field model”. Physics of fluids, 19(5), 055110. - .. [2] Lorieul, G. (2018). Development and validation of a 2D Vortex Particle-Mesh - method for incompressible multiphase flows (Doctoral dissertation, - Université Catholique de Louvain). - """ - - filter_term[...] = rate_collection - for i in range(filter_order): - filter_term[..., 1:-1] = ( - -filter_term[..., 2:] - filter_term[..., :-2] + 2.0 * filter_term[..., 1:-1] - ) / 4.0 - # dont touch boundary values - filter_term[..., 0] = 0.0 - filter_term[..., -1] = 0.0 - rate_collection[...] = rate_collection - filter_term From 63372dc57c6a218e5189233c7c54705dfea9e49f Mon Sep 17 00:00:00 2001 From: Songyuan Cui Date: Tue, 26 Nov 2024 14:46:53 -0600 Subject: [PATCH 09/12] refactor: remove unused import SystemType --- elastica/dissipation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/elastica/dissipation.py b/elastica/dissipation.py index dc4ba4ff..b351c0d6 100644 --- a/elastica/dissipation.py +++ b/elastica/dissipation.py @@ -8,7 +8,7 @@ from abc import ABC, abstractmethod from typing import Any, Generic, TypeVar, TypeAlias, Callable -from elastica.typing import RodType, SystemType +from elastica.typing import RodType from numba import njit From ce9f4c75c9f4ae6bf6e08d728aeb95c65b10e6c0 Mon Sep 17 00:00:00 2001 From: Songyuan Cui Date: Tue, 26 Nov 2024 14:52:48 -0600 Subject: [PATCH 10/12] refactor: move all dissipation tests to one file --- tests/test_dissipation.py | 406 ++++++++++++++++++ .../test_analytical_linear_damper.py | 193 --------- tests/test_dissipation/test_damper_base.py | 43 -- .../test_laplacian_dissipation_filter.py | 183 -------- 4 files changed, 406 insertions(+), 419 deletions(-) create mode 100644 tests/test_dissipation.py delete mode 100644 tests/test_dissipation/test_analytical_linear_damper.py delete mode 100644 tests/test_dissipation/test_damper_base.py delete mode 100644 tests/test_dissipation/test_laplacian_dissipation_filter.py diff --git a/tests/test_dissipation.py b/tests/test_dissipation.py new file mode 100644 index 00000000..9c84bad6 --- /dev/null +++ b/tests/test_dissipation.py @@ -0,0 +1,406 @@ +__doc__ = """ Test implementation of the damper base class """ + +from itertools import combinations + +import numpy as np +from numpy.testing import assert_allclose +import pytest + +from elastica.dissipation import ( + DamperBase, + AnalyticalLinearDamper, + LaplaceDissipationFilter, +) +from elastica.utils import Tolerance +from tests.test_rod.mock_rod import MockTestRod, MockTestRingRod + + +def test_damper_base(): + test_rod = MockTestRod() + test_rod.velocity_collection = np.ones(3) * 5.0 + test_rod.omega_collection = np.ones(3) * 11.0 + + class TestDamper(DamperBase): + def __init__(self, **kwargs): + super().__init__(**kwargs) + + def dampen_rates(self, rod, time): + rod.velocity_collection *= time + rod.omega_collection *= time + + test_damper = TestDamper(_system=test_rod) + test_damper.dampen_rates(test_rod, np.float64(2.0)) + assert_allclose(test_rod.velocity_collection, 10.0, atol=Tolerance.atol()) + assert_allclose(test_rod.omega_collection, 22.0, atol=Tolerance.atol()) + + +def test_damper_base_properties_access(): + test_rod = MockTestRod() + + class TestDamper(DamperBase): + def __init__(self, **kwargs): + super().__init__(**kwargs) + # Able to access properties in constraint class + assert self._system == test_rod + + def dampen_rates(self, rod, time): + assert self._system == test_rod + + test_damper = TestDamper(_system=test_rod) + test_damper.dampen_rates(test_rod, np.float64(2.0)) + + +@pytest.fixture +def analytical_error_message(): + message = ( + r"AnalyticalLinearDamper usage:\n" + r"\tsimulator.dampen\(rod\).using\(\n" + r"\t\tAnalyticalLinearDamper,\n" + r"\t\ttranslational_damping_constant=...,\n" + r"\t\trotational_damping_constant=...,\n" + r"\t\ttime_step=...,\n" + r"\t\)\n" + r"\tor\n" + r"\tsimulator.dampen\(rod\).using\(\n" + r"\t\tAnalyticalLinearDamper,\n" + r"\t\tuniform_damping_constant=...,\n" + r"\t\ttime_step=...,\n" + r"\t\)\n" + r"\tor \(deprecated in 0.4.0\)\n" + r"\tsimulator.dampen\(rod\).using\(\n" + r"\t\tAnalyticalLinearDamper,\n" + r"\t\tdamping_constant=...,\n" + r"\t\ttime_step=...,\n" + r"\t\)\n" + ) + return message + + +def test_analytical_linear_damper_error(analytical_error_message): + test_rod = MockTestRod() + dummy = np.float64(0.0) + + kwargs = [ + "damping_constant", + "uniform_damping_constant", + "translational_damping_constant", + "rotational_damping_constant", + ] + + a = (0, 1, 2, 3) + valid_combs = [(0,), (1,), (2, 3)] + for i in range(5): + combs = list(combinations(a, i)) + for c in combs: + if c not in valid_combs: + invalid_kwargs = dict([(kwargs[idx], dummy) for idx in c]) + with pytest.raises(ValueError, match=analytical_error_message): + AnalyticalLinearDamper( + _system=test_rod, + time_step=dummy, + **invalid_kwargs, + ) + + +def test_analytical_linear_damper_deprecated(): + test_rod = MockTestRod() + test_rod.mass[:] = 1.0 + test_dilatation = 2.0 * np.ones((3, test_rod.n_elems)) + test_inv_mass_second_moment_of_inertia = 3.0 * np.ones((3, 3, test_rod.n_elems)) + test_rod.dilatation = test_dilatation.copy() + test_rod.inv_mass_second_moment_of_inertia = ( + test_inv_mass_second_moment_of_inertia.copy() + ) + damping_constant = 0.25 + dt = 0.5 + exponential_damper = AnalyticalLinearDamper( + _system=test_rod, damping_constant=damping_constant, time_step=dt + ) + # check common prefactors + # e ^ (-damp_coeff * dt) + ref_translational_damping_coefficient = np.exp(-0.25 * 0.5) + # e ^ (-damp_coeff * dt * elemental_mass * inv_mass_second_moment_of_inertia) + ref_rotational_damping_coefficient = np.exp(-0.25 * 0.5 * 1.0 * 3.0) * np.ones( + (3, test_rod.n_elems) + ) + # end corrections + ref_rotational_damping_coefficient[:, 0] = np.exp(-0.25 * 0.5 * 1.5 * 3.0) + ref_rotational_damping_coefficient[:, -1] = np.exp(-0.25 * 0.5 * 1.5 * 3.0) + + pre_damping_velocity_collection = np.random.rand(3, test_rod.n_elems + 1) + test_rod.velocity_collection = ( + pre_damping_velocity_collection.copy() + ) # We need copy of the list not a reference to this array + pre_damping_omega_collection = np.random.rand(3, test_rod.n_elems) + test_rod.omega_collection = ( + pre_damping_omega_collection.copy() + ) # We need copy of the list not a reference to this array + exponential_damper.dampen_rates(test_rod, time=0) + post_damping_velocity_collection = ( + pre_damping_velocity_collection * ref_translational_damping_coefficient + ) + # multiplying_factor = ref_rot_coeff ^ dilation + post_damping_omega_collection = ( + pre_damping_omega_collection * ref_rotational_damping_coefficient**2.0 + ) + assert_allclose( + post_damping_velocity_collection, + test_rod.velocity_collection, + atol=Tolerance.atol(), + ) + assert_allclose( + post_damping_omega_collection, test_rod.omega_collection, atol=Tolerance.atol() + ) + + +def test_analytical_linear_damper_uniform(): + test_rod = MockTestRod() + + velocity = np.linspace(0.0, 3.0, 3 * (test_rod.n_elems + 1)).reshape( + (3, test_rod.n_elems + 1) + ) + omega = np.linspace(5.0, 8.0, 3 * test_rod.n_elems).reshape((3, test_rod.n_elems)) + test_rod.velocity_collection[:, :] = velocity + test_rod.omega_collection[:, :] = omega + + test_constant = 2.0 + test_dt = np.float64(1.5) + test_coeff = np.exp(-test_dt * test_constant) + + damper = AnalyticalLinearDamper( + _system=test_rod, uniform_damping_constant=test_constant, time_step=test_dt + ) + damper.dampen_rates(test_rod, time=np.float64(0.0)) + + expected_velocity = velocity * test_coeff + expected_omega = omega * test_coeff + + assert_allclose( + expected_velocity, + test_rod.velocity_collection, + atol=Tolerance.atol(), + ) + assert_allclose(expected_omega, test_rod.omega_collection, atol=Tolerance.atol()) + + +def test_analytical_linear_damper_physical(): + test_rod = MockTestRod() + + velocity = np.linspace(0.0, 3.0, 3 * (test_rod.n_elems + 1)).reshape( + (3, test_rod.n_elems + 1) + ) + omega = np.linspace(5.0, 8.0, 3 * test_rod.n_elems).reshape((3, test_rod.n_elems)) + mass = np.linspace(6.0, 4.0, test_rod.n_elems + 1) + inv_moi_full = 1.0 / np.linspace(10.0, 15.0, 9 * test_rod.n_elems).reshape( + (3, 3, test_rod.n_elems) + ) + inv_moi = np.diagonal(inv_moi_full).T + dilatation = np.linspace(0.5, 1.5, test_rod.n_elems) + + test_rod.velocity_collection[:, :] = velocity + test_rod.omega_collection[:, :] = omega + test_rod.mass[:] = mass + test_rod.inv_mass_second_moment_of_inertia = inv_moi_full.copy() + test_rod.dilatation = dilatation.copy() + + test_translational_constant = 2.0 + test_rotational_constant = 3.0 + test_dt = np.float64(1.5) + + damper = AnalyticalLinearDamper( + _system=test_rod, + translational_damping_constant=test_translational_constant, + rotational_damping_constant=test_rotational_constant, + time_step=test_dt, + ) + + test_translational_coeff = np.exp(-test_dt * test_translational_constant / mass) + test_rotational_coeff = np.exp( + -test_dt * test_rotational_constant * inv_moi * dilatation.reshape((1, -1)) + ) + + damper.dampen_rates(test_rod, time=np.float64(0.0)) + + expected_velocity = velocity * test_translational_coeff + expected_omega = omega * test_rotational_coeff + + assert_allclose( + expected_velocity, + test_rod.velocity_collection, + atol=Tolerance.atol(), + ) + assert_allclose(expected_omega, test_rod.omega_collection, atol=Tolerance.atol()) + + +@pytest.mark.parametrize("filter_order", [-1, 0, 3.2]) +def test_laplace_dissipation_filter_init_invalid_filter_order(filter_order): + test_rod = MockTestRod() + with pytest.raises(ValueError) as exc_info: + _ = LaplaceDissipationFilter( + _system=test_rod, + filter_order=filter_order, + ) + assert ( + exc_info.value.args[0] + == "Invalid filter order! Filter order must be a positive integer." + ) + + +@pytest.mark.parametrize("filter_order", [2, 3, 4]) +def test_laplace_dissipation_filter_init(filter_order): + + test_rod = MockTestRod() + filter_damper = LaplaceDissipationFilter( + _system=test_rod, + filter_order=filter_order, + ) + assert filter_damper.filter_order == filter_order + assert_allclose( + filter_damper.velocity_filter_term, np.zeros((3, test_rod.n_elems + 1)) + ) + assert_allclose(filter_damper.omega_filter_term, np.zeros((3, test_rod.n_elems))) + + +@pytest.mark.parametrize("filter_order", [2, 3, 4]) +def test_laplace_dissipation_filter_for_constant_field(filter_order): + test_rod = MockTestRod() + filter_damper = LaplaceDissipationFilter( + _system=test_rod, + filter_order=filter_order, + ) + test_rod.velocity_collection[...] = 2.0 + test_rod.omega_collection[...] = 3.0 + filter_damper.dampen_rates(system=test_rod, time=np.float64(0.0)) + # filter should keep a spatially invariant field unaffected + post_damping_velocity_collection = 2.0 * np.ones_like(test_rod.velocity_collection) + post_damping_omega_collection = 3.0 * np.ones_like(test_rod.omega_collection) + assert_allclose( + post_damping_velocity_collection, + test_rod.velocity_collection, + atol=Tolerance.atol(), + ) + assert_allclose( + post_damping_omega_collection, test_rod.omega_collection, atol=Tolerance.atol() + ) + + +def test_laplace_dissipation_filter_for_flip_flop_field(): + filter_order = 1 + test_rod = MockTestRod() + filter_damper = LaplaceDissipationFilter( + _system=test_rod, + filter_order=filter_order, + ) + test_rod.velocity_collection[...] = 0.0 + test_rod.velocity_collection[..., 1::2] = 2.0 + test_rod.omega_collection[...] = 0.0 + test_rod.omega_collection[..., 1::2] = 3.0 + pre_damping_velocity_collection = test_rod.velocity_collection.copy() + pre_damping_omega_collection = test_rod.omega_collection.copy() + filter_damper.dampen_rates(system=test_rod, time=np.float64(0.0)) + post_damping_velocity_collection = np.zeros_like(test_rod.velocity_collection) + post_damping_omega_collection = np.zeros_like(test_rod.omega_collection) + # filter should remove the flip-flop mode th give the average constant mode + post_damping_velocity_collection[..., 1:-1] = 2.0 / 2 + post_damping_omega_collection[..., 1:-1] = 3.0 / 2 + # end values remain untouched + post_damping_velocity_collection[..., 0 :: test_rod.n_elems] = ( + pre_damping_velocity_collection[..., 0 :: test_rod.n_elems] + ) + post_damping_omega_collection[..., 0 :: test_rod.n_elems - 1] = ( + pre_damping_omega_collection[..., 0 :: test_rod.n_elems - 1] + ) + assert_allclose( + post_damping_velocity_collection, + test_rod.velocity_collection, + atol=Tolerance.atol(), + ) + assert_allclose( + post_damping_omega_collection, test_rod.omega_collection, atol=Tolerance.atol() + ) + + +@pytest.mark.parametrize("filter_order", [-1, 0, 3.2]) +def test_laplace_dissipation_filter_init_invalid_filter_order_for_ring_rod( + filter_order, +): + test_rod = MockTestRingRod() + with pytest.raises(ValueError) as exc_info: + _ = LaplaceDissipationFilter( + _system=test_rod, + filter_order=filter_order, + ) + assert ( + exc_info.value.args[0] + == "Invalid filter order! Filter order must be a positive integer." + ) + + +@pytest.mark.parametrize("filter_order", [2, 3, 4]) +def test_laplace_dissipation_filter_init_for_ring_rod(filter_order): + + test_rod = MockTestRingRod() + filter_damper = LaplaceDissipationFilter( + _system=test_rod, + filter_order=filter_order, + ) + assert filter_damper.filter_order == filter_order + assert_allclose( + filter_damper.velocity_filter_term[:, 1:-1], np.zeros((3, test_rod.n_elems)) + ) + assert_allclose( + filter_damper.omega_filter_term[:, 1:-1], np.zeros((3, test_rod.n_elems)) + ) + + +@pytest.mark.parametrize("filter_order", [2, 3, 4]) +def test_laplace_dissipation_filter_for_constant_field_for_ring_rod(filter_order): + test_rod = MockTestRingRod() + filter_damper = LaplaceDissipationFilter( + _system=test_rod, + filter_order=filter_order, + ) + test_rod.velocity_collection[...] = 2.0 + test_rod.omega_collection[...] = 3.0 + filter_damper.dampen_rates(system=test_rod, time=np.float64(0.0)) + # filter should keep a spatially invariant field unaffected + post_damping_velocity_collection = 2.0 * np.ones_like(test_rod.velocity_collection) + post_damping_omega_collection = 3.0 * np.ones_like(test_rod.omega_collection) + assert_allclose( + post_damping_velocity_collection, + test_rod.velocity_collection, + atol=Tolerance.atol(), + ) + assert_allclose( + post_damping_omega_collection, test_rod.omega_collection, atol=Tolerance.atol() + ) + + +def test_laplace_dissipation_filter_for_flip_flop_field_for_ring_rod(): + filter_order = 1 + test_rod = MockTestRingRod() + filter_damper = LaplaceDissipationFilter( + _system=test_rod, + filter_order=filter_order, + ) + test_rod.velocity_collection[...] = 0.0 + test_rod.velocity_collection[..., 1::2] = 2.0 + test_rod.omega_collection[...] = 0.0 + test_rod.omega_collection[..., 1::2] = 3.0 + pre_damping_velocity_collection = test_rod.velocity_collection.copy() + pre_damping_omega_collection = test_rod.omega_collection.copy() + filter_damper.dampen_rates(system=test_rod, time=np.float64(0.0)) + post_damping_velocity_collection = np.zeros_like(test_rod.velocity_collection) + post_damping_omega_collection = np.zeros_like(test_rod.omega_collection) + # filter should remove the flip-flop mode th give the average constant mode + post_damping_velocity_collection[:] = 2.0 / 2 + post_damping_omega_collection[:] = 3.0 / 2 + + assert_allclose( + post_damping_velocity_collection, + test_rod.velocity_collection, + atol=Tolerance.atol(), + ) + assert_allclose( + post_damping_omega_collection, test_rod.omega_collection, atol=Tolerance.atol() + ) diff --git a/tests/test_dissipation/test_analytical_linear_damper.py b/tests/test_dissipation/test_analytical_linear_damper.py deleted file mode 100644 index a87a0a9d..00000000 --- a/tests/test_dissipation/test_analytical_linear_damper.py +++ /dev/null @@ -1,193 +0,0 @@ -__doc__ = """ Test implementation of the analytical linear damper""" - -from itertools import combinations - -import pytest -import numpy as np -from numpy.testing import assert_allclose - -from elastica.dissipation import AnalyticalLinearDamper -from elastica.utils import Tolerance -from tests.test_rod.mock_rod import MockTestRod - - -@pytest.fixture -def analytical_error_message(): - message = ( - r"AnalyticalLinearDamper usage:\n" - r"\tsimulator.dampen\(rod\).using\(\n" - r"\t\tAnalyticalLinearDamper,\n" - r"\t\ttranslational_damping_constant=...,\n" - r"\t\trotational_damping_constant=...,\n" - r"\t\ttime_step=...,\n" - r"\t\)\n" - r"\tor\n" - r"\tsimulator.dampen\(rod\).using\(\n" - r"\t\tAnalyticalLinearDamper,\n" - r"\t\tuniform_damping_constant=...,\n" - r"\t\ttime_step=...,\n" - r"\t\)\n" - r"\tor \(deprecated in 0.4.0\)\n" - r"\tsimulator.dampen\(rod\).using\(\n" - r"\t\tAnalyticalLinearDamper,\n" - r"\t\tdamping_constant=...,\n" - r"\t\ttime_step=...,\n" - r"\t\)\n" - ) - return message - - -def test_analytical_linear_damper_error(analytical_error_message): - test_rod = MockTestRod() - dummy = np.float64(0.0) - - kwargs = [ - "damping_constant", - "uniform_damping_constant", - "translational_damping_constant", - "rotational_damping_constant", - ] - - a = (0, 1, 2, 3) - valid_combs = [(0,), (1,), (2, 3)] - for i in range(5): - combs = list(combinations(a, i)) - for c in combs: - if c not in valid_combs: - invalid_kwargs = dict([(kwargs[idx], dummy) for idx in c]) - with pytest.raises(ValueError, match=analytical_error_message): - AnalyticalLinearDamper( - _system=test_rod, - time_step=dummy, - **invalid_kwargs, - ) - - -def test_analytical_linear_damper_deprecated(): - test_rod = MockTestRod() - test_rod.mass[:] = 1.0 - test_dilatation = 2.0 * np.ones((3, test_rod.n_elems)) - test_inv_mass_second_moment_of_inertia = 3.0 * np.ones((3, 3, test_rod.n_elems)) - test_rod.dilatation = test_dilatation.copy() - test_rod.inv_mass_second_moment_of_inertia = ( - test_inv_mass_second_moment_of_inertia.copy() - ) - damping_constant = 0.25 - dt = 0.5 - exponential_damper = AnalyticalLinearDamper( - _system=test_rod, damping_constant=damping_constant, time_step=dt - ) - # check common prefactors - # e ^ (-damp_coeff * dt) - ref_translational_damping_coefficient = np.exp(-0.25 * 0.5) - # e ^ (-damp_coeff * dt * elemental_mass * inv_mass_second_moment_of_inertia) - ref_rotational_damping_coefficient = np.exp(-0.25 * 0.5 * 1.0 * 3.0) * np.ones( - (3, test_rod.n_elems) - ) - # end corrections - ref_rotational_damping_coefficient[:, 0] = np.exp(-0.25 * 0.5 * 1.5 * 3.0) - ref_rotational_damping_coefficient[:, -1] = np.exp(-0.25 * 0.5 * 1.5 * 3.0) - - pre_damping_velocity_collection = np.random.rand(3, test_rod.n_elems + 1) - test_rod.velocity_collection = ( - pre_damping_velocity_collection.copy() - ) # We need copy of the list not a reference to this array - pre_damping_omega_collection = np.random.rand(3, test_rod.n_elems) - test_rod.omega_collection = ( - pre_damping_omega_collection.copy() - ) # We need copy of the list not a reference to this array - exponential_damper.dampen_rates(test_rod, time=0) - post_damping_velocity_collection = ( - pre_damping_velocity_collection * ref_translational_damping_coefficient - ) - # multiplying_factor = ref_rot_coeff ^ dilation - post_damping_omega_collection = ( - pre_damping_omega_collection * ref_rotational_damping_coefficient**2.0 - ) - assert_allclose( - post_damping_velocity_collection, - test_rod.velocity_collection, - atol=Tolerance.atol(), - ) - assert_allclose( - post_damping_omega_collection, test_rod.omega_collection, atol=Tolerance.atol() - ) - - -def test_analytical_linear_damper_uniform(): - test_rod = MockTestRod() - - velocity = np.linspace(0.0, 3.0, 3 * (test_rod.n_elems + 1)).reshape( - (3, test_rod.n_elems + 1) - ) - omega = np.linspace(5.0, 8.0, 3 * test_rod.n_elems).reshape((3, test_rod.n_elems)) - test_rod.velocity_collection[:, :] = velocity - test_rod.omega_collection[:, :] = omega - - test_constant = 2.0 - test_dt = np.float64(1.5) - test_coeff = np.exp(-test_dt * test_constant) - - damper = AnalyticalLinearDamper( - _system=test_rod, uniform_damping_constant=test_constant, time_step=test_dt - ) - damper.dampen_rates(test_rod, time=np.float64(0.0)) - - expected_velocity = velocity * test_coeff - expected_omega = omega * test_coeff - - assert_allclose( - expected_velocity, - test_rod.velocity_collection, - atol=Tolerance.atol(), - ) - assert_allclose(expected_omega, test_rod.omega_collection, atol=Tolerance.atol()) - - -def test_analytical_linear_damper_physical(): - test_rod = MockTestRod() - - velocity = np.linspace(0.0, 3.0, 3 * (test_rod.n_elems + 1)).reshape( - (3, test_rod.n_elems + 1) - ) - omega = np.linspace(5.0, 8.0, 3 * test_rod.n_elems).reshape((3, test_rod.n_elems)) - mass = np.linspace(6.0, 4.0, test_rod.n_elems + 1) - inv_moi_full = 1.0 / np.linspace(10.0, 15.0, 9 * test_rod.n_elems).reshape( - (3, 3, test_rod.n_elems) - ) - inv_moi = np.diagonal(inv_moi_full).T - dilatation = np.linspace(0.5, 1.5, test_rod.n_elems) - - test_rod.velocity_collection[:, :] = velocity - test_rod.omega_collection[:, :] = omega - test_rod.mass[:] = mass - test_rod.inv_mass_second_moment_of_inertia = inv_moi_full.copy() - test_rod.dilatation = dilatation.copy() - - test_translational_constant = 2.0 - test_rotational_constant = 3.0 - test_dt = np.float64(1.5) - - damper = AnalyticalLinearDamper( - _system=test_rod, - translational_damping_constant=test_translational_constant, - rotational_damping_constant=test_rotational_constant, - time_step=test_dt, - ) - - test_translational_coeff = np.exp(-test_dt * test_translational_constant / mass) - test_rotational_coeff = np.exp( - -test_dt * test_rotational_constant * inv_moi * dilatation.reshape((1, -1)) - ) - - damper.dampen_rates(test_rod, time=np.float64(0.0)) - - expected_velocity = velocity * test_translational_coeff - expected_omega = omega * test_rotational_coeff - - assert_allclose( - expected_velocity, - test_rod.velocity_collection, - atol=Tolerance.atol(), - ) - assert_allclose(expected_omega, test_rod.omega_collection, atol=Tolerance.atol()) diff --git a/tests/test_dissipation/test_damper_base.py b/tests/test_dissipation/test_damper_base.py deleted file mode 100644 index fe2e25c6..00000000 --- a/tests/test_dissipation/test_damper_base.py +++ /dev/null @@ -1,43 +0,0 @@ -__doc__ = """ Test implementation of the damper base class""" - -import numpy as np -from numpy.testing import assert_allclose - -from elastica.dissipation import DamperBase -from elastica.utils import Tolerance -from tests.test_rod.mock_rod import MockTestRod - - -def test_damper_base(): - test_rod = MockTestRod() - test_rod.velocity_collection = np.ones(3) * 5.0 - test_rod.omega_collection = np.ones(3) * 11.0 - - class TestDamper(DamperBase): - def __init__(self, **kwargs): - super().__init__(**kwargs) - - def dampen_rates(self, rod, time): - rod.velocity_collection *= time - rod.omega_collection *= time - - test_damper = TestDamper(_system=test_rod) - test_damper.dampen_rates(test_rod, np.float64(2.0)) - assert_allclose(test_rod.velocity_collection, 10.0, atol=Tolerance.atol()) - assert_allclose(test_rod.omega_collection, 22.0, atol=Tolerance.atol()) - - -def test_damper_base_properties_access(): - test_rod = MockTestRod() - - class TestDamper(DamperBase): - def __init__(self, **kwargs): - super().__init__(**kwargs) - # Able to access properties in constraint class - assert self._system == test_rod - - def dampen_rates(self, rod, time): - assert self._system == test_rod - - test_damper = TestDamper(_system=test_rod) - test_damper.dampen_rates(test_rod, np.float64(2.0)) diff --git a/tests/test_dissipation/test_laplacian_dissipation_filter.py b/tests/test_dissipation/test_laplacian_dissipation_filter.py deleted file mode 100644 index 76d92c56..00000000 --- a/tests/test_dissipation/test_laplacian_dissipation_filter.py +++ /dev/null @@ -1,183 +0,0 @@ -__doc__ = """ Test implementation of the Laplacian dissipation filter""" - -import pytest -import numpy as np -from numpy.testing import assert_allclose - -from elastica.dissipation import LaplaceDissipationFilter -from elastica.utils import Tolerance -from tests.test_rod.mock_rod import MockTestRod, MockTestRingRod - - -@pytest.mark.parametrize("filter_order", [-1, 0, 3.2]) -def test_laplace_dissipation_filter_init_invalid_filter_order(filter_order): - test_rod = MockTestRod() - with pytest.raises(ValueError) as exc_info: - _ = LaplaceDissipationFilter( - _system=test_rod, - filter_order=filter_order, - ) - assert ( - exc_info.value.args[0] - == "Invalid filter order! Filter order must be a positive integer." - ) - - -@pytest.mark.parametrize("filter_order", [2, 3, 4]) -def test_laplace_dissipation_filter_init(filter_order): - - test_rod = MockTestRod() - filter_damper = LaplaceDissipationFilter( - _system=test_rod, - filter_order=filter_order, - ) - assert filter_damper.filter_order == filter_order - assert_allclose( - filter_damper.velocity_filter_term, np.zeros((3, test_rod.n_elems + 1)) - ) - assert_allclose(filter_damper.omega_filter_term, np.zeros((3, test_rod.n_elems))) - - -@pytest.mark.parametrize("filter_order", [2, 3, 4]) -def test_laplace_dissipation_filter_for_constant_field(filter_order): - test_rod = MockTestRod() - filter_damper = LaplaceDissipationFilter( - _system=test_rod, - filter_order=filter_order, - ) - test_rod.velocity_collection[...] = 2.0 - test_rod.omega_collection[...] = 3.0 - filter_damper.dampen_rates(system=test_rod, time=np.float64(0.0)) - # filter should keep a spatially invariant field unaffected - post_damping_velocity_collection = 2.0 * np.ones_like(test_rod.velocity_collection) - post_damping_omega_collection = 3.0 * np.ones_like(test_rod.omega_collection) - assert_allclose( - post_damping_velocity_collection, - test_rod.velocity_collection, - atol=Tolerance.atol(), - ) - assert_allclose( - post_damping_omega_collection, test_rod.omega_collection, atol=Tolerance.atol() - ) - - -def test_laplace_dissipation_filter_for_flip_flop_field(): - filter_order = 1 - test_rod = MockTestRod() - filter_damper = LaplaceDissipationFilter( - _system=test_rod, - filter_order=filter_order, - ) - test_rod.velocity_collection[...] = 0.0 - test_rod.velocity_collection[..., 1::2] = 2.0 - test_rod.omega_collection[...] = 0.0 - test_rod.omega_collection[..., 1::2] = 3.0 - pre_damping_velocity_collection = test_rod.velocity_collection.copy() - pre_damping_omega_collection = test_rod.omega_collection.copy() - filter_damper.dampen_rates(system=test_rod, time=np.float64(0.0)) - post_damping_velocity_collection = np.zeros_like(test_rod.velocity_collection) - post_damping_omega_collection = np.zeros_like(test_rod.omega_collection) - # filter should remove the flip-flop mode th give the average constant mode - post_damping_velocity_collection[..., 1:-1] = 2.0 / 2 - post_damping_omega_collection[..., 1:-1] = 3.0 / 2 - # end values remain untouched - post_damping_velocity_collection[..., 0 :: test_rod.n_elems] = ( - pre_damping_velocity_collection[..., 0 :: test_rod.n_elems] - ) - post_damping_omega_collection[..., 0 :: test_rod.n_elems - 1] = ( - pre_damping_omega_collection[..., 0 :: test_rod.n_elems - 1] - ) - assert_allclose( - post_damping_velocity_collection, - test_rod.velocity_collection, - atol=Tolerance.atol(), - ) - assert_allclose( - post_damping_omega_collection, test_rod.omega_collection, atol=Tolerance.atol() - ) - - -@pytest.mark.parametrize("filter_order", [-1, 0, 3.2]) -def test_laplace_dissipation_filter_init_invalid_filter_order_for_ring_rod( - filter_order, -): - test_rod = MockTestRingRod() - with pytest.raises(ValueError) as exc_info: - _ = LaplaceDissipationFilter( - _system=test_rod, - filter_order=filter_order, - ) - assert ( - exc_info.value.args[0] - == "Invalid filter order! Filter order must be a positive integer." - ) - - -@pytest.mark.parametrize("filter_order", [2, 3, 4]) -def test_laplace_dissipation_filter_init_for_ring_rod(filter_order): - - test_rod = MockTestRingRod() - filter_damper = LaplaceDissipationFilter( - _system=test_rod, - filter_order=filter_order, - ) - assert filter_damper.filter_order == filter_order - assert_allclose( - filter_damper.velocity_filter_term[:, 1:-1], np.zeros((3, test_rod.n_elems)) - ) - assert_allclose( - filter_damper.omega_filter_term[:, 1:-1], np.zeros((3, test_rod.n_elems)) - ) - - -@pytest.mark.parametrize("filter_order", [2, 3, 4]) -def test_laplace_dissipation_filter_for_constant_field_for_ring_rod(filter_order): - test_rod = MockTestRingRod() - filter_damper = LaplaceDissipationFilter( - _system=test_rod, - filter_order=filter_order, - ) - test_rod.velocity_collection[...] = 2.0 - test_rod.omega_collection[...] = 3.0 - filter_damper.dampen_rates(system=test_rod, time=np.float64(0.0)) - # filter should keep a spatially invariant field unaffected - post_damping_velocity_collection = 2.0 * np.ones_like(test_rod.velocity_collection) - post_damping_omega_collection = 3.0 * np.ones_like(test_rod.omega_collection) - assert_allclose( - post_damping_velocity_collection, - test_rod.velocity_collection, - atol=Tolerance.atol(), - ) - assert_allclose( - post_damping_omega_collection, test_rod.omega_collection, atol=Tolerance.atol() - ) - - -def test_laplace_dissipation_filter_for_flip_flop_field_for_ring_rod(): - filter_order = 1 - test_rod = MockTestRingRod() - filter_damper = LaplaceDissipationFilter( - _system=test_rod, - filter_order=filter_order, - ) - test_rod.velocity_collection[...] = 0.0 - test_rod.velocity_collection[..., 1::2] = 2.0 - test_rod.omega_collection[...] = 0.0 - test_rod.omega_collection[..., 1::2] = 3.0 - pre_damping_velocity_collection = test_rod.velocity_collection.copy() - pre_damping_omega_collection = test_rod.omega_collection.copy() - filter_damper.dampen_rates(system=test_rod, time=np.float64(0.0)) - post_damping_velocity_collection = np.zeros_like(test_rod.velocity_collection) - post_damping_omega_collection = np.zeros_like(test_rod.omega_collection) - # filter should remove the flip-flop mode th give the average constant mode - post_damping_velocity_collection[:] = 2.0 / 2 - post_damping_omega_collection[:] = 3.0 / 2 - - assert_allclose( - post_damping_velocity_collection, - test_rod.velocity_collection, - atol=Tolerance.atol(), - ) - assert_allclose( - post_damping_omega_collection, test_rod.omega_collection, atol=Tolerance.atol() - ) From d3c778056ec7ddfee4cbb6f45f2f5a9bbab1ff44 Mon Sep 17 00:00:00 2001 From: Songyuan Cui Date: Tue, 26 Nov 2024 14:54:56 -0600 Subject: [PATCH 11/12] doc: change docstring of test_disspation.py --- tests/test_dissipation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_dissipation.py b/tests/test_dissipation.py index 9c84bad6..3c4d72d8 100644 --- a/tests/test_dissipation.py +++ b/tests/test_dissipation.py @@ -1,4 +1,4 @@ -__doc__ = """ Test implementation of the damper base class """ +__doc__ = """ Test implementation of the dissipation implementation """ from itertools import combinations From 565f6795d638e3fa6b84bf9a5897b0dddb2c42c9 Mon Sep 17 00:00:00 2001 From: Songyuan Cui Date: Tue, 26 Nov 2024 14:56:56 -0600 Subject: [PATCH 12/12] refactor: required uncessary '.view()' --- elastica/dissipation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/elastica/dissipation.py b/elastica/dissipation.py index b351c0d6..83bbf0b7 100644 --- a/elastica/dissipation.py +++ b/elastica/dissipation.py @@ -259,7 +259,7 @@ def _physical_damping_protocol( rotational_damping_constant: np.float64, time_step: np.float64, ) -> DampenType: - nodal_mass = self._system.mass.view() + nodal_mass = self._system.mass self._translational_damping_coefficient = np.exp( -translational_damping_constant / nodal_mass * time_step )