From bf4fc813a3cd6192c56e74ec369c3312e0a693c4 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 13 Jan 2023 09:07:47 -0500 Subject: [PATCH 1/6] Added new runtime checks to the debug compile mode to catch calls to unallocated variables --- cmake/Modules/SetFortranFlags.cmake | 6 +++--- src/fraggle/fraggle_generate.f90 | 32 ++++++++++++++++++++++------- src/swiftest/swiftest_util.f90 | 2 ++ 3 files changed, 30 insertions(+), 10 deletions(-) diff --git a/cmake/Modules/SetFortranFlags.cmake b/cmake/Modules/SetFortranFlags.cmake index 103ff7e45..5ce3c4925 100644 --- a/cmake/Modules/SetFortranFlags.cmake +++ b/cmake/Modules/SetFortranFlags.cmake @@ -110,10 +110,10 @@ SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG}" "-ftrace=full" # GNU (g95) ) -# Check array bounds +# Check everything SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG}" - Fortran "-check bounds" # Intel - "/check:bounds" # Intel Windows + Fortran "-check" # Intel + "/check" # Intel Windows "-fcheck=bounds" # GNU (New style) "-fbounds-check" # GNU (Old style) "-Mbounds" # Portland Group diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 51c4817e9..47e55e51a 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -447,7 +447,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters logical, intent(out) :: lfailure !! Did the velocity computation fail? ! Internals - integer(I4B) :: i, j, loop, try, istart, nfrag + integer(I4B) :: i, j, loop, try, istart, nfrag, nlast logical :: lhitandrun, lsupercat real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual real(DP) :: vmag, vesc, dE, E_residual, ke_avail, ke_remove, dE_best, E_residual_best, fscale, f_spin, f_orbit, dE_metric @@ -456,7 +456,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu ! For the initial "guess" of fragment velocities, this is the minimum and maximum velocity relative to escape velocity that the fragments will have real(DP) :: vmin_guess = 1.5_DP real(DP) :: vmax_guess = 10.0_DP - real(DP) :: delta_v + real(DP) :: delta_v, volume integer(I4B), parameter :: MAXLOOP = 20 integer(I4B), parameter :: MAXTRY = 100 real(DP), parameter :: SUCCESS_METRIC = 1.0e-3_DP @@ -497,7 +497,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu lfailure = .false. dE_metric = huge(1.0_DP) - outer: do try = 1, MAXTRY + outer: do try = 1, nfrag - 1 ! Scale the magnitude of the velocity by the distance from the impact point ! This will reduce the chances of fragments colliding with each other immediately, and is more physically correct do concurrent(i = 1:nfrag) @@ -551,7 +551,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu E_residual_best = E_residual dE_best = dE - do concurrent(i = 1:nfrag) + do concurrent(i = 1:fragments%nbody) fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:) end do @@ -576,7 +576,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu ke_remove = min(f_spin * E_residual, fragments%ke_spin_tot) ke_rot_remove(:) = ke_remove * (fragments%ke_spin(:) / fragments%ke_spin_tot) where(ke_rot_remove(:) > fragments%ke_spin(:)) ke_rot_remove(:) = fragments%ke_spin(:) - do concurrent(i = 1:nfrag, fragments%ke_spin(i) > 10*sqrt(tiny(1.0_DP))) + do concurrent(i = 1:fragments%nbody, fragments%ke_spin(i) > 10*sqrt(tiny(1.0_DP))) fscale = sqrt((fragments%ke_spin(i) - ke_rot_remove(i))/fragments%ke_spin(i)) fragments%rot(:,i) = fscale * fragments%rot(:,i) end do @@ -586,16 +586,34 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu call fragments%set_coordinate_system() end do - ! We didn't converge. Try another configuration and see if we get a better result + ! We didn't converge. Reset the fragment positions and velocities and try a new configuration with some slightly different parameters + ! Reduce the number of fragments by one + nlast = fragments%nbody + fragments%Ip(:,1) = fragments%mass(1) * impactors%Ip(:,1) + fragments%mass(nlast) * impactors%Ip(:,nlast) + fragments%mass(1) = fragments%mass(1) + fragments%mass(nlast) + fragments%Ip(:,1) = fragments%Ip(:,1) / fragments%mass(1) + fragments%Gmass(1) = fragments%Gmass(1) + fragments%mass(nlast) + volume = 4.0_DP / 3.0_DP * PI * (fragments%radius(1) + fragments%radius(nlast)) + fragments%density(1) = fragments%mass(1) / volume + fragments%radius(1) = (3._DP * volume / (4._DP * PI))**(THIRD) + fragments%Ip(:,nlast) = 0.0_DP + fragments%mass(nlast) = 0.0_DP + fragments%Gmass(nlast) = 0.0_DP + fragments%radius(nlast) = 0.0_DP + fragments%nbody = nlast - 1 + call fragments%reset() call fraggle_generate_pos_vec(collider_local) call fraggle_generate_rot_vec(collider_local, nbody_system, param) + + ! Increase the spatial size factor to get a less dense cloud collider_local%fail_scale = collider_local%fail_scale*1.01_DP - ! Bring the minimum and maximum velocities closer together for the next round + ! Bring the minimum and maximum velocities closer together delta_v = 0.125_DP * (vmax_guess - vmin_guess) vmin_guess = vmin_guess + delta_v vmax_guess = vmax_guess - delta_v + end do outer lfailure = dE_best > 0.0_DP diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 48c325b46..d79a9b042 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -1770,6 +1770,8 @@ module subroutine swiftest_util_peri_body(self, nbody_system, param) ! Internals integer(I4B) :: i + if (self%nbody == 0) return + select type(self) class is (swiftest_pl) if (self%lfirst) self%isperi(:) = 0 From 8b8be4faedf3ee4f208aacc1c0930798cb50ca5f Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 13 Jan 2023 09:10:32 -0500 Subject: [PATCH 2/6] More fixes to calls to unallocateds. Also work on a scheme to reduce the number of fragments when the Fraggle velocity function fails to converge. --- src/fraggle/fraggle_generate.f90 | 2 +- src/swiftest/swiftest_discard.f90 | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 47e55e51a..ddb14c261 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -589,7 +589,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu ! We didn't converge. Reset the fragment positions and velocities and try a new configuration with some slightly different parameters ! Reduce the number of fragments by one nlast = fragments%nbody - fragments%Ip(:,1) = fragments%mass(1) * impactors%Ip(:,1) + fragments%mass(nlast) * impactors%Ip(:,nlast) + fragments%Ip(:,1) = fragments%mass(1) * fragments%Ip(:,1) + fragments%mass(nlast) * fragments%Ip(:,nlast) fragments%mass(1) = fragments%mass(1) + fragments%mass(nlast) fragments%Ip(:,1) = fragments%Ip(:,1) / fragments%mass(1) fragments%Gmass(1) = fragments%Gmass(1) + fragments%mass(nlast) diff --git a/src/swiftest/swiftest_discard.f90 b/src/swiftest/swiftest_discard.f90 index 2b69012b7..01a83f5db 100644 --- a/src/swiftest/swiftest_discard.f90 +++ b/src/swiftest/swiftest_discard.f90 @@ -92,6 +92,8 @@ module subroutine swiftest_discard_tp(self, nbody_system, param) logical, dimension(:), allocatable :: ldiscard integer(I4B) :: npl, ntp + if (self%nbody == 0) return + associate(tp => self, cb => nbody_system%cb, pl => nbody_system%pl) ntp = tp%nbody npl = pl%nbody From a5a1b8644b4144b24f03df8d7e84c1eb058b3eef Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 13 Jan 2023 09:50:27 -0500 Subject: [PATCH 3/6] Fixed more bugs and trying new schemes for convergence --- src/collision/collision_util.f90 | 12 ++++++------ src/encounter/encounter_io.f90 | 4 ++-- src/fraggle/fraggle_generate.f90 | 14 +++++++------- 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 728a8ac24..4a089dfa6 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -248,11 +248,11 @@ module subroutine collision_util_get_energy_and_momentum(self, nbody_system, par fragments%L_spin(:,i) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * fragments%rot(:,i) end do call swiftest_util_get_potential_energy(nfrag, [(.true., i = 1, nfrag)], constraint_system%cb%Gmass, fragments%Gmass, fragments%mass, fragments%rb, fragments%pe) - fragments%be = sum(-3*fragments%Gmass(:)*fragments%mass(:)/(5*fragments%radius(:))) - fragments%L_orbit_tot(:) = sum(fragments%L_orbit(:,:),dim=2) - fragments%L_spin_tot(:) = sum(fragments%L_spin(:,:),dim=2) - fragments%ke_orbit_tot = sum(fragments%ke_orbit(:)) - fragments%ke_spin_tot = sum(fragments%ke_spin(:)) + fragments%be = sum(-3*fragments%Gmass(1:nfrag)*fragments%mass(1:nfrag)/(5*fragments%radius(1:nfrag))) + fragments%L_orbit_tot(:) = sum(fragments%L_orbit(:,1:nfrag),dim=2) + fragments%L_spin_tot(:) = sum(fragments%L_spin(:,1:nfrag),dim=2) + fragments%ke_orbit_tot = sum(fragments%ke_orbit(1:nfrag)) + fragments%ke_spin_tot = sum(fragments%ke_spin(1:nfrag)) end if end select @@ -755,7 +755,7 @@ module subroutine collision_util_shift_vector_to_origin(m_frag, vec_frag) mvec_frag(:) = 0.0_DP mtot = sum(m_frag) - nfrag = size(m_frag) + nfrag = count(m_frag > tiny(0.0_DP)) do i = 1, nfrag mvec_frag = mvec_frag(:) + vec_frag(:,i) * m_frag(i) diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index 50180113c..37e784c41 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -40,10 +40,10 @@ module subroutine encounter_io_netcdf_dump(self, param) end if end do + nc%max_tslot = nc%max_tslot + maxval(self%tmap(1:self%iframe)) call nc%close() - call self%reset() ! Update the time slot tracker - nc%max_tslot = nc%max_tslot + maxval(self%tmap(1:self%iframe)) + call self%reset() end if end select diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index ddb14c261..67ad37e82 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -457,9 +457,8 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu real(DP) :: vmin_guess = 1.5_DP real(DP) :: vmax_guess = 10.0_DP real(DP) :: delta_v, volume - integer(I4B), parameter :: MAXLOOP = 20 - integer(I4B), parameter :: MAXTRY = 100 - real(DP), parameter :: SUCCESS_METRIC = 1.0e-3_DP + integer(I4B), parameter :: MAXLOOP = 100 + real(DP), parameter :: SUCCESS_METRIC = 1.0e-2_DP class(collision_fraggle), allocatable :: collider_local character(len=STRMAX) :: message @@ -497,7 +496,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu lfailure = .false. dE_metric = huge(1.0_DP) - outer: do try = 1, nfrag - 1 + outer: do try = 1, nfrag - 2 ! Scale the magnitude of the velocity by the distance from the impact point ! This will reduce the chances of fragments colliding with each other immediately, and is more physically correct do concurrent(i = 1:nfrag) @@ -510,7 +509,6 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu fscale = log(vmax_guess - vmin_guess + 1.0_DP) / log(maxval(vscale(:))) vscale(:) = vscale(:)**fscale + vmin_guess - 1.0_DP - ! Set the velocities of all fragments using all of the scale factors determined above do concurrent(i = 1:nfrag) j = fragments%origin_body(i) @@ -526,7 +524,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu ! Add more velocity dispersion to disruptions vs hit and runs. vmag = vesc * vscale(i) rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) - vimp_unit(:) = .unit. rimp(:) + vimp_unit(:) = .unit. (rimp(:) + 0.5_DP * impactors%bounce_unit(:)) fragments%vc(:,i) = vmag * vscale(i) * vimp_unit(:) + vrot(:) end if end do @@ -541,6 +539,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu L_residual(:) = collider_local%L_total(:,2) - collider_local%L_total(:,1) fragments%L_spin(:,1) = fragments%L_spin(:,1) - L_residual(:) fragments%rot(:,1) = fragments%L_spin(:,1) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(:,1)) + fragments%rotmag(:) = .mag.fragments%rot(:,:) call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) @@ -573,7 +572,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu fragments%vc(:,:) = fscale * fragments%vc(:,:) f_spin = 1.0_DP - f_orbit - ke_remove = min(f_spin * E_residual, fragments%ke_spin_tot) + ke_remove = min(f_spin * E_residual, 0.01_DP*fragments%ke_spin_tot) ke_rot_remove(:) = ke_remove * (fragments%ke_spin(:) / fragments%ke_spin_tot) where(ke_rot_remove(:) > fragments%ke_spin(:)) ke_rot_remove(:) = fragments%ke_spin(:) do concurrent(i = 1:fragments%nbody, fragments%ke_spin(i) > 10*sqrt(tiny(1.0_DP))) @@ -600,6 +599,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu fragments%mass(nlast) = 0.0_DP fragments%Gmass(nlast) = 0.0_DP fragments%radius(nlast) = 0.0_DP + fragments%status(nlast) = INACTIVE fragments%nbody = nlast - 1 call fragments%reset() From 6dc488941a84c3df48750d92e016bdfea547c6ae Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 13 Jan 2023 10:00:10 -0500 Subject: [PATCH 4/6] Trying one more convergence scheme. The behavior still doesn't make sense. It keeps putting all the angular momentum in spin and so never converges on energy --- src/fraggle/fraggle_generate.f90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 67ad37e82..c2de18f73 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -537,8 +537,10 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") ! Check for any residual angular momentum, and if there is any, put it into spin of the largest body L_residual(:) = collider_local%L_total(:,2) - collider_local%L_total(:,1) - fragments%L_spin(:,1) = fragments%L_spin(:,1) - L_residual(:) - fragments%rot(:,1) = fragments%L_spin(:,1) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(:,1)) + do i = 1, fragments%nbody + fragments%L_spin(:,i) = fragments%L_spin(:,i) - L_residual(:) * fragments%mass(i) / fragments%mtot + fragments%rot(:,i) = fragments%L_spin(:,i) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i)) + end do fragments%rotmag(:) = .mag.fragments%rot(:,:) call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") From 71fcd33c787ad412cf8a629a25f9ef19bb1e5866 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 13 Jan 2023 14:29:33 -0500 Subject: [PATCH 5/6] Found a set of criteria that can handle the weird edge case I've been throwing at it. It reduces the number of fragments if there is a failure. --- src/fraggle/fraggle_generate.f90 | 47 +++++++++++++++++++++----------- 1 file changed, 31 insertions(+), 16 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index c2de18f73..846a63b57 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -103,6 +103,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur real(DP), dimension(NDIM) :: dL character(len=STRMAX) :: message real(DP), parameter :: fail_scale_initial = 1.001_DP + integer(I4B) :: nfrag_start ! The minimization and linear solvers can sometimes lead to floating point exceptions. Rather than halting the code entirely if this occurs, we ! can simply fail the attempt and try again. So we need to turn off any floating point exception halting modes temporarily @@ -114,9 +115,10 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur class is (swiftest_nbody_system) select type(param) class is (swiftest_parameters) - associate(impactors => self%impactors, fragments => self%fragments, nfrag => self%fragments%nbody, pl => nbody_system%pl) + associate(impactors => self%impactors, pl => nbody_system%pl) - write(message,*) nfrag + nfrag_start = self%fragments%nbody + write(message,*) nfrag_start call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle generating " // trim(adjustl(message)) // " fragments.") if (param%lflatten_interactions) then @@ -135,6 +137,10 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur call fraggle_generate_rot_vec(self, nbody_system, param) call fraggle_generate_vel_vec(self, nbody_system, param, lfailure_local) + if (self%fragments%nbody /= nfrag_start) then + write(message,*) self%fragments%nbody + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle found a solution with " // trim(adjustl(message)) // " fragments" ) + end if call self%get_energy_and_momentum(nbody_system, param, phase="after") dL = self%L_total(:,2)- self%L_total(:,1) @@ -407,6 +413,7 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param) fragments%rot(:,i) = impactors%rot(:,2) fragments%vc(:,i) = impactors%vc(:,2) end do + fragments%rotmag(:) = .mag.fragments%rot(:,:) call collider%get_energy_and_momentum(nbody_system, param, phase="after") dKE = 0.5_DP * (collider%pe(2) - collider%pe(1) + collider%be(2) - collider%be(1) - impactors%Qloss) KE_final = max(KE_init + dKE,0.0_DP) @@ -416,9 +423,9 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param) Lbefore(:) = mass_init * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (impactors%vb(:,2) - impactors%vb(:,1)) Lafter(:) = mass_final * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (v_final * impactors%bounce_unit(:)) - L_spin(:) = impactors%L_spin(:,1) + (Lbefore(:) - Lafter(:)) - fragments%rot(:,1) = L_spin(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1)) - fragments%rotmag(1) = .mag.fragments%rot(:,1) + L_spin(:) = impactors%L_spin(:,1) + random_scale_factor * (Lbefore(:) - Lafter(:)) + !fragments%rot(:,1) = L_spin(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1)) + !fragments%rotmag(1) = .mag.fragments%rot(:,1) ! Add in some random spin noise. The magnitude will be scaled by the before-after amount and the direction will be random do i = 2,nfrag call random_number(rotdir) @@ -450,7 +457,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu integer(I4B) :: i, j, loop, try, istart, nfrag, nlast logical :: lhitandrun, lsupercat real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual - real(DP) :: vmag, vesc, dE, E_residual, ke_avail, ke_remove, dE_best, E_residual_best, fscale, f_spin, f_orbit, dE_metric + real(DP) :: vmag, vesc, dE, E_residual, ke_min, ke_avail, ke_remove, dE_best, E_residual_best, fscale, f_spin, f_orbit, dE_metric integer(I4B), dimension(collider%fragments%nbody) :: vsign real(DP), dimension(collider%fragments%nbody) :: vscale, ke_rot_remove ! For the initial "guess" of fragment velocities, this is the minimum and maximum velocity relative to escape velocity that the fragments will have @@ -458,6 +465,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu real(DP) :: vmax_guess = 10.0_DP real(DP) :: delta_v, volume integer(I4B), parameter :: MAXLOOP = 100 + integer(I4B), parameter :: MAXTRY = 1000 real(DP), parameter :: SUCCESS_METRIC = 1.0e-2_DP class(collision_fraggle), allocatable :: collider_local character(len=STRMAX) :: message @@ -496,7 +504,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu lfailure = .false. dE_metric = huge(1.0_DP) - outer: do try = 1, nfrag - 2 + outer: do try = 1, maxtry ! Scale the magnitude of the velocity by the distance from the impact point ! This will reduce the chances of fragments colliding with each other immediately, and is more physically correct do concurrent(i = 1:nfrag) @@ -524,7 +532,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu ! Add more velocity dispersion to disruptions vs hit and runs. vmag = vesc * vscale(i) rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) - vimp_unit(:) = .unit. (rimp(:) + 0.5_DP * impactors%bounce_unit(:)) + vimp_unit(:) = .unit. (rimp(:) + vsign(i) * impactors%bounce_unit(:)) fragments%vc(:,i) = vmag * vscale(i) * vimp_unit(:) + vrot(:) end if end do @@ -532,15 +540,22 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu ! Every time the collision-frame velocities are altered, we need to be sure to shift everything back to the center-of-mass frame call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) call fragments%set_coordinate_system() + ke_min = 0.5_DP * fragments%mtot * vesc**2 do loop = 1, MAXLOOP call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") + ke_avail = max(fragments%ke_orbit_tot - ke_min, 0.0_DP) ! Check for any residual angular momentum, and if there is any, put it into spin of the largest body L_residual(:) = collider_local%L_total(:,2) - collider_local%L_total(:,1) - do i = 1, fragments%nbody - fragments%L_spin(:,i) = fragments%L_spin(:,i) - L_residual(:) * fragments%mass(i) / fragments%mtot - fragments%rot(:,i) = fragments%L_spin(:,i) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i)) - end do + if (ke_avail < epsilon(1.0_DP)) then + do i = 1, fragments%nbody + fragments%L_spin(:,i) = fragments%L_spin(:,i) - L_residual(:) * fragments%mass(i) / fragments%mtot + fragments%rot(:,i) = fragments%L_spin(:,i) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i)) + end do + else + fragments%L_spin(:,1) = fragments%L_spin(:,1) - L_residual(:) + fragments%rot(:,1) = fragments%L_spin(:,1) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(:,1)) + end if fragments%rotmag(:) = .mag.fragments%rot(:,:) call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") @@ -562,8 +577,6 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu end if if ((dE_best < 0.0_DP) .and. (dE_metric <= SUCCESS_METRIC * try)) exit outer ! As the tries increase, we relax the success metric. What was once a failure might become a success - ke_avail = fragments%ke_orbit_tot - 0.5_DP * fragments%mtot * vesc**2 - ! Remove a constant amount of velocity from the bodies so we don't shift the center of mass and screw up the momentum f_spin = (fragments%ke_spin_tot )/ (fragments%ke_spin_tot + fragments%ke_orbit_tot) f_orbit = 1.0_DP - f_spin @@ -574,11 +587,12 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu fragments%vc(:,:) = fscale * fragments%vc(:,:) f_spin = 1.0_DP - f_orbit - ke_remove = min(f_spin * E_residual, 0.01_DP*fragments%ke_spin_tot) + ke_remove = min(f_spin * E_residual, 0.9_DP*fragments%ke_spin_tot) ke_rot_remove(:) = ke_remove * (fragments%ke_spin(:) / fragments%ke_spin_tot) where(ke_rot_remove(:) > fragments%ke_spin(:)) ke_rot_remove(:) = fragments%ke_spin(:) do concurrent(i = 1:fragments%nbody, fragments%ke_spin(i) > 10*sqrt(tiny(1.0_DP))) fscale = sqrt((fragments%ke_spin(i) - ke_rot_remove(i))/fragments%ke_spin(i)) + fragments%rotmag(i) = fscale * fragments%rotmag(i) fragments%rot(:,i) = fscale * fragments%rot(:,i) end do @@ -588,13 +602,14 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu end do ! We didn't converge. Reset the fragment positions and velocities and try a new configuration with some slightly different parameters + if (fragments%nbody == 2) exit ! Reduce the number of fragments by one nlast = fragments%nbody fragments%Ip(:,1) = fragments%mass(1) * fragments%Ip(:,1) + fragments%mass(nlast) * fragments%Ip(:,nlast) fragments%mass(1) = fragments%mass(1) + fragments%mass(nlast) fragments%Ip(:,1) = fragments%Ip(:,1) / fragments%mass(1) fragments%Gmass(1) = fragments%Gmass(1) + fragments%mass(nlast) - volume = 4.0_DP / 3.0_DP * PI * (fragments%radius(1) + fragments%radius(nlast)) + volume = 4.0_DP / 3.0_DP * PI * ((fragments%radius(1))**3 + (fragments%radius(nlast))**3) fragments%density(1) = fragments%mass(1) / volume fragments%radius(1) = (3._DP * volume / (4._DP * PI))**(THIRD) fragments%Ip(:,nlast) = 0.0_DP From 420195007a8e2214bfd1be2680269fbe06efd5ea Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 13 Jan 2023 14:33:04 -0500 Subject: [PATCH 6/6] Impvoved the reporting format for Fraggle --- src/fraggle/fraggle_generate.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 846a63b57..e3beb38e1 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -634,12 +634,12 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu end do outer lfailure = dE_best > 0.0_DP + write(message, *) try*loop if (lfailure) then - write(message,*) "Fraggle velocity calculation failed to converge after ",try*loop," steps. This collision may have added energy." + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle velocity calculation failed to converge after " // trim(adjustl(message)) // " steps. This collision would add energy.") else - write(message,*) "Fraggle velocity calculation converged after ",try*loop," steps." + call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Fraggle velocity calculation converged after " // trim(adjustl(message)) // " steps.") end if - call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) end associate end associate