From 8c0c89d0af7a8a50d71bbe1e1de5b0a921d6bda9 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 7 Aug 2018 14:59:42 -0300 Subject: [PATCH] INTERPOLATIONS: Avoid infinite loop in interfaces_tvd() and interfaces_mgp(). In case an interpolated cell, which should be positive, is zero or negative, and the sum of derivatives is zero as well, the correction of derivatives can enter into an infinite loop. It shouldn't normally happen, since such a situation is unphysical. Avoid the infinite loop and print a warning about encountered problems. Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index 617d486..40af9e2 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -695,6 +695,7 @@ module interpolations use coordinates , only : im , jm , km use coordinates , only : ib , jb , kb , ie , je , ke use coordinates , only : ibl, jbl, kbl, ieu, jeu, keu + use error , only : print_warning ! local variables are not implicit by default ! @@ -768,9 +769,15 @@ module interpolations ! variables ! if (positive) then - do while (q(i,j,k) <= sum(abs(dq(1:NDIMS)))) - dq(:) = 0.5d+00 * dq(:) - end do + if (q(i,j,k) > 0.0d+00) then + do while (q(i,j,k) <= sum(abs(dq(1:NDIMS)))) + dq(:) = 0.5d+00 * dq(:) + end do + else + call print_warning("interpolations::interfaces_tvd" & + , "Positive variable is not positive!") + dq(:) = 0.0d+00 + end if end if ! interpolate states @@ -933,6 +940,7 @@ module interpolations use coordinates , only : im , jm , km use coordinates , only : ib , jb , kb , ie , je , ke use coordinates , only : ibl, jbl, kbl, ieu, jeu, keu + use error , only : print_warning ! local variables are not implicit by default ! @@ -1051,9 +1059,15 @@ module interpolations ! variables ! if (positive) then - do while (q(i,j,k) <= sum(abs(dq(1:NDIMS)))) - dq(:) = 0.5d+00 * dq(:) - end do + if (q(i,j,k) > 0.0d+00) then + do while (q(i,j,k) <= sum(abs(dq(1:NDIMS)))) + dq(:) = 0.5d+00 * dq(:) + end do + else + call print_warning("interpolations::interfaces_mgp" & + , "Positive variable is not positive!") + dq(:) = 0.0d+00 + end if end if ! interpolate states