From 6b138aaa4ea146f2a823dff7da32eaeaf9a5a30d Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 17 Feb 2023 15:04:41 -0500 Subject: [PATCH] Added a check to see if bodies are overlapping before shifting them in the collider consolidation routine. Otherwise, a floating point exception can be raised --- src/collision/collision_resolve.f90 | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index 3b5402162..b9cf5663a 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -143,15 +143,17 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa ! Shift the impactors so that they are not overlapping rlim = sum(impactors%radius(1:2)) - vrel = impactors%vb(:,2) - impactors%vb(:,1) rrel = impactors%rb(:,2) - impactors%rb(:,1) - vrel_mag = .mag. vrel rrel_mag = .mag. rrel - vrel_unit = .unit. vrel - r_dot_vunit = dot_product(rrel,vrel_unit) - b = sqrt(rrel_mag**2 - r_dot_vunit**2) - dt = (sqrt(rlim**2 - b**2) + r_dot_vunit)/vrel_mag - impactors%rb(:,1:2) = impactors%rb(:,1:2) - dt * impactors%vb(:,1:2) + if (rrel_mag < rlim) then + vrel = impactors%vb(:,2) - impactors%vb(:,1) + vrel_mag = .mag. vrel + vrel_unit = .unit. vrel + r_dot_vunit = dot_product(rrel,vrel_unit) + b = sqrt(rrel_mag**2 - r_dot_vunit**2) + dt = (sqrt(rlim**2 - b**2) + r_dot_vunit)/vrel_mag + impactors%rb(:,1:2) = impactors%rb(:,1:2) - dt * impactors%vb(:,1:2) + end if xcom(:) = (impactors%mass(1) * impactors%rb(:, 1) + impactors%mass(2) * impactors%rb(:, 2)) / sum(impactors%mass(:)) vcom(:) = (impactors%mass(1) * impactors%vb(:, 1) + impactors%mass(2) * impactors%vb(:, 2)) / sum(impactors%mass(:))