From d3a4aa330203000f2d2f73880778d1ef52b48b41 Mon Sep 17 00:00:00 2001 From: Wolfgang Kastaun Date: Fri, 9 Oct 2020 11:15:39 +0200 Subject: [PATCH] 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 --- docsrc/sphinx/index.rst | 4 +- library/Con2Prim_IMHD/con2prim_imhd.cc | 129 ++++++++++++++++-- .../include/con2prim_imhd_internals.h | 46 ++++++- meson.build | 2 +- meson_options.txt | 1 + tests/benchmarks/scripts/plot_benchmark.py | 29 ++-- .../benchmarks/src/benchmark_con2prim_mhd.cc | 16 ++- tests/meson.build | 5 +- tests/src/test_con2prim_mhd.cc | 9 +- tests/src/test_eos.cc | 1 - tests/src/test_utils.cc | 1 - 11 files changed, 209 insertions(+), 34 deletions(-) diff --git a/docsrc/sphinx/index.rst b/docsrc/sphinx/index.rst index 1e548b2..c8474d1 100644 --- a/docsrc/sphinx/index.rst +++ b/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] `_ +(2020). Prospective users should be sure to check for updated versions that will likely be made available in a public repository. diff --git a/library/Con2Prim_IMHD/con2prim_imhd.cc b/library/Con2Prim_IMHD/con2prim_imhd.cc index 6e8fdf1..7e4fa9f 100644 --- a/library/Con2Prim_IMHD/con2prim_imhd.cc +++ b/library/Con2Prim_IMHD/con2prim_imhd.cc @@ -7,7 +7,6 @@ #include #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 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,10 +342,26 @@ 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) { @@ -354,11 +369,22 @@ void con2prim_mhd::operator()(prim_vars_mhd& pv, cons_vars_mhd& cv, 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) { @@ -366,12 +392,6 @@ void con2prim_mhd::operator()(prim_vars_mhd& pv, cons_vars_mhd& cv, 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 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 ibracket, + const interval 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{muc0, muc1}; +} diff --git a/library/Con2Prim_IMHD/include/con2prim_imhd_internals.h b/library/Con2Prim_IMHD/include/con2prim_imhd_internals.h index 732d17f..5119b9f 100644 --- a/library/Con2Prim_IMHD/include/con2prim_imhd_internals.h +++ b/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; +}; + +///Class for handling rare corner case +class rarecase { + public: + rarecase( + const interval bracket, ///< Initial master root bracket + const interval rgrho, ///< Allowed density range + const froot& f ///< Master root function + ); + + /// Root bracket on which solution is unique + interval 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 +}; + } } diff --git a/meson.build b/meson.build index 5c38f5f..b3a1ea3 100644 --- a/meson.build +++ b/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' diff --git a/meson_options.txt b/meson_options.txt index 7258aea..5f85614 100644 --- a/meson_options.txt +++ b/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) diff --git a/tests/benchmarks/scripts/plot_benchmark.py b/tests/benchmarks/scripts/plot_benchmark.py index 3f7067b..cb59cad 100755 --- a/tests/benchmarks/scripts/plot_benchmark.py +++ b/tests/benchmarks/scripts/plot_benchmark.py @@ -5,7 +5,7 @@ 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 @@ -13,11 +13,11 @@ def plot_perf_zeps(path, eosname, dat_bz, dat_bl): 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) # # diff --git a/tests/benchmarks/src/benchmark_con2prim_mhd.cc b/tests/benchmarks/src/benchmark_con2prim_mhd.cc index 07028c0..f5517c0 100644 --- a/tests/benchmarks/src/benchmark_con2prim_mhd.cc +++ b/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 < @@ -195,7 +194,11 @@ bool test_con2prim_mhd::compare_prims( const real_t f1 = cs1*cs1 * (1.0 + a1) / a1; const real_t f = max(f0,f1); - hope.isclose(pv0.press, pv1.press, f*acc + acc_eps, 0, + const real_t F = (pv0.eps/pv0.press) + * eos.at_rho_eps_ye(pv0.rho, + pv0.eps, pv0.ye).dpress_deps(); + + hope.isclose(pv0.press, pv1.press, f*acc + F*acc_eps, 0, "press"); } hope(dv <= v * (acc/wsqr + acc_max), @@ -738,7 +741,7 @@ BOOST_AUTO_TEST_CASE( test_con2prim_phys_hybr ) const real_t ye_fixed{0.25}; auto iz = log_spacing(min_z, max_z, 100); - auto ith = linear_spacing(1e-4, 0.99, 10); + auto ith = linear_spacing(2e-6, 0.99, 10); //2e-6 = 1e-4 / epsmax auto irho = log_spacing(par.atmo_rho*5, tst.eos.range_rho().max()/1.1, 10); auto ib = linear_spacing(0.0, max_b, 10); diff --git a/tests/src/test_eos.cc b/tests/src/test_eos.cc index 6f071af..6efe0e4 100644 --- a/tests/src/test_eos.cc +++ b/tests/src/test_eos.cc @@ -1,4 +1,3 @@ -#define BOOST_TEST_DYN_LINK #define BOOST_TEST_MODULE EOS #include diff --git a/tests/src/test_utils.cc b/tests/src/test_utils.cc index d690ec8..a2d8176 100644 --- a/tests/src/test_utils.cc +++ b/tests/src/test_utils.cc @@ -1,4 +1,3 @@ -#define BOOST_TEST_DYN_LINK #include #include "test_utils.h"