User: rhaas
Date: 2013/08/13 09:55 AM
Modified:
/trunk/
param.ccl, schedule.ccl
/trunk/src/
GRHydro_Reconstruct.F90, GRHydro_WENOReconstruct.F90
Log:
GRHydro: Implemented WENO-z. Ability to choose between standard WENO5 and "adaptive epsilon"-variant from WHAM.
From: Christian Reisswig <reisswig(a)tapir.caltech.edu>
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_Reconstruct.F90
Delta lines: +1 -1
===================================================================
--- trunk/src/GRHydro_Reconstruct.F90 2013-07-15 22:33:13 UTC (rev 562)
+++ trunk/src/GRHydro_Reconstruct.F90 2013-08-13 14:55:54 UTC (rev 563)
@@ -135,7 +135,7 @@
! 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)
File [modified]: GRHydro_WENOReconstruct.F90
Delta lines: +111 -59
===================================================================
--- trunk/src/GRHydro_WENOReconstruct.F90 2013-07-15 22:33:13 UTC (rev 562)
+++ trunk/src/GRHydro_WENOReconstruct.F90 2013-08-13 14:55:54 UTC (rev 563)
@@ -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 @@
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
+ DECLARE_CCTK_FUNCTIONS
CCTK_INT :: allocstat
@@ -51,24 +53,44 @@
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 @@
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 @@
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 @@
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 @@
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
- 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
-
+ 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
+
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) &
Directory: /trunk/
==================
File [modified]: param.ccl
Delta lines: +10 -5
===================================================================
--- trunk/param.ccl 2013-07-15 22:33:13 UTC (rev 562)
+++ trunk/param.ccl 2013-08-13 14:55:54 UTC (rev 563)
@@ -136,11 +136,12 @@
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 @@
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"
File [modified]: schedule.ccl
Delta lines: +1 -1
===================================================================
--- trunk/schedule.ccl 2013-07-15 22:33:13 UTC (rev 562)
+++ trunk/schedule.ccl 2013-08-13 14:55:54 UTC (rev 563)
@@ -471,7 +471,7 @@
}
-if (CCTK_Equals(recon_method,"weno"))
+if (CCTK_Equals(recon_method,"weno") || CCTK_Equals(recon_method,"weno-z"))
{
schedule GRHydro_WENOSetup AT CCTK_Basegrid
