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

3 files changed, 54 insertions, 21 deletions

diff --git a/src/GRHydro_Eigenproblem.F90 b/src/GRHydro_Eigenproblem.F90
index 03576dd..69a6232 100644
--- a/src/GRHydro_Eigenproblem.F90
+++ b/src/GRHydro_Eigenproblem.F90
@@ -88,7 +88,7 @@ subroutine eigenvalues(handle,rho,velx,vely,velz,eps, &
   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
diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90
index 6a2fc27..6c59629 100644
--- a/src/GRHydro_Reconstruct.F90
+++ b/src/GRHydro_Reconstruct.F90
@@ -49,6 +49,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
   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 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
          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 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
                     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
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