Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Improvements to the energy budget calcs and bug fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Dec 27, 2022
1 parent 112e048 commit 6ed292e
Show file tree
Hide file tree
Showing 6 changed files with 113 additions and 189 deletions.
21 changes: 19 additions & 2 deletions src/collision/collision_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand Down Expand Up @@ -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
Expand Down
139 changes: 85 additions & 54 deletions src/collision/collision_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
!!
Expand All @@ -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
Expand Down
11 changes: 5 additions & 6 deletions src/fraggle/fraggle_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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

Expand Down Expand Up @@ -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
Expand Down
27 changes: 1 addition & 26 deletions src/fraggle/fraggle_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
Loading

0 comments on commit 6ed292e

Please sign in to comment.