Merge branch 'encounter_storage' into debug

examples/Fragmentation/Fragmentation_Movie.py

@@ -36,7 +36,7 @@
 # Define the names and initial conditions of the various fragmentation simulation types
 # ----------------------------------------------------------------------------------------------------------------------
 available_movie_styles = ["disruption_headon", "supercatastrophic_off_axis", "hitandrun"]
-movie_title_list = ["Head-on Disrutption", "Off-axis Supercatastrophic", "Hit and Run"]
+movie_title_list = ["Head-on Disruption", "Off-axis Supercatastrophic", "Hit and Run"]
 movie_titles = dict(zip(available_movie_styles, movie_title_list))
 
 # These initial conditions were generated by trial and error
@@ -135,8 +135,8 @@ def center(Gmass, x, y):
             y_com = np.sum(Gmass * y) / np.sum(Gmass)
             return x_com, y_com
 
-        Gmass, x, y, point_rad = next(self.data_stream(frame))
-        x_com, y_com = center(Gmass, x, y)
+        Gmass, rh, point_rad = next(self.data_stream(frame))
+        x_com, y_com = center(Gmass, rh[:,1], rh[:,2])
         self.scatter_artist.set_offsets(np.c_[x - x_com, y - y_com])
         self.scatter_artist.set_sizes(point_rad)
         return self.scatter_artist,
@@ -147,10 +147,9 @@ def data_stream(self, frame=0):
             ds = ds.where(ds['name'] != "Sun", drop=True)
             radius = ds['radius'].values
             Gmass = ds['Gmass'].values
-            x = ds['xhx'].values
-            y = ds['xhy'].values
+            rh = ds['rh'].values
             point_rad = 2 * radius * self.ax_pt_size
-            yield Gmass, x, y, point_rad
+            yield Gmass, rh, point_rad
 
 if __name__ == "__main__":
 
@@ -168,7 +167,6 @@ def data_stream(self, frame=0):
       movie_styles = available_movie_styles.copy()
 
    for style in movie_styles:
-       param_file = Path(style) / "param.in"
        bin_file = Path(style) / "bin.nc"
        if bin_file.exists():
            user_selection = input(f"An older simulation of {movie_titles[style]} exists. Do you want to re-run it? [y/N] ")
@@ -186,9 +184,9 @@ def data_stream(self, frame=0):
 
        # Pull in the Swiftest output data from the parameter file and store it as a Xarray dataset.
        if run_new:
-           sim = swiftest.Simulation(param_file=param_file, rotation=True, init_cond_format = "XV", compute_conservation_values=True)
+           sim = swiftest.Simulation(simdir=style, rotation=True, init_cond_format = "XV", compute_conservation_values=True)
            sim.add_solar_system_body("Sun")
-           sim.add_body(Gmass=body_Gmass[style], radius=body_radius[style], rh=pos_vectors[style], vh=vel_vectors[style], rot=rot_vectors[style])
+           sim.add_body(Gmass=body_Gmass[style], radius=body_radius[style], rh=pos_vectors[style], vh=vel_vectors[style]) #, rot=rot_vectors[style])
 
            # Set fragmentation parameters
            minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body

python/swiftest/swiftest/init_cond.py

@@ -284,10 +284,9 @@ def vec2xr(param: Dict, **kwargs: Any):
     kwargs = {k:kwargs[k] for k,v in kwargs.items() if v is not None}
     if param['ROTATION']:
         if "rot" not in kwargs and "Gmass" in kwargs:
-            warnings.warn("Rotation vectors must be given when rotation is enabled for massive bodies",stacklevel=2)
-            return
+            kwargs['rot'] = np.zeros((len(kwargs['Gmass']),3))
         if "Ip" not in kwargs and "rot" in kwargs:
-            kwargs['Ip'] = np.full_like(rot, 0.4)
+            kwargs['Ip'] = np.full_like(kwargs['rot'], 0.4)
 
     if "time" not in kwargs:
         kwargs["time"] = np.array([0.0])

python/swiftest/swiftest/io.py

@@ -529,12 +529,8 @@ def swifter_stream(f, param):
 
         tlab = []
         if param['OUT_FORM'] == 'XV' or param['OUT_FORM'] == 'XVEL':
-            tlab.append('rhx')
-            tlab.append('rhy')
-            tlab.append('rhz')
-            tlab.append('vhx')
-            tlab.append('vhy')
-            tlab.append('vhz')
+            tlab.append('rh')
+            tlab.append('vh')
         if param['OUT_FORM'] == 'EL' or param['OUT_FORM'] == 'XVEL':
             tlab.append('a')
             tlab.append('e')
