GRHydro: Fixing Con2PrimM to call the proper point-wise routines

This re-introduces routines that work for hybrid/hot EOS and
corresponding changes in pointwise routines for hot EOS error checking
and temperature treatment by adding old EOSOmni pointwise routine.

From: Philipp Moesta <pmoesta@tapir.caltech.edu>

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@511 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45

1 files changed, 381 insertions, 67 deletions

diff --git a/src/GRHydro_Con2PrimM_pt_EOSOmni.c b/src/GRHydro_Con2PrimM_pt_EOSOmni.c
index 1d7e15c..96e2bd8 100644
--- a/src/GRHydro_Con2PrimM_pt_EOSOmni.c
+++ b/src/GRHydro_Con2PrimM_pt_EOSOmni.c
@@ -61,7 +61,7 @@
 #define NEWT_TOL2       (1.0e-15)  /* TOL of new 1D^*_{v^2} gnr2 method */
 #define MIN_NEWT_TOL2   (1.0e-10)  /* TOL of new 1D^*_{v^2} gnr2 method */
 
-#define W_TOO_BIG        (1.e20)    /* \gamma^2 (\rho_0 + u + p) is assumed
+#define W_TOO_BIG       (1.e20)    /* \gamma^2 (\rho_0 + u + p) is assumed
                                       to always be smaller than this.  This
                                       is used to detect solver failures */
 
@@ -73,13 +73,15 @@
 
 /* Local Globals */
 struct LocGlob {
-  CCTK_REAL Bsq,QdotBsq,Qtsq,Qdotn,D,half_Bsq ;
+  CCTK_REAL Bsq,bsq,QdotBsq,VdotBsq,Qtsq,Qdotn,D,half_Bsq,tau ;
 } ;
 
 struct eosomnivars {
   CCTK_INT eoshandle,eoskeytemp;
-  CCTK_REAL eosprec,eos_y_e[1],eos_temp[1];
+  CCTK_REAL prec,eosprec,eos_y_e[1],eos_temp[1];
   CCTK_INT eoskeyerr[1],eosanyerr[1];
+  CCTK_INT eosrl,eosit,eosi,eosj,eosk;
+  CCTK_REAL eosx,eosy,eosz;
 } ;
 
 // Declarations: 
@@ -122,9 +124,9 @@ static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w,CCTK_REAL gammaeos)
 {
   return((gammaeos-1.)*(w - rho0)/gammaeos) ;
 }
