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

Commit

Permalink
Fixed a number of issues related to saving particle ids in NetCDF fil…
Browse files Browse the repository at this point in the history 
…es and fixing a problem with the merge after Fraggle failure where the mass was not being set properly
  • Loading branch information
@daminton
daminton committed Jan 12, 2023
1 parent b7acdaa commit a794972
Show file tree
Hide file tree
Showing 12 changed files with 3,093 additions and 60 deletions.
2 changes: 0 additions & 2 deletions python/swiftest/swiftest/io.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -813,14 +813,12 @@ def process_netcdf_input(ds, param):
-------
ds : xarray dataset
"""
#

if param['OUT_TYPE'] == "NETCDF_DOUBLE":
ds = fix_types(ds,ftype=np.float64)
elif param['OUT_TYPE'] == "NETCDF_FLOAT":
ds = fix_types(ds,ftype=np.float32)

ds = ds.where(ds['id']>=0, drop=True)
return ds

def swiftest2xr(param, verbose=True):
Expand Down
2 changes: 0 additions & 2 deletions src/base/base_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,8 +22,6 @@ module base
type, abstract :: base_parameters
character(len=:), allocatable :: integrator !! Name of the nbody integrator used
character(len=:), allocatable :: param_file_name !! The name of the parameter file
integer(I4B) :: maxid = -1 !! The current maximum particle id number
integer(I4B) :: maxid_collision = 0 !! The current maximum collision id number
real(DP) :: t0 = 0.0_DP !! Integration reference time
real(DP) :: tstart = -1.0_DP !! Integration start time
real(DP) :: tstop = -1.0_DP !! Integration stop time
Expand Down
11 changes: 5 additions & 6 deletions src/collision/collision_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -177,12 +177,12 @@ module subroutine collision_generate_merge(self, nbody_system, param, t)
! Internals
integer(I4B) :: i, j, k, ibiggest
real(DP), dimension(NDIM) :: L_spin_new
real(DP) :: volume, G
real(DP) :: volume
character(len=STRMAX) :: message

select type(nbody_system)
class is (swiftest_nbody_system)
associate(impactors => nbody_system%collider%impactors)
associate(impactors => self%impactors)
message = "Merging"
call collision_io_collider_message(nbody_system%pl, impactors%id, message)
call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
Expand All @@ -194,7 +194,7 @@ module subroutine collision_generate_merge(self, nbody_system, param, t)

! Generate the merged body as a single fragment
call self%setup_fragments(1)
associate(fragments => nbody_system%collider%fragments)
associate(fragments => self%fragments)

! Calculate the initial energy of the nbody_system without the collisional family
call self%get_energy_and_momentum(nbody_system, param, phase="before")
Expand All @@ -207,9 +207,8 @@ module subroutine collision_generate_merge(self, nbody_system, param, t)

! Compute the physical properties of the new body after the merge.
volume = 4._DP / 3._DP * PI * sum(impactors%radius(:)**3)
G = nbody_system%collider%impactors%Gmass(1) / nbody_system%collider%impactors%mass(1)
fragments%mass(1) = impactors%mass_dist(1)
fragments%Gmass(1) = G * fragments%mass(1)
fragments%mass(1) = sum(impactors%mass(:))
fragments%Gmass(1) =sum(impactors%Gmass(:))
fragments%density(1) = fragments%mass(1) / volume
fragments%radius(1) = (3._DP * volume / (4._DP * PI))**(THIRD)
if (param%lrotation) then
Expand Down
13 changes: 7 additions & 6 deletions src/collision/collision_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -142,12 +142,13 @@ module collision
!! to resolve collision by defining extended types of encounters_impactors and/or encounetr_fragments
!!
!! The generate method for this class is the merge model. This allows any extended type to have access to the merge procedure by selecting the collision_basic parent class
class(collision_fragments(:)), allocatable :: fragments !! Object containing information on the pre-collision system
class(collision_impactors), allocatable :: impactors !! Object containing information on the post-collision system
class(base_nbody_system), allocatable :: before !! A snapshot of the subset of the nbody_system involved in the collision
class(base_nbody_system), allocatable :: after !! A snapshot of the subset of the nbody_system containing products of the collision
integer(I4B) :: status !! Status flag to pass to the collision list once the collision has been resolved
integer(I4B) :: collision_id !! ID number of this collision event
class(collision_fragments(:)), allocatable :: fragments !! Object containing information on the pre-collision system
class(collision_impactors), allocatable :: impactors !! Object containing information on the post-collision system
class(base_nbody_system), allocatable :: before !! A snapshot of the subset of the nbody_system involved in the collision
class(base_nbody_system), allocatable :: after !! A snapshot of the subset of the nbody_system containing products of the collision
integer(I4B) :: status !! Status flag to pass to the collision list once the collision has been resolved
integer(I4B) :: collision_id !! ID number of this collision event
integer(I4B) :: maxid_collision = 0 !! The current maximum collision id number

! Scale factors used to scale dimensioned quantities to a more "natural" system where important quantities (like kinetic energy, momentum) are of order ~1
real(DP) :: dscale = 1.0_DP !! Distance dimension scale factor
Expand Down
18 changes: 9 additions & 9 deletions src/collision/collision_resolve.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -378,7 +378,7 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status)
write(newname, FRAGFMT) fragments%id(i)
call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=t, name=newname, &
origin_rh=plnew%rh(:,i), origin_vh=plnew%vh(:,i), &
collision_id=param%maxid_collision)
collision_id=collider%maxid_collision)
end do
do i = 1, nimpactors
if (impactors%id(i) == ibiggest) then
Expand All @@ -402,7 +402,7 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status)
write(newname, FRAGFMT) fragments%id(i)
call plnew%info(i)%set_value(origin_type=origin_type, origin_time=t, name=newname, &
origin_rh=plnew%rh(:,i), origin_vh=plnew%vh(:,i), &
collision_id=param%maxid_collision)
collision_id=collider%maxid_collision)
end do
do i = 1, nimpactors
if (impactors%id(i) == ibiggest) cycle
Expand All @@ -414,8 +414,8 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status)
case(MERGED)
write(origin_type,*) "Merger"
call plnew%info(1)%copy(pl%info(ibiggest))
param%maxid = param%maxid + 1
plnew%id(1) = param%maxid
nbody_system%maxid = nbody_system%maxid + 1
plnew%id(1) = nbody_system%maxid

! Appends an index number to the end of the original name to make it unique, but still identifiable as the original.
! If there is already an index number appended, replace it
Expand All @@ -427,7 +427,7 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status)
plnew%status(1) = NEW_PARTICLE
call plnew%info(1)%set_value(origin_type=origin_type, origin_time=t, name=newname, &
origin_rh=plnew%rh(:,1), origin_vh=plnew%vh(:,1), &
collision_id=param%maxid_collision)
collision_id=collider%maxid_collision)
do i = 1, nimpactors
if (impactors%id(i) == ibiggest) cycle

Expand Down Expand Up @@ -559,11 +559,11 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec)
if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLIDED)) cycle

! Advance the collision id number and save it
param%maxid_collision = max(param%maxid_collision, maxval(nbody_system%pl%info(:)%collision_id))
param%maxid_collision = param%maxid_collision + 1
collider%collision_id = param%maxid_collision
collider%maxid_collision = max(collider%maxid_collision, maxval(nbody_system%pl%info(:)%collision_id))
collider%maxid_collision = collider%maxid_collision + 1
collider%collision_id = collider%maxid_collision
write(idstr,*) collider%collision_id
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Collision id " // trim(adjustl(idstr)))
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "collision_id " // trim(adjustl(idstr)))

! Get the collision regime
call impactors%get_regime(nbody_system, param)
Expand Down
Loading

0 comments on commit a794972

Please sign in to comment.