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

Commit

Permalink
Improved the way that semi-interacting bodies are identified from input
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed May 25, 2023
1 parent e10f09e commit 660eb48
Show file tree
Hide file tree
Showing 2 changed files with 39 additions and 51 deletions.
10 changes: 5 additions & 5 deletions examples/Basic_Simulation/initial_conditions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,13 +54,13 @@
capom_pl = default_rng().uniform(0.0, 360.0, npl)
omega_pl = default_rng().uniform(0.0, 360.0, npl)
capm_pl = default_rng().uniform(0.0, 360.0, npl)
GM_pl = (np.array([6e23, 8e23, 1e24, 3e24, 5e24]) / sim.param['MU2KG']) * sim.GU
R_pl = np.full(npl, (3 * (GM_pl / sim.GU) / (4 * np.pi * density_pl)) ** (1.0 / 3.0))
Rh_pl = a_pl * ((GM_pl) / (3 * sim.GU)) ** (1.0 / 3.0)
M_pl = np.array([6e23, 8e23, 1e24, 3e24, 5e24]) * sim.KG2MU
R_pl = np.full(npl, (3 * M_pl/ (4 * np.pi * density_pl)) ** (1.0 / 3.0))
Ip_pl = np.full((npl,3),0.4,)
rot_pl = np.zeros((npl,3))
mtiny = 1.01 * np.max(M_pl)

sim.add_body(name=name_pl, a=a_pl, e=e_pl, inc=inc_pl, capom=capom_pl, omega=omega_pl, capm=capm_pl, Gmass=GM_pl, radius=R_pl, rhill=Rh_pl, Ip=Ip_pl, rot=rot_pl)
sim.add_body(name=name_pl, a=a_pl, e=e_pl, inc=inc_pl, capom=capom_pl, omega=omega_pl, capm=capm_pl, mass=M_pl, radius=R_pl, Ip=Ip_pl, rot=rot_pl)

# Add 10 user-defined test particles.
ntp = 10
Expand All @@ -74,7 +74,7 @@
capm_tp = default_rng().uniform(0.0, 360.0, ntp)

sim.add_body(name=name_tp, a=a_tp, e=e_tp, inc=inc_tp, capom=capom_tp, omega=omega_tp, capm=capm_tp)
sim.set_parameter(tstart=0.0, tstop=1.0e3, dt=0.01, istep_out=100, dump_cadence=0)
sim.set_parameter(tstart=0.0, tstop=1.0e3, dt=0.01, istep_out=100, dump_cadence=0, compute_conservation_values=True, mtiny=mtiny)
# Display the run configuration parameters.
sim.get_parameter()
sim.save()
Expand Down
80 changes: 34 additions & 46 deletions src/swiftest/swiftest_io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -123,8 +123,6 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal)
real(DP), dimension(NDIM) :: L_total_now, L_orbit_now, L_spin_now
real(DP) :: ke_orbit_now, ke_spin_now, pe_now, E_orbit_now, be_now, be_cb_now, be_cb_orig, te_now
real(DP) :: GMtot_now
character(len=STRMAX) :: errmsg
integer(I4B), parameter :: EGYIU = 72
character(len=*), parameter :: EGYTERMFMT = '(" DL/L0 = ", ES12.5, "; DE_orbit/|E0| = ", ES12.5, "; DE_total/|E0| = ", ES12.5, "; DM/M0 = ", ES12.5)'

associate(nbody_system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody, display_unit => param%display_unit)
Expand Down Expand Up @@ -204,10 +202,6 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal)
end associate

return

667 continue
write(*,*) "Error writing energy and momentum tracking file: " // trim(adjustl(errmsg))
call base_util_exit(FAILURE)
end subroutine swiftest_io_conservation_report


