From 33819dbffc286765f63ad0520a5a441e8b8eb8b0 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 30 Dec 2022 13:16:10 -0500 Subject: [PATCH] Tweaked initial position calc to get better convergence in some circumstances. It still fails with the off-axis disruption case --- src/fraggle/fraggle_generate.f90 | 23 ++++++++++++++--------- src/fraggle/fraggle_module.f90 | 3 ++- 2 files changed, 16 insertions(+), 10 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 6baa4b3fa..c3dda504c 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 + integer(I4B) :: try, subtry real(DP) :: dEtot, dLmag integer(I4B), parameter :: MAXTRY = 100 logical :: lk_plpl, lfailure_local @@ -153,9 +153,13 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur ! Compute the "before" energy/momentum and the budgets call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) call self%set_budgets() - call fraggle_generate_pos_vec(self) - call fraggle_generate_rot_vec(self) - call fraggle_generate_vel_vec(self) + 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 + 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)) @@ -309,7 +313,7 @@ module subroutine fraggle_generate_pos_vec(collider) logical :: lcat, lhitandrun integer(I4B), parameter :: MAXLOOP = 10000 real(DP) :: rdistance - real(DP), parameter :: fail_scale = 1.01_DP ! Scale factor to apply to cloud radius and distance if cloud generation fails + real(DP), parameter :: fail_scale = 1.001_DP ! Scale factor to apply to cloud radius and distance if cloud generation fails associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) @@ -319,7 +323,7 @@ module subroutine fraggle_generate_pos_vec(collider) ! We will treat the first two fragments of the list as special cases. ! Place the first two bodies at the centers of the two fragment clouds, but be sure they are sufficiently far apart to avoid overlap rdistance = .mag. (impactors%rc(:,2) - impactors%rc(:,1)) - sum(fragments%radius(1:2)) - rdistance = min(0.5_DP*rdistance, 1e-6_DP*impactors%radius(2)) + rdistance = 0.5_DP*rdistance fragment_cloud_radius(:) = impactors%radius(:) @@ -420,7 +424,7 @@ module subroutine fraggle_generate_rot_vec(collider) end subroutine fraggle_generate_rot_vec - module subroutine fraggle_generate_vel_vec(collider) + module subroutine fraggle_generate_vel_vec(collider, lfailure) !! Author: David A. Minton !! !! Generates an initial velocity distribution. For disruptions, the velocity magnitude is set to be @@ -429,15 +433,15 @@ module subroutine fraggle_generate_vel_vec(collider) implicit none ! Arguments 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, iend, n, ndof + integer(I4B) :: i, j, loop, istart, n, ndof logical :: lhitandrun, lnoncat real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, Lresidual real(DP) :: vmag, vesc, rotmag, ke_residual, ke_per_dof, ke_tot integer(I4B), dimension(collider%fragments%nbody) :: vsign real(DP), dimension(collider%fragments%nbody) :: vscale, mass_vscale, ke_avail integer(I4B), parameter :: MAXLOOP = 100 - real(DP), parameter :: TOL = 1e-1 associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) @@ -541,6 +545,7 @@ module subroutine fraggle_generate_vel_vec(collider) fragments%rot(:,i) = fragments%rot(:,i) + Lresidual(:) / (nfrag*fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i)) end do end do + lfailure = ke_residual < 0.0_DP do concurrent(i = 1:nfrag) fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:) diff --git a/src/fraggle/fraggle_module.f90 b/src/fraggle/fraggle_module.f90 index 130b3af5c..9022a0a0c 100644 --- a/src/fraggle/fraggle_module.f90 +++ b/src/fraggle/fraggle_module.f90 @@ -82,9 +82,10 @@ module subroutine fraggle_generate_rot_vec(collider) class(collision_fraggle), intent(inout) :: collider !! Collision system object end subroutine fraggle_generate_rot_vec - module subroutine fraggle_generate_vel_vec(collider) + module subroutine fraggle_generate_vel_vec(collider, lfailure) implicit none class(collision_fraggle), intent(inout) :: collider !! Collision system object + logical, intent(out) :: lfailure !! Did the velocity computation fail? end subroutine fraggle_generate_vel_vec module subroutine fraggle_util_setup_fragments_system(self, nfrag)