From b591ed48992eaf73a7ee15844a21c1407d245c5a Mon Sep 17 00:00:00 2001 From: rhaas Date: Sat, 6 Jul 2013 18:10:53 +0000 Subject: GRHydro: Con2PrimM: Set points to atmosphere if in excised region. From: Christian Reisswig git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@551 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- src/GRHydro_Con2PrimM.F90 | 48 +++++++++++++++++++++++++++++++++++++++++++---- 1 file 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) -- cgit v1.2.3