Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
More bug fixes and improvements to all models except the off-axis dis…
Browse files Browse the repository at this point in the history 
…ruption
  • Loading branch information
@daminton
daminton committed Dec 31, 2022
1 parent c6c76d6 commit fd89184
Show file tree
Hide file tree
Showing 8 changed files with 115 additions and 105 deletions.
2 changes: 1 addition & 1 deletion examples/Fragmentation/Fragmentation_Movie.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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])],
Expand Down
7 changes: 1 addition & 6 deletions src/collision/collision_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
62 changes: 32 additions & 30 deletions src/collision/collision_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
15 changes: 14 additions & 1 deletion src/collision/collision_resolve.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
42 changes: 21 additions & 21 deletions src/collision/collision_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down
Loading

0 comments on commit fd89184

Please sign in to comment.