diff --git a/inputs/torus.pin b/inputs/torus.pin index 82dd178f..99a970fe 100644 --- a/inputs/torus.pin +++ b/inputs/torus.pin @@ -124,6 +124,7 @@ do_fd_on_grid = false derefine_poles = false r_outer = 250 +# hexp_br = 1000. # Initial break radius for hyperexponential coordinates enable_floors = true diff --git a/scripts/bash/format.sh b/scripts/bash/format.sh index 23943ff8..47cc6f89 100755 --- a/scripts/bash/format.sh +++ b/scripts/bash/format.sh @@ -51,6 +51,7 @@ if ! command -v ${PFM} &> /dev/null; then else PFM=$(command -v ${PFM}) echo "black Python formatter found: ${PFM}" + echo "black version: $(${PFM} --version)" fi echo "Formatting Python files..." @@ -59,6 +60,6 @@ for f in $(git grep --untracked -ail res -- :/*.py); do if [ ${VERBOSE} -ge 1 ]; then echo ${f} fi - ${PFM} ${REPO}/${f} + ${PFM} -q ${REPO}/${f} done echo "...Done" diff --git a/scripts/python/plot_torus.py b/scripts/python/plot_torus.py index 44ba73b7..efc596ce 100644 --- a/scripts/python/plot_torus.py +++ b/scripts/python/plot_torus.py @@ -288,6 +288,9 @@ def plot_frame(ifname, fname, savefig, geomfile=None, rlim=40, coords="cartesian parser.add_argument( "--numax", type=float, default=1.0e2, help="Maximum frequency (Hz)" ) + parser.add_argument( + "--rlim", type=float, default=40.0, help="Maximum radius to plot" + ) parser.add_argument( "--nnu", type=int, default=100, help="Number of frequency support points" ) @@ -302,4 +305,4 @@ def plot_frame(ifname, fname, savefig, geomfile=None, rlim=40, coords="cartesian matplotlib.use("Agg") for i, fname in enumerate(args.files): - plot_frame(i, fname, args.savefig, coords=args.coords) + plot_frame(i, fname, args.savefig, rlim=args.rlim, coords=args.coords) diff --git a/src/geometry/fmks.cpp b/src/geometry/fmks.cpp index 85875759..a8bbe80f 100644 --- a/src/geometry/fmks.cpp +++ b/src/geometry/fmks.cpp @@ -46,6 +46,12 @@ void Initialize(ParameterInput *pin, StateDescriptor *geometry) { Real alpha = pin->GetOrAddReal("coordinates", "poly_alpha", 14); Real x0 = pin->GetReal("parthenon/mesh", "x1min"); Real smooth = pin->GetOrAddReal("coordinates", "smooth", 0.5); + const Real hexp_br_large = 1.e30; + PARTHENON_REQUIRE(hexp_br_large < std::numeric_limits::max(), + "Default hyperexponential break radius must be representable!"); + Real hexp_br = pin->GetOrAddReal("coordinates", "hexp_br", hexp_br_large); + Real hexp_nsq = pin->GetOrAddReal("coordinates", "hexp_nsq", 4.0); + Real hexp_csq = pin->GetOrAddReal("coordinates", "hexp_csq", 4.0); Real a = pin->GetReal("geometry", "a"); Real Rh = 1.0 + sqrt(1.0 - a * a); Real xh = log(Rh); @@ -72,6 +78,9 @@ void Initialize(ParameterInput *pin, StateDescriptor *geometry) { x0 = std::min(x0, std::max(0., x1min - nghost * dx1)); params.Add("x0", x0); params.Add("smooth", smooth); + params.Add("hexp_br", hexp_br); + params.Add("hexp_nsq", hexp_nsq); + params.Add("hexp_csq", hexp_csq); params.Add("Rh", Rh); params.Add("xh", xh); params.Add("r_isco_p", r_isco_p); diff --git a/src/geometry/mckinney_gammie_ryan.cpp b/src/geometry/mckinney_gammie_ryan.cpp index 09d362a4..9769cfec 100644 --- a/src/geometry/mckinney_gammie_ryan.cpp +++ b/src/geometry/mckinney_gammie_ryan.cpp @@ -37,7 +37,11 @@ McKinneyGammieRyan GetTransformation(StateDescriptor *pkg) { Real alpha = pkg->Param("alpha"); Real x0 = pkg->Param("x0"); Real smooth = pkg->Param("smooth"); - return McKinneyGammieRyan(derefine_poles, h, xt, alpha, x0, smooth); + Real hexp_br = pkg->Param("hexp_br"); + Real hexp_nsq = pkg->Param("hexp_nsq"); + Real hexp_csq = pkg->Param("hexp_csq"); + return McKinneyGammieRyan(derefine_poles, h, xt, alpha, x0, smooth, hexp_br, hexp_nsq, + hexp_csq); } } // namespace Geometry diff --git a/src/geometry/mckinney_gammie_ryan.hpp b/src/geometry/mckinney_gammie_ryan.hpp index c6481373..10c695e3 100644 --- a/src/geometry/mckinney_gammie_ryan.hpp +++ b/src/geometry/mckinney_gammie_ryan.hpp @@ -46,14 +46,20 @@ class McKinneyGammieRyan { public: McKinneyGammieRyan() : derefine_poles_(true), h_(0.3), xt_(0.82), alpha_(14), x0_(0), smooth_(0.5), - norm_(GetNorm_(alpha_, xt_)) {} + norm_(GetNorm_(alpha_, xt_)), hexp_br_(1000.), hexp_nsq_(4.0), hexp_csq_(4.0) {} McKinneyGammieRyan(Real x0) // this is the most common use-case : derefine_poles_(true), h_(0.3), xt_(0.82), alpha_(14), x0_(x0), smooth_(0.5), - norm_(GetNorm_(alpha_, xt_)) {} + norm_(GetNorm_(alpha_, xt_)), hexp_br_(1000.), hexp_nsq_(4.0), hexp_csq_(4.0) {} + McKinneyGammieRyan(bool derefine_poles, Real h, Real xt, Real alpha, Real x0, + Real smooth, Real hexp_br, Real hexp_nsq, Real hexp_csq) + : derefine_poles_(derefine_poles), h_(h), xt_(xt), alpha_(alpha), x0_(x0), + smooth_(smooth), norm_(GetNorm_(alpha_, xt_)), hexp_br_(hexp_br), + hexp_nsq_(hexp_nsq), hexp_csq_(hexp_csq_) {} McKinneyGammieRyan(bool derefine_poles, Real h, Real xt, Real alpha, Real x0, Real smooth) : derefine_poles_(derefine_poles), h_(h), xt_(xt), alpha_(alpha), x0_(x0), - smooth_(smooth), norm_(GetNorm_(alpha_, xt_)) {} + smooth_(smooth), norm_(GetNorm_(alpha_, xt_)), hexp_br_(1000.), hexp_nsq_(4.0), + hexp_csq_(4.0) {} KOKKOS_INLINE_FUNCTION void operator()(Real X1, Real X2, Real X3, Real C[NDSPACE], Real Jcov[NDSPACE][NDSPACE], Real Jcon[NDSPACE][NDSPACE]) const { @@ -66,7 +72,11 @@ class McKinneyGammieRyan { C[1] = th; C[2] = X3; - const Real drdX1 = std::exp(X1); + // const Real drdX1 = std::exp(X1); + const Real hexp_dr = X1 - hexp_br_; + const Real drdX1 = + std::exp(X1 + (hexp_dr > 0.) * hexp_csq_ * std::pow(hexp_dr, hexp_nsq_)) * + (1. + (hexp_dr > 0.) * hexp_csq_ * hexp_nsq_ * std::pow(hexp_dr, hexp_nsq_ - 1.)); const Real dthGdX2 = M_PI + M_PI * (1 - h_) * std::cos(2 * M_PI * X2); Real dthdX1, dthdX2; if (derefine_poles_) { @@ -142,7 +152,10 @@ class McKinneyGammieRyan { private: KOKKOS_INLINE_FUNCTION - Real r_(const Real X1) const { return std::exp(X1); } + Real r_(const Real X1) const { + const Real hexp_dr = X1 - hexp_br_; + return std::exp(X1 + (hexp_dr > 0.) * hexp_csq_ * std::pow(hexp_dr, hexp_nsq_)); + } KOKKOS_INLINE_FUNCTION void th_(const Real X1, const Real X2, Real &y, Real &th) const { Real thG = thG_(X2); @@ -169,7 +182,6 @@ class McKinneyGammieRyan { else if (th < M_PI) th = M_PI - SINGSMALL; } - // if (std::fabs(M_PI - th) < robust::EPS()) th = } KOKKOS_INLINE_FUNCTION Real thG_(const Real X2) const { @@ -191,6 +203,9 @@ class McKinneyGammieRyan { Real x0_ = 0; // start point of smooth region Real smooth_ = 0.5; Real norm_; + Real hexp_br_ = 1000.; + Real hexp_nsq_ = 4.; + Real hexp_csq_ = 4.; }; template <> diff --git a/tst/unit/geometry/test_geometry.cpp b/tst/unit/geometry/test_geometry.cpp index 2fec197d..48532dbf 100644 --- a/tst/unit/geometry/test_geometry.cpp +++ b/tst/unit/geometry/test_geometry.cpp @@ -549,9 +549,12 @@ TEST_CASE("McKinneyGammieRyan", "[geometry]") { constexpr Real poly_alpha = 14.0; constexpr Real mks_smooth = 0.5; constexpr Real x0 = -0.4409148982097008; + constexpr Real hexp_br = 1000.; + constexpr Real hexp_nsq = 1.; + constexpr Real hexp_csq = 4.; WHEN("We create a modifier object") { McKinneyGammieRyan transformation(derefine_poles, hslope, poly_xt, poly_alpha, x0, - mks_smooth); + mks_smooth, hexp_br, hexp_nsq, hexp_csq); THEN("The modifier can be called on device") { int n_wrong = 100; Kokkos::parallel_reduce( @@ -590,6 +593,36 @@ TEST_CASE("McKinneyGammieRyan", "[geometry]") { if (GetDifference(Jcon[0][0], 0.32222996) > EPS) update += 1; if (GetDifference(Jcon[1][0], 0.00112507) > EPS) update += 1; if (GetDifference(Jcon[1][1], 0.43532481) > EPS) update += 1; + + // Now for a location where hyper-exponential coordinates are active + constexpr Real Xfar[] = {0., 8.5, X[2], X[3]}; + transformation(Xfar[1], Xfar[2], Xfar[3], C, Jcov, Jcon); + // r + if (GetDifference(C[0], 4.914768840299134354e+03) > EPS) update += 1; + // theta + if (GetDifference(C[1], 1.178102621601838651e+00) > EPS) update += 1; + // phi + if (GetDifference(C[2], 3.2724923474893677) > EPS) update += 1; + // Jcov + if (std::abs(Jcov[0][1] > EPS)) update += 1; + if (std::abs(Jcov[0][2] > EPS)) update += 1; + if (std::abs(Jcov[1][2] > EPS)) update += 1; + if (std::abs(Jcov[2][0] > EPS)) update += 1; + if (std::abs(Jcov[2][1] > EPS)) update += 1; + if (GetDifference(Jcov[2][2], 1) > EPS) update += 1; + if (GetDifference(Jcov[0][0], 4.914768840299134354e+03) > EPS) update += 1; + if (GetDifference(Jcov[1][0], -2.015412862105097876e-04) > EPS) update += 1; + if (GetDifference(Jcov[1][1], 2.933523525877780980e+00) > EPS) update += 1; + // Jcon + if (std::abs(Jcon[0][1] > EPS)) update += 1; + if (std::abs(Jcon[0][2] > EPS)) update += 1; + if (std::abs(Jcon[1][2] > EPS)) update += 1; + if (std::abs(Jcon[2][0] > EPS)) update += 1; + if (std::abs(Jcon[2][1] > EPS)) update += 1; + if (GetDifference(Jcon[2][2], 1) > EPS) update += 1; + if (GetDifference(Jcon[0][0], 2.034683690106441685e-04) > EPS) update += 1; + if (GetDifference(Jcon[1][0], 1.397884708672635997e-08) > EPS) update += 1; + if (GetDifference(Jcon[1][1], 3.408869883532895662e-01) > EPS) update += 1; }, n_wrong); REQUIRE(n_wrong == 0);