From 88ba918de0663bd05b739a17686cd35afdb30f2f Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 10 Aug 2021 11:14:24 -0400 Subject: [PATCH] Added energy file to parameter read. Started particle info io methods in SyMBA --- src/io/io.f90 | 18 +++-- src/modules/swiftest_classes.f90 | 2 +- src/modules/symba_classes.f90 | 77 +++++++++----------- src/symba/symba_io.f90 | 116 +++++++++++++++++++++++-------- src/symba/symba_setup.f90 | 2 + src/symba/symba_util.f90 | 5 ++ 6 files changed, 139 insertions(+), 81 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index 52183460c..e19ce2558 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -31,11 +31,11 @@ module subroutine io_conservation_report(self, param, lterminal) Euntracked => self%Euntracked, Eorbit_orig => param%Eorbit_orig, Mtot_orig => param%Mtot_orig, & Ltot_orig => param%Ltot_orig(:), Lmag_orig => param%Lmag_orig, Lorbit_orig => param%Lorbit_orig(:), Lspin_orig => param%Lspin_orig(:), & lfirst => param%lfirstenergy) - if (lfirst) then - if (param%out_stat == "OLD") then - open(unit = EGYIU, file = ENERGY_FILE, form = "formatted", status = "old", action = "write", position = "append") - else - open(unit = EGYIU, file = ENERGY_FILE, form = "formatted", status = "replace", action = "write") + if (param%energy_out /= "") then + if (lfirst .and. (param%out_stat /= "OLD")) then + open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "replace", action = "write") + else + open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "old", action = "write", position = "append") write(EGYIU,EGYHEADER) end if end if @@ -59,8 +59,10 @@ module subroutine io_conservation_report(self, param, lterminal) lfirst = .false. end if - write(EGYIU,EGYFMT) param%t, Eorbit_now, Ecollisions, Ltot_now, Mtot_now - flush(EGYIU) + if (param%energy_out /= "") then + write(EGYIU,EGYFMT) param%t, Eorbit_now, Ecollisions, Ltot_now, Mtot_now + close(EGYIU) + end if if (.not.lfirst .and. lterminal) then Lmag_now = norm2(Ltot_now) Lerror = norm2(Ltot_now - Ltot_orig) / Lmag_orig @@ -422,6 +424,8 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) param%enc_out = param_value case ("DISCARD_OUT") param%discard_out = param_value + case ("ENERGY_OUT") + param%energy_out = param_value case ("EXTRA_FORCE") call io_toupper(param_value) if (param_value == "YES" .or. param_value == 'T') param%lextra_force = .true. diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 25f258d0c..ff32faf80 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -45,7 +45,7 @@ module swiftest_classes real(QP) :: DU2M = -1.0_QP !! Converts distance unit to centimeters real(DP) :: GU = -1.0_DP !! Universal gravitational constant in the system units real(DP) :: inv_c2 = -1.0_DP !! Inverse speed of light squared in the system units - character(STRMAX) :: ennergy_out = "" !! Name of output energy and momentum report file + character(STRMAX) :: energy_out = "" !! Name of output energy and momentum report file ! Logical flags to turn on or off various features of the code logical :: lrhill_present = .false. !! Hill radii are given as an input rather than calculated by the code (can be used to inflate close encounter regions manually) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index b80a9cab8..f11c5d444 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -27,33 +27,17 @@ module symba_classes procedure :: writer => symba_io_param_writer end type symba_parameters - !******************************************************************************************************************************** - ! symba_cb class definitions and method interfaces - !******************************************************************************************************************************* - !> SyMBA central body particle class - type, extends(helio_cb) :: symba_cb - real(DP) :: M0 = 0.0_DP !! Initial mass of the central body - real(DP) :: dM = 0.0_DP !! Change in 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 - contains - end type symba_cb - !******************************************************************************************************************************** ! symba_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, extends(swiftest_base) :: symba_particle_info + 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.) 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 :: dump => symba_io_dump_particle_info !! I/O routine for dumping particle info to file - procedure :: initialize => symba_io_initialize_particle_info !! I/O routine for reading in particle info data - procedure :: read_frame => symba_io_read_frame_info !! I/O routine for reading in a single frame of particle info - procedure :: write_frame => symba_io_write_frame_info !! I/O routine for writing out a single frame of particle info end type symba_particle_info !******************************************************************************************************************************** @@ -66,6 +50,19 @@ module symba_classes integer(I4B), dimension(:), allocatable :: child !! Index of children particles end type symba_kinship + !******************************************************************************************************************************** + ! symba_cb class definitions and method interfaces + !******************************************************************************************************************************* + !> SyMBA central body particle class + type, extends(helio_cb) :: symba_cb + real(DP) :: M0 = 0.0_DP !! Initial mass of the central body + real(DP) :: dM = 0.0_DP !! Change in 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 + !******************************************************************************************************************************** ! symba_pl class definitions and method interfaces !******************************************************************************************************************************* @@ -118,6 +115,7 @@ 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 @@ -327,12 +325,12 @@ 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(self, param, msg) - use swiftest_classes, only : swiftest_parameters + module subroutine symba_io_dump_particle_info(system, param, tpidx, plidx) implicit none - class(symba_particle_info), intent(inout) :: self !! Swiftest base object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - character(*), optional, intent(in) :: msg !! Message to display with dump operation + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA extensions + 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 module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) @@ -357,22 +355,21 @@ 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_initialize_particle_info(self, param) - use swiftest_classes, only : swiftest_parameters + module subroutine symba_io_initialize_particle_info(system, param) implicit none - class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + 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_io_initialize_particle_info - module subroutine symba_io_read_frame_info(self, iu, param, form, ierr) - use swiftest_classes, only : swiftest_parameters - implicit none - class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - character(*), intent(in) :: form !! Input format code ("XV" or "EL") - integer(I4B), intent(out) :: ierr !! Error code - end subroutine symba_io_read_frame_info + !module subroutine symba_io_read_frame_info(self, iu, param, form, ierr) + ! use swiftest_classes, only : swiftest_parameters + ! implicit none + ! class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object + ! integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + ! class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! character(*), intent(in) :: form !! Input format code ("XV" or "EL") + ! integer(I4B), intent(out) :: ierr !! Error code + !end subroutine symba_io_read_frame_info module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters @@ -403,14 +400,6 @@ module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn) integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration end subroutine symba_kick_pltpenc - module subroutine symba_io_write_frame_info(self, iu, param) - use swiftest_classes, only : swiftest_parameters - implicit none - class(symba_particle_info), intent(in) :: self !! SyMBA particle info object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine symba_io_write_frame_info - module subroutine symba_setup_initialize_system(self, param) use swiftest_classes, only : swiftest_parameters implicit none diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index faa3d446b..af7a4f706 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -2,25 +2,91 @@ use swiftest contains - module subroutine symba_io_dump_particle_info(self, param, msg) + module subroutine symba_io_dump_particle_info(system, param, tpidx, plidx) !! author: David A. Minton !! - !! Dumps the particle information data to a file + !! 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 - class(symba_particle_info), intent(inout) :: self !! Swiftest base object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - character(*), optional, intent(in) :: msg !! Message to display with dump operation + ! 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 + 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) :: 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) + case('NEW', 'UNKNOWN', 'REPLACE') + open(unit = iu, file = param%particle_file, status = param%out_stat, form = 'UNFORMATTED', iostat = ierr) + case default + write(*,*) 'Invalid status code',trim(adjustl(param%out_stat)) + call util_exit(FAILURE) + end select + if (ierr /= 0) then + write(*, *) "Swiftest error:" + write(*, *) " particle output file already exists or cannot be accessed" + call util_exit(FAILURE) + end if + + lfirst = .false. + else + open(unit = iu, file = param%particle_file, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr) + if (ierr /= 0) then + write(*, *) "Swiftest error:" + write(*, *) " unable to open binary output file for APPEND" + call util_exit(FAILURE) + 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)) + end do + end select + end if + + if (present(tpidx) .and. (system%tp%nbody > 0) .and. size(tpidx) > 0) then + 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)) + end do + end select + end if + + close(unit = iu, iostat = ierr) + if (ierr /= 0) then + write(*, *) "Swiftest error:" + write(*, *) " unable to close particle output file" + call util_exit(FAILURE) + end if + + return end subroutine symba_io_dump_particle_info - module subroutine symba_io_initialize_particle_info(self, param) + module subroutine symba_io_initialize_particle_info(system, param) !! author: David A. Minton !! !! Initializes a particle info data structure, either starting a new one or reading one in !! from a file if it is a restarted run implicit none - class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! 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 + + return end subroutine symba_io_initialize_particle_info @@ -194,19 +260,19 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms end subroutine symba_io_param_writer - module subroutine symba_io_read_frame_info(self, iu, param, form, ierr) - !! author: David A. Minton - !! - !! Reads a single frame of a particle info data from a file. - implicit none - class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - character(*), intent(in) :: form !! Input format code ("XV" or "EL") - integer(I4B), intent(out) :: ierr !! Error code - - ierr = 0 - end subroutine symba_io_read_frame_info + !module subroutine symba_io_read_frame_info(self, iu, param, form, ierr) + ! !! author: David A. Minton + ! !! + ! !! Reads a single frame of a particle info data from a file. + ! implicit none + ! class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object + ! integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + ! class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! character(*), intent(in) :: form !! Input format code ("XV" or "EL") + ! integer(I4B), intent(out) :: ierr !! Error code +! +! ierr = 0 +! end subroutine symba_io_read_frame_info module subroutine symba_io_write_discard(self, param) @@ -280,13 +346,5 @@ module subroutine symba_io_write_discard(self, param) return end subroutine symba_io_write_discard - - module subroutine symba_io_write_frame_info(self, iu, param) - implicit none - class(symba_particle_info), intent(in) :: self !! SyMBA particle info object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine symba_io_write_frame_info - end submodule s_symba_io diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90 index ab8b5543e..8d727be5c 100644 --- a/src/symba/symba_setup.f90 +++ b/src/symba/symba_setup.f90 @@ -162,10 +162,12 @@ module subroutine symba_setup_tp(self, n, param) if (allocated(self%nplenc)) deallocate(self%nplenc) if (allocated(self%levelg)) deallocate(self%levelg) if (allocated(self%levelm)) deallocate(self%levelm) + if (allocated(self%info)) deallocate(self%info) allocate(self%nplenc(n)) allocate(self%levelg(n)) allocate(self%levelm(n)) + allocate(self%info(n)) self%nplenc(:) = 0 self%levelg(:) = -1 diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 90f5a06e5..efb1832a1 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -144,6 +144,7 @@ module subroutine symba_util_append_tp(self, source, lsource_mask) call util_append(self%nplenc, source%nplenc, nold, nsrc, lsource_mask) call util_append(self%levelg, source%levelg, nold, nsrc, lsource_mask) call util_append(self%levelm, source%levelm, nold, nsrc, lsource_mask) + call util_append(self%info, source%info, nold, nsrc, lsource_mask) call util_append_tp(self, source, lsource_mask) ! Note: helio_tp does not have its own append method, so we skip back to the base class end associate @@ -253,6 +254,7 @@ module subroutine symba_util_fill_tp(self, inserts, lfill_list) call util_fill(keeps%nplenc, inserts%nplenc, lfill_list) call util_fill(keeps%levelg, inserts%levelg, lfill_list) call util_fill(keeps%levelm, inserts%levelm, lfill_list) + call util_fill(keeps%info, inserts%info, lfill_list) call util_fill_tp(keeps, inserts, lfill_list) ! Note: helio_tp does not have its own fill method, so we skip back to the base class class default @@ -520,6 +522,7 @@ module subroutine symba_util_resize_tp(self, nnew) call util_resize(self%nplenc, nnew) call util_resize(self%levelg, nnew) call util_resize(self%levelm, nnew) + call util_resize(self%info, nnew) call util_resize_tp(self, nnew) @@ -679,6 +682,7 @@ module subroutine symba_util_sort_rearrange_tp(self, ind) if (allocated(tp%nplenc)) tp%nplenc(1:ntp) = tp_sorted%nplenc(ind(1:ntp)) if (allocated(tp%levelg)) tp%levelg(1:ntp) = tp_sorted%levelg(ind(1:ntp)) if (allocated(tp%levelm)) tp%levelm(1:ntp) = tp_sorted%levelm(ind(1:ntp)) + if (allocated(tp%info)) tp%info(1:ntp) = tp_sorted%info(ind(1:ntp)) deallocate(tp_sorted) end associate @@ -836,6 +840,7 @@ module subroutine symba_util_spill_tp(self, discards, lspill_list, ldestructive) call util_spill(keeps%nplenc, discards%nplenc, lspill_list, ldestructive) call util_spill(keeps%levelg, discards%levelg, lspill_list, ldestructive) call util_spill(keeps%levelm, discards%levelm, lspill_list, ldestructive) + call util_spill(keeps%info, discards%info, lspill_list, ldestructive) call util_spill_tp(keeps, discards, lspill_list, ldestructive) class default