diff --git a/examples/symba_energy_momentum/param.disruption_headon.in b/examples/symba_energy_momentum/param.disruption_headon.in index 4d187d26b..24e79be5d 100644 --- a/examples/symba_energy_momentum/param.disruption_headon.in +++ b/examples/symba_energy_momentum/param.disruption_headon.in @@ -21,7 +21,8 @@ ENC_OUT enc.disruption_headon.dat EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded RHILL_PRESENT yes ! Hill's sphere radii in input file -GMTINY 1.0e-16 +GMTINY 1.0e-16 +MIN_GMFRAG 1.0e-10 FRAGMENTATION yes MU2KG 1.98908e30 TU2S 3.1556925e7 diff --git a/examples/symba_energy_momentum/param.disruption_off_axis.in b/examples/symba_energy_momentum/param.disruption_off_axis.in index b4ba7f8ba..58922133f 100644 --- a/examples/symba_energy_momentum/param.disruption_off_axis.in +++ b/examples/symba_energy_momentum/param.disruption_off_axis.in @@ -22,7 +22,8 @@ DISCARD_OUT discard.disruption_off_axis.out EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded RHILL_PRESENT yes ! Hill's sphere radii in input file -GMTINY 1.0e-16 +GMTINY 1.0e-16 +MIN_GMFRAG 1.0e-12 FRAGMENTATION yes MU2KG 1.98908e30 TU2S 3.1556925e7 diff --git a/examples/symba_energy_momentum/param.supercatastrophic_headon.in b/examples/symba_energy_momentum/param.supercatastrophic_headon.in index aaa1db99d..1effdd15b 100644 --- a/examples/symba_energy_momentum/param.supercatastrophic_headon.in +++ b/examples/symba_energy_momentum/param.supercatastrophic_headon.in @@ -17,12 +17,11 @@ CHK_RMIN 0.005 CHK_RMAX 1e6 CHK_EJECT -1.0 ! ignore this check CHK_QMIN -1.0 ! ignore this check -DISCARD_OUT discard.supercatastrophic_headon.out -ENC_OUT enc.supercatastrophic_headon.dat EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded RHILL_PRESENT yes ! Hill's sphere radii in input file -GMTINY 1.0e-16 +GMTINY 1.0e-16 +MIN_GMFRAG 1.0e-11 FRAGMENTATION yes MU2KG 1.98908e30 TU2S 3.1556925e7 diff --git a/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in b/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in index 0bf836be5..47bce3775 100644 --- a/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in +++ b/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in @@ -21,7 +21,8 @@ ENC_OUT enc.supercatastrophic_off_axis.dat EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded RHILL_PRESENT yes ! Hill's sphere radii in input file -GMTINY 1.0e-16 +GMTINY 1.0e-16 +MIN_GMFRAG 1.0e-11 FRAGMENTATION yes MU2KG 1.98908e30 TU2S 3.1556925e7 diff --git a/src/fraggle/fraggle_regime.f90 b/src/fraggle/fraggle_regime.f90 index a20b08d29..4ce4699fa 100644 --- a/src/fraggle/fraggle_regime.f90 +++ b/src/fraggle/fraggle_regime.f90 @@ -56,6 +56,7 @@ module subroutine fraggle_regime_colliders(self, frag, system, param) frag%mass_dist(1) = min(max(mlr, 0.0_DP), mtot) frag%mass_dist(2) = min(max(mslr, 0.0_DP), mtot) frag%mass_dist(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot) + frag%mtot = sum(frag%mass_dist(1:3)) ! Convert quantities back to the system units and save them into the fragment system frag%mass_dist(:) = (frag%mass_dist(:) / param%MU2KG) @@ -98,6 +99,7 @@ subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb real(DP), parameter :: DENSITY1 = 1000.0_DP !standard density parameter from LS12 [kg/m3] real(DP), parameter :: MU_BAR = 0.37_DP !0.385#0.37#0.3333# 3.978 # 1/3 material parameter for hydrodynamic planet-size bodies (LS12) real(DP), parameter :: BETA = 2.85_DP !slope of sfd for remnants from LS12 2.85 + real(DP), parameter :: ETA = -1.50_DP !! LS12 eq. (44) real(DP), parameter :: C1 = 2.43_DP !! Kokubo & Genda (2010) eq. (3) real(DP), parameter :: C2 = -0.0408_DP !! Kokubo & Genda (2010) eq. (3) real(DP), parameter :: C3 = 1.86_DP !! Kokubo & Genda (2010) eq. (3) @@ -108,7 +110,7 @@ subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb real(DP) :: a1, alpha, aint, b, bcrit, c_star, egy, zeta, l, lint, mu, phi, theta real(DP) :: Qr, Qrd_pstar, Qr_erosion, Qr_supercat real(DP) :: Vhr, Verosion, Vescp, Vhill, Vimp, Vsupercat - real(DP) :: Mint, Mtot + real(DP) :: Mint, Mtot, Mtmp real(DP) :: Rp, rhill real(DP) :: Mresidual real(DP) :: U_binding @@ -157,7 +159,7 @@ subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb Qr = mu*(Vimp**2) / Mtot / 2.0_DP !calculate mass largest remnant Mlr - Mlr = (1.0_DP - Qr / Qrd_pstar / 2.0_DP) * Mtot ! [kg] # LS12 eq (5) + Mlr = max((1.0_DP - Qr / Qrd_pstar / 2.0_DP) * Mtot, min_mfrag) ! [kg] # LS12 eq (5) !calculate Vsupercat Qr_supercat = SUPERCAT_QRATIO * Qrd_pstar ! See LS12 Section 4.1 @@ -196,7 +198,7 @@ subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb Qloss = 0.0_DP else Mlr = m1 - Mslr = calc_Qrd_rev(m2, m1, Mint, den1, den2, Vimp, c_star) + Mslr = max(calc_Qrd_rev(m2, m1, Mint, den1, den2, Vimp, c_star), min_mfrag) regime = COLLRESOLVE_REGIME_HIT_AND_RUN !hit and run Qloss = (c_star + 1.0_DP) * U_binding ! Qr end if @@ -207,13 +209,13 @@ subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb Mslr = 0.0_DP Qloss = 0.0_DP else - Mslr = Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA) ! LS12 eq (37) + Mslr = max(Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA), min_mfrag) ! LS12 eq (37) regime = COLLRESOLVE_REGIME_DISRUPTION !disruption Qloss = (c_star + 1.0_DP) * U_binding ! Qr - Qr_erosion end if else if (Vimp > Vsupercat) then - Mlr = Mtot * 0.1_DP * (Qr / (Qrd_pstar * SUPERCAT_QRATIO))**(-1.5_DP) !LS12 eq (44) - Mslr = Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA) !LS12 eq (37) + Mlr = max(Mtot * 0.1_DP * (Qr / (Qrd_pstar * SUPERCAT_QRATIO))**(ETA), min_mfrag) !LS12 eq (44) + Mslr = max(Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA), min_mfrag) !LS12 eq (37) regime = COLLRESOLVE_REGIME_SUPERCATASTROPHIC ! supercatastrophic Qloss = (c_star + 1.0_DP) * U_binding ! Qr - Qr_supercat else @@ -221,6 +223,12 @@ subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb end if end if + if (Mslr > Mlr) then ! The second-largest fragment is actually larger than the largest, so we will swap them + Mtmp = Mlr + Mlr = Mslr + Mslr = Mtmp + end if + Mresidual = Mtot - Mlr - Mslr if (Mresidual < 0.0_DP) then ! prevents final masses from going negative Mlr = Mlr + Mresidual diff --git a/src/fraggle/fraggle_set.f90 b/src/fraggle/fraggle_set.f90 index 7ea47e14f..f9238f5f8 100644 --- a/src/fraggle/fraggle_set.f90 +++ b/src/fraggle/fraggle_set.f90 @@ -27,60 +27,106 @@ module subroutine fraggle_set_budgets_fragments(self, colliders) end subroutine fraggle_set_budgets_fragments - module subroutine fraggle_set_mass_dist_fragments(self, colliders) + module subroutine fraggle_set_mass_dist_fragments(self, colliders, param) !! author: David A. Minton !! !! Sets the mass of fragments based on the mass distribution returned by the regime calculation. !! This subroutine must be run after the the setup rourtine has been run on the fragments implicit none ! Arguments - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters ! Internals - integer(I4B) :: i, jproj, jtarg + integer(I4B) :: i, jproj, jtarg, nfrag, istart real(DP), dimension(2) :: volume real(DP), dimension(NDIM) :: Ip_avg + real(DP) :: mfrag, mremaining, min_mfrag + real(DP), parameter :: BETA = 2.85_DP + integer(I4B), parameter :: NFRAGMAX = 100 !! Maximum number of fragments that can be generated + integer(I4B), parameter :: NFRAGMIN = 7 !! Minimum number of fragments that can be generated (set by the fraggle_generate algorithm for constraining momentum and energy) + integer(I4B), parameter :: NFRAG_SIZE_MULTIPLIER = 3 !! Log-space scale factor that scales the number of fragments by the collisional system mass - associate(frag => self, nfrag => self%nbody) + associate(frag => self) ! Get mass weighted mean of Ip and density volume(1:2) = 4._DP / 3._DP * PI * colliders%radius(1:2)**3 - frag%density(1:nfrag) = frag%mtot / sum(volume(:)) Ip_avg(:) = (colliders%mass(1) * colliders%Ip(:,1) + colliders%mass(2) * colliders%Ip(:,2)) / frag%mtot - do i = 1, nfrag - frag%Ip(:, i) = Ip_avg(:) - end do + if (colliders%mass(1) > colliders%mass(2)) then + jtarg = 1 + jproj = 2 + else + jtarg = 2 + jproj = 1 + end if + + select type(param) + class is (symba_parameters) + min_mfrag = (param%min_GMfrag / param%GU) + nfrag = ceiling(NFRAG_SIZE_MULTIPLIER * log(frag%mtot / min_mfrag)) + nfrag = max(min(nfrag, NFRAGMAX), NFRAGMIN) + class default + min_mfrag = 0.0_DP + nfrag = NFRAGMAX + end select select case(frag%regime) case(COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - ! Assign the first and second-largest masses that came out of the regime determination subroutine to the first two fragments - frag%mass(1:2) = frag%mass_dist(1:2) + istart = 2 + case(COLLRESOLVE_REGIME_HIT_AND_RUN) + istart = 1 + case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) + call frag%setup(1, param) + frag%mass(1) = frag%mass_dist(1) + frag%radius(1) = colliders%radius(jtarg) + frag%density(1) = frag%mass_dist(1) / volume(jtarg) + frag%Ip(:, 1) = colliders%Ip(:,1) + return + case default + write(*,*) "fraggle_set_mass_dist_fragments error: Unrecognized regime code",frag%regime + end select - ! Distribute remaining mass among the remaining fragments - do i = 3, nfrag - frag%mass(i) = (frag%mtot - sum(frag%mass(1:2))) / (nfrag - 2) - end do - frag%radius(1:nfrag) = (3 * frag%mass(1:nfrag) / (4 * PI * frag%density(1:nfrag)))**(1.0_DP / 3.0_DP) + i = istart + mremaining = max(frag%mtot - sum(frag%mass_dist(1:istart)), 0.0_DP) + do while (i <= NFRAGMAX) + mfrag = (i - istart + 1)**(-3._DP / BETA) * frag%mass_dist(istart) + if (mremaining - mfrag < 0.0_DP) exit + mremaining = mremaining - mfrag + i = i + 1 + end do + + call frag%setup(nfrag, param) + frag%mass(1:istart) = frag%mass_dist(1:istart) + mremaining = max(frag%mtot - sum(frag%mass_dist(1:istart)), 0.0_DP) + do i = istart + 1, nfrag + mfrag = (i - istart + 1)**(-3._DP / BETA) * frag%mass_dist(istart) + mfrag = min(mfrag, mremaining) + frag%mass(i) = mfrag + mremaining = mremaining - mfrag + end do + select case(frag%regime) case(COLLRESOLVE_REGIME_HIT_AND_RUN) - if (colliders%mass(1) > colliders%mass(2)) then - jtarg = 1 - jproj = 2 - else - jtarg = 2 - jproj = 1 - end if frag%mass(1) = frag%mass_dist(1) frag%radius(1) = colliders%radius(jtarg) frag%density(1) = frag%mass_dist(1) / volume(jtarg) - - frag%mass(2:nfrag) = (frag%mtot - frag%mass(1)) / (nfrag - 1) - frag%mass(nfrag) = frag%mass(nfrag) + (frag%mtot - sum(frag%mass(:))) - frag%radius(2:nfrag) = (3 * frag%mass(2:nfrag) / (4 * PI * frag%density(2:nfrag)))**(1.0_DP / 3.0_DP) - case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - frag%mass(1) = frag%mtot - frag%radius(1) = (3 * sum(volume(:)) / (4 * PI))**(1._DP / 3._DP) + frag%Ip(:, 1) = colliders%Ip(:,1) + istart = 2 case default - write(*,*) "fraggle_set_mass_dist_fragments error: Unrecognized regime code",frag%regime + istart = 1 end select + if (mremaining > 0.0_DP) then + ! Distribute remaining mass among the fragments + mfrag = 1._DP + mremaining / sum(frag%mass(istart:nfrag)) + do i = istart, nfrag + frag%mass(i) = frag%mass(i) * mfrag + end do + mremaining = frag%mtot - sum(frag%mass(1:nfrag)) + frag%mass(nfrag) = frag%mass(nfrag) + mremaining + end if + frag%density(istart:nfrag) = frag%mtot / sum(volume(:)) + frag%radius(istart:nfrag) = (3 * frag%mass(istart:nfrag) / (4 * PI * frag%density(istart:nfrag)))**(1.0_DP / 3.0_DP) + do i = istart, nfrag + frag%Ip(:, i) = Ip_avg(:) + end do end associate @@ -181,7 +227,6 @@ module subroutine fraggle_set_natural_scale_factors(self, colliders) associate(frag => self) ! Find the center of mass of the collisional system - frag%mtot = sum(colliders%mass(:)) frag%xbcom(:) = (colliders%mass(1) * colliders%xb(:,1) + colliders%mass(2) * colliders%xb(:,2)) / frag%mtot frag%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / frag%mtot diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index 8dc81c303..5a8a4a392 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -173,7 +173,7 @@ module subroutine fraggle_regime_colliders(self, frag, system, param) class(fraggle_colliders), intent(inout) :: self !! Fraggle colliders object class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragment system object class(swiftest_nbody_system), intent(in) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameter + class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters end subroutine fraggle_regime_colliders module subroutine fraggle_set_budgets_fragments(self, colliders) @@ -188,10 +188,11 @@ module subroutine fraggle_set_coordinate_system(self, colliders) class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object end subroutine fraggle_set_coordinate_system - module subroutine fraggle_set_mass_dist_fragments(self, colliders) + module subroutine fraggle_set_mass_dist_fragments(self, colliders, param) implicit none - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters end subroutine fraggle_set_mass_dist_fragments module subroutine fraggle_set_natural_scale_factors(self, colliders) diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 20e90909d..6cd2fd093 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -1,8 +1,6 @@ submodule (symba_classes) s_symba_collision use swiftest - integer(I4B), parameter :: NFRAG_DISRUPT = 12 - integer(I4B), parameter :: NFRAG_SUPERCAT = 20 contains module function symba_collision_casedisruption(system, param, colliders, frag) result(status) @@ -26,17 +24,14 @@ module function symba_collision_casedisruption(system, param, colliders, frag) select case(frag%regime) case(COLLRESOLVE_REGIME_DISRUPTION) message = "Disruption between" - nfrag = NFRAG_DISRUPT case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) message = "Supercatastrophic disruption between" - nfrag = NFRAG_SUPERCAT end select call symba_collision_collider_message(system%pl, colliders%idx, message) call fraggle_io_log_one_message(message) ! Collisional fragments will be uniformly distributed around the pre-impact barycenter - call frag%setup(nfrag, param) - call frag%set_mass_dist(colliders) + call frag%set_mass_dist(colliders, param) ! Generate the position and velocity distributions of the fragments call frag%generate_fragments(colliders, system, param, lfailure) @@ -53,6 +48,7 @@ module function symba_collision_casedisruption(system, param, colliders, frag) end select else ! Populate the list of new bodies + nfrag = frag%nbody write(message, *) nfrag call fraggle_io_log_one_message("Generating " // trim(adjustl(message)) // " fragments") select case(frag%regime) @@ -106,10 +102,8 @@ module function symba_collision_casehitandrun(system, param, colliders, frag) r nfrag = 0 lpure = .true. else ! Imperfect hit and run, so we'll keep the largest body and destroy the other - nfrag = NFRAG_DISRUPT - 1 lpure = .false. - call frag%setup(nfrag, param) - call frag%set_mass_dist(colliders) + call frag%set_mass_dist(colliders, param) ! Generate the position and velocity distributions of the fragments call frag%generate_fragments(colliders, system, param, lpure) @@ -118,6 +112,7 @@ module function symba_collision_casehitandrun(system, param, colliders, frag) r call fraggle_io_log_one_message("Should have been a pure hit and run instead") nfrag = 0 else + nfrag = frag%nbody write(message, *) nfrag call fraggle_io_log_one_message("Generating " // trim(adjustl(message)) // " fragments") end if @@ -174,8 +169,7 @@ module function symba_collision_casemerge(system, param, colliders, frag) resul select type(pl => system%pl) class is (symba_pl) - call frag%setup(1, param) - call frag%set_mass_dist(colliders) + call frag%set_mass_dist(colliders, param) ibiggest = colliders%idx(maxloc(pl%Gmass(colliders%idx(:)), dim=1)) frag%id(1) = pl%id(ibiggest) frag%xb(:,1) = frag%xbcom(:) @@ -723,7 +717,7 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) associate(info => pl%info, pl_adds => system%pl_adds, cb => system%cb, npl => pl%nbody) ! Add the colliders%idx bodies to the subtraction list ncolliders = colliders%ncoll - nfrag = size(frag%mass(:)) + nfrag = frag%nbody ! Setup new bodies allocate(plnew, mold=pl)