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

[roland.haas at physics.gatech.edu]

User: rhaas
Date: 2012/05/13 02:01 PM

Modified:
 /trunk/src/
  GRHydro_Con2PrimM.F90

Log:
 Error checking for C2P_failed > 1, not only equals 1
 and check for atmosphere using density after
 Con2PrimM_pt was performed.

File Changes:

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

File [modified]: GRHydro_Con2PrimM.F90
Delta lines: +28 -1
===================================================================
--- trunk/src/GRHydro_Con2PrimM.F90	2012-05-11 22:46:36 UTC (rev 319)
+++ trunk/src/GRHydro_Con2PrimM.F90	2012-05-13 19:01:37 UTC (rev 320)
@@ -243,6 +243,34 @@
 #endif          
               
            end if
+           
+           if ( rho(i,j,k) .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) .or. GRHydro_C2P_failed(i,j,k) .ge. 1) 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+ &
+                   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))
+              
+              dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance)
+              rho(i,j,k) = GRHydro_rho_min
+              scon(i,j,k,:) = 0.d0
+              vel(i,j,k,:) = 0.d0
+              w_lorentz(i,j,k) = 1.d0
+              xtemp = 0.0d0
+   
+              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)
+
+              ! w_lorentz=1, so the expression for tau reduces to:
+
+              !!$ 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)  = sqrt(det) * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k)
+
+           end if
+            
         end do
      end do
   end do
@@ -487,7 +515,6 @@
                  !$OMP END CRITICAL
               endif
            endif
-           
         end do
      end do
   end do