Finish the implementation of con2prim using the entropy equation.

  * Tested successfully pointwisely. It still needs to be tested on evolution.

From: Bruno Coutinho Mundim <bcmsma@astro.rit.edu>

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

4 files changed, 649 insertions, 10 deletions

diff --git a/interface.ccl b/interface.ccl
index 89cae27..db6615d 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -84,6 +84,27 @@ void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, CCTK_R
                            CCTK_INT OUT epsnegative, \
                            CCTK_REAL OUT retval)
 
+void FUNCTION Con2PrimGenMee(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, \
+     CCTK_REAL INOUT prec, CCTK_REAL INOUT gamma_eos, \
+     CCTK_REAL INOUT dens, \
+     CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \
+     CCTK_REAL INOUT tau, \
+     CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \
+     CCTK_REAL INOUT entropycons, \
+     CCTK_REAL INOUT y_e, CCTK_REAL INOUT temp, CCTK_REAL INOUT rho, \
+     CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \
+     CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \
+     CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \
+     CCTK_REAL INOUT Bvecsq, \
+     CCTK_REAL INOUT w_lorentz, \
+     CCTK_REAL INOUT gxx, CCTK_REAL INOUT gxy, CCTK_REAL INOUT gxz, \
+     CCTK_REAL INOUT gyy, CCTK_REAL INOUT gyz, CCTK_REAL INOUT gzz, \
+     CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \
+     CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \
+     CCTK_REAL INOUT det, \
+     CCTK_INT OUT epsnegative, \
+     CCTK_REAL OUT retval)
+
 #void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, CCTK_REAL INOUT prec,CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \
 #                           CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \
 #                           CCTK_REAL INOUT tau, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \
@@ -207,6 +228,7 @@ PROVIDES FUNCTION UpperMet WITH UpperMetric LANGUAGE Fortran
 #PROVIDES FUNCTION Con2Prim WITH Con2Prim_pt LANGUAGE Fortran
 PROVIDES FUNCTION Con2PrimPoly WITH Con2Prim_ptPolytype LANGUAGE Fortran
 PROVIDES FUNCTION Con2PrimGenM WITH GRHydro_Con2PrimM_pt LANGUAGE Fortran
+PROVIDES FUNCTION Con2PrimGenMee WITH GRHydro_Con2PrimM_ptee LANGUAGE Fortran
 PROVIDES FUNCTION Con2PrimGen WITH Con2Prim_pt LANGUAGE Fortran
 PROVIDES FUNCTION Con2PrimPolyM WITH GRHydro_Con2PrimM_Polytype_pt LANGUAGE Fortran
 PROVIDES FUNCTION Prim2ConGen WITH prim2con LANGUAGE Fortran
diff --git a/src/GRHydro_Con2PrimM.F90 b/src/GRHydro_Con2PrimM.F90
index d59cc9c..0960b0d 100644
--- a/src/GRHydro_Con2PrimM.F90
+++ b/src/GRHydro_Con2PrimM.F90
@@ -334,7 +334,11 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
               Bvecz_tmp = Bprim(i,j,k,3)
 
               keytemp = 0
-
+              !Watch out for the values returned to b2. Here b2 is the Bprim^2
+              !while inside the point-wise con2prim routines it is the square
+              !of the comoving B-field, b^{\mu} b_{\mu}. It is overwritten 
+              !in this routine, but we may need to find a better notation 
+              !avoid future confusions.
               call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, &
                    GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), &
                    scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), &
@@ -347,7 +351,13 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
                    epsnegative,GRHydro_C2P_failed(i,j,k))
 
               if(evolve_entropy.ne.0) then
