From 9fba0611fac9c6596c643dd2c49904ef4c4dfbd4 Mon Sep 17 00:00:00 2001 From: tomc271 Date: Mon, 28 Apr 2025 13:17:52 +0100 Subject: [PATCH 01/21] Add TokamakOptions::CylindricalCoordinatesToCartesian() method for converting to cartesian coordinates from cylindrical coordinates. Add new struct Coordinates3D, to encapsulate (x,y,z) coordinates. --- include/bout/tokamak_coordinates.hxx | 12 ++++++++++++ src/mesh/tokamak_coordinates.cxx | 9 ++++++++- 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/include/bout/tokamak_coordinates.hxx b/include/bout/tokamak_coordinates.hxx index f13e718787..4b89d3a4ce 100644 --- a/include/bout/tokamak_coordinates.hxx +++ b/include/bout/tokamak_coordinates.hxx @@ -11,8 +11,18 @@ using FieldMetric = MetricTensor::FieldMetric; namespace bout { + struct Coordinates3D { + + Field3D x; + Field3D y; + Field3D z; + + Coordinates3D(Field3D x, Field3D y, Field3D z) : x(x), y(y), z(z) {} + }; + struct TokamakOptions { Field2D Rxy; + Field2D Zxy; Field2D Bpxy; Field2D Btxy; Field2D Bxy; @@ -23,6 +33,8 @@ namespace bout { TokamakOptions(Mesh &mesh); void normalise(BoutReal Lbar, BoutReal Bbar, BoutReal ShearFactor); + + Coordinates3D CylindricalCoordinatesToCartesian(); }; BoutReal get_sign_of_bp(const Field2D &Bpxy); diff --git a/src/mesh/tokamak_coordinates.cxx b/src/mesh/tokamak_coordinates.cxx index e5a7636b7d..854ebc7982 100644 --- a/src/mesh/tokamak_coordinates.cxx +++ b/src/mesh/tokamak_coordinates.cxx @@ -17,7 +17,7 @@ namespace bout { TokamakOptions::TokamakOptions(Mesh &mesh) { mesh.get(Rxy, "Rxy"); - // mesh->get(Zxy, "Zxy"); + mesh.get(Zxy, "Zxy"); mesh.get(Bpxy, "Bpxy"); mesh.get(Btxy, "Btxy"); mesh.get(Bxy, "Bxy"); @@ -28,6 +28,13 @@ namespace bout { } } + Coordinates3D TokamakOptions::CylindricalCoordinatesToCartesian() { + Field3D x = Rxy * cos(Zxy); + Field3D y = Rxy * sin(Zxy); + Field3D z = Rxy * sin(Zxy); + return Coordinates3D(x, y, z); + } + void TokamakOptions::normalise(BoutReal Lbar, BoutReal Bbar, BoutReal ShearFactor) { Rxy /= Lbar; Bpxy /= Bbar; From 4c50b3da8460837e46b4433359e692413b8c4756 Mon Sep 17 00:00:00 2001 From: tomc271 Date: Mon, 28 Apr 2025 15:00:12 +0100 Subject: [PATCH 02/21] Add test fixture CoordinateTransformTest, and test CylindricalToCartesian --- tests/unit/CMakeLists.txt | 1 + .../mesh/test_change_coordinate_system.cxx | 47 +++++++++++++++++++ 2 files changed, 48 insertions(+) create mode 100644 tests/unit/mesh/test_change_coordinate_system.cxx diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index 47253c508f..e76a15dd57 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -67,6 +67,7 @@ set(serial_tests_source ./mesh/test_boundary_factory.cxx ./mesh/test_boutmesh.cxx ./mesh/test_coordinates.cxx + ./mesh/test_change_coordinate_system.cxx ./mesh/test_coordinates_accessor.cxx ./mesh/test_interpolation.cxx ./mesh/test_invert3x3.cxx diff --git a/tests/unit/mesh/test_change_coordinate_system.cxx b/tests/unit/mesh/test_change_coordinate_system.cxx new file mode 100644 index 0000000000..f96e545fd7 --- /dev/null +++ b/tests/unit/mesh/test_change_coordinate_system.cxx @@ -0,0 +1,47 @@ +#include "gtest/gtest.h" +#include "bout/coordinates.hxx" +#include "bout/mesh.hxx" +#include "fake_mesh_fixture.hxx" +#include + + +using bout::globals::mesh; + +class CoordinateTransformTest : public FakeMeshFixture { +public: + using FieldMetric = Coordinates::FieldMetric; + + CoordinateTransformTest() : FakeMeshFixture() {} +}; + +TEST_F(CoordinateTransformTest, CylindricalToCartesian) { + + auto tokamak_options = bout::TokamakOptions(*mesh); + + for (int i = 0; i < tokamak_options.Rxy.getNx(); i++) { + for (int j = 0; j < tokamak_options.Rxy.getNy(); j++) { + tokamak_options.Rxy(i, j) = ((float) i + 1) / 1000 * ((float) j + 1) / 1000; + } + } + + bout::Coordinates3D cartesian_coords = tokamak_options.CylindricalCoordinatesToCartesian(); + + for (int jx = 0; jx < mesh->xstart; jx++) { + for (int jy = 0; jy < mesh->ystart; jy++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + + auto actual_x = cartesian_coords.x(jx, jy, jz); + auto actual_y = cartesian_coords.y(jx, jy, jz); + auto actual_z = cartesian_coords.z(jx, jy, jz); + + auto expected_x = tokamak_options.Rxy(jx, jy) * cos(tokamak_options.toroidal_angle(jx, jy, jz)); + auto expected_y = tokamak_options.Rxy(jx, jy) * sin(tokamak_options.toroidal_angle(jx, jy, jz)); + auto expected_z = tokamak_options.Zxy(jx, jy); + + EXPECT_EQ(actual_x, expected_x); + EXPECT_EQ(actual_y, expected_y); + EXPECT_EQ(actual_z, expected_z); + } + } + } +} From a63f7d6c96b1cd9cc8d950d78bcb246bef98f27e Mon Sep 17 00:00:00 2001 From: tomc271 Date: Mon, 28 Apr 2025 17:08:44 +0100 Subject: [PATCH 03/21] Fix CylindricalCoordinatesToCartesian() Add toroidal_angle field to TokamakOptions. Iterate over all coordinate points. --- include/bout/tokamak_coordinates.hxx | 1 + src/mesh/tokamak_coordinates.cxx | 16 +++++++++++++--- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/include/bout/tokamak_coordinates.hxx b/include/bout/tokamak_coordinates.hxx index 4b89d3a4ce..f75825ad3a 100644 --- a/include/bout/tokamak_coordinates.hxx +++ b/include/bout/tokamak_coordinates.hxx @@ -29,6 +29,7 @@ namespace bout { Field2D hthe; FieldMetric I; FieldMetric dx; + FieldMetric toroidal_angle; TokamakOptions(Mesh &mesh); diff --git a/src/mesh/tokamak_coordinates.cxx b/src/mesh/tokamak_coordinates.cxx index 854ebc7982..0ac5db1a2b 100644 --- a/src/mesh/tokamak_coordinates.cxx +++ b/src/mesh/tokamak_coordinates.cxx @@ -26,12 +26,22 @@ namespace bout { if (mesh.get(dx, "dpsi")) { dx = mesh.getCoordinates()->dx(); } + mesh.get(toroidal_angle, "z"); } Coordinates3D TokamakOptions::CylindricalCoordinatesToCartesian() { - Field3D x = Rxy * cos(Zxy); - Field3D y = Rxy * sin(Zxy); - Field3D z = Rxy * sin(Zxy); + Field3D x = emptyFrom(Rxy); + Field3D y = emptyFrom(Rxy); + Field3D z = emptyFrom(Zxy); + for (int i = 0; i < toroidal_angle.getNx(); i++) { + for (int j = 0; j < toroidal_angle.getNy(); j++) { + for (int k = 0; k < toroidal_angle.getNz(); k++) { + x(i, j, k) = Rxy(i, j) * cos(toroidal_angle(i, j, k)); + y(i, j, k) = Rxy(i, j) * sin(toroidal_angle(i, j, k)); + z(i, j, k) = Zxy(i, j); + } + } + } return Coordinates3D(x, y, z); } From 71f9429f2e55f35e57baceaf1c096dbcf272c1ae Mon Sep 17 00:00:00 2001 From: tomc271 Date: Mon, 28 Apr 2025 17:26:52 +0100 Subject: [PATCH 04/21] Fix test CylindricalToCartesian --- tests/unit/mesh/test_change_coordinate_system.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit/mesh/test_change_coordinate_system.cxx b/tests/unit/mesh/test_change_coordinate_system.cxx index f96e545fd7..f6075e5485 100644 --- a/tests/unit/mesh/test_change_coordinate_system.cxx +++ b/tests/unit/mesh/test_change_coordinate_system.cxx @@ -26,8 +26,8 @@ TEST_F(CoordinateTransformTest, CylindricalToCartesian) { bout::Coordinates3D cartesian_coords = tokamak_options.CylindricalCoordinatesToCartesian(); - for (int jx = 0; jx < mesh->xstart; jx++) { - for (int jy = 0; jy < mesh->ystart; jy++) { + for (int jx = 0; jx < mesh->xend; jx++) { + for (int jy = 0; jy < mesh->yend; jy++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { auto actual_x = cartesian_coords.x(jx, jy, jz); From 06df2d62bd0df9a813732622a4471a29407d2725 Mon Sep 17 00:00:00 2001 From: tomc271 Date: Mon, 28 Apr 2025 17:55:43 +0100 Subject: [PATCH 05/21] Calculate toroidal angle, as 2*pi/n --- src/mesh/tokamak_coordinates.cxx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/mesh/tokamak_coordinates.cxx b/src/mesh/tokamak_coordinates.cxx index 0ac5db1a2b..8161acb34a 100644 --- a/src/mesh/tokamak_coordinates.cxx +++ b/src/mesh/tokamak_coordinates.cxx @@ -4,6 +4,7 @@ #include "bout/bout_types.hxx" #include "bout/field2d.hxx" #include "bout/utils.hxx" +#include "bout/constants.hxx" namespace bout { @@ -26,7 +27,8 @@ namespace bout { if (mesh.get(dx, "dpsi")) { dx = mesh.getCoordinates()->dx(); } - mesh.get(toroidal_angle, "z"); +// mesh.get(toroidal_angle, "z"); + toroidal_angle = 2 * PI / Rxy.size(); } Coordinates3D TokamakOptions::CylindricalCoordinatesToCartesian() { From 49ae23b3ebd76e30f972855314764d3b860386e1 Mon Sep 17 00:00:00 2001 From: tomc271 Date: Mon, 28 Apr 2025 17:56:38 +0100 Subject: [PATCH 06/21] Modify generation of test values for Rxy --- tests/unit/mesh/test_change_coordinate_system.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit/mesh/test_change_coordinate_system.cxx b/tests/unit/mesh/test_change_coordinate_system.cxx index f6075e5485..f1498a4982 100644 --- a/tests/unit/mesh/test_change_coordinate_system.cxx +++ b/tests/unit/mesh/test_change_coordinate_system.cxx @@ -20,7 +20,7 @@ TEST_F(CoordinateTransformTest, CylindricalToCartesian) { for (int i = 0; i < tokamak_options.Rxy.getNx(); i++) { for (int j = 0; j < tokamak_options.Rxy.getNy(); j++) { - tokamak_options.Rxy(i, j) = ((float) i + 1) / 1000 * ((float) j + 1) / 1000; + tokamak_options.Rxy(i, j) = sqrt(SQ(i) + SQ(j)); } } From 5387b8a1e2f64cbcfb353706daf3e988b893491e Mon Sep 17 00:00:00 2001 From: tomc271 Date: Wed, 14 May 2025 16:13:51 +0100 Subject: [PATCH 07/21] Use a realistic set of test points as input --- tests/unit/mesh/test_change_coordinate_system.cxx | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/tests/unit/mesh/test_change_coordinate_system.cxx b/tests/unit/mesh/test_change_coordinate_system.cxx index f1498a4982..5edab230cf 100644 --- a/tests/unit/mesh/test_change_coordinate_system.cxx +++ b/tests/unit/mesh/test_change_coordinate_system.cxx @@ -2,6 +2,7 @@ #include "bout/coordinates.hxx" #include "bout/mesh.hxx" #include "fake_mesh_fixture.hxx" +#include "bout/constants.hxx" #include @@ -16,11 +17,16 @@ class CoordinateTransformTest : public FakeMeshFixture { TEST_F(CoordinateTransformTest, CylindricalToCartesian) { + double R0 = 2.0; // major radius + double r[9] = {0.10, 0.15, 0.20, 0.25, 0.30}; // minor radius + double theta[3] = {1.07712, 3.17151, 5.26591}; // poloidal angle + auto tokamak_options = bout::TokamakOptions(*mesh); for (int i = 0; i < tokamak_options.Rxy.getNx(); i++) { - for (int j = 0; j < tokamak_options.Rxy.getNy(); j++) { - tokamak_options.Rxy(i, j) = sqrt(SQ(i) + SQ(j)); + for (int j = 0; j < tokamak_options.Zxy.getNy(); j++) { + tokamak_options.Rxy(i, j) = R0 + r[i] * cos(theta[j]); + tokamak_options.Zxy(i, j) = r[i] * sin(theta[j]); } } From c1ed640c4eeb663be1d9083d645ca8f633b2b680 Mon Sep 17 00:00:00 2001 From: tomc271 Date: Wed, 21 May 2025 14:01:56 +0100 Subject: [PATCH 08/21] Use standard library functions for sin() and cos() --- src/mesh/tokamak_coordinates.cxx | 4 ++-- tests/unit/mesh/test_change_coordinate_system.cxx | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/mesh/tokamak_coordinates.cxx b/src/mesh/tokamak_coordinates.cxx index 8161acb34a..5f22b415e3 100644 --- a/src/mesh/tokamak_coordinates.cxx +++ b/src/mesh/tokamak_coordinates.cxx @@ -38,8 +38,8 @@ namespace bout { for (int i = 0; i < toroidal_angle.getNx(); i++) { for (int j = 0; j < toroidal_angle.getNy(); j++) { for (int k = 0; k < toroidal_angle.getNz(); k++) { - x(i, j, k) = Rxy(i, j) * cos(toroidal_angle(i, j, k)); - y(i, j, k) = Rxy(i, j) * sin(toroidal_angle(i, j, k)); + x(i, j, k) = Rxy(i, j) * std::cos(toroidal_angle(i, j, k)); + y(i, j, k) = Rxy(i, j) * std::sin(toroidal_angle(i, j, k)); z(i, j, k) = Zxy(i, j); } } diff --git a/tests/unit/mesh/test_change_coordinate_system.cxx b/tests/unit/mesh/test_change_coordinate_system.cxx index 5edab230cf..8651c17478 100644 --- a/tests/unit/mesh/test_change_coordinate_system.cxx +++ b/tests/unit/mesh/test_change_coordinate_system.cxx @@ -25,8 +25,8 @@ TEST_F(CoordinateTransformTest, CylindricalToCartesian) { for (int i = 0; i < tokamak_options.Rxy.getNx(); i++) { for (int j = 0; j < tokamak_options.Zxy.getNy(); j++) { - tokamak_options.Rxy(i, j) = R0 + r[i] * cos(theta[j]); - tokamak_options.Zxy(i, j) = r[i] * sin(theta[j]); + tokamak_options.Rxy(i, j) = R0 + r[i] * std::cos(theta[j]); + tokamak_options.Zxy(i, j) = r[i] * std::sin(theta[j]); } } @@ -40,8 +40,8 @@ TEST_F(CoordinateTransformTest, CylindricalToCartesian) { auto actual_y = cartesian_coords.y(jx, jy, jz); auto actual_z = cartesian_coords.z(jx, jy, jz); - auto expected_x = tokamak_options.Rxy(jx, jy) * cos(tokamak_options.toroidal_angle(jx, jy, jz)); - auto expected_y = tokamak_options.Rxy(jx, jy) * sin(tokamak_options.toroidal_angle(jx, jy, jz)); + auto expected_x = tokamak_options.Rxy(jx, jy) * std::cos(tokamak_options.toroidal_angle(jx, jy, jz)); + auto expected_y = tokamak_options.Rxy(jx, jy) * std::sin(tokamak_options.toroidal_angle(jx, jy, jz)); auto expected_z = tokamak_options.Zxy(jx, jy); EXPECT_EQ(actual_x, expected_x); From e6e1d6215401fc69012e01502e5c03f94fd8cd79 Mon Sep 17 00:00:00 2001 From: tomc271 Date: Fri, 16 May 2025 10:03:05 +0100 Subject: [PATCH 09/21] Use a vector for the toroidal angle values. Iterate over all r, theta values --- include/bout/tokamak_coordinates.hxx | 2 +- src/mesh/tokamak_coordinates.cxx | 19 ++++++++++++------- .../mesh/test_change_coordinate_system.cxx | 16 ++++++++-------- 3 files changed, 21 insertions(+), 16 deletions(-) diff --git a/include/bout/tokamak_coordinates.hxx b/include/bout/tokamak_coordinates.hxx index f75825ad3a..7075cf8ccf 100644 --- a/include/bout/tokamak_coordinates.hxx +++ b/include/bout/tokamak_coordinates.hxx @@ -29,7 +29,7 @@ namespace bout { Field2D hthe; FieldMetric I; FieldMetric dx; - FieldMetric toroidal_angle; + std::vector toroidal_angles; TokamakOptions(Mesh &mesh); diff --git a/src/mesh/tokamak_coordinates.cxx b/src/mesh/tokamak_coordinates.cxx index 5f22b415e3..e6272035df 100644 --- a/src/mesh/tokamak_coordinates.cxx +++ b/src/mesh/tokamak_coordinates.cxx @@ -28,18 +28,23 @@ namespace bout { dx = mesh.getCoordinates()->dx(); } // mesh.get(toroidal_angle, "z"); - toroidal_angle = 2 * PI / Rxy.size(); + const auto d_phi = TWOPI / mesh.LocalNz; + auto current_phi = 0.0; + for (int k = 0; k < mesh.LocalNz; k++) { + toroidal_angles.push_back(current_phi); + current_phi += d_phi; + } } Coordinates3D TokamakOptions::CylindricalCoordinatesToCartesian() { Field3D x = emptyFrom(Rxy); Field3D y = emptyFrom(Rxy); Field3D z = emptyFrom(Zxy); - for (int i = 0; i < toroidal_angle.getNx(); i++) { - for (int j = 0; j < toroidal_angle.getNy(); j++) { - for (int k = 0; k < toroidal_angle.getNz(); k++) { - x(i, j, k) = Rxy(i, j) * std::cos(toroidal_angle(i, j, k)); - y(i, j, k) = Rxy(i, j) * std::sin(toroidal_angle(i, j, k)); + for (int i = 0; i < Rxy.getNx(); i++) { + for (int j = 0; j < Rxy.getNy(); j++) { + for (uint k = 0; k < toroidal_angles.size(); k++) { + x(i, j, k) = Rxy(i, j) * cos(toroidal_angles[k]); + y(i, j, k) = Rxy(i, j) * sin(toroidal_angles[k]); z(i, j, k) = Zxy(i, j); } } @@ -66,7 +71,7 @@ namespace bout { const BoutReal sign_of_bp = get_sign_of_bp(tokamak_options.Bpxy); - auto *coord = mesh.getCoordinates(); + auto* coord = mesh.getCoordinates(); const auto g11 = SQ(tokamak_options.Rxy * tokamak_options.Bpxy); const auto g22 = 1.0 / SQ(tokamak_options.hthe); diff --git a/tests/unit/mesh/test_change_coordinate_system.cxx b/tests/unit/mesh/test_change_coordinate_system.cxx index 8651c17478..2fedadb3c2 100644 --- a/tests/unit/mesh/test_change_coordinate_system.cxx +++ b/tests/unit/mesh/test_change_coordinate_system.cxx @@ -18,15 +18,15 @@ class CoordinateTransformTest : public FakeMeshFixture { TEST_F(CoordinateTransformTest, CylindricalToCartesian) { double R0 = 2.0; // major radius - double r[9] = {0.10, 0.15, 0.20, 0.25, 0.30}; // minor radius - double theta[3] = {1.07712, 3.17151, 5.26591}; // poloidal angle + std::array r = {0.10, 0.15, 0.20, 0.25, 0.30}; // minor radius + std::array theta = {0.0, 1.07712, 3.17151, 5.26591}; // poloidal angle auto tokamak_options = bout::TokamakOptions(*mesh); - for (int i = 0; i < tokamak_options.Rxy.getNx(); i++) { - for (int j = 0; j < tokamak_options.Zxy.getNy(); j++) { - tokamak_options.Rxy(i, j) = R0 + r[i] * std::cos(theta[j]); - tokamak_options.Zxy(i, j) = r[i] * std::sin(theta[j]); + for (int i = 0; i < r.size(); i++) { + for (int j = 0; j < theta.size(); j++) { + tokamak_options.Rxy(i, j) = R0 + r[i] * cos(theta[j]); + tokamak_options.Zxy(i, j) = r[i] * sin(theta[j]); } } @@ -40,8 +40,8 @@ TEST_F(CoordinateTransformTest, CylindricalToCartesian) { auto actual_y = cartesian_coords.y(jx, jy, jz); auto actual_z = cartesian_coords.z(jx, jy, jz); - auto expected_x = tokamak_options.Rxy(jx, jy) * std::cos(tokamak_options.toroidal_angle(jx, jy, jz)); - auto expected_y = tokamak_options.Rxy(jx, jy) * std::sin(tokamak_options.toroidal_angle(jx, jy, jz)); + auto expected_x = tokamak_options.Rxy(jx, jy) * std::cos(tokamak_options.toroidal_angles[jz]); + auto expected_y = tokamak_options.Rxy(jx, jy) * std::sin(tokamak_options.toroidal_angles[jz]); auto expected_z = tokamak_options.Zxy(jx, jy); EXPECT_EQ(actual_x, expected_x); From bd86979c3199f16a347c22d08fb7358628759abb Mon Sep 17 00:00:00 2001 From: tomc271 Date: Wed, 21 May 2025 13:53:59 +0100 Subject: [PATCH 10/21] Prefer const --- tests/unit/mesh/test_change_coordinate_system.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/unit/mesh/test_change_coordinate_system.cxx b/tests/unit/mesh/test_change_coordinate_system.cxx index 2fedadb3c2..6b84fcac34 100644 --- a/tests/unit/mesh/test_change_coordinate_system.cxx +++ b/tests/unit/mesh/test_change_coordinate_system.cxx @@ -17,9 +17,9 @@ class CoordinateTransformTest : public FakeMeshFixture { TEST_F(CoordinateTransformTest, CylindricalToCartesian) { - double R0 = 2.0; // major radius - std::array r = {0.10, 0.15, 0.20, 0.25, 0.30}; // minor radius - std::array theta = {0.0, 1.07712, 3.17151, 5.26591}; // poloidal angle + const double R0 = 2.0; // major radius + const std::array r = {0.10, 0.15, 0.20, 0.25, 0.30}; // minor radius + const std::array theta = {0.0, 1.07712, 3.17151, 5.26591}; // poloidal angle auto tokamak_options = bout::TokamakOptions(*mesh); From 392ae19537e3ee1f371c4f49f668892109b14420 Mon Sep 17 00:00:00 2001 From: tomc271 Date: Wed, 21 May 2025 13:57:47 +0100 Subject: [PATCH 11/21] Use range-based loops --- src/mesh/tokamak_coordinates.cxx | 8 +++++--- .../unit/mesh/test_change_coordinate_system.cxx | 16 ++++++++++------ 2 files changed, 15 insertions(+), 9 deletions(-) diff --git a/src/mesh/tokamak_coordinates.cxx b/src/mesh/tokamak_coordinates.cxx index e6272035df..a8ceab74f1 100644 --- a/src/mesh/tokamak_coordinates.cxx +++ b/src/mesh/tokamak_coordinates.cxx @@ -42,10 +42,12 @@ namespace bout { Field3D z = emptyFrom(Zxy); for (int i = 0; i < Rxy.getNx(); i++) { for (int j = 0; j < Rxy.getNy(); j++) { - for (uint k = 0; k < toroidal_angles.size(); k++) { - x(i, j, k) = Rxy(i, j) * cos(toroidal_angles[k]); - y(i, j, k) = Rxy(i, j) * sin(toroidal_angles[k]); + int k = 0; + for (int angle : toroidal_angles) { + x(i, j, k) = Rxy(i, j) * std::cos(angle); + y(i, j, k) = Rxy(i, j) * std::sin(angle); z(i, j, k) = Zxy(i, j); + k++; } } } diff --git a/tests/unit/mesh/test_change_coordinate_system.cxx b/tests/unit/mesh/test_change_coordinate_system.cxx index 6b84fcac34..cf954bd51f 100644 --- a/tests/unit/mesh/test_change_coordinate_system.cxx +++ b/tests/unit/mesh/test_change_coordinate_system.cxx @@ -18,16 +18,20 @@ class CoordinateTransformTest : public FakeMeshFixture { TEST_F(CoordinateTransformTest, CylindricalToCartesian) { const double R0 = 2.0; // major radius - const std::array r = {0.10, 0.15, 0.20, 0.25, 0.30}; // minor radius - const std::array theta = {0.0, 1.07712, 3.17151, 5.26591}; // poloidal angle + const std::array r_values = {0.10, 0.15, 0.20, 0.25, 0.30}; // minor radius + const std::array theta_values = {0.0, 1.07712, 3.17151, 5.26591}; // poloidal angle auto tokamak_options = bout::TokamakOptions(*mesh); - for (int i = 0; i < r.size(); i++) { - for (int j = 0; j < theta.size(); j++) { - tokamak_options.Rxy(i, j) = R0 + r[i] * cos(theta[j]); - tokamak_options.Zxy(i, j) = r[i] * sin(theta[j]); + int i = 0; + for (auto r: r_values) { + int j = 0; + for (auto theta: theta_values) { + tokamak_options.Rxy(i, j) = R0 + r * std::cos(theta); + tokamak_options.Zxy(i, j) = r * std::sin(theta); + j++; } + i++; } bout::Coordinates3D cartesian_coords = tokamak_options.CylindricalCoordinatesToCartesian(); From 25c850e7cd77b2a364e16854cd57b0f8edbac7b4 Mon Sep 17 00:00:00 2001 From: tomc271 Date: Wed, 21 May 2025 13:55:20 +0100 Subject: [PATCH 12/21] Pass parameters by reference --- include/bout/tokamak_coordinates.hxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/bout/tokamak_coordinates.hxx b/include/bout/tokamak_coordinates.hxx index 7075cf8ccf..dd1629799c 100644 --- a/include/bout/tokamak_coordinates.hxx +++ b/include/bout/tokamak_coordinates.hxx @@ -17,7 +17,7 @@ namespace bout { Field3D y; Field3D z; - Coordinates3D(Field3D x, Field3D y, Field3D z) : x(x), y(y), z(z) {} + Coordinates3D(Field3D& x, Field3D& y, Field3D& z) : x(x), y(y), z(z) {} }; struct TokamakOptions { @@ -31,7 +31,7 @@ namespace bout { FieldMetric dx; std::vector toroidal_angles; - TokamakOptions(Mesh &mesh); + TokamakOptions(Mesh& mesh); void normalise(BoutReal Lbar, BoutReal Bbar, BoutReal ShearFactor); From bf20998417128b33698f6d4cedfcc757a132630e Mon Sep 17 00:00:00 2001 From: tomc271 Date: Wed, 21 May 2025 20:39:42 +0100 Subject: [PATCH 13/21] FakeMeshFixture has nx=3, ny=5, nz=7 --- tests/unit/mesh/test_change_coordinate_system.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit/mesh/test_change_coordinate_system.cxx b/tests/unit/mesh/test_change_coordinate_system.cxx index cf954bd51f..f84c6db6f0 100644 --- a/tests/unit/mesh/test_change_coordinate_system.cxx +++ b/tests/unit/mesh/test_change_coordinate_system.cxx @@ -18,8 +18,8 @@ class CoordinateTransformTest : public FakeMeshFixture { TEST_F(CoordinateTransformTest, CylindricalToCartesian) { const double R0 = 2.0; // major radius - const std::array r_values = {0.10, 0.15, 0.20, 0.25, 0.30}; // minor radius - const std::array theta_values = {0.0, 1.07712, 3.17151, 5.26591}; // poloidal angle + const std::array r_values = {0.1, 0.2, 0.3}; // minor radius + const std::array theta_values = {0.0, 1.25663, 2.51327, 3.76991, 5.02654}; // poloidal angle auto tokamak_options = bout::TokamakOptions(*mesh); From 7296f4b72e5693476fd9a06b5e30d5adfe574303 Mon Sep 17 00:00:00 2001 From: tomc271 Date: Wed, 21 May 2025 21:10:15 +0100 Subject: [PATCH 14/21] Use nx and ny variables, to avoid magic numbers --- tests/unit/mesh/test_change_coordinate_system.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit/mesh/test_change_coordinate_system.cxx b/tests/unit/mesh/test_change_coordinate_system.cxx index f84c6db6f0..4215ff04eb 100644 --- a/tests/unit/mesh/test_change_coordinate_system.cxx +++ b/tests/unit/mesh/test_change_coordinate_system.cxx @@ -18,8 +18,8 @@ class CoordinateTransformTest : public FakeMeshFixture { TEST_F(CoordinateTransformTest, CylindricalToCartesian) { const double R0 = 2.0; // major radius - const std::array r_values = {0.1, 0.2, 0.3}; // minor radius - const std::array theta_values = {0.0, 1.25663, 2.51327, 3.76991, 5.02654}; // poloidal angle + const std::array r_values = {0.1, 0.2, 0.3}; // minor radius + const std::array theta_values = {0.0, 1.25663, 2.51327, 3.76991, 5.02654}; // poloidal angle auto tokamak_options = bout::TokamakOptions(*mesh); From 19d64545f22812fa78dbe7461bffe0c39ffa4ce2 Mon Sep 17 00:00:00 2001 From: tomc271 Date: Wed, 21 May 2025 21:36:11 +0100 Subject: [PATCH 15/21] Bug fix: toroidal angle is not an int --- src/mesh/tokamak_coordinates.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh/tokamak_coordinates.cxx b/src/mesh/tokamak_coordinates.cxx index a8ceab74f1..cc81358833 100644 --- a/src/mesh/tokamak_coordinates.cxx +++ b/src/mesh/tokamak_coordinates.cxx @@ -43,7 +43,7 @@ namespace bout { for (int i = 0; i < Rxy.getNx(); i++) { for (int j = 0; j < Rxy.getNy(); j++) { int k = 0; - for (int angle : toroidal_angles) { + for (auto angle : toroidal_angles) { x(i, j, k) = Rxy(i, j) * std::cos(angle); y(i, j, k) = Rxy(i, j) * std::sin(angle); z(i, j, k) = Zxy(i, j); From b9757cd3e271139bb1bd49b12595435b4c1b7020 Mon Sep 17 00:00:00 2001 From: tomc271 Date: Thu, 22 May 2025 10:23:25 +0100 Subject: [PATCH 16/21] Add explanatory comments to test --- tests/unit/mesh/test_change_coordinate_system.cxx | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/tests/unit/mesh/test_change_coordinate_system.cxx b/tests/unit/mesh/test_change_coordinate_system.cxx index 4215ff04eb..342b4a3115 100644 --- a/tests/unit/mesh/test_change_coordinate_system.cxx +++ b/tests/unit/mesh/test_change_coordinate_system.cxx @@ -17,6 +17,12 @@ class CoordinateTransformTest : public FakeMeshFixture { TEST_F(CoordinateTransformTest, CylindricalToCartesian) { + // arrange + + // Set up test values + // Calculate cylindrical coordinates (Rxy, Zxy) + // from (2D) orthogonal poloidal coordinates (r, theta) + const double R0 = 2.0; // major radius const std::array r_values = {0.1, 0.2, 0.3}; // minor radius const std::array theta_values = {0.0, 1.25663, 2.51327, 3.76991, 5.02654}; // poloidal angle @@ -34,8 +40,10 @@ TEST_F(CoordinateTransformTest, CylindricalToCartesian) { i++; } + // act bout::Coordinates3D cartesian_coords = tokamak_options.CylindricalCoordinatesToCartesian(); + // assert for (int jx = 0; jx < mesh->xend; jx++) { for (int jy = 0; jy < mesh->yend; jy++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { From 0fb5b3f5ede5d5b48a96469ca0ea88de16fe767c Mon Sep 17 00:00:00 2001 From: tomc271 Date: Thu, 22 May 2025 13:08:51 +0100 Subject: [PATCH 17/21] Don't instantiate Field3D from Field2D --- src/mesh/tokamak_coordinates.cxx | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/mesh/tokamak_coordinates.cxx b/src/mesh/tokamak_coordinates.cxx index cc81358833..740d19e325 100644 --- a/src/mesh/tokamak_coordinates.cxx +++ b/src/mesh/tokamak_coordinates.cxx @@ -37,9 +37,12 @@ namespace bout { } Coordinates3D TokamakOptions::CylindricalCoordinatesToCartesian() { - Field3D x = emptyFrom(Rxy); - Field3D y = emptyFrom(Rxy); - Field3D z = emptyFrom(Zxy); + + auto* mesh = Rxy.getMesh(); + Field3D x = Field3D(0.0, mesh); + Field3D y = Field3D(0.0, mesh); + Field3D z = Field3D(0.0, mesh); + for (int i = 0; i < Rxy.getNx(); i++) { for (int j = 0; j < Rxy.getNy(); j++) { int k = 0; From 02d4cbfe246ed9002ce44d24007e7edf28be54d1 Mon Sep 17 00:00:00 2001 From: tomc271 Date: Thu, 22 May 2025 13:09:29 +0100 Subject: [PATCH 18/21] Just test min and max values --- .../mesh/test_change_coordinate_system.cxx | 43 ++++++++++++------- 1 file changed, 27 insertions(+), 16 deletions(-) diff --git a/tests/unit/mesh/test_change_coordinate_system.cxx b/tests/unit/mesh/test_change_coordinate_system.cxx index 342b4a3115..3777beb3f1 100644 --- a/tests/unit/mesh/test_change_coordinate_system.cxx +++ b/tests/unit/mesh/test_change_coordinate_system.cxx @@ -25,7 +25,12 @@ TEST_F(CoordinateTransformTest, CylindricalToCartesian) { const double R0 = 2.0; // major radius const std::array r_values = {0.1, 0.2, 0.3}; // minor radius - const std::array theta_values = {0.0, 1.25663, 2.51327, 3.76991, 5.02654}; // poloidal angle + const std::array theta_values = { // poloidal angle + 0.0, + PI / 2, + PI, + 3 * PI / 2, + 2 * PI}; auto tokamak_options = bout::TokamakOptions(*mesh); @@ -44,22 +49,28 @@ TEST_F(CoordinateTransformTest, CylindricalToCartesian) { bout::Coordinates3D cartesian_coords = tokamak_options.CylindricalCoordinatesToCartesian(); // assert - for (int jx = 0; jx < mesh->xend; jx++) { - for (int jy = 0; jy < mesh->yend; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + const auto max_r = *std::max_element(begin(r_values), end(r_values)); + const auto expected_max_x = R0 + max_r; + // With nz=7, there is no toroidal coordinate point at exactly pi/2; the nearest point is at 2/7 * 2pi + const auto expected_max_y = (R0 + max_r) * std::sin(TWOPI * 2 / 7); + const auto expected_max_z = max_r; - auto actual_x = cartesian_coords.x(jx, jy, jz); - auto actual_y = cartesian_coords.y(jx, jy, jz); - auto actual_z = cartesian_coords.z(jx, jy, jz); + // With nz=7, there is no toroidal coordinate point at exactly pi; the nearest point is at 3/7 * 2pi + const auto expected_min_x = -1 * (R0 + max_r) * std::cos(TWOPI / 7 / 2); + const auto expected_min_y = -1 * (R0 + max_r) * std::sin(TWOPI * 2 / 7); + const auto expected_min_z = -1 * expected_max_z; - auto expected_x = tokamak_options.Rxy(jx, jy) * std::cos(tokamak_options.toroidal_angles[jz]); - auto expected_y = tokamak_options.Rxy(jx, jy) * std::sin(tokamak_options.toroidal_angles[jz]); - auto expected_z = tokamak_options.Zxy(jx, jy); + const auto actual_max_x = max(cartesian_coords.x, false, "RGN_ALL"); + const auto actual_max_y = max(cartesian_coords.y, false, "RGN_ALL"); + const auto actual_max_z = max(cartesian_coords.z, false, "RGN_ALL"); + const auto actual_min_x = min(cartesian_coords.x, false, "RGN_ALL"); + const auto actual_min_y = min(cartesian_coords.y, false, "RGN_ALL"); + const auto actual_min_z = min(cartesian_coords.z, false, "RGN_ALL"); - EXPECT_EQ(actual_x, expected_x); - EXPECT_EQ(actual_y, expected_y); - EXPECT_EQ(actual_z, expected_z); - } - } - } + EXPECT_EQ(expected_max_x, actual_max_x); + EXPECT_EQ(expected_max_y, actual_max_y); + EXPECT_EQ(expected_max_z, actual_max_z); + EXPECT_EQ(expected_min_x, actual_min_x); + EXPECT_EQ(expected_min_y, actual_min_y); + EXPECT_EQ(expected_min_z, actual_min_z); } From 086a8b1f48493bbcd11807a36c26e25e90d435ab Mon Sep 17 00:00:00 2001 From: tomc271 Date: Tue, 10 Jun 2025 12:49:11 +0100 Subject: [PATCH 19/21] Add nx, ny, nz as optional arguments to FakeMeshFixture constructor --- tests/unit/fake_mesh_fixture.hxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/unit/fake_mesh_fixture.hxx b/tests/unit/fake_mesh_fixture.hxx index 4c33093aee..dcda38de81 100644 --- a/tests/unit/fake_mesh_fixture.hxx +++ b/tests/unit/fake_mesh_fixture.hxx @@ -27,13 +27,13 @@ /// using MyTest = FakeMeshFixture; class FakeMeshFixture : public ::testing::Test { public: - FakeMeshFixture() { + FakeMeshFixture(int nx_ = nx, int ny_ = ny, int nz_ = nz) { WithQuietOutput quiet_info{output_info}; WithQuietOutput quiet_warn{output_warn}; delete bout::globals::mesh; bout::globals::mpi = new MpiWrapper(); - bout::globals::mesh = new FakeMesh(nx, ny, nz); + bout::globals::mesh = new FakeMesh(nx_, ny_, nz_); bout::globals::mesh->createDefaultRegions(); static_cast(bout::globals::mesh)->setCoordinates(nullptr); test_coords = std::make_shared( @@ -71,7 +71,7 @@ public: dynamic_cast(bout::globals::mesh)->createBoundaryRegions(); delete mesh_staggered; - mesh_staggered = new FakeMesh(nx, ny, nz); + mesh_staggered = new FakeMesh(nx_, ny_, nz_); mesh_staggered->StaggerGrids = true; dynamic_cast(mesh_staggered)->setCoordinates(nullptr); dynamic_cast(mesh_staggered)->setCoordinates(nullptr, CELL_XLOW); From 0d432140fd50f07b68f96f0fe36103590d95c5ed Mon Sep 17 00:00:00 2001 From: tomc271 Date: Tue, 10 Jun 2025 13:00:30 +0100 Subject: [PATCH 20/21] Use nz=8 in CoordinateTransformTest so that there will be a toroidal point at pi/2. --- .../unit/mesh/test_change_coordinate_system.cxx | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/tests/unit/mesh/test_change_coordinate_system.cxx b/tests/unit/mesh/test_change_coordinate_system.cxx index 3777beb3f1..c2a36ed18a 100644 --- a/tests/unit/mesh/test_change_coordinate_system.cxx +++ b/tests/unit/mesh/test_change_coordinate_system.cxx @@ -5,6 +5,9 @@ #include "bout/constants.hxx" #include +static constexpr int NX = 3; +static constexpr int NY = 5; +static constexpr int NZ = 8; using bout::globals::mesh; @@ -12,7 +15,7 @@ class CoordinateTransformTest : public FakeMeshFixture { public: using FieldMetric = Coordinates::FieldMetric; - CoordinateTransformTest() : FakeMeshFixture() {} + CoordinateTransformTest() : FakeMeshFixture(NX, NY, NZ) {} }; TEST_F(CoordinateTransformTest, CylindricalToCartesian) { @@ -24,8 +27,8 @@ TEST_F(CoordinateTransformTest, CylindricalToCartesian) { // from (2D) orthogonal poloidal coordinates (r, theta) const double R0 = 2.0; // major radius - const std::array r_values = {0.1, 0.2, 0.3}; // minor radius - const std::array theta_values = { // poloidal angle + const std::array r_values = {0.1, 0.2, 0.3}; // minor radius + const std::array theta_values = { // poloidal angle 0.0, PI / 2, PI, @@ -51,13 +54,11 @@ TEST_F(CoordinateTransformTest, CylindricalToCartesian) { // assert const auto max_r = *std::max_element(begin(r_values), end(r_values)); const auto expected_max_x = R0 + max_r; - // With nz=7, there is no toroidal coordinate point at exactly pi/2; the nearest point is at 2/7 * 2pi - const auto expected_max_y = (R0 + max_r) * std::sin(TWOPI * 2 / 7); + const auto expected_max_y = (R0 + max_r); const auto expected_max_z = max_r; - // With nz=7, there is no toroidal coordinate point at exactly pi; the nearest point is at 3/7 * 2pi - const auto expected_min_x = -1 * (R0 + max_r) * std::cos(TWOPI / 7 / 2); - const auto expected_min_y = -1 * (R0 + max_r) * std::sin(TWOPI * 2 / 7); + const auto expected_min_x = -1 * (R0 + max_r); + const auto expected_min_y = -1 * (R0 + max_r); const auto expected_min_z = -1 * expected_max_z; const auto actual_max_x = max(cartesian_coords.x, false, "RGN_ALL"); From 64b8bd0d37df382b980981b965c6efa2d03d396a Mon Sep 17 00:00:00 2001 From: tomc271 Date: Thu, 12 Jun 2025 12:38:34 +0100 Subject: [PATCH 21/21] Use type alias FakeMeshFixture = FakeMeshFixture_tmpl<3, 5, 7> as a shim to allow FakeMeshFixture to be used with default values for nx, ny, nz --- tests/unit/fake_mesh_fixture.hxx | 24 ++++++++++++------- .../mesh/test_change_coordinate_system.cxx | 4 ++-- 2 files changed, 17 insertions(+), 11 deletions(-) diff --git a/tests/unit/fake_mesh_fixture.hxx b/tests/unit/fake_mesh_fixture.hxx index dcda38de81..d4c10107e3 100644 --- a/tests/unit/fake_mesh_fixture.hxx +++ b/tests/unit/fake_mesh_fixture.hxx @@ -25,15 +25,19 @@ /// alias to make a new test: /// /// using MyTest = FakeMeshFixture; -class FakeMeshFixture : public ::testing::Test { +/// +/// Type alias FakeMeshFixture = FakeMeshFixture_tmpl<3, 5, 7>; +/// is used as a shim to allow FakeMeshFixture to be used with default values for nx, ny, nz +template +class FakeMeshFixture_tmpl : public ::testing::Test { public: - FakeMeshFixture(int nx_ = nx, int ny_ = ny, int nz_ = nz) { + FakeMeshFixture_tmpl() { WithQuietOutput quiet_info{output_info}; WithQuietOutput quiet_warn{output_warn}; delete bout::globals::mesh; bout::globals::mpi = new MpiWrapper(); - bout::globals::mesh = new FakeMesh(nx_, ny_, nz_); + bout::globals::mesh = new FakeMesh(NX, NY, NZ); bout::globals::mesh->createDefaultRegions(); static_cast(bout::globals::mesh)->setCoordinates(nullptr); test_coords = std::make_shared( @@ -71,7 +75,7 @@ public: dynamic_cast(bout::globals::mesh)->createBoundaryRegions(); delete mesh_staggered; - mesh_staggered = new FakeMesh(nx_, ny_, nz_); + mesh_staggered = new FakeMesh(NX, NY, NZ); mesh_staggered->StaggerGrids = true; dynamic_cast(mesh_staggered)->setCoordinates(nullptr); dynamic_cast(mesh_staggered)->setCoordinates(nullptr, CELL_XLOW); @@ -123,7 +127,7 @@ public: ->setCoordinates(test_coords_staggered, CELL_ZLOW); } - ~FakeMeshFixture() override { + ~FakeMeshFixture_tmpl() override { delete bout::globals::mesh; bout::globals::mesh = nullptr; delete mesh_staggered; @@ -134,12 +138,14 @@ public: Options::cleanup(); } - static constexpr int nx = 3; - static constexpr int ny = 5; - static constexpr int nz = 7; + static constexpr int nx = NX; + static constexpr int ny = NY; + static constexpr int nz = NZ; Mesh* mesh_staggered = nullptr; std::shared_ptr test_coords{nullptr}; std::shared_ptr test_coords_staggered{nullptr}; -}; \ No newline at end of file +}; + +using FakeMeshFixture = FakeMeshFixture_tmpl<3, 5, 7>; \ No newline at end of file diff --git a/tests/unit/mesh/test_change_coordinate_system.cxx b/tests/unit/mesh/test_change_coordinate_system.cxx index c2a36ed18a..d9ae71e14c 100644 --- a/tests/unit/mesh/test_change_coordinate_system.cxx +++ b/tests/unit/mesh/test_change_coordinate_system.cxx @@ -11,11 +11,11 @@ static constexpr int NZ = 8; using bout::globals::mesh; -class CoordinateTransformTest : public FakeMeshFixture { +class CoordinateTransformTest : public FakeMeshFixture_tmpl { public: using FieldMetric = Coordinates::FieldMetric; - CoordinateTransformTest() : FakeMeshFixture(NX, NY, NZ) {} + CoordinateTransformTest() : FakeMeshFixture_tmpl() {} }; TEST_F(CoordinateTransformTest, CylindricalToCartesian) {