-
+//CCTK_INT  *handle, CCTK_INT GRHydro_reflevel, CCTK_INT *keytemp, CCTK_REAL *prec, 
 void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) (
-    CCTK_INT  *handle, CCTK_INT *keytemp, CCTK_REAL *prec, 
+    CCTK_INT  *handle, CCTK_INT *rl, CCTK_INT *ii, CCTK_INT *jj, CCTK_INT *kk, CCTK_REAL *x, CCTK_REAL *y, CCTK_REAL *z, CCTK_INT *keytemp, CCTK_REAL *eos_prec, CCTK_REAL *prec,  
     CCTK_REAL *gamma_eos,
     CCTK_REAL *dens_in, 
     CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, 
@@ -172,7 +174,7 @@ INPUT:  (using GRHydro variable defintions)
    tau        =  \alpha^2 \sqrt(\gamma) T^{00} - D 
    g[x,y,z][x,y,x] = spatial metric corresponding to \gamma 
    u[x,y,z][x,y,z] = inverse of the spatial metric, g[x,y,z][x,y,x]
-   det        =  gamma
+   det        =  sqrt(\gamma)
    B[x,y,z]   =  Bvec[0,1,2] 
    bsq        = b^\mu b_\mu  
 
@@ -205,10 +207,8 @@ RETURN VALUE: of retval = (i*100 + j)  where
              ( used to be  5 -> failure: rho,uu <= 0   but now sets epsnegative to non-zero )
 
 **********************************************************************************/
-
-
 void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  ( 
-    CCTK_INT  *handle, CCTK_INT *keytemp, CCTK_REAL *prec,
+    CCTK_INT  *handle, CCTK_INT *rl, CCTK_INT *ii, CCTK_INT *jj, CCTK_INT *kk, CCTK_REAL *x, CCTK_REAL *y, CCTK_REAL *z, CCTK_INT *keytemp, CCTK_REAL *eos_prec, CCTK_REAL *prec,  
     CCTK_REAL *gamma_eos,
     CCTK_REAL *dens_in, 
     CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, 
@@ -234,14 +234,14 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   CCTK_REAL sx, sy, sz;
   CCTK_REAL usx, usy, usz;
   CCTK_REAL tau, dens, gammaeos;
-  CCTK_REAL QdotB;
+  CCTK_REAL QdotB,VdotB;
   CCTK_REAL rho0,u,p,w,gammasq,gamma,gtmp,W_last,W,vsq;
   CCTK_REAL g_o_WBsq, QdB_o_W;
   CCTK_REAL detg = (*det);
   CCTK_REAL sqrt_detg = sqrt(detg);
   CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; 
   CCTK_INT i,j, i_increase ;
-
+  CCTK_INT nf, nfudgemax, failwarnmode, failinfomode;
   DECLARE_CCTK_PARAMETERS;
   
   struct LocGlob lg; 
@@ -250,53 +250,70 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   eosvars.eoshandle = *handle;
   //printf("handle = %i\n",*handle);
   eosvars.eoskeytemp = *keytemp;
-  eosvars.eosprec = *prec;
+  eosvars.eosprec = *eos_prec;
+  eosvars.prec = *prec;
   eosvars.eos_y_e[0] = *y_e_in;
   eosvars.eos_temp[0] = *temp_in;
   eosvars.eoskeyerr[0] = 0;
   eosvars.eosanyerr[0] = 0;
-
+  eosvars.eosrl = *rl;
+  eosvars.eosit = 1;
+  eosvars.eosi = *ii;
+  eosvars.eosj = *jj;
+  eosvars.eosk = *kk;
+  eosvars.eosx = *x;
+  eosvars.eosy = *y;
+  eosvars.eosz = *z;
+ 
   gammaeos = *gamma_eos;
 
+  nf =0;
+  nfudgemax = 30; 
+  failinfomode = 1;
+  failwarnmode = 0;
+  
   /* Assume ok initially: */
   *retval = 0.;
   *epsnegative = 0; 
 
 #if(DEBUG_CON2PRIMM)
+  fprintf(stdout," *x1        = %26.16e \n", *x        );
+  fprintf(stdout," *y1        = %26.16e \n", *y         );
+  fprintf(stdout," *z1        = %26.16e \n", *z         );
   fprintf(stdout," *dens        = %26.16e \n", *dens_in         );
-  fprintf(stdout," *sx                 = %26.16e \n", *sx_in           );    
-  fprintf(stdout," *sy          = %26.16e \n", *sy_in            );    
-  fprintf(stdout," *sz                 = %26.16e \n", *sz_in           );    
-  fprintf(stdout," *tau         = %26.16e \n", *tau_in           );    
-  fprintf(stdout," *Bconsx         = %26.16e \n", *Bconsx_in );    
-  fprintf(stdout," *Bconsy         = %26.16e \n", *Bconsy_in );    
+  fprintf(stdout," *sx          = %26.16e \n", *sx_in     );    
+  fprintf(stdout," *sy          = %26.16e \n", *sy_in     );    
+  fprintf(stdout," *sz          = %26.16e \n", *sz_in     );    
+  fprintf(stdout," *tau         = %26.16e \n", *tau_in    );    
+  fprintf(stdout," *Bconsx      = %26.16e \n", *Bconsx_in );    
+  fprintf(stdout," *Bconsy      = %26.16e \n", *Bconsy_in );    
   fprintf(stdout," *Bconsz      = %26.16e \n", *Bconsz_in );    
-  fprintf(stdout," *rho         = %26.16e \n", *rho           );    
-  fprintf(stdout," *velx        = %26.16e \n", *velx           );    
-  fprintf(stdout," *vely        = %26.16e \n", *vely          );    
-  fprintf(stdout," *velz        = %26.16e \n", *velz          );    
+  fprintf(stdout," *rho         = %26.16e \n", *rho       );    
+  fprintf(stdout," *velx        = %26.16e \n", *velx      );    
+  fprintf(stdout," *vely        = %26.16e \n", *vely      );    
+  fprintf(stdout," *velz        = %26.16e \n", *velz      );    
   fprintf(stdout," *epsilon     = %26.16e \n", *epsilon     );
   fprintf(stdout," *temp_in     = %26.16e \n", *temp_in     );
   fprintf(stdout," *y_e_in     = %26.16e \n", *y_e_in     );
   fprintf(stdout," *pressure    = %26.16e \n", *pressure    );
-  fprintf(stdout," *Bx                 = %26.16e \n", *Bx           );    
-  fprintf(stdout," *By                 = %26.16e \n", *By           );    
-  fprintf(stdout," *Bz                 = %26.16e \n", *Bz           );    
-  fprintf(stdout," *bsq                = %26.16e \n", *bsq          );    
+  fprintf(stdout," *Bx          = %26.16e \n", *Bx        );    
+  fprintf(stdout," *By          = %26.16e \n", *By        );    
+  fprintf(stdout," *Bz          = %26.16e \n", *Bz        );    
+  fprintf(stdout," *bsq         = %26.16e \n", *bsq       );    
   fprintf(stdout," *w_lorentz   = %26.16e \n", *w_lorentz   );
-  fprintf(stdout," *gxx         = %26.16e \n", *gxx           );    
-  fprintf(stdout," *gxy         = %26.16e \n", *gxy           );    
-  fprintf(stdout," *gxz         = %26.16e \n", *gxz           );    
-  fprintf(stdout," *gyy         = %26.16e \n", *gyy           );    
-  fprintf(stdout," *gyz         = %26.16e \n", *gyz           );    
-  fprintf(stdout," *gzz         = %26.16e \n", *gzz           );    
-  fprintf(stdout," *uxx         = %26.16e \n", *uxx           );    
-  fprintf(stdout," *uxy         = %26.16e \n", *uxy           );    
-  fprintf(stdout," *uxz                = %26.16e \n", *uxz          );    
-  fprintf(stdout," *uyy         = %26.16e \n", *uyy           );    
-  fprintf(stdout," *uyz         = %26.16e \n", *uyz           );    
-  fprintf(stdout," *uzz                = %26.16e \n", *uzz          );    
-  fprintf(stdout," *det                = %26.16e \n", *det          );    
+  fprintf(stdout," *gxx         = %26.16e \n", *gxx       );    
+  fprintf(stdout," *gxy         = %26.16e \n", *gxy       );    
+  fprintf(stdout," *gxz         = %26.16e \n", *gxz       );    
+  fprintf(stdout," *gyy         = %26.16e \n", *gyy       );    
+  fprintf(stdout," *gyz         = %26.16e \n", *gyz       );    
+  fprintf(stdout," *gzz         = %26.16e \n", *gzz       );    
+  fprintf(stdout," *uxx         = %26.16e \n", *uxx       );    
+  fprintf(stdout," *uxy         = %26.16e \n", *uxy       );    
+  fprintf(stdout," *uxz         = %26.16e \n", *uxz       );    
+  fprintf(stdout," *uyy         = %26.16e \n", *uyy       );    
+  fprintf(stdout," *uyz         = %26.16e \n", *uyz       );    
+  fprintf(stdout," *uzz         = %26.16e \n", *uzz       );    
+  fprintf(stdout," *det         = %26.16e \n", *det       );    
   fprintf(stdout," *epsnegative = %10d    \n", *epsnegative );
   fprintf(stdout," *retval      = %26.16e \n", *retval      );
   fflush(stdout);
@@ -329,6 +346,13 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
        (*gyz) * (*By) * (*Bz) );
     
   QdotB = (sx * (*Bx) + sy * (*By) + sz * (*Bz)) ;
+
+  VdotB = ((*gxx) * (*velx) + (*gxy) * (*vely) + (*gxz) * (*velz))* (*Bx) + 
+           ((*gxy) * (*velx) + (*gyy) * (*vely) + (*gyz) * (*velz))* (*By) +
+           ((*gxz) * (*velx) + (*gyz) * (*vely) + (*gzz) * (*velz))* (*Bz); 
+
+  lg.VdotBsq = VdotB*VdotB;
+
   lg.QdotBsq = QdotB*QdotB ;
 
   lg.Qdotn = -(tau + dens) ;
@@ -338,6 +362,8 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   lg.D = dens;
 
   lg.half_Bsq = 0.5*lg.Bsq;
+  
+  lg.tau = tau;
 
   /* calculate W from last timestep and use for guess */
   vsq = 
@@ -354,7 +380,8 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   }
   if(vsq < 0. || vsq > 1. ) {
     *retval = 2.;
-    fprintf(stdout," *retval      = %26.16e \n", *retval      );
+    fprintf(stdout," *retval (inside loop)      = %26.16e \n", *retval      );
+    fflush(stdout);
     return;
   }
 
