From de714ff3531a0d0f507a0811e53c326ecb76e584 Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Wed, 19 Apr 2023 15:31:18 -0600 Subject: [PATCH 1/7] First attempt --- src/geometry/mckinney_gammie_ryan.hpp | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/src/geometry/mckinney_gammie_ryan.hpp b/src/geometry/mckinney_gammie_ryan.hpp index c64813739..79e8344e0 100644 --- a/src/geometry/mckinney_gammie_ryan.hpp +++ b/src/geometry/mckinney_gammie_ryan.hpp @@ -46,14 +46,15 @@ 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_(1.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_(1.0), hexp_csq_(4.0) {} McKinneyGammieRyan(bool derefine_poles, Real h, Real xt, Real alpha, Real x0, - Real smooth) + 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_)) {} + smooth_(smooth), norm_(GetNorm_(alpha_, xt_), hexp_br_(hexp_br), + hexp_nsq_(hexp_nsq), hexp_csq_(hexp_csq_) {} 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 +67,9 @@ 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_); const Real dthGdX2 = M_PI + M_PI * (1 - h_) * std::cos(2 * M_PI * X2); Real dthdX1, dthdX2; if (derefine_poles_) { @@ -142,7 +145,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); @@ -191,6 +197,9 @@ class McKinneyGammieRyan { Real x0_ = 0; // start point of smooth region Real smooth_ = 0.5; Real norm_; + Real hexp_br_ = 1000.; + Real hexp_nsq_ = 1.; + Real hexp_cs1_ = 4.; }; template <> From 99f77b1a817b984bb2ebd1cf44816127fa9566ad Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Wed, 19 Apr 2023 16:00:52 -0600 Subject: [PATCH 2/7] Formatting etc. --- src/geometry/fmks.cpp | 6 ++++++ src/geometry/mckinney_gammie_ryan.cpp | 6 +++++- src/geometry/mckinney_gammie_ryan.hpp | 17 +++++++++++------ tst/unit/geometry/test_geometry.cpp | 5 ++++- 4 files changed, 26 insertions(+), 8 deletions(-) diff --git a/src/geometry/fmks.cpp b/src/geometry/fmks.cpp index 858757596..4921b0e8a 100644 --- a/src/geometry/fmks.cpp +++ b/src/geometry/fmks.cpp @@ -46,6 +46,9 @@ 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); + Real hexp_br = pin->GetOrAddReal("coordinates", "hexp_br", 1000.); + Real hexp_nsq = pin->GetOrAddReal("coordinates", "hexp_nsq", 1.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 +75,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 09d362a45..9769cfec1 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 79e8344e0..3131ad07e 100644 --- a/src/geometry/mckinney_gammie_ryan.hpp +++ b/src/geometry/mckinney_gammie_ryan.hpp @@ -46,15 +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_), hexp_br_(1000.), hexp_nsq_(1.0), hexp_csq_(4.0) {} + norm_(GetNorm_(alpha_, xt_)), hexp_br_(1000.), hexp_nsq_(1.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_)), hexp_br_(1000.), hexp_nsq_(1.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), + 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_)), hexp_br_(1000.), hexp_nsq_(1.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 { @@ -67,9 +72,10 @@ 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_); + const Real drdX1 = + std::exp(X1 + (hexp_dr > 0.) * hexp_csq_ * std::pow(hexp_dr, hexp_nsq_)); const Real dthGdX2 = M_PI + M_PI * (1 - h_) * std::cos(2 * M_PI * X2); Real dthdX1, dthdX2; if (derefine_poles_) { @@ -175,7 +181,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 { @@ -199,7 +204,7 @@ class McKinneyGammieRyan { Real norm_; Real hexp_br_ = 1000.; Real hexp_nsq_ = 1.; - Real hexp_cs1_ = 4.; + Real hexp_csq_ = 4.; }; template <> diff --git a/tst/unit/geometry/test_geometry.cpp b/tst/unit/geometry/test_geometry.cpp index 2fec197de..520647033 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( From 487d557ccecbeb3108f1fe2311b007e654666413 Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Thu, 20 Apr 2023 12:25:41 -0600 Subject: [PATCH 3/7] Fix bug in dxdx --- scripts/python/phoedf.py | 2 +- scripts/python/plot_torus.py | 7 ++++--- src/geometry/fmks.cpp | 2 +- src/geometry/mckinney_gammie_ryan.hpp | 9 +++++---- 4 files changed, 11 insertions(+), 9 deletions(-) diff --git a/scripts/python/phoedf.py b/scripts/python/phoedf.py index ae8ee698c..8a0b5e9a6 100644 --- a/scripts/python/phoedf.py +++ b/scripts/python/phoedf.py @@ -75,7 +75,7 @@ def flatten_indices(mu, nu): self.gcov = np.zeros([self.NumBlocks, 4, 4, self.Nx3, self.Nx2, self.Nx1]) for mu in range(4): for nu in range(4): - self.gcov[:,mu,nu,:,:,:] = self.flatgcov[:,flatten_indices(mu,nu),:,:,:] + self.gcov[:,mu,nu,:,:,:] = self.flatgcov[:,flatten_indices(mu,nu),:,4:-4,4:-4] del(self.flatgcov) self.gcon = np.zeros([self.NumBlocks, 4, 4, self.Nx3, self.Nx2, self.Nx1]) for b in range(self.NumBlocks): diff --git a/scripts/python/plot_torus.py b/scripts/python/plot_torus.py index 285c565c7..5b4bbfb5a 100644 --- a/scripts/python/plot_torus.py +++ b/scripts/python/plot_torus.py @@ -107,7 +107,7 @@ def plot_frame(ifname, fname, savefig, geomfile=None, rlim=40, coords='cartesian ax = axes[0,2] Pg = dfile.GetPg() - Pm = np.clip(dfile.GetPm(), 1.e-20, 1.e20) + Pm = np.clip(dfile.GetPm(), 1.e-100, 1.e100) lbeta = np.log10(Pg/Pm) for b in range(nblocks): im = ax.pcolormesh(xplot[b,:,:], yplot[b,:,:], lbeta[b,0,:,:].transpose(), vmin=-2, vmax=2, @@ -182,7 +182,7 @@ def plot_frame(ifname, fname, savefig, geomfile=None, rlim=40, coords='cartesian ax = axes[1,1] Pg = dfile.GetPg() - Pm = np.clip(dfile.GetPm(), 1.e-20, 1.e20) + Pm = np.clip(dfile.GetPm(), 1.e-100, 1.e100) lbeta = np.log10(Pg/Pm) for b in range(nblocks): im = ax.pcolormesh(xplot[b,:,:], yplot[b,:,:], lbeta[b,0,:,:].transpose(), vmin=-2, vmax=2, @@ -217,6 +217,7 @@ def plot_frame(ifname, fname, savefig, geomfile=None, rlim=40, coords='cartesian parser.add_argument('--savefig', type=bool, default=False, help='Whether to save figure') parser.add_argument('--numin', type=float, default=1.e-4, help='Minimum frequency (Hz)') parser.add_argument('--numax', type=float, default=1.e2, help='Maximum frequency (Hz)') + parser.add_argument('--rlim', type=float, default=40., help='Maximum radius to plot') parser.add_argument('--nnu', type=int, default=100, help='Number of frequency support points') parser.add_argument('--coords', type=str, default="cartesian", help='Coordinates to plot in') parser.add_argument('files', type=str, nargs='+', @@ -226,4 +227,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 4921b0e8a..0eae7c6c7 100644 --- a/src/geometry/fmks.cpp +++ b/src/geometry/fmks.cpp @@ -47,7 +47,7 @@ void Initialize(ParameterInput *pin, StateDescriptor *geometry) { Real x0 = pin->GetReal("parthenon/mesh", "x1min"); Real smooth = pin->GetOrAddReal("coordinates", "smooth", 0.5); Real hexp_br = pin->GetOrAddReal("coordinates", "hexp_br", 1000.); - Real hexp_nsq = pin->GetOrAddReal("coordinates", "hexp_nsq", 1.0); + 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); diff --git a/src/geometry/mckinney_gammie_ryan.hpp b/src/geometry/mckinney_gammie_ryan.hpp index 3131ad07e..d6bf79c56 100644 --- a/src/geometry/mckinney_gammie_ryan.hpp +++ b/src/geometry/mckinney_gammie_ryan.hpp @@ -46,10 +46,10 @@ class McKinneyGammieRyan { public: McKinneyGammieRyan() : derefine_poles_(true), h_(0.3), xt_(0.82), alpha_(14), x0_(0), smooth_(0.5), - norm_(GetNorm_(alpha_, xt_)), hexp_br_(1000.), hexp_nsq_(1.0), hexp_csq_(4.0) {} + 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_)), hexp_br_(1000.), hexp_nsq_(1.0), hexp_csq_(4.0) {} + 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), @@ -58,7 +58,7 @@ class McKinneyGammieRyan { 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_)), hexp_br_(1000.), hexp_nsq_(1.0), + 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], @@ -75,7 +75,8 @@ class McKinneyGammieRyan { // 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_)); + 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_) { From df17579e87c797dec763c04523a5b6c122b2acee Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Thu, 20 Apr 2023 12:48:10 -0600 Subject: [PATCH 4/7] Add hyperexp unit test --- tst/unit/geometry/test_geometry.cpp | 30 +++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/tst/unit/geometry/test_geometry.cpp b/tst/unit/geometry/test_geometry.cpp index 520647033..48532dbf2 100644 --- a/tst/unit/geometry/test_geometry.cpp +++ b/tst/unit/geometry/test_geometry.cpp @@ -593,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); From 1da83200b03798aec19dc0bdd56230821051b503 Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Thu, 4 May 2023 11:42:27 -0600 Subject: [PATCH 5/7] dont mess with parthenon version --- external/parthenon | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/parthenon b/external/parthenon index eb16b502d..04eceab86 160000 --- a/external/parthenon +++ b/external/parthenon @@ -1 +1 @@ -Subproject commit eb16b502dc03d5288ef7ea06fb3399b480e17374 +Subproject commit 04eceab86606ae6dc38d1e2fcf1d2b0def501fc7 From 44a84b08c7b1abe94d790c347b8fdccef83ac109 Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Thu, 4 May 2023 12:02:41 -0600 Subject: [PATCH 6/7] Fix defaults and make black less verbose --- scripts/bash/format.sh | 3 ++- scripts/python/plot_torus.py | 2 ++ src/geometry/fmks.cpp | 2 +- src/geometry/mckinney_gammie_ryan.hpp | 2 +- 4 files changed, 6 insertions(+), 3 deletions(-) diff --git a/scripts/bash/format.sh b/scripts/bash/format.sh index 23943ff80..47cc6f89d 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 c57edc572..efc596ceb 100644 --- a/scripts/python/plot_torus.py +++ b/scripts/python/plot_torus.py @@ -42,6 +42,8 @@ def plot_frame(ifname, fname, savefig, geomfile=None, rlim=40, coords="cartesian dfile = phoedf(fname) + rad_active = dfile.Params["radiation/active"] + a = dfile.Params["geometry/a"] hslope = dfile.Params["geometry/h"] diff --git a/src/geometry/fmks.cpp b/src/geometry/fmks.cpp index 0eae7c6c7..4eb7cf377 100644 --- a/src/geometry/fmks.cpp +++ b/src/geometry/fmks.cpp @@ -46,7 +46,7 @@ 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); - Real hexp_br = pin->GetOrAddReal("coordinates", "hexp_br", 1000.); + Real hexp_br = pin->GetOrAddReal("coordinates", "hexp_br", 1.e6); 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"); diff --git a/src/geometry/mckinney_gammie_ryan.hpp b/src/geometry/mckinney_gammie_ryan.hpp index d6bf79c56..10c695e36 100644 --- a/src/geometry/mckinney_gammie_ryan.hpp +++ b/src/geometry/mckinney_gammie_ryan.hpp @@ -204,7 +204,7 @@ class McKinneyGammieRyan { Real smooth_ = 0.5; Real norm_; Real hexp_br_ = 1000.; - Real hexp_nsq_ = 1.; + Real hexp_nsq_ = 4.; Real hexp_csq_ = 4.; }; From dcb59ebbe7d6abca377f8a2baca4a0e12477fe8e Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Thu, 4 May 2023 13:30:39 -0600 Subject: [PATCH 7/7] Larger initial break radius --- inputs/torus.pin | 1 + src/geometry/fmks.cpp | 5 ++++- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/inputs/torus.pin b/inputs/torus.pin index 82dd178fb..99a970fe9 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/src/geometry/fmks.cpp b/src/geometry/fmks.cpp index 4eb7cf377..a8bbe80f7 100644 --- a/src/geometry/fmks.cpp +++ b/src/geometry/fmks.cpp @@ -46,7 +46,10 @@ 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); - Real hexp_br = pin->GetOrAddReal("coordinates", "hexp_br", 1.e6); + 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");