include Y_e_con in atmosphere handling

original commit by Christian Ott (cott)



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

2 files changed, 9 insertions, 6 deletions

diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90
index 8279057..b1e1507 100644
--- a/src/GRHydro_Con2Prim.F90
+++ b/src/GRHydro_Con2Prim.F90
@@ -49,6 +49,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
   
+  integer :: CCTK_MyProc
   integer :: i, j, k, itracer, nx, ny, nz
   CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, pmin, epsmin
   logical :: epsnegative
@@ -62,7 +63,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
 
 ! begin EOS Omni vars
   integer :: n,keytemp,anyerr,keyerr(1)
-  real*8  :: xpress,xeps,xtemp,xye 
+  real*8  :: xpress(1),xeps(1),xtemp(1),xye(1)
   n = 1;keytemp = 0;anyerr = 0;keyerr(1) = 0
   xpress = 0.0d0;xeps = 0.0d0;xtemp = 0.0d0;xye = 0.0d0
 ! end EOS Omni vars
@@ -171,7 +172,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
                  GRHydro_reflevel, GRHydro_C2P_failed(i,j,k))
          else
             call Con2Prim_pt_hot(i,j,k,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),&
+                 scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(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),temperature(i,j,k),y_e(i,j,k),&
                  press(i,j,k),w_lorentz(i,j,k), &
                  uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), &
@@ -184,7 +185,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
                GRHydro_C2P_failed(i,j,k) = 0
                local_perc_ptol = GRHydro_perc_ptol*100.0d0
                call Con2Prim_pt_hot(i,j,k,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),&
+                    scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(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),&
                     temperature(i,j,k),y_e(i,j,k),press(i,j,k),w_lorentz(i,j,k), &
                     uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), &
@@ -539,7 +540,7 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, &
 
 end subroutine Con2Prim_pt
 
-subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, vely, &
+subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, ye_con, rho, velx, vely, &
      velz, epsilon, temp, ye, press, w_lorentz, uxx, uxy, uxz, uyy, &
      uyz, uzz, det, x, y, z, r, epsnegative, GRHydro_rho_min, pmin, epsmin, &
      GRHydro_reflevel, GRHydro_C2P_failed, local_perc_ptol)
@@ -548,7 +549,7 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
   
   DECLARE_CCTK_PARAMETERS
 
-  CCTK_REAL dens, sx, sy, sz, tau, rho, velx, vely, velz, epsilon, &
+  CCTK_REAL dens, sx, sy, sz, tau, ye_con, rho, velx, vely, velz, epsilon, &
        press, uxx, uxy, uxz, uyy, uyz, uzz, det, w_lorentz, x, &
        y, z, r, GRHydro_rho_min
   CCTK_REAL temp, ye
@@ -895,6 +896,7 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve
     dens = sqrt(det) * rho
     temp = GRHydro_hot_atmo_temp
     ye = GRHydro_hot_atmo_Y_e
+    ye_con = dens * ye
     keytemp=1
     call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n, &
          rho,epsilon,temp,ye,xpress,keyerr,anyerr)
diff --git a/src/GRHydro_UpdateMask.F90 b/src/GRHydro_UpdateMask.F90
index d408533..d6bf540 100644
--- a/src/GRHydro_UpdateMask.F90
+++ b/src/GRHydro_UpdateMask.F90
@@ -41,7 +41,7 @@ subroutine GRHydroUpdateAtmosphereMask(CCTK_ARGUMENTS)
 
   frac = CCTK_DELTA_TIME
 
-  if(evolve_temper.ne.1) then
+  if(evolve_temper.ne.1.and.evolve_Y_e.ne.1) then
      !$OMP PARALLEL DO PRIVATE(k,j,i)
      do k = 1, cctk_lsh(3)
         do j = 1, cctk_lsh(2)
@@ -72,6 +72,7 @@ subroutine GRHydroUpdateAtmosphereMask(CCTK_ARGUMENTS)
                    (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, \
                  type_bits, atmosphere)) .or. &
                       (dens(i,j,k) + frac * densrhs(i,j,k) .le. 0.d0) ) then
+                 y_e_con_rhs(i,j,k) = 0.0d0
                  densrhs(i,j,k) = 0.0d0
                  srhs(i,j,k,:)   = 0.0d0
                  taurhs(i,j,k)  = 0.0d0