-                 entropy(i,j,k) = entropycons(i,j,k)/dens(i,j,k)*rho(i,j,k)
+                if(GRHydro_C2P_failed(i,j,k).ne.0) then
+                  !Use previous time step for rho:
+                  entropy(i,j,k) = entropycons(i,j,k)/dens(i,j,k)*rho(i,j,k)
+                else
+                  !Use the current correct value of rho returned by con2prim:
+                  entropy(i,j,k) = entropycons(i,j,k)/dens(i,j,k)*rho_tmp
+                endif
               endif
 
               if(GRHydro_C2P_failed(i,j,k).ne.0) then
@@ -430,14 +440,28 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
                   Bvecy_tmp = Bprim(i,j,k,2)
                   Bvecz_tmp = Bprim(i,j,k,3)
 
-                  call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, &
-                       dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, &
-                       Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3),rho_tmp,&
-                       velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,&
-                       Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,&
-                       g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), &
-                       uxx,uxy,uxz,uyy,uyz,uzz,det, &
-                       epsnegative,GRHydro_C2P_failed(i,j,k))
+                  if(evolve_entropy.ne.0) then
+                    call GRHydro_Con2PrimM_ptee(GRHydro_eos_handle, keytemp, &
+                         GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), &
+                         scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), &
+                         Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3), &
+                         entropycons(i,j,k), xye(1), &
+                         xtemp(1),rho_tmp,velx_tmp,vely_tmp,velz_tmp,&
+                         eps_tmp,press_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,&
+                         w_lorentz_tmp,g11(i,j,k),g12(i,j,k),g13(i,j,k),&
+                         g22(i,j,k),g23(i,j,k),g33(i,j,k), &
+                         uxx,uxy,uxz,uyy,uyz,uzz,det, &
+                         epsnegative,GRHydro_C2P_failed(i,j,k))
+                  else
+                    call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, &
+                         dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, &
+                         Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3),rho_tmp,&
+                         velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,&
+                         Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,&
+                         g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), &
+                         uxx,uxy,uxz,uyy,uyz,uzz,det, &
+                         epsnegative,GRHydro_C2P_failed(i,j,k))
+                  end if 
 
                   rho(i,j,k) = rho_tmp 
                   press(i,j,k) = press_tmp 
