Fixed memory leaks. The main one was due to opening already-open NetC…

…DF files when writing and reading headers

src/encounter/encounter_check.f90

@@ -402,8 +402,6 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, lvdotr,
          npl_last = npl
       end if
 
-      call timer%reset()
-      call timer%start()
       !$omp parallel do default(private) schedule(static) &
       !$omp shared(x, v, renc, boundingbox) &
       !$omp firstprivate(dt, npl, n)
@@ -419,31 +417,17 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, lvdotr,
                                                x(dim,1:npl) + renc(1:npl) + vshift_max(1:npl) * v(dim,1:npl) * dt])
       end do
       !$omp end parallel do 
-      call timer%stop()
-      write(*,*) "plpl  sort  : ",timer%count_stop_step - timer%count_start_step
 
-      call timer%reset()
-      call timer%start()
       call boundingbox%sweep(npl, nenc, index1, index2)
-      call timer%stop()
-      write(*,*) "plpl  sweep : ",timer%count_stop_step - timer%count_start_step
 
       if (nenc > 0) then
          ! Now that we have identified potential pairs, use the narrow-phase process to get the final values
          allocate(lencounter(nenc))
          allocate(lvdotr(nenc))
 
-         call timer%reset()
-         call timer%start()
          call encounter_check_all(nenc, index1, index2, x, v, x, v, renc, renc, dt, lencounter, lvdotr)
-         call timer%stop()
-         write(*,*) "plpl  check : ",timer%count_stop_step - timer%count_start_step
 
-         call timer%reset()
-         call timer%start()
          call encounter_check_reduce_broadphase(npl, nenc, index1, index2, lencounter, lvdotr)
-         call timer%stop()
-         write(*,*) "plpl  reduce: ",timer%count_stop_step - timer%count_start_step
          deallocate(lencounter)
       end if
 
@@ -520,14 +504,8 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt
                                                 xplt(dim,1:nplt) + renct(1:nplt) + vpltshift_max(1:nplt) * vplt(dim,1:nplt) * dt])
       end do
       !$omp end parallel do
-      call timer%stop()
-      write(*,*) "plplm sort  : ",timer%count_stop_step - timer%count_start_step
-
-      call timer%reset()
-      call timer%start()
+
       call boundingbox%sweep(nplm, nplt, nenc, index1, index2)
-      call timer%stop()
-      write(*,*) "plplm sweep : ",timer%count_stop_step - timer%count_start_step
 
       if (nenc > 0) then
          ! Shift tiny body indices back into the range of the input position and velocity arrays
@@ -537,15 +515,9 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt
          allocate(lencounter(nenc))
          allocate(lvdotr(nenc))
 
-         call timer%reset()
-         call timer%start()
          call encounter_check_all(nenc, index1, index2, xplm, vplm, xplt, vplt, rencm, renct, dt, lencounter, lvdotr)
-         call timer%stop()
-         write(*,*) "plplm check : ",timer%count_stop_step - timer%count_start_step
 
          call encounter_check_reduce_broadphase(ntot, nenc, index1, index2, lencounter, lvdotr)
-         call timer%stop()
-         write(*,*) "plplm reduce: ",timer%count_stop_step - timer%count_start_step
       end if
       return
    end subroutine encounter_check_all_sort_and_sweep_plplm
@@ -706,8 +678,6 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, lvdotr, inde
 
       call util_index_array(ind_arr, npl) 
 
-      call timer%reset()
-      call timer%start()
       !$omp parallel do default(private) schedule(static)&
       !$omp shared(x, v, renc, lenc, ind_arr) &
       !$omp firstprivate(npl, dt)
@@ -719,8 +689,6 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, lvdotr, inde
                                                          renc(i), renc(:), dt, ind_arr(:), lenc(i))
       end do
       !$omp end parallel do
-      call timer%stop()
-      write(*,*) "plpl  triang: ",timer%count_stop_step - timer%count_start_step
 
       call encounter_check_collapse_ragged_list(lenc, npl, nenc, index1, index2, lvdotr)
 
@@ -758,8 +726,6 @@ subroutine encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vp
 
       call util_index_array(ind_arr, nplt)
 
-      call timer%reset()
-      call timer%start()
       !$omp parallel do default(private) schedule(static)&
       !$omp shared(xplm, vplm, xplt, vplt, rencm, renct, lenc, ind_arr) &
       !$omp firstprivate(nplm, nplt, dt)
@@ -771,8 +737,6 @@ subroutine encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vp
                                                           rencm(i), renct(:), dt, ind_arr(:), lenc(i))
       end do
       !$omp end parallel do
