Version 1.1 of the ReprimAnd library.

This is the exact version described in the final paper as published in Phys. Rev. D
Changes:
Recovery algorithm now takes into account rare corner case
Improvements and Bugfixes in benchmark/accuracy plots
Improvements of unit tests

Author: Wolfgang Kastaun

docsrc/sphinx/index.rst

@@ -20,7 +20,9 @@ describing the recovery scheme:
 
 Wolfgang Kastaun, Jay Vijay Kalinani, Riccardo Ciolfi, 
 "Robust Recovery of Primitive Variables in Relativistic Ideal 
-Magnetohydrodynamics" (2020).
+Magnetohydrodynamics",  
+`arXiv:2005.01821 [gr-qc] <https://arxiv.org/abs/2005.01821>`_ 
+(2020).
 
 Prospective users should be sure to check for updated versions
 that will likely be made available in a public repository.

library/Con2Prim_IMHD/con2prim_imhd.cc

@@ -7,7 +7,6 @@
 #include <cmath>
 #include "find_roots.h"
 
-
 using namespace EOS_Toolkit;
 using namespace EOS_Toolkit::detail;
 using namespace std;
@@ -252,7 +251,7 @@ auto froot::initial_bracket(report& errs) const -> interval<real_t>
     f_upper g(h0, rsqr, rbsqr, bsqr);
 
     ROOTSTAT status;
