diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 84f8936f0..807f22419 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -72,7 +72,7 @@ rot_vectors = {"disruption_headon" : [np.array([0.0, 0.0, 0.0]), np.array([0.0, 0.0, 0.0])], - "disruption_off_axis": [np.array([0.0, 0.0, 6.0e4]), + "disruption_off_axis": [np.array([0.0, 0.0, -6.0e4]), np.array([0.0, 0.0, 1.0e5])], "supercatastrophic_headon": [np.array([0.0, 0.0, 0.0]), np.array([0.0, 0.0, 0.0])], diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 53900f353..202cc5297 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -155,7 +155,7 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) ! Internals integer(I4B) :: i, j, k, ibiggest real(DP), dimension(NDIM) :: Lspin_new - real(DP) :: volume, dpe + real(DP) :: volume character(len=STRMAX) :: message select type(nbody_system) @@ -208,11 +208,6 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) ! Get the energy of the system after the collision call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) - ! Keep track of the component of potential energy that is now not considered because two bodies became one - dpe = self%pe(2) - self%pe(1) - nbody_system%Ecollisions = nbody_system%Ecollisions - dpe - nbody_system%Euntracked = nbody_system%Euntracked + dpe - ! Update any encounter lists that have the removed bodies in them so that they instead point to the new body do k = 1, nbody_system%plpl_encounter%nenc do j = 1, impactors%ncoll diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 5fe3893cb..23a0d20b4 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -89,36 +89,38 @@ module collision !> Class definition for the variables that describe a collection of fragments in barycentric coordinates type, extends(base_multibody) :: collision_fragments - real(DP) :: mtot !! Total mass of fragments - class(base_particle_info), dimension(:), allocatable :: info !! Particle metadata information - integer(I4B), dimension(nbody) :: status !! An integrator-specific status indicator - real(DP), dimension(NDIM,nbody) :: rh !! Heliocentric position - real(DP), dimension(NDIM,nbody) :: vh !! Heliocentric velocity - real(DP), dimension(NDIM,nbody) :: rb !! Barycentric position - real(DP), dimension(NDIM,nbody) :: vb !! Barycentric velocity - real(DP), dimension(NDIM,nbody) :: rot !! rotation vectors of fragments - real(DP), dimension(NDIM,nbody) :: Ip !! Principal axes moment of inertia for fragments - real(DP), dimension(nbody) :: mass !! masses of fragments - real(DP), dimension(nbody) :: radius !! Radii of fragments - real(DP), dimension(nbody) :: density !! Radii of fragments - real(DP), dimension(NDIM,nbody) :: rc !! Position vectors in the collision coordinate frame - real(DP), dimension(NDIM,nbody) :: vc !! Velocity vectors in the collision coordinate frame - real(DP), dimension(nbody) :: rmag !! Array of radial distance magnitudes of individual fragments in the collisional coordinate frame - real(DP), dimension(nbody) :: vmag !! Array of radial distance magnitudes of individual fragments in the collisional coordinate frame - real(DP), dimension(nbody) :: rotmag !! Array of rotation magnitudes of individual fragments - real(DP), dimension(NDIM,nbody) :: r_unit !! Array of radial direction unit vectors of individual fragments in the collisional coordinate frame - real(DP), dimension(NDIM,nbody) :: v_unit !! Array of velocity direction unit vectors of individual fragments in the collisional coordinate frame - real(DP), dimension(NDIM,nbody) :: t_unit !! Array of tangential direction unit vectors of individual fragments in the collisional coordinate frame - real(DP), dimension(NDIM,nbody) :: n_unit !! Array of normal direction unit vectors of individual fragments in the collisional coordinate frame - integer(I1B), dimension(nbody) :: origin_body !! Array of indices indicating which impactor body (1 or 2) the fragment originates from - real(DP), dimension(NDIM) :: Lorbit !! Orbital angular momentum vector of all fragments - real(DP), dimension(NDIM) :: Lspin !! Spin angular momentum vector of all fragments - real(DP) :: ke_orbit !! Orbital kinetic energy of all fragments - real(DP) :: ke_spin !! Spin kinetic energy of all fragments - real(DP) :: ke_budget !! Kinetic energy budget for computing fragment trajectories - real(DP), dimension(NDIM) :: L_budget !! Angular momentum budget for computing fragment trajectories - real(DP), dimension(nbody) :: ke_orbit_frag !! Orbital kinetic energy of each individual fragment - real(DP), dimension(nbody) :: ke_spin_frag !! Spin kinetic energy of each individual fragment + real(DP) :: mtot !! Total mass of fragments + class(base_particle_info), dimension(:), allocatable :: info !! Particle metadata information + integer(I4B), dimension(nbody) :: status !! An integrator-specific status indicator + real(DP), dimension(NDIM,nbody) :: rh !! Heliocentric position + real(DP), dimension(NDIM,nbody) :: vh !! Heliocentric velocity + real(DP), dimension(NDIM,nbody) :: rb !! Barycentric position + real(DP), dimension(NDIM,nbody) :: vb !! Barycentric velocity + real(DP), dimension(NDIM,nbody) :: rot !! rotation vectors of fragments + real(DP), dimension(NDIM,nbody) :: Ip !! Principal axes moment of inertia for fragments + real(DP), dimension(nbody) :: mass !! masses of fragments + real(DP), dimension(nbody) :: radius !! Radii of fragments + real(DP), dimension(nbody) :: density !! Radii of fragments + real(DP), dimension(NDIM,nbody) :: rc !! Position vectors in the collision coordinate frame + real(DP), dimension(NDIM,nbody) :: vc !! Velocity vectors in the collision coordinate frame + real(DP), dimension(nbody) :: rmag !! Array of radial distance magnitudes of individual fragments in the collisional coordinate frame + real(DP), dimension(nbody) :: vmag !! Array of radial distance magnitudes of individual fragments in the collisional coordinate frame + real(DP), dimension(nbody) :: rotmag !! Array of rotation magnitudes of individual fragments + real(DP), dimension(NDIM,nbody) :: r_unit !! Array of radial direction unit vectors of individual fragments in the collisional coordinate frame + real(DP), dimension(NDIM,nbody) :: v_unit !! Array of velocity direction unit vectors of individual fragments in the collisional coordinate frame + real(DP), dimension(NDIM,nbody) :: t_unit !! Array of tangential direction unit vectors of individual fragments in the collisional coordinate frame + real(DP), dimension(NDIM,nbody) :: n_unit !! Array of normal direction unit vectors of individual fragments in the collisional coordinate frame + integer(I1B), dimension(nbody) :: origin_body !! Array of indices indicating which impactor body (1 or 2) the fragment originates from + real(DP), dimension(NDIM) :: Lorbit_tot !! Orbital angular momentum vector of all fragments + real(DP), dimension(NDIM) :: Lspin_tot !! Spin angular momentum vector of all fragments + real(DP), dimension(NDIM,nbody) :: Lorbit !! Orbital angular momentum vector of each individual fragment + real(DP), dimension(NDIM,nbody) :: Lspin !! Spin angular momentum vector of each individual fragment + real(DP) :: ke_orbit_tot !! Orbital kinetic energy of all fragments + real(DP) :: ke_spin_tot !! Spin kinetic energy of all fragments + real(DP) :: ke_budget !! Kinetic energy budget for computing fragment trajectories + real(DP), dimension(NDIM) :: L_budget !! Angular momentum budget for computing fragment trajectories + real(DP), dimension(nbody) :: ke_orbit !! Orbital kinetic energy of each individual fragment + real(DP), dimension(nbody) :: ke_spin !! Spin kinetic energy of each individual fragment contains procedure :: reset => collision_util_reset_fragments !! Deallocates all allocatable arrays and sets everything else to 0 procedure :: get_angular_momentum => collision_util_get_angular_momentum !! Calcualtes the current angular momentum of the fragments diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index 9e7e1e468..b0fff79dd 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -497,6 +497,7 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) logical :: lgoodcollision integer(I4B) :: i, loop, ncollisions integer(I4B), parameter :: MAXCASCADE = 1000 + real(DP) :: dpe select type (nbody_system) class is (swiftest_nbody_system) @@ -537,11 +538,23 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) idx_parent(2) = pl%kin(idx2(i))%parent call impactors%consolidate(nbody_system, param, idx_parent, lgoodcollision) if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLIDED)) cycle - call impactors%get_regime(nbody_system, param) + call collision_history%take_snapshot(param,nbody_system, t, "before") + + call nbody_system%get_energy_and_momentum(param) + collider%pe(1) = nbody_system%pe + call collider%generate(nbody_system, param, t) + + call nbody_system%get_energy_and_momentum(param) + collider%pe(2) = nbody_system%pe + dpe = collider%pe(2) - collider%pe(1) + nbody_system%Ecollisions = nbody_system%Ecollisions - dpe + nbody_system%Euntracked = nbody_system%Euntracked + dpe + call collision_history%take_snapshot(param,nbody_system, t, "after") + plpl_collision%status(i) = collider%status call impactors%reset() end do diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 9573c1452..2878f95ec 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -186,13 +186,14 @@ module subroutine collision_util_get_angular_momentum(self) integer(I4B) :: i associate(fragments => self, nfrag => self%nbody) - fragments%Lorbit(:) = 0.0_DP - fragments%Lspin(:) = 0.0_DP do i = 1, nfrag - fragments%Lorbit(:) = fragments%Lorbit(:) + fragments%mass(i) * (fragments%rc(:, i) .cross. fragments%vc(:, i)) - fragments%Lspin(:) = fragments%Lspin(:) + fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:, i) * fragments%rot(:, i) + fragments%Lorbit(:,i) = fragments%mass(i) * (fragments%rc(:,i) .cross. fragments%vc(:, i)) + fragments%Lspin(:,i) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i) * fragments%rot(:,i) end do + + fragments%Lorbit_tot(:) = sum(fragments%Lorbit, dim=2) + fragments%Lspin_tot(:) = sum(fragments%Lspin, dim=2) end associate return @@ -210,18 +211,16 @@ module subroutine collision_util_get_kinetic_energy(self) integer(I4B) :: i associate(fragments => self, nfrag => self%nbody) - fragments%ke_orbit_frag(:) = 0.0_DP - fragments%ke_spin_frag(:) = 0.0_DP do concurrent(i = 1:nfrag) - fragments%ke_orbit_frag(i) = fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i)) - fragments%ke_spin_frag(i) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * dot_product(fragments%rot(:,i),fragments%rot(:,i) ) + fragments%ke_orbit(i) = fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i)) + fragments%ke_spin(i) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * dot_product(fragments%rot(:,i),fragments%rot(:,i) ) end do - fragments%ke_orbit_frag(:) = fragments%ke_orbit_frag(:) / 2 - fragments%ke_spin_frag(:) = fragments%ke_spin_frag(:) / 2 - fragments%ke_orbit = sum(fragments%ke_orbit_frag(:)) - fragments%ke_spin = sum(fragments%ke_spin_frag(:)) + fragments%ke_orbit(:) = fragments%ke_orbit(:) / 2 + fragments%ke_spin(:) = fragments%ke_spin(:) / 2 + fragments%ke_orbit_tot = sum(fragments%ke_orbit(:)) + fragments%ke_spin_tot = sum(fragments%ke_spin(:)) end associate @@ -266,14 +265,15 @@ module subroutine collision_util_get_energy_momentum(self, nbody_system, param, end do ke_orbit = ke_orbit / 2 ke_spin = ke_spin / 2 + else call fragments%get_angular_momentum() - Lorbit(:) = fragments%Lorbit(:) - Lspin(:) = fragments%Lspin(:) + Lorbit(:) = fragments%Lorbit_tot(:) + Lspin(:) = fragments%Lspin_tot(:) call fragments%get_kinetic_energy() - ke_orbit = fragments%ke_orbit - ke_spin = fragments%ke_spin + ke_orbit = fragments%ke_orbit_tot + ke_spin = fragments%ke_spin_tot end if ! Calculate the current fragment energy and momentum balances @@ -459,7 +459,7 @@ module subroutine collision_util_set_coordinate_collider(self) fragments%r_unit(:,:) = .unit. fragments%rc(:,:) fragments%v_unit(:,:) = .unit. fragments%vc(:,:) fragments%n_unit(:,:) = .unit. (fragments%rc(:,:) .cross. fragments%vc(:,:)) - fragments%t_unit(:,:) = .unit. (fragments%r_unit(:,:) .cross. fragments%n_unit(:,:)) + fragments%t_unit(:,:) = -.unit. (fragments%r_unit(:,:) .cross. fragments%n_unit(:,:)) end associate @@ -594,11 +594,11 @@ module subroutine collision_util_setup_fragments_collider(self, nfrag) self%fragments%density(:) = 0.0_DP self%fragments%rmag(:) = 0.0_DP self%fragments%vmag(:) = 0.0_DP - self%fragments%Lorbit(:) = 0.0_DP - self%fragments%Lspin(:) = 0.0_DP + self%fragments%Lorbit_tot(:) = 0.0_DP + self%fragments%Lspin_tot(:) = 0.0_DP self%fragments%L_budget(:) = 0.0_DP - self%fragments%ke_orbit = 0.0_DP - self%fragments%ke_spin = 0.0_DP + self%fragments%ke_orbit_tot = 0.0_DP + self%fragments%ke_spin_tot = 0.0_DP self%fragments%ke_budget = 0.0_DP return diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 7e72c8f90..6b93500a5 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -30,7 +30,6 @@ module subroutine fraggle_generate(self, nbody_system, param, t) ! Internals integer(I4B) :: i, ibiggest, nfrag character(len=STRMAX) :: message - real(DP) :: dpe select type(nbody_system) class is (swiftest_nbody_system) @@ -38,6 +37,7 @@ module subroutine fraggle_generate(self, nbody_system, param, t) class is (swiftest_pl) associate(impactors => self%impactors, status => self%status) + select case (impactors%regime) case (COLLRESOLVE_REGIME_HIT_AND_RUN) call self%hitandrun(nbody_system, param, t) @@ -56,9 +56,6 @@ module subroutine fraggle_generate(self, nbody_system, param, t) call self%set_mass_dist(param) call self%disrupt(nbody_system, param, t) - dpe = self%pe(2) - self%pe(1) - nbody_system%Ecollisions = nbody_system%Ecollisions - dpe - nbody_system%Euntracked = nbody_system%Euntracked + dpe associate (fragments => self%fragments) ! Populate the list of new bodies @@ -79,6 +76,7 @@ module subroutine fraggle_generate(self, nbody_system, param, t) end select call collision_resolve_mergeaddsub(nbody_system, param, t, status) + end associate end associate end select @@ -221,17 +219,16 @@ module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t) ! Result integer(I4B) :: status !! Status flag assigned to this outcome ! Internals - integer(I4B) :: i, ibiggest, nfrag, jtarg, jproj + integer(I4B) :: i, ibiggest, jtarg, jproj, nfrag logical :: lpure character(len=STRMAX) :: message - real(DP) :: dpe select type(nbody_system) class is (swiftest_nbody_system) select type(pl => nbody_system%pl) class is (swiftest_pl) - associate(impactors => self%impactors, nfrag => self%fragments%nbody) + associate(impactors => self%impactors) 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))) @@ -244,35 +241,23 @@ module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t) jproj = 1 end if - ! The frag disruption model (and its extended types allow for non-pure hit and run. + ! The Fraggle disruption model (and its extended types allow for non-pure hit and run. if (impactors%mass_dist(2) > 0.9_DP * impactors%mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Pure hit and run. No new fragments generated.") nfrag = 0 call self%collision_basic%hitandrun(nbody_system, param, t) lpure = .true. return - else ! Imperfect hit and run, so we'll keep the largest body and destroy the other - lpure = .false. - call self%set_mass_dist(param) - - ! Generate the position and velocity distributions of the fragments - call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) - call self%disrupt(nbody_system, param, t, lpure) - call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) - - dpe = self%pe(2) - self%pe(1) - nbody_system%Ecollisions = nbody_system%Ecollisions - dpe - nbody_system%Euntracked = nbody_system%Euntracked + dpe - - if (lpure) then - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Should have been a pure hit and run instead") - nfrag = 0 - else - nfrag = self%fragments%nbody - write(message, *) nfrag - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") - end if end if + lpure = .false. + call self%set_mass_dist(param) + + ! Generate the position and velocity distributions of the fragments + call self%disrupt(nbody_system, param, t, lpure) + nfrag = self%fragments%nbody + + write(message, *) nfrag + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") ibiggest = impactors%id(maxloc(pl%Gmass(impactors%id(:)), dim=1)) self%fragments%id(1) = pl%id(ibiggest) @@ -431,11 +416,12 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure) ! Internals integer(I4B) :: i, j, loop, istart, n, ndof logical :: lhitandrun, lnoncat - real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, Lresidual + real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, Lresidual, vshear, vunit 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 + integer(I4B), parameter :: MAXLOOP = 1000 + 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) @@ -467,7 +453,7 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure) ! Give the fragment velocities a random value that is scaled with fragment mass call random_number(mass_vscale) mass_vscale(:) = (mass_vscale(:) + 1.0_DP) / 2 - mass_vscale(:) = mass_vscale(:) * (fragments%mtot / fragments%mass(:))**(0.125_DP) ! The 1/8 power is arbitrary. It just gives the velocity a small mass dependence + mass_vscale(:) = mass_vscale(:) * (fragments%mtot / fragments%mass(:))**(0.125_DP) ! The power is arbitrary. It just gives the velocity a small mass dependence mass_vscale(:) = mass_vscale(:) / minval(mass_vscale(:)) ! Set the velocities of all fragments using all of the scale factors determined above @@ -498,15 +484,15 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure) do loop = 1, MAXLOOP call collider%set_coordinate_system() call fragments%get_kinetic_energy() - ke_residual = fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin) - if (ke_residual > 0.0_DP) exit + ke_residual = fragments%ke_budget - (fragments%ke_orbit_tot + fragments%ke_spin_tot) + if ((ke_residual > 0.0_DP).and.(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%rmag(:) + 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_frag(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 @@ -523,8 +509,8 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure) 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 - ke_per_dof/(fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) + 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) @@ -533,11 +519,22 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure) ! 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(:)) - 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)) + 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 + call fragments%get_angular_momentum() + Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:)) + do concurrent(i = istart:nfrag) + fragments%Lspin(:,i) = fragments%Lspin(:,i) + Lresidual(:) / (nfrag-istart+1) + 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(:)) + end do lfailure = ke_residual < 0.0_DP diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 5e3997ab7..e59ed130b 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -229,12 +229,15 @@ module subroutine fraggle_util_set_natural_scale_factors(self) class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object ! Internals integer(I4B) :: i + real(DP) :: vesc associate(collider => self, fragments => self%fragments, impactors => self%impactors) ! Set primary scale factors (mass, length, and time) based on the impactor properties at the time of collision - collider%mscale = fragments%mtot - collider%dscale = sum(impactors%radius(:)) - collider%tscale = collider%dscale / (.mag.(impactors%vc(:,2) - impactors%vc(:,1))) + collider%mscale = minval(fragments%mass(:)) + collider%dscale = minval(fragments%radius(:)) + + vesc = sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:))) + collider%tscale = collider%dscale / vesc ! Set secondary scale factors for convenience collider%vscale = collider%dscale / collider%tscale diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index a0ac7bc85..a492bdda6 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -211,7 +211,7 @@ module swiftest real(DP) :: dR = 0.0_DP !! Change in the radius of the central body contains procedure :: read_in => swiftest_io_read_in_cb !! Read in central body initial conditions from an ASCII file - procedure :: write_frame => swiftest_io_netcdf_write_frame_cb !! I/O routine for writing out a single frame of time-series data for all bodies in the nbody_system in NetCDF format + procedure :: write_frame => swiftest_io_netcdf_write_frame_cb !! I/O routine for writing out a single frame of time-series data for all bodies in the system in NetCDF format procedure :: write_info => swiftest_io_netcdf_write_info_cb !! Dump contents of particle information metadata to file end type swiftest_cb