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

Commit

Permalink
Made substantials improvements to Fraggle for robustness and ability …
Browse files Browse the repository at this point in the history 
…to converge on the required momentum and energy constraints
  • Loading branch information
@daminton
daminton committed Jan 7, 2023
1 parent 036f2c9 commit 36ec14c
Show file tree
Hide file tree
Showing 13 changed files with 596 additions and 584 deletions.
16 changes: 8 additions & 8 deletions examples/Fragmentation/Fragmentation_Movie.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,16 +73,16 @@

rot_vectors = {"disruption_headon" : [np.array([0.0, 0.0, 0.0]),
np.array([0.0, 0.0, 0.0])],
"disruption_off_axis": [np.array([0.0, 0.0, -6.0e4]),
np.array([0.0, 0.0, 1.0e5])],
"disruption_off_axis": [np.array([0.0, 0.0, -6.0e3]),
np.array([0.0, 0.0, 1.0e4])],
"supercatastrophic_headon": [np.array([0.0, 0.0, 0.0]),
np.array([0.0, 0.0, 0.0])],
"supercatastrophic_off_axis": [np.array([0.0, 0.0, -6.0e4]),
np.array([0.0, 0.0, 1.0e5])],
"hitandrun_disrupt" : [np.array([0.0, 0.0, 6.0e4]),
np.array([0.0, 0.0, 1.0e5])],
"hitandrun_pure" : [np.array([0.0, 0.0, 6.0e4]),
np.array([0.0, 0.0, 1.0e5])]
"supercatastrophic_off_axis": [np.array([0.0, 0.0, -6.0e3]),
np.array([0.0, 0.0, 1.0e4])],
"hitandrun_disrupt" : [np.array([0.0, 0.0, 6.0e3]),
np.array([0.0, 0.0, 1.0e4])],
"hitandrun_pure" : [np.array([0.0, 0.0, 6.0e3]),
np.array([0.0, 0.0, 1.0e4])]
}

body_Gmass = {"disruption_headon" : [1e-7, 1e-10],
Expand Down
2 changes: 1 addition & 1 deletion src/base/base_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,7 +118,7 @@ subroutine abstract_io_dump_param(self, param_file_name)
import base_parameters
implicit none
class(base_parameters),intent(in) :: self !! Output collection of parameters
character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in)
character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in)
end subroutine abstract_io_dump_param

subroutine abstract_io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
Expand Down
10 changes: 5 additions & 5 deletions src/collision/collision_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t)
class(base_parameters), intent(inout) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! The time of the collision
! Internals
integer(I4B) :: i,j,nfrag
integer(I4B) :: i,j,nimp
real(DP), dimension(NDIM) :: vcom, rnorm

select type(nbody_system)
Expand All @@ -62,8 +62,8 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t)
allocate(before%pl, source=pl)
end select

nfrag = size(impactors%id(:))
do i = 1, nfrag
nimp = size(impactors%id(:))
do i = 1, nimp
j = impactors%id(i)
vcom(:) = pl%vb(:,j) - impactors%vbcom(:)
rnorm(:) = .unit. (impactors%rb(:,2) - impactors%rb(:,1))
Expand Down Expand Up @@ -183,7 +183,7 @@ module subroutine collision_generate_merge(self, nbody_system, param, t)
associate(fragments => nbody_system%collider%fragments)

! Calculate the initial energy of the nbody_system without the collisional family
call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.)
call self%get_energy_and_momentum(nbody_system, param, phase="before")

! The new body's metadata will be taken from the largest of the two impactor bodies, so we need
! its index in the main pl structure
Expand Down Expand Up @@ -216,7 +216,7 @@ module subroutine collision_generate_merge(self, nbody_system, param, t)
fragments%vc(:,1) = 0.0_DP

! Get the energy of the system after the collision
call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.)
call self%get_energy_and_momentum(nbody_system, param, phase="after")

! Update any encounter lists that have the removed bodies in them so that they instead point to the new body
do k = 1, nbody_system%plpl_encounter%nenc
Expand Down
78 changes: 46 additions & 32 deletions src/collision/collision_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,14 +121,10 @@ module collision
real(DP) :: ke_spin_tot !! Spin kinetic energy of all fragments
real(DP) :: pe !! Potential energy of all fragments
real(DP) :: be !! Binding energy of all fragments
real(DP) :: E_budget !! Kinetic energy budget for computing fragment trajectories
real(DP), dimension(NDIM) :: L_budget !! Angular momentum budget for computing fragment trajectories
real(DP), dimension(nbody) :: ke_orbit !! Orbital kinetic energy of each individual fragment
real(DP), dimension(nbody) :: ke_spin !! Spin kinetic energy of each individual fragment
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_energy => collision_util_get_energy !! Calcualtes the current kinetic energy of the fragments
procedure :: set_coordinate_system => collision_util_set_coordinate_fragments !! Sets the coordinate system of the fragments
final :: collision_final_fragments !! Finalizer deallocates all allocatables
end type collision_fragments
Expand All @@ -146,27 +142,37 @@ module collision
integer(I4B) :: status !! Status flag to pass to the collision list once the collision has been resolved
integer(I4B) :: collision_id !! ID number of this collision event

