diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index d17191912..84f8936f0 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -152,8 +152,7 @@ def __init__(self, sim, animfile, title, style, nskip=1): self.fig, self.ax = self.setup_plot() # Then setup FuncAnimation. - self.ani = animation.FuncAnimation(self.fig, self.update_plot, interval=1, frames=range(0,nframes,nskip), - blit=True) + self.ani = animation.FuncAnimation(self.fig, self.update_plot, interval=1, frames=range(0,nframes,nskip), blit=True) self.ani.save(animfile, fps=60, dpi=300, extra_args=['-vcodec', 'libx264']) print(f"Finished writing {animfile}") diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 8c025f026..c5b435fb5 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -133,7 +133,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur end if call ieee_set_flag(ieee_all, .false.) ! Set all fpe flags to quiet - !call self%set_natural_scale() + call self%set_natural_scale() call fragments%reset() @@ -157,11 +157,10 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur call fraggle_generate_rot_vec(self) call fraggle_generate_vel_vec(self) call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) - exit dEtot = self%Etot(2) - self%Etot(1) dLmag = .mag. (self%Ltot(:,2) - self%Ltot(:,1)) - lfailure_local = (dEtot > FRAGGLE_ETOL) + 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: " // & @@ -195,7 +194,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur trim(adjustl(message)) // " tries") end if - !call self%set_original_scale() + call self%set_original_scale() ! Restore the big array if (lk_plpl) call pl%flatten(param) @@ -234,7 +233,7 @@ module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t) class is (swiftest_nbody_system) select type(pl => nbody_system%pl) class is (swiftest_pl) - associate(impactors => self%impactors) + associate(impactors => self%impactors, nfrag => self%fragments%nbody) message = "Hit and run between" call collision_io_collider_message(nbody_system%pl, impactors%id, message) call swiftest_io_log_one_message(COLLISION_LOG_OUT, trim(adjustl(message))) @@ -310,7 +309,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.1_DP ! Scale factor to apply to cloud radius and distance if cloud generation fails + real(DP), parameter :: fail_scale = 1.01_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) @@ -431,13 +430,14 @@ module subroutine fraggle_generate_vel_vec(collider) ! Arguments class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object ! Internals - integer(I4B) :: i,j, loop, istart + 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 + 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-4 associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) @@ -501,30 +501,44 @@ module subroutine fraggle_generate_vel_vec(collider) call collider%set_coordinate_system() call fragments%get_kinetic_energy() ke_residual = fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin) + if (abs(ke_residual) <= fragments%ke_budget * TOL) exit ! Make sure we don't take away too much orbital kinetic energy, otherwise the fragment can't escape - ke_avail(:) = fragments%ke_orbit_frag(:) - impactors%Gmass(1)*impactors%mass(2)/fragments%vmag(:) - ke_per_dof = -ke_residual/((nfrag - istart + 1)) + ke_avail(:) = fragments%ke_orbit_frag(:) - 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) + count(fragments%ke_spin_frag(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) - fragments%vmag(i) = sqrt(2 * (fragments%ke_orbit_frag(i) - ke_per_dof)/fragments%mass(i)) - fragments%vc(:,i) = fragments%vmag(i) * fragments%v_unit(:,i) + 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_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 = max(rotmag, 0.0_DP) + fragments%rotmag(i) = sqrt(rotmag) + fragments%rot(:,i) = rotmag * .unit.fragments%rot(:,i) end do - ! do concurrent(i = istart:nfrag, fragments%ke_spin_frag(i) > ke_per_dof) - ! fragments%rotmag(i) = sqrt(2 * (fragments%ke_spin_frag(i) - ke_per_dof)/(fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i))) - ! fragments%rot(:,i) = fragments%rotmag(i) * .unit.fragments%rot(:,i) - ! end do - call fragments%get_kinetic_energy() - ke_residual = fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin) - ! Check for any residual angular momentum, and if there is any, put it into spin of the biggest body + ! Check for any residual angular momentum, and if there is any, put it into spin call collider%set_coordinate_system() call fragments%get_angular_momentum() Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit(:) + fragments%Lspin(:)) - rotmag = .mag. fragments%rot(:,1) - fragments%rot(:,1) = fragments%rot(:,1) + Lresidual(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(:,1)) - rotmag = .mag. fragments%rot(:,1) - - if (ke_residual >= 0.0_DP) exit - + do concurrent(i = 1:nfrag) + rotmag = .mag. fragments%rot(:,i) + fragments%rot(:,i) = fragments%rot(:,i) + Lresidual(:) / (nfrag*fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i)) + end do end do do concurrent(i = 1:nfrag) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 9e519584c..5e3997ab7 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -231,39 +231,40 @@ module subroutine fraggle_util_set_natural_scale_factors(self) integer(I4B) :: i associate(collider => self, fragments => self%fragments, impactors => self%impactors) - ! Set scale factors - collider%Escale = 0.5_DP * ( impactors%mass(1) * dot_product(impactors%vb(:,1), impactors%vb(:,1)) & - + impactors%mass(2) * dot_product(impactors%vb(:,2), impactors%vb(:,2))) - collider%dscale = sum(impactors%radius(:)) + ! Set primary scale factors (mass, length, and time) based on the impactor properties at the time of collision collider%mscale = fragments%mtot - collider%vscale = sqrt(collider%Escale / collider%mscale) - collider%tscale = collider%dscale / collider%vscale + collider%dscale = sum(impactors%radius(:)) + collider%tscale = collider%dscale / (.mag.(impactors%vc(:,2) - impactors%vc(:,1))) + + ! Set secondary scale factors for convenience + collider%vscale = collider%dscale / collider%tscale + collider%Escale = collider%mscale * collider%vscale**2 collider%Lscale = collider%mscale * collider%dscale * collider%vscale ! Scale all dimensioned quantities of impactors and fragments - impactors%rbcom(:) = impactors%rbcom(:) / collider%dscale - impactors%vbcom(:) = impactors%vbcom(:) / collider%vscale - impactors%rbimp(:) = impactors%rbimp(:) / collider%dscale - impactors%bounce_unit(:) = impactors%bounce_unit(:) / collider%vscale - impactors%rb(:,:) = impactors%rb(:,:) / collider%dscale - impactors%vb(:,:) = impactors%vb(:,:) / collider%vscale - impactors%rc(:,:) = impactors%rc(:,:) / collider%dscale - impactors%vc(:,:) = impactors%vc(:,:) / collider%vscale - impactors%mass(:) = impactors%mass(:) / collider%mscale - impactors%Gmass(:) = impactors%Gmass(:) / 0.5_DP * collider%vscale**2 * collider%dscale - impactors%Mcb = impactors%Mcb / collider%mscale - impactors%radius(:) = impactors%radius(:) / collider%dscale - impactors%Lspin(:,:) = impactors%Lspin(:,:) / collider%Lscale - impactors%Lorbit(:,:) = impactors%Lorbit(:,:) / collider%Lscale + impactors%rbcom(:) = impactors%rbcom(:) / collider%dscale + impactors%vbcom(:) = impactors%vbcom(:) / collider%vscale + impactors%rbimp(:) = impactors%rbimp(:) / collider%dscale + impactors%rb(:,:) = impactors%rb(:,:) / collider%dscale + impactors%vb(:,:) = impactors%vb(:,:) / collider%vscale + impactors%rc(:,:) = impactors%rc(:,:) / collider%dscale + impactors%vc(:,:) = impactors%vc(:,:) / collider%vscale + impactors%mass(:) = impactors%mass(:) / collider%mscale + impactors%Gmass(:) = impactors%Gmass(:) / (collider%dscale**3/collider%tscale**2) + impactors%Mcb = impactors%Mcb / collider%mscale + impactors%radius(:) = impactors%radius(:) / collider%dscale + impactors%Lspin(:,:) = impactors%Lspin(:,:) / collider%Lscale + impactors%Lorbit(:,:) = impactors%Lorbit(:,:) / collider%Lscale + impactors%bounce_unit(:) = impactors%bounce_unit(:) / collider%vscale do i = 1, 2 - impactors%rot(:,i) = impactors%Lspin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i)) + impactors%rot(:,i) = impactors%Lspin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i)) end do - fragments%mtot = fragments%mtot / collider%mscale - fragments%mass = fragments%mass / collider%mscale - fragments%radius = fragments%radius / collider%dscale - impactors%Qloss = impactors%Qloss / collider%Escale + fragments%mtot = fragments%mtot / collider%mscale + fragments%mass = fragments%mass / collider%mscale + fragments%radius = fragments%radius / collider%dscale + impactors%Qloss = impactors%Qloss / collider%Escale end associate return @@ -288,13 +289,13 @@ module subroutine fraggle_util_set_original_scale_factors(self) associate(collider => self, fragments => self%fragments, impactors => self%impactors) ! Restore scale factors - impactors%rbcom(:) = impactors%rbcom(:) * collider%dscale - impactors%vbcom(:) = impactors%vbcom(:) * collider%vscale - impactors%rbimp(:) = impactors%rbimp(:) * collider%dscale + impactors%rbcom(:) = impactors%rbcom(:) * collider%dscale + impactors%vbcom(:) = impactors%vbcom(:) * collider%vscale + impactors%rbimp(:) = impactors%rbimp(:) * collider%dscale impactors%bounce_unit(:) = impactors%bounce_unit(:) * collider%vscale - + impactors%mass = impactors%mass * collider%mscale - impactors%Gmass(:) = impactors%Gmass(:) * 0.5_DP * collider%vscale**2 * collider%dscale + impactors%Gmass(:) = impactors%Gmass(:) * (collider%dscale**3/collider%tscale**2) impactors%Mcb = impactors%Mcb * collider%mscale impactors%mass_dist = impactors%mass_dist * collider%mscale impactors%radius = impactors%radius * collider%dscale @@ -305,7 +306,7 @@ module subroutine fraggle_util_set_original_scale_factors(self) impactors%Lspin = impactors%Lspin * collider%Lscale impactors%Lorbit = impactors%Lorbit * collider%Lscale do i = 1, 2 - impactors%rot(:,i) = impactors%Lspin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i)) + impactors%rot(:,i) = impactors%Lspin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i)) end do fragments%mtot = fragments%mtot * collider%mscale