From 58dac0d2a278511529047bbc462433408a0886b2 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 10 Aug 2021 12:33:36 -0400 Subject: [PATCH] Refactored PARTICLE_FILE to PARTICLE_OUT for consistency. Added in particle info io to SyMBA --- .../param.disruption_headon.in | 2 +- .../param.disruption_off_axis.in | 2 +- .../symba_energy_momentum/param.escape.in | 2 +- examples/symba_energy_momentum/param.sun.in | 2 +- .../param.supercatastrophic_headon.in | 2 +- .../param.supercatastrophic_off_axis.in | 2 +- src/io/io.f90 | 2 +- src/modules/symba_classes.f90 | 7 ++-- src/symba/symba_io.f90 | 40 +++++++++++++------ src/symba/symba_setup.f90 | 15 ++++++- src/symba/symba_util.f90 | 5 ++- 11 files changed, 56 insertions(+), 25 deletions(-) diff --git a/examples/symba_energy_momentum/param.disruption_headon.in b/examples/symba_energy_momentum/param.disruption_headon.in index 0f3e88752..4a291535e 100644 --- a/examples/symba_energy_momentum/param.disruption_headon.in +++ b/examples/symba_energy_momentum/param.disruption_headon.in @@ -7,7 +7,7 @@ TP_IN tp.in IN_TYPE ASCII ISTEP_OUT 1 ! output cadence every year BIN_OUT bin.disruption_headon.dat -PARTICLE_FILE particle.disruption_headon.dat +PARTICLE_OUT particle.disruption_headon.dat OUT_TYPE REAL8 ! double precision real output OUT_FORM XV ! osculating element output OUT_STAT REPLACE diff --git a/examples/symba_energy_momentum/param.disruption_off_axis.in b/examples/symba_energy_momentum/param.disruption_off_axis.in index ef32a5c2f..0dfbae80a 100644 --- a/examples/symba_energy_momentum/param.disruption_off_axis.in +++ b/examples/symba_energy_momentum/param.disruption_off_axis.in @@ -7,7 +7,7 @@ TP_IN tp.in IN_TYPE ASCII ISTEP_OUT 1 ! output cadence every year BIN_OUT bin.disruption_off_axis.dat -PARTICLE_FILE particle.disruption_off_axis.dat +PARTICLE_OUT particle.disruption_off_axis.dat OUT_TYPE REAL8 ! double precision real output OUT_FORM XV ! osculating element output OUT_STAT REPLACE diff --git a/examples/symba_energy_momentum/param.escape.in b/examples/symba_energy_momentum/param.escape.in index 5db2c3fe4..90d118017 100644 --- a/examples/symba_energy_momentum/param.escape.in +++ b/examples/symba_energy_momentum/param.escape.in @@ -7,7 +7,7 @@ TP_IN tp.in IN_TYPE ASCII ISTEP_OUT 1 ! output cadence every year BIN_OUT bin.escape.dat -PARTICLE_FILE particle.escape.dat +PARTICLE_OUT particle.escape.dat OUT_TYPE REAL8 ! double precision real output OUT_FORM XV ! osculating element output OUT_STAT REPLACE diff --git a/examples/symba_energy_momentum/param.sun.in b/examples/symba_energy_momentum/param.sun.in index a21b5817b..a7748b19c 100644 --- a/examples/symba_energy_momentum/param.sun.in +++ b/examples/symba_energy_momentum/param.sun.in @@ -9,7 +9,7 @@ IN_TYPE ASCII ISTEP_OUT 1 ISTEP_DUMP 1 BIN_OUT bin.sun.dat -PARTICLE_FILE particle.sun.dat +PARTICLE_OUT particle.sun.dat OUT_TYPE REAL8 OUT_FORM XV ! osculating element output OUT_STAT REPLACE diff --git a/examples/symba_energy_momentum/param.supercatastrophic_headon.in b/examples/symba_energy_momentum/param.supercatastrophic_headon.in index 47c239556..e9b60e7da 100644 --- a/examples/symba_energy_momentum/param.supercatastrophic_headon.in +++ b/examples/symba_energy_momentum/param.supercatastrophic_headon.in @@ -7,7 +7,7 @@ TP_IN tp.in IN_TYPE ASCII ISTEP_OUT 1 ! output cadence every year BIN_OUT bin.supercatastrophic_headon.dat -PARTICLE_FILE particle.supercatastrophic_headon.dat +PARTICLE_OUT particle.supercatastrophic_headon.dat OUT_TYPE REAL8 ! double precision real output OUT_FORM XV ! osculating element output OUT_STAT REPLACE diff --git a/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in b/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in index 64759828c..0bf836be5 100644 --- a/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in +++ b/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in @@ -7,7 +7,7 @@ TP_IN tp.in IN_TYPE ASCII ISTEP_OUT 1 ! output cadence every year BIN_OUT bin.supercatastrophic_off_axis.dat -PARTICLE_FILE particle.supercatastrophic_off_axis.dat +PARTICLE_OUT particle.supercatastrophic_off_axis.dat OUT_TYPE REAL8 ! double precision real output OUT_FORM XV ! osculating element output OUT_STAT REPLACE diff --git a/src/io/io.f90 b/src/io/io.f90 index e2908b957..03fdc2e17 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -498,7 +498,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) read(param_value, *) param%Ecollisions case("EUNTRACKED") read(param_value, *) param%Euntracked - case ("NPLMAX", "NTPMAX", "GMTINY", "PARTICLE_FILE", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters + case ("NPLMAX", "NTPMAX", "GMTINY", "PARTICLE_OUT", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters case default write(iomsg,*) "Unknown parameter -> ",param_name iostat = -1 diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index de747708c..4628202f8 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -18,7 +18,7 @@ module symba_classes integer(I4B), parameter :: PARTICLEUNIT = 44 !! File unit number for the binary particle info output file type, extends(swiftest_parameters) :: symba_parameters - character(STRMAX) :: particle_file = PARTICLE_OUTFILE !! Name of output particle information file + character(STRMAX) :: particle_out = PARTICLE_OUTFILE !! Name of output particle information file real(DP) :: GMTINY = -1.0_DP !! Smallest mass that is fully gravitating integer(I4B), dimension(:), allocatable :: seed !! Random seeds logical :: lfragmentation = .false. !! Do fragmentation modeling instead of simple merger. @@ -325,10 +325,11 @@ module subroutine symba_io_write_discard(self, param) 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, tpidx, plidx) + module subroutine symba_io_dump_particle_info(system, param, lincludecb, tpidx, plidx) 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 + 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 end subroutine symba_io_dump_particle_info diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 906929bef..7751b7b21 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -2,7 +2,7 @@ use swiftest contains - module subroutine symba_io_dump_particle_info(system, param, tpidx, plidx) + module subroutine symba_io_dump_particle_info(system, param, lincludecb, tpidx, plidx) !! author: David A. Minton !! !! Dumps the particle information data to a file. @@ -10,21 +10,22 @@ module subroutine symba_io_dump_particle_info(system, param, tpidx, plidx) implicit none ! Arguments 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_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 ! Internals logical, save :: lfirst = .true. - integer(I4B), parameter :: iu = 22 + integer(I4B), parameter :: LUN = 22 integer(I4B) :: i, ierr if (.not.present(tpidx) .and. .not.present(plidx)) return if (lfirst) then select case(param%out_stat) case('APPEND') - open(unit = iu, file = param%particle_file, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr) + open(unit = LUN, file = param%particle_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr) case('NEW', 'UNKNOWN', 'REPLACE') - open(unit = iu, file = param%particle_file, status = param%out_stat, form = 'UNFORMATTED', iostat = ierr) + open(unit = LUN, file = param%particle_out, status = param%out_stat, form = 'UNFORMATTED', iostat = ierr) case default write(*,*) 'Invalid status code',trim(adjustl(param%out_stat)) call util_exit(FAILURE) @@ -37,7 +38,7 @@ module subroutine symba_io_dump_particle_info(system, param, tpidx, plidx) lfirst = .false. else - open(unit = iu, file = param%particle_file, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr) + open(unit = LUN, file = param%particle_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr) if (ierr /= 0) then write(*, *) "Swiftest error:" write(*, *) " unable to open binary output file for APPEND" @@ -45,12 +46,22 @@ module subroutine symba_io_dump_particle_info(system, param, tpidx, plidx) end if end if + if (present(lincludecb)) then + if (lincludecb) then + select type(cb => system%cb) + class is (symba_cb) + write(LUN) cb%id + write(LUN) cb%info + end select + end if + end if + if (present(plidx) .and. (system%pl%nbody > 0) .and. size(plidx) > 0) then select type(pl => system%pl) class is (symba_pl) do i = 1, size(plidx) - write(iu) pl%id(plidx(i)) - write(iu) pl%info(plidx(i)) + write(LUN) pl%id(plidx(i)) + write(LUN) pl%info(plidx(i)) end do end select end if @@ -59,13 +70,13 @@ module subroutine symba_io_dump_particle_info(system, param, tpidx, plidx) select type(tp => system%tp) class is (symba_tp) do i = 1, size(tpidx) - write(iu) tp%id(tpidx(i)) - write(iu) tp%info(tpidx(i)) + write(LUN) tp%id(tpidx(i)) + write(LUN) tp%info(tpidx(i)) end do end select end if - close(unit = iu, iostat = ierr) + close(unit = LUN, iostat = ierr) if (ierr /= 0) then write(*, *) "Swiftest error:" write(*, *) " unable to close particle output file" @@ -75,6 +86,7 @@ module subroutine symba_io_dump_particle_info(system, param, tpidx, plidx) return end subroutine symba_io_dump_particle_info + module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! @@ -119,6 +131,8 @@ 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. @@ -217,7 +231,7 @@ 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_FILE"; write(param_value, Afmt) trim(adjustl(param%particle_file)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "PARTICLE_OUT"; write(param_value, Afmt) trim(adjustl(param%particle_out)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "GMTINY"; write(param_value, Rfmt) param%Gmtiny; write(unit, Afmt) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "FRAGMENTATION"; write(param_value, Lfmt) param%lfragmentation; write(unit, Afmt) adjustl(param_name), adjustl(param_value) if (param%lfragmentation) then @@ -259,7 +273,7 @@ module subroutine symba_io_read_particle(system, param) logical :: lmatch type(symba_particle_info) :: tmpinfo - open(unit = LUN, file = param%particle_file, status = 'OLD', form = 'UNFORMATTED', iostat = ierr) + open(unit = LUN, file = param%particle_out, status = 'OLD', form = 'UNFORMATTED', iostat = ierr) if (ierr /= 0) then write(*, *) "Swiftest error:" write(*, *) " unable to open binary particle file for reading" diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90 index a5870ed52..e06fb20b5 100644 --- a/src/symba/symba_setup.f90 +++ b/src/symba/symba_setup.f90 @@ -12,6 +12,7 @@ module subroutine symba_setup_initialize_particle_info(system, param) class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA extensions ! Internals integer(I4B) :: i + integer(I4B), dimension(:), allocatable :: idx select type(cb => system%cb) class is (symba_cb) @@ -19,6 +20,7 @@ module subroutine symba_setup_initialize_particle_info(system, param) 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.) end select select type(pl => system%pl) @@ -29,6 +31,11 @@ module subroutine symba_setup_initialize_particle_info(system, param) pl%info(i)%origin_xh(:) = pl%xh(:,i) pl%info(i)%origin_vh(:) = 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) @@ -38,9 +45,15 @@ module subroutine symba_setup_initialize_particle_info(system, param) tp%info(i)%origin_time = param%t0 tp%info(i)%origin_xh(:) = tp%xh(:,i) tp%info(i)%origin_vh(:) = tp%vh(:,i) - end do + 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 + return end subroutine symba_setup_initialize_particle_info diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index efb1832a1..4c0f256e3 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -381,7 +381,10 @@ module subroutine symba_util_rearray_pl(self, system, param) if (allocated(pl%xend)) deallocate(pl%xend) ! Add in any new bodies - call pl%append(pl_adds, lsource_mask=[(.true., i=1, pl_adds%nbody)]) + if (pl_adds%nbody > 0) then + call pl%append(pl_adds, lsource_mask=[(.true., i=1, pl_adds%nbody)]) + call symba_io_dump_particle_info(system, param, plidx=[(i, i = 1, pl%nbody)]) + end if ! If there are still bodies in the system, sort by mass in descending order and re-index if (pl%nbody > 0) then