@@ -368,18 +395,86 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   CCTK_REAL uold = u; 
  
   CCTK_REAL dum1,dum2;
+//  eosvars.eosanyerr[0]=5;
   p = pressure_rho0_eps_eosomni(rho0,*epsilon,&dum1,&dum2,&eosvars) ;  // EOSOMNI
-  CCTK_REAL p2 = pressure_rho0_u(rho0,u,gammaeos) ;  // EOS
+  //  p = pressure_rho0_u(rho0,u,gammaeos) ;  // EOS
+ 
+  w = rho0 + u + p;
 
-  if (fabs(p-p2) > 1.0e-8){
-    printf("p = %20.16g,  p2 = %20.16g, delta pressure = %20.16g \n",p, p2, p-p2);
-  }
+  // we need this for the hot EOS failsafe tau reset later 
 
-  w = rho0 + u + p ;
+  lg.bsq = lg.Bsq/(w*w) + lg.VdotBsq;
+  
+  W_last = w*gammasq;
+
+//  fprintf(stdout," anyerr[0] (before loop) = %d \n", eosvars.eosanyerr[0]);
+//  fprintf(stdout,"vsq: rho0=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,epsilon,eosvars->eos_temp[0],eosvars->eos_y_e[0]);
+//  fprintf(stdout,"rho0=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,*epsilon,eosvars.eos_temp[0],eosvars.eos_y_e[0]);
+//  fprintf(stdout,"keyerr=%d, anyerr=%d \n", eosvars.eoskeyerr[0], eosvars.eosanyerr[0]);
+//  #pragma omp critical 
+  if (eosvars.eosanyerr[0] != 0) {
+     fprintf(stdout,"inside 1: anyerr=%d\n",eosvars.eosanyerr[0]);
+     if (eosvars.eosrl >= GRHydro_c2p_warn_from_reflevel) {
+        fprintf(stdout,"inside 2: anyerr=%d\n",eosvars.eosanyerr[0]);
+        if(eosvars.eoskeyerr[0] == 667 && GRHydro_eos_hot_eps_fix != 0) {
+           fprintf(stdout,"inside 3: anyerr=%d, keyerr=%d\n",eosvars.eosanyerr[0],eosvars.eoskeyerr[0]);
+           // Handling of the case when no new temperature can be
+           // found for a given epsilon. The amount of times
+           // we get to this place should be monitored, as this
+           // 'failsafe' mode creates artificial heat.
+
+           nf = 0;
+           while(eosvars.eosanyerr[0] != 0 && nf<= nfudgemax) {
+              eosvars.eosanyerr[0] = 0;
+              tau = tau + rho0*abs((*epsilon))*0.05*w*w;
+              *epsilon = (*epsilon) + abs((*epsilon))*0.05;
+              u = (*epsilon) * rho0;
+              EOS_Omni_press(eosvars.eoshandle,eosvars.eoskeytemp,eosvars.eosprec,1,
+                 &rho0,epsilon,eosvars.eos_temp,
+                 eosvars.eos_y_e,&p,eosvars.eoskeyerr,eosvars.eosanyerr);
+              nf = nf + 1;
+           } 
+           fprintf(stdout, "Iterations of heat injection: %d \n",nf);
+           fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars.eosit,eosvars.eosrl,eosvars.eosi,eosvars.eosj,eosvars.eosk,eosvars.eosx,eosvars.eosy,eosvars.eosz);
+           if(nf >= nfudgemax) { 
+              if(GRHydro_c2p_reset_eps_tau_hot_eos != 1) {
+                 CCTK_WARN(failinfomode,"EOS error in c2p 0: injected heat too many times");
+                 fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars.eosit,eosvars.eosrl,eosvars.eosi,eosvars.eosj,eosvars.eosk,eosvars.eosx,eosvars.eosy,eosvars.eosz);
+                 //CCTK_WARN(failinfomode,warnline);
+                 fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,*epsilon,eosvars.eos_temp[0],eosvars.eos_y_e[0]);
+                 //call CCTK_WARN(failinfomode,warnline)
+                 fprintf(stdout,"code=%d\n",eosvars.eoskeyerr[0]);
+                 //call CCTK_WARN(failinfomode,warnline)
+                 CCTK_WARN(failwarnmode,"Aborting!!!");
+              } 
+              else {
+                 CCTK_WARN(failinfomode,"EOS error in c2p 0: LAST RESORT -- reset eps and tau based on old temp");
+                 fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars.eosit,eosvars.eosrl,eosvars.eosi,eosvars.eosj,eosvars.eosk,eosvars.eosx,eosvars.eosy,eosvars.eosz);
+                 fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,*epsilon,eosvars.eos_temp[0],eosvars.eos_y_e[0]);
+                 eosvars.eoskeytemp =1; 
+                 EOS_Omni_press(eosvars.eoshandle,eosvars.eoskeytemp,eosvars.eosprec,1,
+                 &rho0,epsilon,eosvars.eos_temp,eosvars.eos_y_e,&p,eosvars.eoskeyerr,eosvars.eosanyerr);
+                 tau = (rho0 + rho0*(*epsilon) + lg.bsq + p)*w*w - (p + lg.bsq/2.0) - w*w*lg.VdotBsq - rho0*w;
+                 u = (*epsilon) * rho0;
+                 //tau = ( (rho + rho*epsilon) + xpress ) * w_lorentz**2 - xpress - rho*w_lorentz;
+                 fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,*epsilon,eosvars.eos_temp[0],eosvars.eos_y_e[0]);
+                 eosvars.eoskeytemp =0; 
+              }
+           }
+        }
+        else {
+          CCTK_WARN(failinfomode,"EOS error in c2p 0");
+          fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars.eosit,eosvars.eosrl,eosvars.eosi,eosvars.eosj,eosvars.eosk,eosvars.eosx,eosvars.eosy,eosvars.eosz);
+          fprintf(stdout,"rho=%f, dens=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,dens,*epsilon,eosvars.eos_temp[0],eosvars.eos_y_e[0]);
+          fprintf(stdout,"code=%d\n",eosvars.eoskeyerr[0]);
+          CCTK_WARN(failwarnmode,"Aborting!!!");
+       }
+    } 
+  }
 
