diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index c3dda504c..7e72c8f90 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -100,7 +100,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur real(DP), intent(in) :: t !! Time of collision logical, optional, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? ! Internals - integer(I4B) :: try, subtry + integer(I4B) :: try real(DP) :: dEtot, dLmag integer(I4B), parameter :: MAXTRY = 100 logical :: lk_plpl, lfailure_local @@ -140,55 +140,49 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur lfailure_local = .false. try = 1 - do while (try < MAXTRY) - write(message,*) try - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle try " // trim(adjustl(message))) - if (lfailure_local) then - call fragments%reset() - try = try + 1 - end if - lfailure_local = .false. - - ! Use the disruption collision model to generate initial conditions - ! Compute the "before" energy/momentum and the budgets - call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) - call self%set_budgets() - do subtry = 1, MAXTRY - call fraggle_generate_pos_vec(self) - call fraggle_generate_rot_vec(self) - call fraggle_generate_vel_vec(self,lfailure_local) - if (.not.lfailure_local) exit - end do + write(message,*) try + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle try " // trim(adjustl(message))) + if (lfailure_local) then + call fragments%reset() + try = try + 1 + end if + lfailure_local = .false. - call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) - dEtot = self%Etot(2) - self%Etot(1) - dLmag = .mag. (self%Ltot(:,2) - self%Ltot(:,1)) - - lfailure_local = (dEtot > -impactors%Qloss) - if (lfailure_local) then - write(message, *) "dEtot: ",dEtot, "dEtot/Qloss", dEtot / impactors%Qloss - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to energy gain: " // & - trim(adjustl(message))) - cycle - lfailure_local = .false. - end if + ! Use the disruption collision model to generate initial conditions + ! Compute the "before" energy/momentum and the budgets + call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) + call self%set_budgets() + do try = 1, MAXTRY + call fraggle_generate_pos_vec(self) + call fraggle_generate_rot_vec(self) + call fraggle_generate_vel_vec(self,lfailure_local) + if (.not.lfailure_local) exit + end do - lfailure_local = ((abs(dLmag) / (.mag.self%Ltot(:,1))) > FRAGGLE_LTOL) - if (lfailure_local) then - write(message,*) dLmag / (.mag.self%Ltot(:,1)) - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to high angular momentum error: " // & - trim(adjustl(message))) - cycle - end if + call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) + dEtot = self%Etot(2) - self%Etot(1) + dLmag = .mag. (self%Ltot(:,2) - self%Ltot(:,1)) - ! Check if any of the usual floating point exceptions happened, and fail the try if so - call ieee_get_flag(ieee_usual, fpe_flag) - lfailure_local = any(fpe_flag) - if (.not.lfailure_local) exit - write(message,*) "Fraggle failed due to a floating point exception: ", fpe_flag - call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + lfailure_local = (dEtot > -impactors%Qloss) + if (lfailure_local) then + write(message, *) "dEtot: ",dEtot, "dEtot/Qloss", dEtot / impactors%Qloss + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to energy gain: " // & + trim(adjustl(message))) + end if + + lfailure_local = ((abs(dLmag) / (.mag.self%Ltot(:,1))) > FRAGGLE_LTOL) + if (lfailure_local) then + write(message,*) dLmag / (.mag.self%Ltot(:,1)) + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to high angular momentum error: " // & + trim(adjustl(message))) + end if + + ! Check if any of the usual floating point exceptions happened, and fail the try if so + call ieee_get_flag(ieee_usual, fpe_flag) + lfailure_local = any(fpe_flag) + write(message,*) "Fraggle failed due to a floating point exception: ", fpe_flag + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) - end do write(message,*) try if (lfailure_local) then call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle fragment generation failed after " // & @@ -530,7 +524,7 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure) fragments%vc(:,i) = fragments%vmag(i) * .unit.fragments%vc(:,i) end do do concurrent(i = istart:nfrag, fragments%ke_spin_frag(i) > ke_per_dof) - rotmag = fragments%rotmag(i)**2 - 2*ke_per_dof/(fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) + rotmag = fragments%rotmag(i)**2 - ke_per_dof/(fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) rotmag = max(rotmag, 0.0_DP) fragments%rotmag(i) = sqrt(rotmag) fragments%rot(:,i) = fragments%rotmag(i) * .unit.fragments%rot(:,i)