Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Fixed issues getting the before/after collision snapshots properly re…
Browse files Browse the repository at this point in the history
…corded
  • Loading branch information
daminton committed Jan 5, 2023
1 parent 8afaf54 commit c25ac25
Show file tree
Hide file tree
Showing 4 changed files with 42 additions and 33 deletions.
2 changes: 0 additions & 2 deletions python/swiftest/swiftest/io.py
Original file line number Diff line number Diff line change
Expand Up @@ -804,8 +804,6 @@ def process_netcdf_input(ds, param):
Performs several tasks to convert raw NetCDF files output by the Fortran program into a form that
is used by the Python side. These include:
- Ensuring all types are correct
- Removing any bad id values (empty id slots)
- Swapping the id and name dimension if the names are unique
Parameters
----------
Expand Down
53 changes: 24 additions & 29 deletions python/swiftest/swiftest/simulation_class.py
Original file line number Diff line number Diff line change
Expand Up @@ -211,12 +211,12 @@ def __init__(self,read_param: bool = False, read_old_output: bool = False, simdi
but not each other (SyMBA only). *Note.* Only mtiny or gmtiny is accepted, not both.
Parameter input file equivalent: `GMTINY`
close_encounter_check : bool, default True
Check for close encounters between bodies. If set to True, then the radii of massive bodies must be included
Check for close encounter between bodies. If set to True, then the radii of massive bodies must be included
in initial conditions.
Parameter input file equivalent: `CHK_CLOSE`
encounter_save : {"NONE","TRAJECTORY","CLOSEST", "BOTH"}, default "NONE"
Indicate if and how encounter data should be saved. If set to "TRAJECTORY", the position and velocity vectors
of all bodies undergoing close encounters are saved at each intermediate step to the encounter files.
of all bodies undergoing close encounter are saved at each intermediate step to the encounter files.
If set to "CLOSEST", the position and velocities at the point of closest approach between pairs of bodies are
computed and stored to the encounter files. If set to "BOTH", then this stores the values that would be computed
in "TRAJECTORY" and "CLOSEST". If set to "NONE" no trajectory information is saved.
Expand Down Expand Up @@ -316,8 +316,8 @@ def __init__(self,read_param: bool = False, read_old_output: bool = False, simdi
self.param = {}
self.data = xr.Dataset()
self.init_cond = xr.Dataset()
self.encounters = xr.Dataset()
self.collisions = xr.Dataset()
self.encounter = xr.Dataset()
self.collision = xr.Dataset()

self.simdir = Path(simdir)
if self.simdir.exists():
Expand Down Expand Up @@ -1035,11 +1035,11 @@ def set_feature(self,
Parameters
----------
close_encounter_check : bool, optional
Check for close encounters between bodies. If set to True, then the radii of massive bodies must be included
Check for close encounter between bodies. If set to True, then the radii of massive bodies must be included
in initial conditions.
encounter_save : {"NONE","TRAJECTORY","CLOSEST","BOTH"}, default "NONE"
Indicate if and how encounter data should be saved. If set to "TRAJECTORY", the position and velocity vectors
of all bodies undergoing close encounters are saved at each intermediate step to the encounter files.
of all bodies undergoing close encounter are saved at each intermediate step to the encounter files.
If set to "CLOSEST", the position and velocities at the point of closest approach between pairs of bodies are
computed and stored to the encounter files. If set to "BOTH", then this stores the values that would be computed
in "TRAJECTORY" and "CLOSEST". If set to "NONE" no trajectory information is saved.
Expand Down Expand Up @@ -2738,8 +2738,8 @@ def read_output_file(self,read_init_cond : bool = True):
else:
self.init_cond = self.data.isel(time=0)

self.read_encounters()
self.read_collisions()
self.read_encounter()
self.read_collision()
if self.verbose:
print("Finished reading Swiftest dataset files.")

Expand All @@ -2752,13 +2752,13 @@ def read_output_file(self,read_init_cond : bool = True):
warnings.warn('Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".',stacklevel=2)
return

def read_encounters(self):
def read_encounter(self):
enc_files = glob(f"{self.simdir}{os.path.sep}encounter_*.nc")
if len(enc_files) == 0:
return

if self.verbose:
print("Reading encounter history file as .encounters")
print("Reading encounter history file as .encounter")

enc_files.sort()

Expand All @@ -2767,33 +2767,28 @@ def _preprocess(ds, param):
return io.process_netcdf_input(ds,param)
partial_func = partial(_preprocess, param=self.param)

self.encounters = xr.open_mfdataset(enc_files,parallel=True,combine="nested",concat_dim="time",join="left",preprocess=partial_func,mask_and_scale=True)
self.encounters = io.process_netcdf_input(self.encounters, self.param)
self.encounter = xr.open_mfdataset(enc_files,parallel=True,combine="nested",concat_dim="time",join="left",preprocess=partial_func,mask_and_scale=True)
self.encounter = io.process_netcdf_input(self.encounter, self.param)
# Remove any overlapping time values
tgood,tid = np.unique(self.encounters.time,return_index=True)
self.encounters = self.encounters.isel(time=tid)
tgood,tid = np.unique(self.encounter.time,return_index=True)
self.encounter = self.encounter.isel(time=tid)
# Remove any NaN values
tgood=self.encounters.time.where(~np.isnan(self.encounters.time),drop=True)
self.encounters = self.encounters.sel(time=tgood)
tgood=self.encounter.time.where(~np.isnan(self.encounter.time),drop=True)
self.encounter = self.encounter.sel(time=tgood)

return

def read_collisions(self):
col_files = glob(f"{self.simdir}{os.path.sep}collision_*.nc")
if len(col_files) == 0:
return
def read_collision(self):

col_file = self.simdir / "collision.nc"
if not os.path.exists(col_file):
return

col_files.sort()
if self.verbose:
print("Reading collision history file as .collisions")

# This is needed in order to pass the param argument down to the io.process_netcdf_input function
def _preprocess(ds, param):
return io.process_netcdf_input(ds,param)
partial_func = partial(_preprocess, param=self.param)
print("Reading collision history file as .collision")

self.collisions = xr.open_mfdataset(col_files,parallel=True, combine="nested", concat_dim="collision", preprocess=partial_func,mask_and_scale=True)
self.collisions = io.process_netcdf_input(self.collisions, self.param)
self.collision = xr.open_dataset(col_file)
self.collision = io.process_netcdf_input(self.collision, self.param)

return

Expand Down
3 changes: 1 addition & 2 deletions src/collision/collision_io.f90
Original file line number Diff line number Diff line change
Expand Up @@ -105,11 +105,10 @@ module subroutine collision_io_netcdf_dump(self, param)

call nc%open(param)

do i = 1, self%nframes
do i = 1, self%iframe
if (allocated(self%frame(i)%item)) then
select type(snapshot => self%frame(i)%item)
class is (collision_snapshot)
param%ioutput = i
call snapshot%write_frame(self,param)
end select
else
Expand Down
17 changes: 17 additions & 0 deletions src/collision/collision_util.f90
Original file line number Diff line number Diff line change
Expand Up @@ -695,6 +695,7 @@ module subroutine collision_util_snapshot(self, param, nbody_system, t, arg)
snapshot%collider%pe(1) = nbody_system%pe
snapshot%collider%be(1) = nbody_system%be
snapshot%collider%te(1) = nbody_system%te

case("after")
! Get record the energy of the sytem after the collision
call nbody_system%get_energy_and_momentum(param)
Expand All @@ -707,6 +708,22 @@ module subroutine collision_util_snapshot(self, param, nbody_system, t, arg)
snapshot%collider%be(2) = nbody_system%be
snapshot%collider%te(2) = nbody_system%te

select type(before_snap => snapshot%collider%before )
class is (swiftest_nbody_system)
select type(before_orig => nbody_system%collider%before)
class is (swiftest_nbody_system)
call move_alloc(before_orig%pl, before_snap%pl)
end select
end select

select type(after_snap => snapshot%collider%after )
class is (swiftest_nbody_system)
select type(after_orig => nbody_system%collider%after)
class is (swiftest_nbody_system)
call move_alloc(after_orig%pl, after_snap%pl)
end select
end select

! Save the snapshot for posterity
call collision_util_save_snapshot(nbody_system%collision_history,snapshot)
deallocate(snapshot)
Expand Down

0 comments on commit c25ac25

Please sign in to comment.