[Users] Proposed modification to the EOS_Omni thorn
Federico Cipolletta
cipo87 at gmail.com
Wed May 6 04:12:43 CDT 2020
We propose the attached patch for the file
EOS_Omni/src/nuc_eos_cxx/helpers.hh of EOS_Omni in order to fix a bug
that causes GRHydro to fail in some situations involving tabulated
EOSs.
We are currently working on the development of our GRMHD code,
"Spritz", which will be capable of dealing with tabulated EOS,
considering the evolution of Temperature, T, and Electron Fraction,
Y_e. In order to test the implementation, we are currently performing
some TOV tests with both Spritz and GRHydro, adopting the LS220
tabulated EOS that we downloaded from the stellarcollapse website. We
realized that the EOS_Omni thorn may present some issues when dealing
with tabulated nuclear EOS. In particular, when T must be found in
relation to other thermodynamic quantities, this is normally done by
the EOS_Omni thorn by computing T through a Newton-Raphson routine
with a "fall-back" on Bisection method in the case of too many
iterations (the procedure is related to the routines
"nuc_eos_findtemp" and "bisection" contained in the file
"EOS_Omni/src/nuc_eos_cxx/helpers.hh" - if the number of iterations
performed by Newton-Raphson is greater that one hard-coded 200 value,
then the Bisection method is called). With our tests for the LS220
EOS, we do see that the calculation of T, in some cases, may not
perform well. Specifically, with the GRHydro code, roughly after 144
iterations the T goes above the prescribed parameter
"GRHydro_max_temp" (which we imagine being chosen in accordance with
the maximum value of T given in the Table), so the simulation is
stopped. With Spritz, as well, we clearly get one wrong T value (well
above T_max given in the table).
This issue is probably related to a weak dependence of T on \epsilon
and can be solved by using the bisection method, see Section IV of
Galeazzi, Kastaun, Rezzolla, and Font, PRD (2013). Lorenzo Sala, one
of our collaborators, already implemented this modification in the
aforementioned routines of the EOS_Omni thorn. In the modified
version, if the energy dependence on the temperature is weak, i.e., if
the derivative of the function is too large and therefore the
Newton-Raphson routine may overshoot (in the routine
"nuc_eos_findtemp"), the code switches to a safer bisection method (in
the routine "bisection"). In addition, the "bisection" routine
presents a slight modification in order to be sure of not going below
the minimum of tabulated temperature value. With these modifications,
both tests with GRHydro and Spritz result to be successful and well in
accordance.
We would like to propose the inclusion of the modified EOS_Omni thorn
in the future releases of the EinsteinToolkit, in order to provide the
community with this more robust version. For your reference, I attach:
1- the file "helpers.hh" containing our proposed modifications;
2- a patch file "diff_STDvsMOD.txt" where I compared the standard
(STD) vs the modified (MOD) helpers.hh file
3 - a file "LS220_RHOmax.pdf" with a plot showing the evolution of
maximum density for the TOV test with the LS220 EOS and comparing the
results obtained with Spritz and GRHydro plus the standard and
modified EOS_Omni thorn;
4- a file "LS220_Tmax.pdf", which is the same as 3 but for maximum temperature.
Best Regards,
Federico Cipolletta.
-------------- next part --------------
--- STD/helpers.hh 2020-04-30 09:19:51.575537041 +0200
+++ MOD/helpers.hh 2020-05-02 14:36:21.035245866 +0200
@@ -1,5 +1,8 @@
+// helpers.hh edited by Lorenzo Sala
+
#include "nuc_eos.hh"
#include <cstdlib>
+#include <iostream>
namespace nuc_eos {
@@ -368,7 +371,7 @@
const double dltp = log(1.2);
const double dltm = log(0.8);
- const double leps0_prec = leps0*prec;
+ double leps0_prec = fabs(leps0*prec);
// temporary local vars
double lt, lt1, lt2;
@@ -380,8 +383,31 @@
double delx,dely,delz;
int idx[8];
+ // LSMOD (Modification made by Lorenzo Sala)
+ // LSMOD: The following lines calculate eps in
+ // f2a = eps(rho,Tmin, Ye) and f1a = eps(rho,Tmax,Ye)
+ get_interp_spots(lr,ltmax,ye,&delx,&dely,&delz,idx);
+ nuc_eos_C_linterp_one(idx,delx,dely,delz,&f1a,iv);
+ get_interp_spots(lr,ltmin,ye,&delx,&dely,&delz,idx);
+ nuc_eos_C_linterp_one(idx,delx,dely,delz,&f2a,iv);
// prepare
+ // check if your energy is actually tabulated at this rho and ye.
+ // f2a is the energy evaluated at ltmin, so it is the minimum energy tabulated
+ // at this rho ad ye.
+ // If leps0 <= f2a, then ltout is likely to be the minimum temperature tabulated.
+ if(leps0 <= f2a) { // + 1.0E-6
+ *ltout = ltmin;
+ return;
+ }
+
+ /* // If leps0 >= f1a, then ltout is likely to be the maximum temperature tabulated.
+ if(leps0 >= f1a) { // + 1.0E-6
+ *ltout = ltmax;
+ return;
+ } */
+
+ // otherwise, proceed finding extrema for applying bisection method.
lt = lt0;
lt1 = MIN(lt0 + dlt0p,ltmax);
lt2 = MAX(lt0 + dlt0m,ltmin);
@@ -416,6 +442,10 @@
bcount++;
if(CCTK_BUILTIN_EXPECT(bcount >= maxbcount, false)) {
+#if DEBUG
+ fprintf(stderr,"bcount out of range it %d, lr: %15.6E, lt1: %15.6E, lt2: %15.6E, f1a: %18.11E, f2a: %18.11E leps0: %18.11E, ye: %15.6E\n",
+ bcount,lr,lt1,lt2,f1a,f2a,leps0,ye);
+#endif
*keyerrt = 667;
return;
}
@@ -510,7 +540,12 @@
const int irho = MIN(MAX(1 + (int)(( lr - logrho[0] - 1.0e-12) * drhoi),1),nrho-1);
const int iye = MIN(MAX(1 + (int)(( ye - yes[0] - 1.0e-12) * dyei),1),nye-1);
- while(it < itmax) {
+ /* ******* if temp low for high density, switch directly to bisection.
+ Verifying Newton-Raphson result evaluating the derivative.
+ The variable shouldgotobisection will be modified accordingly
+ to the value of derivative of eps wrt temp ******* */
+ bool shouldgotobisection = false; // LSMOD
+ while(it < itmax && shouldgotobisection == false) {
it++;
// step 2: check if the two bounding values of the temperature
@@ -572,6 +607,13 @@
dlepsdlti = (logtemp[itemp]-logtemp[itemp-1])/(epst2-epst1);
ldt = -(leps - leps0) * dlepsdlti * fac;
+ //LSMOD: too large a dlt means that the energy dependence on the temperature
+ // is weak ==> We'd better try bisection.
+ // Factor 1/12.0 come from tests by LSMOD
+ // This is done in order to limit the "velocity" of T variation
+ // given by Newton-Raphson.
+ if(ldt > (ltmax-ltmin) / 12.0 ) shouldgotobisection = true;
+
ltn = MIN(MAX(lt + ldt,ltmin),ltmax);
lt = ltn;
@@ -605,9 +647,9 @@
bisection(lr,lt0,ye,leps0,prec,ltout,1,keyerrt);
#if DEBUG
if(*keyerrt==667) {
- fprintf(stderr,"This is worse. Bisection failed!\n");
- abort();
- }
+ fprintf(stderr,"This is worse. Bisection failed!\n");
+ abort();
+ }
#endif
-------------- next part --------------
A non-text attachment was scrubbed...
Name: LS220_RHOmax.pdf
Type: application/pdf
Size: 24045 bytes
Desc: not available
Url : http://lists.einsteintoolkit.org/pipermail/users/attachments/20200506/577bbd68/attachment-0002.pdf
-------------- next part --------------
A non-text attachment was scrubbed...
Name: LS220_Tmax.pdf
Type: application/pdf
Size: 19419 bytes
Desc: not available
Url : http://lists.einsteintoolkit.org/pipermail/users/attachments/20200506/577bbd68/attachment-0003.pdf
-------------- next part --------------
A non-text attachment was scrubbed...
Name: helpers.hh
Type: text/x-c++hdr
Size: 22878 bytes
Desc: not available
Url : http://lists.einsteintoolkit.org/pipermail/users/attachments/20200506/577bbd68/attachment-0001.bin
More information about the Users
mailing list