From b4dabe6674a88dc0bf0c9efcaaa2855372ca4691 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 21 Oct 2022 11:27:17 -0400 Subject: [PATCH 1/2] Fixed issue that occurs when collisions are head-on and there is no initial orbital angular momentum. This caused a divide-by-zero error when normalizing the momentum vector --- src/fraggle/fraggle_set.f90 | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/fraggle/fraggle_set.f90 b/src/fraggle/fraggle_set.f90 index 6902bfb48..baf29f8b0 100644 --- a/src/fraggle/fraggle_set.f90 +++ b/src/fraggle/fraggle_set.f90 @@ -163,7 +163,7 @@ module subroutine fraggle_set_coordinate_system(self, colliders) ! Internals integer(I4B) :: i real(DP), dimension(NDIM) :: x_cross_v, delta_r, delta_v, Ltot - real(DP) :: r_col_norm, v_col_norm + real(DP) :: r_col_norm, v_col_norm, L_mag real(DP), dimension(NDIM, self%nbody) :: L_sigma associate(frag => self, nfrag => self%nbody) @@ -176,7 +176,12 @@ module subroutine fraggle_set_coordinate_system(self, colliders) ! and the y-axis aligned with the pre-impact distance vector. Ltot = colliders%L_orbit(:,1) + colliders%L_orbit(:,2) + colliders%L_spin(:,1) + colliders%L_spin(:,2) frag%y_coll_unit(:) = delta_r(:) / r_col_norm - frag%z_coll_unit(:) = Ltot(:) / (.mag. Ltot(:)) + L_mag = .mag.Ltot(:) + if (L_mag > tiny(L_mag)) then + frag%z_coll_unit(:) = Ltot(:) / L_mag + else + frag%z_coll_unit(:) = 0.0_DP + end if ! The cross product of the y- by z-axis will give us the x-axis frag%x_coll_unit(:) = frag%y_coll_unit(:) .cross. frag%z_coll_unit(:) @@ -228,6 +233,7 @@ module subroutine fraggle_set_natural_scale_factors(self, colliders) colliders%mass(:) = colliders%mass(:) / frag%mscale colliders%radius(:) = colliders%radius(:) / frag%dscale colliders%L_spin(:,:) = colliders%L_spin(:,:) / frag%Lscale + colliders%L_orbit(:,:) = colliders%L_orbit(:,:) / frag%Lscale do i = 1, 2 colliders%rot(:,i) = colliders%L_spin(:,i) / (colliders%mass(i) * colliders%radius(i)**2 * colliders%Ip(3, i)) From cd31bf3329182add4d8dbfaf782441936142bb33 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 21 Oct 2022 11:32:31 -0400 Subject: [PATCH 2/2] Pulled the final momentum calc for the colliders out of the consolidation loop where it didn't belong --- src/symba/symba_collision.f90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index edee2455e..dcb2b674b 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -549,17 +549,17 @@ function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, collid density(j) = colliders%mass(j) / volume(j) colliders%radius(j) = (3 * volume(j) / (4 * PI))**(1.0_DP / 3.0_DP) if (param%lrotation) colliders%Ip(:, j) = colliders%Ip(:, j) / colliders%mass(j) - - xcom(:) = (colliders%mass(1) * colliders%xb(:, 1) + colliders%mass(2) * colliders%xb(:, 2)) / sum(colliders%mass(:)) - vcom(:) = (colliders%mass(1) * colliders%vb(:, 1) + colliders%mass(2) * colliders%vb(:, 2)) / sum(colliders%mass(:)) - mxc(:, 1) = colliders%mass(1) * (colliders%xb(:, 1) - xcom(:)) - mxc(:, 2) = colliders%mass(2) * (colliders%xb(:, 2) - xcom(:)) - vcc(:, 1) = colliders%vb(:, 1) - vcom(:) - vcc(:, 2) = colliders%vb(:, 2) - vcom(:) - colliders%L_orbit(:,:) = mxc(:,:) .cross. vcc(:,:) end do lflag = .true. + xcom(:) = (colliders%mass(1) * colliders%xb(:, 1) + colliders%mass(2) * colliders%xb(:, 2)) / sum(colliders%mass(:)) + vcom(:) = (colliders%mass(1) * colliders%vb(:, 1) + colliders%mass(2) * colliders%vb(:, 2)) / sum(colliders%mass(:)) + mxc(:, 1) = colliders%mass(1) * (colliders%xb(:, 1) - xcom(:)) + mxc(:, 2) = colliders%mass(2) * (colliders%xb(:, 2) - xcom(:)) + vcc(:, 1) = colliders%vb(:, 1) - vcom(:) + vcc(:, 2) = colliders%vb(:, 2) - vcom(:) + colliders%L_orbit(:,:) = mxc(:,:) .cross. vcc(:,:) + ! Destroy the kinship relationships for all members of this colliders%idx call pl%reset_kinship(colliders%idx(:))