Major restructuring of encounter storage classes and methods

examples/Fragmentation/Fragmentation_Movie.py

@@ -41,8 +41,8 @@
 movie_titles = dict(zip(available_movie_styles, movie_title_list))
 
 # These initial conditions were generated by trial and error
-pos_vectors = {"disruption_headon"         : [np.array([1.0, -5.0e-04, 0.0]),
-                                              np.array([1.0,  5.0e-04 ,0.0])],
+pos_vectors = {"disruption_headon"         : [np.array([1.0, -5.0e-05, 0.0]),
+                                              np.array([1.0,  5.0e-05 ,0.0])],
                "supercatastrophic_off_axis": [np.array([1.0, -4.2e-05,      0.0]),
                                               np.array([1.0,  4.2e-05,      0.0])],
                "hitandrun"                 : [np.array([1.0, -2.0e-05,      0.0]),
@@ -203,7 +203,7 @@ def data_stream(self, frame=0):
         minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body
         gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades
         sim.set_parameter(fragmentation=True, fragmentation_save="TRAJECTORY", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False)
-        sim.run(dt=1e-4, tstop=1.0e-4, istep_out=1, dump_cadence=0)
+        sim.run(dt=1e-4, tstop=1.0e-3, istep_out=1, dump_cadence=0)
 
         print("Generating animation")
         anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1)

src/CMakeLists.txt

@@ -30,6 +30,7 @@ SET(FAST_MATH_FILES
             ${SRC}/encounter/encounter_check.f90
             ${SRC}/encounter/encounter_setup.f90
             ${SRC}/encounter/encounter_util.f90
+            ${SRC}/encounter/encounter_io.f90
             ${SRC}/fraggle/fraggle_generate.f90
             ${SRC}/fraggle/fraggle_io.f90
             ${SRC}/fraggle/fraggle_placeholder.f90
@@ -74,6 +75,7 @@ SET(FAST_MATH_FILES
             ${SRC}/util/util_dealloc.f90
             ${SRC}/util/util_exit.f90
             ${SRC}/util/util_fill.f90
+            ${SRC}/util/util_final.f90
             ${SRC}/util/util_flatten.f90
             ${SRC}/util/util_get_energy_momentum.f90
             ${SRC}/util/util_index_array.f90

src/encounter/encounter_io.f90

@@ -0,0 +1,207 @@
+!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh
+!! This file is part of Swiftest.
+!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License 
+!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
+!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty 
+!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
+!! You should have received a copy of the GNU General Public License along with Swiftest. 
+!! If not, see: https://www.gnu.org/licenses. 
+
+submodule (encounter_classes) s_encounter_io
+   use swiftest
+contains
+
+
+   module subroutine encounter_io_dump(self, param)
+      !! author: David A. Minton
+      !!
+      !! Dumps the time history of an encounter to file.
+      implicit none
+      ! Arguments
+      class(encounter_storage(*)),  intent(inout)        :: self   !! Encounter storage object
+      class(swiftest_parameters),   intent(inout)        :: param  !! Current run configuration parameters 
+      ! Internals
+      integer(I4B) :: i
+
+      ! Most of this is just temporary test code just to get something working. Eventually this should get cleaned up.
+
+      do i = 1, self%nframes
+         if (allocated(self%frame(i)%item)) then
+            select type(snapshot => self%frame(i)%item)
+            class is (encounter_snapshot)
+               param%ioutput = self%tslot(i)
+               call snapshot%write_frame(self%nc,param)
+            end select
+         else
+            exit
+         end if
+      end do
+
+
+      return
+   end subroutine encounter_io_dump
+
+
+   module subroutine encounter_io_initialize(self, param)
+      !! author: David A. Minton
+      !!
+      !! Initialize a NetCDF encounter file system. This is a simplified version of the main simulation output NetCDF file, but with fewer variables.
+      use, intrinsic :: ieee_arithmetic
+      use netcdf
+      implicit none
+      ! Arguments
+      class(encounter_io_parameters), intent(inout) :: self    !! Parameters used to identify a particular NetCDF dataset
+      class(swiftest_parameters),     intent(in)    :: param   !! Current run configuration parameters
+      ! Internals
+      integer(I4B) :: nvar, varid, vartype
+      real(DP) :: dfill
+      real(SP) :: sfill
+      logical :: fileExists
+      character(len=STRMAX) :: errmsg
+      integer(I4B) :: ndims
+
+      associate(nc => self)
+         dfill = ieee_value(dfill, IEEE_QUIET_NAN)
+         sfill = ieee_value(sfill, IEEE_QUIET_NAN)
+
+         select case (param%out_type)
+         case("NETCDF_FLOAT")
+            self%out_type = NF90_FLOAT
+         case("NETCDF_DOUBLE")
+            self%out_type = NF90_DOUBLE
+         end select
+
+         ! Check if the file exists, and if it does, delete it
+         inquire(file=nc%enc_file, exist=fileExists)
+         if (fileExists) then
+            open(unit=LUN, file=nc%enc_file, status="old", err=667, iomsg=errmsg)
+            close(unit=LUN, status="delete")
+         end if
+
+         call check( nf90_create(nc%enc_file, NF90_NETCDF4, nc%id), "encounter_io_initialize nf90_create" )
+
+         ! Dimensions
+         call check( nf90_def_dim(nc%id, nc%time_dimname, nc%time_dimsize, nc%time_dimid), "encounter_io_initialize nf90_def_dim time_dimid" ) ! Simulation time dimension
+         call check( nf90_def_dim(nc%id, nc%space_dimname, NDIM , nc%space_dimid), "encounter_io_initialize nf90_def_dim space_dimid" )           ! 3D space dimension
+         call check( nf90_def_dim(nc%id, nc%id_dimname, param%maxid, nc%id_dimid), "encounter_io_initialize nf90_def_dim id_dimid" )       ! dimension to store particle id numbers
+         call check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "encounter_io_initialize nf90_def_dim str_dimid"  )          ! Dimension for string variables (aka character arrays)
+
+         ! Dimension coordinates
+         call check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, nc%time_dimid, nc%time_varid), "encounter_io_initialize nf90_def_var time_varid"  )
+         call check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "encounter_io_initialize nf90_def_var space_varid"  )
+         call check( nf90_def_var(nc%id, nc%id_dimname, NF90_INT, nc%id_dimid, nc%id_varid), "encounter_io_initialize nf90_def_var id_varid"  )
+
+         ! Variables
+         call check( nf90_def_var(nc%id, nc%name_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], nc%name_varid), "encounter_io_initialize nf90_def_var name_varid"  )
+         call check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], nc%ptype_varid), "encounter_io_initialize nf90_def_var ptype_varid"  )
+         call check( nf90_def_var(nc%id, nc%rh_varname,  nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%rh_varid), "encounter_io_initialize nf90_def_var rh_varid"  )
+         call check( nf90_def_var(nc%id, nc%vh_varname,  nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%vh_varid), "encounter_io_initialize nf90_def_var vh_varid"  )
+         call check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%Gmass_varid), "encounter_io_initialize nf90_def_var Gmass_varid"  )
+         call check( nf90_def_var(nc%id, nc%level_varname, NF90_INT, [nc%id_dimid, nc%time_dimid], nc%level_varid), "encounter_io_initialize nf90_def_var level_varid"  )
+         if (param%lclose) then
+            call check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%radius_varid), "encounter_io_initialize nf90_def_var radius_varid"  )
+         end if
+         if (param%lrotation) then
+            call check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%Ip_varid), "encounter_io_initialize nf90_def_var Ip_varid"  )
+            call check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%rot_varid), "encounter_io_initialize nf90_def_var rot_varid"  )
+         end if
+
+         call check( nf90_inquire(nc%id, nVariables=nvar), "encounter_io_initialize nf90_inquire nVariables"  )
+         do varid = 1, nvar
+            call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "encounter_io_initialize nf90_inquire_variable"  )
+            select case(vartype)
+            case(NF90_INT)
+               call check( nf90_def_var_fill(nc%id, varid, 0, NF90_FILL_INT), "encounter_io_initialize nf90_def_var_fill NF90_INT"  )
+            case(NF90_FLOAT)
+               call check( nf90_def_var_fill(nc%id, varid, 0, sfill), "encounter_io_initialize nf90_def_var_fill NF90_FLOAT"  )
+            case(NF90_DOUBLE)
+               call check( nf90_def_var_fill(nc%id, varid, 0, dfill), "encounter_io_initialize nf90_def_var_fill NF90_DOUBLE"  )
+            case(NF90_CHAR)
+               call check( nf90_def_var_fill(nc%id, varid, 0, 0), "encounter_io_initialize nf90_def_var_fill NF90_CHAR"  )
+            end select
+         end do
+
+         ! Take the file out of define mode
+         call check( nf90_enddef(nc%id), "encounter_io_initialize nf90_enddef"  )
+
+         ! Add in the space dimension coordinates
+         call check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "encounter_io_initialize nf90_put_var space"  )
+      end associate
+
+      return
+
+      667 continue
+      write(*,*) "Error creating encounter output file. " // trim(adjustl(errmsg))
+      call util_exit(FAILURE)
+   end subroutine encounter_io_initialize
+
+
+   module subroutine encounter_io_write_frame(self, nc, param)
+      !! author: David A. Minton
+      !!
+      !! Write a frame of output of an encounter trajectory.
+      use netcdf
+      implicit none
+      ! Arguments
+      class(encounter_snapshot),      intent(in)    :: self  !! Swiftest encounter structure
+      class(encounter_io_parameters), intent(inout) :: nc    !! Parameters used to identify a particular encounter io NetCDF dataset
+      class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
+      ! Internals
+      integer(I4B)           :: i, tslot, idslot, old_mode, npl, ntp
+      character(len=NAMELEN) :: charstring
+
+      tslot = param%ioutput
+
+      call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "encounter_io_write_frame nf90_set_fill"  )
+
+      call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "encounter_io_write_frame nf90_put_var time_varid"  )
+
+      associate(pl => self%pl, tp => self%tp)
+         npl = pl%nbody
+         do i = 1, npl
+            idslot = pl%id(i)
+            call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[idslot]), "encounter_io_write_frame nf90_put_var pl id_varid"  )
+            call check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl rh_varid"  )
+            call check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl vh_varid"  )
+            call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl Gmass_varid"  )
+            select type(pl)
+            class is (symba_pl)
+               call check( nf90_put_var(nc%id, nc%level_varid, pl%levelg(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl level_varid"  )
+            end select
+
+            if (param%lclose) call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl radius_varid"  )
+
+            if (param%lrotation) then
+               call check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl Ip_varid"  )
+               call check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl rotx_varid"  )
+            end if
+
+            charstring = trim(adjustl(pl%info(i)%name))
+            call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "encounter_io_write_frame nf90_put_var pl name_varid"  )
+            charstring = trim(adjustl(pl%info(i)%particle_type))
+            call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "encounter_io_write_frame nf90_put_var pl particle_type_varid"  )
+         end do
+
+         ntp = tp%nbody
+         do i = 1, ntp
+            idslot = tp%id(i)
+            call check( nf90_put_var(nc%id, nc%id_varid, tp%id(i), start=[idslot]), "encounter_io_write_frame nf90_put_var tp id_varid"  )
+            call check( nf90_put_var(nc%id, nc%rh_varid, tp%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var tp rh_varid"  )
+            call check( nf90_put_var(nc%id, nc%vh_varid, tp%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var tp vh_varid"  )
+
+            charstring = trim(adjustl(tp%info(i)%name))
+            call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "encounter_io_write_frame nf90_put_var tp name_varid"  )
+            charstring = trim(adjustl(tp%info(i)%particle_type))
+            call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "encounter_io_write_frame nf90_put_var tp particle_type_varid"  )
+         end do
+      end associate
+
+      call check( nf90_set_fill(nc%id, old_mode, old_mode) )
+
+      return
+   end subroutine encounter_io_write_frame
+
+
+
+
+end submodule s_encounter_io

src/encounter/encounter_util.f90

@@ -136,6 +136,36 @@ module subroutine encounter_util_final_list(self)
    end subroutine encounter_util_final_list
 
 
+   module subroutine encounter_util_final_snapshot(self)
+      !! author: David A. Minton
+      !!
+      !! Deallocates allocatable arrays in an encounter snapshot
+      implicit none
+      ! Arguments
+      type(encounter_snapshot),  intent(inout) :: self !! Encounter storage object
+
+      if (allocated(self%pl)) deallocate(self%pl)
+      if (allocated(self%tp)) deallocate(self%tp)
+      self%t = 0.0_DP
+
+      return
+   end subroutine encounter_util_final_snapshot
+
+
+   module subroutine encounter_util_final_storage(self)
+      !! author: David A. Minton
+      !!
+      !! Deallocates allocatable arrays in an encounter snapshot
+      implicit none
+      ! Arguments
+      type(encounter_storage(*)),  intent(inout) :: self !! Encounter storage object
+
+      call util_final_storage(self%swiftest_storage)
+
+      return
+   end subroutine encounter_util_final_storage
+
+
    module subroutine encounter_util_resize_list(self, nnew)
       !! author: David A. Minton
       !!

src/fraggle/fraggle_io.f90

@@ -12,6 +12,26 @@
 
 contains
 
+
+   module subroutine fraggle_io_encounter_dump(self, param)
+      implicit none
+      class(fraggle_storage(*)), intent(inout) :: self   !! Encounter storage object
+      class(swiftest_parameters),          intent(inout) :: param  !! Current run configuration parameters 
+   end subroutine fraggle_io_encounter_dump
+
+   module subroutine fraggle_io_encounter_initialize_output(self, param)
+      implicit none
+      class(fraggle_io_encounter_parameters), intent(inout) :: self    !! Parameters used to identify a particular NetCDF dataset
+      class(swiftest_parameters),             intent(in)    :: param   
+   end subroutine fraggle_io_encounter_initialize_output
+
+   module subroutine fraggle_io_encounter_write_frame(self, nc, param)
+      implicit none
+      class(fraggle_encounter_snapshot), intent(in)    :: self   !! Swiftest encounter structure
+      class(encounter_io_parameters),    intent(inout) :: nc     !! Parameters used to identify a particular encounter io NetCDF dataset
+      class(swiftest_parameters),        intent(inout) :: param  !! Current run configuration parameters
+   end subroutine fraggle_io_encounter_write_frame
+
    module subroutine fraggle_io_log_generate(frag)
       !! author: David A. Minton
       !!

src/fraggle/fraggle_util.f90

@@ -131,6 +131,70 @@ module subroutine fraggle_util_construct_temporary_system(frag, system, param, t
    end subroutine fraggle_util_construct_temporary_system
 
 
+
+   module subroutine fraggle_util_final_colliders(self)
+      !! author: David A. Minton
+      !!
+      !! Finalizer will deallocate all allocatables
+      implicit none
+      ! Arguments
+      type(fraggle_colliders),  intent(inout) :: self !! Fraggle encountar storage object
+
+      if (allocated(self%idx)) deallocate(self%idx)
+
+      return
+   end subroutine fraggle_util_final_colliders
+
+
+   module subroutine fraggle_util_final_fragments(self)
+      !! author: David A. Minton
+      !!
+      !! Finalizer will deallocate all allocatables
+      implicit none
+      ! Arguments
+      type(fraggle_fragments),  intent(inout) :: self !! Fraggle encountar storage object
+
+      if (allocated(self%x_coll)) deallocate(self%x_coll)
+      if (allocated(self%v_coll)) deallocate(self%v_coll)
+      if (allocated(self%v_r_unit)) deallocate(self%v_r_unit)
+      if (allocated(self%v_t_unit)) deallocate(self%v_t_unit)
+      if (allocated(self%v_n_unit)) deallocate(self%v_n_unit)
+      if (allocated(self%rmag)) deallocate(self%rmag)
+      if (allocated(self%rotmag)) deallocate(self%rotmag)
+      if (allocated(self%v_r_mag)) deallocate(self%v_r_mag)
+      if (allocated(self%v_t_mag)) deallocate(self%v_t_mag)
+
+      return
+   end subroutine fraggle_util_final_fragments
+
+
+   module subroutine fraggle_util_final_storage(self)
+      !! author: David A. Minton
+      !!
+      !! Finalizer will deallocate all allocatables
+      implicit none
+      ! Arguments
+      type(fraggle_storage(*)),  intent(inout) :: self !! Fraggle encountar storage object
+
+      call util_final_storage(self%swiftest_storage)
+
+      return
+   end subroutine fraggle_util_final_storage
+
+   module subroutine fraggle_util_final_snapshot(self)
+      !! author: David A. Minton
+      !!
+      !! Finalizer will deallocate all allocatables
+      implicit none
+      ! Arguments
+      type(fraggle_encounter_snapshot),  intent(inout) :: self !! Fraggle encountar storage object
+
+      call encounter_util_final_snapshot(self%encounter_snapshot)
+
+      return
+   end subroutine fraggle_util_final_snapshot
+
+
    module subroutine fraggle_util_get_energy_momentum(self, colliders, system, param, lbefore)
       !! Author: David A. Minton
       !!

src/helio/helio_util.f90

@@ -33,7 +33,7 @@ module subroutine helio_util_final_system(self)
       ! Arguments
       type(helio_nbody_system),  intent(inout) :: self !! Helio nbody system object
 
-      call self%dealloc()
+      call whm_util_final_system(self%whm_nbody_system)
 
       return
    end subroutine helio_util_final_system

src/io/io.f90

@@ -122,9 +122,7 @@ module subroutine io_conservation_report(self, param, lterminal)
       ! Internals
       real(DP), dimension(NDIM)       :: Ltot_now,  Lorbit_now,  Lspin_now
       real(DP)                        :: ke_orbit_now,  ke_spin_now,  pe_now,  Eorbit_now
-      real(DP)                        :: Eorbit_error, Etot_error, Ecoll_error
       real(DP)                        :: GMtot_now
-      real(DP)                        :: Lerror, Merror
       character(len=STRMAX)           :: errmsg
       integer(I4B), parameter         :: EGYIU = 72
       character(len=*), parameter     :: EGYTERMFMT = '(" DL/L0 = ", ES12.5 &