From 6926760bf6547376f90c6fcfe85a3c15c4fd76a5 Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 7 May 2021 18:02:02 -0400 Subject: [PATCH] Fixed some memory issues, and moving rearray back out from loop because of indexing problem --- src/symba/symba_collision.f90 | 25 ++++++++++++------------- src/util/util_resize_pl.f90 | 14 +++++++++++++- 2 files changed, 25 insertions(+), 14 deletions(-) diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 8f43603c8..2c2fccd02 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -45,12 +45,11 @@ subroutine symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, nmer type(family_array), dimension(2) :: parent_child_index_array logical, dimension(nplplenc) :: lplpl_collision logical, dimension(:), allocatable :: lplpl_unique_parent - integer(I4B), dimension(:), pointer :: plparent logical :: ldiscard ! First determine the collisional regime for each colliding pair associate(npl => symba_plA%helio%swiftest%nbody, xbpl => symba_plA%helio%swiftest%xb, statpl => symba_plA%helio%swiftest%status, idpl => symba_plA%helio%swiftest%id, & - idx1 => plplenc_list%index1, idx2 => plplenc_list%index2, plparent => symba_plA%kin%parent) + idx1 => plplenc_list%index1, idx2 => plplenc_list%index2) lplpl_collision(:) = plplenc_list%status(1:nplplenc) == COLLISION ldiscard = any(lplpl_collision) if (.not.ldiscard) return @@ -63,7 +62,7 @@ subroutine symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, nmer ! Get the subset of collisions that involve a unique pair of parents allocate(lplpl_unique_parent(ncollisions)) - lplpl_unique_parent(:) = plparent(idx1(collision_idx(:))) /= plparent(idx2(collision_idx(:))) + lplpl_unique_parent(:) = symba_plA%kin(idx1(collision_idx(:)))%parent /= symba_plA%kin(idx2(collision_idx(:)))%parent nunique_parent = count(lplpl_unique_parent(:)) allocate(unique_parent_idx(nunique_parent)) unique_parent_idx = pack(collision_idx(:), lplpl_unique_parent(:)) @@ -75,12 +74,12 @@ subroutine symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, nmer lplpl_unique_parent(:) = .true. do index_coll = 1, ncollisions index_enc = collision_idx(index_coll) - idx(1) = plparent(idx1(index_enc)) - idx(2) = plparent(idx2(index_enc)) - lplpl_unique_parent(:) = .not. ( any(plparent(idx1(unique_parent_idx(:))) == idx(1)) .or. & - any(plparent(idx2(unique_parent_idx(:))) == idx(1)) .or. & - any(plparent(idx1(unique_parent_idx(:))) == idx(2)) .or. & - any(plparent(idx2(unique_parent_idx(:))) == idx(2)) ) + idx(1) = symba_plA%kin(idx1(index_enc))%parent + idx(2) = symba_plA%kin(idx2(index_enc))%parent + lplpl_unique_parent(:) = .not. ( any(symba_plA%kin(idx1(unique_parent_idx(:)))%parent == idx(1)) .or. & + any(symba_plA%kin(idx2(unique_parent_idx(:)))%parent == idx(1)) .or. & + any(symba_plA%kin(idx1(unique_parent_idx(:)))%parent == idx(2)) .or. & + any(symba_plA%kin(idx2(unique_parent_idx(:)))%parent == idx(2)) ) end do ! Reassemble collision index list to include only those containing the unique pairs of parents, plus all the non-unique pairs that don't @@ -100,7 +99,7 @@ subroutine symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, nmer idx(2) = idx2(index_enc) ! Index values for the parents of this particle pair - idx_parent(:) = plparent(idx(:)) + idx_parent(:) = symba_plA%kin(idx(:))%parent if (any(statpl(idx_parent(:)) /= ACTIVE)) cycle ! One of these two bodies is already gone @@ -123,8 +122,8 @@ subroutine symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, nmer ! Group together the ids and indexes of each collisional parent and its children do j = 1, 2 - allocate(parent_child_index_array(j)%idx(nchild(j)+ 1)) - allocate(parent_child_index_array(j)%id(nchild(j)+ 1)) + allocate(parent_child_index_array(j)%idx(nchild(j) + 1)) + allocate(parent_child_index_array(j)%id(nchild(j) + 1)) associate(idx_arr => parent_child_index_array(j)%idx, & id_arr => parent_child_index_array(j)%id, & ncj => nchild(j), & @@ -304,9 +303,9 @@ subroutine symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, nmer end do deallocate(family) - call symba_rearray_pl(t, symba_plA, nmergeadd, mergeadd_list, discard_plA, param) end do + call symba_rearray_pl(t, symba_plA, nmergeadd, mergeadd_list, discard_plA, param) end associate return diff --git a/src/util/util_resize_pl.f90 b/src/util/util_resize_pl.f90 index 039552c67..2fd8441c3 100644 --- a/src/util/util_resize_pl.f90 +++ b/src/util/util_resize_pl.f90 @@ -18,6 +18,7 @@ subroutine util_resize_pl(symba_plA, npl_new) ! Internals type(symba_pl) :: new_symba_plA + integer(I4B) :: i ! Executable code @@ -35,6 +36,11 @@ subroutine util_resize_pl(symba_plA, npl_new) new_symba_plA%helio%swiftest%vb(:,1:npl_old) = symba_plA%helio%swiftest%vb(:,1:npl_old) new_symba_plA%helio%swiftest%rot(:,1:npl_old) = symba_plA%helio%swiftest%rot(:,1:npl_old) new_symba_plA%helio%swiftest%Ip(:,1:npl_old) = symba_plA%helio%swiftest%Ip(:,1:npl_old) + new_symba_plA%helio%swiftest%info(1:npl_old) = symba_plA%helio%swiftest%info(1:npl_old) + new_symba_plA%kin(1:npl_old) = symba_plA%kin(1:npl_old) + do i = 1, npl_old + if (allocated(symba_plA%kin(i)%child)) allocate(new_symba_plA%kin(i)%child, source = symba_plA%kin(i)%child) + end do else new_symba_plA%helio%swiftest%id(1:npl_new) = symba_plA%helio%swiftest%id(1:npl_new) new_symba_plA%helio%swiftest%status(1:npl_new) = symba_plA%helio%swiftest%status(1:npl_new) @@ -47,6 +53,11 @@ subroutine util_resize_pl(symba_plA, npl_new) new_symba_plA%helio%swiftest%vb(:,1:npl_new) = symba_plA%helio%swiftest%vb(:,1:npl_new) new_symba_plA%helio%swiftest%rot(:,1:npl_new) = symba_plA%helio%swiftest%rot(:,1:npl_new) new_symba_plA%helio%swiftest%Ip(:,1:npl_new) = symba_plA%helio%swiftest%Ip(:,1:npl_new) + new_symba_plA%helio%swiftest%info(1:npl_new) = symba_plA%helio%swiftest%info(1:npl_new) + new_symba_plA%kin(1:npl_new) = symba_plA%kin(1:npl_new) + do i = 1, npl_new + if (allocated(symba_plA%kin(i)%child)) allocate(new_symba_plA%kin(i)%child, source = symba_plA%kin(i)%child) + end do end if call symba_pl_deallocate(symba_plA) call symba_pl_allocate(symba_plA, npl_new) @@ -61,7 +72,8 @@ subroutine util_resize_pl(symba_plA, npl_new) call move_alloc(new_symba_plA%helio%swiftest%vb, symba_plA%helio%swiftest%vb) call move_alloc(new_symba_plA%helio%swiftest%rot, symba_plA%helio%swiftest%rot) call move_alloc(new_symba_plA%helio%swiftest%Ip, symba_plA%helio%swiftest%Ip) - call symba_pl_deallocate(new_symba_plA) + call move_alloc(new_symba_plA%helio%swiftest%info, symba_plA%helio%swiftest%info) + call move_alloc(new_symba_plA%kin, symba_plA%kin) end associate return