@@ -577,12 +573,8 @@ def make_swiftest_labels(param):
     """
     tlab = []
     if param['OUT_FORM'] == 'XV' or param['OUT_FORM'] == 'XVEL':
-        tlab.append('rhx')
-        tlab.append('rhy')
-        tlab.append('rhz')
-        tlab.append('vhx')
-        tlab.append('vhy')
-        tlab.append('vhz')
+        tlab.append('rh')
+        tlab.append('vh')
 
     if param['OUT_FORM'] == 'EL' or param['OUT_FORM'] == 'XVEL':
         tlab.append('a')
@@ -599,18 +591,10 @@ def make_swiftest_labels(param):
         plab.append('rhill')
     clab = ['Gmass', 'radius', 'j2rp2', 'j4rp4']
     if param['ROTATION']:
-        clab.append('Ip1')
-        clab.append('Ip2')
-        clab.append('Ip3')
-        clab.append('rotx')
-        clab.append('roty')
-        clab.append('rotz')
-        plab.append('Ip1')
-        plab.append('Ip2')
-        plab.append('Ip3')
-        plab.append('rotx')
-        plab.append('roty')
-        plab.append('rotz')
+        clab.append('Ip')
+        clab.append('rot')
+        plab.append('Ip')
+        plab.append('rot')
     if param['TIDES']:
         clab.append('k2')
         clab.append('Q')
@@ -619,19 +603,11 @@ def make_swiftest_labels(param):
 
     infolab_float = [
         "origin_time",
-        "origin_rhx",
-        "origin_rhy",
-        "origin_rhz",
-        "origin_vhx",
-        "origin_vhy",
-        "origin_vhz",
+        "origin_rh",
+        "origin_vh",
         "discard_time",
-        "discard_rhx",
-        "discard_rhy",
-        "discard_rhz",
-        "discard_vhx",
-        "discard_vhy",
-        "discard_vhz",
+        "discard_rh",
+        "discard_vh",
         ]
     infolab_int = [
         "collision_id",
@@ -1042,7 +1018,7 @@ def swiftest_particle_2xr(param):
     -------
     infoxr : xarray dataset
     """
-    veclab = ['time_origin', 'rhx_origin', 'py_origin', 'pz_origin', 'vhx_origin', 'vhy_origin', 'vhz_origin']
+    veclab = ['time_origin', 'rh_origin', 'vh_origin']
     id_list = []
     origin_type_list = []
     origin_vec_list = []
@@ -1099,7 +1075,7 @@ def select_active_from_frame(ds, param, framenum=-1):
     # Select only the active particles at this time step
     # Remove the inactive particles
     if param['OUT_FORM'] == 'XV' or param['OUT_FORM'] == 'XVEL':
-        iactive = iframe[count_dim].where((~np.isnan(iframe['Gmass'])) | (~np.isnan(iframe['rhx'])), drop=True)[count_dim]
+        iactive = iframe[count_dim].where((~np.isnan(iframe['Gmass'])) | (~np.isnan(iframe['rh'].isel(space=0))), drop=True)[count_dim]
     else:
         iactive = iframe[count_dim].where((~np.isnan(iframe['Gmass'])) | (~np.isnan(iframe['a'])), drop = True)[count_dim]
     if count_dim == "id":
