diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index f020fa36f..754a65b10 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -987,7 +987,7 @@ end subroutine restructure_failed_fragments end subroutine fragmentation_initialize - module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss) + module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, min_mfrag, Qloss) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Determine the collisional regime of two colliding bodies. @@ -1008,7 +1008,7 @@ module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, v ! 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), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, min_mfrag real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 real(DP), intent(out) :: Qloss !! The residual energy after the collision ! Constants @@ -1082,7 +1082,7 @@ module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, v Qloss = 0.0_DP U_binding = (3.0_DP * Mtot) / (5.0_DP * Rp) ! LS12 eq. 27 - if ((m1 < mtiny).or.(m2 < mtiny)) then + if ((m1 < min_mfrag).or.(m2 < min_mfrag)) then regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime Mlr = Mtot Mslr = 0.0_DP diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 6a85bba3e..bcacef235 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -474,11 +474,11 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? end subroutine fragmentation_initialize - module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss) + module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, min_mfrag, 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), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, min_mfrag real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 real(DP), intent(out) :: Qloss !! Energy lost during the collision end subroutine fragmentation_regime diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 4ed739651..53997f033 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -89,6 +89,7 @@ module symba_classes procedure :: discard => symba_discard_pl !! Process massive body discards procedure :: drift => symba_drift_pl !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level procedure :: encounter_check => symba_encounter_check_pl !! Checks if massive bodies are going through close encounters with each other + procedure :: accel_int => symba_kick_getacch_int_pl !! Compute direct cross (third) term heliocentric accelerations of massive bodiess, with no mutual interactions between bodies below GMTINY procedure :: accel => symba_kick_getacch_pl !! Compute heliocentric accelerations of massive bodies procedure :: setup => symba_setup_pl !! Constructor method - Allocates space for the input number of bodies procedure :: append => symba_util_append_pl !! Appends elements from one structure to another @@ -405,6 +406,11 @@ module subroutine symba_io_read_particle(system, param) class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA extensions end subroutine symba_io_read_particle + module pure subroutine symba_kick_getacch_int_pl(self) + implicit none + class(symba_pl), intent(inout) :: self + end subroutine symba_kick_getacch_int_pl + module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 4ca907fe6..0e37a307e 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -939,7 +939,7 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) logical :: lgoodcollision integer(I4B) :: i, jtarg, jproj, regime real(DP), dimension(2) :: radius_si, mass_si, density_si - real(DP) :: mtiny_si, Mcb_si + real(DP) :: min_mfrag_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 @@ -974,7 +974,7 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) v2_si(:) = plplcollision_list%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%GMTINY / param%GU) * param%MU2KG + min_mfrag_si = (param%min_GMfrag / param%GU) * param%MU2KG mass_res(:) = 0.0_DP @@ -983,7 +983,7 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) !! 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) + v1_si(:), v2_si(:), density_si(jtarg), density_si(jproj), regime, mlr, mslr, min_mfrag_si, Qloss) mass_res(1) = min(max(mlr, 0.0_DP), mtot) mass_res(2) = min(max(mslr, 0.0_DP), mtot) diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 372a40996..ff98fe1ce 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -23,11 +23,11 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc if (self%nbody == 0) return - associate(pl => self, npl => self%nbody, nplpl => self%nplpl) - allocate(lencounter(nplpl), loc_lvdotr(nplpl)) + associate(pl => self, npl => self%nbody, nplplm => self%nplplm) + allocate(lencounter(nplplm), loc_lvdotr(nplplm)) lencounter(:) = .false. - do k = 1, nplpl + do k = 1, nplplm associate(i => pl%k_plpl(1, k), j => pl%k_plpl(2, k)) xr(:) = pl%xh(:, j) - pl%xh(:, i) vr(:) = pl%vh(:, j) - pl%vh(:, i) diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index b3449e04e..bf7adc60c 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -2,6 +2,50 @@ use swiftest contains + module pure subroutine symba_kick_getacch_int_pl(self) + !! author: David A. Minton + !! + !! Compute direct cross (third) term heliocentric accelerations of massive bodies, with no mutual interactions between bodies below GMTINY + !! + !! Adapted from Hal Levison's Swift routine symba5_helio_getacch.f + !! Adapted from David E. Kaufmann's Swifter routine helio_kick_getacch_int.f90 + implicit none + ! Arguments + class(symba_pl), intent(inout) :: self + ! Internals + integer(I4B) :: k + real(DP) :: rji2, irij3, faci, facj, rlim2 + real(DP) :: dx, dy, dz + + associate(pl => self, npl => self%nbody, nplplm => self%nplplm) + do k = 1, nplplm + associate(i => pl%k_plpl(1, k), j => pl%k_plpl(2, k)) + if (pl%lmask(i) .and. pl%lmask(j)) then + dx = pl%xh(1, j) - pl%xh(1, i) + dy = pl%xh(2, j) - pl%xh(2, i) + dz = pl%xh(3, j) - pl%xh(3, i) + rji2 = dx**2 + dy**2 + dz**2 + rlim2 = (pl%radius(i) + pl%radius(j))**2 + if (rji2 > rlim2) then + irij3 = 1.0_DP / (rji2 * sqrt(rji2)) + faci = pl%Gmass(i) * irij3 + facj = pl%Gmass(j) * irij3 + pl%ah(1, i) = pl%ah(1, i) + facj * dx + pl%ah(2, i) = pl%ah(2, i) + facj * dy + pl%ah(3, i) = pl%ah(3, i) + facj * dz + pl%ah(1, j) = pl%ah(1, j) - faci * dx + pl%ah(2, j) = pl%ah(2, j) - faci * dy + pl%ah(3, j) = pl%ah(3, j) - faci * dz + end if + end if + end associate + end do + end associate + + return + end subroutine symba_kick_getacch_int_pl + + module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) !! author: David A. Minton !! diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 05b4d9aa9..c5796ccc8 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -346,10 +346,11 @@ module subroutine symba_util_index_eucl_plpl(self, param) class is (symba_parameters) pl%lmtiny(1:npl) = pl%Gmass(1:npl) < param%GMTINY end select - pl%nplm = count(.not. pl%lmtiny(1:npl)) + nplm = count(.not. pl%lmtiny(1:npl)) + pl%nplm = int(nplm, kind=I4B) nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl, minus first column - nplplm = (nplm * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl, minus first column + nplplm = nplm * npl - nplm * (nplm + 1) / 2 ! number of entries in a strict lower triangle, npl x npl, minus first column including only mutually interacting bodies if (allocated(self%k_plpl)) deallocate(self%k_plpl) ! Reset the index array if it's been set previously allocate(self%k_plpl(2, nplpl)) do i = 1, npl