fix volume

parthenon-hpc-lab · Feb 6, 2025 · 129b36b · 129b36b
1 parent 51ec30e
commit 129b36b
Show file tree

Hide file tree

Showing 3 changed files with 8 additions and 4 deletions.
diff --git a/src/coordinates/uniform_coordinates.hpp b/src/coordinates/uniform_coordinates.hpp
@@ -346,7 +346,7 @@ class UniformCoordinates {
       }
     }
     PARTHENON_FAIL("If you reach this point, someone has added a new value to the the "
-                   "TopologicalElement enum.");
+                   "TopologicalElement or CellLevel enum.");
     return 0.0;
   }
 

diff --git a/src/coordinates/uniform_cylindrical.hpp b/src/coordinates/uniform_cylindrical.hpp
@@ -126,7 +126,9 @@ class UniformCylindrical : public UniformCoordinates<UniformCylindrical> {
   //----------------------------------------
   KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int k, const int j,
                                               const int i) const {
-    return FaceArea<X3DIR>(k, j, i) * Dx<X3DIR>();
+    Real r0 = Xf<X1DIR>(k, j, i);
+    Real r1 = Xf<X1DIR>(k, j, i + 1);
+    return 0.5 * (r1 * r1 - r0 * r0) * Dx<X2DIR>() * Dx<X3DIR>();
   }
 
   KOKKOS_FORCEINLINE_FUNCTION

diff --git a/src/coordinates/uniform_spherical.hpp b/src/coordinates/uniform_spherical.hpp
@@ -147,7 +147,6 @@ class UniformSpherical : public UniformCoordinates<UniformSpherical> {
     }
     Real r0 = Xf<X1DIR>(k, j, i);
     Real r1 = Xf<X1DIR>(k, j, i + 1);
-    Real dcth = std::cos(Xf<X2DIR>(k, j, i)) - std::cos(Xf<X2DIR>(k, j + 1, i));
     return 0.5 * (r1 * r1 - r0 * r0) * Dx<X2DIR>();
   }
 
@@ -156,7 +155,10 @@ class UniformSpherical : public UniformCoordinates<UniformSpherical> {
   //----------------------------------------
   KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int k, const int j,
                                               const int i) const {
-    return FaceArea<X3DIR>(k, j, i) * Dx<X3DIR>();
+    Real dcth = std::cos(Xf<X2DIR>(k, j, i)) - std::cos(Xf<X2DIR>(k, j + 1, i));
+    Real r0 = Xf<X1DIR>(k, j, i);
+    Real r1 = Xf<X1DIR>(k, j, i + 1);
+    return (1.0 / 3.0) * (r1 * r1 * r1 - r0 * r0 * r0) * dcth * Dx<X3DIR>();
   }
 
   KOKKOS_FORCEINLINE_FUNCTION