From 0b806999f04f0faa576ac446c2eef09ccdf6909c Mon Sep 17 00:00:00 2001 From: rhaas Date: Tue, 13 Aug 2013 14:55:54 +0000 Subject: GRHydro: Implemented WENO-z. Ability to choose between standard WENO5 and "adaptive epsilon"-variant from WHAM. From: Christian Reisswig git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@563 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- param.ccl | 15 ++-- schedule.ccl | 2 +- src/GRHydro_Reconstruct.F90 | 2 +- src/GRHydro_WENOReconstruct.F90 | 170 ++++++++++++++++++++++++++-------------- 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) & -- cgit v1.2.3