diff --git a/sources/interpolations.F90 b/sources/interpolations.F90 index 08e637f..de1f420 100644 --- a/sources/interpolations.F90 +++ b/sources/interpolations.F90 @@ -264,85 +264,103 @@ module interpolations interfaces => interfaces_dir reconstruct_states => reconstruct_limo3 order = 3 - eps = max(1.0d-12, eps) + eps = max(1.0d-12, eps) case ("ppm", "PPM") name_rec = "3rd order PPM" interfaces => interfaces_dir reconstruct_states => reconstruct_ppm order = 3 - nghosts = max(nghosts, 4) + nghosts = max(nghosts, 4) case ("weno5z", "weno5-z", "WENO5Z", "WENO5-Z") name_rec = "5th order WENO-Z (Borges et al. 2008)" interfaces => interfaces_dir reconstruct_states => reconstruct_weno5z order = 5 - nghosts = max(nghosts, 4) + nghosts = max(nghosts, 4) case ("weno5yc", "weno5-yc", "WENO5YC", "WENO5-YC") name_rec = "5th order WENO-YC (Yamaleev & Carpenter 2009)" interfaces => interfaces_dir reconstruct_states => reconstruct_weno5yc order = 5 - nghosts = max(nghosts, 4) + nghosts = max(nghosts, 4) 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) + nghosts = max(nghosts, 4) 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) + nghosts = max(nghosts, 4) 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) + nghosts = max(nghosts, 4) 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) + nghosts = max(nghosts, 4) case ("mp5", "MP5") name_rec = "5th order Monotonicity Preserving" interfaces => interfaces_dir reconstruct_states => reconstruct_mp5 order = 5 - nghosts = max(nghosts, 4) + nghosts = max(nghosts, 4) case ("mp7", "MP7") name_rec = "7th order Monotonicity Preserving" interfaces => interfaces_dir reconstruct_states => reconstruct_mp7 order = 7 - nghosts = max(nghosts, 4) + nghosts = max(nghosts, 4) case ("mp9", "MP9") name_rec = "9th order Monotonicity Preserving" interfaces => interfaces_dir reconstruct_states => reconstruct_mp9 order = 9 - nghosts = max(nghosts, 6) + nghosts = max(nghosts, 6) + case ("mp5ld", "MP5LD") + name_rec = "5th order Low-Dissipation Monotonicity Preserving" + interfaces => interfaces_dir + reconstruct_states => reconstruct_mp5ld + order = 5 + nghosts = max(nghosts, 4) + case ("mp7ld", "MP7LD") + name_rec = "7th order Low-Dissipation Monotonicity Preserving" + interfaces => interfaces_dir + reconstruct_states => reconstruct_mp7ld + order = 7 + nghosts = max(nghosts, 4) + case ("mp9ld", "MP9LD") + name_rec = "9th order Low-Dissipation Monotonicity Preserving" + interfaces => interfaces_dir + reconstruct_states => reconstruct_mp9ld + order = 9 + nghosts = max(nghosts, 6) case ("crmp5", "CRMP5") name_rec = "5th order Compact Monotonicity Preserving" interfaces => interfaces_dir reconstruct_states => reconstruct_crmp5 order = 5 - nghosts = max(nghosts, 4) + nghosts = max(nghosts, 4) case ("crmp5l", "crmp5ld", "CRMP5L", "CRMP5LD") name_rec = "5th order Low-Dissipation Compact Monotonicity Preserving" interfaces => interfaces_dir reconstruct_states => reconstruct_crmp5ld order = 5 - nghosts = max(nghosts, 4) + nghosts = max(nghosts, 4) case ("crmp7", "CRMP7") name_rec = "7th order Compact Monotonicity Preserving" interfaces => interfaces_dir reconstruct_states => reconstruct_crmp7 order = 7 - nghosts = max(nghosts, 4) + nghosts = max(nghosts, 4) case ("gp", "GP") write(stmp, '(f16.1)') sgp write(name_rec, '("Gaussian Process (",i1,"-point, δ=",a,")")') ngp & @@ -4165,6 +4183,417 @@ module interpolations ! !=============================================================================== ! +! subroutine RECONSTRUCT_MP5LD: +! ---------------------------- +! +! Subroutine reconstructs the interface states using the fifth order +! low dissipation Monotonicity Preserving (MP) method. +! +! Arguments are described in subroutine reconstruct(). +! +! References: +! +! [1] Suresh, A. & Huynh, H. T., +! "Accurate Monotonicity-Preserving Schemes with Runge-Kutta +! Time Stepping" +! Journal on Computational Physics, +! 1997, vol. 136, pp. 83-99, +! http://dx.doi.org/10.1006/jcph.1997.5745 +! [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen, +! "A 5th order monotonicity-preserving upwind compact difference +! scheme", +! Science China Physics, Mechanics and Astronomy, +! Volume 54, Issue 3, pp. 511-522, +! http://dx.doi.org/10.1007/s11433-010-4220-x +! +!=============================================================================== +! + subroutine reconstruct_mp5ld(h, fc, fl, fr) + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + real(kind=8) , intent(in) :: h + real(kind=8), dimension(:), intent(in) :: fc + real(kind=8), dimension(:), intent(out) :: fl, fr + +! local variables +! + integer :: n, i + +! local arrays for derivatives +! + real(kind=8), dimension(size(fc)) :: fi + real(kind=8), dimension(size(fc)) :: u + +! local parameters +! + real(kind=8), dimension(6), parameter :: & + ce5 = [ 1.20d+01,-8.10d+01, 3.09d+02, & + 2.09d+02,-3.10d+01, 2.00d+00 ] / 4.2d+02 + real(kind=8), dimension(3), parameter :: & + ce3 = [-1.00d+00, 5.00d+00, 2.00d+00 ] / 6.0d+00 + real(kind=8), dimension(2), parameter :: & + ce2 = [ 5.0d-01, 5.0d-01 ] +! +!------------------------------------------------------------------------------- +! +! get the input vector length +! + n = size(fc) + +!! === left-side interpolation === +!! +! reconstruct the interface state using the 5th order interpolation +! + do i = 3, n - 3 + u(i) = sum(ce5(:) * fc(i-2:i+3)) + end do + +! interpolate the interface state of the ghost zones using the interpolations +! of lower orders +! + u( 1) = sum(ce2(:) * fc( 1: 2)) + u( 2) = sum(ce3(:) * fc( 1: 3)) + u(n-2) = sum(ce3(:) * fc(n-3:n-1)) + u(n-1) = sum(ce3(:) * fc(n-2:n )) + u(n ) = fc(n ) + +! apply the monotonicity preserving limiting +! + call mp_limiting(fc(:), u(:)) + +! copy the interpolation to the respective vector +! + fl(1:n) = u(1:n) + +!! === right-side interpolation === +!! +! invert the cell-centered value vector +! + fi(1:n) = fc(n:1:-1) + +! reconstruct the interface state using the 5th order interpolation +! + do i = 3, n - 3 + u(i) = sum(ce5(:) * fi(i-2:i+3)) + end do + +! interpolate the interface state of the ghost zones using the interpolations +! of lower orders +! + u( 1) = sum(ce2(:) * fi( 1: 2)) + u( 2) = sum(ce3(:) * fi( 1: 3)) + u(n-2) = sum(ce3(:) * fi(n-3:n-1)) + u(n-1) = sum(ce3(:) * fi(n-2:n )) + u(n ) = fi(n ) + +! apply the monotonicity preserving limiting +! + call mp_limiting(fi(:), u(:)) + +! copy the interpolation to the respective vector +! + fr(1:n-1) = u(n-1:1:-1) + +! update the interpolation of the first and last points +! + i = n - 1 + fl(1) = 0.5d+00 * (fc(1) + fc(2)) + fr(i) = 0.5d+00 * (fc(i) + fc(n)) + fl(n) = fc(n) + fr(n) = fc(n) + +!------------------------------------------------------------------------------- +! + end subroutine reconstruct_mp5ld +! +!=============================================================================== +! +! subroutine RECONSTRUCT_MP7LD: +! ---------------------------- +! +! Subroutine reconstructs the interface states using the seventh order +! low dissipation Monotonicity Preserving (MP) method. +! +! Arguments are described in subroutine reconstruct(). +! +! References: +! +! [1] Suresh, A. & Huynh, H. T., +! "Accurate Monotonicity-Preserving Schemes with Runge-Kutta +! Time Stepping" +! Journal on Computational Physics, +! 1997, vol. 136, pp. 83-99, +! http://dx.doi.org/10.1006/jcph.1997.5745 +! [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen, +! "A 5th order monotonicity-preserving upwind compact difference +! scheme", +! Science China Physics, Mechanics and Astronomy, +! Volume 54, Issue 3, pp. 511-522, +! http://dx.doi.org/10.1007/s11433-010-4220-x +! +!=============================================================================== +! + subroutine reconstruct_mp7ld(h, fc, fl, fr) + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + real(kind=8) , intent(in) :: h + real(kind=8), dimension(:), intent(in) :: fc + real(kind=8), dimension(:), intent(out) :: fl, fr + +! local variables +! + integer :: n, i + +! local arrays for derivatives +! + real(kind=8), dimension(size(fc)) :: fi + real(kind=8), dimension(size(fc)) :: u + +! local parameters +! + real(kind=8), dimension(8), parameter :: & + ce7 = [-8.00d+00, 6.80d+01,-2.82d+02, 9.22d+02, & + 6.77d+02,-1.35d+02, 1.9d+01, -1.0d+00 ] / 1.26d+03 + real(kind=8), dimension(5), parameter :: & + ce5 = [ 2.0d+00,-1.3d+01, 4.7d+01 & + , 2.7d+01,-3.0d+00 ] / 6.0d+01 + real(kind=8), dimension(3), parameter :: & + ce3 = [-1.0d+00, 5.0d+00, 2.0d+00 ] / 6.0d+00 + real(kind=8), dimension(2), parameter :: & + ce2 = [ 5.0d-01, 5.0d-01 ] +! +!------------------------------------------------------------------------------- +! +! get the input vector length +! + n = size(fc) + +!! === left-side interpolation === +!! +! reconstruct the interface state using the 5th order interpolation +! + do i = 4, n - 4 + u(i) = sum(ce7(:) * fc(i-3:i+4)) + end do + +! interpolate the interface state of the ghost zones using the interpolations +! of lower orders +! + u( 1) = sum(ce2(:) * fc( 1: 2)) + u( 2) = sum(ce3(:) * fc( 1: 3)) + u( 3) = sum(ce5(:) * fc( 1: 5)) + u(n-3) = sum(ce5(:) * fc(n-5:n-1)) + u(n-2) = sum(ce5(:) * fc(n-4:n )) + u(n-1) = sum(ce3(:) * fc(n-2:n )) + u(n ) = fc(n ) + +! apply the monotonicity preserving limiting +! + call mp_limiting(fc(:), u(:)) + +! copy the interpolation to the respective vector +! + fl(1:n) = u(1:n) + +!! === right-side interpolation === +!! +! invert the cell-centered value vector +! + fi(1:n) = fc(n:1:-1) + +! reconstruct the interface state using the 5th order interpolation +! + do i = 4, n - 4 + u(i) = sum(ce7(:) * fi(i-3:i+4)) + end do + +! interpolate the interface state of the ghost zones using the interpolations +! of lower orders +! + u( 1) = sum(ce2(:) * fi( 1: 2)) + u( 2) = sum(ce3(:) * fi( 1: 3)) + u( 3) = sum(ce5(:) * fi( 1: 5)) + u(n-3) = sum(ce5(:) * fi(n-5:n-1)) + u(n-2) = sum(ce5(:) * fi(n-4:n )) + u(n-1) = sum(ce3(:) * fi(n-2:n )) + u(n ) = fi(n ) + +! apply the monotonicity preserving limiting +! + call mp_limiting(fi(:), u(:)) + +! copy the interpolation to the respective vector +! + fr(1:n-1) = u(n-1:1:-1) + +! update the interpolation of the first and last points +! + i = n - 1 + fl(1) = 0.5d+00 * (fc(1) + fc(2)) + fr(i) = 0.5d+00 * (fc(i) + fc(n)) + fl(n) = fc(n) + fr(n) = fc(n) + +!------------------------------------------------------------------------------- +! + end subroutine reconstruct_mp7ld +! +!=============================================================================== +! +! subroutine RECONSTRUCT_MP9LD: +! ---------------------------- +! +! Subroutine reconstructs the interface states using the ninth order +! low dissipation Monotonicity Preserving (MP) method. +! +! Arguments are described in subroutine reconstruct(). +! +! References: +! +! [1] Suresh, A. & Huynh, H. T., +! "Accurate Monotonicity-Preserving Schemes with Runge-Kutta +! Time Stepping" +! Journal on Computational Physics, +! 1997, vol. 136, pp. 83-99, +! http://dx.doi.org/10.1006/jcph.1997.5745 +! [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen, +! "A 5th order monotonicity-preserving upwind compact difference +! scheme", +! Science China Physics, Mechanics and Astronomy, +! Volume 54, Issue 3, pp. 511-522, +! http://dx.doi.org/10.1007/s11433-010-4220-x +! +!=============================================================================== +! + subroutine reconstruct_mp9ld(h, fc, fl, fr) + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + real(kind=8) , intent(in) :: h + real(kind=8), dimension(:), intent(in) :: fc + real(kind=8), dimension(:), intent(out) :: fl, fr + +! local variables +! + integer :: n, i + +! local arrays for derivatives +! + real(kind=8), dimension(size(fc)) :: fi + real(kind=8), dimension(size(fc)) :: u + +! local parameters +! + real(kind=8), dimension(10), parameter :: & + ce9 = [ 4.0000d+01,-4.1500d+02, 2.0450d+03,-6.7150d+03, & + 2.0165d+04, 1.5629d+04,-3.6910d+03, 7.4900d+02, & + -9.1000d+01, 4.0000d+00 ] / 2.7720d+04 + real(kind=8), dimension(7), parameter :: & + ce7 = [-3.000d+00, 2.500d+01,-1.010d+02, 3.190d+02, 2.140d+02, & + -3.800d+01, 4.000d+00 ] / 4.200d+02 + real(kind=8), dimension(5), parameter :: & + ce5 = [ 2.0d+00,-1.3d+01, 4.7d+01, 2.7d+01,-3.0d+00 ] / 6.0d+01 + real(kind=8), dimension(3), parameter :: & + ce3 = [-1.0d+00, 5.0d+00, 2.0d+00 ] / 6.0d+00 + real(kind=8), dimension(2), parameter :: & + ce2 = [ 5.0d-01, 5.0d-01 ] +! +!------------------------------------------------------------------------------- +! +! get the input vector length +! + n = size(fc) + +!! === left-side interpolation === +!! +! reconstruct the interface state using the 9th order interpolation +! + do i = 5, n - 5 + u(i) = sum(ce9(:) * fc(i-4:i+5)) + end do + +! interpolate the interface state of the ghost zones using the interpolations +! of lower orders +! + u( 1) = sum(ce2(:) * fc( 1: 2)) + u( 2) = sum(ce3(:) * fc( 1: 3)) + u( 3) = sum(ce5(:) * fc( 1: 5)) + u( 4) = sum(ce7(:) * fc( 1: 7)) + u(n-4) = sum(ce7(:) * fc(n-7:n-1)) + u(n-3) = sum(ce7(:) * fc(n-6:n )) + u(n-2) = sum(ce5(:) * fc(n-4:n )) + u(n-1) = sum(ce3(:) * fc(n-2:n )) + u(n ) = fc(n ) + +! apply the monotonicity preserving limiting +! + call mp_limiting(fc(:), u(:)) + +! copy the interpolation to the respective vector +! + fl(1:n) = u(1:n) + +!! === right-side interpolation === +!! +! invert the cell-centered value vector +! + fi(1:n) = fc(n:1:-1) + +! reconstruct the interface state using the 9th order interpolation +! + do i = 5, n - 5 + u(i) = sum(ce9(:) * fi(i-4:i+5)) + end do + +! interpolate the interface state of the ghost zones using the interpolations +! of lower orders +! + u( 1) = sum(ce2(:) * fi( 1: 2)) + u( 2) = sum(ce3(:) * fi( 1: 3)) + u( 3) = sum(ce5(:) * fi( 1: 5)) + u( 4) = sum(ce7(:) * fi( 1: 7)) + u(n-4) = sum(ce7(:) * fi(n-7:n-1)) + u(n-3) = sum(ce7(:) * fi(n-6:n )) + u(n-2) = sum(ce5(:) * fi(n-4:n )) + u(n-1) = sum(ce3(:) * fi(n-2:n )) + u(n ) = fi(n ) + +! apply the monotonicity preserving limiting +! + call mp_limiting(fi(:), u(:)) + +! copy the interpolation to the respective vector +! + fr(1:n-1) = u(n-1:1:-1) + +! update the interpolation of the first and last points +! + i = n - 1 + fl(1) = 0.5d+00 * (fc(1) + fc(2)) + fr(i) = 0.5d+00 * (fc(i) + fc(n)) + fl(n) = fc(n) + fr(n) = fc(n) + +!------------------------------------------------------------------------------- +! + end subroutine reconstruct_mp9ld +! +!=============================================================================== +! ! subroutine RECONSTRUCT_CRMP5: ! ---------------------------- !