Skip to content

Added Taylor expansion to polar area calculation #110

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
haykh merged 3 commits into from
Jun 3, 2025
Merged

Added Taylor expansion to polar area calculation #110

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
ce195d3
taylor expand polar area to 4th order when small angle
haykh May 26, 2025
4204cb8
remove override to compile the code
May 30, 2025
b12da67
Merge branch '1.2.0rc' into bug/polararea
haykh Jun 3, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 12 additions & 11 deletions src/global/utils/numeric.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand Down
23 changes: 16 additions & 7 deletions src/metrics/qspherical.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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" };
Expand Down Expand Up @@ -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());
}

Expand Down Expand Up @@ -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));
}
Expand All @@ -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);
}
Expand All @@ -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)));
}
}
}

Expand Down
13 changes: 11 additions & 2 deletions src/metrics/spherical.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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" };
Expand All @@ -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());
}

Expand Down Expand Up @@ -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));
}
}

/**
Expand Down