Merge pull request #166 from lanl/brryan/hyperexp

Hyperexponential coordinates
lanl · May 4, 2023 · ca02a87 · ca02a87
2 parents ca15cb2 + dcb59eb
commit ca02a87
Show file tree

Hide file tree

Showing 7 changed files with 76 additions and 10 deletions.
diff --git a/inputs/torus.pin b/inputs/torus.pin
@@ -124,6 +124,7 @@ do_fd_on_grid = false
 <coordinates>
 derefine_poles = false
 r_outer = 250
+# hexp_br = 1000. # Initial break radius for hyperexponential coordinates
 
 <fixup>
 enable_floors = true

diff --git 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
@@ -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
@@ -46,6 +46,12 @@ void Initialize<FMKSMeshBlock>(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<Real>::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<FMKSMeshBlock>(ParameterInput *pin, StateDescriptor *geometry) {
   x0 = std::min<Real>(x0, std::max<Real>(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
@@ -37,7 +37,11 @@ McKinneyGammieRyan GetTransformation<McKinneyGammieRyan>(StateDescriptor *pkg) {
   Real alpha = pkg->Param<Real>("alpha");
   Real x0 = pkg->Param<Real>("x0");
   Real smooth = pkg->Param<Real>("smooth");
-  return McKinneyGammieRyan(derefine_poles, h, xt, alpha, x0, smooth);
+  Real hexp_br = pkg->Param<Real>("hexp_br");
+  Real hexp_nsq = pkg->Param<Real>("hexp_nsq");
+  Real hexp_csq = pkg->Param<Real>("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
@@ -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
@@ -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);