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

Commit

Permalink
Enabled the disruption and supercat fragmentation and got all of the …
Browse files Browse the repository at this point in the history 
…interfaces straightened out
  • Loading branch information
@daminton
daminton committed Aug 5, 2021
1 parent 8b4d523 commit d0e7fb4
Show file tree
Hide file tree
Showing 3 changed files with 125 additions and 4 deletions.
2 changes: 1 addition & 1 deletion src/fragmentation/fragmentation.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
use swiftest
contains

subroutine fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, &
module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, &
nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure)
!! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
!!
Expand Down
16 changes: 16 additions & 0 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -453,6 +453,22 @@ module subroutine eucl_dist_index_plpl(self)
class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
end subroutine

module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, &
nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure)
implicit none
class(swiftest_nbody_system), intent(inout) :: system
class(swiftest_parameters), intent(in) :: param
integer(I4B), dimension(:), intent(in) :: family
real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip
real(DP), dimension(:), intent(inout) :: mass, radius
integer(I4B), intent(inout) :: nfrag
real(DP), dimension(:), allocatable, intent(inout) :: m_frag, rad_frag
real(DP), dimension(:,:), allocatable, intent(inout) :: Ip_frag
real(DP), dimension(:,:), allocatable, intent(inout) :: xb_frag, vb_frag, rot_frag
logical, intent(out) :: lfailure ! Answers the question: Should this have been a merger instead?
real(DP), intent(inout) :: Qloss
end subroutine fragmentation_initialize

module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss)
implicit none
integer(I4B), intent(out) :: regime
Expand Down
111 changes: 108 additions & 3 deletions src/symba/symba_fragmentation.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,8 +73,8 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v,
Ip_frag(:, i) = Ip_new(:)
end do

!call fragmentation_initialize(pl, param, family, x, v, L_spin, Ip, mass, radius, &
! nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure)
call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, &
nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure)

if (lfailure) then
write(*,*) 'No fragment solution found, so treat as a pure hit-and-run'
Expand All @@ -83,6 +83,7 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v,
else
! Populate the list of new bodies
write(*,'("Generating ",I2.0," fragments")') nfrag
status = DISRUPTION
allocate(plnew, mold=pl)
call plnew%setup(nfrag, param)

Expand Down Expand Up @@ -309,8 +310,112 @@ module function symba_fragmentation_casesupercatastrophic(system, param, family,
! Result
integer(I4B) :: status !! Status flag assigned to this outcome
! Internals
integer(I4B) :: i, nfrag, ibiggest, nstart, nend
real(DP) :: mtot, avg_dens, min_frag_mass
real(DP), dimension(NDIM) :: xcom, vcom
real(DP), dimension(2) :: vol
real(DP), dimension(NDIM) :: Ip_new
real(DP), dimension(:, :), allocatable :: vb_frag, xb_frag, rot_frag, Ip_frag
real(DP), dimension(:), allocatable :: m_frag, rad_frag
logical :: lfailure
class(symba_pl), allocatable :: plnew

select type(pl => system%pl)
class is (symba_pl)
associate(mergeadd_list => system%mergeadd_list, mergesub_list => system%mergesub_list, cb => system%cb)
! Collisional fragments will be uniformly distributed around the pre-impact barycenter
nfrag = NFRAG_SUPERCAT
allocate(m_frag(nfrag))
allocate(rad_frag(nfrag))
allocate(xb_frag(NDIM, nfrag))
allocate(vb_frag(NDIM, nfrag))
allocate(rot_frag(NDIM, nfrag))
allocate(Ip_frag(NDIM, nfrag))

mtot = sum(mass(:))
xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot
vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot

! Get mass weighted mean of Ip and average density
Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mtot
vol(:) = 4._DP / 3._DP * pi * radius(:)**3
avg_dens = mtot / sum(vol(:))

! If we are adding the first and largest fragment (lr), check to see if its mass is SMALLER than an equal distribution of
! mass between all fragments. If so, we will just distribute the mass equally between the fragments
min_frag_mass = mtot / nfrag
if (mass_res(1) < min_frag_mass) then
m_frag(:) = min_frag_mass
else
m_frag(1) = mass_res(1)
m_frag(2:nfrag) = (mtot - mass_res(1)) / (nfrag - 1)
end if
! Distribute any residual mass if there is any and set the radius
m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:)))
rad_frag(:) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP)

do i = 1, nfrag
Ip_frag(:, i) = Ip_new(:)
end do

status = SUPERCATASTROPHIC
call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, &
nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure)

if (lfailure) then
write(*,*) 'No fragment solution found, so treat as a pure hit-and-run'
status = ACTIVE
nfrag = 0
else
! Populate the list of new bodies
write(*,'("Generating ",I2.0," fragments")') nfrag
status = SUPERCATASTROPHIC
allocate(plnew, mold=pl)
call plnew%setup(nfrag, param)

plnew%id(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)]
system%maxid = system%maxid + nfrag
plnew%status(:) = ACTIVE
plnew%lcollision(:) = .false.
plnew%ldiscard(:) = .false.
plnew%xb(:,:) = xb_frag(:, :)
plnew%vb(:,:) = vb_frag(:, :)
do i = 1, nfrag
plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:)
plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:)
end do
plnew%mass(:) = m_frag(:)
plnew%Gmass(:) = param%GU * m_frag(:)
plnew%density(:) = avg_dens
plnew%radius(:) = rad_frag(:)
plnew%info(:)%origin_type = "Disruption"
plnew%info(:)%origin_time = param%t
do i = 1, nfrag
plnew%info(i)%origin_xh(:) = plnew%xh(:,i)
plnew%info(i)%origin_vh(:) = plnew%vh(:,i)
end do
if (param%lrotation) then
plnew%Ip(:,:) = Ip_frag(:,:)
plnew%rot(:,:) = rot_frag(:,:)
end if
if (param%ltides) then
ibiggest = maxloc(pl%Gmass(family(:)), dim=1)
plnew%Q = pl%Q(ibiggest)
plnew%k2 = pl%k2(ibiggest)
plnew%tlag = pl%tlag(ibiggest)
end if

! Append the new merged body to the list and record how many we made
nstart = mergeadd_list%nbody + 1
nend = mergeadd_list%nbody + plnew%nbody
call mergeadd_list%append(plnew)
mergeadd_list%ncomp(nstart:nend) = plnew%nbody

call plnew%setup(0, param)
deallocate(plnew)
end if

end associate
end select

return
end function symba_fragmentation_casesupercatastrophic
Expand Down

0 comments on commit d0e7fb4

Please sign in to comment.