diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 55a8b23a1..2e63029ce 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -324,6 +324,66 @@ module subroutine collision_util_reset_system(self) end subroutine collision_util_reset_system + + + module subroutine collision_util_set_coordinate_system(self) + !! author: David A. Minton + !! + !! Defines the collisional coordinate system, including the unit vectors of both the system and individual fragments. + implicit none + ! Arguments + class(collision_system), intent(inout) :: self !! Collisional system + ! Internals + integer(I4B) :: i + real(DP), dimension(NDIM) :: delta_r, delta_v, Ltot + real(DP) :: L_mag + real(DP), dimension(NDIM, self%fragments%nbody) :: L_sigma + + associate(fragments => self%fragments, impactors => self%impactors, nfrag => self%fragments%nbody) + delta_v(:) = impactors%vb(:, 2) - impactors%vb(:, 1) + delta_r(:) = impactors%rb(:, 2) - impactors%rb(:, 1) + + ! We will initialize fragments on a plane defined by the pre-impact system, with the z-axis aligned with the angular momentum vector + ! and the y-axis aligned with the pre-impact distance vector. + + ! y-axis is the separation distance + fragments%y_coll_unit(:) = .unit.delta_r(:) + Ltot = impactors%Lorbit(:,1) + impactors%Lorbit(:,2) + impactors%Lspin(:,1) + impactors%Lspin(:,2) + + L_mag = .mag.Ltot(:) + if (L_mag > sqrt(tiny(L_mag))) then + fragments%z_coll_unit(:) = .unit.Ltot(:) + else ! Not enough angular momentum to determine a z-axis direction. We'll just pick a random direction + call random_number(fragments%z_coll_unit(:)) + fragments%z_coll_unit(:) = .unit.fragments%z_coll_unit(:) + end if + + ! The cross product of the y- by z-axis will give us the x-axis + fragments%x_coll_unit(:) = fragments%y_coll_unit(:) .cross. fragments%z_coll_unit(:) + + fragments%v_coll_unit(:) = .unit.delta_v(:) + + if (.not.any(fragments%r_coll(:,:) > 0.0_DP)) return + fragments%rmag(:) = .mag. fragments%r_coll(:,:) + + ! 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 system + call random_number(L_sigma(:,:)) + do concurrent(i = 1:nfrag, fragments%rmag(i) > 0.0_DP) + fragments%v_n_unit(:, i) = fragments%z_coll_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%r_coll(:,:) + 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 + + 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 7c9cfa358..e603b5b9d 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -45,7 +45,9 @@ module subroutine fraggle_generate_fragments(self, system, param, lfailure) fpe_quiet_modes(:) = .false. call ieee_set_halting_mode(IEEE_ALL,fpe_quiet_modes) - associate(fragments => self, nfrag => self%nbody, pl => system%pl) + select type(fragments => self%fragments) + class is (fraggle_fragments) + associate(collision => self, impactors => self%impactors, nfrag => fragments%nbody, pl => system%pl) write(message,*) nfrag call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle generating " // trim(adjustl(message)) // " fragments.") @@ -72,7 +74,7 @@ module subroutine fraggle_generate_fragments(self, system, param, lfailure) call fragments%reset() ! Calculate the initial energy of the system without the collisional family - call fragments%get_energy_and_momentum(impactors, system, param, lbefore=.true.) + call collision%get_energy_and_momentum(system, param, lbefore=.true.) ! Start out the fragments close to the initial separation distance. This will be increased if there is any overlap or we fail to find a solution r_max_start = 1.2_DP * .mag.(impactors%rb(:,2) - impactors%rb(:,1)) @@ -162,6 +164,7 @@ module subroutine fraggle_generate_fragments(self, system, param, lfailure) ! Restore the big array if (lk_plpl) call pl%flatten(param) end associate + end select call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily return diff --git a/src/fraggle/fraggle_set.f90 b/src/fraggle/fraggle_set.f90 index fa0e28aeb..a69a69260 100644 --- a/src/fraggle/fraggle_set.f90 +++ b/src/fraggle/fraggle_set.f90 @@ -160,64 +160,6 @@ module subroutine fraggle_set_mass_dist_fragments(self, impactors, param) end subroutine fraggle_set_mass_dist_fragments - module subroutine encounter_set_coordinate_system(self, impactors) - !! author: David A. Minton - !! - !! Defines the collisional coordinate system, including the unit vectors of both the system and individual fragments. - implicit none - ! Arguments - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object - ! Internals - integer(I4B) :: i - real(DP), dimension(NDIM) :: delta_r, delta_v, Ltot - real(DP) :: L_mag - real(DP), dimension(NDIM, self%nbody) :: L_sigma - - associate(fragments => self, nfrag => self%nbody) - delta_v(:) = impactors%vb(:, 2) - impactors%vb(:, 1) - delta_r(:) = impactors%rb(:, 2) - impactors%rb(:, 1) - - ! We will initialize fragments on a plane defined by the pre-impact system, with the z-axis aligned with the angular momentum vector - ! and the y-axis aligned with the pre-impact distance vector. - - ! y-axis is the separation distance - fragments%y_coll_unit(:) = .unit.delta_r(:) - Ltot = impactors%Lorbit(:,1) + impactors%Lorbit(:,2) + impactors%Lspin(:,1) + impactors%Lspin(:,2) - - L_mag = .mag.Ltot(:) - if (L_mag > sqrt(tiny(L_mag))) then - fragments%z_coll_unit(:) = .unit.Ltot(:) - else ! Not enough angular momentum to determine a z-axis direction. We'll just pick a random direction - call random_number(fragments%z_coll_unit(:)) - fragments%z_coll_unit(:) = .unit.fragments%z_coll_unit(:) - end if - - ! The cross product of the y- by z-axis will give us the x-axis - fragments%x_coll_unit(:) = fragments%y_coll_unit(:) .cross. fragments%z_coll_unit(:) - - fragments%v_coll_unit(:) = .unit.delta_v(:) - - if (.not.any(fragments%r_coll(:,:) > 0.0_DP)) return - fragments%rmag(:) = .mag. fragments%r_coll(:,:) - - ! 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 system - call random_number(L_sigma(:,:)) - do concurrent(i = 1:nfrag, fragments%rmag(i) > 0.0_DP) - fragments%v_n_unit(:, i) = fragments%z_coll_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%r_coll(:,:) - 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 encounter_set_coordinate_system - module subroutine fraggle_set_natural_scale_factors(self, impactors) !! author: David A. Minton diff --git a/src/modules/collision_classes.f90 b/src/modules/collision_classes.f90 index d2d2bd9b1..262920269 100644 --- a/src/modules/collision_classes.f90 +++ b/src/modules/collision_classes.f90 @@ -102,11 +102,12 @@ module collision_classes real(DP), dimension(2) :: pe !! Before/after potential energy real(DP), dimension(2) :: Etot !! Before/after total system energy contains - procedure :: regime => collision_regime_system !! Determine which fragmentation regime the set of impactors will be - procedure :: setup => collision_setup_system !! Initializer for the encounter collision system. Allocates the collider and fragments classes and the before/after snapshots - procedure :: get_energy_and_momentum => collision_util_get_energy_momentum !! Calculates total 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 - final :: collision_util_final_system !! Finalizer will deallocate all allocatables + procedure :: regime => collision_regime_system !! Determine which fragmentation regime the set of impactors will be + procedure :: setup => collision_setup_system !! Initializer for the encounter collision system. Allocates the collider and fragments classes and the before/after snapshots + procedure :: get_energy_and_momentum => collision_util_get_energy_momentum !! Calculates total 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 system of the collisional system + final :: collision_util_final_system !! Finalizer will deallocate all allocatables end type collision_system !! NetCDF dimension and variable names for the enounter save object @@ -217,6 +218,12 @@ module subroutine collision_set_coordinate_fragments(self) class(collision_fragments), intent(inout) :: self !! Fragment system object end subroutine collision_set_coordinate_fragments + module subroutine collision_util_set_coordinate_system(self) + implicit none + class(collision_system), intent(inout) :: self !! Collisional system + end subroutine collision_util_set_coordinate_system + + module subroutine collision_setup_fragments(self, n, param) implicit none class(collision_fragments), intent(inout) :: self !! Fragment system object