! For the following variables, index 1 refers to the *entire* n-body nbody_system in its pre-collisional state and index 2 refers to the nbody_system in its post-collisional state
real(DP), dimension(NDIM,2) :: L_orbit !! Before/after orbital angular momentum
real(DP), dimension(NDIM,2) :: L_spin !! Before/after spin angular momentum
real(DP), dimension(NDIM,2) :: L_total !! Before/after total nbody_system angular momentum
! Scale factors used to scale dimensioned quantities to a more "natural" system where important quantities (like kinetic energy, momentum) are of order ~1
real(DP) :: dscale = 1.0_DP !! Distance dimension scale factor
real(DP) :: mscale = 1.0_DP !! Mass scale factor
real(DP) :: tscale = 1.0_DP !! Time scale factor
real(DP) :: vscale = 1.0_DP !! Velocity scale factor (a convenience unit that is derived from dscale and tscale)
real(DP) :: Escale = 1.0_DP !! Energy scale factor (a convenience unit that is derived from dscale, tscale, and mscale)
real(DP) :: Lscale = 1.0_DP !! Angular momentum scale factor (a convenience unit that is derived from dscale, tscale, and mscale)

! For the following variables, index 1 refers to the *entire* n-body system in its pre-collisional state and index 2 refers to the system in its post-collisional state
real(DP), dimension(NDIM,2) :: L_orbit !! Before/after orbital angular momentum
real(DP), dimension(NDIM,2) :: L_spin !! Before/after spin angular momentum
real(DP), dimension(NDIM,2) :: L_total !! Before/after total nbody_system angular momentum
real(DP), dimension(2) :: ke_orbit !! Before/after orbital kinetic energy
real(DP), dimension(2) :: ke_spin !! Before/after spin kinetic energy
real(DP), dimension(2) :: pe !! Before/after potential energy
real(DP), dimension(2) :: be !! Before/after binding energy
real(DP), dimension(2) :: te !! Before/after total system energy

contains
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
procedure :: merge => collision_generate_merge !! Merges the impactors to make a single final body
procedure :: add_fragments => collision_util_add_fragments_to_collider !! Add fragments to nbody_system
procedure :: get_energy_and_momentum => collision_util_get_energy_and_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 :: setup => collision_util_setup_collider !! Initializer for the encounter collision system and the before/after snapshots
procedure :: setup_impactors => collision_util_setup_impactors_collider !! Initializer for the impactors for the encounter collision system. Deallocates old impactors before creating new ones
procedure :: setup_fragments => collision_util_setup_fragments_collider !! Initializer for the fragments of the collision system.
procedure :: add_fragments => collision_util_add_fragments_to_collider !! Add fragments to nbody_system
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
procedure :: merge => collision_generate_merge !! Merges the impactors to make a single final body
procedure :: set_natural_scale => collision_util_set_natural_scale_factors !! Scales dimenional quantities to ~O(1) with respect to the collisional system.
procedure :: set_original_scale => collision_util_set_original_scale_factors !! Restores dimenional quantities back to the original system units
end type collision_basic


