From a518bb79cea6f5d125814ad812a9a89b8303209b Mon Sep 17 00:00:00 2001 From: jeffyujianfu <157439044+jeffyujianfu@users.noreply.github.com> Date: Mon, 1 Apr 2024 19:17:09 -0400 Subject: [PATCH] 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 --- data/fetch_exogcm_opacity.sh | 3 +- src/exo3/cs_velocity_rotation.cpp | 196 ++++++++++++++++++++++++++++ src/exo3/cs_velocity_rotation.hpp | 65 +++++++++ tests/CMakeLists.txt | 1 + tests/test_cs_velocity_rotation.cpp | 86 ++++++++++++ 5 files changed, 350 insertions(+), 1 deletion(-) create mode 100644 src/exo3/cs_velocity_rotation.cpp create mode 100644 src/exo3/cs_velocity_rotation.hpp create mode 100644 tests/test_cs_velocity_rotation.cpp diff --git a/data/fetch_exogcm_opacity.sh b/data/fetch_exogcm_opacity.sh index 5033c06b..7b7047f1 100755 --- 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 new file mode 100644 index 00000000..630fb4bb --- /dev/null +++ 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 \ No newline at end of file diff --git a/src/exo3/cs_velocity_rotation.hpp b/src/exo3/cs_velocity_rotation.hpp new file mode 100644 index 00000000..9df5b4b6 --- /dev/null +++ b/src/exo3/cs_velocity_rotation.hpp @@ -0,0 +1,65 @@ +#ifndef SRC_EXO3_VELOCITY_ROTATION_HPP_ +#define SRC_EXO3_VELOCITY_ROTATION_HPP_ + +#include // 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_ \ No newline at end of file diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 7b8c7a88..5c8eb489 100644 --- 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 new file mode 100644 index 00000000..c64f02fd --- /dev/null +++ b/tests/test_cs_velocity_rotation.cpp @@ -0,0 +1,86 @@ +// C/C++ +#include +#include +#include + +// external +#include + +#include + +// athena +#include + +// exo3 +#include + +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(); +} \ No newline at end of file