From b41b5a000b9833d389cccb57e11c490a456c1df2 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 5 Apr 2017 10:23:15 -0300 Subject: [PATCH] INTERPOLATIONS: Fix MLP limiting. We should take half of the TVD limited derivatives in order to compare properly with the high order interpolated derivative. Fix it. Also, TVD limit both states if the high order interpolation of any of them overshoot. Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 47 ++++++++++++++++++++++-------------------- 1 file changed, 25 insertions(+), 22 deletions(-) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index 02adebb..80fa8c1 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -1136,16 +1136,16 @@ module interpolations ! dql(1) = q(i ,j,k) - q(im1,j,k) dqr(1) = q(ip1,j,k) - q(i ,j,k) - dq (1) = minmod(dql(1), dqr(1)) + dq (1) = limiter_minmod(0.5d+00, dql(1), dqr(1)) dql(2) = q(i,j ,k) - q(i,jm1,k) dqr(2) = q(i,jp1,k) - q(i,j ,k) - dq (2) = minmod(dql(2), dqr(2)) + dq (2) = limiter_minmod(0.5d+00, dql(2), dqr(2)) #if NDIMS == 3 dql(3) = q(i,j,k ) - q(i,j,km1) dqr(3) = q(i,j,kp1) - q(i,j,k ) - dq (3) = minmod(dql(3), dqr(3)) + dq (3) = limiter_minmod(0.5d+00, dql(3), dqr(3)) #endif /* NDIMS == 3 */ ! calculate dqc @@ -1197,33 +1197,36 @@ module interpolations #endif /* NDIMS == 3 */ do m = 1, NDIMS - dq(m) = sign(ap(m), dq(m)) + dq(m) = sign(ap(m), dq(m)) end do - end if - -! calculate the limited variable increments -! - dql(1) = minmod(dq(1), qi(i ,j,k,1,1) - q(i,j,k)) - dqr(1) = minmod(dq(1), - qi(im1,j,k,2,1) + q(i,j,k)) - dql(2) = minmod(dq(2), qi(i,j ,k,1,2) - q(i,j,k)) - dqr(2) = minmod(dq(2), - qi(i,jm1,k,2,2) + q(i,j,k)) -#if NDIMS == 3 - dql(3) = minmod(dq(3), qi(i,j,k ,1,3) - q(i,j,k)) - dqr(3) = minmod(dq(3), - qi(i,j,km1,2,3) + q(i,j,k)) -#endif /* NDIMS == 3 */ ! update the interpolated states ! - qi(i ,j,k,1,1) = q(i,j,k) + dql(1) - qi(im1,j,k,2,1) = q(i,j,k) - dqr(1) + dql(1) = qi(i ,j,k,1,1) - q(i,j,k) + dqr(1) = qi(im1,j,k,2,1) - q(i,j,k) + if (max(abs(dql(1)), abs(dqr(1))) > abs(dq(1))) then + qi(i ,j,k,1,1) = q(i,j,k) + dq(1) + qi(im1,j,k,2,1) = q(i,j,k) - dq(1) + end if + + dql(2) = qi(i,j ,k,1,2) - q(i,j,k) + dqr(2) = qi(i,jm1,k,2,2) - q(i,j,k) + if (max(abs(dql(2)), abs(dqr(2))) > abs(dq(2))) then + qi(i,j ,k,1,2) = q(i,j,k) + dq(2) + qi(i,jm1,k,2,2) = q(i,j,k) - dq(2) + end if - qi(i,j ,k,1,2) = q(i,j,k) + dql(2) - qi(i,jm1,k,2,2) = q(i,j,k) - dqr(2) #if NDIMS == 3 - qi(i,j,k ,1,3) = q(i,j,k) + dql(3) - qi(i,j,km1,2,3) = q(i,j,k) - dqr(3) + dql(3) = qi(i,j,k ,1,3) - q(i,j,k)) + dqr(3) = qi(i,j,km1,2,3) - q(i,j,k)) + if (max(abs(dql(3)), abs(dqr(3))) > abs(dq(3))) then + qi(i,j,k ,1,3) = q(i,j,k) + dq(3) + qi(i,j,km1,2,3) = q(i,j,k) - dq(3) + end if #endif /* NDIMS == 3 */ + end if + end do ! i = ibl, ieu end do ! j = jbl, jeu end do ! k = kbl, keu