* improve atmosphere handling so that it is now possible to use

  hot EOS with an atmosphere



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

1 files changed, 27 insertions, 12 deletions

diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90
index 9a3d2bc..9d69ebd 100644
--- a/src/GRHydro_Con2Prim.F90
+++ b/src/GRHydro_Con2Prim.F90
@@ -88,7 +88,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
        GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
 
   !$OMP PARALLEL DO PRIVATE(i,j,itracer,&
-  !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative)
+  !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, anyerr, keyerr, keytemp)
   do k = 1, nz 
     do j = 1, ny 
       do i = 1, nx
@@ -128,7 +128,6 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
            Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
         endif
 
-
          if ( dens(i,j,k) .le. sqrt(det)*GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then
 
            dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance)
@@ -137,15 +136,26 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
            vel(i,j,k,:) = 0.d0
            w_lorentz(i,j,k) = 1.d0
 
-           call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
-                rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
-
-           call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
-                rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
+           if(evolve_temper.ne.0) then
+              ! set hot atmosphere values
+              temperature(i,j,k) = grhydro_hot_atmo_temp
+              y_e(i,j,k) = grhydro_hot_atmo_Y_e
+              y_e_con(i,j,k) = y_e(i,j,k) * dens(i,j,k)
+              keytemp = 1
+              call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),&
+                   press(i,j,k),keyerr,anyerr)
+           else
+              ! use polytropic EOS
+              call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
+
+              call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
+           endif
 
-           ! w_lorentz=1, so the expression for tau reduces to:
            tau(i,j,k)  = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k)
-
+           
            cycle
 
          end if
@@ -196,6 +206,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
             endif
          endif
 
+
         if (epsnegative) then  
 
 #if 0
@@ -848,10 +859,14 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
     rho = GRHydro_rho_min
     udens = rho
     dens = sqrt(det) * rho
-!    epsilon = (sqrt( (utau + pnew + udens)**2) - pnew -  udens)/udens
-    epsilon = epsminl
+    temp = GRHydro_hot_atmo_temp
+    ye = GRHydro_hot_atmo_Y_e
+    keytemp=1
+    call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n, &
+         rho,epsilon,temp,ye,xpress,keyerr,anyerr)
+    keytemp=0
     ! w_lorentz=1, so the expression for utau reduces to:
-    utau  = rho + rho*epsminl - udens
+    utau  = rho + rho*epsilon - udens
     sx = 0.d0
     sy = 0.d0
     sz = 0.d0