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

Commit

Permalink
More restructruing anc getting close approach saving infrastructure i…
Browse files Browse the repository at this point in the history 
…n place
  • Loading branch information
@daminton
daminton committed Dec 13, 2022
1 parent c0adfd8 commit 64909af
Show file tree
Hide file tree
Showing 11 changed files with 237 additions and 223 deletions.
11 changes: 5 additions & 6 deletions python/swiftest/swiftest/io.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -900,15 +900,14 @@ def string_converter(da):
-------
da : xarray dataset with the strings cleaned up
"""
if da.dtype == np.dtype(object):
da = da.astype('<U32')
da = xstrip_str(da)
elif type(da.values[0]) != np.str_:
da = xstrip_nonstr(da)

da = da.astype('<U32')
da = xstrip_str(da)

return da

def char_converter(da):
"""
"""`
Converts a string to a unicode string
Parameters
Expand Down
14 changes: 3 additions & 11 deletions src/encounter/encounter_setup.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,23 +66,14 @@ module subroutine encounter_setup_list(self, n)
integer(I8B) :: i

if (n < 0) return

if (allocated(self%lvdotr)) deallocate(self%lvdotr)
if (allocated(self%status)) deallocate(self%status)
if (allocated(self%index1)) deallocate(self%index1)
if (allocated(self%index2)) deallocate(self%index2)
if (allocated(self%id1)) deallocate(self%id1)
if (allocated(self%id2)) deallocate(self%id2)
if (allocated(self%x1)) deallocate(self%x1)
if (allocated(self%x2)) deallocate(self%x2)
if (allocated(self%v1)) deallocate(self%v1)
if (allocated(self%v2)) deallocate(self%v2)
call self%dealloc()

self%nenc = n
if (n == 0_I8B) return
self%t = 0.0_DP

allocate(self%lvdotr(n))
allocate(self%lclosest(n))
allocate(self%status(n))
allocate(self%index1(n))
allocate(self%index2(n))
Expand All @@ -94,6 +85,7 @@ module subroutine encounter_setup_list(self, n)
allocate(self%v2(NDIM,n))

self%lvdotr(:) = .false.
self%lclosest(:) = .false.
self%status(:) = INACTIVE
self%index1(:) = 0
self%index2(:) = 0
Expand Down
219 changes: 117 additions & 102 deletions src/encounter/encounter_util.f90
View file

Large diffs are not rendered by default.

3 changes: 2 additions & 1 deletion src/modules/encounter_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@ module encounter_classes
integer(I8B) :: nenc = 0 !! Total number of encounters
logical :: lcollision !! Indicates if the encounter resulted in at least one collision
real(DP) :: t !! Time of encounter
logical, dimension(:), allocatable :: lclosest !! indicates that thie pair of bodies is in currently at its closest approach point
logical, dimension(:), allocatable :: lvdotr !! relative vdotr flag
integer(I4B), dimension(:), allocatable :: status !! status of the interaction
integer(I4B), dimension(:), allocatable :: index1 !! position of the first body in the encounter
Expand Down Expand Up @@ -329,7 +330,7 @@ module subroutine encounter_util_snapshot_collision(self, param, system, t, arg)
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object to store
real(DP), intent(in), optional :: t !! Time of snapshot if different from system time
character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots)
character(*), intent(in), optional :: arg !! "before": takes a snapshot just before the collision. "after" takes the snapshot just after the collision.
end subroutine encounter_util_snapshot_collision

module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg)
Expand Down
2 changes: 1 addition & 1 deletion src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1705,7 +1705,7 @@ module subroutine util_snapshot_system(self, param, system, t, arg)
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object to store
real(DP), intent(in), optional :: t !! Time of snapshot if different from system time
character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots)
character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in encounter snapshots)
end subroutine util_snapshot_system
end interface

Expand Down
8 changes: 4 additions & 4 deletions src/modules/symba_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ module symba_classes
use swiftest_classes, only : swiftest_parameters, swiftest_base, swiftest_particle_info, swiftest_storage, netcdf_parameters
use helio_classes, only : helio_cb, helio_pl, helio_tp, helio_nbody_system
use fraggle_classes, only : fraggle_colliders, fraggle_fragments
use encounter_classes, only : encounter_list, encounter_snapshot, encounter_storage, collision_storage
use encounter_classes, only : encounter_list, encounter_storage, collision_storage
implicit none
public

Expand All @@ -31,8 +31,8 @@ module symba_classes
integer(I4B), dimension(:), allocatable :: seed !! Random seeds for fragmentation modeling
logical :: lfragmentation = .false. !! Do fragmentation modeling instead of simple merger.
character(STRMAX) :: encounter_save = "NONE" !! Indicate if and how encounter data should be saved
logical :: lenc_trajectory_save = .false. !! Indicates that when encounters are saved, the full trajectory through recursion steps are saved
logical :: lenc_closest_save = .false. !! Indicates that when encounters are saved, the closest approach distance between pairs of bodies is saved
logical :: lenc_save_trajectory = .false. !! Indicates that when encounters are saved, the full trajectory through recursion steps are saved
logical :: lenc_save_closest = .false. !! Indicates that when encounters are saved, the closest approach distance between pairs of bodies is saved
type(encounter_storage(nframes=:)), allocatable :: encounter_history !! Stores encounter history for later retrieval and saving to file
type(collision_storage(nframes=:)), allocatable :: collision_history !! Stores encounter history for later retrieval and saving to file
contains
Expand Down Expand Up @@ -212,7 +212,7 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir
implicit none
class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter list object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! current time
real(DP), intent(in) :: dt !! step size
integer(I4B), intent(in) :: irec !! Current recursion level
Expand Down
2 changes: 1 addition & 1 deletion src/setup/setup.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,7 @@ module subroutine setup_construct_system(system, param)

select type(param)
class is (symba_parameters)
if (param%lenc_trajectory_save .or. param%lenc_closest_save) then
if (param%lenc_save_trajectory .or. param%lenc_save_closest) then
allocate(encounter_storage :: param%encounter_history)
associate (encounter_history => param%encounter_history)
allocate(encounter_io_parameters :: encounter_history%nc)
Expand Down
180 changes: 91 additions & 89 deletions src/symba/symba_collision.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -274,7 +274,7 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir
! Arguments
class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter list object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! current time
real(DP), intent(in) :: dt !! step size
integer(I4B), intent(in) :: irec !! Current recursion level
Expand Down Expand Up @@ -302,97 +302,99 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir
class is (symba_pl)
select type(tp => system%tp)
class is (symba_tp)
nenc = self%nenc
allocate(lmask(nenc))
lmask(:) = ((self%status(1:nenc) == ACTIVE) .and. (pl%levelg(self%index1(1:nenc)) >= irec))
if (isplpl) then
lmask(:) = lmask(:) .and. (pl%levelg(self%index2(1:nenc)) >= irec)
else
lmask(:) = lmask(:) .and. (tp%levelg(self%index2(1:nenc)) >= irec)
end if
if (.not.any(lmask(:))) return

allocate(lcollision(nenc))
lcollision(:) = .false.
allocate(lclosest(nenc))
lclosest(:) = .false.

if (isplpl) then
do concurrent(k = 1:nenc, lmask(k))
i = self%index1(k)
j = self%index2(k)
xr(:) = pl%rh(:, i) - pl%rh(:, j)
vr(:) = pl%vb(:, i) - pl%vb(:, j)
rlim = pl%radius(i) + pl%radius(j)
Gmtot = pl%Gmass(i) + pl%Gmass(j)
call symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), Gmtot, rlim, dt, self%lvdotr(k), lcollision(k), lclosest(k))
end do
else
do concurrent(k = 1:nenc, lmask(k))
i = self%index1(k)
j = self%index2(k)
xr(:) = pl%rh(:, i) - tp%rh(:, j)
vr(:) = pl%vb(:, i) - tp%vb(:, j)
call symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%Gmass(i), pl%radius(i), dt, self%lvdotr(k), lcollision(k), lclosest(k))
end do
end if
select type (param)
class is (symba_parameters)
nenc = self%nenc
allocate(lmask(nenc))
lmask(:) = ((self%status(1:nenc) == ACTIVE) .and. (pl%levelg(self%index1(1:nenc)) >= irec))
if (isplpl) then
lmask(:) = lmask(:) .and. (pl%levelg(self%index2(1:nenc)) >= irec)
else
lmask(:) = lmask(:) .and. (tp%levelg(self%index2(1:nenc)) >= irec)
end if
if (.not.any(lmask(:))) return

allocate(lcollision(nenc))
lcollision(:) = .false.
allocate(lclosest(nenc))
lclosest(:) = .false.

if (isplpl) then
do concurrent(k = 1:nenc, lmask(k))
i = self%index1(k)
j = self%index2(k)
xr(:) = pl%rh(:, i) - pl%rh(:, j)
vr(:) = pl%vb(:, i) - pl%vb(:, j)
rlim = pl%radius(i) + pl%radius(j)
Gmtot = pl%Gmass(i) + pl%Gmass(j)
call symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), Gmtot, rlim, dt, self%lvdotr(k), lcollision(k), lclosest(k))
end do
else
do concurrent(k = 1:nenc, lmask(k))
i = self%index1(k)
j = self%index2(k)
xr(:) = pl%rh(:, i) - tp%rh(:, j)
vr(:) = pl%vb(:, i) - tp%vb(:, j)
call symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%Gmass(i), pl%radius(i), dt, self%lvdotr(k), lcollision(k), lclosest(k))
end do
end if

lany_collision = any(lcollision(:))

if (lany_collision) then
call pl%rh2rb(system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary
do k = 1, nenc
i = self%index1(k)
j = self%index2(k)
if (lcollision(k)) self%status(k) = COLLISION
self%tcollision(k) = t
self%x1(:,k) = pl%rh(:,i) + system%cb%rb(:)
self%v1(:,k) = pl%vb(:,i)
if (isplpl) then
self%x2(:,k) = pl%rh(:,j) + system%cb%rb(:)
self%v2(:,k) = pl%vb(:,j)
if (lcollision(k)) then
! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional colliders%idx
if (pl%lcollision(i) .or. pl%lcollision(j)) call pl%make_colliders([i,j])

! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step
pl%lcollision([i, j]) = .true.
pl%status([i, j]) = COLLISION
call pl%info(i)%set_value(status="COLLISION", discard_time=t, discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i))
call pl%info(j)%set_value(status="COLLISION", discard_time=t, discard_rh=pl%rh(:,j), discard_vh=pl%vh(:,j))
end if
else
self%x2(:,k) = tp%rh(:,j) + system%cb%rb(:)
self%v2(:,k) = tp%vb(:,j)
if (lcollision(k)) then
tp%status(j) = DISCARDED_PLR
tp%ldiscard(j) = .true.
write(idstri, *) pl%id(i)
write(idstrj, *) tp%id(j)
write(timestr, *) t
call tp%info(j)%set_value(status="DISCARDED_PLR", discard_time=t, discard_rh=tp%rh(:,j), discard_vh=tp%vh(:,j))
write(message, *) "Particle " // trim(adjustl(tp%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" &
// " collided with massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" &
// " at t = " // trim(adjustl(timestr))
call io_log_one_message(FRAGGLE_LOG_OUT, message)
lany_collision = any(lcollision(:))

if (lany_collision) then
call pl%rh2rb(system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary
do k = 1, nenc
i = self%index1(k)
j = self%index2(k)
if (lcollision(k)) self%status(k) = COLLISION
self%tcollision(k) = t
self%x1(:,k) = pl%rh(:,i) + system%cb%rb(:)
self%v1(:,k) = pl%vb(:,i)
if (isplpl) then
self%x2(:,k) = pl%rh(:,j) + system%cb%rb(:)
self%v2(:,k) = pl%vb(:,j)
if (lcollision(k)) then
! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collider pair
if (pl%lcollision(i) .or. pl%lcollision(j)) call pl%make_colliders([i,j])

! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step
pl%lcollision([i, j]) = .true.
pl%status([i, j]) = COLLISION
call pl%info(i)%set_value(status="COLLISION", discard_time=t, discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i))
call pl%info(j)%set_value(status="COLLISION", discard_time=t, discard_rh=pl%rh(:,j), discard_vh=pl%vh(:,j))
end if
else
self%x2(:,k) = tp%rh(:,j) + system%cb%rb(:)
self%v2(:,k) = tp%vb(:,j)
if (lcollision(k)) then
tp%status(j) = DISCARDED_PLR
tp%ldiscard(j) = .true.
write(idstri, *) pl%id(i)
write(idstrj, *) tp%id(j)
write(timestr, *) t
call tp%info(j)%set_value(status="DISCARDED_PLR", discard_time=t, discard_rh=tp%rh(:,j), discard_vh=tp%vh(:,j))
write(message, *) "Particle " // trim(adjustl(tp%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" &
// " collided with massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" &
// " at t = " // trim(adjustl(timestr))
call io_log_one_message(FRAGGLE_LOG_OUT, message)
end if
end if
end if
end do
end do

! Extract the pl-pl or pl-tp encounter list and return the pl-pl or pl-tp collision_list
select type(self)
class is (symba_plplenc)
call self%extract_collisions(system, param)
class is (symba_pltpenc)
allocate(tmp, mold=self)
call self%spill(tmp, lcollision, ldestructive=.true.) ! Remove this encounter pair from the encounter list
end select
end if
! Extract the pl-pl or pl-tp encounter list and return the pl-pl or pl-tp collision_list
select type(self)
class is (symba_plplenc)
call self%extract_collisions(system, param)
class is (symba_pltpenc)
allocate(tmp, mold=self)
call self%spill(tmp, lcollision, ldestructive=.true.) ! Remove this encounter pair from the encounter list
end select
end if

! Take snapshots of pairs of bodies at close approach (but not collision) if requested
if (any(lclosest(:))) then
! Take snapshots of pairs of bodies at close approach (but not collision) if requested
if (param%lenc_save_closest .and. any(lclosest(:))) call param%encounter_history%take_snapshot(param, system, t, "closest")

end if
end select

end select
end select
Expand Down Expand Up @@ -905,7 +907,7 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param)

call system%colliders%regime(system%fragments, system, param)

if (param%lenc_trajectory_save) call collision_history%take_snapshot(param,system, t, "before")
if (param%lenc_save_trajectory) call collision_history%take_snapshot(param,system, t, "before")
select case (system%fragments%regime)
case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC)
plplcollision_list%status(i) = symba_collision_casedisruption(system, param)
Expand All @@ -917,7 +919,7 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param)
write(*,*) "Error in symba_collision, unrecognized collision regime"
call util_exit(FAILURE)
end select
if (param%lenc_trajectory_save) call collision_history%take_snapshot(param,system, t, "after")
if (param%lenc_save_trajectory) call collision_history%take_snapshot(param,system, t, "after")
deallocate(system%colliders,system%fragments)
end do
end select
Expand Down
Loading

0 comments on commit 64909af

Please sign in to comment.