[Commits] [svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 551)
rhaas at tapir.caltech.edu
rhaas at tapir.caltech.edu
Sat Jul 6 13:10:54 CDT 2013
User: rhaas
Date: 2013/07/06 01:10 PM
Modified:
/trunk/src/
GRHydro_Con2PrimM.F90
Log:
GRHydro: Con2PrimM: Set points to atmosphere if in excised region.
From: Christian Reisswig <reisswig at tapir.caltech.edu>
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_Con2PrimM.F90
Delta lines: +44 -4
===================================================================
--- trunk/src/GRHydro_Con2PrimM.F90 2013-07-06 18:10:49 UTC (rev 550)
+++ trunk/src/GRHydro_Con2PrimM.F90 2013-07-06 18:10:53 UTC (rev 551)
@@ -63,7 +63,7 @@
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 @@
!$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 @@
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)
More information about the Commits
mailing list