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

Commit

Permalink
Added collisional regime determination code
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Aug 5, 2021
1 parent b4f1cb7 commit 4fc17fd
Show file tree
Hide file tree
Showing 4 changed files with 109 additions and 0 deletions.
6 changes: 6 additions & 0 deletions Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,6 +96,11 @@ lib:
ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
ln -s $(SWIFTEST_HOME)/Makefile .; \
make libdir
cd $(SWIFTEST_HOME)/src/fragmentation; \
rm -f Makefile.Defines Makefile; \
ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
ln -s $(SWIFTEST_HOME)/Makefile .; \
make libdir
cd $(SWIFTEST_HOME)/src/gr; \
rm -f Makefile.Defines Makefile; \
ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
Expand Down Expand Up @@ -189,6 +194,7 @@ clean:
cd $(SWIFTEST_HOME)/src/discard; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/drift; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/eucl; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/fragmentation; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/gr; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/helio; rm -f Makefile.Defines Makefile *.gc*
cd $(SWIFTEST_HOME)/src/io; rm -f Makefile.Defines Makefile *.gc*
Expand Down
32 changes: 32 additions & 0 deletions src/fragmentation/fragmentation.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,32 @@
submodule(swiftest_classes) s_fragmentation
use swiftest
contains
module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss)
!! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
!!
!! Determine the collisional regime of two colliding bodies.
!! Current version requires all values to be converted to SI units prior to calling the function
!! References:
!! Kokubo, E., Genda, H., 2010. Formation of Terrestrial Planets from Protoplanets Under a Realistic Accretion
!! Condition. ApJL 714, L21. https://doi.org/10.1088/2041-8205/714/1/L21
!! Leinhardt, Z.M., Stewart, S.T., 2012. Collisions between Gravity-dominated Bodies. I. Outcome Regimes and Scaling
!! Laws 745, 79. https://doi.org/10.1088/0004-637X/745/1/79
!! Mustill, A.J., Davies, M.B., Johansen, A., 2018. The dynamical evolution of transiting planetary systems including
!! a realistic collision prescription. Mon Not R Astron Soc 478, 2896–2908. https://doi.org/10.1093/mnras/sty1273
!! Rufu, R., Aharonson, O., 2019. Impact Dynamics of Moons Within a Planetary Potential. J. Geophys. Res. Planets 124,
!! 1008–1019. https://doi.org/10.1029/2018JE005798
!! Stewart, S.T., Leinhardt, Z.M., 2012. Collisions between Gravity-dominated Bodies. II. The Diversity of Impact
!! Outcomes during the End Stage of Planet Formation. ApJ 751, 32. https://doi.org/10.1088/0004-637X/751/1/32
!!
implicit none
! Arguments
integer(I4B), intent(out) :: regime
real(DP), intent(out) :: Mlr, Mslr
real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny
real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2
real(DP), intent(out) :: Qloss !! The residual energy after the collision

return
end subroutine fragmentation_regime

end submodule s_fragmentation
9 changes: 9 additions & 0 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -452,6 +452,15 @@ module subroutine eucl_dist_index_plpl(self)
class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
end subroutine

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
real(DP), intent(out) :: Mlr, Mslr
real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny
real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2
real(DP), intent(out) :: Qloss !! The residual energy after the collision
end subroutine fragmentation_regime

module pure subroutine gr_kick_getaccb_ns_body(self, system, param)
implicit none
class(swiftest_body), intent(inout) :: self !! Swiftest generic body object
Expand Down
62 changes: 62 additions & 0 deletions src/symba/symba_collision.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -408,6 +408,68 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param)
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions
! Internals
! Internals
integer(I4B), dimension(:), allocatable :: family !! List of indices of all bodies inovlved in the collision
integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision
real(DP), dimension(NDIM,2) :: x, v, L_spin, Ip !! Output values that represent a 2-body equivalent of a possibly 2+ body collision
real(DP), dimension(2) :: mass, radius !! Output values that represent a 2-body equivalent of a possibly 2+ body collision
logical :: lgoodcollision
integer(I4B) :: i, status, jtarg, jproj, regime
real(DP), dimension(2) :: radius_si, mass_si, density_si
real(DP) :: mtiny_si, Mcb_si
real(DP), dimension(NDIM) :: x1_si, v1_si, x2_si, v2_si
real(DP) :: mlr, mslr, mtot, dentot, msys, msys_new, Qloss, impact_parameter
integer(I4B), parameter :: NRES = 3 !! Number of collisional product results
real(DP), dimension(NRES) :: mass_res

associate(plpl_collisions => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, cb => system%cb)
select type(pl => system%pl)
class is (symba_pl)
do i = 1, ncollisions
idx_parent(1) = pl%kin(idx1(i))%parent
idx_parent(2) = pl%kin(idx2(i))%parent
lgoodcollision = symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v, mass, radius, L_spin, Ip)
if (.not. lgoodcollision) cycle
if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved

! Convert all quantities to SI units and determine which of the pair is the projectile vs. target before sending them
! to symba_regime
if (mass(1) > mass(2)) then
jtarg = 1
jproj = 2
else
jtarg = 2
jproj = 1
end if
mass_si(:) = (mass(:)) * param%MU2KG !! The collective mass of the parent and its children
radius_si(:) = radius(:) * param%DU2M !! The collective radius of the parent and its children
x1_si(:) = plpl_collisions%x1(:,i) * param%DU2M !! The position of the parent from inside the step (at collision)
v1_si(:) = plpl_collisions%v1(:,i) * param%DU2M / param%TU2S !! The velocity of the parent from inside the step (at collision)
x2_si(:) = plpl_collisions%x2(:,i) * param%DU2M !! The position of the parent from inside the step (at collision)
v2_si(:) = plpl_collisions%v2(:,i) * param%DU2M / param%TU2S !! The velocity of the parent from inside the step (at collision)
density_si(:) = mass_si(:) / (4.0_DP / 3._DP * PI * radius_si(:)**3) !! The collective density of the parent and its children
Mcb_si = cb%mass * param%MU2KG
mtiny_si = (param%MTINY / param%GU) * param%MU2KG

mass_res(:) = 0.0_DP

mtot = sum(mass_si(:))
dentot = sum(mass_si(:) * density_si(:)) / mtot

!! Use the positions and velocities of the parents from indside the step (at collision) to calculate the collisional regime
call fragmentation_regime(Mcb_si, mass_si(jtarg), mass_si(jproj), radius_si(jtarg), radius_si(jproj), x1_si(:), x2_si(:),&
v1_si(:), v2_si(:), density_si(jtarg), density_si(jproj), regime, mlr, mslr, mtiny_si, Qloss)

mass_res(1) = min(max(mlr, 0.0_DP), mtot)
mass_res(2) = min(max(mslr, 0.0_DP), mtot)
mass_res(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot)
mass_res(:) = (mass_res(:) / param%MU2KG) * param%GU
Qloss = Qloss * (param%GU / param%MU2KG) * (param%TU2S / param%DU2M)**2

!status = symba_fragmentation_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip)
end do
end select
end associate

return
end subroutine symba_collision_resolve_fragmentations
Expand Down

0 comments on commit 4fc17fd

Please sign in to comment.