diff --git a/sources/interpolations.F90 b/sources/interpolations.F90 index edb5a3a..6b937bf 100644 --- a/sources/interpolations.F90 +++ b/sources/interpolations.F90 @@ -52,6 +52,16 @@ module interpolations real(kind=8), intent(in) :: x, a, b real(kind=8) :: c end function + subroutine weno_smoothness_iface(q, betas) + real(kind=8), dimension(:) , intent(in ) :: q + real(kind=8), dimension(:,:), intent( out) :: betas + end subroutine + subroutine weno_weights_iface(h, q, betas, alphas) + real(kind=8) , intent(in ) :: h + real(kind=8), dimension(:) , intent(in ) :: q + real(kind=8), dimension(:,:), intent(in ) :: betas + real(kind=8), dimension(:,:), intent( out) :: alphas + end subroutine end interface ! pointers to the reconstruction and limiter procedures @@ -62,6 +72,11 @@ module interpolations procedure(limiter_iface) , pointer, save :: limiter_prol => null() procedure(limiter_iface) , pointer, save :: limiter_clip => null() +! pointer to WENO scheme subroutines +! + procedure(weno_smoothness_iface), pointer :: weno_smoothness_measurements + procedure(weno_weights_iface) , pointer :: weno_stencil_weights + ! method names ! character(len=255), save :: name_rec = "" @@ -100,6 +115,12 @@ module interpolations ! real(kind=8), save :: ppm_const = 1.25d+00 +! WENO parameters +! + logical :: weno_scheme = .false. + integer :: weno_power = 1 + real(kind=8) :: weno_sensitivity = 1.0d-40 + ! weight toward central scheme for low-dissipation methods ! real(kind=8), save :: cweight = 0.0d+00 @@ -269,42 +290,131 @@ module interpolations reconstruct_states => reconstruct_ppm order = 3 nghosts = max(nghosts, 4) + case ("weno5js", "weno5-js", "WENO5JS", "WENO5-JS") + name_rec = "5th order WENO-JS (Jiang & Shu, 1996)" + interfaces => interfaces_dir + reconstruct_states => reconstruct_weno5 + order = 5 + nghosts = max(nghosts, 4) + weno_scheme = .true. + weno_power = 2 + weno_sensitivity = 1.0d-06 + weno_smoothness_measurements => weno5_smoothness_measurements_js + weno_stencil_weights => weno5_stencil_weights_js case ("weno5z", "weno5-z", "WENO5Z", "WENO5-Z") name_rec = "5th order WENO-Z (Borges et al. 2008)" interfaces => interfaces_dir - reconstruct_states => reconstruct_weno5z + reconstruct_states => reconstruct_weno5 order = 5 nghosts = max(nghosts, 4) + weno_scheme = .true. + weno_power = 2 + weno_sensitivity = 1.0d-40 + weno_smoothness_measurements => weno5_smoothness_measurements_js + weno_stencil_weights => weno5_stencil_weights_z case ("weno5yc", "weno5-yc", "WENO5YC", "WENO5-YC") name_rec = "5th order WENO-YC (Yamaleev & Carpenter 2009)" interfaces => interfaces_dir - reconstruct_states => reconstruct_weno5yc + reconstruct_states => reconstruct_weno5 order = 5 nghosts = max(nghosts, 4) + weno_scheme = .true. + weno_power = 2 + weno_sensitivity = 1.0d-40 + weno_smoothness_measurements => weno5_smoothness_measurements_js + weno_stencil_weights => weno5_stencil_weights_yc case ("weno5ns", "weno5-ns", "WENO5NS", "WENO5-NS") name_rec = "5th order WENO-NS (Ha et al. 2013)" interfaces => interfaces_dir reconstruct_states => reconstruct_weno5ns order = 5 nghosts = max(nghosts, 4) + weno_scheme = .true. + weno_power = 2 + weno_sensitivity = 1.0d-40 + weno_smoothness_measurements => weno5_smoothness_measurements_js + weno_stencil_weights => weno5_stencil_weights_js + case ("weno5z+", "weno5-z+", "WENO5Z+", "WENO5-Z+") + name_rec = "5th order WENO-Z+ (Acker et al. 2016)" + interfaces => interfaces_dir + reconstruct_states => reconstruct_weno5 + order = 5 + nghosts = max(nghosts, 4) + weno_scheme = .true. + weno_power = 2 + weno_sensitivity = 1.0d-40 + weno_smoothness_measurements => weno5_smoothness_measurements_js + weno_stencil_weights => weno5_stencil_weights_zp + case ("weno5at", "weno5-at", "WENO5AT", "WENO5-AT") + name_rec = "5th order WENO-AT (Huang et al. 2022)" + interfaces => interfaces_dir + reconstruct_states => reconstruct_weno5 + order = 5 + nghosts = max(nghosts, 4) + weno_scheme = .true. + weno_power = 2 + weno_sensitivity = 1.0d-40 + weno_smoothness_measurements => weno5_smoothness_measurements_js + weno_stencil_weights => weno5_stencil_weights_at + case ("weno5zc", "weno5-zc", "WENO5ZC", "WENO5-ZC") + name_rec = "5th order WENO-ZC (Baretto et al. 2023)" + interfaces => interfaces_dir + reconstruct_states => reconstruct_weno5 + order = 5 + nghosts = max(nghosts, 4) + weno_scheme = .true. + weno_power = 2 + weno_sensitivity = 1.0d-40 + weno_smoothness_measurements => weno5_smoothness_measurements_js + weno_stencil_weights => weno5_stencil_weights_zc + case ("weno5zc+", "weno5-zc+", "WENO5ZC+", "WENO5-ZC+") + name_rec = "5th order WENO-ZC+ (Baretto et al. 2023)" + interfaces => interfaces_dir + reconstruct_states => reconstruct_weno5 + order = 5 + nghosts = max(nghosts, 4) + weno_scheme = .true. + weno_power = 2 + weno_sensitivity = 1.0d-40 + weno_smoothness_measurements => weno5_smoothness_measurements_js + weno_stencil_weights => weno5_stencil_weights_zcp + case ("weno5z+m", "weno5-z+m", "WENO5Z+M", "WENO5-Z+M") + name_rec = "5th order WENO-Z+M (Luo & Wu 2021)" + interfaces => interfaces_dir + reconstruct_states => reconstruct_weno5zpm + order = 5 + nghosts = max(nghosts, 4) + weno_scheme = .true. + weno_power = 2 + weno_sensitivity = 1.0d-40 + weno_smoothness_measurements => weno5_smoothness_measurements_js case ("crweno5z", "crweno5-z", "CRWENO5Z", "CRWENO5-Z") name_rec = "5th order Compact WENO-Z" interfaces => interfaces_dir reconstruct_states => reconstruct_crweno5z order = 5 nghosts = max(nghosts, 4) + weno_scheme = .true. + weno_power = 2 + weno_sensitivity = 1.0d-40 case ("crweno5yc", "crweno5-yc", "CRWENO5YC", "CRWENO5-YC") name_rec = "5th order Compact WENO-YC" interfaces => interfaces_dir reconstruct_states => reconstruct_crweno5yc order = 5 nghosts = max(nghosts, 4) + weno_scheme = .true. + weno_power = 2 + weno_sensitivity = 1.0d-40 case ("crweno5ns", "crweno5-ns", "CRWENO5NS", "CRWENO5-NS") name_rec = "5th order Compact WENO-NS" interfaces => interfaces_dir reconstruct_states => reconstruct_crweno5ns order = 5 nghosts = max(nghosts, 4) + weno_scheme = .true. + weno_power = 2 + weno_sensitivity = 1.0d-40 case ("mp5", "MP5") name_rec = "5th order Monotonicity Preserving" interfaces => interfaces_dir @@ -426,8 +536,9 @@ module interpolations write(msg,fmt) "'" // uppercase(trim(sreconstruction)) // "'", & new_line('a'), & "'TVD' 'LIMO3', 'PPM', " // & - "'WENO3', 'WENO5Z', 'WENO5YC', 'WENO5NS', " // & - "'CRWENO5Z', 'CRWENO5YC', 'CRWENO5NS', " // & + "'WENO3', 'WENO5JS', 'WENO5Z', 'WENO5YC', " // & + "'WENO5NS', 'WENO5Z+', 'WENO5Z+M', 'WENO5AT', " // & + "'CRWENO5Z', 'CRWENO5YC', 'CRWENO5NS', " // & "'MP5', 'MP7', 'MP9', " // & "'CRMP5', 'CRMP7', 'CRMP9', " // & "'OCMP5', 'OCMP7', 'OCMP9', " // & @@ -439,6 +550,36 @@ module interpolations end select +! process WENO related parameters +! + if (weno_scheme) then + +! get WENO related parameters +! + call get_parameter("weno_power" , weno_power ) + call get_parameter("weno_sensitivity", weno_sensitivity) + +! correct wrong values +! + if (weno_power < 1) then + if (verbose) & + call print_message(loc, & + "The parameter 'weno_power' has to be an integer number " // & + "equal to or greater than 1. Setting it to the default " // & + "value of 1.") + weno_power = 1 + end if + if (weno_sensitivity <= 0.0d+00) then + if (verbose) & + call print_message(loc, & + "The parameter 'weno_sensitivity' has to be a real " // & + "number greater than 0. Setting it to the default " // & + "value of 1e-40.") + weno_sensitivity = 1.0d-40 + end if + + end if ! weno_scheme = .true. + ! select the TVD limiter ! select case(trim(tlimiter)) @@ -665,6 +806,10 @@ module interpolations call print_section(verbose, "Interpolations") call print_parameter(verbose, "reconstruction" , name_rec ) + if (weno_scheme) then + call print_parameter(verbose, "WENO power" , weno_power) + call print_parameter(verbose, "WENO sensitivity", weno_sensitivity) + end if call print_parameter(verbose, "TVD limiter" , name_tlim) call print_parameter(verbose, "prolongation limiter", name_plim) if (mlp) then @@ -2003,320 +2148,739 @@ module interpolations ! !=============================================================================== ! -! subroutine RECONSTRUCT_WENO5Z: -! ----------------------------- +! subroutine RECONSTRUCT_WENO5: +! ---------------------------- ! -! Subroutine reconstructs the interface states using the fifth order +! Subroutine reconstructs the interface states using the 5th order ! Explicit Weighted Essentially Non-Oscillatory (WENO5) method with -! stencil weights by Borges et al. (2008). +! the selected stencil weights. ! ! Arguments are described in subroutine reconstruct(). ! ! References: ! -! [1] Borges, R., Carmona, M., Costa, B., & Don, W.-S., +! [1] Jiang, Guang-Shan & Shu, Chi-Wang, +! "Efficient Implementation of Weighted ENO Schemes", +! Journal of Computational Physics, +! 1996, Volume 126, Issue 1, p. 202-228, +! https://doi.org/10.1006/jcph.1996.0130 +! +!=============================================================================== +! + subroutine reconstruct_weno5(p, h, f, fl, fr) + + implicit none + + logical , intent(in) :: p + real(kind=8) , intent(in) :: h + real(kind=8), dimension(:), intent(in) :: f + real(kind=8), dimension(:), intent(out) :: fl, fr + + integer :: n, i + real(kind=8) :: df + + real(kind=8), dimension(5) :: ff, fb + real(kind=8), dimension(3) :: af, ab, qf, qb, wf, wb + + real(kind=8), dimension(3, size(f)) :: alphas, betas + + real(kind=8), dimension(3), parameter :: d = [1.0d-01, 6.0d-01, 3.0d-01] + real(kind=8), dimension(3), parameter :: & + c1 = [ 2.0d+00,-7.0d+00, 1.1d+01] / 6.0d+00, & + c2 = [-1.0d+00, 5.0d+00, 2.0d+00] / 6.0d+00, & + c3 = [ 2.0d+00, 5.0d+00,-1.0d+00] / 6.0d+00 + +!------------------------------------------------------------------------------- +! +! calculate the smoothness measurements βₖ (eqs. 3.2-3.2 in [1]) +! + call weno_smoothness_measurements(f, betas) + +! calculate the stencil weights αₖ (eqs. 2.16 in [1]) +! + call weno_stencil_weights(h, f, betas, alphas) + +! apply the WENO5 interpolation with the precomputed weights +! + n = size(f) + + do i = 3, n - 2 + + ff = f(i-2:i+2) + fb = f(i+2:i-2:-1) + + af = alphas(1:3 ,i) + ab = alphas(3:1:-1,i) + + qf = [ sum(c1 * ff(1:3)), sum(c2 * ff(2:4)), sum(c3 * ff(3:5)) ] + qb = [ sum(c1 * fb(1:3)), sum(c2 * fb(2:4)), sum(c3 * fb(3:5)) ] + + wf = d * af + wf = wf / sum(wf) + + wb = d * ab + wb = wb / sum(wb) + + fl(i ) = sum(wf * qf) + fr(i-1) = sum(wb * qb) + + end do ! i = 3, n - 2 + +! update the interpolation of the first and last two points +! + fl(1) = 5.0d-01 * (f(1) + f(2)) + df = limiter_tvd(5.0d-01, f(2)-f(1), f(3)-f(2)) + fr(1) = f(2) - df + fl(2) = f(2) + df + df = limiter_tvd(5.0d-01, f(n)-f(n-1), f(n-2)-f(n-3)) + fr(n-2) = f(n-1) - df + fl(n-1) = f(n-1) + df + fr(n-1) = 5.0d-01 * (f(n-1) + f(n)) + fl(n) = f(n) + fr(n) = f(n) + +!------------------------------------------------------------------------------- +! + end subroutine reconstruct_weno5 +! +!=============================================================================== +! +! subroutine WENO5_SMOOTHNESS_MEASUREMENTS_JS: +! ------------------------------------------- +! +! Subroutine calculates the Jiang & Shu (1996) smoothness measurements for +! the 5th order WENO schemes. +! +! Arguments: +! +! q - the vector of cell-integrated values; +! betas - the array of smoothing measurements; +! +! References: +! +! [1] Jiang, Guang-Shan & Shu, Chi-Wang, +! "Efficient Implementation of Weighted ENO Schemes", +! Journal of Computational Physics, +! 1996, Volume 126, Issue 1, p. 202-228, +! https://doi.org/10.1006/jcph.1996.0130 +! +!=============================================================================== +! + subroutine weno5_smoothness_measurements_js(f, betas) + + implicit none + + real(kind=8), dimension(:) , intent(in) :: f + real(kind=8), dimension(:,:), intent( out) :: betas + + integer :: n + + real(kind=8), dimension(size(f)) :: df, df2 + +!------------------------------------------------------------------------------- +! + n = size(f) + + df(1:n-1) = f(2:n) - f(1:n-1) + df(n) = df(n-1) + df2(2:n) = 1.3d+01 * (df(2:n) - df(1:n-1))**2 / 1.2d+01 + + betas(1,3:n ) = df2(2:n-1) + 2.5d-01 * (3.0d+00 * df(2:n-1) - df(1:n-2))**2 + betas(2,2:n ) = df2(2:n ) + 2.5d-01 * ( df(2:n ) + df(1:n-1))**2 + betas(3,1:n-1) = df2(2:n ) + 2.5d-01 * (3.0d+00 * df(1:n-1) - df(2:n ))**2 + +!------------------------------------------------------------------------------- +! + end subroutine weno5_smoothness_measurements_js +! +!=============================================================================== +! +! subroutine WENO5_SMOOTHNESS_MEASUREMENTS_NS: +! ------------------------------------------- +! +! Subroutine calculates the Ha et al. (2013) smoothness measurements for +! the 5th order WENO schemes. +! +! Arguments: +! +! q - the vector of cell-integrated values; +! betas - the array of smoothing measurements; +! +! References: +! +! [1] Youngsoo Ha, Chang Ho Kim, Yeon Ju Lee, & Jungho Yoon, +! "An improved weighted essentially non-oscillatory scheme with a new +! smoothness indicator", +! Journal of Computational Physics, +! 2013, Volume 232, p. 68-86, +! https://doi.org/10.1016/j.jcp.2012.06.016 +! +!=============================================================================== +! + subroutine weno5_smoothness_measurements_ns(f, betas) + + implicit none + + real(kind=8), dimension(:) , intent(in) :: f + real(kind=8), dimension(:,:), intent( out) :: betas + + integer :: n + + real(kind=8), dimension(size(f)) :: df, df2 + +!------------------------------------------------------------------------------- +! + n = size(f) + + df(1:n-1) = f(2:n) - f(1:n-1) + df(n) = df(n-1) + df2(2:n) = 1.3d+01 * (df(2:n) - df(1:n-1))**2 / 1.2d+01 + + betas(1,3:n ) = df2(2:n-1) + 2.5d-01 * (3.0d+00 * df(2:n-1) - df(1:n-2))**2 + betas(2,2:n ) = df2(2:n ) + 2.5d-01 * ( df(2:n ) + df(1:n-1))**2 + betas(3,1:n-1) = df2(2:n ) + 2.5d-01 * (3.0d+00 * df(1:n-1) - df(2:n ))**2 + +!------------------------------------------------------------------------------- +! + end subroutine weno5_smoothness_measurements_ns +! +!=============================================================================== +! +! subroutine WENO5_STENCIL_WEIGHTS_JS: +! ----------------------------------- +! +! Subroutine calculates the Jiang & Shu (1996) stencil weigths for +! the 5th order WENO schemes. +! +! Arguments: +! +! h - the cell size; +! q - the vector of cell-integrated values; +! betas - the array of smoothing measurements; +! alphas - the array of stencil weights; +! +! References: +! +! [1] Jiang, Guang-Shan & Shu, Chi-Wang, +! "Efficient Implementation of Weighted ENO Schemes", +! Journal of Computational Physics, +! 1996, Volume 126, Issue 1, p. 202-228, +! https://doi.org/10.1006/jcph.1996.0130 +! +!=============================================================================== +! + subroutine weno5_stencil_weights_js(h, f, betas, alphas) + + implicit none + + real(kind=8) , intent(in) :: h + real(kind=8), dimension(:) , intent(in) :: f + real(kind=8), dimension(:,:), intent(in ) :: betas + real(kind=8), dimension(:,:), intent( out) :: alphas + +!------------------------------------------------------------------------------- +! + alphas(:,:) = 1.0d+00 / (betas(:,:) + weno_sensitivity)**weno_power + +!------------------------------------------------------------------------------- +! + end subroutine weno5_stencil_weights_js +! +!=============================================================================== +! +! subroutine WENO5_STENCIL_WEIGHTS_Z: +! ---------------------------------- +! +! Subroutine calculates the Borges eta al. (2008) stencil weigths for +! the 5th order WENO schemes. +! +! Arguments: +! +! h - the cell size; +! q - the vector of cell-integrated values; +! betas - the array of smoothing measurements; +! alphas - the array of stencil weights; +! +! References: +! +! [1] Rafael Borges, Monique Carmona, Bruno Costa, & Wai Sun Don, Borges, ! "An improved weighted essentially non-oscillatory scheme for -! hyperbolic conservation laws" +! hyperbolic conservation laws", ! Journal of Computational Physics, -! 2008, vol. 227, pp. 3191-3211, -! http://dx.doi.org/10.1016/j.jcp.2007.11.038 +! 2008, Volume 227, p. 3191-3211, +! https://doi.org/10.1016/j.jcp.2007.11.038 ! !=============================================================================== ! - subroutine reconstruct_weno5z(p, h, fc, fl, fr) + subroutine weno5_stencil_weights_z(h, f, betas, alphas) implicit none - logical , intent(in) :: p - real(kind=8) , intent(in) :: h - real(kind=8), dimension(:), intent(in) :: fc - real(kind=8), dimension(:), intent(out) :: fl, fr + real(kind=8) , intent(in) :: h + real(kind=8), dimension(:) , intent(in) :: f + real(kind=8), dimension(:,:), intent(in ) :: betas + real(kind=8), dimension(:,:), intent( out) :: alphas - integer :: n, i, im1, ip1, im2, ip2 - real(kind=8) :: bl, bc, br, tt, df - real(kind=8) :: al, ac, ar - real(kind=8) :: wl, wc, wr, ww - real(kind=8) :: ql, qc, qr - - real(kind=8), dimension(size(fc)) :: dfm, dfp, df2 - - real(kind=8), parameter :: c1 = 1.3d+01 / 1.2d+01, c2 = 2.5d-01 - - real(kind=8), parameter :: dl = 1.0d-01, dc = 6.0d-01, dr = 3.0d-01 - - real(kind=8), parameter :: a11 = 2.0d+00 / 6.0d+00 & - , a12 = - 7.0d+00 / 6.0d+00 & - , a13 = 1.1d+01 / 6.0d+00 - real(kind=8), parameter :: a21 = - 1.0d+00 / 6.0d+00 & - , a22 = 5.0d+00 / 6.0d+00 & - , a23 = 2.0d+00 / 6.0d+00 - real(kind=8), parameter :: a31 = 2.0d+00 / 6.0d+00 & - , a32 = 5.0d+00 / 6.0d+00 & - , a33 = - 1.0d+00 / 6.0d+00 + integer :: i + real(kind=8) :: tau !------------------------------------------------------------------------------- ! - n = size(fc) + do i = 1, size(f) + tau = abs(betas(1,i) - betas(3,i)) -! calculate the left and right derivatives -! - do i = 1, n - 1 - ip1 = i + 1 - dfp(i ) = fc(ip1) - fc(i) - dfm(ip1) = dfp(i) + alphas(:,i) = 1.0d+00 & + + (tau / (betas(:,i) + weno_sensitivity))**weno_power end do - dfm(1) = dfp(1) - dfp(n) = dfm(n) - -! calculate the absolute value of the second derivative -! - df2(:) = c1 * (dfp(:) - dfm(:))**2 - -! iterate along the vector -! - do i = 3, n - 2 - -! prepare neighbour indices -! - im2 = i - 2 - im1 = i - 1 - ip1 = i + 1 - ip2 = i + 2 - -! calculate βₖ (eqs. 9-11 in [1]) -! - bl = df2(im1) + c2 * (3.0d+00 * dfm(i ) - dfm(im1))**2 - bc = df2(i ) + c2 * ( dfp(i ) + dfm(i ))**2 - br = df2(ip1) + c2 * (3.0d+00 * dfp(i ) - dfp(ip1))**2 - -! calculate τ (below eq. 25 in [1]) -! - tt = abs(br - bl) - -! calculate αₖ (eq. 28 in [1]) -! - al = 1.0d+00 + tt / (bl + eps) - ac = 1.0d+00 + tt / (bc + eps) - ar = 1.0d+00 + tt / (br + eps) - -! calculate weights -! - wl = dl * al - wc = dc * ac - wr = dr * ar - ww = (wl + wr) + wc - wl = wl / ww - wr = wr / ww - wc = 1.0d+00 - (wl + wr) - -! calculate the interpolations of the left state -! - ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i ) - qc = a21 * fc(im1) + a22 * fc(i ) + a23 * fc(ip1) - qr = a31 * fc(i ) + a32 * fc(ip1) + a33 * fc(ip2) - -! calculate the left state -! - fl(i ) = (wl * ql + wr * qr) + wc * qc - -! normalize weights -! - wl = dl * ar - wc = dc * ac - wr = dr * al - ww = (wl + wr) + wc - wl = wl / ww - wr = wr / ww - wc = 1.0d+00 - (wl + wr) - -! calculate the interpolations of the right state -! - ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i ) - qc = a21 * fc(ip1) + a22 * fc(i ) + a23 * fc(im1) - qr = a31 * fc(i ) + a32 * fc(im1) + a33 * fc(im2) - -! calculate the right state -! - fr(im1) = (wl * ql + wr * qr) + wc * qc - - end do ! i = 3, n - 2 - -! update the interpolation of the first and last two points -! - fl(1) = 0.5d+00 * (fc(1) + fc(2)) - df = limiter_tvd(0.5d+00, dfm(2), dfp(2)) - fr(1) = fc(2) - df - fl(2) = fc(2) + df - i = n - 1 - df = limiter_tvd(0.5d+00, dfm(i), dfp(i)) - fr(i-1) = fc(i) - df - fl(i) = fc(i) + df - fr(i) = 0.5d+00 * (fc(i) + fc(n)) - fl(n) = fc(n) - fr(n) = fc(n) !------------------------------------------------------------------------------- ! - end subroutine reconstruct_weno5z + end subroutine weno5_stencil_weights_z ! !=============================================================================== ! -! subroutine RECONSTRUCT_WENO5YC: -! ------------------------------ +! subroutine WENO5_STENCIL_WEIGHTS_YC: +! ----------------------------------- ! -! Subroutine reconstructs the interface states using the fifth order +! Subroutine calculates the Yamaleev & Carpenter (2009) stencil weigths for +! the 5th order WENO schemes. +! +! Arguments: +! +! h - the cell size; +! q - the vector of cell-integrated values; +! betas - the array of smoothing measurements; +! alphas - the array of stencil weights; +! +! References: +! +! [1] Nail K. Yamaleev & Mark H. Carpenter, +! "A systematic methodology for constructing high-order energy stable +! WENO schemes", +! Journal of Computational Physics, +! 2009, Volume 228, p. 4248-4272, +! https://doi.org/10.1016/j.jcp.2009.03.002 +! +!=============================================================================== +! + subroutine weno5_stencil_weights_yc(h, f, betas, alphas) + + implicit none + + real(kind=8) , intent(in) :: h + real(kind=8), dimension(:) , intent(in) :: f + real(kind=8), dimension(:,:), intent(in ) :: betas + real(kind=8), dimension(:,:), intent( out) :: alphas + + integer :: i + real(kind=8) :: tau + +!------------------------------------------------------------------------------- +! + + do i = 3, size(f) - 2 + tau = (6.0d+00 * f(i) - 4.0d+00 * (f(i-1) + f(i+1)) & + + (f(i-2) + f(i+2)))**2 + + alphas(:,i) = 1.0d+00 & + + (tau / (betas(:,i) + weno_sensitivity))**weno_power + end do + +!------------------------------------------------------------------------------- +! + end subroutine weno5_stencil_weights_yc +! +!=============================================================================== +! +! subroutine WENO5_STENCIL_WEIGHTS_ZP: +! ----------------------------------- +! +! Subroutine calculates the Acker eta al. (2016) stencil weigths for +! the 5th order WENO schemes. +! +! Arguments: +! +! h - the cell size; +! q - the vector of cell-integrated values; +! betas - the array of smoothing measurements; +! alphas - the array of stencil weights; +! +! References: +! +! [1] F. Acker, R. B. de R. Borges, & B. Costa, +! "An improved WENO-Z scheme", +! Journal of Computational Physics, +! 2016, Volume 313, p. 726-753, +! https://doi.org/10.1016/j.jcp.2016.01.038 +! +!=============================================================================== +! + subroutine weno5_stencil_weights_zp(h, f, betas, alphas) + + implicit none + + real(kind=8) , intent(in) :: h + real(kind=8), dimension(:) , intent(in) :: f + real(kind=8), dimension(:,:), intent(in ) :: betas + real(kind=8), dimension(:,:), intent( out) :: alphas + + integer :: i + real(kind=8) :: tau, lambda + + real(kind=8), dimension(3) :: xi + + real(kind=8), parameter :: two_thirds = 2.0d+00 / 3.0d+00 + +!------------------------------------------------------------------------------- +! + lambda = h**two_thirds + + do i = 1, size(f) + tau = abs(betas(1,i) - betas(3,i)) + + xi(:) = (tau + weno_sensitivity) / (betas(:,i) + weno_sensitivity) + + alphas(:,i) = 1.0d+00 + xi(:)**weno_power + lambda / xi(:) + end do + +!------------------------------------------------------------------------------- +! + end subroutine weno5_stencil_weights_zp +! +!=============================================================================== +! +! subroutine WENO5_STENCIL_WEIGHTS_AT: +! ----------------------------------- +! +! Subroutine calculates the Huang et al (2022) stencil weigths for +! the 5th order WENO schemes. +! +! Arguments: +! +! h - the cell size; +! q - the vector of cell-integrated values; +! betas - the array of smoothing measurements; +! alphas - the array of stencil weights; +! +! References: +! +! [1] Ziquan Huang, Shichao Zheng, Dongfang Wang, & Xiaogang Deng, +! "A New e-Adaptive Algorithm for Improving Weighted Compact +! Nonlinear Scheme with Applications", +! Aerospace, +! 2022, Volume 9, Issue 369, +! https://doi.org/10.3390/aerospace9070369 +! +!=============================================================================== +! + subroutine weno5_stencil_weights_at(h, f, betas, alphas) + + implicit none + + real(kind=8) , intent(in) :: h + real(kind=8), dimension(:) , intent(in) :: f + real(kind=8), dimension(:,:), intent(in ) :: betas + real(kind=8), dimension(:,:), intent( out) :: alphas + + integer :: i + real(kind=8) :: tau, rho, eps_adp + + real(kind=8), parameter :: small = sqrt(tiny(1.0d+00)) + +!------------------------------------------------------------------------------- +! + do i = 1, size(f) + tau = abs(betas(1,i) - betas(3,i)) + rho = (tau + weno_sensitivity) / (maxval(betas(:,i)) + weno_sensitivity) + + eps_adp = minval(betas(:,i)) / rho**weno_power + + alphas(:,i) = 1.0d+00 / ((betas(:,i) + eps_adp)**2 + weno_sensitivity) + end do + +!------------------------------------------------------------------------------- +! + end subroutine weno5_stencil_weights_at +! +!=============================================================================== +! +! subroutine WENO5_STENCIL_WEIGHTS_ZC: +! ----------------------------------- +! +! Subroutine calculates the Baretto et al. (2023) stencil weigths for +! the 5th order WENO schemes. +! +! Arguments: +! +! h - the cell size; +! q - the vector of cell-integrated values; +! betas - the array of smoothing measurements; +! alphas - the array of stencil weights; +! +! References: +! +! [1] Daniel Barreto, Rafael B. de R. Borges, Bruno Costa, +! & Silvaneo dos Santos, +! "New Weighting Strategies for WENO Schemes", +! arXiv, +! 2023, arXiv:2311.09332, +! https://doi.org/10.48550/arXiv.2311.09332 +! +!=============================================================================== +! + subroutine weno5_stencil_weights_zc(h, f, betas, alphas) + + implicit none + + real(kind=8) , intent(in) :: h + real(kind=8), dimension(:) , intent(in) :: f + real(kind=8), dimension(:,:), intent(in ) :: betas + real(kind=8), dimension(:,:), intent( out) :: alphas + + integer :: i + real(kind=8) :: tau, avg, den + + real(kind=8), dimension(3) :: xi + + real(kind=8), dimension(3), parameter :: c = [7.5d-01, 1.5d+00, 7.5d-01] + +!------------------------------------------------------------------------------- +! + do i = 1, size(f) + tau = abs(betas(1,i) - betas(3,i)) + avg = sum(betas(:,i)) / 3.0d+00 + den = tau + avg + weno_sensitivity + + xi(:) = tau**2 / ((betas(:,i) + weno_sensitivity) * den) + + alphas(:,i) = 1.0d+00 + c * xi(:)**weno_power + end do + +!------------------------------------------------------------------------------- +! + end subroutine weno5_stencil_weights_zc +! +!=============================================================================== +! +! subroutine WENO5_STENCIL_WEIGHTS_ZCP: +! ------------------------------------ +! +! Subroutine calculates the Baretto et al. (2023) stencil weigths for +! the 5th order WENO schemes. +! +! Arguments: +! +! h - the cell size; +! q - the vector of cell-integrated values; +! betas - the array of smoothing measurements; +! alphas - the array of stencil weights; +! +! References: +! +! [1] Daniel Barreto, Rafael B. de R. Borges, Bruno Costa, +! & Silvaneo dos Santos, +! "New Weighting Strategies for WENO Schemes", +! arXiv, +! 2023, arXiv:2311.09332, +! https://doi.org/10.48550/arXiv.2311.09332 +! +!=============================================================================== +! + subroutine weno5_stencil_weights_zcp(h, f, betas, alphas) + + implicit none + + real(kind=8) , intent(in) :: h + real(kind=8), dimension(:) , intent(in) :: f + real(kind=8), dimension(:,:), intent(in ) :: betas + real(kind=8), dimension(:,:), intent( out) :: alphas + + integer :: i + real(kind=8) :: tau, avg, den + + real(kind=8), dimension(3) :: xi + + real(kind=8), dimension(3), parameter :: c = [7.5d-01, 1.5d+00, 7.5d-01] + +!------------------------------------------------------------------------------- +! + do i = 1, size(f) + tau = abs(betas(1,i) - betas(3,i)) + avg = sum(betas(:,i)) / 3.0d+00 + den = tau + avg + weno_sensitivity + + xi(:) = tau**2 / ((betas(:,i) + weno_sensitivity) * den) + + alphas(:,i) = 1.0d+00 + c * xi(:)**weno_power + betas(:,i) / den + end do + +!------------------------------------------------------------------------------- +! + end subroutine weno5_stencil_weights_zcp +! +!=============================================================================== +! +! subroutine RECONSTRUCT_WENO5ZPM: +! ------------------------------- +! +! Subroutine reconstructs the interface states using the 5th order ! Explicit Weighted Essentially Non-Oscillatory (WENO5) method with -! stencil weights by Yamaleev & Carpenter (2009). +! the Luo & Wu (2021) stencil weigths. ! ! Arguments are described in subroutine reconstruct(). ! ! References: ! -! [1] Yamaleev, N. K. & Carpenter, H. C., -! "A Systematic Methodology for Constructing High-Order Energy Stable -! WENO Schemes" +! [1] Jiang, Guang-Shan & Shu, Chi-Wang, +! "Efficient Implementation of Weighted ENO Schemes", ! Journal of Computational Physics, -! 2009, vol. 228, pp. 4248-4272, -! http://dx.doi.org/10.1016/j.jcp.2009.03.002 +! 1996, Volume 126, Issue 1, p. 202-228, +! https://doi.org/10.1006/jcph.1996.0130 +! [2] Xin Luo & Song-ping Wu, +! "An improved WENO-Z+ scheme for solving hyperbolic conservation laws", +! Journal of Computational Physics, +! 2021, Volume 445, Issue 110608, +! https://doi.org/10.1016/j.jcp.2021.110608 ! !=============================================================================== ! - subroutine reconstruct_weno5yc(p, h, fc, fl, fr) + subroutine reconstruct_weno5zpm(p, h, f, fl, fr) implicit none logical , intent(in) :: p real(kind=8) , intent(in) :: h - real(kind=8), dimension(:), intent(in) :: fc + real(kind=8), dimension(:), intent(in) :: f real(kind=8), dimension(:), intent(out) :: fl, fr - integer :: n, i, im1, ip1, im2, ip2 - real(kind=8) :: bl, bc, br, tt, df - real(kind=8) :: al, ac, ar - real(kind=8) :: wl, wc, wr, ww - real(kind=8) :: ql, qc, qr + integer :: n, i + real(kind=8) :: df - real(kind=8), dimension(size(fc)) :: dfm, dfp, df2 + real(kind=8), dimension(5) :: ff, fb + real(kind=8), dimension(3) :: af, ab, qf, qb, wf, wb - real(kind=8), parameter :: c1 = 1.3d+01 / 1.2d+01, c2 = 2.5d-01 + real(kind=8), dimension(3, size(f)) :: alphas_f, alphas_b, betas - real(kind=8), parameter :: dl = 1.0d-01, dc = 6.0d-01, dr = 3.0d-01 - - real(kind=8), parameter :: a11 = 2.0d+00 / 6.0d+00 & - , a12 = - 7.0d+00 / 6.0d+00 & - , a13 = 1.1d+01 / 6.0d+00 - real(kind=8), parameter :: a21 = - 1.0d+00 / 6.0d+00 & - , a22 = 5.0d+00 / 6.0d+00 & - , a23 = 2.0d+00 / 6.0d+00 - real(kind=8), parameter :: a31 = 2.0d+00 / 6.0d+00 & - , a32 = 5.0d+00 / 6.0d+00 & - , a33 = - 1.0d+00 / 6.0d+00 + real(kind=8), dimension(3), parameter :: d = [1.0d-01, 6.0d-01, 3.0d-01] + real(kind=8), dimension(3), parameter :: & + c1 = [ 2.0d+00,-7.0d+00, 1.1d+01] / 6.0d+00, & + c2 = [-1.0d+00, 5.0d+00, 2.0d+00] / 6.0d+00, & + c3 = [ 2.0d+00, 5.0d+00,-1.0d+00] / 6.0d+00 !------------------------------------------------------------------------------- ! - n = size(fc) - -! calculate the left and right derivatives +! calculate the smoothness measurements βₖ (eqs. 3.2-3.2 in [1]) ! - do i = 1, n - 1 - ip1 = i + 1 - dfp(i ) = fc(ip1) - fc(i) - dfm(ip1) = dfp(i) - end do - dfm(1) = dfp(1) - dfp(n) = dfm(n) + call weno_smoothness_measurements(f, betas) -! calculate the absolute value of the second derivative +! calculate the stencil weights αₖ (eqs. 2.16 in [1]) ! - df2(:) = c1 * (dfp(:) - dfm(:))**2 + call weno5_stencil_weights_zpm(h, f, d(1:3: 1), betas, alphas_f) + call weno5_stencil_weights_zpm(h, f, d(3:1:-1), betas, alphas_b) -! iterate along the vector +! apply the WENO5 interpolation with the precomputed weights ! + n = size(f) + do i = 3, n - 2 -! prepare neighbour indices -! - im2 = i - 2 - im1 = i - 1 - ip1 = i + 1 - ip2 = i + 2 + ff = f(i-2:i+2) + fb = f(i+2:i-2:-1) -! calculate βₖ (eq. 19 in [1]) -! - bl = df2(im1) + c2 * (3.0d+00 * dfm(i ) - dfm(im1))**2 - bc = df2(i ) + c2 * ( dfp(i ) + dfm(i ))**2 - br = df2(ip1) + c2 * (3.0d+00 * dfp(i ) - dfp(ip1))**2 + af = alphas_f(1:3 ,i) + ab = alphas_b(3:1:-1,i) -! calculate τ (below eq. 20 in [1]) -! - tt = (6.0d+00 * fc(i) - 4.0d+00 * (fc(im1) + fc(ip1)) & - + (fc(im2) + fc(ip2)))**2 + qf = [ sum(c1 * ff(1:3)), sum(c2 * ff(2:4)), sum(c3 * ff(3:5)) ] + qb = [ sum(c1 * fb(1:3)), sum(c2 * fb(2:4)), sum(c3 * fb(3:5)) ] -! calculate αₖ (eqs. 18 or 58 in [1]) -! - al = 1.0d+00 + tt / (bl + eps) - ac = 1.0d+00 + tt / (bc + eps) - ar = 1.0d+00 + tt / (br + eps) + wf = d * af + wf = wf / sum(wf) -! calculate weights -! - wl = dl * al - wc = dc * ac - wr = dr * ar - ww = (wl + wr) + wc - wl = wl / ww - wr = wr / ww - wc = 1.0d+00 - (wl + wr) + wb = d * ab + wb = wb / sum(wb) -! calculate the interpolations of the left state -! - ql = a11 * fc(im2) + a12 * fc(im1) + a13 * fc(i ) - qc = a21 * fc(im1) + a22 * fc(i ) + a23 * fc(ip1) - qr = a31 * fc(i ) + a32 * fc(ip1) + a33 * fc(ip2) - -! calculate the left state -! - fl(i ) = (wl * ql + wr * qr) + wc * qc - -! normalize weights -! - wl = dl * ar - wc = dc * ac - wr = dr * al - ww = (wl + wr) + wc - wl = wl / ww - wr = wr / ww - wc = 1.0d+00 - (wl + wr) - -! calculate the interpolations of the right state -! - ql = a11 * fc(ip2) + a12 * fc(ip1) + a13 * fc(i ) - qc = a21 * fc(ip1) + a22 * fc(i ) + a23 * fc(im1) - qr = a31 * fc(i ) + a32 * fc(im1) + a33 * fc(im2) - -! calculate the right state -! - fr(im1) = (wl * ql + wr * qr) + wc * qc + fl(i ) = sum(wf * qf) + fr(i-1) = sum(wb * qb) end do ! i = 3, n - 2 ! update the interpolation of the first and last two points ! - fl(1) = 0.5d+00 * (fc(1) + fc(2)) - df = limiter_tvd(0.5d+00, dfm(2), dfp(2)) - fr(1) = fc(2) - df - fl(2) = fc(2) + df - i = n - 1 - df = limiter_tvd(0.5d+00, dfm(i), dfp(i)) - fr(i-1) = fc(i) - df - fl(i) = fc(i) + df - fr(i) = 0.5d+00 * (fc(i) + fc(n)) - fl(n) = fc(n) - fr(n) = fc(n) + fl(1) = 5.0d-01 * (f(1) + f(2)) + df = limiter_tvd(5.0d-01, f(2)-f(1), f(3)-f(2)) + fr(1) = f(2) - df + fl(2) = f(2) + df + df = limiter_tvd(5.0d-01, f(n)-f(n-1), f(n-2)-f(n-3)) + fr(n-2) = f(n-1) - df + fl(n-1) = f(n-1) + df + fr(n-1) = 5.0d-01 * (f(n-1) + f(n)) + fl(n) = f(n) + fr(n) = f(n) !------------------------------------------------------------------------------- ! - end subroutine reconstruct_weno5yc + end subroutine reconstruct_weno5zpm +! +!=============================================================================== +! +! subroutine WENO5_STENCIL_WEIGHTS_ZPM: +! ------------------------------------ +! +! Subroutine calculates the Luo & Wu (2021) stencil weigths for +! the 5th order WENO schemes. +! +! Arguments: +! +! h - the cell size; +! q - the vector of cell-integrated values; +! d - the weights; +! betas - the array of smoothing measurements; +! alphas - the array of stencil weights; +! +! References: +! +! [1] Xin Luo & Song-ping Wu, +! "An improved WENO-Z+ scheme for solving hyperbolic conservation laws", +! Journal of Computational Physics, +! 2021, Volume 445, Issue 110608, +! https://doi.org/10.1016/j.jcp.2021.110608 +! +!=============================================================================== +! + subroutine weno5_stencil_weights_zpm(h, f, d, betas, alphas) + + implicit none + + real(kind=8) , intent(in) :: h + real(kind=8), dimension(:) , intent(in) :: f + real(kind=8), dimension(:) , intent(in) :: d + real(kind=8), dimension(:,:), intent(in ) :: betas + real(kind=8), dimension(:,:), intent( out) :: alphas + + integer :: i + real(kind=8) :: tau, lambda, z + + real(kind=8), dimension(3) :: xi + + real(kind=8), parameter :: two_thirds = 2.0d+00 / 3.0d+00 + +!------------------------------------------------------------------------------- +! + lambda = h**two_thirds + + do i = 1, size(f) + tau = abs(betas(1,i) - betas(3,i)) + + xi(:) = (tau + weno_sensitivity) / (betas(:,i) + weno_sensitivity) + + z = (1.0d+00 + minval(xi)**2) / (1.0d+00 + sum(d * xi**2)) + + alphas(:,i) = 1.0d+00 + xi(:)**weno_power + lambda * z / sqrt(xi(:)) + end do + +!------------------------------------------------------------------------------- +! + end subroutine weno5_stencil_weights_zpm ! !=============================================================================== ! @@ -2331,12 +2895,12 @@ module interpolations ! ! References: ! -! [1] Ha, Y., Kim, C. H., Lee, Y. J., & Yoon, J., +! [1] Youngsoo Ha, Chang Ho Kim, Yeon Ju Lee, & Jungho Yoon, ! "An improved weighted essentially non-oscillatory scheme with a new ! smoothness indicator", ! Journal of Computational Physics, -! 2013, vol. 232, pp. 68-86 -! http://dx.doi.org/10.1016/j.jcp.2012.06.016 +! 2013, Volume 232, p. 68-86, +! https://doi.org/10.1016/j.jcp.2012.06.016 ! !=============================================================================== !