diff --git a/include/bout/tokamak_coordinates.hxx b/include/bout/tokamak_coordinates.hxx index f13e718787..dd1629799c 100644 --- a/include/bout/tokamak_coordinates.hxx +++ b/include/bout/tokamak_coordinates.hxx @@ -11,18 +11,31 @@ 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; Field2D hthe; FieldMetric I; FieldMetric dx; + std::vector toroidal_angles; - TokamakOptions(Mesh &mesh); + 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..740d19e325 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 { @@ -17,7 +18,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"); @@ -26,6 +27,34 @@ namespace bout { if (mesh.get(dx, "dpsi")) { dx = mesh.getCoordinates()->dx(); } +// mesh.get(toroidal_angle, "z"); + 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() { + + 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; + 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); + k++; + } + } + } + return Coordinates3D(x, y, z); } void TokamakOptions::normalise(BoutReal Lbar, BoutReal Bbar, BoutReal ShearFactor) { @@ -47,7 +76,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/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/fake_mesh_fixture.hxx b/tests/unit/fake_mesh_fixture.hxx index 4c33093aee..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() { + 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 new file mode 100644 index 0000000000..d9ae71e14c --- /dev/null +++ b/tests/unit/mesh/test_change_coordinate_system.cxx @@ -0,0 +1,77 @@ +#include "gtest/gtest.h" +#include "bout/coordinates.hxx" +#include "bout/mesh.hxx" +#include "fake_mesh_fixture.hxx" +#include "bout/constants.hxx" +#include + +static constexpr int NX = 3; +static constexpr int NY = 5; +static constexpr int NZ = 8; + +using bout::globals::mesh; + +class CoordinateTransformTest : public FakeMeshFixture_tmpl { +public: + using FieldMetric = Coordinates::FieldMetric; + + CoordinateTransformTest() : FakeMeshFixture_tmpl() {} +}; + +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 = { // poloidal angle + 0.0, + PI / 2, + PI, + 3 * PI / 2, + 2 * PI}; + + auto tokamak_options = bout::TokamakOptions(*mesh); + + 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++; + } + + // act + bout::Coordinates3D cartesian_coords = tokamak_options.CylindricalCoordinatesToCartesian(); + + // assert + const auto max_r = *std::max_element(begin(r_values), end(r_values)); + const auto expected_max_x = R0 + max_r; + const auto expected_max_y = (R0 + max_r); + const auto expected_max_z = max_r; + + 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"); + 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(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); +}