GRHydro: Con2PrimM: Set points to atmosphere if in excised region.

From: Christian Reisswig <reisswig@tapir.caltech.edu>

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

1 files changed, 44 insertions, 4 deletions

diff --git a/src/GRHydro_Con2PrimM.F90 b/src/GRHydro_Con2PrimM.F90
index d587d02..6517d41 100644
--- a/src/GRHydro_Con2PrimM.F90
+++ b/src/GRHydro_Con2PrimM.F90
@@ -63,7 +63,7 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
   CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, sdet, pmin(1), epsmin(1)
   CCTK_REAL :: oob, b2, d2, s2, bscon, bxhat, byhat, bzhat, bhatscon
   CCTK_REAL :: Wm, Wm0, Wm_resid, Wmold
-  CCTK_REAL :: s2m, s2m0, s2m_resid, s2mold, s2max, taum
+  CCTK_REAL :: s2m, s2m0, s2m_resid, s2mold, s2max, taum, dummy1, dummy2
   CCTK_INT :: niter
   CCTK_INT :: epsnegative
   character(len=100) warnline
@@ -167,14 +167,17 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
   !$OMP sdet,d2,s2,oob,bscon,bxhat,byhat,bzhat, &
   !$OMP bhatscon,Wm,Wm0,Wm_resid,Wmold,s2m,s2m0,s2m_resid,s2mold,s2max, &
   !$OMP taum,niter,rho_tmp,press_tmp,eps_tmp,velx_tmp,vely_tmp,velz_tmp, &
-  !$OMP w_lorentz_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,bdotv,magpress)
+  !$OMP w_lorentz_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,bdotv,magpress, dummy1, dummy2)
   do k = 1, nz 
      do j = 1, ny 
         do i = 1, nx
            
+           !do not compute if in atmosphere region!
+           if (atmosphere_mask(i,j,k) .gt. 0) cycle
+           
            !do not compute if in atmosphere or in excised region
-           if ((atmosphere_mask(i,j,k) .ne. 0) .or. &
-                (hydro_excision_mask(i,j,k) .ne. 0)) cycle
+           !if ((atmosphere_mask(i,j,k) .ne. 0) .or. &
+           !     (hydro_excision_mask(i,j,k) .ne. 0)) cycle
            
            epsnegative = 0
            
@@ -201,6 +204,43 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
               
            endif
            
+           ! In excised region, set to atmosphere!    
+           if (GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .gt. 0)) then
+              SET_ATMO_MIN(dens(i,j,k), sqrt(det)*GRHydro_rho_min, r(i,j,k)) !sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance)
+              SET_ATMO_MIN(rho(i,j,k), GRHydro_rho_min, r(i,j,k)) !GRHydro_rho_min
+              scon(i,j,k,:) = 0.d0
+              vup(i,j,k,:) = 0.d0
+              w_lorentz(i,j,k) = 1.d0
+              
+              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
+                keytemp = 0
+                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
+
+              !!$ tau does need to take into account the existing B-field
+              !!$ with w_lorentz=1, we find tau = sqrtdet*(rho (1+eps+b^2/2)) - dens  [Press drops out]
+              tau(i,j,k)  = sdet * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k)
+              
+              if(tau(i,j,k).le.sdet*b2*0.5d0)then
+                tau(i,j,k) = GRHydro_tau_min + sdet*b2*0.5d0
+              endif
+               
+              cycle
+           
+           endif
+           
            if(evolve_Y_e.ne.0) then
               Y_e(i,j,k) = max(min(Y_e_con(i,j,k) / dens(i,j,k),GRHydro_Y_e_max),&
                  GRHydro_Y_e_min)