-      call timer%stop()
-      write(*,*) "plplm triang: ",timer%count_stop_step - timer%count_start_step
 
       call encounter_check_collapse_ragged_list(lenc, nplm, nenc, index1, index2, lvdotr)
 
@@ -991,7 +955,6 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, nenc, ind
       type(encounter_list), dimension(n1+n2) :: lenc         !! Array of encounter lists (one encounter list per body)
       integer(I4B), dimension(:), allocatable, save :: ind_arr
       integer(I4B), dimension(:), allocatable :: ibeg, iend
-      type(walltimer) :: timer 
 
       ntot = n1 + n2
       call util_index_array(ind_arr, ntot)
@@ -1005,12 +968,8 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, nenc, ind
       end do
       ! Sweep the intervals for each of the massive bodies along one dimension
       ! This will build a ragged pair of index lists inside of the lenc data structure
-      call timer%reset()
-      call timer%start()
       call encounter_check_sweep_aabb_all_double_list(n1, n2, self%aabb(1)%ind(:), reshape(ibeg(:), [SWEEPDIM, ntot]), reshape(iend(:), [SWEEPDIM, ntot]), ind_arr(:), lenc(:))
-      call timer%stop()
-      write(*,*) "sweep double: ",timer%count_stop_step - timer%count_start_step
-
+
       call encounter_check_collapse_ragged_list(lenc, ntot, nenc, index1, index2)
 
       ! Reorder the pairs and sort the first index in order to remove any duplicates

src/main/swiftest_driver.f90

@@ -123,6 +123,8 @@ program swiftest_driver
       end do
    end associate
 
+   call nbody_system%dealloc()
+
    call util_exit(SUCCESS)
 
    stop

src/netcdf/netcdf.f90

@@ -650,8 +650,6 @@ module subroutine netcdf_read_hdr_system(self, iu, param)
 
       tslot = int(param%ioutput, kind=I4B) + 1
 
-      call check( nf90_open(param%outfile, NF90_NOWRITE, iu%ncid) )
-
       call check( nf90_get_var(iu%ncid, iu%time_varid, param%t,       start=[tslot]) )
       call check( nf90_get_var(iu%ncid, iu%npl_varid,  self%pl%nbody, start=[tslot]) )
       call check( nf90_get_var(iu%ncid, iu%ntp_varid,  self%tp%nbody, start=[tslot]) )
@@ -1117,8 +1115,6 @@ module subroutine netcdf_write_hdr_system(self, iu, param)
 
       tslot = int(param%ioutput, kind=I4B) + 1
 
-      call check( nf90_open(param%outfile, nf90_write, iu%ncid) )
-
       call check( nf90_put_var(iu%ncid, iu%time_varid, param%t, start=[tslot]) )
       call check( nf90_put_var(iu%ncid, iu%npl_varid, self%pl%nbody, start=[tslot]) )
       call check( nf90_put_var(iu%ncid, iu%ntp_varid, self%tp%nbody, start=[tslot]) )

src/symba/symba_util.f90

@@ -251,6 +251,8 @@ module subroutine symba_util_dealloc_pl(self)
       implicit none
       ! Arguments
       class(symba_pl),  intent(inout) :: self !! SyMBA massive body object
+      ! Internals
+      integer(I4B) :: i
 
       if (allocated(self%lcollision)) deallocate(self%lcollision)
       if (allocated(self%lencounter)) deallocate(self%lencounter)
@@ -262,6 +264,11 @@ module subroutine symba_util_dealloc_pl(self)
       if (allocated(self%isperi)) deallocate(self%isperi)
       if (allocated(self%peri)) deallocate(self%peri)
       if (allocated(self%atp)) deallocate(self%atp)
+
+      do i = 1, self%nbody
+         call self%kin(i)%dealloc()
+      end do
+
       if (allocated(self%kin)) deallocate(self%kin)
 
       call util_dealloc_pl(self)
@@ -408,7 +415,7 @@ module subroutine symba_util_flatten_eucl_plpl(self, param)
          if (param%lflatten_interactions) then
             if (allocated(self%k_plpl)) deallocate(self%k_plpl) ! Reset the index array if it's been set previously
             allocate(self%k_plpl(2, nplpl), stat=err)
-            if (err /=0) then ! An error occurred trying to allocate this big array. This probably means it's too big to fit in memory, and so we will force the run back into triangular mode
+            if (err /= 0) then ! An error occurred trying to allocate this big array. This probably means it's too big to fit in memory, and so we will force the run back into triangular mode
                param%lflatten_interactions = .false.
             else
                do concurrent (i=1:npl, j=1:npl, j>i)