-  W_last = w*gammasq ;
+  lg.tau = tau;
 
-  //fprintf(stdout," p                      = %26.16e \n", p                );
+  //fprintf(stdout," p                  = %26.16e \n", p                );
 
   // Make sure that W is large enough so that v^2 < 1 : 
   i_increase = 0;
@@ -410,16 +505,51 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   /* Problem with solver, so return denoting error before doing anything further */
   if( ((*retval) != 0.) || (W == FAIL_VAL) ) {
     *retval = *retval*100.+1.;
-    fprintf(stdout," *retval      = %26.16e \n", *retval      );
-    fprintf(stdout," x_3d[%d] = %26.16e \n", 0, x_3d[0]     );
-    fprintf(stdout," x_3d[%d] = %26.16e \n", 1, x_3d[1]     );
-    fprintf(stdout," x_3d[%d] = %26.16e \n", 2, x_3d[2]     );
+  fprintf(stdout," *retval      = %26.16e \n", *retval      );
+  fprintf(stdout," *dens        = %26.16e \n", *dens_in         );
+  fprintf(stdout," *sx          = %26.16e \n", *sx_in     );    
+  fprintf(stdout," *sy          = %26.16e \n", *sy_in     );    
+  fprintf(stdout," *sz          = %26.16e \n", *sz_in     );    
+  fprintf(stdout," *tau         = %26.16e \n", *tau_in    );    
+  fprintf(stdout," *Bconsx      = %26.16e \n", *Bconsx_in );    
+  fprintf(stdout," *Bconsy      = %26.16e \n", *Bconsy_in );    
+  fprintf(stdout," *Bconsz      = %26.16e \n", *Bconsz_in );    
+  fprintf(stdout," *rho         = %26.16e \n", *rho       );    
+  fprintf(stdout," *velx        = %26.16e \n", *velx      );    
+  fprintf(stdout," *vely        = %26.16e \n", *vely      );    
+  fprintf(stdout," *velz        = %26.16e \n", *velz      );    
+  fprintf(stdout," *epsilon     = %26.16e \n", *epsilon     );
+  fprintf(stdout," *temp_in     = %26.16e \n", *temp_in     );
+  fprintf(stdout," *y_e_in     = %26.16e \n", *y_e_in     );
+  fprintf(stdout," *pressure    = %26.16e \n", *pressure    );
+  fprintf(stdout," *Bx          = %26.16e \n", *Bx        );    
+  fprintf(stdout," *By          = %26.16e \n", *By        );    
+  fprintf(stdout," *Bz          = %26.16e \n", *Bz        );    
+  fprintf(stdout," *bsq         = %26.16e \n", *bsq       );    
+  fprintf(stdout," *w_lorentz   = %26.16e \n", *w_lorentz   );
+  fprintf(stdout," *gxx         = %26.16e \n", *gxx       );    
+  fprintf(stdout," *gxy         = %26.16e \n", *gxy       );    
+  fprintf(stdout," *gxz         = %26.16e \n", *gxz       );    
+  fprintf(stdout," *gyy         = %26.16e \n", *gyy       );    
+  fprintf(stdout," *gyz         = %26.16e \n", *gyz       );    
+  fprintf(stdout," *gzz         = %26.16e \n", *gzz       );    
+  fprintf(stdout," *uxx         = %26.16e \n", *uxx       );    
+  fprintf(stdout," *uxy         = %26.16e \n", *uxy       );    
+  fprintf(stdout," *uxz         = %26.16e \n", *uxz       );    
+  fprintf(stdout," *uyy         = %26.16e \n", *uyy       );    
+  fprintf(stdout," *uyz         = %26.16e \n", *uyz       );    
+  fprintf(stdout," *uzz         = %26.16e \n", *uzz       );    
+  fprintf(stdout," *det         = %26.16e \n", *det       );    
+  fprintf(stdout," *epsnegative = %10d    \n", *epsnegative );
+  fprintf(stdout," *retval      = %26.16e \n", *retval      );
+  fflush(stdout);
     return;
   }
   else{
     if(W <= 0. || W > W_TOO_BIG) {
       *retval = 3.;
       fprintf(stdout," *retval      = %26.16e \n", *retval      );
+      fflush(stdout);
       return;
     }
   }
