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

Commit

Permalink
Updated the swiftest_discard_system method to save all discard bodies…
Browse files Browse the repository at this point in the history 
… to the collisions.nc file
  • Loading branch information
@daminton
daminton committed Feb 22, 2024
1 parent 7d5a7df commit 9eee641
Show file tree
Hide file tree
Showing 10 changed files with 227 additions and 158 deletions.
111 changes: 59 additions & 52 deletions src/collision/collision_io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -467,62 +467,69 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param)
if (allocated(tp)) deallocate(tp)
select case(stage)
case(1)
if (.not. allocated(before%pl)) cycle
allocate(pl, source=before%pl)
if (allocated(before%pl)) allocate(pl, source=before%pl)
if (allocated(before%tp)) allocate(tp, source=before%tp)
case(2)
if (.not. allocated(after%pl)) cycle
allocate(pl, source=after%pl)
if (allocated(after%pl)) allocate(pl, source=after%pl)
if (allocated(after%tp)) allocate(tp, source=after%tp)
end select
npl = pl%nbody

! This ensures that there first idslot will have the first body in it, not id 0 which is the default for a new
! idvals array
if (.not.allocated(nc%idvals)) allocate(nc%idvals, source=pl%id)
do i = 1, npl
call nc%find_idslot(pl%id(i), idslot)
call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[ idslot ]), &
"collision_io_netcdf_write_frame_snapshot nf90_put_var id_varid: pl")
charstring = trim(adjustl(pl%info(i)%name))
call netcdf_io_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot ], &
count=[NAMELEN, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var name_varid: pl")
charstring = trim(adjustl(pl%info(i)%particle_type))
call netcdf_io_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot, stage, eslot], &
count=[NAMELEN, 1, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var particle_type_varid: pl")
call netcdf_io_check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1, idslot, stage, eslot], &
count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var rh_varid: pl")
call netcdf_io_check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1, idslot, stage, eslot], &
count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var vh_varid: pl")
call netcdf_io_check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[ idslot, stage, eslot]), &
"collision_io_netcdf_write_frame_snapshot nf90_put_var Gmass_varid: pl")
call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[ idslot, stage, eslot]), &
"collision_io_netcdf_write_frame_snapshot nf90_put_var radius_varid: pl")
if (param%lrotation) then
call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, stage, eslot], &
count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var Ip_varid: pl")
call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i)*RAD2DEG, &
start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), &
"collision_io_netcdf_write_frame_snapshot nf90_put_var rotx_varid: pl")
end if
end do

ntp = tp%nbody
do i = 1, ntp
call nc%find_idslot(tp%id(i), idslot)
call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, tp%id(i), start=[ idslot ]), &
"collision_io_netcdf_write_frame_snapshot nf90_put_var id_varid: tp" )
charstring = trim(adjustl(tp%info(i)%name))
call netcdf_io_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot ], &
count=[NAMELEN, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var name_varid: tp" )
charstring = trim(adjustl(tp%info(i)%particle_type))
call netcdf_io_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot, stage, eslot], &
count=[NAMELEN, 1, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var particle_type_varid: tp" )
call netcdf_io_check( nf90_put_var(nc%id, nc%rh_varid, tp%rh(:,i), start=[1, idslot, stage, eslot], &
count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var rh_varid: tp" )
call netcdf_io_check( nf90_put_var(nc%id, nc%vh_varid, tp%vh(:,i), start=[1, idslot, stage, eslot], &
count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var vh_varid: tp" )
end do

if (.not. (allocated(pl) .or. allocated(tp))) cycle

if (allocated(pl)) then
npl = pl%nbody
! This ensures that there first idslot will have the first body in it, not id 0 which is the default for a new
! idvals array
if (.not.allocated(nc%idvals)) allocate(nc%idvals, source=pl%id)
do i = 1, npl
call nc%find_idslot(pl%id(i), idslot)
call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[ idslot ]), &
"collision_io_netcdf_write_frame_snapshot nf90_put_var id_varid: pl")
charstring = trim(adjustl(pl%info(i)%name))
call netcdf_io_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot ], &
count=[NAMELEN, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var name_varid: pl")
charstring = trim(adjustl(pl%info(i)%particle_type))
call netcdf_io_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot, stage, eslot], &
count=[NAMELEN, 1, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var particle_type_varid: pl")
call netcdf_io_check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1, idslot, stage, eslot], &
count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var rh_varid: pl")
call netcdf_io_check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1, idslot, stage, eslot], &
count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var vh_varid: pl")
call netcdf_io_check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[ idslot, stage, eslot]), &
"collision_io_netcdf_write_frame_snapshot nf90_put_var Gmass_varid: pl")
call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[ idslot, stage, eslot]), &
"collision_io_netcdf_write_frame_snapshot nf90_put_var radius_varid: pl")
if (param%lrotation) then
call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, stage, eslot], &
count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var Ip_varid: pl")
call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i)*RAD2DEG, &
start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), &
"collision_io_netcdf_write_frame_snapshot nf90_put_var rotx_varid: pl")
end if
end do
end if

