merge branch hot_and_MHD_temp_dev into branch at revision r185

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

1 files changed, 42 insertions, 48 deletions

diff --git a/src/GRHydro_Con2PrimM.F90 b/src/GRHydro_Con2PrimM.F90
index 401fa4d..8b595d9 100644
--- a/src/GRHydro_Con2PrimM.F90
+++ b/src/GRHydro_Con2PrimM.F90
@@ -49,13 +49,6 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
   
-#if !USE_EOS_OMNI
-#ifdef _EOS_BASE_INC_
-#undef _EOS_BASE_INC_
-#endif
-#include "EOS_Base.inc"
-#endif
-
   integer :: i, j, k, itracer, nx, ny, nz
   CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, pmin, epsmin
   CCTK_REAL :: b2
@@ -67,7 +60,6 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
 
   CCTK_REAL :: local_min_tracer
 
-#if USE_EOS_OMNI
 ! begin EOS Omni vars
   integer :: n = 1
   integer :: keytemp = 0
@@ -78,7 +70,6 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
   real*8  :: xtemp = 0.0d0
   real*8  :: xye = 0.0d0
 ! end EOS Omni vars
-#endif
 
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
   call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
@@ -94,16 +85,11 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
      local_min_tracer = 0d0
   end if
 
-#if USE_EOS_OMNI
   ! this is a poly call
   call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
          GRHydro_rho_min,xeps,xtemp,xye,pmin,keyerr,anyerr)
   call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
        GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
-#else           
-  pmin   = EOS_Pressure(GRHydro_polytrope_handle, GRHydro_rho_min, 1.0d0)
-  epsmin = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, GRHydro_rho_min, pmin)
-#endif
 
   !$OMP PARALLEL DO PRIVATE(i,j,itracer,&
   !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative)
@@ -138,6 +124,11 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
               
            endif
            
+           if(evolve_Y_e.ne.0) then
+              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
               
               b2=gxx(i,j,k)*Bvec(i,j,k,1)**2+gyy(i,j,k)*Bvec(i,j,k,2)**2+gzz(i,j,k)*Bvec(i,j,k,3)**2+ &
@@ -149,16 +140,13 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
               scon(i,j,k,:) = 0.d0
               vel(i,j,k,:) = 0.d0
               w_lorentz(i,j,k) = 1.d0
-#if USE_EOS_OMNI
+
               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)
-#else
-              press(i,j,k) = EOS_Pressure(GRHydro_polytrope_handle, rho(i,j,k), eps(i,j,k))
-              eps(i,j,k) = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, rho(i,j,k), press(i,j,k))
-#endif
+
               ! w_lorentz=1, so the expression for tau reduces to:
               
 !!$ tau does need to take into account the existing B-field
@@ -169,16 +157,19 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
               
            end if
            
+           if(evolve_temper.eq.0) then
+              call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, dens(i,j,k), &
+                   scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), &
+                   rho(i,j,k),&
+                   vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),&
+                   gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
+                   uxx,uxy,uxz,uyy,uyz,uzz,det, &
+                   Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,&
+                   epsnegative,GRHydro_C2P_failed(i,j,k))
+           else
+              stop "Please implement finite T con2prim routine in MHD part!"
+           endif
 
-           call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, dens(i,j,k), &
-                scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), &
-                rho(i,j,k),&
-                vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),&
-                gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
-                uxx,uxy,uxz,uyy,uyz,uzz,det, &
-                Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,&
-                epsnegative,GRHydro_C2P_failed(i,j,k))
-           
            
            if (epsnegative .ne. 0) then  
               
@@ -209,16 +200,13 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
               scon(i,j,k,:) = 0.d0
               vel(i,j,k,:) = 0.d0
               w_lorentz(i,j,k) = 1.d0
-#if USE_EOS_OMNI
+
               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)
-#else
-              press(i,j,k) = EOS_Pressure(GRHydro_polytrope_handle, rho(i,j,k), eps(i,j,k))
-              eps(i,j,k) = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, rho(i,j,k), press(i,j,k))
-#endif
+
               b2=gxx(i,j,k)*Bvec(i,j,k,1)**2+gyy(i,j,k)*Bvec(i,j,k,2)**2+gzz(i,j,k)*Bvec(i,j,k,3)**2+ &
                    2.0*(gxy(i,j,k)*Bvec(i,j,k,1)*Bvec(i,j,k,2)+gxz(i,j,k)*Bvec(i,j,k,1)*Bvec(i,j,k,3)+ &
                    gyz(i,j,k)*Bvec(i,j,k,2)*Bvec(i,j,k,3))
@@ -301,15 +289,7 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS)
 
   CCTK_REAL :: local_min_tracer
 
-#if !USE_EOS_OMNI
-#ifdef _EOS_BASE_INC_
-#undef _EOS_BASE_INC_
-#endif
-#include "EOS_Base.inc"
-#endif
-
 
-#if USE_EOS_OMNI
 ! begin EOS omni
   CCTK_INT  :: keyerr(1) = 0
   CCTK_INT  :: anyerr = 0
@@ -319,21 +299,15 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS)
   CCTK_REAL :: xeps = 0.0d0
   CCTK_REAL :: xtemp = 0.0d0
 ! end EOS omni
-#endif
+
   
 
-#if USE_EOS_OMNI
   ! this is a poly call
   call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
          GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr)
 
   call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
        GRHydro_rho_min,epsmin,xtemp,xye,pmin,epsmin,keyerr,anyerr)
-#else
-  pmin=EOS_Pressure(GRHydro_polytrope_handle, GRHydro_rho_min, 1.0d0) 
-  epsmin = EOS_SpecificIntEnergy(GRHydro_polytrope_handle, GRHydro_rho_min, pmin)
-#endif
-
                                                                               
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
   call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
@@ -393,6 +367,11 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS)
               end if
               
            endif
+
+           if(evolve_Y_e.ne.0) then
+              Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
+           endif
+           
            
            call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, densminus(i,j,k), &
                 sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), tauminus(i,j,k), &
@@ -574,6 +553,11 @@ subroutine Conservative2PrimitivePolytypeM(CCTK_ARGUMENTS)
               
            endif
 
+           if(evolve_Y_e.ne.0) then
+              Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
+           endif
+           
+
 !!$           call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, dens(i,j,k), &
 !!$                scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k),&
 !!$                rho(i,j,k),&
@@ -768,6 +752,11 @@ subroutine Conservative2PrimitiveGeneralM(CCTK_ARGUMENTS)
                  end if
               enddo
            endif
+
+           if(evolve_Y_e.ne.0) then
+              Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
+           endif
+           
            
            if ( (dens(i,j,k).le.sqrt(GRHydro_Det(i,j,k)) * &
                 GRHydro_rho_min*(1.0d0+GRHydro_atmo_tolerance)) &
@@ -883,6 +872,11 @@ subroutine Con2PrimPolytypeGeneralM(CCTK_ARGUMENTS)
             end if
          enddo
         end if
+
+        if(evolve_Y_e.ne.0) then
+           Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
+        endif
+
   
         if (dens(i,j,k).le.sqrt(GRHydro_Det(i,j,k)) * &
              GRHydro_rho_min*(1.0d0+GRHydro_atmo_tolerance)) then