From 8bb85e39f8c1605c32ed35e094e25681902885ec Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 27 Dec 2022 21:52:34 -0500 Subject: [PATCH] Rearranged some of the methods in the collision system to keep consistency when computing spins --- src/collision/collision_generate.f90 | 19 ++--------- src/collision/collision_module.f90 | 34 +++++++++++++------ src/collision/collision_util.f90 | 51 +++++++++++++++++++++++++++- src/fraggle/fraggle_generate.f90 | 34 ++----------------- src/fraggle/fraggle_module.f90 | 8 ----- src/fraggle/fraggle_util.f90 | 28 --------------- 6 files changed, 78 insertions(+), 96 deletions(-) diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 165ee42a8..684970bef 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -382,6 +382,7 @@ module subroutine collision_generate_disrupt(self, nbody_system, param, t, lfail logical, optional, intent(out) :: lfailure ! Internals + call self%set_budgets() call collision_generate_simple_pos_vec(self) call self%set_coordinate_system() call collision_generate_simple_vel_vec(self) @@ -487,9 +488,7 @@ module subroutine collision_generate_simple_rot_vec(collider) class(collision_basic), intent(inout) :: collider !! Fraggle collision system object ! Internals real(DP), dimension(NDIM) :: Lresidual, Lbefore, Lafter, Lspin, rot - real(DP) :: ke_residual, rotmag_old, rotmag_new, vmag_old, vmag_new integer(I4B) :: i, loop - integer(I4B) :: MAXLOOP = 10 associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) @@ -505,19 +504,10 @@ module subroutine collision_generate_simple_rot_vec(collider) fragments%rot(:,1) = Lspin(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1)) - Lresidual(:) = sum(impactors%Lspin(:,:) + impactors%Lorbit(:,:), dim=2) - Lspin(:) - ! Randomize the rotational vector direction of the n>1 fragments and distribute the residual momentum amongst them call random_number(fragments%rot(:,2:nfrag)) - rot(:) = Lresidual(:) / sum(fragments%mass(2:nfrag) * fragments%radius(2:nfrag)**2 * fragments%Ip(3,2:nfrag)) - - fragments%rot(:,2:nfrag) = .unit. fragments%rot(:,2:nfrag) * .mag. rot(:) - - if (.mag.(Lresidual(:)) > tiny(1.0_DP)) then - do i = 2,nfrag - fragments%rot(:,i) = fragments%rot(:,i) + Lresidual(:) / ((nfrag - 1) * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) - end do - end if + fragments%rot(:,2:nfrag) = .unit. fragments%rot(:,2:nfrag) * .mag. fragments%rot(:,1) + call fragments%set_spins() end associate @@ -576,12 +566,9 @@ module subroutine collision_generate_simple_vel_vec(collider) mass_vscale(:) = sqrt(2*mass_vscale(:) / maxval(mass_vscale(:))) fragments%vc(:,1) = .mag.impactors%vc(:,1) * impactors%bounce_unit(:) do concurrent(i = 2:nfrag) - j = fragments%origin_body(i) vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rb(:,j) + impactors%rbcom(:)) - vmag = .mag.impactors%vc(:,j) * vscale(i) * mass_vscale(i) - if (lhitandrun) then fragments%vc(:,i) = vmag * 0.5_DP * impactors%bounce_unit(:) * vsign(i) + vrot(:) else diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index ed94b6be1..4e33889c9 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -114,10 +114,13 @@ module collision real(DP), dimension(NDIM) :: Lspin !! Spin angular momentum vector of all fragments real(DP) :: ke_orbit !! Orbital kinetic energy of all fragments 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 :: reset => collision_util_reset_fragments !! Deallocates all allocatable arrays and sets everything else to 0 procedure :: get_angular_momentum => collision_util_get_angular_momentum !! Calcualtes the current angular momentum of the fragments procedure :: get_kinetic_energy => collision_util_get_kinetic_energy !! Calcualtes the current kinetic energy of the fragments + procedure :: set_spins => collision_util_set_spins !! Calcualtes the spins of all fragments from the angular momentum budget and residual final :: collision_final_fragments !! Finalizer deallocates all allocatables end type collision_fragments @@ -149,6 +152,7 @@ module collision 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_budgets => collision_util_set_budgets !! Sets the energy and momentum budgets of the fragments based on the collider value procedure :: set_coordinate_system => collision_util_set_coordinate_collider !! Sets the coordinate system of the collisional system procedure :: generate => collision_generate_basic !! Merges the impactors to make a single final body procedure :: hitandrun => collision_generate_hitandrun !! Merges the impactors to make a single final body @@ -393,7 +397,6 @@ 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_add_fragments_to_collider(self, nbody_system, param) implicit none class(collision_basic), intent(in) :: self !! Collision system object @@ -410,15 +413,14 @@ module subroutine collision_util_construct_temporary_system(self, nbody_system, class(base_parameters), allocatable, intent(out) :: tmpparam !! Output temporary configuration run parameters end subroutine collision_util_construct_temporary_system - module subroutine collision_util_get_angular_momentum(self) implicit none - class(collision_fragments(*)), intent(inout) :: self !! Fraggle fragment system object + class(collision_fragments(*)), intent(inout) :: self !! Fragment system object end subroutine collision_util_get_angular_momentum module subroutine collision_util_get_kinetic_energy(self) implicit none - class(collision_fragments(*)), intent(inout) :: self !! Fraggle fragment system object + class(collision_fragments(*)), intent(inout) :: self !! Fragment system object end subroutine collision_util_get_kinetic_energy module subroutine collision_util_reset_fragments(self) @@ -426,11 +428,10 @@ 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) + module subroutine collision_util_set_budgets(self) 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 + class(collision_basic), intent(inout) :: self !! Collision system object + end subroutine collision_util_set_budgets module subroutine collision_util_set_coordinate_collider(self) implicit none @@ -442,6 +443,17 @@ module subroutine collision_util_set_coordinate_impactors(self) class(collision_impactors), intent(inout) :: self !! collisional system end subroutine collision_util_set_coordinate_impactors + 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_set_spins(self) + implicit none + class(collision_fragments(*)), intent(inout) :: self !! Collision fragment system object + end subroutine collision_util_set_spins + module subroutine collision_util_setup_collider(self, nbody_system) implicit none class(collision_basic), intent(inout) :: self !! Encounter collision system object @@ -554,7 +566,7 @@ subroutine collision_final_plpl(self) !! Finalizer will deallocate all allocatables implicit none ! Arguments - type(collision_list_plpl), intent(inout) :: self !! Fraggle encountar storage object + type(collision_list_plpl), intent(inout) :: self !! PL-PL collision list object call self%dealloc() @@ -567,7 +579,7 @@ subroutine collision_final_pltp(self) !! Finalizer will deallocate all allocatables implicit none ! Arguments - type(collision_list_pltp), intent(inout) :: self !! Fraggle encountar storage object + type(collision_list_pltp), intent(inout) :: self !! PL-TP collision list object call self%dealloc() @@ -581,7 +593,7 @@ subroutine collision_final_snapshot(self) !! Finalizer will deallocate all allocatables implicit none ! Arguments - type(collision_snapshot), intent(inout) :: self !! Fraggle encountar storage object + type(collision_snapshot), intent(inout) :: self !! Collsion snapshot object call encounter_final_snapshot(self%encounter_snapshot) diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 94dce4f6d..6deacbf3d 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -175,7 +175,6 @@ module subroutine collision_util_get_idvalues_snapshot(self, idvals) end subroutine collision_util_get_idvalues_snapshot - module subroutine collision_util_get_angular_momentum(self) !! Author: David A. Minton !! @@ -417,6 +416,31 @@ module subroutine collision_util_reset_system(self) end subroutine collision_util_reset_system + module subroutine collision_util_set_budgets(self) + !! author: David A. Minton + !! + !! Sets the energy and momentum budgets of the fragments based on the collider values and the before/after values of energy and momentum + implicit none + ! Arguments + class(collision_basic), intent(inout) :: self !! Fraggle collision system object + ! Internals + real(DP) :: dEtot + real(DP), dimension(NDIM) :: dL + + associate(impactors => self%impactors, fragments => self%fragments) + + dEtot = self%Etot(1) + dL(:) = self%Ltot(:,1) + + fragments%L_budget(:) = -dL(:) + fragments%ke_budget = -(dEtot - impactors%Qloss) + + end associate + + return + end subroutine collision_util_set_budgets + + module subroutine collision_util_set_coordinate_collider(self) !! author: David A. Minton !! @@ -675,6 +699,31 @@ module subroutine collision_util_set_mass_dist(self, param) end subroutine collision_util_set_mass_dist + module subroutine collision_util_set_spins(self) + !! Author: David A. Minton + !! + !! Distributes any residual angular momentum into the spins of the n>1 fragments + implicit none + ! Arguments + class(collision_fragments(*)), intent(inout) :: self !! Collision fragment system object + ! Internals + integer(I4B) :: i + real(DP), dimension(NDIM) :: Lresidual + + call self%get_angular_momentum() + Lresidual(:) = self%L_budget(:) - (self%Lorbit(:) + self%Lspin(:)) + + ! Distribute residual angular momentum amongst the fragments + if (.mag.(Lresidual(:)) > tiny(1.0_DP)) then + do i = 2,self%nbody + self%rot(:,i) = self%rot(:,i) + Lresidual(:) / ((self%nbody - 1) * self%mass(i) * self%radius(i)**2 * self%Ip(3,i)) + end do + end if + + return + end subroutine collision_util_set_spins + + module subroutine collision_util_setup_collider(self, nbody_system) !! author: David A. Minton !! diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 8a22b41e4..958c1febf 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -197,7 +197,7 @@ subroutine fraggle_generate_minimize(collider, lfailure) end do ! Set spins in order to conserve angular momentum - call fraggle_generate_set_spins(fragments) + call collision_util_set_spins(fragments) if (.not.lfailure) exit tol = tol * 2 ! Keep increasing the tolerance until we converge on a solution @@ -238,7 +238,7 @@ function E_objective_function(val_input) result(fval) call collision_util_shift_vector_to_origin(tmp_frag%mass, tmp_frag%vc) ! Set spins in order to conserve angular momentum - call fraggle_generate_set_spins(tmp_frag) + call collision_util_set_spins(tmp_frag) ! Get the current kinetic energy of the system call tmp_frag%get_kinetic_energy() @@ -261,34 +261,4 @@ end function E_objective_function end subroutine fraggle_generate_minimize - - subroutine fraggle_generate_set_spins(fragments) - !! Author: David A. Minton - !! - !! Distributes any residual angular momentum into the spins of the n>1 fragments - implicit none - ! Arguments - class(fraggle_fragments(*)), intent(inout) :: fragments !! Fraggle fragment system object - ! Internals - integer(I4B) :: i - real(DP), dimension(NDIM) :: Lresidual - - call fragments%get_angular_momentum() - Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit(:) + fragments%Lspin(:)) - - ! Distribute residual angular momentum amongst the fragments - if (.mag.(Lresidual(:)) > tiny(1.0_DP)) then - do i = 2,fragments%nbody - fragments%rot(:,i) = fragments%rot(:,i) + Lresidual(:) / ((fragments%nbody - 1) * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) - end do - end if - - return - end subroutine fraggle_generate_set_spins - - - - - - end submodule s_fraggle_generate diff --git a/src/fraggle/fraggle_module.f90 b/src/fraggle/fraggle_module.f90 index a40e196b4..d54afa791 100644 --- a/src/fraggle/fraggle_module.f90 +++ b/src/fraggle/fraggle_module.f90 @@ -20,8 +20,6 @@ module fraggle 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 - 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 :: 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) @@ -39,7 +37,6 @@ module fraggle real(DP) :: Lscale = 1.0_DP !! Angular momentum scale factor (a convenience unit that is derived from dscale, tscale, and mscale) contains procedure :: disrupt => fraggle_generate_disrupt !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. - procedure :: set_budgets => fraggle_util_set_budgets !! Sets the energy and momentum budgets of the fragments based on the collider value procedure :: set_natural_scale => fraggle_util_set_natural_scale_factors !! Scales dimenional quantities to ~O(1) with respect to the collisional system. procedure :: set_original_scale => fraggle_util_set_original_scale_factors !! Restores dimenional quantities back to the original system units procedure :: setup_fragments => fraggle_util_setup_fragments_system !! Initializer for the fragments of the collision system. @@ -84,11 +81,6 @@ module subroutine fraggle_util_reset_system(self) class(collision_fraggle), intent(inout) :: self !! Collision system object end subroutine fraggle_util_reset_system - module subroutine fraggle_util_set_budgets(self) - implicit none - class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object - end subroutine fraggle_util_set_budgets - module subroutine fraggle_util_set_natural_scale_factors(self) implicit none class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index f85eed056..94e245c0b 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -87,34 +87,6 @@ module subroutine fraggle_util_reset_system(self) end subroutine fraggle_util_reset_system - module subroutine fraggle_util_set_budgets(self) - !! author: David A. Minton - !! - !! Sets the energy and momentum budgets of the fragments based on the collider values and the before/after values of energy and momentum - implicit none - ! Arguments - class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object - ! Internals - real(DP) :: dEtot - real(DP), dimension(NDIM) :: dL - - associate(impactors => self%impactors) - select type(fragments => self%fragments) - class is (fraggle_fragments(*)) - - dEtot = self%Etot(1) - dL(:) = self%Ltot(:,1) - - fragments%L_budget(:) = -dL(:) - fragments%ke_budget = -(dEtot - impactors%Qloss) - - end select - end associate - - return - end subroutine fraggle_util_set_budgets - - module subroutine fraggle_util_set_natural_scale_factors(self) !! author: David A. Minton !!