Restructured code to promote encounter subroutines to its own module.…

… Refactored encounter list type

Makefile

@@ -48,6 +48,7 @@ SWIFTEST_MODULES =   swiftest_globals.f90 \
                      swiftest_operators.f90 \
 							lambda_function.f90\
                      swiftest_classes.f90 \
+							encounter_classes.f90 \
                      fraggle_classes.f90 \
                      whm_classes.f90 \
                      rmvs_classes.f90 \

src/encounter/encounter_io.f90

@@ -0,0 +1,90 @@
+submodule (encounter_classes) s_encounter_io
+   use swiftest
+contains
+
+
+   module subroutine encounter_io_write_frame(iu, t, id1, id2, Gmass1, Gmass2, radius1, radius2, xh1, xh2, vh1, vh2)
+      !! author: David A. Minton
+      !!
+      !! Write a single frame of close encounter data to output binary files
+      !!
+      !! Adapted from David E. Kaufmann's Swifter routine: io_write_encounter.f90
+      !! Adapted from Hal Levison's Swift routine io_write_encounter.f
+      implicit none
+      ! Arguments
+      integer(I4B),           intent(in) :: iu               !! Open file unit number
+      real(DP),               intent(in) :: t                !! Time of encounter
+      integer(I4B),           intent(in) :: id1, id2         !! ids of the two encountering bodies
+      real(DP),               intent(in) :: Gmass1, Gmass2   !! G*mass of the two encountering bodies
+      real(DP),               intent(in) :: radius1, radius2 !! Radii of the two encountering bodies
+      real(DP), dimension(:), intent(in) :: xh1, xh2         !! Heliocentric position vectors of the two encountering bodies 
+      real(DP), dimension(:), intent(in) :: vh1, vh2         !! Heliocentric velocity vectors of the two encountering bodies  
+      ! Internals
+      character(len=STRMAX)   :: errmsg
+
+      write(iu, err = 667, iomsg = errmsg) t
+      write(iu, err = 667, iomsg = errmsg) id1, xh1(1), xh1(2), xh1(3), vh1(1), vh1(2), Gmass1, radius1
+      write(iu, err = 667, iomsg = errmsg) id2, xh2(1), xh2(2), xh2(3), vh2(1), vh2(2), Gmass2, radius2
+
+      return
+      667 continue
+      write(*,*) "Error writing encounter file: " // trim(adjustl(errmsg))
+      call util_exit(FAILURE)
+   end subroutine
+
+   module subroutine encounter_io_write_list(self, pl, encbody, param)
+      implicit none
+      ! Arguments
+      class(encounter_list),  intent(in) :: self    !! Swiftest encounter list object
+      class(swiftest_pl),         intent(in) :: pl      !! Swiftest massive body object
+      class(swiftest_body),       intent(in) :: encbody !! Encountering body - Swiftest generic body object (pl or tp) 
+      class(swiftest_parameters), intent(in) :: param   !! Current run configuration parameters 
+      ! Internals
+      logical , save          :: lfirst = .true.
+      integer(I4B)            :: k, ierr
+      character(len=STRMAX)   :: errmsg
+
+      if (param%enc_out == "" .or. self%nenc == 0) return
+
+      open(unit = LUN, file = param%enc_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr, iomsg = errmsg)
+      if (ierr /= 0) then
+         if (lfirst) then
+            open(unit = LUN, file = param%enc_out, status = 'NEW', form = 'UNFORMATTED', err = 667, iomsg = errmsg)
+         else
+            goto 667
+         end if
+      end if
+      lfirst = .false.
+
+      associate(ind1 => self%index1, ind2 => self%index2)
+         select type(encbody)
+         class is (swiftest_pl)
+            do k = 1, self%nenc
+               call encounter_io_write_frame(LUN, self%t(k), &
+                                             pl%id(ind1(k)),     encbody%id(ind2(k)), &
+                                             pl%Gmass(ind1(k)),  encbody%Gmass(ind2(k)), &
+                                             pl%radius(ind1(k)), encbody%radius(ind2(k)), &
+                                             self%x1(:,k),       self%x2(:,k), &
+                                             self%v1(:,k),       self%v2(:,k))
+            end do
+         class is (swiftest_tp)
+            do k = 1, self%nenc
+               call encounter_io_write_frame(LUN, self%t(k), &
+                                             pl%id(ind1(k)),     encbody%id(ind2(k)), &
+                                             pl%Gmass(ind1(k)),  0.0_DP, &
+                                             pl%radius(ind1(k)), 0.0_DP, &
+                                             self%x1(:,k),       self%x2(:,k), &
+                                             self%v1(:,k),       self%v2(:,k))
+            end do 
+         end select
+      end associate
+
+      close(unit = LUN, err = 667, iomsg = errmsg)
+
+      return
+      667 continue
+      write(*,*) "Error writing encounter file: " // trim(adjustl(errmsg))
+      call util_exit(FAILURE)
+   end subroutine encounter_io_write_list
+
+end submodule s_encounter_io

