From d54cb62247dd0fd4109550db9d647ee5ecec15d6 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 23 Dec 2022 06:10:56 -0500 Subject: [PATCH] Restructured a bit so that the fragment SFD generator is in the main collision module, not Fraggle --- src/collision/collision_generate.f90 | 2 +- src/collision/collision_module.f90 | 19 ++-- src/collision/collision_setup.f90 | 2 +- src/collision/collision_util.f90 | 136 +++++++++++++++++++++++++++ src/fraggle/fraggle_module.f90 | 23 ++--- src/fraggle/fraggle_set.f90 | 136 --------------------------- src/fraggle/fraggle_setup.f90 | 1 + src/swiftest/swiftest_setup.f90 | 2 +- 8 files changed, 159 insertions(+), 162 deletions(-) diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 421cf0558..59b78fb26 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -219,7 +219,7 @@ end subroutine collision_generate_bounce module subroutine collision_generate_simple(self, nbody_system, param, t) implicit none - class(collision_simple), intent(inout) :: self !! Simple fragment system object + class(collision_simple_disruption), intent(inout) :: self !! Simple fragment system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! The time of the collision diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 7b8331d1b..98463fb93 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -147,11 +147,12 @@ module collision final :: collision_final_bounce !! Finalizer will deallocate all allocatables end type collision_bounce - type, extends(collision_merge) :: collision_simple + type, extends(collision_merge) :: collision_simple_disruption contains - procedure :: generate => collision_generate_simple !! If a collision would result in a disruption [TODO: SOMETHING LIKE CHAMBERS 2012] - final :: collision_final_simple !! Finalizer will deallocate all allocatables - end type collision_simple + procedure :: generate => collision_generate_simple !! If a collision would result in a disruption [TODO: SOMETHING LIKE CHAMBERS 2012] + procedure :: set_mass_dist => collision_util_set_mass_dist !! Sets the distribution of mass among the fragments depending on the regime type + final :: collision_final_simple !! Finalizer will deallocate all allocatables + end type collision_simple_disruption !! NetCDF dimension and variable names for the enounter save object @@ -224,7 +225,7 @@ end subroutine collision_generate_bounce module subroutine collision_generate_simple(self, nbody_system, param, t) implicit none - class(collision_simple), intent(inout) :: self !! Simple fragment nbody_system object + class(collision_simple_disruption), intent(inout) :: self !! Simple fragment nbody_system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! The time of the collision @@ -389,6 +390,12 @@ module subroutine collision_util_reset_fragments(self) class(collision_fragments(*)), intent(inout) :: self end subroutine collision_util_reset_fragments + module subroutine collision_util_set_mass_dist(self, param) + implicit none + class(collision_simple_disruption), intent(inout) :: self !! Simple disruption collision object + class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters + end subroutine collision_util_set_mass_dist + module subroutine collision_util_get_idvalues_snapshot(self, idvals) implicit none class(collision_snapshot), intent(in) :: self !! Collision snapshot object @@ -554,7 +561,7 @@ subroutine collision_final_simple(self) !! Finalizer will deallocate all allocatables implicit none ! Arguments - type(collision_simple), intent(inout) :: self !! Collision system object + type(collision_simple_disruption), intent(inout) :: self !! Collision system object call self%reset() if (allocated(self%impactors)) deallocate(self%impactors) diff --git a/src/collision/collision_setup.f90 b/src/collision/collision_setup.f90 index 5b621a44f..2c4880cb0 100644 --- a/src/collision/collision_setup.f90 +++ b/src/collision/collision_setup.f90 @@ -58,11 +58,11 @@ module subroutine collision_setup_fragments_system(self, nfrag) if (allocated(self%fragments)) deallocate(self%fragments) allocate(collision_fragments(nfrag) :: self%fragments) + self%fragments%nbody = nfrag return end subroutine collision_setup_fragments_system - end submodule s_collision_setup \ No newline at end of file diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 9338597a3..410090107 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -435,6 +435,142 @@ module subroutine collision_util_set_coordinate_system(self) end subroutine collision_util_set_coordinate_system + module subroutine collision_util_set_mass_dist(self, 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 routine has been run on the fragments + !! + implicit none + ! Arguments + class(collision_simple_disruption), intent(inout) :: self !! Fraggle collision system object + class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters + ! Internals + integer(I4B) :: i, jproj, jtarg, nfrag, istart + real(DP), dimension(2) :: volume + real(DP), dimension(NDIM) :: Ip_avg + real(DP) :: mfrag, mremaining, min_mfrag, mtot + 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 + integer(I4B), parameter :: iMlr = 1 + integer(I4B), parameter :: iMslr = 2 + integer(I4B), parameter :: iMrem = 3 + + associate(impactors => self%impactors) + ! Get mass weighted mean of Ip and density + volume(1:2) = 4._DP / 3._DP * PI * impactors%radius(1:2)**3 + mtot = sum(impactors%mass(:)) + Ip_avg(:) = (impactors%mass(1) * impactors%Ip(:,1) + impactors%mass(2) * impactors%Ip(:,2)) / mtot + + if (impactors%mass(1) > impactors%mass(2)) then + jtarg = 1 + jproj = 2 + else + jtarg = 2 + jproj = 1 + end if + + select case(impactors%regime) + case(COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC, COLLRESOLVE_REGIME_HIT_AND_RUN) + ! The first two bins of the mass_dist are the largest and second-largest fragments that came out of collision_regime. + ! The remainder from the third bin will be distributed among nfrag-2 bodies. The following code will determine nfrag based on + ! the limits bracketed above and the model size distribution of fragments. + ! Check to see if our size distribution would give us a smaller number of fragments than the maximum number + + select type(param) + class is (swiftest_parameters) + min_mfrag = (param%min_GMfrag / param%GU) + ! The number of fragments we generate is bracked by the minimum required by fraggle_generate (7) and the + ! maximum set by the NFRAG_SIZE_MULTIPLIER which limits the total number of fragments to prevent the nbody + ! code from getting an overwhelmingly large number of fragments + nfrag = ceiling(NFRAG_SIZE_MULTIPLIER * log(mtot / min_mfrag)) + nfrag = max(min(nfrag, NFRAGMAX), NFRAGMIN) + class default + min_mfrag = 0.0_DP + nfrag = NFRAGMAX + end select + + i = iMrem + mremaining = impactors%mass_dist(iMrem) + do while (i <= nfrag) + mfrag = (1 + i - iMslr)**(-3._DP / BETA) * impactors%mass_dist(iMslr) + if (mremaining - mfrag < 0.0_DP) exit + mremaining = mremaining - mfrag + i = i + 1 + end do + if (i < nfrag) nfrag = max(i, NFRAGMIN) ! The sfd would actually give us fewer fragments than our maximum + call self%setup_fragments(nfrag) + + case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) + + call self%setup_fragments(1) + select type(fragments => self%fragments) + class is (collision_fragments(*)) + fragments%mass(1) = impactors%mass_dist(1) + fragments%radius(1) = impactors%radius(jtarg) + fragments%density(1) = impactors%mass_dist(1) / volume(jtarg) + if (param%lrotation) fragments%Ip(:, 1) = impactors%Ip(:,1) + end select + return + case default + write(*,*) "collision_util_set_mass_dist_fragments error: Unrecognized regime code",impactors%regime + end select + + select type(fragments => self%fragments) + class is (collision_fragments(*)) + fragments%mtot = mtot + + ! Make the first two bins the same as the Mlr and Mslr values that came from collision_regime + fragments%mass(1) = impactors%mass_dist(iMlr) + fragments%mass(2) = impactors%mass_dist(iMslr) + + ! Distribute the remaining mass the 3:nfrag bodies following the model SFD given by slope BETA + mremaining = impactors%mass_dist(iMrem) + do i = iMrem, nfrag + mfrag = (1 + i - iMslr)**(-3._DP / BETA) * impactors%mass_dist(iMslr) + fragments%mass(i) = mfrag + mremaining = mremaining - mfrag + end do + + ! If there is any residual mass (either positive or negative) we will distribute remaining mass proportionally among the the fragments + if (mremaining < 0.0_DP) then ! If the remainder is negative, this means that that the number of fragments required by the SFD is smaller than our lower limit set by fraggle_generate. + istart = iMrem ! We will reduce the mass of the 3:nfrag bodies to prevent the second-largest fragment from going smaller + else ! If the remainder is postiive, this means that the number of fragments required by the SFD is larger than our upper limit set by computational expediency. + istart = iMslr ! We will increase the mass of the 2:nfrag bodies to compensate, which ensures that the second largest fragment remains the second largest + end if + mfrag = 1._DP + mremaining / sum(fragments%mass(istart:nfrag)) + fragments%mass(istart:nfrag) = fragments%mass(istart:nfrag) * mfrag + + ! There may still be some small residual due to round-off error. If so, simply add it to the last bin of the mass distribution. + mremaining = fragments%mtot - sum(fragments%mass(1:nfrag)) + fragments%mass(nfrag) = fragments%mass(nfrag) + mremaining + + ! Compute physical properties of the new fragments + select case(impactors%regime) + case(COLLRESOLVE_REGIME_HIT_AND_RUN) ! The hit and run case always preserves the largest body intact, so there is no need to recompute the physical properties of the first fragment + fragments%radius(1) = impactors%radius(jtarg) + fragments%density(1) = impactors%mass_dist(iMlr) / volume(jtarg) + fragments%Ip(:, 1) = impactors%Ip(:,1) + istart = 2 + case default + istart = 1 + end select + + fragments%density(istart:nfrag) = fragments%mtot / sum(volume(:)) + fragments%radius(istart:nfrag) = (3 * fragments%mass(istart:nfrag) / (4 * PI * fragments%density(istart:nfrag)))**(1.0_DP / 3.0_DP) + do i = istart, nfrag + fragments%Ip(:, i) = Ip_avg(:) + end do + + end select + end associate + + return + end subroutine collision_util_set_mass_dist + + subroutine collision_util_save_snapshot(collision_history, snapshot) !! author: David A. Minton !! diff --git a/src/fraggle/fraggle_module.f90 b/src/fraggle/fraggle_module.f90 index 83aaee06b..2471b9975 100644 --- a/src/fraggle/fraggle_module.f90 +++ b/src/fraggle/fraggle_module.f90 @@ -15,10 +15,8 @@ module fraggle implicit none public - !> Class definition for the variables that describe a collection of fragments by Fraggle barycentric coordinates type, extends(collision_fragments) :: fraggle_fragments - real(DP), dimension(nbody) :: v_r_mag !! Array of radial direction velocity magnitudes of individual fragments real(DP), dimension(nbody) :: v_t_mag !! Array of tangential direction velocity magnitudes of individual fragments real(DP), dimension(nbody) :: v_n_mag !! Array of normal direction velocity magnitudes of individual fragments @@ -28,16 +26,15 @@ module fraggle real(DP) :: ke_spin !! Spin kinetic energy of all fragments real(DP) :: ke_budget !! Kinetic energy budget for computing fragment trajectories real(DP), dimension(NDIM) :: L_budget !! Angular momentum budget for computing fragment trajectories - contains - procedure :: get_angular_momentum => fraggle_util_get_angular_momentum !! Calcualtes the current angular momentum of the fragments - procedure :: reset => fraggle_util_reset_fragments !! Resets all position and velocity-dependent fragment quantities in order to do a fresh calculation (does not reset mass, radius, or other values that get set prior to the call to fraggle_generate) - procedure :: restructure => fraggle_util_restructure !! Restructure the inputs after a failed attempt failed to find a set of positions and velocities that satisfy the energy and momentum constraints - final :: fraggle_final_fragments !! Finalizer will deallocate all allocatables + procedure :: get_angular_momentum => fraggle_util_get_angular_momentum !! Calcualtes the current angular momentum of the fragments + procedure :: reset => fraggle_util_reset_fragments !! Resets all position and velocity-dependent fragment quantities in order to do a fresh calculation (does not reset mass, radius, or other values that get set prior to the call to fraggle_generate) + procedure :: restructure => fraggle_util_restructure !! Restructure the inputs after a failed attempt failed to find a set of positions and velocities that satisfy the energy and momentum constraints + final :: fraggle_final_fragments !! Finalizer will deallocate all allocatables end type fraggle_fragments - type, extends(collision_merge) :: collision_fraggle + type, extends(collision_simple_disruption) :: collision_fraggle ! Scale factors used to scale dimensioned quantities to a more "natural" system where important quantities (like kinetic energy, momentum) are of order ~1 real(DP) :: dscale = 1.0_DP !! Distance dimension scale factor real(DP) :: mscale = 1.0_DP !! Mass scale factor @@ -48,13 +45,12 @@ module fraggle contains procedure :: generate => fraggle_generate_system !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. procedure :: set_budgets => fraggle_set_budgets !! Sets the energy and momentum budgets of the fragments based on the collider value - procedure :: set_mass_dist => fraggle_set_mass_dist !! Sets the distribution of mass among the fragments depending on the regime type procedure :: set_natural_scale => fraggle_set_natural_scale_factors !! Scales dimenional quantities to ~O(1) with respect to the collisional system. procedure :: set_original_scale => fraggle_set_original_scale_factors !! Restores dimenional quantities back to the original system units procedure :: setup_fragments => fraggle_setup_fragments_system !! Initializer for the fragments of the collision system. procedure :: construct_temporary_system => fraggle_util_construct_temporary_system !! Constructs temporary n-body system in order to compute pre- and post-impact energy and momentum procedure :: reset => fraggle_util_reset_system !! Deallocates all allocatables - final :: fraggle_final_system !! Finalizer will deallocate all allocatables + final :: fraggle_final_system !! Finalizer will deallocate all allocatables end type collision_fraggle @@ -98,12 +94,6 @@ module subroutine fraggle_set_budgets(self) class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object end subroutine fraggle_set_budgets - module subroutine fraggle_set_mass_dist(self, param) - implicit none - class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object - class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters - end subroutine fraggle_set_mass_dist - module subroutine fraggle_set_natural_scale_factors(self) implicit none class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object @@ -176,7 +166,6 @@ end function fraggle_util_vmag_to_vb end interface contains - subroutine fraggle_final_fragments(self) !! author: David A. Minton diff --git a/src/fraggle/fraggle_set.f90 b/src/fraggle/fraggle_set.f90 index 21ba46588..03f2e5e88 100644 --- a/src/fraggle/fraggle_set.f90 +++ b/src/fraggle/fraggle_set.f90 @@ -40,142 +40,6 @@ module subroutine fraggle_set_budgets(self) end subroutine fraggle_set_budgets - module subroutine fraggle_set_mass_dist(self, 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 routine has been run on the fragments - !! - implicit none - ! Arguments - class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object - class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters - ! Internals - integer(I4B) :: i, jproj, jtarg, nfrag, istart - real(DP), dimension(2) :: volume - real(DP), dimension(NDIM) :: Ip_avg - real(DP) :: mfrag, mremaining, min_mfrag, mtot - 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 - integer(I4B), parameter :: iMlr = 1 - integer(I4B), parameter :: iMslr = 2 - integer(I4B), parameter :: iMrem = 3 - - associate(impactors => self%impactors) - ! Get mass weighted mean of Ip and density - volume(1:2) = 4._DP / 3._DP * PI * impactors%radius(1:2)**3 - mtot = sum(impactors%mass(:)) - Ip_avg(:) = (impactors%mass(1) * impactors%Ip(:,1) + impactors%mass(2) * impactors%Ip(:,2)) / mtot - - if (impactors%mass(1) > impactors%mass(2)) then - jtarg = 1 - jproj = 2 - else - jtarg = 2 - jproj = 1 - end if - - select case(impactors%regime) - case(COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC, COLLRESOLVE_REGIME_HIT_AND_RUN) - ! The first two bins of the mass_dist are the largest and second-largest fragments that came out of collision_regime. - ! The remainder from the third bin will be distributed among nfrag-2 bodies. The following code will determine nfrag based on - ! the limits bracketed above and the model size distribution of fragments. - ! Check to see if our size distribution would give us a smaller number of fragments than the maximum number - - select type(param) - class is (swiftest_parameters) - min_mfrag = (param%min_GMfrag / param%GU) - ! The number of fragments we generate is bracked by the minimum required by fraggle_generate (7) and the - ! maximum set by the NFRAG_SIZE_MULTIPLIER which limits the total number of fragments to prevent the nbody - ! code from getting an overwhelmingly large number of fragments - nfrag = ceiling(NFRAG_SIZE_MULTIPLIER * log(mtot / min_mfrag)) - nfrag = max(min(nfrag, NFRAGMAX), NFRAGMIN) - class default - min_mfrag = 0.0_DP - nfrag = NFRAGMAX - end select - - i = iMrem - mremaining = impactors%mass_dist(iMrem) - do while (i <= nfrag) - mfrag = (1 + i - iMslr)**(-3._DP / BETA) * impactors%mass_dist(iMslr) - if (mremaining - mfrag < 0.0_DP) exit - mremaining = mremaining - mfrag - i = i + 1 - end do - if (i < nfrag) nfrag = max(i, NFRAGMIN) ! The sfd would actually give us fewer fragments than our maximum - call self%setup_fragments(nfrag) - - case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - - call self%setup_fragments(1) - select type(fragments => self%fragments) - class is (fraggle_fragments(*)) - fragments%mass(1) = impactors%mass_dist(1) - fragments%radius(1) = impactors%radius(jtarg) - fragments%density(1) = impactors%mass_dist(1) / volume(jtarg) - if (param%lrotation) fragments%Ip(:, 1) = impactors%Ip(:,1) - end select - return - case default - write(*,*) "fraggle_set_mass_dist_fragments error: Unrecognized regime code",impactors%regime - end select - - select type(fragments => self%fragments) - class is (fraggle_fragments(*)) - fragments%mtot = mtot - - ! Make the first two bins the same as the Mlr and Mslr values that came from collision_regime - fragments%mass(1) = impactors%mass_dist(iMlr) - fragments%mass(2) = impactors%mass_dist(iMslr) - - ! Distribute the remaining mass the 3:nfrag bodies following the model SFD given by slope BETA - mremaining = impactors%mass_dist(iMrem) - do i = iMrem, nfrag - mfrag = (1 + i - iMslr)**(-3._DP / BETA) * impactors%mass_dist(iMslr) - fragments%mass(i) = mfrag - mremaining = mremaining - mfrag - end do - - ! If there is any residual mass (either positive or negative) we will distribute remaining mass proportionally among the the fragments - if (mremaining < 0.0_DP) then ! If the remainder is negative, this means that that the number of fragments required by the SFD is smaller than our lower limit set by fraggle_generate. - istart = iMrem ! We will reduce the mass of the 3:nfrag bodies to prevent the second-largest fragment from going smaller - else ! If the remainder is postiive, this means that the number of fragments required by the SFD is larger than our upper limit set by computational expediency. - istart = iMslr ! We will increase the mass of the 2:nfrag bodies to compensate, which ensures that the second largest fragment remains the second largest - end if - mfrag = 1._DP + mremaining / sum(fragments%mass(istart:nfrag)) - fragments%mass(istart:nfrag) = fragments%mass(istart:nfrag) * mfrag - - ! There may still be some small residual due to round-off error. If so, simply add it to the last bin of the mass distribution. - mremaining = fragments%mtot - sum(fragments%mass(1:nfrag)) - fragments%mass(nfrag) = fragments%mass(nfrag) + mremaining - - ! Compute physical properties of the new fragments - select case(impactors%regime) - case(COLLRESOLVE_REGIME_HIT_AND_RUN) ! The hit and run case always preserves the largest body intact, so there is no need to recompute the physical properties of the first fragment - fragments%radius(1) = impactors%radius(jtarg) - fragments%density(1) = impactors%mass_dist(iMlr) / volume(jtarg) - fragments%Ip(:, 1) = impactors%Ip(:,1) - istart = 2 - case default - istart = 1 - end select - - fragments%density(istart:nfrag) = fragments%mtot / sum(volume(:)) - fragments%radius(istart:nfrag) = (3 * fragments%mass(istart:nfrag) / (4 * PI * fragments%density(istart:nfrag)))**(1.0_DP / 3.0_DP) - do i = istart, nfrag - fragments%Ip(:, i) = Ip_avg(:) - end do - - end select - end associate - - return - end subroutine fraggle_set_mass_dist - - module subroutine fraggle_set_natural_scale_factors(self) !! author: David A. Minton !! diff --git a/src/fraggle/fraggle_setup.f90 b/src/fraggle/fraggle_setup.f90 index 534924c74..e54ca46fd 100644 --- a/src/fraggle/fraggle_setup.f90 +++ b/src/fraggle/fraggle_setup.f90 @@ -24,6 +24,7 @@ module subroutine fraggle_setup_fragments_system(self, nfrag) if (allocated(self%fragments)) deallocate(self%fragments) allocate(fraggle_fragments(nbody=nfrag) :: self%fragments) self%fragments%nbody = nfrag + return end subroutine fraggle_setup_fragments_system diff --git a/src/swiftest/swiftest_setup.f90 b/src/swiftest/swiftest_setup.f90 index ae7f0dc40..c693c96ed 100644 --- a/src/swiftest/swiftest_setup.f90 +++ b/src/swiftest/swiftest_setup.f90 @@ -121,7 +121,7 @@ module subroutine swiftest_setup_construct_system(nbody_system, param) case("BOUNCE") allocate(collision_bounce :: nbody_system%collider) case("SIMPLE") - allocate(collision_simple :: nbody_system%collider) + allocate(collision_simple_disruption :: nbody_system%collider) case("FRAGGLE") allocate(collision_fraggle :: nbody_system%collider) end select