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

Commit

Permalink
Added start of merger resolving code
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Aug 4, 2021
1 parent 8c8c8a0 commit 96e3170
Show file tree
Hide file tree
Showing 2 changed files with 44 additions and 21 deletions.
4 changes: 2 additions & 2 deletions src/modules/symba_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -283,7 +283,7 @@ end subroutine symba_io_dump_particle_info

module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
implicit none
class(symba_parameters), intent(inout) :: self !! Collection of SyMBA parameters
class(symba_parameters), intent(inout) :: self !! Current run configuration parameters with SyMBA additionss
integer, intent(in) :: unit !! File unit number
character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT.
!! If you do not include a char-literal-constant, the iotype argument contains only DT.
Expand All @@ -294,7 +294,7 @@ end subroutine symba_io_param_reader

module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, iomsg)
implicit none
class(symba_parameters),intent(in) :: self !! Collection of SyMBA parameters
class(symba_parameters),intent(in) :: self !! Current run configuration parameters with SyMBA additions
integer, intent(in) :: unit !! File unit number
character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT.
!! If you do not include a char-literal-constant, the iotype argument contains only DT.
Expand Down
61 changes: 42 additions & 19 deletions src/symba/symba_collision.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,14 +143,15 @@ pure elemental function symba_collision_check_one(xr, yr, zr, vxr, vyr, vzr, Gmt
end function symba_collision_check_one


function symba_collision_consolidate_familes(pl, idx_parent, family, x, v, mass, radius, L_spin, Ip) result(lflag)
function symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v, mass, radius, L_spin, Ip) result(lflag)
!! author: David A. Minton
!!
!! Loops through the pl-pl collision list and groups families together by index. Outputs the indices of all family members,
!! and pairs of quantities (x and v vectors, mass, radius, L_spin, and Ip) that can be used to resolve the collisional outcome.
implicit none
! Arguments
class(symba_pl), intent(inout) :: pl !! SyMBA massive body object
class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions
integer(I4B), dimension(2), intent(inout) :: idx_parent !! Index of the two bodies considered the "parents" of the collision
integer(I4B), dimension(:), allocatable, intent(out) :: family !! List of indices of all bodies inovlved in the collision
real(DP), dimension(NDIM,2), intent(out) :: x, v, L_spin, Ip !! Output values that represent a 2-body equivalent of a possibly 2+ body collision
Expand Down Expand Up @@ -209,17 +210,17 @@ function symba_collision_consolidate_familes(pl, idx_parent, family, x, v, mass,
family = [parent_child_index_array(1)%idx(:),parent_child_index_array(2)%idx(:)]
fam_size = count(pl%status(family(:)) == ACTIVE)
family = pack(family(:), pl%status(family(:)) == ACTIVE)
L_spin(:,:) = 0.0_DP
Ip(:,:) = 0.0_DP

! Find the barycenter of each body along with its children, if it has any
do j = 1, 2
x(:, j) = pl%xb(:, idx_parent(j))
v(:, j) = pl%vb(:, idx_parent(j))
Ip(:, j) = mass(j) * pl%Ip(:, idx_parent(j))
! Assume principal axis rotation about axis corresponding to highest moment of inertia (3rd Ip)
L_spin(:, j) = Ip(3, j) * radius(j)**2 * pl%rot(:, idx_parent(j))
if (param%lrotation) L_spin(:, j) = Ip(3, j) * radius(j)**2 * pl%rot(:, idx_parent(j))

! Use conservation of energy and angular momentum to adjust velocities and distances to be those equivalent to where they would be when the mutual radii are touching
!call vector_adjust(mass(j), x(:,j), v(:,j))
if (nchild(j) > 0) then
do i = 1, nchild(j) ! Loop over all children and take the mass weighted mean of the properties
idx_child = parent_child_index_array(j)%idx(i + 1)
Expand All @@ -230,31 +231,33 @@ function symba_collision_consolidate_familes(pl, idx_parent, family, x, v, mass,
volchild = (4.0_DP / 3.0_DP) * PI * pl%radius(idx_child)**3
volume(j) = volume(j) + volchild
! Get angular momentum of the child-parent pair and add that to the spin
xcom(:) = (mass(j) * x(:,j) + mchild * xchild(:)) / (mass(j) + mchild)
vcom(:) = (mass(j) * v(:,j) + mchild * vchild(:)) / (mass(j) + mchild)
xc(:) = x(:, j) - xcom(:)
vc(:) = v(:, j) - vcom(:)
xcrossv(:) = xc(:) .cross. vc(:)
L_spin(:, j) = L_spin(:, j) + mass(j) * xcrossv(:)

xc(:) = xchild(:) - xcom(:)
vc(:) = vchild(:) - vcom(:)
xcrossv(:) = xc(:) .cross. vc(:)
L_spin(:, j) = L_spin(:, j) + mchild * xcrossv(:)

! Add the child's spin
L_spin(:, j) = L_spin(:, j) + mchild * pl%Ip(3, idx_child) * pl%radius(idx_child)**2 * pl%rot(:, idx_child)
if (param%lrotation) then
xcom(:) = (mass(j) * x(:,j) + mchild * xchild(:)) / (mass(j) + mchild)
vcom(:) = (mass(j) * v(:,j) + mchild * vchild(:)) / (mass(j) + mchild)
xc(:) = x(:, j) - xcom(:)
vc(:) = v(:, j) - vcom(:)
xcrossv(:) = xc(:) .cross. vc(:)
L_spin(:, j) = L_spin(:, j) + mass(j) * xcrossv(:)

xc(:) = xchild(:) - xcom(:)
vc(:) = vchild(:) - vcom(:)
xcrossv(:) = xc(:) .cross. vc(:)
L_spin(:, j) = L_spin(:, j) + mchild * xcrossv(:)

L_spin(:, j) = L_spin(:, j) + mchild * pl%Ip(3, idx_child) * pl%radius(idx_child)**2 * pl%rot(:, idx_child)
Ip(:, j) = Ip(:, j) + mchild * pl%Ip(:, idx_child)
end if

! Merge the child and parent
mass(j) = mass(j) + mchild
x(:, j) = xcom(:)
v(:, j) = vcom(:)
Ip(:, j) = Ip(:, j) + mchild * pl%Ip(:, idx_child)
end do
end if
density(j) = mass(j) / volume(j)
radius(j) = ((3 * mass(j)) / (density(j) * 4 * pi))**(1.0_DP / 3.0_DP)
Ip(:, j) = Ip(:, j) / mass(j)
if (param%lrotation) Ip(:, j) = Ip(:, j) / mass(j)
end do
lflag = .true.

Expand Down Expand Up @@ -418,6 +421,26 @@ module subroutine symba_collision_resolve_mergers(self, system, param)
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions
! Internals
integer(I4B) :: i
logical :: lgoodcollision
integer(I4B), dimension(:), allocatable :: family !! List of indices of all bodies inovlved in the collision
integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision
real(DP), dimension(NDIM,2) :: x, v, L_spin, Ip !! Output values that represent a 2-body equivalent of a possibly 2+ body collision
real(DP), dimension(2) :: mass, radius !! Output values that represent a 2-body equivalent of a possibly 2+ body collision

associate(plpl_collisions => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2)
select type(pl => system%pl)
class is (symba_pl)
do i = 1, ncollisions
idx_parent(1) = pl%kin(idx1(i))%parent
idx_parent(2) = pl%kin(idx2(i))%parent
lgoodcollision = symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v, mass, radius, L_spin, Ip)
if (.not. lgoodcollision) cycle
if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved
!call symba_collision_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip)
end do
end select
end associate

return
end subroutine symba_collision_resolve_mergers
Expand Down

0 comments on commit 96e3170

Please sign in to comment.