Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Floating point exception catching fixes.
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Jan 6, 2023
1 parent 5edb71d commit f0dca45
Show file tree
Hide file tree
Showing 5 changed files with 28 additions and 25 deletions.
2 changes: 0 additions & 2 deletions src/collision/collision_check.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,7 +114,6 @@ module subroutine collision_check_plpl(self, nbody_system, param, t, dt, irec, l
lany_collision = any(lcollision(:))
lany_closest = (param%lenc_save_closest .and. any(self%lclosest(:)))


if (lany_collision .or. lany_closest) then
call pl%rh2rb(nbody_system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary
do k = 1, nenc
Expand Down Expand Up @@ -176,7 +175,6 @@ module subroutine collision_check_pltp(self, nbody_system, param, t, dt, irec, l
logical, dimension(:), allocatable :: lcollision, lmask
real(DP), dimension(NDIM) :: xr, vr
integer(I4B) :: i, j, k, nenc
real(DP) :: rlim
logical :: lany_closest
character(len=STRMAX) :: timestr, idstri, idstrj, message
class(collision_list_pltp), allocatable :: tmp
Expand Down
14 changes: 9 additions & 5 deletions src/encounter/encounter_check.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -706,12 +706,16 @@ elemental module subroutine encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc,
r2min = r2
else
v2 = vxr**2 + vyr**2 + vzr**2
tmin = -vdotr / v2

if (tmin < dt) then
r2min = r2 - vdotr**2 / v2
if (v2 <= VSMALL) then
r2min = r2
else
r2min = r2 + 2 * vdotr * dt + v2 * dt**2
tmin = -vdotr / v2

if (tmin < dt) then
r2min = r2 - vdotr**2 / v2
else
r2min = r2 + 2 * vdotr * dt + v2 * dt**2
end if
end if
end if
else
Expand Down
2 changes: 1 addition & 1 deletion src/globals/globals_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,6 +118,6 @@ module globals
!> Miscellaneous constants:
integer(I4B), parameter :: NDIM = 3 !! Number of dimensions in our reality
integer(I4B), parameter :: NDIM2 = 2 * NDIM !! 2x the number of dimensions
real(DP), parameter :: VSMALL = 2 * epsilon(1._DP) !! Very small number used to prevent floating underflow
real(DP), parameter :: VSMALL = sqrt(TINY(1._DP)) !! Very small number used to prevent floating underflow

end module globals
33 changes: 17 additions & 16 deletions src/swiftest/swiftest_orbel.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
!! If not, see: https://www.gnu.org/licenses.

submodule (swiftest) s_swiftest_orbel
real(DP), parameter :: TINYVALUE = 4.0e-15_DP !! Tiny value used to prevent floating point errors. Value set based on the Swifter TINY parameter.
contains

module subroutine swiftest_orbel_el2xv_vec(self, cb)
Expand Down Expand Up @@ -68,7 +69,7 @@ pure subroutine swiftest_orbel_el2xv(mu, a, ie, inc, capom, omega, capm, x, v)
else
e = ie
em1 = e - 1._DP
if (abs(em1) < VSMALL) then
if (abs(em1) < TINYVALUE) then
iorbit_type = PARABOLA
else if (e > 1.0_DP) then
iorbit_type = HYPERBOLA
Expand Down Expand Up @@ -258,17 +259,17 @@ real(DP) pure function swiftest_orbel_flon(e,icapn)
bigb = -(+0.5_DP * b + sq)**(1.0_DP / 3.0_DP)
x = biga + bigb
swiftest_orbel_flon = x
! If capn is VSMALL (or zero) no need to go further than cubic even for
! If capn is TINYVALUE (or zero) no need to go further than cubic even for
! e =1.
if( capn >= VSMALL) then
if( capn >= TINYVALUE) then
do i = 1,IMAX
x2 = x**2
f = a0 + x * (a1 + x2 * (a3 + x2 * (a5 + x2 * (a7 + x2 * (a9 + x2 * (a11 + x2))))))
fp = b1 + x2 * (b3 + x2 * (b5 + x2 * (b7 + x2 * (b9 + x2 * (b11 + 13 * x2)))))
dx = -f / fp
swiftest_orbel_flon = x + dx
! if we have converged here there's no point in going on
if(abs(dx) <= VSMALL) exit
if(abs(dx) <= TINYVALUE) exit
x = swiftest_orbel_flon
end do
end if
Expand Down Expand Up @@ -344,7 +345,7 @@ real(DP) pure function swiftest_orbel_fget(e,capn)
dx = -f / (fp + dx * fpp / 2._DP + dx**2 * fppp / 6._DP)
swiftest_orbel_fget = x + dx
! if we have converged here there's no point in going on
if(abs(dx) <= VSMALL) return
if(abs(dx) <= TINYVALUE) return
x = swiftest_orbel_fget
end do

Expand Down Expand Up @@ -720,13 +721,13 @@ pure elemental module subroutine swiftest_orbel_xv2aeq(mu, px, py, pz, vx, vy, v
h2 = dot_product(hvec(:), hvec(:))
if (h2 == 0.0_DP) return
energy = 0.5_DP * v2 - mu / r
if (abs(energy * r / mu) < sqrt(VSMALL)) then
if (abs(energy * r / mu) < sqrt(TINYVALUE)) then
iorbit_type = PARABOLA
else
a = -0.5_DP * mu / energy
if (a < 0.0_DP) then
fac = -h2 / (mu * a)
if (fac > VSMALL) then
if (fac > TINYVALUE) then
iorbit_type = HYPERBOLA
else
iorbit_type = PARABOLA
Expand All @@ -738,7 +739,7 @@ pure elemental module subroutine swiftest_orbel_xv2aeq(mu, px, py, pz, vx, vy, v
select case (iorbit_type)
case (ELLIPSE)
fac = 1.0_DP - h2 / (mu * a)
if (fac > VSMALL) e = sqrt(fac)
if (fac > TINYVALUE) e = sqrt(fac)
q = a * (1.0_DP - e)
case (PARABOLA)
a = 0.5_DP * h2 / mu
Expand Down Expand Up @@ -790,13 +791,13 @@ pure module subroutine swiftest_orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, c
if (h2 == 0.0_DP) return
rdotv = dot_product(x(:), v(:))
energy = 0.5_DP * v2 - mu / r
if (abs(energy * r / mu) < sqrt(VSMALL)) then
if (abs(energy * r / mu) < sqrt(TINYVALUE)) then
iorbit_type = PARABOLA
else
a = -0.5_DP * mu / energy
if (a < 0.0_DP) then
fac = -h2 / (mu * a)
if (fac > VSMALL) then
if (fac > TINYVALUE) then
iorbit_type = HYPERBOLA
else
iorbit_type = PARABOLA
Expand All @@ -808,7 +809,7 @@ pure module subroutine swiftest_orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, c
select case (iorbit_type)
case (ELLIPSE)
fac = 1.0_DP - h2 / (mu * a)
if (fac > VSMALL) then
if (fac > TINYVALUE) then
e = sqrt(fac)
cape = 0.0_DP
face = (a - r) / (a * e)
Expand Down Expand Up @@ -955,7 +956,7 @@ pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, in
inc = acos(fac)
end if
fac = sqrt(hvec(1)**2 + hvec(2)**2) / h
if (fac**2 < VSMALL) then
if (fac**2 < TINYVALUE) then
u = atan2(py, px)
if (hvec(3) < 0.0_DP) u = -u
else
Expand All @@ -964,13 +965,13 @@ pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, in
end if
if (capom < 0.0_DP) capom = capom + TWOPI
if (u < 0.0_DP) u = u + TWOPI
if (abs(energy * r / mu) < sqrt(VSMALL)) then
if (abs(energy * r / mu) < sqrt(TINYVALUE)) then
iorbit_type = parabola
else
a = -0.5_DP * mu / energy
if (a < 0.0_DP) then
fac = -h2 / (mu * a)
if (fac > VSMALL) then
if (fac > TINYVALUE) then
iorbit_type = HYPERBOLA
else
iorbit_type = PARABOLA
Expand All @@ -982,7 +983,7 @@ pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, in
select case (iorbit_type)
case (ELLIPSE)
fac = 1.0_DP - h2 / (mu * a)
if (fac > VSMALL) then
if (fac > TINYVALUE) then
e = sqrt(fac)
cape = 0.0_DP
face = (a - r) / (a * e)
Expand Down Expand Up @@ -1031,7 +1032,7 @@ pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, in
if (omega < 0.0_DP) omega = omega + TWOPI
varpi = mod(omega + capom, TWOPI)
lam = mod(capm + varpi, TWOPI)
if (e > VSMALL) then
if (e > TINYVALUE) then
cf = 1.0_DP / e * (a * (1.0_DP - e**2)/r - 1.0_DP)
rdot = sign(sqrt(v2 - (h / r)**2),rdotv)
sf = a * (1.0_DP - e**2) / (h * e) * rdot
Expand Down
2 changes: 1 addition & 1 deletion src/symba/symba_step.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -194,7 +194,7 @@ recursive module subroutine symba_step_recur_system(self, param, t, ireci)
nbody_system%irec = ireci
dtl = param%dt / (NTENC**ireci)
dth = 0.5_DP * dtl
IF (dtl / param%dt < VSMALL) THEN
IF (dtl / param%dt < epsilon(1.0_DP)) THEN
write(*, *) "SWIFTEST Warning:"
write(*, *) " In symba_step_recur_system, local time step is too small"
write(*, *) " Roundoff error will be important!"
Expand Down

0 comments on commit f0dca45

Please sign in to comment.