@@ -428,6 +558,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   if( vsq >= 1. ) {
     *retval = 4.;
     fprintf(stdout," *retval      = %26.16e \n", *retval      );
+    fflush(stdout);
     return;
   }
 
@@ -445,9 +576,80 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   } else {
     u=x_3d[2];
     *epsilon = u/rho0;
+
+    // Update tau in case we have changed it during the NR iterations
+
+    tau = lg.tau;
+
+//    eosvars.eosanyerr[0]=5;
+
     p = pressure_rho0_eps_eosomni(rho0,*epsilon,&dum1,&dum2,&eosvars) ;  // EOSOMNI
-   // printf("%g %g %g %g\n",rho0,u,*epsilon,p);
-  }
+
+//    #pragma omp critical 
+    if (eosvars.eosanyerr[0] != 0) {
+     fprintf(stdout,"inside 4: anyerr=%d\n",eosvars.eosanyerr[0]);
+       if (eosvars.eosrl >= GRHydro_c2p_warn_from_reflevel) {
+          if(eosvars.eoskeyerr[0] == 667 && GRHydro_eos_hot_eps_fix != 0) {
+             // Handling of the case when no new temperature can be
+             // found for a given epsilon. The amount of times
+             // we get to this place should be monitored, as this
+             // 'failsafe' mode creates artificial heat.
+             nf = 0;
+             while(eosvars.eosanyerr[0] != 0 && nf<= nfudgemax) {
+                eosvars.eosanyerr[0] = 0;
+                tau = tau + rho0*abs((*epsilon))*0.05*w*w;
+                *epsilon = (*epsilon) + abs((*epsilon))*0.05;
+                u = (*epsilon) * rho0;
+                EOS_Omni_press(eosvars.eoshandle,eosvars.eoskeytemp,eosvars.eosprec,1,
+                   &rho0,epsilon,eosvars.eos_temp,
+                   eosvars.eos_y_e,&p,eosvars.eoskeyerr,eosvars.eosanyerr);
+                nf = nf + 1;
+             } 
+             fprintf(stdout, "Iterations of heat injection: %d \n",nf);
+             fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars.eosit,eosvars.eosrl,eosvars.eosi,eosvars.eosj,eosvars.eosk,eosvars.eosx,eosvars.eosy,eosvars.eosz);
+             if(nf >= nfudgemax) { 
+                if(GRHydro_c2p_reset_eps_tau_hot_eos != 1) {
+                   CCTK_WARN(failinfomode,"EOS error in c2p 0: injected heat too many times");
+                   fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars.eosit,eosvars.eosrl,eosvars.eosi,eosvars.eosj,eosvars.eosk,eosvars.eosx,eosvars.eosy,eosvars.eosz);
+                   //CCTK_WARN(failinfomode,warnline);
+                   fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,*epsilon,eosvars.eos_temp[0],eosvars.eos_y_e[0]);
+                   //call CCTK_WARN(failinfomode,warnline)
+                   fprintf(stdout,"code=%d\n",eosvars.eoskeyerr[0]);
+                   //call CCTK_WARN(failinfomode,warnline)
+                   CCTK_WARN(failwarnmode,"Aborting!!!");
+                } 
+                else {
+                   CCTK_WARN(failinfomode,"EOS error in c2p 0: LAST RESORT -- reset eps and tau based on old temp");
+                   fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars.eosit,eosvars.eosrl,eosvars.eosi,eosvars.eosj,eosvars.eosk,eosvars.eosx,eosvars.eosy,eosvars.eosz);
+                   fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,*epsilon,eosvars.eos_temp[0],eosvars.eos_y_e[0]);
+                   eosvars.eoskeytemp =1; 
+                   EOS_Omni_press(eosvars.eoshandle,eosvars.eoskeytemp,eosvars.eosprec,1,
+                   &rho0,epsilon,eosvars.eos_temp,eosvars.eos_y_e,&p,eosvars.eoskeyerr,eosvars.eosanyerr);
+                 //  tau = (rho0 + rho0*(*epsilon) + lg.bsq + p)*w*w - (p + lg.bsq/2.0) - rho0*w;
+                   tau = (rho0 + rho0*(*epsilon) + lg.bsq + p)*w*w - (p + lg.bsq/2.0) - w*w*lg.VdotBsq - rho0*w;
+                   u = (*epsilon) * rho0;
+                   //tau = ( (rho + rho*epsilon) + xpress ) * w_lorentz**2 - xpress - rho*w_lorentz;
+                   //fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,*epsilon,eosvars.eos_temp[0],eosvars.eos_y_e[0]);
+                   eosvars.eoskeytemp =0; 
+                }
+             }
+          }
+          else {
+            CCTK_WARN(failinfomode,"EOS error in c2p 0");
+            fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars.eosit,eosvars.eosrl,eosvars.eosi,eosvars.eosj,eosvars.eosk,eosvars.eosx,eosvars.eosy,eosvars.eosz);
+            fprintf(stdout,"rho=%f, dens=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,dens,*epsilon,eosvars.eos_temp[0],eosvars.eos_y_e[0]);
+            fprintf(stdout,"code=%d\n",eosvars.eoskeyerr[0]);
+            CCTK_WARN(0,"Aborting!!!");
+         }
+       } 
+     }
+
+    // printf("%g %g %g %g\n",rho0,u,*epsilon,p);
+   }
+
+   // update tau in case we have modified it in C2P
+
+   *tau_in = tau*sqrt_detg;
 
   // User may want to handle this case differently, e.g. do NOT return upon 
   // a negative rho/u, calculate v^i so that rho/u can be floored by other routine:
