GRHydro: Implemented WENO-z. Ability to choose between standard WENO5 and "adaptive epsilon"-variant from WHAM.

From: Christian Reisswig <reisswig@tapir.caltech.edu>

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

4 files changed, 123 insertions, 66 deletions

diff --git a/param.ccl b/param.ccl
index 726ee35..e7427cc 100644
--- a/param.ccl
+++ b/param.ccl
@@ -136,11 +136,12 @@ CCTK_INT GRHydro_MaxNumSandRVars "The maximum number of save and restore variabl
 
 keyword recon_method "Which reconstruction method to use" STEERABLE=RECOVER
 {
-  "tvd"	 :: "Slope limited TVD"
-  "ppm"  :: "PPM reconstruction"
-  "eno"  :: "ENO reconstruction"
-  "weno" :: "WENO reconstruction"
-  "mp5"  :: "MP5 reconstruction"
+  "tvd"	   :: "Slope limited TVD"
+  "ppm"    :: "PPM reconstruction"
+  "eno"    :: "ENO reconstruction"
+  "weno"   :: "WENO reconstruction"
+  "weno-z" :: "WENO-Z reconstruction"
+  "mp5"    :: "MP5 reconstruction"
 } "tvd"
 
 keyword method_type "Which type of method to use"
@@ -286,6 +287,10 @@ real weno_eps "WENO epsilon parameter"
   0:*  :: "small and positive"
 } 1e-26
 
+boolean weno_adaptive_epsilon "use modified smoothness indicators that take into account scale of solution (adaptive epsilon)"
+{
+} yes
+
 keyword riemann_solver "Which Riemann solver to use" STEERABLE=always
 {
   "Roe"		:: "Standard Roe solver"
diff --git a/schedule.ccl b/schedule.ccl
index 6bcf4d4..a3fbefb 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -471,7 +471,7 @@ if (CCTK_Equals(recon_method,"eno"))
 
 }
 
