Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Filter Damper class #123

Merged
merged 12 commits into from
Jul 5, 2022
7 changes: 6 additions & 1 deletion docs/api/damping.rst
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ Damping is used to numerically stabilize the simulations.

DamperBase
ExponentialDamper
LaplaceDissipationFilter

Compatibility
~~~~~~~~~~~~~
Expand All @@ -23,7 +24,8 @@ Compatibility
Damping/Numerical Dissipation Rod Rigid Body
=============================== ==== ===========
ExponentialDamper ✅ ✅
=============================== ==== ===========
LaplaceDissipationFilter ✅ ❌
=============================== ==== ===========


Built-in Constraints
Expand All @@ -37,3 +39,6 @@ Built-in Constraints

.. autoclass:: ExponentialDamper
:special-members: __init__

.. autoclass:: LaplaceDissipationFilter
:special-members: __init__
132 changes: 130 additions & 2 deletions elastica/dissipation.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,10 +6,13 @@
__all__ = [
"DamperBase",
"ExponentialDamper",
"LaplaceDissipationFilter",
]
from abc import ABC, abstractmethod

from elastica.typing import SystemType
from elastica.typing import RodType, SystemType

from numba import njit

import numpy as np

Expand Down Expand Up @@ -85,7 +88,7 @@ class ExponentialDamper(DamperBase):
--------
How to set exponential damper for rod or rigid body:

>>> simulator.dampin(rod).using(
>>> simulator.dampen(rod).using(
... ExponentialDamper,
... damping_constant=0.1,
... time_step = 1E-4, # Simulation time-step
Expand Down Expand Up @@ -151,3 +154,128 @@ def dampen_rates(self, rod: SystemType, time: float):
rod.omega_collection[:] = rod.omega_collection * np.power(
self.rotational_exponential_damping_coefficient, rod.dilatation
)


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] 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
bhosale2 marked this conversation as resolved.
Show resolved Hide resolved
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):
"""
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
self.velocity_filter_term = np.zeros_like(self._system.velocity_collection)
self.omega_filter_term = np.zeros_like(self._system.omega_collection)

def dampen_rates(self, rod: RodType, time: float) -> None:
nb_filter_rate(
rate_collection=rod.velocity_collection,
filter_term=self.velocity_filter_term,
filter_order=self.filter_order,
)
nb_filter_rate(
rate_collection=rod.omega_collection,
filter_term=self.omega_filter_term,
filter_order=self.filter_order,
)


@njit(cache=True)
def nb_filter_rate(
rate_collection: np.ndarray, filter_term: np.ndarray, 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
-----
bhosale2 marked this conversation as resolved.
Show resolved Hide resolved
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
3 changes: 2 additions & 1 deletion elastica/typing.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,4 +3,5 @@

from typing import Type, Union

SystemType = Union[Type[RodBase], Type[RigidBodyBase]]
RodType = Type[RodBase]
SystemType = Union[RodType, Type[RigidBodyBase]]
98 changes: 95 additions & 3 deletions tests/test_dissipation.py
Original file line number Diff line number Diff line change
@@ -1,11 +1,15 @@
__doc__ = """ Test Dissipation module for in Elastica implementation"""

# System imports
from elastica.dissipation import DamperBase, ExponentialDamper, LaplaceDissipationFilter
from elastica.utils import Tolerance

import numpy as np
from test_rod.test_rods import MockTestRod
from elastica.dissipation import DamperBase, ExponentialDamper
from numpy.testing import assert_allclose
from elastica.utils import Tolerance

import pytest

from test_rod.test_rods import MockTestRod


def test_damper_base():
Expand Down Expand Up @@ -114,3 +118,91 @@ def test_exponential_damper():
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()
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_elem + 1))
)
assert_allclose(filter_damper.omega_filter_term, np.zeros((3, test_rod.n_elem)))


@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(rod=test_rod, time=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(rod=test_rod, time=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_elem
] = pre_damping_velocity_collection[..., 0 :: test_rod.n_elem]
post_damping_omega_collection[
..., 0 :: test_rod.n_elem - 1
] = pre_damping_omega_collection[..., 0 :: test_rod.n_elem - 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()
)