diff --git a/src/interpolations.F90 b/src/interpolations.F90 index 8d5eba4..f75fd03 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -184,6 +184,7 @@ module interpolations if (verbose .and. ng < 2) & call print_warning("interpolations:initialize_interpolation" & , "Increase the number of ghost cells (at least 2).") + eps = max(1.0d-12, eps) case ("weno5z", "weno5-z", "WENO5Z", "WENO5-Z") name_rec = "5th order WENO-Z (Borges et al. 2008)" reconstruct_states => reconstruct_weno5z @@ -484,12 +485,12 @@ module interpolations ! ! calculate the left- and right-side interface interpolations ! - do i = 1, n + do i = 2, n - 1 ! calculate left and right indices ! - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 ! calculate left and right side derivatives ! @@ -505,11 +506,14 @@ module interpolations fl(i ) = f(i) + df fr(im1) = f(i) - df - end do ! i = 1, n + end do ! i = 2, n - 1 ! update the interpolation of the first and last points ! - fl(1) = f(1) + i = n - 1 + fl(1) = 0.5d+00 * (f(1) + f(2)) + fr(i) = 0.5d+00 * (f(i) + f(n)) + fl(n) = f(n) fr(n) = f(n) !------------------------------------------------------------------------------- @@ -563,12 +567,12 @@ module interpolations ! iterate along the vector ! - do i = 1, n + do i = 2, n - 1 ! prepare neighbour indices ! - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 ! calculate the left and right derivatives ! @@ -625,12 +629,15 @@ module interpolations ! fr(im1) = (wp * fp + wm * fm) / ww - end do ! i = 1, n + end do ! i = 2, n - 1 ! update the interpolation of the first and last points ! - fl(1) = f (1) - fr(n) = fl(n) + i = n - 1 + fl(1) = 0.5d+00 * (f(1) + f(2)) + fr(i) = 0.5d+00 * (f(i) + f(n)) + fl(n) = f(n) + fr(n) = f(n) !------------------------------------------------------------------------------- ! @@ -686,12 +693,12 @@ module interpolations ! iterate over positions and interpolate states ! - do i = 1, n + do i = 2, n - 1 ! prepare neighbour indices ! - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 ! prepare left and right differences ! @@ -709,11 +716,7 @@ module interpolations ! calculate values at i + ½ ! - if (dfr == 0.0d+00) then - - fl(i) = f(i) - - else + if (abs(dfr) > eps) then ! calculate the slope ratio (eq. 2.8 in [1]) ! @@ -735,15 +738,15 @@ module interpolations ! fl(i) = f(i) + dfr * (xl * f1 + xi * f2) + else + + fl(i) = f(i) + end if ! calculate values at i - ½ ! - if (dfl == 0.0d+00) then - - fr(im1) = f(i) - - else + if (abs(dfl) > eps) then ! calculate the slope ratio (eq. 2.8 in [1]) ! @@ -765,14 +768,21 @@ module interpolations ! fr(im1) = f(i) - dfl * (xl * f1 + xi * f2) + else + + fr(im1) = f(i) + end if - end do ! i = 1, n + end do ! i = 2, n - 1 ! update the interpolation of the first and last points ! - fl(1) = f (1) - fr(n) = fl(n) + i = n - 1 + fl(1) = 0.5d+00 * (f(1) + f(2)) + fr(i) = 0.5d+00 * (f(i) + f(n)) + fl(n) = f(n) + fr(n) = f(n) !------------------------------------------------------------------------------- ! @@ -863,13 +873,13 @@ module interpolations ! iterate along the vector ! - do i = 1, n + do i = 2, n - 1 ! prepare neighbour indices ! - im1 = max(1, i - 1) im2 = max(1, i - 2) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 ip2 = min(n, i + 2) ! calculate βₖ (eqs. 9-11 in [1]) @@ -928,12 +938,15 @@ module interpolations ! fr(im1) = (wl * ql + wr * qr) + wc * qc - end do ! i = 1, n + end do ! i = 2, n - 1 ! update the interpolation of the first and last points ! - fl(1) = fr(1) - fr(n) = fl(n) + i = n - 1 + fl(1) = 0.5d+00 * (f(1) + f(2)) + fr(i) = 0.5d+00 * (f(i) + f(n)) + fl(n) = f(n) + fr(n) = f(n) !------------------------------------------------------------------------------- ! @@ -2489,11 +2502,11 @@ module interpolations ! obtain the face values using high order interpolation ! - do i = 1, n + do i = 2, n - 1 im2 = max(1, i - 2) - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 ip2 = min(n, i + 2) fl(i) = (4.7d+01 * f(i ) + (2.7d+01 * f(ip1) - 1.3d+01 * f(im1)) & @@ -2503,14 +2516,14 @@ module interpolations - (3.0d+00 * f(im2) - 2.0d+00 * f(ip2))) & / 6.0d+01 - end do ! i = 1, n + end do ! i = 2, n - 1 ! apply monotonicity preserving limiting ! - do i = 1, n + do i = 2, n - 1 - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 if (dfm(i) * dfp(i) >= 0.0d+00) then sigma = kappa @@ -2580,12 +2593,15 @@ module interpolations ! fr(im1) = fr(i) - end do + end do ! n = 2, n - 1 ! update the interpolation of the first and last points ! - fl(1) = fr(1) - fr(n) = fl(n) + i = n - 1 + fl(1) = 0.5d+00 * (f(1) + f(2)) + fr(i) = 0.5d+00 * (f(i) + f(n)) + fl(n) = f(n) + fr(n) = f(n) !------------------------------------------------------------------------------- ! @@ -2662,10 +2678,10 @@ module interpolations ! prepare the tridiagonal system coefficients for the interior ! - do i = 1, n + do i = ng, n - ng + 1 - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 a(i,1) = 3.0d-01 b(i,1) = 6.0d-01 @@ -2678,15 +2694,36 @@ module interpolations r(i,1) = (f(im1) + 1.9d+01 * f(i ) + 1.0d+01 * f(ip1)) / 3.0d+01 r(i,2) = (f(ip1) + 1.9d+01 * f(i ) + 1.0d+01 * f(im1)) / 3.0d+01 - end do ! i = 1, n + end do ! i = ng, n - ng + 1 ! interpolate ghost zones using explicit method (left-side reconstruction) ! - do i = 1, ng + do i = 2, ng im2 = max(1, i - 2) - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 + ip2 = i + 2 + + a(i,1) = 0.0d+00 + b(i,1) = 1.0d+00 + c(i,1) = 0.0d+00 + + r(i,1) = (4.7d+01 * f(i ) + (2.7d+01 * f(ip1) - 1.3d+01 * f(im1)) & + - (3.0d+00 * f(ip2) - 2.0d+00 * f(im2))) & + / 6.0d+01 + + end do ! i = 2, ng + a(1,1) = 0.0d+00 + b(1,1) = 1.0d+00 + c(1,1) = 0.0d+00 + r(1,1) = 0.5d+00 * (f(1) + f(2)) + + do i = n - ng, n - 1 + + im2 = i - 2 + im1 = i - 1 + ip1 = i + 1 ip2 = min(n, i + 2) a(i,1) = 0.0d+00 @@ -2697,32 +2734,40 @@ module interpolations - (3.0d+00 * f(ip2) - 2.0d+00 * f(im2))) & / 6.0d+01 - end do ! i = 1, ng - - do i = n - ng, n - - im2 = max(1, i - 2) - im1 = max(1, i - 1) - ip1 = min(n, i + 1) - ip2 = min(n, i + 2) - - a(i,1) = 0.0d+00 - b(i,1) = 1.0d+00 - c(i,1) = 0.0d+00 - - r(i,1) = (4.7d+01 * f(i ) + (2.7d+01 * f(ip1) - 1.3d+01 * f(im1)) & - - (3.0d+00 * f(ip2) - 2.0d+00 * f(im2))) & - / 6.0d+01 - - end do ! i = n - ng, n + end do ! i = n - ng, n - 1 + a(n,1) = 0.0d+00 + b(n,1) = 1.0d+00 + c(n,1) = 0.0d+00 + r(n,1) = f(n) ! interpolate ghost zones using explicit method (right-side reconstruction) ! - do i = 1, ng + 1 + do i = 2, ng + 1 im2 = max(1, i - 2) - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 + ip2 = i + 2 + + a(i,2) = 0.0d+00 + b(i,2) = 1.0d+00 + c(i,2) = 0.0d+00 + + r(i,2) = (4.7d+01 * f(i ) + (2.7d+01 * f(im1) - 1.3d+01 * f(ip1)) & + - (3.0d+00 * f(im2) - 2.0d+00 * f(ip2))) & + / 6.0d+01 + + end do ! i = 2, ng + 1 + a(1,2) = 0.0d+00 + b(1,2) = 1.0d+00 + c(1,2) = 0.0d+00 + r(1,2) = f(1) + + do i = n - ng + 1, n - 1 + + im2 = i - 2 + im1 = i - 1 + ip1 = i + 1 ip2 = min(n, i + 2) a(i,2) = 0.0d+00 @@ -2733,24 +2778,11 @@ module interpolations - (3.0d+00 * f(im2) - 2.0d+00 * f(ip2))) & / 6.0d+01 - end do ! i = 1, ng + 1 - - do i = n - ng + 1, n - - im2 = max(1, i - 2) - im1 = max(1, i - 1) - ip1 = min(n, i + 1) - ip2 = min(n, i + 2) - - a(i,2) = 0.0d+00 - b(i,2) = 1.0d+00 - c(i,2) = 0.0d+00 - - r(i,2) = (4.7d+01 * f(i ) + (2.7d+01 * f(im1) - 1.3d+01 * f(ip1)) & - - (3.0d+00 * f(im2) - 2.0d+00 * f(ip2))) & - / 6.0d+01 - - end do ! i = n - ng + 1, n + end do ! i = n - ng + 1, n - 1 + a(n,2) = 0.0d+00 + b(n,2) = 1.0d+00 + c(n,2) = 0.0d+00 + r(n,2) = 0.5d+00 * (f(n-1) + f(n)) ! solve the tridiagonal system of equations for the left-side interpolation ! @@ -2758,10 +2790,10 @@ module interpolations ! apply the monotonicity preserving limiting ! - do i = 1, n + do i = 2, n - 1 - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 if (dfm(i) * dfp(i) >= 0.0d+00) then sigma = kappa @@ -2798,7 +2830,7 @@ module interpolations end if - end do + end do ! i = 2, n - 1 ! solve the tridiagonal system of equations for the right-side interpolation ! @@ -2806,10 +2838,10 @@ module interpolations ! apply the monotonicity preserving limiting ! - do i = 1, n + do i = 2, n - 1 - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 if (dfm(i) * dfp(i) >= 0.0d+00) then sigma = kappa @@ -2851,12 +2883,15 @@ module interpolations ! fr(im1) = fr(i) - end do + end do ! i = 2, n - 1 ! update the interpolation of the first and last points ! - fl(1) = fr(1) - fr(n) = fl(n) + i = n - 1 + fl(1) = 0.5d+00 * (f(1) + f(2)) + fr(i) = 0.5d+00 * (f(i) + f(n)) + fl(n) = f(n) + fr(n) = f(n) !------------------------------------------------------------------------------- !