Expand Down Expand Up @@ -380,26 +386,23 @@ module subroutine collision_util_add_fragments_to_collider(self, nbody_system, p
class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters
end subroutine collision_util_add_fragments_to_collider

module subroutine collision_util_get_angular_momentum(self)
implicit none
class(collision_fragments(*)), intent(inout) :: self !! Fragment system object
end subroutine collision_util_get_angular_momentum

module subroutine collision_util_get_energy(self)
module subroutine collision_util_construct_after_system(collider, nbody_system, param, after_system)
!! Author: David A. Minton
!!
!! Constructs a temporary internal system consisting of active bodies and additional fragments. This internal temporary system is used to calculate system energy with and without fragments
implicit none
class(collision_fragments(*)), intent(inout) :: self !! Fragment system object
end subroutine collision_util_get_energy
! Arguments
class(collision_basic), intent(inout) :: collider !! Collision system object
class(base_nbody_system), intent(in) :: nbody_system !! Original swiftest nbody system object
class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters
class(base_nbody_system), allocatable, intent(out) :: after_system !! Output temporary swiftest nbody system object
end subroutine collision_util_construct_after_system

module subroutine collision_util_reset_fragments(self)
implicit none
class(collision_fragments(*)), intent(inout) :: self
end subroutine collision_util_reset_fragments

module subroutine collision_util_set_budgets(self)
implicit none
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
class(collision_basic), intent(inout) :: self !! collisional system
Expand Down Expand Up @@ -444,14 +447,14 @@ module subroutine collision_util_get_idvalues_snapshot(self, idvals)
integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values saved in this snapshot
end subroutine collision_util_get_idvalues_snapshot

module subroutine collision_util_get_energy_momentum(self, nbody_system, param, lbefore)
module subroutine collision_util_get_energy_and_momentum(self, nbody_system, param, phase)
use base, only : base_nbody_system, base_parameters
implicit none
class(collision_basic), intent(inout) :: self !! Encounter collision system object
class(collision_basic), intent(inout) :: self !! Encounter collision system object
class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
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
end subroutine collision_util_get_energy_momentum
character(len=*), intent(in) :: phase !! One of "before" or "after", indicating which phase of the calculation this needs to be done
end subroutine collision_util_get_energy_and_momentum

module subroutine collision_util_index_map(self)
implicit none
Expand All @@ -476,6 +479,17 @@ module subroutine collision_util_snapshot(self, param, nbody_system, t, arg)
real(DP), intent(in), optional :: t !! Time of snapshot if different from nbody_system time
character(*), intent(in), optional :: arg !! "before": takes a snapshot just before the collision. "after" takes the snapshot just after the collision.
end subroutine collision_util_snapshot

module subroutine collision_util_set_natural_scale_factors(self)
implicit none
class(collision_basic), intent(inout) :: self !! collision system object
end subroutine collision_util_set_natural_scale_factors

module subroutine collision_util_set_original_scale_factors(self)
implicit none
class(collision_basic), intent(inout) :: self !! collision system object
end subroutine collision_util_set_original_scale_factors

end interface

contains
Expand Down
31 changes: 22 additions & 9 deletions src/collision/collision_resolve.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -312,13 +312,13 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status)
implicit none
! Arguments
class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(base_parameters), intent(inout) :: param !! Current run configuration parameters with Swiftest additions
real(DP), intent(in) :: t !! Time of collision
integer(I4B), intent(in) :: status !! Status flag to assign to adds
class(base_parameters), intent(inout) :: param !! Current run configuration parameters with Swiftest additions
real(DP), intent(in) :: t !! Time of collision
integer(I4B), intent(in) :: status !! Status flag to assign to adds
! Internals
integer(I4B) :: i, ibiggest, ismallest, iother, nimpactors, nfrag
logical, dimension(:), allocatable :: lmask
class(swiftest_pl), allocatable :: plnew, plsub
logical, dimension(:), allocatable :: lmask
class(swiftest_pl), allocatable :: plnew, plsub
character(*), parameter :: FRAGFMT = '("Newbody",I0.7)'
character(len=NAMELEN) :: newname, origin_type
real(DP) :: volume
Expand Down Expand Up @@ -488,7 +488,8 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec)
real(DP), intent(in) :: dt !! Current simulation step size
integer(I4B), intent(in) :: irec !! Current recursion level
! Internals
real(DP) :: E_orbit_before, E_orbit_after
real(DP) :: E_before, E_after, dLmag
real(DP), dimension(NDIM) :: L_orbit_before, L_orbit_after, L_spin_before, L_spin_after, L_before, L_after, dL_orbit, dL_spin, dL
logical :: lplpl_collision
character(len=STRMAX) :: timestr
integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision
Expand All @@ -515,7 +516,10 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec)
! Get the energy before the collision is resolved
if (param%lenergy) then
call nbody_system%get_energy_and_momentum(param)
E_orbit_before = nbody_system%te
E_before = nbody_system%te
L_orbit_before(:) = nbody_system%L_orbit(:)
L_spin_before(:) = nbody_system%L_spin(:)
L_before(:) = nbody_system%L_total(:)
end if

do loop = 1, MAXCASCADE
Expand Down Expand Up @@ -584,8 +588,17 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec)

if (param%lenergy) then
call nbody_system%get_energy_and_momentum(param)
E_orbit_after = nbody_system%te
nbody_system%E_collisions = nbody_system%E_collisions + (E_orbit_after - E_orbit_before)
E_after = nbody_system%te
nbody_system%E_collisions = nbody_system%E_collisions + (E_after - E_before)

L_orbit_after(:) = nbody_system%L_orbit(:)
L_spin_after(:) = nbody_system%L_spin(:)
L_after(:) = nbody_system%L_total(:)

dL_orbit = L_orbit_after(:) - L_orbit_before(:)
dL_spin = L_spin_after(:) - L_spin_before(:)
dL = L_after(:) - L_before(:)
dLmag = .mag.dL
end if

end associate
Expand Down
Loading

0 comments on commit 36ec14c

Please sign in to comment.