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

Commit

Permalink
Added discard info and rearranged code to properly record discard inf…
Browse files Browse the repository at this point in the history 
…ormation. Still hunting for mass loss bug
  • Loading branch information
@daminton
daminton committed Aug 30, 2021
1 parent a2b8b4a commit 71cc587
Show file tree
Hide file tree
Showing 15 changed files with 460 additions and 283 deletions.
8 changes: 6 additions & 2 deletions python/swiftest/swiftest/io.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -718,8 +718,12 @@ def clean_string_values(param, ds):
-------
ds : xarray dataset with the strings cleaned up
"""
ds['name'] = ds['name'].str.decode(encoding='utf-8')
ds['particle_type'] = ds['particle_type'].str.decode(encoding='utf-8')
if 'name' in ds:
ds['name'] = ds['name'].str.decode(encoding='utf-8')
if 'particle_type' in ds:
ds['particle_type'] = ds['particle_type'].str.decode(encoding='utf-8')
if 'status' in ds:
ds['status'] = ds['status'].str.decode(encoding='utf-8')
if 'origin_type' in ds:
ds['origin_type'] = ds['origin_type'].str.decode(encoding='utf-8')
return ds
Expand Down
5 changes: 5 additions & 0 deletions src/discard/discard.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,13 +128,15 @@ subroutine discard_cb_tp(tp, system, param)
write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too far from the central body at t = " // trim(adjustl(timestr))
tp%ldiscard(i) = .true.
tp%lmask(i) = .false.
call tp%info(i)%set_value(status="DISCARDED_RMAX", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i))
else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then
tp%status(i) = DISCARDED_RMIN
write(idstr, *) tp%id(i)
write(timestr, *) param%t
write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too close to the central body at t = " // trim(adjustl(timestr))
tp%ldiscard(i) = .true.
tp%lmask(i) = .false.
call tp%info(i)%set_value(status="DISCARDED_RMIN", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i), discard_body_id=cb%id)
else if (param%rmaxu >= 0.0_DP) then
rb2 = dot_product(tp%xb(:, i), tp%xb(:, i))
vb2 = dot_product(tp%vb(:, i), tp%vb(:, i))
Expand All @@ -146,6 +148,7 @@ subroutine discard_cb_tp(tp, system, param)
write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr))
tp%ldiscard(i) = .true.
tp%lmask(i) = .false.
call tp%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i))
end if
end if
end if
Expand Down Expand Up @@ -195,6 +198,7 @@ subroutine discard_peri_tp(tp, system, param)
write(timestr, *) param%t
write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " perihelion distance too small at t = " // trim(adjustl(timestr))
tp%ldiscard(i) = .true.
call tp%info(i)%set_value(status="DISCARDED_PERI", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i), discard_body_id=pl%id(j))
end if
end if
end if
Expand Down Expand Up @@ -243,6 +247,7 @@ subroutine discard_pl_tp(tp, system, param)
// " too close to massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstrj)) &
// " at t = " // trim(adjustl(timestr))
tp%ldiscard(i) = .true.
call tp%info(i)%set_value(status="DISCARDED_PLR", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i), discard_body_id=pl%id(j))
exit
end if
end do
Expand Down
56 changes: 30 additions & 26 deletions src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -938,7 +938,7 @@ module subroutine io_read_in_body(self, param)
! Internals
integer(I4B), parameter :: LUN = 7 !! Unit number of input file
integer(I4B) :: iu = LUN
integer(I4B) :: nbody
integer(I4B) :: i, nbody
logical :: is_ascii, is_pl
character(len=:), allocatable :: infile
real(DP) :: t
Expand Down Expand Up @@ -974,6 +974,9 @@ module subroutine io_read_in_body(self, param)
ierr = self%read_frame(iu, param)
self%status(:) = ACTIVE
self%lmask(:) = .true.
do i = 1, nbody
call self%info(i)%set_value(status="ACTIVE")
end do
end if
close(iu, err = 667, iomsg = errmsg)

Expand Down Expand Up @@ -1547,7 +1550,7 @@ module subroutine io_write_discard(self, param)
implicit none
! Arguments
class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
! Internals
integer(I4B), parameter :: LUN = 40
integer(I4B) :: i
Expand All @@ -1561,9 +1564,17 @@ module subroutine io_write_discard(self, param)
class(swiftest_body), allocatable :: pltemp
character(len=STRMAX) :: errmsg, out_stat

if (param%discard_out == "") return

associate(tp_discards => self%tp_discards, nsp => self%tp_discards%nbody, pl => self%pl, npl => self%pl%nbody)

! Record the discarded body metadata information to file
if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then
call param%nciu%open(param)
call tp_discards%write_particle_info(param%nciu)
call param%nciu%close(param)
end if

if (param%discard_out == "") return

if (nsp == 0) return
if (lfirst) then
out_stat = param%out_stat
Expand All @@ -1589,22 +1600,22 @@ module subroutine io_write_discard(self, param)
write(LUN, VECFMT, err = 667, iomsg = errmsg) tp_discards%vh(1, i), tp_discards%vh(2, i), tp_discards%vh(3, i)
end do
if (param%lbig_discard) then
if (param%lgr) then
allocate(pltemp, source = pl)
call pltemp%pv2v(param)
allocate(vh, source = pltemp%vh)
deallocate(pltemp)
else
allocate(vh, source = pl%vh)
end if
if (param%lgr) then
allocate(pltemp, source = pl)
call pltemp%pv2v(param)
allocate(vh, source = pltemp%vh)
deallocate(pltemp)
else
allocate(vh, source = pl%vh)
end if

write(LUN, NPLFMT) npl
do i = 1, npl
write(LUN, PLNAMEFMT, err = 667, iomsg = errmsg) pl%id(i), pl%Gmass(i), pl%radius(i)
write(LUN, VECFMT, err = 667, iomsg = errmsg) pl%xh(1, i), pl%xh(2, i), pl%xh(3, i)
write(LUN, VECFMT, err = 667, iomsg = errmsg) vh(1, i), vh(2, i), vh(3, i)
end do
deallocate(vh)
write(LUN, NPLFMT) npl
do i = 1, npl
write(LUN, PLNAMEFMT, err = 667, iomsg = errmsg) pl%id(i), pl%Gmass(i), pl%radius(i)
write(LUN, VECFMT, err = 667, iomsg = errmsg) pl%xh(1, i), pl%xh(2, i), pl%xh(3, i)
write(LUN, VECFMT, err = 667, iomsg = errmsg) vh(1, i), vh(2, i), vh(3, i)
end do
deallocate(vh)
end if
close(LUN)
end associate
Expand Down Expand Up @@ -1878,13 +1889,6 @@ module subroutine io_write_frame_system(self, param)
end if
end select

select type(param)
class is (symba_parameters)
param%nciu%ltrack_origin = param%lfragmentation
class default
param%nciu%ltrack_origin = .false.
end select

select case(param%out_stat)
case('APPEND')
call param%nciu%open(param)
Expand Down
Loading

0 comments on commit 71cc587

Please sign in to comment.