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@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