diff --git a/Makefile b/Makefile index 63cfb0ee0..162129149 100644 --- a/Makefile +++ b/Makefile @@ -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 .; \ @@ -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* diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 new file mode 100644 index 000000000..d809bc03a --- /dev/null +++ b/src/fragmentation/fragmentation.f90 @@ -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 \ No newline at end of file diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 10e60a491..af751fdcf 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -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 diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index ad0e64079..843de7dd6 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -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