From 0d836c5717ed9a3f8a0a33f85bbc1524ac2cda14 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 27 Aug 2021 17:32:29 -0400 Subject: [PATCH] Major restructuring to bring the particle information to all integrators, not just SyMBA. Name is now contains therein --- src/io/io.f90 | 199 +++++++++++++++++++++++++++++-- src/modules/swiftest_classes.f90 | 176 ++++++++++++++++++++------- src/modules/swiftest_globals.f90 | 17 ++- src/modules/symba_classes.f90 | 119 ++++++------------ src/netcdf/netcdf.f90 | 4 +- src/setup/setup.f90 | 51 +++++++- src/symba/symba_collision.f90 | 196 ++++++++++++++++-------------- src/symba/symba_io.f90 | 142 ++++------------------ src/symba/symba_setup.f90 | 76 +++++------- src/symba/symba_util.f90 | 142 +--------------------- src/util/util_append.f90 | 30 ++++- src/util/util_fill.f90 | 23 +++- src/util/util_resize.f90 | 37 +++++- src/util/util_sort.f90 | 23 +++- src/util/util_spill.f90 | 43 ++++++- 15 files changed, 730 insertions(+), 548 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index 17c75aeea..8dcb37053 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -111,7 +111,90 @@ module subroutine io_dump_param(self, param_file_name) end subroutine io_dump_param - module subroutine io_dump_swiftest(self, param) + module subroutine io_dump_particle_info(self, iu) + !! author: David A. Minton + !! + !! Reads in particle information object information from an open file unformatted file + implicit none + ! Arguments + class(swiftest_particle_info), intent(in) :: self !! Particle metadata information object + integer(I4B), intent(in) :: iu !! Open file unit number + ! Internals + character(STRMAX) :: errmsg + + write(iu, err = 667, iomsg = errmsg) self%name + write(iu, err = 667, iomsg = errmsg) self%particle_type + + return + + 667 continue + write(*,*) "Error writing particle metadata information from file: " // trim(adjustl(errmsg)) + call util_exit(FAILURE) + end subroutine io_dump_particle_info + + + module subroutine io_dump_particle_info_base(self, param, idx) + !! author: David A. Minton + !! + !! Dumps the particle information data to a file. + !! Pass a list of array indices for test particles (tpidx) and/or massive bodies (plidx) to append + implicit none + ! Arguments + class(swiftest_base), intent(inout) :: self !! Swiftest base object (can be cb, pl, or tp) + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + integer(I4B), dimension(:), optional, intent(in) :: idx !! Array of test particle indices to append to the particle file + + ! Internals + logical, save :: lfirst = .true. + integer(I4B), parameter :: LUN = 22 + integer(I4B) :: i + character(STRMAX) :: errmsg + + if (lfirst) then + select case(param%out_stat) + case('APPEND') + open(unit = LUN, file = param%particle_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', err = 667, iomsg = errmsg) + case('NEW', 'UNKNOWN', 'REPLACE') + open(unit = LUN, file = param%particle_out, status = param%out_stat, form = 'UNFORMATTED', err = 667, iomsg = errmsg) + case default + write(*,*) 'Invalid status code',trim(adjustl(param%out_stat)) + call util_exit(FAILURE) + end select + + lfirst = .false. + else + open(unit = LUN, file = param%particle_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', err = 667, iomsg = errmsg) + end if + + select type(self) + class is (swiftest_cb) + write(LUN, err = 667, iomsg = errmsg) self%id + call self%info%dump(LUN) + class is (swiftest_body) + if (present(idx)) then + do i = 1, size(idx) + write(LUN, err = 667, iomsg = errmsg) self%id(idx(i)) + call self%info(idx(i))%dump(LUN) + end do + else + do i = 1, self%nbody + write(LUN, err = 667, iomsg = errmsg) self%id(i) + call self%info(i)%dump(LUN) + end do + end if + end select + + close(unit = LUN, err = 667, iomsg = errmsg) + + return + + 667 continue + write(*,*) "Error reading central body file: " // trim(adjustl(errmsg)) + call util_exit(FAILURE) + end subroutine io_dump_particle_info_base + + + module subroutine io_dump_base(self, param) !! author: David A. Minton !! !! Dump massive body data to files @@ -152,7 +235,7 @@ module subroutine io_dump_swiftest(self, param) 667 continue write(*,*) "Error dumping body data to file " // trim(adjustl(errmsg)) call util_exit(FAILURE) - end subroutine io_dump_swiftest + end subroutine io_dump_base module subroutine io_dump_system(self, param) @@ -554,7 +637,9 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) read(param_value, *, err = 667, iomsg = iomsg) param%Euntracked case ("MAXID") read(param_value, *, err = 667, iomsg = iomsg) param%maxid - case ("NPLMAX", "NTPMAX", "GMTINY", "MIN_GMFRAG", "PARTICLE_OUT", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters + case ("PARTICLE_OUT") + param%particle_out = param_value + case ("NPLMAX", "NTPMAX", "GMTINY", "MIN_GMFRAG", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters case default write(*,*) "Unknown parameter -> ",param_name iostat = -1 @@ -771,6 +856,7 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) write(param_name, Afmt) "OUT_FORM"; write(param_value, Afmt) trim(adjustl(param%out_form)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "OUT_STAT"; write(param_value, Afmt) "APPEND"; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) end if + write(param_name, Afmt) "PARTICLE_OUT"; write(param_value, Afmt) trim(adjustl(param%particle_out)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "ENC_OUT"; write(param_value, Afmt) trim(adjustl(param%enc_out)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "CHK_RMIN"; write(param_value, Rfmt) param%rmin; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "CHK_RMAX"; write(param_value, Rfmt) param%rmax; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) @@ -906,7 +992,7 @@ module subroutine io_read_in_cb(self, param) self%id = 0 param%maxid = 0 open(unit = iu, file = param%incbfile, status = 'old', form = 'FORMATTED', err = 667, iomsg = errmsg) - read(iu, *, err = 667, iomsg = errmsg) self%name + read(iu, *, err = 667, iomsg = errmsg) self%info%name read(iu, *, err = 667, iomsg = errmsg) self%Gmass self%mass = real(self%Gmass / param%GU, kind=DP) read(iu, *, err = 667, iomsg = errmsg) self%radius @@ -1036,7 +1122,9 @@ module function io_read_frame_body(self, iu, param) result(ierr) select case(param%in_type) case (REAL4_TYPE, REAL8_TYPE) read(iu, err = 667, iomsg = errmsg) self%id(:) - read(iu, err = 667, iomsg = errmsg) self%name(:) + associate(name => self%info%name) + read(iu, err = 667, iomsg = errmsg) name(:) + end associate select case (param%in_form) case (XV) @@ -1082,9 +1170,9 @@ module function io_read_frame_body(self, iu, param) result(ierr) select type(self) class is (swiftest_pl) if (param%lrhill_present) then - read(iu, *, err = 667, iomsg = errmsg) self%name(i), val, self%rhill(i) + read(iu, *, err = 667, iomsg = errmsg) self%info(i)%name, val, self%rhill(i) else - read(iu, *, err = 667, iomsg = errmsg) self%name(i), val + read(iu, *, err = 667, iomsg = errmsg) self%info(i)%name, val end if self%Gmass(i) = real(val, kind=DP) self%mass(i) = real(val / param%GU, kind=DP) @@ -1160,7 +1248,7 @@ module function io_read_frame_cb(self, iu, param) result(ierr) character(len=STRMAX) :: errmsg read(iu, err = 667, iomsg = errmsg) self%id - read(iu, err = 667, iomsg = errmsg) self%name + read(iu, err = 667, iomsg = errmsg) self%info%name read(iu, err = 667, iomsg = errmsg) self%Gmass self%mass = self%Gmass / param%GU read(iu, err = 667, iomsg = errmsg) self%radius @@ -1279,6 +1367,7 @@ function io_read_hdr(iu, t, npl, ntp, out_form, out_type) result(ierr) return end function io_read_hdr + module subroutine io_read_in_param(self, param_file_name) !! author: David A. Minton !! @@ -1314,6 +1403,96 @@ module subroutine io_read_in_param(self, param_file_name) end subroutine io_read_in_param + module subroutine io_read_in_particle_info(self, iu) + !! author: David A. Minton + !! + !! Reads in particle information object information from an open file unformatted file + implicit none + ! Arguments + class(swiftest_particle_info), intent(inout) :: self !! Particle metadata information object + integer(I4B), intent(in) :: iu !! Open file unit number + ! Internals + character(STRMAX) :: errmsg + + read(iu, err = 667, iomsg = errmsg) self%name + read(iu, err = 667, iomsg = errmsg) self%particle_type + + return + + 667 continue + write(*,*) "Error reading particle metadata information from file: " // trim(adjustl(errmsg)) + call util_exit(FAILURE) + end subroutine io_read_in_particle_info + + + module subroutine io_read_particle_info_system(self, param) + !! author: David A. Minton + !! + !! Reads an old particle information file for a restartd run + implicit none + ! Arguments + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B), parameter :: LUN = 22 + integer(I4B) :: i, id, idx + logical :: lmatch + character(STRMAX) :: errmsg + class(swiftest_particle_info), allocatable :: tmpinfo + + open(unit = LUN, file = param%particle_out, status = 'OLD', form = 'UNFORMATTED', err = 667, iomsg = errmsg) + + allocate(tmpinfo, mold=self%cb%info) + + select type(cb => self%cb) + class is (swiftest_cb) + select type(pl => self%pl) + class is (swiftest_pl) + select type(tp => self%tp) + class is (swiftest_tp) + associate(npl => pl%nbody, ntp => tp%nbody) + do + lmatch = .false. + read(LUN, err = 667, iomsg = errmsg, end = 333) id + + if (id == cb%id) then + call cb%info%read_in(LUN) + lmatch = .true. + else + if (npl > 0) then + idx = findloc(pl%id(1:npl), id, dim=1) + if (idx /= 0) then + call pl%info(idx)%read_in(LUN) + lmatch = .true. + end if + end if + if (.not.lmatch .and. ntp > 0) then + idx = findloc(tp%id(1:ntp), id, dim=1) + if (idx /= 0) then + call tp%info(idx)%read_in(LUN) + lmatch = .true. + end if + end if + end if + if (.not.lmatch) then + call tmpinfo%read_in(LUN) + end if + end do + end associate + close(unit = LUN, err = 667, iomsg = errmsg) + end select + end select + end select + + 333 continue + return + + 667 continue + write(*,*) "Error reading particle information file: " // trim(adjustl(errmsg)) + call util_exit(FAILURE) + end subroutine io_read_particle_info_system + + module subroutine io_toupper(string) !! author: David A. Minton !! @@ -1495,7 +1674,7 @@ module subroutine io_write_frame_body(self, iu, param) associate(n => self%nbody) if (n == 0) return write(iu, err = 667, iomsg = errmsg) self%id(1:n) - write(iu, err = 667, iomsg = errmsg) self%name(1:n) + write(iu, err = 667, iomsg = errmsg) self%info(1:n)%name if ((param%out_form == XV) .or. (param%out_form == XVEL)) then write(iu, err = 667, iomsg = errmsg) self%xh(1, 1:n) write(iu, err = 667, iomsg = errmsg) self%xh(2, 1:n) @@ -1556,7 +1735,7 @@ module subroutine io_write_frame_cb(self, iu, param) associate(cb => self) write(iu, err = 667, iomsg = errmsg) cb%id - write(iu, err = 667, iomsg = errmsg) cb%name + write(iu, err = 667, iomsg = errmsg) cb%info%name write(iu, err = 667, iomsg = errmsg) cb%Gmass write(iu, err = 667, iomsg = errmsg) cb%radius write(iu, err = 667, iomsg = errmsg) cb%j2rp2 diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 06435cbe7..34e530bb9 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -31,6 +31,7 @@ module swiftest_classes character(STRMAX) :: out_type = NETCDF_DOUBLE_TYPE !! Binary format of output file character(STRMAX) :: out_form = XVEL !! Data to write to output file character(STRMAX) :: out_stat = 'NEW' !! Open status for output binary file + character(STRMAX) :: particle_out = PARTICLE_OUTFILE !! Name of output particle information file integer(I4B) :: istep_dump = -1 !! Number of time steps between dumps real(DP) :: rmin = -1.0_DP !! Minimum heliocentric radius for test particle real(DP) :: rmax = -1.0_DP !! Maximum heliocentric radius for test particle @@ -139,15 +140,28 @@ module swiftest_classes procedure :: open => netcdf_open !! Opens a NetCDF file end type netcdf_parameters + !******************************************************************************************************************************** + ! swiftest_swiftest_particle_info class definitions and method interfaces + !******************************************************************************************************************************* + !> Class definition for the particle origin information object. This object is used to track time, location, and collisional regime + !> of fragments produced in collisional events. + type :: swiftest_particle_info + character(len=NAMELEN) :: name !! Non-unique name + character(len=NAMELEN) :: particle_type !! String containing a description of the particle type (e.g. CentralBody, MassiveBody, TestParticle) + contains + procedure :: dump => io_dump_particle_info !! Dumps contents of particle information to file + procedure :: read_in => io_read_in_particle_info !! Read in a particle information object from an open file + end type swiftest_particle_info + !******************************************************************************************************************************** ! swiftest_base class definitions and methods !******************************************************************************************************************************** type, abstract :: swiftest_base - !! An superclass for a generic Swiftest object - logical :: lintegrate = .false. !! Flag indicating that this object should be integrated in the current step + !! An abstract superclass for a generic Swiftest object contains !! The minimal methods that all systems must have - procedure :: dump => io_dump_swiftest !! Dump contents to file + procedure :: dump => io_dump_base !! Dump contents to file + procedure :: dump_particle_info => io_dump_particle_info_base !! Dump contents of particle information metadata to file procedure :: write_frame_netcdf => netcdf_write_frame_base !! I/O routine for writing out a single frame of time-series data for all bodies in the system in NetCDF format generic :: write_frame => write_frame_netcdf !! Set up generic procedure that will switch between NetCDF or Fortran binary depending on arguments end type swiftest_base @@ -157,30 +171,30 @@ module swiftest_classes !******************************************************************************************************************************** !> A concrete lass for the central body in a Swiftest simulation type, abstract, extends(swiftest_base) :: swiftest_cb - character(len=NAMELEN) :: name !! Non-unique name - integer(I4B) :: id = 0 !! External identifier (unique) - real(DP) :: mass = 0.0_DP !! Central body mass (units MU) - real(DP) :: Gmass = 0.0_DP !! Central mass gravitational term G * mass (units GU * MU) - real(DP) :: radius = 0.0_DP !! Central body radius (units DU) - real(DP) :: density = 1.0_DP !! Central body mass density - calculated internally (units MU / DU**3) - real(DP) :: j2rp2 = 0.0_DP !! J2*R^2 term for central body - real(DP) :: j4rp4 = 0.0_DP !! J4*R^2 term for central body - real(DP), dimension(NDIM) :: aobl = 0.0_DP !! Barycentric acceleration due to central body oblatenes - real(DP), dimension(NDIM) :: atide = 0.0_DP !! Barycentric acceleration due to central body oblatenes - real(DP), dimension(NDIM) :: aoblbeg = 0.0_DP !! Barycentric acceleration due to central body oblatenes at beginning of step - real(DP), dimension(NDIM) :: aoblend = 0.0_DP !! Barycentric acceleration due to central body oblatenes at end of step - real(DP), dimension(NDIM) :: atidebeg = 0.0_DP !! Barycentric acceleration due to central body oblatenes at beginning of step - real(DP), dimension(NDIM) :: atideend = 0.0_DP !! Barycentric acceleration due to central body oblatenes at end of step - real(DP), dimension(NDIM) :: xb = 0.0_DP !! Barycentric position (units DU) - real(DP), dimension(NDIM) :: vb = 0.0_DP !! Barycentric velocity (units DU / TU) - real(DP), dimension(NDIM) :: agr = 0.0_DP !! Acceleration due to post-Newtonian correction - real(DP), dimension(NDIM) :: Ip = 0.0_DP !! Unitless principal moments of inertia (I1, I2, I3) / (MR**2). Principal axis rotation assumed. - real(DP), dimension(NDIM) :: rot = 0.0_DP !! Body rotation vector in inertial coordinate frame (units rad / TU) - real(DP) :: k2 = 0.0_DP !! Tidal Love number - real(DP) :: Q = 0.0_DP !! Tidal quality factor - real(DP) :: tlag = 0.0_DP !! Tidal phase lag angle - real(DP), dimension(NDIM) :: L0 = 0.0_DP !! Initial angular momentum of the central body - real(DP), dimension(NDIM) :: dL = 0.0_DP !! Change in angular momentum of the central body + class(swiftest_particle_info), allocatable :: info !! Particle metadata information + integer(I4B) :: id = 0 !! External identifier (unique) + real(DP) :: mass = 0.0_DP !! Central body mass (units MU) + real(DP) :: Gmass = 0.0_DP !! Central mass gravitational term G * mass (units GU * MU) + real(DP) :: radius = 0.0_DP !! Central body radius (units DU) + real(DP) :: density = 1.0_DP !! Central body mass density - calculated internally (units MU / DU**3) + real(DP) :: j2rp2 = 0.0_DP !! J2*R^2 term for central body + real(DP) :: j4rp4 = 0.0_DP !! J4*R^2 term for central body + real(DP), dimension(NDIM) :: aobl = 0.0_DP !! Barycentric acceleration due to central body oblatenes + real(DP), dimension(NDIM) :: atide = 0.0_DP !! Barycentric acceleration due to central body oblatenes + real(DP), dimension(NDIM) :: aoblbeg = 0.0_DP !! Barycentric acceleration due to central body oblatenes at beginning of step + real(DP), dimension(NDIM) :: aoblend = 0.0_DP !! Barycentric acceleration due to central body oblatenes at end of step + real(DP), dimension(NDIM) :: atidebeg = 0.0_DP !! Barycentric acceleration due to central body oblatenes at beginning of step + real(DP), dimension(NDIM) :: atideend = 0.0_DP !! Barycentric acceleration due to central body oblatenes at end of step + real(DP), dimension(NDIM) :: xb = 0.0_DP !! Barycentric position (units DU) + real(DP), dimension(NDIM) :: vb = 0.0_DP !! Barycentric velocity (units DU / TU) + real(DP), dimension(NDIM) :: agr = 0.0_DP !! Acceleration due to post-Newtonian correction + real(DP), dimension(NDIM) :: Ip = 0.0_DP !! Unitless principal moments of inertia (I1, I2, I3) / (MR**2). Principal axis rotation assumed. + real(DP), dimension(NDIM) :: rot = 0.0_DP !! Body rotation vector in inertial coordinate frame (units rad / TU) + real(DP) :: k2 = 0.0_DP !! Tidal Love number + real(DP) :: Q = 0.0_DP !! Tidal quality factor + real(DP) :: tlag = 0.0_DP !! Tidal phase lag angle + real(DP), dimension(NDIM) :: L0 = 0.0_DP !! Initial angular momentum of the central body + real(DP), dimension(NDIM) :: dL = 0.0_DP !! Change in angular momentum of the central body contains procedure :: read_in => io_read_in_cb !! I/O routine for reading in central body data procedure :: read_frame => io_read_frame_cb !! I/O routine for reading out a single frame of time-series data for the central body @@ -196,7 +210,7 @@ module swiftest_classes !! Superclass that defines the generic elements of a Swiftest particle logical :: lfirst = .true. !! Run the current step as a first integer(I4B) :: nbody = 0 !! Number of bodies - character(len=NAMELEN), dimension(:), allocatable :: name !! Non-unique name + class(swiftest_particle_info), dimension(:), allocatable :: info !! Particle metadata information integer(I4B), dimension(:), allocatable :: id !! External identifier (unique) integer(I4B), dimension(:), allocatable :: status !! An integrator-specific status indicator logical, dimension(:), allocatable :: ldiscard !! Body should be discarded @@ -357,22 +371,24 @@ module swiftest_classes procedure(abstract_step_system), deferred :: step ! Concrete classes that are common to the basic integrator (only test particles considered for discard) - procedure :: discard => discard_system !! Perform a discard step on the system - procedure :: conservation_report => io_conservation_report !! Compute energy and momentum and print out the change with time - procedure :: dump => io_dump_system !! Dump the state of the system to a file - procedure :: get_old_t_final => io_get_old_t_final_system !! Validates the dump file to check whether the dump file initial conditions duplicate the last frame of the binary output. - procedure :: read_frame => io_read_frame_system !! Read in a frame of input data from file - procedure :: write_discard => io_write_discard !! Write out information about discarded test particles - procedure :: write_frame => io_write_frame_system !! Append a frame of output data to file - procedure :: write_hdr_bin => io_write_hdr_system !! Write a header for an output frame in Fortran binary format - procedure :: write_hdr_netcdf => netcdf_write_hdr_system !! Write a header for an output frame in NetCDF format - procedure :: initialize => setup_initialize_system !! Initialize the system from input files - procedure :: step_spin => tides_step_spin_system !! Steps the spins of the massive & central bodies due to tides. - procedure :: set_msys => util_set_msys !! Sets the value of msys from the masses of system bodies. - procedure :: get_energy_and_momentum => util_get_energy_momentum_system !! Calculates the total system energy and momentum - procedure :: rescale => util_rescale_system !! Rescales the system into a new set of units - procedure :: validate_ids => util_valid_id_system !! Validate the numerical ids passed to the system and save the maximum value - generic :: write_hdr => write_hdr_bin, write_hdr_netcdf !! Generic method call for writing headers + procedure :: discard => discard_system !! Perform a discard step on the system + procedure :: conservation_report => io_conservation_report !! Compute energy and momentum and print out the change with time + procedure :: dump => io_dump_system !! Dump the state of the system to a file + procedure :: get_old_t_final => io_get_old_t_final_system !! Validates the dump file to check whether the dump file initial conditions duplicate the last frame of the binary output. + procedure :: read_frame => io_read_frame_system !! Read in a frame of input data from file + procedure :: read_particle_info => io_read_particle_info_system !! Read in particle metadata from file + procedure :: write_discard => io_write_discard !! Write out information about discarded test particles + procedure :: write_frame => io_write_frame_system !! Append a frame of output data to file + procedure :: write_hdr_bin => io_write_hdr_system !! Write a header for an output frame in Fortran binary format + procedure :: write_hdr_netcdf => netcdf_write_hdr_system !! Write a header for an output frame in NetCDF format + procedure :: initialize => setup_initialize_system !! Initialize the system from input files + procedure :: init_particle_info => setup_initialize_particle_info_system !! Initialize the system from input files + procedure :: step_spin => tides_step_spin_system !! Steps the spins of the massive & central bodies due to tides. + procedure :: set_msys => util_set_msys !! Sets the value of msys from the masses of system bodies. + procedure :: get_energy_and_momentum => util_get_energy_momentum_system !! Calculates the total system energy and momentum + procedure :: rescale => util_rescale_system !! Rescales the system into a new set of units + procedure :: validate_ids => util_valid_id_system !! Validate the numerical ids passed to the system and save the maximum value + generic :: write_hdr => write_hdr_bin, write_hdr_netcdf !! Generic method call for writing headers end type swiftest_nbody_system type :: swiftest_encounter @@ -614,11 +630,24 @@ module subroutine io_dump_param(self, param_file_name) character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) end subroutine io_dump_param - module subroutine io_dump_swiftest(self, param) + module subroutine io_dump_particle_info_base(self, param, idx) + implicit none + class(swiftest_base), intent(inout) :: self !! Swiftest base object (can be cb, pl, or tp) + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + integer(I4B), dimension(:), optional, intent(in) :: idx !! Array of test particle indices to append to the particle file + end subroutine io_dump_particle_info_base + + module subroutine io_dump_particle_info(self, iu) + implicit none + class(swiftest_particle_info), intent(in) :: self !! Swiftest particle info metadata object + integer(I4B), intent(in) :: iu !! Open unformatted file unit number + end subroutine io_dump_particle_info + + module subroutine io_dump_base(self, param) implicit none class(swiftest_base), intent(inout) :: self !! Swiftest base object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine io_dump_swiftest + end subroutine io_dump_base module subroutine io_dump_system(self, param) implicit none @@ -689,6 +718,12 @@ module subroutine io_read_in_param(self, param_file_name) character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) end subroutine io_read_in_param + module subroutine io_read_in_particle_info(self, iu) + implicit none + class(swiftest_particle_info), intent(inout) :: self !! Particle metadata information object + integer(I4B), intent(in) :: iu !! Open file unit number + end subroutine io_read_in_particle_info + module function io_read_frame_body(self, iu, param) result(ierr) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest body object @@ -713,6 +748,12 @@ module function io_read_frame_system(self, iu, param) result(ierr) integer(I4B) :: ierr !! Error code: returns 0 if the read is successful end function io_read_frame_system + module subroutine io_read_particle_info_system(self, param) + implicit none + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine io_read_particle_info_system + module subroutine io_write_discard(self, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object @@ -953,6 +994,12 @@ module subroutine setup_encounter(self, n) integer(I4B), intent(in) :: n !! Number of encounters to allocate space for end subroutine setup_encounter + module subroutine setup_initialize_particle_info_system(self, param) + implicit none + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + end subroutine setup_initialize_particle_info_system + module subroutine setup_initialize_system(self, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object @@ -1030,6 +1077,14 @@ module subroutine util_append_arr_I4B(arr, source, nold, nsrc, lsource_mask) logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_arr_I4B + module subroutine util_append_arr_info(arr, source, nold, nsrc, lsource_mask) + implicit none + class(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: arr !! Destination array + class(swiftest_particle_info), dimension(:), allocatable, intent(in) :: source !! Array to append + integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + end subroutine util_append_arr_info + module subroutine util_append_arr_logical(arr, source, nold, nsrc, lsource_mask) implicit none logical, dimension(:), allocatable, intent(inout) :: arr !! Destination array @@ -1037,6 +1092,7 @@ module subroutine util_append_arr_logical(arr, source, nold, nsrc, lsource_mask) integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_arr_logical + end interface interface @@ -1190,6 +1246,13 @@ module subroutine util_fill_arr_I4B(keeps, inserts, lfill_list) logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps end subroutine util_fill_arr_I4B + module subroutine util_fill_arr_info(keeps, inserts, lfill_list) + implicit none + class(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: keeps !! Array of values to keep + class(swiftest_particle_info), dimension(:), allocatable, intent(in) :: inserts !! Array of values to insert into keep + logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps + end subroutine util_fill_arr_info + module subroutine util_fill_arr_logical(keeps, inserts, lfill_list) implicit none logical, dimension(:), allocatable, intent(inout) :: keeps !! Array of values to keep @@ -1252,6 +1315,12 @@ module subroutine util_resize_arr_I4B(arr, nnew) integer(I4B), intent(in) :: nnew !! New size end subroutine util_resize_arr_I4B + module subroutine util_resize_arr_info(arr, nnew) + implicit none + class(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: arr !! Array to resize + integer(I4B), intent(in) :: nnew !! New size + end subroutine util_resize_arr_info + module subroutine util_resize_arr_logical(arr, nnew) implicit none logical, dimension(:), allocatable, intent(inout) :: arr !! Array to resize @@ -1430,6 +1499,13 @@ module subroutine util_sort_rearrange_arr_I4B(arr, ind, n) integer(I4B), intent(in) :: n !! Number of elements in arr and ind to rearrange end subroutine util_sort_rearrange_arr_I4B + module subroutine util_sort_rearrange_arr_info(arr, ind, n) + implicit none + class(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: arr !! Destination array + integer(I4B), dimension(:), intent(in) :: ind !! Index to rearrange against + integer(I4B), intent(in) :: n !! Number of elements in arr and ind to rearrange + end subroutine util_sort_rearrange_arr_info + module subroutine util_sort_rearrange_arr_logical(arr, ind, n) implicit none logical, dimension(:), allocatable, intent(inout) :: arr !! Destination array @@ -1520,6 +1596,14 @@ module subroutine util_spill_arr_I8B(keeps, discards, lspill_list, ldestructive) logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not end subroutine util_spill_arr_I8B + module subroutine util_spill_arr_info(keeps, discards, lspill_list, ldestructive) + implicit none + class(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: keeps !! Array of values to keep + class(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: discards !! Array of discards + logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discardss + logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not + end subroutine util_spill_arr_info + module subroutine util_spill_arr_logical(keeps, discards, lspill_list, ldestructive) implicit none logical, dimension(:), allocatable, intent(inout) :: keeps !! Array of values to keep diff --git a/src/modules/swiftest_globals.f90 b/src/modules/swiftest_globals.f90 index 3df99a8d1..dd6752b8f 100644 --- a/src/modules/swiftest_globals.f90 +++ b/src/modules/swiftest_globals.f90 @@ -55,6 +55,10 @@ module swiftest_globals character(*), parameter :: XV = 'XV' !! Symbolic name for binary output file contents for cartesian position and velocity vectors character(*), parameter :: XVEL = 'XVEL' !! Symbolic name for binary output file contents for both cartesian position and velocity and orbital elements + character(*), parameter :: CB_TYPE_NAME = "CentralBody" + character(*), parameter :: PL_TYPE_NAME = "MassiveBody" + character(*), parameter :: TP_TYPE_NAME = "TestParticle" + ! OpenMP Parameters integer(I4B) :: nthreads = 1 !! Number of OpenMP threads integer(I4B), parameter :: NTHERSHOLD = 1000 !! Threshold value for OpenMP loop parallelization @@ -113,11 +117,14 @@ module swiftest_globals character(*), dimension(2), parameter :: DUMP_PARAM_FILE = ['dump_param1.in', 'dump_param2.in'] !> Default file names that can be changed by the user in the parameters file - character(*), parameter :: CB_INFILE = 'cb.in' - character(*), parameter :: PL_INFILE = 'pl.in' - character(*), parameter :: TP_INFILE = 'tp.in' - character(*), parameter :: BIN_OUTFILE = 'bin.dat' - integer(I4B), parameter :: BINUNIT = 20 !! File unit number for the binary output file + character(*), parameter :: CB_INFILE = 'cb.in' + character(*), parameter :: PL_INFILE = 'pl.in' + character(*), parameter :: TP_INFILE = 'tp.in' + character(*), parameter :: BIN_OUTFILE = 'bin.dat' + integer(I4B), parameter :: BINUNIT = 20 !! File unit number for the binary output file + character(*), parameter :: PARTICLE_OUTFILE = 'particle.dat' + integer(I4B), parameter :: PARTICLEUNIT = 44 !! File unit number for the binary particle info output file + !> Miscellaneous constants: integer(I4B), parameter :: NDIM = 3 !! Number of dimensions in our reality diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 25f9b5940..5294efa64 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -4,7 +4,7 @@ module symba_classes !! Definition of classes and methods specific to the Democratic SyMBAcentric Method !! Adapted from David E. Kaufmann's Swifter routine: helio.f90 use swiftest_globals - use swiftest_classes, only : swiftest_parameters, swiftest_base, swiftest_encounter + use swiftest_classes, only : swiftest_parameters, swiftest_base, swiftest_encounter, swiftest_particle_info use helio_classes, only : helio_cb, helio_pl, helio_tp, helio_nbody_system use rmvs_classes, only : rmvs_chk_ind implicit none @@ -14,11 +14,8 @@ module symba_classes integer(I4B), private, parameter :: NTENC = 3 real(DP), private, parameter :: RHSCALE = 6.5_DP real(DP), private, parameter :: RSHELL = 0.48075_DP - character(*), parameter :: PARTICLE_OUTFILE = 'particle.dat' - integer(I4B), parameter :: PARTICLEUNIT = 44 !! File unit number for the binary particle info output file type, extends(swiftest_parameters) :: symba_parameters - character(STRMAX) :: particle_out = PARTICLE_OUTFILE !! Name of output particle information file real(DP) :: GMTINY = -1.0_DP !! Smallest G*mass that is fully gravitating real(DP) :: min_GMfrag = -1.0_DP !! Smallest G*mass that can be produced in a fragmentation event integer(I4B), dimension(:), allocatable :: seed !! Random seeds @@ -29,16 +26,17 @@ module symba_classes end type symba_parameters !******************************************************************************************************************************** - ! symba_particle_info class definitions and method interfaces + ! symba_swiftest_particle_info class definitions and method interfaces !******************************************************************************************************************************* !> Class definition for the particle origin information object. This object is used to track time, location, and collisional regime !> of fragments produced in collisional events. - type :: symba_particle_info - sequence - character(len=32) :: origin_type !! String containing a description of the origin of the particle (e.g. Initial Conditions, Supercatastrophic, Disruption, etc.) + type, extends(swiftest_particle_info) :: symba_particle_info + character(len=NAMELEN) :: origin_type !! String containing a description of the origin of the particle (e.g. Initial Conditions, Supercatastrophic, Disruption, etc.) real(DP) :: origin_time !! The time of the particle's formation real(DP), dimension(NDIM) :: origin_xh !! The heliocentric distance vector at the time of the particle's formation real(DP), dimension(NDIM) :: origin_vh !! The heliocentric velocity vector at the time of the particle's formation + contains + procedure :: read_in => symba_io_read_in_particle_info !! Reads in SyMBA particle information metadata from an open unformatted file end type symba_particle_info !******************************************************************************************************************************** @@ -60,7 +58,6 @@ module symba_classes real(DP) :: dGM = 0.0_DP !! Change in G*mass of the central body real(DP) :: R0 = 0.0_DP !! Initial radius of the central body real(DP) :: dR = 0.0_DP !! Change in the radius of the central body - type(symba_particle_info) :: info contains end type symba_cb @@ -82,7 +79,6 @@ module symba_classes real(DP), dimension(:), allocatable :: peri !! perihelion distance real(DP), dimension(:), allocatable :: atp !! semimajor axis following perihelion passage type(symba_kinship), dimension(:), allocatable :: kin !! Array of merger relationship structures that can account for multiple pairwise mergers in a single step - type(symba_particle_info), dimension(:), allocatable :: info contains procedure :: make_family => symba_collision_make_family_pl !! When a single body is involved in more than one collision in a single step, it becomes part of a family procedure :: index => symba_util_index_eucl_plpl !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix @@ -118,7 +114,6 @@ module symba_classes integer(I4B), dimension(:), allocatable :: nplenc !! number of encounters with planets this time step integer(I4B), dimension(:), allocatable :: levelg !! level at which this particle should be moved integer(I4B), dimension(:), allocatable :: levelm !! deepest encounter level achieved this time step - type(symba_particle_info), dimension(:), allocatable :: info contains procedure :: drift => symba_drift_tp !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level procedure :: encounter_check => symba_encounter_check_tp !! Checks if any test particles are undergoing a close encounter with a massive body @@ -173,19 +168,20 @@ module symba_classes ! symba_nbody_system class definitions and method interfaces !******************************************************************************************************************************** type, extends(helio_nbody_system) :: symba_nbody_system - class(symba_merger), allocatable :: pl_adds !! List of added bodies in mergers or collisions - class(symba_pltpenc), allocatable :: pltpenc_list !! List of massive body-test particle encounters in a single step - class(symba_plplenc), allocatable :: plplenc_list !! List of massive body-massive body encounters in a single step - class(symba_plplenc), allocatable :: plplcollision_list !! List of massive body-massive body collisions in a single step - integer(I4B) :: irec !! System recursion level + class(symba_merger), allocatable :: pl_adds !! List of added bodies in mergers or collisions + class(symba_pltpenc), allocatable :: pltpenc_list !! List of massive body-test particle encounters in a single step + class(symba_plplenc), allocatable :: plplenc_list !! List of massive body-massive body encounters in a single step + class(symba_plplenc), allocatable :: plplcollision_list !! List of massive body-massive body collisions in a single step + integer(I4B) :: irec !! System recursion level contains - procedure :: write_discard => symba_io_write_discard !! Write out information about discarded and merged planets and test particles in SyMBA - procedure :: initialize => symba_setup_initialize_system !! Performs SyMBA-specific initilization steps - procedure :: step => symba_step_system !! Advance the SyMBA nbody system forward in time by one step - procedure :: interp => symba_step_interp_system !! Perform an interpolation step on the SymBA nbody system - procedure :: set_recur_levels => symba_step_set_recur_levels_system !! Sets recursion levels of bodies and encounter lists to the current system level - procedure :: recursive_step => symba_step_recur_system !! Step interacting planets and active test particles ahead in democratic heliocentric coordinates at the current recursion level, if applicable, and descend to the next deeper level if necessary - procedure :: reset => symba_step_reset_system !! Resets pl, tp,and encounter structures at the start of a new step + procedure :: write_discard => symba_io_write_discard !! Write out information about discarded and merged planets and test particles in SyMBA + procedure :: initialize => symba_setup_initialize_system !! Performs SyMBA-specific initilization steps + procedure :: init_particle_info => symba_setup_initialize_particle_info_system !! Initialize the system from input files + procedure :: step => symba_step_system !! Advance the SyMBA nbody system forward in time by one step + procedure :: interp => symba_step_interp_system !! Perform an interpolation step on the SymBA nbody system + procedure :: set_recur_levels => symba_step_set_recur_levels_system !! Sets recursion levels of bodies and encounter lists to the current system level + procedure :: recursive_step => symba_step_recur_system !! Step interacting planets and active test particles ahead in democratic heliocentric coordinates at the current recursion level, if applicable, and descend to the next deeper level if necessary + procedure :: reset => symba_step_reset_system !! Resets pl, tp,and encounter structures at the start of a new step end type symba_nbody_system interface @@ -362,21 +358,12 @@ module subroutine symba_util_index_eucl_plpl(self, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine symba_util_index_eucl_plpl - module subroutine symba_io_write_discard(self, param) - use swiftest_classes, only : swiftest_parameters - implicit none - class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine symba_io_write_discard - - module subroutine symba_io_dump_particle_info(system, param, lincludecb, tpidx, plidx) + module subroutine symba_io_dump_particle_info(self, iu) implicit none - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA extensions - logical, optional, intent(in) :: lincludecb !! Set to true to include the central body (default is false) - integer(I4B), dimension(:), optional, intent(in) :: tpidx !! Array of test particle indices to append to the particle file - integer(I4B), dimension(:), optional, intent(in) :: plidx !! Array of massive body indices to append to the particle file + class(symba_particle_info), intent(in) :: self !! Particle metadata information object + integer(I4B), intent(in) :: iu !! Open file unit number end subroutine symba_io_dump_particle_info + module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) implicit none @@ -400,11 +387,18 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 end subroutine symba_io_param_writer - module subroutine symba_io_read_particle(system, param) + module subroutine symba_io_read_in_particle_info(self, iu) implicit none - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system file - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA extensions - end subroutine symba_io_read_particle + class(symba_particle_info), intent(inout) :: self !! Particle metadata information object + integer(I4B), intent(in) :: iu !! Open file unit number + end subroutine symba_io_read_in_particle_info + + module subroutine symba_io_write_discard(self, param) + use swiftest_classes, only : swiftest_parameters + implicit none + class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + end subroutine symba_io_write_discard module subroutine symba_kick_getacch_int_pl(self) implicit none @@ -440,11 +434,12 @@ module subroutine symba_kick_encounter(self, system, dt, irec, sgn) integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration end subroutine symba_kick_encounter - module subroutine symba_setup_initialize_particle_info(system, param) + module subroutine symba_setup_initialize_particle_info_system(self, param) + use swiftest_classes, only : swiftest_parameters implicit none - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA extensions - end subroutine symba_setup_initialize_particle_info + class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA extensions + end subroutine symba_setup_initialize_particle_info_system module subroutine symba_setup_initialize_system(self, param) use swiftest_classes, only : swiftest_parameters @@ -524,14 +519,6 @@ end subroutine symba_step_reset_system end interface interface util_append - module subroutine symba_util_append_arr_info(arr, source, nold, nsrc, lsource_mask) - implicit none - type(symba_particle_info), dimension(:), allocatable, intent(inout) :: arr !! Destination array - type(symba_particle_info), dimension(:), allocatable, intent(in) :: source !! Array to append - integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively - logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to - end subroutine symba_util_append_arr_info - module subroutine symba_util_append_arr_kin(arr, source, nold, nsrc, lsource_mask) implicit none type(symba_kinship), dimension(:), allocatable, intent(inout) :: arr !! Destination array @@ -582,13 +569,6 @@ end subroutine symba_util_copy_encounter end interface interface util_fill - module subroutine symba_util_fill_arr_info(keeps, inserts, lfill_list) - implicit none - type(symba_particle_info), dimension(:), allocatable, intent(inout) :: keeps !! Array of values to keep - type(symba_particle_info), dimension(:), allocatable, intent(in) :: inserts !! Array of values to insert into keep - logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps - end subroutine symba_util_fill_arr_info - module subroutine symba_util_fill_arr_kin(keeps, inserts, lfill_list) implicit none type(symba_kinship), dimension(:), allocatable, intent(inout) :: keeps !! Array of values to keep @@ -631,12 +611,6 @@ end subroutine symba_util_rearray_pl end interface interface util_resize - module subroutine symba_util_resize_arr_info(arr, nnew) - implicit none - type(symba_particle_info), dimension(:), allocatable, intent(inout) :: arr !! Array to resize - integer(I4B), intent(in) :: nnew !! New size - end subroutine symba_util_resize_arr_info - module subroutine symba_util_resize_arr_kin(arr, nnew) implicit none type(symba_kinship), dimension(:), allocatable, intent(inout) :: arr !! Array to resize @@ -679,13 +653,6 @@ end subroutine symba_util_sort_tp end interface interface util_sort_rearrange - module subroutine symba_util_sort_rearrange_arr_info(arr, ind, n) - implicit none - type(symba_particle_info), dimension(:), allocatable, intent(inout) :: arr !! Destination array - integer(I4B), dimension(:), intent(in) :: ind !! Index to rearrange against - integer(I4B), intent(in) :: n !! Number of elements in arr and ind to rearrange - end subroutine symba_util_sort_rearrange_arr_info - module subroutine symba_util_sort_rearrange_arr_kin(arr, ind, n) implicit none type(symba_kinship), dimension(:), allocatable, intent(inout) :: arr !! Destination array @@ -709,14 +676,6 @@ end subroutine symba_util_sort_rearrange_tp end interface interface util_spill - module subroutine symba_util_spill_arr_info(keeps, discards, lspill_list, ldestructive) - implicit none - type(symba_particle_info), dimension(:), allocatable, intent(inout) :: keeps !! Array of values to keep - type(symba_particle_info), dimension(:), allocatable, intent(inout) :: discards !! Array of discards - logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discardss - logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not - end subroutine symba_util_spill_arr_info - module subroutine symba_util_spill_arr_kin(keeps, discards, lspill_list, ldestructive) implicit none type(symba_kinship), dimension(:), allocatable, intent(inout) :: keeps !! Array of values to keep diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index 2abcbead1..785e10a92 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -231,7 +231,7 @@ module subroutine netcdf_write_frame_base(self, iu, param) j = ind(i) idslot = self%id(j) + 1 call check( nf90_put_var(iu%ncid, iu%id_varid, self%id(j), start=[idslot]) ) - name = trim(adjustl(self%name(j))) + name = trim(adjustl(self%info(j)%name)) strlen = len(name) call check( nf90_put_var(iu%ncid, iu%name_varid, name, start=[1, idslot], count=[strlen, 1]) ) if ((param%out_form == XV) .or. (param%out_form == XVEL)) then @@ -277,7 +277,7 @@ module subroutine netcdf_write_frame_base(self, iu, param) class is (swiftest_cb) idslot = self%id + 1 call check( nf90_put_var(iu%ncid, iu%id_varid, self%id, start=[idslot]) ) - name = trim(adjustl(self%name)) + name = trim(adjustl(self%info%name)) strlen = len(name) call check( nf90_put_var(iu%ncid, iu%name_varid, name, start=[1, idslot], count=[strlen, 1]) ) call check( nf90_put_var(iu%ncid, iu%Gmass_varid, self%Gmass, start=[idslot, tslot]) ) diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 070308838..fe78ee9e3 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -67,6 +67,13 @@ module subroutine setup_construct_system(system, param) call util_exit(FAILURE) end select + select type(system) + class is (symba_nbody_system) + allocate(symba_particle_info :: system%cb%info) + class default + allocate(swiftest_particle_info :: system%cb%info) + end select + return end subroutine setup_construct_system @@ -127,6 +134,33 @@ module subroutine setup_encounter(self, n) end subroutine setup_encounter + module subroutine setup_initialize_particle_info_system(self, param) + !! author: David A. Minton + !! + !! Setup up particle information metadata from initial conditions + ! + implicit none + ! Arguments + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + + associate(cb => self%cb, pl => self%pl, npl => self%pl%nbody, tp => self%tp, ntp => self%tp%nbody) + cb%info%particle_type = CB_TYPE_NAME + call cb%dump_particle_info(param) + if (npl > 0) then + pl%info(1:npl)%particle_type = PL_TYPE_NAME + call pl%dump_particle_info(param) + end if + if (ntp > 0) then + tp%info(1:ntp)%particle_type = TP_TYPE_NAME + call tp%dump_particle_info(param) + end if + end associate + + return + end subroutine setup_initialize_particle_info_system + + module subroutine setup_initialize_system(self, param) !! author: David A. Minton !! @@ -152,6 +186,12 @@ module subroutine setup_initialize_system(self, param) if (.not.param%lrhill_present) call self%pl%set_rhill(self%cb) self%pl%lfirst = param%lfirstkick self%tp%lfirst = param%lfirstkick + + if (param%lrestart) then + call self%read_particle_info(param) + else + call self%init_particle_info(param) + end if return end subroutine setup_initialize_system @@ -173,7 +213,7 @@ module subroutine setup_body(self, n, param) self%lfirst = .true. if (allocated(self%id)) deallocate(self%id) - if (allocated(self%name)) deallocate(self%name) + if (allocated(self%info)) deallocate(self%info) if (allocated(self%status)) deallocate(self%status) if (allocated(self%ldiscard)) deallocate(self%ldiscard) if (allocated(self%xh)) deallocate(self%xh) @@ -186,7 +226,7 @@ module subroutine setup_body(self, n, param) if (allocated(self%lmask)) deallocate(self%lmask) allocate(self%id(n)) - allocate(self%name(n)) + allocate(self%info(n)) allocate(self%status(n)) allocate(self%ldiscard(n)) allocate(self%xh(NDIM, n)) @@ -199,7 +239,7 @@ module subroutine setup_body(self, n, param) allocate(self%lmask(n)) self%id(:) = 0 - self%name(:) = "UNNAMED" + self%info(:)%name = "UNNAMED" self%status(:) = INACTIVE self%lmask(:) = .false. self%ldiscard(:) = .false. @@ -247,10 +287,12 @@ module subroutine setup_pl(self, n, param) call setup_body(self, n, param) if (n <= 0) return + if (allocated(self%info)) deallocate(self%info) if (allocated(self%mass)) deallocate(self%mass) if (allocated(self%Gmass)) deallocate(self%Gmass) if (allocated(self%rhill)) deallocate(self%rhill) + allocate(swiftest_particle_info :: self%info(n)) allocate(self%mass(n)) allocate(self%Gmass(n)) allocate(self%rhill(n)) @@ -311,14 +353,17 @@ module subroutine setup_tp(self, n, param) call setup_body(self, n, param) if (n <= 0) return + if (allocated(self%info)) deallocate(self%info) if (allocated(self%isperi)) deallocate(self%isperi) if (allocated(self%peri)) deallocate(self%peri) if (allocated(self%atp)) deallocate(self%atp) + allocate(swiftest_particle_info :: self%info(n)) allocate(self%isperi(n)) allocate(self%peri(n)) allocate(self%atp(n)) + self%info(:)%particle_type = TP_TYPE_NAME self%isperi(:) = 0 self%peri(:) = 0.0_DP self%atp(:) = 0.0_DP diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 034038ded..4e21b9411 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -838,101 +838,117 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, class is (symba_pl) select type(pl_discards => system%pl_discards) class is (symba_merger) - associate(pl_adds => system%pl_adds, cb => system%cb) - - ! Add the family bodies to the subtraction list - nfamily = size(family(:)) - nfrag = size(m_frag(:)) - lmask(:) = .false. - lmask(family(:)) = .true. - pl%status(family(:)) = MERGED - nstart = pl_discards%nbody + 1 - nend = pl_discards%nbody + nfamily - call pl_discards%append(pl, lmask) - pl%ldiscard(family(:)) = .true. - pl%lcollision(family(:)) = .true. - - ! Record how many bodies were subtracted in this event - pl_discards%ncomp(nstart:nend) = nfamily - - ! Setup new bodies - allocate(plnew, mold=pl) - call plnew%setup(nfrag, param) - ibiggest = family(maxloc(pl%Gmass(family(:)), dim=1)) - - ! Copy over identification, information, and physical properties of the new bodies from the fragment list - plnew%id(1:nfrag) = id_frag(1:nfrag) - param%maxid = param%maxid + nfrag - plnew%xb(:, 1:nfrag) = xb_frag(:, 1:nfrag) - plnew%vb(:, 1:nfrag) = vb_frag(:, 1:nfrag) - do i = 1, nfrag - plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) - plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) - end do - plnew%mass(1:nfrag) = m_frag(1:nfrag) - plnew%Gmass(1:nfrag) = param%GU * m_frag(1:nfrag) - plnew%radius(1:nfrag) = rad_frag(1:nfrag) - plnew%density(1:nfrag) = m_frag(1:nfrag) / rad_frag(1:nfrag) + select type(info => pl%info) + class is (symba_particle_info) + associate(pl_adds => system%pl_adds, cb => system%cb) + + ! Add the family bodies to the subtraction list + nfamily = size(family(:)) + nfrag = size(m_frag(:)) + lmask(:) = .false. + lmask(family(:)) = .true. + pl%status(family(:)) = MERGED + nstart = pl_discards%nbody + 1 + nend = pl_discards%nbody + nfamily + call pl_discards%append(pl, lmask) + pl%ldiscard(family(:)) = .true. + pl%lcollision(family(:)) = .true. + + ! Record how many bodies were subtracted in this event + pl_discards%ncomp(nstart:nend) = nfamily + + ! Setup new bodies + allocate(plnew, mold=pl) + call plnew%setup(nfrag, param) + ibiggest = family(maxloc(pl%Gmass(family(:)), dim=1)) - select case(status) - case(DISRUPTION) - plnew%info(1:nfrag)%origin_type = "Disruption" - plnew%status(1:nfrag) = NEW_PARTICLE - plnew%info(1:nfrag)%origin_time = param%t + ! Copy over identification, information, and physical properties of the new bodies from the fragment list + plnew%id(1:nfrag) = id_frag(1:nfrag) + param%maxid = param%maxid + nfrag + plnew%xb(:, 1:nfrag) = xb_frag(:, 1:nfrag) + plnew%vb(:, 1:nfrag) = vb_frag(:, 1:nfrag) do i = 1, nfrag - plnew%info(i)%origin_xh(:) = plnew%xh(:,i) - plnew%info(i)%origin_vh(:) = plnew%vh(:,i) + plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) + plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) end do - case(SUPERCATASTROPHIC) - plnew%info(1:nfrag)%origin_type = "Supercatastrophic" - plnew%status(1:nfrag) = NEW_PARTICLE - plnew%info(1:nfrag)%origin_time = param%t - do i = 1, nfrag - plnew%info(i)%origin_xh(:) = plnew%xh(:,i) - plnew%info(i)%origin_vh(:) = plnew%vh(:,i) - end do - case(HIT_AND_RUN_DISRUPT) - plnew%info(1) = pl%info(ibiggest) - plnew%status(1) = OLD_PARTICLE - plnew%status(2:nfrag) = NEW_PARTICLE - plnew%info(2:nfrag)%origin_type = "Hit and run fragment" - plnew%info(2:nfrag)%origin_time = param%t - do i = 2, nfrag - plnew%info(i)%origin_xh(:) = plnew%xh(:,i) - plnew%info(i)%origin_vh(:) = plnew%vh(:,i) - end do - case(MERGED) - plnew%info(1) = pl%info(ibiggest) - plnew%status(1) = OLD_PARTICLE - end select - - if (param%lrotation) then - plnew%Ip(:, 1:nfrag) = Ip_frag(:, 1:nfrag) - plnew%rot(:, 1:nfrag) = rot_frag(:, 1:nfrag) - end if + plnew%mass(1:nfrag) = m_frag(1:nfrag) + plnew%Gmass(1:nfrag) = param%GU * m_frag(1:nfrag) + plnew%radius(1:nfrag) = rad_frag(1:nfrag) + plnew%density(1:nfrag) = m_frag(1:nfrag) / rad_frag(1:nfrag) - if (param%ltides) then - plnew%Q = pl%Q(ibiggest) - plnew%k2 = pl%k2(ibiggest) - plnew%tlag = pl%tlag(ibiggest) - end if + select type(info => plnew%info) + class is (symba_particle_info) + select case(status) + case(DISRUPTION) + info(1:nfrag)%origin_type = "Disruption" + plnew%status(1:nfrag) = NEW_PARTICLE + info(1:nfrag)%origin_time = param%t + do i = 1, nfrag + info(i)%origin_xh(:) = plnew%xh(:,i) + info(i)%origin_vh(:) = plnew%vh(:,i) + end do + case(SUPERCATASTROPHIC) + info(1:nfrag)%origin_type = "Supercatastrophic" + plnew%status(1:nfrag) = NEW_PARTICLE + info(1:nfrag)%origin_time = param%t + do i = 1, nfrag + info(i)%origin_xh(:) = plnew%xh(:,i) + info(i)%origin_vh(:) = plnew%vh(:,i) + end do + case(HIT_AND_RUN_DISRUPT) + select type(plinfo => pl%info) + class is (symba_particle_info) + info(1)%name = plinfo(ibiggest)%name + info(1)%origin_xh(:) = plinfo(ibiggest)%origin_xh(:) + info(1)%origin_vh(:) = plinfo(ibiggest)%origin_vh(:) + end select + plnew%status(1) = OLD_PARTICLE + plnew%status(2:nfrag) = NEW_PARTICLE + info(2:nfrag)%origin_type = "Hit and run fragment" + info(2:nfrag)%origin_time = param%t + do i = 2, nfrag + info(i)%origin_xh(:) = plnew%xh(:,i) + info(i)%origin_vh(:) = plnew%vh(:,i) + end do + case(MERGED) + select type(plinfo => pl%info) + class is (symba_particle_info) + info(1)%name = plinfo(ibiggest)%name + info(1)%origin_xh(:) = plinfo(ibiggest)%origin_xh(:) + info(1)%origin_vh(:) = plinfo(ibiggest)%origin_vh(:) + end select + plnew%status(1) = OLD_PARTICLE + end select + end select + + if (param%lrotation) then + plnew%Ip(:, 1:nfrag) = Ip_frag(:, 1:nfrag) + plnew%rot(:, 1:nfrag) = rot_frag(:, 1:nfrag) + end if + + if (param%ltides) then + plnew%Q = pl%Q(ibiggest) + plnew%k2 = pl%k2(ibiggest) + plnew%tlag = pl%tlag(ibiggest) + end if - call plnew%set_mu(cb) - !Copy over or set integration parameters for new bodies - plnew%lcollision(1:nfrag) = .false. - plnew%ldiscard(1:nfrag) = .false. - plnew%levelg(1:nfrag) = pl%levelg(ibiggest) - plnew%levelm(1:nfrag) = pl%levelm(ibiggest) - - ! Append the new merged body to the list and record how many we made - nstart = pl_adds%nbody + 1 - nend = pl_adds%nbody + plnew%nbody - call pl_adds%append(plnew, lsource_mask=[(.true., i=1, plnew%nbody)]) - pl_adds%ncomp(nstart:nend) = plnew%nbody - - call plnew%setup(0, param) - deallocate(plnew) - end associate + call plnew%set_mu(cb) + !Copy over or set integration parameters for new bodies + plnew%lcollision(1:nfrag) = .false. + plnew%ldiscard(1:nfrag) = .false. + plnew%levelg(1:nfrag) = pl%levelg(ibiggest) + plnew%levelm(1:nfrag) = pl%levelm(ibiggest) + + ! Append the new merged body to the list and record how many we made + nstart = pl_adds%nbody + 1 + nend = pl_adds%nbody + plnew%nbody + call pl_adds%append(plnew, lsource_mask=[(.true., i=1, plnew%nbody)]) + pl_adds%ncomp(nstart:nend) = plnew%nbody + + call plnew%setup(0, param) + deallocate(plnew) + end associate + end select end select end select diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 542db70c9..281fe7cdb 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -2,76 +2,26 @@ use swiftest contains - module subroutine symba_io_dump_particle_info(system, param, lincludecb, tpidx, plidx) + module subroutine symba_io_dump_particle_info(self, iu) !! author: David A. Minton !! - !! Dumps the particle information data to a file. - !! Pass a list of array indices for test particles (tpidx) and/or massive bodies (plidx) to append + !! Reads in particle information object information from an open file unformatted file implicit none ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA extensions - logical, optional, intent(in) :: lincludecb !! Set to true to include the central body (default is false) - integer(I4B), dimension(:), optional, intent(in) :: tpidx !! Array of test particle indices to append to the particle file - integer(I4B), dimension(:), optional, intent(in) :: plidx !! Array of massive body indices to append to the particle file + class(symba_particle_info), intent(in) :: self !! Particle metadata information object + integer(I4B), intent(in) :: iu !! Open file unit number ! Internals - logical, save :: lfirst = .true. - integer(I4B), parameter :: LUN = 22 - integer(I4B) :: i - character(STRMAX) :: errmsg - - if (lfirst) then - select case(param%out_stat) - case('APPEND') - open(unit = LUN, file = param%particle_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', err = 667, iomsg = errmsg) - case('NEW', 'UNKNOWN', 'REPLACE') - open(unit = LUN, file = param%particle_out, status = param%out_stat, form = 'UNFORMATTED', err = 667, iomsg = errmsg) - case default - write(*,*) 'Invalid status code',trim(adjustl(param%out_stat)) - call util_exit(FAILURE) - end select + character(STRMAX) :: errmsg - lfirst = .false. - else - open(unit = LUN, file = param%particle_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', err = 667, iomsg = errmsg) - end if - - if (present(lincludecb)) then - if (lincludecb) then - select type(cb => system%cb) - class is (symba_cb) - write(LUN, err = 667, iomsg = errmsg) cb%id - write(LUN, err = 667, iomsg = errmsg) cb%info - end select - end if - end if - - if (present(plidx) .and. (system%pl%nbody > 0)) then - select type(pl => system%pl) - class is (symba_pl) - do i = 1, size(plidx) - write(LUN, err = 667, iomsg = errmsg) pl%id(plidx(i)) - write(LUN, err = 667, iomsg = errmsg) pl%info(plidx(i)) - end do - end select - end if - - if (present(tpidx) .and. (system%tp%nbody > 0)) then - select type(tp => system%tp) - class is (symba_tp) - do i = 1, size(tpidx) - write(LUN, err = 667, iomsg = errmsg) tp%id(tpidx(i)) - write(LUN, err = 667, iomsg = errmsg) tp%info(tpidx(i)) - end do - end select - end if - - close(unit = LUN, err = 667, iomsg = errmsg) + write(iu, err = 667, iomsg = errmsg) self%origin_type + write(iu, err = 667, iomsg = errmsg) self%origin_time + write(iu, err = 667, iomsg = errmsg) self%origin_xh(:) + write(iu, err = 667, iomsg = errmsg) self%origin_vh(:) return 667 continue - write(*,*) "Error reading central body file: " // trim(adjustl(errmsg)) + write(*,*) "Error writing particle metadata information from file: " // trim(adjustl(errmsg)) call util_exit(FAILURE) end subroutine symba_io_dump_particle_info @@ -120,8 +70,6 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms ifirst = ilast + 1 param_value = io_get_token(line_trim, ifirst, ilast, iostat) select case (param_name) - case ("PARTICLE_OUT") - param%particle_out = param_value case ("FRAGMENTATION") call io_toupper(param_value) if (param_value == "YES" .or. param_value == "T") self%lfragmentation = .true. @@ -226,7 +174,6 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms ! Special handling is required for writing the random number seed array as its size is not known until runtime ! For the "SEED" parameter line, the first value will be the size of the seed array and the rest will be the seed array elements - write(param_name, Afmt) "PARTICLE_OUT"; write(param_value, Afmt) trim(adjustl(param%particle_out)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "GMTINY"; write(param_value, Rfmt) param%GMTINY; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "MIN_GMFRAG"; write(param_value, Rfmt) param%min_GMfrag; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "FRAGMENTATION"; write(param_value, Lfmt) param%lfragmentation; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) @@ -256,70 +203,29 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms end subroutine symba_io_param_writer - module subroutine symba_io_read_particle(system, param) + module subroutine symba_io_read_in_particle_info(self, iu) !! author: David A. Minton !! - !! Reads an old particle information file for a restartd run + !! Reads in particle information object information from an open file unformatted file implicit none - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system file - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA extensions - + ! Arguments + class(symba_particle_info), intent(inout) :: self !! Particle metadata information object + integer(I4B), intent(in) :: iu !! Open file unit number ! Internals - integer(I4B), parameter :: LUN = 22 - integer(I4B) :: i, id, idx - logical :: lmatch - type(symba_particle_info) :: tmpinfo - character(STRMAX) :: errmsg - - open(unit = LUN, file = param%particle_out, status = 'OLD', form = 'UNFORMATTED', err = 667, iomsg = errmsg) - - select type(cb => system%cb) - class is (symba_cb) - select type(pl => system%pl) - class is (symba_pl) - select type(tp => system%tp) - class is (symba_tp) - associate(npl => pl%nbody, ntp => tp%nbody) - do - lmatch = .false. - read(LUN, err = 667, iomsg = errmsg, end = 333) id - - if (id == cb%id) then - read(LUN, err = 667, iomsg = errmsg) cb%info - lmatch = .true. - else - if (npl > 0) then - idx = findloc(pl%id(1:npl), id, dim=1) - if (idx /= 0) then - read(LUN, err = 667, iomsg = errmsg) pl%info(idx) - lmatch = .true. - end if - end if - if (.not.lmatch .and. ntp > 0) then - idx = findloc(tp%id(1:ntp), id, dim=1) - if (idx /= 0) then - read(LUN, err = 667, iomsg = errmsg) tp%info(idx) - lmatch = .true. - end if - end if - end if - if (.not.lmatch) then - read(LUN, err = 667, iomsg = errmsg) tmpinfo - end if - end do - end associate - close(unit = LUN, err = 667, iomsg = errmsg) - end select - end select - end select + character(STRMAX) :: errmsg + + call io_read_in_particle_info(self, iu) + read(iu, err = 667, iomsg = errmsg) self%origin_type + read(iu, err = 667, iomsg = errmsg) self%origin_time + read(iu, err = 667, iomsg = errmsg) self%origin_xh(:) + read(iu, err = 667, iomsg = errmsg) self%origin_vh(:) - 333 continue return 667 continue - write(*,*) "Error reading particle information file: " // trim(adjustl(errmsg)) + write(*,*) "Error reading particle metadata information from file: " // trim(adjustl(errmsg)) call util_exit(FAILURE) - end subroutine symba_io_read_particle + end subroutine symba_io_read_in_particle_info module subroutine symba_io_write_discard(self, param) diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90 index 166f6e804..f7cc82024 100644 --- a/src/symba/symba_setup.f90 +++ b/src/symba/symba_setup.f90 @@ -2,60 +2,48 @@ use swiftest contains - module subroutine symba_setup_initialize_particle_info(system, param) + module subroutine symba_setup_initialize_particle_info_system(self, param) !! author: David A. Minton !! !! Initializes a new particle information data structure with initial conditions recorded implicit none ! Argumets - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA extensions + class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i - integer(I4B), dimension(:), allocatable :: idx - - select type(cb => system%cb) - class is (symba_cb) - cb%info%origin_type = "Central body" - cb%info%origin_time = param%t0 - cb%info%origin_xh(:) = 0.0_DP - cb%info%origin_vh(:) = 0.0_DP - call symba_io_dump_particle_info(system, param, lincludecb=.true.) + + select type(cbinfo => self%cb%info) + class is (symba_particle_info) + cbinfo%origin_type = "Initial conditions" + cbinfo%origin_time = param%t0 + cbinfo%origin_xh(:) = 0.0_DP + cbinfo%origin_vh(:) = 0.0_DP end select - select type(pl => system%pl) - class is (symba_pl) - do i = 1, pl%nbody - pl%info(i)%origin_type = "Initial conditions" - pl%info(i)%origin_time = param%t0 - pl%info(i)%origin_xh(:) = pl%xh(:,i) - pl%info(i)%origin_vh(:) = pl%vh(:,i) + select type(plinfo => self%pl%info) + class is (symba_particle_info) + do i = 1, self%pl%nbody + plinfo(i)%origin_type = "Initial conditions" + plinfo(i)%origin_time = param%t0 + plinfo(i)%origin_xh(:) = self%pl%xh(:,i) + plinfo(i)%origin_vh(:) = self%pl%vh(:,i) end do - if (pl%nbody > 0) then - allocate(idx(pl%nbody)) - call symba_io_dump_particle_info(system, param, plidx=[(i, i=1, pl%nbody)]) - deallocate(idx) - end if end select - select type(tp => system%tp) - class is (symba_tp) - do i = 1, tp%nbody - tp%info(i)%origin_type = "Initial conditions" - tp%info(i)%origin_time = param%t0 - tp%info(i)%origin_xh(:) = tp%xh(:,i) - tp%info(i)%origin_vh(:) = tp%vh(:,i) + select type(tpinfo => self%tp%info) + class is (symba_particle_info) + do i = 1, self%tp%nbody + tpinfo(i)%origin_type = "Initial conditions" + tpinfo(i)%origin_time = param%t0 + tpinfo(i)%origin_xh(:) = self%tp%xh(:,i) + tpinfo(i)%origin_vh(:) = self%tp%vh(:,i) end do - if (tp%nbody > 0) then - allocate(idx(tp%nbody)) - call symba_io_dump_particle_info(system, param, tpidx=[(i, i=1, tp%nbody)]) - deallocate(idx) - end if end select - + call setup_initialize_particle_info_system(self, param) return - end subroutine symba_setup_initialize_particle_info + end subroutine symba_setup_initialize_particle_info_system module subroutine symba_setup_initialize_system(self, param) @@ -77,14 +65,6 @@ module subroutine symba_setup_initialize_system(self, param) call system%pltpenc_list%setup(0) call system%plplenc_list%setup(0) call system%plplcollision_list%setup(0) - select type(param) - class is (symba_parameters) - if (param%lrestart) then - call symba_io_read_particle(system, param) - else - call symba_setup_initialize_particle_info(system, param) - end if - end select end associate return @@ -135,6 +115,8 @@ module subroutine symba_setup_pl(self, n, param) call setup_pl(self, n, param) if (n <= 0) return + + if (allocated(self%info)) deallocate(self%info) if (allocated(self%lcollision)) deallocate(self%lcollision) if (allocated(self%lencounter)) deallocate(self%lencounter) if (allocated(self%lmtiny)) deallocate(self%lmtiny) @@ -148,6 +130,7 @@ module subroutine symba_setup_pl(self, n, param) if (allocated(self%kin)) deallocate(self%kin) if (allocated(self%info)) deallocate(self%info) + allocate(symba_particle_info :: self%info(n)) allocate(self%lcollision(n)) allocate(self%lencounter(n)) allocate(self%lmtiny(n)) @@ -159,7 +142,6 @@ module subroutine symba_setup_pl(self, n, param) allocate(self%peri(n)) allocate(self%atp(n)) allocate(self%kin(n)) - allocate(self%info(n)) self%lcollision(:) = .false. self%lencounter(:) = .false. diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index aa2fe0234..f5201ac6c 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -2,33 +2,6 @@ use swiftest contains - module subroutine symba_util_append_arr_info(arr, source, nold, nsrc, lsource_mask) - !! author: David A. Minton - !! - !! Append a single array of particle information type onto another. If the destination array is not allocated, or is not big enough, this will allocate space for it. - implicit none - ! Arguments - type(symba_particle_info), dimension(:), allocatable, intent(inout) :: arr !! Destination array - type(symba_particle_info), dimension(:), allocatable, intent(in) :: source !! Array to append - integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively - logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to - ! Internals - integer(I4B) :: nnew - - if (.not. allocated(source)) return - - nnew = count(lsource_mask(1:nsrc)) - if (.not.allocated(arr)) then - allocate(arr(nold+nnew)) - else - call util_resize(arr, nold + nnew) - end if - - arr(nold + 1:nold + nnew) = pack(source(1:nsrc), lsource_mask(1:nsrc)) - - return - end subroutine symba_util_append_arr_info - module subroutine symba_util_append_arr_kin(arr, source, nold, nsrc, lsource_mask) !! author: David A. Minton @@ -210,26 +183,6 @@ module subroutine symba_util_copy_encounter(self, source) end subroutine symba_util_copy_encounter - module subroutine symba_util_fill_arr_info(keeps, inserts, lfill_list) - !! author: David A. Minton - !! - !! Performs a fill operation on a single array of particle origin information types - !! This is the inverse of a spill operation - implicit none - ! Arguments - type(symba_particle_info), dimension(:), allocatable, intent(inout) :: keeps !! Array of values to keep - type(symba_particle_info), dimension(:), allocatable, intent(in) :: inserts !! Array of values to insert into keep - logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps - - if (.not.allocated(keeps) .or. .not.allocated(inserts)) return - - keeps(:) = unpack(keeps(:), .not.lfill_list(:), keeps(:)) - keeps(:) = unpack(inserts(:), lfill_list(:), keeps(:)) - - return - end subroutine symba_util_fill_arr_info - - module subroutine symba_util_fill_arr_kin(keeps, inserts, lfill_list) !! author: David A. Minton !! @@ -501,7 +454,7 @@ module subroutine symba_util_rearray_pl(self, system, param) lmask(npl+1:npl+nadd) = pl%status(npl+1:npl+nadd) == NEW_PARTICLE npl = pl%nbody - call symba_io_dump_particle_info(system, param, plidx=pack([(i, i=1, npl)], lmask)) + call pl%dump_particle_info(param, idx=pack([(i, i=1, npl)], lmask)) deallocate(lmask) end if @@ -578,40 +531,6 @@ module subroutine symba_util_rearray_pl(self, system, param) end subroutine symba_util_rearray_pl - module subroutine symba_util_resize_arr_info(arr, nnew) - !! author: David A. Minton - !! - !! Resizes an array component of type character string. Array will only be resized if has previously been allocated. Passing nnew = 0 will deallocate. - implicit none - ! Arguments - type(symba_particle_info), dimension(:), allocatable, intent(inout) :: arr !! Array to resize - integer(I4B), intent(in) :: nnew !! New size - ! Internals - type(symba_particle_info), dimension(:), allocatable :: tmp !! Temporary storage array in case the input array is already allocated - integer(I4B) :: nold !! Old size - - if (.not. allocated(arr) .or. nnew < 0) return - - nold = size(arr) - if (nnew == nold) return - - if (nnew == 0) then - deallocate(arr) - return - end if - - allocate(tmp(nnew)) - if (nnew > nold) then - tmp(1:nold) = arr(1:nold) - else - tmp(1:nnew) = arr(1:nnew) - end if - call move_alloc(tmp, arr) - - return - end subroutine symba_util_resize_arr_info - - module subroutine symba_util_resize_arr_kin(arr, nnew) !! author: David A. Minton !! @@ -803,26 +722,6 @@ module subroutine symba_util_sort_tp(self, sortby, ascending) end subroutine symba_util_sort_tp - module subroutine symba_util_sort_rearrange_arr_info(arr, ind, n) - !! author: David A. Minton - !! - !! Rearrange a single array of particle information type in-place from an index list. - implicit none - ! Arguments - type(symba_particle_info), dimension(:), allocatable, intent(inout) :: arr !! Destination array - integer(I4B), dimension(:), intent(in) :: ind !! Index to rearrange against - integer(I4B), intent(in) :: n !! Number of elements in arr and ind to rearrange - ! Internals - type(symba_particle_info), dimension(:), allocatable :: tmp !! Temporary copy of array used during rearrange operation - - if (.not. allocated(arr) .or. n <= 0) return - allocate(tmp, source=arr) - tmp(1:n) = arr(ind(1:n)) - call move_alloc(tmp, arr) - - return - end subroutine symba_util_sort_rearrange_arr_info - module subroutine symba_util_sort_rearrange_arr_kin(arr, ind, n) !! author: David A. Minton @@ -908,45 +807,6 @@ module subroutine symba_util_sort_rearrange_tp(self, ind) end subroutine symba_util_sort_rearrange_tp - module subroutine symba_util_spill_arr_info(keeps, discards, lspill_list, ldestructive) - !! author: David A. Minton - !! - !! Performs a spill operation on a single array of particle origin information types - !! This is the inverse of a spill operation - implicit none - ! Arguments - type(symba_particle_info), dimension(:), allocatable, intent(inout) :: keeps !! Array of values to keep - type(symba_particle_info), dimension(:), allocatable, intent(inout) :: discards !! Array of discards - logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discardss - logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not - ! Internals - integer(I4B) :: nspill, nkeep, nlist - - nkeep = count(.not.lspill_list(:)) - nspill = count(lspill_list(:)) - nlist = size(lspill_list(:)) - - if (.not.allocated(keeps) .or. nspill == 0) return - if (.not.allocated(discards)) then - allocate(discards(nspill)) - else if (size(discards) /= nspill) then - deallocate(discards) - allocate(discards(nspill)) - end if - - discards(:) = pack(keeps(1:nlist), lspill_list(1:nlist)) - if (ldestructive) then - if (nkeep > 0) then - keeps(:) = pack(keeps(1:nlist), .not. lspill_list(1:nlist)) - else - deallocate(keeps) - end if - end if - - return - end subroutine symba_util_spill_arr_info - - module subroutine symba_util_spill_arr_kin(keeps, discards, lspill_list, ldestructive) !! author: David A. Minton !! diff --git a/src/util/util_append.f90 b/src/util/util_append.f90 index 971e3c950..11664827e 100644 --- a/src/util/util_append.f90 +++ b/src/util/util_append.f90 @@ -116,6 +116,34 @@ module subroutine util_append_arr_I4B(arr, source, nold, nsrc, lsource_mask) end subroutine util_append_arr_I4B + module subroutine util_append_arr_info(arr, source, nold, nsrc, lsource_mask) + !! author: David A. Minton + !! + !! Append a single array of particle information type onto another. If the destination array is not allocated, or is not big enough, this will allocate space for it. + implicit none + ! Arguments + class(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: arr !! Destination array + class(swiftest_particle_info), dimension(:), allocatable, intent(in) :: source !! Array to append + integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + ! Internals + integer(I4B) :: nnew + + if (.not. allocated(source)) return + + nnew = count(lsource_mask(1:nsrc)) + if (.not.allocated(arr)) then + allocate(arr(nold+nnew)) + else + call util_resize(arr, nold + nnew) + end if + + arr(nold + 1:nold + nnew) = pack(source(1:nsrc), lsource_mask(1:nsrc)) + + return + end subroutine util_append_arr_info + + module subroutine util_append_arr_logical(arr, source, nold, nsrc, lsource_mask) !! author: David A. Minton !! @@ -156,7 +184,7 @@ module subroutine util_append_body(self, source, lsource_mask) logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to associate(nold => self%nbody, nsrc => source%nbody) - call util_append(self%name, source%name, nold, nsrc, lsource_mask) + call util_append(self%info, source%info, nold, nsrc, lsource_mask) call util_append(self%id, source%id, nold, nsrc, lsource_mask) call util_append(self%status, source%status, nold, nsrc, lsource_mask) call util_append(self%ldiscard, source%ldiscard, nold, nsrc, lsource_mask) diff --git a/src/util/util_fill.f90 b/src/util/util_fill.f90 index 4a5a70311..b1bf951d0 100644 --- a/src/util/util_fill.f90 +++ b/src/util/util_fill.f90 @@ -82,6 +82,27 @@ module subroutine util_fill_arr_I4B(keeps, inserts, lfill_list) return end subroutine util_fill_arr_I4B + + module subroutine util_fill_arr_info(keeps, inserts, lfill_list) + !! author: David A. Minton + !! + !! Performs a fill operation on a single array of particle origin information types + !! This is the inverse of a spill operation + implicit none + ! Arguments + class(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: keeps !! Array of values to keep + class(swiftest_particle_info), dimension(:), allocatable, intent(in) :: inserts !! Array of values to insert into keep + logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps + + if (.not.allocated(keeps) .or. .not.allocated(inserts)) return + + keeps(:) = unpack(keeps(:), .not.lfill_list(:), keeps(:)) + keeps(:) = unpack(inserts(:), lfill_list(:), keeps(:)) + + return + end subroutine util_fill_arr_info + + module subroutine util_fill_arr_logical(keeps, inserts, lfill_list) !! author: David A. Minton !! @@ -119,7 +140,7 @@ module subroutine util_fill_body(self, inserts, lfill_list) !> Fill all the common components associate(keeps => self) call util_fill(keeps%id, inserts%id, lfill_list) - call util_fill(keeps%name, inserts%name, lfill_list) + call util_fill(keeps%info, inserts%info, lfill_list) call util_fill(keeps%status, inserts%status, lfill_list) call util_fill(keeps%ldiscard, inserts%ldiscard, lfill_list) call util_fill(keeps%lmask, inserts%lmask, lfill_list) diff --git a/src/util/util_resize.f90 b/src/util/util_resize.f90 index 1dee5fdb3..e76d53a8d 100644 --- a/src/util/util_resize.f90 +++ b/src/util/util_resize.f90 @@ -138,6 +138,41 @@ module subroutine util_resize_arr_I4B(arr, nnew) end subroutine util_resize_arr_I4B + + module subroutine util_resize_arr_info(arr, nnew) + !! author: David A. Minton + !! + !! Resizes an array component of type character string. Array will only be resized if has previously been allocated. Passing nnew = 0 will deallocate. + implicit none + ! Arguments + class(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: arr !! Array to resize + integer(I4B), intent(in) :: nnew !! New size + ! Internals + class(swiftest_particle_info), dimension(:), allocatable :: tmp !! Temporary storage array in case the input array is already allocated + integer(I4B) :: nold !! Old size + + if (.not. allocated(arr) .or. nnew < 0) return + + nold = size(arr) + if (nnew == nold) return + + if (nnew == 0) then + deallocate(arr) + return + end if + + allocate(tmp(nnew)) + if (nnew > nold) then + tmp(1:nold) = arr(1:nold) + else + tmp(1:nnew) = arr(1:nnew) + end if + call move_alloc(tmp, arr) + + return + end subroutine util_resize_arr_info + + module subroutine util_resize_arr_logical(arr, nnew) !! author: David A. Minton !! @@ -181,7 +216,7 @@ module subroutine util_resize_body(self, nnew) class(swiftest_body), intent(inout) :: self !! Swiftest body object integer(I4B), intent(in) :: nnew !! New size neded - call util_resize(self%name, nnew) + call util_resize(self%info, nnew) call util_resize(self%id, nnew) call util_resize(self%status, nnew) call util_resize(self%ldiscard, nnew) diff --git a/src/util/util_sort.f90 b/src/util/util_sort.f90 index b2a5464aa..4b96dd9d1 100644 --- a/src/util/util_sort.f90 +++ b/src/util/util_sort.f90 @@ -328,7 +328,7 @@ module subroutine util_sort_rearrange_body(self, ind) associate(n => self%nbody) call util_sort_rearrange(self%id, ind, n) - call util_sort_rearrange(self%name, ind, n) + call util_sort_rearrange(self%info, ind, n) call util_sort_rearrange(self%status, ind, n) call util_sort_rearrange(self%ldiscard, ind, n) call util_sort_rearrange(self%xh, ind, n) @@ -459,6 +459,27 @@ module subroutine util_sort_rearrange_arr_logical(arr, ind, n) end subroutine util_sort_rearrange_arr_logical + module subroutine util_sort_rearrange_arr_info(arr, ind, n) + !! author: David A. Minton + !! + !! Rearrange a single array of particle information type in-place from an index list. + implicit none + ! Arguments + class(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: arr !! Destination array + integer(I4B), dimension(:), intent(in) :: ind !! Index to rearrange against + integer(I4B), intent(in) :: n !! Number of elements in arr and ind to rearrange + ! Internals + class(swiftest_particle_info), dimension(:), allocatable :: tmp !! Temporary copy of array used during rearrange operation + + if (.not. allocated(arr) .or. n <= 0) return + allocate(tmp, source=arr) + tmp(1:n) = arr(ind(1:n)) + call move_alloc(tmp, arr) + + return + end subroutine util_sort_rearrange_arr_info + + module subroutine util_sort_rearrange_pl(self, ind) !! author: David A. Minton !! diff --git a/src/util/util_spill.f90 b/src/util/util_spill.f90 index 63cce317a..66e5d22d8 100644 --- a/src/util/util_spill.f90 +++ b/src/util/util_spill.f90 @@ -199,7 +199,46 @@ module subroutine util_spill_arr_I8B(keeps, discards, lspill_list, ldestructive) return end subroutine util_spill_arr_I8B - + + + module subroutine util_spill_arr_info(keeps, discards, lspill_list, ldestructive) + !! author: David A. Minton + !! + !! Performs a spill operation on a single array of particle origin information types + !! This is the inverse of a spill operation + implicit none + ! Arguments + class(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: keeps !! Array of values to keep + class(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: discards !! Array of discards + logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discardss + logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not + ! Internals + integer(I4B) :: nspill, nkeep, nlist + + nkeep = count(.not.lspill_list(:)) + nspill = count(lspill_list(:)) + nlist = size(lspill_list(:)) + + if (.not.allocated(keeps) .or. nspill == 0) return + if (.not.allocated(discards)) then + allocate(discards(nspill)) + else if (size(discards) /= nspill) then + deallocate(discards) + allocate(discards(nspill)) + end if + + discards(:) = pack(keeps(1:nlist), lspill_list(1:nlist)) + if (ldestructive) then + if (nkeep > 0) then + keeps(:) = pack(keeps(1:nlist), .not. lspill_list(1:nlist)) + else + deallocate(keeps) + end if + end if + + return + end subroutine util_spill_arr_info + module subroutine util_spill_arr_logical(keeps, discards, lspill_list, ldestructive) !! author: David A. Minton @@ -257,7 +296,7 @@ module subroutine util_spill_body(self, discards, lspill_list, ldestructive) !> Spill all the common components associate(keeps => self) call util_spill(keeps%id, discards%id, lspill_list, ldestructive) - call util_spill(keeps%name, discards%name, lspill_list, ldestructive) + call util_spill(keeps%info, discards%info, lspill_list, ldestructive) call util_spill(keeps%status, discards%status, lspill_list, ldestructive) call util_spill(keeps%lmask, discards%lmask, lspill_list, ldestructive) call util_spill(keeps%ldiscard, discards%ldiscard, lspill_list, ldestructive)