From 3015359e93e1db4bddc84af4388ea7a2a59c9950 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 28 Dec 2022 14:32:58 -0500 Subject: [PATCH] Put the spin/radial/tangential velocitiy solvers back into Fraggle. Now that I have a much better handle on the initial guess, I will attempt to resurrect these solvers to get better constraints on energy and momentum --- src/fraggle/fraggle_generate.f90 | 376 +++++++++++++++++++++++++++++++ 1 file changed, 376 insertions(+) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 4ed5efac1..98ac872b5 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -255,4 +255,380 @@ end function E_objective_function end subroutine fraggle_generate_minimize + + subroutine fraggle_generate_spins(collider, f_spin, lfailure) + !! 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 without blowing the fragment kinetic energy budget. + !! + !! A failure will trigger a restructuring of the fragments so we will try new values of the radial position distribution. + implicit none + ! Arguments + class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object + real(DP), intent(in) :: f_spin !! Fraction of energy or momentum that goes into spin (whichever gives the lowest kinetic energy) + logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds! + ! Internals + real(DP), dimension(NDIM) :: L_remainder, rot_L, rot_ke + integer(I4B) :: i + character(len=STRMAX) :: message + + associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) + select type(fragments => collider%fragments) + class is (fraggle_fragments(*)) + lfailure = .false. + + ! Start the first two bodies with the same rotation as the original two impactors, then distribute the remaining angular momentum among the rest + L_remainder(:) = fragments%L_budget(:) + fragments%rot(:,:) = 0.0_DP + + fragments%ke_spin = 0.0_DP + if (norm2(L_remainder(:)) > FRAGGLE_LTOL) then + do i = 1, nfrag + ! Convert a fraction (f_spin) of either the remaining angular momentum or kinetic energy budget into spin, whichever gives the smaller rotation so as not to blow any budgets + rot_ke(:) = sqrt(2 * f_spin * fragments%ke_budget / (nfrag * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3, i))) & + * L_remainder(:) / norm2(L_remainder(:)) + rot_L(:) = f_spin * L_remainder(:) / (nfrag * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3, i)) + if (norm2(rot_ke) < norm2(rot_L)) then + fragments%rot(:,i) = rot_ke(:) + else + fragments%rot(:, i) = rot_L(:) + end if + fragments%ke_spin = fragments%ke_spin + fragments%mass(i) * fragments%Ip(3, i) * fragments%radius(i)**2 & + * dot_product(fragments%rot(:, i), fragments%rot(:, i)) + end do + end if + fragments%ke_spin = 0.5_DP * fragments%ke_spin + + lfailure = ((fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit) < 0.0_DP) + + if (lfailure) then + call swiftest_io_log_one_message(COLLISION_LOG_OUT, " ") + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Spin failure diagnostics") + write(message, *) fragments%ke_budget + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_budget : " // trim(adjustl(message))) + write(message, *) fragments%ke_spin + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_spin : " // trim(adjustl(message))) + write(message, *) fragments%ke_orbit + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_orbit : " // trim(adjustl(message))) + write(message, *) fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_remainder : " // trim(adjustl(message))) + end if + + end select + end associate + + return + end subroutine fraggle_generate_spins + + + subroutine fraggle_generate_tan_vel(collider, lfailure) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Adjusts the tangential velocities and spins of a collection of fragments such that they conserve angular momentum without blowing the fragment kinetic energy budget. + !! This procedure works in several stages, with a goal to solve the angular and linear momentum constraints on the fragments, while still leaving a positive balance of + !! our fragment kinetic energy (fragments%ke_budget) that we can put into the radial velocity distribution. + !! + !! The first thing we'll try to do is solve for the tangential velocities of the first 6 fragments, using angular and linear momentum as constraints and an initial + !! tangential velocity distribution for the remaining bodies (if there are any) that distributes their angular momentum equally between them. + !! If that doesn't work and we blow our kinetic energy budget, we will attempt to find a tangential velocity distribution that minimizes the kinetic energy while + !! conserving momentum. + !! + !! A failure will trigger a restructuring of the fragments so we will try new values of the radial position distribution. + implicit none + ! Arguments + class(collision_fraggle), intent(inout) :: collider !! Fraggle fragment system object + logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds + ! Internals + integer(I4B) :: i + real(DP), parameter :: TOL_MIN = 1e-1_DP ! This doesn't have to be very accurate, as we really just want a tangential velocity distribution with less kinetic energy than our initial guess. + real(DP), parameter :: TOL_INIT = 1e-14_DP + real(DP), parameter :: VNOISE_MAG = 1e-3_DP !! Magnitude of the noise to apply to initial conditions to help minimizer find a solution in case of failure + integer(I4B), parameter :: MAXLOOP = 10 + real(DP) :: tol + real(DP), dimension(:), allocatable :: v_t_initial + real(DP), dimension(collider%fragments%nbody) :: kefrag, vnoise + 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 + + associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) + select type(fragments => collider%fragments) + class is (fraggle_fragments(*)) + lfailure = .false. + + allocate(v_t_initial, mold=fragments%v_t_mag) + v_t_initial(:) = 0.0_DP + fragments%v_t_unit(:,:) = 0.0_DP + + ! Next we will solve for the tangential component of the velocities that both conserves linear momentum and uses the remaining angular momentum not used in spin. + ! This will be done using a linear solver that solves for the tangential velocities of the first 6 fragments, constrained by the linear and angular momentum vectors, + ! which is embedded in a non-linear minimizer that will adjust the tangential velocities of the remaining i>6 fragments to minimize kinetic energy for a given momentum solution + ! The initial conditions fed to the minimizer for the fragments will be the remaining angular momentum distributed between the fragments. + call fragments%get_angular_momentum() + L_remainder(:) = fragments%L_budget(:) - fragments%Lspin(:) + do i = 1, nfrag + v_t_initial(i) = norm2(L_remainder(:)) / ((nfrag - i + 1) * fragments%mass(i) * norm2(fragments%v_r_unit(:,i))) + Li(:) = fragments%mass(i) * (fragments%v_r_unit(:,i) .cross. (v_t_initial(i) * fragments%v_t_unit(:, i))) + L_remainder(:) = L_remainder(:) - Li(:) + end do + + ! Find the local kinetic energy minimum for the system that conserves linear and angular momentum + objective_function = lambda_obj(tangential_objective_function, lfailure) + + tol = TOL_INIT + do while(tol < TOL_MIN) + call minimize_bfgs(objective_function, nfrag-6, v_t_initial(7:nfrag), tol, MAXLOOP, lfailure, output_v) + fragments%v_t_mag(7:nfrag) = output_v(:) + ! 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) + if (.not.lfailure) exit + tol = tol * 2_DP ! Keep increasing the tolerance until we converge on a solution + call random_number(vnoise(1:nfrag)) ! Adding a bit of noise to the initial conditions helps it find a solution more often + vnoise(:) = 1.0_DP + VNOISE_MAG * (2 * vnoise(:) - 1._DP) + v_t_initial(:) = v_t_initial(:) * vnoise(:) + 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 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) + fragments%v_t_unit(:,i) = fragments%vb(:,i) - impactors%vbcom(:) + end do + + ! Now do a kinetic energy budget check to make sure we are still within the budget. + kefrag = 0.0_DP + do concurrent(i = 1:nfrag) + kefrag(i) = fragments%mass(i) * dot_product(fragments%vb(:, i), fragments%vb(:, i)) + end do + fragments%ke_orbit = 0.5_DP * sum(kefrag(:)) + + ! If we are over the energy budget, flag this as a failure so we can try again + lfailure = ((fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit) < 0.0_DP) + if (lfailure) then + call swiftest_io_log_one_message(COLLISION_LOG_OUT, " ") + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Tangential velocity failure diagnostics") + call fragments%get_angular_momentum() + L_frag_tot = fragments%Lspin(:) + fragments%Lorbit(:) + write(message, *) .mag.(L_frag_tot(:)) / (.mag.fragments%L_budget(:)) + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "|L_remainder| : " // trim(adjustl(message))) + write(message, *) fragments%ke_budget + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_budget : " // trim(adjustl(message))) + write(message, *) fragments%ke_spin + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_spin : " // trim(adjustl(message))) + write(message, *) fragments%ke_orbit + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_tangential : " // trim(adjustl(message))) + write(message, *) fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_radial : " // trim(adjustl(message))) + end if + end select + end associate + + return + + contains + function solve_fragment_tan_vel(lfailure, v_t_mag_input) result(v_t_mag_output) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Adjusts the positions, velocities, and spins of a collection of fragments such that they conserve angular momentum + implicit none + ! Arguments + logical, intent(out) :: lfailure !! Error flag + real(DP), dimension(:), optional, intent(in) :: v_t_mag_input !! Unknown tangential velocities for fragments 7:nfrag + ! Internals + integer(I4B) :: i + ! Result + real(DP), dimension(:), allocatable :: v_t_mag_output + + real(DP), dimension(2 * NDIM, 2 * NDIM) :: A ! LHS of linear equation used to solve for momentum constraint in Gauss elimination code + real(DP), dimension(2 * NDIM) :: b ! RHS of linear equation used to solve for momentum constraint in Gauss elimination code + real(DP), dimension(NDIM) :: L_lin_others, L_orb_others, L, vtmp + + associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) + select type(fragments => collider%fragments) + class is (fraggle_fragments(*)) + lfailure = .false. + ! We have 6 constraint equations (2 vector constraints in 3 dimensions each) + ! The first 3 are that the linear momentum of the fragments is zero with respect to the collisional barycenter + ! The second 3 are that the sum of the angular momentum of the fragments is conserved from the pre-impact state + L_lin_others(:) = 0.0_DP + L_orb_others(:) = 0.0_DP + do i = 1, nfrag + if (i <= 2 * NDIM) then ! The tangential velocities of the first set of bodies will be the unknowns we will solve for to satisfy the constraints + A(1:3, i) = fragments%mass(i) * fragments%v_t_unit(:, i) + A(4:6, i) = fragments%mass(i) * fragments%rmag(i) * (fragments%v_r_unit(:, i) .cross. fragments%v_t_unit(:, i)) + else if (present(v_t_mag_input)) then + vtmp(:) = v_t_mag_input(i - 6) * fragments%v_t_unit(:, i) + L_lin_others(:) = L_lin_others(:) + fragments%mass(i) * vtmp(:) + L(:) = fragments%mass(i) * (fragments%v_r_unit(:, i) .cross. vtmp(:)) + L_orb_others(:) = L_orb_others(:) + L(:) + end if + end do + b(1:3) = -L_lin_others(:) + b(4:6) = fragments%L_budget(:) - fragments%Lspin(:) - L_orb_others(:) + allocate(v_t_mag_output(nfrag)) + v_t_mag_output(1:6) = solve_linear_system(A, b, 6, lfailure) + if (present(v_t_mag_input)) v_t_mag_output(7:nfrag) = v_t_mag_input(:) + end select + end associate + return + end function solve_fragment_tan_vel + + + function tangential_objective_function(v_t_mag_input, lfailure) result(fval) + !! Author: David A. Minton + !! + !! Objective function for evaluating how close our fragment velocities get to minimizing KE error from our required value + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: v_t_mag_input !! Unknown tangential component of velocity vector set previously by angular momentum constraint + logical, intent(out) :: lfailure !! Error flag + ! Result + real(DP) :: fval + ! Internals + integer(I4B) :: i + real(DP), dimension(NDIM,collider%fragments%nbody) :: v_shift + real(DP), dimension(collider%fragments%nbody) :: v_t_new, kearr + real(DP) :: keo + + associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) + select type(fragments => collider%fragments) + class is (fraggle_fragments(*)) + lfailure = .false. + + v_t_new(:) = solve_fragment_tan_vel(v_t_mag_input=v_t_mag_input(:), lfailure=lfailure) + v_shift(:,:) = fraggle_util_vmag_to_vb(fragments%v_r_mag, fragments%v_r_unit, v_t_new, fragments%v_t_unit, fragments%mass, impactors%vbcom) + + kearr = 0.0_DP + do concurrent(i = 1:nfrag) + kearr(i) = fragments%mass(i) * dot_product(v_shift(:, i), v_shift(:, i)) + end do + keo = 0.5_DP * sum(kearr(:)) + fval = keo + lfailure = .false. + end select + end associate + + return + end function tangential_objective_function + + end subroutine fraggle_generate_tan_vel + + + subroutine fraggle_generate_rad_vel(collider, lfailure) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! + !! Adjust the fragment velocities to set the fragment orbital kinetic energy. This will minimize the difference between the fragment kinetic energy and the energy budget + implicit none + ! Arguments + class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object + logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds! + ! Internals + real(DP), parameter :: TOL_MIN = FRAGGLE_ETOL ! This needs to be more accurate than the tangential step, as we are trying to minimize the total residual energy + 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 + 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 + + 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 + + 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 + if (.not.lfailure) exit + tol = tol * 2 ! Keep increasing the tolerance until we converge on a solution + v_r_initial(:) = fragments%v_r_mag(:) + call random_number(vnoise(1:nfrag)) ! Adding a bit of noise to the initial conditions helps it find a solution more often + vnoise(:) = 1.0_DP + VNOISE_MAG * (2 * vnoise(:) - 1._DP) + v_r_initial(:) = v_r_initial(:) * vnoise(:) + end do + + ! 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%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 i = 1, nfrag + fragments%v_t_unit(:, i) = fragments%vb(:, i) - impactors%vbcom(:) + fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * dot_product(fragments%vb(:, i), fragments%vb(:, i)) + end do + fragments%ke_orbit = 0.5_DP * fragments%ke_orbit + + lfailure = abs((fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin)) / fragments%ke_budget) > FRAGGLE_ETOL + if (lfailure) then + call swiftest_io_log_one_message(COLLISION_LOG_OUT, " ") + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Radial velocity failure diagnostics") + write(message, *) fragments%ke_budget + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_budget : " // trim(adjustl(message))) + write(message, *) fragments%ke_spin + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_spin : " // trim(adjustl(message))) + write(message, *) fragments%ke_orbit + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_orbit : " // trim(adjustl(message))) + write(message, *) fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin) + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_remainder : " // trim(adjustl(message))) + end if + end select + end associate + return + + contains + function radial_objective_function(v_r_mag_input) result(fval) + !! Author: David A. Minton + !! + !! Objective function for evaluating how close our fragment velocities get to minimizing KE error from our required value + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: v_r_mag_input !! Unknown radial component of fragment velocity vector + ! Result + real(DP) :: fval !! The objective function result, which is the square of the difference between the calculated fragment kinetic energy and our target + !! Minimizing this brings us closer to our objective + ! Internals + integer(I4B) :: i + real(DP), dimension(:,:), allocatable :: v_shift + real(DP), dimension(collider%fragments%nbody) :: kearr + real(DP) :: keo, ke_radial, rotmag2, vmag2 + + associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) + select type(fragments => collider%fragments) + class is (fraggle_fragments(*)) + 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 + 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) + end do + keo = 2 * fragments%ke_budget - sum(kearr(:)) + ke_radial = fragments%ke_budget - fragments%ke_orbit - fragments%ke_spin + ! The following ensures that fval = 0 is a local minimum, which is what the BFGS method is searching for + fval = (keo / (2 * ke_radial))**2 + end select + end associate + + return + end function radial_objective_function + + end subroutine fraggle_generate_rad_vel + end submodule s_fraggle_generate