From 9627d1f8f6d96abeeaf3b92571c607b7dc0ddcb6 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 23 Dec 2022 11:00:07 -0500 Subject: [PATCH] Refactoring and cleaning, aiming to get collisional coordinate system computed earlier --- src/collision/collision_module.f90 | 60 ++++++++++++----------- src/collision/collision_resolve.f90 | 3 ++ src/collision/collision_util.f90 | 75 ++++++++++++++++++----------- src/fraggle/fraggle_generate.f90 | 4 +- 4 files changed, 84 insertions(+), 58 deletions(-) diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index ca0002f24..570e1cccc 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -65,21 +65,22 @@ module collision real(DP), dimension(:), allocatable :: mass_dist !! Distribution of fragment mass determined by the regime calculation (largest fragment, second largest, and remainder) real(DP) :: Mcb !! Mass of central body (used to compute potential energy in regime determination) - ! Values in a coordinate frame centered on the collider barycenter and collisional nbody_system unit vectors - real(DP), dimension(NDIM) :: x_unit !! x-direction unit vector of collisional nbody_system - real(DP), dimension(NDIM) :: y_unit !! y-direction unit vector of collisional nbody_system - real(DP), dimension(NDIM) :: z_unit !! z-direction unit vector of collisional nbody_system - real(DP), dimension(NDIM) :: v_unit !! velocity direction unit vector of collisional nbody_system + ! Values in a coordinate frame centered on the collider barycenter and collisional system unit vectors + real(DP), dimension(NDIM) :: x_unit !! x-direction unit vector of collisional system + real(DP), dimension(NDIM) :: y_unit !! y-direction unit vector of collisional system + real(DP), dimension(NDIM) :: z_unit !! z-direction unit vector of collisional system + real(DP), dimension(NDIM) :: v_unit !! velocity direction unit vector of collisional system real(DP), dimension(NDIM) :: rbcom !! Center of mass position vector of the collider nbody_system in nbody_system barycentric coordinates real(DP), dimension(NDIM) :: vbcom !! Velocity vector of the center of mass of the collider nbody_system in nbody_system barycentric coordinates real(DP), dimension(NDIM) :: rbimp !! Impact point position vector of the collider nbody_system in nbody_system barycentric coordinates real(DP), dimension(NDIM) :: vbimp !! The impact point velocity vector is the component of the velocity in the distance vector direction contains - procedure :: consolidate => collision_resolve_consolidate_impactors !! Consolidates a multi-body collision into an equivalent 2-body collision - procedure :: get_regime => collision_regime_impactors !! Determine which fragmentation regime the set of impactors will be - procedure :: reset => collision_util_reset_impactors !! Resets the collider object variables to 0 and deallocates the index and mass distributions - final :: collision_final_impactors !! Finalizer will deallocate all allocatables + procedure :: consolidate => collision_resolve_consolidate_impactors !! Consolidates a multi-body collision into an equivalent 2-body collision + procedure :: get_regime => collision_regime_impactors !! Determine which fragmentation regime the set of impactors will be + procedure :: reset => collision_util_reset_impactors !! Resets the collider object variables to 0 and deallocates the index and mass distributions + procedure :: set_coordinate_system => collision_util_set_coordinate_impactors !! Sets the coordinate system of the impactors + final :: collision_final_impactors !! Finalizer will deallocate all allocatables end type collision_impactors @@ -112,7 +113,7 @@ module collision type :: collision_merge - !! This class defines a collisional nbody_system that stores impactors and fragments. This is written so that various collision models (i.e. Fraggle) could potentially be used + !! This class defines a collisional system that stores impactors and fragments. This is written so that various collision models (i.e. Fraggle) could potentially be used !! to resolve collision by defining extended types of encounters_impactors and/or encounetr_fragments !! !! The generate method for this class is the merge model. This allows any extended type to have access to the merge procedure by selecting the collision_merge parent class @@ -131,15 +132,15 @@ module collision real(DP), dimension(2) :: pe !! Before/after potential energy real(DP), dimension(2) :: Etot !! Before/after total nbody_system energy contains - procedure :: setup => collision_util_setup_system !! Initializer for the encounter collision system and the before/after snapshots - procedure :: setup_impactors => collision_util_setup_impactors_system !! Initializer for the impactors for the encounter collision system. Deallocates old impactors before creating new ones - procedure :: setup_fragments => collision_util_setup_fragments_system !! Initializer for the fragments of the collision system. - procedure :: add_fragments => collision_util_add_fragments_to_system !! Add fragments to nbody_system + procedure :: setup => collision_util_setup_collider !! Initializer for the encounter collision system and the before/after snapshots + procedure :: setup_impactors => collision_util_setup_impactors_collider !! Initializer for the impactors for the encounter collision system. Deallocates old impactors before creating new ones + procedure :: setup_fragments => collision_util_setup_fragments_collider !! Initializer for the fragments of the collision system. + procedure :: add_fragments => collision_util_add_fragments_to_collider !! Add fragments to nbody_system procedure :: construct_temporary_system => collision_util_construct_temporary_system !! Constructs temporary n-body nbody_system in order to compute pre- and post-impact energy and momentum procedure :: get_energy_and_momentum => collision_util_get_energy_momentum !! Calculates total nbody_system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) procedure :: reset => collision_util_reset_system !! Deallocates all allocatables - procedure :: set_coordinate_system => collision_util_set_coordinate_system !! Sets the coordinate nbody_system of the collisional nbody_system - procedure :: generate => collision_generate_merge !! Merges the impactors to make a single final body + procedure :: set_coordinate_system => collision_util_set_coordinate_collider !! Sets the coordinate system of the collisional system + procedure :: generate => collision_generate_merge !! Merges the impactors to make a single final body end type collision_merge type, extends(collision_merge) :: collision_bounce @@ -348,34 +349,39 @@ module subroutine collision_resolve_pltp(self, nbody_system, param, t, dt, irec) integer(I4B), intent(in) :: irec !! Current recursion level end subroutine collision_resolve_pltp - module subroutine collision_util_set_coordinate_system(self) + module subroutine collision_util_set_coordinate_collider(self) implicit none - class(collision_merge), intent(inout) :: self !! Collisional nbody_system - end subroutine collision_util_set_coordinate_system + class(collision_merge), intent(inout) :: self !! collisional system + end subroutine collision_util_set_coordinate_collider - module subroutine collision_util_setup_system(self, nbody_system) + module subroutine collision_util_set_coordinate_impactors(self) + implicit none + class(collision_impactors), intent(inout) :: self !! collisional system + end subroutine collision_util_set_coordinate_impactors + + module subroutine collision_util_setup_collider(self, nbody_system) implicit none class(collision_merge), intent(inout) :: self !! Encounter collision system object class(base_nbody_system), intent(in) :: nbody_system !! Current nbody system. Used as a mold for the before/after snapshots - end subroutine collision_util_setup_system + end subroutine collision_util_setup_collider - module subroutine collision_util_setup_impactors_system(self) + module subroutine collision_util_setup_impactors_collider(self) implicit none class(collision_merge), intent(inout) :: self !! Encounter collision system object - end subroutine collision_util_setup_impactors_system + end subroutine collision_util_setup_impactors_collider - module subroutine collision_util_setup_fragments_system(self, nfrag) + module subroutine collision_util_setup_fragments_collider(self, nfrag) implicit none class(collision_merge), intent(inout) :: self !! Encounter collision system object integer(I4B), intent(in) :: nfrag !! Number of fragments to create - end subroutine collision_util_setup_fragments_system + end subroutine collision_util_setup_fragments_collider - module subroutine collision_util_add_fragments_to_system(self, nbody_system, param) + module subroutine collision_util_add_fragments_to_collider(self, nbody_system, param) implicit none class(collision_merge), intent(in) :: self !! Collision system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters - end subroutine collision_util_add_fragments_to_system + end subroutine collision_util_add_fragments_to_collider module subroutine collision_util_construct_temporary_system(self, nbody_system, param, tmpsys, tmpparam) implicit none diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index 4dbfc4ea1..dae44b9cc 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -148,6 +148,9 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa ! Destroy the kinship relationships for all members of this impactors%id call pl%reset_kinship(impactors%id(:)) + ! Set the coordinate system of the impactors + call impactors%set_coordinate_system() + end associate end select return diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index fb0cbee06..06a39f4ed 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -11,7 +11,7 @@ use swiftest contains - module subroutine collision_util_add_fragments_to_system(self, nbody_system, param) + module subroutine collision_util_add_fragments_to_collider(self, nbody_system, param) !! Author: David A. Minton !! !! Adds fragments to the temporary system pl object @@ -58,7 +58,7 @@ module subroutine collision_util_add_fragments_to_system(self, nbody_system, par end select return - end subroutine collision_util_add_fragments_to_system + end subroutine collision_util_add_fragments_to_collider module subroutine collision_util_construct_temporary_system(self, nbody_system, param, tmpsys, tmpparam) @@ -381,7 +381,7 @@ module subroutine collision_util_reset_system(self) end subroutine collision_util_reset_system - module subroutine collision_util_set_coordinate_system(self) + module subroutine collision_util_set_coordinate_collider(self) !! author: David A. Minton !! !! Defines the collisional coordinate nbody_system, including the unit vectors of both the nbody_system and individual fragments. @@ -390,11 +390,45 @@ module subroutine collision_util_set_coordinate_system(self) class(collision_merge), intent(inout) :: self !! Collisional nbody_system ! Internals integer(I4B) :: i - real(DP), dimension(NDIM) :: delta_r, delta_v, Ltot - real(DP) :: L_mag, mtot real(DP), dimension(NDIM, self%fragments%nbody) :: L_sigma associate(fragments => self%fragments, impactors => self%impactors, nfrag => self%fragments%nbody) + call impactors%set_coordinate_system() + + if ((.not.allocated(self%fragments)) .or. (nfrag == 0) .or. (.not.any(fragments%rc(:,:) > 0.0_DP))) return + + fragments%rmag(:) = .mag. fragments%rc(:,:) + + ! Randomize the tangential velocity direction. + ! This helps to ensure that the tangential velocity doesn't completely line up with the angular momentum vector, otherwise we can get an ill-conditioned nbody_system + call random_number(L_sigma(:,:)) + do concurrent(i = 1:nfrag, fragments%rmag(i) > 0.0_DP) + fragments%v_n_unit(:, i) = impactors%z_unit(:) + 2e-1_DP * (L_sigma(:,i) - 0.5_DP) + end do + + ! Define the radial, normal, and tangential unit vectors for each individual fragment + fragments%v_r_unit(:,:) = .unit. fragments%rc(:,:) + fragments%v_n_unit(:,:) = .unit. fragments%v_n_unit(:,:) + fragments%v_t_unit(:,:) = .unit. (fragments%v_n_unit(:,:) .cross. fragments%v_r_unit(:,:)) + + end associate + + return + end subroutine collision_util_set_coordinate_collider + + + module subroutine collision_util_set_coordinate_impactors(self) + !! author: David A. Minton + !! + !! Defines the collisional coordinate nbody_system, including the unit vectors of both the nbody_system and individual fragments. + implicit none + ! Arguments + class(collision_impactors), intent(inout) :: self !! Collisional nbody_system + ! Internals + real(DP), dimension(NDIM) :: delta_r, delta_v, Ltot + real(DP) :: L_mag, mtot + + associate(impactors => self) delta_v(:) = impactors%vb(:, 2) - impactors%vb(:, 1) delta_r(:) = impactors%rb(:, 2) - impactors%rb(:, 1) @@ -427,27 +461,10 @@ module subroutine collision_util_set_coordinate_system(self) ! The "bounce" unit vector is the projection of impactors%vbimp(:) = dot_product(delta_v(:),impactors%x_unit(:)) * impactors%x_unit(:) - - if ((.not.allocated(self%fragments)) .or. (nfrag == 0) .or. (.not.any(fragments%rc(:,:) > 0.0_DP))) return - - fragments%rmag(:) = .mag. fragments%rc(:,:) - - ! Randomize the tangential velocity direction. - ! This helps to ensure that the tangential velocity doesn't completely line up with the angular momentum vector, otherwise we can get an ill-conditioned nbody_system - call random_number(L_sigma(:,:)) - do concurrent(i = 1:nfrag, fragments%rmag(i) > 0.0_DP) - fragments%v_n_unit(:, i) = impactors%z_unit(:) + 2e-1_DP * (L_sigma(:,i) - 0.5_DP) - end do - - ! Define the radial, normal, and tangential unit vectors for each individual fragment - fragments%v_r_unit(:,:) = .unit. fragments%rc(:,:) - fragments%v_n_unit(:,:) = .unit. fragments%v_n_unit(:,:) - fragments%v_t_unit(:,:) = .unit. (fragments%v_n_unit(:,:) .cross. fragments%v_r_unit(:,:)) - end associate return - end subroutine collision_util_set_coordinate_system + end subroutine collision_util_set_coordinate_impactors module subroutine collision_util_set_mass_dist(self, param) @@ -586,7 +603,7 @@ module subroutine collision_util_set_mass_dist(self, param) end subroutine collision_util_set_mass_dist - module subroutine collision_util_setup_system(self, nbody_system) + module subroutine collision_util_setup_collider(self, nbody_system) !! author: David A. Minton !! !! Initializer for the encounter collision system. Sets up impactors and the before/after snapshots, @@ -604,10 +621,10 @@ module subroutine collision_util_setup_system(self, nbody_system) allocate(self%after, mold=nbody_system) return - end subroutine collision_util_setup_system + end subroutine collision_util_setup_collider - module subroutine collision_util_setup_impactors_system(self) + module subroutine collision_util_setup_impactors_collider(self) !! author: David A. Minton !! !! Initializer for the impactors for the encounter collision system. Deallocates old impactors before creating new ones @@ -619,10 +636,10 @@ module subroutine collision_util_setup_impactors_system(self) allocate(collision_impactors :: self%impactors) return - end subroutine collision_util_setup_impactors_system + end subroutine collision_util_setup_impactors_collider - module subroutine collision_util_setup_fragments_system(self, nfrag) + module subroutine collision_util_setup_fragments_collider(self, nfrag) !! author: David A. Minton !! !! Initializer for the fragments of the collision system. @@ -636,7 +653,7 @@ module subroutine collision_util_setup_fragments_system(self, nfrag) self%fragments%nbody = nfrag return - end subroutine collision_util_setup_fragments_system + end subroutine collision_util_setup_fragments_collider subroutine collision_util_save_snapshot(collision_history, snapshot) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 3cfb3ee82..b809e45f4 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -632,7 +632,7 @@ subroutine fraggle_generate_tan_vel(collider, lfailure) end do fragments%v_t_mag(1:nfrag) = solve_fragment_tan_vel(v_t_mag_input=v_t_initial(7:nfrag), lfailure=lfailure) - ! Perform one final shift of the radial velocity vectors to align with the center of mass of the collisional nbody_system (the origin) + ! Perform one final shift of the radial velocity vectors to align with the center of mass of the collisional system (the origin) fragments%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%v_r_unit(:,1:nfrag), fragments%v_t_mag(1:nfrag), & fragments%v_t_unit(:,1:nfrag), fragments%mass(1:nfrag), impactors%vbcom(:)) do concurrent (i = 1:nfrag) @@ -813,7 +813,7 @@ subroutine fraggle_generate_rad_vel(collider, lfailure) v_r_initial(:) = v_r_initial(:) * vnoise(:) end do - ! Shift the radial velocity vectors to align with the center of mass of the collisional nbody_system (the origin) + ! Shift the radial velocity vectors to align with the center of mass of the collisional system (the origin) fragments%ke_orbit = 0.0_DP fragments%ke_spin = 0.0_DP fragments%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%v_r_unit(:,1:nfrag), &