From 40232d24554e2782f64e398747d7f67c477bc9b7 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 21 Oct 2022 09:45:40 -0400 Subject: [PATCH] Fixed NaN issue in potential energy computation. Inconsistent array bounds on mask (this may be the core reason why fraggle fails sometimes) --- src/util/util_get_energy_momentum.f90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index ae59c158a..65e9b90d0 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -196,10 +196,12 @@ subroutine util_get_energy_potential_triangular(npl, lmask, GMcb, Gmass, mass, x !$omp firstprivate(npl) & !$omp reduction(+:pe) do i = 1, npl - do concurrent(j = i+1:npl, lmask(i) .and. lmask(j)) - pepl(j) = - (Gmass(i) * mass(j)) / norm2(xb(:, i) - xb(:, j)) - end do - pe = pe + sum(pepl(i+1:npl), lmask(i)) + if (lmask(i)) then + do concurrent(j = i+1:npl, lmask(i) .and. lmask(j)) + pepl(j) = - (Gmass(i) * mass(j)) / norm2(xb(:, i) - xb(:, j)) + end do + pe = pe + sum(pepl(i+1:npl), lmask(i+1:npl)) + end if end do !$omp end parallel do pe = pe + sum(pecb(1:npl), lmask(1:npl))