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()