[Commits] [svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 542)

[rhaas at tapir.caltech.edu]

User: rhaas
Date: 2013/07/06 01:10 PM

Modified:
 /trunk/src/
  GRHydro_Con2PrimHot.F90

Log:
 GRHydro: Critical bugfix: Con2Prim was probably never executed in post_recover_variables or for new points in post_regrid since
 excision mask does typically not get initialized before MoL_PostStep!
 
 Also: Instead of aborting Con2Prim, set points to atmopshere for a point where hydro_excision_mask > 0.
 
 From: Christian Reisswig <reisswig at scriwalker.(none)>

File Changes:

Directory: /trunk/src/
======================

File [modified]: GRHydro_Con2PrimHot.F90
Delta lines: +7 -4
===================================================================
--- trunk/src/GRHydro_Con2PrimHot.F90	2013-07-06 18:10:12 UTC (rev 541)
+++ trunk/src/GRHydro_Con2PrimHot.F90	2013-07-06 18:10:16 UTC (rev 542)
@@ -35,7 +35,7 @@
   CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup
 
 ! begin EOS Omni vars
-  integer :: n,keytemp,anyerr,keyerr(1)
+  integer :: n,keytemp,anyerr,keyerr(1), reset_to_atmo
   real*8  :: xpress(1),xeps(1),xtemp(1),xye(1),xrho(1)
   n = 1;keytemp = 0;anyerr = 0;keyerr(1) = 0
   xpress = 0.0d0;xeps = 0.0d0;xtemp = 0.0d0;xye = 0.0d0
@@ -89,9 +89,8 @@
       do i = 1, nx
 
 
-        !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
+        !do not compute if in atmosphere
+        if (atmosphere_mask(i,j,k) .gt. 0) cycle
          
         epsnegative = .false.
 
@@ -121,8 +120,12 @@
                 GRHydro_Y_e_min)
         endif
 
+        reset_to_atmo = 0
         IF_BELOW_ATMO(dens(i,j,k), sqrt(det)*GRHydro_rho_min, GRHydro_atmo_tolerance, r(i,j,k)) then
+           reset_to_atmo = 1
+        endif
 
+        if (reset_to_atmo .gt. 0 .or. 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)) 
            SET_ATMO_MIN(rho(i,j,k), GRHydro_rho_min, r(i,j,k)) 
            scon(i,j,k,:) = 0.d0