@@ -455,11 +657,12 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   if(rho0 <= 0.) { 
 //    *epsnegative = 1; 
     fprintf(stdout," *epsnegative = %10d    \n", *epsnegative );
-    fprintf(stdout," rho0              = %26.16e \n", rho0                ); 
-    fprintf(stdout," u                      = %26.16e \n", u                );
-    fprintf(stdout," W                      = %26.16e \n", W                );
-    fprintf(stdout," vsq                      = %26.16e \n", vsq                );
-    fprintf(stdout," uold                      = %26.16e \n", uold                );
+    fprintf(stdout," rho0       = %26.16e \n", rho0             ); 
+    fprintf(stdout," u                  = %26.16e \n", u                );
+    fprintf(stdout," W                  = %26.16e \n", W                );
+    fprintf(stdout," vsq                = %26.16e \n", vsq              );
+    fprintf(stdout," uold               = %26.16e \n", uold             );
+    fflush(stdout);
     return;
   }
 
@@ -482,8 +685,48 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
     *velz = 0.0;
     *w_lorentz = 1.0;
   }
-
 #if(DEBUG_CON2PRIMM)
+  fprintf(stdout," *retvalo      = %26.16e \n", *retval      );
+  fprintf(stdout," *x2        = %26.16e \n", *x        );
+  fprintf(stdout," *y2        = %26.16e \n", *y         );
+  fprintf(stdout," *z2        = %26.16e \n", *z         );
+  fprintf(stdout," *denso        = %26.16e \n", *dens_in         );
+  fprintf(stdout," *sxo                 = %26.16e \n", *sx_in     );    
+  fprintf(stdout," *syo         = %26.16e \n", *sy_in     );    
+  fprintf(stdout," *szo                 = %26.16e \n", *sz_in     );    
+  fprintf(stdout," *tauo        = %26.16e \n", *tau_in    );    
+  fprintf(stdout," *Bconsxo     = %26.16e \n", *Bconsx_in );    
+  fprintf(stdout," *Bconsyo     = %26.16e \n", *Bconsy_in );    
+  fprintf(stdout," *Bconszo      = %26.16e \n", *Bconsz_in );    
+  fprintf(stdout," *rhoo         = %26.16e \n", *rho      );    
+  fprintf(stdout," *velxo        = %26.16e \n", *velx     );    
+  fprintf(stdout," *velyo        = %26.16e \n", *vely     );    
+  fprintf(stdout," *velzo        = %26.16e \n", *velz     );    
+  fprintf(stdout," *epsilono     = %26.16e \n", *epsilon     );
+  fprintf(stdout," *temp_ino     = %26.16e \n", *temp_in     );
+  fprintf(stdout," *y_e_ino     = %26.16e \n", *y_e_in     );
+  fprintf(stdout," *pressureo    = %26.16e \n", *pressure    );
+  fprintf(stdout," *Bxo                 = %26.16e \n", *Bx        );    
+  fprintf(stdout," *Byo                 = %26.16e \n", *By        );    
+  fprintf(stdout," *Bzo                 = %26.16e \n", *Bz        );    
+  fprintf(stdout," *bsqo                = %26.16e \n", *bsq       );    
+  fprintf(stdout," *w_lorentzo   = %26.16e \n", *w_lorentz   );
+  fprintf(stdout," *gxxo        = %26.16e \n", *gxx       );    
+  fprintf(stdout," *gxyo        = %26.16e \n", *gxy       );    
+  fprintf(stdout," *gxzo        = %26.16e \n", *gxz       );    
+  fprintf(stdout," *gyyo        = %26.16e \n", *gyy       );    
+  fprintf(stdout," *gyzo        = %26.16e \n", *gyz       );    
+  fprintf(stdout," *gzzo        = %26.16e \n", *gzz       );    
+  fprintf(stdout," *uxxo        = %26.16e \n", *uxx       );    
+  fprintf(stdout," *uxyo        = %26.16e \n", *uxy       );    
+  fprintf(stdout," *uxzo                = %26.16e \n", *uxz       );    
+  fprintf(stdout," *uyyo        = %26.16e \n", *uyy       );    
+  fprintf(stdout," *uyzo        = %26.16e \n", *uyz       );    
+  fprintf(stdout," *uzzo                = %26.16e \n", *uzz       );    
+  fprintf(stdout," *deto                = %26.16e \n", *det       );    
+  fprintf(stdout," *epsnegativeo = %10d    \n", *epsnegative );
+
+
   fprintf(stdout,"rho          = %26.16e \n",*rho      );
   fprintf(stdout,"epsilon      = %26.16e \n",*epsilon  );
   fprintf(stdout,"pressure     = %26.16e \n",*pressure );
