[Commits] [svn:einsteintoolkit] GRHydro/trunk/ (Rev. 466)
rhaas at tapir.caltech.edu
rhaas at tapir.caltech.edu
Mon Jan 14 08:23:54 CST 2013
User: rhaas
Date: 2013/01/14 08:23 AM
Modified:
/trunk/
param.ccl
/trunk/src/
GRHydro_Eigenproblem.F90, GRHydro_Reconstruct.F90, GRHydro_WENOReconstruct.F90
Log:
GRHydro: Improved WENO implementation: (i) Pick adaptive epsilon-parameter according to WHAM code paper, (ii) Catch epsilon<0 for non-microphysical EOS. Tested with shock tube and TOV.
From: Christian Reisswig <reisswig at scriwalker.(none)>
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_Eigenproblem.F90
Delta lines: +1 -1
===================================================================
--- trunk/src/GRHydro_Eigenproblem.F90 2013-01-14 14:23:52 UTC (rev 465)
+++ trunk/src/GRHydro_Eigenproblem.F90 2013-01-14 14:23:54 UTC (rev 466)
@@ -88,7 +88,7 @@
if(cs2.lt.0.0d0) then
!$OMP CRITICAL
if (abs(cs2) .gt. 1.0d-4) then
- write(warnline,'(a50,6g16.7)') 'abs(cs2), rho, dpdrho, press*dpdeps/rho**2, eps, press/rho: ', abs(cs2), rho, dpdrho, press * dpdeps / (rho**2), eps, press/rho
+ write(warnline,'(a60,6g16.7)') 'abs(cs2), rho, dpdrho, press*dpdeps/rho**2, eps, press/rho: ', abs(cs2), rho, dpdrho, press * dpdeps / (rho**2), eps, press/rho
call CCTK_WARN(1,warnline)
call CCTK_WARN(1,"cs2 < 0! Check speed of sound calculation!")
cs2 = 0.0d0
File [modified]: GRHydro_Reconstruct.F90
Delta lines: +14 -1
===================================================================
--- trunk/src/GRHydro_Reconstruct.F90 2013-01-14 14:23:52 UTC (rev 465)
+++ trunk/src/GRHydro_Reconstruct.F90 2013-01-14 14:23:54 UTC (rev 466)
@@ -49,6 +49,7 @@
CCTK_REAL :: local_min_tracer, dummy1, dummy2
CCTK_INT :: type_bits, not_trivial
CCTK_REAL :: agxx, agxy, agxz, agyy, agyz, agzz, w
+ character*256 :: warnline
! save memory when MP is not used
CCTK_INT :: GRHydro_UseGeneralCoordinates
@@ -130,7 +131,7 @@
SET_ATMO_MIN(dummy2, GRHydro_rho_min, r(i,j,k))
if(rhoplus(i,j,k).lt.dummy2 .or. &
rhominus(i,j,k).lt.dummy2) then
-
+ !.or. epsplus(i,j,k) .lt. 0.0d0 .or. epsminus(i,j,k) .lt. 0.0d0) then
rhoplus(i,j,k) = rho(i,j,k)
rhominus(i,j,k) = rho(i,j,k)
epsplus(i,j,k) = eps(i,j,k)
@@ -198,9 +199,21 @@
end where
end if
end if
+
+ ! Check if epsilon became negative for ideal gas EOS.
+ if (evolve_temper.eq.0) then
+ if (epsplus(i,j,k) .lt. 0.0d0) then
+ epsplus(i,j,k) = eps(i,j,k)
+ endif
+ if (epsminus(i,j,k) .lt. 0.0d0) then
+ epsminus(i,j,k) = eps(i,j,k)
+ endif
+ endif
+
! Riemann problem might not be trivial anymore!!!!
SpaceMask_SetStateBitsF90(space_mask, i-xoffset, j-yoffset, k-zoffset, type_bits, not_trivial)
SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, not_trivial)
+
enddo
enddo
enddo
File [modified]: GRHydro_WENOReconstruct.F90
Delta lines: +39 -19
===================================================================
--- trunk/src/GRHydro_WENOReconstruct.F90 2013-01-14 14:23:52 UTC (rev 465)
+++ trunk/src/GRHydro_WENOReconstruct.F90 2013-01-14 14:23:54 UTC (rev 466)
@@ -158,14 +158,15 @@
CCTK_INT :: order, nx, i, j, k, r
CCTK_REAL, dimension(nx) :: v, vplus, vminus
- CCTK_REAL, dimension(order, 1-order:nx+order) :: vdiff
CCTK_INT, dimension(nx) :: hydro_excision_mask
logical, dimension(nx) :: trivial_rp
logical, dimension(nx) :: excise
logical :: normal_weno
- CCTK_REAL :: large = 1.d10, gamma1, gamma2, gamma3, beta1, beta2, beta3, w1, w2, w3, wbar1, wbar2, wbar3
+ CCTK_REAL :: large = 1.d10, gamma1, gamma2, gamma3, beta1, beta2, beta3, vnorm, betanorm
+ CCTK_REAL :: wplus1, wplus2, wplus3, wbarplus1, wbarplus2, wbarplus3
+ CCTK_REAL :: wminus1, wminus2, wminus3, wbarminus1, wbarminus2, wbarminus3
vminus = 0.d0
vplus = 0.d0
@@ -173,12 +174,9 @@
excise = .false.
trivial_rp = .false.
-
-
!!$ Initialize excision
do i = 1, nx
if (GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i) .ne. 0)) then
- vdiff(1, i) = large * (-1.d0*(1.d0+mod(i,10)))**i
trivial_rp(i) = .true.
excise(i) = .true.
if (i > 1) then
@@ -208,6 +206,7 @@
if (normal_weno) then
!!$ Compute smoothness indicators
+!!$ This is from Tchekhovskoy et al 2007 (WHAM code paper).
beta1 = beta_shu(1,1)*v(i-2)**2 &
+ beta_shu(1,2)*v(i-2)*v(i-1) &
+ beta_shu(1,3)*v(i-1)**2 &
@@ -230,25 +229,46 @@
+ beta_shu(3,6)*v(i+2)**2
- wbar1 = 1.0d0/16.d0 / (weno_eps + beta1)**2
- wbar2 = 5.0d0/8.d0 / (weno_eps + beta2)**2
- wbar3 = 5.0d0/16.d0 / (weno_eps + beta3)**2
+ vnorm = (v(i-2)**2 + v(i-1)**2 + v(i)**2 + v(i+1)**2 + v(i+2)**2)
+
+ beta1 = beta1 + 100.0d0*weno_eps*(vnorm + 1.0d0)
+ beta2 = beta2 + 100.0d0*weno_eps*(vnorm + 1.0d0)
+ beta3 = beta3 + 100.0d0*weno_eps*(vnorm + 1.0d0)
+
+ betanorm = beta1 + beta2 + beta3
+
+ beta1 = beta1 / betanorm
+ beta2 = beta2 / betanorm
+ beta3 = beta3 / betanorm
+
+ wbarplus1 = 1.0d0/16.0d0 / (weno_eps + beta1)**2
+ wbarplus2 = 5.0d0/8.0d0 / (weno_eps + beta2)**2
+ wbarplus3 = 5.0d0/16.0d0 / (weno_eps + beta3)**2
- w1 = wbar1 / (wbar1 + wbar2 + wbar3)
- w2 = wbar2 / (wbar1 + wbar2 + wbar3)
- w3 = wbar3 / (wbar1 + wbar2 + wbar3)
+ wplus1 = wbarplus1 / (wbarplus1 + wbarplus2 + wbarplus3)
+ wplus2 = wbarplus2 / (wbarplus1 + wbarplus2 + wbarplus3)
+ wplus3 = wbarplus3 / (wbarplus1 + wbarplus2 + wbarplus3)
+
+ wbarminus1 = 5.0d0/16.0d0 / (weno_eps + beta1)**2
+ wbarminus2 = 5.0d0/8.0d0 / (weno_eps + beta2)**2
+ wbarminus3 = 1.0d0/16.0d0 / (weno_eps + beta3)**2
+ wminus1 = wbarminus1 / (wbarminus1 + wbarminus2 + wbarminus3)
+ wminus2 = wbarminus2 / (wbarminus1 + wbarminus2 + wbarminus3)
+ wminus3 = wbarminus3 / (wbarminus1 + wbarminus2 + wbarminus3)
+
!!$ Calculate the reconstruction
do j = 1, 5
- vplus(i) = vplus(i) + w1 * weno_coeffs(1,j)*v(i-3+j) &
- + w2 * weno_coeffs(2,j)*v(i-3+j) &
- + w3 * weno_coeffs(3,j)*v(i-3+j)
- vminus(i) = vminus(i) + w1 * weno_coeffs(3,6-j)*v(i-3+j) &
- + w2 * weno_coeffs(2,6-j)*v(i-3+j) &
- + w3 * weno_coeffs(1,6-j)*v(i-3+j)
-
+ vplus(i) = vplus(i) + wplus1 * weno_coeffs(1,j)*v(i-3+j) &
+ + wplus2 * weno_coeffs(2,j)*v(i-3+j) &
+ + wplus3 * weno_coeffs(3,j)*v(i-3+j)
+ vminus(i) = vminus(i) + wminus1 * weno_coeffs(3,6-j)*v(i-3+j) &
+ + wminus2 * weno_coeffs(2,6-j)*v(i-3+j) &
+ + wminus3 * weno_coeffs(1,6-j)*v(i-3+j)
end do
-
+ !vminus(i) = v(i)
+ !vplus(i) = v(i)
+
end if
end do
Directory: /trunk/
==================
File [modified]: param.ccl
Delta lines: +1 -1
===================================================================
--- trunk/param.ccl 2013-01-14 14:23:52 UTC (rev 465)
+++ trunk/param.ccl 2013-01-14 14:23:54 UTC (rev 466)
@@ -265,7 +265,7 @@
real weno_eps "WENO epsilon parameter"
{
0:* :: "small and positive"
-} 1e-6
+} 1e-26
keyword riemann_solver "Which Riemann solver to use"
{
More information about the Commits
mailing list