src/encounter/encounter_setup.f90

@@ -0,0 +1,61 @@
+submodule (encounter_classes) s_encounter_setup
+   use swiftest
+contains
+
+   module subroutine encounter_setup_list(self, n)
+      !! author: David A. Minton
+      !!
+      !! A constructor that sets the number of encounters and allocates and initializes all arrays  
+      !!
+      implicit none
+      ! Arguments
+      class(encounter_list), intent(inout) :: self !! Swiftest encounter structure
+      integer(I4B),              intent(in)    :: n    !! Number of encounters to allocate space for
+
+      if (n < 0) return
+
+      if (allocated(self%lvdotr)) deallocate(self%lvdotr)
+      if (allocated(self%status)) deallocate(self%status)
+      if (allocated(self%index1)) deallocate(self%index1)
+      if (allocated(self%index2)) deallocate(self%index2)
+      if (allocated(self%id1)) deallocate(self%id1)
+      if (allocated(self%id2)) deallocate(self%id2)
+      if (allocated(self%x1)) deallocate(self%x1)
+      if (allocated(self%x2)) deallocate(self%x2)
+      if (allocated(self%v1)) deallocate(self%v1)
+      if (allocated(self%v2)) deallocate(self%v2)
+      if (allocated(self%t)) deallocate(self%t)
+
+      self%nenc = n
+      if (n == 0) return
+
+      allocate(self%lvdotr(n))
+      allocate(self%status(n))
+      allocate(self%index1(n))
+      allocate(self%index2(n))
+      allocate(self%id1(n))
+      allocate(self%id2(n))
+      allocate(self%x1(NDIM,n))
+      allocate(self%x2(NDIM,n))
+      allocate(self%v1(NDIM,n))
+      allocate(self%v2(NDIM,n))
+      allocate(self%t(n))
+
+      self%lvdotr(:) = .false.
+      self%status(:) = INACTIVE
+      self%index1(:) = 0
+      self%index2(:) = 0
+      self%id1(:) = 0
+      self%id2(:) = 0
+      self%x1(:,:) = 0.0_DP
+      self%x2(:,:) = 0.0_DP
+      self%v1(:,:) = 0.0_DP
+      self%v2(:,:) = 0.0_DP
+      self%t(:) = 0.0_DP
+
+      return
+   end subroutine encounter_setup_list
+
+end submodule s_encounter_setup
+
+

src/encounter/encounter_util.f90