Expand All @@ -226,7 +220,6 @@ module subroutine swiftest_io_display_run_information(self, param, integration_t
real(DP) :: tfrac !! Fraction of total simulation time completed
type(progress_bar), save :: pbar !! Object used to print out a progress bar
character(len=64) :: pbarmessage
character(*), parameter :: symbacompactfmt = '(";NPLM",ES22.15,$)'
#ifdef COARRAY
character(*), parameter :: co_statusfmt = '("Image: ",I4, "; Time = ", ES12.5, "; fraction done = ", F6.3, ' // &
'"; Number of active pl, tp = ", I6, ", ", I6)'
Expand Down Expand Up @@ -465,7 +458,6 @@ module subroutine swiftest_io_get_args(integrator, param_file_name, display_styl
character(len=STRMAX), dimension(:), allocatable :: arg
integer(I4B), dimension(:), allocatable :: ierr
integer :: i,narg
character(len=*),parameter :: linefmt = '(A)'

narg = command_argument_count()
if (narg > 0) then
Expand Down Expand Up @@ -653,6 +645,7 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, nc, param)
! Internals
integer(I4B) :: itmax, idmax, tslot
real(DP), dimension(:), allocatable :: vals
logical, dimension(:), allocatable :: plmask, tpmask
real(DP), dimension(1) :: rtemp
real(DP), dimension(NDIM) :: rot0, Ip0, L
real(DP) :: mass0
Expand Down Expand Up @@ -689,7 +682,8 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, nc, param)
self%L_total_orig(:) = self%L_orbit_orig(:) + self%L_spin_orig(:)

call netcdf_io_check( nf90_get_var(nc%id, nc%Gmass_varid, vals, start=[1,tslot], count=[idmax,1]), "netcdf_io_get_t0_values_system Gmass_varid" )
self%GMtot_orig = vals(1) + sum(vals(2:idmax), vals(2:idmax) == vals(2:idmax))
call nc%get_valid_masks(plmask,tpmask)
self%GMtot_orig = vals(1) + sum(vals(2:idmax), plmask(:))

cb%GM0 = vals(1)
cb%dGM = cb%Gmass - cb%GM0
Expand Down Expand Up @@ -1030,7 +1024,7 @@ module subroutine swiftest_io_netcdf_open(self, param, readonly)
end subroutine swiftest_io_netcdf_open


module subroutine swiftest_io_netcdf_get_valid_masks(self, plmask, tpmask)
module subroutine swiftest_io_netcdf_get_valid_masks(self, plmask, tpmask, plmmask, Gmtiny)
!! author: David A. Minton
!!
!! Given an open NetCDF, returns logical masks indicating which bodies in the body arrays are active pl and tp type at the current time.
Expand All @@ -1039,15 +1033,17 @@ module subroutine swiftest_io_netcdf_get_valid_masks(self, plmask, tpmask)
use, intrinsic :: ieee_arithmetic
implicit none
! Arguments
class(swiftest_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset
logical, dimension(:), allocatable, intent(out) :: plmask !! Logical mask indicating which bodies are massive bodies
logical, dimension(:), allocatable, intent(out) :: tpmask !! Logical mask indicating which bodies are test particles
class(swiftest_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset
logical, dimension(:), allocatable, intent(out) :: plmask !! Logical mask indicating which bodies are massive bodies
logical, dimension(:), allocatable, intent(out) :: tpmask !! Logical mask indicating which bodies are test particles
logical, dimension(:), allocatable, intent(out), optional :: plmmask !! Logical mask indicating which bodies are fully interacting massive bodies
real(DP), intent(in), optional :: Gmtiny !! The cutoff G*mass between semi-interacting and fully interacting massive bodies
! Internals
real(DP), dimension(:), allocatable :: Gmass, a
real(DP), dimension(:,:), allocatable :: rh
integer(I4B), dimension(:), allocatable :: body_status
logical, dimension(:), allocatable :: lvalid
integer(I4B) :: idmax, status, i
integer(I4B) :: idmax, status
logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes

call ieee_get_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily
Expand Down Expand Up @@ -1094,6 +1090,14 @@ module subroutine swiftest_io_netcdf_get_valid_masks(self, plmask, tpmask)
! Select only active bodies
plmask(:) = plmask(:) .and. lvalid(:)
tpmask(:) = tpmask(:) .and. lvalid(:)

if (present(plmmask) .and. present(Gmtiny)) then
allocate(plmmask, source=plmask)
where(plmask(:))
plmmask = Gmass(:) > Gmtiny
endwhere
end if

call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes)

end associate
Expand Down Expand Up @@ -1316,10 +1320,6 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier

ierr = 0
return

667 continue
write(*,*) "Error reading nbody_system frame in netcdf_io_read_frame_system"

end function swiftest_io_netcdf_read_frame_system


Expand All @@ -1335,21 +1335,15 @@ module subroutine swiftest_io_netcdf_read_hdr_system(self, nc, param)
class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for reading a NetCDF dataset to file
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
! Internals
integer(I4B) :: status, idmax
integer(I4B) :: status
logical, dimension(:), allocatable :: plmask, tpmask, plmmask
real(DP), dimension(:), allocatable :: Gmtemp

associate(tslot => nc%tslot)
call netcdf_io_check( nf90_get_var(nc%id, nc%time_varid, self%t, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar time_varid" )
call nc%get_valid_masks(plmask, tpmask)

if (param%lmtiny_pl) then
idmax = size(plmask)
allocate(plmmask(idmax))
allocate(Gmtemp(idmax))
call netcdf_io_check( nf90_get_var(nc%id, nc%Gmass_varid, Gmtemp, start=[1,tslot], count=[idmax,1]), "netcdf_io_read_hdr_system nf90_getvar Gmass_varid" )
where(Gmtemp(:) /= Gmtemp(:)) Gmtemp(:) = 0.0_DP
plmmask(:) = plmask(:) .and. Gmtemp(:) > param%GMTINY
call nc%get_valid_masks(plmask, tpmask, plmmask, param%GMTINY)
else
call nc%get_valid_masks(plmask, tpmask)
end if

status = nf90_inq_varid(nc%id, nc%npl_varname, nc%npl_varid)
Expand Down Expand Up @@ -1624,7 +1618,7 @@ module subroutine swiftest_io_netcdf_write_frame_body(self, nc, param)
class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
! Internals
integer(I4B) :: i, j, idslot, old_mode
integer(I4B) :: i, j, idslot, old_mode, tmp
integer(I4B), dimension(:), allocatable :: ind
real(DP), dimension(NDIM) :: vh !! Temporary variable to store heliocentric velocity values when converting from pseudovelocity in GR-enabled runs
real(DP) :: a, e, inc, omega, capom, capm, varpi, lam, f, cape, capf
Expand Down Expand Up @@ -1729,7 +1723,7 @@ module subroutine swiftest_io_netcdf_write_frame_body(self, nc, param)
end select
#endif

call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode), "netcdf_io_write_frame_body nf90_set_fill old_mode" )
call netcdf_io_check( nf90_set_fill(nc%id, old_mode, tmp), "netcdf_io_write_frame_body nf90_set_fill old_mode" )

return
end subroutine swiftest_io_netcdf_write_frame_body
Expand All @@ -1745,7 +1739,7 @@ module subroutine swiftest_io_netcdf_write_frame_cb(self, nc, param)
class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
! Internals
integer(I4B) :: idslot, old_mode
integer(I4B) :: idslot, old_mode, tmp

associate(tslot => nc%tslot)
call self%write_info(nc, param)
Expand All @@ -1766,7 +1760,7 @@ module subroutine swiftest_io_netcdf_write_frame_cb(self, nc, param)
call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, self%rot(:) * RAD2DEG, start=[1, idslot, tslot], count=[NDIM,1,1]), "swiftest_io_netcdf_write_frame_cby nf90_put_var cb rot_varid" )
end if

