From ce195d3eb79872a181a1ece0f7c135bd00263c71 Mon Sep 17 00:00:00 2001 From: hayk Date: Mon, 26 May 2025 11:11:42 -0400 Subject: [PATCH 1/2] taylor expand polar area to 4th order when small angle --- src/global/utils/numeric.h | 23 +++++++++++----------- src/metrics/qspherical.h | 40 ++++++++++++++++++++++++++++++-------- src/metrics/spherical.h | 30 +++++++++++++++++++++++++--- 3 files changed, 71 insertions(+), 22 deletions(-) diff --git a/src/global/utils/numeric.h b/src/global/utils/numeric.h index cc1191b62..afd662335 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 4df5db866..4a5a9d380 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()); } @@ -97,6 +99,20 @@ namespace metric { return min_dx; } + /** + * total volume of the region described by the metric (in physical units) + */ + [[nodiscard]] + auto totVolume() const -> real_t override { + if constexpr (D == Dim::_1D) { + raise::Error("1D spherical metric not applicable", HERE); + } else if constexpr (D == Dim::_2D) { + return (SQR(x1_max) - SQR(x1_min)) * (x2_max - x2_min); + } else { + return (SQR(x1_max) - SQR(x1_min)) * (x2_max - x2_min) * (x3_max - x3_min); + } + } + /** * metric component with lower indices: h_ij * @param x coordinate array in code units @@ -155,11 +171,11 @@ namespace metric { */ Inline auto sqrt_det_h(const coord_t& 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)); } @@ -171,7 +187,7 @@ namespace metric { */ Inline auto sqrt_det_h_tilde(const coord_t& 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); } @@ -183,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(48) - SQR(dtheta)) * SQR(dtheta) / + static_cast(384); + } else { + return dchi * exp_chi * SQR(r0 + exp_chi) * + (ONE - math::cos(eta2theta(HALF * deta))); + } } } @@ -284,7 +307,8 @@ namespace metric { * @note tetrad/sph <-> cntrv <-> cov */ template - Inline auto transform(const coord_t& xi, const real_t& v_in) const -> real_t { + Inline auto transform(const coord_t& xi, const real_t& v_in) const + -> real_t { static_assert(i > 0 && i <= 3, "Invalid index i"); static_assert(in != out, "Invalid vector transformation"); if constexpr ((in == Idx::T && out == Idx::Sph) || diff --git a/src/metrics/spherical.h b/src/metrics/spherical.h index 54e12367d..05ec6b018 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()); } @@ -76,6 +78,20 @@ namespace metric { return ONE / math::sqrt(ONE / SQR(dx1) + ONE / SQR(dx2)); } + /** + * total volume of the region described by the metric (in physical units) + */ + [[nodiscard]] + auto totVolume() const -> real_t override { + if constexpr (D == Dim::_1D) { + raise::Error("1D spherical metric not applicable", HERE); + } else if constexpr (D == Dim::_2D) { + return (SQR(x1_max) - SQR(x1_min)) * (x2_max - x2_min); + } else { + return (SQR(x1_max) - SQR(x1_min)) * (x2_max - x2_min) * (x3_max - x3_min); + } + } + /** * metric component with lower indices: h_ij * @param x coordinate array in code units @@ -152,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(48) - SQR(dtheta)) * SQR(dtheta) / + static_cast(384); + } else { + return dr * SQR(x1 * dr + x1_min) * (ONE - math::cos(HALF * dtheta)); + } } /** @@ -252,7 +275,8 @@ namespace metric { * @note tetrad/sph <-> cntrv <-> cov */ template - Inline auto transform(const coord_t& xi, const real_t& v_in) const -> real_t { + Inline auto transform(const coord_t& xi, const real_t& v_in) const + -> real_t { static_assert(i > 0 && i <= 3, "Invalid index i"); static_assert(in != out, "Invalid vector transformation"); if constexpr ((in == Idx::T && out == Idx::Sph) || From 4204cb8c7201a0471826880ba81fed9c8cedec5f Mon Sep 17 00:00:00 2001 From: Sasha Chernoglazov Date: Fri, 30 May 2025 12:59:39 -0400 Subject: [PATCH 2/2] remove override to compile the code --- src/metrics/qspherical.h | 2 +- src/metrics/spherical.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/metrics/qspherical.h b/src/metrics/qspherical.h index 4a5a9d380..3b1c8264a 100644 --- a/src/metrics/qspherical.h +++ b/src/metrics/qspherical.h @@ -103,7 +103,7 @@ namespace metric { * total volume of the region described by the metric (in physical units) */ [[nodiscard]] - auto totVolume() const -> real_t override { + auto totVolume() const -> real_t { if constexpr (D == Dim::_1D) { raise::Error("1D spherical metric not applicable", HERE); } else if constexpr (D == Dim::_2D) { diff --git a/src/metrics/spherical.h b/src/metrics/spherical.h index 05ec6b018..ee8018c52 100644 --- a/src/metrics/spherical.h +++ b/src/metrics/spherical.h @@ -82,7 +82,7 @@ namespace metric { * total volume of the region described by the metric (in physical units) */ [[nodiscard]] - auto totVolume() const -> real_t override { + auto totVolume() const -> real_t { if constexpr (D == Dim::_1D) { raise::Error("1D spherical metric not applicable", HERE); } else if constexpr (D == Dim::_2D) {