@@ -759,13 +1002,13 @@ static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *e
     /*  before stopping                                                          */
     /*****************************************************************************/
     
-    if( (fabs(errx) <= NEWT_TOL) && (doing_extra == 0) && (EXTRA_NEWT_ITER > 0) ) {
+    if( (fabs(errx) <= eosvars->prec) && (doing_extra == 0) && (EXTRA_NEWT_ITER > 0) ) {
       doing_extra = 1;
     }
 
     if( doing_extra == 1 ) i_extra++ ;
 
-    if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) 
+    if( ((fabs(errx) <= eosvars->prec)&&(doing_extra == 0)) 
         || (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) {
       keep_iterating = 0;
     }
@@ -781,10 +1024,10 @@ static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *e
   }
 
 
-  if( fabs(errx) <= NEWT_TOL ){
+  if( fabs(errx) <= eosvars->prec ){
     return(0);
   }
-  else if( (fabs(errx) <= MIN_NEWT_TOL) && (fabs(errx) > NEWT_TOL) ){
+  else if( (fabs(errx) <= eosvars->prec) && (fabs(errx) > eosvars->prec) ){
     return(0);
   }
   else {
@@ -904,23 +1147,94 @@ static void func_vsq_eosomni(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[],
                              struct eosomnivars *eosvars, struct LocGlob *lgp)
 {
 
-  CCTK_REAL  W, vsq, u, p_tmp,epsilon,Wsq, dPdvsq, dPdW;
-  CCTK_REAL res0, LorInv, rho0, Winv, QB2Winv2, drho0_dv, res1, j11, detJ,detJinv;
+  CCTK_REAL p_tmp;
+  CCTK_REAL  W, vsq, u, epsilon,Wsq, dPdvsq, dPdW;
+  CCTK_REAL res0, Lor, LorInv, rho0, Winv, QB2Winv2, drho0_dv, res1, j11, detJ,detJinv;
   CCTK_REAL QB2Winv3, j10,B2plusW,detJ_gcf,t36,B2plusW_sq, j00,j01,j12,j20,j21,j22;
   CCTK_REAL dpress_dv,dpress_du,dpdrho,dpdeps,c00,c01,c02,c10,c11,c12,c20,c21,c22;
 
+  CCTK_INT nf, nfudgemax, failwarnmode, failinfomode;
+  DECLARE_CCTK_PARAMETERS;
+
+  nf =0;
+  nfudgemax =30;
+  failinfomode =2;
+  failwarnmode =1;
 
   W = x[0];
   vsq = x[1];
   u = x[2];
   
   Wsq = W*W;
-  
+  Lor = 1.0/sqrt(1.0-vsq);
   LorInv = sqrt(1.0-vsq);
   rho0 = lgp->D * LorInv;
   epsilon = u/rho0;
   p_tmp = pressure_rho0_eps_eosomni(rho0,epsilon,&dpdrho,&dpdeps,eosvars);
+  
+//  fprintf(stdout,"vsq: keyerr=%d, eosrl = %d, anyerr = %d \n",eosvars->eoskeyerr[0],eosvars->eosrl, eosvars->eosanyerr[0]);
+//  fprintf(stdout,"vsq: rho0=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,epsilon,eosvars->eos_temp[0],eosvars->eos_y_e[0]);
+  #pragma omp critical 
+  if (eosvars->eosanyerr[0] != 0) {
+     fprintf(stdout,"inside vsq_omni: anyerr=%d rl=%d\n",eosvars->eosanyerr[0], eosvars->eosrl);
+     if (eosvars->eosrl >= GRHydro_c2p_warn_from_reflevel) {
+        if(eosvars->eoskeyerr[0] == 667 && GRHydro_eos_hot_eps_fix != 0) {
+           // Handling of the case when no new temperature can be
+           // found for a given epsilon. The amount of times
+           // we get to this place should be monitored, as this
+           // 'failsafe' mode creates artificial heat.
+           nf = 0;
+           while(eosvars->eosanyerr[0] != 0 && nf<= nfudgemax) {
+              eosvars->eosanyerr[0] = 0;
+              lgp->tau = lgp->tau + rho0*abs(epsilon)*0.05*Lor*Lor;
+              epsilon = epsilon + abs(epsilon)*0.05;
+              u = epsilon*rho0;
+              CCTK_WARN(failwarnmode,"Injecting heat inside vsq_omni without resetting tau");
+              fprintf(stdout,"Injecting heat inside vsq_omni without resetting tau: anyerr=%d\n",eosvars->eosanyerr[0]);
+              EOS_Omni_press(eosvars->eoshandle,eosvars->eoskeytemp,eosvars->eosprec,1,
+                 &rho0,&epsilon,eosvars->eos_temp,
+                 eosvars->eos_y_e,&p_tmp,eosvars->eoskeyerr,eosvars->eosanyerr);
+              nf = nf + 1;
+           } 
+           fprintf(stdout, "Iterations of heat injection: %d \n",nf);
+           fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars->eosit,eosvars->eosrl,eosvars->eosi,eosvars->eosj,eosvars->eosk,eosvars->eosx,eosvars->eosy,eosvars->eosz);
+           if(nf >= nfudgemax) { 
+              if(GRHydro_c2p_reset_eps_tau_hot_eos != 1) {
+                 CCTK_WARN(failinfomode,"EOS error in c2p 0: injected heat too many times");
+                 fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars->eosit,eosvars->eosrl,eosvars->eosi,eosvars->eosj,eosvars->eosk,eosvars->eosx,eosvars->eosy,eosvars->eosz);
+                 //CCTK_WARN(failinfomode,warnline);
+                 fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,epsilon,eosvars->eos_temp[0],eosvars->eos_y_e[0]);
+                 //call CCTK_WARN(failinfomode,warnline)
+                 fprintf(stdout,"code=%d\n",eosvars->eoskeyerr[0]);
+                 //call CCTK_WARN(failinfomode,warnline)
+                 CCTK_WARN(failwarnmode,"Aborting!!!");
+              } 
+              else {
+                 CCTK_WARN(failinfomode,"EOS error in c2p 0: LAST RESORT -- reset eps and tau based on old temp");
+                 fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars->eosit,eosvars->eosrl,eosvars->eosi,eosvars->eosj,eosvars->eosk,eosvars->eosx,eosvars->eosy,eosvars->eosz);
+                 fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,epsilon,eosvars->eos_temp[0],eosvars->eos_y_e[0]);
+                 eosvars->eoskeytemp =1; 
+                 EOS_Omni_press(eosvars->eoshandle,eosvars->eoskeytemp,eosvars->eosprec,1,
+                 &rho0,&epsilon,eosvars->eos_temp,eosvars->eos_y_e,&p_tmp,eosvars->eoskeyerr,eosvars->eosanyerr);
+                 //tau = (rho0 + rho0*(*epsilon) + lg.bsq + p)*w*w - (p + lg.bsq/2.0) - w*w*lg.VdotBsq - rho0*w;
+                 lgp->tau = (rho0 + rho0*epsilon + lgp->bsq + p_tmp)*Lor*Lor - (p_tmp + lgp->bsq/2.0) - Lor*Lor*lgp->VdotBsq - rho0*Lor;
+                 u = epsilon*rho0;
+                 eosvars->eoskeytemp = 0;
+                 fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,epsilon,eosvars->eos_temp[0],eosvars->eos_y_e[0]);
+              }
+           }
+        }
+        else {
+          CCTK_WARN(failinfomode,"EOS error in c2p 0");
+          fprintf(stdout, "it=%d, rl=%d, i=%d, j=%d, k=%d, x=%f,y=%f,z=%f\n",eosvars->eosit,eosvars->eosrl,eosvars->eosi,eosvars->eosj,eosvars->eosk,eosvars->eosx,eosvars->eosy,eosvars->eosz);
+          fprintf(stdout,"rho=%f, epsilon=%f, temp=%f, ye=%f\n", rho0,epsilon,eosvars->eos_temp[0],eosvars->eos_y_e[0]);
+          fprintf(stdout,"code=%d\n",eosvars->eoskeyerr[0]);
+          CCTK_WARN(failwarnmode,"Aborting!!!");
+       }
+    } 
+  }
 
+  lgp->Qdotn = -(lgp->tau + lgp->D);
   B2plusW = lgp->Bsq+W;
   B2plusW_sq = B2plusW*B2plusW;
   Winv = 1/W;