From 82c93299426345e493675d0b2d98324450868c30 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 2 Dec 2022 12:27:05 -0500 Subject: [PATCH] Added a resize method and more operations leading up to saving the encounter list history --- src/encounter/encounter_util.f90 | 38 +++++++++++++++++++++++++++++ src/modules/encounter_classes.f90 | 8 +++++- src/modules/symba_classes.f90 | 1 + src/symba/symba_encounter_check.f90 | 10 ++++++-- src/symba/symba_step.f90 | 1 + 5 files changed, 55 insertions(+), 3 deletions(-) diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index 76ac0e492..8cdc5c94f 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -176,6 +176,44 @@ module subroutine encounter_util_resize_list(self, nnew) end subroutine encounter_util_resize_list + module subroutine encounter_util_resize_storage(self, nnew) + !! author: David A. Minton + !! + !! Checks the current size of the encounter storage against the required size and extends it by a factor of 2 more than requested if it is too small. + !! Note: The reason to extend it by a factor of 2 is for performance. When there are many enounters per step, resizing every time you want to add an + !! encounter takes significant computational effort. Resizing by a factor of 2 is a tradeoff between performance (fewer resize calls) and memory managment + !! Memory usage grows by a factor of 2 each time it fills up, but no more. + implicit none + ! Arguments + class(encounter_storage(*)), allocatable, intent(inout) :: self !! Swiftest encounter list + integer(I4B), intent(in) :: nnew !! New size of list needed + ! Internals + type(encounter_storage(nframes=:)), allocatable :: tmp + integer(I4B) :: i, nold + logical :: lmalloc + + + lmalloc = allocated(self) + if (lmalloc) then + nold = self%nframes + else + nold = 0 + end if + + if (nnew > nold) then + allocate(encounter_storage(nnew) :: tmp) + if (lmalloc) then + do i = 1, nold + if (allocated(self%frame(i)%item)) tmp%frame(i) = self%frame(i)%item + end do + deallocate(self) + end if + call move_alloc(tmp,self) + end if + + return + end subroutine encounter_util_resize_storage + module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestructive) !! author: David A. Minton !! diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index 17bb81d9d..10c9cceec 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -43,7 +43,7 @@ module encounter_classes end type encounter_list type, extends(swiftest_storage) :: encounter_storage - !! A class that that is used to store simulation history data between file output + !! A class that that is used to store simulation history data between file output contains procedure :: dump => encounter_io_dump_storage_list end type encounter_storage @@ -279,6 +279,12 @@ module subroutine encounter_util_resize_list(self, nnew) integer(I8B), intent(in) :: nnew !! New size of list needed end subroutine encounter_util_resize_list + module subroutine encounter_util_resize_storage(self, nnew) + implicit none + class(encounter_storage(*)), allocatable, intent(inout) :: self !! Swiftest encounter list + integer(I4B), intent(in) :: nnew !! New size of list needed + end subroutine encounter_util_resize_storage + module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestructive) implicit none class(encounter_list), intent(inout) :: self !! Swiftest encounter list diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 988d75d95..f303d6143 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -186,6 +186,7 @@ module symba_classes class(symba_plplenc), allocatable :: plplcollision_list !! List of massive body-massive body collisions in a single step integer(I4B) :: irec !! System recursion level type(encounter_storage(nframes=:)), allocatable :: encounter_history + integer(I4B) :: iframe = 0 !! Encounter history frame number contains procedure :: write_discard => symba_io_write_discard !! Write out information about discarded and merged planets and test particles in SyMBA procedure :: initialize => symba_setup_initialize_system !! Performs SyMBA-specific initilization steps diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 26996302f..8983dfa97 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -34,7 +34,7 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l lany_encounter = .false. if (self%nbody == 0) return - associate(pl => self, plplenc_list => system%plplenc_list) + associate(pl => self, plplenc_list => system%plplenc_list, cb => system%cb, iframe => system%iframe, encounter_history => system%encounter_history) npl = pl%nbody nplm = pl%nplm @@ -57,7 +57,7 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l call move_alloc(index2, plplenc_list%index2) end if - if (lany_encounter) then + if (lany_encounter) then do k = 1_I8B, nenc i = plplenc_list%index1(k) j = plplenc_list%index2(k) @@ -65,6 +65,10 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l plplenc_list%id2(k) = pl%id(j) plplenc_list%status(k) = ACTIVE plplenc_list%level(k) = irec + plplenc_list%x1(:,k) = pl%xh(:,i) + plplenc_list%x2(:,k) = pl%xh(:,j) + plplenc_list%v1(:,k) = pl%vb(:,i) - cb%vb(:) + plplenc_list%v2(:,k) = pl%vb(:,j) - cb%vb(:) pl%lencounter(i) = .true. pl%lencounter(j) = .true. pl%levelg(i) = irec @@ -74,6 +78,8 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l pl%nplenc(i) = pl%nplenc(i) + 1 pl%nplenc(j) = pl%nplenc(j) + 1 end do + iframe = iframe + 1 + encounter_history%frame(iframe) = plplenc_list end if end associate diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index e24eeec31..d936ebcea 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -275,6 +275,7 @@ module subroutine symba_step_reset_system(self, param) nenc_old = system%plplenc_list%nenc call system%plplenc_list%setup(0_I8B) call system%plplcollision_list%setup(0_I8B) + system%iframe = 0 if (npl > 0) then pl%lcollision(1:npl) = .false. call pl%reset_kinship([(i, i=1, npl)])