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

Commit

Permalink
More cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Dec 16, 2022
1 parent 9e4dd98 commit 265483b
Show file tree
Hide file tree
Showing 3 changed files with 30 additions and 30 deletions.
20 changes: 10 additions & 10 deletions src/collision/collision_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -347,34 +347,34 @@ module subroutine collision_util_set_coordinate_system(self)
! and the y-axis aligned with the pre-impact distance vector.

! y-axis is the separation distance
fragments%y_coll_unit(:) = .unit.delta_r(:)
impactors%y_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(:)
impactors%z_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(:)
call random_number(impactors%z_unit(:))
impactors%z_unit(:) = .unit.impactors%z_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(:)
impactors%x_unit(:) = impactors%y_unit(:) .cross. impactors%z_unit(:)

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

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

! 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)
fragments%v_n_unit(:, i) = impactors%z_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_r_unit(:,:) = .unit. fragments%rc(:,:)
fragments%v_n_unit(:,:) = .unit. fragments%v_n_unit(:,:)
fragments%v_t_unit(:,:) = .unit. (fragments%v_n_unit(:,:) .cross. fragments%v_r_unit(:,:))

Expand Down
32 changes: 16 additions & 16 deletions src/fraggle/fraggle_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,35 +213,35 @@ subroutine fraggle_generate_pos_vec(fragments, impactors, r_max_start)

! We will treat the first two fragments of the list as special cases. They get initialized the maximum distances apart along the original impactor distance vector.
! This is done because in a regular disruption, the first body is the largest, the second the second largest, and the rest are smaller equal-mass fragments.
call random_number(fragments%r_coll(:,3:nfrag))
call random_number(fragments%rc(:,3:nfrag))
loverlap(:) = .true.
do while (any(loverlap(3:nfrag)))
fragments%r_coll(:, 1) = impactors%rb(:, 1) - fragments%rbcom(:)
fragments%r_coll(:, 2) = impactors%rb(:, 2) - fragments%rbcom(:)
fragments%rc(:, 1) = impactors%rb(:, 1) - fragments%rbcom(:)
fragments%rc(:, 2) = impactors%rb(:, 2) - fragments%rbcom(:)
r_max = r_max + 0.1_DP * rad
do i = 3, nfrag
if (loverlap(i)) then
call random_number(fragments%r_coll(:,i))
fragments%r_coll(:,i) = 2 * (fragments%r_coll(:, i) - 0.5_DP)
fragments%r_coll(:, i) = fragments%r_coll(:,i) + fdistort * vunit(:)
fragments%r_coll(:, i) = r_max * fragments%r_coll(:, i)
fragments%r_coll(:, i) = fragments%r_coll(:, i) + (fragments%rbimp(:) - fragments%rbcom(:)) ! Shift the center of the fragment cloud to the impact point rather than the CoM
!if (lnoncat .and. dot_product(fragments%r_coll(:,i), runit(:)) < 0.0_DP) fragments%r_coll(:, i) = -fragments%r_coll(:, i) ! Make sure the fragment cloud points away from the impact point
call random_number(fragments%rc(:,i))
fragments%rc(:,i) = 2 * (fragments%rc(:, i) - 0.5_DP)
fragments%rc(:, i) = fragments%rc(:,i) + fdistort * vunit(:)
fragments%rc(:, i) = r_max * fragments%rc(:, i)
fragments%rc(:, i) = fragments%rc(:, i) + (fragments%rbimp(:) - fragments%rbcom(:)) ! Shift the center of the fragment cloud to the impact point rather than the CoM
!if (lnoncat .and. dot_product(fragments%rc(:,i), runit(:)) < 0.0_DP) fragments%rc(:, i) = -fragments%rc(:, i) ! Make sure the fragment cloud points away from the impact point
end if
end do
loverlap(:) = .false.
do j = 1, nfrag
do i = j + 1, nfrag
dis = .mag.(fragments%r_coll(:,j) - fragments%r_coll(:,i))
dis = .mag.(fragments%rc(:,j) - fragments%rc(:,i))
loverlap(i) = loverlap(i) .or. (dis <= (fragments%radius(i) + fragments%radius(j)))
end do
end do
end do
call fraggle_util_shift_vector_to_origin(fragments%mass, fragments%r_coll)
call fraggle_util_shift_vector_to_origin(fragments%mass, fragments%rc)
call fragments%set_coordinate_system(impactors)

