From 28e5cde369f9653c67ab19a5c646a37dc5b13535 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 31 Dec 2022 12:23:19 -0500 Subject: [PATCH] More tweaks to find the best fragment distribution for the tricky cases --- src/fraggle/fraggle_generate.f90 | 121 ++++++++++++++++--------------- 1 file changed, 64 insertions(+), 57 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index c0afb2393..8e434ed6a 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -370,13 +370,14 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure) class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object logical, intent(out) :: lfailure !! Did the velocity computation fail? ! Internals - integer(I4B) :: i, j, loop, istart, n, ndof + integer(I4B) :: i, j, loop, try, istart, n, ndof logical :: lhitandrun, lsupercat real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, Lresidual, vshear, vunit real(DP) :: vmag, vesc, rotmag, ke_residual, ke_per_dof, ke_tot, ke_residual_min integer(I4B), dimension(collider%fragments%nbody) :: vsign real(DP), dimension(collider%fragments%nbody) :: vscale, mass_vscale, ke_avail - integer(I4B), parameter :: MAXLOOP = 1000 + integer(I4B), parameter :: MAXLOOP = 100 + integer(I4B), parameter :: MAXTRY = 100 real(DP), parameter :: TOL = 1e-6 class(collision_fragments(:)), allocatable :: fragments @@ -442,68 +443,74 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure) end if call fragments%set_coordinate_system() ke_residual_min = -huge(1.0_DP) - do loop = 1, MAXLOOP - call fragments%get_kinetic_energy() - ke_residual = fragments%ke_budget - (fragments%ke_orbit_tot + fragments%ke_spin_tot) - if ((abs(ke_residual) < abs(ke_residual_min)) .or. ((ke_residual >= 0.0_DP) .and. (ke_residual_min < 0.0_DP))) then ! This is our best case so far. Save it for posterity - if (allocated(collider%fragments)) deallocate(collider%fragments) - allocate(collider%fragments, source=fragments) - ke_residual_min = ke_residual - if ((ke_residual > 0.0_DP) .and. (ke_residual < TOL * fragments%ke_budget)) exit - end if - ! Make sure we don't take away too much orbital kinetic energy, otherwise the fragment can't escape - ke_avail(:) = fragments%ke_orbit(:) - impactors%Gmass(1)*impactors%mass(2)/fragments%rmag(:) - ke_tot = 0.0_DP - ke_per_dof = -ke_residual - do i = 1, 2*(nfrag - istart + 1) - n = count(ke_avail(istart:nfrag) > -ke_residual/i) - if (ke_residual < 0.0_DP) n = n + count(fragments%ke_spin(istart:nfrag) > -ke_residual/i) - if (abs(n * ke_per_dof) > ke_tot) then - ke_per_dof = -ke_residual/i - ke_tot = n * ke_per_dof - ndof = i - if (abs(ke_tot) > abs(ke_residual)) then - ke_tot = -ke_residual - ke_per_dof = ke_tot/n - exit - end if + outer: do try = 1, MAXTRY + do loop = 1, MAXLOOP + call fragments%get_kinetic_energy() + ke_residual = fragments%ke_budget - (fragments%ke_orbit_tot + fragments%ke_spin_tot) + if ((abs(ke_residual) < abs(ke_residual_min)) .or. ((ke_residual >= 0.0_DP) .and. (ke_residual_min < 0.0_DP))) then ! This is our best case so far. Save it for posterity + if (allocated(collider%fragments)) deallocate(collider%fragments) + allocate(collider%fragments, source=fragments) + ke_residual_min = ke_residual + if ((ke_residual > 0.0_DP) .and. (ke_residual < TOL * fragments%ke_budget)) exit outer end if - end do - do concurrent(i = istart:nfrag, ke_avail(i) > ke_per_dof) - vmag = max(fragments%vmag(i)**2 - 2*ke_per_dof/fragments%mass(i),vesc**2) - fragments%vmag(i) = sqrt(vmag) - fragments%vc(:,i) = fragments%vmag(i) * .unit.fragments%vc(:,i) - end do - do concurrent(i = istart:nfrag, fragments%ke_spin(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 = max(rotmag, 0.0_DP) - fragments%rotmag(i) = sqrt(rotmag) - fragments%rot(:,i) = fragments%rotmag(i) * .unit.fragments%rot(:,i) - end do + ! Make sure we don't take away too much orbital kinetic energy, otherwise the fragment can't escape + ke_avail(:) = fragments%ke_orbit(:) - impactors%Gmass(1)*impactors%mass(2)/fragments%rmag(:) + ke_tot = 0.0_DP + ke_per_dof = -ke_residual + do i = 1, 2*(nfrag - istart + 1) + n = count(ke_avail(istart:nfrag) > -ke_residual/i) + if (ke_residual < 0.0_DP) n = n + count(fragments%ke_spin(istart:nfrag) > -ke_residual/i) + if (abs(n * ke_per_dof) > ke_tot) then + ke_per_dof = -ke_residual/i + ke_tot = n * ke_per_dof + ndof = i + if (abs(ke_tot) > abs(ke_residual)) then + ke_tot = -ke_residual + ke_per_dof = ke_tot/n + exit + end if + end if + end do + do concurrent(i = istart:nfrag, ke_avail(i) > ke_per_dof) + vmag = max(fragments%vmag(i)**2 - 2*ke_per_dof/fragments%mass(i),vesc**2) + fragments%vmag(i) = sqrt(vmag) + fragments%vc(:,i) = fragments%vmag(i) * .unit.fragments%vc(:,i) + end do + do concurrent(i = istart:nfrag, fragments%ke_spin(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 = max(rotmag, 0.0_DP) + fragments%rotmag(i) = sqrt(rotmag) + fragments%rot(:,i) = fragments%rotmag(i) * .unit.fragments%rot(:,i) + end do - call fragments%set_coordinate_system() - if (lsupercat) then - ! Put some of the residual angular momentum into velocity shear. Not too much, or we get some weird trajectories + call fragments%set_coordinate_system() + if (lsupercat) then + ! Put some of the residual angular momentum into velocity shear. Not too much, or we get some weird trajectories + call fragments%get_angular_momentum() + Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:)) + do concurrent(i = istart:nfrag) + vunit(:) = .unit. (Lresidual(:) .cross. fragments%r_unit(:,i)) + vshear(:) = vunit(:) * (.mag.Lresidual(:) / ((nfrag-istart+1)*fragments%mass(i) * fragments%rmag(i))) + fragments%vc(:,i) = fragments%vc(:,i) + vshear(:) + end do + end if + ! Check for any residual angular momentum, and if there is any, put it into spin call fragments%get_angular_momentum() Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:)) - do concurrent(i = istart:nfrag) - vunit(:) = .unit. (Lresidual(:) .cross. fragments%r_unit(:,i)) - vshear(:) = vunit(:) * (.mag.Lresidual(:) / ((nfrag-istart+1)*fragments%mass(i) * fragments%rmag(i))) - fragments%vc(:,i) = fragments%vc(:,i) + vshear(:) + do concurrent(i = 1:nfrag) + fragments%Lspin(:,i) = fragments%Lspin(:,i) + Lresidual(:) / nfrag + fragments%rot(:,i) = fragments%Lspin(:,i) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i)) end do - end if - ! Check for any residual angular momentum, and if there is any, put it into spin - call fragments%get_angular_momentum() - Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:)) - do concurrent(i = 1:nfrag) - fragments%Lspin(:,i) = fragments%Lspin(:,i) + Lresidual(:) / nfrag - fragments%rot(:,i) = fragments%Lspin(:,i) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i)) - end do - call fragments%get_angular_momentum() - Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:)) + call fragments%get_angular_momentum() + Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:)) - end do + end do + ! We didn't converge. Try another configuration and see if we get a better result + call fraggle_generate_pos_vec(collider) + call fraggle_generate_rot_vec(collider) + collider%fail_scale = collider%fail_scale*1.01_DP + end do outer lfailure = ke_residual < 0.0_DP do concurrent(i = 1:nfrag)