From 6ed292e69d924a1baf549b8bd61be3bf8b21ef30 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 27 Dec 2022 11:50:43 -0500 Subject: [PATCH] Improvements to the energy budget calcs and bug fixes --- src/collision/collision_module.f90 | 21 ++++- src/collision/collision_util.f90 | 139 ++++++++++++++++++----------- src/fraggle/fraggle_generate.f90 | 11 ++- src/fraggle/fraggle_module.f90 | 27 +----- src/fraggle/fraggle_util.f90 | 102 +-------------------- src/misc/minimizer_module.f90 | 2 +- 6 files changed, 113 insertions(+), 189 deletions(-) diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index cfac7b49f..ed94b6be1 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -110,9 +110,15 @@ module collision real(DP), dimension(NDIM,nbody) :: v_t_unit !! Array of tangential direction unit vectors of individual fragments in the collisional coordinate frame real(DP), dimension(NDIM,nbody) :: v_n_unit !! Array of normal direction unit vectors of individual fragments in the collisional coordinate frame integer(I1B), dimension(nbody) :: origin_body !! Array of indices indicating which impactor body (1 or 2) the fragment originates from + real(DP), dimension(NDIM) :: Lorbit !! Orbital angular momentum vector of all fragments + 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 contains - procedure :: reset => collision_util_reset_fragments !! Deallocates all allocatable arrays and sets everything else to 0 - final :: collision_final_fragments !! Finalizer deallocates all allocatables + 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 + final :: collision_final_fragments !! Finalizer deallocates all allocatables end type collision_fragments @@ -404,6 +410,17 @@ 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 + 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 + end subroutine collision_util_get_kinetic_energy + module subroutine collision_util_reset_fragments(self) implicit none class(collision_fragments(*)), intent(inout) :: self diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 4a8ee6b74..f182b6dfb 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -175,6 +175,59 @@ 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 + !! + !! Calculates the current angular momentum of the fragments + implicit none + ! Arguments + class(collision_fragments(*)), intent(inout) :: self !! Fraggle fragment system object + ! Internals + integer(I4B) :: i + + associate(fragments => self, nfrag => self%nbody) + fragments%Lorbit(:) = 0.0_DP + fragments%Lspin(:) = 0.0_DP + + do i = 1, nfrag + fragments%Lorbit(:) = fragments%Lorbit(:) + fragments%mass(i) * (fragments%rc(:, i) .cross. fragments%vc(:, i)) + fragments%Lspin(:) = fragments%Lspin(:) + fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:, i) * fragments%rot(:, i) + end do + end associate + + return + end subroutine collision_util_get_angular_momentum + + + module subroutine collision_util_get_kinetic_energy(self) + !! Author: David A. Minton + !! + !! Calculates the current kinetic energy of the fragments + implicit none + ! Argument + class(collision_fragments(*)), intent(inout) :: self !! Fraggle fragment system object + ! Internals + integer(I4B) :: i + + associate(fragments => self, nfrag => self%nbody) + fragments%ke_orbit = 0.0_DP + fragments%ke_spin = 0.0_DP + + do i = 1, nfrag + fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i)) + fragments%ke_spin = fragments%ke_spin + fragments%mass(i) * fragments%Ip(3,i) * dot_product(fragments%rot(:,i),fragments%rot(:,i) ) + end do + + fragments%ke_orbit = fragments%ke_orbit / 2 + fragments%ke_spin = fragments%ke_spin / 2 + + end associate + + return + end subroutine collision_util_get_kinetic_energy + + module subroutine collision_util_get_energy_momentum(self, nbody_system, param, lbefore) !! Author: David A. Minton !! @@ -189,71 +242,49 @@ module subroutine collision_util_get_energy_momentum(self, nbody_system, param, class(base_parameters), intent(inout) :: param !! Current swiftest run configuration parameters logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the nbody_system, with impactors included and fragments excluded or vice versa ! Internals - class(base_nbody_system), allocatable, save :: tmpsys - class(base_parameters), allocatable, save :: tmpparam - integer(I4B) :: npl_before, npl_after, stage,i + integer(I4B) :: stage,i + real(DP), dimension(NDIM) :: Lorbit, Lspin + real(DP) :: ke_orbit, ke_spin + select type(nbody_system) class is (swiftest_nbody_system) select type(param) class is (swiftest_parameters) associate(fragments => self%fragments, impactors => self%impactors, nfrag => self%fragments%nbody, pl => nbody_system%pl, cb => nbody_system%cb) - ! Because we're making a copy of the massive body object with the excludes/fragments appended, we need to deallocate the - ! big k_plpl array and recreate it when we're done, otherwise we run the risk of blowing up the memory by - ! allocating two of these ginormous arrays simulteouously. This is not particularly efficient, but as this - ! subroutine should be called relatively infrequently, it shouldn't matter too much. - - npl_before = pl%nbody - npl_after = npl_before + nfrag if (lbefore) then - call self%construct_temporary_system(nbody_system, param, tmpsys, tmpparam) - select type(tmpsys) - class is (swiftest_nbody_system) - ! Build the exluded body logical mask for the *before* case: Only the original bodies are used to compute energy and momentum - tmpsys%pl%status(impactors%id(1:impactors%ncoll)) = ACTIVE - tmpsys%pl%status(npl_before+1:npl_after) = INACTIVE - end select - else - if (.not.allocated(tmpsys)) then - write(*,*) "Error in collision_util_get_energy_momentum. " // & - " This must be called with lbefore=.true. at least once before calling it with lbefore=.false." - call util_exit(FAILURE) - end if - select type(tmpsys) - class is (swiftest_nbody_system) - ! Build the exluded body logical mask for the *after* case: Only the new bodies are used to compute energy and momentum - call self%add_fragments(tmpsys, tmpparam) - tmpsys%pl%status(impactors%id(1:impactors%ncoll)) = INACTIVE - do concurrent(i = npl_before+1:npl_after) - tmpsys%pl%status(i) = ACTIVE - tmpsys%pl%rot(:,i) = 0.0_DP - tmpsys%pl%vb(:,i) = impactors%vbcom(:) - end do - end select - end if - select type(tmpsys) - class is (swiftest_nbody_system) - if (param%lflatten_interactions) call tmpsys%pl%flatten(param) - - call tmpsys%get_energy_and_momentum(param) + Lorbit(:) = sum(impactors%Lorbit(:,:), dim=2) + Lspin(:) = sum(impactors%Lspin(:,:), dim=2) + ke_orbit = 0.0_DP + ke_spin = 0.0_DP + do concurrent(i = 1:2) + ke_orbit = ke_orbit + impactors%mass(i) * dot_product(impactors%vc(:,i), impactors%vc(:,i)) + ke_spin = ke_spin + impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i) * dot_product(impactors%rot(:,i), impactors%rot(:,i)) + end do + else + call fragments%get_angular_momentum() + Lorbit(:) = fragments%Lorbit(:) + Lspin(:) = fragments%Lspin(:) + call fragments%get_kinetic_energy() + ke_orbit = fragments%ke_orbit + ke_spin = fragments%ke_spin - ! Calculate the current fragment energy and momentum balances - if (lbefore) then - stage = 1 - else - stage = 2 - end if - self%Lorbit(:,stage) = tmpsys%Lorbit(:) - self%Lspin(:,stage) = tmpsys%Lspin(:) - self%Ltot(:,stage) = tmpsys%Ltot(:) - self%ke_orbit(stage) = tmpsys%ke_orbit - self%ke_spin(stage) = tmpsys%ke_spin - self%pe(stage) = tmpsys%pe - self%Etot(stage) = tmpsys%ke_orbit + tmpsys%ke_spin - end select + end if + ! Calculate the current fragment energy and momentum balances + if (lbefore) then + stage = 1 + else + stage = 2 + end if + self%Lorbit(:,stage) = Lorbit(:) + self%Lspin(:,stage) = Lspin(:) + self%Ltot(:,stage) = Lorbit(:) + Lspin(:) + self%ke_orbit(stage) = ke_orbit + self%ke_spin(stage) = ke_spin + self%Etot(stage) = ke_orbit + ke_spin end associate end select end select diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 7ca7cac06..2738039e3 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -80,7 +80,6 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur write(message,*) try call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle try " // trim(adjustl(message))) if (lfailure_local) then - !call fragments%restructure(impactors, try, f_spin, r_max_start) call fragments%reset() try = try + 1 end if @@ -163,7 +162,7 @@ subroutine fraggle_generate_minimize(collider, lfailure) ! Internals real(DP), parameter :: TOL_MIN = 1.0e-6_DP real(DP), parameter :: TOL_INIT = 1e-8_DP - integer(I4B), parameter :: MAXLOOP = 30 + integer(I4B), parameter :: MAXLOOP = 50 real(DP), dimension(collider%fragments%nbody) :: input_v real(DP), dimension(:), allocatable :: output_v type(lambda_obj) :: Efunc @@ -174,7 +173,7 @@ subroutine fraggle_generate_minimize(collider, lfailure) select type(fragments => collider%fragments) class is (fraggle_fragments(*)) - nelem = 2 * (nfrag - 1) + nelem = nfrag lfailure = .false. ! Find the local kinetic energy minimum for the nbody_system that conserves linear and angular momentum Efunc = lambda_obj(E_objective_function) @@ -193,7 +192,7 @@ subroutine fraggle_generate_minimize(collider, lfailure) fragments%vmag(1:nfrag) = output_v(1:nfrag) - do concurrent(i=2:nfrag) + do concurrent(i=1:nfrag) fragments%vc(:,i) = abs(fragments%vmag(i)) * fragments%v_r_unit(:,i) end do @@ -240,8 +239,8 @@ function E_objective_function(val_input) result(fval) call fraggle_generate_set_spins(tmp_frag) ! Get the current kinetic energy of the system - call fragments%get_kinetic_energy() - deltaE = fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin) + call tmp_frag%get_kinetic_energy() + deltaE = tmp_frag%ke_budget - (tmp_frag%ke_orbit + tmp_frag%ke_spin) ! Use the deltaE as the basis of our objective function, with a higher penalty for having excess kinetic energy compared with having a deficit if (deltaE < 0.0_DP) then diff --git a/src/fraggle/fraggle_module.f90 b/src/fraggle/fraggle_module.f90 index f2ad578c4..a40e196b4 100644 --- a/src/fraggle/fraggle_module.f90 +++ b/src/fraggle/fraggle_module.f90 @@ -20,17 +20,11 @@ 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), dimension(NDIM) :: Lorbit !! Orbital angular momentum vector of all fragments - 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 :: get_angular_momentum => fraggle_util_get_angular_momentum !! Calcualtes the current angular momentum of the fragments - procedure :: get_kinetic_energy => fraggle_util_get_kinetic_energy !! Calcualtes the current kinetic energy of the fragments + 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) - procedure :: restructure => fraggle_util_restructure !! Restructure the inputs after a failed attempt failed to find a set of positions and velocities that satisfy the energy and momentum constraints final :: fraggle_final_fragments !! Finalizer will deallocate all allocatables end type fraggle_fragments @@ -80,16 +74,6 @@ module subroutine fraggle_util_construct_temporary_system(self, nbody_system, pa class(base_parameters), allocatable, intent(out) :: tmpparam !! Output temporary configuration run parameters end subroutine fraggle_util_construct_temporary_system - module subroutine fraggle_util_get_angular_momentum(self) - implicit none - class(fraggle_fragments(*)), intent(inout) :: self !! Fraggle fragment system object - end subroutine fraggle_util_get_angular_momentum - - module subroutine fraggle_util_get_kinetic_energy(self) - implicit none - class(fraggle_fragments(*)), intent(inout) :: self !! Fraggle fragment system object - end subroutine fraggle_util_get_kinetic_energy - module subroutine fraggle_util_reset_fragments(self) implicit none class(fraggle_fragments(*)), intent(inout) :: self @@ -100,15 +84,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_restructure(self, impactors, try, f_spin, r_max_start) - implicit none - class(fraggle_fragments(*)), intent(inout) :: self !! Fraggle fragment system object - class(collision_impactors), intent(in) :: impactors !! Fraggle collider system object - integer(I4B), intent(in) :: try !! The current number of times Fraggle has tried to find a solution - real(DP), intent(inout) :: f_spin !! Fraction of energy/momentum that goes into spin. This decreases ater a failed attempt - real(DP), intent(inout) :: r_max_start !! The maximum radial distance that the position calculation starts with. This increases after a failed attempt - end subroutine fraggle_util_restructure - module subroutine fraggle_util_set_budgets(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 07d845053..f85eed056 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -37,57 +37,6 @@ module subroutine fraggle_util_construct_temporary_system(self, nbody_system, pa end subroutine fraggle_util_construct_temporary_system - module subroutine fraggle_util_get_angular_momentum(self) - !! Author: David A. Minton - !! - !! Calculates the current angular momentum of the fragments - implicit none - ! Arguments - class(fraggle_fragments(*)), intent(inout) :: self !! Fraggle fragment system object - ! Internals - integer(I4B) :: i - - associate(fragments => self, nfrag => self%nbody) - fragments%Lorbit(:) = 0.0_DP - fragments%Lspin(:) = 0.0_DP - - do i = 1, nfrag - fragments%Lorbit(:) = fragments%Lorbit(:) + fragments%mass(i) * (fragments%rc(:, i) .cross. fragments%vc(:, i)) - fragments%Lspin(:) = fragments%Lspin(:) + fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:, i) * fragments%rot(:, i) - end do - end associate - - return - end subroutine fraggle_util_get_angular_momentum - - - module subroutine fraggle_util_get_kinetic_energy(self) - !! Author: David A. Minton - !! - !! Calculates the current kinetic energy of the fragments - implicit none - ! Argument - class(fraggle_fragments(*)), intent(inout) :: self !! Fraggle fragment system object - ! Internals - integer(I4B) :: i - - associate(fragments => self, nfrag => self%nbody) - fragments%ke_orbit = 0.0_DP - fragments%ke_spin = 0.0_DP - - do i = 1, nfrag - fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i)) - fragments%ke_spin = fragments%ke_spin + fragments%mass(i) * fragments%Ip(3,i) * dot_product(fragments%rot(:,i),fragments%rot(:,i) ) - end do - - fragments%ke_orbit = fragments%ke_orbit / 2 - fragments%ke_spin = fragments%ke_spin / 2 - - end associate - - return - end subroutine fraggle_util_get_kinetic_energy - module subroutine fraggle_util_reset_fragments(self) !! author: David A. Minton !! @@ -138,53 +87,6 @@ module subroutine fraggle_util_reset_system(self) end subroutine fraggle_util_reset_system - module subroutine fraggle_util_restructure(self, impactors, try, f_spin, r_max_start) - !! Author: David A. Minton - !! - !! Restructure the inputs after a failed attempt failed to find a set of positions and velocities that satisfy the energy and momentum constraints - implicit none - ! Arguments - class(fraggle_fragments(*)), intent(inout) :: self !! Fraggle fragment system object - class(collision_impactors), intent(in) :: impactors !! Fraggle collider system object - integer(I4B), intent(in) :: try !! The current number of times Fraggle has tried to find a solution - real(DP), intent(inout) :: f_spin !! Fraction of energy/momentum that goes into spin. This decreases ater a failed attempt - real(DP), intent(inout) :: r_max_start !! The maximum radial distance that the position calculation starts with. This increases after a failed attempt - ! Internals - real(DP), save :: ke_tot_deficit, r_max_start_old, ke_avg_deficit_old - real(DP) :: delta_r, delta_r_max, ke_avg_deficit - real(DP), parameter :: ke_avg_deficit_target = 0.0_DP - - ! Introduce a bit of noise in the radius determination so we don't just flip flop between similar failed positions - associate(fragments => self) - call random_number(delta_r_max) - delta_r_max = sum(impactors%radius(:)) * (1.0_DP + 2e-1_DP * (delta_r_max - 0.5_DP)) - if (try == 1) then - ke_tot_deficit = - (fragments%ke_budget - fragments%ke_orbit - fragments%ke_spin) - ke_avg_deficit = ke_tot_deficit - delta_r = delta_r_max - else - ! Linearly interpolate the last two failed solution ke deficits to find a new distance value to try - ke_tot_deficit = ke_tot_deficit - (fragments%ke_budget - fragments%ke_orbit - fragments%ke_spin) - ke_avg_deficit = ke_tot_deficit / try - delta_r = (r_max_start - r_max_start_old) * (ke_avg_deficit_target - ke_avg_deficit_old) & - / (ke_avg_deficit - ke_avg_deficit_old) - if (abs(delta_r) > delta_r_max) delta_r = sign(delta_r_max, delta_r) - end if - r_max_start_old = r_max_start - r_max_start = r_max_start + delta_r ! The larger lever arm can help if the problem is in the angular momentum step - ke_avg_deficit_old = ke_avg_deficit - - if (f_spin > epsilon(1.0_DP)) then ! Try reducing the fraction in spin - f_spin = f_spin / 2 - else - f_spin = 0.0_DP - end if - end associate - - return - end subroutine fraggle_util_restructure - - module subroutine fraggle_util_set_budgets(self) !! author: David A. Minton !! @@ -200,8 +102,8 @@ module subroutine fraggle_util_set_budgets(self) select type(fragments => self%fragments) class is (fraggle_fragments(*)) - dEtot = self%Etot(2) - self%Etot(1) - dL(:) = self%Ltot(:,2) - self%Ltot(:,1) + dEtot = self%Etot(1) + dL(:) = self%Ltot(:,1) fragments%L_budget(:) = -dL(:) fragments%ke_budget = -(dEtot - impactors%Qloss) diff --git a/src/misc/minimizer_module.f90 b/src/misc/minimizer_module.f90 index 667f061db..8beab3879 100644 --- a/src/misc/minimizer_module.f90 +++ b/src/misc/minimizer_module.f90 @@ -128,7 +128,7 @@ module subroutine minimize_bfgs(f, N, x0, eps, maxloop, lerr, x1) real(DP), dimension(:), intent(out), allocatable :: x1 ! Internals integer(I4B) :: i, j, k, l, conv - real(DP), parameter :: graddelta = 1e-6_DP !! Delta x for gradient calculations + real(DP), parameter :: graddelta = 1e-4_DP !! Delta x for gradient calculations real(DP), dimension(N) :: S !! Direction vectors real(DP), dimension(N,N) :: H !! Approximated inverse Hessian matrix real(DP), dimension(N) :: grad1 !! gradient of f