diff --git a/src/coordinates/uniform_cartesian.hpp b/src/coordinates/uniform_cartesian.hpp index 32f38bd6a9b5..c1ef9542b96a 100644 --- a/src/coordinates/uniform_cartesian.hpp +++ b/src/coordinates/uniform_cartesian.hpp @@ -20,12 +20,13 @@ namespace parthenon { class UniformCartesian : public UniformCoordinates { public: UniformCartesian() = default; - UniformCartesian(const RegionSize &rs, ParameterInput *pin) - : UniformCoordinates(rs, pin) {} + UniformCartesian(const RegionSize &rs, ParameterInput *pin) + : UniformCoordinates(rs, pin) {} UniformCartesian(const UniformCartesian &src, int coarsen) : UniformCoordinates(src, coarsen) {} constexpr static const char *name_ = "UniformCartesian"; + private: }; diff --git a/src/coordinates/uniform_coordinates.hpp b/src/coordinates/uniform_coordinates.hpp index 93dd6a1f5ebe..b617c84d23c5 100644 --- a/src/coordinates/uniform_coordinates.hpp +++ b/src/coordinates/uniform_coordinates.hpp @@ -75,17 +75,24 @@ class UniformCoordinates { template KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int k, const int j, const int i) const { static_assert(dir > 0 && dir < 4); - if constexpr (dir == X1DIR) return static_cast(this)->template Dxc(i); - else if constexpr (dir == X2DIR) return static_cast(this)->template Dxc(j); - else if constexpr (dir == X3DIR) return static_cast(this)->template Dxc(k); + if constexpr (dir == X1DIR) + return static_cast(this)->template Dxc(i); + else if constexpr (dir == X2DIR) + return static_cast(this)->template Dxc(j); + else if constexpr (dir == X3DIR) + return static_cast(this)->template Dxc(k); PARTHENON_FAIL("Unknown dir."); return 0.0; } - KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int dir, const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int dir, const int k, const int j, + const int i) const { assert(dir > 0 && dir < 4); - if (dir == X1DIR) return static_cast(this)->template Dxc(k, j, i); - else if (dir == X2DIR) return static_cast(this)->template Dxc(k, j, i); - else if (dir == X3DIR) return static_cast(this)->template Dxc(k, j, i); + if (dir == X1DIR) + return static_cast(this)->template Dxc(k, j, i); + else if (dir == X2DIR) + return static_cast(this)->template Dxc(k, j, i); + else if (dir == X3DIR) + return static_cast(this)->template Dxc(k, j, i); PARTHENON_FAIL("Unknown dir."); return 0.0; } @@ -115,9 +122,12 @@ class UniformCoordinates { template KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int k, const int j, const int i) const { static_assert(dir > 0 && dir < 4); - if constexpr (dir == X1DIR) return static_cast(this)->template Xc(i); - else if constexpr (dir == X2DIR) return static_cast(this)->template Xc(j); - else if constexpr (dir == X3DIR) return static_cast(this)->template Xc(k); + if constexpr (dir == X1DIR) + return static_cast(this)->template Xc(i); + else if constexpr (dir == X2DIR) + return static_cast(this)->template Xc(j); + else if constexpr (dir == X3DIR) + return static_cast(this)->template Xc(k); return 0; // To appease compiler } @@ -171,27 +181,32 @@ class UniformCoordinates { template KOKKOS_FORCEINLINE_FUNCTION Real X(const int k, const int j, const int i) const { static_assert(dir > 0 && dir < 4); - if constexpr (dir == X1DIR) return X(i); - else if constexpr (dir == X2DIR) return X(j); - else if constexpr (dir == X3DIR) return X(k); + if constexpr (dir == X1DIR) + return X(i); + else if constexpr (dir == X2DIR) + return X(j); + else if constexpr (dir == X3DIR) + return X(k); return 0.0; } template - KOKKOS_FORCEINLINE_FUNCTION - Real Scale(const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real Scale(const int k, const int j, const int i) const { if constexpr (dir > 0 && dir < 4) return 1.0; PARTHENON_FAIL("Unknown dir"); return 0.0; } template - KOKKOS_FORCEINLINE_FUNCTION - Real Scale(const int dir, const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real Scale(const int dir, const int k, const int j, + const int i) const { assert(dir > 0 && dir < 4); - if (dir == X1DIR) return static_cast(this)->template Scale(k, j, i); - else if (dir == X2DIR) return static_cast(this)->template Scale(k, j, i); - else if (dir == X3DIR) return static_cast(this)->template Scale(k, j, i); + if (dir == X1DIR) + return static_cast(this)->template Scale(k, j, i); + else if (dir == X2DIR) + return static_cast(this)->template Scale(k, j, i); + else if (dir == X3DIR) + return static_cast(this)->template Scale(k, j, i); return 0.0; } @@ -199,15 +214,20 @@ class UniformCoordinates { // CellWidth: Width of cells at cell centers //---------------------------------------- template - KOKKOS_FORCEINLINE_FUNCTION Real CellWidth(const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real CellWidth(const int k, const int j, + const int i) const { assert(dir > 0 && dir < 4); return dx_[dir - 1]; } - KOKKOS_FORCEINLINE_FUNCTION Real CellWidth(const int dir, const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real CellWidth(const int dir, const int k, const int j, + const int i) const { assert(dir > 0 && dir < 4); - if (dir == X1DIR) return static_cast(this)->template CellWidth(k, j, i); - else if (dir == X2DIR) return static_cast(this)->template CellWidth(k, j, i); - else if (dir == X3DIR) return static_cast(this)->template CellWidth(k, j, i); + if (dir == X1DIR) + return static_cast(this)->template CellWidth(k, j, i); + else if (dir == X2DIR) + return static_cast(this)->template CellWidth(k, j, i); + else if (dir == X3DIR) + return static_cast(this)->template CellWidth(k, j, i); return 0.0; } @@ -215,11 +235,13 @@ class UniformCoordinates { // EdgeLength: Length of cell edges //---------------------------------------- template - KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int k, const int j, + const int i) const { assert(dir > 0 && dir < 4); return CellWidth(k, j, i); } - KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int dir, const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int dir, const int k, const int j, + const int i) const { return CellWidth(dir, k, j, i); } @@ -231,9 +253,10 @@ class UniformCoordinates { static_assert(dir > 0 && dir < 4); return area_[dir - 1]; } - KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(const int dir, const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(const int dir, const int k, const int j, + const int i) const { assert(dir > 0 && dir < 4); - switch(dir) { + switch (dir) { case X1DIR: return static_cast(this)->template FaceArea(k, j, i); case X2DIR: @@ -249,7 +272,8 @@ class UniformCoordinates { //---------------------------------------- // CellVolume //---------------------------------------- - KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int k, const int j, + const int i) const { return cell_volume_; } @@ -330,6 +354,7 @@ class UniformCoordinates { const std::array &GetXmin() const { return xmin_; } const std::array &GetStartIndex() const { return istart_; } const char *Name() const { return System::name_; } + private: std::array istart_; std::array xmin_, dx_, area_; diff --git a/src/coordinates/uniform_cylindrical.hpp b/src/coordinates/uniform_cylindrical.hpp index 877507058c8f..908095a9b9e4 100644 --- a/src/coordinates/uniform_cylindrical.hpp +++ b/src/coordinates/uniform_cylindrical.hpp @@ -18,16 +18,17 @@ namespace parthenon { class UniformCylindrical : public UniformCoordinates { - using base_t = UniformCoordinates; + using base_t = UniformCoordinates; + public: + using base_t::CellWidth; using base_t::Dxc; - using base_t::Xc; using base_t::Scale; - using base_t::CellWidth; using base_t::Volume; + using base_t::Xc; UniformCylindrical() = default; - UniformCylindrical(const RegionSize &rs, ParameterInput *pin) - : UniformCoordinates(rs, pin) {} + UniformCylindrical(const RegionSize &rs, ParameterInput *pin) + : UniformCoordinates(rs, pin) {} UniformCylindrical(const UniformCylindrical &src, int coarsen) : UniformCoordinates(src, coarsen) {} constexpr static const char *name_ = "UniformCylindrical"; @@ -38,7 +39,7 @@ class UniformCylindrical : public UniformCoordinates { template KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int idx) const { static_assert(dir > 0 && dir < 4); - return Xc(idx) - Xc(idx-1); + return Xc(idx) - Xc(idx - 1); } //---------------------------------------- @@ -58,8 +59,7 @@ class UniformCylindrical : public UniformCoordinates { } template - KOKKOS_FORCEINLINE_FUNCTION - Real Scale(const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real Scale(const int k, const int j, const int i) const { static_assert(dir > 0 && dir < 4); using TE = TopologicalElement; if constexpr (dir == X1DIR || dir == X2DIR) return 1.0; @@ -71,7 +71,8 @@ class UniformCylindrical : public UniformCoordinates { // CellWidth: width of cell through the centroid //---------------------------------------- template - KOKKOS_FORCEINLINE_FUNCTION Real CellWidth(const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real CellWidth(const int k, const int j, + const int i) const { using TE = TopologicalElement; static_assert(dir > 0 && dir < 4); if constexpr (dir == X1DIR || dir == X2DIR) return Dx(); @@ -83,7 +84,8 @@ class UniformCylindrical : public UniformCoordinates { // EdgeLength: Length of cell edges //---------------------------------------- template - KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int k, const int j, + const int i) const { static_assert(dir > 0 && dir < 4); if constexpr (dir == X1DIR || dir == X2DIR) { // radial and z directions are trivial @@ -92,10 +94,13 @@ class UniformCylindrical : public UniformCoordinates { // phi direction return Xf(k, j, i) * Dx(); } - KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int dir, const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int dir, const int k, const int j, + const int i) const { assert(dir > 0 && dir < 4); - if (dir == X1DIR) return EdgeLength(k, j, i); - else if (dir == X2DIR) return EdgeLength(k, j, i); + if (dir == X1DIR) + return EdgeLength(k, j, i); + else if (dir == X2DIR) + return EdgeLength(k, j, i); return EdgeLength(k, j, i); } @@ -109,35 +114,46 @@ class UniformCylindrical : public UniformCoordinates { return Xf(k, j, i) * Dx() * Dx(); } else if constexpr (dir == X2DIR) { Real r0 = Xf(k, j, i); - Real r1 = Xf(k, j, i+1); - return 0.5 * (r1*r1 - r0*r0) * Dx(); + Real r1 = Xf(k, j, i + 1); + return 0.5 * (r1 * r1 - r0 * r0) * Dx(); } Real r0 = Xf(k, j, i); - Real r1 = Xf(k, j, i+1); - return 0.5 * (r1*r1 - r0*r0) * Dx(); + Real r1 = Xf(k, j, i + 1); + return 0.5 * (r1 * r1 - r0 * r0) * Dx(); } //---------------------------------------- // CellVolume //---------------------------------------- - KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int k, const int j, + const int i) const { return FaceArea(k, j, i) * Dx(); } KOKKOS_FORCEINLINE_FUNCTION - Real Volume(CellLevel cl, TopologicalElement el, const int k, const int j, const int i) { + Real Volume(CellLevel cl, TopologicalElement el, const int k, const int j, + const int i) { using TE = TopologicalElement; if (cl == CellLevel::same) { - if (el == TE::CC) return CellVolume(k, j, i); - else if (el == TE::F1) return FaceArea(k, j, i); - else if (el == TE::F2) return FaceArea(k, j, i); - else if (el == TE::F3) return FaceArea(k, j, i); - else if (el == TE::E1) return EdgeLength(k, j, i); - else if (el == TE::E2) return EdgeLength(k, j, i); - else if (el == TE::E3) return EdgeLength(k, j, i); - else if (el == TE::NN) return 1.0; + if (el == TE::CC) + return CellVolume(k, j, i); + else if (el == TE::F1) + return FaceArea(k, j, i); + else if (el == TE::F2) + return FaceArea(k, j, i); + else if (el == TE::F3) + return FaceArea(k, j, i); + else if (el == TE::E1) + return EdgeLength(k, j, i); + else if (el == TE::E2) + return EdgeLength(k, j, i); + else if (el == TE::E3) + return EdgeLength(k, j, i); + else if (el == TE::NN) + return 1.0; } else { - PARTHENON_FAIL("Have not yet implemented fine fields for UniformCylindrical coordinates."); + PARTHENON_FAIL( + "Have not yet implemented fine fields for UniformCylindrical coordinates."); } PARTHENON_FAIL("If you reach this point, someone has added a new value to the the " "TopologicalElement enum."); diff --git a/src/coordinates/uniform_spherical.hpp b/src/coordinates/uniform_spherical.hpp index 6f8cae6dc938..b67a8783f03c 100644 --- a/src/coordinates/uniform_spherical.hpp +++ b/src/coordinates/uniform_spherical.hpp @@ -18,16 +18,17 @@ namespace parthenon { class UniformSpherical : public UniformCoordinates { - using base_t = UniformCoordinates; + using base_t = UniformCoordinates; + public: + using base_t::CellWidth; using base_t::Dxc; - using base_t::Xc; using base_t::Scale; - using base_t::CellWidth; using base_t::Volume; + using base_t::Xc; UniformSpherical() = default; - UniformSpherical(const RegionSize &rs, ParameterInput *pin) - : UniformCoordinates(rs, pin) {} + UniformSpherical(const RegionSize &rs, ParameterInput *pin) + : UniformCoordinates(rs, pin) {} UniformSpherical(const UniformSpherical &src, int coarsen) : UniformCoordinates(src, coarsen) {} constexpr static const char *name_ = "UniformSpherical"; @@ -38,7 +39,7 @@ class UniformSpherical : public UniformCoordinates { template KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int idx) const { static_assert(dir > 0 && dir < 4); - return Xc(idx) - Xc(idx-1); + return Xc(idx) - Xc(idx - 1); } //---------------------------------------- @@ -66,14 +67,15 @@ class UniformSpherical : public UniformCoordinates { } template - KOKKOS_FORCEINLINE_FUNCTION - Real Scale(const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real Scale(const int k, const int j, const int i) const { static_assert(dir > 0 && dir < 4); using TE = TopologicalElement; - if constexpr (dir == X1DIR) return 1.0; + if constexpr (dir == X1DIR) + return 1.0; else { const Real r = X(k, j, i); - if constexpr (dir == X2DIR) return r; + if constexpr (dir == X2DIR) + return r; else { const Real th = X(k, j, i); return r * std::sin(th); @@ -86,12 +88,16 @@ class UniformSpherical : public UniformCoordinates { // CellWidth: width of cell through the centroid //---------------------------------------- template - KOKKOS_FORCEINLINE_FUNCTION Real CellWidth(const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real CellWidth(const int k, const int j, + const int i) const { using TE = TopologicalElement; static_assert(dir > 0 && dir < 4); - if constexpr (dir == X1DIR) return Dx(); - else if constexpr (dir == X2DIR) return Xc(i) * Dx(); - else if constexpr (dir == X3DIR) return Xc(i) * std::sin(Xc(j)) * Dx(); + if constexpr (dir == X1DIR) + return Dx(); + else if constexpr (dir == X2DIR) + return Xc(i) * Dx(); + else if constexpr (dir == X3DIR) + return Xc(i) * std::sin(Xc(j)) * Dx(); return 0.0; } @@ -99,7 +105,8 @@ class UniformSpherical : public UniformCoordinates { // EdgeLength: Length of cell edges //---------------------------------------- template - KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int k, const int j, + const int i) const { static_assert(dir > 0 && dir < 4); if constexpr (dir == X1DIR) { // radial direction is trivial @@ -111,10 +118,13 @@ class UniformSpherical : public UniformCoordinates { // phi direction return Xf(k, j, i) * std::sin(Xf(k, j, i)) * Dx(); } - KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int dir, const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int dir, const int k, const int j, + const int i) const { assert(dir > 0 && dir < 4); - if (dir == X1DIR) return EdgeLength(k, j, i); - else if (dir == X2DIR) return EdgeLength(k, j, i); + if (dir == X1DIR) + return EdgeLength(k, j, i); + else if (dir == X2DIR) + return EdgeLength(k, j, i); return EdgeLength(k, j, i); } @@ -125,41 +135,54 @@ class UniformSpherical : public UniformCoordinates { KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(const int k, const int j, const int i) const { static_assert(dir > 0 && dir < 4); if constexpr (dir == X1DIR) { - Real dOmega = Dx() * (std::cos(Xf(k, j, i)) - std::cos(Xf(k, j+1, i))); + Real dOmega = + Dx() * (std::cos(Xf(k, j, i)) - std::cos(Xf(k, j + 1, i))); Real r = Xf(k, j, i); return r * r * dOmega; } else if constexpr (dir == X2DIR) { Real r0 = Xf(k, j, i); - Real r1 = Xf(k, j, i+1); - return (r1*r1*r1 - r0*r0*r0) * std::sin(Xf(k, j, i)) * Dx() / 3.0; + Real r1 = Xf(k, j, i + 1); + return (r1 * r1 * r1 - r0 * r0 * r0) * std::sin(Xf(k, j, i)) * Dx() / + 3.0; } Real r0 = Xf(k, j, i); - Real r1 = Xf(k, j, i+1); - Real dcth = std::cos(Xf(k, j, i)) - std::cos(Xf(k, j+1, i)); - return (r1*r1*r1 - r0*r0*r0) * dcth / 3.0; + Real r1 = Xf(k, j, i + 1); + Real dcth = std::cos(Xf(k, j, i)) - std::cos(Xf(k, j + 1, i)); + return (r1 * r1 * r1 - r0 * r0 * r0) * dcth / 3.0; } //---------------------------------------- // CellVolume //---------------------------------------- - KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int k, const int j, + const int i) const { return FaceArea(k, j, i) * Dx(); } KOKKOS_FORCEINLINE_FUNCTION - Real Volume(CellLevel cl, TopologicalElement el, const int k, const int j, const int i) { + Real Volume(CellLevel cl, TopologicalElement el, const int k, const int j, + const int i) { using TE = TopologicalElement; if (cl == CellLevel::same) { - if (el == TE::CC) return CellVolume(k, j, i); - else if (el == TE::F1) return FaceArea(k, j, i); - else if (el == TE::F2) return FaceArea(k, j, i); - else if (el == TE::F3) return FaceArea(k, j, i); - else if (el == TE::E1) return EdgeLength(k, j, i); - else if (el == TE::E2) return EdgeLength(k, j, i); - else if (el == TE::E3) return EdgeLength(k, j, i); - else if (el == TE::NN) return 1.0; + if (el == TE::CC) + return CellVolume(k, j, i); + else if (el == TE::F1) + return FaceArea(k, j, i); + else if (el == TE::F2) + return FaceArea(k, j, i); + else if (el == TE::F3) + return FaceArea(k, j, i); + else if (el == TE::E1) + return EdgeLength(k, j, i); + else if (el == TE::E2) + return EdgeLength(k, j, i); + else if (el == TE::E3) + return EdgeLength(k, j, i); + else if (el == TE::NN) + return 1.0; } else { - PARTHENON_FAIL("Have not yet implemented fine fields for UniformSpherical coordinates."); + PARTHENON_FAIL( + "Have not yet implemented fine fields for UniformSpherical coordinates."); } PARTHENON_FAIL("If you reach this point, someone has added a new value to the the " "TopologicalElement enum."); diff --git a/src/utils/index_split.hpp b/src/utils/index_split.hpp index 4767bcff8aee..84b71210766f 100644 --- a/src/utils/index_split.hpp +++ b/src/utils/index_split.hpp @@ -67,13 +67,9 @@ class IndexSplit { } KOKKOS_FORCEINLINE_FUNCTION - int get_i(const int idx) const { - return idx % (ibe_entire_ + 1); - } + int get_i(const int idx) const { return idx % (ibe_entire_ + 1); } KOKKOS_FORCEINLINE_FUNCTION - int get_deltaj(const int idx) const { - return idx / (ibe_entire_ + 1); - } + int get_deltaj(const int idx) const { return idx / (ibe_entire_ + 1); } KOKKOS_INLINE_FUNCTION bool is_i_ghost(const int idx) const { diff --git a/tst/unit/test_coordinates.cpp b/tst/unit/test_coordinates.cpp index e2a98dbd8a25..990019f493b7 100644 --- a/tst/unit/test_coordinates.cpp +++ b/tst/unit/test_coordinates.cpp @@ -25,12 +25,12 @@ using Real = double; using parthenon::ParameterInput; using parthenon::RegionSize; -using parthenon::X1DIR; -using parthenon::X2DIR; -using parthenon::X3DIR; using parthenon::UniformCartesian; using parthenon::UniformCylindrical; using parthenon::UniformSpherical; +using parthenon::X1DIR; +using parthenon::X2DIR; +using parthenon::X3DIR; #include @@ -58,7 +58,7 @@ TEST_CASE("Checking UniformSpherical") { ParameterInput pin; parthenon::Globals::nghost = 2; GIVEN("A coordinate object") { - const Real rout = 3.5; + const Real rout = 3.5; const Real rin = 0.0; std::array xmin{rin, 0.0, 0.0}; std::array xmax{rout, M_PI, 2 * M_PI}; @@ -68,7 +68,7 @@ TEST_CASE("Checking UniformSpherical") { const auto istart = c.GetStartIndex(); Real dr = (xmax[0] - xmin[0]) / nx[0]; const auto cxmin = c.GetXmin(); - REQUIRE(std::abs(cxmin[0] - (xmin[0] - parthenon::Globals::nghost*dr)) < 1.e-14); + REQUIRE(std::abs(cxmin[0] - (xmin[0] - parthenon::Globals::nghost * dr)) < 1.e-14); int i0 = 6 + istart[0]; Real r0 = xmin[0] + (i0 - istart[0]) * dr; REQUIRE(c.Xf(i0) == r0); @@ -89,7 +89,10 @@ TEST_CASE("Checking UniformSpherical") { } } } - REQUIRE(std::abs(volume - (4.0 / 3.0 * M_PI * (std::pow(rout,3) - std::pow(rin,3)))) / volume < 1.e-13); + REQUIRE( + std::abs(volume - (4.0 / 3.0 * M_PI * (std::pow(rout, 3) - std::pow(rin, 3)))) / + volume < + 1.e-13); } parthenon::Globals::nghost = nghost_save; } @@ -99,7 +102,7 @@ TEST_CASE("Checking UniformCylindrical") { ParameterInput pin; parthenon::Globals::nghost = 2; GIVEN("A coordinate object") { - const Real rout = 3.5; + const Real rout = 3.5; const Real rin = 0.0; const Real zlo = 1.3; const Real zhi = 2.78; @@ -111,7 +114,7 @@ TEST_CASE("Checking UniformCylindrical") { const auto istart = c.GetStartIndex(); Real dr = (xmax[0] - xmin[0]) / nx[0]; const auto cxmin = c.GetXmin(); - REQUIRE(std::abs(cxmin[0] - (xmin[0] - parthenon::Globals::nghost*dr)) < 1.e-14); + REQUIRE(std::abs(cxmin[0] - (xmin[0] - parthenon::Globals::nghost * dr)) < 1.e-14); int i0 = 6 + istart[0]; Real r0 = xmin[0] + (i0 - istart[0]) * dr; REQUIRE(c.Xf(i0) == r0); @@ -132,7 +135,8 @@ TEST_CASE("Checking UniformCylindrical") { } } } - REQUIRE(std::abs(volume - M_PI * (rout * rout - rin * rin) * (zhi - zlo)) / volume < 1.e-13); + REQUIRE(std::abs(volume - M_PI * (rout * rout - rin * rin) * (zhi - zlo)) / volume < + 1.e-13); } parthenon::Globals::nghost = nghost_save; } \ No newline at end of file