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@scriwalker.(none)>

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@466 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45

1 files changed, 39 insertions, 19 deletions

diff --git a/src/GRHydro_WENOReconstruct.F90 b/src/GRHydro_WENOReconstruct.F90
index e7273a0..fb03f66 100644
--- a/src/GRHydro_WENOReconstruct.F90
+++ b/src/GRHydro_WENOReconstruct.F90
@@ -158,14 +158,15 @@ subroutine GRHydro_WENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, &
 
   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 @@ subroutine GRHydro_WENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, &
   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 @@ subroutine GRHydro_WENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, &
     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 @@ subroutine GRHydro_WENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, &
              + 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