if (allocated(tp)) then
ntp = tp%nbody
! This ensures that there first idslot will have the first body in it, not id 0 which is the default for a new
! idvals array
if (.not.allocated(nc%idvals)) allocate(nc%idvals, source=tp%id)
do i = 1, ntp
call nc%find_idslot(tp%id(i), idslot)
call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, tp%id(i), start=[ idslot ]), &
"collision_io_netcdf_write_frame_snapshot nf90_put_var id_varid: tp" )
charstring = trim(adjustl(tp%info(i)%name))
call netcdf_io_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot ], &
count=[NAMELEN, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var name_varid: tp" )
charstring = trim(adjustl(tp%info(i)%particle_type))
call netcdf_io_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot, stage, eslot], &
count=[NAMELEN, 1, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var particle_type_varid: tp" )
call netcdf_io_check( nf90_put_var(nc%id, nc%rh_varid, tp%rh(:,i), start=[1, idslot, stage, eslot], &
count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var rh_varid: tp" )
call netcdf_io_check( nf90_put_var(nc%id, nc%vh_varid, tp%vh(:,i), start=[1, idslot, stage, eslot], &
count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var vh_varid: tp" )
end do
end if
end do
end select
end select
Expand Down
6 changes: 3 additions & 3 deletions src/collision/collision_resolve.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -532,7 +532,7 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status)
allocate(plsub, mold=pl)
call pl%spill(plsub, lmask, ldestructive=.false.)

call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, nimpactors)])
! call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, nimpactors)])

! Save the before/after snapshots
select type(before => collider%before)
Expand Down Expand Up @@ -644,7 +644,7 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec)
! Destroy the collision list now that the collisions are resolved
call plpl_collision%setup(0_I8B)

if ((nbody_system%pl_adds%nbody == 0) .and. (nbody_system%pl_discards%nbody == 0)) exit
if ((nbody_system%pl_adds%nbody == 0) .and. (.not.any(pl%ldiscard(:)))) exit
if (allocated(idnew)) deallocate(idnew)
nnew = nbody_system%pl_adds%nbody
allocate(idnew, source=nbody_system%pl_adds%id)
Expand Down Expand Up @@ -743,7 +743,7 @@ module subroutine collision_resolve_pltp(self, nbody_system, param, t, dt, irec)
! Restructure the massive bodies based on the outcome of the collision
call tp%rearray(nbody_system, param)

! Discard the collider
! Check for discards
call nbody_system%tp%discard(nbody_system, param)