call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode), "swiftest_io_netcdf_write_frame_cb nf90_set_fill old_mode" )
call netcdf_io_check( nf90_set_fill(nc%id, old_mode, tmp), "swiftest_io_netcdf_write_frame_cb nf90_set_fill old_mode" )
end associate

return
Expand Down Expand Up @@ -1848,7 +1842,7 @@ module subroutine swiftest_io_netcdf_write_info_body(self, nc, param)
class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
! Internals
integer(I4B) :: i, j, idslot, old_mode
integer(I4B) :: i, j, idslot, old_mode, tmp
integer(I4B), dimension(:), allocatable :: ind
character(len=NAMELEN) :: charstring
character(len=NAMELEN), dimension(self%nbody) :: origin_type
Expand Down Expand Up @@ -1892,7 +1886,7 @@ module subroutine swiftest_io_netcdf_write_info_body(self, nc, param)
end associate
end select

call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode) )
call netcdf_io_check( nf90_set_fill(nc%id, old_mode, tmp) )
return
end subroutine swiftest_io_netcdf_write_info_body

Expand All @@ -1906,7 +1900,7 @@ module subroutine swiftest_io_netcdf_write_info_cb(self, nc, param)
class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
! Internals
integer(I4B) :: idslot, old_mode
integer(I4B) :: idslot, old_mode, tmp
character(len=NAMELEN) :: charstring

