fix: omega damping rescaled to mass prefactor

GazzolaLab · May 21, 2022 · 89b7bc8 · 89b7bc8
1 parent 4e7321e
commit 89b7bc8
Show file tree

Hide file tree

Showing 4 changed files with 27 additions and 14 deletions.
diff --git a/elastica/rod/cosserat_rod.py b/elastica/rod/cosserat_rod.py
@@ -817,16 +817,16 @@ def _compute_internal_forces(
 
 @numba.njit(cache=True)
 def _compute_damping_torques(
-    damping_torques, J_omega_upon_e, dissipation_constant_for_torques
+    damping_torques, omega_collection, dissipation_constant_for_torques
 ):
     """
-    Update <damping torque> given <angular velocity and mass second moment of inertia>.
+    Update <damping torque> given <angular velocity>.
     """
-    blocksize = damping_torques.shape[1]
+    blocksize = omega_collection.shape[1]
     for i in range(3):
         for k in range(blocksize):
             damping_torques[i, k] = (
-                dissipation_constant_for_torques[k] * J_omega_upon_e[i, k]
+                dissipation_constant_for_torques[k] * omega_collection[i, k]
             )
 
 
@@ -910,7 +910,7 @@ def _compute_internal_torques(
     unsteady_dilatation = J_omega_upon_e * dilatation_rate / dilatation
 
     _compute_damping_torques(
-        damping_torques, J_omega_upon_e, dissipation_constant_for_torques
+        damping_torques, omega_collection, dissipation_constant_for_torques
     )
 
     blocksize = internal_torques.shape[1]

diff --git a/elastica/rod/factory_function.py b/elastica/rod/factory_function.py
@@ -204,12 +204,15 @@ def allocate(
 
     # Custom nu for torques
     dissipation_constant_for_torques = np.zeros((n_elements))
+    elemental_mass = (mass[1:] + mass[:-1]) / 2.0
     if nu_for_torques is None:
         dissipation_constant_for_torques[:] = (dissipation_constant_for_forces / mass)[
             :n_elements
-        ]
+        ] * elemental_mass
     else:
-        dissipation_constant_for_torques[:] = np.asarray(nu_for_torques)
+        dissipation_constant_for_torques[:] = (
+            np.asarray(nu_for_torques) * elemental_mass
+        )
     _assert_dim(
         dissipation_constant_for_torques, 2, "dissipation constant (nu) for torque)"
     )

diff --git a/tests/test_math/test_governing_equations.py b/tests/test_math/test_governing_equations.py
@@ -409,15 +409,17 @@ def test_compute_damping_forces_torques(self, n_elem, nu):
         # Construct velocity and omega
         test_rod.velocity_collection[:] = 1.0
         test_rod.omega_collection[:] = 1.0
-        J_omega_upon_e = test_rod.omega_collection.copy()
         # Compute damping forces and torques
         damping_forces = (
             np.repeat(np.array([1.0, 1.0, 1.0])[:, np.newaxis], n_elem + 1, axis=1)
             * nu
             * test_rod.mass
         )
+        elemental_mass = (test_rod.mass[1:] + test_rod.mass[:-1]) / 2.0
         damping_torques = (
-            np.repeat(np.array([1.0, 1.0, 1.0])[:, np.newaxis], n_elem, axis=1) * nu
+            np.repeat(np.array([1.0, 1.0, 1.0])[:, np.newaxis], n_elem, axis=1)
+            * nu
+            * elemental_mass
         )
         # Compute geometry from state
         _compute_geometry_from_state(
@@ -435,7 +437,7 @@ def test_compute_damping_forces_torques(self, n_elem, nu):
         )
         _compute_damping_torques(
             test_rod.damping_torques,
-            J_omega_upon_e,
+            test_rod.omega_collection,
             test_rod.dissipation_constant_for_torques,
         )
 

diff --git a/tests/test_rod_initialisation.py b/tests/test_rod_initialisation.py
@@ -1191,7 +1191,8 @@ def test_constant_nu_for_torques(n_elems):
         shear_modulus=shear_modulus,
         nu_for_torques=nu_for_torques,
     )
-    correct_nu = nu_for_torques
+    elemental_mass = (mockrod.mass[1:] + mockrod.mass[:-1]) / 2.0
+    correct_nu = nu_for_torques * elemental_mass
     test_nu = mockrod.dissipation_constant_for_torques
     assert_allclose(correct_nu, test_nu, atol=Tolerance.atol())
 
@@ -1236,7 +1237,8 @@ def test_varying_nu_for_torques(n_elems):
         shear_modulus=shear_modulus,
         nu_for_torques=nu_for_torques,
     )
-    correct_nu = nu_for_torques
+    elemental_mass = (mockrod.mass[1:] + mockrod.mass[:-1]) / 2.0
+    correct_nu = nu_for_torques * elemental_mass
     test_nu = mockrod.dissipation_constant_for_torques
     assert_allclose(correct_nu, test_nu, atol=Tolerance.atol())
 
@@ -1319,7 +1321,8 @@ def test_constant_nu_for_torques_if_not_input(n_elems):
         youngs_modulus,
         shear_modulus=shear_modulus,
     )
-    correct_nu = nu
+    elemental_mass = (mockrod.mass[1:] + mockrod.mass[:-1]) / 2.0
+    correct_nu = nu * elemental_mass
     test_nu = mockrod.dissipation_constant_for_torques
     assert_allclose(correct_nu, test_nu, atol=Tolerance.atol())
 
@@ -1593,7 +1596,12 @@ def test_straight_rod(n_elems):
         nu * mockrod.mass,
         atol=Tolerance.atol(),
     )
-    assert_allclose(mockrod.dissipation_constant_for_torques, nu, atol=Tolerance.atol())
+    elemental_mass = (mockrod.mass[1:] + mockrod.mass[:-1]) / 2.0
+    assert_allclose(
+        mockrod.dissipation_constant_for_torques,
+        nu * elemental_mass,
+        atol=Tolerance.atol(),
+    )
     assert_allclose(
         rest_voronoi_lengths, mockrod.rest_voronoi_lengths, atol=Tolerance.atol()
     )