do concurrent(i = 1:nfrag)
fragments%rb(:,i) = fragments%r_coll(:,i) + fragments%rbcom(:)
fragments%rb(:,i) = fragments%rc(:,i) + fragments%rbcom(:)
end do

fragments%rbcom(:) = 0.0_DP
Expand Down Expand Up @@ -408,7 +408,7 @@ subroutine fraggle_generate_tan_vel(fragments, impactors, lfailure)
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), fragments%vbcom(:))
do concurrent (i = 1:nfrag)
fragments%v_coll(:,i) = fragments%vb(:,i) - fragments%vbcom(:)
fragments%vc(:,i) = fragments%vb(:,i) - fragments%vbcom(:)
end do

! Now do a kinetic energy budget check to make sure we are still within the budget.
Expand All @@ -422,7 +422,7 @@ subroutine fraggle_generate_tan_vel(fragments, impactors, lfailure)
ke_diff = fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit
lfailure = ke_diff < 0.0_DP
if (.not.lfailure) exit
fragments%r_coll(:,:) = fragments%r_coll(:,:) * 1.1_DP
fragments%rc(:,:) = fragments%rc(:,:) * 1.1_DP
end do
if (lfailure) then
call io_log_one_message(FRAGGLE_LOG_OUT, " ")
Expand Down Expand Up @@ -476,7 +476,7 @@ function solve_fragment_tan_vel(lfailure, v_t_mag_input) result(v_t_mag_output)
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%r_coll(:, i) .cross. vtmp(:))
L(:) = fragments%mass(i) * (fragments%rc(:, i) .cross. vtmp(:))
L_orb_others(:) = L_orb_others(:) + L(:)
end if
end do
Expand Down Expand Up @@ -582,7 +582,7 @@ subroutine fraggle_generate_rad_vel(fragments, impactors, lfailure)
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), fragments%vbcom(:))
do i = 1, nfrag
fragments%v_coll(:, i) = fragments%vb(:, i) - fragments%vbcom(:)
fragments%vc(:, i) = fragments%vb(:, i) - fragments%vbcom(:)
fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * norm2(fragments%vb(:, i))
fragments%ke_spin = fragments%ke_spin + fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * norm2(fragments%rot(:,i))
end do
Expand Down
8 changes: 4 additions & 4 deletions src/fraggle/fraggle_set.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -244,12 +244,12 @@ module subroutine fraggle_set_original_scale_factors(self)
fragments%mass = fragments%mass * fragments%mscale
fragments%radius = fragments%radius * fragments%dscale
fragments%rot = fragments%rot / fragments%tscale
fragments%r_coll = fragments%r_coll * fragments%dscale
fragments%v_coll = fragments%v_coll * fragments%vscale
fragments%rc = fragments%rc * fragments%dscale
fragments%vc = fragments%vc * fragments%vscale

do i = 1, fragments%nbody
fragments%rb(:, i) = fragments%r_coll(:, i) + fragments%rbcom(:)
fragments%vb(:, i) = fragments%v_coll(:, i) + fragments%vbcom(:)
fragments%rb(:, i) = fragments%rc(:, i) + fragments%rbcom(:)
fragments%vb(:, i) = fragments%vc(:, i) + fragments%vbcom(:)
end do

impactors%Qloss = impactors%Qloss * fragments%Escale
Expand Down

0 comments on commit 265483b

Please sign in to comment.