-if (CCTK_Equals(recon_method,"weno")) 
+if (CCTK_Equals(recon_method,"weno") || CCTK_Equals(recon_method,"weno-z"))
 {
 
   schedule GRHydro_WENOSetup AT CCTK_Basegrid
diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90
index cf6d103..14a5f31 100644
--- a/src/GRHydro_Reconstruct.F90
+++ b/src/GRHydro_Reconstruct.F90
@@ -135,7 +135,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
     ! this handles MHD and non-MHD
     call GRHydro_ENOReconstruct_drv(CCTK_PASS_FTOF)
 
-  else if (CCTK_EQUALS(recon_method,"weno")) then
+  else if (CCTK_EQUALS(recon_method,"weno") .or. CCTK_EQUALS(recon_method,"weno-z")) then
     ! this handles MHD and non-MHD
     call GRHydro_WENOReconstruct_drv(CCTK_PASS_FTOF)
   
diff --git a/src/GRHydro_WENOReconstruct.F90 b/src/GRHydro_WENOReconstruct.F90
index 64b5d18..635cdf0 100644
--- a/src/GRHydro_WENOReconstruct.F90
+++ b/src/GRHydro_WENOReconstruct.F90
@@ -11,6 +11,7 @@
 #include "cctk.h"
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
+#include "cctk_Functions.h"
 
  /*@@
    @routine    GRHydro_WENOSetup
@@ -39,6 +40,7 @@ subroutine GRHydro_WENOSetup(CCTK_ARGUMENTS)
 
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
 
   CCTK_INT :: allocstat
 
@@ -51,24 +53,44 @@ subroutine GRHydro_WENOSetup(CCTK_ARGUMENTS)
   endif
 
   ! Set stencils
-  weno_coeffs(1,1) =   3.0d0/8.0d0
-  weno_coeffs(1,2) =  -5.0d0/4.0d0
-  weno_coeffs(1,3) =  15.0d0/8.0d0
-  weno_coeffs(1,4) =   0.0d0
-  weno_coeffs(1,5) =   0.0d0
-  
-  weno_coeffs(2,1) =   0.0d0
-  weno_coeffs(2,2) =  -1.0d0/8.0d0
-  weno_coeffs(2,3) =   3.0d0/4.0d0
-  weno_coeffs(2,4) =   3.0d0/8.0d0
-  weno_coeffs(2,5) =   0.0d0
-  
-  weno_coeffs(3,1) =   0.0d0
-  weno_coeffs(3,2) =   0.0d0
-  weno_coeffs(3,3) =   3.0d0/8.0d0
-  weno_coeffs(3,4) =   3.0d0/4.0d0
-  weno_coeffs(3,5) =  -1.0d0/8.0d0
-
+  if (CCTK_EQUALS(recon_method,"weno-z")) then
+      weno_coeffs(1,1) =   2.0d0/6.0d0
+      weno_coeffs(1,2) =  -7.0d0/6.0d0
+      weno_coeffs(1,3) =  11.0d0/6.0d0
+      weno_coeffs(1,4) =   0.0d0
+      weno_coeffs(1,5) =   0.0d0
+      
+      weno_coeffs(2,1) =   0.0d0
+      weno_coeffs(2,2) =  -1.0d0/6.0d0
+      weno_coeffs(2,3) =   5.0d0/6.0d0
+      weno_coeffs(2,4) =   2.0d0/6.0d0
+      weno_coeffs(2,5) =   0.0d0
+      
+      weno_coeffs(3,1) =   0.0d0
+      weno_coeffs(3,2) =   0.0d0
+      weno_coeffs(3,3) =   2.0d0/6.0d0
+      weno_coeffs(3,4) =   5.0d0/6.0d0
+      weno_coeffs(3,5) =  -1.0d0/6.0d0
+  else
+      weno_coeffs(1,1) =   3.0d0/8.0d0
+      weno_coeffs(1,2) =  -5.0d0/4.0d0
+      weno_coeffs(1,3) =  15.0d0/8.0d0
+      weno_coeffs(1,4) =   0.0d0
+      weno_coeffs(1,5) =   0.0d0
+      
+      weno_coeffs(2,1) =   0.0d0
+      weno_coeffs(2,2) =  -1.0d0/8.0d0
+      weno_coeffs(2,3) =   3.0d0/4.0d0
+      weno_coeffs(2,4) =   3.0d0/8.0d0
+      weno_coeffs(2,5) =   0.0d0
+      
+      weno_coeffs(3,1) =   0.0d0
+      weno_coeffs(3,2) =   0.0d0
+      weno_coeffs(3,3) =   3.0d0/8.0d0
+      weno_coeffs(3,4) =   3.0d0/4.0d0
+      weno_coeffs(3,5) =  -1.0d0/8.0d0
+  endif
+   
   ! Shu smoothness indicator stencil coefficients
   beta_shu(1,1) =   4.0d0/3.0d0
   beta_shu(1,2) = -19.0d0/3.0d0
@@ -152,7 +174,8 @@ subroutine GRHydro_WENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, &
   implicit none
 
   DECLARE_CCTK_PARAMETERS
-
+  DECLARE_CCTK_FUNCTIONS
+  
   CCTK_INT :: order, nx, i, j
   CCTK_REAL, dimension(nx) :: v, vplus, vminus
 
@@ -160,6 +183,7 @@ subroutine GRHydro_WENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, &
   logical, dimension(nx) :: trivial_rp
   logical, dimension(nx) :: excise
   logical :: normal_weno
+  logical :: wenoZ
 
   CCTK_REAL :: beta1, beta2, beta3, vnorm, betanorm
   CCTK_REAL :: wplus1, wplus2, wplus3, wbarplus1, wbarplus2, wbarplus3
@@ -171,6 +195,12 @@ subroutine GRHydro_WENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, &
   excise = .false.
   trivial_rp = .false.
 
+  if (CCTK_EQUALS(recon_method,"weno-z")) then
+     wenoZ = .true.
+  else
+     wenoZ = .false.
+  endif
+  
 !!$    Initialize excision
   do i = 1, nx
     if (GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i) .ne. 0)) then
@@ -202,58 +232,80 @@ 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      &
-             + beta_shu(1,4)*v(i-2)*v(i)    &
-             + beta_shu(1,5)*v(i-1)*v(i)    &
-             + beta_shu(1,6)*v(i)**2
-    
-      beta2  = beta_shu(2,1)*v(i-1)**2      &
-             + beta_shu(2,2)*v(i-1)*v(i)    &
-             + beta_shu(2,3)*v(i)**2        &
-             + beta_shu(2,4)*v(i-1)*v(i+1)  &
-             + beta_shu(2,5)*v(i)*v(i+1)    &
-             + beta_shu(2,6)*v(i+1)**2
-    
-      beta3  = beta_shu(3,1)*v(i)**2       &
-             + beta_shu(3,2)*v(i)*v(i+1)   &
-             + beta_shu(3,3)*v(i+1)**2     &
-             + beta_shu(3,4)*v(i)*v(i+2)   &
-             + beta_shu(3,5)*v(i+1)*v(i+2) &
-             + beta_shu(3,6)*v(i+2)**2
-    
+      if (wenoZ) then
+
+         beta1 = 13.0d0/12.d0*(v(i-2)-2.0d0*v(i-1)+v(i))**2 + 1.0d0/4.d0*(v(i-2)-4.0d0*v(i-1)+3.0d0*v(i))**2
+         beta2 = 13.0d0/12.d0*(v(i-1)-2.0d0*v(i)+v(i+1))**2 + 1.0d0/4.d0*(v(i-1)-v(i+1))**2
+         beta3 = 13.0d0/12.d0*(v(i)-2.0d0*v(i+1)+v(i+2))**2 + 1.0d0/4.d0*(3.0d0*v(i)-4.0d0*v(i+1)+v(i+2))**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)
+      !!$    Compute weights according to weno-z alorithm
+         wbarplus1 = 1.0d0/10.0d0 * (1.0d0 + abs(beta1-beta3) / (weno_eps + beta1))
+         wbarplus2 = 3.0d0/5.0d0 * (1.0d0 + abs(beta1-beta3) / (weno_eps + beta2))
+         wbarplus3 = 3.0d0/10.0d0 * (1.0d0 + abs(beta1-beta3) / (weno_eps + beta3))
+
+         wbarminus1 = 3.0d0/10.0d0 * (1.0d0 + abs(beta1-beta3) / (weno_eps + beta1))
+         wbarminus2 = 3.0d0/5.0d0 * (1.0d0 + abs(beta1-beta3) / (weno_eps + beta2))
+         wbarminus3 = 1.0d0/10.0d0 * (1.0d0 + abs(beta1-beta3) / (weno_eps + beta3))
+
       
-      betanorm = beta1 + beta2 + beta3
+      else
       
-      beta1 = beta1 / betanorm
-      beta2 = beta2 / betanorm
-      beta3 = beta3 / betanorm
+   !!$    Compute smoothness indicators
+         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      &
+               + beta_shu(1,4)*v(i-2)*v(i)    &
+               + beta_shu(1,5)*v(i-1)*v(i)    &
+               + beta_shu(1,6)*v(i)**2
+      
+         beta2  = beta_shu(2,1)*v(i-1)**2      &
+               + beta_shu(2,2)*v(i-1)*v(i)    &
+               + beta_shu(2,3)*v(i)**2        &
+               + beta_shu(2,4)*v(i-1)*v(i+1)  &
+               + beta_shu(2,5)*v(i)*v(i+1)    &
+               + beta_shu(2,6)*v(i+1)**2
+      
+         beta3  = beta_shu(3,1)*v(i)**2       &
+               + beta_shu(3,2)*v(i)*v(i+1)   &
+               + beta_shu(3,3)*v(i+1)**2     &
+               + beta_shu(3,4)*v(i)*v(i+2)   &
+               + beta_shu(3,5)*v(i+1)*v(i+2) &
+               + beta_shu(3,6)*v(i+2)**2
+      
+   !!$    This is modification is suggested by Tchekhovskoy et al 2007 (WHAM code paper).
+         if (weno_adaptive_epsilon.ne.0) then
+            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
+         endif
+         
+         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
+      
+         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
+
+      endif
       
-      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
-    
       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) + wplus1 * weno_coeffs(1,j)*v(i-3+j) &