@@ -0,0 +1,141 @@
+submodule (encounter_classes) s_encounter_util
+   use swiftest
+contains
+
+   module subroutine encounter_util_append_list(self, source, lsource_mask)
+      !! author: David A. Minton
+      !!
+      !! Append components from one Swiftest body object to another. 
+      !! This method will automatically resize the destination body if it is too small
+      implicit none
+      ! Arguments
+      class(encounter_list), intent(inout) :: self         !! Swiftest encounter list object
+      class(encounter_list), intent(in)    :: source       !! Source object to append
+      logical, dimension(:),     intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+      ! Internals
+      integer(I4B) :: nold, nsrc
+
+      nold = self%nenc
+      nsrc = source%nenc
+      call util_append(self%lvdotr, source%lvdotr, nold, nsrc, lsource_mask)
+      call util_append(self%status, source%status, nold, nsrc, lsource_mask)
+      call util_append(self%index1, source%index1, nold, nsrc, lsource_mask)
+      call util_append(self%index2, source%index2, nold, nsrc, lsource_mask)
+      call util_append(self%id1, source%id1, nold, nsrc, lsource_mask)
+      call util_append(self%id2, source%id2, nold, nsrc, lsource_mask)
+      call util_append(self%x1, source%x1, nold, nsrc, lsource_mask)
+      call util_append(self%x2, source%x2, nold, nsrc, lsource_mask)
+      call util_append(self%v1, source%v1, nold, nsrc, lsource_mask)
+      call util_append(self%v2, source%v2, nold, nsrc, lsource_mask)
+      call util_append(self%t, source%t, nold, nsrc, lsource_mask)
+      self%nenc = nold + count(lsource_mask(1:nsrc))
+
+      return
+   end subroutine encounter_util_append_list
+
+
+   module subroutine encounter_util_copy_list(self, source)
+      !! author: David A. Minton
+      !!
+      !! Copies elements from the source encounter list into self.
+      implicit none
+      ! Arguments
+      class(encounter_list), intent(inout) :: self   !! Encounter list 
+      class(encounter_list), intent(in)    :: source !! Source object to copy into
+
+      associate(n => source%nenc)
+         self%nenc = n
+         self%lvdotr(1:n) = source%lvdotr(1:n) 
+         self%status(1:n) = source%status(1:n) 
+         self%index1(1:n) = source%index1(1:n)
+         self%index2(1:n) = source%index2(1:n)
+         self%id1(1:n) = source%id1(1:n)
+         self%id2(1:n) = source%id2(1:n)
+         self%x1(:,1:n) = source%x1(:,1:n)
+         self%x2(:,1:n) = source%x2(:,1:n)
+         self%v1(:,1:n) = source%v1(:,1:n)
+         self%v2(:,1:n) = source%v2(:,1:n)
+         self%t(1:n) = source%t(1:n)
+      end associate
+
+      return
+   end subroutine encounter_util_copy_list
+
+
+   module subroutine encounter_util_resize_list(self, nnew)
+      !! author: David A. Minton
+      !!
+      !! Checks the current size of the encounter list 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_list), intent(inout) :: self !! Swiftest encounter list 
+      integer(I4B),              intent(in)    :: nnew !! New size of list needed
+      ! Internals
+      class(encounter_list), allocatable :: enc_temp
+      integer(I4B)                           :: nold
+      logical                                :: lmalloc
+
+      lmalloc = allocated(self%status)
+      if (lmalloc) then
+         nold = size(self%status)
+      else
+         nold = 0
+      end if
+      if (nnew > nold) then
+         if (lmalloc) allocate(enc_temp, source=self)
+         call self%setup(2 * nnew)
+         if (lmalloc) then
+            call self%copy(enc_temp)
+            deallocate(enc_temp)
+         end if
+      else
+         self%status(nnew+1:nold) = INACTIVE
+      end if
+      self%nenc = nnew
+
+      return
+   end subroutine encounter_util_resize_list
+
+
+   module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestructive)
+      !! author: David A. Minton
+      !!
+      !! Move spilled (discarded) Swiftest encounter structure from active list to discard list
+      implicit none
+      ! Arguments
+      class(encounter_list), intent(inout) :: self         !! Swiftest encounter list 
+      class(encounter_list), intent(inout) :: discards     !! Discarded object 
+      logical, dimension(:),     intent(in)    :: lspill_list  !! Logical array of bodies to spill into the discards
+      logical,                   intent(in)    :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list
+      ! Internals
+      integer(I4B) :: nenc_old
+
+      associate(keeps => self)
+         call util_spill(keeps%lvdotr, discards%lvdotr, lspill_list, ldestructive)
+         call util_spill(keeps%status, discards%status, lspill_list, ldestructive)
+         call util_spill(keeps%index1, discards%index1, lspill_list, ldestructive)
+         call util_spill(keeps%index2, discards%index2, lspill_list, ldestructive)
+         call util_spill(keeps%id1, discards%id1, lspill_list, ldestructive)
+         call util_spill(keeps%id2, discards%id2, lspill_list, ldestructive)
+         call util_spill(keeps%x1, discards%x1, lspill_list, ldestructive)
+         call util_spill(keeps%x2, discards%x2, lspill_list, ldestructive)
+         call util_spill(keeps%v1, discards%v1, lspill_list, ldestructive)
+         call util_spill(keeps%v2, discards%v2, lspill_list, ldestructive)
+         call util_spill(keeps%t, discards%t, lspill_list, ldestructive)
+
+         nenc_old = keeps%nenc
+
+         ! This is the base class, so will be the last to be called in the cascade. 
+         ! Therefore we need to set the nenc values for both the keeps and discareds
+         discards%nenc = count(lspill_list(1:nenc_old))
+         if (ldestructive) keeps%nenc = nenc_old - discards%nenc
+      end associate
+
+      return
+   end subroutine encounter_util_spill_list
+
+
+end submodule s_encounter_util