@@ -1164,12 +1140,12 @@ def swiftest_xr2infile(ds, param, in_type="NETCDF_DOUBLE", infile_name=None,fram
     J4 = np.double(cb['j4rp4'])
     cbname = cb['name'].values[0]
     if param['ROTATION']:
-        Ip1cb = np.double(cb['Ip1'])
-        Ip2cb = np.double(cb['Ip2'])
-        Ip3cb = np.double(cb['Ip3'])
-        rotxcb = np.double(cb['rotx'])
-        rotycb = np.double(cb['roty'])
-        rotzcb = np.double(cb['rotz'])
+        Ip1cb = np.double(cb['Ip'].values[0])
+        Ip2cb = np.double(cb['Ip'].values[1])
+        Ip3cb = np.double(cb['Ip'].values[2])
+        rotxcb = np.double(cb['rot'].values[0])
+        rotycb = np.double(cb['rot'].values[1])
+        rotzcb = np.double(cb['rot'].values[2])
     cbid = int(0)
 
     if in_type == 'ASCII':
@@ -1196,16 +1172,16 @@ def swiftest_xr2infile(ds, param, in_type="NETCDF_DOUBLE", infile_name=None,fram
             if param['CHK_CLOSE']:
                print(pli['radius'].values[0], file=plfile)
             if param['IN_FORM'] == 'XV':
-                print(pli['rhx'].values[0], pli['rhy'].values[0], pli['rhz'].values[0], file=plfile)
-                print(pli['vhx'].values[0], pli['vhy'].values[0], pli['vhz'].values[0], file=plfile)
+                print(pli['rh'].values[0,0], pli['rh'].values[0,1], pli['rh'].values[0,2], file=plfile)
+                print(pli['vh'].values[0,0], pli['vh'].values[0,1], pli['vh'].values[0,2], file=plfile)
             elif param['IN_FORM'] == 'EL':
                 print(pli['a'].values[0], pli['e'].values[0], pli['inc'].values[0], file=plfile)
                 print(pli['capom'].values[0], pli['omega'].values[0], pli['capm'].values[0], file=plfile)
             else:
                 print(f"{param['IN_FORM']} is not a valid input format type.")
             if param['ROTATION']:
-                print(pli['Ip1'].values[0], pli['Ip2'].values[0], pli['Ip3'].values[0], file=plfile)
-                print(pli['rotx'].values[0], pli['roty'].values[0], pli['rotz'].values[0], file=plfile)
+                print(pli['Ip'].values[0,0], pli['Ip'].values[0,1], pli['Ip'].values[0,2], file=plfile)
+                print(pli['rot'].values[0,0], pli['rot'].values[0,1], pli['rot'].values[0,2], file=plfile)
         plfile.close()
 
         # TP file
@@ -1215,8 +1191,8 @@ def swiftest_xr2infile(ds, param, in_type="NETCDF_DOUBLE", infile_name=None,fram
             tpi = tp.sel(id=i)
             print(tpi['name'].values[0], file=tpfile)
             if param['IN_FORM'] == 'XV':
-                print(tpi['rhx'].values[0], tpi['rhy'].values[0], tpi['rhz'].values[0], file=tpfile)
-                print(tpi['vhx'].values[0], tpi['vhy'].values[0], tpi['vhz'].values[0], file=tpfile)
+                print(tpi['rh'].values[0,0], tpi['rh'].values[0,1], tpi['rh'].values[0,2], file=tpfile)
+                print(tpi['vh'].values[0,0], tpi['vh'].values[0,1], tpi['vh'].values[0,2], file=tpfile)
             elif param['IN_FORM'] == 'EL':
                 print(tpi['a'].values[0], tpi['e'].values[0], tpi['inc'].values[0], file=tpfile)
                 print(tpi['capom'].values[0], tpi['omega'].values[0], tpi['capm'].values[0], file=tpfile)
@@ -1279,8 +1255,8 @@ def swifter_xr2infile(ds, param, framenum=-1):
                 print(i.values, pli['Gmass'].values, file=plfile)
             if param['CHK_CLOSE']:
                 print(pli['radius'].values, file=plfile)
-            print(pli['rhx'].values, pli['rhy'].values, pli['rhz'].values, file=plfile)
-            print(pli['vhx'].values, pli['vhy'].values, pli['vhz'].values, file=plfile)
+            print(pli['rh'].values[0,0], pli['ry'].values[0,1], pli['rh'].values[0,2], file=plfile)
+            print(pli['vh'].values[0,0], pli['vh'].values[0,1], pli['vh'].values[0,2], file=plfile)
         plfile.close()
 
         # TP file
@@ -1289,8 +1265,8 @@ def swifter_xr2infile(ds, param, framenum=-1):
         for i in tp.id:
             tpi = tp.sel(id=i)
             print(i.values, file=tpfile)
-            print(tpi['rhx'].values, tpi['rhy'].values, tpi['rhz'].values, file=tpfile)
-            print(tpi['vhx'].values, tpi['vhy'].values, tpi['vhz'].values, file=tpfile)
+            print(tpi['rh'].values[0,0], tpi['ry'].values[0,1], tpi['rh'].values[0,2], file=tpfile)
+            print(tpi['vh'].values[0,0], tpi['vh'].values[0,1], tpi['vh'].values[0,2], file=tpfile)
         tpfile.close()
     else:
         # Now make Swiftest files

python/swiftest/swiftest/simulation_class.py

@@ -2811,6 +2811,7 @@ def initial_conditions_from_bin(self, framenum=-1, new_param=None, new_param_fil
             A dataset containing the extracted initial condition data.
         """
 
+        frame = None
         if codename != "Swiftest":
             self.save(new_param_file, framenum, codename)
             return

src/encounter/encounter_io.f90

@@ -105,6 +105,7 @@ module subroutine encounter_io_initialize_output(self, param)
             call check( nf90_def_var(nciu%id, nciu%Ip_varname, nciu%out_type, [nciu%space_dimid, nciu%id_dimid, nciu%time_dimid], nciu%Ip_varid), "encounter_io_initialize_output nf90_def_var Ip_varid"  )
             call check( nf90_def_var(nciu%id, nciu%rot_varname, nciu%out_type, [nciu%space_dimid, nciu%id_dimid, nciu%time_dimid], nciu%rot_varid), "encounter_io_initialize_output nf90_def_var rot_varid"  )
          end if
+
          call check( nf90_inquire(nciu%id, nVariables=nvar), "encounter_io_initialize_output nf90_inquire nVariables"  )
          do varid = 1, nvar
             call check( nf90_inquire_variable(nciu%id, varid, xtype=vartype, ndims=ndims), "encounter_io_initialize_output nf90_inquire_variable"  )
@@ -132,23 +133,38 @@ module subroutine encounter_io_initialize_output(self, param)
    end subroutine encounter_io_initialize_output
 
 
-   module subroutine encounter_io_write_frame(self, iu, param)
+   module subroutine encounter_io_write_frame(self, nciu, param)
       !! author: David A. Minton
       !!
       !! Write a frame of output of an encounter list structure.
       implicit none
       ! Arguments
       class(encounter_list),          intent(in)    :: self   !! Swiftest encounter structure
-      class(encounter_io_parameters), intent(inout) :: iu     !! Parameters used to identify a particular encounter io NetCDF dataset
+      class(encounter_io_parameters), intent(inout) :: nciu   !! Parameters used to identify a particular encounter io NetCDF dataset
       class(swiftest_parameters),     intent(inout) :: param  !! Current run configuration parameters
       ! Internals
-      integer(I4B)                                  :: i,old_mode, n
+      integer(I4B)                             :: tslot,i,old_mode, n
+      character(len=NAMELEN)                   :: charstring
+
+      tslot = nciu%ienc_frame
+      call check( nf90_set_fill(nciu%id, nf90_nofill, old_mode), "encounter_io_write_frame_base nf90_set_fill"  )
+      call check( nf90_put_var(nciu%id, nciu%time_varid, self%t, start=[tslot]), "encounter_io_write_frame nf90_put_var time_varid"  )
+
+      ! charstring = trim(adjustl(self%info(j)%name))
+      ! call check( nf90_put_var(nciu%id, nciu%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "encounter_io_write_info_base nf90_put_var name_varid"  )
+
+      ! charstring = trim(adjustl(self%info(j)%particle_type))
+      ! call check( nf90_put_var(nciu%id, nciu%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "encounter_io_write_info_base nf90_put_var particle_type_varid"  )
 
-      i = iu%ienc_frame
-      n = int(self%nenc, kind=I4B)
-      call check( nf90_set_fill(iu%id, nf90_nofill, old_mode), "encounter_io_write_frame_base nf90_set_fill"  )
-      call check( nf90_put_var(iu%id, iu%time_varid, self%t, start=[i]), "encounter_io_write_frame nf90_put_var time_varid"  )
 
+      ! call check( nf90_put_var(nciu%id, nciu%rh_varid, self%rh(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame_base nf90_put_var rh_varid"  )
+      ! call check( nf90_put_var(nciu%id, nciu%vh_varid, self%vh(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame_base nf90_put_var vh_varid"  )
+      ! call check( nf90_put_var(nciu%id, nciu%Gmass_varid, self%Gmass(j), start=[idslot, tslot]), "encounter_io_write_frame_base nf90_put_var body Gmass_varid"  )
+      ! if (param%lclose) call check( nf90_put_var(nciu%id, nciu%radius_varid, self%radius(j), start=[idslot, tslot]), "encounter_io_write_frame_base nf90_put_var body radius_varid"  )
+      ! if (param%lrotation) then
+      !    call check( nf90_put_var(nciu%id, nciu%Ip_varid, self%Ip(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame_base nf90_put_var body Ip_varid"  )
+      !    call check( nf90_put_var(nciu%id, nciu%rot_varid, self%rot(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame_base nf90_put_var body rotx_varid"  )
+      ! end if
 
       return
    end subroutine encounter_io_write_frame

src/encounter/encounter_setup.f90

@@ -77,12 +77,6 @@ module subroutine encounter_setup_list(self, n)
       if (allocated(self%x2)) deallocate(self%x2)
       if (allocated(self%v1)) deallocate(self%v1)
       if (allocated(self%v2)) deallocate(self%v2)
-      if (allocated(self%Gmass1)) deallocate(self%Gmass1)
-      if (allocated(self%Gmass2)) deallocate(self%Gmass2)
-      if (allocated(self%radius1)) deallocate(self%radius1)
-      if (allocated(self%radius2)) deallocate(self%radius2)
-      if (allocated(self%name1)) deallocate(self%name1)
-      if (allocated(self%name2)) deallocate(self%name2)
 
       self%nenc = n
       if (n == 0_I8B) return
@@ -98,12 +92,6 @@ module subroutine encounter_setup_list(self, n)
       allocate(self%x2(NDIM,n))
       allocate(self%v1(NDIM,n))
       allocate(self%v2(NDIM,n))
-      allocate(self%Gmass1(n))
-      allocate(self%Gmass2(n))
-      allocate(self%radius1(n))
-      allocate(self%radius2(n))
-      allocate(self%name1(n))
-      allocate(self%name2(n))
 
       self%lvdotr(:) = .false.
       self%status(:) = INACTIVE
@@ -115,14 +103,6 @@ module subroutine encounter_setup_list(self, n)
       self%x2(:,:) = 0.0_DP
       self%v1(:,:) = 0.0_DP
       self%v2(:,:) = 0.0_DP
-      self%Gmass1(:) = 0.0_DP
-      self%Gmass2(:) = 0.0_DP
-      self%radius1(:) = 0.0_DP
-      self%radius2(:) = 0.0_DP
-      do i = 1_I8B, n
-         self%name1(i) = "UNNAMED"
-         self%name2(i) = "UNNAMED"
-      end do
 
       return
    end subroutine encounter_setup_list

src/encounter/encounter_util.f90

@@ -36,12 +36,6 @@ module subroutine encounter_util_append_list(self, source, 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%Gmass1, source%Gmass1, nold, nsrc, lsource_mask)
-      call util_append(self%Gmass2, source%Gmass2, nold, nsrc, lsource_mask)
-      call util_append(self%radius1, source%radius1, nold, nsrc, lsource_mask)
-      call util_append(self%radius2, source%radius2, nold, nsrc, lsource_mask)
-      call util_append(self%name1, source%name1, nold, nsrc, lsource_mask)
-      call util_append(self%name2, source%name2, nold, nsrc, lsource_mask)
       self%nenc = nold + count(lsource_mask(1:nsrc))
 
       return
@@ -70,12 +64,6 @@ module subroutine encounter_util_copy_list(self, source)
          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%Gmass1(1:n) = source%Gmass1(1:n)
-         self%Gmass2(1:n) = source%Gmass2(1:n)
-         self%radius1(1:n) = source%radius1(1:n)
-         self%radius2(1:n) = source%radius2(1:n)
-         self%name1(1:n) = source%name1(1:n)
-         self%name2(1:n) = source%name2(1:n)
       end associate
 
       return
@@ -115,12 +103,6 @@ module subroutine encounter_util_dealloc_list(self)
       if (allocated(self%x2)) deallocate(self%x2)
       if (allocated(self%v1)) deallocate(self%v1)
       if (allocated(self%v2)) deallocate(self%v2)
-      if (allocated(self%Gmass1)) deallocate(self%Gmass1)
-      if (allocated(self%Gmass2)) deallocate(self%Gmass2)
-      if (allocated(self%radius1)) deallocate(self%radius1)
-      if (allocated(self%radius2)) deallocate(self%radius2)
-      if (allocated(self%name1)) deallocate(self%name1)
-      if (allocated(self%name2)) deallocate(self%name2)
 
       return
    end subroutine encounter_util_dealloc_list
@@ -216,12 +198,6 @@ module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestru
          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%Gmass1, discards%Gmass1, lspill_list, ldestructive)
-         call util_spill(keeps%Gmass2, discards%Gmass2, lspill_list, ldestructive)
-         call util_spill(keeps%radius1, discards%radius1, lspill_list, ldestructive)
-         call util_spill(keeps%radius2, discards%radius2, lspill_list, ldestructive)
-         call util_spill(keeps%name1, discards%name1, lspill_list, ldestructive)
-         call util_spill(keeps%name2, discards%name2, lspill_list, ldestructive)
 
          nenc_old = keeps%nenc