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

Commit

Permalink
More restructuring and cleaning (not there yet)
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Dec 16, 2022
1 parent e031a2e commit 9e4dd98
Show file tree
Hide file tree
Showing 4 changed files with 77 additions and 65 deletions.
60 changes: 60 additions & 0 deletions src/collision/collision_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -324,6 +324,66 @@ module subroutine collision_util_reset_system(self)
end subroutine collision_util_reset_system




module subroutine collision_util_set_coordinate_system(self)
!! author: David A. Minton
!!
!! Defines the collisional coordinate system, including the unit vectors of both the system and individual fragments.
implicit none
! Arguments
class(collision_system), intent(inout) :: self !! Collisional system
! Internals
integer(I4B) :: i
real(DP), dimension(NDIM) :: delta_r, delta_v, Ltot
real(DP) :: L_mag
real(DP), dimension(NDIM, self%fragments%nbody) :: L_sigma

associate(fragments => self%fragments, impactors => self%impactors, nfrag => self%fragments%nbody)
delta_v(:) = impactors%vb(:, 2) - impactors%vb(:, 1)
delta_r(:) = impactors%rb(:, 2) - impactors%rb(:, 1)

! We will initialize fragments on a plane defined by the pre-impact system, with the z-axis aligned with the angular momentum vector
! and the y-axis aligned with the pre-impact distance vector.

! y-axis is the separation distance
fragments%y_coll_unit(:) = .unit.delta_r(:)
Ltot = impactors%Lorbit(:,1) + impactors%Lorbit(:,2) + impactors%Lspin(:,1) + impactors%Lspin(:,2)

L_mag = .mag.Ltot(:)
if (L_mag > sqrt(tiny(L_mag))) then
fragments%z_coll_unit(:) = .unit.Ltot(:)
else ! Not enough angular momentum to determine a z-axis direction. We'll just pick a random direction
call random_number(fragments%z_coll_unit(:))
fragments%z_coll_unit(:) = .unit.fragments%z_coll_unit(:)
end if

! The cross product of the y- by z-axis will give us the x-axis
fragments%x_coll_unit(:) = fragments%y_coll_unit(:) .cross. fragments%z_coll_unit(:)

fragments%v_coll_unit(:) = .unit.delta_v(:)

if (.not.any(fragments%r_coll(:,:) > 0.0_DP)) return
fragments%rmag(:) = .mag. fragments%r_coll(:,:)

! Randomize the tangential velocity direction.
! This helps to ensure that the tangential velocity doesn't completely line up with the angular momentum vector, otherwise we can get an ill-conditioned system
call random_number(L_sigma(:,:))
do concurrent(i = 1:nfrag, fragments%rmag(i) > 0.0_DP)
fragments%v_n_unit(:, i) = fragments%z_coll_unit(:) + 2e-1_DP * (L_sigma(:,i) - 0.5_DP)
end do

! Define the radial, normal, and tangential unit vectors for each individual fragment
fragments%v_r_unit(:,:) = .unit. fragments%r_coll(:,:)
fragments%v_n_unit(:,:) = .unit. fragments%v_n_unit(:,:)
fragments%v_t_unit(:,:) = .unit. (fragments%v_n_unit(:,:) .cross. fragments%v_r_unit(:,:))

end associate

return
end subroutine collision_util_set_coordinate_system


subroutine collision_util_save_snapshot(collision_history, snapshot)
!! author: David A. Minton
!!
Expand Down
7 changes: 5 additions & 2 deletions src/fraggle/fraggle_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,9 @@ module subroutine fraggle_generate_fragments(self, system, param, lfailure)
fpe_quiet_modes(:) = .false.
call ieee_set_halting_mode(IEEE_ALL,fpe_quiet_modes)

associate(fragments => self, nfrag => self%nbody, pl => system%pl)
select type(fragments => self%fragments)
class is (fraggle_fragments)
associate(collision => self, impactors => self%impactors, nfrag => fragments%nbody, pl => system%pl)

write(message,*) nfrag
call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle generating " // trim(adjustl(message)) // " fragments.")
Expand All @@ -72,7 +74,7 @@ module subroutine fraggle_generate_fragments(self, system, param, lfailure)
call fragments%reset()

! Calculate the initial energy of the system without the collisional family
call fragments%get_energy_and_momentum(impactors, system, param, lbefore=.true.)
call collision%get_energy_and_momentum(system, param, lbefore=.true.)

! Start out the fragments close to the initial separation distance. This will be increased if there is any overlap or we fail to find a solution
r_max_start = 1.2_DP * .mag.(impactors%rb(:,2) - impactors%rb(:,1))
Expand Down Expand Up @@ -162,6 +164,7 @@ module subroutine fraggle_generate_fragments(self, system, param, lfailure)
! Restore the big array
if (lk_plpl) call pl%flatten(param)
end associate
end select
call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily

return
Expand Down
58 changes: 0 additions & 58 deletions src/fraggle/fraggle_set.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -160,64 +160,6 @@ module subroutine fraggle_set_mass_dist_fragments(self, impactors, param)
end subroutine fraggle_set_mass_dist_fragments


module subroutine encounter_set_coordinate_system(self, impactors)
!! author: David A. Minton
!!
!! Defines the collisional coordinate system, including the unit vectors of both the system and individual fragments.
implicit none
! Arguments
class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object
class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object
! Internals
integer(I4B) :: i
real(DP), dimension(NDIM) :: delta_r, delta_v, Ltot
real(DP) :: L_mag
real(DP), dimension(NDIM, self%nbody) :: L_sigma

associate(fragments => self, nfrag => self%nbody)
delta_v(:) = impactors%vb(:, 2) - impactors%vb(:, 1)
delta_r(:) = impactors%rb(:, 2) - impactors%rb(:, 1)

! We will initialize fragments on a plane defined by the pre-impact system, with the z-axis aligned with the angular momentum vector
! and the y-axis aligned with the pre-impact distance vector.

! y-axis is the separation distance
fragments%y_coll_unit(:) = .unit.delta_r(:)
Ltot = impactors%Lorbit(:,1) + impactors%Lorbit(:,2) + impactors%Lspin(:,1) + impactors%Lspin(:,2)

L_mag = .mag.Ltot(:)
if (L_mag > sqrt(tiny(L_mag))) then
fragments%z_coll_unit(:) = .unit.Ltot(:)
else ! Not enough angular momentum to determine a z-axis direction. We'll just pick a random direction
call random_number(fragments%z_coll_unit(:))
fragments%z_coll_unit(:) = .unit.fragments%z_coll_unit(:)
end if

! The cross product of the y- by z-axis will give us the x-axis
fragments%x_coll_unit(:) = fragments%y_coll_unit(:) .cross. fragments%z_coll_unit(:)

fragments%v_coll_unit(:) = .unit.delta_v(:)

if (.not.any(fragments%r_coll(:,:) > 0.0_DP)) return
fragments%rmag(:) = .mag. fragments%r_coll(:,:)

! Randomize the tangential velocity direction.
! This helps to ensure that the tangential velocity doesn't completely line up with the angular momentum vector, otherwise we can get an ill-conditioned system
call random_number(L_sigma(:,:))
do concurrent(i = 1:nfrag, fragments%rmag(i) > 0.0_DP)
fragments%v_n_unit(:, i) = fragments%z_coll_unit(:) + 2e-1_DP * (L_sigma(:,i) - 0.5_DP)
end do

! Define the radial, normal, and tangential unit vectors for each individual fragment
fragments%v_r_unit(:,:) = .unit. fragments%r_coll(:,:)
fragments%v_n_unit(:,:) = .unit. fragments%v_n_unit(:,:)
fragments%v_t_unit(:,:) = .unit. (fragments%v_n_unit(:,:) .cross. fragments%v_r_unit(:,:))

end associate

return
end subroutine encounter_set_coordinate_system


module subroutine fraggle_set_natural_scale_factors(self, impactors)
!! author: David A. Minton
Expand Down
17 changes: 12 additions & 5 deletions src/modules/collision_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,11 +102,12 @@ module collision_classes
real(DP), dimension(2) :: pe !! Before/after potential energy
real(DP), dimension(2) :: Etot !! Before/after total system energy
contains
procedure :: regime => collision_regime_system !! Determine which fragmentation regime the set of impactors will be
procedure :: setup => collision_setup_system !! Initializer for the encounter collision system. Allocates the collider and fragments classes and the before/after snapshots
procedure :: get_energy_and_momentum => collision_util_get_energy_momentum !! Calculates total 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
final :: collision_util_final_system !! Finalizer will deallocate all allocatables
procedure :: regime => collision_regime_system !! Determine which fragmentation regime the set of impactors will be
procedure :: setup => collision_setup_system !! Initializer for the encounter collision system. Allocates the collider and fragments classes and the before/after snapshots
procedure :: get_energy_and_momentum => collision_util_get_energy_momentum !! Calculates total 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_coordinate_system => collision_util_set_coordinate_system !! Sets the coordinate system of the collisional system
final :: collision_util_final_system !! Finalizer will deallocate all allocatables
end type collision_system

!! NetCDF dimension and variable names for the enounter save object
Expand Down Expand Up @@ -217,6 +218,12 @@ module subroutine collision_set_coordinate_fragments(self)
class(collision_fragments), intent(inout) :: self !! Fragment system object
end subroutine collision_set_coordinate_fragments

module subroutine collision_util_set_coordinate_system(self)
implicit none
class(collision_system), intent(inout) :: self !! Collisional system
end subroutine collision_util_set_coordinate_system


module subroutine collision_setup_fragments(self, n, param)
implicit none
class(collision_fragments), intent(inout) :: self !! Fragment system object
Expand Down

0 comments on commit 9e4dd98

Please sign in to comment.