! This string of spaces of length NAMELEN is used to clear out any old data left behind inside the string variables
Expand Down Expand Up @@ -1935,7 +1929,7 @@ module subroutine swiftest_io_netcdf_write_info_cb(self, nc, param)
call netcdf_io_check( nf90_put_var(nc%id, nc%discard_rh_varid, self%info%discard_rh(:), start=[1, idslot], count=[NDIM,1]), "netcdf_io_write_info_body nf90_put_var cb discard_rh_varid" )
call netcdf_io_check( nf90_put_var(nc%id, nc%discard_vh_varid, self%info%discard_vh(:), start=[1, idslot], count=[NDIM,1]), "netcdf_io_write_info_body nf90_put_var cb discard_vh_varid" )
end if
call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode) )
call netcdf_io_check( nf90_set_fill(nc%id, old_mode, tmp) )

return
end subroutine swiftest_io_netcdf_write_info_cb
Expand Down Expand Up @@ -1989,7 +1983,6 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i
character(*),parameter :: linefmt = '(A)' !! Format code for simple text string
integer(I4B) :: nseeds, nseeds_from_file
logical :: seed_set = .false. !! Is the random seed set in the input file?
character(len=:), allocatable :: integrator
real(DP) :: tratio, y
#ifdef COARRAY
type(swiftest_parameters), codimension[*], save :: coparam
Expand Down Expand Up @@ -2517,9 +2510,6 @@ module subroutine swiftest_io_param_writer(self, unit, iotype, v_list, iostat, i
integer, intent(out) :: iostat !! IO status code
character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0
! Internals
character(*),parameter :: Ifmt = '(I0)' !! Format label for integer values
character(*),parameter :: Rfmt = '(ES25.17)' !! Format label for real values
character(*),parameter :: Lfmt = '(L1)' !! Format label for logical values
integer(I4B) :: nseeds

associate(param => self)
Expand Down Expand Up @@ -2588,7 +2578,6 @@ module subroutine swiftest_io_param_writer(self, unit, iotype, v_list, iostat, i
iomsg = "UDIO not implemented"
end associate

667 continue
return
end subroutine swiftest_io_param_writer

Expand Down Expand Up @@ -2935,7 +2924,7 @@ module subroutine swiftest_io_read_in_system(self, nc, param)
if (ierr /=0) call base_util_exit(FAILURE)
end if

param%loblatecb = ((self%cb%j2rp2 /= 0.0_DP) .or. (self%cb%j4rp4 /= 0.0_DP))
param%loblatecb = ((abs(self%cb%j2rp2) > 0.0_DP) .or. (abs(self%cb%j4rp4) > 0.0_DP))
if (.not.param%loblatecb) then
if (allocated(self%pl%aobl)) deallocate(self%pl%aobl)
if (allocated(self%tp%aobl)) deallocate(self%tp%aobl)
Expand Down Expand Up @@ -3083,7 +3072,6 @@ module subroutine swiftest_io_read_in_param(self, param_file_name)
call self%reader(LUN, iotype= "none", v_list=[""], iostat = ierr, iomsg = errmsg)
if (ierr == 0) return

667 continue
write(self%display_unit,*) "Error reading parameter file: " // trim(adjustl(errmsg))
call base_util_exit(FAILURE)
end subroutine swiftest_io_read_in_param
Expand Down Expand Up @@ -3177,7 +3165,7 @@ module subroutine swiftest_io_initialize_output_file_system(self, nc, param)
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
! Internals
logical, save :: lfirst = .true. !! Flag to determine if this is the first call of this method
character(len=STRMAX) :: errmsg
character(len=2*STRMAX) :: errmsg
logical :: fileExists

associate (pl => self%pl, tp => self%tp, npl => self%pl%nbody, ntp => self%tp%nbody)
Expand All @@ -3191,13 +3179,13 @@ module subroutine swiftest_io_initialize_output_file_system(self, nc, param)
select case(param%out_stat)
case('APPEND')
if (.not.fileExists) then
errmsg = param%outfile // " not found! You must specify OUT_STAT = NEW, REPLACE, or UNKNOWN"
errmsg = trim(adjustl(param%outfile)) // " not found! You must specify OUT_STAT = NEW, REPLACE, or UNKNOWN"
goto 667
end if
call nc%open(param)
case('NEW')
if (fileExists) then
errmsg = param%outfile // " Alread Exists! You must specify OUT_STAT = APPEND, REPLACE, or UNKNOWN"
errmsg = trim(adjustl(param%outfile)) // " Alread Exists! You must specify OUT_STAT = APPEND, REPLACE, or UNKNOWN"
goto 667
end if
call nc%initialize(param)
Expand Down

0 comments on commit 660eb48

Please sign in to comment.