GRHydro: Implemented WENO5. Tested with a shock tube.

From: Christian Reisswig <reisswig@scriwalker.(none)>

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

1 files changed, 276 insertions, 0 deletions

diff --git a/src/GRHydro_WENOReconstruct.F90 b/src/GRHydro_WENOReconstruct.F90
new file mode 100644
index 0000000..e7273a0
--- /dev/null
+++ b/src/GRHydro_WENOReconstruct.F90
@@ -0,0 +1,276 @@
+ /*@@
+   @file      GRHydro_WENOReconstruct.F90
+   @date      Fri Jan 3 2013
+   @author    Ian Hawke, Christian Reisswig
+   @desc 
+   Routines to set up the coefficient array and to perform one dimensional 
+   ENO reconstruction of arbitrary order.
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+ /*@@
+   @routine    GRHydro_WENOSetup
+   @date       Fri Jan 3 2013
+   @author     Christian Reisswig
+   @desc 
+   Sets up the coefficient array for WENO reconstruction. 
+   Uses the notation of Shu, equation (2.21), in 
+   ''High Order ENO and WENO Schemes for CFD''
+   (see ThornGuide for full reference).
+   One exception: (Shu) r -> (Here) i: avoiding name clash.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_WENOSetup(CCTK_ARGUMENTS)
+
+  USE GRHydro_WENOScalars
+
+  implicit none
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT :: i, j, k, l, m, q, allocstat
+  CCTK_REAL :: denominator, numerator, numerator_product
+
+  if(.not.coeffs_allocated) then
+     ! Right now we hardcode to 5th order
+     allocate(weno_coeffs(3, 5), STAT=allocstat)
+     if (allocstat .ne. 0) call CCTK_WARN(0, "Failed to allocate WENO coefficient arrays!")
+     allocate(beta_shu(3, 6), STAT=allocstat)
+     if (allocstat .ne. 0) call CCTK_WARN(0, "Failed to allocate smoothness indicator stencil coefficient arrays!")
+     coeffs_allocated = .true.
+  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
+
+  ! Shu smoothness indicator stencil coefficients
+  beta_shu(1,1) =   4.0d0/3.0d0
+  beta_shu(1,2) = -19.0d0/3.0d0
+  beta_shu(1,3) =  25.0d0/3.0d0
+  beta_shu(1,4) =  11.0d0/3.0d0
+  beta_shu(1,5) = -31.0d0/3.0d0
+  beta_shu(1,6) =  10.0d0/3.0d0
+  
+  beta_shu(2,1) =   4.0d0/3.0d0
+  beta_shu(2,2) = -13.0d0/3.0d0
+  beta_shu(2,3) =  13.0d0/3.0d0
+  beta_shu(2,4) =   5.0d0/3.0d0
+  beta_shu(2,5) = -13.0d0/3.0d0
+  beta_shu(2,6) =   4.0d0/3.0d0
+  
+  beta_shu(3,1) =  10.0d0/3.0d0
+  beta_shu(3,2) = -31.0d0/3.0d0
+  beta_shu(3,3) =  25.0d0/3.0d0
+  beta_shu(3,4) =  11.0d0/3.0d0
+  beta_shu(3,5) = -19.0d0/3.0d0
+  beta_shu(3,6) =   4.0d0/3.0d0
+  
+end subroutine GRHydro_WENOSetup
+
+ /*@@
+   @routine    GRHydro_ENOShutdown
+   @date       Fri Jan  3 2013
+   @author     Christian Reisswig
+   @desc 
+   Deallocates the coefficient arrays
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_WENOShutdown(CCTK_ARGUMENTS)
+
+  USE GRHydro_WENOScalars
+
+  implicit none
+
+  DECLARE_CCTK_ARGUMENTS
+
+  CCTK_INT :: deallocstat
+
+  if(coeffs_allocated) then
+     deallocate(weno_coeffs, STAT = deallocstat)
+     if (deallocstat .ne. 0) call CCTK_WARN(0, "Failed to deallocate WENO coefficients.")
+     deallocate(beta_shu, STAT = deallocstat)
+     if (deallocstat .ne. 0) call CCTK_WARN(0, "Failed to deallocate shu smoothness indicator coefficients.")
+     coeffs_allocated = .false.
+  endif
+
+end subroutine GRHydro_WENOShutdown
+
+ /*@@
+   @routine    GRHydro_ENOReconstruct1d
+   @date       Fri Jan 3 2013
+   @author     Christian Reisswig
+   @desc 
+   Perform a one dimensional reconstruction of a given array using WENO.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+#define SpaceMask_CheckStateBitsF90_1D(mask,i,type_bits,state_bits) \
+  (iand(mask((i)),(type_bits)).eq.(state_bits))
+
+subroutine GRHydro_WENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, &
+     hydro_excision_mask)
+  USE GRHydro_WENOScalars
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+
+  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
+
+  vminus = 0.d0
+  vplus = 0.d0
+
+  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
+        trivial_rp(i-1) = .true.
+      end if
+    end if
+  end do
+
+  do i = 3, nx-2
+!!$    Handle excision
+    normal_weno = .true.
+    if (i < nx) then
+     if (excise(i+1)) then
+      vminus(i) = v(i)
+      vplus(i) = v(i)
+      normal_weno = .false.
+     end if
+    end if
+    if (i > 1) then
+     if (excise(i-1)) then
+      vminus(i) = v(i)
+      vplus(i) = v(i)
+      normal_weno = .false.
+     end if
+    end if
+
+    if (normal_weno) then
+
+!!$    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
+    
+      
+      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
+    
+      w1 = wbar1 / (wbar1 + wbar2 + wbar3)
+      w2 = wbar2 / (wbar1 + wbar2 + wbar3)
+      w3 = wbar3 / (wbar1 + wbar2 + wbar3)
+    
+!!$    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)
+      
+      end do
+      
+    end if
+  end do
+
+  do i = 1, nx
+    if (excise(i)) then
+      if (i > 1) then
+        if (.not. excise(i-1)) then
+          vminus(i) = vplus(i-1)
+        end if
+      end if
+      vplus(i) = vminus(i)
+    end if
+  end do
+  do i = nx, 1, -1
+    if (excise(i)) then
+      if (i < nx) then
+        if (.not. excise(i+1)) then
+          vplus(i) = vminus(i+1)
+        end if
+      end if
+      vminus(i) = vplus(i)
+    end if
+  end do
+
+end subroutine GRHydro_WENOReconstruct1d