fix and test transform_vector_cartesian_to_cylinder (#4094)

Description of changes:
- Fixes missing ``sqrt`` in calculation of vector transformation
- Unit test the function

src/utils/include/utils/math/coordinate_transformation.hpp

@@ -77,15 +77,14 @@ inline Vector3d transform_vector_cartesian_to_cylinder(Vector3d const &vec,
   std::tie(theta, rotation_axis) = rotation_params(axis, z_axis);
   auto const rotated_pos = vec_rotate(rotation_axis, theta, pos);
   auto const rotated_vec = vec_rotate(rotation_axis, theta, vec);
+  auto const r = std::sqrt(rotated_pos[0] * rotated_pos[0] +
+                           rotated_pos[1] * rotated_pos[1]);
   // v_r = (x * v_x + y * v_y) / sqrt(x^2 + y^2)
   auto const v_r =
-      (rotated_pos[0] * rotated_vec[0] + rotated_pos[1] * rotated_vec[1]) /
-      std::sqrt(rotated_pos[0] * rotated_pos[0] +
-                rotated_pos[1] * rotated_pos[1]);
-  // v_phi = (x * v_y - y * v_x ) / (x^2 + y^2)
+      (rotated_pos[0] * rotated_vec[0] + rotated_pos[1] * rotated_vec[1]) / r;
+  // v_phi = (x * v_y - y * v_x ) / sqrt(x^2 + y^2)
   auto const v_phi =
-      (rotated_pos[0] * rotated_vec[1] - rotated_pos[1] * rotated_vec[0]) /
-      (rotated_pos[0] * rotated_pos[0] + rotated_pos[1] * rotated_pos[1]);
+      (rotated_pos[0] * rotated_vec[1] - rotated_pos[1] * rotated_vec[0]) / r;
   return Vector3d{v_r, v_phi, rotated_vec[2]};
 }
 

src/utils/tests/coordinate_transformation.cpp

@@ -85,3 +85,21 @@ BOOST_AUTO_TEST_CASE(cylinder_to_cartesian_test) {
     BOOST_CHECK(transformed_z[i] == expected_z[i]);
   }
 }
+
+BOOST_AUTO_TEST_CASE(vector_cart_to_cyl_test) {
+  constexpr auto eps = 1e-13;
+  Vector3d const pos{{1.1, 2.2, 3.3}};
+  auto const axis = (Vector3d{{4.4, 5.5, 6.6}}).normalized();
+  Vector3d const vec{{7.7, 8.8, 9.9}};
+
+  auto const vec_cyl = transform_vector_cartesian_to_cylinder(vec, axis, pos);
+
+  // cylindrical basis vectors at pos
+  auto const e_z = axis;
+  auto const e_r = (pos - (pos * axis) * axis).normalized();
+  auto const e_phi = Utils::vector_product(e_z, e_r);
+
+  BOOST_CHECK_SMALL(vec_cyl[0] - vec * e_r, eps);
+  BOOST_CHECK_SMALL(vec_cyl[1] - vec * e_phi, eps);
+  BOOST_CHECK_SMALL(vec_cyl[2] - vec * e_z, eps);
+}

testsuite/python/observable_cylindrical.py

@@ -92,12 +92,13 @@ def set_particles(self):
                  b * np.sin(i * 2.0 * np.pi / (len(self.params['ids']) + 1)),
                  i * (self.params['max_z'] - self.params['min_z']) /
                  (len(self.params['ids']) + 1) - self.params['center'][2]])
-            v_y = (position[0] * np.sqrt(position[0]**2 + position[1]**2) *
-                   self.v_phi + position[1] * self.v_r) / np.sqrt(
-                       position[0]**2 + position[1]**2)
-            v_x = (self.v_r * np.sqrt(position[0]**2 + position[1]**2)
-                   - position[1] * v_y) / position[0]
-            velocity = np.array([v_x, v_y, self.v_z])
+
+            e_z = np.array([0, 0, 1])
+            e_r = position - (position * e_z) * e_z
+            e_r /= np.linalg.norm(e_r)
+            e_phi = np.cross(e_z, e_r)
+            velocity = e_r * self.v_r + e_phi * self.v_phi + e_z * self.v_z 
+
             velocity = self.swap_axis(velocity, self.params['axis'])
             position = self.swap_axis(position, self.params['axis'])
             position += np.array(self.params['center'])

testsuite/python/observable_cylindricalLB.py

@@ -106,11 +106,14 @@ def set_fluid_velocity(self):
         for i, _ in enumerate(node_positions):
             position = np.array(
                 [node_positions[i], node_positions[i], node_positions[i]])
-            v_y = (position[0] * np.sqrt(position[0]**2 + position[1]**2) * self.v_phi +
-                   position[1] * self.v_r) / np.sqrt(position[0]**2 + position[1]**2)
-            v_x = (self.v_r * np.sqrt(position[0]**2 + position[1]**2) -
-                   position[1] * v_y) / position[0]
-            velocity = np.array([v_x, v_y, self.v_z])
+
+            e_z = np.array([0, 0, 1])
+            e_r = position - (position * e_z) * e_z
+            e_r /= np.linalg.norm(e_r)
+            e_phi = np.cross(e_z, e_r)
+
+            velocity = e_r * self.v_r + e_phi * self.v_phi + e_z * self.v_z
+
             velocity = self.swap_axis(velocity, self.params['axis'])
             position = self.swap_axis(position, self.params['axis'])
             position += np.array(self.params['center'])

testsuite/python/tests_common.py

@@ -164,7 +164,7 @@ def transform_vel_from_cartesian_to_polar_coordinates(pos, vel):
     """
     return np.array([
         (pos[0] * vel[0] + pos[1] * vel[1]) / np.sqrt(pos[0]**2 + pos[1]**2),
-        (pos[0] * vel[1] - pos[1] * vel[0]) / (pos[0]**2 + pos[1]**2), vel[2]])
+        (pos[0] * vel[1] - pos[1] * vel[0]) / np.sqrt(pos[0]**2 + pos[1]**2), vel[2]])
 
 
 def convert_vec_body_to_space(system, part, vec):