Update on coordinate transform between cubed sphere panels (#138)

summary -add vs_velocity_rotation.hpp Defined transformation of velocity from panel 1 into neighboring panels and back. -add test_vs_velocity_rotation.cpp Test cases for the transformation, ensure that all inverses are valid. --------- Co-authored-by: mac/cli <[email protected]>
chengcli · Apr 1, 2024 · a518bb7 · a518bb7
1 parent dccef2f
commit a518bb7
Show file tree

Hide file tree

Showing 5 changed files with 350 additions and 1 deletion.
diff --git a/data/fetch_exogcm_opacity.sh b/data/fetch_exogcm_opacity.sh
@@ -5,7 +5,8 @@ FILE="ck_data_01242024"
 DATA_DIR=$1/
 
 # Dropbox link
-DROPBOX_LINK="https://www.dropbox.com/scl/fi/tt2vez9jadfcq3j0wdw8x/ck_data_01242024.tar.gz?rlkey=dv7tmqpfy2uvs0u5lt2wz3c4k&dl=0"
+#DROPBOX_LINK="https://www.dropbox.com/scl/fi/tt2vez9jadfcq3j0wdw8x/ck_data_01242024.tar.gz?rlkey=dv7tmqpfy2uvs0u5lt2wz3c4k&dl=0"
+DROPBOX_LINK="https://www.dropbox.com/scl/fi/6nohk2uv89pupe8tjj9ib/ck_data_01242024.tar.gz?rlkey=58gjv5re6fusz1s9jk9fgpx95&dl=0"
 
 # Read expected SHA256 from the file
 EXPECTED_SHA256=$(grep "$FILE" ${DATA_DIR}checksums.txt | awk '{print $1}')

diff --git a/src/exo3/cs_velocity_rotation.cpp b/src/exo3/cs_velocity_rotation.cpp
@@ -0,0 +1,196 @@
+#include "cs_velocity_rotation.hpp"
+
+namespace CubedSphereUtility {
+
+//! Transform cubed sphere velocity to local cartesian velocity
+void vel_zab_to_zxy(Real *v1, Real *v2, Real *v3, Real a, Real b) {
+  Real x = tan(a);
+  Real y = tan(b);
+
+  Real vx = *v2;
+  Real vy = *v3;
+  Real vz = *v1;
+
+  Real delta = sqrt(x*x + y*y + 1);
+  Real C = sqrt(1 + x*x);
+  Real D = sqrt(1 + y*y);
+
+  *v1 = (vz - D * x * vx - C * y * vy) / delta;
+  *v2 = 
+      (x * vz + D * vx) / delta;
+  *v3 =
+      (y * vz + C * vy) / delta;
+}
+
+//! Transform local cartesian velocity to cubed sphere velocity
+void vel_zxy_to_zab(Real *v1, Real *v2, Real *v3, Real a, Real b) {
+  Real x = tan(a);
+  Real y = tan(b);
+
+  Real vx = *v2;
+  Real vy = *v3;
+  Real vz = *v1;
+
+  Real delta = sqrt(x*x + y*y + 1);
+  Real C = sqrt(1 + x*x);
+  Real D = sqrt(1 + y*y);
+
+  *v1 = (vz + x * vx + y * vy) / delta;
+  *v2 =
+      (-x * vz / D + vx * (1 + y*y) / D - vy * x * y / D) / delta;
+  *v3 =
+      (-y * vz / C - x * y * vx / C + (1 + x*x) * vy / C) / delta;
+}
+
+//! Transform cubed sphere velocity from panel 1 to panel 2
+//! \param a $x = \tan(\xi)$ coordinates
+//! \param b $y = \tan(\eta)$ coordinat
+void vel_zab_from_p1(Real *vz, Real *vx, Real *vy, Real a, Real b,
+                            int panel) {
+  vel_zab_to_zxy(vz, vx, vy, a, b);
+  Real v1 = *vz;
+  Real v2 = *vx;
+  Real v3 = *vy;
+
+  switch (panel) {
+    case 2:
+      // z->y, x->-x, y->z
+      *vz = v3;
+      *vx = -v2;
+      *vy = v1;
+      vel_zxy_to_zab(vz, vx, vy, a, b);
+      break;
+    case 3:
+      // z->-x, x->z, y->y
+      *vz = v2;
+      *vx = -v1;
+      *vy = v3;
+      vel_zxy_to_zab(vz, vx, vy, a, b);
+      break;
+    case 4:
+      // z->-x, x->-y, y->z
+      *vz = -v2;
+      *vx = -v3;
+      *vy = v1;
+      vel_zxy_to_zab(vz, vx, vy, a, b);
+      break;
+    case 6:
+      // z->-y, x->x, y->z
+      *vz = -v3;
+      *vx = v2;
+      *vy = v1;
+      vel_zxy_to_zab(vz, vx, vy, a, b);
+      break;
+  }
+}
+
+void vel_zab_from_p2(Real *vz, Real *vx, Real *vy, Real a, Real b,
+                            int panel) {
+  vel_zab_to_zxy(vz, vx, vy, a, b);
+  Real v1 = *vz;
+  Real v2 = *vx;
+  Real v3 = *vy;
+  switch (panel) {
+    case 1:
+      // z->y, x->-x, y->z
+      *vz = v3;
+      *vx = -v2;
+      *vy = v1;
+      vel_zxy_to_zab(vz, vx, vy, a, b);
+      break;
+    case 2:
+      break;
+    case 4:
+      break;
+    case 5:
+      break;
+  }
+}
+
+void vel_zab_from_p3(Real *vz, Real *vx, Real *vy, Real a, Real b,
+                            int panel) {
+  vel_zab_to_zxy(vz, vx, vy, a, b);
+  Real v1 = *vz;
+  Real v2 = *vx;
+  Real v3 = *vy;
+  switch (panel) {
+    case 1:
+      // z->x, x->-z, y->y
+      *vz = -v2;
+      *vx = v1;
+      *vy = v3;
+      vel_zxy_to_zab(vz, vx, vy, a, b);
+      break;
+    case 2:
+      break;
+    case 5:
+      break;
+    case 6:
+      break;
+  }
+}
+
+void vel_zab_from_p4(Real *vz, Real *vx, Real *vy, Real a, Real b,
+                            int panel) {
+  vel_zab_to_zxy(vz, vx, vy, a, b);
+  Real v1 = *vz;
+  Real v2 = *vx;
+  Real v3 = *vy;
+  switch (panel) {
+    case 1:
+      // z->y, x->-z, y->-x
+      *vz = v3;
+      *vx = -v1;
+      *vy = -v2;
+      vel_zxy_to_zab(vz, vx, vy, a, b);
+      break;
+    case 2:
+      break;
+    case 5:
+      break;
+    case 6:
+      break;
+  }
+}
+
+void vel_zab_from_p5(Real *vz, Real *vx, Real *vy, Real a, Real b,
+                            int panel) {
+  vel_zab_to_zxy(vz, vx, vy, a, b);
+  Real v1 = *vz;
+  Real v2 = *vx;
+  Real v3 = *vy;
+  switch (panel) {
+    case 2:
+      break;
+    case 3:
+      break;
+    case 4:
+      break;
+    case 6:
+      break;
+  }
+}
+
+void vel_zab_from_p6(Real *vz, Real *vx, Real *vy, Real a, Real b,
+                            int panel) {
+  vel_zab_to_zxy(vz, vx, vy, a, b);
+  Real v1 = *vz;
+  Real v2 = *vx;
+  Real v3 = *vy;
+  switch (panel) {
+    case 1:
+      // z->y, x->x, y->-z
+      *vz = v3;
+      *vx = v2;
+      *vy = -v1;
+      vel_zxy_to_zab(vz, vx, vy, a, b);
+      break;
+    case 3:
+      break;
+    case 4:
+      break;
+    case 5:
+      break;
+  }
+}
+}  // namespace CubedSphereUtility
diff --git a/src/exo3/cs_velocity_rotation.hpp b/src/exo3/cs_velocity_rotation.hpp
@@ -0,0 +1,65 @@
+#ifndef SRC_EXO3_VELOCITY_ROTATION_HPP_
+#define SRC_EXO3_VELOCITY_ROTATION_HPP_
+
+#include <athena/athena.hpp>  // Real
+
+//! Panel number assiangments
+//! (1) Top
+//! (2) Front
+//! (3) Left
+//!
+//!                                           z
+//!                ___________                |
+//!                |\        .\           x___| (1)
+//!                | \   1   . \              /
+//!                |  \_________\           y/
+//!                | 3 |     .  |
+//!                \. .|......  |
+//! z___ (3)        \  |    2 . |
+//!    /|            \ |       .|             y
+//!   / |             \|________|             |
+//!  y  x                                 (2) |___x
+//!                                          /
+//!                                        z/
+
+//! Panel number assiangments
+//! (4) Right
+//! (5) Bottom
+//! (6) Back
+//!                                        y  x
+//!                __________              | /
+//!                |\       .\             |/___z
+//!                | \      . \           (4)
+//!     y  z       |  \________\
+//!     | /        |  |  6  .  |
+//! x___|/         |..|......  |
+//!      (6)       \  |     . 4|       (5) ___ x
+//!                 \ |  5   . |          /|
+//!                  \|_______.|         / |
+//!                                     y  z
+
+namespace CubedSphereUtility {
+
+//! Transform cubed sphere velocity to local cartesian velocity
+void vel_zab_to_zxy(Real *v1, Real *v2, Real *v3, Real a, Real b);
+//! Transform local cartesian velocity to cubed sphere velocity
+void vel_zxy_to_zab(Real *v1, Real *v2, Real *v3, Real a, Real b);
+
+//! Transform cubed sphere velocity from panel 1 to panel 2
+//! \param a $x = \tan(\xi)$ coordinates
+//! \param b $y = \tan(\eta)$ coordinat
+void vel_zab_from_p1(Real *vz, Real *vx, Real *vy, Real a, Real b,
+                            int panel);
+void vel_zab_from_p2(Real *vz, Real *vx, Real *vy, Real a, Real b,
+                            int panel);
+void vel_zab_from_p3(Real *vz, Real *vx, Real *vy, Real a, Real b,
+                            int panel);
+void vel_zab_from_p4(Real *vz, Real *vx, Real *vy, Real a, Real b,
+                            int panel);
+void vel_zab_from_p5(Real *vz, Real *vx, Real *vy, Real a, Real b,
+                            int panel);
+void vel_zab_from_p6(Real *vz, Real *vx, Real *vy, Real a, Real b,
+                            int panel);
+}  // namespace CubedSphereUtility
+
+#endif  // SRC_EXO3_VELOCITY_ROTATION_HPP_
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
@@ -17,6 +17,7 @@ setup_test(test_read_cia_ff)
 setup_test(test_read_rayleigh)
 setup_test(test_read_yaml)
 setup_test(test_helios_ck_read)
+setup_test(test_cs_velocity_rotation)
 
 if (PVFMM_OPTION STREQUAL "ENABLE_PVFMM")
   setup_test(test_nbody)

diff --git a/tests/test_cs_velocity_rotation.cpp b/tests/test_cs_velocity_rotation.cpp
@@ -0,0 +1,86 @@
+// C/C++
+#include <algorithm>
+#include <cmath>
+#include <vector>
+
+// external
+#include <gtest/gtest.h>
+
+#include <application/application.hpp>
+
+// athena
+#include <athena/athena_arrays.hpp>
+
+// exo3
+#include <exo3/cs_velocity_rotation.hpp>
+
+TEST(vel_zab_to_zxy, test_ab_to_xy_to_ab) {
+  Real result[3] = {1,2,3};
+  Real *vz = result;
+  Real *vx = result + 1;
+  Real *vy = result + 2;
+  Real expected_result[3] = {1,2,3};
+  CubedSphereUtility::vel_zxy_to_zab(vz, vx, vy, PI/3, PI/4);
+  std::cout << *vz << " " << *vx << " " << *vy;
+  CubedSphereUtility::vel_zab_to_zxy(vz, vx, vy, PI/3, PI/4);
+  for (int i = 0; i < 3; i++) {
+    EXPECT_NEAR(result[i], expected_result[i], 1e-5);
+  }
+}
+
+TEST(vel_zab_from_test, test_case_1_to_2_to_1) {
+  Real result[3] = {1,2,3};
+  Real *vz = result;
+  Real *vx = result + 1;
+  Real *vy = result + 2;
+  Real expected_result[3] = {1,2,3};
+  CubedSphereUtility::vel_zab_from_p1(vz, vx, vy, PI/5*2, PI/8*3, 2);
+  CubedSphereUtility::vel_zab_from_p2(vz, vx, vy, PI/5*2, PI/8*3, 1);
+  for (int i = 0; i < 3; i++) {
+    EXPECT_NEAR(result[i], expected_result[i], 1e-5);
+  }
+}
+
+TEST(vel_zab_from_test, test_case_1_to_3_to_1) {
+  Real result[3] = {1,2,3};
+  Real *vz = result;
+  Real *vx = result + 1;
+  Real *vy = result + 2;
+  Real expected_result[3] = {1,2,3};
+  CubedSphereUtility::vel_zab_from_p1(vz, vx, vy, PI/5*2, PI/8*3, 3);
+  CubedSphereUtility::vel_zab_from_p3(vz, vx, vy, PI/5*2, PI/8*3, 1);
+  for (int i = 0; i < 3; i++) {
+    EXPECT_NEAR(result[i], expected_result[i], 1e-5);
+  }
+}
+
+TEST(vel_zab_from_test, test_case_1_to_4_to_1) {
+  Real result[3] = {1,2,3};
+  Real *vz = result;
+  Real *vx = result + 1;
+  Real *vy = result + 2;
+  Real expected_result[3] = {1,2,3};
+  CubedSphereUtility::vel_zab_from_p1(vz, vx, vy, PI/5*2, PI/8*3, 4);
+  CubedSphereUtility::vel_zab_from_p4(vz, vx, vy, PI/5*2, PI/8*3, 1);
+  for (int i = 0; i < 3; i++) {
+    EXPECT_NEAR(result[i], expected_result[i], 1e-5);
+  }
+}
+
+TEST(vel_zab_from_test, test_case_1_to_6_to_1) {
+  Real result[3] = {1,2,3};
+  Real *vz = result;
+  Real *vx = result + 1;
+  Real *vy = result + 2;
+  Real expected_result[3] = {1,2,3};
+  CubedSphereUtility::vel_zab_from_p1(vz, vx, vy, PI/5*2, PI/8*3, 6);
+  CubedSphereUtility::vel_zab_from_p6(vz, vx, vy, PI/5*2, PI/8*3, 1);
+  for (int i = 0; i < 3; i++) {
+    EXPECT_NEAR(result[i], expected_result[i], 1e-5);
+  }
+}
+
+int main(int argc, char **argv) {
+  testing::InitGoogleTest(&argc, argv);
+  return RUN_ALL_TESTS();
+}