diff --git a/src/GRHydro_Con2PrimM_ptee.c b/src/GRHydro_Con2PrimM_ptee.c
new file mode 100644
index 0000000..de58d4c
--- /dev/null
+++ b/src/GRHydro_Con2PrimM_ptee.c
@@ -0,0 +1,592 @@
+/***********************************************************************************
+    Copyright 2006 Scott C. Noble, Charles F. Gammie, Jonathan C. McKinney, 
+                   and Luca Del Zanna.
+
+                        PVS_GRMHD
+
+    This file was derived from PVS_GRMHD.  The authors of PVS_GRMHD include 
+    Scott C. Noble, Charles F. Gammie, Jonathan C. McKinney, and Luca Del Zanna.
+    PVS_GRMHD is available under the GPL from:
+    http://rainman.astro.uiuc.edu/codelib/  
+
+    You are morally obligated to cite the following  paper in his/her 
+    scientific literature that results from use of this file:
+
+    [1] Noble, S. C., Gammie, C. F., McKinney, J. C., \& Del Zanna, L. \ 2006, 
+        Astrophysical Journal, 641, 626.
+
+    PVS_GRMHD is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    PVS_GRMHD is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with PVS_GRMHD; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+
+    If the user has any questions, please direct them to Scott C. Noble at 
+    scn@astro.rit.edu  . 
+
+***********************************************************************************/
+
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdarg.h>
+#include <string.h>
+#include <math.h>
+#include <float.h>
+#include <complex.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+/* Set this to be 1 if you want debug output */
+#define DEBUG_CON2PRIMM (0)
+
+
+/* Adiabatic index used for the state equation */
+
+#define MAX_NEWT_ITER   (30)       /* Max. # of Newton-Raphson iterations for find_root_2D(); */
+#define NEWT_TOL        (1.0e-10)  /* Min. of tolerance allowed for Newton-Raphson iterations */
+#define MIN_NEWT_TOL    (1.0e-10)  /* Max. of tolerance allowed for Newton-Raphson iterations */
+#define EXTRA_NEWT_ITER (2)
+
+#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
+                                      to always be smaller than this.  This
+				      is used to detect solver failures */
+
+#define FAIL_VAL        (1.e30)    /* Generic value to which we set variables when a problem arises */
+
+/**************************************************
+  The following functions assume a Gamma-law EOS:
+***************************************************/
+
+/* Local Globals */
+struct LocGlob {
+  CCTK_REAL Bsq, QdotBsq, Qtsq, Qdotn, D, half_Bsq, Sc, g_o_gm1,
+           W_for_gnr2, rho_for_gnr2, W_for_gnr2_old, rho_for_gnr2_old, drho_dW ;
+} ;
+
+// Declarations: 
+
+
+void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptee)  ( 
+    CCTK_INT  *handle, CCTK_INT *keytemp, CCTK_REAL *prec,
+    CCTK_REAL *gamma_eos,
+    CCTK_REAL *dens_in, 
+    CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, 
+    CCTK_REAL *tau_in,     
+    CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_in, 
+    CCTK_REAL *entropycons_in,
+    CCTK_REAL *y_e_in, CCTK_REAL* temp_in,
+    CCTK_REAL *rho, 
+    CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz,
+    CCTK_REAL *epsilon, CCTK_REAL *pressure,
+    CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, 
+    CCTK_REAL *bsq,
+    CCTK_REAL *w_lorentz, 
+    CCTK_REAL *gxx, CCTK_REAL *gxy, CCTK_REAL *gxz, 
+    CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, 
+    CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz,
+    CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz,
+    CCTK_REAL *det,
+    CCTK_INT  *epsnegative,
+    CCTK_REAL *retval) ;
+
+static CCTK_INT general_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos,
+               struct LocGlob *lgp,
+               void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [],
+                    CCTK_REAL [][1], CCTK_REAL *,
+                    CCTK_REAL *, CCTK_REAL, struct LocGlob *) );
+
+static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[],
+          CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, 
+          CCTK_REAL gammaeos, struct LocGlob *lgp);
+
+/**********************************************************************/
+/**********************************************************************************
+
+  Con2PrimM_pt():
+  -----------------------------
+
+   -- Attempts an inversion from GRMHD conserved variables to primitive variables assuming a guess.
+
+   -- Uses the 2D method of Noble et al. (2006): 
+       -- Solves for two independent variables (W,v^2) via a 2D
+          Newton-Raphson method 
+       -- Can be used (in principle) with a general equation of state. 
+
+   -- Minimizes two residual functions using a homemade Newton-Raphson routine. 
+       -- It is homemade so that it can catch exceptions and handle them correctly, plus it is 
+          optimized for this particular  problem. 
+
+  -- Note that the notation used herein is that of Noble et al. (2006) except for the argument 
+     list. 
+
+
+INPUT:  (using GRHydro variable defintions)
+
+   s[x,y,z]   =  scons[0,1,2]  = \alpha \sqrt(\gamma) T^0_i 
+   dens, tau  =  as defined in GRHydro and are assumed to be densitized (i.e. with sqrt(\gamma))   
+   dens       =  D = \sqrt(\gamma) W \rho    
+   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        =  sqrt(\gamma)
+   B[x,y,z]   =  Bvec[0,1,2] 
+   bsq        = b^\mu b_\mu  
+
+   epsnegative = (integer)
+               = 0  if rho and epsilon are positive
+              != 0  otherwise 
+
+
+  --  (users should set B[x,y,z] = 0  for hydrodynamic runs) 
+
+
+OUTPUT:  (using GRHydro variable defintions)
+   rho, eps   =  as defined in GRHydro, primitive variables
+   vel[x,y,z] =  as defined in GRHydro, primitive variables
+
+
+RETURN VALUE: of retval = (i*100 + j)  where 
+         i = 0 ->  Newton-Raphson solver either was not called (yet or not used) 
+                   or returned successfully;
+             1 ->  Newton-Raphson solver did not converge to a solution with the 
+                   given tolerances;
+             2 ->  Newton-Raphson procedure encountered a numerical divergence 
+                   (occurrence of "nan" or "+/-inf" ;
+	     
+         j = 0 -> success 
+             1 -> failure: some sort of failure in Newton-Raphson; 
+             2 -> failure: unphysical vsq = v^2  value at initial guess;
+	     3 -> failure: W<0 or W>W_TOO_BIG
+             4 -> failure: v^2 > 1 
+             ( used to be  5 -> failure: rho,uu <= 0   but now sets epsnegative to non-zero )
+
+**********************************************************************************/
+
+void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptee)  ( 
+    CCTK_INT  *handle, CCTK_INT *keytemp, CCTK_REAL *prec,
+    CCTK_REAL *gamma_eos,
+    CCTK_REAL *dens_in, 
+    CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, 
+    CCTK_REAL *tau_in,     
+    CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_in, 
+    CCTK_REAL *entropycons_in,
+    CCTK_REAL *y_e_in, CCTK_REAL* temp_in,
+    CCTK_REAL *rho, 
+    CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz,
+    CCTK_REAL *epsilon, CCTK_REAL *pressure,
+    CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, 
+    CCTK_REAL *bsq,
+    CCTK_REAL *w_lorentz, 
+    CCTK_REAL *gxx, CCTK_REAL *gxy, CCTK_REAL *gxz, 
+    CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, 
+    CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz,
+    CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz,
+    CCTK_REAL *det,
+    CCTK_INT  *epsnegative,
+    CCTK_REAL *retval)
+
+{
+  CCTK_REAL x_1d[1];
+  CCTK_REAL sx, sy, sz;
+  CCTK_REAL usx, usy, usz;
+  CCTK_REAL tau, dens, gammaeos;
+  CCTK_REAL QdotB;
+  CCTK_REAL rho0,u,p,w,gammasq,gamma,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_REAL rho_gm1;
+  CCTK_REAL utsq;
+ 
+  DECLARE_CCTK_PARAMETERS;
+ 
+  struct LocGlob lg; 
+
+  gammaeos = *gamma_eos;
+
+  /* Assume ok initially: */
+  // BCM: But let the driver function take care of its initialization
+  //*retval = 0.; 
+  *epsnegative = 0; 
+
+#if(DEBUG_CON2PRIMM)
+  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," *entropycons      = %26.16e \n", *entropycons_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," *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);
+#endif
+
+  /* First undensitize all conserved variables : */
+  sx   = (  *sx_in)   * inv_sqrt_detg;
+  sy   = (  *sy_in)   * inv_sqrt_detg;
+  sz   = (  *sz_in)   * inv_sqrt_detg;
+  tau  = ( *tau_in)   * inv_sqrt_detg;
+  dens = (*dens_in)   * inv_sqrt_detg;
+
+  usx  =  (*uxx)*sx + (*uxy)*sy + (*uxz)*sz;
+  usy  =  (*uxy)*sx + (*uyy)*sy + (*uyz)*sz;
+  usz  =  (*uxz)*sx + (*uyz)*sy + (*uzz)*sz;
+
+  *Bx = (*Bconsx_in) * inv_sqrt_detg;
+  *By = (*Bconsy_in) * inv_sqrt_detg;
+  *Bz = (*Bconsz_in) * inv_sqrt_detg;
+
+  lg.Sc = (*entropycons_in) * inv_sqrt_detg;
+  // Calculate various scalars (Q.B, Q^2, etc)  from the conserved variables:
+
+  lg.g_o_gm1 = gammaeos/(gammaeos-1.0);
+  
+  lg.Bsq = 
+    (*gxx) * (*Bx) * (*Bx) + 
+    (*gyy) * (*By) * (*By) + 
+    (*gzz) * (*Bz) * (*Bz) + 
+    2*( 
+       (*gxy) * (*Bx) * (*By) +
+       (*gxz) * (*Bx) * (*Bz) +
+       (*gyz) * (*By) * (*Bz) );
+    
+  QdotB = (sx * (*Bx) + sy * (*By) + sz * (*Bz)) ;
+  lg.QdotBsq = QdotB*QdotB ;
+
+  lg.Qdotn = -(tau + dens) ;
+
+  lg.Qtsq = (usx * sx  +  usy * sy  +  usz * sz) ;
+
+  lg.D = dens;
+
+  lg.half_Bsq = 0.5*lg.Bsq;
+
+  /* calculate W from last timestep and use for guess */
+  vsq = 
+    (*gxx) * (*velx) * (*velx) + 
+    (*gyy) * (*vely) * (*vely) + 
+    (*gzz) * (*velz) * (*velz) + 
+    2*( 
+       (*gxy) * (*velx) * (*vely) +
+       (*gxz) * (*velx) * (*velz) +
+       (*gyz) * (*vely) * (*velz) );
+
+  if( (vsq < 0.) && (fabs(vsq) < 1.0e-13) ) { 
+    vsq = fabs(vsq);
+  }
+  if(vsq < 0. || vsq > 1. ) {
+    *retval = 2.;
+    return;
+  }
+
+  gammasq = 1. / (1. - vsq);
+  gamma  = sqrt(gammasq);
+	
+  // Always calculate rho from D and gamma so that using D in EOS remains consistent
+  //   i.e. you don't get positive values for dP/d(vsq) . 
+  rho0 = lg.D / gamma ;
+  rho_gm1 = pow(rho0,(gammaeos-1.));
+  p = lg.Sc * rho_gm1 / gamma;
+  u = p / (gammaeos-1.);
+  w = rho0 + u + p ;
+
+//  W_last = w*gammasq ;
+
+  // Calculate W and vsq: 
+  x_1d[0] = rho0;
+//  *retval = 1.0*twod_newton_raphson( x_2d, gammaeos, &lg, func_vsq ) ;  
+  *retval = general_newton_raphson( x_1d, gammaeos, &lg, func_rho ) ;  
+  rho0 = x_1d[0];
+	
+  /* Problem with solver, so return denoting error before doing anything further */
+  if( ((*retval) != 0.) || (rho0 == FAIL_VAL) ) {
+    *retval = *retval*100.+1.;
+    return;
+  }
+  else{
+    if( rho0 > W_TOO_BIG) {
+      *retval = 3.;
+      return;
+    }
+  }
+
+  // Calculate v^2:
+  utsq = (lg.D-rho0)*(lg.D+rho0)/(rho0*rho0);
+  gammasq = 1.+utsq;
+  gamma = sqrt(gammasq);
+
+  if( utsq < 0. ) {
+    *retval = 4.;
+    return;
+  }
+
+
+  // Recover the primitive variables from the scalars and conserved variables:
+  rho0 = lg.D / gamma;
+  rho_gm1 = pow(rho0,(gammaeos-1.));
+  p = lg.Sc * rho_gm1 / gamma;
+  u = p / (gammaeos-1.);
+  w = rho0 + u + p ;
+  W = w * gammasq;
+
+  // 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:
+  if( (rho0 <= 0.) || (u <= 0.) ) { 
+    *epsnegative = 1; 
+    return;
+  }
+
+  *rho = rho0;
+  *epsilon = u / rho0;
+  *w_lorentz = gamma; 
+  *pressure = p ; 
+
+  g_o_WBsq = 1./(W+lg.Bsq);
+  QdB_o_W  = QdotB / W; 
+  *bsq = lg.Bsq / gammasq + QdB_o_W*QdB_o_W;
+
+  *velx = g_o_WBsq * ( usx + QdB_o_W*(*Bx) ) ;
+  *vely = g_o_WBsq * ( usy + QdB_o_W*(*By) ) ;
+  *velz = g_o_WBsq * ( usz + QdB_o_W*(*Bz) ) ;
+
+  if (*rho <= rho_abs_min*(1.0+GRHydro_atmo_tolerance) ) {
+    *rho = rho_abs_min;
+    *velx = 0.0;
+    *vely = 0.0;
+    *velz = 0.0;
+    *w_lorentz = 1.0;
+  }
+
+
+#if(DEBUG_CON2PRIMM)
+  fprintf(stdout,"rho          = %26.16e \n",*rho      );
+  fprintf(stdout,"epsilon      = %26.16e \n",*epsilon  );
+  fprintf(stdout,"pressure     = %26.16e \n",*pressure );
+  fprintf(stdout,"w_lorentz    = %26.16e \n",*w_lorentz);
+  fprintf(stdout,"bsq          = %26.16e \n",*bsq      );
+  fprintf(stdout,"velx         = %26.16e \n",*velx     );
+  fprintf(stdout,"vely         = %26.16e \n",*vely     );
+  fprintf(stdout,"velz         = %26.16e \n",*velz     );
+  fprintf(stdout,"gammaeos     = %26.16e \n",gammaeos  );
+  fflush(stdout);
+#endif
+
+  /* done! */
+  return;
+
+}
+
+/**********************************************************************/
+/************************************************************
+
+  general_newton_raphson(): 
+
+    -- performs Newton-Rapshon method on an arbitrary system.
+
+    -- inspired in part by Num. Rec.'s routine newt();
+
+    Arguements: 
+
+       -- x[]   = set of independent variables to solve for;
+       -- n     = number of independent variables and residuals;
+       -- funcd = name of function that calculates residuals, etc.;
+
+*****************************************************************/
+static CCTK_INT general_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos,
+               struct LocGlob *lgp,
+               void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [],
+                    CCTK_REAL [][1], CCTK_REAL *,
+                    CCTK_REAL *, CCTK_REAL, struct LocGlob *) )
+{
+  CCTK_REAL f, df, dx[1], x_old[1], resid[1],
+    jac[1][1];
+  CCTK_REAL errx, x_orig[1];
+  CCTK_INT  n_iter, i_extra, doing_extra;
+  CCTK_REAL W,W_old;
+
+  int   keep_iterating;
+
+
+  // Initialize various parameters and variables:
+  errx = 1. ;
+  df = f = 1.;
+  i_extra = doing_extra = 0;
+  x_old[0] = x_orig[0] = x[0] ;
+
+  W = W_old = 0.;
+
+  n_iter = 0;
+
+
+  /* Start the Newton-Raphson iterations : */
+  keep_iterating = 1;
+  while( keep_iterating ) {
+
+    (*funcd) (x, dx, resid, jac, &f, &df, gammaeos, lgp);  /* returns with new dx, f, df */
+
+    /* Save old values before calculating the new: */
+    errx = 0.;
+    x_old[0] = x[0] ;
+    /* don't use line search : */
+    x[0] += dx[0]  ;
+
+    /****************************************/
+    /* Calculate the convergence criterion */
+    /****************************************/
+
+    /* For the new criterion, always look at error in "W" : */
+    // METHOD specific:
+    errx  = (x[0]==0.) ?  fabs(dx[0]) : fabs(dx[0]/x[0]);
+
+
+    /****************************************/
+    /* Make sure that the new x[] is physical : */
+    /****************************************/
+    x[0] = fabs(x[0]);
+
+
+    /*****************************************************************************/
+    /* If we've reached the tolerance level, then just do a few extra iterations */
+    /*  before stopping                                                          */
+    /*****************************************************************************/
+
+    if( (fabs(errx) <= NEWT_TOL) && (doing_extra == 0) && (EXTRA_NEWT_ITER > 0) ) {
+      doing_extra = 1;
+    }
+
+    if( doing_extra == 1 ) i_extra++ ;
+
+    if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) ||
+   (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) {
+      keep_iterating = 0;
+    }
+
+    n_iter++;
+
+  }   // END of while(keep_iterating)
+
+
+  /*  Check for bad untrapped divergences : */
+  if( (!finite(f)) || (!finite(df)) || (!finite(x[0]))  ) {
+    return(2);
+  }
+
+
+  if( fabs(errx) > MIN_NEWT_TOL){
+    return(1);
+  }
+  if( (fabs(errx) <= MIN_NEWT_TOL) && (fabs(errx) > NEWT_TOL) ){
+    return(0);
+  }
+  if( fabs(errx) <= NEWT_TOL ){
+    return(0);
+  }
+
+  return(0);
+
+}
+
+/**********************************************************************/
+/*********************************************************************************
+   func_rho():
+
+        -- residual/jacobian routine to calculate rho Qtsq equation with 
+            the definition of W 
+        W  =  ( 1 + GAMMA * K_atm * rho^(GAMMA-1)/(GAMMA-1) ) D^2 / rho
+              substituted in. 
+
+     Arguments:
+          x   = current value of independent var's (on input & output);
+         dx   = Newton-Raphson step (on output);
+        resid = residuals based on x (on output);
+         jac  = Jacobian matrix based on x (on output);
+         f    =  resid.resid/2  (on output)
+        df    = -2*f;  (on output)
+         n    = dimension of x[];
+ *********************************************************************************/
+static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[],
+          CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, 
+          CCTK_REAL gammaeos, struct LocGlob *lgp)
+{
+
+  CCTK_REAL t1,t2;
+  CCTK_REAL t12;
+  CCTK_REAL t17;
+  CCTK_REAL t3, t100,rhosq,t200,rho,W,dWdrho,dvsqdrho,vsq;
+
+  rho = x[0];
+  rhosq = rho*rho;
+  t200 = 1./(lgp->D*lgp->D);
+  t100 = lgp->g_o_gm1*lgp->Sc*pow(rho,(gammaeos-1.));
+  W = lgp->D*( lgp->D + t100 ) / rho ;
+  dWdrho = lgp->D * ( -lgp->D + t100*(gammaeos-2.) ) / rhosq;
+  t1 = W*W;
+  t2 = lgp->Bsq+W;
+  //    t3 = pow(Bsq+W,2.0);
+  t3 = t2*t2;
+  vsq = (lgp->D-rho)*(lgp->D+rho)*t200;
+  dvsqdrho = -2*rho*t200;
+  resid[0] = t1*(lgp->Qtsq-vsq*t3)+lgp->QdotBsq*(t2+W);
+  t12 = lgp->Bsq*lgp->Bsq;
+  t17 = dWdrho*vsq;
+  jac[0][0] = 2*lgp->QdotBsq*dWdrho
+    +((lgp->Qtsq-vsq*t12)*2*dWdrho+(-6*t17*lgp->Bsq-dvsqdrho*t12
+                +(-2*dvsqdrho*lgp->Bsq-4*t17-dvsqdrho*W)*W)*W)*W;
+
+  dx[0] = -resid[0]/jac[0][0];
+  *f = 0.5*resid[0]*resid[0];
+  *df = -2. * (*f);
+
+  return;
+
+}
+
+
+/****************************************************************************** 
+             END   
+ ******************************************************************************/
+
+
+#undef DEBUG_CON2PRIMM
diff --git a/src/make.code.defn b/src/make.code.defn
index c16764e..7bdd9f6 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -46,6 +46,7 @@ SRCS = 	Utils.F90 \
 	GRHydro_Tmunu.F90 \
 	GRHydro_Con2PrimM.F90 \
 	GRHydro_Con2PrimM_pt.c \
+	GRHydro_Con2PrimM_ptee.c \
 	GRHydro_Con2PrimM_polytype_pt.c \
 	GRHydro_EigenproblemM.F90 \
 	GRHydro_FluxM.F90 \