-    mu_max = findroot_using_deriv(g, status, ndigits2);
+    mu_max = findroot_using_deriv(g, status, ndigits2, ndigits2+4);
 
     if (status != ROOTSTAT::SUCCESS) {
       if (status == ROOTSTAT::NOT_CONVERGED) {
@@ -343,35 +342,56 @@ void con2prim_mhd::operator()(prim_vars_mhd& pv, cons_vars_mhd& cv,
     set_to_nan(pv, cv);
     return;
   }
+  
+  rarecase nc(bracket, eos.range_rho(), f);
+
+  if (nc.rho_too_big) 
+  {
+    errs.set_range_rho(d, d);
+    set_to_nan(pv, cv);
+    return;
+  }
 
+  if (nc.rho_too_small) 
+  {
+    errs.set_atmo_set();
+    atmo.set(pv, cv, g);
+    return;
+  }
+
+  
   ROOTSTAT status;  
-  bracket = findroot_no_deriv(f, bracket, acc, max_iter, status);
-  assert(bracket.contains(sol.lmu));
+  bracket = findroot_no_deriv(f, nc.bracket, acc, max_iter, status);
 
   errs.iters = sol.calls;
   if (status != ROOTSTAT::SUCCESS) {
     if (status == ROOTSTAT::NOT_CONVERGED) {
       errs.set_root_conv();
     }
     else if (status == ROOTSTAT::NOT_BRACKETED) {
+      if (nc.rho_big) { //That's why
+        errs.set_range_rho(d, d);
+        set_to_nan(pv, cv);
+        return;
+      }
+      if (nc.rho_small) {
+        errs.set_atmo_set();
+        atmo.set(pv, cv, g);
+        return;
+      }
       errs.set_root_bracket();
     }
     set_to_nan(pv, cv);
     return;
   }
+  assert(bracket.contains(sol.lmu));
 
 
   if (sol.rho < atmo.rho_cut) {
     errs.set_atmo_set();
     atmo.set(pv, cv, g);
     return;
   }
-  
-  if (sol.rho_raw > eos.range_rho()) {
-    errs.set_range_rho(d, sol.rho);
-    set_to_nan(pv, cv);
-    return;
-  }
 
   auto rgeps = eos.range_eps(sol.rho, sol.ye);
   if (sol.eps_raw > rgeps) {
@@ -430,3 +450,92 @@ void con2prim_mhd::operator()(prim_vars_mhd& pv, cons_vars_mhd& cv,
 }
 
 
+
+
+
+f_rare::f_rare(real_t wtarg_, const froot& f_)
+: v2targ(1.0 - 1.0/(wtarg_*wtarg_)), f(f_) {}
+
+
+/**
+This implements a root function for finding mu from W_hat
+**/
+auto f_rare::operator()(const real_t mu) const 
+-> std::pair<real_t, real_t>
+{
+  real_t x     = f.x_from_mu(mu);
+  real_t xsqr  = x*x;
+  real_t rfsqr = f.rfsqr_from_mu_x(mu,x);
+  real_t vsqr  = mu*mu * rfsqr;  
+  
+  real_t y     = vsqr - v2targ;
+  real_t dy = 2.0*mu*x*(xsqr * f.rsqr + mu * (xsqr+x+1.0) * f.rbsqr);
+  
+  return {y, dy};
+}
+
+
+/**
+This checks for the corner case where the density might cross the 
+allowed bounds while finding the root of the master function. In
+that case, a tighter initial root finding interval is constructed
+to guarantee uniqueness.
+**/
+rarecase::rarecase(const interval<real_t> ibracket,
+                   const interval<real_t> rgrho, const froot& f)
+{
+  
+  real_t muc0 = ibracket.min();
+  real_t muc1 = ibracket.max();
+  const int ndigits = 30;
+          
+  
+  if (f.d > rgrho.max()) {
+    real_t wc = f.d / rgrho.max(); 
+    if (wc  > f.winf) {
+      rho_too_big = true; 
+    }
+    else {
+      f_rare g(wc, f);
+      
+      if (g(muc1).first <= 0) {
+        rho_too_big = true; 
+      }
+      else {
+        if (g(muc0).first < 0) {
+          ROOTSTAT status;
+          real_t mucc = findroot_using_deriv(g, ibracket, 
+                                             status, ndigits, ndigits+2);
+          assert(status == ROOTSTAT::SUCCESS); 
+          muc0 = max(muc0, mucc);
+          rho_big = true;
+        }
+      }
+    }
+  }
+
+  if (f.d <  f.winf * rgrho.min()) {
+    real_t wc = f.d / rgrho.min(); 
+    if (wc < 1) {
+      rho_too_small = true; 
+    }
+    else {
+      f_rare g(wc, f);
+      if (g(muc0).first >= 0) {
+        rho_too_small = true; 
+      }
+      else {
+        if (g(muc1).first > 0) {
+          ROOTSTAT status;
+          real_t mucc = findroot_using_deriv(g, ibracket, 
+                                             status, ndigits, ndigits+2);
+          assert(status == ROOTSTAT::SUCCESS); 
+          muc1 = min(muc1, mucc);
+          rho_small = true;
+        }
+      }
+    }
+  }
+
+  bracket = interval<real_t>{muc0, muc1};
+}

library/Con2Prim_IMHD/include/con2prim_imhd_internals.h

@@ -10,6 +10,9 @@
 namespace EOS_Toolkit {
 namespace detail {
 
+class rarecase;
+class f_rare;
+
 /// Function object representing the root function.
 /** This contains all the fixed parameters defining the function.
     It also remembers intermediate results from the last evaluation,
@@ -45,10 +48,14 @@ class froot {
   static real_t get_eps_raw(real_t mu,  real_t qf, 
                             real_t rfsqr, real_t w);  
 
-   
+  friend class rarecase;
+  friend class f_rare;
+     
   public:
 
   using value_t = real_t; 
+
+
 
   ///Store intermediate results obtained when evaluating function.
   struct cache {
@@ -134,6 +141,43 @@ class f_upper {
   const real_t bsqr;    ///< Fixed parameter \f$ b^2 = \frac{B^2}{D} \f$
 };
 
+
+
+
+///Root function used by rarecase class   
+class f_rare {
+  const real_t v2targ;  ///< Target squared velocity   
+  const froot& f;
+
+  public:
+  using value_t = real_t;
+  
+  f_rare(
+    real_t wtarg_,      ///< Target Lorentz factor
+    const froot& f_
+  );
+
+  auto operator()(real_t mu) const -> std::pair<real_t,real_t>;  
+};
+
+///Class for handling rare corner case 
+class rarecase {  
+  public:
+  rarecase(
+    const interval<real_t> bracket,   ///< Initial master root bracket 
+    const interval<real_t> rgrho,     ///< Allowed density range
+    const froot& f                    ///< Master root function
+  );
+
+  /// Root bracket on which solution is unique
+  interval<real_t> bracket;  
+  
+  bool rho_too_big { false };    ///< Density definitely too large
+  bool rho_big { false };        ///< Possibly too large
+  bool rho_too_small { false };  ///< Density definitely too small
+  bool rho_small { false };      ///< Possibly too small
+};
+
 }
 }
 

meson.build

@@ -1,6 +1,6 @@
 
 project('RePrimAnd', ['cpp','c'], default_options : ['cpp_std=c++11'],
-        version : '1.0', license : 'CC BY-NC-SA 4.0')
+        version : '1.1', license : 'CC BY-NC-SA 4.0')
 
 project_headers_dest = 'reprimand'
 

meson_options.txt

@@ -1,3 +1,4 @@
 
 option('build_documentation', type : 'boolean', value : true)
 option('build_benchmarks', type : 'boolean', value : true)
+option('build_tests', type : 'boolean', value : true)

tests/benchmarks/scripts/plot_benchmark.py

@@ -5,19 +5,19 @@
 import sys
 from stdplt import *
 
-
+GAUSS_GU = 1.1973228161170991e-20
 
 def plot_perf_zeps(path, eosname, dat_bz, dat_bl):
   figname = "perf_eos%s_z_eps" % eosname
 
   with Figure(os.path.join(path, figname)) as fig:
     grid = twopanel(fig)
 
-    maxit     = max(int(dat_bz[-1].max()), int(dat_bl[-1].max()))
+    maxit     = max(int(dat_bz[-1][2].max()), int(dat_bl[-1][2].max()))
     bticks    = np.arange(1,maxit, max(1,int(maxit/10)))
     clrmap    = get_color_map("viridis", over='r')
 
-    for (x0,x1,it),ax in zip([dat_bz, dat_bl], grid):
+    for (x0,x1,(z,eps,it,p,B)),ax in zip([dat_bz, dat_bl], grid):
       im = plot_color(it, x0, x1, vmin=0.999, vmax=maxit,
                  cmap=clrmap, bar=False, axes=ax)
       ax.set_aspect('auto')
@@ -33,12 +33,23 @@ def plot_perf_zb(path, eosname, dat_c, dat_h):
   with Figure(os.path.join(path, figname)) as fig:
     grid = twopanel(fig)
 
-    maxit    = max(int(dat_c[-1].max()), int(dat_h[-1].max()))
+    maxit    = max(int(dat_c[-1][2].max()), int(dat_h[-1][2].max()))
     bticks   = np.arange(1,maxit, max(1,int(maxit/10)))
     clrmap    = get_color_map("viridis", over='r')
-    for (x0,x1,it),ax in zip([dat_c, dat_h], grid):
+    for (x0,x1,(z,b,it,p,B)),ax in zip([dat_c, dat_h], grid):
       im = plot_color(it, x0, x1, vmin=0.999, vmax=maxit, 
                       cmap=clrmap, bar=False, axes=ax)
+      
+      z1d = z[:,0]
+      b1d = b[0,:]
+      lgbeta = np.log10(B**2 / p)
+      cs = ax.contour(np.log10(z1d), np.log10(b1d), lgbeta.T, 
+                   levels=np.arange(-2,11,4), colors='y')
+
+      ax.clabel(cs, inline=True)
+      ax.contour(np.log10(z1d), np.log10(b1d), B.T, 
+                   levels=[1e16*GAUSS_GU], colors=['r'])
+
       ax.set_aspect('auto')
       ax.set_ylabel(r'$\log_{10}(B / \sqrt{D})$')
     grid[1].set_xlabel(r'$\log_{10}(z[c])$')
@@ -49,8 +60,8 @@ def plot_perf_zb(path, eosname, dat_c, dat_h):
 
 def load_data(path, eos):
     def read(t):
-      c0,c1,d = load_grid(os.path.join(path, t % eos), 3)
-      return c0,c1,d[2]
+      c0,c1,d = load_grid(os.path.join(path, t % eos), 5)
+      return c0,c1,d
     #
     files = ["perf_eos%s_z_b_cold.dat", "perf_eos%s_z_b_hot.dat",
              "perf_eos%s_z_eps_Bzero.dat", 
@@ -63,8 +74,8 @@ def main():
   path = sys.argv[1]
   for eos in ['ig', 'hyb']:
     dat_c, dat_h, dat_bz, dat_bl = load_data(path, eos)
-    plot_perf_zeps(path, eos, dat_c, dat_h)
-    plot_perf_zb(path, eos, dat_bz, dat_bl)
+    plot_perf_zeps(path, eos, dat_bz, dat_bl)
+    plot_perf_zb(path, eos, dat_c, dat_h)
   #
 #
 

tests/benchmarks/src/benchmark_con2prim_mhd.cc

@@ -75,7 +75,10 @@ void map_z_eps(eos_thermal eos, con2prim_mhd& cv2pv, real_t b, real_t rho,
       }
       assert(!rep.failed());
       os << setw(20) << z << setw(20) << epsth 
-         << setw(20) << rep.iters << endl;
+         << setw(20) << rep.iters 
+         << setw(20) << pv.press 
+         << setw(20) << b*sqrt(cv.dens)
+         << endl;
     }
     os << endl;
   }
@@ -85,7 +88,7 @@ void map_z_b(eos_thermal eos, con2prim_mhd& cv2pv, real_t epsth,
              real_t rho, real_t ye, sm_metric3& g, string fn) 
 {
   auto iz     = log_spacing(1e-2, 1e3, 200);
-  auto ib     = log_spacing(1e-4, 1e1, 200);
+  auto ib     = log_spacing(1e-4, 1e4, 200);
 
   const real_t eps0  = eos.range_eps(rho, ye).min();
   const real_t eps   = eps0 + epsth;
@@ -104,7 +107,10 @@ void map_z_b(eos_thermal eos, con2prim_mhd& cv2pv, real_t epsth,
         assert(false);
       }
       assert(!rep.failed());
-      os << setw(20) << z << setw(20) << b <<setw(20) << rep.iters << endl;
+      os << setw(20) << z << setw(20) << b <<setw(20) << rep.iters 
+         << setw(20) << pv.press 
+         << setw(20) << b*sqrt(cv.dens)
+         << endl;
     }
     os << endl;
   }
@@ -144,14 +150,14 @@ int main(int argc, char *argv[])
   eos_thermal eos_ig = get_eos_ig();
   atmosphere atmo_ig = get_atmo(eos_ig, 1e-6);
   const real_t acc          = 1e-8;  
-  con2prim_mhd cv2pv_ig(eos_ig, atmo_ig.rho, false, 2e3, 1e2, 
+  con2prim_mhd cv2pv_ig(eos_ig, atmo_ig.rho, false, 2e3, 5e4, 
                         atmo_ig, acc, 100);
 
 
   auto eos_hyb = load_eos_thermal(PATH_EOS_HYB, units::geom_solar());
 
   atmosphere atmo_hyb = get_atmo(eos_hyb, 0.0);
-  con2prim_mhd cv2pv_hyb(eos_hyb, atmo_hyb.rho, false, 2e3, 1e2, 
+  con2prim_mhd cv2pv_hyb(eos_hyb, atmo_hyb.rho, false, 2e3, 5e4, 
                          atmo_hyb, acc, 100);