associate(idx1 => pltp_collision%index1, idx2 => pltp_collision%index2)
Expand Down
1 change: 1 addition & 0 deletions src/helio/helio_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ module subroutine helio_util_setup_initialize_system(self, system_history, param
call self%tp%h2b(self%cb)

! Make sure that the discard list gets allocated initially
call self%pl_discards%setup(0, param)
call self%tp_discards%setup(0, param)
call self%pl%set_mu(self%cb)
call self%tp%set_mu(self%cb)
Expand Down
1 change: 1 addition & 0 deletions src/rmvs/rmvs_discard.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ module subroutine rmvs_discard_tp(self, nbody_system, param)
call swiftest_io_log_one_message(COLLISION_LOG_OUT,message)
tp%ldiscard(i) = .true.
tp%lmask(i) = .false.
pl%ldiscard(iplperP) = .true.
call tp%info(i)%set_value(status="DISCARDED_PLQ", discard_time=t, discard_rh=tp%rh(:,i), &
discard_vh=tp%vh(:,i), discard_body_id=pl%id(iplperP))
end if
Expand Down
78 changes: 56 additions & 22 deletions src/swiftest/swiftest_discard.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,33 +21,80 @@ module subroutine swiftest_discard_system(self, param)
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
! Internals
logical :: lpl_discards, ltp_discards, lpl_check, ltp_check
logical, dimension(:), allocatable :: ldiscard
integer(I4B) :: i, nstart, nend, nsub
character(len=STRMAX) :: idstr
class(swiftest_pl), allocatable :: plsub
class(swiftest_tp), allocatable :: tpsub

lpl_check = allocated(self%pl_discards)
ltp_check = allocated(self%tp_discards)

associate(nbody_system => self,tp => self%tp,pl => self%pl,tp_discards => self%tp_discards,pl_discards => self%pl_discards)
associate(nbody_system => self,tp => self%tp,pl => self%pl,tp_discards => self%tp_discards,pl_discards => self%pl_discards, &
npl => self%pl%nbody, ntp => self%tp%nbody, t => self%t, collision_history => self%collision_history, &
collider => self%collider)
lpl_discards = .false.
ltp_discards = .false.
if (lpl_check .and. pl%nbody > 0) then
pl%ldiscard = pl%status(:) /= ACTIVE
call pl%discard(nbody_system, param)
lpl_discards = (pl_discards%nbody > 0)
lpl_discards = any(pl%ldiscard(1:npl))
end if

if (ltp_check .and. tp%nbody > 0) then
tp%ldiscard = tp%status(:) /= ACTIVE
call tp%discard(nbody_system, param)
ltp_discards = (tp_discards%nbody > 0)
ltp_discards = any(tp%ldiscard(1:ntp))
lpl_discards = any(pl%ldiscard(1:npl))
end if

if (ltp_discards.or.lpl_discards) then
if (lpl_discards) then
if (param%lenergy) call self%conservation_report(param, lterminal=.false.)
call pl_discards%setup(0,param)
end if
if (ltp_discards) then
allocate(ldiscard, source=tp%ldiscard(:))
allocate(tpsub, mold=tp)
call tp%spill(tpsub, ldiscard, ldestructive=.true.)
nsub = tpsub%nbody
nstart = tp_discards%nbody + 1
nend = tp_discards%nbody + nsub
call tp_discards%append(tpsub, lsource_mask=[(.true., i = 1, nsub)])
deallocate(ldiscard)
select type(before => collider%before)
class is (swiftest_nbody_system)
if (allocated(before%tp)) deallocate(before%tp)
allocate(before%tp, source=tp_discards)
end select
call tp_discards%setup(0,param)
end if

if (lpl_discards) then ! In the base integrators, massive bodies are not true discards. The discard is
! simply used to trigger a snapshot.
if (param%lenergy) call self%conservation_report(param, lterminal=.false.)
allocate(ldiscard, source=pl%ldiscard(:))
allocate(plsub, mold=pl)
call pl%spill(plsub, ldiscard, ldestructive=.false.)
nsub = plsub%nbody
nstart = pl_discards%nbody + 1
nend = pl_discards%nbody + nsub
call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, nsub)])
deallocate(ldiscard)
pl%ldiscard(1:npl) = .false.
! Save the before snapshots
select type(before => collider%before)
class is (swiftest_nbody_system)
if (allocated(before%pl)) deallocate(before%pl)
allocate(before%pl, source=pl_discards)
end select
call pl_discards%setup(0,param)
end if
! Advance the collision id number and save it
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
collider%impactors%regime = COLLRESOLVE_REGIME_MERGE
write(idstr,*) collider%collision_id
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "collision_id " // trim(adjustl(idstr)))

call collision_history%take_snapshot(param,nbody_system, t, "particle")
end if

end associate
Expand Down Expand Up @@ -86,15 +133,11 @@ module subroutine swiftest_discard_tp(self, nbody_system, param)
class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object
class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameter
! Internals
logical, dimension(self%nbody) :: ldiscard
integer(I4B) :: i, nstart, nend, nsub
class(swiftest_tp), allocatable :: tpsub

if (self%nbody == 0) return

associate(tp => self, ntp => self%nbody, cb => nbody_system%cb, pl => nbody_system%pl, npl => nbody_system%pl%nbody, &
tp_discards => nbody_system%tp_discards)
tp_discards => nbody_system%tp_discards, pl_discards => nbody_system%pl_discards)

if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. &
(param%rmaxu >= 0.0_DP) .or. ((param%qmin >= 0.0_DP) .and. (param%qmin_coord == "BARY"))) then
Expand All @@ -113,15 +156,6 @@ module subroutine swiftest_discard_tp(self, nbody_system, param)
if (param%lclose) then
call swiftest_discard_pl_tp(tp, nbody_system, param)
end if
if (any(tp%ldiscard(1:ntp))) then
ldiscard(1:ntp) = tp%ldiscard(1:ntp)
allocate(tpsub, mold=tp)
call tp%spill(tpsub, ldiscard, ldestructive=.false.)
nsub = tpsub%nbody
nstart = tp_discards%nbody + 1
nend = tp_discards%nbody + nsub
call tp_discards%append(tpsub, lsource_mask=[(.true., i = 1, nsub)])
end if
end associate

return
Expand Down Expand Up @@ -242,7 +276,7 @@ subroutine swiftest_discard_peri_tp(tp, nbody_system, param)
call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
tp%ldiscard(i) = .true.
call tp%info(i)%set_value(status="DISCARDED_PERI", discard_time=nbody_system%t, discard_rh=tp%rh(:,i), &
discard_vh=tp%vh(:,i), discard_body_id=pl%id(j))
discard_vh=tp%vh(:,i), discard_body_id=cb%id)
end if
end if
end if
Expand Down
Loading

0 comments on commit 9eee641

Please sign in to comment.