From bed6564e0bef980b714a963fdf631f94781ea6c3 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 23 Dec 2022 07:14:18 -0500 Subject: [PATCH] Cleanup and fixed up some of the coordinate system vector calcs --- src/collision/collision_generate.f90 | 120 +++++++++++++-------------- src/collision/collision_module.f90 | 21 ++--- src/collision/collision_util.f90 | 11 ++- src/fraggle/fraggle_generate.f90 | 6 +- 4 files changed, 81 insertions(+), 77 deletions(-) diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 59b78fb26..2d4a006a8 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -12,6 +12,64 @@ use swiftest contains + module subroutine collision_generate_bounce(self, nbody_system, param, t) + !! author: David A. Minton + !! + !! In this collision model, if the collision would result in a disruption, the bodies are instead "bounced" off + !! of the center of mass. This is done as a reflection in the 2-body equivalent distance vector direction. + implicit none + ! Arguments + class(collision_bounce), intent(inout) :: self !! Bounce 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 + ! Internals + integer(I4B) :: i,j,nfrag + real(DP), dimension(NDIM) :: vcom, rnorm + + select type(nbody_system) + class is (swiftest_nbody_system) + select type (pl => nbody_system%pl) + class is (swiftest_pl) + associate(impactors => nbody_system%collider%impactors, fragments => nbody_system%collider%fragments) + select case (impactors%regime) + case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + nfrag = size(impactors%id(:)) + do i = 1, nfrag + j = impactors%id(i) + vcom(:) = pl%vb(:,j) - impactors%vbcom(:) + rnorm(:) = .unit. (impactors%rb(:,2) - impactors%rb(:,1)) + ! Do the reflection + vcom(:) = vcom(:) - 2 * dot_product(vcom(:),rnorm(:)) * rnorm(:) + pl%vb(:,j) = impactors%vbcom(:) + vcom(:) + self%status = DISRUPTED + pl%status(j) = ACTIVE + pl%ldiscard(j) = .false. + pl%lcollision(j) = .false. + end do + select type(before => self%before) + class is (swiftest_nbody_system) + select type(after => self%after) + class is (swiftest_nbody_system) + allocate(after%pl, source=before%pl) ! Be sure to save the pl so that snapshots still work + end select + end select + case (COLLRESOLVE_REGIME_HIT_AND_RUN) + call collision_generate_hitandrun(self, nbody_system, param, t) + case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) + call self%collision_merge%generate(nbody_system, param, t) ! Use the default collision model, which is merge + case default + write(*,*) "Error in swiftest_collision, unrecognized collision regime" + call util_exit(FAILURE) + end select + end associate + end select + end select + + return + end subroutine collision_generate_bounce + + module subroutine collision_generate_hitandrun(collider, nbody_system, param, t) !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! @@ -159,70 +217,12 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) end subroutine collision_generate_merge - module subroutine collision_generate_bounce(self, nbody_system, param, t) - !! author: David A. Minton - !! - !! In this collision model, if the collision would result in a disruption, the bodies are instead "bounced" off - !! of the center of mass. This is done as a reflection in the 2-body equivalent distance vector direction. - implicit none - ! Arguments - class(collision_bounce), intent(inout) :: self !! Bounce 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 - ! Internals - integer(I4B) :: i,j,nfrag - real(DP), dimension(NDIM) :: vcom, rnorm - - select type(nbody_system) - class is (swiftest_nbody_system) - select type (pl => nbody_system%pl) - class is (swiftest_pl) - associate(impactors => nbody_system%collider%impactors, fragments => nbody_system%collider%fragments) - select case (impactors%regime) - case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - nfrag = size(impactors%id(:)) - do i = 1, nfrag - j = impactors%id(i) - vcom(:) = pl%vb(:,j) - impactors%vbcom(:) - rnorm(:) = .unit. (impactors%rb(:,2) - impactors%rb(:,1)) - ! Do the reflection - vcom(:) = vcom(:) - 2 * dot_product(vcom(:),rnorm(:)) * rnorm(:) - pl%vb(:,j) = impactors%vbcom(:) + vcom(:) - self%status = DISRUPTED - pl%status(j) = ACTIVE - pl%ldiscard(j) = .false. - pl%lcollision(j) = .false. - end do - select type(before => self%before) - class is (swiftest_nbody_system) - select type(after => self%after) - class is (swiftest_nbody_system) - allocate(after%pl, source=before%pl) ! Be sure to save the pl so that snapshots still work - end select - end select - case (COLLRESOLVE_REGIME_HIT_AND_RUN) - call collision_generate_hitandrun(self, nbody_system, param, t) - case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - call self%collision_merge%generate(nbody_system, param, t) ! Use the default collision model, which is merge - case default - write(*,*) "Error in swiftest_collision, unrecognized collision regime" - call util_exit(FAILURE) - end select - end associate - end select - end select - - return - end subroutine collision_generate_bounce - - - module subroutine collision_generate_simple(self, nbody_system, param, t) + module subroutine collision_generate_simple_disruption(self, nbody_system, param, t) implicit none 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 - end subroutine collision_generate_simple + end subroutine collision_generate_simple_disruption end submodule s_collision_generate \ No newline at end of file diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 40ef8b012..ca0002f24 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -66,13 +66,14 @@ module collision 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 !! z-direction unit vector of collisional nbody_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) :: 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 + 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 @@ -149,7 +150,7 @@ module collision 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] + procedure :: generate => collision_generate_simple_disruption !! 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 @@ -223,13 +224,13 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t) real(DP), intent(in) :: t !! The time of the collision end subroutine collision_generate_bounce - module subroutine collision_generate_simple(self, nbody_system, param, t) + module subroutine collision_generate_simple_disruption(self, nbody_system, param, t) implicit none 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 - end subroutine collision_generate_simple + end subroutine collision_generate_simple_disruption module subroutine collision_io_collider_message(pl, collidx, collider_message) implicit none diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 62e275d24..784802dba 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -87,7 +87,7 @@ module subroutine collision_util_construct_temporary_system(self, nbody_system, ! Set up a new system based on the original if (allocated(tmpparam)) deallocate(tmpparam) if (allocated(tmpsys)) deallocate(tmpsys) - allocate(tmpparam, source=param) + allocate(tmpparam_local, source=param) call swiftest_util_setup_construct_system(tmpsys_local, tmpparam_local) ! No test particles necessary for energy/momentum calcs @@ -411,10 +411,13 @@ module subroutine collision_util_set_coordinate_system(self) ! The cross product of the y- by z-axis will give us the x-axis impactors%x_unit(:) = impactors%y_unit(:) .cross. impactors%z_unit(:) - impactors%v_unit(:) = .unit.delta_v(:) - if (.not.any(fragments%rc(:,:) > 0.0_DP)) return + ! 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. @@ -623,7 +626,7 @@ module subroutine collision_util_setup_fragments_system(self, nfrag) return end subroutine collision_util_setup_fragments_system - + subroutine collision_util_save_snapshot(collision_history, snapshot) !! author: David A. Minton !! diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 73a2c0255..e555b11f2 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -252,9 +252,6 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai return end if - ! This is a factor that will "distort" the shape of the frgment cloud in the direction of the impact velocity - f_spin= .mag. (runit(:) .cross. vunit(:)) - if (param%lflatten_interactions) then lk_plpl = allocated(pl%k_plpl) if (lk_plpl) deallocate(pl%k_plpl) @@ -288,6 +285,9 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai call fraggle_generate_pos_vec(collider, r_max_start) call collider%set_coordinate_system() + ! This is a factor that will "distort" the shape of the fragment cloud in the direction of the impact velocity + f_spin= .mag. (impactors%y_unit(:) .cross. impactors%v_unit(:)) + ! Initial velocity guess will be the barycentric velocity of the colliding nbody_system so that the budgets are based on the much smaller collisional-frame velocities do concurrent (i = 1:nfrag) fragments%vb(:, i) = impactors%vbcom(:)