From aeef4fe87c7d220c38ca1272d5f9febb75d3f1e0 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 29 Dec 2022 11:42:42 -0500 Subject: [PATCH] Refactored name of disruption model --- src/collision/collision_generate.f90 | 22 ++++++++--------- src/collision/collision_module.f90 | 18 +++++++------- src/collision/collision_util.f90 | 7 ++++-- src/fraggle/fraggle_generate.f90 | 37 ++++++++++++++++------------ 4 files changed, 46 insertions(+), 38 deletions(-) diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index c3348841d..e91a19b7a 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -191,7 +191,7 @@ module subroutine collision_generate_hitandrun(self, nbody_system, param, t) end subroutine collision_generate_hitandrun - module subroutine collision_generate_simple(self, nbody_system, param, t) + module subroutine collision_generate_disruption(self, nbody_system, param, t) !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Create the fragments resulting from a non-catastrophic disruption collision @@ -263,7 +263,7 @@ module subroutine collision_generate_simple(self, nbody_system, param, t) end select end select return - end subroutine collision_generate_simple + end subroutine collision_generate_disruption module subroutine collision_generate_merge(self, nbody_system, param, t) @@ -384,15 +384,15 @@ module subroutine collision_generate_disrupt(self, nbody_system, param, t, lfail logical, optional, intent(out) :: lfailure ! Internals - call collision_generate_simple_pos_vec(self) - call collision_generate_simple_rot_vec(self) - call collision_generate_simple_vel_vec(self) + call collision_generate_disruption_pos_vec(self) + call collision_generate_disruption_rot_vec(self) + call collision_generate_disruption_vel_vec(self) return end subroutine collision_generate_disrupt - module subroutine collision_generate_simple_pos_vec(collider) + module subroutine collision_generate_disruption_pos_vec(collider) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Initializes the position vectors of the fragments around the center of mass based on the collision style. @@ -477,10 +477,10 @@ module subroutine collision_generate_simple_pos_vec(collider) end associate return - end subroutine collision_generate_simple_pos_vec + end subroutine collision_generate_disruption_pos_vec - module subroutine collision_generate_simple_rot_vec(collider) + module subroutine collision_generate_disruption_rot_vec(collider) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Calculates the spins of a collection of fragments such that they conserve angular momentum @@ -519,10 +519,10 @@ module subroutine collision_generate_simple_rot_vec(collider) end associate return - end subroutine collision_generate_simple_rot_vec + end subroutine collision_generate_disruption_rot_vec - module subroutine collision_generate_simple_vel_vec(collider) + module subroutine collision_generate_disruption_vel_vec(collider) !! Author: David A. Minton !! !! Generates an initial velocity distribution. For disruptions, the velocity magnitude is set to be @@ -644,7 +644,7 @@ module subroutine collision_generate_simple_vel_vec(collider) end associate return - end subroutine collision_generate_simple_vel_vec + end subroutine collision_generate_disruption_vel_vec 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 88698e360..217e9a74e 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -167,7 +167,7 @@ module collision type, extends(collision_basic) :: collision_disruption contains - procedure :: generate => collision_generate_simple !! A simple disruption models that does not constrain energy loss in collisions + procedure :: generate => collision_generate_disruption !! A simple disruption models that does not constrain energy loss in collisions procedure :: disrupt => collision_generate_disrupt !! Disrupt the colliders into the fragments 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 @@ -258,28 +258,28 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) real(DP), intent(in) :: t !! The time of the collision end subroutine collision_generate_merge - module subroutine collision_generate_simple(self, nbody_system, param, t) + module subroutine collision_generate_disruption(self, nbody_system, param, t) implicit none class(collision_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_disruption - module subroutine collision_generate_simple_pos_vec(collider) + module subroutine collision_generate_disruption_pos_vec(collider) implicit none class(collision_disruption), intent(inout) :: collider !! Collision system object - end subroutine collision_generate_simple_pos_vec + end subroutine collision_generate_disruption_pos_vec - module subroutine collision_generate_simple_rot_vec(collider) + module subroutine collision_generate_disruption_rot_vec(collider) implicit none class(collision_basic), intent(inout) :: collider !! Collision system object - end subroutine collision_generate_simple_rot_vec + end subroutine collision_generate_disruption_rot_vec - module subroutine collision_generate_simple_vel_vec(collider) + module subroutine collision_generate_disruption_vel_vec(collider) implicit none class(collision_disruption), intent(inout) :: collider !! Collision system object - end subroutine collision_generate_simple_vel_vec + end subroutine collision_generate_disruption_vel_vec 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 13d97b4e1..d7ae591f1 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -700,11 +700,14 @@ module subroutine collision_util_set_spins(self) 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)) + do i = 1,self%nbody + self%rot(:,i) = self%rot(:,i) + Lresidual(:) / (self%nbody * self%mass(i) * self%radius(i)**2 * self%Ip(3,i)) end do end if + call self%get_angular_momentum() + Lresidual(:) = self%L_budget(:) - (self%Lorbit(:) + self%Lspin(:)) + return end subroutine collision_util_set_spins diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 6fe15633b..cff097b4b 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -93,9 +93,9 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur call self%set_budgets() ! Minimize difference between energy/momentum and budgets - call fraggle_generate_minimize(self, lfailure_local) - call fraggle_generate_tan_vel(self, lfailure_local) - ! call fraggle_generate_rad_vel(self, lfailure_local) + ! call fraggle_generate_minimize(self, lfailure_local) + ! call fraggle_generate_tan_vel(self, lfailure_local) + call fraggle_generate_rad_vel(self, lfailure_local) call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) dEtot = self%Etot(2) - self%Etot(1) @@ -188,6 +188,7 @@ subroutine fraggle_generate_minimize(collider, lfailure) fval = E_objective_function(input_v) call minimize_bfgs(Efunc, nelem, input_v, tol, MAXLOOP, lfailure, output_v) + fval = E_objective_function(output_v) fragments%vmag(1:nfrag) = output_v(1:nfrag) do concurrent(i=1:nfrag) fragments%vc(:,i) = abs(fragments%vmag(i)) * fragments%v_r_unit(:,i) @@ -287,7 +288,7 @@ subroutine fraggle_generate_tan_vel(collider, lfailure) type(lambda_obj_err) :: objective_function real(DP), dimension(NDIM) :: Li, L_remainder, L_frag_tot character(len=STRMAX) :: message - real(DP), dimension(:), allocatable :: output_v + real(DP), dimension(:), allocatable :: v_t_output logical :: lfirst_func associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) @@ -310,10 +311,11 @@ subroutine fraggle_generate_tan_vel(collider, lfailure) ! ! Find the local kinetic energy minimum for the system that conserves linear and angular momentum objective_function = lambda_obj(tangential_objective_function, lfailure) + fval = tangential_objective_function(v_t_output(:), lfailure) - call minimize_bfgs(objective_function, nfrag-6, v_t_initial(7:nfrag), tol, MAXLOOP, lfailure, output_v) - fval = tangential_objective_function(output_v(:), lfailure) - fragments%v_t_mag(7:nfrag) = output_v(:) + call minimize_bfgs(objective_function, nfrag-6, v_t_initial(7:nfrag), tol, MAXLOOP, lfailure, v_t_output) + fval = tangential_objective_function(v_t_initial(:), lfailure) + fragments%v_t_mag(7:nfrag) = v_t_output(:) ! Now that the KE-minimized values of the i>6 fragments are found, calculate the momentum-conserving solution for tangential velociteis v_t_initial(7:nfrag) = fragments%v_t_mag(7:nfrag) @@ -474,33 +476,35 @@ subroutine fraggle_generate_rad_vel(collider, lfailure) real(DP), parameter :: TOL_INIT = 1e-14_DP real(DP), parameter :: VNOISE_MAG = 1e-10_DP !! Magnitude of the noise to apply to initial conditions to help minimizer find a solution in case of failure integer(I4B), parameter :: MAXLOOP = 100 - real(DP) :: ke_radial, tol + real(DP) :: ke_radial, tol, fval integer(I4B) :: i real(DP), dimension(:), allocatable :: v_r_initial real(DP), dimension(collider%fragments%nbody) :: vnoise type(lambda_obj) :: objective_function character(len=STRMAX) :: message - real(DP), dimension(:), allocatable :: output_v + real(DP), dimension(:), allocatable :: v_r_output associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) select type(fragments => collider%fragments) class is (fraggle_fragments(*)) ! Set the "target" ke for the radial component ke_radial = fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit + + do i = 1, nfrag + fragments%v_t_mag(i) = dot_product(fragments%vc(:,i), fragments%v_t_unit(:,i)) + fragments%v_r_mag(i) = dot_product(fragments%vc(:,i), fragments%v_r_unit(:,i)) + end do allocate(v_r_initial, source=fragments%v_r_mag) - ! Initialize radial velocity magnitudes with a random value that related to equipartition of kinetic energy with some noise - call random_number(vnoise(1:nfrag)) - vnoise(:) = 1.0_DP + VNOISE_MAG * (2 * vnoise(:) - 1.0_DP) - v_r_initial(1:nfrag) = sqrt(abs(2 * ke_radial) / (fragments%mass(1:nfrag) * nfrag)) * vnoise(1:nfrag) ! Initialize the lambda function using a structure constructor that calls the init method ! Minimize the ke objective function using the BFGS optimizer objective_function = lambda_obj(radial_objective_function) tol = TOL_INIT do while(tol < TOL_MIN) - call minimize_bfgs(objective_function, nfrag, v_r_initial, tol, MAXLOOP, lfailure, output_v) - fragments%v_r_mag = output_v + fval = radial_objective_function(v_r_initial) + call minimize_bfgs(objective_function, nfrag, v_r_initial, tol, MAXLOOP, lfailure, v_r_output) + fragments%v_r_mag = v_r_output if (.not.lfailure) exit tol = tol * 2 ! Keep increasing the tolerance until we converge on a solution v_r_initial(:) = fragments%v_r_mag(:) @@ -515,7 +519,7 @@ subroutine fraggle_generate_rad_vel(collider, lfailure) fragments%v_t_mag(1:nfrag), fragments%v_t_unit(:,1:nfrag), fragments%mass(1:nfrag), impactors%vbcom(:)) do i = 1, nfrag fragments%vc(:, i) = fragments%vb(:, i) - impactors%vbcom(:) - fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * dot_product(fragments%vb(:, i), fragments%vb(:, i)) + fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * dot_product(fragments%vc(:, i), fragments%vc(:, i)) end do fragments%ke_orbit = 0.5_DP * fragments%ke_orbit @@ -559,6 +563,7 @@ function radial_objective_function(v_r_mag_input) result(fval) allocate(v_shift, mold=fragments%vb) v_shift(:,:) = fraggle_util_vmag_to_vb(v_r_mag_input, fragments%v_r_unit, fragments%v_t_mag, fragments%v_t_unit, fragments%mass, impactors%vbcom) do i = 1,fragments%nbody + v_shift(:,i) = v_shift(:,i) - impactors%vbcom(:) rotmag2 = fragments%rot(1,i)**2 + fragments%rot(2,i)**2 + fragments%rot(3,i)**2 vmag2 = v_shift(1,i)**2 + v_shift(2,i)**2 + v_shift(3,i)**2 kearr(i) = fragments%mass(i) * (fragments%Ip(3, i) * fragments%radius(i)**2 * rotmag2 + vmag2)