diff --git a/src/global/utils/numeric.h b/src/global/utils/numeric.h
index cc1191b6..afd66233 100644
--- a/src/global/utils/numeric.h
+++ b/src/global/utils/numeric.h
@@ -80,17 +80,18 @@ inline constexpr double INV_64 = 0.015625;
 #define CROSS_x3(ax1, ax2, ax3, bx1, bx2, bx3) ((ax1) * (bx2) - (ax2) * (bx1))
 
 namespace constant {
-  inline constexpr std::uint64_t RandomSeed = 0x123456789abcdef0;
-  inline constexpr double        HALF_PI    = 1.57079632679489661923;
-  inline constexpr double        PI         = 3.14159265358979323846;
-  inline constexpr double        INV_PI     = 0.31830988618379067154;
-  inline constexpr double        PI_SQR     = 9.86960440108935861882;
-  inline constexpr double        INV_PI_SQR = 0.10132118364233777144;
-  inline constexpr double        TWO_PI     = 6.28318530717958647692;
-  inline constexpr double        E          = 2.71828182845904523536;
-  inline constexpr double        SQRT2      = 1.41421356237309504880;
-  inline constexpr double        INV_SQRT2  = 0.70710678118654752440;
-  inline constexpr double        SQRT3      = 1.73205080756887729352;
+  inline constexpr std::uint64_t RandomSeed  = 0x123456789abcdef0;
+  inline constexpr double        HALF_PI     = 1.57079632679489661923;
+  inline constexpr double        PI          = 3.14159265358979323846;
+  inline constexpr double        INV_PI      = 0.31830988618379067154;
+  inline constexpr double        PI_SQR      = 9.86960440108935861882;
+  inline constexpr double        INV_PI_SQR  = 0.10132118364233777144;
+  inline constexpr double        TWO_PI      = 6.28318530717958647692;
+  inline constexpr double        E           = 2.71828182845904523536;
+  inline constexpr double        SQRT2       = 1.41421356237309504880;
+  inline constexpr double        INV_SQRT2   = 0.70710678118654752440;
+  inline constexpr double        SQRT3       = 1.73205080756887729352;
+  inline constexpr double        SMALL_ANGLE = 1e-3;
 } // namespace constant
 
 namespace convert {
diff --git a/src/metrics/qspherical.h b/src/metrics/qspherical.h
index 2af1f02a..4a5a9d38 100644
--- a/src/metrics/qspherical.h
+++ b/src/metrics/qspherical.h
@@ -38,6 +38,7 @@ namespace metric {
     const real_t r0, h, chi_min, eta_min, phi_min;
     const real_t dchi, deta, dphi;
     const real_t dchi_inv, deta_inv, dphi_inv;
+    const bool   small_angle;
 
   public:
     static constexpr const char*       Label { "qspherical" };
@@ -69,7 +70,8 @@ namespace metric {
       , dphi { (x3_max - x3_min) / nx3 }
       , dchi_inv { ONE / dchi }
       , deta_inv { ONE / deta }
-      , dphi_inv { ONE / dphi } {
+      , dphi_inv { ONE / dphi }
+      , small_angle { eta2theta(HALF * deta) < constant::SMALL_ANGLE } {
       set_dxMin(find_dxMin());
     }
 
@@ -169,11 +171,11 @@ namespace metric {
      */
     Inline auto sqrt_det_h(const coord_t<D>& x) const -> real_t {
       if constexpr (D == Dim::_2D) {
-        real_t exp_chi { math::exp(x[0] * dchi + chi_min) };
+        const real_t exp_chi { math::exp(x[0] * dchi + chi_min) };
         return dchi * deta * exp_chi * dtheta_deta(x[1] * deta + eta_min) *
                SQR(r0 + exp_chi) * math::sin(eta2theta(x[1] * deta + eta_min));
       } else if constexpr (D == Dim::_3D) {
-        real_t exp_chi { math::exp(x[0] * dchi + chi_min) };
+        const real_t exp_chi { math::exp(x[0] * dchi + chi_min) };
         return dchi * deta * dphi * exp_chi * dtheta_deta(x[1] * deta + eta_min) *
                SQR(r0 + exp_chi) * math::sin(eta2theta(x[1] * deta + eta_min));
       }
@@ -185,7 +187,7 @@ namespace metric {
      */
     Inline auto sqrt_det_h_tilde(const coord_t<D>& x) const -> real_t {
       if constexpr (D != Dim::_1D) {
-        real_t exp_chi { math::exp(x[0] * dchi + chi_min) };
+        const real_t exp_chi { math::exp(x[0] * dchi + chi_min) };
         return dchi * deta * exp_chi * dtheta_deta(x[1] * deta + eta_min) *
                SQR(r0 + exp_chi);
       }
@@ -197,9 +199,16 @@ namespace metric {
      */
     Inline auto polar_area(const real_t& x1) const -> real_t {
       if constexpr (D != Dim::_1D) {
-        real_t exp_chi { math::exp(x1 * dchi + chi_min) };
-        return dchi * exp_chi * SQR(r0 + exp_chi) *
-               (ONE - math::cos(eta2theta(HALF * deta)));
+        const real_t exp_chi { math::exp(x1 * dchi + chi_min) };
+        if (small_angle) {
+          const real_t dtheta = eta2theta(HALF * deta);
+          return dchi * exp_chi * SQR(r0 + exp_chi) *
+                 (static_cast<real_t>(48) - SQR(dtheta)) * SQR(dtheta) /
+                 static_cast<real_t>(384);
+        } else {
+          return dchi * exp_chi * SQR(r0 + exp_chi) *
+                 (ONE - math::cos(eta2theta(HALF * deta)));
+        }
       }
     }
 
diff --git a/src/metrics/spherical.h b/src/metrics/spherical.h
index bf4c3453..05ec6b01 100644
--- a/src/metrics/spherical.h
+++ b/src/metrics/spherical.h
@@ -33,6 +33,7 @@ namespace metric {
 
     const real_t dr, dtheta, dphi;
     const real_t dr_inv, dtheta_inv, dphi_inv;
+    const bool   small_angle;
 
   public:
     static constexpr const char*       Label { "spherical" };
@@ -59,7 +60,8 @@ namespace metric {
       , dphi((x3_max - x3_min) / nx3)
       , dr_inv { ONE / dr }
       , dtheta_inv { ONE / dtheta }
-      , dphi_inv { ONE / dphi } {
+      , dphi_inv { ONE / dphi }
+      , small_angle { HALF * dtheta < constant::SMALL_ANGLE } {
       set_dxMin(find_dxMin());
     }
 
@@ -166,9 +168,16 @@ namespace metric {
     /**
      * differential area at the pole (used in axisymmetric solvers)
      * @param x1 radial coordinate along the axis (code units)
+     * @note uses small-angle approximation when the resolution is too high
      */
     Inline auto polar_area(const real_t& x1) const -> real_t {
-      return dr * SQR(x1 * dr + x1_min) * (ONE - math::cos(HALF * dtheta));
+      if (small_angle) {
+        return dr * SQR(x1 * dr + x1_min) *
+               (static_cast<real_t>(48) - SQR(dtheta)) * SQR(dtheta) /
+               static_cast<real_t>(384);
+      } else {
+        return dr * SQR(x1 * dr + x1_min) * (ONE - math::cos(HALF * dtheta));
+      }
     }
 
     /**